主成分分析(PCA)原理与实现

时间  2019-12-10
标签 成分 分析 pca 原理 实现 繁體版
原文   原文链接

主成分分析原理与实现

  主成分分析是一种矩阵的压缩算法，在减小矩阵维数的同时尽量的保留原矩阵的信息，简单来讲就是将 \(n×m\)的矩阵转换成\(n×k\)的矩阵，仅保留矩阵中所存在的主要特性，从而能够大大节省空间和数据量。最近课上学到这个知识，感受颇有意思，就在网上找一些博客进行学习，发现网上关于这方面的介绍不少，可是感受都不太全面，单靠某一个介绍仍是没法理解，固然这可能也跟我的基础有关。因此我在这里根据本身的理解写一个总结性的帖子，与你们分享同时也方便本身复习。对于主成分分析，能够参照如下几篇博客：html

  1. PCA的数学原理该博客介绍了主成分中的数学原理，给出了比较清晰的数学解释。简单易懂，可是有一些细节并无涉及到，因此仍是不能彻底理解。ios

  2. PCA 原理：为何用协方差矩阵介绍了为何在降维的时候采用协方差矩阵，可是对于协方差矩阵的解释不详细。算法

  3. 关于协方差矩阵的理解对协方差矩阵的进行了详细的推导，解释了为何能够经过\(A A^T\)来计算协方差矩阵。网络

  4. 矩阵求导、几种重要的矩阵及经常使用的矩阵求导公式对矩阵求导进行了介绍。提到了可能会用到的一些求导公式。性能

  5. UFLDL 教程学习笔记（四）主成分分析对主成分的原理和使用进行了介绍。学习

1. 数学原理

  数学原理的介绍部分能够参考文献1，该博客对主成分分析的数学原理进行了很直观的介绍。这里我根据本身的理解进行简单介绍。网站

基变换

图一（图片来源于文献 1）

  对于一个坐标点\((3,2)\)，咱们知道其表明的意思是在二维坐标里其横坐标为3，纵坐标为2。其实这隐含了一个假设，即其横纵坐标的基为\((1,0)和(0,1)\)。对于通常的二维向量，这彷佛是你们的默认状况，就像随便给出一个数字\(10\)，你们会认为这是\(10\)进制表示，除非特殊标明，不会把它看成其余进制来理解。对于任意一个坐标点\((x,y)\)，咱们能够将其表示为：
\[ \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix} \cdot \begin{pmatrix} x \\ y \\ \end{pmatrix} = \begin{pmatrix} x \\ y \\ \end{pmatrix} \]
其中\(\begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix}\)的每个行向量表明一个基向量。spa

若是我想更换基向量怎么办呢，如上图所示，若是我想知道\((3,2)\)在\((\sqrt{2}/2,\sqrt{2}/2)与(-\sqrt{2}/2,\sqrt{2}/2)\)基下的坐标值，该如何计算呢？回顾基本的数学知识，咱们发现对于一个向量在一个基上的值其实就是该向量在该基向量上的投影。因此，已知基向量，咱们能够很容易求得，对于一个向量，如\((3,2)\)，其在基\((\sqrt{2}/2,\sqrt{2}/2)与(-\sqrt{2}/2,\sqrt{2}/2)\)上的投影为：
\[ \begin{pmatrix} \sqrt{2}/2 & \sqrt{2}/2 \\ -\sqrt{2}/2 & \sqrt{2}/2 \\ \end{pmatrix} \cdot \begin{pmatrix} 3 \\ 2 \\ \end{pmatrix} = \begin{pmatrix} 5\sqrt{2}/2 \\ -\sqrt{2}/2 \\ \end{pmatrix} \]
直观的图表示如上图所示。.net

  再回到主成分分析上来，若是咱们想对一个矩阵\(A\)进行降维，其中\(A\)为：
\[ \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \\ \end{pmatrix} \]
行向量表明样本，列向量表明特征，因此其矩阵含义为m个具备n个特征的样本值。对于每个样本具备的n个特征值，其特征值之间可能会存在很大的耦合，就如文献1中所列举的那样，特征M表明是否为男性，特征F表明是否为女性，由于一我的的性别只能为其中的一个（不考虑特殊状况）。因此这两个特征只留一个就好了，因此就能够省下一半的空间。这个例子有些极端，可是并不影响理解。code

