MDS（multidimensional scaling）多维尺度分析

1. MDS

Multidimensional scaling，简称MDS，中文翻译成多维尺度分析（暂时这样翻译吧。。。）。其原理是利用成对样本间的相似性，去构建合适的低维空间，使得样本在此空间的距离和在高维空间中的样本间的相似性尽可能的保持一致。我们可以用这种方式来可视化数据分布。

2. MDS基本原理

首先，给定 $n$ 个 $d$ 维的样本数据，

$$\mathbf{X} = \left[ \begin{matrix} x_{11} & x_{12} & \cdots & x_{1d} \\ x_{21} & x_{22} & \cdots & x_{2d} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \cdots & x_{nd} \\ \end{matrix} \right]\tag{1}$$

其中这 $n$ 个样本之间的欧式距离可以表示成如下的矩阵

$$d(\mathbf{X})= \left[ \begin{matrix} \delta_{11} &\delta_{12} &\cdots &\delta_{1n} \\ \delta_{21} &\delta_{22} &\cdots &\delta_{2n} \\ \vdots &\vdots &\ddots &\vdots \\ \delta_{n1} &\delta_{n2} &\cdots &\delta_{nn} \\ \end{matrix} \right] \tag{2}$$
其中， $\mathbf{\sigma}_{ij}\mathbf{(X)}=\delta_{ij}$ 是样本 $i$ 个 样本 $j$ 之间的欧式距离，有

$$\delta_{ij}=||\mathbf{x}_{i} - \mathbf{x}_{j}|| \\ = \sqrt{\sum_{s=1}^{d}(x_{is}-x_{js})^{2}} \tag{3}$$
MDS的基本原理是在一个低维的空间中，将 $n$个物体嵌入到该空间中，使得彼此间的相似度尽可能地保留。如果这个低维的空间是2维或者3维，则可以画出该$n$个物体的可视化结果。假设将这 $n$个样本映射到新的$p$维空间中，映射成的矩阵形式如下
$$\mathbf{Y} = \left[ \begin{matrix} y_{11} & y_{12} & \cdots & y_{1p} \\ y_{21} & y_{22} & \cdots & y_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ y_{n1} & y_{n2} & \cdots & y_{np} \\ \end{matrix} \right]\tag{4}$$
其中 $p<<d$。一般 $p$取$2$或者 $3$。通常，计量多元尺度法（Metric MDS）是最小化如下的目标函数
$$\sigma(\mathbf{X}) =\sum_{i=1}^{n}{\sum_{j=1}^{n}{(\delta_{ij}-||\mathbf{y}_{i} - \mathbf{y}_{j}||)^{2}}} \tag{5}$$
对于欧式距离，可以被任意旋转和变换，因为这些变换不会改变样本间的距离。所以，对于 $\mathbf{Y}$，解是不唯一的。

3. Basic majorization theory (基本的最优化理论)

上述的MDS的目标函数可以用如下的优化理论来优化。其基本原理是如下描述。
首先，假设待优化的函数 $f(x)$ 有最优的最小解，但是直接求解 $f(x)$ 的解析解是非常困难的。因此，我们可以找一个简单的，更容易求解最优值的函数 $g_{y}(x)$ 来代替求解 $f(x)$。 其中，$g_{y}(x)$ 满足如下条件
对于任意的 $x$，有

$$g_{y}(x) \ge f(x) \tag{6}$$
其中, $y$ 是一个固定点，并且 $g_{y}(x)$ 必须“着陆”于 $f(x)$，即 $f_{y}(x)$ 必须满足如下条件

$$g_{y}(y) = f(y)$$
则，假设 $x^{*}$ 是 $g_{y}(x)$ 的最小值，则根据(6)，有 $f(x^{*}) \le g_{y}(x^{*})$，又因为 $g_{y}(x^{*})$ 是函数 $g_{y}(x)$ 的最小值，则有 $g_{y}(x^{*}) \le g_{y}(y)$。所以，我们可以得到如下的一个不等式链：
$$f(x^{*}) \le g_{y}(x^{*}) \le g_{y}(y) = f(y) \tag{7}$$

通过(7)我们发现，经过 $g_{y}(x)$ 函数，$f(x)$ “超着” 下降的”方向”走去($f(x^{*}) \le f(y)$)。
若我们反复迭代，则会找到函数 $f(x)$ 的一个局部最小值，若$f(x)$ 是凸函数，则会找到 $f(x)$ 的最小值。
该优化方法的迭代流程如下：

