Python機器學習筆記：主成分分析（PCA）算法

一：引入問題

　　首先看一個表格，下表是某些學生的語文，數學，物理，化學成績統計：

　　首先，假設這些科目成績不相關，也就是說某一科目考多少分與其他科目沒有關系，那么如何判斷三個學生的優秀程度呢？首先我們一眼就能看出來，數學，物理，化學這三門課的成績構成了這組數據的主成分（很顯然，數學作為第一主成分，因為數據成績拉的最開）。

　　那么為什么我們能一眼看出來呢？

　　當然是我們的坐標軸選對了！！

　　下面，我們繼續看一個表格，下標是一組學生的數學，物理，化學，語文，歷史，英語成績統計：

　　那么這個表我們能一眼看出來嗎？

　　數據太多了，以至於看起來有些凌亂，無法直接看出這組數據的主成分，因為在坐標系下這組數據分布的很散亂。究其原因，是因為無法撥開遮住肉眼的迷霧，如果把這些數據在相應的空間中表示出來，也許你就能換一個觀察角度找出主成分，如下圖所示：

　　簡單的二維或者三維我們可以想象出來其分布狀態，那么對於更高維的數據，能想象出來其分布嗎？還有，就算能描述分布，如何精確地找到這些主成分的軸？如何衡量你提取的主成分到底占了整個數據的多少信息？所以，我們就要用到主成分分析的處理方法。

　　為了說明什么是數據的主成分，我們首先得了解數據降維，數據降維時怎么回事？

二，數據降維

　　假設三維空間中有一系列點，這些點分布在一個過原點的斜面上，如果你用自然坐標x,y,z這三個軸來表示這組數據的話，需要使用三個維度，而事實上，這些點的分布僅僅是在一個二維的平面上，那么問題出在哪里？如果你仔細想想，能不能把x,y,z坐標系旋轉一下，使數據所在平面與x，y平面重合？這就對了！如果把旋轉后的坐標記為x'，y'，z'，那么這組數據的表示只用x'和y'兩個維度表示即可！

　　當然了，如果想恢復原來的表示方式，那就得把這兩個坐標之間的變換矩陣存下來。這樣就能把數據維度降下來了！但是，我們要看到這個過程的本質，如果把這些數據按行或者按類排成一個矩陣，那么這個矩陣的秩就是2！這些數據之間是有相關性的，這些數據構成的過原點的向量的最大線性無關組包含2個向量，這就是為什么一開始就假設平面過原點的原因！

　　那么如果不過原點呢？這就是數據中心化的緣故！將坐標原點平移到數據中心，這樣原本不相關的數據在這個新坐標系中就有相關性了！有趣的是，三點一定共面，也就是三維空間中任意三點中心化后都是線性相關的，一般來講n維空間中n個點一定能在一個n-1維子空間中分析！

　　總結一下這個例子，數據降維后並沒有丟棄任何東西，因為這些數據在平面以外的第三個維度的分量都為0。現在，假設這些數據在z'軸有一個很小的抖動，那么我們仍然用上述的二維表示這些數據，理由是我們可以認為這兩個軸的信息是數據的主成分，而這些信息對於我們的分析已經足夠了，z'軸上的抖動很有可能是噪音，也就是說本來這組數據是有相關性的，噪聲的引入，導致了數據不完全相關，但是，這些數據在z'軸上的分布與原點構成的夾角非常小，也就是說在z'軸上有很大的相關性，綜合考慮，就可以認為數據在x'，y'軸上的投影構成了數據的主成分！

　　所以說，降維肯定意味着信息的丟失，不過鑒於實際數據本身常常存在的相關性，我們可以想辦法在降維的同時將信息的損失盡量降低。

　　下面在說一個極端的情況，也許在現實中不會出現，但是 類似的情況還是很常見的。

　　假設某學籍數據有兩列M和F，其中M列的取值是如果此學生為男性，則取值為1，為女性則取0；而F列是學生為女性，則取值為0，男性則為1.此時如果我們統計全部學籍數據，會發現對於任何一條記錄來說，當M為1時F必定為0，反之當M為0時F必定為1，在這種情況下，我們將M或者F去掉實際上沒有任何信息的損失，因為只要保留一列就可以完全還原另一列。

　　那么降維我們差不多說清楚了，現在我們將自己面對的數據抽象為一組向量，那么下面我們有必要研究一些向量的數學性質，而這些數學性質將成為后續推導出PCA的理論基礎。

 三，PCA基本數學原理

3.1 內積與投影

　　下面先看一個向量運算：內積。

　　兩個維數相同的向量的內積被定義為：

　　內積運算將兩個向量映射為一個實數，其計算方式非常容易理解，但是其意義並不明顯，下面我們分析內積的幾何意義。假設A和B是兩個n維向量，我們知道n維向量可以等價表示為n維空間中的一條從原點發射的有向線段，為了簡單起見，我們假設A和B均為二維向量，則

