古典多維標度法 (MDS)

本文的閱讀等級:中級

下圖顯示一份手寫數字的樣本,其中每一數字以大小為 16\times 16 像素 (pixel) 的灰階圖片儲存。讀者不妨想像樣本所含的200張數字圖片對應於 \mathbb{R}^{256} 空間的200個數據點。我們提出下面的問題:給定這份樣本資料,如何「目視」數據點於高維空間的散佈?主成分分析 (principal components analysis) 是當今最常採行的一種降維技術 (見“主成分分析”)。在保留數據集的最大變異前提下,將高維數據點正交投影至一個特定的二維空間,此空間由對應樣本共變異數矩陣的最大兩個特徵值的特徵向量擴張而成。如此一來,我們可在平面上觀察所有數據點的投影位置 (稱為主成分係數)。

Digits

手寫數字樣本

 
不過,在某些應用場合,我們僅知道任兩數據點的相異性 (dissimilarity)。舉例來說,手寫數字包含許多變異,如位移、旋轉、伸縮與形變,直接計算兩數字圖片於同一像素位置的灰階差距並不能反映實際的型態差異,我們必須先把兩圖放在可供比較的基準上。為了降低上述變異造成的影響,在比對圖片之前,我們容許一圖 (或兩圖) 些微調整轉變 (見“最小平方法於圖形比對的應用”),並採用各種複雜的圖片相異性算法。因為這些緣故,主成分分析不適用於手寫數字圖片的降維。本文介紹一個建立於數據點的相異性的降維方法,稱為多維標度法 (multidimensional scaling,簡稱 MDS)。下圖顯示手寫數字集經多維標度法處理後得到的二維標度散佈圖。根據相異性的定義,多維標度法可區分為公制 (metric) 與非公制 (nonmetric),前者採用歐幾里得距離 (簡稱歐氏距離),後者則泛指任何非歐氏距離[1]。本文將介紹公制,也稱古典多維標度法,並解說古典多維標度法與主成分分析的關係。

MDS-digits

二維標度散佈圖

 

考慮 p 維歐幾里得空間的 n 的點 \mathbf{x}_1,\ldots,\mathbf{x}_n,其中 \mathbf{x}_r=(x_{r1},\ldots,x_{rp})^T。兩點 \mathbf{x}_r\mathbf{x}_s 的歐氏距離可由下列算式求得:

\displaystyle d_{rs}^2=\sum_{i=1}^p(x_{ri}-x_{si})^2=\Vert\mathbf{x}_r-\mathbf{x}_s\Vert^2=(\mathbf{x}_r-\mathbf{x}_s)^T(\mathbf{x}_r-\mathbf{x}_s)

B=[b_{rs}]n\times n 階內積矩陣,其中 b_{rs} 等於 \mathbf{x}_r\mathbf{x}_s 的內積,b_{rs}=\mathbf{x}_r^T\mathbf{x}_s。因為 \mathbf{x}_r^T\mathbf{x}_s=\mathbf{x}_s^T\mathbf{x}_rB 是對稱矩陣。給定任兩數據點的距離平方 \{d_{rs}^2,r,s=1,\ldots,n\},下面我先解說如何計算內積矩陣 B,隨後推導滿足所有兩數據點之距離的數據集座標 \{\mathbf{x}_r, r=1,\ldots,n\}

 
為了排除平移的影響,我們將數據點的中心置於座標系統的原點,即有

\displaystyle \sum_{r=1}^n\mathbf{x}_r=\mathbf{0}

將歐氏距離平方乘開,

\displaystyle d_{rs}^2=\mathbf{x}_r^T\mathbf{x}_r+\mathbf{x}_s^T\mathbf{x}_s-2\mathbf{x}_r^T\mathbf{x}_s

由上式可推得下列結果 (詳見[2]):

\displaystyle \begin{aligned} \frac{1}{n}\sum_{r=1}^nd_{rs}^2&=\frac{1}{n}\sum_{r=1}^n\mathbf{x}_r^T\mathbf{x}_r+\mathbf{x}_s^T\mathbf{x}_s\\ \frac{1}{n}\sum_{s=1}^nd_{rs}^2&=\mathbf{x}_r^T\mathbf{x}_r+\frac{1}{n}\sum_{s=1}^n\mathbf{x}_s^T\mathbf{x}_s\\ \frac{1}{n^2}\sum_{r=1}^n\sum_{s=1}^nd_{rs}^2&=\frac{2}{n}\sum_{r=1}^n\mathbf{x}_r^T\mathbf{x}_r. \end{aligned}

