機器學習筆記----四大降維方法之PCA（內帶python及matlab實現）

大家看了之后，可以點一波關注或者推薦一下，以后我也會盡心盡力地寫出好的文章和大家分享。

本文先導：在我們平時看NBA的時候，可能我們只關心球員是否能把球打進，而不太關心這個球的顏色，品牌，只要有3D效果，看到球員扣籃的動作就可以了，比如下圖：

如果我們直接對籃球照片進行幾百萬像素的處理，會有幾千維甚至幾萬維的數據要計算，計算量很大。而往往我們只需要大概勾勒出籃球的大概形狀就可以描述問題，所以必須對此類數據降維，這樣會使處理數據更加輕松。這個在人臉識別中必須要降維，因為我們在做特征提取的時候幾萬維的數據很難全部有效計算，PCA技術作為四大降維技術之一對於此類問題有很好的處理效果。

一 前瞻知識

向量的內積，矩陣的意義，矩陣特征值，正交基，方差，協方差。

 

二 降維的概念

降維通俗地說，就是把高維數據降成低維數據。因為對於機器學習算法，處理上萬級別維數的數據是家常便飯的事情，但是處理是能處理，這會帶來巨大的時間開銷。影響算法的時間復雜度，降維就是在減少維數的前提下，盡量保證數據的完整性。這里需要說明的是，降維不是單純的刪除掉某些特征值，而是把高維數據通過一定的矩陣變換映射到低維空間，現在我們舉一個例子。對於A(3,4),B(5，6),C(3,5)D(-1，2)四個點，如下圖所示：

 

 

我們現在要把他投影在一維坐標上，我想大部分人都是想投影在x軸上，如圖：

等等，不是應該四個點嗎，怎么映射變成了三個點了，如果現實是這樣，那豈不是丟失很多數據？那么我們的目標現在就是找到一條直線，讓這四個點都能在投影在這個直線上，那么數據就不會缺少，而且讓各個點之間的距離最大化。這個有點像我之前說的LDA，鏈接：

http://www.cnblogs.com/xiaohuahua108/p/5985826.html

思路如下：

要想讓各個點映射到直線之后各個點之間的距離最大，我們不妨用方差來表示。

各個點的方差之和為：

    

其中a_i表示各個點的向量表示，u表示各個點的平均坐標的向量表示，一般我們在處理pca的時候，會先減去平均值，得到下面式子：

現在的優化目標就是找到一條直線，使二維上的點投影到這條直線上，使這些點在這條直線上的方差最大。

但是在三維甚至超維的空間中，我們先找到一條直線，讓方差最大，這是第一個方向，接着我們找第二個方向，但是如果還是找方差最大的話，那么這個方向是不是和第一個方向重合了？這樣一直找下去，所有點都集中在一起了，豈不是亂套了？所以我們希望這兩個方向正交，這樣每次找一條直線，都跟原來的方向正交這樣就可以讓數據最大化的顯示在低維空間了。。

三PCA降維方法詳解

1、主成分分析的數學模型及幾何意義

設有n個樣品，每個樣品觀測p項指標(變量)：X₁,X₂, ….X_p，得到原始

數據資料陣：

 

其中

用數據矩陣X的p個向量(即p個指標向量)X₁，X₂，…X_p作線性組合(即

綜合指標向量)為：

簡寫成

其中，Xi是n維向量，所以Fi也是n維向量。上述方程組要求：

且系數a_ij由下列原則決定：

(1) Fi與Fj(i≠j,i,j=1,…p)不相關；

(2) F1是X1 ,X2，…,Xp的一切線性組合(系數滿足上述方程組)中方差最大的，F2與F1不相關的X1 ,X2，…,Xp一切線性組合中方差最大的，…，Fp是與F1，F2，…，Fp-1都不相關的X1 ,X2，…,Xp的一切線性組合中方差最大的。

