如果無法理解一種算法的原理，那么就整理出它的流程，再在直接的計算過程中慢慢體會它的思想。

– 我說的。

許多問題在一開始都很難理解其原理。像PCA這種經典的算法，其實從頭至尾的每個步驟都非常嚴密，細節的知識特別多。但最近看到不少的書、網上的資料都是在通篇扯原理，講道理，到最后就只列一個調包的程序，讀完之后還是總覺得不清楚它到底是怎么玩的，這就沒多大意思了。

理解原理是進行深入研究的必要前提，但要做到這一點的確很難。但這個其實不是我們的問題，而是因為經典的方法通常是被許多優秀的科學家和工程師不斷完善才有了今天這副完美的樣子。作為凡人的我們一時半會理解不透是非常正常的事情。

在這種時候最好的辦法就是不去糾結過深的推導和證明（因為已經被證明是對的了），從復雜的理論中暫時抽離出來單純地梳理其操作流程，在這些具體的計算過程中再去體會那些理論本身的意義。這是一個屢試不爽的方法。在這個過程中，我們需要自覺地避開那些與我們主要目的無關的信息，梳理出最需要的技術脈絡。這自然需要的是閱讀理解能力、知識轉換能力（動手能力），長期堅持下去又可以不斷提升這些方面的本事，形成一個非常良性又高效的循環。

本文就以PCA為例，展示一遍上述思路的具體操作步驟。

1、問題描述

注意，為了和程序數據結構匹配這里使用行向量。

原特征
$$\mathbf{x}=\left(x_1,x_2,\dots ,x_d\right)$$
新特征

$$\mathbf{w}=\left(w_1,w_2,\dots ,w_d\right)$$
記樣本為：
$$\mathbf{X}=?\begin{array}{c}{\mathbf{x}}^1\\ {\mathbf{x}}^2\\ ?\\ {\mathbf{x}}^n\end{array}?$$
問題目標：

求出特征變換矩陣$A$使得：
$${\mathbf{X}}_{n\times d}{A}_{d\times d}={\mathbf{W}}_{n\times d}$$

2、計算過程

Step 1: 計算原特征協方差矩陣
$$C=(\mathbf{X}?\mathbf{u}{)}^T(\mathbf{X}?\mathbf{u})$$
Step 2: 計算協方差矩陣特征分解:
$$C\xi =\lambda \xi$$
Step 3: 特征值排序
$${\lambda }_1\ge {\lambda }_2\ge {\lambda }_d$$
Step 4: 對應特征向量排序
$${\xi }_1,{\xi }_2,\dots ,{\xi }_d$$
Step 5: 獲得變換矩陣：
$$A=?\begin{array}{c}{\xi }_1\\ {\xi }_2\\ ?\\ {\xi }_d\end{array}?$$
Step 6: 計算特征貢獻率(第$k$個新特征的貢獻率):
$$\mathbf{r}=\frac{{\lambda }_k}{{\sum }_{i=1}^d{\lambda }_i}$$
回顧：若要計算前 $k$ 個主成分，則只需要取 $\mathbf{W}=\mathbf{X}A$ 的前 $k$ 列即可。

3、Python代碼

先上代碼，直接對照公式一步步來：

import numpy as np import scipy.linalg as lina x = np.random.rand(10,5) #隨機生成一組樣本 x -= x.mean(axis=0) # 見詳注1 C = x.T.dot(x) # 計算自協方差矩陣 lam,v= lina.eig(C) # 特征分解，v是 new_index = np.argsort(lam)[::-1] # 特征值排序，見詳注2 A = -v[:,new_index] # 得到A w = x.dot(A) # 計算變換后的特征 r = lam[new_index]/lam.sum() # 計算所有特征對應的貢獻率

測試結果（只看關鍵結果）：

w[:,:2] # 新特征的前2個>>> array([[-0.3939524518, -0.4184678305],[-0.5907434013, 0.2033346207],[-0.4585388051, -0.111367225 ],[ 0.4552495673, -0.0405062598],[-0.2335902798, -0.4260334862],[ 0.4523182771, 0.039755097 ],[ 0.0902288594, 0.1869543779],[ 0.089419155 , 0.7656098218],[ 0.7645053936, -0.3353675658],[-0.1748963144, 0.1360884499]])r # 各個特征值對應的貢獻率>>> array([0.4026073116, 0.2589988934, 0.2088275432, 0.0902665298,0.0392997221])

4、與sklearn結果對比

pca = PCA(n_components=2)pca.fit(x) # x還是最開始那個xpca.explained_variance_ratio_>>> array([0.4026073116, 0.2589988934]) # 前2個特征對應的貢獻率（完全一致）pca.transform(x) # 降維變換（完全一致）>>> array([[-0.3939524518, -0.4184678305],[-0.5907434013, 0.2033346207],[-0.4585388051, -0.111367225 ],[ 0.4552495673, -0.0405062598],[-0.2335902798, -0.4260334862],[ 0.4523182771, 0.039755097 ],[ 0.0902288594, 0.1869543779],[ 0.089419155 , 0.7656098218],[ 0.7645053936, -0.3353675658],[-0.1748963144, 0.1360884499]])

再仔細計算一下其相差的值：

w1 = w[:,0]w2 = pca.transform(c)[:,0] ((w1-w2)**2).sum() >>> 2.2980196096428498e-30

相差只有1e-30的量級，注意這里是平方和，那么完全可以視作只是一些系統誤差（包括最后一位的舍入誤差）造成的。

再看r的值：

pca.explained_variance_ratio_ - r[:2] >>> array([ 1.1102230246e-16, -1.6653345369e-16])

仍然只有1e-16的量級。

因此上述方法和sklearn中的方法完全一致。

5、詳注

• 詳注1：x -= x.mean(axis=0)

這里x.mean(axis=0) 表示求出x中每列的平均值，返回一個一維數組。這里之所以可以讓不同形狀的數組做減法是用到了python自帶的broadcasting機制（廣播機制），它會自動將一維數組擴充至二維，使其變成每一行都完全相同的矩陣，行數與x相同。因此不需要單獨再做其它處理。

• 詳注2:

這里首先要明白一點的是，利用eig函數得到的特征值和特征向量是未排序的，同時特征值與特征向量的位置對應，也就是說特征值的下標與特征向量第二維的下標（也可以理解成列標）是完全一樣的。因此只需要對特征值進行排序，同時獲得其降序排列后新的下標位置即可。

其次，np.argsort只返回按升序排列得到的新下標集。那么如果對新下標集進行反向處理，即可得到降序排列的新下標。（有點繞，但這個邏輯是層次清楚的）。而符號[::-1]正好是一維數組倒序，因此所有操作放在一起就有了這段代碼new_index = np.argsort(lam)[::-1]。

總結

以上是生活随笔為你收集整理的PCA计算流程详解与实现（Python详细编码，全部测试正确，与sklearn完全一致，只有7行代码）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。