多维缩放(MDS)算法的详细推导及Python实现

  最近对流形学习比较感兴趣，因为流形学习用到了MDS算法，所以写一篇博客记录下来，主要是对西瓜书里面的公式进行详细的推导，并给出Python代码实现。MDS算法即多维缩放(Multiple Dimensional Scaling)算法，是一种常见的降维算法，核心思想是要求原始空间样本间的距离在低维空间得到保持。
  现给定训练样本 χ = { x 1 , x 2 , . . . , x m } \chi=\{x_1, x_2,...,x_m\} χ={x1​,x2​,...,xm​}，其中$x_i\in R^d,i=1,...,m$，共$m$个样本，每个样本的长度为$d$。我们的目标是获得在低维$d^{\prime }$空间的表示 Z = { z 1 , z 2 , . . . , z m } ∈ R d ′ × m Z=\{z_1, z_2,...,z_m\}\in R^{d'\times{m}} Z={z1​,z2​,...,zm​}∈Rd′×m，且$d^{\prime }<d$。记$m$个样本在原始空间的距离矩阵为 D = { d i j } ∈ R m × m D=\{d_{ij}\}\in R^{m\times{m}} D={dij​}∈Rm×m，其中$d_{ij}$为样本$x_i$与$x_j$的距离，且$d_{ij}=||x_i-x_j|{|}_2^2$，我们的目标是使得样本在低维空间的距离等于其在高维空间的距离，即$d_{ij}=||x_i-x_j|{|}_2^2=||z_i-z_j|{|}_2^2$。
  令$B=Z^{T}Z\in R^{m\times m}$，$B$为降维后的内积矩阵，即$b_{ij}=z_i^T{z}_j$，则$$d_{ij}=||z_i-z_j|{|}_2^2=(z_i{)}^2+(z_j{)}^2-2z_i^T{z}_j=b_{ii}+b_{jj}-2b_{ij}$$此时令降维后的样本$Z$被中心化${\sum }_{i=1}^m{z}_i=0$，这里的中心化要说明一下，第一遍看的时候迷迷糊糊，跑完实验才彻底理解，中心化指的是所有的低维向量的某一个维度之和等于零，而非低维向量自身的所有维度相加等于零。那么$$\begin{array}{rl}\sum _{i=1}^m{b}_{ij}& =b_{1j}+b_{2j}+...+b_{mj}\\ & =z_1^T{z}_j+z_2^T{z}_j+...+z_m^T{z}_j\\ & =(z_1^T+z_2^T+...+z_m^T)z_j\\ & =0\end{array}$$同理，${\sum }_{j=1}^m{b}_{ij}=0$。接下来可以推导得到：$$\begin{array}{rl}\sum _{i=1}^m{d}_{ij}^2& =d_{1j}^2+d_{2j}^2+...+d_{mj}^2\\ & =b_{11}+b_{jj}-2b_{1j}+b_{22}+b_{jj}-2b_{2j}+...+b_{mm}+b_{jj}-2b_{mj}\\ & =b_{11}+b_{22}+...+b_{mm}+m{b}_{jj}-2(b_{1j}+b_{2j}+...+b_{mj})\\ & =tr(B)+m{b}_{jj}\end{array}$$同理，${\sum }_{j=1}^m{d}_{ij}^2=tr(B)+m{b}_{ii}$，最后：$$\begin{array}{rl}\sum _{i=1}^m\sum _{j=1}^m{d}_{ij}^2& =\sum _{j=1}^m{d}_{1j}^2+\sum _{j=1}^m{d}_{2j}^2+...+\sum _{j=1}^m{d}_{mj}^2\\ & =tr(B)+m{b}_{11}+tr(B)+m{b}_{22}+...+tr(B)+m{b}_{mm}\\ & =m\times tr(B)+m(b_{11}+b_{22}+...+b_{mm})\\ & =2m\times tr(B)\end{array}$$在得到原始空间样本的距离矩阵$D$和${\sum }_{i=1}^m{d}_{ij}^2$、${\sum }_{j=1}^m{d}_{ij}^2$以及${\sum }_{i=1}^m{\sum }_{j=1}^m{d}_{ij}^2$之后，就可以求得降维后的内积矩阵$B$：$$\begin{array}{rl}{b}_{ij}& =\frac{1}{2}(b_{ii}+b_{jj}-d_{ij}^2)\\ & =\frac{1}{2}(\frac{1}{m}[tr(B)+m{b}_{jj}]+\frac{1}{m}[tr(B)+m{b}_{ii}]-\frac{2}{m}tr(B)-d_{ij}^2)\\ & =\frac{1}{2}(\frac{1}{m}\sum _{i=1}^m{d}_{ij}^2+\frac{1}{m}\sum _{j=1}^m{d}_{ij}^2-\frac{1}{m^2}\sum _{i=1}^m\sum _{j=1}^m{d}_{ij}^2-d_{ij}^2)\end{array}$$  下一步就是对内积矩阵$B$进行特征值分解：$B=V\Lambda V^T$，其中$\Lambda =diag({\lambda }_1,{\lambda }_2,...,{\lambda }_d)$特征值构成的对角矩阵,且${\lambda }_1\ge {\lambda }_2\ge ...\ge {\lambda }_d$。现取前$d^{\prime }$个特征值构成对角矩阵${\Lambda }^{\prime }=diag({\lambda }_1,{\lambda }_2,...,{\lambda }_{d^{\prime }})\in R^{d^{\prime }\times d^{\prime }}$，对应的特征向量矩阵$V^{\prime }\in R^{m\times d^{\prime }}$。因为$B=Z^{T}Z$，所以$Z={\Lambda }^{\prime 1/2}{V}^{\prime T}\in R^{d^{\prime }\times m}$，这就得到了最后的低维嵌入向量$Z$。公式一大堆，看似很复杂，最后总结一下就三个步骤：首先求原始空间的距离矩阵$D$，然后根据$D$求降维后的内积矩阵$B$，最后根据$B$求低维嵌入$Z$。
  下面通过代码来说明一下，这里只给出部分核心代码，完整代码见我的GitHub，代码部分参考了博客MDS算法，在此表示感谢。首先是距离矩阵$D$的求解：

def get_distance_matrix(data):
	expand_ = data[:, np.newaxis, :]
	repeat1 = np.repeat(expand_, data.shape[0], axis=1)
	repeat2 = np.swapaxes(repeat1, 0, 1)
	D = np.linalg.norm(repeat1 - repeat2, ord=2, axis=-1, keepdims=True).squeeze(-1)
	return D

然后是求解内积矩阵$B$：

def get_matrix_B(D):
	assert D.shape[0] == D.shape[1]
	DD = np.square(D)
	sum_ = np.sum(DD, axis=1) / D.shape[0]
	Di = np.repeat(sum_[:, np.newaxis], D.shape[0], axis=1)
	Dj = np.repeat(sum_[np.newaxis, :], D.shape[0], axis=0)
	Dij = np.sum(DD) / ((D.shape[0])**2) * np.ones([D.shape[0], D.shape[0]])
	B = (Di + Dj - DD- Dij) / 2
	return B

最后根据$B$求解$Z$：

def MDS(data, n=2):
	D = get_distance_matrix(data)
	B = get_matrix_B(D)
	B_value, B_vector = np.linalg.eigh(B)
	Be_sort = np.argsort(-B_value)
	B_value = B_value[Be_sort]               # 降序排列的特征值
	B_vector = B_vector[:,Be_sort]           # 降序排列的特征值对应的特征向量
	Bez = np.diag(B_value[0:n])
	Bvz = B_vector[:, 0:n]
	Z = np.dot(np.sqrt(Bez), Bvz.T).T
	return Z

最后贴一张实验结果图，将大家都爱的瑞士卷(swiss roll)降到二维，本博客到此结束啦！