主成分分析

Posted on 04/15/2013 by ccjou

本文的閱讀等級：高級

美國作家梭羅 (Henry D. Thoreau) 在《湖濱散記》談到他的幽居生活時，說道^[1]：

我們的生活消耗在瑣碎之中。一個老實的人除了十指之外，便不必有更大的數字了，頂多加上十個足趾，其餘不妨勉強一下。簡單，簡單，簡單啊！我說，最好你的事祇兩三件，不要一百件或一千件；不必一百萬一百萬地計算，半打不夠計算嗎？總之，賬目可以記在大拇指甲上就好了。

我們也許不能複製梭羅在瓦爾登湖 (Walden) 的簡單生活，但是我們永遠可以通過化繁為簡來改善現況。處於資訊爆炸的時代，我們不免要面對變數很多且樣本數很大的資料。在分析高維度 (變數很多) 數據時，降維 (dimension reduction) 常是一個必要的前處理工作。主成分分析 (principal components analysis，簡稱 PCA) 由英國統計學家皮爾生 (Karl Pearson) 於1901年提出^[2]，是一種降低數據維度的有效技術。主成分分析的主要構想是分析共變異數矩陣 (covariance matrix) 的特徵性質 (見“共變異數矩陣與常態分布”)，以得出數據的主成分 (即特徵向量) 與它們的權值 (即特徵值)；透過保留低階主成分 (對應大特徵值)，捨棄高階主成分 (對應小特徵值)，達到減少數據集維度，同時保留最大數據集變異的目的。本文從線性代數觀點介紹主成分分析，並討論實際應用時可能遭遇的一些問題。

假設我們有一筆維數等於 $p$ ，樣本大小是 $n$ 的數據 $\{\mathbf{x}_1,\ldots,\mathbf{x}_n\}$ ，即每一數據點 $\mathbf{x}_k\in\mathbb{R}^p$ 包含 $p$ 個變數的量測值。如果我們想要以單一向量 $\mathbf{a}$ 來代表整組數據，可用均方誤差作為誤差函數：

$\displaystyle E_0(\mathbf{a})=\frac{1}{n-1}\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{a}\Vert^2$ 。

上式中，誤差平方和除以 $n-1$ 而非 $n$ 的緣由請見“樣本平均數、變異數和共變異數”。理想的中心向量 $\mathbf{a}$ 必須有最小的均方誤差，滿足此條件的向量是

$\displaystyle \mathbf{m}=\frac{1}{n}\sum_{k=1}^n\mathbf{x}_k$ ，

稱為樣本平均數向量。證明如下：

$\displaystyle\begin{aligned} E_0(\mathbf{a})&=\frac{1}{n-1}\sum_{k=1}^n\Vert(\mathbf{x}_k-\mathbf{m})+(\mathbf{m}-\mathbf{a})\Vert^2\\ &=\frac{1}{n-1}\left(\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2+\sum_{k=1}^n\Vert\mathbf{m}-\mathbf{a}\Vert^2+2\sum_{k=1}^n(\mathbf{x}_k-\mathbf{m})^T(\mathbf{m}-\mathbf{a})\right)\\ &=\frac{1}{n-1}\left(\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2+n\Vert\mathbf{m}-\mathbf{a}\Vert^2+2\left(\sum_{k=1}^n\mathbf{x}_k-n\mathbf{m}\right)^T(\mathbf{m}-\mathbf{a})\right). \end{aligned}$

根據樣本平均數向量 $\mathbf{m}$ 的表達式，上式最後一項等於零。因為 $\Vert\mathbf{m}-\mathbf{a}\Vert^2\ge 0$ ，可知 $E_0(\mathbf{a})\ge\frac{1}{n-1}\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2$ ，等號於 $\mathbf{a}=\mathbf{m}$ 時成立。

數據點 $\mathbf{x}_k$ 與其分布平均的差， $\mathbf{x}_k-\mathbf{m}$ ，稱為離差 (deviation)。在變數數目 $p$ 很大時，如何「目視」離差的散布情況？最簡單的作法是降維，譬如，將散布於高維空間 $\mathbb{R}^p$ 的離差「壓縮」至一直線上。令直線 $L$ 穿越 $\mathbf{m}$ ， $\mathbf{w}$ 為其指向向量，且 $\Vert\mathbf{w}\Vert=1$ 。直線 $L$ 上的任一點可表示如下：

$\mathbf{x}=\mathbf{m}+c\mathbf{w}$ ，

其中 $c$ 為一純量， $\vert c\vert$ 代表 $\mathbf{x}$ 至 $\mathbf{m}$ 的距離。一旦 $\mathbf{w}$ 給定，我們可以用 $\mathbf{m}+c_k\mathbf{w}$ 來近似 $\mathbf{x}_k$ ，也就是以 $c_k\mathbf{w}$ 近似離差 $\mathbf{x}_k-\mathbf{m}$ 。如同樣本平均數向量的設定方式，最佳的近似係數 $c_1,\ldots,c_n$ 必須最小化均方誤差：

