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

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

首先给出PCA算法的过程：

对每一维度均值化为0，计算出每一维度的均值，每一个数均减去这个均值
对于均值化为零的矩阵，求其维度的协方差矩阵C
求C矩阵的特征值V和特征向量矩阵P
用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;//样本均值化为0X.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;//读取数据//initdouble 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;
}