主成分分析與低秩矩陣近似

本文的閱讀等級:高級

假設我們有一筆維數等於 p,樣本大小為 n 的數據 \{\mathbf{x}_1,\ldots,\mathbf{x}_n\},也就是說每一個數據點 \mathbf{x}_i=(x_{i1},\ldots,x_{ip})^T\in\mathbb{R}^p 包含 p 個變數的量測值。沿用統計學與數據科學的慣例 (見“數據矩陣的列與行”),定義 n\times p 階數據矩陣

X=\begin{bmatrix}  \mathbf{x}_1^T\\  \vdots\\  \mathbf{x}_n^T  \end{bmatrix}=\begin{bmatrix}  x_{11}&\cdots&x_{1p}\\  \vdots&\ddots&\vdots\\  x_{n1}&\cdots&x_{np}  \end{bmatrix}

其中 x_{ij} 代表第 j 個變數的第 i 個量測值,i=1,\ldots,nj=1,\ldots,p。在不造成混淆的情況下,以下用 x_j 表示第 j 個變數。如果數據包含大量的變數 (p 很大) 或者變數之間存在顯著的共線性關係[1],你可以設計一個從向量空間 \mathbb{R}^p 映至 \mathbb{R}^k 的線性映射,1\le k<p,數據點 \mathbf{x}_1,\ldots,\mathbf{x}_n 經映射後的像 (image) 構築另一筆變數較少且兩兩變數不存在線性相關性的新數據,這個方法稱為主成分分析 (principal components analysis)。從統計學的觀點,主成分分析的目的是找到少量的新變數,稱為降維 (dimension reduction),同時盡可能地保留變數的總變異量。從線性代數的觀點,主成分分析其實是一種矩陣近似法,我們希望得到一個最近似於原數據矩陣 X 的低秩 (low rank) 同尺寸矩陣。本文證明證明主成分分析與低秩矩陣近似在本質上是相同的問題。

 
先回顧主成分分析的問題陳述與主要結果 (見“主成分分析”)。因為我們不關心變數的平均數,假設數據矩陣 X 的每一行 (column,對應一個變數) 的平均數都等於零 (如果不成立,將該行所有元減去平均數),即 n 個數據點的中心位於原點,\sum_{i=1}^n\mathbf{x}_i=\mathbf{0}。定義 p\times p 階樣本變異數矩陣 (sample covariance matrix)

\displaystyle  S=[s_{ij}]=\frac{1}{n-1}\sum_{i=1}^n\mathbf{x}_i\mathbf{x}_i^T=\frac{1}{n-1}X^TX

其中 s_{ij} 為變數 x_ix_j 的樣本共變異數 (見“樣本平均數、變異數和共變異數”)。明顯地,S 是一個實對稱半正定矩陣,設正交對角化形式為

\displaystyle  S=Q\Lambda Q^T

其中 \Lambda=\mathrm{diag}(\lambda_1,\ldots,\lambda_p) 是特徵值矩陣,\lambda_1\ge\cdots\ge\lambda_p\ge 0 代表主成分的權值,Q=\begin{bmatrix}  \mathbf{q}_1&\cdots&\mathbf{q}_p  \end{bmatrix} 是單範正交 (orthonormal) 的特徵向量, \mathbf{q}_i^T\mathbf{q}_j=\delta_{ij}[2],組成的 p\times p 階正交主成分矩陣,Q^T=Q^{-1}。數據樣本 X 包含 p 個變數的總變異量為

\displaystyle\begin{aligned}  \sum_{j=1}^ps_{jj}&=\hbox{trace}S=\hbox{trace}(Q\Lambda Q^T)=\hbox{trace}(\Lambda Q^TQ)\\  &=\hbox{trace}\Lambda=\lambda_1+\cdots+\lambda_p,\end{aligned}

上面使用了跡數的循環不變性。主成分分析要找到一組 p 維單範正交向量 \{\mathbf{w}_1,\ldots,\mathbf{w}_k\},稱為主軸,使得 \sum_{j=1}^k\mathbf{w}_j^TS\mathbf{w}_j 有最大值,我們用下列約束優化 (constrained optimization) 問題表示:

