主成分分析與奇異值分解

本文的閱讀等級：高級

給定一份樣本大小為 $n$ 的數據 $\{\mathbf{x}_1,\ldots,\mathbf{x}_n\}$ ，其中 $\mathbf{x}_k$ 是 $p$ 維實向量，記錄 $p$ 個變數的觀測值。所有的數據點 $\mathbf{x}_k$ 扣除平均數向量 $\mathbf{m}=\frac{1}{n}\sum_{k=1}^n\mathbf{x}_k$ 可得 $n\times p$ 階離差矩陣 (deviation matrix) $X=[x_{kj}]$ ，表示如下：

$X=\begin{bmatrix} (\mathbf{x}_1-\mathbf{m})^T\\ \vdots\\ (\mathbf{x}_n-\mathbf{m})^T \end{bmatrix}$ ，

其中 $x_{kj}$ 是第 $k$ 個數據點的第 $j$ 個變數值，也就是說， $X$ 的每一列 (row) 對應一個數據點，每一行 (column) 對應一個變數。假設 $X$ 不存在常數行，即每個變數總是存在若干變異。如欲將數據予以標準化 (每一變數的平均數等於 $0$ ，變異數等於 $1$ )，將 $X$ 的每一行的所有元除以該變數的樣本標準差 (樣本變異數的平方根)，即有

$\tilde{X}=\begin{bmatrix} (\mathbf{x}_1-\mathbf{m})^T\\ \vdots\\ (\mathbf{x}_n-\mathbf{m})^T \end{bmatrix}\begin{bmatrix} 1/s_1&&\\ &\ddots&\\ &&1/s_p \end{bmatrix}$ ，

其中 $s_j^2=\frac{1}{n-1}\sum_{k=1}^n(x_{kj}-m_j)^2$ 是第 $j$ 個變數的樣本變異數， $m_j=\frac{1}{n}\sum_{k=1}^nx_{kj}$ 是第 $j$ 個變數的樣本平均數 (見“樣本平均數、變異數和共變異數”)。令 $D=\mathrm{diag}(s_1,\ldots,s_p)$ 。標準化後的離差矩陣可表示為 $\tilde{X}=XD^{-1}$ 。當數據集的變數總數 $p$ 很大或變數具有相關性時，主成分分析 (principal components analysis) 可產生一組由 $r\le p$ 個不具相關性的新變數構成的數據。主成分分析的詳細推導過程請見“主成分分析”，本文僅列舉主要結果與性質，接著介紹如何利用奇異值分解計算，最後解釋主成分分析所隱含的資料生成模型。

定義 $p\times p$ 階樣本共變異數矩陣

$\displaystyle S=\frac{1}{n-1}X^TX$ ，

和樣本相關矩陣

$\displaystyle R=\frac{1}{n-1}\tilde{X}^T\tilde{X}=\frac{1}{n-1}D^{-1}X^TXD^{-1}=D^{-1}SD^{-1}$ 。

以下討論限定為 $S$ ，只要將 $S$ 替換成 $R$ ，便可得到相應結果。實對稱矩陣 $S$ 可正交對角化為

$S=W\Lambda W^T$ ，

其中 $W=\begin{bmatrix} \mathbf{w}_1&\cdots&\mathbf{w}_p \end{bmatrix}$ 是由單範正交 (orthonormal) 主成分 (特徵向量) 構成的 $p\times p$ 階正交矩陣，滿足 $W^TW=WW^T=I$ ； $\Lambda=\mathrm{diag}(\lambda_1,\ldots,\lambda_p)$ 是特徵值矩陣， $\lambda_1\ge\cdots\ge\lambda_p\ge 0$ 。令 $z_{kj}$ 是數據點 $\mathbf{x}_k-\mathbf{m}$ 所含主成分 $\mathbf{w}_j$ 的組合係數，稱為主成分係數，可由正交投影算得：

$z_{kj}=(\mathbf{x}_k-\mathbf{m})^T\mathbf{w}_j,~~1\le k\le n,~~1\le j\le p$ ，

或合併成 $n\times p$ 階主成分係數矩陣

$Z=\begin{bmatrix} (\mathbf{x}_1-\mathbf{m})^T\\ \vdots\\ (\mathbf{x}_n-\mathbf{m})^T \end{bmatrix}\begin{bmatrix} \mathbf{w}_1&\cdots&\mathbf{w}_p \end{bmatrix}=XW$ 。

特徵值 $\lambda_j$ 是第 $j$ 個主成分係數 (即 $Z$ 的第 $j$ 行) 的樣本變異數，因此代表主成分 $\mathbf{w}_j$ 的權值。

在數值計算上，矩陣乘法運算可能因引進捨入或截斷誤差而造成破壞性的影響。例如，

$A=\begin{bmatrix} 1&1&1\\ \epsilon&0&0\\ 0&\epsilon&0\\ 0&0&\epsilon \end{bmatrix}$ ，

其中 $\epsilon$ 為極小數。若數值計算誤差使得 $1+\epsilon^2\approx 1$ ，交互乘積