如何求滿足上述要求的方程組的系數a_ij呢？下一節將會看到每個方程式中的系數向量(a_1i,a_2i, …,a_pi)，i=1,2, …,p不是別的而恰好是X的協差陣∑的特征值所對應的特征向量，也就是說，數學上可以證明使Var(F1)達到最大，這個最大值是在協方差陣∑的第一個特征值所對應特征向量處達到。依此類推使Var(Fp)達到最大值是在協方差陣∑的第p個特征值所對應特征向量處達到。

主成分的幾何意義

從代數學觀點看主成分就是p個變量X₁ ,X₂，…,X_p的一些特殊的線性組合，而在幾何上這些線性組合正是把X₁ ,X₂，…,X_p構成的坐標系旋轉產生的新坐標系，新坐標軸使之通過樣品變差最大的方向(或說具有最大的樣品方差)。下面以最簡單的二元正態變量來說明主成分的幾何意義。

設有n個樣品，每個樣品有p個變量記為X₁ ,X₂，…,X_p，它們的綜合變量記為F₁，F₂，…，F_p 。當p=2時，原變量是X₁，X₂，它們有下圖的相關關系：

主成分的意義

對於二元正態分布變量，n個分散的點大致形成為一個橢圓，若在橢圓長軸方向取坐標軸F₁，在短軸方向聚F₂，這相當於在平面上作一個坐標變換，即按逆時針方向旋轉θ角度，根據旋轉軸變換公式新老坐標之間有關系：

矩陣表示為：

顯然U^T=U^-1且是正交矩陣，即U^TU=I。

從上圖還容易看出二維平面上的n個點的波動(可用方差表示)大部分可以歸結為在F₁軸上的波動，而在F₂軸上的波動是較小的。如果上圖的橢圓是相當扁平的，那么我們可以只考慮F₁方向上的波動，忽略F₂方向的波動。這樣一來，二維可以降為一維了，只取第一個綜合變量F₁即可。而F₁是橢圓的長軸。一般情況，p個變量組成p維空間，n個樣品就是p維空間的n個點，對p元正態分布變量來說，找主成分的問題就是找P維空間中橢球體的主軸問題。

2、主成分分析的推導

在下面推導過程中，要用到線性代數中的兩個定理：

定理一若A是p*p階實對稱陣，則一定可以找到正交陣U使

，其中λ1，λ2，…，λp是A的特征根。

定理二若上述矩陣A的特征根所對應的單位特征向量為u1，u2，…，up令

則實對稱A 屬於不同特征根所對應的特征向量是正交的，即

設，其中a=(a₁，a₂，…，a_p)^T，X=(X₁，X₂，…，X_p)^T，求主成分就是尋找X的線性函數aTX使相應得方差盡可能地大，即使

達到最大值，且a^Ta=1。

設協方差矩陣∑的特征根為l₁，l₂，…，l_p，不妨假設l₁³l₂ ³ …³l_p>0，相應的單位特征向量為u₁， u₂，…， u_p。令

由前面線性代數定理可知：U^TU=UU^T=I，且

因此

所以

而且，當a=u₁時有

因此，a=u₁使Var(a^TX)=a^T∑a達到最大值，且

同理

而且

上述推導表明：X1，X2，…，Xp的主成分就是以E的特征向量為系數的線性組合，它們互不相關，其方差為∑的特征根。

由於∑的特征根l₁³l₂ ³ …³l_p>0，所以有VarF₁³ VarF₂ ³ …³VarF_p>0。了解這一點也就可以明白為什么主成分的名次是按特征根取值大小的順序排列的。

在解決實際問題時，一般不是取p個主成分，而是根據累計貢獻率的大小取前k個。稱第一主成分的貢獻率為，由於有，所以。因此第一主成分的貢獻率就是第一主成分的方差在全部方差中的比值。這個值越大，表明第一主成分綜合X₁，X₂，…，X_p信息的力越強。

前兩個主成分的累計貢獻率定義為，前k個主成分的累計貢獻率定義為。如果前k個主成分的貢獻率達到85%，表明取前 k個主成分包含了全部測量指標所具有的信息，這樣既減少了變量的個數又便於對實際問題進行分析和研究。

