7.1-PCA

概述

  PCA有兩種通俗易懂的解:
1. 是最大化投影后數據的方差(讓數據更分散)；
2. 是最小化投影造成的損失。

  這兩個思路最後都能推導出同樣的結果。
  下圖應該是對PCA第二種解釋展示得最好的一張圖片了.

原理

  在理解特徵提取與處理時，涉及高維特徵向量的問題往往容易陷入維度災難。隨着數據集維度的增加，算法學習需要的樣本數量呈指數級增加。有些應用中，遇到這樣的大數據是非常不利的，而且從大數據集中學習需要更多的內存和處理能力。另外，隨着維度的增加，數據的稀疏性會越來越高。在高維向量空間中探索同樣的數據集比在同樣稀疏的數據集中探索更加困難。
  PCA將數據投射到一個低維子空間實現降維。例如，二維數據集降維就是把點投射成一條線，數據集的每個樣本都可以用一個值表示，不需要兩個值。三維數據集可以降成二維，就是把變量映射成一個平面。一般情況下，n維數據集可以通過映射降成k維子空間，其中k≤n。
假如你是一本養花工具宣傳冊的攝影師，你正在拍攝一個水壺。水壺是三維的，但是照片是二維的，爲了更全面的把水壺展示給客戶，你需要從不同角度拍幾張圖片。下圖是你從四個方向拍的照片：

  第一張圖裏水壺的背面可以看到，但是看不到前面。第二張圖是拍前面，可以看到壺嘴，這張圖可以提供了第一張圖缺失的信息，但是壺把看不到了。從第三張俯視圖裏無法看出壺的高度。第四張圖是你真正想要的，水壺的高度，頂部，壺嘴和壺把都清晰可見。
PCA的設計理念與此類似，它可以將高維數據集映射到低維空間的同時，儘可能的保留更多信息。PCA旋轉數據集與其主成分對齊，將最多的信息保留到第一主成分中。假設我們有下圖所示的數據集：

  數據集看起來像一個從原點到右上角延伸的細長扁平的橢圓。要降低整個數據集的維度，我們必須把點映射成一條線。下圖中的兩條線都是數據集可以映射的，映射到哪條線樣本變化最大？

  顯然，樣本映射到黑色虛線的變化比映射到紅色點線的變化要大的多。實際上，這條黑色虛線就是第一主成分。第二主成分必須與第一主成分正交，也就是說第二主成分必須是在統計學上獨立的，會出現在與第一主成分垂直的方向，如下圖所示：

  後面的每個主成分也會盡量多的保留剩下的信息，唯一的要求就是每一個主成分需要和前面的主成分正交。
現在假設數據集是三維的，散點圖看起來像是沿着一個軸旋轉的圓盤。

  這些點可以通過旋轉和變換使圓盤完全變成二維的。現在這些點看着像一個橢圓，第三維上基本沒有信息，可以被忽略。
當數據集不同維度上的方差分佈不均勻的時候，PCA最有用。（如果是一個球殼行數據集，PCA不能有效的發揮作用，因爲各個方向上的方差都相等；沒有丟失大量的信息維度一個都不能忽略）。

理論依據

  首先看一個例子：這裏給出5個同學的成績

header 1	學生1	學生2	學生3	學生4	學生5
語文成績	50	60	70	80	90
數學成績	80	82	84	86	88

  此時如果將二維特徵【語文成績，數學成績】壓縮成1維，僅僅用數學成績表達似乎不是最佳方案，即使沿着投影方向2，只用語文成績表達，其分散程度似乎也比不上投影方向1的降維結果。
  如何找到最佳的（方差最大，即最分散。這也容易理解，畢竟分散的數據更容易區分開來）降維方式呢？顯然最佳方案不總是用語文/數學成績單一替代，這便需要本文的主角——PCA算法。簡而言之：PCA算法其表現形式是降維，同時也是一種特徵融合算法。
  推導之前先介紹幾個操作的幾何意義：

1. 內積和投影

  下面先來看一個高中就學過的向量運算：內積。兩個維數相同的向量的內積被定義爲：

