【转】浅谈对主成分分析（PCA）算法的理解

以前对PCA算法有过一段时间的研究，但没整理成文章，最近项目又打算用到PCA算法，故趁热打铁整理下PCA算法的知识。本文观点旨在抛砖引玉，不是权威，更不能尽信，只是本人的一点体会。

主成分分析（PCA）是多元统计分析中用来分析数据的一种方法，它是用一种较少数量的特征对样本进行描述以达到降低特征空间维数的方法，它的本质实际上是K-L变换。PCA方法最著名的应用应该是在人脸识别中特征提取及数据维，我们知道输入200*200大小的人脸图像，单单提取它的灰度值作为原始特征，则这个原始特征将达到40000维，这给后面分类器的处理将带来极大的难度。著名的人脸识别Eigenface算法就是采用PCA算法，用一个低维子空间描述人脸图像，同时用保存了识别所需要的信息。下面先介绍下PCA算法的本质K-L变换。

1、K-L变换（卡洛南-洛伊（Karhunen-Loeve）变换）：最优正交变换  

• 一种常用的特征提取方法；
• 最小均方误差意义下的最优正交变换；
• 在消除模式特征之间的相关性、突出差异性方面有最优的效果。

离散K-L变换：对向量x（可以想象成 M维=width*height 的人脸图像原始特征）用确定的完备正交归一向量系u_j展开：

 

 

 

这个公式由来我想应该是任一n维欧式空间V均存在正交基，利用施密特正交化过程即可构建这个正交基。

现在我们希望用d个有限项来估计向量x，公式如下：

 

计算该估计的均方误差如下：

 

要使用均方误差最小，我们采用Langrange乘子法进行求解：

 

                                                           

 

因此，当满足上式时，

取得最小值。

 

 

即相关矩阵R的d个特征向量（对应d个特征值从大到小排列）为基向量来展开向量x时，其均方误差最小，为：

 

因此，K-L变换定义：当取矩阵R的d个最大特征值对应的特征向量来展开x时，其截断均方误差最小。这d个特征向量组成的正交坐标系称作x所在的D维空间的d维K-L变换坐标系， x在K-L坐标系上的展开系数向量y称作x的K-L变换。

 

---

总结下，K-L变换的方法：对相关矩阵R的特征值由大到小进行排队，

则均方误差最小的x近似于：

                                                               

矩阵形式：

 

上式两边乘以U的转置，得

                                             

 

向量y就是变换（降维）后的系数向量，在人脸识别Eigenface算法中就是用系数向量y代替原始特征向量x进行识别。

---

 

下面，我们来看看相关矩阵R到底是什么样子。

 

因此，我们可以看出相关矩阵R是一个实对称矩阵（或者严谨的讲叫正规矩阵），正规矩阵有什么特点呢？？学过《矩阵分析》的朋友应该知道：

若矩阵R是一个实对称矩阵，则必定存在正交矩阵U，使得R相似于对角形矩阵，即：

 

 因此，我们可以得出这样一个结论：

                                                     

 

 降维后的系数向量y的相关矩阵是对角矩阵，即通过K-L变换消除原有向量x的各分量间的相关性，从而有可能去掉那些带有较少信息的分量以达到降低特征维数的目的。

 

2、主成分分析（PCA)

 

主成分分析（PCA）的原理就是将一个高维向量x,通过一个特殊的特征向量矩阵U，投影到一个低维的向量空间中，表征为一个低维向量y，并且仅仅损失了一些次要信息。也就是说，通过低维表征的向量和特征向量矩阵，可以基本重构出所对应的原始高维向量。

在人脸识别中，特征向量矩阵U称为特征脸（eigenface）空间，因此其中的特征向量u_i进行量化后可以看出人脸轮廓，在下面的实验中可以看出。

以人脸识别为例，说明下PCA的应用。

设有N个人脸训练样本，每个样本由其像素灰度值组成一个向量x_i，则样本图像的像素点数即为x_i的维数，M=width*height ，由向量构成的训练样本集为。

该样本集的平均向量为：

平均向量又叫平均脸。

 

样本集的协方差矩阵为：

 

求出协方差矩阵的特征向量u_i和对应的特征值，这些特征向量组成的矩阵U就是人脸空间的正交基底，用它们的线性组合可以重构出样本中任意的人脸图像，（如果有朋友不太理解这句话的意思，请看下面的总结2。）并且图像信息集中在特征值大的特征向量中，即使丢弃特征值小的向量也不会影响图像质量。