　　那么在二維平面上A和B可以用兩條發自原點的有向線段表示，如下圖：

　　現在我們從A點向B所在直線引入一條垂線，我們知道垂線與B的交點叫做A在B上的投影，再假設A與B的夾角為a，則投影的矢量長度為（這里假設向量B的模為1）

　　其中是向量A的模，也就是A線段的標量長度。

　　注意這里區分了矢量長度和標量長度，標量長度總是大於等於0，值就是線段的長度；而矢量長度可能為負，其絕對值是線性長度，而符號取決於其方向與標准方向相同或者相反。

　　可能大家還沒明白內積和這個東西有什么關系，不過我們將內積表示為另一種我們熟悉的方式：

　　現在明白了點吧，A與B的內積等於A到B的投影長度乘以B的模，在進一步，如果我們假設B的模為1，即讓，那么就變成了下面公式（也就是上面我們說的投影的矢量長度）：

　　也就是說，設向量B的模為1，則A與B的內積值等於A向B所在直線投影的矢量長度！這就是內積的一種幾何解釋，也就是我們得到的第一個重要結論。

3.2 基

　　下面我們繼續在二維空間討論向量，上文說過，一個二維向量可以對應二維笛卡爾直角坐標系中從原點出發的一條有向線段，例如下面這個向量：

 　　在代數表示方面，我們經常使用線段終點的點的坐標表示向量，例如上面的向量可以表示為（3,2），這里我們再熟悉不過的向量表示。

　　不過我們常常忽略，只有一個（3,2）本身是不能夠精確表示一個向量的，我們仔細看一下，這里的坐標（3,2）實際上表示的是向量在x軸上的投影值3，在y軸上的投影值為2。也就是說我們使其隱式引入一個定義：以x軸和y軸上正方向長度為1的向量為標准，那么一個向量（3,2）實際上是說在x軸投影為3而y軸投影為2.注意投影是一個矢量，所以可以為負。

　　更正式的說，向量（x,y）實際上表示線性組合：

　　不難證明所有二維向量都可以表示為這樣的線性組合，此處（1,0）和（0，1）叫做二維空間的一組基。

　　所以，要准確描述向量，首先要確定一組基，然后給出基所在的各個直線上的投影值，就可以了，只不過我們經常省略第一步，而默認以（1,0）和（0,1）為基。

　　我們之所以默認選擇（1,0）和（0,1）為基，當然是比較方便，因為他們分別是x和y軸正方向上的單位向量，因此就使得二維平面上點坐標和向量一一對應，非常方便。但實際上任何兩個線性無關的二維向量都可以成為一組基，所謂線性無關在二維平面內可以直觀認為是兩個不再一條直線上的向量。

　　例如，（1，1）和（1，-1）也可以成為一組基。一般來說，我們希望基的模是1，因為從內積的意義可以看到，如果基的模式1，那么就可以方便的用向量點乘基而直接獲得其在新基上的坐標了！實際上，對應於任何一個向量我們總可以找到其同方向上模為1的向量，只要讓兩個分量分別除以模就好了，例如上面的基就可以變為：。

　　現在我們想獲得（3,2）在新基上的坐標，即在兩個方向上的投影矢量值，那么根據內積的幾何意義，我們只要分別計算（3,2）和兩個基的內積，不難得到新的坐標為。

　　下圖給出了新的基以及（3,2）在新基上坐標值的示意圖：

　　另外這里要注意的是，我們列舉的例子中基是正交的（即內積為0，或者說相互垂直），但是可以成為一組基的唯一要求就是線性無關，非正交的基也是可以的，不過因為正交基有較好的性質，所以一般使用的基都是正交的。

3.3 基變換的矩陣表示

　　下面我們找一種簡單的方式來表示基變換，還是拿上面的例子，想一下，將（3,2）變換為新基上的坐標，就是用（3,2）與第一個基做內積運算，作為第一個新的坐標分量，然后用（3,2）與第二個基做內積運算，作為第二個新坐標的分量。實際上，我們可以用矩陣想成的形式簡潔的表示這個變換：

　　那么其中矩陣的兩行分別為兩個基，乘以原向量，其結果剛好為新基的坐標，可以稍微擴展一下，如果我們有m個二維向量，只要將二維向量按照列排成一個兩行m列的矩陣，然后用“基矩陣”乘以這個矩陣，就得到了所有這些向量在新基下的值，例如（1,1），（2，2），（3,3）想變換到剛才那組基上，則可以變為這樣：

　　於是一組向量的基變換被表示為矩陣的相乘。

　　一般地，如果我們有M個N維向量，想將其變換為由R個N維向量表示的新空間中，那么首先將R個基按照行組成矩陣A,，然后將向量按照列組成矩陣B，那么兩個矩陣的乘積AB就是變換結果，其中AB的第m列為A中的第M列變換后的結果。