(a1,a2,⋯,an)T⋅(b1,b2,⋯,bn)T=a1b1+a2b2+⋯+anbn $(a_1,a_2,\cdots ,a_n{)}^T\cdot (b_1,b_2,\cdots ,b_n{)}^T=a_1{b}_1+a_2{b}_2+\cdots +a_n{b}_n$

  內積運算將兩個向量映射爲一個實數。其計算方式非常容易理解，但是其意義並不明顯。下面我們分析內積的幾何意義。假設A和B是兩個n維向量，我們知道n維向量可以等價表示爲n維空間中的一條從原點發射的有向線段，爲了簡單起見我們假設A和B均爲二維向量，則 A=(x1,y1),B=(x2,y2) $A=(x_1,y_1),B=(x_2,y_2)$。則在二維平面上A和B可以用兩條發自原點的有向線段表示，見下圖：

  現在我們從A點向B所在直線引一條垂線。我們知道垂線與B的交點叫做A在B上的投影，再設A與B的夾角是a，則投影的矢量長度爲 |A|cos(a) $|A|cos(a)$，其中 |A|=x21+y21‾‾‾‾‾‾‾√ $|A|=\sqrt{x_1^2+y_1^2}$是向量A的模，也就是A線段的標量長度。

注意這裏我們專門區分了矢量長度和標量長度，標量長度總是大於等於0，值就是線段的長度；而矢量長度可能爲負，其絕對值是線段長度，而符號取決於其方向與標準方向相同或相反。
  到這裏還是看不出內積和這東西有什麼關係，不過如果我們將內積表示爲另一種我們熟悉的形式：

A⋅B=|A||B|cos(a) $A\cdot B=|A||B|cos(a)$

  現在事情似乎是有點眉目了：A與B的內積等於A到B的投影長度乘以B的模。再進一步，如果我們假設B的模爲1，即讓 |B|=1 $|B|=1$，那麼就變成了：

A⋅B=|A|cos(a) $A\cdot B=|A|cos(a)$

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

基的概念

給出下面這個向量：

  在代數表示方面，我們經常用線段終點的點座標表示向量，例如上面的向量可以表示爲(3,2)，這是我們再熟悉不過的向量表示。
  不過我們常常忽略，只有一個(3,2)本身是不能夠精確表示一個向量的。我們仔細看一下，這裏的3實際表示的是向量在x軸上的投影值是3，在y軸上的投影值是2。也就是說我們其實隱式引入了一個定義：以x軸和y軸上正方向長度爲1的向量爲標準。那麼一個向量(3,2)實際是說在x軸投影爲3而y軸的投影爲2。注意投影是一個矢量，所以可以爲負。更正式的說，向量(x,y)實際上表示線性組合：

x(1,0)T+y(0,1)T $x(1,0{)}^T+y(0,1{)}^T$

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

  所以，要準確描述向量，首先要確定一組基，然後給出在基所在的各個直線上的投影值，就可以了。只不過我們經常省略第一步，而默認以(1,0)和(0,1)爲基。基就是骨架，基的個數即爲撐起給定維度空間的最小向量個數，通常設定基德模爲1。

