PCA降维的计算原理与应用
在科研、工程应用、生活中,我们所获取的数据往往包含着很多冗余信息,这些冗余信息往往对数据分析造成干扰,增加数据分析的复杂度。此时我们则需要对这些数据进行预处理,预处理的原则是:既能抓住其主要特征,又能剔除冗余信息,从而减少数据量。PCA降维就是这样的一种数据预处理算法。
本文首先讲解PCA降维的计算原理,再使用C++与Opencv来实现该算法,并与Opencv现有的PCA函数接口进行降维结果的对比。看到这里,可能有人会问,Opencv都有现成的函数可以调用了,为什么还要自己去写呢?我想说的是,对于学习者来说,重复造轮子并不是坏事,它可以让我们更加深刻的理解造轮子的过程,从而才有改进和创新的机会,相反,如果只是使用别人造好的轮子,不深刻理解其构造原理,压根就没有改进创新的空间了。
1. 计算原理
假设有m行n列的数据,计为矩阵X0,其每一行数据看作一个一维行向量,那么该数据本来有m个一维行向量,我们要使用PCA降维算法把其降为k个一维行向量(k < m),计算过程如下:
(1) 求出每行数据的平均值。
(2) 去平均处理,把每行数据都减去本行数据的平均值。
(3) 去平均处理之后,同样得到m行n列的数据,计为矩阵X1。
(4) 按以下公式计算X0的协方差矩阵,其中“*”表示矩阵乘法,得到的协方差矩阵Cov为m行m列矩阵:
(5) 计算协方差矩阵Cov的特征值与对应的特征向量。得到m个特征值,对应m个特征向量,其中每个特征向量的长度又为m,也即所有特征向量组成m行m列的矩阵。
(6) 将特征值按照从大到小排列,并根据排列后特征值的顺序来按行从上到下排列特征向量(每个特征向量为一行),使特征值与特征向量仍保持对应。
比如本来特征值依次为a1,a2,a3,a4,a5,对应的特征向量为:
v1
v2
v3
v4
v5
排序之后,a3>a1>a4>a5>a2,那么特征向量的顺序也作对应的调整:
v3
v1
v4
v5
v2
此时计特征向量组成的矩阵为V。
(7) 取矩阵V的前k行数据,得到k行m列的矩阵P,计算Y=P*X1,Y矩阵即为最终得到的k行n列的降维数据,从而实现把数据从m维降为n维。
上代码:
void do_PCA(Mat src, Mat &pca, int k)
{
if (src.rows <= 1 || src.cols <= 1)
{
return;
}
k = k > src.rows ? src.rows : k;
Mat src_float;
src.convertTo(src_float, CV_32F);
//求每行的平均值
Mat col_mean = Mat::zeros(1, src.rows, CV_32FC1);
float *col_mean_p = col_mean.ptr<float>(0);
for (int i = 0; i < src_float.rows; i++)
{
float *p = src_float.ptr<float>(i);
for (int j = 0; j < src_float.cols; j++)
{
col_mean_p[i] += p[j];
}
}
col_mean /= src_float.cols;
//去平均值
for (int i = 0; i < src_float.rows; i++)
{
float *p = src_float.ptr<float>(i);
for (int j = 0; j < src_float.cols; j++)
{
p[j] = p[j] - col_mean_p[i];
}
}
//计算协方差矩阵
Mat X = src_float*src_float.t();
X = X / src_float.cols;
Mat eValuesMat, eVectorsMat;
//调用opencv接口计算特征值与特征向量
eigen(X, eValuesMat, eVectorsMat); //这里得到的特征值已经按照从大到小排序了,特征向量也与特征值相对应
Mat Vectors_k;
eVectorsMat(Rect(0, 0, eVectorsMat.cols, k)).copyTo(Vectors_k); //取特征向量的前k行, k*r
pca = Vectors_k*src_float; //k*r * r*c = k*c
pca = pca.clone(); //确保矩阵连续,拷贝一份
printf("pca.rows=%d, pca.cols=%d\n", pca.rows, pca.cols);
}
测试代码:
void pca_test(void)
{
Mat img = imread("lena.jpg", CV_LOAD_IMAGE_GRAYSCALE);
int k = 100; //从图像原有的m行降为100行,列不变
Mat pca_m;
do_PCA(img, pca_m, k);
imshow("原图", img);
imshow("本文实现算法 PCA降维后", pca_m);
//使用Opencv的PCA函数接口
PCA pca(img, Mat(), CV_PCA_DATA_AS_COL, k);
Mat dst = pca.project(img); //dst则为最终降维数据
imshow("Opencv PCA降维后", dst);
waitKey(0);
}
运行上述代码,得到的结果如下。可以看到,本文实现的算法与Opencv函数的计算结果使一致的,说明我们的计算过程没错,鼓掌~再接再厉~
原图
本文实现的PCA算法的降维结果
Opencv现有的PCA算法的降维结果
微信公众号如下,欢迎关注,欢迎交流: