多为缩放（MDS）

多维缩放（Multiple Dimensional Scaling，MDS）

1 绪论

假设任意测试样本𝑥附近任意小的𝛿距离范围内总能找到一个训练样本，即训练样本的采集密度足够大，或称为“密采样”（dense sample）。

然而，这个假设在现实任务中通常很难满足，例如𝛿=0.001，仅考虑单个属性，则仅需1000个样本点平均分布在归一化后的属性取值范围内，即可使得任意测试样本在其附近0.001距离范围内总能找到一个训练样本，此时最近邻分类器（1NN）的错误率不超过贝叶斯最优分类器的错误率的两倍。然而，这仅是属性维数为1的情形，若有更多的属性，则情况会发生显著变化。例如，假定属性维数为20，若要求样本满足密采样条件，则至少需$(10^3{)}^{20}=10^{60}$个样本。现实应用中属性维数经常成千上万，要满足密采样条件所需的样本数目是无法达到的天文数字。此外，许多学习方法都涉及距离计算，而高维空间会给距离计算带来很大的麻烦。例如，当维数很高时，甚至连计算内积都不再容易。

在高维情形下出现的数据样本稀疏、距离计算困难等问题，是所有机器学习方法共同面临的严重障碍，被称为“维数灾难”（curse of dimensionality）。

以一个简单的例子说明，我们考虑设计一个用于区分猫狗照片的分类器，我们提取了毛发以及形状作为特征，但是仅仅两个特征似乎并不能很好的区分猫和狗，那么考虑增加更多特征效果是否更好呢？

如图，我们嵌入到三维空间中就可以使用一个超平面将猫狗进行很好的分类。若将这个超平面映射到二维平面上，即如第三幅图所示。

如此看来，增加更多的特征效果似乎是越来越好了。

事实上，从下图中可以看出，分类器的性能随着特征个数的变化不断增加，但过了某一个值后，性能不升反降。这种情况被称为“过拟合”，即在训练集上拟合很好，但缺乏足够的泛化能力。“过拟合”通常也是“维数灾难”的直接体现。

那么，是否有办法将高维空间中的数据映射到一个低维“子空间”（subspace）中，且数据还能保持高维空间中的特性，可以线性区分呢？

“降维”（dimension reduction），是缓解维数灾难的一个重要途径，亦称“维数约简”。它一般通过某种数学变换将原始高维属性空间转变为一个低维“子空间” ，在这个子空间中样本密度大幅提高，距离计算也变得更容易。

2 多维缩放

“多维缩放”（Multiple Dimensional Scaling，MDS）最初的许多理论都是由数学心理学家发展起来的，他们在《心理测量学》杂志上发表了他们的许多发现。尽管MDS起源于行为科学，但现在已经变得越来越流行，并被广泛应用于各种应用领域。这种受欢迎程度反映在许多基于计算机的统计数据包中。

多维缩放，即MDS是对一组对象之间的距离或差异的视觉表示。 “对象”可以是颜色，面孔，地图坐标，各种真实存在的事物。除了将不同之处解释为图上的距离之外，MDS还可以作为高维数据的降维技术。简而言之，MDS的主要目的是将这些差异保持在降维中。
MDS算法思想很简单，一句话就是保持样本在原空间和低维空间的距离不变。

3 MDS算法解析

假设𝑁个样本在原始空间中的距离矩阵为$D=(d_{ij}{)}_{N\times N}$:
$$D=\left[\begin{array}{cccc}{d}_{1,1}& d_{1,2}& \cdots & d_{1,n}\\ d_{2,1}& d_{2,2}& \cdots & d_{2,n}\\ ⋮& ⋮& \ddots & ⋮\\ d_{n,1}& d_{n,2}& \cdots & d_{n,n}\end{array}\right]$$
其中$d_{ij}=\Vert {\vec{x}}_i-{\vec{x}}_j\Vert$为样本${\vec{x}}_i$到样本$𝑥_𝑗 $的距离。我们的目标是获得样本在$𝑛^{′$维空间的表示$𝑍∈ℝ}{𝑁×𝑛^′}$， $n^{\prime }\le n$，且任意两个样本在$n^{\prime }$维空间中的欧式距离等于原始空间中的距离，即$\Vert z_i-z_j\Vert =d_{ij}$。

令$B=Z{Z}^T\in R^{N\times N}$，即：
$$B=\left[\begin{array}{cccc}{b}_{1,1}& b_{1,2}& \cdots & b_{1,N}\\ b_{2,1}& b_{2,2}& \cdots & b_{2,N}\\ ⋮& ⋮& \ddots & ⋮\\ b_{N,1}& b_{N,2}& \cdots & b_{N,N}\end{array}\right]$$
其中$b_{ij}=z_i\cdot z_j$为降维后样本的内积。
则根据降维后样本的欧氏距离保持不变有：
$d_{ij}^2=\Vert z_i-z_j{\Vert }^2=\Vert z_i{\Vert }^2+\Vert z_i{\Vert }^2-2z_i^T{z}_j=b_{ii}+b_{jj}-2b_{ij}$

