PCA降維在MATLAB上的實(shí)現(xiàn)

PCA降維在MATLAB上的實(shí)現(xiàn)學(xué)院計(jì)算機(jī)科學(xué)與技術(shù)專業(yè)計(jì)算機(jī)科學(xué)與技術(shù)年級(jí)2011級(jí)姓名周忠儒/王云標(biāo)學(xué)號(hào)30111216058/051指導(dǎo)教師魏建國(guó)2014年5月28日PCA降維在MATLAB上的實(shí)現(xiàn)TOC\o"1-3"\h\u32538一實(shí)驗(yàn)?zāi)康?23986二實(shí)驗(yàn)環(huán)境 221415三實(shí)驗(yàn)原理 2105891、PCA降維方法原理 2314952、MATLAB 391093、PCA降維方法詳解 3203501）、原始數(shù)據(jù)： 383142）、協(xié)方差矩陣的求法： 479513）、計(jì)算協(xié)方差矩陣的特征向量和特征值： 6162964）、選擇成分組成模式矢量： 6302075）、得到降維后的數(shù)據(jù)： 7769四實(shí)驗(yàn)代碼詳解 78857五實(shí)驗(yàn)結(jié)果 811239六實(shí)驗(yàn)總結(jié) 9PCA降維在MATLAB上的實(shí)現(xiàn)一實(shí)驗(yàn)?zāi)康?掌握PCA降維的基本內(nèi)容2了解MATLAB的基本用法3用PCA降維算法處理圖像數(shù)據(jù)二實(shí)驗(yàn)環(huán)境M三實(shí)驗(yàn)原理1、PCA降維方法原理PCA的原理就是將原來的樣本數(shù)據(jù)投影到一個(gè)新的空間中，相當(dāng)于我們?cè)诰仃嚪治隼锩鎸W(xué)習(xí)的將一組矩陣映射到另外的坐標(biāo)系下。通過一個(gè)轉(zhuǎn)換坐標(biāo)，也可以理解成把一組坐標(biāo)轉(zhuǎn)換到另外一組坐標(biāo)系下，但是在新的坐標(biāo)系下，表示原來的原本不需要那么多的變量，只需要原來樣本的最大的一個(gè)線性無關(guān)組的特征值對(duì)應(yīng)的空間的坐標(biāo)即可。PCA即主成分分析，是圖像處理中經(jīng)常用到的降維方法，大家知道，我們?cè)谔幚碛嘘P(guān)數(shù)字圖像處理方面的問題時(shí)，比如經(jīng)常用的圖像的查詢問題，在一個(gè)幾萬(wàn)或者幾百萬(wàn)甚至更大的數(shù)據(jù)庫(kù)中查詢一幅相近的圖像。這時(shí)，我們通常的方法是對(duì)圖像庫(kù)中的圖片提取響應(yīng)的特征，如顏色，紋理，sift，surf，vlad等等特征，然后將其保存，建立響應(yīng)的數(shù)據(jù)索引，然后對(duì)要查詢的圖像提取相應(yīng)的特征，與數(shù)據(jù)庫(kù)中的圖像特征對(duì)比，找出與之最近的圖片。2、MATLABMATLAB（矩陣實(shí)驗(yàn)室）是MATrixLABoratory的縮寫，是一款由美國(guó)TheMathWorks公司出品的商業(yè)數(shù)學(xué)軟件。MATLAB是一種用于算法開發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計(jì)算的高級(jí)技術(shù)計(jì)算語(yǔ)言和交互式環(huán)境。除了矩陣運(yùn)算、繪制函數(shù)/數(shù)據(jù)圖像等常用功能外，MATLAB還可以用來創(chuàng)建用戶界面及與調(diào)用其它語(yǔ)言（包括C，C++和FORTRAN）編寫的程序。3、PCA降維方法詳解1）、原始數(shù)據(jù)：由于本實(shí)驗(yàn)數(shù)據(jù)過于龐大，為了方便，我們假定數(shù)據(jù)是二維的，借助網(wǎng)絡(luò)上的一組數(shù)據(jù)，如下：x=[2.5,0.5,2.2,1.9,3.1,2.3,2,1,1.5,1.1]Ty=[2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9]T2）、計(jì)算協(xié)方差矩陣（1）協(xié)方差矩陣：以下是含有n個(gè)樣本的集合中，一些數(shù)理統(tǒng)計(jì)的相關(guān)概念：均值：標(biāo)準(zhǔn)差：方差：在這里，標(biāo)準(zhǔn)差和方差一般是用來描述一維數(shù)據(jù)的，但現(xiàn)實(shí)生活我們常常遇到含有多維數(shù)據(jù)的數(shù)據(jù)集，例如上學(xué)時(shí)免不了要統(tǒng)計(jì)多個(gè)學(xué)科的考試成績(jī)。面對(duì)這樣的數(shù)據(jù)集，我們當(dāng)然可以按照每一維獨(dú)立的計(jì)算其方差，但是通常我們還想了解這幾科成績(jī)之間的關(guān)系，這時(shí)，我們就要用協(xié)方差，協(xié)方差就是一種用來度量?jī)蓚€(gè)隨機(jī)變量關(guān)系的統(tǒng)計(jì)量，其定義為：從協(xié)方差的定義上我們也可以看出一些顯而易見的性質(zhì)，如：cov(X,Y)=var(X);cov(X,Y)=var(Y,X);需要注意的是，協(xié)方差也只能處理二維問題，那維數(shù)多了自然就需要計(jì)算多個(gè)協(xié)方差，比如n維的數(shù)據(jù)集就需要計(jì)算個(gè)協(xié)方差，那自然而然的我們會(huì)想到使用矩陣來組織這些數(shù)據(jù)。給出協(xié)方差矩陣的定義：這個(gè)定義還是很容易理解的，我們可以舉一個(gè)簡(jiǎn)單的三維的例子，假設(shè)數(shù)據(jù)集有三個(gè)維度，則協(xié)方差矩陣為：可見，協(xié)方差矩陣是一個(gè)對(duì)稱的矩陣，而且對(duì)角線是各個(gè)維度上的方差。2）、協(xié)方差矩陣的求法：協(xié)方差矩陣計(jì)算的是不同維度之間的協(xié)方差，而不是不同樣本之間的。下面我們將在matlab中用一個(gè)例子進(jìn)行詳細(xì)說明：首先，隨機(jī)產(chǎn)生一個(gè)10*3維的整數(shù)矩陣作為樣本集，10為樣本的個(gè)數(shù)，3為樣本的維數(shù)。MySample=fix(rand(10,3)*50)根據(jù)公式，計(jì)算協(xié)方差需要計(jì)算均值，那是按行計(jì)算均值還是按列呢，我一開始就老是困擾這個(gè)問題。前面我們也特別強(qiáng)調(diào)了，協(xié)方差矩陣是計(jì)算不同維度間的協(xié)方差，要時(shí)刻牢記這一點(diǎn)。樣本矩陣的每行是一個(gè)樣本，每列為一個(gè)維度，所以我們要按列計(jì)算均值。為了描述方便，我們先將三個(gè)維度的數(shù)據(jù)分別賦值：dim1=MySample(:,1);dim2=MySample(:,2);dim3=MySample(:,3);計(jì)算dim1與dim2，dim1與dim3，dim2與dim3的協(xié)方差：sum((dim1-mean(dim1)).*(dim2-mean(dim2)))/(size(MySample,1)-1)sum((dim1-mean(dim1)).*(dim3-mean(dim3)))/(size(MySample,1)-1)