(1) 选择一个初始值 $y = y_{0}$，找到一个函数,使得 $g_{y}(x)$ 满足公式(6)和公式(7)
(2) 第 $i+1$ 步，求解函数 $g_{y}(x)$ ，得到最小值 $y_{i+1}$。
(3) 如果 $f(y_{i}) - f(y_{i+1}) < \epsilon$，则结束，$f(x)$ 的最优值是 $f(y_{i+1})$；否则，另 $y=y_{i+1}$，转到步骤(2)

4. Cauchy-Schwartz inequality(柯西-施瓦茨不等式)

在优化MDS的目标函数的时候，会用到第3节中的基本的最优化理论，而在寻找 $g_{y}(x)$ 函数的时候，就会用到本小节的柯西-施瓦茨不等式。
首先，先给出柯西-施瓦茨不等式，然后证明其不等式成立。

$$(\sum_{i=1}^{n}{x_{i}y_{i}})^{2} \le (\sum_{i=1}^{n}{x_{i}^{2}})(\sum_{i=1}^{n}{y_{i}^{2}})\tag{8}$$
不等式成立的条件是
$$\frac{x_{1}}{y_{1}} = \frac{x_{2}}{y_{2}} = \cdots = \frac{x_{n}}{y_{n}} \tag{9}$$
证明：
考虑一个一元二次方程
$$(x_{1}t+y_{1})^{2}+(x_{2}t+y_{2})^{2}+\cdots+(x_{n}t+y_{n})^{2}=0\tag{10}$$
将(10)转化形式如下
$$(\sum_{i=1}^{n}{x_{i}^{2}})t^{2} + (2\sum_{i=1}^{n}{x_{i}y_{i}})t + \sum_{i=1}^{n}{y_{i}^{2}}=0 \tag{11}$$
很显然，该方程要不没有实数解，要不就是有重根，所以，其判别式 $\Delta \le 0$，所以有
$$\Delta = (2\sum_{i=1}^{n}{x_{i}y_{i}})^{2} - 4(\sum_{i=1}^{n}{x_{i}^{2}})(\sum_{i=1}^{n}{y_{i}^{2}}) \\ =4(\sum_{i=1}^{n}{x_{i}y_{i}})^{2} - 4(\sum_{i=1}^{n}{x_{i}^{2}})(\sum_{i=1}^{n}{y_{i}^{2}}) \le 0 \tag{12}$$
所以有公式(8)成立。当该方程是重根时，有 $\Delta=0$，即 $(x_{1}t+y_{1})=(x_{2}t+y_{2})=\cdots=(x_{n}t+y_{n})=0$，
有公式(9)成立。

证明完毕。
其中，柯西-施瓦茨不等式可以化成如下的形式：

$$\sum_{i=1}^{n}{x_{i}y_{i}} \le \sqrt{(\sum_{i=1}^{n}{x_{i}^{2}})(\sum_{i=1}^{n}{y_{i}^{2}})}\tag{13}$$
我们在优化MDS目标函数的时候，会用到公式(13)。

5. Basic SMACOF methodology

本小节将采用Basic SMACOF methodology对MDS进行优化，其中会用到3、4小节中的内容。
首先，我们将公式(5)展开，拆成3部分。

$$\sigma(\mathbf{Y}) =\sum_{i=1}^{n}{\sum_{j=1}^{n}{(\delta_{ij}-||\mathbf{y}_{i} - \mathbf{y}_{j}||)^{2}}} \\ =\sum_{i=1}^{n}{\sum_{j=1}^{n}{\delta_{ij}^{2}}}+\sum_{i=1}^{n}{\sum_{j=1}^{n}{||\mathbf{y}_{i} - \mathbf{y}_{j}||^{2}}} - 2\sum_{i=1}^{n}{\sum_{j=1}^{n}{\delta_{ij}||\mathbf{y}_{i} - \mathbf{y}_{j}||}} \tag{14}$$

观察公式(14)有，
$\sum_{i=1}^{n}{\delta_{ij}^{2}}$ 是一个常数项，对目标函数优化没有影响。
$\sum_{i=1}^{n}{\sum_{j=1}^{n}{||\mathbf{y}_{i} - \mathbf{y}_{j}||^{2}}}$ 是一个凸的二次函数，有最小值。

而第三项，根据柯西-施瓦茨不等式，可以得知其有界。如下所示：

