主成分分析 | Principal Components Analysis | PCA
理論
僅僅使用基本的線性代數知識，就可以推匯出一種簡單的機器學習演算法，主成分分析（Principal Components Analysis, PCA）。
假設有 $m$ 個點的集合：$\left\{\boldsymbol{x}^{(1)}, \ldots, \boldsymbol{x}^{(m)}\right\}$ in $\mathbb{R}^{n}$，我們希望對這些點進行有失真壓縮（lossy compression）。有失真壓縮是指，失去一些精度作為代價，用更少的儲存空間來儲存這些點。我們當然希望損失的精度儘可能少。
一種方法是，我們對這些點進行編碼：以低維的方式，來表示這些點。
那就是說，對於任意點：$\boldsymbol{x}^{(i)} \in \mathbb{R}^{n}$，我們會計算相應的碼向量（code vector）：$\boldsymbol{c}^{(i)} \in \mathbb{R}^{l}$。
如果 $l$ 小於 $m$，那麼我們就是用一個低維向量來表示一個高維向量。需要的儲存空間變小，當然同時可能會有精度損失。
我們希望知道某個編碼函式（encoding function），根據輸入，計算編碼（code），即滿足：
\begin{equation}
f(\boldsymbol{x})=\boldsymbol{c}
\end{equation}
並且還希望知道某個解碼函式（decoding function），根據編碼，計算重構的輸入（reconstructed input），即滿足：
\begin{equation}
\boldsymbol{x} \approx g(f(\boldsymbol{x}))
\end{equation}
PCA 就是選擇特定的解碼函式。具體說來，我們選擇矩陣乘法，作為從編碼空間到 $\mathbb{R}^{n}$ 的對映。令 $g(\boldsymbol{c})=\boldsymbol{D} \boldsymbol{c}$，其中 $\boldsymbol{D} \in \mathbb{R}^{n \times l}$ 就是解碼函式的矩陣定義。
PCA 還限制了 $\boldsymbol{D}$ 的每一列互為正交（注意：除非 $l=n$，否則 $\boldsymbol{D}$ 不是嚴格意義上的正交矩陣）。
此外，為了不受尺度的影響，我們還限定了 $\boldsymbol{D}$ 的每一列都是單位向量。
為了把基本思想轉換為可實現的演算法，首先我們需要知道：如何生成每個輸入點 $\boldsymbol{x}$ 的最優編碼點 $C^{*}$。其中一個方法就是，最小化輸入點 $\boldsymbol{x}$ 和重構點 $g\left(\boldsymbol{c}^{*}\right)$ 的距離。在主成分演算法中，我們使用 $L^{2}$ 範數：
\begin{equation}
\boldsymbol{c}^{*}=\underset{\boldsymbol{c}}{\arg \min }\|\boldsymbol{x}-g(\boldsymbol{c})\|_{2}
\end{equation}
我們可以切換為計算平方的 $L^{2}$ 範數，它們優化的 $\boldsymbol{c}$ 值相同：
\begin{equation}
\boldsymbol{c}^{*}=\underset{\boldsymbol{c}}{\arg \min }\|\boldsymbol{x}-g(\boldsymbol{c})\|_{2}^{2}
\end{equation}
該函式可以化簡為：
\begin{equation}
\begin{array}{l}
(\boldsymbol{x}-g(\boldsymbol{c}))^{\top}(\boldsymbol{x}-g(\boldsymbol{c})) \\
=\boldsymbol{x}^{\top} \boldsymbol{x}-\boldsymbol{x}^{\top} g(\boldsymbol{c})-g(\boldsymbol{c})^{\top} \boldsymbol{x}+g(\boldsymbol{c})^{\top} g(\boldsymbol{c}) \\
=\boldsymbol{x}^{\top} \boldsymbol{x}-2 \boldsymbol{x}^{\top} g(\boldsymbol{c})+g(\boldsymbol{c})^{\top} g(\boldsymbol{c})
\end{array}
\end{equation}
我們可以刪除第一項，因為它和 $\boldsymbol{x}$ 無關，於是方程最終可以化簡為：
\begin{equation}
\boldsymbol{c}^{*}=\underset{\boldsymbol{c}}{\arg \min }-2 \boldsymbol{x}^{\top} g(\boldsymbol{c})+g(\boldsymbol{c})^{\top} g(\boldsymbol{c})
\end{equation}
根據 $g(\boldsymbol{c})$ 的定義，代入可得：
\begin{equation}
\begin{aligned}
\boldsymbol{c}^{*} &=\underset{\boldsymbol{c}}{\arg \min }-2 \boldsymbol{x}^{\top} \boldsymbol{D} \boldsymbol{c}+\boldsymbol{c}^{\top} \boldsymbol{D}^{\top} \boldsymbol{D} \boldsymbol{c} \\
&=\arg \min -2 \boldsymbol{x}^{\top} \boldsymbol{D} \boldsymbol{c}+\boldsymbol{c}^{\top} \boldsymbol{I}_{l} \boldsymbol{c} \\
&=\arg \min -2 \boldsymbol{x}^{\top} \boldsymbol{D} \boldsymbol{c}+\boldsymbol{c}^{\top} \boldsymbol{c}
\end{aligned}
\end{equation}
使用向量積分，我們可以求解這個優化問題：
\begin{equation}
\begin{array}{c}
\nabla_{\boldsymbol{c}}\left(-2 \boldsymbol{x}^{\top} \boldsymbol{D} \boldsymbol{c}+\boldsymbol{c}^{\top} \boldsymbol{c}\right)=\mathbf{0} \\
-2 \boldsymbol{D}^{\top} \boldsymbol{x}+2 \boldsymbol{c}=\mathbf{0} \\
\boldsymbol{c}=\boldsymbol{D}^{\top} \boldsymbol{x}
\end{array}
\end{equation}
這使得該演算法很高效，我們可以用矩陣-向量運算，來最優地對 $\boldsymbol{x}$ 進行編碼。要對一個向量進行編碼，我們使用編碼函式：
\begin{equation}
f(\boldsymbol{x})=\boldsymbol{D}^{\top} \boldsymbol{x}
\end{equation}
我們同樣可以定義 PCA 的重構運算：
\begin{equation}
r(\boldsymbol{x})=g(f(\boldsymbol{x}))=\boldsymbol{D} \boldsymbol{D}^{\top} \boldsymbol{x}
\end{equation}
下一步，我們需要選擇解碼矩陣 $\boldsymbol{D}$，與前面最小化一個點上的重構誤差不同，由於我們需要對所有點都使用相同的矩陣 $\boldsymbol{D}$，因此我們必須針對所有點，最小化誤差矩陣的 Frobenius 範數（又稱 F-範數）：
\begin{equation}
\boldsymbol{D}^{*}=\underset{\boldsymbol{D}}{\arg \min } \sqrt{\sum_{i, j}\left(x_{j}^{(i)}-r\left(\boldsymbol{x}^{(i)}\right)_{j}\right)^{2}} \text { subject to } \boldsymbol{D}^{\top} \boldsymbol{D}=\boldsymbol{I}_{l}
\end{equation}
為了推導計算 $\boldsymbol{D}^{*}$ 的演算法，我們首先考慮 $l=1$ 的情況。此時，$\boldsymbol{D}^{*}$ 就是一個向量：$d$。將式（10）代入到式（11），可得：
\begin{equation}
\boldsymbol{d}^{*}=\underset{\boldsymbol{d}}{\arg \min } \sum_{i}\left\|\boldsymbol{x}^{(i)}-\boldsymbol{d} \boldsymbol{d}^{\top} \boldsymbol{x}^{(i)}\right\|_{2}^{2} \text { subject to }\|\boldsymbol{d}\|_{2}=1
\end{equation}
$\boldsymbol{d}^{\top} \boldsymbol{x}^{(i)}$ 是一個標量，移動到式子左邊，可得：
\begin{equation}
\boldsymbol{d}^{*}=\underset{\boldsymbol{d}}{\arg \min } \sum_{i}\left\|\boldsymbol{x}^{(i)}-\boldsymbol{d}^{\top} \boldsymbol{x}^{(i)} \boldsymbol{d}\right\|_{2}^{2} \text { subject to }\|\boldsymbol{d}\|_{2}=1
\end{equation}
標量的轉置等於標量本身，可得：
\begin{equation}
\boldsymbol{d}^{*}=\underset{\boldsymbol{d}}{\arg \min } \sum_{i}\left\|\boldsymbol{x}^{(i)}-\boldsymbol{x}^{(i) \top} \boldsymbol{d} \boldsymbol{d}\right\|_{2}^{2} \text { subject to }\|\boldsymbol{d}\|_{2}=1
\end{equation}
令 $\boldsymbol{X} \in \mathbb{R}^{m \times n}$，其中 $\boldsymbol{X}_{i,:}=\boldsymbol{x}^{(i)^{\top}}$，方程整理為：
\begin{equation}
\boldsymbol{d}^{*}=\underset{\boldsymbol{d}}{\arg \min }\left\|\boldsymbol{X}-\boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right\|_{F}^{2} \text { subject to } \boldsymbol{d}^{\top} \boldsymbol{d}=1
\end{equation}
先不考慮約束條件，我們可以簡化 Frobenius 範數：
\begin{equation}
\begin{array}{l}
\underset{\boldsymbol{d}}{\arg \min }\left\|\boldsymbol{X}-\boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right\|_{F}^{2} \\
=\underset{\boldsymbol{d}}{\arg \min } \operatorname{Tr}\left(\left(\boldsymbol{X}-\boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right)^{\top}\left(\boldsymbol{X}-\boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right)\right) \\
=\underset{\boldsymbol{d}}{\arg \min } \operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X}-\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}-\boldsymbol{d} \boldsymbol{d}^{\top} \boldsymbol{X}^{\top} \boldsymbol{X}+\boldsymbol{d} \boldsymbol{d}^{\top} \boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right) \\
=\underset{\boldsymbol{d}}{\arg \min } \operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X}\right)-\operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right)-\operatorname{Tr}\left(\boldsymbol{d} \boldsymbol{d}^{\top} \boldsymbol{X}^{\top} \boldsymbol{X}\right)+\operatorname{Tr}\left(\boldsymbol{d} \boldsymbol{d}^{\top} \boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right) \\
=\underset{\boldsymbol{d}}{\arg \min }-\operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right)-\operatorname{Tr}\left(\boldsymbol{d} \boldsymbol{d}^{\top} \boldsymbol{X}^{\top} \boldsymbol{X}\right)+\operatorname{Tr}\left(\boldsymbol{d} \boldsymbol{d}^{\top} \boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right) \\
=\underset{\boldsymbol{d}}{\arg \min }-2 \operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right)+\operatorname{Tr}\left(\boldsymbol{d} \boldsymbol{d}^{\top} \boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right) \\
=\underset{\boldsymbol{d}}{\arg \min }-2 \operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right)+\operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top} \boldsymbol{d} \boldsymbol{d}^{\top}\right)
\end{array}
\end{equation}
現在，我們重新引入約束：
\begin{equation}
\begin{array}{l}
\underset{\boldsymbol{d}}{\arg \min }-2 \operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right)+\operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top} \boldsymbol{d} \boldsymbol{d}^{\top}\right) \text { subject to } \boldsymbol{d}^{\top} \boldsymbol{d}=1 \\
=\underset{\boldsymbol{d}}{\arg \min }-2 \operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right)+\operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right) \text { subject to } \boldsymbol{d}^{\top} \boldsymbol{d}=1 \\
=\underset{\boldsymbol{d}}{\arg \min }-\operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right) \text { subject to } \boldsymbol{d}^{\top} \boldsymbol{d}=1 \\
=\underset{\boldsymbol{d}}{\arg \max } \operatorname{Tr}\left(\boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d} \boldsymbol{d}^{\top}\right) \text { subject to } \boldsymbol{d}^{\top} \boldsymbol{d}=1 \\
=\underset{\boldsymbol{d}}{\arg \max } \operatorname{Tr}\left(\boldsymbol{d}^{\top} \boldsymbol{X}^{\top} \boldsymbol{X} \boldsymbol{d}\right) \text { subject to } \boldsymbol{d}^{\top} \boldsymbol{d}=1
\end{array}
\end{equation}
我們可以使用特徵值分解，來求解這個優化問題。具體地說，通過計算 $\boldsymbol{X}^{\top} \boldsymbol{X}$ 最大特徵值的特徵向量，可以求出最優值 $d$ 。
這個推導式針對 $l=1$ 的情況而言的，它只能夠求出第一個主成分。在更一般的情況，如果我們想要求出主成分的基，矩陣 $\boldsymbol{D}$ 由最大特徵值的 $l$ 個特徵向量確定。這個可以通過歸納法證明。
實踐
PCA 通過資料降維，常見的應用有加速機器學習演算法和資料視覺化。
資料視覺化
載入 Iris 資料集
1 import pandas as pd 2 3 url = "https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data" 4 # load dataset into Pandas DataFrame 5 df = pd.read_csv(url, names=['sepal length', 'sepal width', 'petal length', 'petal width', 'target']) 6 print(df.head())
輸出：
sepal length sepal width petal length petal width target 0 5.1 3.5 1.4 0.2 Iris-setosa 1 4.9 3.0 1.4 0.2 Iris-setosa 2 4.7 3.2 1.3 0.2 Iris-setosa 3 4.6 3.1 1.5 0.2 Iris-setosa 4 5.0 3.6 1.4 0.2 Iris-setosa
資料標準化
PCA 會受到資料尺度的影響，因此在使用 PCA 之前最好進行標準化。使用 StandardScaler 可以標準化資料集特徵到一個單位尺度（均值為 0，方差為 1）：這是很多機器學習演算法效能優化的條件。
1 from sklearn.preprocessing import StandardScaler 2 3 features = ['sepal length', 'sepal width', 'petal length', 'petal width'] 4 # Separating out the features 5 x = df.loc[:, features].values 6 # Separating out the target 7 y = df.loc[:, ['target']].values 8 # Standardizing the features 9 x = StandardScaler().fit_transform(x)
使用 PCA 投影到 2D
下面程式碼將原始資料從 4 維投影到 2 維。降維後的每個主成分通常不會有特定的語義。
1 from sklearn.decomposition import PCA 2 3 pca = PCA(n_components=2) 4 principalComponents = pca.fit_transform(x) 5 principalDf = pd.DataFrame(data=principalComponents, 6 columns=['principal component 1', 'principal component 2'])
最後把特徵和標籤放在一起：
finalDf = pd.concat([principalDf, df[['target']]], axis=1)
2D 投影的視覺化
下面程式碼繪製二維資料：
1 import matplotlib.pyplot as plt
2
3 fig = plt.figure(figsize=(8, 8))
4 ax = fig.add_subplot(1, 1, 1)
5 ax.set_xlabel('Principal Component 1', fontsize=15)
6 ax.set_ylabel('Principal Component 2', fontsize=15)
7 ax.set_title('2 component PCA', fontsize=20)
8 targets = ['Iris-setosa', 'Iris-versicolor', 'Iris-virginica']
9 colors = ['r', 'g', 'b']
10 for target, color in zip(targets, colors):
11 indicesToKeep = finalDf['target'] == target
12 ax.scatter(finalDf.loc[indicesToKeep, 'principal component 1']
13 , finalDf.loc[indicesToKeep, 'principal component 2']
14 , c=color
15 , s=50)
16 ax.legend(targets)
17 ax.grid()
18 plt.show()
輸出：
可解釋方差
可解釋方差（explained variance）可以體現出：有多少資訊（方差）來自每個主成分。當你把 4 維降為 2 維時，你損失了一些資訊。通過 explained_variance_ratio_，你可以看到主成分包含了多少資訊。
print(pca.explained_variance_ratio_)
輸出：
[0.72770452 0.23030523]
加速機器學習演算法
如果你的學習演算法，由於輸入維度太高的原因，變得很慢，那麼使用 PCA 來加速可能是一個不錯的選擇。
這裡我們使用 MNIST 資料集。
下載並載入資料
1 from sklearn.datasets import fetch_openml
2
3 mnist = fetch_openml('mnist_784')
這裡有 70000 張影象，每張影象維度是 784。
將資料劃分為訓練集和測試集
1 from sklearn.model_selection import train_test_split 2 3 # test_size: what proportion of original data is used for test set 4 train_img, test_img, train_lbl, test_lbl = train_test_split(mnist.data, mnist.target, test_size=1 / 7.0, random_state=0)
資料標準化
1 from sklearn.preprocessing import StandardScaler 2 3 scaler = StandardScaler() 4 # Fit on training set only. 5 scaler.fit(train_img) 6 # Apply transform to both the training set and the test set. 7 train_img = scaler.transform(train_img) 8 test_img = scaler.transform(test_img)
匯入並使用 PCA
與硬編碼選擇多少主成分不同，下面程式碼選擇剛好可以保留 95% 以上資訊的主成分個數，因此主成分個數是可變的。
1 from sklearn.decomposition import PCA 2 3 # Make an instance of the Model 4 pca = PCA(0.95)
再訓練集進行 PCA 擬合：
pca.fit(train_img)
檢視有多少個主成分：
print(pca.n_components_)
輸出：
327
對訓練集和測試集使用 transform
train_img = pca.transform(train_img) test_img = pca.transform(test_img)
對變換後的資料使用邏輯迴歸
1 import time
2 from sklearn.linear_model import LogisticRegression
3
4 # lbfgs 很快但是可能警報不能收斂
5 logisticRegr = LogisticRegression(solver='lbfgs')
6 t0 = time.time()
7 logisticRegr.fit(train_img, train_lbl)
8 print(f'time elapsed: {time.time() - t0}')
9
10 print(logisticRegr.score(test_img, test_lbl))
注：如果使用 lbfgs 報錯：
AttributeError: 'str' object has no attribute 'decode'
scikit-learn 需要升級到 0.24 以上。
使用 PCA 後的擬合時間
以下是某一組測試效能圖表，可以作為參考：
壓縮表徵後的影象重構
PCA 當然也可以對影象進行壓縮：
如果對這段程式碼感興趣，可以參考：github。
參考
- Goodfellow I, Bengio Y, Courville A, et al. Deep learning[M]. Cambridge: MIT press, 2016.
- PCA using Python (scikit-learn)