PCA openCA 中使用 PCA 降维

PCA（principal component analysis）翻译过来就是主分量分析，是一种流行的数据降维方法。通过数据降维可以实现数据的压缩，同时方便数据分析和提高算法的处理速度。PCA的原理就是通过正交变换，最大化样本协方差阵的对角元素，最小化非对角元素。具体的介绍可以参考Shlens, J., A tutorial on principal component analysis. Systems Neurobiology Laboratory, University of California at San Diego, 2005.这个介绍的很不错，从信号处理的角度出发分析了PCA。由于PCA的计算简单，效果显著所以得到很多人的青睐，但是PCA应用本身是基于一定假设的：1.线性。即特征的变换是线性变换，作用有限，目前也有非线性的特征变换kernel PCA。2.处理的数据分布式服从指数族概率密度函数，即能通过均值和协方差来表征数据的分布，因为只有在这个情况下信噪比和协方差矩阵才能表示噪声和数据冗余。好在实际应用时，我们碰到的大多数数据是服从高斯分布或近似高斯分布的（这得益与概率论理的中心极限定理，这个概念经常会碰到，多理解理解有好处）。
PCA是一个无参数的数据分析技术，所以变换结果只受样本影响，结果不受经验影响，与用户独立。就是说PCA无法判断哪个特征比较重要，虽然这个特征的方差比较小。所有特征都平等对待，不能发挥出先验知识的作用。而且经过变换后特征的含义也无法与实际相联系，仅仅是个数据上的表示。毕竟天下没有免费的午餐，有得必有失，PCA是得到的便捷多。
在这里顺便提下PCA、KL变换、奇异值分解的关系，其实这三个变换殊途同归，都实现相类似的功能。在信号处理中都得到广泛应用，真要说有什么区别就是中间的过程细微区别。奇异值分解个人感觉是矩阵论上的一种矩阵分界方法，并不和实际工程性的技术挂钩，只是很多问题如PCA可以通过奇异值分解的方法得到相同的结果，同时奇异值分解也反映了PCA的不同维度的内在转换关系。
在OPENCV的PCA类里提供了PCA的功能，还提供向量投影，和重构的功能，具体可以参考opencv的函数手册。其实这部份是把平均数,共变量矩阵,以及特征值及特征向量都计算出来了,全部都包在cvCalcPCA( const CvArr* data, CvArr* mean,CvArr* eigenvals, CvArr* eigenvects, int flags )的函式里面,因此可以不必特地的用cvCalcCovarMatrix()求得共变量矩阵,也不需要再由共变量矩阵套用cvEigenVV()求得它的EigenValue以及EigenVector了cvCalcPCA()=cvCalcCovarMatrix()+cvEigenVV(),不仅如此,cvCalcPCA()使用上更是灵活,当向量的维度数目比输入的数据大的时候(例如Eigenface),它的共变量矩阵cvCalcCovarMatrix就会自动转成CV_COVAR_SCRAMBLED,而当输入数据量比向量维度大的时候,它亦会自动转成CV_COVAR_NORMAL的形态,而且OpenCV也提供了计算投影量cvProjectPCA(),以及反向投影的函式cvBackProjectPCA()。输入样本可以是一行一个样本也可以试一列一个样本，感觉行方便点。里面用到了OPENCV C++接口里面MAT结构，这个挺好用的，大家可以看看，类似与matlab的矩阵，操作更灵活了。

　opencv中使用步骤：

在OPENCV中使用PCA非常简单，只要几条语句就可以了。

1、初始化数据

//每一行表示一个样本

CvMat* pData = cvCreateMat( 总的样本数, 每个样本的维数, CV_32FC1 );

CvMat* pMean = cvCreateMat(1, 样本的维数, CV_32FC1);

//pEigVals中的每个数表示一个特征值

CvMat* pEigVals = cvCreateMat(1, min(总的样本数,样本的维数), CV_32FC1);

//每一行表示一个特征向量

CvMat* pEigVecs = cvCreateMat( min(总的样本数,样本的维数), 样本的维数, CV_32FC1);

2、PCA处理,计算出平均向量pMean,特征值pEigVals和特征向量pEigVecs

cvCalcPCA( pData, pMean, pEigVals, pEigVecs, CV_PCA_DATA_AS_ROW );

3、选出前P个特征向量(主成份),然后投影,结果保存在pResult中，pResult中包含了P个系数

