理解PCA原理与C++\Matlab实现

1. PCA原理

图像处理等领域经常会用到主成分分析，这样做的好处是使要分析的数据的维度降低了，但是数据的主要信息还能保留下来。它的原理是这样的，对于给定一组数据（列向量）：



将其中心化后表示为： 



其中u为输入列向量的均值。中心化后的数据在第一主轴u1（既是数据的主方向，这里假设为单位向量）方向上分布散的最开，也就是说在u1方向上的投影的绝对值之和最大（也可以说方差最大），计算投影的方法就是将x与u1做内积，由于只需要求u1的方向，所以设u1是单位向量。也就是最大化下式：



也即最大化：



采用平方可以把绝对值符号拿掉，光滑曲线处理起来方便。两个向量做内积可以转化成矩阵乘法的形式，所以目标函数括号里面就是矩阵乘法表示内积，转置以后的行向量乘以列向量得到一个数。因为一个数的转置还是其本身，所以又可以将目标函数化为：



经过上述的过程已经将绝对值符号去掉，将原来的问题转换成为了对等的问题，之后便可以用矩阵之间的乘法关系将上述用矩阵乘法定义展开的式子进行进一步合并，这样就可以把括号去掉,去掉以后由于u1（之前定义的最大主方向）和i无关是的变量，可以把它拿到求和符外面：



其实括号里面是一个矩阵乘以自身的转置。X矩阵的第i列就是xi，于是目标函数最后化为：



在没有前面的1/n，那就是一个标准的二次型！并且XX'(为了方便，用'表示转置)得到的矩阵是一个半正定的对称阵（具体的证明过程可以参考矩阵教材资料）。由于是条件极值，便直接向导拉格朗日算子：则最优化问题可以转换为：
目标函数和约束条件构成了一个最大化问题（之前假设了u1是单位向量）：



构造拉格朗日函数：



对u1求导数就可以得到函数的极值，由于最优化函数是严格的凸函数，则得到：


从上式中可以直观看出u1即为XX'特征值λ对应的特征向量！XX'的所有特征值和特征向量都满足上式。所以，取最大的那个特征值，那么得到的目标值就最大。因而在图像进行压缩的过程中选取最大的几个特征值，便可以对图像进行压缩，裁剪掉相对来说“冗余”的部分。
至于PCA的运用，应该是PCA方法在人脸识别中特征提取及数据维，对于输入200*200大小的人脸图像，单单提取它的灰度值作为原始特征，则这个原始特征将达到40000维，这给后面分类器的处理将带来极大的难度。著名的人脸识别Eigenface算法就是采用PCA算法，用一个低维子空间描述人脸图像，同时用保存了识别所需要的信息。基于K-L变换的高斯人脸检测算法在多姿态、有遮挡物情景下检测率很高，实时性好，鲁棒性强，对监控系统的智能化发展具有重要的实际意义。
在PCA主成成分分析中有两个运用到的理论：最大方差理论、最小错误理论。在信号处理中认为信号具有较大的方差，噪声有较小的方差，信噪比就是信号与噪声的方差比，越大越好。因此认为，最好的k维特征是将n维样本点转换为k维后，每一维上的样本方差都很大。最小错误理论是选取合适的特征向量使得恢复的图像数据尽量与原始数据之间的差别最小化，需要在数据降维程度和数据的损失之间选取平衡。

2. 实现

2.1 基于Opencv的实现

Opencv中涉及到PCA的函数主要有：

//Opencv中PCA类的主要函数有：

//参数: InputArray data			原始矩阵数据
//参数: InputArray mean			原始矩阵数据均值，输入空函数会自己计算
//参数: int flags				每行(CV_PCA_DATA_AS_ROW)/列(CV_PCA_DATA_AS_COL)代表一个样本
//参数: int maxComponents = 0	保留多少特征值，默认全保留
PCA::PCA(InputArray data, InputArray mean, int flags, int maxComponents = 0)

//参数: InputArray data			原始矩阵数据
//参数: InputArray mean			原始矩阵数据均值，输入空函数会自己计算
//参数: int flags				每行(CV_PCA_DATA_AS_ROW)/列(CV_PCA_DATA_AS_COL)代表一个样本
//参数: double retainedVariance	保留多少特征值，百分比，默认全保留
PCA::PCA(InputArray data, InputArray mean, int flags, double retainedVariance)

