一、主成分分析法的思想

我們在研究某些問題時，需要處理帶有很多變量的數據，比如研究房價的影響因素，需要考慮的變量有物價水平、土地價格、利率、就業率、城市化率等。變量和數據很多，但是可能存在噪音和冗余，因為這些變量中有些是相關的，那麽就可以從相關的變量中選擇一個，或者將幾個變量綜合為一個變量，作為代表。用少數變量來代表所有的變量，用來解釋所要研究的問題，就能從化繁為簡，抓住關鍵，這也就是降維的思想。

主成分分析法（Principal Component Analysis，PCA）就是一種運用線性代數的知識來進行數據降維的方法，它將多個變量轉換為少數幾個不相關的綜合變量來比較全面地反映整個數據集。這是因為數據集中的原始變量之間存在一定的相關關系，可用較少的綜合變量來綜合各原始變量之間的信息。這些綜合變量稱為主成分，各主成分之間彼此不相關，即所代表的的信息不重疊。

那麽主成分分析法是如何降維的呢？我們從坐標變換的角度來獲得一個感性的認識。

我們先從最簡單的情形開始，假定數據集中的原始變量只有兩個，即數據是二維的，每個觀測值都用標準的X-y坐標軸來表示。如果每一個維度都服從正態分布（這比較常見），那麽這些數據就會形成橢圓形狀的點陣。如下圖所示，橢圓有一個長軸和一個短軸，二者是垂直的。

在短軸上，觀測點數據的變化比較小，如果把這些點垂直地投影到短軸上，那麽有很多點的投影會重合，這相當於很多數據點的信息都沒有被充分利用到；而在長軸上，觀測點的數據變化比較大。因此，如果坐標軸和橢圓的長短軸平行，那麽代表長軸的變量直接可以從數據集的原始變量中找到，它描述了數據的主要變化，而另一個原始變量就代表短軸的變量，描述的是數據的次要變化。在極端情況下，短軸退化成一個點，那麽就只能用長軸的變量來解釋數據點的所有變化，就可以把二維數據降至一維。

但是，坐標軸通常並不和橢圓的長短軸平行，就像上圖所展示的那樣。因此，需要構建新的坐標系，使得新坐標系的坐標軸與橢圓的長短軸重合或平行。這需要用到坐標變換，把觀測點在原坐標軸的坐標轉換到新坐標系下，同時也把原始變量轉換為了長軸的變量和短軸的變量，這種轉換是通過對原始變量進行線性組合的方式而完成的。

比如一個觀測點在原X-y坐標系中的坐標為(4,5)，坐標基為（1,0）和（0,1），如果長軸為斜率是1的線，短軸為斜率是-1的線，新坐標系以長軸和短軸作為坐標軸，那麽新坐標基可以取為（1/√2, 1/√2）和（-1/√2, 1/√2）。我們把兩個坐標基按行放置，作為變換矩陣，乘以原坐標，也就是對原坐標進行線性組合，可以得到該點在新坐標系下的坐標。可以看到變換後長軸變量的值遠大於短軸變量的值。

如果長軸變量解釋了數據集中的大部分變化，那麽就可以用長軸變量來代表原來的兩個變量，從而把二維數據降至一維。橢圓的長軸和短軸的長度相差越大，這種做法的效果也就越好。

接著我們把二維變量推廣到多維變量，具有多維變量的數據集其觀測點的形狀類似於一個高維橢球，同樣的，把高維橢球的軸都找出來，再把代表數據大部分信息的k個最長的軸作為新變量（相互垂直），也就是k個主成分，那麽主成分分析就完成了。

選擇的主成分越少，越能體現降維二字的內涵，可是不可避免會舍棄越多的信息。因此以什麽標準來決定我們應該選幾個主成分呢？

二、主成分的選擇問題

到這裏，我們應該有三個問題需要思考：一是進行坐標變換的矩陣是怎麽得到的呢？二是用什麽指標來衡量一個主成分所能解釋的數據變化的大小？三是以什麽標準來決定選多少個主成分呢？

首先來解決第二和第三個問題。

假定我們有m個觀測值，每個觀測值有n個特征（變量），那麽將其按列排成n行m列的矩陣，並且每一行都減去該行的均值，得到矩陣X（減去均值是為了下面方便求方差和協方差）。並按行把X整理成n個行向量的形式，即用X₁, X₂, ..., X_n來表示n個原始變量。

