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

主成分分析原理与实现

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

  1. PCA的数学原理该博客介绍了主成分中的数学原理，给出了比较清晰的数学解释。简单易懂，但是有一些细节并没有涉及到，所以还是不能完全理解。
  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}\)的每一个行向量代表一个基向量。

如果我想更换基向量怎么办呢，如上图所示，如果我想知道\((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}
\]
直观的图表示如上图所示。

??再回到主成分分析上来，如果我们想对一个矩阵\(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代表是否为女性，因为一个人的性别只能为其中的一个（不考虑特殊情况）。所以这两个特征只留一个就行了，所以就可以省下一半的空间。这个例子有些极端，但是并不影响理解。

图二（图片来源于网络）

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


参考文献【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]https://blog.csdn.net/a10767891/article/details/80288463

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

[4]https://blog.csdn.net/daaikuaichuan/article/details/80620518

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

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

原文地址：https://www.cnblogs.com/xinyuyang/p/11178676.html

时间： 2024-10-11 06:14:47

主成分分析(PCA)原理与实现的相关文章

主成分分析法原理及其python实现

主成分分析法原理及其python实现 前言: 这片文章主要参考了Andrew Ng的Machine Learning课程讲义,我进行了翻译,并配上了一个python演示demo加深理解. 本文主要介绍一种降维算法,主成分分析法,Principal Components Analysis,简称PCA,这种方法的目标是找到一个数据近似集中的子空间,至于如何找到这个子空间,下文会给出详细的介绍,PCA比其他降维算法更加直接,只需要进行一次特征向量的计算即可.(在Matlab,python,R中这个可以

PCA原理分析

动机 在机器学习领域中,我们常常会遇到维数很高的数据,有些数据的特征维度高达上百万维,很显然这样的数据是无法直接计算的,而且维度这么高,其中包含的信息一定有冗余,这时就需要进行降维,总的来说,我们降维的主要目的有如下几条: 在原始的高维空间中,包含有冗余信息以及噪音信息,在实际应用例如图像识别中造成了误差,降低了准确率:而通过降维,我们希望减少冗余信息所造成的误差,提高识别(或其他应用)的精度. 希望通过降维算法来寻找数据内部的本质结构特征. 通过降维来加速后续计算的速度 还有其他很多目的,如解

机器学习中的数学(4)-线性判别分析（LDA）, 主成分分析(PCA)

机器学习中的数学(4)-线性判别分析(LDA), 主成分分析(PCA) 版权声明: 本文由LeftNotEasy发布于http://leftnoteasy.cnblogs.com, 本文可以被全部的转载或者部分使用,但请注明出处,如果有问题,请联系[email protected] 前言: 第二篇的文章中谈到,和部门老大一宁出去outing的时候,他给了我相当多的机器学习的建议,里面涉及到很多的算法的意义.学习方法等等.一宁上次给我提到,如果学习分类算法,最好从线性的入手,线性分类器最简单的就是

主成分分析(PCA)

相对与网上很多人分享的有关PCA的经历,我第一次接触PCA却不是从人脸表情识别开始的,但我所在的实验室方向之一是人脸的研究,最后也会回到这个方向上来吧. PCA(principal components analysis)是一种非常有用的统计技术,它已经应用于人脸识别和图像压缩领域中,并且是高维数据计算模型的常用技术.简单说是把高维数据将成低维数据,比如100000x100000的矩阵降成100000x100的. 从例子中也看得出在数学模型中直观看到的是对矩阵进行的各种各样的变形最终达到我们所需

线性判别分析（LDA）, 主成分分析(PCA)及其推导【转】

前言: 如果学习分类算法,最好从线性的入手,线性分类器最简单的就是LDA,它可以看做是简化版的SVM,如果想理解SVM这种分类器,那理解LDA就是很有必要的了. 谈到LDA,就不得不谈谈PCA,PCA是一个和LDA非常相关的算法,从推导.求解.到算法最终的结果,都有着相当的相似. 本次的内容主要是以推导数学公式为主,都是从算法的物理意义出发,然后一步一步最终推导到最终的式子,LDA和PCA最终的表现都是解一个矩阵特征值的问题,但是理解了如何推导,才能更深刻的理解其中的含义.本次内容要求读者有一些

机器学习——降维（主成分分析PCA、线性判别分析LDA、奇异值分解SVD、局部线性嵌入LLE）

机器学习--降维(主成分分析PCA.线性判别分析LDA.奇异值分解SVD.局部线性嵌入LLE) 以下资料并非本人原创,因为觉得石头写的好,所以才转发备忘 (主成分分析(PCA)原理总结)[https://mp.weixin.qq.com/s/XuXK4inb9Yi-4ELCe_i0EA] 来源:?石头?机器学习算法那些事?3月1日 主成分分析(Principal components analysis,以下简称PCA)是最常用的降维方法之一,在数据压缩和消除冗余方面具有广泛的应用,本文由浅入深的

降维（一）----说说主成分分析(PCA)的源头

降维(一)----说说主成分分析(PCA)的源头 降维系列: 降维(一)----说说主成分分析(PCA)的源头 降维(二)----Laplacian Eigenmaps --------------------- 主成分分析(PCA) 在很多教程中做了介绍,但是为何通过协方差矩阵的特征值分解能够得到数据的主成分?协方差矩阵和特征值为何如此神奇,我却一直没弄清.今天终于把整个过程整理出来,方便自己学习,也和大家交流. 提出背景 以二维特征为例,两个特征之间可能存在线性关系的(例如这两个特征分别是运

K-L变换和 主成分分析PCA

一.K-L变换 说PCA的话,必须先介绍一下K-L变换了. K-L变换是Karhunen-Loeve变换的简称,是一种特殊的正交变换.它是建立在统计特性基础上的一种变换,有的文献也称其为霍特林(Hotelling)变换,因为他在1933年最先给出将离散信号变换成一串不相关系数的方法.K-L变换的突出优点是它能去相关性,而且是均方误差(Mean Square Error,MSE)意义下的最佳变换. 下面就简单的介绍一下K-L变换了. 设,随机向量X ∈Rn(n阶列向量),它的均值向量为mX,则其协

深度学习入门教程UFLDL学习实验笔记三：主成分分析PCA与白化whitening

主成分分析与白化是在做深度学习训练时最常见的两种预处理的方法,主成分分析是一种我们用的很多的降维的一种手段,通过PCA降维,我们能够有效的降低数据的维度,加快运算速度.而白化就是为了使得每个特征能有同样的方差,降低相邻像素的相关性. 主成分分析PCA PCA算法可以将输入向量转换为一个维数低很多的近似向量.我们在这里首先用2D的数据进行试验,其数据集可以在UFLDL网站的相应页面http://ufldl.stanford.edu/wiki/index.php/Exercise:PCA_in_2D