（十一）MDS算法

1.用MDS的场景

 从降维的层面来说，由于MDS是一种降维方法，那么它和PCA等其他降维方式有什么不同呢，什么样的场景适用于使用MDS呢？与PCA不同的是，MDS保证了原始数据点之间的距离与降维后数据点的距离一致。
 另一方面，假如我们只知道一组点之间的距离，我们如何反演出它的相对坐标？这也是MDS可以做到的。
 MDS的这种降维是继承了原始数据空间中的欧式距离度量。使得数据点能在低维空间重构其相对位置。
 举个例子，比如如果用欧式距离度量两种商品的相似度，那么，商品特征繁多，我们就需要在不破坏商品之间的这种相似度的情况下对数据进行降维。
 再比如我们手头只有一个邻接矩阵，代表飞机坠机后的残骸之间的距离，现在新来了一个残骸，多个声呐探测后能够计算出这个残骸与其他所有残骸的距离，但是我们无法得知其具体位置，那么我们可以通过MDS反演出其相对距离，那么我们只要知道其中一个残骸的位置，这个残骸的位置也就很明朗了。根据这个相对距离缩放到真实距离就行。

2.MDS算法的数学推导

定理：在p维欧式空间中有一组点集，这组点集之间的欧式距离由矩阵D给出，即\forall \vec{x_{r}},\vec{x_{s}}\in X_{[n\times p]}\in R^{p};D=[d^{2}_{rs}]_{n\times n};d^{2}_{rs}=||\vec{x_{r}}-\vec{x_{s}}||_{2}^{2}=(\vec{x_{r}}-\vec{x_{s}})^{T}(\vec{x_{r}}-\vec{x_{s}})令A=[a_{rs}]_{n\times n}且a_{rs}=-\frac{1}{2}d^{2}_{rs}定义一个B=[b_{rs}]_{n\times n}=HAH;H=I_{n}-\frac{1}{n}\vec{1}_{n}\vec{1}_{n}^{T};
在如上定义之下，当且仅当B是半正定矩阵，D是欧式距离矩阵。

 我们现在给出构造性证明：
 首先是必要性：
即：

如果D是点集Z=(z_{1},z_{2},...,z_{n})^{T}的欧氏距离矩阵；那么：
b_{rs}=(z_{r}-\bar{z})^{T}(z_{s}-\bar{z});\bar{z}=\vec{1}_{n}\cdot Z=\frac{1}{n}\sum_{i=1}^{n}z_{i};B=HZZ^{T}H;

 证明：
 B=HAH=(I_{n}-\frac{1}{n}\vec{1}_{n}\vec{1}_{n}^{T})A(I_{n}-\frac{1}{n}\vec{1}_{n}\vec{1}_{n}^{T})\\=A+\frac{1}{n^{2}}\vec{1}_{n}\vec{1}_{n}^{T}A\vec{1}_{n}\vec{1}_{n}^{T}-\frac{1}{n}\vec{1}_{n}\vec{1}_{n}^{T}A-\frac{1}{n}A\vec{1}_{n}\vec{1}_{n}^{T}
 矩阵A的列平均:\bar a_{s}=\frac {1}{n}\vec{1}_{n}^{T}A=\frac {1}{n} \sum^{n}_{s=1}a_{rs}
 矩阵A的行平均:\bar a_{r}=\frac{1}{n}A\vec{1}_{n}=\frac {1}{n} \sum^{n}_{r=1}a_{rs}
 矩阵A的整体平均:\bar a_{r}=\frac{1}{n^2}\vec{1}_{n}^{T} A\vec{1}_{n}=\frac {1}{n^2}\sum^{n}_{s=1} \sum^{n}_{r=1}a_{rs}
 那么：
