PCA（主成分分析）原理、步骤及代码实现_如何对普通矩阵进行pca

1、PCA处理步骤：

假设一组数据：m条n维数据

例如点云：
x1, y1, z1            x1, x2, ...xn
x2, y2, z2      ->    y1, y2, ...yn
...                   z1, z2, ...zn
xn, yn zn
1
2
3
4
5

1. 将原始数据组成n行m列矩阵X（每一行为同一字段）
   2.将X的每一行进行零均值化（即每一行数据均减去改行均值）
   3.求出均值化后矩阵的协方差矩阵
   4.求出协方差矩阵对应的特征值、特征向量
   5.将特征向量依据特征值的大小从上到下，按行排列成矩阵，取前K行组成矩阵P
   6.Y = PX，Y即为降维到K维后的数据。

2、数学原理可以参考此处

3、C++借助Eigen库代码实现

//基于Eigen库实现PCA，此代码经本人调试完成，未出现错误
//示例数据：5个2维数据   
/* 
1 1
1 3
2 3
4 4
2 4
*/

#include <iostream>
#include <algorithm>
#include<cstdlib>
#include <fstream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

MatrixXd featurnormal(MatrixXd &X)
{
	// Eigen:基本功能
	// VectorXd    等价于Matrix<double,Dynamic,1>;       //列向量，一列
	// RowVectorXd 等价于Matrix<double,1,Dynamic>;       //行向量，一行
	// MatrixXd    等价于Matrix<double,Dynamic,Dynamic>;

	//计算每一维的均值
	MatrixXd X1 = X.transpose();  //先来个普通转置【行列变换】

	VectorXd meanval = X1.rowwise().mean();//colwise按列求均值。//rowwise按行求均值。 //X1.rowwise()不可用
	cout << "meanval" << endl << meanval << endl;
	
	//样本均值化为0	
	X1.colwise() -= meanval;   //按列减，每一列与均值做差
	
	return X1;
}

void ComputeCov(MatrixXd &X1, MatrixXd &C)
{
	//计算协方差矩阵设我们有m个n维数据记录，将其按列排成n乘m的矩阵X，设C =X*XT/m
    //此处除以m，具体原因你深思[可参考此文](https://blog.csdn.net/hustqb/article/details/78394058)
		
	C = X1*X1.adjoint();//相当于原始数据XT*X    adjiont()求伴随矩阵 相当于求矩阵的转置
	/*cout << "C:" << endl << C<< endl;*/
	C = C.array() / X1.cols();//C.array()矩阵的数组形式  ///逐元素除以m

	cout << "C:" << endl << C << endl;
}

void ComputEig(MatrixXd &C, MatrixXd &vec, MatrixXd &val)
{
	//计算特征向量和特征值 SelfAdjointEigenSolver自动将计算得到的特征向量和特征值排序
	SelfAdjointEigenSolver<MatrixXd> eig(C);
	//所求特征向量依据对应的特征值从小到大排列。
	vec = eig.eigenvectors();
	val = eig.eigenvalues();
	
	//故应该依靠翻转将特征向量依据特征值从大到小排列。【如下】
	vec.colwise().reverse();
	cout << "***************************************" << endl;
	cout << "vec.colwise().reverse():" << endl << vec.colwise().reverse() << endl;
	cout << "***************************************" << endl;
	val.colwise().reverse();
	cout << "val.colwise().reverse():" << endl << val.colwise().reverse() << endl;
	cout << "***************************************" << endl;
}

// 计算维度
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;
	}
	//cout << val.rows() - dim << endl;
	return val.rows() - dim;
}

int main()
{
	//文件操作   文件路径使用： \\ 或者 /
	ifstream fin("F:\\PCLprincipal\\pca070808\\Project1\\m5n2.txt");
	if (fin.is_open())
	{
		cout << "载入成功" << endl;
	}
	ofstream fout("X1X2outpot1.txt");

	//定义所需要的量
	const int m = 5, n = 2, z = 2;
	MatrixXd X(m, n), C(z, z);  //C为协方差阵
	MatrixXd vec, val;

	//读取数据
	double in[z];
	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 << endl;
	cout << "X为原始的数据：" << endl << X << endl;
	cout << endl;
	
	// 1 均值0标准化
	MatrixXd X1 = featurnormal(X);
	cout << "X1均值0标准化的数据：" << endl << X1.transpose() << endl;   ///.transpose()方便为了查看 N*3
	cout << endl;

	// 2 计算协方差
	ComputeCov(X1, C);

	// 3 计算特征值和特征向量
	ComputEig(C, vec, val);

	// 4 计算损失率，确定降低的维数
	int dim = ComputDim(val);
	//int dim = 1;
	// 5计算结果

	//但当涉猎
	//MatrixXd res001 = vec.transpose().rightCols(dim).transpose();
	//MatrixXd res001 = vec.bottomRows(dim);
	//MatrixXd res001 = vec.transpose().leftCols(dim).transpose();
	/*cout <<"X.colwise().reverse()"<<endl<< X.colwise().reverse() << endl;*/ //垂直翻转，优秀的很，

	MatrixXd res001 = vec.colwise().reverse().topRows(dim);  
	//按列翻转取前dim行 特征向量，实现依据特征值从大到小排列获取特征向量。
	
	cout << "res001:" << endl << res001 << endl;  //2*3
	MatrixXd res = res001*X1; //2*3  *  N*3  Y=PX 此Y即为降维之后的结果。

	//方便查看
	cout << "res为用PCA算法之后的数据：" << endl << res.transpose() << endl;
	
	// 6输出结果
	fout << "the result is" << res.rows() << "x" << res.cols() << "after pca algorithm" << endl;
	fout << res.transpose();
	fout.close();
	
	system("pause");
	return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

• 矩阵：依据矩阵有实数矩阵与复数矩阵数据类型分为：
• 转置：共轭转置、普通转置
  1.普通转置：实数矩阵的共轭转置矩阵就是转置矩阵【行与列对换】
  2.共轭转置：复数矩阵的共轭转置矩阵就是将形如a+bi的数变成a-bi,实数的共轭是它本身.
• Eigen库详解1
• Eigen库详解2
• Eigen库详解3
• Eigen库详解4