PCA算法的原理以及c++实现(Eigen库实现)源码

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/Qer_computerscience/article/details/71246634

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 = C.array() / X.rows() - 1;
}
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;
	}
	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);
	MatrixXd vec, val;
	
	//读取数据
	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];
	}
	//pca

	//零均值化
	featurenormalize(X);
	//计算协方差
	computeCov(X, C);
	//计算特征值和特征向量
	computeEig(C, vec, val);
	//计算损失率,确定降低维数
	int dim = computeDim(val);
	//计算结果
	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;
}






展开阅读全文

没有更多推荐了,返回首页