基變換的矩陣表示

  但基不是唯一的，例如上一個例子中選(1,0)和(0,1)爲基，同樣，上面的基可以變爲 (12√,12√) $(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$和 (−12√,12√) $(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$.
  現在我們想獲得(3,2)在新基上的座標，即在兩個方向上的投影矢量值，就需要利用到基變換的矩陣表示。
  根據下圖，我們只要分別計算(3,2)和兩個基的內積，不難得到新的座標爲 (52√,−12√) $(\frac{5}{\sqrt{2}},-\frac{1}{\sqrt{2}})$.

  同樣我們可以利用基變換的矩陣表示進行求解：

  再舉一例：(1,1)，(2,2)，(3,3)，想變換到剛纔那組基上，則可以這樣表示：

  於是一組向量的基變換被幹淨的表示爲矩陣的相乘。
  推廣成一般性結論：
  一般的，如果我們有M個N維向量，想將其變換爲由R個N維向量表示的新空間中，那麼首先將R個基按行組成矩陣A，然後將向量按列組成矩陣B，那麼兩矩陣的乘積AB就是變換結果，其中AB的第m列爲A中第m列變換後的結果。
  對應數學表達：

  其中 pi $p_i$是一個行向量，表示第i個基， aj $a_j$是一個列向量，表示第j個原始數據記錄。
  上述分析同時給矩陣相乘找到了一種物理解釋：兩個矩陣相乘的意義是將右邊矩陣中的每一列列向量變換到左邊矩陣中每一行行向量爲基所表示的空間中去。這種操作便是降維。

最優基的選取方法

方差

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

Var(a)=1m∑mi=1(ai−μ)2 $Var(a)=\frac{1}{m}\underset{i=1}{\overset{m}{\sum }}(a_i-\mu {)}^2$

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

Var(a)=1m∑mi=1a2i $Var(a)=\frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{a}_i^2$

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

協方差

  對於上面二維降成一維的問題來說，找到那個使得方差最大的方向就可以了。不過對於更高維，還有一個問題需要解決。考慮三維降到二維問題。與之前相同，首先我們希望找到一個方向使得投影后方差最大，這樣就完成了第一個方向的選擇，繼而我們選擇第二個投影方向。
  如果我們還是單純只選擇方差最大的方向，很明顯，這個方向與第一個方向應該是「幾乎重合在一起」，顯然這樣的維度是沒有用的，因此，應該有其他約束條件。從直觀上說，讓兩個字段儘可能表示更多的原始信息，我們是不希望它們之間存在（線性）相關性的，因爲相關性意味着兩個字段不是完全獨立，必然存在重複表示的信息。
  數學上可以用兩個字段的協方差表示其相關性，由於已經讓每個字段均值爲0，則：

Cov(a,b)=1m∑mi=1aibi $Cov(a,b)=\frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{a}_i{b}_i$

  可以看到，在字段均值爲0的情況下，兩個字段的協方差簡潔的表示爲其內積除以元素數m。
  當協方差爲0時，表示兩個字段完全獨立。爲了讓協方差爲0，我們選擇第二個基時只能在與第一個基正交的方向上選擇。因此最終選擇的兩個方向一定是正交的。
  至此，我們得到了降維問題的優化目標：將一組N維向量降爲K維（K大於0，小於N），其目標是選擇K個單位（模爲1）正交基，使得原始數據變換到這組基上後，各字段兩兩間協方差爲0，而字段的方差則儘可能大（在正交的約束下，取最大的K個方差）。

協方差矩陣

  上面我們導出了優化目標，但是這個目標似乎不能直接作爲操作指南（或者說算法），因爲它只說要什麼，但根本沒有說怎麼做。所以我們要繼續在數學上研究計算方案。
  我們看到，最終要達到的目的與字段內方差及字段間協方差有密切關係。因此我們希望能將兩者統一表示，仔細觀察發現，兩者均可以表示爲內積的形式，而內積又與矩陣相乘密切相關。於是我們來了靈感：
  假設我們只有a和b兩個字段，那麼我們將它們按行組成矩陣X：

X=(a1b1a2b2⋯⋯ambm) $X=\left(\begin{array}{cccc}{a}_1& a_2& \cdots & a_m\\ b_1& b_2& \cdots & b_m\end{array}\right)$

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

1mXXT=(1m∑mi=1a2i1m∑mi=1aibi1m∑mi=1aibi1m∑mi=1b2i) $\frac{1}{m}X{X}^T=\left(\begin{array}{cc}\frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{a}_i^2& \frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{a}_i{b}_i\\ \frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{a}_i{b}_i& \frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{b}_i^2\end{array}\right)$

  奇蹟出現了！這個矩陣對角線上的兩個元素分別是兩個字段的方差，而其它元素是a和b的協方差。兩者被統一到了一個矩陣的。
  根據矩陣相乘的運算法則，這個結論很容易被推廣到一般情況：
  設我們有m個n維數據記錄，將其按列排成n乘m的矩陣X，設 C=1mXXT $C=\frac{1}{m}X{X}^T$，則C是一個對稱陣，其對角線分別爲各個字段的方差，而第i行和j列的元素和第j行i列的元素相同，表示i和j兩個字段的協方差。

協方差矩陣對角化

  根據上述推導，我們發現要達到優化目前，等價於將協方差矩陣對角化：即除對角線外的其它元素化爲0，並且在對角線上將元素按大小從上到下排列，這樣我們就達到了優化目的。這樣說可能還不是很明晰，我們進一步看下原矩陣與基變換後矩陣協方差矩陣的關係：
  設原始數據矩陣X對應的協方差矩陣爲C，而P是一組基按行組成的矩陣，設Y=PX，則Y爲X對P做基變換後的數據。設Y的協方差矩陣爲D，我們推導一下D與C的關係：

D=====1mYYT1m(PX)(PX)T1mPXXTPTP(1mXXT)PTPCPT $\begin{array}{lll}D& =& \frac{1}{m}Y{Y}^T\\ & =& \frac{1}{m}(PX)(PX{)}^T\\ & =& \frac{1}{m}PX{X}^T{P}^T\\ & =& P(\frac{1}{m}X{X}^T)P^T\\ & =& PC{P}^T\end{array}$

  現在事情很明白了！我們要找的P不是別的，而是能讓原始協方差矩陣對角化的P。換句話說，優化目標變成了尋找一個矩陣P，滿足 PCPT $PC{P}^T$是一個對角矩陣，並且對角元素按從大到小依次排列，那麼P的前K行就是要尋找的基，用P的前K行組成的矩陣乘以X就使得X從N維降到了K維並滿足上述優化條件。
  至此，我們離「發明」PCA還有僅一步之遙！
  現在所有焦點都聚焦在了協方差矩陣對角化問題上，有時，我們真應該感謝數學家的先行，因爲矩陣對角化在線性代數領域已經屬於被玩爛了的東西，所以這在數學上根本不是問題。
  由上文知道，協方差矩陣C是一個是對稱矩陣，在線性代數上，實對稱矩陣有一系列非常好的性質：
1. 實對稱矩陣不同特徵值對應的特徵向量必然正交。
2. 設特徵向量 λ $\lambda$重數爲r，則必然存在r個線性無關的特徵向量對應於 λ $\lambda$，因此可以將這r個特徵向量單位正交化.

  由上面兩條可知，一個n行n列的實對稱矩陣一定可以找到n個單位正交特徵向量，設這n個特徵向量爲 e1,e2,⋯,en $e_1,e_2,\cdots ,e_n$，我們將其按列組成矩陣：

E=(e1e2⋯en) $E=\left(\begin{array}{cccc}{e}_1& e_2& \cdots & e_n\end{array}\right)$

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

ETCE=Λ=⎛⎝⎜⎜⎜⎜λ1λ2⋱λn⎞⎠⎟⎟⎟⎟ $E^{T}CE=\mathrm{\Lambda }=\left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_n\end{array}\right)$

  其中 Λ $\mathrm{\Lambda }$爲對角矩陣，其對角元素爲各特徵向量對應的特徵值（可能有重複）。
  以上結論不再給出嚴格的數學證明，對證明感興趣的朋友可以參考線性代數書籍關於「實對稱矩陣對角化」的內容。
  到這裏，我們發現我們已經找到了需要的矩陣P：

P=ET $P=E^T$

  P是協方差矩陣的特徵向量單位化後按行排列出的矩陣，其中每一行都是C的一個特徵向量。如果設P按照Λ中特徵值的從大到小，將特徵向量從上到下排列，則用P的前K行組成的矩陣乘以原始數據矩陣X，就得到了我們需要的降維後的數據矩陣Y。
  至此我們完成了整個PCA的數學原理討論。在下面的一節，我們將給出PCA的一個實例。

算法

  總結一下PCA的算法步驟：
  設有m條n維數據。
1. 將原始數據按列組成n行m列矩陣X.
2. 將X的每一行（代表一個屬性字段）進行零均值化，即減去這一行的均值.
3. 求出協方差矩陣 C=1mXXT $C=\frac{1}{m}X{X}^T$
4. 求出協方差矩陣的特徵值及對應的特徵向量.
5. 將特徵向量按對應特徵值大小從上到下按行排列成矩陣，取前k行組成矩陣P.
6. Y=PX $Y=PX$即爲降維到k維後的數據.

實例

A-示例1：二維數據降維

  這裏以上文提到的

(−1−2−10002101) $\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right)$

  爲例，我們用PCA方法將這組二維數據其降到一維。
  因爲這個矩陣的每行已經是零均值，這裏我們直接求協方差矩陣：

