主成分分析原理與實現
主成分分析是一種矩陣的壓縮算法，在減少矩陣維數的同時盡可能的保留原矩陣的信息，簡單來說就是將 \(n×m\)的矩陣轉換成\(n×k\)的矩陣，僅保留矩陣中所存在的主要特性，從而可以大大節省空間和數據量。最近課上學到這個知識，感覺很有意思，就在網上找一些博客進行學習，發現網上關於這方面的介紹很多，但是感覺都不太全面，單靠某一個介紹還是無法理解，當然這可能也跟個人基礎有關。所以我在這里根據自己的理解寫一個總結性的帖子，與大家分享同時也方便自己復習。對於主成分分析，可以參照以下幾篇博客：
-
PCA的數學原理該博客介紹了主成分中的數學原理，給出了比較清晰的數學解釋。簡單易懂，但是有一些細節並沒有涉及到，所以還是不能完全理解。
-
PCA 原理：為什么用協方差矩陣介紹了為什么在降維的時候采用協方差矩陣，但是對於協方差矩陣的解釋不詳細。
-
關於協方差矩陣的理解對協方差矩陣的進行了詳細的推導，解釋了為什么可以通過\(A A^T\)來計算協方差矩陣。
-
矩陣求導、幾種重要的矩陣及常用的矩陣求導公式對矩陣求導進行了介紹。提到了可能會用到的一些求導公式。
-
UFLDL 教程學習筆記（四）主成分分析對主成分的原理和使用進行了介紹。
1. 數學原理
數學原理的介紹部分可以參考文獻1，該博客對主成分分析的數學原理進行了很直觀的介紹。這里我根據自己的理解進行簡單介紹。
對於一個坐標點\((3,2)\)，我們知道其代表的意思是在二維坐標里其橫坐標為3，縱坐標為2。其實這隱含了一個假設，即其橫縱坐標的基為\((1,0)和(0,1)\)。對於一般的二維向量，這似乎是大家的默認情況，就像隨便給出一個數字\(10\)，大家會認為這是\(10\)進制表示，除非特殊標明，不會把它當作其他進制來理解。對於任意一個坐標點\((x,y)\)，我們可以將其表示為：
其中\(\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)\)上的投影為：
直觀的圖表示如上圖所示。
再回到主成分分析上來，如果我們想對一個矩陣\(A\)進行降維，其中\(A\)為：
行向量代表樣本，列向量代表特征，所以其矩陣含義為m個具有n個特征的樣本值。對於每一個樣本具有的n個特征值，其特征值之間可能會存在很大的耦合，就如文獻1中所列舉的那樣，特征M代表是否為男性，特征F代表是否為女性，因為一個人的性別只能為其中的一個（不考慮特殊情況）。所以這兩個特征只留一個就行了，所以就可以省下一半的空間。這個例子有些極端，但是並不影響理解。
同樣對於一個具有n個特征的集合來說，很難說這n個特征都是完全有必要的，所以我們就想辦法來精簡一些特征。選取少於n個的基向量組，將數據投影在這個向量組上，減少空間的同時又能保證信息量。首先需要明確的一點是什么才算好的基向量？首先舉一個將二維空間的數據投影到一維空間的情況。如上圖所示，對於空間中的這些點，我們應該怎么投影才能夠盡可能的保持數據的信息量呢？通過上圖中可以看出，如果將數據投影到PC1上，那么所有的數據點較為分散，與之相反，如果投影到PC2上，則數據較為集中。考慮一個極端的情況，假如所有的點在投影之后全部集中在一個點上，這樣好嗎？當然不！如果所有的點都集中到一個點上，那就說明所有的點都沒有差別，信息全部丟失了。所以我們希望當數據點投影到某個坐標軸之上以后，數據越分散越好，而衡量一組數據是否發散恰好有一個統計名詞“方差”，也就是說投影過后的點值方差越大越好。同時，如果數據被投影到多個基向量上，那么我們希望這些基向量之間的耦合程度越小越好，也就說基向量之間應該是正交的，如圖三所示（建議點擊鏈接去相應網站查看3D演示）。因為如果不考慮基向量之間的正交性，只考慮方差最大的話，那么所求得的值其實都是一樣的。關於在不同的基向量上的投影的線性相關度也有一個度量標准--協方差。那么我們的目標明確了，使得相同特征之間方差越大越好，不同特征之間協方差越小越好。
那么這些方差，協方差什么的怎么計算呢？這里可以先給出一個結論，將\(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\)，我們希望將它投影在一組新的基空間上，使之矩陣大小得到壓縮。即：
我們要做的就是將n個特征壓縮為N個特征。對於壓縮過的數據投影，根據上面的敘述可知，我們希望對於相同特征之間方差越大越好，不同特征之間協方差越小越好，並且我們已經知道該如何計算方差和協方差了。
所以現在的目標很明確，我們要做的就是求得\(P\)，使得\(Cov(D)\)的對角線元素盡可能大，非對角線元素盡可能小。學過線性代數的應該都知道，對於\(A^T A\)矩陣來說，其特征向量就滿足這一條件。因為已知\(A^T A\)矩陣為對稱矩陣，所以可知：
其中\(\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中。
下面我通過散點圖的方式直觀的展示分析的效果：
參考文獻
[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