CvMat* pResult = cvCreateMat( 总的样本数, PCA变换后的样本维数(即主成份的数目), CV_32FC1 );

cvProjectPCA( pData, pMean, pEigVecs, pResult );

4、 重构,结果保存在pRecon中

CvMat* pRecon = cvCreateMat( 总的样本数, 每个样本的维数, CV_32FC1 );

cvBackProjectPCA( pResult, pMean, pEigVecs, pRecon );

5、重构误差的计算

计算pRecon和pData的"差"就可以了.

使用时如果是想用PCA判断“是非”问题，则可以先用正样本计算主成分，判断时，对需要判断得数据进行投影，然后重构，计算重构出的数据与原数据的差异，如果差异在给定范围内，可以认为“是”。

如果相用PCA进行分类，例如对数字进行分类，则先用所有数据（0-9的所有样本）计算主成分，然后对每一类数据进行投影，计算投影的系数，可简单得求平均。即对每一类求出平均系数。分类时，将需要分类得数据进行投影，得到系数，与先前计算出得每一类得平均系数进行比较，可判为最接近得一类。当然这只是最简单得使用方法。

网友chinasun84在vC6.0下使用PCA的经验，正好适合我啊，mark！

最近在研究人脸识别，打算用PCA（主原分析）对图像数据进行降维后用神经网络训练的方法实现，在网上找了一下PCA的C++算法，发现很难用，而且速度奇慢，后来知道opencv上有实现PCA算法的函数，于是下载了一个2.0版，发现原来已经不支持VC6了，由于之前的代码都是在VC6下实现的，现在要移植也不太可能，无奈之下，用了1.0，但是却不知为什么只要维数大于数据量时就出错，真的要放弃？最后就孤注一掷，在VS2008下作一个封装2.0版下PCA类相关函数的库，再在VC6上链接。

　　首先做成了静态库，居然链接不成功，看来VC6与VS2008的静态库还是不兼容，我想如果做动态库链接时应该还是会有这个问题啊，就想放弃了，不过经过高人指点，可以直接用LoadLibrary的方法调用动态库的函数啊。好，这个方法可以完全避开编译器和链接器，是最后的机会了。



　　首先声明一个struct

　　struct PCA

{void* pca;}

这个pca就是指向opencv2.0中PCA对象的指针。



　　接着声明创建和释放该对象的函数。



　　最后再声明其它PCA类方法的函数，如Project。

　　extern "C" void Project(PCA* pca, float* vec, int veclen, float*& res, int reslen);



这里记着最好要加上extern "C"，好让编译器用C函数的方法对该函数进行编译，这样在dll库的函数名就会是Project，否则用C++编译的就会显示加上参数的一串很长的函数名，这时就需要VC6的一个工具Depends打开dll文件才能看到正确的函数名。



　　在VS2008下写好了动态库，就要在VC6下调用了，这个网上就有很多资料了，大概步骤如下：

　　// 声明这个函数指针

　　typedef bool (__cdecl *ProjectPtr)(CPCA*, const data_t*, int, data_t*&, int);



　　// 加载这个动态库

　　HMODULE hModule = LoadLibrary("pcalib.dll");



　　// 取得这个函数的地址

　　ProjectPtr projectPtr = (ProjectPtr)GetProcAddress(hModule, "Project");



　　// 创建对象

　　PCA* pca = (createPCAPtr)(orl_input, TrainNum, InputUnits, MaxComponents);



// 调用这个函数

　 (projectPtr)(pca, orl_input[i], InputUnits, train_input[i], MaxComponents);



　　该方法的确行，而且opencv中的PCA算法速度超快，opevcv果然强大，里面还有很多矩阵运算等算法，一点也不比Matlab差。

#include "cxcore.h"

#include <iostream>

void PrintMatrix(CvMat *Matrix, int Rows, int Cols)
{
for (int i = 0; i < Rows; i++)
{
for (int j = 0; j < Cols; j++)
{
printf("%.2f\t", cvGet2D(Matrix, i, j).val[0]);
}
printf("\n");
}
}

