資料降維（PCA、KPCA、PPCA）及C++實作

1、何為資料降維

1.1維數災難：往往滿足采樣條件所需的樣本數目巨大、樣本稀疏、距離計算困難。

1.2降維：利用數學變換将原始高維屬性空間轉變為低維“子空間”，即在高維采樣資料中提取能夠表達原始資料的特征。

1.3 降維優點：資料集更易懂、使用和顯示；降低算法計算開銷；去除噪聲。

2、一些降維算法

Principal Component Analysis （PCA）

Linear Discriminant Analysis（LDA）

Locally linear embedding（LLE）

Laplacian Eigenmaps

本文主要針對以下三種算法：

2.1 PCA：PCA算法是一種線性投影技術，利用降維後使資料的方差最大原則保留盡可能多的資訊；

2.2 KPCA：PCA僅考慮了資料的二階統計資訊，而沒有利用高階統計資訊，忽略了資料的非線性相關性，而KPCA，通過非線性變換将資料映射到了高維，在高維空間中進行特征提取，獲得了更好的特征提取性能；

2.3 PPCA：PCA沒有将資料的機率分布考慮，PPCA對PCA做了機率上的解釋，延伸了PCA算法。

總之：PPCA和KPCA都是針對PCA算法的缺陷，做出了不同方向上的改進。

3 PCA、KPCA、PPCA算法步驟

3.1 PCA： 資料在低維線性空間上的正交投影，這個線性空間被稱為主子空間，使得投影資料的方差被最大化。

》将原始資料按列組成n行m列矩陣X；

》将X的每一行（代表一個屬性字段）進行零均值化，即減去這一行的均值；

》求出協方差矩陣C；

》求出協方差矩陣的特征值及對應的特征向量；

》将特征向量按對應特征值大小從上到下按行排列成矩陣，取前k行組成矩陣P；

》Y=PX即為降維到k維後的資料。

3.2 KPCA：通過使用一個非線性核替換線性的方式來對高維資料向低維投影。

》将原始資料按列組成m行n列矩陣X；

》計算核矩陣，先標明高斯徑向核函數中的參數，計算核矩陣K，修正核矩陣得到KL；

》求出協方差矩陣C，運用雅可比疊代方法計算KL特征值與特征向量；

》将特征向量按對應特征值大小從上到下按行排列成矩陣，取前k行組成矩陣；

》通過施密特正交化方法機關正交化特征向量得到P；

》Y=PX即為降維到k維後的資料。

3.3 PPCA：PPCA是一種考慮每個變量機率分布的方法，在确定主元和誤差的機率函數後，通過期望最大（EM）算法建立模型。

》将原始資料按列組成n行m列矩陣；

》對原始訓練樣本資料進行标準中心化處理得到X；

》在隐含變量x 的條件下得到觀測資料的機率分布；

》通過EM 算法獲得機率PCA 的模型參數W(因子矩陣)和方差；

》舍去不滿足因子矩陣與方差特定關系的歸一化資料；

》剩餘滿足條件資料即為降維到k維後的資料。

4、C++實作PCA、KPCA

（環境：Visual Studio 2013）