C=15(−1−2−10002101)⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜−1−1020−20011⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜65454565⎞⎠⎟⎟⎟ $C=\frac{1}{5}\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right)\left(\begin{array}{cc}-1& -2\\ \\ -1& 0\\ \\ 0& 0\\ \\ 2& 1\\ \\ 0& 1\end{array}\right)=\left(\begin{array}{cc}\frac{6}{5}& \frac{4}{5}\\ \\ \frac{4}{5}& \frac{6}{5}\end{array}\right)$

  然後求其特徵值和特徵向量，具體求解方法不再詳述，可以參考相關資料。求解後特徵值爲：

λ1=2,λ2=2/5 ${\lambda }_1=2,{\lambda }_2=2/5$

  其對應的特徵向量分別是：

c1(11),c2⎛⎝⎜⎜−1111),c2⎛⎝⎜⎜−11⎞⎠⎟⎟ $c_1\left(\begin{array}{c}1\\ 1\end{array}\right),c_2\left(\begin{array}{c}-1\\ \\ 1\end{array}\right)$

  其中對應的特徵向量分別是一個通解，c1和c2可取任意實數。那麼標準化後的特徵向量爲：

(1/⎛⎝⎜⎜−11⎞⎠⎟⎟ $c_1\left(\begin{array}{c}1\\ 1\end{array}\right),c_2\left(\begin{array}{c}-1\\ \\ 1\end{array}\right)$

  其中對應的特徵向量分別是一個通解，c1和c2可取任意實數。那麼標準化後的特徵向量爲：

(1/2‾√1⎞⎠⎟⎟