$\displaystyle \begin{aligned} E_1(\{c_k\},\mathbf{w})&=\frac{1}{n-1}\sum_{k=1}^n\Vert(\mathbf{m}+c_k\mathbf{w})-\mathbf{x}_k\Vert^2\\ &=\frac{1}{n-1}\sum_{k=1}^n\Vert c_k\mathbf{w}-(\mathbf{x}_k-\mathbf{m}) \Vert^2\\ &=\frac{1}{n-1}\left(\sum_{k=1}^nc_k^2\Vert\mathbf{w}\Vert^2-2\sum_{k=1}^nc_k\mathbf{w}^T(\mathbf{x}_k-\mathbf{m})+\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2\right). \end{aligned}$

因為 $\Vert\mathbf{w}\Vert=1$ ，計算偏導數 $\partial E_1/\partial c_k$ 並設為零，

$\displaystyle \frac{\partial E_1}{\partial c_k}=\frac{1}{n-1}\left(2c_k-2\mathbf{w}^T(\mathbf{x}_k-\mathbf{m})\right)=0$ 。

由此解得

$c_k=\mathbf{w}^T(\mathbf{x}_k-\mathbf{m}),~~k=1,\ldots,n$ 。

從幾何觀點解釋， $c_k\mathbf{w}$ 即是離差 $\mathbf{x}_k-\mathbf{m}$ 在直線 $L$ 的正交投影，以下稱為正交原則 (見圖1，參閱“正交投影──威力強大的線代工具”)。

圖1 正交投影即為最佳近似

接下來我們尋找使 $E_1$ 最小化的直線方向 $\mathbf{w}$ 。將先前解出的最佳係數代回 $E_1$ ，整理可得

$\displaystyle\begin{aligned} E_1(\mathbf{w})&=\frac{1}{n-1}\left(\sum_{k=1}^nc_k^2-2\sum_{k=1}^nc_k^2+\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2\right)\\ &=-\frac{1}{n-1}\sum_{k=1}^n\left(\mathbf{w}^T(\mathbf{x}_k-\mathbf{m})\right)^2+\frac{1}{n-1}\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2\\ &=-\frac{1}{n-1}\sum_{k=1}^n\mathbf{w}^T(\mathbf{x}_k-\mathbf{m})(\mathbf{x}_k-\mathbf{m})^T\mathbf{w}+\frac{1}{n-1}\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2\\ &=-\mathbf{w}^T\left(\frac{1}{n-1}\sum_{k=1}^n(\mathbf{x}_k-\mathbf{m})(\mathbf{x}_k-\mathbf{m})^T\right)\mathbf{w}+\frac{1}{n-1}\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2, \end{aligned}$

上面使用了 $\mathbf{w}^T(\mathbf{x}_k-\mathbf{m})=(\mathbf{x}_k-\mathbf{m})^T\mathbf{w}$ 。令

$\displaystyle S=\frac{1}{n-1}\sum_{k=1}^n(\mathbf{x}_k-\mathbf{m})(\mathbf{x}_k-\mathbf{m})^T$ ，

稱為樣本共變異數矩陣。乘開上式，確認 $S$ 的 $(i,j)$ 元 $s_{ij}$ 即為第 $i$ 個變數和第 $j$ 個變數的樣本共變異數：

$\displaystyle s_{ij}=\frac{1}{n-1}\sum_{k=1}^n(x_{ki}-m_i)(x_{kj}-m_j)$ ，

其中 $x_{ki}$ 是第 $k$ 個數據點的第 $i$ 個變數 (即 $\mathbf{x}_k$ 的第 $i$ 元)， $m_i=\frac{1}{n}\sum_{k=1}^nx_{ki}$ 是第 $i$ 個變數的樣本平均數 (即 $\mathbf{m}$ 的第 $i$ 元)。明顯地， $s_{ij}=s_{ji}$ ，故 $S$ 是一 $p\times p$ 階對稱矩陣。對於任一 $p$ 維向量 $\mathbf{y}$ ，

$\displaystyle \mathbf{y}^TS\mathbf{y}=\frac{1}{n-1}\mathbf{y}^T\sum_{k=1}^n(\mathbf{x}_k-\mathbf{m})(\mathbf{x}_k-\mathbf{m})^T\mathbf{y}=\frac{1}{n-1}\sum_{k=1}^n\left(\mathbf{y}^T(\mathbf{x}_k-\mathbf{m})\right)^2\ge 0$ ，

即知 $S$ 是半正定矩陣。數據集的總變異量 (即離差平方和) $\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2$ 是一常數，所以最小化 $E_1(\mathbf{w})$ 等價於最大化 $\mathbf{w}^TS\mathbf{w}$ 。我們的問題變成求解下列約束二次型：

$\displaystyle \max_{\Vert\mathbf{w}\Vert=1}\mathbf{w}^TS\mathbf{w}$ 。

解法如下 (另一個解法見“Hermitian 矩陣特徵值的變化界定”)。使用 Lagrange 乘數法 (見“Lagrange 乘數法”)，定義

$L(\mathbf{w},\mu)=\mathbf{w}^TS\mathbf{w}-\mu(\mathbf{w}^T\mathbf{w}-1)$ ，

其中 $\mu$ 是未定的 Lagrange 乘數。計算偏導數 (即梯度) 並設為零，

$\displaystyle \frac{\partial L}{\partial \mathbf{w}}=2S\mathbf{w}-2\mu\mathbf{w}=\mathbf{0}$ ，

可得

$S\mathbf{w}=\mu\mathbf{w}$ 。