\displaystyle  \begin{array}{ll}  \hbox{maximize}&\displaystyle\sum_{j=1}^k\mathbf{w}_j^TS\mathbf{w}_j\\    \hbox{subject to}&\mathbf{w}_i^T\mathbf{w}_j=\delta_{ij}.\end{array}

z_{ij}=\mathbf{w}_j^T\mathbf{x}_i 為新變數,稱為主成分,z_j 的第 i 個量測值,稱為主成分係數,\mathbf{w}_j^TS\mathbf{w}_jz_j 的樣本變異數 (稍後會說明)。主成分分析給出的主軸就是樣本共變異數矩陣 S 的單範正交特徵向量 \mathbf{w}_j=\mathbf{q}_j,新變數 z_j 的樣本變異數為特徵值 \lambda_jj=1,\ldots,k

 
考慮這個秩受限的矩陣近似問題 (見“奇異值分解於矩陣近似的應用”):給定一個 n\times p 階實矩陣 X,設 \hat{X}n\times p 階實矩陣,

\displaystyle  \begin{array}{ll}  \hbox{minimize}&\Vert X-\hat{X}\Vert_F^2\\    \hbox{subject to}&\hbox{rank}\hat{X}=k,\end{array}

其中 \Vert\cdot\Vert_F 是 Frobenius 範數。我們一邊引進符號一邊說明解法。欲產生矩陣 \hat{X} 滿足 \hbox{rank}\hat{X}=k,最簡單的線性代數方法是將 X 所包含的數據點 \mathbf{x}_i\in\mathbb{R}^p 正交投影至一個維數等於 k 的子空間,並以投影 \hat{\mathbf{x}}_i 近似數據點 \mathbf{x}_i。令 \mathcal{W}\mathbb{R}^p 的一個子空間,\dim\mathcal{W}=k1\le k<p,並令 \{\mathbf{w}_1,\ldots,\mathbf{w}_k\}\mathcal{W} 的一組單範正交基底,合併為 p\times k 階基向量矩陣 W=\begin{bmatrix}  \mathbf{w}_1&\cdots&\mathbf{w}_k  \end{bmatrix}。因此,(W^TW)_{ij}=\mathbf{w}_i^T\mathbf{w}_j=\delta_{ij},也就是說 W^TW=I_k (但 WW^T\neq I_p,除非 k=p)。我們採用單範正交基底的用意一方面在於簡化正交投影計算,另一方面是為了使新變數之間不存在線性相關性。寫出值域等於子空間 \mathcal{W} 的正交投影矩陣 (見“正交投影──威力強大的線代工具”)

\displaystyle  P=W(W^TW)^{-1}W^T=WW^T

數據點 \mathbf{x}_k 至子空間 \mathcal{W} 的正交投影為 \hat{\mathbf{x}}_i=P\mathbf{x}_i=WW^T\mathbf{x}_i。接下來,定義新數據點 \mathbf{z}_i=(z_{i1},\ldots,z_{ik})^T=W^T\mathbf{x}_i\in\mathbb{R}^k,其中 z_{ij}=\mathbf{w}_j^T\mathbf{x}_i。基向量矩陣的轉置 W^T:\mathbb{R}^p\to\mathbb{R}^k 就是實現降維的線性映射。在子空間 \mathcal{W} 的投影 \hat{\mathbf{x}}_i 可表示為基向量 \mathbf{w}_1,\ldots,\mathbf{w}_k 的唯一線性組合:

\hat{\mathbf{x}}_i=W\mathbf{z}_i=\begin{bmatrix}  \mathbf{w}_1&\cdots&\mathbf{w}_k  \end{bmatrix}\begin{bmatrix}  z_{i1}\\  \vdots\\  z_{ik}  \end{bmatrix}=z_{i1}\mathbf{w}_1+\cdots+z_{ik}\mathbf{w}_k

