聲明：版權所有，轉載請聯系作者并注明出處:
http://blog.csdn.net/u013719780?viewmode=contents

知乎專欄：
https://www.zhihu.com/people/feng-xue-ye-gui-zi

前面兩篇文章詳細講解了線性判別分析LDA，說到LDA，就不能不提到主成份分析，簡稱為PCA，是一種非監督學習算法，經常被用來進行數據降維、有損數據壓縮、特征抽取、數據可視化（Jolliffe, 2002）。它也被稱為Karhunen-Loève變換。

1. PCA原理

PCA的思想是將n

維特征映射到k維空間上k<n，這k維特征是全新的正交特征，是重新構造出來的k維特征，而不是簡單地從n維特征中去除其余n?k

維特征。那么如何衡量投影向量的優劣呢？在數學上有三種方法衡量投影的優劣！PCA可以被定義為數據在低維線性空間上的正交投影，這個線性空間被稱為主?空間（principal subspace），使得投影數據的?差被最?化（Hotelling, 1933），即最大方差理論。等價地，它也可以被定義為使得平均投影代價最?的線性投影，即最小誤差理論。平均投影代價是指數據點和它們的投影之間的平均平?距離（Pearson, 1901）。還有另一個理論也可以解釋PCA原理，即坐標軸相關度理論。這里簡單探討前兩種，最后一種在討論PCA意義時簡單概述。

1.1 最大方差理論

在信號處理中認為信號具有較大的方差，噪聲有較小的方差，信噪比就是信號與噪聲的方差比，越大越好。因此我們認為，最好的k

維特征是將n維樣本點變換為k

維后，每一維上的樣本方差都盡可能的大。

首先，考慮在一維空間 (k=1) 上的投影。我們可以使用 n

維向量u定義這個空間的方
向。為了方便(并且不失一般性)，我們假定選擇一個單位向量，從而 uTu=1 (注意，我們只對u的方向感興趣，而對 u

本身的大小不感興趣)。

如上圖所示，紅色點表示原樣本點x(i)

，u是藍色直線的斜率也是直線的方向向量，而且是單位向量，直線上的藍色點表示原樣本點x(i)在u上的投影。容易知道投影點離原點的距離是x(i)Tu

，由于這些原始樣本點的每一維特征均值都為0，因此投影到u上的樣本點的均值仍然是0。

假設原始數據集為Xmxn

，我們的目標是找到最佳的投影空間Wnxk=(w1,w2,…,wk)，其中wi是單位向量且wi與wj(i≠j)正交， 何為最佳的W？就是原始樣本點投影到W

上之后，使得投影后的樣本點方差最大。

由于投影后均值為0，因此投影后的總方差為：

?

1m∑i=1m(x(i)Tw)2=1m∑i=1mwTx(i)x(i)Tw=∑i=1mwT(1mx(i)x(i)T)w1

?

1mx(i)x(i)T

是不是似曾相識，沒錯，它就是原始數據集X的協方差矩陣(因為x(i)的均值為0，因為無偏估計的原因，一般協方差矩陣除以m?1

，這里用m)。

記λ=1m∑mi=1(x(i)Tw)2

, ∑=1mx(i)x(i)T

, 則有

?

λ=wT∑w

.

?

上式兩邊同時左乘w

，注意到wwT=1

(單位向量)，則有

?

λw=∑w

.

?

所以w

是矩陣∑

的特征值所對應的特征向量。

欲使投影后的總方差最大，即λ

最大，因此最佳的投影向量w是特征值λ最大時對應的特征向量，因此，當我們將w設置為與具有最大的特征值λ

的特征向量相等時，方差會達到最大值。這個特征向量被稱為第一主成分。

我們可以用一種增量的方式定義額外的主成分，方法為：在所有與那些已經考慮過的方向正交的所有可能的方向中，將新的方向選擇為最大化投影方差的方向。如果我們考慮 k

維投影空間的一般情形，那么最大化投影數據方差的最優線性投影由數據協方差矩陣 ∑ 的 k 個特征向量 w1,...,wk 定義，對應于 k 個最大的特征值 λ1,...,λk

。可以通過歸納法很容易地證明出來。

因此，我們只需要對協方差矩陣進行特征值分解，得到的前k

大特征值對應的特征向量就是最佳的k維新特征，而且這k維新特征是正交的。得到前k個u以后，原始數據集X

通過變換可以得到新的樣本。

PCA算法流程

算法輸入：數據集Xmxn

* 按列計算數據集X的均值Xmean，然后令Xnew=X?Xmean；
* 求解矩陣Xnew的協方差矩陣，并將其記為Cov；
* 計算協方差矩陣COv的特征值和相應的特征向量；
* 將特征值按照從大到小的排序，選擇其中最大的k個，然后將其對應的k個特征向量分別作為列向量組成特征向量矩陣Wnxk;
* 計算XnewW，即將數據集Xnew投影到選取的特征向量上，這樣就得到了我們需要的已經降維的數據集XnewW

。

注意，計算一個nxn

矩陣的完整的特征向量分解的時間復雜度為 O(n3) 。如果我們將數據集投影到前 k 個主成分中，那么我們只需尋找前 k 個特征值和特征向量。這可以使用更高效的方法得到,例如冪方法(power method) (Golub and Van Loan, 1996)，它的時間復雜度為 O(kn2)

，或者我們也可以使用 EM 算法。

1.2 最小平方誤差理論

如上圖所示，假設有這樣的二維樣本點（紅色點），按照前文我們講解的最大方差理論，我們的目標是是求一條直線，使得樣本點投影到直線或者平面上的點的方差最大。本質是求直線或者平面，那么度量直線求的好不好，不僅僅只有方差最大化的方法。再回想我們最開始學習的線性回歸等，目的也是求一個線性函數使得直線能夠最佳擬合樣本點，那么我們能不能認為最佳的直線就是回歸后的直線呢？回歸時我們的最小二乘法度量的是樣本點到直線的坐標軸距離。比如這個問題中，特征是x，類標簽是y。回歸時最小二乘法度量的是距離d。如果使用回歸方法來度量最佳直線，那么就是直接在原始樣本上做回歸了，跟特征選擇就沒什么關系了。

因此，我們打算選用另外一種評價直線好壞的方法，使用點到直線的距離d′

來度量。

現在有m

個樣本點x(1),...,x(m)，每個樣本點為n維。將樣本點x(i)在直線上的投影記為x(1)′

，那么我們就是要最小化

?

∑i=1m(x(i)′?x(i))2

?

這個公式稱作最小平方誤差（Least Squared Error）。

初中我們就已經知道確定一條直線，只需要知道直線經過某一個點和其方向即可。

首先，我們確定直線經過的點，假設要在空間中找一點x0

來代表這m個樣本點，“代表”這個詞不是量化的，因此要量化的話，我們就是要找一個n維的點x0

，使得

?

J0(x0)=∑i=1m(x0?x(i))2

?

最小。其中J0(x0)

是平方錯誤評價函數(squared-error criterion function)，假設xˉ為m

個樣本點的均值，即

?

xˉ=1m∑i=1mx(i)

?

則

?

J0(x0)=∑i=1m(x0?x(i))2=∑i=1m((x0?xˉ)?(x(i)?xˉ))2=∑i=1m(x0?xˉ)2?2∑i=1m(x0?xˉ)T(x(i)?xˉ)+∑i=1m(x(i)?xˉ)2=∑i=1m(x0?xˉ)2?2(x0?xˉ)T∑i=1m(x(i)?xˉ)+∑i=1m(x(i)?xˉ)2=∑i=1m(x0?xˉ)2+∑i=1m(x(i)?xˉ)2

?

顯然，上式的第二項與x0

無關，因此，J0(x0)在xˉ

處有最小值。

接下來，我們確定直線的方向向量。我們已經知道直線經過點xˉ

，假設直線的方向是單位向量e??。那么直線上任意一點x(i)′

有：

?

x(i)′=xˉ+aie??

?

其中，ai

是x(i)′到點xˉ

的距離。

