相關文章

PCA的原理及MATLAB實現

UFLDL教程：Exercise:PCA in 2D & PCA and Whitening

python-A comparison of various Robust PCA implementations

－－－－－－－－－－－－－－－－－－－－－－－－－－－－

本文參考：
http://blog.csdn.net/gwh111/article/details/11742735
http://www.mathworks.com/matlabcentral/newsreader/view_thread/326056
https://www.zhihu.com/question/37094362
http://blog.csdn.net/zhaoxinfan/article/details/9138165
http://blog.csdn.net/stdcoutzyx/article/details/37568225
http://blog.csdn.net/ice110956/article/details/14250745
http://blog.csdn.net/xl890727/article/details/16898315
http://www.cnblogs.com/tornadomeet/archive/2012/09/06/2673104.html
http://blog.csdn.net/llp1992/article/details/45065609
http://www.cnblogs.com/tornadomeet/archive/2012/12/30/2839615.html
http://blog.csdn.net/ice110956/article/details/20936351?utm_source=tuicool&utm_medium=referral

關于PCA原理可以直接參考下面的文章
深入理解PCA

PCA是經常用來減少數據集的維數，同時保留數據集中對方差貢獻最大的特征來達到簡化數據集的目的。
PCA的原理就是將原來的樣本數據投影到一個新的空間中，相當于我們在矩陣分析里面學習的將一組矩陣映射到另外的坐標系下。通過一個轉換坐標，也可以理解成把一組坐標轉換到另外一組坐標系下，但是在新的坐標系下，表示原來的原本不需要那么多的變量，只需要原來樣本的最大的一個線性無關組的特征值對應的空間的坐標即可。
PCA算法核心思想就是將 n 維特征映射到 k 維上（k < n），這 k 維是全新的正交特征。我們將這 k 維成為主元，是重新構造出來的 k 維特征，而不是簡單地從 n 維特征中取出其余 n-k 維特征。

PCA實現的步驟
（1）把原始數據中每個樣本用一個向量表示，然后把所有樣本組合起來構成一個矩陣。為了避免樣本的單位的影響，樣本集需要標準化（一般都是去均值化）。
（2）求該矩陣的協方差矩陣。
（3）求步驟2中得到的協方差矩陣的特征值和特征向量。
（4）將求出的特征向量按照特征值的大小進行組合形成一個映射矩陣，并根據指定的PCA保留的特征個數取出映射矩陣的前n行或者前n列作為最終的映射矩陣。
（5）用步驟4的映射矩陣對原始數據進行映射，達到數據降維的目的。

PCA理論的matlab實現

function [Y,V,E,D] = newpca(X) %其中X為輸入數據，X的每一列是一個輸入樣本。返回值Y是對X進行PCA分析后的投影矩陣。 %V是與X有關的協方差矩陣特征向量的白化矩陣，E是對應的特征向量（列）構成的矩陣， %D是對應的特征值構成的對角矩陣（特征值處于對角線上）。 %返回值中的白化矩陣，特征向量和特征值都是按照對應特征值大小進行排序后了的。 % do PCA on image patches % % INPUT variables: % X matrix with image patches as columns % % OUTPUT variables: % Y the project matrix of the input data X without whiting % V whitening matrix % E principal component transformation (orthogonal) % D variances of the principal components%去除直流成分 X = X-ones(size(X,1),1)*mean(X);% Calculate the eigenvalues and eigenvectors of the new covariance matrix. covarianceMatrix = X*X'/size(X,2); %求出其協方差矩陣 %E是特征向量構成，它的每一列是特征向量，D是特征值構成的對角矩陣 %這些特征值和特征向量都沒有經過排序 [E, D] = eig(covarianceMatrix); % Sort the eigenvalues and recompute matrices % 因為sort函數是升序排列，而需要的是降序排列，所以先取負號,diag(a)是取出a的對角元素構成 % 一個列向量，這里的dummy是降序排列后的向量，order是其排列順序 [dummy,order] = sort(diag(-D)); E = E(:,order);%將特征向量按照特征值大小進行降序排列，每一列是一個特征向量 Y = E'*X; d = diag(D); %d是一個列向量 %dsqrtinv是列向量，特征值開根號后取倒,仍然是與特征值有關的列向量 %其實就是求開根號后的逆矩陣 dsqrtinv = real(d.^(-0.5)); Dsqrtinv = diag(dsqrtinv(order));%是一個對角矩陣，矩陣中的元素時按降序排列好了的特征值（經過取根號倒后） D = diag(d(order));%D是一個對角矩陣，其對角元素由特征值從大到小構成 V = Dsqrtinv*E';%特征值矩陣乘以特征向量矩陣

示例：二維數據的實現過程

