用opencv實現的PCA演算法，非API呼叫

理論參考文獻：但此文沒有程式碼實現。這裡自己實現一下，讓理解更為深刻

問題：如果在IR中我們建立的文件-詞項矩陣中，有兩個詞項為“learn”和“study”，在傳統的向量空間模型中，覺得兩者獨立。

然而從語義的角度來講。兩者是相似的。並且兩者出現頻率也類似，是不是能夠合成為一個特徵呢？

       《模型選擇和規則化》談到的特徵選擇的問題。就是要剔除的特徵主要是和類標籤無關的特徵。比方“學生的名字”就和他的“成績”無關，使用的是互資訊的方法。

       而這裡的特徵非常多是和類標籤有關的，但裡面存在噪聲或者冗餘。

在這樣的情況下，須要一種特徵降維的方法來降低特徵數。降低噪音和冗餘，降低過度擬合的可能性。

        PCA的思想是將n維特徵對映到k維上（k<n）。這k維是全新的正交特徵。

這k維特徵稱為主元，是又一次構造出來的k維特徵，而不是簡單地從n維特徵中去除其餘n-k維特徵。

         ＰＣＡ計算過程：

  如果我們得到的2維資料例如以下：

     

  行代表了例子。列代表特徵，這裡有10個例子，每一個例子兩個特徵。能夠這樣覺得。有10篇文件，x是10篇文件中“learn”出現的TF-IDF。y是10篇文件中“study”出現的TF-IDF。

第一步 分別求x和y的平均值，然後對於全部的例子。都減去相應的均值。這裡x的均值是1.81。y的均值是1.91，那麼一個例子減去均值後即為（0.69,0.49），得到

     

      第二步 ，求特徵協方差矩陣。假設資料是3維。那麼協方差矩陣是

     

     這裡僅僅有x和y，求解得

     

     對角線上各自是x和y的方差，非對角線上是協方差。

協方差是衡量兩個變數同一時候變化的變化程度。協方差大於0表示x和y若一個增。還有一個也增；小於0表示一個增，一個減。假設ｘ和ｙ是統計獨立的。那麼二者之間的協方差就是０；可是協方差是０，並不能說明ｘ和ｙ是獨立的。協方差絕對值越大，兩者對彼此的影響越大。反之越小。協方差是沒有單位的量，因此，假設相同的兩個變數所採用的量綱發生變化，它們的協方差也會產生樹枝上的變化。

    第三步 。求協方差的特徵值和特徵向量。得到

     

     上面是兩個特徵值，以下是相應的特徵向量，特徵值0.0490833989相應特徵向量為 。這裡的特徵向量都歸一化為單位向量。

     第四步 ，將特徵值依照從大到小的順序排序。選擇當中最大的k個。然後將其相應的k個特徵向量分別作為列向量組成特徵向量矩陣。

     這裡特徵值僅僅有兩個，我們選擇當中最大的那個，這裡是1.28402771，相應的特徵向量是 。

      第五步 ，將樣本點投影到選取的特徵向量上。如果例子數為m，特徵數為n，減去均值後的樣本矩陣為DataAdjust(m*n)。協方差矩陣是n*n。選取的k個特徵向量組成的矩陣為EigenVectors(n*k)。那麼投影后的資料FinalData為

     

     這裡是

     FinalData(10*1) = DataAdjust(10*2矩陣)×特徵向量

     得到結果是

     

     這樣，就將原始例子的n維特徵變成了k維，這k維就是原始特徵在k維上的投影。

詳細程式碼例如以下：

#include<cv.h>
#include "opencv2/core/core.hpp"
#include "opencv2/contrib/contrib.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/objdetect/objdetect.hpp"
using namespace std;
using namespace cv;