我們重新定義最小平方誤差：

?

J1(a1,a2,...,am,e??)=∑i=1m(x(i)′?x(i))2=∑i=1m((xˉ+aie??)?x(i))2=∑i=1m(aie???(x(i)?xˉ))2=∑i=1ma2ie??2?2∑i=1maie??T(x(i)?xˉ)+∑i=1m(x(i)?xˉ)2

?

我們首先固定e??

，將其看做是常量，然后令J1 關于 ai

的導數等于0，則有：

?

ai=e??T(x(i)?xˉ),

?

這個結果意思是說，如果知道了e??

，那么將(x(i)?xˉ)與e??做內積，就可以知道了x(i)在e??上的投影離xˉ

的長度距離，不過這個結果不用求都知道。

然后是固定ai

，對e??求偏導數，我們先將 ai代入J1

，得

?

J1(e??)=∑i=1ma2ie??2?2∑i=1maie??T(x(i)?xˉ)+∑i=1m(x(i)?xˉ)2=∑i=1ma2ie??2?2∑i=1ma2i+∑i=1m(x(i)?xˉ)2=?∑i=1m(e??T(x(i)?xˉ))2+∑i=1m(x(i)?xˉ)2=?∑i=1me??T(x(i)?xˉ)(x(i)?xˉ)Te??+∑i=1m(x(i)?xˉ)2=?e??TSe??+∑i=1m(x(i)?xˉ)2

?

其中S=∑mi=1(x(i)?xˉ)(x(i)?xˉ)T

，與協方差矩陣類似，只是缺少個分母n?1

，我們稱之為散列矩陣(scatter matrix)。

現在我們就可以用拉格朗日乘數法求解方向向量e??

了。令

?

f(e??)=?e??TSe??+∑i=1m(x(i)?xˉ)2+λ(e??Te???1)

?

令上式關于e??

的偏導數等于0，則可得

?

Se??=λe??

?

兩邊除以n?1

就變成了對協方差矩陣求特征值向量了。

從不同的思路出發，最后得到同一個結果，對協方差矩陣求特征向量，求得后特征向量上就成為了新的坐標，如下圖：

這時候點都聚集在新的坐標軸周圍，因為我們使用的最小平方誤差的意義就在此。另外，PRML書上從線性子空間的角度進行了詳細的闡述，有興趣的讀者可以看看。

2. PCA算法優缺點：

優點：

• 它是無監督學習，完全無參數限制的。在PCA的計算過程中完全不需要人為的設定參數或是根據任何經驗模型對計算進行干預，最后的結果只與數據相關，與用戶是獨立的。

• 用PCA技術可以對數據進行降維，同時對新求出的“主元”向量的重要性進行排序，根據需要取前面最重要的部分，將后面的維數省去，可以達到降維從而簡化模型或是對數據進行壓縮的效果。同時最大程度的保持了原有數據的信息。

• 各主成分之間正交，可消除原始數據成分間的相互影響。

• 計算方法簡單，易于在計算機上實現。

缺點：

• 如果用戶對觀測對象有一定的先驗知識，掌握了數據的一些特征，卻無法通過參數化等方法對處理過程進行干預，可能會得不到預期的效果，效率也不高。

• 貢獻率小的主成分往往可能含有對樣本差異的重要信息。

• 特征值矩陣的正交向量空間是否唯一有待討論。

• 在非高斯分布的情況下，PCA方法得出的主元可能并不是最優的，此時在尋找主元時不能將方差作為衡量重要性的標準。

  3. 代碼實現