将协方差矩阵的特征值按大到小排序：。由大于的对应的特征向量构成主成分，主成分构成的变换矩阵为：

                               

这样每一幅人脸图像都可以投影到构成的特征脸子空间中，U的维数为M×d。有了这样一个降维的子空间，任何一幅人脸图像都可以向其作投影，即并获得一组坐标系数，即低维向量y，维数d×1,为称为KL分解系数。这组系数表明了图像在子空间的位置，从而可以作为人脸识别的依据。

---

有朋友可能不太理解，第一部分讲K-L变换的时候，求的是相关矩阵的特征向量和特征值，这里怎么求的是协方差矩阵?

其实协方差矩阵也是：

，可以看出其实用代替x就成了相关矩阵R，相当于原始样本向量都减去个平均向量，实质上还是一样的，协方差矩阵也是实对称矩阵。

---

总结下：

1、在人脸识别过程中，对输入的一个测试样本x，求出它与平均脸的偏差，则在特征脸空间U的投影，可以表示为系数向量y：

                                     

U的维数为M×d，的维数为M×1，y的维数d×1。若M为200*200=40000维，取200个主成分，即200个特征向量，则最后投影的系数向量y维数降维200维。

2、根据1中的式子，可以得出：

                                

这里的x就是根据投影系数向量y重构出的人脸图像，丢失了部分图像信息，但不会影响图像质量。

 

 3、PCA算法实验  

 

在计算机视觉库OpenCV中较新的版本中，封装了PCA算法的类。下面是对PCA算法做的一些实验，有助于加深对PCA算法的理解。代码来自于tornadomeet，我并没有作太多修改，加多一些说明。

链接在此：http://www.cnblogs.com/tornadomeet/archive/2012/09/06/2673104.html

运行环境为：WindowsXP+QT+OpenCV2.3.1。

 

部分函数说明如下：

Mat Mat::reshape(int cn, int rows=0) const

　　该函数是改变Mat的尺寸，即保持尺寸大小=行数*列数*通道数 不变。其中第一个参数为变换后Mat的通道数，如果为0，代表变换前后通道数不变。第二个参数为变换后Mat的行数，如果为0也是代表变换前后通道数不变。但是该函数本身不复制数据。

 

　　void Mat::convertTo(OutputArray m, int rtype, double alpha=1, double beta=0 ) const

　　该函数其实是对原Mat的每一个值做一个线性变换。参数1为目的矩阵，参数2为目d矩阵的类型，参数3和4变换的系数，看完下面的公式就明白了：

　　

 

　　PCA::PCA(InputArray data, InputArray mean, int flags, int maxComponents=0)

　　该构造函数的第一个参数为要进行PCA变换的输入Mat；参数2为该Mat的均值向量；参数3为输入矩阵数据的存储方式，如果其值为CV_PCA_DATA_AS_ROW则说明输入Mat的每一行代表一个样本，同理当其值为CV_PCA_DATA_AS_COL时，代表输入矩阵的每一列为一个样本；最后一个参数为该PCA计算时保留的最大主成分的个数。如果是缺省值，则表示所有的成分都保留。

 

　　Mat PCA::project(InputArray vec) const

　　该函数的作用是将输入数据vec(该数据是用来提取PCA特征的原始数据)投影到PCA主成分空间中去，返回每一个样本主成分特征组成的矩阵。因为经过PCA处理后，原始数据的维数降低了，因此原始数据集中的每一个样本的维数都变了，由改变后的样本集就组成了本函数的返回值。下面由一个图说明：

                                                              

 

　　Mat PCA::backProject(InputArray vec) const

　　一般调用backProject（）函数前需调用project()函数，因为backProject()函数的参数vec就是经过PCA投影降维过后的矩阵dst。 因此backProject()函数的作用就是用vec来重构原始数据集（关于该函数的本质就是上面总结2的公式）。由一个图说明如下：

                                                                 

　　另外PCA类中还有几个成员变量，mean,eigenvectors, eigenvalues等分别对应着原始数据的均值，协方差矩阵的特征值和特征向量。

实验结果：