　　數學表示為：

　　其中P_{i是一個行向量，表示第i個基，aj是一個列向量，表示第j個原始數據記錄。}

　　特別要注意的是，這里R可以小於N，而R決定了變換后數據的維數，也就是說，我們可以將一個N維數據變換到更低維度的空間中去，變換后的維度取決於基的數量，因此這種矩陣相乘的表示也可以表示為降維變換。 

 　　最后，上述分析同時給矩陣相乘找到了一種物理解釋：兩個矩陣相乘的意義是將右邊矩陣中每一列列向量變換到左邊矩陣中每一行行向量為基所表示的空間中去。更抽象的說：一個矩陣可以表示為一種線性變換。

3.4 協方差矩陣及優化目標

　　上面我們討論了選擇不同的基可以對同樣一組數據給出不同的表示，而且如果基的數量少於向量的本身的維數，則可以達到降維的效果，但是我們還沒有回答最關鍵的一個問題：如何選擇基才是最優的，或者說，如果我們有一組N維向量，現在要將其降到K維（K小於N），那么我們應該如何選擇K個基才能最大程度保留原有的信息？

　　要完全數學化這個問題非常繁雜，這里我們用一個非形式化的直觀方法來看這個問題。

　　為了避免過於抽象的討論，我們仍然以一個具體的例子展開，假設我們的數據由五條記錄組成，將它們表示為矩陣形式：

　　其中每一列為一條數據記錄，而一行為一個字段，為了后續處理方便，我們首先將每個字段內所有值都減去字段均值，其結果是將每個字段都變為均值為0（這樣做的好處后面可以看到）。

　　我們看上面的數據，第一個字段的均值為2，第二個字段的均值為3，所以變換后：

　　我們可以看到五條數據在平面直角坐標系內的樣子：

　　現在問題來了：如果我們必須使用一維來表示這些數據，又希望盡量保留原始的信息，你要如何選擇？

　　通過上一節對及變換的討論我們知道，這個問題實際上是要在二維平面中選擇一個方向，將所有數據都投影到這個方向所在的直線上，用投影值表示原始記錄，這是一個實際的二維降到一維的問題。

　　那么如何選擇這個方向（或者說是基）才能盡量保留最多的原始信息呢？一種直觀的看法是：希望投影后的投影值盡可能分散。

　　以上圖為例，可以看出如果向x軸投影，那么最左邊的兩個點會重疊在一起，中間的兩個點也會重疊在一起，於是本身四個各不相同的二維點投影后只剩下兩個不同的值了，這是一種嚴重的信息丟失，同理，如果向y軸投影最上面的兩個點和分布在x軸上的兩個點也會重疊，所以看來x和y軸都不是最好的投影選擇。我們直觀目測，如果向通過第一象限和第三象限的斜線投影，則五個點在投影后還是可以區分的。

　　下面我們用數學方法表述這個問題。

3.5 方差

　　上文說道，我們希望投影后投影值盡可能分散，而這種分散程度，可以用數學上的方差來表述，此處，一個字段的方差可以看做事每個元素與字段均值的差的平方和的均值，即：

　　由於上面我們已經將每個字段的均值都化0 了，因此方差可以直接用每個元素的平方和除以元素個數表示：

　　於是上面的問題被形式化表示為：尋找一個一維基，使得所有數據變換為這個基上的坐標表示后，方差值最大。

3.6 協方差

　　對於上面二維降成一維的問題來說，找到那個使得方差最大的方向就可以了，不過對於更高維，還有一個問題需要解決，考慮三維降到二維問題，與之前相同，首先我們希望找到一個方向使得投影后方差最大，這樣就完成了第一個方向的選擇，繼而我們選擇第二個投影方向。

　　如果我們還是單純的只選擇方差最大的方向，很顯然，這個方向與第一個方向應該是“幾乎重合在一起”，顯然這樣的維度是沒有用的，因此應該有其他約束條件。從直觀上講，讓兩個字段盡可能表示更多的原始信息，我們是不希望他們之間存在線性相關性，因為相關性意味着兩個字段不是完全獨立，必然存在重復表示的信息。

 　　數字上可以用兩個字段的協方差表示其相關性，由於已經讓每個字段均值為0，則：

　　可以看出，在字段均值為0的情況下，兩個字段的協方差簡潔的表示為其內積除以元素數m。

　　當協方差為0時，表示兩個字段完全獨立，為了讓協方差為0，我們選擇第二個即時只能在與第一個基正交的方向上選擇。因此最終選擇的兩個方向一定是正交的。

　　至此，我們得到了降維問題的優化目標：將一組N維向量降維k維（K大於0，小於N），其目標是選擇K個單位（模為1）正交基，使得原始數據變換到這組基上后，各字段兩兩間協方差為0，而字段的方差則盡可能大（在正交的約束下，取最大的k個方差）。