$$-2\sum_{i=1}^{n}{\sum_{j=1}^{n}{\delta_{ij}||\mathbf{y}_{i} - \mathbf{y}_{j}||}} \\ = -2\sum_{i=1}^{n}{\sum_{j=1}^{n}{\delta_{ij}\sqrt{\sum_{s=1}^{p}(y_{is}-y_{js})^{2}}}}\\ \le -2\sum_{i=1}^{n}{\sum_{j=1}^{n}{\delta_{ij}\frac{\sum_{s=1}^{p}(y_{is}-y_{js})(z_{is}-z_{js})}{\sqrt{\sum_{s=1}^{p}{(z_{is}-z_{js})^{2}}}}}}\tag{15}$$
其中， $z_{ij}$ 是常数。当 $y_{ij}=z_{ij}$ 时，上述不等式变成等式，符合第3小节的优化理论。根据第3小节，可得到：
$$\sigma(\mathbf{Y})=\sum_{i=1}^{n}{\sum_{j=1}^{n}{\delta_{ij}^{2}}}+\sum_{i=1}^{n}{\sum_{j=1}^{n}{||\mathbf{y}_{i} - \mathbf{y}_{j}||^{2}}} - 2\sum_{i=1}^{n}{\sum_{j=1}^{n}{\delta_{ij}||\mathbf{y}_{i} - \mathbf{y}_{j}||}} \\ \le \sum_{i=1}^{n}{\sum_{j=1}^{n}{\delta_{ij}^{2}}}+\sum_{i=1}^{n}{\sum_{j=1}^{n}{||\mathbf{y}_{i} - \mathbf{y}_{j}||^{2}}} -2\sum_{i=1}^{n}{\sum_{j=1}^{n}{\delta_{ij}\frac{\sum_{s=1}^{p}(y_{is}-y_{js})(z_{is}-z_{js})}{\sqrt{\sum_{s=1}^{p}{(z_{is}-z_{js})^{2}}}}}} \\ =\sigma_{\mathbf{Z}}(\mathbf{Y}) \tag{16}$$
则，我们根据第3小节的优化理论，迭代可求解得到 $\sigma(\mathbf{Y})$ 的最优解。

6. Basic SMACOF methodology 矩阵形式

为了便于优化，可以将公式(16)转化成矩阵的形式。对于公式(16)，我们设

$$\sigma(\mathbf{Y})=\sum_{i=1}^{n}{\sum_{j=1}^{n}{\delta_{ij}^{2}}}+\sum_{i=1}^{n}{\sum_{j=1}^{n}{||\mathbf{y}_{i} - \mathbf{y}_{j}||^{2}}} - 2\sum_{i=1}^{n}{\sum_{j=1}^{n}{\delta_{ij}||\mathbf{y}_{i} - \mathbf{y}_{j}||}} \\ = \eta_{\sigma}^{2} + \eta^{2}(\mathbf{Y})-2\rho(\mathbf{Y}) \tag{17}$$
对于公式(17)中的 $\eta^{2}(\mathbf{Y})$，有

$$||\mathbf{y}_{i} - \mathbf{y}_{j}||^{2} \\ =\sum_{s=1}^{p}(y_{is}-y_{js})^{2} \\ =tr(\mathbf{Y^{T}}\mathbf{A}_{ij}\mathbf{Y})\tag{18}$$
其中， $\mathbf{A}_{ij}=(\mathbf{e}_{i}-\mathbf{e}_{j})(\mathbf{e}_{i}-\mathbf{e}_{j})^{\mathbf{T}}$， $a_{ii}=a_{jj}=1$， $a_{ij}=a_{ji}=-1$，其他地方等于0。则
$$\eta^{2}(\mathbf{Y})=\sum_{i=1}^{n}{\sum_{j=1}^{n}{\sum_{s=1}^{p}(y_{is}-y_{js})^{2}}} \\ =\sum_{i=1}^{n}{\sum_{j=1}^{n}{tr(\mathbf{Y^{T}}\mathbf{A}_{ij}\mathbf{Y})}}\\ =tr(\mathbf{Y^{T}}(\sum_{i=1}^{n}{\sum_{j=1}^{n}{\mathbf{A}_{ij}}})\mathbf{Y})\\ =tr(\mathbf{Y^{T}}\mathbf{V}\mathbf{Y}) \tag{19}$$

其中，tr是矩阵的迹。矩阵 $\mathbf{V}$ 是函数 $\eta^{2}(\mathbf{Y})$ 的海森矩阵，其形式如下:

$$\mathbf{V}= \left[ \begin{matrix} 2n-2 &-2 &\cdots &-2 \\ -2 &2n-2 &\cdots &-2 \\ \vdots &\vdots &\ddots &\vdots \\ -2 &-2 &\cdots &2n-2 \\ \end{matrix} \right] \tag{20}$$

对于公式(17)中的 $\rho(\mathbf{Y})$ ，有