換句話說,新數據點 \mathbf{z}_i 是投影 \hat{\mathbf{x}}_i 參考基底 (座標系統) \{\mathbf{w}_1,\ldots,\mathbf{w}_k\} 的座標向量,\mathbf{z}_i\in\mathbb{R}^k\hat{\mathbf{x}}_i\in\mathbb{R}^p 具有一對一的對應關係。在實際應用上,我們用下面的 n\times k 階座標向量矩陣 Z 取代 n\times p 階數據矩陣 X 以達到降維的目的:

Z=\begin{bmatrix}  \mathbf{z}_1^T\\  \vdots\\  \mathbf{z}_n^T  \end{bmatrix}=\begin{bmatrix}  \mathbf{x}_1^TW\\  \vdots\\  \mathbf{x}_n^TW  \end{bmatrix}=XW

投影 \hat{\mathbf{x}}_i 構成的近似數據矩陣則為

\hat{X}=\begin{bmatrix}  \hat{\mathbf{x}}_1^T\\  \vdots\\  \hat{\mathbf{x}}_n^T  \end{bmatrix}=\begin{bmatrix}  \mathbf{z}_1^TW^T\\  \vdots\\  \mathbf{z}_n^TW^T  \end{bmatrix}=XWW^T

請注意,新數據矩陣 Z 的每一行的平均數都等於零,

\displaystyle  \sum_{i=1}^nz_{ij}=\sum_{i=1}^n\mathbf{w}_j^T\mathbf{x}_i=\mathbf{w}_j^T\left(\sum_{i=1}^n\mathbf{x}_i\right)=\mathbf{w}_j^T\mathbf{0}=0,~~j=1,\ldots,k,

便有 \sum_{i=1}^n\mathbf{z}_i=\mathbf{0}

 
在秩受限的條件下,前述矩陣近似問題可以表示成底下的約束優化問題:

\displaystyle  \begin{array}{ll}  \hbox{minimize}&\displaystyle\sum_{i=1}^n\left\| \mathbf{x}_i-WW^T\mathbf{x}_i\right\|^2=\left\| X-XWW^T\right\|^2_F\\    \hbox{subject to}&W^TW=I_k.  \end{array}

使用 Frobenius 範數的跡數表達式 \Vert A\Vert_F^2=\hbox{trace}(A^TA)=\hbox{trace}(AA^T) (見“矩陣範數”),可得

\displaystyle  \begin{aligned}  \left\| X-XWW^T\right\|^2_F&=\hbox{trace}\left(\left(X-XWW^T\right)\left(X-XWW^T\right)^T\right)\\  &=\hbox{trace}\left(XX^T\right)-2\,\hbox{trace}\left(XWW^TX^T\right)+\hbox{trace}\left(XWW^TWW^TX^T\right)\\  &=\hbox{trace}\left(XX^T\right)-\hbox{trace}\left(XWW^TX^T\right)\\  &=\Vert X\Vert_F^2-\Vert XW\Vert_F^2.  \end{aligned}

因為 \Vert X\Vert_F 是常數,我們有形式較簡化的等價矩陣近似問題:

\displaystyle  \begin{array}{ll}  \hbox{maximize}&\Vert XW\Vert_F^2\\  \hbox{subject to}&W^TW=I_k.  \end{array}

 
接著證明秩受限的矩陣近似與主成分分析是相同的問題。從主成分分析的角度來看,因為 \sum_{i=1}^nz_{ij}=0

\displaystyle  \begin{aligned}  \mathbf{w}_j^TS\mathbf{w}_j&=\frac{1}{n-1}\mathbf{w}_j^TX^TX\mathbf{w}_j=\frac{1}{n-1}\left\|X\mathbf{w}_j\right\|^2\\  &=\frac{1}{n-1}\left\|\begin{bmatrix}  \mathbf{x}_1^T\mathbf{w}_j\\  \vdots\\  \mathbf{x}_n^T\mathbf{w}_j  \end{bmatrix}\right\|^2=\frac{1}{n-1}\sum_{i=1}^nz_{ij}^2,\end{aligned}