直線 $L$ 的指向向量 $\mathbf{w}$ 就是樣本共變異數矩陣 $S$ 的特徵向量。因為 $\mathbf{w}^TS\mathbf{w}=\mu\mathbf{w}^T\mathbf{w}=\mu$ ，欲使 $\mathbf{w}^TS\mathbf{w}$ 有最大值，我們必須選擇對應最大特徵值的特徵向量。實對稱半正定矩陣有非負的特徵值 (見“半正定矩陣的判別方法”)，故可設 $S$ 的特徵值為 $\lambda_1\ge\lambda_2\cdots\ge\lambda_p\ge 0$ 。此外， $S$ 的特徵向量構成一個單範正交集 (orthonormal set)，見“實對稱矩陣可正交對角化的證明”，這個性質將於稍後使用。我們選擇對應最大特徵值 $\lambda_1$ 的特徵向量作為 $\mathbf{w}$ ，透過離差 $\mathbf{x}_k-\mathbf{m}$ 在直線 $L$ 的正交投影 $c_k=\mathbf{w}^T(\mathbf{x}_k-\mathbf{m})$ 即可粗估數據集的散布情況。

如果我們想獲得較為精確的近似結果，以上過程可推廣至更高維度。考慮

$\mathbf{x}=\mathbf{m}+z_{1}\mathbf{w}_1+\cdots+z_r\mathbf{w}_r$ ，

其中 $r\le p$ 。我們的目標是尋找一組單範正交集 $\{\mathbf{w}_1,\ldots,\mathbf{w}_r\}$ ，即 $\mathbf{w}_i^T\mathbf{w}_j=\delta_{ij}$ (稱為 Kronecker 記號， $\delta_{ij}=1$ 若 $i=j$ ，否則 $\delta_{ij}=0$ )，使最小化

$\displaystyle E_r\left(\{\mathbf{w}_j\}\right)=\sum_{k=1}^n\left\|\left(\mathbf{m}+\sum_{j=1}^rz_{kj}\mathbf{w}_j\right)-\mathbf{x}_k\right\|^2$ 。

根據正交原則，組合係數 $z_{kj}$ 就是離差 $\mathbf{x}_k-\mathbf{m}$ 至單位向量 $\mathbf{w}_j$ 的正交投影量，即 $z_{kj}=\mathbf{w}_j^T(\mathbf{x}_k-\mathbf{m})$ 。重複前述推導方式，可以證明^[3]

$\displaystyle E_r\left(\{\mathbf{w}_j\}\right)=-(n-1)\sum_{j=1}^r\mathbf{w}_j^TS\mathbf{w}_j+\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2$ 。

緊接的任務要解出約束二次型：

$\displaystyle \max_{\mathbf{w}_i^T\mathbf{w}_j=\delta_{ij}}\sum_{j=1}^r\mathbf{w}_j^TS\mathbf{w}_j$ 。

仍然使用 Lagrange 乘數法，定義

$\displaystyle L\left(\{\mathbf{w}_j\},\{\mu_{ij}\}\right)=\sum_{j=1}^r\mathbf{w}_j^TS\mathbf{w}_j-\sum_{i=1}^r\sum_{j=1}^r\mu_{ij}(\mathbf{w}_i^T\mathbf{w}_j-\delta_{ij})$ ，

其中 $\{\mu_{ij}\}$ 是未定的 Lagrange 乘數。計算偏導數，並設結果為零向量：

$\displaystyle \frac{\partial L}{\partial \mathbf{w}_j}=2S\mathbf{w}_j-2\mu_{jj}\mathbf{w}_j-\sum_{i\neq j}(\mu_{ij}+\mu_{ji})\mathbf{w}_i=\mathbf{0},~~j=1,\ldots,r$ 。

對於 $i\neq j$ ，設 $\mu_{ij}+\mu_{ji}=0$ ，則

$S\mathbf{w}_j=\mu_{jj}\mathbf{w}_j,~~j=1,\ldots,r$ 。

換句話說，當 $\mathbf{w}_1,\ldots,\mathbf{w}_r$ 是樣本共變數矩陣 $S$ 的最大 $r$ 個特徵值 $\lambda_1,\ldots,\lambda_r$ 的對應 (單範正交) 特徵向量時，目標函數 $\sum_{j=1}^r\mathbf{w}_j^TS\mathbf{w}_j=\sum_{j=1}^r\lambda_j\mathbf{w}_j^T\mathbf{w}_j=\sum_{j=1}^r\lambda_j$ 有最大值。當然，這並不是唯一解，但可以證明所有的解都有相同的目標函數值 $\sum_{j=1}^r\mathbf{w}_j^TS\mathbf{w}_j$ ^[4]。不過，如果要求 $\{\mathbf{w}_1,\ldots,\mathbf{w}_r\}$ 滿足 $\max_{\mathbf{w}_i^T\mathbf{w}_j=\delta_{ij}}\sum_{j=1}^r\mathbf{w}_j^TS\mathbf{w}_j$ ， $1\le r\le p$ ，則樣本共變數矩陣 $S$ 的特徵向量集是唯一的選擇^[5]。

我們以 $\mathbf{m}+z_{k1}\mathbf{w}_1+\cdots+z_{kr}\mathbf{w}_r$ 近似數據點 $\mathbf{x}_k$ ，樣本共變異數矩陣 $S$ 的特徵向量 $\{\mathbf{w}_1,\ldots,\mathbf{w}_r\}$ 描述了數據集的主要成分，因此稱為主成分。那麼 $S$ 的特徵值又有甚麼涵義呢？特徵值 $\lambda_j$ 就是主成分 $\mathbf{w}_j$ 的係數 $z_j$ (稱為主成分係數) 的變異數。首先證明 $z_j$ 的平均數是零：對於 $j=1,\ldots,r$ ，

