協(xié)方差矩陣與主成成分分析(整理后)

1、統(tǒng)計學(xué)的基本概念統(tǒng)計里最基本的概念就是樣本的均值，方差，或者再加個標準差。首先我們給你一個含有n個樣本的集合，依次給出這些概念的公式描述，這些高中學(xué)過數(shù)學(xué)的孩子都應(yīng)該知道吧，一帶而過。均值：標準差：方差：很顯然，均值描述的是樣本集合的中間點，它告訴我們的信息是很有限的，而標準差給我們描述的則是樣本集合的各個樣本點到均值的距離之平均。以這兩個集合為例，0，8，12，20和8，9，11，12，兩個集合的均值都是10，但顯然兩個集合差別是很大的，計算兩者的標準差，前者是8.3，后者是1.8，顯然后者較為集中，故其標準差小一些，標準差描述的就是這種“散布度”。之所以除以n-1而不是除以n，是因為這樣能

2、使我們以較小的樣本集更好的逼近總體的標準差，即統(tǒng)計上所謂的“無偏估計”。而方差則僅僅是標準差的平方。為什么需要協(xié)方差？上面幾個統(tǒng)計量看似已經(jīng)描述的差不多了，但我們應(yīng)該注意到，標準差和方差一般是用來描述一維數(shù)據(jù)的，但現(xiàn)實生活我們常常遇到含有多維數(shù)據(jù)的數(shù)據(jù)集，最簡單的大家上學(xué)時免不了要統(tǒng)計多個學(xué)科的考試成績。面對這樣的數(shù)據(jù)集，我們當然可以按照每一維獨立的計算其方差，但是通常我們還想了解更多，協(xié)方差就是這樣一種用來度量兩個隨機變量關(guān)系的統(tǒng)計量，我們可以仿照方差的定義：來度量各個維度偏離其均值的程度，標準差可以這么來定義：協(xié)方差的結(jié)果有什么意義呢？如果結(jié)果為正值，則說明兩者是正相關(guān)的(從協(xié)方差可以引出

3、“相關(guān)系數(shù)”的定義)，那必須的結(jié)果為負值就說明負相關(guān)的，如果為0，也是就是統(tǒng)計上說的“相互獨立”。從協(xié)方差的定義上我們也可以看出一些顯而易見的性質(zhì)，如：協(xié)方差多了就是協(xié)方差矩陣上一節(jié)提到的猥瑣和受歡迎的問題是典型二維問題，而協(xié)方差也只能處理二維問題，那維數(shù)多了自然就需要計算多個協(xié)方差，比如n維的數(shù)據(jù)集就需要計算個協(xié)方差，那自然而然的我們會想到使用矩陣來組織這些數(shù)據(jù)。給出協(xié)方差矩陣的定義：這個定義還是很容易理解的，我們可以舉一個簡單的三維的例子，假設(shè)數(shù)據(jù)集有三個維度，則協(xié)方差矩陣為可見，協(xié)方差矩陣是一個對稱的矩陣，而且對角線是各個維度上的方差。Matlab協(xié)方差實戰(zhàn)上面涉及的內(nèi)容都比較容易，協(xié)方

4、差矩陣似乎也很簡單，但實戰(zhàn)起來就很容易讓人迷茫了。必須要明確一點，協(xié)方差矩陣計算的是不同維度之間的協(xié)方差，而不是不同樣本之間的。這個我將結(jié)合下面的例子說明，以下的演示將使用Matlab，為了說明計算原理，不直接調(diào)用Matlab的cov函數(shù)(藍色部分為Matlab代碼)。首先，隨機產(chǎn)生一個10*3維的整數(shù)矩陣作為樣本集，10為樣本的個數(shù)，3為樣本的維數(shù)。1MySample = fix(rand(10,3)*50)根據(jù)公式，計算協(xié)方差需要計算均值，那是按行計算均值還是按列呢，我一開始就老是困擾這個問題。前面我們也特別強調(diào)了，協(xié)方差矩陣是計算不同維度間的協(xié)方差，要時刻牢記這一點。樣本矩陣的每行是一個