图二（图片来源于网络）

  一样对于一个具备n个特征的集合来讲，很难说这n个特征都是彻底有必要的，因此咱们就想办法来精简一些特征。选取少于n个的基向量组，将数据投影在这个向量组上，减小空间的同时又能保证信息量。首先须要明确的一点是什么才算好的基向量？首先举一个将二维空间的数据投影到一维空间的状况。如上图所示，对于空间中的这些点，咱们应该怎么投影才可以尽量的保持数据的信息量呢？经过上图中能够看出，若是将数据投影到PC1上，那么全部的数据点较为分散，与之相反，若是投影到PC2上，则数据较为集中。考虑一个极端的状况，假如全部的点在投影以后所有集中在一个点上，这样好吗？固然不！若是全部的点都集中到一个点上，那就说明全部的点都没有差异，信息所有丢失了。因此咱们但愿当数据点投影到某个坐标轴之上之后，数据越分散越好，而衡量一组数据是否发散刚好有一个统计名词“方差”，也就是说投影事后的点值方差越大越好。同时，若是数据被投影到多个基向量上，那么咱们但愿这些基向量之间的耦合程度越小越好，也就说基向量之间应该是正交的，如图三所示（建议点击连接去相应网站查看3D演示）。由于若是不考虑基向量之间的正交性，只考虑方差最大的话，那么所求得的值其实都是同样的。关于在不一样的基向量上的投影的线性相关度也有一个度量标准--协方差。那么咱们的目标明确了，使得相同特征之间方差越大越好，不一样特征之间协方差越小越好。

PCA图示

图三（参考文献【6】）

  那么这些方差，协方差什么的怎么计算呢？这里能够先给出一个结论，将\(A\)向量的每一列减去该列的平均值获得一个新的\(A\)矩阵。而后计算\(Cov=1/m \cdot A^T\cdot A\)，获得一个\(n\times n\)的矩阵\(Cov\)，那么\(Cov\)的对角线上的元素\(c_{ii}\)即为第i个特征的方差，对于其余元素\(c_{ij}\)表示第i个和第j个特征的协方差，很明显该矩阵是对称矩阵。关于该矩阵的求解方式能够参考文献3，其介绍的很详细，这里就再也不重复。须要注意的一点是这里\(Cov=1/m \cdot A^T\cdot A\)是由于A矩阵的列向量为特征，因此才这样计算。若是A矩阵的行列向量所表达的含义相反则\(Cov=1/m \cdot A\cdot A^T\)。

  已经知道了计算协方差矩阵的方法，下面看一下怎么跟咱们要作的结合在一块儿。再次总结一下咱们要作的是什么，对于一个已有的矩阵\(A\)，咱们但愿将它投影在一组新的基空间上，使之矩阵大小获得压缩。即：
\[ D_{m,N} = A_{mn} \cdot P_{nN}, \,\,\,\,given (N<n) \]

咱们要作的就是将n个特征压缩为N个特征。对于压缩过的数据投影，根据上面的叙述可知，咱们但愿对于相同特征之间方差越大越好，不一样特征之间协方差越小越好，而且咱们已经知道该如何计算方差和协方差了。
\[ Cov(D)_{NN} = D^T \cdot D = P^T A^T A P. \]
因此如今的目标很明确，咱们要作的就是求得\(P\)，使得\(Cov(D)\)的对角线元素尽量大，非对角线元素尽量小。学过线性代数的应该都知道，对于\(A^T A\)矩阵来讲，其特征向量就知足这一条件。由于已知\(A^T A\)矩阵为对称矩阵，因此可知：
\[ P^T (A^T A) P = P^{-1} (A^T A) P = \Lambda \]
其中\(\Lambda\)为\(A^T A\)的特征值组成的对角阵，\(P\)为相应的的特征向量组。

  至此，咱们就找到了进行主成分分析的方法：

  1. 首先对矩阵A进行处理，使得其每一列（或者行）减去其相应列的平均值，使得每一列的平均值都为0，而后计算\(B = A^TA\)。
  2. 求B矩阵的特征值和特征向量，将特征值进行排序，并选取前N大的特征值，选取其对应的特征向量组成特征向量组\(P_{nN}\)。
  3. \(D_{m,N} = A_{mn} \cdot P_{nN}\)即为最终想要获得的值。