void printMat(CvMat* matric)
{
  int row = matric->rows,col = matric->cols,i,j;
  for(i = 0;i < row;i++)
  {
    for(j = 0;j < col;j++)
    {
      float num = cvGet2D(matric,i,j).val[0];
      cout << num << " ";
    }
    cout << endl;
  }
}
int main()
{
  printf("-------輸入原矩陣的行和列-----------\n");
  int row,col,i,j;
  cin >> row >> col;
  float* data = new float [row*col];
  printf("-------輸入原矩陣的元素-------------\n");
  for(i = 0;i < col;i++)
  {
    for(j = 0;j < row;j++)cin >> data[j * col + i]; 
  }
  CvMat* inputs = cvCreateMat(row,col,CV_32FC1);
  cvSetData(inputs,data,inputs->step);
  printf("------輸入的矩陣為------------------\n");
  printMat(inputs);
  CvMat* means = cvCreateMat(1,col,CV_32FC1);//均值
  CvMat* EigenValue = cvCreateMat(col,1,CV_32FC1);//特徵值
  CvMat* EigenVector = cvCreateMat(col,col,CV_32FC1);//特徵向量


//	eigen(inputs,EigenValue,EigenVector);

  cvCalcPCA(inputs,means,EigenValue,EigenVector,CV_PCA_DATA_AS_ROW);//以行作為特徵維數計算PCA的引數
  printf("-------------均值為-----------------\n");
  printMat(means);
  printf("-------------特徵值為---------------\n");
  printMat(EigenValue);
  printf("-------------特徵向量為-------------\n");
  printMat(EigenVector);
  printf("-----輸入PCA的降維繫數。按特徵值總和的比例確定PCA的維度，大小在0~1之間------\n");
  float rate;
  cin >> rate;
  float EigenValueSum = cvSum(EigenValue).val[0],curSum = 0.0;
  for(i = 0;i < col;i++)
  {
    curSum += cvGet2D(EigenValue,i,0).val[0];//cvCalcPCA求出的特徵值和特徵向量都是排好序的。能夠直接使用
    if(curSum > EigenValueSum * rate)break;
  }
  int postDim = i;
  cout << "----從原來的 " << col <<" 維降到 " << postDim << " 維"<<endl;
  CvMat* postEigenVector = cvCreateMat(col,postDim,CV_32FC1);//取前k個特徵向量組成的投影矩陣
  CvMat sub;
  cvGetSubRect(EigenVector,postEigenVector,cvRect(0,0,postDim,col));
  for(i = 0;i < row;i++)
  {
    cvGetSubRect(inputs,⊂,cvRect(0,i,col,1));
    cvSub(⊂,means,⊂);
  }
  printf("------------投影矩陣為------------\n");
  printMat(postEigenVector);
  printf("-----------原矩陣單位化後為----------\n");
  printMat(inputs);
  CvMat* output = cvCreateMat(row,postDim,CV_32FC1);
  cvmMul(inputs,postEigenVector,output);//cvMul 是多維向量的點乘,cvmMul 是矩陣乘
  printf("-----------PCA處理後的結果為---------\n");
  printMat(output);

  delete[] data;
  cvReleaseMat(&inputs);
  cvReleaseMat(&means);
  cvReleaseMat(&EigenValue);
  cvReleaseMat(&EigenVector);
  cvReleaseMat(&postEigenVector);
  cvReleaseMat(&output);
  return 0;
}

例子測試：

fangjian@fangjian:~/study/code$ ./opencv 

-------輸入原矩陣的行和列-----------

10 2

-------輸入原矩陣的元素-------------

2.5 0.5 2.2 1.9 3.1 2.3 2 1 1.5 1.1 2.4 0.7 2.9 2.2 3.0 2.7 1.6 1.1 1.6 0.9

------輸入的矩陣為------------------

2.5 2.4 

0.5 0.7 

2.2 2.9 

1.9 2.2 

3.1 3 

2.3 2.7 

2 1.6 

1 1.1 

1.5 1.6 

1.1 0.9 

-------------均值為-----------------

1.81 1.91 

-------------特徵值為---------------

1.15562 

0.044175 

-------------特徵向量為-------------

0.677873 0.735179 

0.735179 -0.677873 

-----輸入PCA的降維繫數，按特徵值總和的比例確定PCA的維度。大小在0~1之間------

0.98

----從原來的 2 維降到 1 維

------------投影矩陣為------------

0.677873 

0.735179 

-----------原矩陣單位化後為----------

0.69 0.49 

-1.31 -1.21 

0.39 0.99 

0.0899999 0.29 

1.29 1.09 

0.49 0.79 

0.19 -0.31 

-0.81 -0.81 

-0.31 -0.31 

-0.71 -1.01 