　　然后我們用X乘以X的轉置，並乘上系數1/m：

　　這時候我們會發現，這個矩陣對角線上的兩個元素分別是兩個字段的方差，而其他元素是a和b的協方差，兩者被統一到了一個矩陣的。

　　根據矩陣相乘的運算法則，這個結論很容易被推廣到一般情況：

　　設我們有m個n維數據記錄，將其按列排成n乘m的矩陣X，設，則C是一個對稱矩陣，其對角線分別是各個字段的方差，而第l行j列和j行i列元素相同，表示i和j兩個字段的協方差。

3.7 協方差矩陣

　　上面我們導出了優化目標，但是這個目標似乎不能直接作為操作指南（或者說算法），因為它只說要什么，但是根本沒有說怎么做，所以我們要在數學上繼續研究計算方案。

　　我們看到，最終要達到的目標與字段內方差及字段間協方差有密切關系。因此我們希望能將兩者統一表示，仔細觀察發現，兩者均可以表示為內積的形式，而內積又與矩陣相乘密切相關。於是，我們來了靈感：

　　假設我們只有a和b 兩個字段，那么我們將他們按行組成矩陣X：

 

3.8 協方差矩陣對角化

　　根據上述推導，我們發現要達到優化目前等價於將協方差矩陣對角化：即除對角線外的其他元素化為0，並且在對角線上將元素按照大小從上到下排列，這樣我們就達到了優化目的，這樣說可能還不清晰，我們進一步看下原矩陣與基變換后矩陣協方差矩陣的關系：

　　設原始數據矩陣X對於的協方差矩陣為C，而P是一組基按行組成的矩陣，設Y=PX，則Y為X對P做基變換后的數據，設Y的協方差矩陣為D，我們推導一下D與C的關系：

　　現在事情很明白，我們要找的P不是別的，而是能讓原始協方差矩陣對角化的P，換句話說，優化目標變成了尋找一個矩陣P，滿足PCP^T是一個對角矩陣，並且對角元素按照從大到小依次排列，那么P的前K行就是要尋找的基，用P的前K行就是要尋找的基，用P的前K行組成的矩陣乘以X就使得X從N維降到了K維並滿足上述優化條件。

　　至此，我們離“發明”PCA還有一步之遙！

　　現在所有的焦點都聚集在了協方差矩陣對角化問題上，有時，我們真應該感謝數學家的先行，因為矩陣對角化在線性代數領域已經屬於被玩爛的東西，所以這在數學上根本不是問題。

　　由上文知道，協方差矩陣C是一個對稱矩陣，在線性代數上，實對稱矩陣有一系列非常好的性質：

　　1）實對稱矩陣不同特征值對應的特征向量必然正交。

　　2）設特征向量重數為r，則必然存在r個線性無關的特征向量對應於，因此可以將這r個特征向量單位正交化。

　　有上面兩條可知，一個n行n列的實對稱矩陣一定可以找到n個單位正交特征向量，設這n個特征向量為，我們將其按照列組成矩陣：

　　則對協方差矩陣C有如下結論：

　　其中為對稱矩陣，其對角元素為各特征向量對應的特征值（可能有重復）。

　　到這里，我們發現我們已經找到了需要的矩陣P：

　　P是協方差矩陣的特征向量單位化后按照行排列出的矩陣，其中每一行都是C的一個特征向量，如果設P按照中特征值從大到小，將特征向量從上到下排列，則用P的前K行組成的矩陣乘以原始數據矩陣X，就可以得到我們需要的降維后的數據矩陣Y。

　　至此，我們完成了整個PCA的數學原理討論。

3.9 對上面例子整合

　　1，原始數據集矩陣X：

　　2，均值為（2， 3），求均值后：

　　3，再求協方差矩陣：

　　4，特征值：

　　5，對應的特征向量：

　　6，標准化：

　　7，選擇較大特征值對應的特征向量：

　　8，執行PCA變換：Y=PX 得到的Y就是PCA降維后的值 數據集矩陣：

　　9，計算代碼：

from  numpy  import *

li = [[1,1],[1,3],[2,3],[4,4],[2,4]]
matrix = mat(li)
# 求均值
mean_matrix = mean(matrix,axis=0)
# print(mean_matrix.shape)
# 減去平均值
Dataadjust = matrix - mean_matrix
# print(Dataadjust.shape)
#計算特征值和特征向量
covMatrix = cov(Dataadjust,rowvar=0)
eigValues , eigVectors = linalg.eig(covMatrix)
# print(eigValues.shape)
# print(eigVectors.shape)
# 對特征值進行排序
eigValuesIndex = argsort(eigValues)
# print(eigValuesIndex)
# 保留前K個最大的特征值
eigValuesIndex = eigValuesIndex[:-(1000000):-1]
# print(eigValuesIndex)
# 計算出對應的特征向量
trueEigVectors = eigVectors[:,eigValuesIndex]
# print(trueEigVectors)
# 選擇較大特征值對應的特征向量
maxvector_eigval = trueEigVectors[：，0]
# print(maxvector_eigval)
# 執行PCA變換：Y=PX 得到的Y就是PCA降維后的值 數據集矩陣
pca_result = maxvector_eigval * Dataadjust.T
print(pca_result)

　　