实验是用4个人人脸图像，其中每个人分别有5张，共计20张人脸图片。用这些图片组成原始数据集来提取他们的PCA主特征脸。该20张图片如下所示：

                                                   

软件运行结果：

实验中保留4个特征向量作为人脸图像的正交基底，运行结果如下：

                                                              

   其中第一行的3张人脸分别为20张原图中的3张，这里取的是3个不同人的。

　　第二行中显示的3张人脸重构的人脸图像，可以看出由于只取了4个特征向量作为正交基底，因此重构后的人脸图像一些细节会丢失。如果增加保留的特征向量个数，则能较好的重构出人脸图像。

　　第3行的人脸图为取的原始数据协方差矩阵特征向量的最前面3个，因此这3个人脸为最具代表人脸特征的3个PCA人脸特征。

实验主要部分代码：

pcaface.h

 1 #ifndef PCAFACE_H
 2 #define PCAFACE_H
 3 #include <opencv2/core/core.hpp>
 4 #include <opencv2/highgui/highgui.hpp>
 5 #include <opencv2/imgproc/imgproc.hpp>
 6 
 7 using namespace cv;
 8 
 9 #include <QDialog>
10 
11 namespace Ui {
12 class PCAFace;
13 }
14 
15 class PCAFace : public QDialog
16 {
17     Q_OBJECT
18     
19 public:
20     explicit PCAFace(QWidget *parent = 0);
21     ~PCAFace();
22 
23     Mat normalize(const Mat& src);
24 
25     
26 protected:
27     void changeEvent(QEvent *e);
28     
29 private slots:
30     void on_startButton_clicked();
31 
32     void on_closeButton_clicked();
33 
34 private:
35     Ui::PCAFace *ui;
36     Mat src_face1, src_face2, src_face3;
37     Mat project_face1, project_face2, project_face3;
38     Mat dst;
39     Mat pca_face1, pca_face2, pca_face3;
40     vector<Mat> src;
41     int total;
42 };
43 
44 #endif // PCAFACE_H

