PCA和協方差的理解

學過概率統計的孩子都知道，統計裡最基本的概念就是樣本的均值，方差，或者再加個標準差。首先我們給你一個含有n個樣本的集合，依次給出這些概念的公式描述，這些高中學過數學的孩子都應該知道吧，一帶而過。

很顯然，均值描述的是樣本集合的中間點，它告訴我們的資訊是很有限的，而標準差給我們描述的則是樣本集合的各個樣本點到均值的距離之平均。以這兩個集合為例，[0，8，12，20]和[8，9，11，12]，兩個集合的均值都是10，但顯然兩個集合差別是很大的，計算兩者的標準差，前者是8.3，後者是1.8，顯然後者較為集中，故其標準差小一些，標準差描述的就是這種“散佈度”。之所以除以n-1而不是除以n，是因為這樣能使我們以較小的樣本集更好的逼近總體的標準差，即統計上所謂的“無偏估計”。而方差則僅僅是標準差的平方。

為什麼需要協方差？

上面幾個統計量看似已經描述的差不多了，但我們應該注意到，標準差和方差一般是用來描述一維資料的，但現實生活我們常常遇到含有多維資料的資料集，最簡單的大家上學時免不了要統計多個學科的考試成績。面對這樣的資料集，我們當然可以按照每一維獨立的計算其方差，但是通常我們還想了解更多，比如，一個男孩子的猥瑣程度跟他受女孩子歡迎程度是否存在一些聯絡啊，嘿嘿~協方差就是這樣一種用來度量兩個隨機變數關係的統計量，我們可以仿照方差的定義：

來度量各個維度偏離其均值的程度，標準差可以這麼來定義：

協方差的結果有什麼意義呢？如果結果為正值，則說明兩者是正相關的(從協方差可以引出“相關係數”的定義)，也就是說一個人越猥瑣就越受女孩子歡迎，嘿嘿，那必須的~結果為負值就說明負相關的，越猥瑣女孩子越討厭，可能嗎？如果為0，也是就是統計上說的“相互獨立”。

從協方差的定義上我們也可以看出一些顯而易見的性質，如：

協方差多了就是協方差矩陣

上一節提到的猥瑣和受歡迎的問題是典型二維問題，而協方差也只能處理二維問題，那維數多了自然就需要計算多個協方差，比如n維的資料集就需要計算 n! / ((n-2)!*2) 個協方差，那自然而然的我們會想到使用矩陣來組織這些資料。給出協方差矩陣的定義：

這個定義還是很容易理解的，我們可以舉一個簡單的三維的例子，假設資料集有三個維度，則協方差矩陣為

可見，協方差矩陣是一個對稱的矩陣，而且對角線是各個維度上的方差。

Matlab協方差實戰

上面涉及的內容都比較容易，協方差矩陣似乎也很簡單，但實戰起來就很容易讓人迷茫了。必須要明確一點，協方差矩陣計算的是不同維度之間的協方差，而不是不同樣本之間的。

這個我將結合下面的例子說明，以下的演示將使用Matlab，為了說明計算原理，不直接呼叫Matlab的cov函式(藍色部分為Matlab程式碼)。

首先，隨機產生一個10*3維的整數矩陣作為樣本集，10為樣本的個數，3為樣本的維數。

mysample = fix(rand(10,3)*50)

根據公式，計算協方差需要計算均值，那是按行計算均值還是按列呢，我一開始就老是困擾這個問題。前面我們也特別強調了，協方差矩陣是計算不同維度間的協方差，要時刻牢記這一點。樣本矩陣的每行是一個樣本，每列為一個維度，所以我們要按列計算均值。為了描述方便，我們先將三個維度的資料分別賦值：

>> dim1 = mysample(:,1);
>> dim2 = mysample(:,2);
>> dim3 = mysample(:,3);

計算dim1與dim2，dim1與dim3，dim2與dim3的協方差：

>> sum((dim1 - mean(dim1)) .* (dim2 - mean(dim2))) / (size(mysample, 1) - 1)  %得到 -147.0667
>> sum((dim1 - mean(dim1)) .* (dim3 - mean(dim3))) / (size(mysample, 1) - 1)  %得到  -82.2667
>> sum((dim2 - mean(dim2)) .* (dim3 - mean(dim3))) / (size(mysample, 1) - 1)  %得到   76.5111

