Multidimensional scaling

Multidimensional scaling, 簡稱 MDS, 是個資料分析或是資料降維的工具.

這裡我們要談一下從數學角度來說 MDS 的原理及做法, 更精確的說, 這裡講的是 classical MDS.


假設我們有 $n$ 筆 $p$ 維的資料, 記成 $$ \{x_1, x_2, \cdots, x_n\} \in \mathbb{R}^p $$ 那我們可以據此構造出一個距離平方矩陣 $D$, Euclidian distance matrix, 簡稱 EDM, 其中 $D_{ij} = \|x_i-x_j\|^2$, 也就是 $(x_i)$ 及 $x_j$ 這兩個點距離的平方.


舉例來說

以下我們用 Matlab 隨機生出 $R^2$ 空間中的 $5$ 個點, 並求出他的 EDM:

% Matlab program
X = rand(2,5);                                  % R^2 中的 5 個點
D = squareform(pdist(X', 'squaredeuclidean'));  % 求出其 EDM
D =
     0    0.3839    0.2333    0.2675    0.6373
0.3839         0    0.1571    0.1272    0.0470
0.2333    0.1571         0    0.0028    0.3725
0.2675    0.1272    0.0028         0    0.3222
0.6373    0.0470    0.3725    0.3222         0

可看出 EDM 的幾個性質

  1. 一定是個對稱矩陣
  2. 所有元素都非負
  3. 對角線元素必為零

而 MDS 要做的事是以上的反問題:

  • 假設我們拿到一個 EDM 矩陣, 我想要把它原始生成的點 $x_i$ 找出來.

這樣子的問題稱為 classical MDS.

classical MDS 處理這些距離用 Euclidean distance 來量出來的矩陣.


要講 MDS 做法之前我們先定義一個矩陣稱為置中矩陣 $H$, centering matrix, 也就是把一組資料平移使得其中心為坐標原點: $$ H = I_n - \frac{1}{n}{\bf 1}{\bf 1}^T, $$ 其中 $I_n$ 是 $n\times n$ 的單位矩陣, ${\bf 1}$ 則是元素全為 $1$ 的 $n\times 1$ 向量.

可以輕易看出來 $H$ 是一個對稱矩陣, $H^T=H$.

若我們將原始資料收集一起成一個 $p\times n$ 矩陣 $X = [x_1, \cdots, x_n]$, 也就是把原始資料每一筆當成一個 column 排隊排好, 則有 $$ X H = X (I_n - \frac{1}{n}{\bf 1}{\bf 1}^T) = X - \frac{1}{n}\left(\sum_i x_i\right) {\bf 1}^T = X - \mu {\bf 1}^T = Y, $$ 其中 $\mu=\frac{1}{n}\sum_i x_i \in \mathbb{R}^p$ 就是原始資料的平均, 而 $Y=[y_1, \cdots, y_n]$, $y_i = x_i-\mu$, 也就是把每筆資料減去平均之後記成 $y_i$ 再排排站好. 所以的確 $XH$ 就是把資料平移使得其中心為坐標原點.

推導一下後也可以發現, 若將每筆資料當成一個 row 擺好變成 $\hat{X}$, 則要用 $H\hat{X}$ 將資料做平移使資料中心為原點.


我們有以下這個定理

Theorem

給定原始資料 $X = [x_1, \cdots, x_n]$, 且其相對應的 EDM 為 $D$, 則我們有 $$Y^T Y = -\frac{1}{2}HDH.$$

這個定理證明很簡單.

Proof

(Sketch, not complete, please full-in the details by yourself)

我們先重新整理一下這個 EDM 矩陣: $$ \begin{aligned} D_{ij} &= \|x_i-x_j\|^2 = <x_i-x_j, x_i-x_j> \\ &= <x_i, x_i> + <x_j, x_j> - 2<x_i, x_j> \end{aligned} $$

接著我們定義一個 $n\times 1$ 向量 ${\bf k}$, 其中 ${\bf k}_i = <x_i, x_i>$, 則可以將 $D$ 改寫為 $$D = {\bf k}{\bf 1}^T + {\bf 1}{\bf k}^T - 2 X^TX.$$ 接著兩邊乘上 $H$, 並利用 ${\bf 1}^T H = 0$ 以及 $H {\bf 1} = 0$ (分別將一個常數列向量以及常數行向量置中都會得到零向量), 我們可得 $$HDH = -2HX^TXH = -2(XH)^T(XH) = -2Y^TY.$$ 故得證.


這個定理告訴我們的是, 如果我們只知道一組資料的 EDM 矩陣 $D$, 那要怎樣把資料給還原回來.

Remark

當然, 不可能把原始資料完全還原回來! 如同定理中所述我們所能算出的等號左邊是 $Y$, 也就是置中後的原始資料. 事實上我們也很容易想像, 若將一組資料平移或是旋轉後其 EDM 應該是完全不會變的. 也就是說我們若只知道 EDM, 則其原始資料應該有無限多解. 這裡所謂的還原回來是找到其中一組解, 其他所有可能的解則都可以將之平移旋轉後得到.


定理等號右邊由於是個對稱矩陣, 所以可以對角化為 $$ -\frac{1}{2}HDH = V\Lambda V^T = V\sqrt{\Lambda}\sqrt{\Lambda} V^T, $$ 其中 $\Lambda$ 是個對角矩陣包含所有特徵值, 而 $\sqrt{\Lambda}$ 則是將 $\Lambda$ 對角線元素都開根號.

跟定理對照一下可以輕易地看出來, 平移後的原始資料點可以被還原出來: $Y = \sqrt{\Lambda}V^T$.

這就是 classical MDS 的作法!! 非常簡單.

Remark:

若原始資料 $x_i\in R^p$ (並假設排排站之後的 $X$ 其 rank 為 $p$), 則 $Y^TY$ 的 rank 為 $p$, 必有至少 $n-p$ 個為零的特徵值. 所以, 若我們將 $-\frac{1}{2}HDH$ 對角化後發現有 $m$ 個為零的特徵值, 表示原資料的 $p=n-m$, 那我們就把這些零特徵值都拿掉, 使 $\sqrt{\Lambda}$ 為一個 $p\times n$ 的矩陣, 這樣我們就有 $Y = \sqrt{\Lambda}V^T \in R^{p\times n}$.

Remark:

若將 $Y$ 做奇異值分解 (Singular Value Decomposition, SVD), 得到 $$ Y = \hat{U}\hat{\Sigma}\hat{V}^T, \quad \hat{U}\in R^{p\times p}, \quad \hat{\Sigma}\in R^{p\times n}, \quad \hat{V}\in R^{n\times n}, $$ 並且 $\hat{U}^T\hat{U}=I_p$ 以及 $\hat{V}^T\hat{V}=I_n$. 則 $Y^TY = \hat{V}\hat{\Sigma}^T\hat{\Sigma}\hat{V}^T$. 對照一下定理可以看出, 如果 $D$ 是個 EDM 矩陣, 那他的 eigenvalues 一定都是非負.

依照上面兩個 remark 稍微整理一下 MDS 的做法:

  1. 將 EDM 矩陣 $D$ 做 double centering 並乘以 $-1/2$ 求出 $-\frac{1}{2}HDH$.
  2. 做對角化得到 $\Lambda$ 以及 $V$
  3. 拿掉所有零特徵值及其相對應的特徵向量, 求出 $Y = \sqrt{\Lambda}V^T$

一般來說我們會將 MDS 當成一個降維的工具, 希望在低維度(二或三, 或 $k$)找到一組資料使得其 EDM 與原始資料的 EDM 最像. 而作法就是保留前 $k$ 個特徵值及其特徵向量.

給定一個 $n\times n$ 的 EDM 以及欲投影的維度 $k$, MDS 簡單的程式如下:

% Matlab program
function Y = multidimensional_scaling(D, k)

%   Input: D, n*n EDM 矩陣, 元素為距離平方
%          k, 要降到的維度, k 為正整數, k<=p
%   Output: Y: k*n data matrix, k 個 features 以及 n 個 samples

    n = size(D, 1);                         % n 個 samples
    mu = sum(D, 1)/n; D = D - mu;           %
    mu = sum(D, 2)/n; B = D - mu;           %
    B = -0.5*B;                             % B = -0.5*H*D*H
    [V, D] = eig(B, 'vector');              % 求出特徵值及特徵向量
    [sqD, ind] = sort(sqrt(D), 'descend');  % 將特徵值按大小排列
    sqD = sqD(1:k);                         % 取前 k 個
    V = V(:, ind(1:k));                     % 取相對應的特徵向量
    Y = V'.*(sqD*ones(1,n));                % 求出 k 維資料點
end

Extension

廣義一點的 MDS 可以想像是, 我們拿到的距離也許不是用歐式距離量的. 比如說手裡有許多照片, 照片與照片兩兩之間的距離就有非常多種測量的方式. 而不管是用什麼方式量的, 只要他們是 距離(metric) , 我們就可以定出一個 $D$ 矩陣, 其中 $d_{ij}$ 就是 sample points $x_i$ 與 $x_j$ 之間的距離.

接著我們可以定義一個 cost function, 稱之為 stress: $$ Stress_D(x_1, \cdots, x_n) = \left(\sum_i\sum_j\left(d_{ij} - \|x_i-x_j\|\right)^2\right)^{1/2}. $$ Metric Multidimensional scaling (mMDS) 要做的事就是要找到一組資料點 $\{x_1, \cdots, x_n\}$ 使得上式 stress 有最小值.

Remark:

若 $D$ 是用廣義距離造出來的, 就無法保證 $-\frac{1}{2}HDH$ 這矩陣的特徵值都是正的. 不過我們依然可以用 classical MDS 的做法來做, 只是這時候我們會將所有負的特徵值全都丟掉.


Extension

原始點資料若是在某個 weighted inner product space 裡, 內積定義為 $$ <x, y>_Q = x^TQy. $$ 則我們可以將原本的定理推廣如下

Theorem

$$Y^T QY = -\frac{1}{2}HDH.$$


Extension: classical MDS v.s. PCA

給定 $n$ 筆 $p$ 維的資料, 記成 $$ \{x_1, x_2, \cdots, x_n\} \in R^p, $$

  1. 我們可以造出其 EDM, 再用 classical MDS 投影到 $k$ 維.
  2. 我們也可以直接用 PCA 投影到 $k$ 維.

有趣的是以上兩個做法得到完全相同的結果.

不管 PCA 或 classical MDS 最後都是考慮置中後的資料, 所以我們只需看 $Y$ 即可, $y_i = x_i-\mu$, 而 $\mu$ 是平均數.

  1. PCA 是對 $\Sigma = YY^T$ 做 spectral decomposition
  2. classical MDS 是對 $Y^TY$ 做 spectral decomposition

若將 $Y$ 做 SVD 得到 $Y = \hat{U}\hat{\Sigma}\hat{V}^T$, 則

  1. PCA 我們知道最後投影的結果為 $B = \hat{\Sigma}\hat{V}^T$
  2. classical MDS 我們也是投影到 $\hat{\Sigma}\hat{V}^T$

所以雖然出發點不同, 不過結果真的一模一樣.

PCA = classical MDS.


References:

  1. NCTS mini-course on manifold learning
  2. What is principal component analysis?
  3. GeostatsGuy Lectures - Multidimensional Scaling

Avatar
Te-Sheng Lin (林得勝)
Associate Professor

The focus of my research is concerned with the development of analytical and computational tools, and further to communicate with scientists from other disciplines to solve engineering problems in practice.

Related