　　10，代碼執行結果：

[[-2.12132034 -0.70710678  0.          2.12132034  0.70710678]]

　　

四，主成分分析（PCA）算法步驟

　　介紹一個PCA的教程：A tutorial on Principal Components Analysis ——Lindsay I Smith

 

　　PCA（Principal Components Analysis）即主成分分析，是一種常用的數據分析手段，是圖像處理中經常用到的降維方法。對於一組不同維度之間可能存在線性相關關系的數據，PCA能夠把這組數據通過正交變換變成各個維度之間線性無關的數據，經過PCA處理的數據中的各個樣本之間的關系往往更直觀，所以它是一種非常常用的數據分析和預處理工具。PCA處理之后的數據各個維度之間是線性無關的，通過剔除方差較小的那些維度上的數據，我們可以達到數據降維的目的。

　　PCA從原始變量出發，通過旋轉變化（即原始變量的線性組合）構建出一組新的，互不相關的新變量，這些變量盡可能多的解釋原始數據之間的差異性（即數據內在的結構），他們就成為原始數據的主成分。由於這些變量不相關，因此他們無重疊的各自解釋一部分差異性。依照每個變量解釋時差異性大小排序，他們成為第一主成分，第二主成分，以此類推。

　　主成分分析（PCA）是一種基於變量協方差矩陣對數據進行壓縮降維，去噪的有效方法，PCA的思想是將n維特征映射到k維上（k<n），這k維特征稱為主元（主成分），是舊特征的線性組合，這些線性組合最大化樣本方差，盡量使用新的k個特征互不相關。這k維是全新的正交特征，是重新構造出來的k維特征，而不是簡單地從n維特征中取出其余n-k維特征。

　　說了這么多，下面說一下PCA降維的算法步驟。

　　設有m條n維數據：

1） 將原始數據按列組成n行m列矩陣X

2）將X的每一行（代表一個屬性字段）進行零均值化（去平均值），即減去這一行的均值

3）求出協方差矩陣  

4）求出協方差矩陣的特征值及對應的特征向量

5）將特征向量按對應特征值大小從上到下按行排列成矩陣，取前k行組成矩陣P（保留最大的k各特征向量）

6）Y=PX 即為降維到K維后的數據

 

五，核主成分分析KPCA介紹

　　在上面的PCA算法中，我們假設存在一個線性的超平面，可以讓我們對數據進行投影。但是有些時候，數據不是線性的，不能直接進行PCA降維。這里就需要用到和支持向量機一樣的核函數的思想，先把數據集從 n 維映射到線性可分的高維 N>n，然后再從N維降維到一個低維度 n'，這里的維度之間滿足 n' < n< N。

　　使用了核函數的主成分分析一般稱之為核主成分分析（Kernelized  PCA，以下簡稱 KPCA。假設高維空間的數據是由 n 維空間的數據通過映射 Φ 產生）。

　　則對於 n 維空間的特征分解：

 

 　映射為：

 

 　　通過在高維空間進行協方差矩陣的特征值分解，然后用和PCA一樣的方法進行降維。一般來說，映射 Φ 不用顯式的計算，而是在需要計算的時候通過核函數完成。由於KPCA需要核函數的運算，因此它的計算量要比PCA大很多。

六，PCA算法總結

　　這里對PCA算法做一個總結，作為一個非監督學習的降維方法，它只需要特征值分解，就可以對數據進行壓縮，去噪。因此在實際場景應用很廣泛。為了克服PCA一些缺點，出現了很多PCA的變種，比如上面為解決非線性降維的KPCA，還有解決內存限制的增量PCA方法 Incremental  PCA，以及解決稀疏數據降維的PCA方法Sparse  PCA等。

　　PCA算法的主要優點：

• 1）僅僅需要以方差衡量信息量，不受數據集以外的因素影響
• 2）各主成分之間正交，可消除原始數據成分間的互相影響的因素
• 3）計算方法簡單，主要運算是特征值分解，易於實現

　　PCA算法的主要缺點：

• 1）主成分各個特征維度的含義具有一定的模糊性，不如原始樣本特征的解釋性強
• 2）方差小的非主成分也可能含有對樣本差異的重要信息，因降維丟棄可能對后續數據處理有影響

七，實例PCA計算過程

　　現在假設有一組數據如下：