將前面兩式相加再減去第三式,

\displaystyle \mathbf{x}_r^T\mathbf{x}_r+\mathbf{x}_s^T\mathbf{x}_s=\frac{1}{n}\sum_{r=1}^nd_{rs}^2+\frac{1}{n}\sum_{s=1}^nd_{rs}^2-\frac{1}{n^2}\sum_{r=1}^n\sum_{s=1}^nd_{rs}^2

將上式代入 d_{rs}^2 的展開式,可得

\displaystyle\begin{aligned} b_{rs}&=\mathbf{x}_r^T\mathbf{x}_s\\ &=-\frac{1}{2}\left(d_{rs}^2-\mathbf{x}_r^T\mathbf{x}_r-\mathbf{x}_s^T\mathbf{x}_s\right)\\ &=-\frac{1}{2}\left(d_{rs}^2-\frac{1}{n}\sum_{r=1}^nd_{rs}^2-\frac{1}{n}\sum_{s=1}^nd_{rs}^2+\frac{1}{n^2}\sum_{r=1}^n\sum_{s=1}^nd_{rs}^2\right). \end{aligned}

這個結果看似繁雜,我們將它改寫為較簡明的表達式。令 a_{rs}=-\frac{1}{2}d_{rs}^2。統計學家常採用下列記號:

\displaystyle\begin{aligned} a_{r.}&=\frac{1}{n}\sum_{s=1}^na_{rs}\\ a_{.s}&=\frac{1}{n}\sum_{r=1}^na_{rs}\\ a_{..}&=\frac{1}{n^2}\sum_{s=1}^n\sum_{r=1}^na_{rs}.\end{aligned}

據此,b_{rs}=a_{rs}-a_{r.}-a_{.s}+a_{..}。定義 n\times n 階矩陣 A=[a_{rs}],內積矩陣 B 可表示為

\displaystyle B=CAC

其中 C=I-\frac{1}{n}E,這裡 E 的每一元都等於 1 (證明見[3])。

 
接著討論如何重建數據點的座標。考慮內積矩陣 B 的另一個表達式:

B=\begin{bmatrix} \mathbf{x}_1^T\\ \vdots\\ \mathbf{x}_n^T \end{bmatrix}\begin{bmatrix} \mathbf{x}_1&\cdots&\mathbf{x}_n \end{bmatrix}=XX^T

上面我們令 X=\begin{bmatrix} \mathbf{x}_1&\cdots&\mathbf{x}_n \end{bmatrix}^Tn\times p 階座標矩陣。因為 B 是實對稱半正定矩陣,故可正交對角化為 (見“實對稱矩陣可正交對角化的證明”)

B=QMQ^T

其中 M=\hbox{diag}(\mu_1,\ldots,\mu_n)\mu_1\ge\cdots\ge\mu_{n}\ge 0B 的特徵值,Q=\begin{bmatrix} \mathbf{q}_1&\cdots&\mathbf{q}_n \end{bmatrix}B 的單範正交 (orthonormal) 特徵向量構成的正交矩陣 (orthogonal matrix),Q^T=Q^{-1}。交互乘積不改變矩陣秩,

\displaystyle \hbox{rank}B=\hbox{rank}(XX^T)=\hbox{rank}X\le\min\{n,p\}

在一般情況下,n>p,即有 \hbox{rank}B\le p。所以,B 至多有 p 個非零特徵值 \mu_1,\ldots,\mu_p (即至少有 n-p 個零特徵值),可知

B=\begin{bmatrix} \mathbf{q}_1&\cdots&\mathbf{q}_p \end{bmatrix}\begin{bmatrix} \mu_1&&\\ &\ddots&\\ &&\mu_p \end{bmatrix}\begin{bmatrix} \mathbf{q}_1^T\\ \vdots\\ \mathbf{q}_p^T \end{bmatrix}

因為 B=XX^T,即得座標矩陣

X=\begin{bmatrix} \mathbf{q}_1&\cdots&\mathbf{q}_p \end{bmatrix}\begin{bmatrix} \sqrt{\mu_1}&&\\ &\ddots&\\ &&\sqrt{\mu_p} \end{bmatrix}=\begin{bmatrix} \sqrt{\mu_1}\mathbf{q}_1&\cdots&\sqrt{\mu_p}\mathbf{q}_p \end{bmatrix}

注意,重建的座標矩陣 X 並不唯一存在。設 X'=XP,其中 P 是任意 p\times p 階正交矩陣,P^T=P^{-1},則

X'X'^T=(XP)(XP)^T=XPP^TX^T=XX^T=B

