PCA代码实现详解

#include
#include
#include“Eigen/Eigen“  
#include
using namespace Eigen;  
using namespace std;
using namespace cv;
void BubbleSort（float *pData int Count）;//冒泡排序  
const int num = 5;    // 样本数量
const int dim = 10;    // 样本的维度（一张图片的维度）  
const int pca_dim = 4;// pca dimension of each sample  

int main（）
{
	// 原始数据：每一行代表一个样本（这里为了测试时随机生成的）
	MatrixXf data =MatrixXf::Random（num dim）;
	// 求每一维度上的平均值，即每一列的平均值，最后求出来是一个行向量
	MatrixXf mean =data.colwise（）.mean（）;     
	cout<<“------------------Orignal data-------------------------\n“<	//cout<<“----------------Mean------------------\n“<
	MatrixXf C（num dim）;  // data-mean
	MatrixXf C_T（dim num）;// C的转置
	for（int i=0; i	{
		C.row（i） =data.row（i）-mean;  //C= data-mean
	}
	C_T =C.transpose（）; 
	/*
	协方差矩阵A的维度为dim*dim，这里除不除dim最后结果都一样，
	在书中的定义中是要除的，估计是有数学含义的。
	*/
	MatrixXf A =C_T*C;
	EigenSolver sol_A;
    sol_A.compute（A）;
	//求协方差矩阵的特征值和特征向量
	MatrixXf eigenvals_A =sol_A.eigenvalues（）.real（）;
	MatrixXf eigenvecs_A =sol_A.eigenvectors（）.real（）;
	//排序
	float a[dim]={0.0};
	for（int i=0; i	{
		a[i] =eigenvals_A（i 0）;
	}
	BubbleSort（a dim）; 
	MatrixXf sorted_eigenvals_A（dim 1）;
	MatrixXf sorted_eigenvecs_A（dim dim）;
	//对特征值从大到小排序
	for（int i=0; i	{
		sorted_eigenvals_A（i 0） =a[i];
	}
	//对应的改变特征向量的排序
	for（int i=0; i	{
		for（int j=0; j		{
			if（eigenvals_A（i 0） == a[j]）
			{
				sorted_eigenvecs_A.col（j） =eigenvecs_A.col（i）;
			}
		}
	}
	
	//求PCA的投影矩阵及降维结果
	MatrixXf project_mat（dim pca_dim）; //投影矩阵
	MatrixXf res（num pca_dim）;         //最终输出，维度为num*pca_dim
	for（int i=0; i	{
		project_mat.col（i） =sorted_eigenvecs_A.col（i）;
	}
	
	/*
	res就是原始数据降维后的矩阵
	特别注意，这里是C*project_mat即（data-mean）*project_mat
	在网上看到很多人这里都是直接用的data*project_mat那样是不对的
	*/
	res =C* project_mat;    

	//-----------------------------OpenCV中的PCA--------------------------------//
	CvMat *pca_input         =cvCreateMat（ num dim CV_32FC1）;               // 输入
	CvMat *pca_avg           =cvCreateMat（ 1 dim CV_32FC1）;		          // 平均值 
	CvMat *pca_eigenvalue    =cvCreateMat（ 1 std::min（num dim） CV_32FC1）;  // 特征值  
	CvMat *pca_eigenvector   =cvCreateMat（ std::min（num dim） dim CV_32FC1）;// 特征向量
	CvMat *pca_eigenvector_T =cvCreateMat（dim std::min（num dim） CV_32FC1）; // 特征向量
	CvMat *pca_output        =cvCreateMat（num pca_dim CV_32FC1）;            // 输出
	// 数据导入
	for（int i=0; i	{
		for（int j=0; j		{
			cvSetReal2D（pca_input i j data（i j））;
		}
	}
	cvCalcPCA（pca_input pca_avg pca_eigenvalue pca_eigenvector CV_PCA_DATA_AS_ROW）;
	

	cout<<“\n---------------Eeigen Tool降维结果-----------------“<	cout<
	cout<<“---------------OpenCV Tool降维结果---------------“<	cvProjectPCA（pca_input pca_avg pca_eigenvector pca_output）;
	// 数据输出
	for（i