Matlab自帶PCA分析
阿新 • • 發佈：2019-01-06
本文轉載自博主：blog.csdn.net/watkinsong/article
ＰＣＡ原理：
PCA的原理就是將原來的樣本資料投影到一個新的空間中，相當於我們在矩陣分析裡面學習的將一組矩陣對映到另外的座標系下。通過一個轉換座標，也可以理解成把一組座標轉換到另外一組座標系下，但是在新的座標系下，表示原來的原本不需要那麼多的變數，只需要原來樣本的最大的一個線性無關組的特徵值對應的空間的座標即可。
比如，原來的樣本是30*1000000的維數，就是說我們有30個樣本，每個樣本有1000000個特徵點，這個特徵點太多了，我們需要對這些樣本的特徵點進行降維。那麼在降維的時候會計算一個原來樣本矩陣的協方差矩陣，這裡就是1000000*1000000，當然，這個矩陣太大了，計算的時候有其他的方式進行處理，這裡只是講解基本的原理，然後通過這個1000000*1000000的協方差矩陣計算它的特徵值和特徵向量，最後獲得具有最大特徵值的特徵向量構成轉換矩陣。比如我們的前29個特徵值已經能夠佔到所有特徵值的99%以上，那麼我們只需要提取前29個特徵值對應的特徵向量即可。這樣就構成了一個1000000*29的轉換矩陣，然後用原來的樣本乘以這個轉換矩陣，就可以得到原來的樣本資料在新的特徵空間的對應的座標。30*1000000 * 1000000*29 = 30 *29， 這樣原來的訓練樣本每個樣本的特徵值的個數就降到了29個。
一般來說，PCA降維後的每個樣本的特徵的維數，不會超過訓練樣本的個數，因為超出的特徵是沒有意義的。
下面是百度百科中對pca降維的一段解釋，還是挺清晰的：
“對於一個訓練集，100個物件模板，特徵是10維，那麼它可以建立一個100*10的矩陣，作為樣本。求這個樣本的協方差矩陣，得到一個10*10的協方差矩陣，然後求出這個協方差矩陣的特徵值和特徵向量，應該有10個特徵值和特徵向量，我們根據特徵值的大小，取前四個特徵值所對應的特徵向量，構成一個10*4的矩陣，這個矩陣就是我們要求的特徵矩陣，100*10的樣本矩陣乘以這個10*4的特徵矩陣，就得到了一個100*4的新的降維之後的樣本矩陣，每個特徵的維數下降了。
當給定一個測試的特徵集之後，比如1*10維的特徵，乘以上面得到的10*4的特徵矩陣，便可以得到一個1*4的特徵，用這個特徵去分類。我對 ＰＣＡ的一些瞭解
我的pca迷惑
迷惑一
剛開始接觸ＰＣＡ的時候，諮詢了一個浙大的博士朋友，這朋友告訴我，如果對訓練樣本進行降維，那麼樣本的數量必須大於特徵的維數，然後我當時就迷惑了，那我怎麼辦啊，我的人臉表情影象頂多有幾百張就算多的了，但是每個影象提取的特徵的維數將近有幾十萬，我不可能找那麼多樣本去啊。當時有這個迷惑也是因為matlab給出的一個實現在pca降維的函式的說明，就是princomp，這個函式的說明也是用的樣本的個數多餘特徵的維數。後來經過試驗是證實，證實了那個浙大的博士的認識是錯誤的，pca降維肯定不需要樣本的個數大於特徵的維數，要不然還降維個什麼意思。比如我有30*1000000的特徵矩陣，那麼降維後肯定是每個樣本在新的空間中的表示的特徵維數不超過30.迷惑二
迷惑三
網上有不同的pca降維的程式碼，每個程式碼也實現的不一樣，那麼對於同一個資料是否是pca降維以後都是獲得相同的資料呢，也就是說不管你用哪種方式進行pca降維，不管你是從哪裡下載到的或者自己根據演算法實現的pca降維，同樣的矩陣降維以後的資料是否一致？這個我個人認為，不同的演算法最後導致的pca降維的資料肯定不一致。因為pca降維以後，只是把原來的資料對映到新的特徵空間，所以如果你的演算法不同，那麼選擇的協方差矩陣肯定就不同，最後獲得的轉換矩陣肯定也不一樣。那麼訓練樣本和測試樣本和不同的轉換矩陣相乘以後最終肯定會獲得不同的降維座標。所以使用不同的演算法應該最後不會有相同的座標結果，這個也是我一直實驗的結果，我也使用了matlab自帶的princomp降維，並且使用相同的資料使用網上下載的一些降維方法進行降維，得到的資料都不一致。 比如說princomp這個matlab自帶的函式，在降維之前就將每一個樣本減去了一個所有樣本的平均值，也可能有很多樣本沒有減去平均值。princomp這裡使用一行表示一個樣本，每行包括這個樣本的所有的特徵值。而網上大部分都是每一列表示一個樣本，這樣這一列的所有行都表示這個樣本的特徵值。網上的程式使用列表示樣本是有一定好處的，比如我的樣本是1000000*30，總共有30個訓練樣本，每個樣本的特徵值個數是1000000，那麼這個矩陣獲得的協方差矩陣是30*30，計算起來非常的方便，不想30*1000000這樣的矩陣獲得到的協方差矩陣式1000000*1000000，直接就記憶體溢位了，不過matlab有自己的實現方式，巧妙的解決了這個問題。pca的實現（matlab）
我在網上看了很多pca降維的例子，都大同小異，原理差不多，都是活的原來矩陣的協方差矩陣，然後計算協方差矩陣的特徵值和特徵向量，最後通過特徵向量的根據特徵值由大到小的排序進行KL變換神馬的獲得一個轉換矩陣。1. matlab自帶的實現方式
PCA在matlab中的實現舉例 以下資料來自matlab的help，翻譯和註解部分由筆者新增：(重點部分添加了翻譯！) princomp-----函式名稱 Principal component analysis (PCA) on data Syntax------函式呼叫語法 [COEFF,SCORE] = princomp(X) [COEFF,SCORE,latent] = princomp(X) [COEFF,SCORE,latent,tsquare] = princomp(X) [...] = princomp(X,'econ') Description -----函式描述 COEFF = princomp(X) performs principal components analysis (PCA) on the n-by-p data matrix X, and returns the principal component coefficients, also known as loadings. Rows of X correspond to observations, columns to variables. COEFF is a p-by-p matrix, each column containing coefficients for one principal component. The columns are in order of decreasing component variance. 在n行p列的資料集X上做主成分分析。返回主成分系數。X的每行表示一個樣本的觀測值，每一列表示特徵變數。COEFF是一個p行p列的矩陣，每一列包含一個主成分的係數，列是按主成分變數遞減順序排列。(按照這個翻譯很難理解，其實COEFF是X矩陣所對應的協方差陣V的所有特徵向量組成的矩陣，即變換矩陣或稱投影矩陣，COEFF每列對應一個特徵值的特徵向量，列的排列順序是按特徵值的大小遞減排序，後面有具體例子解釋，見說明1) princomp centers X by subtracting off column means, but does not rescale the columns of X. To perform principal components analysis with standardized variables, that is, based on correlations, use princomp(zscore(X)). To perform principal components analysis directly on a covariance or correlation matrix, use pcacov. 計算PCA的時候，MATLAB自動對列進行了去均值的操作，但是並不對資料進行規格化，如果要規格化的話，用princomp(zscore(X))。另外，如果直接有現成的協方差陣，用函式pcacov來計算。 [COEFF,SCORE] = princomp(X) returns SCORE, the principal component scores; that is, the representation of X in the principal component space. Rows of SCORE correspond to observations, columns to components. 返回的SCORE是對主分的打分，也就是說原X矩陣在主成分空間的表示。SCORE每行對應樣本觀測值，每列對應一個主成份(變數)，它的行和列的數目和X的行列數目相同。 [COEFF,SCORE,latent] = princomp(X) returns latent, a vector containing the eigenvalues of the covariance matrix of X. 返回的latent是一個向量，它是X所對應的協方差矩陣的特徵值向量。 [COEFF,SCORE,latent,tsquare] = princomp(X) returns tsquare, which contains Hotelling's T2 statistic for each data point. 返回的tsquare，是表示對每個樣本點Hotelling的T方統計量(我也不很清楚是什麼東東)。 The scores are the data formed by transforming the original data into the space of the principal components. The values of the vector latent are the variance of the columns of SCORE. Hotelling's T2 is a measure of the multivariate distance of each observation from the center of the data set. 所得的分(scores)表示由原資料X轉變到主成分空間所得到的資料。latent向量的值表示SCORE矩陣每列的方差(見說明2)。Hotelling的T方是用來衡量多變數間的距離，這個距離是指樣本觀測值到資料集中心的距離。 When n <= p, SCORE(:,n:p) and latent(n:p) are necessarily zero, and the columns of COEFF(:,n:p) define directions that are orthogonal to X. [...] = princomp(X,'econ') returns only the elements of latent that are not necessarily zero, and the corresponding columns of COEFF and SCORE, that is, when n <= p, only the first n-1. This can be significantly faster when p is much larger than n. 當維數p超過樣本個數n的時候，用[...] = princomp(X,'econ')來計算，這樣會顯著提高計算速度 Examples--舉例 (上面說了那麼多廢話，看了還不一定懂，還不如舉例容易理解，下面樣本資料集為ingredients，這個資料集是matlab自帶的) Compute principal components for the ingredients data in the Hald data set, and the variance accounted for by each component. load hald; %載入matlab內部資料 [pc,score,latent,tsquare] = princomp(ingredients); %呼叫pca分析函式 ingredients,score,pc,latent,tsquare %顯示得到的結果 ingredients = 7 26 6 60 1 29 15 52 11 56 8 20 11 31 8 47 7 52 6 33 11 55 9 22 3 71 17 6 1 31 22 44 2 54 18 22 21 47 4 26 1 40 23 34 11 66 9 12 10 68 8 12 score = 36.8218 -6.8709 -4.5909 0.3967 29.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 pc = -0.0678 -0.6460 0.5673 0.5062 -0.6785 -0.0200 -0.5440 0.4933 0.0290 0.7553 0.4036 0.5156 0.7309 -0.1085 -0.4684 0.4844 latent = 517.7969 67.4964 12.4054 0.2372 tsquare = 5.6803 3.0758 6.0002 2.6198 3.3681 0.5668 3.4818 3.9794 2.6086 7.4818 4.1830 2.2327 2.7216 %下面我們來做一個驗證 %下面為計算ingredients協方差矩陣： cov_ingredients=cov(ingredients) cov_ingredients = 34.6026 20.9231 -31.0513 -24.1667 20.9231 242.1410 -13.8782 -253.4167 -31.0513 -13.8782 41.0256 3.1667 -24.1667 -253.4167 3.1667 280.1667 %下面為計算ingredients所對應的協方差矩陣(也就是cov_ingredients矩陣)的特徵值和特徵 %向量，下面的矩陣V為特徵向量，D為特徵值(對比上面的latent)組成的對角線矩陣 [V,D] = eig(cov_ingredients) V = 0.5062 0.5673 0.6460 -0.0678 0.4933 -0.5440 0.0200 -0.6785 0.5156 0.4036 -0.7553 0.0290 0.4844 -0.4684 0.1085 0.7309 D = 0.2372 0 0 0 0 12.4054 0 0 0 0 67.4964 0 0 0 0 517.7969 %說明1：對比一下矩陣V和矩陣pc，現在很容易明白為什麼COEFF是按列遞減順序排列的 % 了！(V中第三列與pc中倒數第三列差個負號，學過線性代數的人都知道這沒問題) %下面再驗證一下說明2 diag(cov(score)) ans = 517.7969 67.4964 12.4054 0.2372 %說明2：以上結果顯示latent確實表示SCORE矩陣每列的方差，517.7969表示第一列方差 下面做圖表示結果： 上面說了半天還沒有達到我們終極想要的，其實我們要的是由函式[pc,score,latent,tsquare] = princomp(ingredients)它所產生的pc和latent。由latent可以算出降維後的空間所能表示原空間的程度，只要這個累積的值大於95%就行了。 The following command and plot show that two components account for 98% of the variance: cumsum(latent)./sum(latent) ans = 0.86597 0.97886 0.9996 1 %由以上ans值可以看出前兩個主成分就能表示原空間的97.886%,所以取pc中的前兩列可 %做主成分變換矩陣tranMatrix = pc(:,1:2)。則從原來的4維空間降到2維空間。對任意一個 %原空間樣本,例如a=(7 ,26 ,6 ,60)變到低維空間的表示式為a1 = a*tranMatrix。(當然你也可 %以取pc中的前三列，由原來的4維空間變到3維空間) biplot(pc(:,1:2),'Scores',score(:,1:2),'VarLabels',... {'X1' 'X2' 'X3' 'X4'})這位博主（watkinsong）的理解： 上面這個matlab函式的說明呢，只是引用百度百科，也可以看看matlab的函式說明，但是多少還是有點難懂。 我把我的理解簡單的說說。 [COEFF, SCORE, LATENT, TSQUARED] = PRINCOMP(X)
上面這個函式，coeff矩陣是返回的轉換矩陣，也就是把樣本轉換到新的空間中的準換矩陣，這個準換矩陣式比較大的，比如你的降維矩陣式30*100000，那麼這個準換矩陣一般都是10000*29的維數。 score是原來的樣本矩陣在新的座標系中的表示，也就是原來的樣本乘上轉換矩陣，但是還不是直接乘，要減去一個樣本的均值。將原來的資料轉換到新的樣本空間中的演算法是這樣實現的： x0 = bsxfun(@minus,x,mean(x,1)); score = x0 * coeff; 然後就會得到和[COEFF, SCORE, LATENT, TSQUARED] = PRINCOMP(X) 輸出一樣的score資料。 同時這個也是原來的樣本矩陣降維後的結果，如果使用降維後的資料就使用這個資料。一般情況下，如果你的每個樣本的特徵維數遠遠大於樣本數，比如30*1000000的維數,princomp要加上'econ', 就是princomp(x,'econ')這樣使用，可以很大程度的加快計算速度，而且不會記憶體溢位，否則會經常報記憶體溢位。 [...] = PRINCOMP(X,'econ') returns only the elements of LATENT that are
not necessarily zero, i.e., when N <= P, only the first N-1, and the
corresponding columns of COEFF and SCORE. This can be significantly
faster when P >> N.
latent是返回的按降序排列的特徵值，根據這個你可以手動的選擇降維以後的資料要選擇前多少列。
cumsum(latent)./sum(latent)，通過這樣計算特徵值的累計貢獻率，一般來說都選擇前95%的特徵值對應的特徵向量，還是原來的矩陣30*1000000，如果你計算得到前25個特徵值的累計貢獻率已經超過99.9%，那麼就完全可以只要降維後的資料的前25列。 tsquared是個什麼東西我也不知道。。。不過貌似很少有人能用到，網路上也沒有神馬資料，各位如果需要用的再查閱吧，一般情況下也用不到。 如果你需要對測試樣本降維，一般情況下，使用matlab自帶的方式，肯定需要對測試樣本減去一個訓練樣本均值，因為你在給訓練樣本降維的時候減去了均值，所以測試樣本也要減去均值，然後乘以coeff這個矩陣，就獲得了測試樣本降維後的資料。比如說你的測試樣本是1*1000000，那麼乘上一個1000000*29的降維矩陣，就獲得了1*29的降維後的測試樣本的降維資料。 princomp(x)使用的行表示一個樣本，每行的所有的列資料都是這個樣本的特徵值。降維以後比如是30*29，那麼每一行就是降維以後的資料。每個樣本有29個特徵值。
關於網路上的一些解釋個人理解（僅供大家參考理解）
1.
原文地址：http://www.cnblogs.com/sunwufan/archive/2011/08/31/2159952.html 原文：最近看了些主成分分析，混跡Matlab論壇，翻了n多帖子，對princomp函式有了些瞭解。
在此只講一些個人理解，並沒有用術語，只求通俗。
貢獻率：每一維資料對於區分整個資料的貢獻，貢獻率最大的顯然是主成分，第二大的是次主成分......
[coef,score,latent,t2] = princomp(x);（個人觀點）：
x：為要輸入的n維原始資料。帶入這個matlab自帶函式，將會生成新的n維加工後的資料（即score）。此資料與之前的n維原始資料一一對應。
score：生成的n維加工後的資料存在score裡。它是對原始資料進行的分析，進而在新的座標系下獲得的資料。他將這n維資料按貢獻率由大到小排列。（即在改變座標系的情況下，又對n維資料排序）
latent：是一維列向量，每一個數據是對應score裡相應維的貢獻率，因為資料有n維所以列向量有n個數據。由大到小排列（因為score也是按貢獻率由大到小排列）。
coef：是係數矩陣。通過cofe可以知道x是怎樣轉換成score的。
則模型為從原始資料出發：score= bsxfun(@minus,x,mean(x,1))*coef;(作用：可以把測試資料通過此方法轉變為新的座標系)
逆變換：x= bsxfun(@plus,score*inv(coef),mean(x,1))
例子：
View
Code
%% %清屏 clear %% %初始化資料 a=[-14.8271317103068,-3.00108550936016,1.52090778549498,3.95534842970601;-16.2288612441648,-2.80187433749996,-0.410815700402130,1.47546694457079;-15.1242838039605,-2.59871263957451,-0.359965674446737,1.34583763509479;-15.7031424565913,-2.53005662064257,0.255003254103276,-0.179334985754377;-17.7892158910100,-3.32842422986555,0.255791146332054,1.65118282449042;-17.8126324036279,-4.09719527953407,-0.879821957489877,-0.196675865428539;-14.9958877514765,-3.90753364293621,-0.418298866141441,-0.278063876667954;-15.5246706309866,-2.08905845264568,-1.16425848541704,-1.16976057326753;]; x=a; %% %呼叫princomp函式 [coef,score,latent,t2] = princomp(x); score %測試score是否和score_test一樣 score_test=bsxfun(@minus,x,mean(x,1))*coef; score_test latent=100*latent/sum(latent)%將latent總和統一為100，便於觀察貢獻率 pareto(latent);%呼叫matla畫圖
上圖是通過自帶函式繪製，當貢獻率累加至95%，以後的維數會不在顯示，最多隻顯示10維。
下面用自己編寫的表示：
之前的錯誤認識：
1.認為主成分分析中latent顯示的貢獻值是原始資料的，其實是加工後的資料的。解釋:對原始資料既然選擇PCA方法，那麼計算機認為原始資料每維之間可能存在關聯，你想去掉關聯、降低維數。所以採用這種方法的。所以計算機並不關心原始資料的貢獻值，因為你不會去用了，用的是加工後的資料（這也是為什麼當把輸入資料每一維的順序改變後，score、latent不受影響的原因）。
2.認為PCA分析後自動降維，不對。PCA後會有貢獻值，是輸入者根據自己想要的貢獻值進行維數的改變，進而生成資料。（一般大家會取貢獻值在85%以上，要求高一點95%）。
3.PCA分析，只根據輸入資料的特徵進行主成分分析，與輸出有多少型別，每個資料對應哪個型別無關。如果樣本已經分好型別，那PCA後勢必對結果的準確性有一定影響，我認為對於此類資料的PCA，就是在降維與準確性間找一個平衡點的問題，讓資料即不會維數多而使運算複雜，又有較高的解析度。
我的個人見解：這篇文章中的解釋挺靠譜的，可以用來參考。第二點其實matlab的輸出結果score這個資料已經是降維後的資料，不過大家可以根據自己的需要取前多少列的資料。
2。
原文地址：http://www.ilovematlab.cn/thread-54600-1-1.html 部分原文： 回覆 8# 5342245 的帖子 設原始資料為X，先不做任何預處理。[coef,score,latent,t2] = princomp(X);
則那些引數的底層演算法大體過程如下：
x0 = bsxfun(@minus,X,mean(X,1)); %x0為將X去均值後的資料。
[coef,ignore] = eig(x0'*x0); 這就是coef的由來。 【當然最終的還有排序什麼亂七八糟的。。】
scroe = x0*coef % 這就是score的由來，就是一個簡單的線性變換，將原來的X的座標轉換到主成分空間中的座標。僅此而已
則模型為從原始資料出發：
score = bsxfun(@minus,X,mean(X,1))*coef;
逆變換：
X = bsxfun(@plus,score*inv(coef),mean(X,1))
以上這些你可以自己驗證，看是否正確。
關於你的第三問。對於每一個主成分，就看coef的相應的列就能知道原始的變數那個對該主成分貢獻大了啊。。
上面是沒有預處理的。如果加了可逆的預處理。則原始資料亦可從預處理後的資料表示出。進而 bla bla....
===============這回夠通俗易懂吧。。O(∩_∩)O
PS：pca演算法流程，你熟悉嗎？只要知道那個演算法過程。這些都不難理解啊。。
建議您看看書把pca演算法流程再過一遍。。否則別人再怎麼說也沒用。。。
我的個人見解： 這裡我想說的是，再對測試樣本進行降維的時候，一定要減去訓練樣本的均值，使用訓練樣本得到的轉換矩陣，保證訓練樣本和測試樣本轉換到相同的樣本空間中，這樣才有意思。大家有時間可以去看看英文的資料，說的都比較詳細。再用測試樣本減去均值以後，就可以進行轉換了。 很多同學可能在開始的時候和我一樣，都是不知道如果對測試樣本進行降維，很多人就選擇了還是使用princomp這個函式處理測試樣本，那麼這樣測試樣本被對映到一個新的空間中，和原來的訓練樣本完全不是在一個空間，一點意義都沒有，還是要使用測試樣本減去均值，然後乘上訓練樣本降維的時候獲得降維矩陣，轉換到相同的空間中。 基本的對pca的認識就都說完了，比較亂，沒有條理，不過如果認真看下來的話，應該還是可以理解的。目前網上沒有關於pca的綜合的介紹個注意事項，說以我就把我的經驗和大家分享一下，還望文明轉載，轉載宣告出處。我也沒有對pca進行詳細的學習，肯定有不正確的地方，還請大家多多指教，共同探討。