主成分分析PCA（Principal Component Analysis）在sklearn中的應用及部分源碼分析

最近太忙，又有一段時間沒寫東西了。

pca是機器學習中一個重要的降維技術，是特征提取的代表。關於pca的實現原理，在此不做過多贅述，相關參考書和各大神牛的博客都已經有各種各樣的詳細介紹。 如需學習相關數學理論，請移駕。T_T

簡單說一下pca的實現，首先對於一個矩陣X，我們計算X·X^T，顯然這個一個半正定矩陣，可以做特征值分解，然后取出k個最大的特征值及其對應的特征向量就可以表達整個原矩陣。若X·X^T=p^-1Λp，因為p是單位矩陣，所以p^-1=p^T，即X·X^T=p^-1·Λ^1/2·(p^-1·Λ^1/2)^T，也就是降維的X后用來p^-1·Λ^1/2表示。

其實從SVD的角度來理解也是一樣的，若X=UΣV^T，則X·X^T=UΣ²U^T，同樣我們用來UΣ來表示原X。

當我看sklearn的文檔時，文檔並沒有具體解釋它的方法得到的結果在數學上的表示什么，鑽研了半天，看了源碼后才知道。

sklearn的方法是通過SVD來實現的。這里着重介紹sklearn的pca類中的一個屬性(components_)和兩個方法(fit，transform)。

首先，給定一個矩陣，設置參數后，通過調用fit方法得到降維模型，也就是一個基矩陣。我們看一下fit中的部分關鍵代碼。

 

...
self.mean_ = np.mean(X, axis=0)
X -= self.mean_
U, S, V = linalg.svd(X, full_matrices=False)
U, V = svd_flip(U, V)
components_ = V
...
self.components_ = components_[:n_components]
...

 

首先做的工作就是對數據進行按列中心化，然后做svd分解，然后把V的前k個向量保存為模型，模型的關鍵內容就是components_。（看似幾百行代碼，T_T）

接下來看看transfrom的部分關鍵源碼。

 

...
X = np.dot(X, self.components_.T)
return X

 

我只貼這兩句，T_T，因為真的只有這兩句特別重要。也就是說我們要把一個新矩陣用到訓練好的pca模型中時，其實只是做了一次矩陣乘法而已。

 

怎么來理解作者的做法呢？

其實就是我們上面的提到的，用svd方式來實現pca時，我們實際上用UΣ來表示降維后的數據。綜合svd公式，可以看成XV=UΣ，也就是把原矩陣與V做乘法，實際上這里理解為投影，把X投影到V的單位正交基所表示的子空間中去得到X的低維表示。對於一個新的矩陣Y，同樣用YV來表示其在V子空間降維后的結果，這也就是為什么transform方法為什么最關鍵的步驟只有一步乘法了。回過頭看，fit訓練模型就是要得到V，然后在transform降維時只需要一步乘法就可以了。

我們用下面的代碼做個小實驗。

 

import numpy as np
from sklearn.decomposition import PCA
from sklearn import preprocessing
from sklearn.utils.extmath import svd_flip
 
svd = np.linalg.svd
 

A = np.random.randint(0,5, (5,3))

print('A=')
print(A)

pca = PCA(n_components=2, svd_solver='full')

# print(A)
model = pca.fit(A)
print('model.components_=')
print(model.components_)
X = model.transform(A)
print('pca(A)=')
print(X)

A = A - np.mean(A, axis=0)
u,s,v = svd(A, full_matrices=False)

print('V=')
print(v[:2])
print('pva_by_svd=')
print(np.dot(A, v[:2].T))

 

運行得到的結果如下

A=
[[3 1 4]
 [4 1 1]
 [1 2 0]
 [1 4 2]
 [1 3 2]]
model.components_=
[[ 0.70734192 -0.61231721  0.35317848]
 [ 0.19273774 -0.31363722 -0.92977624]]
pca(A)=
[[ 2.21911523 -1.47640531]
 [ 1.86692171  1.50566115]
 [-1.22059974  1.54358693]
 [-1.73887721 -0.94323999]
 [-1.12656    -0.62960277]]
V=
[[-0.70734192  0.61231721 -0.35317848]
 [ 0.19273774 -0.31363722 -0.92977624]]
pva_by_svd=
[[-2.21911523 -1.47640531]
 [-1.86692171  1.50566115]
 [ 1.22059974  1.54358693]
 [ 1.73887721 -0.94323999]
 [ 1.12656    -0.62960277]]

咦，結果跟之前的分析有點小小的差別，通過svd和sklearn的結果對比，V的第一行和降維結果的第一列正負相反了。

多運行幾次，這個現象並不一定出現。其實是通過上面svd_flip函數來實現的。

sklearn對奇異分解結果進行了一個處理，因為u_i*σ_i*v_i=(-u_i)*σ_i*(-v_i)，也就是u和v同時取反得到的結果是一樣的，而這會導致通過pca降維得到不一樣的結果（雖然都是正確的）。為了追求唯一的表示，首先定位u_i向量中絕對值最大的元素位置，如果它為負數，則u_i和v_i取反，否則不變。這部分的源碼在這里。