$\displaystyle \frac{1}{n}\sum_{k=1}^nz_{kj}=\frac{1}{n}\sum_{k=1}^n\mathbf{w}_j^T(\mathbf{x}_k-\mathbf{m})=\frac{1}{n}\mathbf{w}_j^T\left(\sum_{k=1}^n\mathbf{x}_k-n\mathbf{m}\right)=0$ 。

再來計算 $z_j$ 的樣本變異數：

$\displaystyle\begin{aligned} s^2_{z_j}&=\frac{1}{n-1}\sum_{k=1}^nz_{kj}^2=\frac{1}{n-1}\sum_{k=1}^n(\mathbf{w}_j^T(\mathbf{x}_k-\mathbf{m}))((\mathbf{x}_k-\mathbf{m})^T\mathbf{w}_j)\\ &=\mathbf{w}_j^T\left(\frac{1}{n-1}\sum_{k=1}^T(\mathbf{x}_k-\mathbf{m})(\mathbf{x}_k-\mathbf{m})^T\right)\mathbf{w}_j=\mathbf{w}_j^TS\mathbf{w}_j=\lambda_j\mathbf{w}_j^T\mathbf{w}_j=\lambda_j. \end{aligned}$

所以說，特徵值 $\lambda_j$ 表示主成分 $\mathbf{w}_j$ 的權值。另外必須一提的是，主成分係數 $z_i$ 和 $z_j$ 沒有相關性，即樣本共變異數 $s_{z_iz_j}$ 為零。當 $i\neq j$ ，

$\displaystyle\begin{aligned} s_{z_iz_j}&=\frac{1}{n-1}\sum_{k=1}z_{ki}z_{kj}=\frac{1}{n-1}\sum_{k=1}^n(\mathbf{w}_i^T(\mathbf{x}_k-\mathbf{m}))((\mathbf{x}_k-\mathbf{m})^T\mathbf{w}_j)\\ &=\mathbf{w}_i^TS\mathbf{w}_j=\lambda_j\mathbf{w}_i^T\mathbf{w}_j=0.\end{aligned}$

給定維度等於 $p$ 的數據集 $\{\mathbf{x}_1,\ldots,\mathbf{x}_n\}$ ，主成分分析的計算程序與結果可整理成簡明的矩陣形式。

計算樣本平均 $\mathbf{m}=\frac{1}{n}\sum_{k=1}^n\mathbf{x}_k$ ，定義 $n\times p$ 階離差矩陣

$X=\begin{bmatrix} (\mathbf{x}_1-\mathbf{m})^T\\ (\mathbf{x}_2-\mathbf{m})^T\\ \vdots\\ (\mathbf{x}_n-\mathbf{m})^T \end{bmatrix}=\begin{bmatrix} x_{11}-m_1&x_{12}-m_2&\cdots&x_{1p}-m_p\\ x_{21}-m_1&x_{22}-m_2&\cdots&x_{2p}-m_p\\ \vdots&\vdots&\ddots&\vdots\\ x_{n1}-m_1&x_{n2}-m_2&\cdots&x_{np}-m_p \end{bmatrix}$ ，

$p\times p$ 階樣本共變異數矩陣則是

$\displaystyle S=\frac{1}{n-1}\sum_{k=1}^n(\mathbf{x}_k-\mathbf{m})(\mathbf{x}_k-\mathbf{m})^T=\frac{1}{n-1}X^TX$ 。
將 $S$ 正交對角化為

$S=W\Lambda W^T$ ，

其中 $\Lambda=\mathrm{diag}(\lambda_1,\ldots,\lambda_p)$ 是特徵值矩陣， $\lambda_1\ge\cdots\ge\lambda_p\ge 0$ 代表主成分的權值， $W=\begin{bmatrix} \mathbf{w}_1&\cdots&\mathbf{w}_p \end{bmatrix}$ 是單範正交特徵向量構成的 $p\times p$ 階正交主成分矩陣， $W^TW=WW^T=I_p$ 。圖2顯示 $p=2$ 的資料散布圖，樣本平均數向量 $\mathbf{m}$ ，以及主成分 $\mathbf{w}_1$ 和 $\mathbf{w}_2$ 。圖中橢圓的長軸平方與短軸平方之比等於主成分係數 $z_1$ 的變異數與 $z_2$ 的變異數之比，即 $\lambda_1:\lambda_2$ 。

圖2 資料散布圖與主成分
定義 $n\times p$ 階主成分係數矩陣 $Z=[z_{kj}]$ ，其中 $z_{kj}=(\mathbf{x}_k-\mathbf{m})^T\mathbf{w}_j$ ，因此

$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$ 。

上式等號兩邊右乘 $W^T$ ，可得 $X=ZW^T$ 。換一個說法，數據點 $\mathbf{x}_k$ 的主成分分解式為

$\displaystyle \mathbf{x}_k=\mathbf{m}+\sum_{j=1}^pz_{kj}\mathbf{w}_j,~~k=1,\ldots,n$ 。