function [lowData,reconMat] = PCA(data,K) [row , col] = size(data); meanValue = mean(data); %varData = var(data,1,1); normData = data - repmat(meanValue,[row,1]); covMat = cov(normData(:,1),normData(:,2));%求取協方差矩陣 [eigVect,eigVal] = eig(covMat);%求取特征值和特征向量 [sortMat, sortIX] = sort(eigVal,'descend'); [B,IX] = sort(sortMat(1,:),'descend'); len = min(K,length(IX)); eigVect(:,IX(1:1:len)); lowData = normData * eigVect(:,IX(1:1:len)); reconMat = (lowData * eigVect(:,IX(1:1:len))') + repmat(meanValue,[row,1]); % 將降維后的數據轉換到新空間 end

下面給出PCA的Matlab實現部分：Matlab自帶的函數
Matlab中已經包含了實現了的PCA算法，可以通過princomp函數調用。其形式為：
[COEFF,SCORE, latent]=princomp(X);
X：為要輸入的n維原始數據。帶入這個matlab自帶函數，將會生成新的n維加工后的數據（即SCORE）。此數據與之前的n維原始數據一一對應。
COEFF為主成分分量，即變換空間中的那些基向量，是系數矩陣。通過COEFF可以知道x是怎樣轉換成SCORE的。
在n行p列的數據集X上做主成分分析。返回主成分系數。X的每行表示一個樣本的觀測值，每一列表示特征變量。COEFF是一個p行p列的矩陣，是X矩陣所對應的協方差陣V的所有特征向量組成的矩陣，即變換矩陣或稱投影矩陣，COEFF每列對應一個特征值的特征向量，列的排列順序是按特征值的大小遞減排序。
SCORE為主成分，即X的低維表示。生成的n維加工后的數據存在SCORE里。它是對原始數據進行的分析，進而在新的坐標系下獲得的數據。他將這n維數據按貢獻率由大到小排列。（即在改變坐標系的情況下，又對n維數據排序）
latent為一個包含樣本協方差矩陣的本征值的向量。是一維列向量，每一個數據是對應SCORE里相應維的貢獻率，因為數據有n維所以列向量有n個數據。由大到小排列（因為SCORE也是按貢獻率由大到小排列）。
princomp這個matlab自帶的函數，在降維之前就將每一個樣本減去了一個所有樣本的平均值。princomp這里使用一行表示一個樣本，每行包括這個樣本的所有的特征值。

其實我們要的是由函數[COEFF,SCORE, latent] = princomp(X)它所產生的pc和latent。由latent可以算出降維后的空間所能表示原空間的程度，只要這個累積的值大于95%就行了
cumsum(latent)./sum(latent)

舉例說明：

load hald; %載入matlab內部數據 [COEFF,SCORE,latent,tsquare] = princomp(ingredients); %調用pca分析函數 % 下面我們具體說明COEFF，SCORE，latent %下面為計算ingredients所對應的協方差矩陣(也就是cov_ingredients矩陣)的特征值和特征向量，下面的矩陣V為特征向量，D為特征值(對比上面的latent)組成的對角線矩陣 [V,D] = eig(cov_ingredients)V =0.5062 0.5673 0.6460 -0.06780.4933 -0.5440 0.0200 -0.67850.5156 0.4036 -0.7553 0.02900.4844 -0.4684 0.1085 0.7309D =0.2372 0 0 00 12.4054 0 00 0 67.4964 00 0 0 517.7969 %說明1：對比一下矩陣V和矩陣COEFF，現在很容易明白為什么COEFF是按列遞減順序排列的了！(V中第三列與COEFF中倒數第三列差個負號，學過線性代數的人都知道這沒問題) %下面再驗證一下說明2 diag(cov(SCORE))ans =517.796967.496412.40540.2372 %說明2：以上結果顯示latent確實表示SCORE矩陣每列的方差，517.7969表示第一列方差[COEFF,SCORE,latent,tsquare] = princomp(ingredients); % 上面這個函數，COEFF矩陣是返回的轉換矩陣，也就是把樣本轉換到新的空間中的準換矩陣，這個準換矩陣式比較大的。 % SCORE是原來的樣本矩陣在新的坐標系中的表示，也就是原來的樣本乘上轉換矩陣，但是,還不是直接乘，要減去一個樣本的均值。將原來的數據轉換到新的樣本空間中的算法是這樣實現的： x0 = bsxfun(@minus,ingredients,mean(ingredients,1)); NewSCORE1 = x0 * COEFF;NewSCORE1 =36.8218 -6.8709 -4.5909 0.396729.6073 4.6109 -2.2476 -0.3958-12.9818 -4.2049 0.9022 -1.1261 23.7147 -6.6341 1.8547 -0.3786-0.5532 -4.4617 -6.0874 0.1424 -10.8125 -3.6466 0.9130 -0.1350 -32.5882 8.9798 -1.6063 0.0818 22.6064 10.7259 3.2365 0.3243-9.2626 8.9854 -0.0169 -0.5437 -3.2840 -14.1573 7.0465 0.3405 9.2200 12.3861 3.4283 0.4352-25.5849 -2.7817 -0.3867 0.4468 -26.9032 -2.9310 -2.4455 0.4116 NewSCORE2 = ingredients * COEFF;NewSCORE2 =25.9105 -7.0189 -35.8545 48.526618.6960 4.4628 -33.5112 47.7341-23.8931 -4.3530 -30.3613 47.0039 12.8034 -6.7821 -29.4088 47.7514-11.4645 -4.6098 -37.3510 48.2723 -21.7238 -3.7946 -30.3506 47.9950 -43.4995 8.8318 -32.8699 48.2117 11.6951 10.5779 -28.0270 48.4543-20.1739 8.8373 -31.2805 47.5862 -14.1953 -14.3053 -24.2171 48.4705 -1.6913 12.2380 -27.8352 48.5651 -36.4962 -2.9297 -31.6503 48.5768 -37.8145 -3.0790 -33.7091 48.5416 %對比NewSCORE1和NewSCORE2 與 SCORE，可以看成 SCORE與NewSCORE1是一樣。 % 下面我們就來確定最后的轉換矩陣 cumsum(latent)./sum(latent) ans =0.86600.97890.99961.0000SelectNum = cumsum(latent)./sum(latent)index = find(SelectNum >= 0.95);ForwardNum = index(1); %選擇由latent可以算出降維后的空間所能表示原空間的程度，只要這個累積的值大于95%的數目ForwardNum，就選擇前ForwardNum個主成分作為變換矩陣。 %做主成分變換矩陣，也就是降維矩陣 tranMatrix = COEFF(:,1:ForwardNum)tranMatrix =-0.0678 -0.6460 -0.6785 -0.0200 0.0290 0.75530.7309 -0.1085 % Newingredients = x0 * tranMatrix Newingredients =36.8218 -6.870929.6073 4.6109-12.9818 -4.2049 23.7147 -6.6341-0.5532 -4.4617 -10.8125 -3.6466 -32.5882 8.9798 22.6064 10.7259-9.2626 8.9854 -3.2840 -14.1573 9.2200 12.3861-25.5849 -2.7817 -26.9032 -2.9310 % 這樣Newingredients就是ingredients經過PCA降維的最終的結果。

測試樣本降維
如果你需要對測試樣本降維，一般情況下，使用matlab自帶的方式，肯定需要對測試樣本減去一個訓練樣本均值，因為你在給訓練樣本降維的時候減去了均值，所以測試樣本也要減去均值，然后乘以coeff這個矩陣，就獲得了測試樣本降維后的數據。使用訓練樣本得到的轉換矩陣，保證訓練樣本和測試樣本轉換到相同的樣本空間中。比如說你的測試樣本是1*1000000，那么乘上一個1000000*29的降維矩陣，就獲得了1*29的降維后的測試樣本的降維數據。

princomp(x)使用的行表示一個樣本，每行的所有的列數據都是這個樣本的特征值。降維以后比如是30*29，那么每一行就是降維以后的數據。每個樣本有29個特征值。

%X為訓練集，Xtest為測試集 [COEFF,SCORE, latent] = princomp(X);SelectNum = cumsum(latent)./sum(latent)index = find(SelectNum >= 0.95);ForwardNum = index(1); %選擇由latent可以算出降維后的空間所能表示原空間的程度，只要這個累積的值大于95%的數目ForwardNum，就選擇前ForwardNum個主成分作為變換矩陣。 %做主成分變換矩陣 NewtrainScores = SCORE(:,1:ForwardNum); **NewtrainScores就是X經過PCA降維后的數據**% tranMatrix = COEFF(:,1:ForwardNum) [row , col] = size(Xtest); meanValue = mean(X); normXtest = Xtest - repmat(meanValue,[row,1]); % 第二種方法 %normXtest = bsxfun(@minus,Xtest,meanValue); % center the test data NewtestScores = normXtest*tranMatrix; %**NewtestScores 就是Xtest經過PCA降維后的數據。**

備注說明：
對測試樣本進行降維的時候，一定要減去訓練樣本的均值，使用訓練樣本得到的轉換矩陣，保證訓練樣本和測試樣本轉換到相同的樣本空間中，這樣才有意義。
很多人就選擇還是使用princomp這個函數處理測試樣本，那么這樣測試樣本被映射到一個新的空間中，和原來的訓練樣本完全不是在一個空間，一點意義都沒有，還是要使用測試樣本減去訓練樣本的均值，然后乘上訓練樣本降維的時候獲得降維矩陣，轉換到相同的空間中。

project the test data on the reduced set of the principal components. Note that princomp centers data before computing the principal components. You need to center the test data as well. It is best to do that using the means obtained from the training data. Something like this would do:mu = mean(Xtrain); [coefs,scores,variances] = princomp(Xtrain,'econ'); % do your thing to select the components Xtest = bsxfun(@minus,Xtest,mu); % center the test data testScores = Xtest*coefs(:,1:d);

具體詳情大家可以參考
http://blog.csdn.net/gwh111/article/details/11742735
http://www.mathworks.com/matlabcentral/newsreader/view_thread/326056

總結

以上是生活随笔為你收集整理的PCA的原理及MATLAB实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。