　　行代表了樣例，列代表了特征，這里有10個樣例，每個樣例兩個特征，可以這樣認為，有10篇文檔，x是10篇文檔中“learn”出現的TF-IDF，y是10篇文檔中“study”出現的IF_IDF。

　　第一步，分別求x和y的平均值，然后對所有的樣例，都減去對應的均值。這里x的均值為1.81，y的均值為1.91，那么第一個一個樣例減去均值后為（0.69,0.49），以此類推得到：

 

　　第二步，計算特征協方差矩陣，如果數據是三維的，那么協方差矩陣為：

　　這里是2維的，只有x和y，求解得：

　　

　　對角線上分別是x和y的方差，非對角線上是協方差。協方差是衡量兩個變量同時變化的變化程度。協方差大於0表示x和y若一個增加，另一個也增加；協方差小於0寶石一個增加，則另一個減少。如果x和y是統計獨立的，那么二者之間的協方差就是0；但是協方差是0 ，並不能說明x和y是獨立的。協方差絕對值越大，兩者對彼此的影響越大，反之越小。協方差是沒有單位的量，因此，如果同樣的兩個變量所采用的量綱發生變化，他們的協方差也會產生數值上的變化。

 

　　第三步，計算協方差矩陣的特征向量和特征值，選取特征向量

 

　　 　　上面兩個特征值，下面是對應的特征向量，特征值0.490833989對應的特征向量是（-0.735178656, 0.677873399），這里的特征向量是正交的、歸一化的，即長度為1。

 

　　第四步，將特征值按照從大到小的順序排序，選擇其中最大的k個，然后將其對應的k各特征向量分別作為列向量組成特征向量矩陣。

 

　　如果數據中有n維，計算出n個特征向量和特征值，選擇前k個特征向量，然后最終的數據集合只有k維，取的特征向量命名為FeatureVector。 

  　  

　　這里的特征值只有兩個，我們選擇其中最大的那個，這里是1.28402771，對應的特征向量是(-0.677873399, -0.735178656)T。 

　　第五步，將樣本點投影到選取的特征向量上，得到新的數據集

　　假設樣例數為m，特征數為n，減去均值后的樣本矩陣為DataAdjust(m*n)，協方差矩陣是n*n，選取的k個特征向量組成的矩陣為EigenVectors(n*k)。那么投影后的數據FinalData為

　　得到結果為

　　

　　這樣，就將原始樣例的n維特征變成了k維，這k維就是原始特征在k維上的投影。

　　代碼如下：

from  numpy  import *

li = [[2.5,2.4],[0.5,0.7],[2.2,2.9],[1.9,2.2],[3.1,3.0],[2.3,2.7],[2.0,1.6],[1.0,1.1],[1.5,1.6],[1.1,0.9]]
matrix = mat(li)
# 求均值
mean_matrix = mean(matrix,axis=0)
# print(mean_matrix)    #[[1.81 1.91]]
# 減去平均值
Dataadjust = matrix - mean_matrix
# print(Dataadjust)
'''
[[ 0.69  0.49]
 [-1.31 -1.21]
 [ 0.39  0.99]
 [ 0.09  0.29]
 [ 1.29  1.09]
 [ 0.49  0.79]
 [ 0.19 -0.31]
 [-0.81 -0.81]
 [-0.31 -0.31]
 [-0.71 -1.01]]
'''
#計算特征值和特征向量
covMatrix = cov(Dataadjust,rowvar=0)
# print(covMatrix)
'''
[[0.61655556 0.61544444]
 [0.61544444 0.71655556]]
'''
eigValues , eigVectors = linalg.eig(covMatrix)
# print(eigValues)
# print(eigVectors)
'''
[0.0490834  1.28402771]

[[-0.73517866 -0.6778734 ]
 [ 0.6778734  -0.73517866]]'''
# 對特征值進行排序
eigValuesIndex = argsort(eigValues)
# print(eigValuesIndex)
# 保留前K個最大的特征值
eigValuesIndex = eigValuesIndex[:-(1000000):-1]
# print(eigValuesIndex)
# 計算出對應的特征向量
trueEigVectors = eigVectors[:,eigValuesIndex]
# print(trueEigVectors)
'''
[[-0.6778734  -0.73517866]
 [-0.73517866  0.6778734 ]]
 '''
# # 選擇較大特征值對應的特征向量
maxvector_eigval = trueEigVectors[：，0]
# print(maxvector_eigval)
'''
[-0.6778734  -0.73517866]
'''
# # 執行PCA變換：Y=PX 得到的Y就是PCA降維后的值 數據集矩陣
pca_result = maxvector_eigval * Dataadjust.T
# print(pca_result)
'''
[[-0.82797019  1.77758033 -0.99219749 -0.27421042 -1.67580142 -0.9129491
   0.09910944  1.14457216  0.43804614  1.22382056]]
'''

　　

