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