一、簡介
PCA（Principal Components Analysis）即主成分分析，是圖像處理中經常用到的降維方法，大家知道，我們在處理有關數字圖像處理方面的問題時，比如經常用的圖像的查詢問題，在一個幾萬或者幾百萬甚至更大的數據庫中查詢一幅相近的圖像。這時，我們通常的方法是對圖像庫中的圖片提取響應的特征，如顏色，紋理，sift，surf，vlad等等特征，然后將其保存，建立響應的數據索引，然后對要查詢的圖像提取相應的特征，與數據庫中的圖像特征對比，找出與之最近的圖片。這里，如果我們為了提高查詢的准確率，通常會提取一些較為復雜的特征，如sift，surf等，一幅圖像有很多個這種特征點，每個特征點又有一個相應的描述該特征點的128維的向量，設想如果一幅圖像有300個這種特征點，那么該幅圖像就有300*vector（128維）個，如果我們數據庫中有一百萬張圖片，這個存儲量是相當大的，建立索引也很耗時，如果我們對每個向量進行PCA處理，將其降維為64維，是不是很節約存儲空間啊？對於學習圖像處理的人來說，都知道PCA是降維的，但是，很多人不知道具體的原理，為此，我寫這篇文章，來詳細闡述一下PCA及其具體計算過程：
二、PCA原理
1、原始數據：
為了方便，我們假定數據是二維的，借助網絡上的一組數據，如下：
x=[2.5, 0.5, 2.2, 1.9, 3.1, 2.3, 2, 1,1.5, 1.1]T
y=[2.4, 0.7, 2.9, 2.2, 3.0, 2.7, 1.6, 1.1, 1.6, 0.9]T
2、計算協方差矩陣
什么是協方差矩陣？相信看這篇文章的人都學過數理統計，一些基本的常識都知道，但是，也許你很長時間不看了，都忘差不多了，為了方便大家更好的理解，這里先簡單的回顧一下數理統計的相關知識，當然如果你知道協方差矩陣的求法你可以跳過這里。
（1）協方差矩陣：
首先我們給你一個含有n個樣本的集合，依次給出數理統計中的一些相關概念：
均值： 
標准差： 
方差： 
既然我們都有這么多描述數據之間關系的統計量，為什么我們還要用協方差呢？我們應該注意到，標准差和方差一般是用來描述一維數據的，但現實生活我們常常遇到含有多維數據的數據集，最簡單的大家上學時免不了要統計多個學科的考試成績。面對這樣的數據集，我們當然可以按照每一維獨立的計算其方差，但是通常我們還想了解這幾科成績之間的關系，這時，我們就要用協方差，協方差就是一種用來度量兩個隨機變量關系的統計量，其定義為：
從協方差的定義上我們也可以看出一些顯而易見的性質，如：
需要注意的是，協方差也只能處理二維問題，那維數多了自然就需要計算多個協方差，比如n維的數據集就需要計算CN2【此乃組合數基本公式】個協方差，那自然而然的我們會想到使用矩陣來組織這些數據。給出協方差矩陣的定義：
這個定義還是很容易理解的，我們可以舉一個簡單的三維的例子，假設數據集有三個維度{x,y,z}，則協方差矩陣為
可見，協方差矩陣是一個對稱的矩陣，而且對角線是各個維度上的方差。
（2）協方差矩陣的求法：
協方差矩陣計算的是不同維度之間的協方差，而不是不同樣本之間的。下面我們將在matlab中用一個例子進行詳細說明：
首先，隨機產生一個10*3維的整數矩陣作為樣本集，10為樣本的個數，3為樣本的維數。
- MySample = fix(rand(10,3)*50)
根據公式，計算協方差需要計算均值，那是按行計算均值還是按列呢，我一開始就老是困擾這個問題。前面我們也特別強調了，協方差矩陣是計算不同維度間的協方差，要時刻牢記這一點。樣本矩陣的每行是一個樣本，每列為一個維度，所以我們要按列計算均值。為了描述方便，我們先將三個維度的數據分別賦值：
dim1 = MySample(:,1); dim2 = MySample(:,2); dim3 = MySample(:,3); %計算dim1與dim2，dim1與dim3，dim2與dim3的協方差： sum( (dim1-mean(dim1)) .*(dim2-mean(dim2)) ) / ( size(MySample,1)-1 ) %得到 74.5333 sum( (dim1-mean(dim1)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到 -10.0889 sum( (dim2-mean(dim2)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到 -10***000 %搞清楚了這個后面就容易多了，協方差矩陣的對角線就是各個維度上的方差，下面我們依次計算： std(dim1)^2 % 得到 108.3222 std(dim2)^2 % 得到 260.6222 std(dim3)^2 % 得到 94.1778 %這樣，我們就得到了計算協方差矩陣所需要的所有數據，調用Matlab自帶的cov函數進行驗證： cov(MySample)
可以看到跟我們計算的結果是一樣的，說明我們的計算是正確的。但是通常我們不用這種方法，而是用下面簡化的方法進行計算：
先讓樣本矩陣中心化，即每一維度減去該維度的均值，然后直接用新的到的樣本矩陣乘上它的轉置，然后除以(N-1)即可。其實這種方法也是由前面的公式通道而來，只不過理解起來不是很直觀而已。大家可以自己寫個小的矩陣看一下就明白了。其Matlab代碼實現如下：
1 X = MySample –repmat(mean(MySample),10,1); %中心化樣本矩陣 2 C = (X’*X)./(size(X,1)-1) 3 %為方便對matlab不太明白的人，小小說明一下各個函數，同樣，對matlab有一定基礎的人直接跳過： 4 %B = repmat(A,m,n ) %%將矩陣 A復制 m×n塊，即把 A 作為 B的元素，B由 m×n個 A平鋪而成。B的維數是 [size(A,1)*m, (size(A,2)*n] 5 %B = mean(A)的說明： 6 %如果你有這樣一個矩陣：A = [1 2 3; 3 36; 4 6 8; 4 7 7]; 7 %用mean(A)（默認dim=1）就會求每一列的均值 8 % ans = 9 % 3.0000 4.5000 6.0000 10 % 用mean(A,2)就會求每一行的均值 11 % ans = 12 % 2.0000 13 % 4.0000 14 % 6.0000 15 % 6.0000 16 size(A,n)%% 如果在size函數的輸入參數中再添加一項n，並用1或2為n賦值，則 size將返回矩陣的行數或列數。其中r=size(A,1)該語句返回的是矩陣A的行數， %c=size(A,2)該語句返回的是矩陣A的列數
上面我們簡單說了一下協方差矩陣及其求法，言歸正傳，我們用上面簡化求法，求出樣本的協方差矩陣為：
3、計算協方差矩陣的特征向量和特征值
因為協方差矩陣為方陣，我們可以計算它的特征向量和特征值，如下：
我們可以看到這些矢量都是單位矢量，也就是它們的長度為1，這對PCA來說是很重要的。
4、選擇成分組成模式矢量
求出協方差矩陣的特征值及特征向量之后，按照特征值由大到小進行排列，這將給出成分的重要性級別。現在，如果你喜歡，可以忽略那些重要性很小的成分，當然這會丟失一些信息，但是如果對應的特征值很小，你不會丟失很多信息。如果你已經忽略了一些成分，那么最后的數據集將有更少的維數，精確地說，如果你的原始數據是n維的，你選擇了前p個主要成分，那么你現在的數據將僅有p維。現在我們要做的是組成一個模式矢量，這只是幾個矢量組成的矩陣的一個有意思的名字而已，它由你保持的所有特征矢量構成，每一個特征矢量是這個矩陣的一列。
對於我們的數據集，因為有兩個特征矢量，因此我們有兩個選擇。我們可以用兩個特征矢量組成模式矢量：
我們也可以忽略其中較小特征值的一個特征矢量，從而得到如下模式矢量：
5、得到降維后的數據
其中rowFeatureVector是由模式矢量作為列組成的矩陣的轉置，因此它的行就是原來的模式矢量，而且對應最大特征值的特征矢量在該矩陣的最上一行。rowdataAdjust是每一維數據減去均值后，所組成矩陣的轉置，即數據項目在每一列中，每一行是一維，對我們的樣本來說即是，第一行為x維上數據，第二行為y維上的數據。FinalData是最后得到的數據，數據項目在它的列中，維數沿着行。
這將給我們什么結果呢？這將僅僅給出我們選擇的數據。我們的原始數據有兩個軸（x和y），所以我們的原始數據按這兩個軸分布。我們可以按任何兩個我們喜歡的軸表示我們的數據。如果這些軸是正交的，這種表達將是最有效的，這就是特征矢量總是正交的重要性。我們已經將我們的數據從原來的xy軸表達變換為現在的單個特征矢量表達。
說明：如果要恢復原始數據，只需逆過程計算即可，即：
到此為止，相信你已經掌握了PCA的原理了。
三 . PCA的應用
PCA及其改進算法主要應用的人臉識別領域，是人臉識別的經典算法之一。OpenCv2.4以后的版本實現了三種經典的人臉識別算法，其中就包括PCA。對openCv比較老的版本也可以調用PCA的算法去做，只是稍顯復雜而已，網上有一篇博文如下：
http://www.cognotics.com/opencv/servo_2007_series/part_5/index.html
該代碼運行在openCv2.1之前的版本當中，但是該代碼有個重要的bug就是特征數K被設置為固定的值，而選擇更小的值的時候，代碼會crash。
PCA另外一個主要的用途是作為其他算法的預處理，術語叫做數據的白化。由於PCA具有壓縮數據的作用，所以可以認為經過PCA處理過之后的數據是不相關的，但一般未必是獨立的。實際可用的PCA算法一般不是以解析解的形式給出的，而是在線學習算法。有很多的原因決定了只能使用在線學習算法。在線學習算法主要有基於神經網絡學習的算法和遞歸最小二乘法，相關的文獻如下：
http://wenku.baidu.com/view/c91f31c058f5f61fb73666f8.html
要注意的是openCv的實現不是在線學習算法。
四. PCA 的實現
1 #ifndef _FCE__PCA__H__ 2 #define _FCE__PCA__H__ 3 4 #define HIGH_PRECISON 5 6 #ifdef HIGH_PRECISON 7 #define real float 8 #else 9 #define real double 10 #endif 11 12 13 #ifdef _cplusplus 14 { 15 extern "C" 16 #endif 17 18 19 typedef struct _FCE_PCA{ 20 int count; //the number of sample 21 int n; // the number of features 22 real *covariance; 23 real *mean; 24 real *z; 25 }FCE_PCA; 26 27 28 FCE_PCA *fce_pca_init(int n); 29 30 void fce_pca_push_add(FCE_PCA *pca, real *v); 31 32 int fce_pca_solve_eig(FCE_PCA *pca, real *eigenvector, real *eigenvalue); 33 34 void fce_pca_free(FCE_PCA *pca); 35 36 #ifdef _cplusplus 37 } 38 #endif 39 40 #endif
1 #include "fce_pca.h" 2 3 4 #define FCE_MIN(i,j) (((i) > (j)) ? (j) : (i)) 5 #define FCE_MAX(i,j) (((i) > (j)) ? (i) : (j)) 6 7 FCE_PCA *fce_pca_init(int n){ 8 FCE_PCA *pca; 9 real zero = 0.0; 10 if(n <= 1) 11 return NULL; 12 13 pca = (FCE_PCA* )malloc(sizeof(FCE_PCA)); 14 if (pca == NULL){ 15 return NULL; 16 } 17 18 pca->n = n; 19 pca->z = (real* )malloc(sizeof(*pca->z) * n); 20 if (pca->z == NULL){ 21 free(pca); 22 return NULL; 23 } 24 25 memset(pca->z, zero, sizeof(*pca->z) * n); 26 27 pca->count=0; 28 pca->covariance= (real* )malloc(sizeof(real) * n * n); 29 if (pca->covariance == NULL){ 30 free(pca->z); 31 free(pca); 32 return NULL; 33 } 34 35 memset(pca->covariance, zero, sizeof(real) * n * n); 36 37 pca->mean = (real* )malloc(sizeof(real) * n); 38 if (pca->mean == NULL){ 39 free(pca->covariance); 40 free(pca->z); 41 free(pca); 42 return NULL; 43 } 44 45 memset(pca->mean, zero, sizeof(real) * n); 46 47 return pca; 48 } 49 50 void fce_pca_free(FCE_PCA *pca){ 51 free(pca->covariance); 52 free(pca->mean); 53 free(pca->z); 54 free(pca); 55 } 56 57 void fce_pca_push_add(FCE_PCA *pca, real *v){ 58 int i, j; 59 const int n = pca->n; 60 for(i = 0; i < n; i++){ 61 pca->mean[i] += v[i]; 62 for(j = i; j < n; j++) 63 pca->covariance[j + i * n] += v[i]*v[j]; 64 } 65 pca->count++; 66 } 67 68 int fce_pca_solve_eig(FCE_PCA *pca, real *eigenvector, real *eigenvalue){ 69 int i, j, pass; 70 int k = 0; 71 const int n = pca->n; 72 real *z = pca->z; 73 real zero = 0.0; 74 75 memset(eigenvector, zero, sizeof(real)*n*n); 76 77 for(j = 0; j < n; j++){ 78 pca->mean[j] /= pca->count; 79 eigenvector[j + j * n] = 1.0; 80 for(i = 0; i <= j; i++){ 81 pca->covariance[j + i * n] /= pca->count; 82 pca->covariance[j + i * n] -= pca->mean[i] * pca->mean[j]; 83 pca->covariance[i + j * n] = pca->covariance[j + i * n]; 84 } 85 eigenvalue[j] = pca->covariance[j + j*n]; 86 z[j] = 0; 87 } 88 89 for(pass=0; pass < 50; pass++){ 90 real sum = 0; 91 for(i = 0; i < n; i++) 92 for(j = i+1; j < n; j++) 93 sum += fabs(pca->covariance[j + i * n]); 94 95 if(sum == 0){ 96 for(i = 0; i < n; i++){ 97 real maxvalue = -1; 98 for(j = i; j < n; j++){ 99 if(eigenvalue[j] > maxvalue){ 100 maxvalue = eigenvalue[j]; 101 k= j; 102 } 103 } 104 eigenvalue[k] = eigenvalue[i]; 105 eigenvalue[i] = maxvalue; 106 for(j = 0; j < n; j++){ 107 real tmp = eigenvector[k + j * n]; 108 eigenvector[k + j * n] = eigenvector[i + j * n]; 109 eigenvector[i + j * n] = tmp; 110 } 111 } 112 return pass; 113 } 114 115 for(i = 0; i < n; i++){ 116 for(j = i + 1; j < n; j++){ 117 real covar = pca->covariance[j + i * n]; 118 real t,c,s,tau,theta, h; 119 120 if(pass < 3 && fabs(covar) < sum / (5*n*n)) 121 continue; 122 if(fabs(covar) <= 0.00000000001) 123 continue; 124 if(pass >=3 && fabs((eigenvalue[j]+z[j])/covar) > (1LL<<32) && fabs((eigenvalue[i]+z[i])/covar) > (1LL<<32)){ 125 pca->covariance[j + i * n]=0.0; 126 continue; 127 } 128 129 h = (eigenvalue[j] + z[j]) - (eigenvalue[i] + z[i]); 130 theta = 0.5 * h/covar; 131 t = 1.0 /(fabs(theta) + sqrt(1.0 + theta * theta)); 132 if(theta < 0.0) t = -t; 133 134 c = 1.0 /sqrt(1 + t * t); 135 s = t * c; 136 tau = s /(1.0 + c); 137 z[i] -= t * covar; 138 z[j] += t * covar; 139 140 #define ROTATE(a,i,j,k,l) {\ 141 real g =a[j + i*n];\ 142 real h =a[l + k*n];\ 143 a[j + i*n] = g - s * (h + g * tau);\ 144 a[l + k*n] = h + s * (g - h * tau); } 145 for(k = 0; k < n; k++) { 146 if(k != i && k != j){ 147 ROTATE(pca->covariance,FCE_MIN(k,i),FCE_MAX(k,i),FCE_MIN(k,j),FCE_MAX(k,j)) 148 } 149 ROTATE(eigenvector,k,i,k,j) 150 } 151 pca->covariance[j + i * n]=0.0; 152 } 153 } 154 for (i = 0; i < n; i++) { 155 eigenvalue[i] += z[i]; 156 z[i]=0.0; 157 } 158 } 159 160 return 0; 161 }