主成分分析(PCA)原理與實現

時間 2019-12-10

標籤成分分析 pca 原理實現简体版

原文原文鏈接

主成分分析原理與實現

主成分分析是一種矩陣的壓縮算法，在減小矩陣維數的同時儘量的保留原矩陣的信息，簡單來講就是將 \(n×m\)的矩陣轉換成\(n×k\)的矩陣，僅保留矩陣中所存在的主要特性，從而能夠大大節省空間和數據量。最近課上學到這個知識，感受頗有意思，就在網上找一些博客進行學習，發現網上關於這方面的介紹不少，可是感受都不太全面，單靠某一個介紹仍是沒法理解，固然這可能也跟我的基礎有關。因此我在這裏根據本身的理解寫一個總結性的帖子，與你們分享同時也方便本身複習。對於主成分分析，能夠參照如下幾篇博客：html

PCA的數學原理該博客介紹了主成分中的數學原理，給出了比較清晰的數學解釋。簡單易懂，可是有一些細節並無涉及到，因此仍是不能徹底理解。ios
PCA 原理：爲何用協方差矩陣介紹了爲何在降維的時候採用協方差矩陣，可是對於協方差矩陣的解釋不詳細。算法
關於協方差矩陣的理解對協方差矩陣的進行了詳細的推導，解釋了爲何能夠經過\(A A^T\)來計算協方差矩陣。網絡
矩陣求導、幾種重要的矩陣及經常使用的矩陣求導公式對矩陣求導進行了介紹。提到了可能會用到的一些求導公式。性能
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演示）。由於若是不考慮基向量之間的正交性，只考慮方差最大的話，那麼所求得的值其實都是同樣的。關於在不一樣的基向量上的投影的線性相關度也有一個度量標準--協方差。那麼咱們的目標明確了，使得相同特徵之間方差越大越好，不一樣特徵之間協方差越小越好。

圖三（參考文獻【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\)爲相應的的特徵向量組。

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

首先對矩陣A進行處理，使得其每一列（或者行）減去其相應列的平均值，使得每一列的平均值都爲0，而後計算\(B = A^TA\)。
求B矩陣的特徵值和特徵向量，將特徵值進行排序，並選取前N大的特徵值，選取其對應的特徵向量組成特徵向量組\(P_{nN}\)。
\(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/