第一部分的例子說明了通過一個n×n的轉換矩陣對數據集中的原始變量進行線性組合，就可以得到n個新的變量。轉換矩陣可以有很多個，也就是變換的坐標系有很多個，但是只有一個可以由原始變量得到主成分。我們先不管這個獨特的矩陣是怎麽得到的，假定我們已經得到了這個轉換矩陣P，那麽把轉換後的n個主成分記為Y₁, Y₂, ..., Y_n，那麽由Y=PX，就可以得到主成分矩陣：

這n個行向量都是主成分，彼此之間是線性無關的，按照對數據變化解釋力的強度降序排列（並非被挑出來的前k個行向量才叫做主成分）。

那麽如何衡量每一個主成分所能解釋的數據變化的大小呢？

我們先看n=2時，主成分為Y₁和Y₂，原變量為X₁和X₂。從下圖可見Y₁為長軸變量，數據沿著這條軸的分布比較分散，數據的變化比較大，因此可以用Y₁作為第一主成分來替代X₁和X₂。那用什麽指標來量化數據的變化和分散程度呢？用方差！

我們把向量X₁和X₂的元素記為x_1t、x_2t（t=1,2,...,m），把主成分Y₁和Y₂的元素記為y_1t、y_2t（t=1,2,...,m），那麽整個數據集上的方差可以如下表示（數據早已經減去均值，所以行向量的均值為0）。第一主成分Y₁所能解釋的數據的變化，可以用在主成分的方差來衡量，也就是。也可以用主成分的方差占總體方差的比重來衡量，這裏假設為85%，這個比例越大，則反映的信息越多。

我們回到有n個原始變量和n個主成分的例子，那麽選擇合適的轉換矩陣P來計算得到主成分矩陣Y時，要讓單個主成分在數據集上的方差盡可能大。那麽選擇主成分的第一個一般標準是少數k個主成分（1≤k<n）的方差占數據集總體方差的比例超過85%。

於是我們初步解決了第二個問題和第三個問題，也就是如果已知轉換矩陣P和主成分矩陣Y，那麽就用一個主成分的方差占數據集總體方差的比例，來衡量該主成分能解釋的數據集方差的大小，然後按這個比例從大到小進行排序，並進行累加，如果到第k個主成分時，累加的比例恰好等於或者超過85%，那麽就選擇這k個主成分作為新變量，對數據集進行降維。

接下來問題倒回至第一個問題，也就是求解第二個問題和第三個問題的前提：轉換矩陣P怎麽算出來？

三、求解轉換矩陣和主成分矩陣

前面我們說了主成分矩陣Y的一個特點是，單個主成分向量Y_i的方差占總體方差的比例盡可能大，而且按照方差占比的大小，對所有的主成分進行降序排列。另外還有一個特性是所有的主成分都是線性無關的，或者說是正交的，那麽所有主成分中，任意兩個主成分Y_i和Y_j的協方差都是0。

第一個特點涉及到主成分的方差，第二個特點涉及到主成分之間的協方差，這自然而然讓我們想到協方差矩陣的概念，因為主成分矩陣Y的協方差矩陣的對角元素，就是每個主成分的方差，而非對角元素就是協方差。由於協方差為0，那麽主成分矩陣的協方差矩陣為一個對角矩陣，且對角元素是降序排列的！

由於數據集已經減去了均值，那麽同樣，主成分矩陣中的行向量也是0均值的，於是某兩個主成分的協方差為;

進一步得到主成分矩陣Y的協方差矩陣為：

那知道了主成分矩陣Y的協方差矩陣是對角矩陣，對於我們求出轉換矩陣P和主成分矩陣有什麽用呢？

有的，我們把Y=PX這個等式代入協方差矩陣中進行變換，就把已知的數據X和需要求的P都放到了協方差矩陣中：

比較神奇的是，主成分矩陣Y的協方差矩陣可以由數據集X的協方差矩陣得到。

數據集X的協方差矩陣顯然是一個實對稱矩陣，實對稱矩陣有一系列好用的性質：

1、n階實對稱矩陣A必然可以對角化，而且相似對角陣的對角元素都是矩陣的特征值；

2、n階實對稱矩陣A的不同特征值對應的特征向量是正交的（必然線性無關）；

3、n階實對稱矩陣A的某一特征值λ_k如果是k重特征根，那麽必有k個線性無關的特征向量與之對應。

因此數據集X的協方差矩陣作為n階實對稱矩陣，一定可以找到n個單位正交特征向量將其相似對角化。設這n個單位特征向量為e₁, e₂, ..., e_n，並按列組成一個矩陣：

那麽數據集X的協方差矩陣可以對角化為：

相似對角陣上的元素λ₁、λ₂、... 、λ_n是協方差矩陣的特征值（可能存在多重特征值），E中對應位置的列向量是特征值對應的單位特征向量。

