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