5、樣本，每列為一個維度，所以我們要按列計算均值。為了描述方便，我們先將三個維度的數(shù)據(jù)分別賦值：123dim1 = MySample(:,1); dim2 = MySample(:,2); dim3 = MySample(:,3);計算dim1與dim2，dim1與dim3，dim2與dim3的協(xié)方差：123sum( (dim1-mean(dim1) .* (dim2-mean(dim2) ) / ( size(MySample,1)-1 ) % 得到 74.5333 sum( (dim1-mean(dim1) .* (dim3-mean(dim3) ) / ( size(MySample,1)-

6、1 ) % 得到 -10.0889 sum( (dim2-mean(dim2) .* (dim3-mean(dim3) ) / ( size(MySample,1)-1 ) % 得到 -106.4000搞清楚了這個后面就容易多了，協(xié)方差矩陣的對角線就是各個維度上的方差，下面我們依次計算：123std(dim1)2 % 得到 108.3222 std(dim2)2 % 得到 260.6222 std(dim3)2 % 得到 94.1778這樣，我們就得到了計算協(xié)方差矩陣所需要的所有數(shù)據(jù)，調(diào)用Matlab自帶的cov函數(shù)進行驗證：1cov(MySample)把我們計算的數(shù)據(jù)對號入座，是不是一摸一樣

7、？Update：今天突然發(fā)現(xiàn)，原來協(xié)方差矩陣還可以這樣計算，先讓樣本矩陣中心化，即每一維度減去該維度的均值，使每一維度上的均值為0，然后直接用新的到的樣本矩陣乘上它的轉(zhuǎn)置，然后除以(N-1)即可。其實這種方法也是由前面的公式通道而來，只不過理解起來不是很直觀，但在抽象的公式推導(dǎo)時還是很常用的！同樣給出Matlab代碼實現(xiàn)：12X = MySample - repmat(mean(MySample),10,1); % 中心化樣本矩陣，使各維度均值為0 C = (X'*X)./(size(X,1)-1);總結(jié)理解協(xié)方差矩陣的關(guān)鍵就在于牢記它計算的是不同維度之間的協(xié)方差，而不是不同樣本之間，

8、拿到一個樣本矩陣，我們最先要明確的就是一行是一個樣本還是一個維度，心中明確這個整個計算過程就會順流而下，這么一來就不會迷茫了P.S.寫論文要選Latex，在wordpress里編輯公式還得用Latex，用Latex還真對得起咱學(xué)計算機這張臉自從上次談了協(xié)方差矩陣之后，感覺寫這種科普性文章還不錯，那我就再談一把協(xié)方差矩陣吧。上次那篇文章在理論層次介紹了下協(xié)方差矩陣，沒準很多人覺得這東西用處不大，其實協(xié)方差矩陣在好多學(xué)科里都有很重要的作用，比如多維的正態(tài)分布，再比如今天我們今天的主角主成分分析(Principal Component Analysis，簡稱PCA)。結(jié)合PCA相信能對協(xié)方差矩陣有個

9、更深入的認識PCA的緣起PCA大概是198x年提出來的吧，簡單的說，它是一種通用的降維工具。在我們處理高維數(shù)據(jù)的時候，為了能降低后續(xù)計算的復(fù)雜度，在“預(yù)處理”階段通常要先對原始數(shù)據(jù)進行降維，而PCA就是干這個事的。本質(zhì)上講，PCA就是將高維的數(shù)據(jù)通過線性變換投影到低維空間上去，但這個投影可不是隨便投投，要遵循一個指導(dǎo)思想，那就是：找出最能夠代表原始數(shù)據(jù)的投影方法。這里怎么理解這個思想呢？“最能代表原始數(shù)據(jù)”希望降維后的數(shù)據(jù)不能失真，也就是說，被PCA降掉的那些維度只能是那些噪聲或是冗余的數(shù)據(jù)。這里的噪聲和冗余我認為可以這樣認識：· 噪聲：我們常說“噪音污染”，意思就是“噪聲”干擾我們

10、想聽到的真正聲音。同樣，假設(shè)樣本中某個主要的維度A，它能代表原始數(shù)據(jù)，是“我們真正想聽到的東西”，它本身含有的“能量”(即該維度的方差，為啥？別急，后文該解釋的時候就有啦)本來應(yīng)該是很大的，但由于它與其他維度有那么一些千絲萬縷的相關(guān)性，受到這些個相關(guān)維度的干擾，它的能量被削弱了，我們就希望通過PCA處理后，使維度A與其他維度的相關(guān)性盡可能減弱，進而恢復(fù)維度A應(yīng)有的能量，讓我們“聽的更清楚”！ · 冗余：冗余也就是多余的意思，就是有它沒它都一樣，放著就是占地方。同樣，假如樣本中有些個維度，在所有的樣本上變化不明顯(極端情況：在所有的樣本中該維度都等于同一個數(shù))，也就是說該維度上的方差接

11、近于零，那么顯然它對區(qū)分不同的樣本絲毫起不到任何作用，這個維度即是冗余的，有它沒它一個樣，所以PCA應(yīng)該去掉這些維度。 這么一分析，那么PCA的最終目的就是“降噪”和消滅這些“冗余”的維度，以使降低維度的同時保存數(shù)據(jù)原有的特征不失真。后面我們將結(jié)合例子繼續(xù)討論。協(xié)方差矩陣PCA實現(xiàn)的關(guān)鍵前面我們說了，PCA的目的就是“降噪”和“去冗余”?！敖翟搿钡哪康木褪鞘贡Ａ粝聛淼木S度間的相關(guān)性盡可能小，而“去冗余”的目的就是使保留下來的維度含有的“能量”即方差盡可能大。那首先，我們得需要知道各維度間的相關(guān)性以及個維度上的方差！那有什么數(shù)據(jù)結(jié)構(gòu)能同時表現(xiàn)不同維度間的相關(guān)性以及各個維度上的方差呢？自然是非協(xié)方

12、差矩陣莫屬。回憶下淺談協(xié)方差矩陣的內(nèi)容，協(xié)方差矩陣度量的是維度與維度之間的關(guān)系，而非樣本與樣本之間。協(xié)方差矩陣的主對角線上的元素是各個維度上的方差(即能量)，其他元素是兩兩維度間的協(xié)方差(即相關(guān)性)。我們要的東西協(xié)方差矩陣都有了，先來看“降噪”，讓保留下的不同維度間的相關(guān)性盡可能小，也就是說讓協(xié)方差矩陣中非對角線元素都基本為零。達到這個目的的方式自然不用說，線代中講的很明確矩陣對角化。而對角化后得到的矩陣，其對角線上是協(xié)方差矩陣的特征值，它還有兩個身份：首先，它還是各個維度上的新方差；其次，它是各個維度本身應(yīng)該擁有的能量(能量的概念伴隨特征值而來)。這也就是我們?yōu)楹卧谇懊娣Q“方差”為“能量”的

13、原因。也許第二點可能存在疑問，但我們應(yīng)該注意到這個事實，通過對角化后，剩余維度間的相關(guān)性已經(jīng)減到最弱，已經(jīng)不會再受“噪聲”的影響了，故此時擁有的能量應(yīng)該比先前大了。看完了“降噪”，我們的“去冗余”還沒完呢。對角化后的協(xié)方差矩陣，對角線上較小的新方差對應(yīng)的就是那些該去掉的維度。所以我們只取那些含有較大能量(特征值)的維度，其余的就舍掉即可。PCA的本質(zhì)其實就是對角化協(xié)方差矩陣。下面就讓我們跟著上面的感覺來推推公式吧。假設(shè)我們有一個樣本集X，里面有N個樣本，每個樣本的維度為d。即：將這些樣本組織成樣本矩陣的形式，即每行為一個樣本，每一列為一個維度，得到樣本矩陣S：。我們先將樣本進行中心化(將樣本中

14、每一個元素都減去樣本的均值。 中心化的意思其實是使樣本數(shù)據(jù)在零點附近波動。)，即保證每個維度的均值為零，只需讓矩陣的每一列除以對應(yīng)的均值即可。很多算法都會先將樣本中心化，以保證所有維度上的偏移都是以零為基點的。然后，對樣本矩陣計算其協(xié)方差矩陣，按照淺談協(xié)方差矩陣里末尾的update，我們知道，協(xié)方差矩陣可以簡單的按下式計算得到：下面，根據(jù)我們上文的推理，將協(xié)方差矩陣C對角化。注意到，這里的矩陣C是是對稱矩陣，對稱矩陣對角化就是找到一個正交矩陣P，滿足：。具體操作是：先對C進行特征值分解，得到特征值矩陣(對角陣)即為，得到特征向量矩陣并正交化即為。顯然，。假如我們?nèi)∽畲蟮那皃(p<d)個特

15、征值對應(yīng)的維度，那么這個p個特征值組成了新的對角陣，對應(yīng)的p個特征向量組成了新的特征向量矩陣。實際上，這個新的特征向量矩陣就是投影矩陣，為什么這么說呢？假設(shè)PCA降維后的樣本矩陣為，顯然，根據(jù)PCA的目的，中的各個維度間的協(xié)方差基本為零，也就是說，的協(xié)方差矩陣應(yīng)該為。即滿足：而我們又有公式：代入可得：由于樣本矩陣的每一行是一個樣本，特征向量矩陣的每一列是一個特征向量。右乘相當于每個樣本以的特征向量為基進行線性變換，得到的新樣本矩陣中每個樣本的維數(shù)變?yōu)榱藀，完成了降維操作。實際上，中的特征向量就是低維空間新的坐標系，稱之為“主成分”。這就是“主成分分析”的名稱由來。同時，的協(xié)方差矩陣為近對角陣，

16、說明不同維度間已經(jīng)基本獨立，噪聲和冗余的數(shù)據(jù)已經(jīng)不見了。至此，整個PCA的過程已經(jīng)結(jié)束，小小總結(jié)一下：1. 形成樣本矩陣，樣本中心化 2. 計算樣本矩陣的協(xié)方差矩陣 3. 對協(xié)方差矩陣進行特征值分解，選取最大的p個特征值對應(yīng)的特征向量組成投影矩陣 4. 對原始樣本矩陣進行投影，得到降維后的新樣本矩陣 Matlab中PCA實戰(zhàn)首先，隨機產(chǎn)生一個10*3維的整數(shù)矩陣作為樣本集，10為樣本的個數(shù)，3為樣本的維數(shù)。1S = fix(rand(10,3)*50);計算協(xié)方差矩陣：1234S = S - repmat(mean(S),10,1); C = (S'*S)./(size(S,1)-1)

17、; or C = cov(S);對協(xié)方差矩陣進行特征值分解：1P,Lambda = eig(C);這里由于三個方差沒有明顯特別小的，所以我們都保留下來，雖然維度沒有降，但觀察Lambda(即PCA后的樣本協(xié)方差矩陣)和C(即原始的協(xié)方差矩陣)，可以發(fā)現(xiàn)，3個維度上的方差有所變化，但對角線之和沒有變，能量重新得到了分配，這就是“降噪”的功勞。最后我們得到降維后的樣本矩陣：1S1 = S*P;為了驗證，我們調(diào)用matlab自帶的主成分分析函數(shù)princomp：1COEFF,SCORE = princomp(S) % COEFF表示投影矩陣，SCORE表示投影后新樣本矩陣對比，可以發(fā)現(xiàn)，SCORE和S1在不考慮維度順序和正負的情況下是完全吻合的，之所以我們計算的S1的維度順序不同，是因為通常都是將投影矩陣P按能量(特征值)的降序排列的，而剛才我們用eig函數(shù)得到的結(jié)