用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还能够用于预測矩阵中缺失的元素。