主成分係數 $(z_{k1},\ldots,z_{kp})$ 是離差 $\mathbf{x}_k-\mathbf{m}$ 參考單範正交基底 $\mathfrak{B}=\{\mathbf{w}_1,\ldots,\mathbf{w}_p\}$ 的座標向量。

最後還有一些與實際應用面相關的細節需要釐清。

我們應當保留多少低階主成分 (對應大特徵值的特徵向量)？也就是說，如何選擇 $r$ ？常用的一種方式是設定近似數據 $\mathbf{m}+\sum_{j=1}^rz_{kj}\mathbf{w}_j$ 的變異與原始數據 $\mathbf{x}_k=\mathbf{m}+\sum_{j=1}^pz_{kj}\mathbf{w}_j$ 的變異的比例。譬如，選擇最小的 $r$ 使得

$\displaystyle \frac{\sum_{j=1}^r\lambda_j}{\sum_{j=1}^p\lambda_j}\ge 0.8$ ，

表示我們保留了 80% 的數據變異。
若數據的變數具有不同的變異，主成分方向會受到變異大的變數所決定。如欲排除這個影響，我們可以用樣本相關矩陣取代樣本共變異數矩陣。在套用主成分分析之前，預先將每一變數予以標準化 (standardized)，如下：

$\displaystyle \begin{aligned} \tilde{X}&=\begin{bmatrix} (\mathbf{x}_1-\mathbf{m})^T\\ (\mathbf{x}_2-\mathbf{m})^T\\ \vdots\\ (\mathbf{x}_n-\mathbf{m})^T \end{bmatrix}\begin{bmatrix} 1/s_1&&&\\ &1/s_2&&\\ &&\ddots&\\ &&&1/s_p \end{bmatrix}\\ &=\begin{bmatrix} (x_{11}-m_1)/s_1&(x_{12}-m_2)/s_2&\cdots&(x_{1p}-m_p)/s_p\\ (x_{21}-m_1)/s_1&(x_{22}-m_2)/s_2&\cdots&(x_{2p}-m_p)/s_p\\ \vdots&\vdots&\ddots&\vdots\\ (x_{n1}-m_1)/s_1&(x_{n2}-m_2)/s_2&\cdots&(x_{np}-m_p)/s_p \end{bmatrix},\end{aligned}$

其中 $s_i^2$ 是第 $i$ 個變數的樣本變異數，即 $s_i^2=\frac{1}{n-1}\sum_{k=1}^n(x_{ki}-m_i)^2$ 。標準化後的數據集的樣本共變異數矩陣即為樣本相關矩陣

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

請讀者自行驗證 $R$ 的 $(i,j)$ 元就是第 $i$ 個變數與第 $j$ 個變數的相關係數 (見“相關係數”)。這時候，數據集的總變異等於維數 $p$ ，原因如下：

$\displaystyle \sum_{j=1}^p\lambda_j=\mathrm{trace}\Lambda=\mathrm{trace}(W^TRW)=\mathrm{trace}(RWW^T)=\mathrm{trace}R=p$ 。

上面使用了跡數循環不變性 $\mathrm{trace}(AB)=\mathrm{trace}(BA)$ (見“跡數的性質與應用”)，最後一個等式係因 $R$ 的主對角元皆為 $1$ 。
如何得到數值穩定的主成分 $\mathbf{w}_1,\ldots,\mathbf{w}_p$ ，權值 $\lambda_1,\ldots,\lambda_p$ ，以及主成分係數 $z_{kj}$ ， $k=1,\ldots,n$ ， $j=1,\ldots,p$ ？答案是奇異值分解 (singular value decomposition)。通過主成分分析與奇異值分解的關係可以顯現主成分分析隱含的其他訊息 (見“主成分分析與奇異值分解”)。

註解
[1] Henry D. Thoreau 的 Walden，原文如下：“Our life is frittered away by detail. An honest man has hardly need to count more than his ten fingers, or in extreme cases he may add his ten toes, and lump the rest. Simplicity, simplicity, simplicity! I say, let your affairs be as two or three, and not a hundred or a thousand; instead of a million count half a dozen, and keep your accounts on your thumb-nail.” 中譯取自《湖濱散記》，吳明實譯，今日世界出版社，1978年，83頁。

[2] 維基百科：主成分分析

[3]

$\displaystyle \begin{aligned} E_r\left(\{\mathbf{w}_j\}\right)&=\sum_{k=1}^n\left\|\left(\mathbf{m}+\sum_{j=1}^rz_{kj}\mathbf{w}_j\right)-\mathbf{x}_k\right\|^2\\ &=\sum_{k=1}^n\left\|\sum_{j=1}^rz_{kj}\mathbf{w}_j-(\mathbf{x}_k-\mathbf{m})\right\|^2\\ &=\sum_{k=1}^n\left\|\sum_{j=1}^rz_{kj}\mathbf{w}_j\right\|^2-2\sum_{k=1}^n\sum_{j=1}^rz_{kj}\mathbf{w}_j^T(\mathbf{x}_k-\mathbf{m})+\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2\\ &=\sum_{k=1}^n\sum_{j=1}^rz_{kj}^2-2\sum_{k=1}^n\sum_{j=1}^rz_{kj}^2+\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2\\ &=-\sum_{k=1}^n\sum_{j=1}^r\left(\mathbf{w}_j^T(\mathbf{x}_k-\mathbf{m})\right)^2+\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2\\ &=-\sum_{j=1}^r\mathbf{w}_j^T\left(\sum_{k=1}^n(\mathbf{x}_k-\mathbf{m})(\mathbf{x}_k-\mathbf{m})^T\right)\mathbf{w}_j+\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2\\ &=-(n-1)\sum_{j=1}^r\mathbf{w}_j^TS\mathbf{w}_j+\sum_{k=1}^n\Vert\mathbf{x}_k-\mathbf{m}\Vert^2 \end{aligned}$

[4] 考慮

$\displaystyle S\mathbf{w}_j=\mu_{jj}\mathbf{w}_j+\frac{1}{2}\sum_{i\neq j}(\mu_{ij}+\mu_{ji})\mathbf{w}_i,~~j=1,\ldots,r$ 。

令 $W=\begin{bmatrix} \mathbf{w}_1&\cdots&\mathbf{w}_r \end{bmatrix}$ 且

$\displaystyle M=\begin{bmatrix} \mu_{11}&\cdots&\frac{\mu_{1r}+\mu_{r1}}{2}\\ \vdots&\ddots&\vdots\\ \frac{\mu_{1r}+\mu_{r1}}{2}&\cdots&\mu_{rr} \end{bmatrix}$ 。

將條件方程組寫成矩陣形式：

$SW=WM$ 。

因為 $W$ 有單範正交行向量， $W^TW=I_r$ 。上式左乘 $W^T$ ，即有 $W^TSW=M$ 。下面證明滿足條件式的任何解 $W$ 和 $M$ 都有相同的目標函數值：

$\displaystyle \sum_{j=1}^r\mathbf{w}_j^TS\mathbf{w}_j=\mathrm{trace}(W^TSW)$ 。

將實對稱矩陣 $M$ 正交對角化為 $D=Q^TMQ$ ，其中 $Q$ 是一正交矩陣， $Q^TQ=QQ^T=I_r$ ， $D=\mathrm{diag}(\lambda_1,\ldots,\lambda_r)$ 。條件式 $W^TSW=M$ 左乘 $Q^T$ ，右乘 $Q$ ，可得

$Q^TW^TSWQ=(WQ)^TS(WQ)=\tilde{W}^TS\tilde{W}=D$ ，

上面我們令 $\tilde{W}=WQ$ 。因為 $\tilde{W}^T\tilde{W}=(WQ)^T(WQ)=Q^TW^TWQ=Q^TQ=I_r$ ， $\tilde{W}$ 也有單範正交行向量，故知 $\tilde{W}$ 和 $D$ 也是一組解。利用跡數循環不變性，

$\mathrm{trace}(W^TSW)=\mathrm{trace}(Q\tilde{W}^TS\tilde{W}Q^T)=\mathrm{trace}(\tilde{W}^TS\tilde{W}Q^TQ)=\mathrm{trace}(\tilde{W}^TS\tilde{W})$ ，

證明 $W$ 和 $\tilde{W}$ 所含的兩組單範正交集有相同的目標函數值 $\mathrm{trace}D=\sum_{j=1}^r\lambda_j$ 。

[5] 假設我們已經求得最大化 $\mathbf{w}_1^TS\mathbf{w}_1$ 的單位向量 $\mathbf{w}_1$ ，滿足 $S\mathbf{w}_1=\lambda_1\mathbf{w}_1$ 。下一步要找出單位向量 $\mathbf{w}_2$ 使最大化 $\mathbf{w}_2^TS\mathbf{w}_2$ ，並滿足 $\mathbf{w}_2^T\mathbf{w}_1=0$ 。寫出

$\displaystyle L(\mathbf{w}_2,\alpha,\beta)=\mathbf{w}^T_2S\mathbf{w}_2-\alpha(\mathbf{w}_2^T\mathbf{w}_2-1)-\beta\mathbf{w}_2^T\mathbf{w}_1$ ，

其中 $\alpha,\beta$ 是 Lagrange 乘數。計算偏導數，設結果為零：

$\displaystyle \frac{\partial L}{\partial \mathbf{w}_2}=2S\mathbf{w}_2-2\alpha\mathbf{w}_2-\beta\mathbf{w}_1=\mathbf{0}$ 。

上式左乘 $\mathbf{w}_1^T$ ，可得 $\mathbf{w}_1^TS\mathbf{w}_2=\alpha\mathbf{w}_1^T\mathbf{w}_2+\frac{\beta}{2}\mathbf{w}_1^T\mathbf{w}_1=\frac{\beta}{2}$ 。另一方面， $\mathbf{w}_1^TS\mathbf{w}_2=(S\mathbf{w}_1)^T\mathbf{w}_2=(\lambda_1\mathbf{w}_1)^T\mathbf{w}_2=0$ 。所以， $\beta=0$ ，即得 $S\mathbf{w}_2=\alpha\mathbf{w}_2$ 。明顯地，我們選擇 $\alpha=\lambda_2$ ， $\mathbf{w}_2$ 是對應的單位特徵向量。同樣的推導步驟可以推廣至 $r=3,4,\ldots$ 的情況。

This entry was posted in 機器學習 and tagged 相關係數, Lagrange 乘數法, PCA, SVD, 共變異數矩陣, 半正定矩陣, 奇異值分解, 樣本共變異數, 樣本平均數, 主成分分析. Bookmark the permalink.

20 Responses to 主成分分析

小葉 says:

04/15/2013 at 9:34 pm

最近又去聽了微積分的課，剛好聽到了講述simpison數值積分方法。數值積分方法特別可以處理非平滑的函數數據。其中可達到一定精確度的simpison方法，用了拋物線來近似分散的數據，其手法就是揉合上述的正交投影以及黎曼和的概念。因為正交投影，所以可以控制誤差到一定範圍內，黎曼積分也就可以派上用場。只是我好奇的是，因為訊號處理常常用到傅立葉分析來搜尋出有用的訊號，那麼是否可以用傅立葉分析來做數值積分呢?

Reply
- ccjou says:
  
  04/16/2013 at 8:33 am
  
  Simpson 積分法則的原理是Newton-Cotes公式： $\int_a^bf(x)dx\approx\sum_{i=0}^nw_if(x_i)$ ， $x_i=hi+x_0$ ， $h=(b-a)/n$ 。通常我們使用一內插多項式來近似 $f(x)$ 以得到 $w_i$ 的計算公式。這個式子看起來很像是Hilbert空間的函數近似問題：最小化 $\Vert f-\sum_{i=1}^nw_i\phi_i\Vert^2$ ， $\{\phi_i\}_{i=1}^n$ 是我們喜愛的基底，最佳近似確實是 $f$ 在這個基底所擴張的子空間的正交投影。
  
  傅立葉變換是一種積分變換，它可以用來計算函數近似的三角級數(trigonometric series)，或證明一些定理，但我不知道是否亦可作數值積分用。
  
  Reply
- TUK says:
  
  08/10/2020 at 2:52 pm
  
  傅立葉積分轉換是把時域的資訊轉換成頻率域的資訊，你可以理解成他把時域中某個頻率的能量積分然後放到頻率域中去。所以沒錯，傅立葉轉換本身是有積分的效果在的，只是他把要拿來積分的東西分離成正交的獨立函數分別處理。
  
  https://zh.wikipedia.org/wiki/%E8%B0%B1%E5%AF%86%E5%BA%A6
  https://zh.wikipedia.org/wiki/%E5%B8%95%E5%A1%9E%E7%93%A6%E5%B0%94%E5%AE%9A%E7%90%86
  
  Reply
Hermione says:

08/27/2013 at 5:06 pm

請問主成分分析是否一定要進行標準化?

Reply
- ccjou says:
  
  08/27/2013 at 6:33 pm
  
  如果變數的尺規(scale)相似，譬如，全部都是EEG訊號，那麼使用共變異數矩陣即可。如果變數的尺規差異很大，譬如，年所得，年齡，體重，那麼應該使用相關係數矩陣(即標準化後的共變異數矩陣)。
  
  Reply
張盛東 says:

12/25/2013 at 10:33 am

周老師，我有兩個問題：
1）對E1求ck的偏導數是否第一個2前面的負號應該去掉？
2）對L({Wj},{u[i,j]})求偏導後，為什麼可以假設u[i,j]+u[j,i]=0從而將wj消去？

Reply
- ccjou says:
  
  12/25/2013 at 11:08 am
  
  謝謝指正。
  
  你可以試一下 $r=2$ 的例子，未知數有4個： $\mu_{11},\mu_{12},\mu_{21},\mu_{22}$ 。
  線性方程組為
  $\begin{bmatrix} 2w_1&w_2&w_2&0\\ 0&w_1&w_1&2w_2 \end{bmatrix}\begin{bmatrix} \mu_{11}\\ \mu_{12}\\ \mu_{21}\\ \mu_{22} \end{bmatrix}=\begin{bmatrix} 2Sw_1\\ 2Sw_2 \end{bmatrix}$
  係數矩陣A只有三個線性獨立的行向量(因為 $w_1$ 和 $w_2$ 不為零向量)，故 $\text{rank} A=3$ ，也就是說有無限多組解，故可設 $\mu_{12}+\mu_{21}=0$ 。同理，對於 $r^2$ 個未知數 $\{\mu_{ij}\}$ ， $\text{rank}A=r+\binom{r}{2}=(r^2+r)/2$ ，即有 $\dim N(A)=r^2-(r^2+r)/2=(r^2-r)/2=\binom{r}{2}$ 。所以我們可以放心地設 $\mu_{ij}+\mu_{ji}=0$ ， $i\neq j$ 。
  
  突然想起在高中生時，我的數學老師常告誡我要將算式寫清楚。
  
  Reply
  - 張盛東 says:
    
    12/25/2013 at 3:32 pm
    
    老師，我還是有一些不明白。我自己推導一次後覺得u[i,j]+u[j,i]=0是必然發生的事，為什麼老師您要“設”呢？
    
    Reply
    - ccjou says:
      
      12/25/2013 at 3:37 pm
      
      啊，真的嗎？你有空時方便將過程貼上網或寄給我看看嗎？
      
      $S\mathbf{w}_j=\mu_{jj}\mathbf{w}_j+\frac{1}{2}\sum_{i\neq j}(\mu_{ij}+\mu_{ji})\mathbf{w}_i,~~j=1,\ldots,r$
      
      左乘 $\mathbf{w}_i^T$ ， $i\neq j$ ，可得
      
      $\mathbf{w}_i^TS\mathbf{w}_j=\frac{1}{2}(\mu_{ij}+\mu_{ji})$
      
      上式不能推斷 $\mu_{ij}+\mu_{ji}=0$ ，雖然 $\mathbf{w}_i^T\mathbf{w}_j=0$ ，但 $\mathbf{w}_i^TS\mathbf{w}_j$ 未必等於零。
      
      Reply
      - 張盛東 says:
        
        12/26/2013 at 11:01 pm
        
        我說“推導”有些誇張了，其實我就寫了幾個步驟然後直覺地覺得，u[i,j]和u[j,i]（i not equal to j）所對應的在A中的行向量必然一樣,所以所有形如
        c * [0 0 … 1 … -1 0 …]’ 必然在A的零空間中，其中c為某一非零實數，1位於u[i,j]的位置，-1位於u[j,i]的位置。所有這些互為正交向量構成了A的零空間的基底，因此所有齊次解滿足u[i,j]+u[j,i]=0 if i is not equal to j。要令方程組成立，[Sw[1],…Sw[r]]’必須在X=span{W[1],…,w[r]}中，也就是說X是S的不變子空間。後來發現在推導的時候想當然地認為X就是S的特徵空間(eigenspace)所以必然有Sw[k]=lamda*w[k]從而得到方程組的特解中也必須有u[i,j]+u[j,i]=0 if i is not equal to j這個結論。我曾嘗試過證明X是滿足所有條件的唯一不變子空間但失敗了。
        老師，不知道為什麼，我對這個假設沒有“安全感”，萬一該假設不成立那麼之後的證明就只有或然性而不是必然性了。
        
        Reply
        
        ccjou says:
        
        12/27/2013 at 8:57 am
        
        我明白了。沒有安全感是好事。朱熹說：「讀書，始讀未知有疑。其次則漸漸有疑，中則節節是疑。過了這一番后，疑漸漸解，以至融會貫通，都無所疑，方始是學。」