假设降维后的样本集𝑍被中心化，即${\sum }_{i=1}^N{z}_i=0$，则矩阵𝐵的每行之和均为零，每列之和均为零。即：
$$\sum _{i=1}^N{b}_{ij}=0,j=1,2,\dots ,N$$
$$\sum _{j=1}^N{b}_{ij}=0,i=1,2,\dots ,N$$
于是有：
$$\sum _{i=1}^N{d}_{ij}^2=\sum _{i=1}^N{b}_{ii}+N{b}_{jj}=tr(B)+N{b}_{jj}$$
$$\sum _{j=1}^N{d}_{ij}^2=\sum _{j=1}^N{b}_{jj}+N{b}_{ii}=tr(B)+N{b}_{ii}$$
$$\sum _{i=1}^N\sum _{j=1}^N{d}_{ij}^2=\sum _{i=1}^{N}tr(B)+N{b}_{ii}=2Ntr(B)$$

令：
$$d_{i\ast }^2=\frac{1}{N}\sum _{j=1}^N{d}_{ij}^2=\frac{tr(B)}{N}+b_{ii}$$
$$d_{j\ast }^2=\frac{1}{N}\sum _{i=1}^N{d}_{ij}^2=\frac{tr(B)}{N}+b_{jj}$$
$$d_{\ast \ast }^2=\frac{1}{N^2}\sum _{i=1}^N\sum _{j=1}^N{d}_{ij}^2=\frac{2tr(B)}{N}$$
代入$d_{ij}^2=b_{ii}+b_{jj}-2b_{ij}$，有：
$b_{ij}=\frac{b_{ii}+b_{jj}-d_{ij}^2}{2}=\frac{d_{i\ast }^2+d_{j\ast }^2-d_{\ast \ast }^2-d_{ij}^2}{2}$
右式根据$d_{ij}$给出了$b_{ij}$，因此可以根据原始空间中的距离矩阵𝐷求出在降维后空间的内积矩阵𝐵。

我们得到了内积矩阵𝐵，那么如何求解矩阵𝑍呢？

对矩阵𝐵做特征值分解，设$B=V\Lambda V^T$，其中$\Lambda =diag({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_N)$为特征值构成的对角矩阵，${\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_N$，𝑉为特征向量矩阵。
假定特征值中有$n^{\ast }$个非零特征值，它们构成对角矩阵${\Lambda }^{\ast }=diag({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{n^{\ast }})$。令$V^{\ast }$为对应的特征向量矩阵，则$Z=V^{\ast }{\Lambda }^{\ast 1/2}\in R^{N\times n^{\ast }}，n^{\prime }=n^{\ast }$。

在现实应用中为了有效降维，往往仅需要降维后的距离与原始空间中的距离尽可能相等，而不必严格相等。 此时可取$n^{\prime }\ll n^{\ast }\le n$个最大特征值构成对角矩阵$\tilde{\Lambda }=diag({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{(}{n}^{\prime }))$。令$\tilde{V}$表示对应的特征向量矩阵，则$Z=\tilde{V}{\tilde{\Lambda }}^{1/2}\in R^{N\times n^{\prime }}$。

4 MDS伪码说明

应该是Markdown表格有问题，这里就直接使用图片来展示。

5 以瑞士卷例子来展示MDS算法

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as p3d
from sklearn import datasets, manifold
from sklearn.datasets import make_swiss_roll


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

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

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

def test_iris():
	iris = datasets.load_iris()
	data = iris.data                         # [150, 4]
	target = iris.target                     # [150]
	# print(np.shape(data))
	Z = MDS(data)
	figure1 = plt.figure()
	plt.subplot(1, 2, 1)
	plt.plot(Z[target==0, 0], Z[target==0, 1], 'r*', markersize=20)
	plt.plot(Z[target==1, 0], Z[target==1, 1], 'bo', markersize=20)
	plt.plot(Z[target==2, 0], Z[target==2, 1], 'gx', markersize=20)
	plt.title('CUSTOM')
	plt.subplot(1, 2, 2)
	Z1 = manifold.MDS(n_components=2).fit_transform(data)
	plt.plot(Z1[target==0,0], Z1[target==0,1], 'r*', markersize=20)
	plt.plot(Z1[target==1,0], Z1[target==1,1], 'bo', markersize=20)
	plt.plot(Z1[target==2,0], Z1[target==2,1], 'gx', markersize=20)
	plt.title('SKLEARN')
	plt.show()

def test_swiss_roll(n_samples=5000):
	X, t = make_swiss_roll(n_samples=n_samples, noise=0.2, random_state=42) # X为坐标 t为颜色
	Z = MDS(X)
	# figure1=plt.figure()
	axes = [-11.5, 14, -2, 23, -12, 15]	
	fig = plt.figure()
	ax = fig.add_subplot(131, projection='3d')
	ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=t, cmap=plt.cm.hot)
	ax.view_init(10, 60)
	ax.set_xlabel("$x$", fontsize=18)
	ax.set_ylabel("$y$", fontsize=18)
	ax.set_zlabel("$z$", fontsize=18)
	ax.set_xlim(axes[0:2])
	ax.set_ylim(axes[2:4])
	ax.set_zlim(axes[4:6])
	plt.title('3D swiss roll')
	ax2 = fig.add_subplot(132)
	ax2.scatter(Z[:, 0], Z[:, 1], c=t, cmap=plt.cm.hot)
	ax2.set_xlabel("$x$", fontsize=18)
	ax2.set_ylabel("$y$", fontsize=18)
	plt.title('after MDS')
	ax3 = fig.add_subplot(133)
	Z = manifold.MDS(n_components=2).fit_transform(X)
	ax3.scatter(Z[:, 0], Z[:, 1], c=t, cmap=plt.cm.hot)
	plt.title('sklearn MDS')
	plt.show()

if __name__ == '__main__':
	test_iris()
	test_swiss_roll()

程序运行可见：https://aistudio.baidu.com/aistudio/projectdetail/2507328