二次型 \mathbf{w}_j^TS\mathbf{w}_j 是新變數 z_j 的樣本變異數。再者,

\displaystyle   \Vert XW\Vert_F^2=\Vert Z\Vert_F^2=\sum_{j=1}^k\left(\sum_{i=1}^nz_{ij}^2\right)=(n-1)\sum_{j=1}^k\mathbf{w}_j^TS\mathbf{w}_j

因此,最小化矩陣誤差平方 \Vert X-\hat{X}\Vert_F^2 等同於最大化新變數 z_1,\ldots,z_k 的總變異量 \sum_{j=1}^k\mathbf{w}_j^TS\mathbf{w}_j=\frac{1}{n-1}\Vert Z\Vert_F^2 (常數 n-1 不影響 W 的最佳值)。

 
下面我們使用奇異值分解推導 p\times k 階基向量矩陣 W=\begin{bmatrix}  \mathbf{w}_1&\cdots&\mathbf{w}_k  \end{bmatrix} (見“奇異值分解 (SVD)”)。令 n\times p 階數據矩陣 X 的奇異值分解為 X=U\Sigma V^T,其中 Un\times n 階正交矩陣,U^T=U^{-1}Vp\times p 階正交矩陣,V^T=V^{-1},且 \Sigma=\begin{bmatrix}  D&0\\  0&0  \end{bmatrix},分塊 D=\hbox{diag}(\sigma_1,\ldots,\sigma_r) 是由非零奇異值 \sigma_1\ge\cdots\ge\sigma_r>0 構成的 r\times r 階對角矩陣,r=\hbox{rank}X\le\min\{n,p\}。將奇異值分解代入目標函數,

\displaystyle  \begin{aligned}  \left\|XW\right\|_F^2&=\left\|U\Sigma V^TW\right\|_F^2=\hbox{trace}\left((U\Sigma V^TW)^T(U\Sigma V^TW)\right)\\  &=\hbox{trace}\left(W^TV\Sigma^TU^TU\Sigma V^TW\right)=\hbox{trace}\left(P^T\Sigma^T\Sigma P\right)\\  &=\Vert \Sigma P\Vert_F^2,  \end{aligned}

上面令 p\times k 階矩陣 P=V^TW。請注意,P^TP=W^TVV^TW=W^TW=I_k。我們得到一個更簡單的等價問題:

\displaystyle  \begin{array}{ll}  \hbox{maximize}&\Vert \Sigma P\Vert_F^2\\  \hbox{subject to}&P^TP=I_k.\end{array}

你只要找到 P=[p_{ij}],即得 W=VP。寫出明確的算式,

\displaystyle  \Vert \Sigma P\Vert_F^2=\sum_{i=1}^r\sigma_i^2\sum_{j=1}^kp_{ij}^2

矩陣 P 的每一行都是單位向量,\sum_{i=1}^pp_{ij}^2=1。再者,k\le p,我們可以推論 \sum_{j=1}^kp_{ij}^2\le 1。為便利說明,假設 k\le r (這與實際應用情況吻合)。上式的最大值發生於 \sum_{j=1}^kp_{ij}^2=1i=1,\ldots,k,也就是 P=\begin{bmatrix}  I_k\\  0  \end{bmatrix}。令 V=\begin{bmatrix}  V_1&V_2  \end{bmatrix},其中 V_1p\times k 階分塊,V_2p\times (p-k) 階分塊。所求的最佳基向量矩陣 W 即為 V 的領先 k 個行向量構成的分塊,W=V_1;最佳近似矩陣則為

\begin{aligned}  \hat{X}&=XWW^T=U\Sigma V^TV_1V_1^T=U\Sigma \begin{bmatrix}  V_1^T\\  V_2^T  \end{bmatrix}V_1V_1^T\\  &=U\Sigma\begin{bmatrix}  I_k\\  0  \end{bmatrix}V_1^T=U\begin{bmatrix}  D&0\\  0&0  \end{bmatrix}\begin{bmatrix}  I_k&0\\  0&0  \end{bmatrix}V^T\\  &=U\begin{bmatrix}  D'&0\\  0&0  \end{bmatrix} V^T,  \end{aligned}

