多维标度法MDS古典解的证明与R语言实现

文章目录

写在前面

最近忙于作业一直没有更新，后续会慢慢补上。这次写一下前一阶段学习的多维标度 (Multi-Dimensional Scaling, MDS) 分析的古典解部分，包括欧式型距离矩阵判定定理的证明以及利用R语言进行MDS的古典解实现。

MDS概念与基本思想

MDS是以空间分布的形式表现对象之间相似性或亲疏关系的一种多元数据分析方法，其基本思想是将高维坐标中的点投影到低维空间中，保持点彼此之间的相似性尽可能不变。其结果主要是偏好图(又称多维标度图)等。

一些基本概念与定义

距离阵

一个$n\times n$矩阵$D=(d_{ij})$，若满足 $D^T=D$, $\,d_{ii}=0$, $\,d_{ij}\geqslant0$, $\,(i,\,j=1,\,2,\,\cdots,\,n;\,i\neq j)$，则称 $\!D\!$ 为距离阵。

欧式型距离阵

一个距离阵$D=(d_{ij})$称为欧式型的，若存在某个正整数$p$及$p$维空间$\mathbb{R}^p$中的$n$个点 $x_1,\,\cdots,\,x_n$, 使得
$\begin{aligned} d^2_{ij} &= (x_i-x_j)'(x_i-x_j), \\ \text{其中，令}\\ A &= (a_{ij}),\ \ \,\ a_{ij}=-\frac12d_{ij}^2, \\ B &= H'AH,\ H=I_n-\frac1n1_n1'_n. \\ \end{aligned}$

欧式型距离阵判定定理

一个$n\times n$的距离阵$D$是欧式型的充要条件是:$B\geqslant0$.

证明 ★ \bigstar ★

必要性

设距离阵$D$是欧式型的，则由欧式型距离阵的定义可知，存在$x_1,\,\cdots,\,x_n\in\mathbb{R}^p$，使得
$d_{ij}^2=-2a_{ij}=(x_i-x_j)'(x_i-x_j),$
由$B = H'AH,\ H=I_n-\frac1n\mathbf{1}_n\mathbf{1}'_n$可得：
$B=H'AH=A-\dfrac1nAJ-\dfrac1nJA+\dfrac1{n^2}JAJ,$
上式中，$J=\mathbf{1}_n\mathbf{1}_n'$。并注意到
​
$\dfrac1nAJ=\begin{bmatrix} \bar{a}_{1.}\\ \bar{a}_{2.}\\ \vdots\\ \bar{a}_{n.}\\ \end{bmatrix} \mathbf{1}_n',\ \dfrac1nJA=\mathbf{1}_n(\bar{a}_{.1},\,\bar{a}_{.2},\,\cdots,\,\bar{a}_{.n}),\ \dfrac1{n^2}JAJ=\bar{a}{..}\mathbf{1}_n\mathbf{1}_n',$

其中，
$\bar{a}_{i.}=\dfrac1n\sum\limits_{j=1}^na_{ij},\ \bar{a}_{.j}=\dfrac1n\sum\limits_{i=1}^na{ij},\ \bar{a}_{..}=\dfrac1{n^2}\sum\limits_{i=1}^n\sum\limits_{j=1}^na_{ij}.$

将上述各式代入 $B=H'AH=A-\dfrac1nAJ-\dfrac1nJA+\dfrac1{n^2}JAJ$

$b_{ij}=a_{ij}-\bar{a}_{i.}-\bar{a}_{.j}+\bar{a}_{..},$

再由$d_{ij}^2=-2a_{ij}=(x_i-x_j)'(x_i-x_j)$式可分别求得$a_{ij},\,\bar{a}_{i.},\,\bar{a}_{.j},\,\bar{a}_{..}$，将其代入式$b_{ij}=a_{ij}-\bar{a}_{i.}-\bar{a}_{.j}+\bar{a}_{..}$，得到
$\begin{aligned} b_{ij} =&a_{ij}-\bar{a}_{i.}-\bar{a}_{.j}+\bar{a}_{..}\\ =&a_{ij}-\frac1n\sum_{j=1}^na_{ij}-\frac1n\sum_{i=1}^na_{ij}+\frac1{n^2}\sum_{i=1}^n\sum_{j=1}^na_{ij} \\ =&-\frac12(x_i^2+x_j^2-2x_ix_j)\\ &+\frac{1}{2n}\left[\sum_{j=1}^n(x_i^2+x_j^2-2x_ix_j)+\sum_{i=1}^n(x_i^2+x_j^2-2x_ix_j)\right]\\ &-\frac{1}{2n^2}\sum_{i=1}^n\sum_{j=1}^n(x_i^2+x_j^2-2x_ix_j)\\ =&-\frac12(x_i^2+x_j^2-2x_ix_j)\\ &+\frac{1}{2n}\left(nx_i^2+\sum_{j=1}^nx_j^2-2x_i\sum_{j=1}^nx_j+\sum_{i=1}^nx_i^2+nx_j^2-2x_j\sum_{i=1}^nx_i\right)\\ &-\frac{1}{2n^2}\sum_{i=1}^n\left(nx_i^2+\sum_{j=1}^nx_j^2-2x_i\sum_{j=1}^nx_j\right) \\ =&x_ix_j+\frac{1}{2n}\left(\sum_{j=1}^nx_j^2-2x_i\sum_{j=1}^nx_j+\sum_{i=1}^nx_i^2-2x_j\sum_{i=1}^nx_i\right)\\ &-\frac{1}{2n}\left(\sum_{i=1}^nx_i^2+\sum_{j=1}^nx_j^2-\frac2n\sum_{i=1}^nx_i\sum_{j=1}^nx_j\right) \\ =&x_ix_j-\frac{1}{n}\left(x_i\sum_{j=1}^nx_j+x_j\sum_{i=1}^nx_i\right)+\frac{1}{n^2}\sum_{i=1}^nx_i\sum_{j=1}^nx_j \\ =&x_ix_j-x_i\bar{x}-x_j\bar{x}+\bar{x}^2=(x_i-\bar{x})'(x_j-\bar{x}) \end{aligned}$

