多维尺度变换

转载链接：点击打开链接

多维尺度变换简介

         图形化数据的一个重要目标就是了解数据点相互之间离得多远或多近。通常，你可以使用散点图来更明显的查看。然而，在一些分析中，你使用的数据可能不是以点的形式给出的，而是观察数据之间的相似性或不相似性，这时将无点可绘制。

         即使你的数据是以点的形式给出，散点图有时未必是有用的。比如对一些数据，衡量两个点多近的方式可能并不是Euclidean距离，然而原始数据的散点图在比较Euclidean距离时是比较直观的，但比较其他类型的距离时并不总是有效，如街区距离等。而且，对于变量比较多的情况也是很难把距离直观的展现出来的，除非数据只有很小的维数，所以降维常常是必要的。

         多维尺度变换MDS就是一类解决上面问题的方法，它能直观的展现数据点之间在多种距离上离得多近，还可以产生一个低维的数据点表示。MDS不需要原始数据，只需要衡量点间距离的距离矩阵。

经典多维尺度变换

经典多维尺度变换简介

         函数cmdscale可以用来实现经典的多维尺度变换，也称主坐标分析。它以内点距离矩阵为输入，输出一系列点。理想情况下，这些点是二维或三维的，它们之间的Euclidean距离产生与原来相同的距离矩阵，因此这些点的散点图可以直观的呈现原始的距离。

         作为一个非常简单的例子，你可以从距离重构原始点集。首先创建一组四维点集，其中的第四维上是较小的成分，然后计算点之间的距离：

X = [ normrnd(0,1,10,3), normrnd(0,.1,10,1) ];

D = pdist(X,'euclidean');

         接着，使用cmdscale产生相同距离的点集，它接受距离方阵或距离上三角向量为输入，上面pdist的结果D就是距离上三角向量：

[Y,eigvals] = cmdscale(D);

        函数cmdscale产生两个输出，第一个Y是包含重构点的矩阵，第二个eigvals是被称为标量积矩阵（最简单的形式就是Y*Y’）的顺序特征值向量。这些特征值的相对大小表示了使用这些重构点再现原始距离矩阵D时，Y的相应各列的相对贡献。

format short g

[eigvals eigvals/max(abs(eigvals))]

ans =

       12.623            1

       4.3699      0.34618

       1.9307      0.15295

     0.025884    0.0020505

  1.7192e-015  1.3619e-016

  6.8727e-016  5.4445e-017

  4.4367e-017  3.5147e-018

 -9.2731e-016 -7.3461e-017

  -1.327e-015 -1.0513e-016

 -1.9232e-015 -1.5236e-016

         如果eigvals只包含正或0（包括舍入误差在内）特征值，那么与正特征值对应的Y的各列可以在Euclidean距离意义下对D进行准确重构，如使用pdist函数计算的距离将和D近似相等

maxerr4 = max(abs(D - pdist(Y))) % exact reconstruction

maxerr4 =

  2.6645e-015

         如果eigvals中2或3个特征值比其他都大很多，那么基于Y相应列产生的距离矩阵将与D几乎一样。在这样意义下，对应的列将形成一个低维的近似表示，然而并不是总能找到一个好的低维重构。

% good reconstruction in 3D

maxerr3 = max(abs(D - pdist(Y(:,1:3))))

maxerr3 =

     0.029728

 

% poor reconstruction in 2D

maxerr2 = max(abs(D - pdist(Y(:,1:2))))

maxerr2 =

      0.91641

         可以看到三维重构产生的D很好，但二维重构的误差与D中的最大值等量级：

max(max(D))

ans =

       3.4686

         有时，eigvals中包含一些负的特征值，这表明矩阵D中的距离不能准确重构，即无法配置点，使得它们的Euclidean距离矩阵是给定的D矩阵。但如果绝对值最大的负特征值相对于最大正特征值在幅值上很小，那么cmdscale函数返回的Y还是可以重构D的。

多维尺度变换的例子

         假设只给定10个US城市间的距离，cmdscale可以重建这些城市间的地图。首先，创建距离矩阵，并把他传给cmdscale。在这个例子中，D是距离矩阵：它是方阵、对称阵，对角线上元素为0，其他元素为正。

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);

         然后看cmdscale返回的特征值，一些特征值是负的，说明原始的距离不是Euclidean距离，这是因为地球的表面不是平面。

         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

         然而，在这种情况下，最大的两个正特征值在幅度上比其他特征值都要大很多，所以尽管有负特征值的存在，Y的前两维坐标还是足以可靠的重构D

Dtriu = D(find(tril(ones(10),-1)))';

maxrelerr = max(abs(Dtriu-pdist(Y(:,1:2))))./max(Dtriu)

maxrelerr =

    0.0075371

         下面是重构的城市位置地图，重构的原点是任意的。可以看到虽然结果不是完全准确的，但与真实结果很接近。