int main()
{
// 样本数少于维数
float a[] = {
1.5,2.3,
3.0,1.7,
1.2,2.9,
2.1,2.2,
3.1,3.1,
1.3,2.7,
2.0,1.7,
1.0,2.0,
0.5,0.6,
1.0,0.9
};

const int rows = 10, cols = 2;

CvMat* mat = cvCreateMat(rows, cols, CV_32FC1);
cvSetData(mat, a, mat->step);

std::cout << "original matrix: " << std::endl;
PrintMatrix(mat, rows, cols);
std::cout << "================================" << std::endl;

const int dim = 2; // dimension
CvMat* avg2 = cvCreateMat(1, cols, CV_32FC1);
CvMat* eigenVector2 = cvCreateMat(dim, cols, CV_32FC1);
CvMat* eigenValue2 = cvCreateMat(dim, 1, CV_32FC1);
//CvMat* vector_pca=cvCreateMat(dim, 3, CV_32FC1);

cvCalcPCA(mat, avg2, eigenValue2, eigenVector2, CV_PCA_DATA_AS_ROW);
//cvProjectPCA(mat,avg2,eigenVector2,vector_pca);

std::cout << "average2: " << std::endl;
PrintMatrix(avg2, 1, cols);
std::cout << "eigenVector2: " << std::endl;
PrintMatrix(eigenVector2, dim, cols);
std::cout << "eigenValues2: " << std::endl;
PrintMatrix(eigenValue2, dim, 1);
std::cout << "================================" << std::endl;

//PrintMatrix(vector_pca, dim, 3);

return 0;
}


我们再来看看另一位作者的博客（写得相当清楚）：







PCA的缘起

PCA大概是198x年提出来的吧，简单的说，它是一种通用的降维工具。在我们处理高维数据的时候，为了能降低后续计算的复杂度，在“预处理”阶段通常要先对原始数据进行降维，而PCA就是干这个事的。本质上讲，PCA就是将高维的数据通过线性变换投影到低维空间上去，但这个投影可不是随便投投，要遵循一个指导思想，那就是：找出最能够代表原始数据的投影方法。这里怎么理解这个思想呢？“最能代表原始数据”希望降维后的数据不能失真，也就是说，被PCA降掉的那些维度只能是那些噪声或是冗余的数据。这里的噪声和冗余我认为可以这样认识：

噪声：我们常说“噪音污染”，意思就是“噪声”干扰我们想听到的真正声音。同样，假设样本中某个主要的维度A，它能代表原始数据，是“我们真正想听到的东西”，它本身含有的“能量”(即该维度的方差，为啥？别急，后文该解释的时候就有啦~)本来应该是很大的，但由于它与其他维度有那么一些千丝万缕的相关性，受到这些个相关维度的干扰，它的能量被削弱了，我们就希望通过PCA处理后，使维度A与其他维度的相关性尽可能减弱，进而恢复维度A应有的能量，让我们“听的更清楚”！

冗余：冗余也就是多余的意思，就是有它没它都一样，放着就是占地方。同样，假如样本中有些个维度，在所有的样本上变化不明显(极端情况：在所有的样本中该维度都等于同一个数)，也就是说该维度上的方差接近于零，那么显然它对区分不同的样本丝毫起不到任何作用，这个维度即是冗余的，有它没它一个样，所以PCA应该去掉这些维度。

这么一分析，那么PCA的最终目的就是“降噪”和消灭这些“冗余”的维度，以使降低维度的同时保存数据原有的特征不失真。后面我们将结合例子继续讨论。

协方差矩阵——PCA实现的关键

前面我们说了，PCA的目的就是“降噪”和“去冗余”。“降噪”的目的就是使保留下来的维度间的相关性尽可能小，而“去冗余”的目的就是使保留下来的维度含有的“能量”即方差尽可能大。那首先的首先，我们得需要知道各维度间的相关性以及个维度上的方差啊！那有什么数据结构能同时表现不同维度间的相关性以及各个维度上的方差呢？自然是非协方差矩阵莫属。协方差矩阵度量的是维度与维度之间的关系，而非样本与样本之间。协方差矩阵的主对角线上的元素是各个维度上的方差(即能量)，其他元素是两两维度间的协方差(即相关性)。我们要的东西协方差矩阵都有了，先来看“降噪”，让保留下的不同维度间的相关性尽可能小，也就是说让协方差矩阵中非对角线元素都基本为零。达到这个目的的方式自然不用说，线代中讲的很明确——矩阵对角化。而对角化后得到的矩阵，其对角线上是协方差矩阵的特征值，它还有两个身份：首先，它还是各个维度上的新方差；其次，它是各个维度本身应该拥有的能量(能量的概念伴随特征值而来)。这也就是我们为何在前面称“方差”为“能量”的原因。也许第二点可能存在疑问，但我们应该注意到这个事实，通过对角化后，剩余维度间的相关性已经减到最弱，已经不会再受“噪声”的影响了，故此时拥有的能量应该比先前大了。看完了“降噪”，我们的“去冗余”还没完呢。对角化后的协方差矩阵，对角线上较小的新方差对应的就是那些该去掉的维度。所以我们只取那些含有较大能量(特征值)的维度，其余的就舍掉即可。PCA的本质其实就是对角化协方差矩阵。