$$\delta_{ij}||\mathbf{y}_{i} - \mathbf{y}_{j}|| \\ =\delta_{ij} \sqrt{\sum_{s=1}^{p}(y_{is}-y_{js})^{2}}\\ =\delta_{ij}\frac{\sum_{s=1}^{p}(y_{is}-y_{js})^{2}}{\sqrt{\sum_{s=1}^{p}(y_{is}-y_{js})^{2}}}\\ = tr(\mathbf{Y^{T}}s_{ij}(\mathbf{Y})\mathbf{A}_{ij}\mathbf{Y})\tag{21}$$

其中

$$\begin{equation} s_{ij}(\mathbf{Y})= \begin{cases} \delta_{ij}/d_{ij}(\mathbf{Y}) & \mbox{if } d_{ij}(\mathbf{Y}) > 0,\\ 0 & \mbox{if } d_{ij}(\mathbf{Y}) = 0, \end{cases} \end{equation}$$

则

$$\rho(\mathbf{Y})=\sum_{i=1}^{n}{\sum_{j=1}^{n}{tr(\mathbf{Y^{T}}s_{ij}(\mathbf{Y})\mathbf{A}_{ij}\mathbf{Y})}}\\ =tr(\mathbf{Y^{T}}(\sum_{i=1}^{n}{\sum_{j=1}^{n}{s_{ij}(\mathbf{Y})\mathbf{A}_{ij}}})\mathbf{Y})\\ = tr(\mathbf{Y^{T}}B_{ij}(\mathbf{Y})\mathbf{Y}) \tag{22}$$

综上，

$$\sigma(\mathbf{Y})=\\ 1+tr(\mathbf{Y^{T}VY})-2tr(\mathbf{Y^{T}}B({\mathbf{Y}})\mathbf{Y})\tag{23}$$

根据公式(16)，有

$$\sigma(\mathbf{Y})=\\ 1+tr(\mathbf{Y^{T}VY})-2tr(\mathbf{Y^{T}}B({\mathbf{Y}})\mathbf{Y})\\ \le 1+tr(\mathbf{Y^{T}VY})-2tr(\mathbf{Y^{T}}B({\mathbf{Z}})\mathbf{Z}) = \tau_{\mathbf{Z}}(\mathbf{Y}) \tag{24}$$

至此，我们可以优化 $\tau_{\mathbf{Z}}(\mathbf{Y})$ 来迭代求解 $\sigma(\mathbf{Y})$ 。

$$\frac{\partial{\tau_{\mathbf{Z}}(\mathbf{Y})}}{\partial{\mathbf{Y}}} \\ = 2 \mathbf{VY} - 2B(\mathbf{Z})\mathbf{Z} = \mathbf{0} \tag{25}$$

求解(25)，有

$$\mathbf{Y}^{*}=\mathbf{V}^{-1}B(\mathbf{Z})\mathbf{Z} \tag{26}$$
其中， $\mathbf{V}^{-1}$ 是 $\mathbf{V}$的广义逆。

7. MDS python代码

在sklearn机器学习包中，有MDS的相关代码，示例如下：

import numpy as np
import pandas as pd
from sklearn.manifold import MDS
import matplotlib.pyplot as plt

data = np.array([(5,0,0,0,0,1,0,2,1,0),
                 (0,0,3,0,3,0,1,0,0,1),
                 (2,0,0,0,0,0,1,0,0,0),
                 (1,0,1,0,2,0,0,0,0,1),
                 (5,0,2,0,0,4,2,2,3,7),
                 (0,3,0,1,0,0,0,0,0,0),
                 (0,0,0,6,0,0,0,0,0,1),
                 (0,5,0,0,0,0,0,0,0,0),
                 (0,1,0,0,0,0,0,0,0,0),
                 (0,2,0,0,0,0,0,0,0,0)
                ]
               )
index = ['auto1','auto2','auto3','auto4','auto5','moto1','moto2','moto3','moto4','moto5']
columns = ['car','bike','cars','his','tires','she','ive','her','#k','are']
Word = pd.DataFrame(data,index,columns)

mds = MDS()
mds.fit(data)
a = mds.embedding_
plt.scatter(a[0:5,0],a[0:5,1],color='turquoise')
plt.scatter(a[5:10,0],a[5:10,1],color='red')

8 参考文献

[1] Multidimensional Scaling Using Majorization: SMACOF in R
http://gifi.stat.ucla.edu/janspubs/2008/reports/deleeuw_mair_R_08.pdf

[2] https://en.wikipedia.org/wiki/Stress_majorization