$A^TA=\begin{bmatrix} 1+\epsilon^2&1&1\\ 1&1+\epsilon^2&1\\ 1&1&1+\epsilon^2 \end{bmatrix}$

變成秩-1矩陣，違反 $\mathrm{rank}(A^TA)=\mathrm{rank}A$ 。因為這個緣故， $\frac{1}{n-1}X^TX$ 的正交譜分解 (對角化) 僅可用於理論推導，不宜用於數值計算。奇異值分解提供了一個不涉及交互乘積的解決辦法。考慮離差矩陣 $X$ 的「瘦」奇異值分解 (見“奇異值分解 (SVD)”)

$X=U\Sigma V^T$ ，

其中 $n\times p$ 階矩陣 $U$ 的行向量是單範正交左奇異向量 $\{\mathbf{u}_1,\ldots,\mathbf{u}_p\}$ ， $U^TU=I_p$ ，但 $UU^T$ 不一定等於 $I_n$ ； $\Sigma=\mathrm{diag}(\sigma_1,\ldots,\sigma_p)$ ， $\sigma_1\ge\cdots\ge\sigma_p\ge 0$ 是奇異值； $p\times p$ 階矩陣 $V$ 的行向量是單範正交右奇異向量 $\{\mathbf{v}_1,\ldots,\mathbf{v}_p\}$ ， $VV^T=V^TV=I_p$ 。提醒讀者， $X$ 和 $\tilde{X}$ 的奇異值分解並不相同。將奇異值分解代入交互乘積，

$X^TX=(U\Sigma V^T)^T(U\Sigma V^T)=V\Sigma^TU^TU\Sigma V^T=V\Sigma^2V^T$ ，

故樣本共變異數矩陣可表示為

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

比較 $S$ 的兩種正交對角化表達式，可得特徵值矩陣

$\displaystyle \Lambda=\frac{1}{n-1}\Sigma^2$ ，

主成分矩陣

$W=V$ ，

以及主成分係數矩陣

$Z=XV=(U\Sigma V^T)V=U\Sigma$ 。

奇異值分解的正交矩陣 $V$ 即為主成分矩陣 $W$ ，奇異值 $\sigma_j$ 決定特徵值 $\lambda_j=\frac{1}{n-1}\sigma_j^2$ ，矩陣 $U=[u_{kj}]$ 和奇異值 $\sigma_j$ 共同決定主成分係數 $z_{kj}=u_{kj}\sigma_j$ ，那麼 $U$ 本身具有甚麼意義呢？若 $X$ 有線性獨立的行向量，即 $\mathrm{rank}X=p$ ，則 $\sigma_1\ge\cdots\ge\sigma_p>0$ ， $\Sigma=\sqrt{n-1}\Lambda^{1/2}$ 可逆，其中 $\Lambda^{1/2}=\mathrm{diag}(\sqrt{\lambda_1},\ldots,\sqrt{\lambda_p})$ 。等式 $Z=U\Sigma$ 右乘 $\Sigma^{-1}$ ，可得

$\displaystyle\begin{aligned} U&=Z\Sigma^{-1}=\frac{1}{\sqrt{n-1}}Z\Lambda^{-1/2}\\ &=\frac{1}{\sqrt{n-1}}\begin{bmatrix} z_{11}&z_{12}&\cdots&z_{1p}\\ z_{21}&z_{22}&\cdots&z_{2P}\\ \vdots&\vdots&\ddots&\vdots\\ z_{n-1}&z_{n2}&\cdots&z_{np} \end{bmatrix}\begin{bmatrix} \frac{1}{\sqrt{\lambda_1}}&&&\\ &\frac{1}{\sqrt{\lambda_2}}&&\\ &&\ddots&\\ &&&\frac{1}{\sqrt{\lambda_p}} \end{bmatrix}\\ &=\frac{1}{\sqrt{n-1}}\begin{bmatrix} \frac{z_{11}}{\sqrt{\lambda_1}}&\frac{z_{12}}{\sqrt{\lambda_2}}&\cdots&\frac{z_{1p}}{\sqrt{\lambda_p}}\\ \frac{z_{21}}{\sqrt{\lambda_1}}&\frac{z_{22}}{\sqrt{\lambda_2}}&\cdots&\frac{z_{2p}}{\sqrt{\lambda_p}}\\ \vdots&\vdots&\ddots&\vdots\\ \frac{z_{n1}}{\sqrt{\lambda_1}}&\frac{z_{n2}}{\sqrt{\lambda_2}}&\cdots&\frac{z_{np}}{\sqrt{\lambda_p}} \end{bmatrix}.\end{aligned}$

因為第 $j$ 個主成分係數的樣本變異數等於 $\lambda_j$ ，

$\tilde{Z}=Z\Lambda^{-1/2}=\sqrt{n-1}U$

就是 $Z$ 經過標準化的結果，稱為因素分數 (factor score) 矩陣，它的第 $j$ 行 $\sqrt{n-1}\mathbf{u}_j$ 表示數據集在第 $j$ 個主成分的得點，其平均數等於 $0$ ，變異數等於 $1$ 。