非经典多维尺度变换

         函数mdscale可以实现非经典多维尺度变换。和cmdscale相似，你也可以使用mdscale来直观展现数据差异或通过降维来呈现高维数据。两个函数都是以差异性矩阵（距离矩阵）作为输入，输出一系列点配置。然而mdscale提供一些额外的重构数据的准则供选择，而且允许有遗失数据存在、权重存在。

         例如，cereal数据包含10个描述早餐谷类食物的变量的测量值，这里使用mdscale函数在二维直观展示这些数据。首先，加载数据，这里为了清晰，选择观察数据中的22个来做分析

load cereal.mat

X = [Calories Protein Fat Sodium Fiber ...

    Carbo Sugars Shelf Potass Vitamins];

% Take a subset from a single manufacturer

mfg1 = strcmp('G',cellstr(Mfg));

X = X(mfg1,:);

size(X)

ans =

    22 10

         然后使用pdist函数将这10维数据转化为差异性描述。Pdist的输出是一个对称差异矩阵，以向量形式存储在为一个包含23*22/2个元素的向量里

dissimilarities = pdist(zscore(X),'cityblock');

size(dissimilarities)

ans =

     1   231

         这个例子代码中首先对数据做标准化，然后使用街区距离来衡量差异性。差异性的衡量方式是依赖于具体应用的，这里为了简单使用街区距离。在一些应用中，原始数据就是以差异性的形式给出。

         然后，使用mdscale函数执行非经典多维尺度变换。和cmdscale不同，这里你需要指定期望维数和重建输出配置的方法。如对这个例子，使用两维和常用的STRESS准则输出计算方法，更多方法参加mdscale函数界面。Mdscale的第二个输出是一个评价输出配置的值，它衡量了输出配置的距离与原始输入差异性的吻合程度：

[Y,stress] =...

mdscale(dissimilarities,2,'criterion','metricstress');

stress

stress =

    0.1856

         下面的散点图在二维中绘制了原十维数据，你可以使用gname函数给选择的点做标记：

plot(Y(:,1),Y(:,2),'o','LineWidth',2);

gname(Name(mfg1))

非度量多维尺度变换

         上面的度量多维尺度变换根据差异矩阵创建了一系列点，这些点的距离矩阵与原始距离矩阵近似，这个要求有时是比较严格的，非度量多维尺度变换就是为了放松这一要求。与前面近似差异矩阵不同，非度量多维尺度变换是近似一个非线性但单调的变换。由于单调性，输出图中，更大或更小的距离与更大或更小的差异对应。然而，非线性暗示了mdscale只是试图保存差异性的次序。因此，在不同的尺度可能有距离的压缩或展开。

         你可以使用mdscale函数以度量MDS相同的方式执行非度量MDS。非度量STRESS准则是一种常用的方法，更多选型参见mdscale函数帮助页面。度量MDS中，mdscale的第二个输出是一个评价输出效果的值。而对于非度量MDS，它衡量了输出配置的距离近似不一致性的吻合程度。不一致性由第三个参数返回，它们是原始差异性的变换值。

[Y,stress,disparities] = ...

mdscale(dissimilarities,2,'criterion','stress');

stress

stress =

    0.1562

        为了检验输出配置对差异性拟合的效果，理解不一致性，可以绘制Shepard图

distances = pdist(Y);

[dum,ord] = sortrows([disparities(:) dissimilarities(:)]);

plot(dissimilarities,distances,'bo', ...

     dissimilarities(ord),disparities(ord),'r.-', ...

     [0 25],[0 25],'k-')

xlabel('Dissimilarities')

ylabel('Distances/Disparities')

legend({'Distances' 'Disparities' '1:1 Line'},...

       'Location','NorthWest');

         上图可以看出，mdscale在二维找到了一组点的配置，它们的点间距离近似了不一致性，而此不一致性正是原始差异性的非线性变换。不一致性对差异性的凹函数曲线表示拟合相对于差异性有小的收缩趋向，这在实际中是确实显而易见的。

         Midscale函数使用迭代算法来找到输出配置，所以结果常常依赖于初始点。默认情况下，mdscale以cmdscale重建的配置作为初始值，这样常常可以得到一个全局最优解。然而，mdscale仍有可能收敛到一个局部极小处。这种情况可以通过使用不同初始点多次执行mdscale来诊断和解决。你可以使用额外参数’start’和’replicates’参数来控制。下面的代码执行了5次MDS，每次从随机的初始配置开始。每次执行的准则值显示如下，最终mdscale返回最优拟合的结果：

opts = statset('Display','final');

[Y,stress] =...

mdscale(dissimilarities,2,'criterion','stress',...

'start','random','replicates',5,'Options',opts);

 

35 iterations, Final stress criterion = 0.156209

31 iterations, Final stress criterion = 0.156209

48 iterations, Final stress criterion = 0.171209

33 iterations, Final stress criterion = 0.175341

32 iterations, Final stress criterion = 0.185881

         注意到mdscale找到了几个不同的局部解，有些解没有以cmdscale为初始点的解低。