from __future__ import print_function from sklearn import datasets import matplotlib.pyplot as plt import matplotlib.cm as cmx import matplotlib.colors as colors import numpy as np %matplotlib inlinedef shuffle_data(X, y, seed=None):if seed:np.random.seed(seed)idx = np.arange(X.shape[0])np.random.shuffle(idx)return X[idx], y[idx]# 正規化數據集 X def normalize(X, axis=-1, p=2):lp_norm = np.atleast_1d(np.linalg.norm(X, p, axis))lp_norm[lp_norm == 0] = 1return X / np.expand_dims(lp_norm, axis)# 標準化數據集 X def standardize(X):X_std = np.zeros(X.shape)mean = X.mean(axis=0)std = X.std(axis=0)# 做除法運算時請永遠記住分母不能等于0的情形# X_std = (X - X.mean(axis=0)) / X.std(axis=0) for col in range(np.shape(X)[1]):if std[col]:X_std[:, col] = (X_std[:, col] - mean[col]) / std[col]return X_std# 劃分數據集為訓練集和測試集 def train_test_split(X, y, test_size=0.2, shuffle=True, seed=None):if shuffle:X, y = shuffle_data(X, y, seed)n_train_samples = int(X.shape[0] * (1-test_size))x_train, x_test = X[:n_train_samples], X[n_train_samples:]y_train, y_test = y[:n_train_samples], y[n_train_samples:]return x_train, x_test, y_train, y_test# 計算矩陣X的協方差矩陣 def calculate_covariance_matrix(X, Y=np.empty((0,0))):if not Y.any():Y = Xn_samples = np.shape(X)[0]covariance_matrix = (1 / (n_samples-1)) * (X - X.mean(axis=0)).T.dot(Y - Y.mean(axis=0))return np.array(covariance_matrix, dtype=float)# 計算數據集X每列的方差 def calculate_variance(X):n_samples = np.shape(X)[0]variance = (1 / n_samples) * np.diag((X - X.mean(axis=0)).T.dot(X - X.mean(axis=0)))return variance# 計算數據集X每列的標準差 def calculate_std_dev(X):std_dev = np.sqrt(calculate_variance(X))return std_dev# 計算相關系數矩陣 def calculate_correlation_matrix(X, Y=np.empty([0])):# 先計算協方差矩陣covariance_matrix = calculate_covariance_matrix(X, Y)# 計算X, Y的標準差std_dev_X = np.expand_dims(calculate_std_dev(X), 1)std_dev_y = np.expand_dims(calculate_std_dev(Y), 1)correlation_matrix = np.divide(covariance_matrix, std_dev_X.dot(std_dev_y.T))return np.array(correlation_matrix, dtype=float)class PCA():"""主成份分析算法PCA，非監督學習算法."""def __init__(self):self.eigen_values = Noneself.eigen_vectors = Noneself.k = 2def transform(self, X):""" 將原始數據集X通過PCA進行降維"""covariance = calculate_covariance_matrix(X)# 求解特征值和特征向量self.eigen_values, self.eigen_vectors = np.linalg.eig(covariance)# 將特征值從大到小進行排序，注意特征向量是按列排的，即self.eigen_vectors第k列是self.eigen_values中第k個特征值對應的特征向量idx = self.eigen_values.argsort()[::-1]eigenvalues = self.eigen_values[idx][:self.k]eigenvectors = self.eigen_vectors[:, idx][:, :self.k]# 將原始數據集X映射到低維空間X_transformed = X.dot(eigenvectors)return X_transformeddef main():# Load the datasetdata = datasets.load_iris()X = data.datay = data.target# 將數據集X映射到低維空間X_trans = PCA().transform(X)x1 = X_trans[:, 0]x2 = X_trans[:, 1]cmap = plt.get_cmap('viridis')colors = [cmap(i) for i in np.linspace(0, 1, len(np.unique(y)))]class_distr = []# Plot the different class distributionsfor i, l in enumerate(np.unique(y)):_x1 = x1[y == l]_x2 = x2[y == l]_y = y[y == l]class_distr.append(plt.scatter(_x1, _x2, color=colors[i]))# Add a legendplt.legend(class_distr, y, loc=1)# Axis labelsplt.xlabel('Principal Component 1')plt.ylabel('Principal Component 2')plt.show()if __name__ == "__main__":main()

版權聲明：本文為博主原創文章，未經博主允許不得轉載。 https://blog.csdn.net/u013719780/article/details/78352262

總結

以上是生活随笔為你收集整理的从零开始实现主成分分析(PCA)算法的全部內容，希望文章能夠幫你解決所遇到的問題。

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