2.实验验证

  下面咱们对该算法进行实际的实现，为了更好的了解PCA的工做原理，同时又保证程序的计算速度，我才用了C语言进行实现，并借助OpenBLAS库进行高效的矩阵运算。OpenBLAS是BLAS标准的一个开源实现，听说也是目前性能和维护的最好的一个。BLAS是Basic Linear Algebra Subprograms的简称，是一个矩阵运算的接口标准。既然是接口标准，那么全部根据该标准的实现都具备相同的使用方式和功能。类似的实现还有BLAS、MKL、ACML等，我使用OpenBLAS进行实现，由于其实现不依赖于任何平台，具备良好的性能，并且亲测易于安装。下面将附上个人实现代码：

//矩阵运算部分 Matrix.cpp
#include<iostream>
#include<stdio.h>
#include<stdlib.h>
//#include "mkl.h"
#include"OpenBLAS/cblas.h"
class Matrix
{
    public:
    //Print matrix;
    bool printMatrix() const;
    //get r.
    int getr() {return r;}
    //get l.
    int getc() {return c;}
    //get a.
    float *geta() {return a;}
    //normalization.
    void nmlt();
    //Compute Coevariance of a, aTxa
    void coev(Matrix &c);
    //Default constructor.
    Matrix():a(NULL), r(0), c(0) {}
    //Constructor with matrix pointer and dimension.
    Matrix(float *aa, int rr, int cc): a(aa), r(rr), c(cc) {}
    //Constructor with only dimension, should allocate space.
    Matrix(int rr, int cc): r(rr), c(cc)
    {
        a = new float[rr*cc];
    }
    //Destructor.
    ~Matrix() {delete []a; a=NULL;}

    protected:
    //Matrix pointer.
    float *a;
    //Dimension n, order lda
    int r,c;
};

extern bool printArray(float *p, int n);

class SquareMatrix:public Matrix
{
    public:
    //Default constructor.
    SquareMatrix(float *aa, int nn):Matrix(aa, nn, nn), n(nn) {}
    SquareMatrix(int nn): Matrix(nn, nn), n(nn){}
    //Destructor.
    ~SquareMatrix() {}
    //Get eigenvalue and eigenvector;
    int ssyevd(float *w);

    private:
    int n;
};
bool Matrix::printMatrix() const
{
    int i=0, j=0;
    float temp(0);
    for(i=0; i<r; i++)
    {
        for(j=0; j<c; j++)
        {
            temp = *(a+c*i+j);
            printf("%7.3f\t", temp);
        }
        std::cout<<std::endl;
    }
}


int SquareMatrix::ssyevd(float *w)
{
    lapack_int res = 0;
    res = LAPACKE_ssyevd(LAPACK_ROW_MAJOR, 'V', 'U', n, a, n, w);
    if(res == 0)
    {
        return res;
    }
    else
    {
        std::cout<<"ERROR:"<<res<<std::endl;
        exit(-1);
    }
}

void Matrix::coev(Matrix &cc)
{
    nmlt();
    cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, c, c, r, 1.0/r, a, c, a, c, 0.0, cc.geta(), c);
}

void Matrix::nmlt()
{
    int i=0,j=0;
    float av = 0.0;
    for(i=0;i<c;i++)
    {
        av = 0.0;
        for(j=0;j<r;j++)
        {
            av+=*(a+i+j*c);
        }
        av = av/r;
        for(j=0;j<r;j++)
        {
            *(a+i+j*c) -= av;
        }
    }
}

bool printArray(float *p, int n)
{
    for(int i=0; i<n; i++)
    {
        printf("%7.3f\t", p[i]);
    }
    std::cout<<std::endl;
    return true;
}
//PCA部分 PCA.cpp
#include<iostream>
#include<stdio.h>
#include<stdlib.h>
//#include "mkl.h"
#include"OpenBLAS/cblas.h"
#include"Matrix.h"
#include"PCA.h"

#define N 5
#define T 0.8f
const char SEP = ',';

static unsigned int R = 5;
static unsigned int C = 5;