另外，我們想知道抽取出的主成分係數與原變數之間的相關性。令 $f_{ij}$ 是第 $i$ 個變數與第 $j$ 個主成分係數的相關係數 (見“相關係數”)，稱為因素負荷 (factor loading)，並令 $F=[f_{ij}]$ 為 $p\times p$ 階因素負荷矩陣。因為 $\tilde{X}$ 和 $\tilde{Z}$ 包含標準化後的數據，立知

$\displaystyle F=\frac{1}{n-1}\tilde{X}^T\tilde{Z}$ 。

將 $\tilde{X}=XD^{-1}$ ， $\tilde{Z}=XV\Lambda^{-1/2}$ ， $S=\frac{1}{n-1}X^TX=V\Lambda V^T$ ， $\Lambda^{1/2}=\frac{1}{\sqrt{n-1}}\Sigma$ 代入上式，可得

$\displaystyle\begin{aligned} F&=\frac{1}{n-1}(D^{-1}X^T)(XV\Lambda^{-1/2})=D^{-1}SV\Lambda^{-1/2}\\ &=D^{-1}(V\Lambda V^T)V\Lambda^{-1/2}=D^{-1}V\Lambda\Lambda^{-1/2}\\ &=D^{-1}V\Lambda^{1/2}=\frac{1}{\sqrt{n-1}}D^{-1}V\Sigma. \end{aligned}$

如果以 $R$ 取代 $S$ ，因素負荷矩陣則為 $F=\frac{1}{\sqrt{n-1}}V\Sigma$ ，但 $V$ 和 $\Sigma$ 也隨之改變，因為 $\tilde{X}$ 和 $X$ 有不同的奇異值分解。

主成分分析提供了一個資料生成的描述模型。利用奇異值分解，離差矩陣 $X$ 可用因素分數 $\tilde{Z}$ 與因素負荷 $F$ 表達如下：

$\displaystyle X=U\Sigma V^T=U\Sigma V^TD^{-1}D=\sqrt{n-1}U\left(\frac{1}{\sqrt{n-1}}D^{-1}V\Sigma\right)^TD=\tilde{Z}F^TD$ 。

取轉置可得

$\displaystyle X^T=\begin{bmatrix} \mathbf{x}_1-\mathbf{m}&\cdots&\mathbf{x}_n-\mathbf{m} \end{bmatrix}=DF\tilde{Z}^T=DF\begin{bmatrix} \tilde{\mathbf{z}}_1&\cdots&\tilde{\mathbf{z}}_n \end{bmatrix}$ 。

上式可以這麼解讀：想像存在 $p$ 個獨立 (無相關) 且強度相同的隨機訊號源，稱為因素 (factor)。每一因素產生平均數是 $0$ ，變異數是 $1$ 的訊號。因素分數 (轉置) 矩陣 $\tilde{Z}^T$ 儲存了對於 $p$ 個因素的 $n$ 次觀測結果，其中向量 $\tilde{\mathbf{z}}_k=(\tilde{z}_{k1},\ldots,\tilde{z}_{kp})^T$ 記錄 $p$ 個因素的第 $k$ 次觀測值。訊號 $\tilde{\mathbf{z}}_k$ 先經過一個叫做因素負荷 $F$ 的混合過程，其中 $f_{ij}$ 決定從第 $j$ 個因素傳遞至第 $i$ 個變數的載入權重 (即相關係數)；接著 $p$ 個變數通過伸縮過程 $D$ 以放大或縮小其數值範圍 (即標準差)；最後再平移 $\mathbf{m}$ ，產生我們觀測到的數據點 $\mathbf{x}_k$ (見下圖，以 $p=4$ 為例)。主成分分析的作用即在從接收的數據集 $X$ 反推資料生成模型中獨立訊號的大小 (即因素分數 $\tilde{Z}$ )，並揭示各個變數所含的混合訊號權重 (即因素負荷 $F$ )。

主成分分析的資料生成模型

4 Responses to 主成分分析與奇異值分解

w says:

09/11/2017 at 11:40 am

老师您好，请问 UU^T 不一定等於 I_n

- ccjou says:
  
  09/11/2017 at 12:06 pm
  
  舉個例子可以幫助理解，比如
  $U=\begin{bmatrix} 0&1\\ 1&0\\ 0&0 \end{bmatrix}$ 。
  
  - w says:
    
    09/11/2017 at 1:37 pm
    
    SVD分解中的U不是m*m的吗？
    
    - w says:
      
      09/11/2017 at 1:50 pm
      
      哦谢谢老师没看仔细一个是In 一个是Ip

	陳倍恩 on 線性代數的第一堂課──矩陣乘法的定義
	輕鬆談如何教學二項式定理？… on 牛頓的二項式定理 (上)
	madhouse on 高斯消去法
	WishMobile on 翻轉 LU 分解
	周子傑 on Karush-Kuhn-Tucker (KKT) 條件
	Cloud Huang on 線性泛函與伴隨