其中 D'=\hbox{diag}(\sigma_1,\ldots,\sigma_k)。因為 \sigma_1\ge\cdots\ge\sigma_k>0\hbox{rank}\hat{X}=\hbox{rank}D'=k

 
最後補充一些延伸討論。

  1. 樣本共變異數矩陣的正交對角化 S=Q\Lambda Q^T 可由數據矩陣的奇異值分解 X=U\Sigma V^T 求出 (見“主成分分析與奇異值分解”)。使用定義,

    \displaystyle  S=\frac{1}{n-1}X^TX=\frac{1}{n-1}V\Sigma^T U^TU\Sigma V^T=V\left(\frac{1}{n-1}\Sigma^T\Sigma\right)V^T

    比較 S 的兩種正交對角化表達式,立得 Q=V\Lambda=\frac{1}{n-1}\hbox{diag}(\sigma_1^2,\ldots,\sigma_p^2)

  2. 定義 k\times k 階新樣本變異數矩陣

    \displaystyle  S_z=\frac{1}{n-1}\sum_{i=1}^n\mathbf{z}_i\mathbf{z}_i^T=\frac{1}{n-1}Z^TZ

    使用 Z=XW=XV_1S=Q\Lambda Q^T=V\Lambda V^T,可得

    \displaystyle  \begin{aligned}  S_z&=\frac{1}{n-1}V_1^TX^TXV_1=V_1^TSV_1=V_1^TV\Lambda V^TV_1\\  &=V_1^T\begin{bmatrix}  V_1&V_2  \end{bmatrix}\Lambda \begin{bmatrix}  V_1^T\\  V_2^T  \end{bmatrix}V_1=\begin{bmatrix}  I_k&0  \end{bmatrix}\Lambda\begin{bmatrix}  I_k\\  0  \end{bmatrix}\\  &=\begin{bmatrix}  \lambda_1&&\\  &\ddots&\\  &&\lambda_k  \end{bmatrix},  \end{aligned}

    證明任兩個新變數的相關係數等於零 (見“相關係數”),變數 z_j 的樣本變異數則為 \lambda_j=\frac{1}{n-1}\sigma_j^2j=1,\ldots,k

  3. 最佳近似數據矩陣 \hat{X} 的樣本共變異數矩陣為

    \displaystyle  \begin{aligned}  \hat{S}&=\frac{1}{n-1}\hat{X}^T\hat{X}=\frac{1}{n-1}V\begin{bmatrix}  D'&0\\  0&0  \end{bmatrix}U^TU\begin{bmatrix}  D'&0\\  0&0  \end{bmatrix}V^T\\  &=\frac{1}{n-1}V\begin{bmatrix}  \begin{matrix}  \sigma_1^2&&\\  &\ddots&\\  &&\sigma_k^2  \end{matrix}&0\\  0&0  \end{bmatrix}V^T=V\begin{bmatrix}  \begin{matrix}  \lambda_1&&\\  &\ddots&\\  &&\lambda_k  \end{matrix}&0\\  0&0  \end{bmatrix}V^T.  \end{aligned}

    主成分分析保留的總變異量比為

    \displaystyle  \frac{\hbox{trace}{S}_z}{\hbox{trace} S}=\frac{\hbox{trace}\hat{S}}{\hbox{trace} S}=\frac{\sum_{j=1}^k\lambda_j}{\sum_{j=1}^p\lambda_j}

 
註解
[1] 若樣本共變異數矩陣 S 的條件數 (condition number) 很大,譬如大於 1,000,則變數之間的共線性關係十分嚴重 (見“條件數”)。
[2] Kronecker 函數 \delta_{ij} 表示單位矩陣 I(i,j) 元。

廣告
本篇發表於 線性代數專欄, 應用之道 並標籤為 , , , , 。將永久鏈結加入書籤。

發表迴響

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

WordPress.com Logo

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

Twitter picture

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

Facebook照片

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

Google+ photo

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

連結到 %s