接下來是高能時刻。我們把這個對角陣Λ上的元素從大到小降序排列，相應的把單位特征向量矩陣E裏的特征向量也進行排列。我們假設上面已經是排列好之後的形式了，那麽由於主成分矩陣的協方差矩陣也是元素從大到小降序排列的對角矩陣，那麽就可以得到：

也就是取X的協方差矩陣的單位特征向量矩陣E，用它的轉置E^T來作為轉換矩陣P，而X的協方差矩陣的特征值λ就是各主成分的方差！有了轉換矩陣P，那麽由PX我們自然就可以得到主成分矩陣Y。如果我們想把數據從n維降至k維，那麽從P中挑出前k個行向量，去乘以數據集X就行，就可以得到前k個主成分。

至此第一個問題，也就是轉換矩陣P和主成分矩陣的求解就可以完成了。

四、主成分的方差貢獻率和累計方差貢獻率

我們來拆細了看各主成分是怎麽得到的。主成分可以由協方差矩陣的單位特征向量和原始變量進行線性組合得到。

P₁就是由，X的協方差矩陣最大特征根λ₁的單位特征向量e₁轉置而成（列向量變為行向量），於是第一主成分就是：

第一主成分的方差是最大的。然後第二主成分滿足：(1)和第一主成分正交，(2)在剩余的其他主成分中，方差最大，表達式為：

同理，第k個主成分的表達式為：

我們知道用主成分的方差來衡量其所能解釋的數據集的方差，而主成分的方差就是X的協方差矩陣的特征值λ，所以第k個主成分的方差就是λ_k。我們來定義一個指標，叫做主成分Y_k的方差貢獻率，它是第k個主成分的方差占總方差的比例：

那麽前k個主成分的方差累計貢獻率為：

如果前k個主成分的方差累計貢獻率超過了85%，那麽說明用前k個主成分去代替原來的n個變量後，不能解釋的方差不足15%，沒有損失太多信息，於是我們可以把n個變量減少為k個變量，達到降維的目的。

五、主成分分析法的流程總結

我們為了推導出主成分分析法的線性代數解法，鋪墊了很多，但推導出的結果卻是相當簡潔漂亮。現在我們省略中間的過程，看主成分分析法的計算流程。

假設我們拿到了一份數據集，有m個樣本，每個樣本由n個特征（變量）來描述，那麽我們可以按照以下的步驟進行降維：

1、將數據集中的每個樣本作為列向量，按列排列構成一個n行m列的矩陣；

2、將矩陣的每一個行向量（每個變量）都減去該行向量的均值，從而使得新行向量的均值為0，得到新的數據集矩陣X；

3、求X的協方差矩陣，並求出協方差矩陣的特征值λ和單位特征向量e；

4、按照特征值從大到小的順序，將單位特征向量排列成矩陣，得到轉換矩陣P，並按PX計算出主成分矩陣；

5、用特征值計算方差貢獻率和方差累計貢獻率，取方差累計貢獻率超過85%的前k個主成分，或者想降至特定的k維，直接取前k個主成分。

六、主成分分析法計算的案例

為了更好地掌握主成分分析法的計算過程，我們來看一個例子。

假設我們想研究上海、北京房地產指數與其他價格指數之間的關系，設定了4個變量，如下表所示。

樣本數據取自1997年1月～2000年6月的統計資料，時間跨度為42個月，因此樣本容量為m=42，為了簡單起見，數據就不展示了。

第一步：計算數據集的協方差矩陣

將每個樣本作為列向量構成一個矩陣，並對矩陣的每一個行向量進行0均值化，得到了4行42列的數據集矩陣X。我們直接由X得到其協方差矩陣：

第二步：計算協方差矩陣的特征值和單位特征向量

我們用numpy來計算，代碼如下：

import numpy as np
from numpy import linalg 

# 協方差矩陣
C = [[1,-0.339,0.444,0.525],
     [-0.339,1,0.076,-0.374],
     [0.444,0.076,1,0.853],
     [0.525,-0.374,0.853,1]]
# 計算特征值和特征向量
value,vector = linalg.eig(C)
print(‘特征值為：‘,np.round(value,4),‘\n‘)
for i in range(4):
    print(‘特征值‘,np.round(value[i],4),‘對應的特征向量為:\n‘,np.round(vector[:,i].T,4),‘\n‘)

# 求每一列的L2範數,如果都是1，則已經單位化了。   
print(‘特征向量已經是單位特征向量了：‘,linalg.norm(vector,ord=2,axis=0))

特征值為： [2.3326 1.0899 0.5399 0.0376] 