所有 X'=XP 皆為符合要求的座標矩陣。換句話說,旋轉、鏡射 (反射) 與平移不改變數據點之間的距離。

 
看一個例子。考慮 \mathbb{R}^2 的四個點,其座標矩陣為

\hat{X}=\left[\!\!\begin{array}{rr} 2&1\\ 1&4\\ -3&-2\\ 0&-3 \end{array}\!\!\right]

假設四點的座標未知,但我們有它們之間的歐氏距離平方,以矩陣 D_2=[d_{rs}^2] 表示為

D_2=\left[\!\!\begin{array}{rrrr} 0&10&34&20\\ 10&0&52&50\\ 34&52&0&10\\ 20&50&10&0 \end{array}\!\!\right]

則有

\displaystyle A=-\frac{1}{2}D_2=\left[\!\!\begin{array}{rrrr} 0&-5&-17&-10\\ -5&0&-26&-25\\ -17&-26&0&-5\\ -10&-25&-5&0 \end{array}\!\!\right]

內積矩陣即為

\displaystyle\begin{aligned} B&=CAC\\ &=\frac{1}{4}\left[\!\!\begin{array}{rrrr} 3&-1&-1&-1\\ -1&3&-1&-1\\ -1&-1&3&-1\\ -1&-1&-1&3 \end{array}\!\!\right]\left[\!\!\begin{array}{rrrr} 0&-5&-17&-10\\ -5&0&-26&-25\\ -17&-26&0&-5\\ -10&-25&-5&0 \end{array}\!\!\right]\frac{1}{4}\left[\!\!\begin{array}{rrrr} 3&-1&-1&-1\\ -1&3&-1&-1\\ -1&-1&3&-1\\ -1&-1&-1&3 \end{array}\!\!\right]\\ &=\left[\!\!\begin{array}{rrrr} 5&6&-8&-3\\ 6&17&-11&-12\\ -8&-11&13&6\\ -3&-12&6&9 \end{array}\!\!\right]. \end{aligned}

計算正交對角化,B=QMQ^T,結果如下:

\displaystyle Q=\left[\!\!\begin{array}{rrrr} -0.302&0.469&-0.815&-0.154\\ -0.663&-0.365&-0.087&0.648\\ 0.527&-0.618&-0.572&0.114\\ 0.438&0.514&-0.006&0.737 \end{array}\!\!\right],~~M=\left[\!\!\begin{array}{cccc} 36.422&0&0&0\\ 0&7.578&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}\!\!\right]

所以重建的座標矩陣為

\displaystyle X=\left[\!\!\begin{array}{rr} -0.302&0.469\\ -0.663&-0.365\\ 0.527&-0.618\\ 0.438&0.514 \end{array}\!\!\right]\left[\!\!\begin{array}{cc} \sqrt{36.422}&0\\ 0&\sqrt{7.578} \end{array}\!\!\right]=\left[\!\!\begin{array}{rr} -1.823&1.291\\ -4.001&-1.005\\ 3.180&-1.701\\ 2.643&1.415 \end{array}\!\!\right]

原始座標與重建座標之間的差異僅在於相對原點的旋轉。

 
最後我們討論古典多維標度法與主成分分析的關係。令 X 表示一 n\times p 階數據矩陣,X=\begin{bmatrix} \mathbf{x}_1&\cdots&\mathbf{x}_n \end{bmatrix}^T。假設樣本平均數向量為零,則 p\times p 階樣本共變異數矩陣為

\displaystyle S=\frac{1}{n-1}\sum_{r=1}^n\mathbf{x}_r^T\mathbf{x}_r=\frac{1}{n-1}X^TX

因為 S 是實對稱半正定矩陣,故可正交對角化為 S=W\Lambda W^T,其中 \Lambda=\hbox{diag}(\lambda_1,\ldots,\lambda_p)\lambda_1\ge\cdots\ge\lambda_p\ge 0S 的特徵值 (即主成分權值),W=\begin{bmatrix} \mathbf{w}_1&\cdots&\mathbf{w}_p \end{bmatrix} 是單範正交的特徵向量構成的主成分矩陣。內積矩陣 B=XX^T 與樣本共變異數矩陣 S=\frac{1}{n-1}X^TX 有甚麼關係呢?首先,XX^TX^TX 有相同的特徵值,前者再加上 n-p 個零特徵值。理由如下。寫出 BQ=QM 的第 i 行,\displaystyle XX^T\mathbf{q}_i=\mu_i\mathbf{q}_i。上式等號兩邊左乘 X^T

