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算法的降维结果