張盛東 says:

12/25/2013 at 12:07 pm

謝謝老師。

Reply
ccjou says:

12/27/2013 at 12:14 pm

我們可以從另一個角度切入：根據註解二，待解的方程式為 $SW=WM$ ， $W^TW=I$ 。因為 ~~W^T=W~~ $M^T=M$ ，故可正交對角化為 $M=QDQ^T$ ，其中 $Q^TQ=I$ ， $D=\text{diag}(d_1,\ldots,d_r)$ 。代入可得 $SW=WQDQ^T$ ，右乘 $Q$ ， $SWQ=WQD$ 。令 $V=WQ$ ， $V^TV=Q^TW^TWQ=I$ 。問題變成 $SV=VD$ ，也就是 $S\mathbf{v}_j=d_{j}\mathbf{v}_j$ ， $1\le j\le r$ 。

Reply
- 張盛東 says:
  
  12/27/2013 at 2:11 pm
  
  經過老師一點明果然豁然開朗。另外，老師的意思是因為M^T=M所以M可正交對角話吧。
  
  非常感謝老師。
  
  Reply
張盛東 says:

12/31/2013 at 2:03 pm

周老師，請問一下主成份分析和主成份回歸的關係為何？主成份回歸是否只是主成份分析的一個應用而已？

Reply
- ccjou says:
  
  12/31/2013 at 3:07 pm
  
  PCR (principal components regression) 就是 PCA+multiple regression。俗話說：蟑螂怕拖鞋，烏龜怕鐵鎚。在資料分析，線性模型怕共線性(collinearity)，非線性模型怕過度擬合(over-fitting)。如果數據輸入矩陣X的維數很大(變數很多)且變數之間相關，則multiple regression的估計係數的變異很大(表示模型不可信賴)。為了解決這個問題，先對X進行主成分分析，取得低維數的主成分係數矩陣Z(各變數無關)，透過適當的變數選擇方法(譬如計算輸出y和z_i的相關係數)挑選出一部分的z變數當作新的數據輸入，然後建模，再將估計出的z_i係數轉換回變數x_j的係數，此法稱為PCR estimator。詳見
  
  http://en.wikipedia.org/wiki/Principal_component_regression
  
  Reply
冠亨 says:

06/09/2017 at 9:38 am

周老師您好：
https://stats.stackexchange.com/questions/66926/what-are-the-four-axes-on-pca-biplot
目前資料分析使用主成分分析方法通常會畫一張叫做 biplot的圖來看看各筆資料在第一集第二主成分所組成的平面上到底比較受到哪一個原始參數的影響，在上述網頁中有提到這張圖到底是怎麼畫的，但是為什麼在biplot中的紅色箭頭向量並非只是主成分中各個變數佔的比例，而是還得乘上(主成分的解釋變異*資料筆數^(1/2))？

Reply
- ccjou says:
  
  06/09/2017 at 10:20 am
  
  請參考
  
  主成分分析與奇異值分解
  
  裡面的factor loading matrix
  
  Reply
duby says:

07/25/2017 at 9:14 am

周老師，注釋3最後一個式子的第一項應該有個n-1的系數

Reply
- ccjou says:
  
  07/31/2017 at 3:52 pm
  
  非常謝謝，已訂正。
  
  Reply

	xmj on 內積的定義
	Ning ChingSan on 線性代數的第一堂課──矩陣乘法的定義
	momo on 兩岸線性代數用詞參照
	訪客 on 克拉瑪公式的簡易幾何證明
	悟 on 條件機率與貝氏定理
	jeremy on 內積與外積是怎麼來的？