作者：韓信子@ShowMeAI
教程地址：http://www.showmeai.tech/tutorials/34
本文地址：http://www.showmeai.tech/article-detail/198
宣告：版權所有，轉載請聯絡平臺與作者並註明出處
引言
在網際網路大資料場景下，我們經常需要面對高維資料，在對這些資料做分析和視覺化的時候，我們通常會面對「高維」這個障礙。在資料探勘和建模的過程中，高維資料也同樣帶來大的計算量，佔據更多的資源，而且許多變數之間可能存在相關性，從而增加了分析與建模的複雜性。
我們希望找到一種方法，在對資料完成降維「壓縮」的同時，儘量減少資訊損失。由於各變數之間存在一定的相關關係，因此可以考慮將關係緊密的變數變成儘可能少的新變數，使這些新變數是兩兩不相關的，那麼就可以用較少的綜合指標分別代表存在於各個變數中的各類資訊。機器學習中的降維演算法就是這樣的一類演算法。
主成分分析（Principal Components Analysis，簡稱PCA）是最重要的資料降維方法之一。在資料壓縮消除冗餘和資料噪音消除等領域都有廣泛的應用。本篇我們來展開講解一下這個演算法。
（本篇降維演算法部分內容涉及到機器學習基礎知識，沒有先序知識儲備的寶寶可以檢視ShowMeAI的文章 圖解機器學習 | 機器學習基礎知識。
1.PCA與最大可分性
對於 \(X = \begin {bmatrix} x_1 \\ x_2 \\ ... \\ x_n \end{bmatrix}\) ， \(X \in R^n\) 。我們希望 \(X\) 從 \(n\) 維降到 \(n^{'}\) 維，同時希望資訊損失最少。比如，從 \(n = 2\) 維降到 \(n^{'} = 1\) 。
左圖為一個典型的例子，假如我們要對一系列人的樣本進行資料降維（每個樣本包含「身高」「體重」兩個維度）。右圖我們既可以降維到第一主成分軸，也可以降維到第二主成分軸。
哪個主成分軸更優呢？從直觀感覺上，我們會認為「第一主成分軸」優於「第二主成分軸」，因為它比較大程度保留了資料之間的區分性（保留大部分資訊）。
對PCA演算法而言，我們希望找到小於原資料維度的若干個投影座標方向，把資料投影在這些方向，獲得壓縮的資訊表示。下面我們就一步一步來推導一下PCA演算法原理。
2.基變換
先來複習一點點數學知識。我們知道要獲得原始資料 \(X\) 新的表示空間 \(Y\) ，最簡單的方法是對原始資料進行線性變換（也叫做基變換） \(Y = PX\) 。其中， \(X\) 是原始樣本， \(P\) 是基向量， \(Y\) 是新表達。
數學表達為：
-
其中 \(p_i\) 是行向量，表示第 \(i\) 個基；
-
\(x_j\) 是一個列向量，表示第 \(j\) 個原始資料記錄。
當 \(r < n\) 時，即「基的維度<資料維度」時，可達到降維的目的，即 \(X \in R^{n \times m} \rightarrow Y \in R^{r \times m }\) 。
以直角座標系下的點 \((3,2)\) 為例，要把點 \((3,2)\) 變換為新基上的座標，就是用 \((3,2)\) 與第一個基做內積運算，作為第一個新的座標分量，然後用 \((3,2)\) 與第二個基做內積運算，作為第二個新座標的分量。
上述變化，線上性代數裡，我們可以用矩陣相乘的形式簡潔的來表示：
再稍微推廣一下，假如我們有m個二維向量，只要將二維向量按列排成一個兩行m列矩陣，然後用「基矩陣」乘以這個矩陣，就得到了所有這些向量在新基下的值。例如(1,1)、(2,2)、(3,3)，想變換到剛才那組基上，可以如下這樣表示：
3.方差
在本文的開始部分，我們提到了，降維的目的是希望壓縮資料但資訊損失最少，也就是說，我們希望投影后的資料儘可能分散開。在數學上，這種分散程度我們用「方差」來表達，方差越大，資料越分散。
-
定義方差 \(Var\) ：對於單一隨機變數 \(a\) ， \(Var(a) = \frac{1}{m} \sum_{i = 1}^m (a_i - \mu)^2\)
-
對資料做去中心化（方便後面操作）： \(Var(a) = \frac{1}{m} \sum_{i = 1}^m a_i ^2\)
\(Var(a)\) 表示 \(a\) 的取值與其數學期望之間的偏離程度。若 \(Var(a)\) 較小，意味著 \(a\) 的取值主要集中在期望 \(\mu\) 也就是 \(E(a)\) )的附近；反之，若 \(Var(a)\) 較大，意味著 \(a\) 的取值比較分散。
我們來看一個具體的例子。假設我們5個樣本資料，分別是 \(x_1 = \begin{bmatrix} 1 \\ 1 \end{bmatrix}\) 、 \(x_2 = \begin{bmatrix} 1 \\ 3\end{bmatrix}\) 、 \(x_3 = \begin{bmatrix} 2 \\ 3\end{bmatrix}\) 、 \(x_4 = \begin{bmatrix} 4 \\ 4\end{bmatrix}\) 、 \(x_5 = \begin{bmatrix} 2 \\ 4 \end{bmatrix}\) ，將它們表示成矩陣形式： \(X = \begin{bmatrix} 1 & 1 & 2 & 4 & 2 \\ 1 & 3 & 3 & 4 & 4 \end {bmatrix}\) 。
為了後續處理方便，我們首先將每個欄位內所有值都減去欄位均值，其結果是將每個欄位都變為均值為0。
我們看上面的資料，設第一個特徵為 \(a\) ，第二個特徵為 \(b\) ，則某個樣本可以寫作 \(x_i = \begin{bmatrix} a \\ b \end {bmatrix}\)
且特徵 \(a\) 的均值為2，特徵 \(b\) 的均值為3。所以，變換後
4.協方差
協方差（Covariance）在概率和統計學中用於衡量兩個變數的總體誤差。比如對於二維隨機變數 \(x_i = \begin{bmatrix} a \\ b \end{bmatrix}\) ，特徵 \(a、b\) 除了自身的數學期望和方差，還需要討論 \(a、b\) 之間互相關係的數學特徵。
定義協方差 \(Cov\) ：
當 \(Cov(a, b) = 0\) 時，變數 \(a、b\) 完全獨立，這也是我們希望達到的優化目標。方差是協方差的一種特殊情況，即當兩個變數是相同的情況 \(Cov(a, a) = Var(a)\) 。
5.協方差矩陣
對於二維隨機變數 \(x_i = \begin{bmatrix} a \\ b \end {bmatrix}\) ，定義協方差矩陣 \(C = \begin{bmatrix} Var(a) & Cov(a, b) \\ Cov(b, a) &Var(b)\end{bmatrix}\) 。
對於 \(n\) 維隨機變數
我們可以看到，協方差矩陣是 \(n\) 行 \(n\) 列的對稱矩陣，主對角線上是方差，而協對角線上是協方差。
我們再來用一個示例對應講解一下。還是同樣的5個樣本資料
-
\(x_1 = \begin{bmatrix} 1 \\ 1 \end{bmatrix}\)
-
\(x_2 = \begin{bmatrix} 1 \\ 3\end{bmatrix}\)
-
\(x_3 = \begin{bmatrix} 2 \\ 3\end{bmatrix}\)
-
\(x_4 = \begin{bmatrix} 4 \\ 4\end{bmatrix}\)
-
\(x_5 = \begin{bmatrix} 2 \\ 4 \end{bmatrix}\)
去中心化後表示成矩陣
那如果有 \(m\) 個樣本的話， \(X =\begin{bmatrix} a_1 & a_2 & \cdots &a_m \\ b_1 & b_2 & \cdots & b_m\end{bmatrix}\) 。對 \(X\) 做一些變換，用 \(X\) 乘以 \(X\) 的轉置，並乘上係數 \(1/m\) ：
這正是協方差矩陣！我們歸納得到：設我們有 \(m\) 個 \(n\) 維資料記錄，將其按列排成 \(n\) 乘 \(m\) 的矩陣 \(X\) ，設 \(C = \frac{1}{m}XX^T\) ，則 \(C\) 是一個對稱矩陣，其對角線分別個各個特徵的方差，而第 \(i\) 行 \(j\) 列和 \(j\) 行 \(i\) 列元素相同，表示 \(i\) 和 \(j\) 兩個特徵之間的協方差。
6.協方差矩陣對角化
再回到我們的場景和目標：
-
現在我們有 \(m\) 個樣本資料，每個樣本有 \(n\) 個特徵，那麼設這些原始資料為 \(X\) ， \(X\) 為 \(n\) 行 \(m\) 列的矩陣。
-
想要找到一個基 \(P\) ，使 \(Y_{r \times m} = P_{r \times n}X_{n \times m}\) ，其中 $r<n $，達到降維的目的。
設 \(X\) 的協方差矩陣為 \(C\) ， \(Y\) 的協方差矩陣為 \(D\) ，且 \(Y = PX\) 。
- 我們的目的變為：對原始資料 \(X\) 做PCA後，得到的 \(Y\) 的協方差矩陣 \(D\) 的各個方向方差最大，協方差為0。
那麼 \(C\) 與 \(D\) 是什麼關係呢？
我們發現，要找的 \(P\) 不是別的，而是能讓原始協方差矩陣對角化的 \(P\) 。
換句話說，優化目標變成了尋找一個矩陣 \(P\) ，滿足 \(PCP^T\) 是一個對角矩陣，並且對角元素按從大到小依次排列，那麼 $P的前 $K行就是要尋找的基，用 \(P\) 的前 \(K\) 行組成的矩陣乘以 \(X\) 就使得 \(X\) 從 \(N\) 維降到了 \(K\) 維並滿足上述優化條件。
最終我們聚焦在協方差矩陣對角化這個問題上。
由上文知道，協方差矩陣 \(C\) 是一個是對稱矩陣，線上性代數上，實對稱矩陣有一系列非常好的性質：
1）實對稱矩陣不同特徵值對應的特徵向量必然正交。
2）設特徵向量 \(\lambda\) 重數為 \(r\) ，則必然存在 \(r\) 個線性無關的特徵向量對應於 \(\lambda\) ，因此可以將這 \(r\) 個特徵向量單位正交化。
由上面兩條可知，一個 \(n\) 行 \(n\) 列的實對稱矩陣一定可以找到 \(n\) 個單位正交特徵向量，設這 \(n\) 個特徵向量為 \(e_1,e_2,⋯,e_n\) ，我們將其按列組成矩陣：
則對協方差矩陣 \(C\) 有如下結論：
其中 \(\Lambda\) 為對角矩陣，其對角元素為各特徵向量對應的特徵值（可能有重複）。
結合上面的公式：
其中， \(D\) 為對角矩陣，我們可以得到：
\(P\) 是協方差矩陣\(C\)的特徵向量單位化後按行排列出的矩陣，其中每一行都是 \(C\) 的一個特徵向量。如果設 \(P\) 按照 \(\Lambda\) 中特徵值的從大到小，將特徵向量從上到下排列，則用 \(P\) 的前 \(K\) \(K\)行組成的矩陣乘以原始資料矩陣 \(X\) ，就得到了我們需要的降維後的資料矩陣 \(Y\) 。
7.PCA演算法
總結一下PCA的演算法步驟：
設有 \(m\) 條 \(n\) 維資料。
1）將原始資料按列組成 \(n\) 行 \(m\) 列矩陣 \(X\)
2）將 \(X\) 的每一行（代表一個特徵）進行零均值化，即減去這一行的均值
3）求出協方差矩陣 \(C=\frac{1}{m}XX^T\)
4）求出協方差矩陣 \(C\) 的特徵值及對應的特徵向量
5）將特徵向量按對應特徵值大小從上到下按行排列成矩陣，取前 \(k\) 行組成矩陣 \(P\)
6） \(Y=PX\) 即為降維到 \(k\) 維後的資料
8.PCA程式碼實踐
我們這裡直接使用python機器學習工具庫scikit-learn來給大家演示PCA演算法應用（相關知識速查可以檢視ShowMeAI文章AI建模工具速查|Scikit-learn使用指南），sklearn工具庫中與PCA相關的類都在sklearn.decomposition包裡，最常用的PCA類就是sklearn.decomposition.PCA。
1）引數介紹
sklearn中的PCA類使用簡單，基本無需調參，一般只需要指定需要降維到的維度，或者降維後的主成分的方差和佔原始維度所有特徵方差和的比例閾值就可以了。
下面是sklearn.decomposition.PCA的主要引數介紹：
-
n_components：PCA降維後的特徵維度數目。
-
whiten：是否進行白化。所謂白化，就是對降維後的資料的每個特徵進行歸一化，讓方差都為1，預設值是False，即不進行白化。
-
svd_solver：奇異值分解SVD的方法，有4個可以選擇的值：{‘auto’,‘full’,‘arpack’,‘randomized’}。
除上述輸入引數，還有兩個PCA類的成員屬性也很重要：
-
① explained_variance_，它代表降維後的各主成分的方差值。
-
② explained_variance_ratio_，它代表降維後的各主成分的方差值佔總方差值的比例。
2）程式碼例項
# 構建資料樣本並視覺化
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
from sklearn.datasets import make_blobs
# X為樣本特徵，Y為樣本簇類別， 共1000個樣本，每個樣本3個特徵，共4個簇
X, y = make_blobs(n_samples=10000, n_features=3, centers=[[3,3, 3], [0,0,0], [1,1,1], [2,2,2]], cluster_std=[0.2, 0.1, 0.2, 0.2],
random_state =9)
fig = plt.figure()
ax = Axes3D(fig, rect=[0, 0, 1, 1], elev=30, azim=20)
plt.scatter(X[:, 0], X[:, 1], X[:, 2],marker='x')
先不降維，只對資料進行投影，看看投影后的三個維度的方差分佈，程式碼如下：
from sklearn.decomposition import PCA
pca = PCA(n_components=3)
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
輸出如下：
[0.98318212 0.00850037 0.00831751]
[3.78521638 0.03272613 0.03202212]
可以看出投影后三個特徵維度的方差比例大約為98.3%：0.8%：0.8%。投影后第一個特徵佔了絕大多數的主成分比例。現在我們來進行降維，從三維降到2維，程式碼如下：
pca = PCA(n_components=2)
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
輸出如下：
[0.98318212 0.00850037]
[3.78521638 0.03272613]
這個結果其實可以預料，因為上面三個投影后的特徵維度的方差分別為：[ 3.78521638 0.03272613]，投影到二維後選擇的肯定是前兩個特徵，而拋棄第三個特徵。為了有個直觀的認識，我們看看此時轉化後的資料分佈，程式碼如下：
X_new = pca.transform(X)
plt.scatter(X_new[:, 0], X_new[:, 1],marker='x')
plt.show()
從上圖可以看出，降維後的資料依然清楚可見之前三維圖中的4個簇。現在我們不直接指定降維的維度，而指定降維後的主成分方差和比例，來試驗一下。
pca = PCA(n_components=0.9)
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
print(pca.n_components_)
我們指定了主成分至少佔90%，輸出如下：
[0.98318212]
[3.78521638]
1
可見只有第一個投影特徵被保留。這也很好理解，我們的第一個主成分佔投影特徵的方差比例高達98%。只選擇這一個特徵維度便可以滿足90%的閾值。我們現在選擇閾值99%看看，程式碼如下：
pca = PCA(n_components=0.99)
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
print(pca.n_components_)
此時的輸出如下：
[0.98318212 0.00850037]
[3.78521638 0.03272613]
2
這個結果也很好理解，因為我們第一個主成分佔了98.3%的方差比例，第二個主成分佔了0.8%的方差比例，兩者一起可以滿足我們的閾值。最後我們看看讓MLE演算法自己選擇降維維度的效果，程式碼如下：
pca = PCA(n_components= 'mle',svd_solver='full')
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
print(pca.n_components_)
輸出結果如下：
[0.98318212]
[3.78521638]
1
可見由於我們的資料的第一個投影特徵的方差佔比高達98.3%，MLE演算法只保留了我們的第一個特徵。
更多無監督學習的演算法模型總結可以檢視ShowMeAI的文章 AI知識技能速查 | 機器學習-無監督學習。
參考連結
ShowMeAI相關文章推薦
- 1.機器學習基礎知識
- 2.模型評估方法與準則
- 3.KNN演算法及其應用
- 4.邏輯迴歸演算法詳解
- 5.樸素貝葉斯演算法詳解
- 6.決策樹模型詳解
- 7.隨機森林分類模型詳解
- 8.迴歸樹模型詳解
- 9.GBDT模型詳解
- 10.XGBoost模型最全解析
- 11.LightGBM模型詳解
- 12.支援向量機模型詳解
- 13.聚類演算法詳解
- 14.PCA降維演算法詳解
ShowMeAI系列教程推薦
- 圖解Python程式設計：從入門到精通系列教程
- 圖解資料分析：從入門到精通系列教程
- 圖解AI數學基礎：從入門到精通系列教程
- 圖解大資料技術：從入門到精通系列教程
- 圖解機器學習演算法：從入門到精通系列教程