　　上面的數據可以認為是learn和study特征融合為一個新的特征叫LS特征，該特征基本上代表了這兩個特征，該過程如下圖所示：

　　正號表示預處理后的樣本點，斜着的兩條線就分別是正交的特征向量（由於協方差矩陣是對稱的，因此其特征向量正交），最后一句矩陣乘法就是將原始樣本點分別往特征向量對應的軸上做投影，下圖是FinalData根據最大特征值對應的特征向量轉化回去后的數據集形式，可看出是將DataAdjust樣本點分別往特征向量對應的軸上做投影：

　　

　　如果取的k=2，那么結果是

     

　　可見，若使用了所有特征向量得到的新的數據集，轉化回去之后，與原來的數據集完全一樣（只是坐標軸旋轉）。

八，python實現主成分（PCA）降維

from numpy import *

def loadDataSet(fileName, delim='\t'):
    fr = open(fileName)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    datArr = [map(float,line) for line in stringArr]
    return mat(datArr)

def pca(dataMat, topNfeat=999999):
    meanVals = mean(dataMat, axis=0)
    DataAdjust = dataMat - meanVals           #減去平均值
    covMat = cov(DataAdjust, rowvar=0)
    eigVals,eigVects = linalg.eig(mat(covMat)) #計算特征值和特征向量
    #print eigVals
    eigValInd = argsort(eigVals)
    eigValInd = eigValInd[:-(topNfeat+1):-1]   #保留最大的前K個特征值
    redEigVects = eigVects[:,eigValInd]        #對應的特征向量
    lowDDataMat = DataAdjust * redEigVects     #將數據轉換到低維新空間
    reconMat = (lowDDataMat * redEigVects.T) + meanVals   #重構數據，用於調試
    return lowDDataMat, reconMat

　　測試數據testSet.txt由1000個數據點組成。下面對數據進行降維，並用matplotlib模塊將降維后的數據和原始數據一起繪制出來。

（鏈接：https://pan.baidu.com/s/19k4ND3ISUjhfWZhj4Fhcbw 提取碼：r3n7 ）

數據：（此數據直接復制可能無法使用，會報錯， could not convert string to float，建議最好下載，或者去我的GitHub拿：https://github.com/LeBron-Jian/MachineLearningNote）

import matplotlib
import matplotlib.pyplot as plt

dataMat = loadDataSet('testSet.txt')
lowDMat, reconMat = pca(dataMat,1)
print "shape(lowDMat): ",shape(lowDMat)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(dataMat[:,0].flatten().A[0],dataMat[:,1].flatten().A[0],marker='^',s=90)
ax.scatter(reconMat[:,0].flatten().A[0],reconMat[:,1].flatten().A[0],marker='o',s=50,c='red')
plt.show()

　　

結果如下圖：

 九，sklearn實現PCA降維

　　　　在scikit-learn中，與PCA相關的類都在sklearn.decomposition包中。最常用的PCA類就是sklearn.decomposition.PCA，我們下面主要也會講解基於這個類的使用的方法。

　　　　除了PCA類以外，最常用的PCA相關類還有KernelPCA類，在上面我們也講到了，它主要用於非線性數據的降維，需要用到核技巧。因此在使用的時候需要選擇合適的核函數並對核函數的參數進行調參。

　　　　另外一個常用的PCA相關類是IncrementalPCA類，它主要是為了解決單機內存限制的。有時候我們的樣本量可能是上百萬+，維度可能也是上千，直接去擬合數據可能會讓內存爆掉， 此時我們可以用IncrementalPCA類來解決這個問題。IncrementalPCA先將數據分成多個batch，然后對每個batch依次遞增調用partial_fit函數，這樣一步步的得到最終的樣本最優降維。

　　　　此外還有SparsePCA和MiniBatchSparsePCA。他們和上面講到的PCA類的區別主要是使用了L1的正則化，這樣可以將很多非主要成分的影響度降為0，這樣在PCA降維的時候我們僅僅需要對那些相對比較主要的成分進行PCA降維，避免了一些噪聲之類的因素對我們PCA降維的影響。SparsePCA和MiniBatchSparsePCA之間的區別則是MiniBatchSparsePCA通過使用一部分樣本特征和給定的迭代次數來進行PCA降維，以解決在大樣本時特征分解過慢的問題，當然，代價就是PCA降維的精確度可能會降低。使用SparsePCA和MiniBatchSparsePCA需要對L1正則化參數進行調參。

9.1 sklearn.decomposition.PCA參數介紹

　　　　下面我們主要基於sklearn.decomposition.PCA來講解如何使用scikit-learn進行PCA降維。PCA類基本不需要調參，一般來說，我們只需要指定我們需要降維到的維度，或者我們希望降維后的主成分的方差和占原始維度所有特征方差和的比例閾值就可以了。

　　　　現在我們對sklearn.decomposition.PCA的主要參數做一個介紹：