搞清楚了這個後面就容易多了，協方差矩陣的對角線就是各個維度上的方差，下面我們依次計算：

>> var(dim1)  %得到 227.8778
>> var(dim2)  %得到 179.8222
>> var(dim3)  %得到 156.7111

這樣，我們就得到了計算協方差矩陣所需要的所有資料，呼叫Matlab自帶的cov函式進行驗證：

>> cov(mysample)

把我們計算的資料對號入座，是不是一摸一樣？

Update

今天突然發現，原來協方差矩陣還可以這樣計算，先讓樣本矩陣中心化，即每一維度減去該維度的均值，使每一維度上的均值為0，然後直接用新的到的樣本矩陣乘上它的轉置，然後除以(N-1)即可。其實這種方法也是由前面的公式推導而來，只不過理解起來不是很直觀，但在抽象的公式推導時還是很常用的！同樣給出Matlab程式碼實現：

>> temp = mysample - repmat(mean(mysample), 10, 1);
>> result = temp' * temp ./ (size(mysample, 1) - 1)

總結

理解協方差矩陣的關鍵就在於牢記它計算的是不同維度之間的協方差，而不是不同樣本之間，拿到一個樣本矩陣，我們最先要明確的就是一行是一個樣本還是一個維度，心中明確這個整個計算過程就會順流而下，這麼一來就不會迷茫了~

再來看看PCA:

問題：假設在IR中我們建立的文件-詞項矩陣中，有兩個詞項為“learn”和“study”，在傳統的向量空間模型中，認為兩者獨立。然而從語義的角度來講，兩者是相似的，而且兩者出現頻率也類似，是不是可以合成為一個特徵呢？

《模型選擇和規則化》談到的特徵選擇的問題，就是要剔除的特徵主要是和類標籤無關的特徵。比如“學生的名字”就和他的“成績”無關，使用的是互資訊的方法。

而這裡的特徵很多是和類標籤有關的，但裡面存在噪聲或者冗餘。在這種情況下，需要一種特徵降維的方法來減少特徵數，減少噪音和冗餘，減少過度擬合的可能性。

PCA的思想是將n維特徵對映到k維上（k<n），這k維是全新的正交特徵。這k維特徵稱為主元，是重新構造出來的k維特徵，而不是簡單地從n維特徵中去除其餘n-k維特徵。

ＰＣＡ計算過程：

　　假設我們得到的2維資料如下：

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

　　第一步分別求x和y的平均值，然後對於所有的樣例，都減去對應的均值。這裡x的均值是1.81，y的均值是1.91，那麼一個樣例減去均值後即為（0.69,0.49），得到

第二步，求特徵協方差矩陣，如果資料是3維，那麼協方差矩陣是

這裡只有x和y，求解得

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

　　第三步，求協方差的特徵值和特徵向量，得到

上面是兩個特徵值，下面是對應的特徵向量，特徵值0.0490833989對應特徵向量為，這裡的特徵向量都歸一化為單位向量。

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

這裡特徵值只有兩個，我們選擇其中最大的那個，這裡是1.28402771，對應的特徵向量是。

第五步，將樣本點投影到選取的特徵向量上。假設樣例數為m，特徵數為n，減去均值後的樣本矩陣為DataAdjust(m*n)，協方差矩陣是n*n，選取的k個特徵向量組成的矩陣為EigenVectors(n*k)。那麼投影后的資料FinalData為

這裡是

FinalData(10*1) = DataAdjust(10*2矩陣)×特徵向量

得到結果是

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

上面的資料可以認為是learn和study特徵融合為一個新的特徵叫做LS特徵，該特徵基本上代表了這兩個特徵。

　　上述過程有個圖描述：

正號表示預處理後的樣本點，斜著的兩條線就分別是正交的特徵向量（由於協方差矩陣是對稱的，因此其特徵向量正交），最後一步的矩陣乘法就是將原始樣本點分別往特徵向量對應的軸上做投影。

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

這就是經過PCA處理後的樣本資料，水平軸（上面舉例為LS特徵）基本上可以代表全部樣本點。整個過程看起來就像將座標系做了旋轉，當然二維可以圖形化表示，高維就不行了。上面的如果k=1，那麼只會留下這裡的水平軸，軸上是所有點在該軸的投影。