b_{rs}=a_{rs}-\bar{a_{s}}-\bar{a_{r}}+\bar a\\=-\frac {1}{2}d^{2}_{rs}+\frac {1}{2n}\sum^{n}_{s=1}d^{2}_{rs}+\frac {1}{2n}\sum^{n}_{r=1}d^{2}_{rs}-\frac {1}{n^2}\sum^{n}_{s=1} \sum^{n}_{r=1}d^{2}_{rs}\\ =-\frac{1}{2} (\vec{z_{r}}-\vec{z_{s}})^{T}(\vec{z_{r}}-\vec{z_{s}})+\frac{1}{2n} \sum^{n}_{s=1}(\vec{z_{r}}-\vec{z_{s}})^{T}(\vec{z_{r}}-\vec{z_{s}})+ \frac{1}{2n} \sum^{n}_{r=1}(\vec{z_{r}}-\vec{z_{s}})^{T}(\vec{z_{r}}-\vec{z_{s}})+\frac{1}{n^2}\sum^{n}_{s=1} \sum^{n}_{r=1}(\vec{z_{r}}-\vec{z_{s}})^{T}(\vec{z_{r}}-\vec{z_{s}})\\=-\frac {1}{2}z_{r}^{T}z_{r}-\frac {1}{2}z_{s}^{T}z_{s}+z_{r}^{T}z_{s}+\frac {1}{2}z_{r}^{T}z_{r}+\frac{1}{2n} \sum^{n}_{j=1}z_{j}^{T}z_{j}-z_{r}^{T}\bar z+ \frac {1}{2}z_{s}^{T}z_{s}+\frac{1}{2n} \sum^{n}_{j=1}z_{j}^{T}z_{j}-z_{s}^{T}\bar z-\frac{1}{2n} \sum^{n}_{j=1}z_{j}^{T}z_{j}-\frac{1}{2n} \sum^{n}_{i=1}z_{i}^{T}z_{i}+\bar z^{T}\bar z\\=z_{r}^{T}z_{s}-z_{r}^{T}\bar z-z_{s}^{T}\bar z+\bar z^{T}\bar z\\=(z_{r}-\bar{z})^{T}(z_{s}-\bar{z})
 所以：
b_{rs}=(z_{r}-\bar{z})^{T}(z_{s}-\bar{z})
B=HZZ^{T}H
由此可得，B一定是一个对称半正定矩阵，必要性证毕。

 接着证明充分性：
即：

如果B是秩为p的对称半正定矩阵，那么可以将B按照如下方式构造：
令\lambda_{1}>\lambda_{2}>...>\lambda_{p}表示B矩阵的所有大于0的特征值，且对应的特征向量为：Z=(Z_{(1)},Z_{(2)},...,Z_{(p)})且满足Z_{(i)}^{T}Z_{(i)}=\lambda_{i};i=1,2,...,p;
那么在p维欧式空间中,点集Z之间的距离可由距离矩阵D给出。

证明：
假设特征值对应的矩阵为：\Lambda
令U=Z\Lambda^{-\frac{1}{2}};因为Z^{T}Z=\Lambda
那么U^{T}U=\Lambda^{-\frac{1}{2}}Z^{T}Z\Lambda^{-\frac{1}{2}}=\Lambda^{-\frac{1}{2}}\Lambda\Lambda^{-\frac{1}{2}}=I_{p};
所以可以将B特征分解为：B=U\Lambda U^{T}=Z\Lambda^{-\frac{1}{2}}\Lambda\Lambda^{-\frac{1}{2}}Z^{T}=Z Z^{T}
即：b_{rs}=z_{r}^{T}z_{s};
我们算一下第r个特征向量和第s个特征向量的欧式距离，以确定其确实可以由距离矩阵给出。
(z_{r}-z_{s})^{T}(z_{r}-z_{s})=z_{r}^{T}z_{r}-2z_{r}^{T}z_{s}+z_{s}^{T}z_{s}\\=b_{rr}+b_{ss}-2b_{rs}
因为b_{rs}=a_{rs}-\bar a_{r}-\bar a_{s}+\bar a
那么b_{ss}=a_{ss}-\bar a_{s}-\bar a_{s}+\bar a
那么b_{rr}=a_{rr}-\bar a_{r}-\bar a_{r}+\bar a
且：a_{rr}=-\frac{1}{2}d_{rr}^{2}=0=a_{ss}
所以：
(z_{r}-z_{s})^{T}(z_{r}-z_{s})=z_{r}^{T}z_{r}-2z_{r}^{T}z_{s}+z_{s}^{T}z_{s}\\=b_{rr}+b_{ss}-2b_{rs}\\=a_{ss}-\bar a_{s}-\bar a_{s}+\bar a+a_{rr}-\bar a_{r}-\bar a_{r}+\bar a-2a_{rs}+2\bar a_{r}+2\bar a_{s}-2\bar a\\=-2a_{rs}\\=d_{rs}^{2}
这恰好证明了那么在p维欧式空间中,点集Z之间的距离可由距离矩阵D给出。充分性证毕。