int main(int argc, char *argv[])
{
    // float *A = new float [N*N]
    // {
    //  1.96f,  -6.49f,  -0.47f,  -7.20f,  -0.65f,
    // -6.49f,   3.80f,  -6.39f,   1.50f,  -6.34f,
    // -0.47f,  -6.39f,   4.17f,  -1.51f,   2.67f,
    // -7.20f,   1.50f,  -1.51f,   5.70f,   1.80f,
    // -0.65f,  -6.34f,   2.67f,   1.80f,  -7.10f
    // };
    if(argc <= 1)
    {
        printf("Usage: PCA [INPUT FILE] [OUTPUT FILE] [ROW] [COLUM]\n");
        printf("INPUT FILE: input file path.\n");
        printf("OUTPUT FILE: output file path.\n");
        printf("ROW: Row of matrix.\n");
        printf("COLUM: Colum of matrix.\n");
        exit(0);
    }
    FILE *input = fopen(argv[1], "r");
    FILE *output = fopen(argv[2], "w+");
    R = atof(argv[3]);
    C = atof(argv[4]);
    printf("Input:%s\nOutput:%s\nR:%d\nC:%d\n",argv[1], argv[2], R, C);
    float *I = new float[R*C]();
    //float *O = new float[R*C]();
    char *label = new char[R];
    //read matrix.
    readMtx(input, I, label);

    SquareMatrix cov = SquareMatrix(C);
    float *eValue = new float[C]();
    Matrix m = Matrix(I, R, C);
    Matrix n = Matrix(R, C);
    // m.printMatrix();
    //compute coveriance matrix.
    m.coev(cov);
    //compute eigenvalue and eigenvector of coveriance matrix.
    cov.ssyevd(eValue);
    //Compute compressed matrix.
    eMtx(m, cov, n);
    //n.printMatrix();
    saveMtx(output, n.geta(), label);

    fclose(input);
    fclose(output);
    delete []label;
    delete []eValue;
    return 0;
}

//eigen matrix
void eMtx(Matrix&a, Matrix&b, Matrix&r)
{
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, a.getr(), b.getc(), a.getc(), 1.0, a.geta(), a.getc(), b.geta(), b.getc(), 0.0, r.geta(), b.getc());
}

bool readUtl(FILE *f, char sep)
{
    char c;
    if((c=fgetc(f))!=EOF && c==sep)
    {
        return true;
    }
    return false;
}

void readMtx(FILE *f, float *m, char *la)
{
    float ft(0.0);
    char ch;
    int i(0),j(0),index(0);
    while(i<R)
    {
        while(!readUtl(f, SEP));
        la[i++] = fgetc(f);
        readUtl(f, SEP);
        for(j=0;j<C-1;j++)
        {
            fscanf(f, "%f,", &m[index++]);
        }
        fscanf(f, "%f", &m[index++]);
        while(!readUtl(f, '\n') && i<R);
    }
}

void saveMtx(FILE *f, float *m, char *la)
{
    int i(0),j(0);
    for(i=0;i<R;i++)
    {
        fprintf(f, "%c,", la[i]);
        for(j=0;j<C-1;j++)
        {
            fprintf(f, "%.4f,", m[i*C+j]);
        }
        fprintf(f, "%.4f", m[i*C+j]);
        fprintf(f, "\n");
    }
}

编译运行：

./PCA wdbc.data wdbc.out 569 30

本文所采用的实验数据为开源数据集，该数据集是有关于乳腺癌诊断的相关数据，共有569条记录，每个记录有30个特征，而且每一条记录都有一个标签，标签为'B'意味着良性，'M'意味着恶性。上述代码对该数据集继续主成分分析，最后将输出矩阵保存在wdbc.out中。
下面我经过散点图的方式直观的展现分析的效果：
一维图

PCA一维映射

其中绿色表明良性，红色表明恶性。从图中能够看出，即便仅映射到一维，不一样类别的数据彷佛就已经很容易分离开了，这是由于咱们选取的这个一维空间正是最大的那个特征值对应的空间，因此包含最多的信息。接下来咱们将数据映射到二维和三维空间：

二维图

PCA二维映射

PCA三维映射

参考文献

[1]http://blog.codinglabs.org/articles/pca-tutorial.html

[2]http://www.javashuo.com/article/p-rsmsmslm-nt.html

[3]https://blog.csdn.net/itplus/article/details/11452743#commentsedit

[4]http://www.javashuo.com/article/p-sioenphr-ep.html

[5]https://blog.csdn.net/itplus/article/details/11451327

[6]http://setosa.io/ev/principal-component-analysis/

相关文章
相关标签/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公众号
   欢迎关注本站公众号,获取更多信息
联系我们 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程