這樣PCA的過程基本結束。在第一步減均值之後，其實應該還有一步對特徵做方差歸一化。比如一個特徵是汽車速度（0到100），一個是汽車的座位數（2到6），顯然第二個的方差比第一個小。因此，如果樣本特徵中存在這種情況，那麼在第一步之後，求每個特徵的標準差，然後對每個樣例在該特徵下的資料除以。

歸納一下，使用我們之前熟悉的表示方法，在求協方差之前的步驟是：

其中是樣例，共m個，每個樣例n個特徵，也就是說是n維向量。是第i個樣例的第j個特徵。是樣例均值。是第j個特徵的標準差。（個人認為這兒有問題，第一步應該是每個維度的均值，第二步應該是每個維度減去其對應的均值，至於三四步是縮放同一個範圍，對結果沒有影響。）

整個PCA過程貌似及其簡單，就是求協方差的特徵值和特徵向量，然後做資料轉換。但是有沒有覺得很神奇，為什麼求協方差的特徵向量就是最理想的k維向量？其背後隱藏的意義是什麼？整個PCA的意義是什麼？

PCA理論基礎

要解釋為什麼協方差矩陣的特徵向量就是k維理想特徵，我看到的有三個理論：分別是最大方差理論、最小錯誤理論和座標軸相關度理論。這裡簡單探討前兩種，最後一種在討論PCA意義時簡單概述。

最大方差理論

在訊號處理中認為訊號具有較大的方差，噪聲有較小的方差，信噪比就是訊號與噪聲的方差比，越大越好。如前面的圖，樣本在橫軸上的投影方差較大，在縱軸上的投影方差較小，那麼認為縱軸上的投影是由噪聲引起的。

因此我們認為，最好的k維特徵是將n維樣本點轉換為k維後，每一維上的樣本方差都很大。

比如下圖有5個樣本點：（已經做過預處理，均值為0，特徵方差歸一）

下面將樣本投影到某一維上，這裡用一條過原點的直線表示（前處理的過程實質是將原點移到樣本點的中心點）。

假設我們選擇兩條不同的直線做投影，那麼左右兩條中哪個好呢？根據我們之前的方差最大化理論，左邊的好，因為投影后的樣本點之間方差最大。

這裡先解釋一下投影的概念：

紅色點表示樣例，藍色點表示在u上的投影，u是直線的斜率也是直線的方向向量，而且是單位向量。藍色點是在u上的投影點，離原點的距離是（即或者）由於這些樣本點（樣例）的每一維特徵均值都為0，因此投影到u上的樣本點（只有一個到原點的距離值）的均值仍然是0。

回到上面左右圖中的左圖，我們要求的是最佳的u，使得投影后的樣本點方差最大。

由於投影后均值為0，因此方差為：

中間那部分很熟悉啊，不就是樣本特徵的協方差矩陣麼（的均值為0，一般協方差矩陣都除以m-1，這裡用m）。

用來表示，表示，那麼上式寫作

由於u是單位向量，即，上式兩邊都左乘u得，

即

We got it！就是的特徵值，u是特徵向量。最佳的投影直線是特徵值最大時對應的特徵向量，其次是第二大對應的特徵向量，依次類推。

因此，我們只需要對協方差矩陣進行特徵值分解，得到的前k大特徵值對應的特徵向量就是最佳的k維新特徵，而且這k維新特徵是正交的。得到前k個u以後，樣例通過以下變換可以得到新的樣本。

其中的第j維就是在上的投影。

通過選取最大的k個u，使得方差較小的特徵（如噪聲）被丟棄。

最小平方誤差理論：

假設有這樣的二維樣本點（紅色點），回顧我們前面探討的是求一條直線，使得樣本點投影到直線上的點的方差最大。本質是求直線，那麼度量直線求的好不好，不僅僅只有方差最大化的方法。再回想我們最開始學習的線性迴歸等，目的也是求一個線性函式使得直線能夠最佳擬合樣本點，那麼我們能不能認為最佳的直線就是迴歸後的直線呢？迴歸時我們的最小二乘法度量的是樣本點到直線的座標軸距離。比如這個問題中，特徵是x，類標籤是y。迴歸時最小二乘法度量的是距離d。如果使用迴歸方法來度量最佳直線，那麼就是直接在原始樣本上做迴歸了，跟特徵選擇就沒什麼關係了。

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

現在有n個樣本點，每個樣本點為m維（這節內容中使用的符號與上面的不太一致，需要重新理解符號的意義）。將樣本點在直線上的投影記為，那麼我們就是要最小化