那么也就是说只要B是半正定的矩阵，那么我们总能找到它的特征向量U=Z\Lambda^{-\frac{1}{2}}，将对应的坐标Z=U\Lambda^{\frac{1}{2}}求出。
我们还能发现一些好玩的性质，比如因为H\vec{1_{n}}=0所以：B\vec{1_{n}}=HA(H\vec{1_{n}})=0
由于不同特征值对应的特征向量是正交的，那么这个式子就告诉我们\vec{1_{n}}是其中一个特征向量。且对应特征值为0。
也就是说：
\vec{1_{n}}Z=0\to\sum_{r=1}^{n}Z_{r}=0
这表明实际上，我们在保持其欧式距离结构的前提下，将坐标原点定在了所有点的均值处。
 回顾一下我们做了什么，我们有一个距离矩阵D，这个距离矩阵的来源可以是我们手上有一堆已有数据产生的，或者是直接给出的；我们试图在低维或者同维度下重构它的数据点，我们首先通过D矩阵找到A矩阵，A矩阵构成B=HAH矩阵，由于距离矩阵一定是对称半正定的，B矩阵一定也是对称半正定的。由于B是对称的，所以我们一定可以对它进行特征分解:B=U\Lambda U^{T};如果此时需要在低维下重构，那么对特征值排序后剔除小的，如果在同维度下重构，那么保持原样即可。我们最后就直接还原出其相对坐标：Z=U\Lambda^{\frac{1}{2}}。

3.MDS中加入新数据点的问题

 我们要是新来一个数据点，或者说新来一点，我只是知道它与其他所有点的距离，这样的情况怎么办。
 我们可能会不假思索的回答说，这个很简单啊，只要把点加入矩阵D中，由n\times n变成(n+1)\times (n+1)就好了，其他重头再来一遍得到一个新的坐标点集Z就好了。
 在数据量不是很大的情况下，这确实不失为一种好方法，但是万一数据量很大，重头再来的代价就让人有点头疼了，这就是为什么我们要看看有什么办法可以绕过这样重头再来的过程。
假设我们新的数据点是Z_{n+1}=(z^{(1)}_{n+1},z^{(2)}_{n+1},...,z^{(p)}_{n+1})
我们得计算一下它和已有的第i个点之间的距离：
d^{2}_{(i,n+1)}=\sum_{k=1}^{p}(z^{(k)}_{n+1}-z^{(k)}_{i})^{2}\\ =\sum_{k=1}^{p}(z^{(k)}_{n+1})^{2}+\sum_{k=1}^{p}(z^{(k)}_{i})^{2}-2\sum_{k=1}^{p}z^{(k)}_{n+1}z^{(k)}_{i}\\ =d^{2}_{n+1}+d^{2}_{i}-2\sum_{k=1}^{p}z^{(k)}_{n+1}z^{(k)}_{i}
这里面我们什么是已知的，什么是未知的呢？
d^{2}_{n+1}和-2\sum_{k=1}^{p}z^{(k)}_{n+1}z^{(k)}_{i};是我们未知的，其他都是已知的。
我们对d^{2}_{(i,n+1)}求和，看看是否能将一个未知量给表示成另一个未知量。
\sum_{i=1}^{n}d^{2}_{(i,n+1)}=\sum_{i=1}^{n}d^{2}_{n+1}+\sum_{i=1}^{n}d^{2}_{i}-2\sum_{k=1}^{p}z^{(k)}_{n+1}\sum_{i=1}^{n}z^{(k)}_{i}
由于上面我们有说过：\sum_{r=1}^{n}Z_{r}=0
所以：\sum_{i=1}^{n}z^{(k)}_{i}=0
因此：
\sum_{i=1}^{n}d^{2}_{(i,n+1)}=\sum_{i=1}^{n}d^{2}_{n+1}+\sum_{i=1}^{n}d^{2}_{i}=nd^{2}_{n+1}+\sum_{i=1}^{n}d^{2}_{i}
d^{2}_{n+1}=\frac{1}{n}\cdot \sum_{i=1}^{n}(d^{2}_{(i,n+1)}-d^{2}_{i})
将这个结果代入原式：
d^{2}_{(i,n+1)}=\frac{1}{n}\cdot \sum_{i=1}^{n}(d^{2}_{i}) +d^{2}_{i}-2\sum_{k=1}^{p}z^{(k)}_{n+1}z^{(k)}_{i}
2\sum_{k=1}^{p}z^{(k)}_{n+1}z^{(k)}_{i}=d^{2}_{i}-d^{2}_{(i,n+1)}-\frac{1}{n}\cdot \sum_{i=1}^{n}(d^{2}_{i}-d^{2}_{(i,n+1)})
我们令\alpha_{i}=d^{2}_{i}-d^{2}_{(i,n+1)}
2\sum_{k=1}^{p}z^{(k)}_{n+1}z^{(k)}_{i}= \alpha_{i}- \frac{1}{n}\cdot \sum_{i=1}^{n}\alpha_{i}
写成矩阵形式：
Z\cdot \vec{z}_{n+1}=\frac{1}{2}\cdot (\vec{\alpha}-\frac{1}{n}\vec{1_{n}}\vec{1_{n}}^{T}\vec{\alpha })

