Multidimensional Scaling -- Reconstructing a Map from Intercity DistancesGiven only the distances between 10 US cities, cmdscale can construct a map of those cities. First, create the distance matrix and pass it to cmdscale. In this example, D is a full distance matrix: it is square and symmetric, has positive entries off the diagonal, and has zeros on the diagonal. cities =
{'Atl','Chi','Den','Hou','LA','Mia','NYC','SF','Sea','WDC'};
D = [ 0 587 1212 701 1936 604 748 2139 2182 543;
587 0 920 940 1745 1188 713 1858 1737 597;
1212 920 0 879 831 1726 1631 949 1021 1494;
701 940 879 0 1374 968 1420 1645 1891 1220;
1936 1745 831 1374 0 2339 2451 347 959 2300;
604 1188 1726 968 2339 0 1092 2594 2734 923;
748 713 1631 1420 2451 1092 0 2571 2408 205;
2139 1858 949 1645 347 2594 2571 0 678 2442;
2182 1737 1021 1891 959 2734 2408 678 0 2329;
543 597 1494 1220 2300 923 205 2442 2329 0];
[Y,eigvals] = cmdscale(D);
Next, look at the eigenvalues returned by cmdscale. Some of these are negative, indicating that the original distances are not Euclidean. This is because of the curvature of the earth. format short g
[eigvals eigvals/max(abs(eigvals))]
ans =
9.5821e+006 1
1.6868e+006 0.17604
8157.3 0.0008513
1432.9 0.00014954
508.67 5.3085e-005
25.143 2.624e-006
5.3394e-010 5.5722e-017
-897.7 -9.3685e-005
-5467.6 -0.0005706
-35479 -0.0037026
However, in this case, the two largest positive eigenvalues are much larger in magnitude than the remaining eigenvalues. So, despite the negative eigenvalues, the first two coordinates of Y are sufficient for a reasonable reproduction of D. Dtriu = D(find(tril(ones(10),-1)))';
maxrelerr = max(abs(Dtriu - pdist(Y(:,1:2)))) ./ max(Dtriu)
maxrelerr =
0.0075371
Here is a plot of the reconstructed city locations as a map. The orientation of the reconstruction is arbitrary: in this case, it happens to be close to, although not exactly, the correct orientation. plot(Y(:,1),Y(:,2),'.');
text(Y(:,1)+25,Y(:,2),cities);
xlabel('Miles'); ylabel('Miles)
|