//参数: InputArray data			原始矩阵数据
//参数: InputArray mean			原始矩阵数据均值，输入空函数会自己计算
//参数: int flags				每行(CV_PCA_DATA_AS_ROW)/列(CV_PCA_DATA_AS_COL)代表一个样本
//参数: double retainedVariance	保留多少特征值，百分比，默认全保留
PCA& PCA::computeVar(InputArray data, InputArray mean, int flags, double retainedVariance)
//参数依次为：原始数据；原始数据均值，输入空会自己计算；每行/列代表一个样本；保留多少特征值，百分比，默认全保留

//原图像，投影到新的空间
Mat PCA::project(InputArray vec) const

//先进行project之后的数据，反映摄到原始图像
C++: Mat PCA::backProject(InputArray vec) const

Opencv测试例子：

const int row_num = 4;
const int col_num = 4;
float array[] =
{1.2, 3.2, 4.1, 6.1,
3.2, 0.23, 4.1, 3.1,
5.1, 0.9, 0.12, 1.9,
3.2, 5.2, 9.0, 2.4};

cv::Mat pcaSet(row_num, col_num, CV_32FC1, cv::Scalar::all(0));
float* data = nullptr;
for (int i = 0; i < row_num; i++)
{
data = pcaSet.ptr<float>(i);
for (int j = 0; j < col_num; j++)
{
*data++ = array[i*col_num + j];
}
}

cout << "origin Matrix:\n" << pcaSet << endl;

cv::PCA pca(pcaSet, cv::Mat(), CV_PCA_DATA_AS_ROW, 3);//参数依次为：原始数据；原始数据均值，输入空会自己计算；每行/列代表一个样本；保留多少特征值，默认全保留

cout << "对特征矩阵进行PCA变换之后的到的矩阵均值：" << pca.mean << endl;//均值

cout << "对特征矩阵进行PCA变换之后的到的特征值：\n" << pca.eigenvalues << endl;//特征值

cout << "对特征矩阵进行PCA变换之后的到的特征值向量：\n" << pca.eigenvectors << endl;//特征向量

cv::Mat dst = pca.project(pcaSet);//映射新空间
cout << "对特征矩阵进行PCA变换之后的到的映射空间：\n" << dst << endl;//特征向量

2.2 Matalb实现

function [Pro_Matrix,Mean_Image]=my_pca(Train_SET,Eigen_NUM)
%输入：
%Train_SET：训练样本集，每列是一个样本，每行一类特征，Dim*Train_Num
%Eigen_NUM：投影维数

%输出：
%Pro_Matrix：投影矩阵
%Mean_Image：均值图像

[Dim,Train_Num]=size(Train_SET);

%当训练样本数大于样本维数时，直接分解协方差矩阵
if Dim<=Train_Num
Mean_Image=mean(Train_SET,2);
Train_SET=bsxfun(@minus,Train_SET,Mean_Image);
R=Train_SET*Train_SET'/(Train_Num-1);

[eig_vec,eig_val]=eig(R);
eig_val=diag(eig_val);
[~,ind]=sort(eig_val,'descend');
W=eig_vec(:,ind);
Pro_Matrix=W(:,1:Eigen_NUM);

else
%构造小矩阵，计算其特征值和特征向量，然后映射到大矩阵
Mean_Image=mean(Train_SET,2);%求取每一行的平均值
Train_SET=bsxfun(@minus,Train_SET,Mean_Image);%每一行减去平均值
R=Train_SET'*Train_SET/(Train_Num-1);   %转换为200*200的方阵 之后矩阵中元素每个除以199

[eig_vec,eig_val]=eig(R);   %计算特征向量 特征值
eig_val=diag(eig_val);  %获取特征值向量
[val,ind]=sort(eig_val,'descend');  %对特征值进行排序 获取特征值所在的位置大小的排序
W=eig_vec(:,ind);   %获取特征值对应的特征向量
Pro_Matrix=Train_SET*W(:,1:Eigen_NUM)*diag(val(1:Eigen_NUM).^(-1/2));
Pro = Train_SET*W(:,1:Eigen_NUM)*diag(val(1:Eigen_NUM).^(-1/2));
end

end

3. 结果