\displaystyle X^TXX^T\mathbf{q}_i=X^TX(X^T\mathbf{q}_i)=\mu_i(X^T\mathbf{q}_i)

另一方面,SW=W\Lambda 的第 i 行是

\displaystyle \frac{1}{n-1}X^TX\mathbf{w}_i=\lambda_i\mathbf{w}_i

比較上面兩式,可得 \lambda_i=\mu_i/(n-1)\mathbf{w}_i=X^T\mathbf{q}_i/\Vert X^T\mathbf{q}_i\Vert。另計算

\displaystyle \Vert X^T\mathbf{q}_i\Vert^2=\mathbf{q}_i^TXX^T\mathbf{q}_i=\mu_i\mathbf{q}_i^T\mathbf{q}_i=\mu_i

使用以上結果,主成分係數矩陣 Z 即為古典多維標度法重建的座標矩陣:

\displaystyle \begin{aligned} Z&=XW=X\begin{bmatrix} \mathbf{w}_1&\cdots&\mathbf{w}_p \end{bmatrix}\\ &=X\begin{bmatrix} X^T\mathbf{q}_1/\Vert X^T\mathbf{q}_1\Vert&\cdots&X^T\mathbf{q}_p/\Vert X^T\mathbf{q}_p\Vert \end{bmatrix}\\ &=XX^T\begin{bmatrix} \mathbf{q}_1/\sqrt{\mu_1}&\cdots&\mathbf{q}_p/\sqrt{\mu_p} \end{bmatrix}\\ &=\begin{bmatrix} \sqrt{\mu_1}\mathbf{q}_1&\cdots&\sqrt{\mu_p}\mathbf{q}_p \end{bmatrix}.\end{aligned}

註解 

[1] 參見IBM SPSS 手冊說明
[2] 第一式推導如下:

\displaystyle  \begin{aligned} \frac{1}{n}\sum_{r=1}^nd_{rs}^2&=\frac{1}{n}\sum_{r=1}^n(\mathbf{x}_r^T\mathbf{x}_r+\mathbf{x}_s^T\mathbf{x}_s-2\mathbf{x}_r^T\mathbf{x}_s)\\ &=\frac{1}{n}\left(\sum_{r=1}^n\mathbf{x}_r^T\mathbf{x}_r+n\mathbf{x}_s^T\mathbf{x}_s-2\left(\sum_{r=1}^n\mathbf{x}_r\right)^T\mathbf{x}_s\right)\\ &=\frac{1}{n}\sum_{r=1}^n\mathbf{x}_r^T\mathbf{x}_r+\mathbf{x}_s^T\mathbf{x}_s. \end{aligned}

將指標 rs 交換,可得第二式:

\displaystyle  \frac{1}{n}\sum_{s=1}^nd_{rs}^2=\mathbf{x}_r^T\mathbf{x}_r+\frac{1}{n}\sum_{s=1}^n\mathbf{x}_r^T\mathbf{x}_s

第三式推導如下:

\displaystyle \begin{aligned} \frac{1}{n^2}\sum_{r=1}^n\sum_{s=1}^nd_{rs}^2&=\frac{1}{n^2}\sum_{r=1}^n\sum_{s=1}^n\left(\mathbf{x}_r^T\mathbf{x}_r+\mathbf{x}_s^T\mathbf{x}_s-2\mathbf{x}_r^T\mathbf{x}_s\right)\\ &=\frac{1}{n^2}\left(n\sum_{r=1}^n\mathbf{x}_r^T\mathbf{x}_r+n\sum_{s=1}^n\mathbf{x}_s^T\mathbf{x}_s-2\sum_{s=1}^n\left(\sum_{r=1}^n\mathbf{x}_r\right)^T\mathbf{x}_s\right)\\ &=\frac{2}{n}\sum_{r=1}^n\mathbf{x}_r^T\mathbf{x}_r. \end{aligned}

[3] 直接將 C=I-\frac{1}{n}E 代入 B=CAC

\displaystyle B=\left(I-n^{-1}E\right)A\left(I-n^{-1}E\right)=A-\frac{1}{n}EA-\frac{1}{n}AE+\frac{1}{n^2}EAE

乘開上式,B(r,s) 元等於

\displaystyle b_{rs}=a_{rs}-\frac{1}{n}\sum_{r=1}^na_{rs}-\frac{1}{n}\sum_{s=1}^na_{rs}+\frac{1}{n^2}\sum_{r=1}^n\sum_{s=1}^na_{rs}

廣告
本篇發表於 機器學習 並標籤為 , , , , 。將永久鏈結加入書籤。

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s