-----------PCA處理後的結果為---------

0.82797 

-1.77758 

0.992197 

0.27421 

1.6758 

0.912949 

-0.0991095 

-1.14457 

-0.438046 

-1.22382 

注意，輸入的降維比例要足夠大，至少保證有1維 ，這裡使用的是0.98

上面的資料能夠覺得是learn和study特徵融合為一個新的特徵叫做LS特徵，該特徵基本上代表了這兩個特徵。

  上述過程有個圖描寫敘述：

     

     正號表示預處理後的樣本點，斜著的兩條線就各自是正交的特徵向量（因為協方差矩陣是對稱的。因此其特徵向量正交）。最後一步的矩陣乘法就是將原始樣本點分別往特徵向量相應的軸上做投影。

     假設取的k=2。那麼結果是

     

     這就是經過PCA處理後的樣本資料，水平軸（上面舉例為LS特徵）基本上能夠代表所有樣本點。

整個過程看起來就像將座標系做了旋轉，當然二維能夠圖形化表示，高維就不行了。

上面的假設k=1，那麼僅僅會留下這裡的水平軸，軸上是所有點在該軸的投影。

     這樣PCA的過程基本結束。

在第一步減均值之後，事實上應該另一步對特徵做方差歸一化。比方一個特徵是汽車速度（0到100）。一個是汽車的座位數（2到6），顯然第二個的方差比第一個小。因此，假設樣本特徵中存在這樣的情況，那麼在第一步之後，求每一個特徵的標準差  。然後對每一個例子在該特徵下的資料除以  。

     歸納一下，使用我們之前熟悉的表示方法，在求協方差之前的步驟是：

     

     當中  是例子，共m個。每一個例子n個特徵，也就是說  是n維向量。  是第i個例子的第j個特徵。  是例子均值。  是第j個特徵的標準差。

     整個PCA過程貌似及其簡單。就是求協方差的特徵值和特徵向量，然後做資料轉換。可是有沒有認為非常奇妙。為什麼求協方差的特徵向量就是最理想的k維向量？其背後隱藏的意義是什麼？整個PCA的意義是什麼？

 PCA理論基礎

     要解釋為什麼協方差矩陣的特徵向量就是k維理想特徵。我看到的有三個理論：各自是最慷慨差理論、最小錯誤理論和座標軸相關度理論。這裡簡單探討前兩種，最後一種在討論PCA意義時簡單概述。

 最慷慨差理論

     在訊號處理中覺得訊號具有較大的方差。噪聲有較小的方差。訊雜比就是訊號與噪聲的方差比。越大越好。如前面的圖，樣本在橫軸上的投影方差較大，在縱軸上的投影方差較小，那麼覺得縱軸上的投影是由噪聲引起的。

因此我們覺得，最好的k維特徵是將n維樣本點轉換為k維後。每一維上的樣本方差都非常大。

     比方下圖有5個樣本點：（已經做過預處理，均值為0。特徵方差歸一）

     

     以下將樣本投影到某一維上，這裡用一條過原點的直線表示（前處理的過程實質是將原點移到樣本點的中心點）。

     

     如果我們選擇兩條不同的直線做投影，那麼左右兩條中哪個好呢？依據我們之前的方差最大化理論。左邊的好。由於投影后的樣本點之間方差最大。

     這裡先解釋一下投影的概念：

     

     紅色點表示例子  。藍色點表示  在u上的投影，u是直線的斜率也是直線的方向向量，並且是單位向量。藍色點是  在u上的投影點。離原點的距離是 （即  或者  ）因為這些樣本點（例子）的每一維特徵均值都為0。因此投影到u上的樣本點（僅僅有一個到原點的距離值）的均值仍然是0。

     回到上面左右圖中的左圖，我們要求的是最佳的u。使得投影后的樣本點方差最大。

     因為投影后均值為0，因此方差為：

     

     中間那部分非常熟悉啊，不就是樣本特徵的協方差矩陣麼（  的均值為0，一般協方差矩陣都除以m-1，這裡用m）。

     用  來表示  ，  表示  ，那麼上式寫作

        

     因為u是單位向量，即  ，上式兩邊都左乘u得，

     即 

     We got it！

  就是  的特徵值，u是特徵向量。