pcaface.cpp

  1 #include "pcaface.h"
  2 #include "ui_pcaface.h"
  3 #include <QString>
  4 #include <iostream>
  5 #include <stdio.h>
  6 
  7 using namespace std;
  8 
  9 PCAFace::PCAFace(QWidget *parent) :
 10     QDialog(parent),
 11     ui(new Ui::PCAFace)
 12 {
 13     ui->setupUi(this);
 14     src_face1 = imread("./images/1.pgm", 0);
 15 
 16     //下面的代码为设置图片显示区域自适应图片的大小
 17     ui->face1Browser->setFixedHeight(src_face1.rows+1);
 18     ui->face1Browser->setFixedWidth(src_face1.cols+1);
 19     ui->face2Browser->setFixedHeight(src_face1.rows+1);
 20     ui->face2Browser->setFixedWidth(src_face1.cols+1);
 21     ui->face3Browser->setFixedHeight(src_face1.rows+1);
 22     ui->face3Browser->setFixedWidth(src_face1.cols+1);
 23 
 24     ui->face4Browser->setFixedHeight(src_face1.rows+1);
 25     ui->face4Browser->setFixedWidth(src_face1.cols+1);
 26     ui->face5Browser->setFixedHeight(src_face1.rows+1);
 27     ui->face5Browser->setFixedWidth(src_face1.cols+1);
 28     ui->face6Browser->setFixedHeight(src_face1.rows+1);
 29     ui->face6Browser->setFixedWidth(src_face1.cols+1);
 30 
 31     ui->face7Browser->setFixedHeight(src_face1.rows+1);
 32     ui->face7Browser->setFixedWidth(src_face1.cols+1);
 33     ui->face8Browser->setFixedHeight(src_face1.rows+1);
 34     ui->face8Browser->setFixedWidth(src_face1.cols+1);
 35     ui->face9Browser->setFixedHeight(src_face1.rows+1);
 36     ui->face9Browser->setFixedWidth(src_face1.cols+1);
 37 
 38     for(int i = 1; i <= 15; i++)
 39     {
 40         stringstream ss;
 41         string num;
 42         ss<<i;        //将整数i读入字符串流
 43         ss>>num;      //将字符串流中的数据传入num，这2句代码即把数字转换成字符
 44         string image_name = ("./images/" + num + ".pgm");//需要读取的图片全名
 45         src.push_back(imread(image_name, 0));
 46     }
 47     total= src[0].rows*src[0].cols;
 48 }
 49 
 50 PCAFace::~PCAFace()
 51 {
 52     delete ui;
 53 }
 54 
 55 void PCAFace::changeEvent(QEvent *e)
 56 {
 57     QDialog::changeEvent(e);
 58     switch (e->type()) {
 59     case QEvent::LanguageChange:
 60         ui->retranslateUi(this);
 61         break;
 62     default:
 63         break;
 64     }
 65 }
 66 
 67 //将Mat内的内容归一化到0～255，归一化后的类型为但通道整型
 68 Mat PCAFace::normalize(const Mat& src) {
 69     Mat srcnorm;
 70     cv::normalize(src, srcnorm, 0, 255, NORM_MINMAX, CV_8UC1);
 71     return srcnorm;
 72 }
 73 
 74 
 75 void PCAFace::on_startButton_clicked()
 76 {
 77     //先显示3张原图
 78     ui->face1Browser->append("<img src=./images/1.pgm>");
 79     ui->face2Browser->append("<img src=./images/7.pgm>");
 80     ui->face3Browser->append("<img src=./images/14.pgm>");
 81 
 82     //mat数组用来存放读取进来的所有图片的数据，其中mat的每一列对应1张图片，该实现在下面的for函数中
 83     Mat mat(total, src.size(), CV_32FC1);
 84     for(int i = 0; i < src.size(); i++)
 85     {
 86         Mat col_tmp = mat.col(i);
 87         src[i].reshape(1, total).col(0).convertTo(col_tmp, CV_32FC1, 1/255.);
 88     }
 89     int number_principal_compent = 4;   //保留最大的主成分数
 90 
 91     //构造pca数据结构
 92     PCA pca(mat, Mat(), CV_PCA_DATA_AS_COL, number_principal_compent);
 93 
 94 
 95     //pca.eigenvectors中的每一行代表输入数据协方差矩阵一个特征向量，且是按照该协方差矩阵的特征值进行排序的
 96     pca_face1 = normalize(pca.eigenvectors.row(0)).reshape(1, src[0].rows);  //第一个主成分脸
 97     imwrite("./result/pca_face1.jpg", pca_face1);//显示主成分特征脸1
 98     ui->face7Browser->append("<img src=./result/pca_face1.jpg>");
 99 
100     pca_face2 = normalize(pca.eigenvectors.row(1)).reshape(1, src[0].rows);  //第二个主成分脸
101     imwrite("./result/pca_face2.jpg", pca_face2);//显示主成分特征脸2
102     ui->face8Browser->append("<img src=./result/pca_face2.jpg>");
103 
104     pca_face3 = normalize(pca.eigenvectors.row(2)).reshape(1, src[0].rows);  //第三个主成分脸
105     imwrite("./result/pca_face3.jpg", pca_face3);//显示主成分特征脸3
106     ui->face9Browser->append("<img src=./result/pca_face3.jpg>");
107 
108     //将原始数据通过PCA方向投影，即通过特征向量的前面几个作用后的数据，因此这里的dst的尺寸变小了
109     dst = pca.project(mat);
110 
111     //通过方向投影重构原始人脸图像
112     project_face1 = normalize(pca.backProject(dst).col(1)).reshape(1, src[0].rows);
113     imwrite("./result/project_face1.jpg", project_face1);
114     ui->face4Browser->append("<img src=./result/project_face1.jpg>");
115 
116     project_face2 = normalize(pca.backProject(dst).col(7)).reshape(1, src[0].rows);
117     imwrite("./result/project_face2.jpg", project_face2);
118     ui->face5Browser->append("<img src=./result/project_face2.jpg>");
119 
120     project_face3 = normalize(pca.backProject(dst).col(14)).reshape(1, src[0].rows);
121     imwrite("./result/project_face3.jpg", project_face3);
122     ui->face6Browser->append("<img src=./result/project_face3.jpg>");
123 }
124 
125 void PCAFace::on_closeButton_clicked()
126 {
127     close();
128 }

实验的工程代码可以在上面的那个链接下载，环境搭建可以参考我之前写个一个文章http://www.cnblogs.com/liu-jun/archive/2012/09/26/Jacky_Liu.html。

作者：Jacky_Liu，转载或分享请注明出处。