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;
}