最佳的投影直線是特徵值  最大時相應的特徵向量。其次是  第二大相應的特徵向量，依次類推。

     因此，我們僅僅須要對協方差矩陣進行特徵值分解。得到的前k大特徵值相應的特徵向量就是最佳的k維新特徵，並且這k維新特徵是正交的。得到前k個u以後，例子 通過下面變換能夠得到新的樣本。

     

     當中的第j維就是  在  上的投影。

     通過選取最大的k個u，使得方差較小的特徵（如噪聲）被丟棄。

最小平方誤差理論：

     如果有這種二維樣本點（紅色點）。回想我們前面探討的是求一條直線。使得樣本點投影到直線上的點的方差最大。本質是求直線，那麼度量直線求的好不好。不只唯獨方差最大化的方法。再回想我們最開始學習的線性迴歸等，目的也是求一個線性函式使得直線可以最佳擬合樣本點。那麼我們能不能覺得最佳的直線就是迴歸後的直線呢？迴歸時我們的最小二乘法度量的是樣本點到直線的座標軸距離。

比方這個問題中，特徵是x。類標籤是y。迴歸時最小二乘法度量的是距離d。如果使用迴歸方法來度量最佳直線。那麼就是直接在原始樣本上做迴歸了，跟特徵選擇就沒什麼關係了。

     因此，我們打算選用第二種評價直線好壞的方法，使用點到直線的距離d’來度量。

     如今有n個樣本點  。每一個樣本點為m維（這節內容中使用的符號與上面的不太一致。須要又一次理解符號的意義）。將樣本點  在直線上的投影記為 ，那麼我們就是要最小化

     

     這個公式稱作最小平方誤差（Least Squared Error）。

     而確定一條直線，一般僅僅須要確定一個點，而且確定方向就可以。

      第一步確定點：

     如果要在空間中找一點  來代表這n個樣本點，“代表”這個詞不是量化的，因此要量化的話。我們就是要找一個m維的點  ，使得

     

     最小。

當中  是平方錯誤評價函式（squared-error criterion function），如果 m 為n個樣本點的均值：

     

     那麼平方錯誤能夠寫作：

     

     後項與  無關。看做常量。而  ，因此最小化  時，

        

       是樣本點均值。

      第二步確定方向：

     我們從  拉出要求的直線（這條直線要過點 m ），如果直線的方向是單位向量 e 。

那麼直線上隨意一點，比方  就能夠用點 m 和 e 來表示

        

     當中  是  到點 m 的距離。

     我們又一次定義最小平方誤差：

     

     這裡的k僅僅是相當於 i 。  就是最小平方誤差函式，當中的未知引數是 和 e 。

     實際上是求  的最小值。首先將上式展開：

     

     我們首先固定 e ，將其看做是常量。  ，然後對  進行求導，得

     

     這個結果意思是說，假設知道了 e 。那麼將  與 e 做內積，就能夠知道了 在 e 上的投影離 m 的長度距離。只是這個結果不用求都知道。

     然後是固定  ，對 e 求偏導數。我們先將公式（8）代入  。得  

     

     當中  與協方差矩陣類似。僅僅是缺少個分母n-1，我們稱之為 雜湊矩陣 （scatter matrix）。

     然後能夠對 e 求偏導數。可是 e 須要首先滿足  。引入拉格朗日乘子 。來使  最大（  最小）。令

     

     求偏導

     

     這裡存在對向量求導數的技巧。方法這裡不多做介紹。

能夠去看一些關於矩陣微積分的資料。這裡求導時能夠將  看作是  ，將  看做是  。

     導數等於0時。得

     

     兩邊除以n-1就變成了。對協方差矩陣求特徵值向量了。

     從不同的思路出發，最後得到同一個結果。對協方差矩陣求特徵向量，求得後特徵向量上就成為了新的座標，例如以下圖：

     

     這時候點都聚集在新的座標軸周圍，由於我們使用的最小平方誤差的意義就在此。

PCA理論意義：

   PCA將n個特徵降維到k個，能夠用來進行資料壓縮，假設100維的向量最後能夠用10維來表示，那麼壓縮率為90%。相同影象處理領域的KL變換使用PCA做影象壓縮。但PCA要保證降維後。還要保證資料的特性損失最小。