特征值 2.3326 對應的特征向量為:
 [ 0.4947 -0.2687  0.5464  0.6201] 

特征值 1.0899 對應的特征向量為:
 [-0.2019  0.8378  0.5004  0.0832] 

特征值 0.5399 對應的特征向量為:
 [-0.844  -0.3399  0.1652  0.3805] 

特征值 0.0376 對應的特征向量為:
 [ 0.0458  0.3322 -0.651   0.681 ] 

特征向量已經是單位特征向量了： [1. 1. 1. 1.]

得到特征值是λ₁=2.3326 ，λ₂=1.0899 ，λ₃=0.5399 ，λ₄=0.0376，已經是從大到小排列好的了。而且特征向量已經是單位特征向量了。

第三步：得到轉換矩陣P和主成分矩陣Y

我們得到第一個主成分如下，也就是用最大特征值的特征向量對原始變量進行線性組合。

其他三個主成分同樣可以得到。

第四步：計算主成分的方差貢獻率和累計方差貢獻率，選擇k個主成分

有了協方差矩陣的特征值，計算就非常簡單了。

# 方差貢獻率
contrib_rate = value/sum(value)
print(‘方差貢獻率為：‘,np.round(contrib_rate,4))

# 累計方差貢獻率
cum_contrib_rate = np.cumsum(contrib_rate)
print(‘\n累計方差貢獻率為：‘,np.round(cum_contrib_rate,4))

方差貢獻率為： [0.5831 0.2725 0.135  0.0094]

累計方差貢獻率為： [0.5831 0.8556 0.9906 1.    ]

得到的結果整理如下。可以看到第一主成分Y₁和第二主成分Y₂的累積方差貢獻率已經達到了85.56%，可以認為用來代替4個原始變量，也不會造成太多信息損失。

七、基於投影方差最大化的數學推導

不要不耐煩，數學還是很有意思的，哈哈。我們下面用其他的方法來推導轉換矩陣和主成分的計算公式，可以把主成分的求解問題轉換為一個約束條件下的求最大值的問題。

假設數據集X有m個樣本，每個樣本是一個n維的列向量，我們把X整理成n行m列的矩陣：X=(X₁, X₂,..., X_n)^T，且已經對行向量進行了0均值化。現在我們希望用主成分分析法將n維變量降至k維。首先進行坐標變換，經過坐標變換後的新坐標系為W=｛w₁, w₂, ..., w_n｝，其中w_i是標準正交基，W^TW是單位向量。如果第一主成分Y₁的方向就是w₁這條坐標軸的方向，那麽樣本投影到w₁上之後會被廣泛散布，使得樣本之間的差別變得特別明顯，也就是投影的方差最大。

設數據集X在w₁上的投影為z₁=w₁^TX，那麽問題就成了希望在w₁的L2範數為1的約束條件下，尋找向量w₁，使得投影的方差最大化。記數據集X的協方差矩陣為Cov(X)=Σ，則投影方差最大化問題為：

寫成拉格朗日問題：

對w₁求導並令其為0，得到如下表達式。Σ是數據集X的協方差矩陣，所以w₁可以看做是協方差矩陣的一個特征值λ的特征向量。

對於以上的式子，等式左右兩邊都左乘一個w₁^T，得到數據集X在w₁上投影的方差，也就是特征值λ。由於w₁是第一主成分Y₁所在的坐標軸，那麽由方差最大得到λ是最大的特征值。

第一主成分怎麽求出來呢？很簡單，Y₁=w₁^TX（這裏的w₁是特指方差最大化的解）就是了。

求出了第一主成分，我們可以再求第二主成分。假設第二主成分Y₂在新坐標軸w₂的方向上，那麽Y₂應該是剩余的主成分中，使數據集在w₂上投影的方差最大的那個。

數據集X在w₂上的投影為z₂=w₂^TX，除了滿足w₂的L2範數為1的條件外，還需要滿足w₂^Tw₁=0，其中w₁是我們已經求出來的。於是投影方差最大化問題寫成拉格朗日的格式為：

對w₂求導，經過一系列的推導，我們最終可以得到Σw₂=λw_2。

那麽w₂可以看做是協方差矩陣Σ的特征向量，對應的特征值為第二大特征值λ₂，第二主成分Y₂=w₂X。

類似的，其他主成分所在的坐標軸的標準正交基w_i是依次遞減的特征值所對應的單位特征向量。

看到最後的人，祝你五一勞動節快樂！！！身體健康，工作順利！

參考資料

1、《PCA數學原理》：http://www.360doc.com/content/13/1124/02/9482_331688889.shtml

2、主成分分析（PCA）原理總結

降維之主成分分析法（PCA）