Z^{T}Z\cdot \vec{z}_{n+1}=\frac{1}{2}\cdot Z^{T}(\vec{\alpha}-\frac{1}{n}\vec{1_{n}}\vec{1_{n}}^{T}\vec{\alpha })

\Lambda\cdot \vec{z}_{n+1}=\frac{1}{2}\cdot Z^{T}(\vec{\alpha}-\frac{1}{n}\vec{1_{n}}\vec{1_{n}}^{T}\vec{\alpha })
因为Z^{T}\vec{1_{n}}=0
所以新加入的数据点可以由下面公式导出：
\vec{z}_{n+1}=\frac{1}{2}\cdot \Lambda^{-1} Z^{T}\vec{\alpha }

4.代码

import numpy as np
class MDS:
    def __init__(self,Distance_Matrix,q):
        self.D=Distance_Matrix
        self.q=q
        self.n=np.shape(self.D)[0]
        self.A=(-1/2)*self.D
        self.H=np.identity(self.n)-(1/self.n)*np.ones([self.n,self.n])
        self.B=self.H.dot(self.A).dot(self.H)
    def rebuild_coordinates(self):
        self.lambda_,self.U=np.linalg.eig(self.B)
        Largest_q_EigenValues=self.lambda_[np.where(self.lambda_>0)]
        Largest_q_EigenVector=self.U[np.where(self.lambda_>0)]
        Largest_q_EigenValues=Largest_q_EigenValues[np.where(np.argsort(Largest_q_EigenValues)<self.q)]
        Largest_q_EigenVector=Largest_q_EigenVector[np.where(np.argsort(Largest_q_EigenValues)<self.q)]
        self.Lambda=np.diag(Largest_q_EigenValues.real)
        self.U=Largest_q_EigenVector.T
        self.Z=self.U.dot(self.Lambda**(1/2)).real
        return self.Z
    def add_new(self,distance_to_everypoint):
        Z=self.rebuild_coordinates()
        d_i_square=np.array([np.sum(Z**2,axis=1)]).T
        alpha=d_i_square-distance_to_everypoint.T
        new_z=(1/2)*np.linalg.inv(self.Lambda).dot(Z.T).dot(alpha)
        return new_z
if __name__ == '__main__':
    from sklearn.datasets import load_iris
    import seaborn as sns
    import matplotlib.pyplot as plt
    def calculate_distance_matrix(X):
        D = []
        for every_point in X:
            dis = np.sum((X - np.array([every_point])) ** 2, axis=1)
            D.append(dis)
        D = np.array(D)
        return D
    def plot(dfData):
        plt.subplots(figsize=(9, 9))
        sns.heatmap(dfData, annot=True, vmax=1, square=True, cmap="Blues")
        plt.show()
    data=load_iris()
    X=np.array([[1,2],[3,4],[2,2],[0,1],[1,3],[5,4]])
    new_point=np.array([[1,4]])
    new_distance=np.array([np.sum((X - new_point) ** 2, axis=1)])
    D=calculate_distance_matrix(X)
    mds=MDS(D,2)
    Z=mds.rebuild_coordinates()
    new_z=mds.add_new(new_distance)
    print(Z)
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 211,290评论 6赞 491
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,107评论 2赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 156,872评论 0赞 347
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,415评论 1赞 283
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,453评论 6赞 385
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,784评论 1赞 290
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,927评论 3赞 406
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,691评论 0赞 266
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,137评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,472评论 2赞 326
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,622评论 1赞 340
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,289评论 4赞 329
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,887评论 3赞 312
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,741评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,977评论 1赞 265
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,316评论 2赞 360
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,490评论 2赞 348