zoukankan      html  css  js  c++  java

  • PCA方法从原理到实现

    转自：http://blog.csdn.net/celerychen2009/article/details/9048033

    深度神经网路已经在语音识别，图像识别等领域取得前所未有的成功。本人在多年之前也曾接触过神经网络。本系列文章主要记录自己对深度神经网络的一些学习心得。


    第五篇，谈谈PCA模型。本来PCA模型与深度学习是没有任何联系的。通常我们只是用PCA来对机器学习的数据做预处理。




    本来想详细记录一下PCA的原理，但网上已经有一篇不错的文章，链接如下：

    http://hi.baidu.com/ifengzh/item/8851b6387aebefc4382ffa60


    本文前面部分内容引用了这篇文章的内容。


    一、简介 

           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为样本的维数。


    [cpp] view plaincopy
    1. MySample = fix(rand(10,3)*50)  

     

    根据公式，计算协方差需要计算均值，那是按行计算均值还是按列呢，我一开始就老是困扰这个问题。前面我们也特别强调了，协方差矩阵是计算不同维度间的协方差，要时刻牢记这一点。样本矩阵的每行是一个样本，每列为一个维度，所以我们要按列计算均值。为了描述方便，我们先将三个维度的数据分别赋值：

    [cpp] view plaincopy
    1. dim1 = MySample(:,1);  
    2. dim2 = MySample(:,2);  
    3. dim3 = MySample(:,3);  
    4. %计算dim1与dim2，dim1与dim3，dim2与dim3的协方差：  
    5. sum( (dim1-mean(dim1)) .*(dim2-mean(dim2)) ) / ( size(MySample,1)-1 ) %得到  74.5333  
    6. sum( (dim1-mean(dim1)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到  -10.0889  
    7. sum( (dim2-mean(dim2)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到  -10***000  
    8. %搞清楚了这个后面就容易多了，协方差矩阵的对角线就是各个维度上的方差，下面我们依次计算：  
    9. std(dim1)^2 % 得到   108.3222  
    10. std(dim2)^2 % 得到   260.6222  
    11. std(dim3)^2 % 得到  94.1778  
    12. %这样，我们就得到了计算协方差矩阵所需要的所有数据，调用Matlab自带的cov函数进行验证：  
    13. cov(MySample)  

     

    可以看到跟我们计算的结果是一样的，说明我们的计算是正确的。但是通常我们不用这种方法，而是用下面简化的方法进行计算：

    先让样本矩阵中心化，即每一维度减去该维度的均值，然后直接用新的到的样本矩阵乘上它的转置，然后除以(N-1)即可。其实这种方法也是由前面的公式通道而来，只不过理解起来不是很直观而已。大家可以自己写个小的矩阵看一下就明白了。其Matlab代码实现如下：

    [cpp] view plaincopy
    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、计算协方差矩阵的特征向量和特征值

    因为协方差矩阵为方阵，我们可以计算它的特征向量和特征值，如下：

    [cpp] view plaincopy
    1. [eigenvectors,eigenvalues] = eig(cov)  

               

    我们可以看到这些矢量都是单位矢量，也就是它们的长度为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 的实现

               前面已经谈到了PCA的实现分为解析解和在线学习算法。解析解适合于数据量小并且数据完全已知的情况下，这里给出一种高效的解析解的实现代码。

      4.1 数据结构定义及API说明如下：
    [cpp] view plaincopy
    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  

               函数 fce_pca_push_add 用于把每一个样本点添加到PCA模型之中，例如，一个人脸的样本数据。
               函数 fce_pca_solve_eig 采用雅克比迭代法快速求解对称矩阵的特征值和特征向量，其它两个函数分别用以创建PCA模型和释放PCA模型。

    4.2  各函数的实现

    [cpp] view plaincopy
    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. }  
  • 相关阅读:
    start tag, end tag issues in IE7, particularly in xslt transformation
    用SandCastle为注释生成chm文档
    Firebug
    架构的重点
    Linux Shell常用技巧(十) 管道组合
    Linux JDK升级
    Linux Shell常用技巧(十二) Shell编程
    Packet Tracer 5.0实验(一) 交换机的基本配置与管理
    Linux Shell常用技巧(六) sort uniq tar split
    Linux Shell常用技巧(二) grep
  • 原文地址：https://www.cnblogs.com/walccott/p/4957058.html
Copyright © 2011-2022 走看看