PCA算法的原理C++ Eigen库实现（附源码下载）

PCA的目的：  pca算法，也叫主成分分析法，能够对一个多样本的多维特征向量构成的矩阵进行分析，分析主成分并去除维度之间的相关性，使线性相关的向量组变成线性无关的向量组。  并且可以对样本进行降维，降高维向量映射到低维度空间，同时确保纬度之间的信息损失尽可能小。

首先给出PCA算法的过程：

1. 对每一维度均值化为0，计算出每一维度的均值，每一个数均减去这个均值
2. 对于均值化为零的矩阵，求其维度的协方差矩阵C
3. 求C矩阵的特征值V和特征向量矩阵P
4. 用P作基底，将原来的特征矩阵映射到这一正交基上

基础知识：

首先理解正交基分解：

我们默认一个n维向量v=（x1,x2,x3……xn）是一个n维空间的以n个单位基底（1,0,0,……），（0,1,0,0……）……所表示的坐标形式。

例如，向量（1,2）可以表示为以(1,0),(0,1)为基底的一个向量。以下称这种基底为默认基底。

同时，在一个n维向量空间中，理论上是有无数多对基底的，要从默认基底变化到任意基底的坐标表示情况，首先要明确，向量v在n维空间上i维的坐标可以表示为v与第i个基底的内积（数量积），即可表示为v=(vi, vk, vj……)其中i，j，k……为选取的基底集合。

那么由m个样本构成的n维特征，表示为n行m列的矩阵X。则将在默认基底下的矩阵X变换到以n个k维为基底的表示，其矩阵表示为Y=PX ，其中P是kxn的n个变换基向量构成的矩阵，Y是变换结果。

那么一个n维的特征矩阵就可以通过正交基P变换到k维空间。

---

接下来写一下为什么要这样做，为什么通过这一系列操作，就能达到PCA的目的呢？

首先理解一下协方差矩阵：  协方差矩阵的对角线上的值代表这一维度的方差，而方差能够代表能量，方差越大，这一维度的能量越大（代表特征的能力越大）。非对角线上的元素代表了两个维度之间的协方差，协方差越小，相关性越小，协方差为0时，两个维度线性无关。

所以假设我们最后经过PCA处理要得到矩阵Y，Y是通过P这一正交基映得到的，即Y = PX， 那么我们要尽可能的使Y的协方差矩阵主对角线上的元素值尽可能大，而非对角线上的元素值尽可能为0。

假设Y的协方差矩阵为D，则：

我们要找的P恰好就是能让C对角化的P，根据线性代数知识，我们知道要使一个矩阵对角化，只需要求其特征值和特征向量，使用特征向量就能通过矩阵乘法将其对角化。

所以P就是C的特征向量矩阵。

这也就是我们所做的要求协方差矩阵，求其特征值和特征向量的原因。

但为了达到降维的目的，我们不能单纯的取全部特征向量来做，这样做的话维度是不变的，只是去除了一些相关性。所以我们要选取一部分特征向量来用。

衡量特征向量好坏的标准就是特征值，特征值越大，其对应的特征向量的维度分量上能量越大。我们对特征值进行从大到小排序，选取前K大的特征值对应的向量作为基底。

K的选取取决于我们想保留多少信息，这个计算方式为 sum（前k特征值）/ sum(所有特征值)

接下来附c++实现代码并带有注释：

#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<fstream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
void featurenormalize(MatrixXd &X)
{
	//计算每一维度均值
	MatrixXd meanval = X.colwise().mean();//每一列的矩阵，列降维
	RowVectorXd meanvecRow = meanval;
	//样本均值化为0
	X.rowwise() -= meanvecRow;
	
}
void computeCov(MatrixXd &X, MatrixXd &C)
{
	//计算协方差矩阵C = XTX / n-1;
	C = X.adjoint() * X;
	//C = X * X.adjoint() ;
	MatrixXd Y;
	Y = X.adjoint();//转置矩阵+
	for (int i = 0; i < Y.rows(); i++)
	{
		for (int j = 0; j < Y.cols(); j++)
		{
			cout << " " << Y(i, j);
		}
		cout << " " << endl;
	}
	//C = C.array() / X.rows() - 1;
	//C = C.array() / (X.rows() - 1);
	C = C.array() / (X.cols() - 1);
	for (int i = 0; i < C.rows(); i++)
	{
		for (int j = 0; j < C.cols(); j++)
		{
			cout << " " << C(i, j);
		}
		cout << " " << endl;
	}
}
void computeEig(MatrixXd &C, MatrixXd &vec, MatrixXd &val)
{
	//计算特征值和特征向量，使用selfadjont按照对阵矩阵的算法去计算，可以让产生的vec和val按照有序排列
	SelfAdjointEigenSolver<MatrixXd> eig(C);

	vec = eig.eigenvectors();
	val = eig.eigenvalues();
}
int computeDim(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.95)
			break;
	}
	//std::cout << "sum: " << sum << " val.sum(): " << val.sum() << "  rows: " << val.rows() << "  dim: " << dim << std::endl;
	return val.rows() - dim;
}
int main()
{
	ifstream fin("siftsmallD.txt");
	ofstream fout("output.txt");
	/*const int m = 10000, n = 128;
	MatrixXd X(10000, 128), C(128, 128);*/
	const int m = 8, n = 6;
	MatrixXd X(m, n), C(3, 3);
	MatrixXd vec, val;

	//读取数据
	//init
	double in[200];
	for (int i = 0; i < m; ++i)
	{
		for (int j = 0; j < n; ++j)
			fin >> in[j];
		for (int j = 1; j <= n; ++j)
		{
			X(i, j - 1) = in[j - 1];
			//cout << " " << X(i, j - 1);
		}
		//cout << " " << endl;
	}
	
	//cout << "X.rows(): " << X.rows() << "  X.cols(): " << X.cols() << endl;

	//pca
	//零均值化
	featurenormalize(X);
	//计算协方差
	computeCov(X, C);
	//计算特征值和特征向量
	computeEig(C, vec, val);
	//计算损失率，确定降低维数
	int dim = computeDim(val);
	std::cout << "dim: " << dim << std::endl;
	//计算结果
	MatrixXd res = X * vec.rightCols(dim);
	//输出结果
	fout << "the result is " << res.rows() << "x" << res.cols() << " after pca algorithm." << endl;
	fout << res;
	system("pause");
	return 0;
}