　　　　1）n_components：這個參數可以幫我們指定希望PCA降維后的特征維度數目。最常用的做法是直接指定降維到的維度數目，此時n_components是一個大於等於1的整數。當然，我們也可以指定主成分的方差和所占的最小比例閾值，讓PCA類自己去根據樣本特征方差來決定降維到的維度數，此時n_components是一個（0，1]之間的數。當然，我們還可以將參數設置為"mle", 此時PCA類會用MLE算法根據特征的方差分布情況自己去選擇一定數量的主成分特征來降維。我們也可以用默認值，即不輸入n_components，此時n_components=min(樣本數，特征數)。

　　　　2）whiten ：判斷是否進行白化。所謂白化，就是對降維后的數據的每個特征進行歸一化，讓方差都為1.對於PCA降維本身來說，一般不需要白化。如果你PCA降維后有后續的數據處理動作，可以考慮白化。默認值是False，即不進行白化。

　　　　3）svd_solver：即指定奇異值分解SVD的方法，由於特征分解是奇異值分解SVD的一個特例，一般的PCA庫都是基於SVD實現的。有4個可以選擇的值：{‘auto’, ‘full’, ‘arpack’, ‘randomized’}。randomized一般適用於數據量大，數據維度多同時主成分數目比例又較低的PCA降維，它使用了一些加快SVD的隨機算法。 full則是傳統意義上的SVD，使用了scipy庫對應的實現。arpack和randomized的適用場景類似，區別是randomized使用的是scikit-learn自己的SVD實現，而arpack直接使用了scipy庫的sparse SVD實現。默認是auto，即PCA類會自己去在前面講到的三種算法里面去權衡，選擇一個合適的SVD算法來降維。一般來說，使用默認值就夠了。

　　　　除了這些輸入參數外，有兩個PCA類的成員值得關注。第一個是explained_variance_，它代表降維后的各主成分的方差值。方差值越大，則說明越是重要的主成分。第二個是explained_variance_ratio_，它代表降維后的各主成分的方差值占總方差值的比例，這個比例越大，則越是重要的主成分。

9.2 PCA實例

　　下面學習一個示例。

　　首先我們生成隨機數據並可視化，代碼如下：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn.datasets.samples_generator 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='o')

 　　生成的圖如下：

 　　我們先不降維，只對數據進行投影，看看投影后的三個維度的方差分析，代碼如下：

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]

 　　可以看出投影后三個特征維度的方差比例大約為 0.983,  0.08， 0.08。投影后第一個特征占了絕大多數的主成分比例。

　　現在我們開始進行降維，從三維降到二維，代碼只需要換參數即可，我們看結果：

[0.98318212 0.00850037]
[3.78521638 0.03272613]

 　　這個結果其實是可以語料的，因為上面三個投影后的特征維度的方差分別為 [3.78521638 0.03272613 0.03202212] ，投影到二維后選擇的肯定是前兩個特征，而拋棄第三個特征。

　　為了有個直觀的認識，我們看看此時轉換后的數據分布，代碼如下：

X_new = pca.transform(X)
plt.scatter(X_new[:, 0], X_new[:, 1], marker='o')
plt.show()

 　　輸出圖如下：

 　　可見降維后的數據依然可以很清楚的看到我們之前三維圖中的四個簇。

　　現在我們看看不直接指定降維的維度，而指定降維后的主成分方差和比例。

pca = PCA(n_components=0.95)
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
print(pca.n_components_)
# [0.98318212]
# [3.78521638]
# 1

 　　我們指定了主成分至少占 95% ，輸出如上，可見只有第一個投影特征被保留，其實也很好理解，我們的第一個主成分占投影特征的方差比例高達 98 %，只選擇這一個特征便可以滿足 95%的閾值。我們現在選擇閾值 99%看看：

[0.98318212 0.00850037]
[3.78521638 0.03272613]
2

 　　其實這個也很好理解，因為我們第一個主成分占了 98.3% 的方差比例，第二個主成分占了0.08的方差比例，兩者一起就可以滿足我們的閾值。

　　最后我們看看讓 MLE算法自己選擇降維維度的效果，代碼如下：

pca = PCA(n_components='mle')
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
print(pca.n_components_)
# [0.98318212]
# [3.78521638]
# 1

 　　可見，由於我們數據的第一個投影特征的方差占比高達 98.3%，MLE算法只保留了我們第一個特征。

 

 

 

完整代碼及其數據，請移步小編的GitHub

　　傳送門：請點擊我

　　如果點擊有誤：https://github.com/LeBron-Jian/MachineLearningNote

 參考文獻：

http://www.360doc.com/content/13/1124/02/9482_331688889.shtml

 http://www.docin.com/p-160874280.html

https://www.zhihu.com/question/54100505/answer/432025686

https://www.cnblogs.com/pinard/p/6243025.html