下面就让我们跟着上面的感觉来推推公式吧。假设我们有一个样本集X，里面有N个样本，每个样本的维度为d。即：

PCA的系统分析（包括PCA原理+opencv的PCA函数内部原理和应用+matlab程序）

将这些样本组织成样本矩阵的形式，即每行为一个样本，每一列为一个维度，得到样本矩阵S（N*d）。我们先将样本进行中心化，即保证每个维度的均值为零，只需让矩阵的每一列除以对应的均值即可。很多算法都会先将样本中心化，以保证所有维度上的偏移都是以零为基点的。然后，对样本矩阵计算其协方差矩阵

PCA的系统分析（包括PCA原理+opencv的PCA函数内部原理和应用+matlab程序）
下面，根据我们上文的推理，将协方差矩阵C对角化。注意到，这里的矩阵C是是对称矩阵，对称矩阵对角化就是找到一个正交矩阵P，满足：PCA的系统分析（包括PCA原理+opencv的PCA函数内部原理和应用+matlab程序）具体操作是：先对C进行特征值分解，得到特征值矩阵(对角阵)，，得到特征向量矩阵并正交化即为P.假如我们取最大的前p(p<d)个特征值对应的维度，那么这个p个特征值组成了新的对角阵PCA的系统分析（包括PCA原理+opencv的PCA函数内部原理和应用+matlab程序）对应的p个特征向量组成了新的特征向量矩阵PCA的系统分析（包括PCA原理+opencv的PCA函数内部原理和应用+matlab程序）
实际上，这个新的特征向量矩阵就是投影矩阵.PCA的系统分析（包括PCA原理+opencv的PCA函数内部原理和应用+matlab程序）


样本矩阵S右乘P1相当于每个样本以P1的特征向量为基进行线性变换,得到的新样本矩阵S1中每个样本的维数变为
了p，完成了降维操作。
Matlab中PCA实战:
S = fix(rand(10,3)*50);
S = S - repmat(mean(S),10,1);
C = (S'*S)./(size(S,1)-1);%或者直接： C = cov(S);(matlab的cov(S)用的是(size(S,1)-1）,而我们在opencv编程时候,直接用的是size(S,1)，导致算出来的特征值可能存在差异！！ )
[P,Lambda] = eig(C);

为了验证，我们调用matlab自带的主成分分析函数princomp：
[COEFF,SCORE] = princomp(S) % COEFF表示投影矩阵，SCORE表示投影后新样本矩阵

经过了一番推公式加敲代码的过程，相信大家对主成分分析应该不陌生了吧，同时对协方差矩阵也有了更深层次的认识了吧，它可不只是花花枪啊。我个人觉得PCA在数学上的理论还是很完备的，相必这也是它能在多种应用中博得鳌头的原因吧。
PS：
计算A的协方差矩阵AA'的特征值和特征向量，但是AA'有可能比较大，所以根据AA'的大小，可以计算AA'或者A'A的特征值，原矩阵和其转置矩阵的特征值是一样的，只是特征向量不一样。推导如下：

假如我们的数据按行存放，A是N*n的矩阵（经过样本中心化之后的矩阵，即保证每个维度的均值为零），n>>N，N是样本个数，n是维数，则协方差矩阵应该是A'A，A'A是n*n维的一个矩阵，这个矩阵非常大，不利于求特征值和特征向量，所以先求AA'的特征值，它是一个N*N维的矩阵。
由矩阵性质，AA'的特征值就是A'A的特征值。下面推导A'A的特征向量和AA'的特征向量的关系。

B = A'A; B*x=b*x; C = AA';
C*y=c*y -->
AA'*y=c*y -->
A'A*(A'*y)=c*(A'*y) -->
c = b, x=A'*y

所以根据AA'的特征向量y可以算出A'A的特征向量x。