PCA算法的数学原理和C++语言（Eigen库）实现

PCA算法的数学原理

    最近在学习图像处理相关方面的知识，在图像压缩时用到主成分分析算法（Principal Component Analysis PCA）。数学理论主要参考了这篇博客点击打开链接，博主写的非常好，通俗易懂。这里总结了一下PCA算法的实现步骤如下：

设有m条n维数据。

1）将原始数据按列组成n行m列矩阵X；

2）将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值；

3）求出协方差矩阵C=1mXXT

4）求出协方差矩阵的特征值及对应的特征向量；

5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P；

6）Y=PXY=PX即为降维到k维后的数据；

C++语言（Eigen库）实现

这里用C++语言和Eigen库实现了点击打开链接这篇博客的例子，实验的数据和结果可以从下面的图中看出，与原文的结果是一样的。

#include <iostream>
#include <algorithm>
#include <fstream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
 

MatrixXd featurnormail(MatrixXd &X)
{
	//计算每一维的均值
	MatrixXd X1 = X.transpose();
	MatrixXd meanval = X1.colwise().mean();

//	cout << "meanval" << endl << meanval <<endl;
	//样本均值化为0
	RowVectorXd meanvecRow = meanval;
	X1.rowwise() -= meanvecRow;
	return X1.transpose();
}
void ComComputeCov(MatrixXd &X, MatrixXd &C)
{
	//计算协方差矩阵
	C = X*X.adjoint();//相当于XT*X adjiont()求伴随矩阵 相当于求矩阵的转置
	C = C.array() / X.cols();//C.array()矩阵的数组形式
}

void ComputEig(MatrixXd &C, MatrixXd &vec, MatrixXd &val)
{
	//计算特征向量和特征值 SelfAdjointEigenSolver自动将计算得到的特征向量和特征值排序
	SelfAdjointEigenSolver<MatrixXd> eig(C);
	vec = eig.eigenvectors();
	val = eig.eigenvalues();
}

// 计算维度
int ComputDim(MatrixXd &val)
{
	int dim;
	double sum = 0;
	for (int i = val.rows() - 1; i >= 0;--i)
	{
		sum += val(i, 0);
		dim = i;
		if (sum / val.sum()>=0.8)//达到所需要的要求时
			break;
	}
	return val.rows() - dim;
}

int main()
{
	//测试
	Eigen::VectorXf v(3);
	v(0) = 2;
	v(1) = 9;
	v(2) = 7;
	Eigen::Matrix2f M;
	M(0, 0) = 2;
	M(0, 1) = 9;
	M(1, 0) = 6;
	M(1, 1) = 12;

	cout << "V:" << endl << v << endl;
	cout << endl;
	cout << "V:" << endl << v.transpose() << endl;
	cout << "M:" << endl << M << endl;
	cout << endl;
	cout << "M:" << endl << M.transpose() << endl;


	//正式开始程序
	//文件操作
	ifstream fin("test.txt");
	ofstream fout("outpot1.txt");
	ofstream fout1("X1.txt");

	//定义所需要的量
	const int m = 2, n = 5, z = 2;
	MatrixXd X(m, n), C(z, z);
	MatrixXd vec, val;
	//读取数据
	double in[200];
	for (int i = 0; i < m; ++i)
	{
		for (int j = 0; j < n; ++j)
		{
			fin >> in[j];
			//cout << "in[" << j << "]" << in[j] << endl;
		}
		for (int j = 1; j <= n; ++j)
			X(i, j - 1) = in[j - 1];
	}
	cout << endl;
	cout << "X为原始的数据：" << endl << X << endl;
	cout << endl;
	// 1 均值0标准化
	MatrixXd X1=featurnormail(X);
	cout << "X1均值0标准化的数据：" << endl << X1 << endl;
	cout << endl;
	// 2 计算协方差
	ComComputeCov(X1, C);

	// 3 计算特征值和特征向量
	ComputEig(C, vec, val);

	// 4 计算损失率，确定降低的维数
	int dim = ComputDim(val);

	// 5计算结果
	MatrixXd res = vec.rightCols(dim).transpose()*X1;
	cout << "res为用PCA算法之后的数据：" << endl << res << endl;
	cout << endl;
	// 6输出结果
	fout << "the result is" << res.rows() << "x" << res.cols() << "after pca algorithm" << endl;
	fout << res;

	cout << "I Love n(Yaner)......" << endl;
	system("pause");
	return 0;
}