sum((dim2-mean(dim2)).*(dim3-mean(dim3)))/(size(MySample,1)-1)

搞清楚了這個(gè)后面就容易多了，協(xié)方差矩陣的對(duì)角線就是各個(gè)維度上的方差，下面我們依次計(jì)算：std(dim1)^2std(dim2)^2

std(dim3)^2

這樣，我們就得到了計(jì)算協(xié)方差矩陣所需要的所有數(shù)據(jù)，調(diào)用Matlab自帶的cov函數(shù)進(jìn)行驗(yàn)證：cov(MySample)可以看到跟我們計(jì)算的結(jié)果是一樣的，說明我們的計(jì)算是正確的。但是通常我們不用這種方法，而是用下面簡(jiǎn)化的方法進(jìn)行計(jì)算：1，先讓樣本矩陣中心化，即每一維度減去該維度的均值。2，然后直接用新的樣本矩陣乘上它的轉(zhuǎn)置。3，然后除以(N-1)即可。其實(shí)這種方法也是由前面的公式通道而來，只不過理解起來不是很直觀而已。其Matlab代碼實(shí)現(xiàn)如下：X=MySample–repmat(mean(MySample),10,1);%中心化樣本矩陣C=(X’*X)./(size(X,1)-1)其中B=repmat(A,m,n)%將矩陣A復(fù)制m×n塊。即把A作為B的元素，B由m×n個(gè)A平鋪而成。B的維數(shù)是[size(A,1)*m,(size(A,2)*n]

B=mean(A)的說明：如果你有這樣一個(gè)矩陣：A=[123;336;468;477];用mean(A)（默認(rèn)dim=1）就會(huì)求每一列的均值ans=

6.0000用mean(A,2)就會(huì)求每一行的均值

ans=

size(A,n)%如果在size函數(shù)的輸入?yún)?shù)中再添加一項(xiàng)n，并用1或2為n賦值，則size將返回矩陣的行數(shù)或列數(shù)。其中r=size(A,1)該語(yǔ)句返回的是矩陣A的行數(shù)，c=size(A,2)該語(yǔ)句返回的是矩陣A的列數(shù)）上面我們簡(jiǎn)單說了一下協(xié)方差矩陣及其求法，言歸正傳，我們用上面簡(jiǎn)化求法，求出樣本的協(xié)方差矩陣為：3）、計(jì)算協(xié)方差矩陣的特征向量和特征值因?yàn)閰f(xié)方差矩陣為方陣，我們可以計(jì)算它的特征向量和特征值，如下：[eigenvectors,eigenvalues]=eig(cov)我們可以看到這些矢量都是單位矢量，也就是它們的長(zhǎng)度為1，這對(duì)PCA來說是很重要的。4）、選擇成分組成模式矢量求出協(xié)方差矩陣的特征值及特征向量之后，按照特征值由大到小進(jìn)行排列，這將給出成分的重要性級(jí)別。現(xiàn)在，如果你喜歡，可以忽略那些重要性很小的成分，當(dāng)然這會(huì)丟失一些信息，但是如果對(duì)應(yīng)的特征值很小，你不會(huì)丟失很多信息。如果你已經(jīng)忽略了一些成分，那么最后的數(shù)據(jù)集將有更少的維數(shù)，精確地說，如果你的原始數(shù)據(jù)是n維的，你選擇了前p個(gè)主要成分，那么你現(xiàn)在的數(shù)據(jù)將僅有p維。現(xiàn)在我們要做的是組成一個(gè)模式矢量，這只是幾個(gè)矢量組成的矩陣的一個(gè)有意思的名字而已，它由你保持的所有特征矢量構(gòu)成，每一個(gè)特征矢量是這個(gè)矩陣的一列。對(duì)于我們的數(shù)據(jù)集，因?yàn)橛袃蓚€(gè)特征矢量，因此我們有兩個(gè)選擇。我們可以用兩個(gè)特征矢量組成模式矢量：我們也可以忽略其中較小特征值的一個(gè)特征矢量，從而得到如下模式矢量：、得到降維后的數(shù)據(jù)其中rowFeatureVector是由模式矢量作為列組成的矩陣的轉(zhuǎn)置，因此它的行就是原來的模式矢量，而且對(duì)應(yīng)最大特征值的特征矢量在該矩陣的最上一行。rowdataAdjust是每一維數(shù)據(jù)減去均值后，所組成矩陣的轉(zhuǎn)置，即數(shù)據(jù)項(xiàng)目在每一列中，每一行是一維，對(duì)我們的樣本來說即是，第一行為x維上數(shù)據(jù)，第二行為y維上的數(shù)據(jù)。FinalData是最后得到的數(shù)據(jù)，數(shù)據(jù)項(xiàng)目在它的列中，維數(shù)沿著行。這將給我們什么結(jié)果呢？這將僅僅給出我們選擇的數(shù)據(jù)。我們的原始數(shù)據(jù)有兩個(gè)軸（x和y），所以我們的原始數(shù)據(jù)按這兩個(gè)軸分布。我們可以按任何兩個(gè)我們喜歡的軸表示我們的數(shù)據(jù)。如果這些軸是正交的，這種表達(dá)將是最有效的，這就是特征矢量總是正交的重要性。我們已經(jīng)將我們的數(shù)據(jù)從原來的xy軸表達(dá)變換為現(xiàn)在的單個(gè)特征矢量表達(dá)。說明：如果要恢復(fù)原始數(shù)據(jù)，只需逆過程計(jì)算即可，即：四實(shí)驗(yàn)代碼詳解按照上述PCA降維算法的原理編寫代碼如下：將im中的數(shù)據(jù)導(dǎo)入MATLAB；m=mean(im,2);求im數(shù)據(jù)每一行的平均值，即得數(shù)為一列向量，每一元素為im每一行的平均值。Train_Number=size(im,2);求im數(shù)據(jù)每一行數(shù)據(jù)的元素個(gè)數(shù)A=[];定義一個(gè)空矩陣Afori=1:Train_Numbertemp=im(:,i)-uint8(m);定義一個(gè)中間矩陣變量，其值為im中每一列減去m的得數(shù)A=[Atemp];將中間矩陣變量temp的值存入A中endfor循環(huán)的作用是將im每一列都減去m并存入A中L=A'*A;L=A矩陣與其轉(zhuǎn)置矩陣的乘積[VD]=eig(L);求矩陣L的特征值與特征向量L_eig_vec=[];定義一個(gè)空矩陣fori=1:15L_eig_vec=[L_eig_vecV(:,i)];End將特征矩陣存入L中PCA=A*L_eig_vec;PCA為A矩