3、主成分分析的計算步驟

設x = ( x₁ , x₂ , ⋯, x_n) ^T 為n 維隨機矢量,則PCA具體計算步驟如下：

(1) 將原始觀察數據組成樣本矩陣X ,每一列為一個觀察樣本x ,每一行代表一維數據。

(2) 計算樣本的協方差矩陣：

(3) 計算協方差矩陣C_x 的特征值λ_i 及相應特征向量u_i ,其中i = 1 ,2 , ⋯, n。

(4) 將特征值按由大到小順序排列,並按照下式計算前m 個主元的累積貢獻率：

累積貢獻率用於衡量新生成分量對原始數據的信息保存程度,通常要求其大於85%即可。

(5) 取前m 個較大特征值對應的特征向量構成變換矩陣T^T：

T^T= ( u₁ , u₂ , ⋯, u_m ) 　m < n

(6) 通過Y= TX 計算前m 個主成分,達到降低維數的目的。

四 代碼實現

1、python 2.7代碼實現

from numpy import *
import matplotlib
import matplotlib.pyplot as plt
def loadDataSet(fileName,delim='\t'):
fr=open(fileName)
stringArr=[line.strip().split(delim)for line in fr.readlines()]
dataArr=[map(float,line )for line in stringArr]
return mat(dataArr)
def pca(dataMat,topNfeat=9999999):
meanVals=mean(dataMat,axis=0)
meanRemoved=dataMat-meanVals
covMat=cov(meanRemoved,rowvar=0)
eigvals,eigVects=linalg.eig(mat(covMat))
eigValInd=argsort(eigvals)
eigValInd=eigValInd[:-(topNfeat+1):-1]
redEigVects=eigVects[:,eigValInd]
lowDDataMat=meanRemoved*redEigVects
reconMat=(lowDDataMat*redEigVects.T)+meanVals
return lowDDataMat,reconMat
def display():
import PCA
dataMat=PCA.loadDataSet('testSet.txt')
lowDmat,reconMat=PCA.pca(dataMat,1)
print shape(lowDmat)
print dataMat
print reconMat

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()
if __name__=='__main__':
d=loadDataSet("testSet.txt")
f=pca(d)
print f

上面是python 的代碼實現，經過上面的推導說明大概大家都能看得懂。

實驗效果圖如下：

實驗結果如下：

 

2、matlab代碼實現

E=load('H:\operate\PCA\testSet.txt');

covE=cov(E'); %EΪԭʼ¾ØÕ󣬶ÔE½øÐÐÖ÷³É·Ö·ÖÎö·¨

[v,d]=eig(covE);

%%%%%%%%%%%%%%%%%%%%5

V_eig = v;

temp = sum(d);

D_eig = temp;

[m,n]=size(v);

for i=1:n

V_eig(:,i) = v(:,n-i+1);

end

[m,n]=size(temp);

for i=1:n

D_eig(:,i) = temp(:,n-i+1);

end

D_eig = D_eig';

V=V_eig(:,1:4);%È¡Ç°Èý¸öÖ÷³É·Ö

V=V';

T=V;

結果：

五 PCA的優缺點

   PCA技術的一個很大的優點是，它是完全無參數限制的。在PCA的計算過程中完全不需要人為的設定參數或是根據任何經驗模型對計算進行干預，最后的結果只與數據相關，與用戶是獨立的。 
但是，這一點同時也可以看作是缺點。如果用戶對觀測對象有一定的先驗知識，掌握了數據的一些特征，卻無法通過參數化等方法對處理過程進行干預，可能會得不到預期的效果，效率也不高。

六 說明

由於今天和女朋友吵架，把女朋友惹生氣了。我心情也不是很好，在匆忙之中完成這篇博客，有什么錯誤地方還請大家指正。么么噠，愛你們的小花。

只想在這里和我漂亮可愛善良的女朋友說一句，我愛你。