PCA降维的计算原理与应用(2)
学习PCA降维算法的时候，在网上看到过两个不同版本的计算过程，一直有点迷糊，到底哪个版本才是对的。后来发现，两个版本的计算方法都没错，区别主要在于把每行看作一维向量，还是把每列看作一维向量。所以本文的主要目的就是总结和对比一下这两种过程略有不同的计算方法。
1. 把每行看作一个一维向量
该计算方法就是我们在之前一篇讲PCA降维算法的文章中所讲述的方法，其对应Opencv接口中的CV_PCA_DATA_AS_COL模式，把每行看作一个一维向量，通过PCA降维减少一维向量的个数，所以从整体来看，矩阵数据的列不变，但是行减少了。如下图所示，矩阵的原本尺寸是m行n列，降维之后就变成了k行n列的矩阵的(k < m)：
在上篇文章中我们已经详细讲了这种计算方法，此处再重复讲解一下，方便读者对比。假设有m行n列的数据，计为矩阵X0，则有：
(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列的降维数据，从而实现数据降维。
1. 把每列看作一个一维向量
该计算方法对应Opencv接口中的CV_PCA_DATA_AS_ROW模式，把每列看作一个一维向量，通过PCA降维减少一维向量的个数，所以从整体来看，矩阵数据的行不变，但是列减少了。如下图所示，矩阵的原本尺寸是m行n列，降维之后就变成了m行k列的矩阵(k < n)：
同样假设有m行n列的数据，计为矩阵X0，则该方法的计算过程如下：
(1) 求出每列数据的平均值。
(2) 去平均处理，把每列数据都减去本列数据的平均值。
(3) 去平均处理之后，同样得到m行n列的数据，计为矩阵X1。
(4) 按以下公式计算X0的协方差矩阵，其中“*”表示矩阵乘法，得到的协方差矩阵Cov为n行n列矩阵：
(5) 计算协方差矩阵Cov的特征值与对应的特征向量。得到n个特征值，对应n个特征向量，其中每个特征向量的长度又为n，也即所有特征向量组成n行n列的矩阵。
(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行n列的矩阵P，计算Y=X1*PT，Y矩阵即为最终得到的m行k列的降维数据，从而实现数据降维。
以上计算过程的实现代码如下：
void do_PCA_col(Mat src, Mat &pca, int k)
{
if (src.rows <= 1 || src.cols <= 1)
{
return;
}
k = k > src.cols ? src.cols : k;
Mat src_float;
src.convertTo(src_float, CV_32F);
//求每列的平均值
Mat col_mean = Mat::zeros(1, src.cols, 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[j] += p[j];
}
}
col_mean /= src_float.rows;
//每列去平均值
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[j];
}
}
//计算协方差矩阵
Mat X = src_float.t()*src_float; //n*m * m*n = n*n
X = X / src_float.rows;
Mat eValuesMat, eVectorsMat;
eigen(X, eValuesMat, eVectorsMat); //这里得到的特征值已经按照从大到小排序了，特征向量也与特征值相对应
Mat Vectors_k;
eVectorsMat(Rect(0, 0, eVectorsMat.cols, k)).copyTo(Vectors_k); //取特征向量的前k行, k*n
//printf("Vectors_k.rows=%d, Vectors_k.cols=%d\n", Vectors_k.rows, Vectors_k.cols);
pca = src_float*Vectors_k.t(); //m*n * n*k = m*k
pca = pca.clone();
}
测试代码如下，分别调用以上实现的函数，以及Opencv现有的PCA接口函数，对472行496列的Lena图像进行列降维。
void pca_test(void)
{
Mat img = imread("lena.jpg", CV_LOAD_IMAGE_GRAYSCALE);
int k = 200; //把列降维到200列
Mat pca_m;
do_PCA_col(img, pca_m, k);
imshow("img", img);
imshow("pca_m", pca_m);
//使用Opencv的PCA函数接口，模式为CV_PCA_DATA_AS_ROW
PCA pca(img, Mat(), CV_PCA_DATA_AS_ROW, k);
Mat dst = pca.project(img);//映射新空间
imshow("PCA降维后", dst);
waitKey(0);
}
运行结果如下。可以看到本文实现的算法与Opencv函数在CV_PCA_DATA_AS_ROW模式下的计算结果是一致的，均得到472行200列的降维数据。
原图
本文实现的PCA算法的降维结果
Opencv现有的PCA算法的降维结果
微信公众号如下，欢迎扫码关注：