pca.h
    #include<math.h>    
     #include<fstream>    
     #include<iostream>    
     #include<stdio.h>    
     #include<iomanip>    
     #include<math.h>
     using namespace std;
      
     typedef struct sourcedata      //聲明了一個原始資料類型
     {
         int m;
         int n;
         double **data;
     }SourceData;
     class PCA
     {
     public:
         PCA(int m, int n);                                //m為行數，n為列數    
         SourceData getdata(const  char *file);           //擷取外部資料
         void standarddata(double **a);                  //資料标準化
         double  product(double *a, double *b);         //向量乘積
         void  swap(double &x, double &y);             //資料交換
         double  **matrixproduct(double **a);            //求解協方差矩陣
         void  selectionsort(double *A, double **v);     //特征值排序
         void  zhengjiao(double **v);                   //向量正交化
         int jcb(double **a, double **v, double eps, int jt);         //求解特征值和特征向量
         int selectcharactor(double *A, double getratio, double *B); //提取主分量
         double  **getProject(int t, double **x, double **v);        //計算降維後特征點  
         void   saveProject(const char *projectfile, double **project, int t); //儲存
         ~PCA(){}
     private:
         int  rows;
         int columns;
     };kpca.h
    #include<math.h>
     #include<fstream>
     #include<iostream>
     #include<stdio.h>
     #include<iomanip>
     #include<math.h>
     #include<stdlib.h>
     #include<time.h>
     using namespace std;
      
     class KPCA
     {
     public:
         KPCA(int m, int n);
         SourceData getdata(const char *file); //擷取外部資料
         int randdef(int n1, int n2); //生成n1到n2随機整數
         double getvar(double **testdata, int m, int n, int l, double left, double right);//通過對随機樣本的最大特征提取效率擷取高斯徑向基函數的參數
         double product(double *a, double *b, int size);  //向量乘積
         double kernel(double var, double *x, double *y, int sign);  //核函數定義
         double **getkernelmatrix(double **a, double var, int sign); //擷取核矩陣
         double **modifykernelmatrix(double **K);  //修正核矩陣
         int jcb(double **a, double **v, double eps, int jt);    //求解矩陣的特征值和特征向量    
         void zhengjiao(double **v); //正交化特征向量
         void swap(double &x, double &y); //交換元素
         void selectionsort(double *A, double **v);  //特征值和特征向量選擇排序
         void saveeigenvectors(double A[], double **v, const char *vectorfile);//儲存特征值和特征向量
         int selectcharactor(double *A, double getratio, double *B);  //提取特征
         double  **getProject(int t, double **x, double **v);  //獲得投影
         void saveProject(const char *projectfile, double **project, int t);  //儲存投影
         ~KPCA(){}
     private:
         int  rows;
         int columns;
     };pca.cpp與kpca.cpp由于篇幅問題未分享，
main.cpp
    #include"pca.h"    
     #include"kpca.h"
     void main()
     {
         //pca
      
         cout << "-----------------------pca------------------------" << endl;
         int i, j, t;                                                        //i,j循環用；t降維後維數
         int m, n;                                                           //m行n列
         double **x, **c, **v, **Project;                                      
         double *A, *B;                                                         //A特征值B貢獻率
         sourcedata pp;
         double eps = 0.000001;                                              //雅克比方法的終止精度            
         double getratio = 0.9;                                              //特征值的提取率      
         const  char *File = "test1.txt";                                    //原始資料檔案名稱    
         const char *projectfile = "pcaproject.txt";                         //處理後的資料檔案名稱    
         PCA pca(2, 3);                                                      //聲明一個臨時對象調用成員函數來擷取資料    
         pp = pca.getdata(File);                                             //擷取外部資料    
         x = pp.data;
         m = pp.m;
         n = pp.n;
      
         cout << "資料的行數為" << m << "，資料的列數為 " << n << endl;
         A = new double[n];                                           //存放特征值
         B = new double[n];                                          //存放特值貢獻率
         v = new double*[n];                                          //存放特征向量
         for (i = 0; i < n; i++)
             v[i] = new double[n];
      
      
         PCA  testpca(m, n);                                        //聲明一個對象并初始化           
         testpca.standarddata(x);                                   //對資料進行标準化處理   X是原始資料        
         c = testpca.matrixproduct(x);                              //    求協方差矩陣  
         i = testpca.jcb(c, v, eps, 100);                           //    求特征值和特征向量    
         for (int k = 0; k < n; k++)
             A[k] = c[k][k];                                        //存特征值    
         testpca.zhengjiao(v);                                     //正交化特征向量    
         testpca.selectionsort(A, v);                              //特征值和特征向量排序        
         t = testpca.selectcharactor(A, getratio, B);              //提取特征值 t為降維後維數   
      
      
         cout << "PCA降維後的維數：" << t << endl;
         cout << "排序後提取的特征值及對應的特征向量" << endl;
         for (i = 0; i <= t - 1; i++)                              //輸出特征值
             printf("%13.7e  ", A[i]);
         printf("\n\n");
      
         for (i = 0; i < n; i++)                                    //輸出特征向量
         {
             for (j = 0; j < t; j++)
                 printf("%13.7e  ", v[i][j]);
             printf("\n");
      
         }
         cout << "特征值的累計貢獻率為" << endl;
         for (i = 0; i < n; i++)
             cout << B[i] << "  ";
         cout << endl;
         cout << "當提取效率是" << getratio << "時提取了前" << t << "個分量" << endl;     //getratio特征提取率
         if (t >= 1 && t <= n)
             Project = testpca.getProject(t, x, v);                //求降維後特征點       
         else
             cout << "error" << endl;
         testpca.saveProject(projectfile, Project, t);             //儲存特征點到TXT檔案  
      
      
      
         //kpca
      
         cout << endl<< "----------------------kpca------------------------" << endl;
         int a;                                                      //循環用
         int l = 50;                                                 //随機提取樣本的數目    
         const char *File2 = "test2.txt";
         const char*eigenvectors = "eigen.txt";                     //特征值和特征向量存儲檔案名稱
         const char *projectfile2 = "kpcaproject.txt";              //降維後特征點檔案存儲名稱
         SourceData pdata;
         double  gaussparameter;                                    //高斯核參數
         double **K, **KL;                                           //高斯核矩陣k及修正核矩陣    
         KPCA kpca(3, 2);
         pdata = kpca.getdata(File2);                               //擷取外部資料
         x = pdata.data;
         m = pdata.m;
         n = pdata.n;
         A = new double[m];
         B = new double[m];
      
         KPCA  testkpca(m, n);                                        //對象    
         gaussparameter = testkpca.getvar(x, m, n, l, 100, 800);       //求高斯核參數    通過對随機樣本的最大特征提取效率擷取高斯徑向基函數的參數
         cout << "高斯核參數： " << gaussparameter << endl;
         K = testkpca.getkernelmatrix(x, gaussparameter, 1);          //求核矩陣
         KL = testkpca.modifykernelmatrix(K);                        //求修正核矩陣    
         c = new double*[m];                                         //定義c、v二維數組
         for (a = 0; a<m; a++)
             c[a] = new double[m];
         v = new double*[m];
         for (a = 0; a<m; a++)
             v[a] = new double[m];
         for (a = 0; a<m; a++)                                        //修正核矩陣放入c
             for (j = 0; j<m; j++)
                 c[a][j] = KL[a][j];
         a = testkpca.jcb(c, v, eps, 10000);                            //求取特征值和特征向量
         cout << "計算特征值的疊代次數為" << a << endl;
         if (a != -1)
         {
             for (a = 0; a<m; a++)
                 A[a] = c[a][a];                                     //特征值存入A
         }
         else
             cout << "不能求得特征值和特征向量" << endl;
         testkpca.zhengjiao(v);                                      //正交化特征向量
         testkpca.saveeigenvectors(A, v, eigenvectors);
         testkpca.selectionsort(A, v);                             //特征值和特征向量排序    
         t = testkpca.selectcharactor(A, getratio, B);            //提取特征值    
      
         cout << "特征值的累計貢獻率是" << endl;
         for (a = 0; a<m; a++)
             cout << B[a] << "  ";
         cout << endl;
         cout << "當提取效率為" << getratio << "時提取了前" << t << "個分量" << endl;
         if (t >= 1 && t <= m)
             Project = testkpca.getProject(t, KL, v);                //求降維後特征點
         else
             cout << "error" << endl;
         testkpca.saveProject(projectfile2, Project, t);               //存入TXT檔案
     }      

結果分析

工程目錄

4.1、pca原始線性資料集，即工程目錄下test1.txt檔案，共6維資料，每個次元150個特征點，部分資料截圖如下：

經過pca算法降維後，生成工程目錄下pcaproject.txt檔案，資料變為4維：