主成份分析（PCA）基本原理/步骤及其C++ 实现与优化（结合Eigen矩阵库）

主成份分析是常用的降维方法，其他降维方法还有线性判别分析LDA，二者的区别见：https://www.cnblogs.com/pinard/p/6244265.html   简要说就是：

1.PCA将原始数据投影到方差最大的方向，LDA将数据投影到不同样本的中心点距离最大的方向。

2. PCA是无监督降维，LDA是有监督降维。

3. 若分类主要依赖均值而非方差，则LDA效果好，反之PCA效果好

 

PCA 的主要步骤:

1. 使用PCA之前必须进行特征缩放！feature scaling

2. 计算特征矩阵X的协方差矩阵Sigma

sigma = 1/m × X^T × X;

相关文献大部分公式都是要求计算协方差矩阵之前先将特征矩阵每一个维度减去平均值，这样是让数据分布以原点为中心，但并非必要，并不影响对数据分布方差的分析。因为协方差矩阵和PCA主要考虑方差而非均值，与LDA正好相反。

协方差矩阵描述了样本的分布形状。

m 是样本数，即特征矩阵X的行数。X 的维度是 m×n，n 是特征向量的维度，即降维之前原始特征数。

得到的协方差矩阵sigma 是 n×n 的矩阵

3. 对协方差矩阵进行奇异值分解。

奇异值分解的几何意义这篇文章讲的特别好：  https://blog.csdn.net/jinshengtao/article/details/18448355

[U, S, V]  = svd(sigma);

U，S，V 都是n×n的矩阵

奇异值（特征值）描述了数据分布的形状。最大特征值（奇异值）对应的特征向量指向数据主要分布方向，即方差最大的方向！

->协方差矩阵特征值从小到大排列对应的特征向量指向数据分布的方差从大到小的方向。协方差矩阵特征值不受刚性变换的影响，而特征向量受刚性变换的影响！

其中 U 是 包含左奇异向量的矩阵，V 是包含右奇异向量的矩阵。S 是一个对角阵，对角线上的元素都是奇异值：s11, s22, s33, ..., snn，奇异值在S中从大到小排列.  特征向量即PCA需要将数据投影的方向！为什么PCA要将数据投影到特征向量的方向即方差最大的方向呢？因为数据的分布无非是用均值和方差来表征，PCA主要考虑方差，投影后保留大部分的方差就意味着保留数据分布的大部分特征！使得样本数据往低维投影后，能尽可能表征原始的数据。

下面这张图很关键：

这样按照上图，就可以取U的前k列，作为Uredue，降维后的特征矩阵 Xreduce = X × Ureduce

将特征向量矩阵取前k列，与原矩阵相乘，这样的几何意义是将原矩阵投影到k个特征向量上，因为矩阵乘法的意义就是一个变换矩阵作用于另一个矩阵X。

协方差矩阵的几何意义详见我这篇博客：https://blog.csdn.net/shaozhenghan/article/details/81291988

C++代码（结合Eigen矩阵库）

do_pca.cpp

对之前的不等式等价变换，如上图。变换后的不等式右边项是固定值，代码实现时放在for循环外面。左边项是累加，每次循环都比上次循环多加一个数。因此把这个累加和定义在循环外面，每次在原来的基础上加一个数。这样就不用每次从头加起。

float sum_sing_part = 0.0;

unsigned int k = 0;

while (k < S.rows())

{      

       sum_sing_part  +=  S.row(k).sum();

       ..........

 

#include "do_pca.h"

using namespace std;
using namespace Eigen;

bool pca (const MatrixXf & X, MatrixXf & X_reduced, const float variance_remain)
{
    // m: number of rows of original data set
    unsigned int m = X.rows();
    // Covariance Matrix Sigma 
    MatrixXf Sigma = 1.0 / m * (X.transpose() * X); 
    // SVD decomposition: [U, S, V] = svd(Sigma);
    JacobiSVD<MatrixXf> svd(Sigma, ComputeFullU);
    // left_singular_matrix
    MatrixXf U = svd.matrixU();
    // singular values vector
    MatrixXf S = svd.singularValues();
    cout << "\n S = \n" << S << endl; // debug
    // (variance_remain*100)% of variance should be retained
    if(variance_remain < 0 || variance_remain > 1.0)
    {
        cout << "\n variance_remain should in [0.0, 1.0]! \n" << endl;
        return(false);
    }
    float sum_sing_remained = variance_remain * S.sum();

    cout << "\n S.sum() = " << S.sum() << endl;  // debug

    float sum_sing_part = 0.0;
    unsigned int k = 0;
    while (k < S.rows())
    {
        sum_sing_part += S.row(k).sum();
        cout << "\n sum_sing_part = " << sum_sing_part << "for k = " << k << endl; // debug 

        if (sum_sing_part >= sum_sing_remained)
        {
            cout << "\n" 
                << " more than " 
                << 100*variance_remain 
                << "% of variance is retained for k = " 
                << k 
                << endl; 
            break;
        }
        ++k;
    }
    // Uk: n*(k+1)
    MatrixXf Uk = U.leftCols(k + 1); 
    // X_reduced: m * (k+1)
    X_reduced = X * Uk; 

    return (true);
}

 

do_pca.h

#ifndef DO_PCA_H
#define DO_PCA_H

#include <pcl/common/eigen.h>

bool pca (const Eigen::MatrixXf & X, Eigen::MatrixXf & X_reduced, const float variance_remain);


#endif

 

写一个测试代码：用随机数矩阵测试一下。

test_pca.cpp

#include "do_pca.h"
#include <ctime>
#include <iostream>

using namespace std;
using namespace Eigen;

int main(int argc, char const *argv[])
{
    srand((unsigned)time(NULL));
    MatrixXf X = (MatrixXf::Random(10,10));
    cout << "\n X before pca\n" << X << endl;
    
    MatrixXf X_reduced;

    if(pca(X, X_reduced, 0.99))
    {
        cout << "\n X after pca \n" << X_reduced << endl;
    }
    return 0;
}

 

cmake make 之后，运行结果为：

因为k从0计数，所以k=6 对应7列。X after pca 之后是10 行 7 列。