其中$\bar{x}=\dfrac1n\sum\limits_{i=1}^nx_i$，即对矩阵B的每个元素，有
$b_{ij}=(x_i-\bar{x})'(x_j-\bar{x}),$

将上式写成矩阵形式，并注意到
$\begin{aligned} H=I_n-\frac1n\mathbf{1}_n\mathbf{1}'_n &=\begin{bmatrix} 1&&\\ &1&\\ &&\ddots&\\ &&&1 \end{bmatrix}_{n\times n}- \frac1n \begin{bmatrix} 1&1&\cdots&1\\ 1&1&\cdots&1\\ \vdots&\vdots&\ddots&\vdots\\ 1&1&\cdots&1 \end{bmatrix}_{n\times n}\\ &= \begin{bmatrix} 1-\frac1n&-\frac1n&\cdots&-\frac1n\\ -\frac1n&1-\frac1n&\cdots&-\frac1n\\ \vdots&\vdots&\ddots&\vdots\\ -\frac1n&-\frac1n&\cdots&1-\frac1n \end{bmatrix}_{n\times n} \end{aligned}$

$X=\begin{bmatrix} x_{11}&x_{12}&\cdots&x_{1p}\\ x_{21}&x_{22}&\cdots&x_{2p}\\ \vdots&\vdots&\ddots&\vdots\\ x_{n1}&x_{n2}&\cdots&x_{np} \end{bmatrix}_{n\times p}$

$\begin{aligned} HX&=\begin{bmatrix} 1-\frac1n&-\frac1n&\cdots&-\frac1n\\ -\frac1n&1-\frac1n&\cdots&-\frac1n\\ \vdots&\vdots&\ddots&\vdots\\ -\frac1n&-\frac1n&\cdots&1-\frac1n \end{bmatrix} \begin{bmatrix} x_{11}&x_{12}&\cdots&x_{1p}\\ x_{21}&x_{22}&\cdots&x_{2p}\\ \vdots&\vdots&\ddots&\vdots\\ x_{n1}&x_{n2}&\cdots&x_{np} \end{bmatrix}_{n\times p} \\\\ &=\begin{bmatrix} x_{11}-\frac{\sum\limits_{i=1}^n x_{i1}}{n}&x_{12}-\frac{\sum\limits_{i=1}^n x_{i2}}{n}&\cdots&x_{1p}-\frac{\sum\limits_{i=1}^n x_{ip}}{n}\\ x_{21}-\frac{\sum\limits_{i=1}^n x_{i1}}{n}&x_{22}-\frac{\sum\limits_{i=1}^n x_{i2}}{n}&\cdots&x_{2p}-\frac{\sum\limits_{i=1}^n x_{ip}}{n}\\ \vdots&\vdots&\ddots&\vdots\\ x_{n1}-\frac{\sum\limits_{i=1}^n x_{i1}}{n}&x_{n2}-\frac{\sum\limits_{i=1}^n x_{i2}}{n}&\cdots&x_{np}-\frac{\sum\limits_{i=1}^n x_{ip}}{n} \end{bmatrix}_{n\times p}\\ &=\begin{bmatrix} x_{11}-\bar{x}&x_{12}-\bar{x}&\cdots&x_{1p}-\bar{x}\\ x_{21}-\bar{x}&x_{22}-\bar{x}&\cdots&x_{2p}-\bar{x}\\ \vdots&\vdots&\ddots&\vdots\\ x_{n1}-\bar{x}&x_{n2}-\bar{x}&\cdots&x_{np}-\bar{x} \end{bmatrix}_{n\times p}\\ &=(X_1-\overline{X},\,X_2-\overline{X},\,\cdots,\,X_n-\overline{X})'\\ \end{aligned}$

其中$X_i=(x_{i1},\,x_{i2},\,\cdots,\,x_{ip}),\ \ i=1,\,2,\,\cdots,\,n$.

于是得到$B=(HX)(HX)'$，用任一$n$维列向量$\alpha$构造二次型可知，
$\alpha' B\alpha=\alpha'(HX)(HX)'\alpha=\alpha'(HX)[\alpha'(HX)]'=\|\alpha' HX\|^2\geqslant0,$
所以$B$为半正定阵，必要性得证，下证充分性。

充分性

记$p=\mathrm{rank}(B),\,\lambda_1,\,\lambda_2,\,\cdots,\,\lambda_p$为$B$的正特征根，$x_{(1)},\,\cdots,\,x_{(p)}$为对应的特征向量。

由已知条件，$B\geqslant0$，根据施密特正交化方法得

$B=H'AH=\varGamma\Lambda\varGamma'$

其中$\Lambda=\mathrm{diag}(\lambda_1,\,\lambda_2,\,\cdots,\,\lambda_p),\,\lambda_1\geqslant\cdots\geqslant\lambda_p$为$B$的$p$个正特征根$\varGamma=X\Lambda^{-\frac12},\,\varGamma$的$p$个列为对应的$p$个标准正交化的特征向量。

取$X=\varGamma\Lambda^{\frac12}$，为一$n\times p$阶矩阵。将$X$写成$X=(x_1,\,x_2,\,\cdots,\,x_n)'=(x_{(1)},\,x_{(2)},\,\cdots,\,x_{(p)})$，于是有
$X'X=(\varGamma\Lambda^{\frac12})'(\varGamma\Lambda^{\frac12})=\Lambda,\ B=XX',$
即$b_{ij}=x_i'x_j$。由此得到$x_i$与$x_j$两点的距离平方
$\begin{aligned} (x_i-x_j)'(x_i-x_j) &= x_i'x_i-2x_i'x_j+x_j'x_j=b_{ii}-2b_{ij}+b_{jj}\\ &= a_{ii}-2a_{ij}+a_{jj}= -2a_{ij}=d_{ij}^2 \end{aligned}$
根据上述推导可知：存在正整数$p$和一个$n\times p$阶矩阵$X=x_1,\,x_2,\,\cdots,\,x_n)=\varGamma\Lambda^{\frac12}$，使得$X$是$D$的构造点，所以$D$为欧式型。

□

MDS古典解计算步骤

1. 由距离阵$D=(d_{ij})$ 构造矩阵$A=(a_{ij})=-\dfrac12d_{ij}^2$;
2. 令矩阵$B=(b_{ij})$, 其中$b_{ij}=a_{ij}-\overline{a}_{i.}-\overline{a}_{.j}+\overline{a}_{..}$;
3. 求矩阵$B$的特征根$\lambda_1\geqslant\lambda_2\geqslant\cdots\geqslant\lambda_n$, 若无负特征根，表明$B\geqslant 0$, 从而距离阵$D$是欧式型的; 若有负特征根，距离阵$D$一定不是欧式型的。此时令
   $a_{1,\,k}=\dfrac{\sum\limits_{i=1}^k\lambda_i}{\sum\limits_{i=1}^n|\lambda_i|},\quad a_{2,\,k}=\dfrac{\sum\limits_{i=1}^k\lambda_i^2}{\sum\limits_{i=1}^n\lambda_i^2},$
   上面两个量相当于主成分分析(PCA)中的累积贡献率;
4. 令$\hat{X}=\big(\hat{x}_{(1)},\,\cdots,\,\hat{x}_{(k)}\big)$, 则$\hat{X}$的行向量$x_1,\,\cdots,\,x_n$即为待求的古典解。

R语言实现

方法一：使用内置的cmdscale()函数

直接输入距离阵即可进行求解

方法二：使用自定义函数

根据定义及计算步骤编写自定义函数

# 输入距离阵D
calc_mds <- function(D) {
    dim <- nrow(D)
    A <- -1/2*D^2
    onen <- matrix(rep(1, dim), nrow=dim)
    H <- diag(dim) - 1/dim * onen %*% t(onen)
    B <- t(H) %*% A %*% H
    val <- eigen(B, symmetric=F)$values
    vec <- eigen(B, symmetric=F)$vectors
    if (all(val>=0)) {
        MDS <- vec[,1:2] %*% diag(sqrt(val[1:2])); 
        data.frame(points = round(MDS, 3))
        }
    else 
        print("距离阵非欧式型！")
    
}

MDS <- calc_mds(D)
plot(MDS)

参考文献

[1] 王斌会. 多元统计分析及R语言建模(第四版). 暨南大学出版社. 2016.