再看回想一下PCA的效果。經過PCA處理後。二維資料投影到一維上能夠有下面幾種情況：

     

     我們覺得左圖好，一方面是投影后方差最大，一方面是點到直線的距離平方和最小，並且直線過樣本點的中心點。為什麼右邊的投影效果比較差？直覺是由於座標軸之間相關，以至於去掉一個座標軸。就會使得座標點無法被單獨一個座標軸確定。

     PCA得到的k個座標軸實際上是k個特徵向量，因為協方差矩陣對稱。因此k個特徵向量正交。

看以下的計算過程。

     如果我們還是用 來表示例子。m個例子，n個特徵。

特徵向量為 e 。  表示第i個特徵向量的第1維。

那麼原始樣本特徵方程能夠用以下式子來表示：

     前面兩個矩陣乘積就是協方差矩陣  （除以m後）。原始的樣本矩陣A是第二個矩陣m*n。

     

     上式能夠簡寫為 

     我們最後得到的投影結果是  。E是k個特徵向量組成的矩陣，展開例如以下：

     

     得到的新的例子矩陣就是m個例子到k個特徵向量的投影，也是這k個特徵向量的線性組合。e之間是正交的。

從矩陣乘法中能夠看出，PCA所做的變換是將原始樣本點（n維），投影到k個正交的座標系中去，丟棄其它維度的資訊。舉個例子，如果宇宙是n維的（霍金說是11維的），我們得到銀河系中每一個星星的座標（相對於銀河系中心的n維向量）。然而我們想用二維座標去逼近這些樣本點。如果算出來的協方差矩陣的特徵向量各自是圖中的水平和豎直方向。那麼我們建議以銀河系中心為原點的x和y座標軸。全部的星星都投影到x和y上。得到以下的圖片。然而我們丟棄了每一個星星離我們的遠近距離等資訊。

     

總結與討論：

PCA技術的一大優點是對資料進行降維的處理。

我們能夠對新求出的“主元”向量的重要性進行排序。依據須要取前面最重要的部分，將後面的維數省去，能夠達到降維從而簡化模型或是對資料進行壓縮的效果。

同一時候最大程度的保持了原有資料的資訊。

     PCA技術的一個非常大的長處是，它是全然無引數限制的。

在PCA的計算過程中全然不須要人為的設定引數或是依據不論什麼經驗模型對計算進行干預，最後的結果僅僅與資料相關。與使用者是獨立的。 

可是。這一點同一時候也能夠看作是缺點。

假設使用者對觀測物件有一定的先驗知識，掌握了資料的一些特徵，卻無法通過引數化等方法對處理過程進行干預。可能會得不到預期的效果，效率也不高。

     

     

圖表 4：黑色點表示取樣資料，排列成轉盤的形狀。 

     easy想象，該資料的主元是  或是旋轉角  。

     如圖表 4中的樣例，PCA找出的主元將是  。可是這顯然不是最優和最簡化的主元。  之間存在著非線性的關係。

依據先驗的知識可知旋轉角  是最優的主元（類比極座標）。則在這樣的情況下，PCA就會失效。

可是，假設增加先驗的知識，對資料進行某種劃歸。就能夠將資料轉化為以  為線性的空間中。這類依據先驗知識對資料預先進行非線性轉換的方法就成為 kernel -PCA，它擴充套件了PCA能夠處理的問題的範圍。又能夠結合一些先驗約束，是比較流行的方法。

     有時資料的分佈並非滿足高斯分佈。如圖表 5所看到的。在非高斯分佈的情況下，PCA方法得出的主元可能並非最優的。在尋找主元時不能將方差作為衡量重要性的標準。要依據資料的分佈情況選擇合適的描寫敘述全然分佈的變數。然後依據概率分散式

     

     來計算兩個向量上資料分佈的相關性。等價的，保持主元間的正交如果。尋找的主元相同要使  。

這一類方法被稱為獨立主元分解(ICA)。

     

     

圖表 5：資料的分佈並不滿足高斯分佈，呈明顯的十字星狀。 

     這樣的情況下。方差最大的方向並非最優主元方向。

     另外PCA還能夠用於預測矩陣中缺失的元素。