从QR分解到PCA,再到人脸识别

        PCA(Principal Component Analysis,主成分分析)是一种很常用的根据变量协方差对数据进行降维、压缩的方法。它的精髓在于尽量用最少数量的维度,尽可能精确地描述数据。

        PCA对数据进行降维的过程可以用下面这个动图来解释(图片摘自http://stats.stackexchange.com/a/140579/93946):

lNHqt        在上图中,一组位于直角坐标系的二维样本集,沿着斜线的方向有很强的相关性。所以如果我们将直角坐标系转换到斜向,也就是让横轴沿斜线方向,纵轴垂直于斜线方向。于是,在这个新的坐标系下,数据点在横轴上分布很分散,但是在纵轴方向比较集中。如果在误差允许范围内,我们完全可以将数据点在新纵轴上的坐标全部置为0,只用新横轴上的坐标来表示点的位置。这样,就完成了对数据降维的过程(即将原始直角坐标系的2个维度减小到新坐标系的1个维度)。对更高维的情况,处理过程与之类似。


PCA人脸识别

        将PCA用于人脸识别的过程如下:

        1.假设有400幅尺寸为100*100的图像,构成10000*400的矩阵X=[x_1,\dots,x_n]

        2.计算均值\mu=\frac{1}{n}\sum_{j=1}^n x_j,令H=\frac{1}{\sqrt{n-1}}[x_1-\mu,\dots,x_n-\mu]

        3.根据定义,计算协方差矩阵\Sigma=HH^T

        4.计算\Sigma的特征值与特征向量,取前h个最大特征值所对应的特征向量,构成矩阵\Phi

        5.矩阵\Phi可对数据降维:\Phi^T X=Y,Y是h行400列的矩阵,也就是将数据从10000维降为h维。

        这种做法一个明显的缺陷在于,\Sigma的维度为10000×10000,直接进行奇异值分解计算量非常大。利用QR分解,作间接的奇异值分解,可以减小计算量。


利用QR分解减小计算量

        基于QR分解的PCA算法步骤如下:

        1.已知\Sigma=HH^T,其中\Sigma为d*d,H为d*n,d代表原始数据的维数,n代表样本数,d远大于n;

        2.H作QR分解,h=QR,其中Q为d*t,R为t*n,1\leq t \leq n

        3.\Sigma=QRR^T Q^T,对R^T作奇异值分解R^T=UDV^T,其中U为n*t,V为t*t,D=diag(\sigma_1,\dots,\sigma_t)

        4.于是\Sigma=QVDU^T UDV^T Q^T=QVD^2 V^T Q^T=QV\Lambda V^T Q^T,其中\Lambda=D^2

        5.由于(QV)^T (QV)=V^T Q^T QV=V^T V=I,所以QV可将\Sigma对角化,QV为\Sigma的特征向量矩阵,\Lambda\Sigma的特征值矩阵;

        6.选取D前h个最大对角元所对应于V中的h个列,构成t*h的矩阵V_h,则降维矩阵\Phi=QV_h

        该过程涉及对一个很大的矩阵的QR分解,和对一个较小矩阵的奇异值分解。计算量与传统方法相比较的结果如下(图片摘自[1]):

computation_comparision

        进一步,进行人脸识别的过程如下:

        1.假设有c个类别,每类包含s个样本,则n=c∗s;

        2.对X计算Y=\Phi^T X,将Y(也称特征脸)按类别计算均值,得到c个长度为h的列向量v_1,\dots,v_c

        3.对于未知类别的新样本x,计算y=\Phi^T x,y的长度为h;

        4.计算距离d(y,v_i),i=1,\dots,c,取距离最小的i作为x的类标号。


距离度量d

        1.欧式距离

Euclidean_Distance

        2.曼哈顿距离

Manhattan_Distance

        3.马氏距离

Mahalanobis_Distance

        在马氏距离中,x与y分布相同,且协方差矩阵为S。加入协方差矩阵的逆矩阵的作用是,将如下图(图片部分取自http://stats.stackexchange.com/a/62147/93946)中呈椭圆分布的数据归一化到圆形分布中,再来比较距离,可以抵消不同样本集在特征空间中的分布差异。

Mahalanobis_Distance_example


C++实现

        环境要求:OpenCV(样本图像的读取),Armadillo(高性能线性代数C++函数库),Intel MKL(替代LAPACK为Armadillo提供矩阵分解计算)。

        项目使用Visual Studio Ultimate 2012建立,不过核心代码只有一个cpp文件。

        全部代码托管在github.com/johnhany/QR-PCA-FaceRec


分类测试

        测试样本采用The AT&T Database of Faces,原称The ORL Database of Faces,包含取自40个人的样本,每人10幅,共400幅图像,图像尺寸92*112。

at_t_datasets        样本库的链接为http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html

        1.欧式距离降维及分类效果:

Euclidean_Distance_result        即使将h设为最大的400(样本数),其分类正确率也只能达到99%。

        2.曼哈顿距离降维及分类效果:

Manhattan_Distance_result        在h为128时,分类正确率可以达到100%,降维能力略好于欧式距离。

        3.马氏距离降维及分类效果:

Mahalanobis_Distance_result        在h仅仅为7的时候,其分类正确率就已经达到100%了。采用马氏距离的PCA方法可以将人脸数据的维度从10000左右降到7,降维效果非常优秀。在h超过9时,分类过程中计算的最大马氏距离超出了双精度浮点数double的上限。


参考文献

[1] Sharma A, Paliwal K K, Imoto S, et al. Principal component analysis using QR decomposition[J]. International Journal of Machine Learning and Cybernetics, 2013, 4(6): 679-683.

[2] Raj D. A Realtime Face Recognition system using PCA and various Distance Classifiers[J]. 2011.

[3] Turk M, Pentland A. Eigenfaces for recognition[J]. Journal of Cognitive Neuroscience, 1991, 3(1): 71-86.

共有31条评论

  1. 你好,我想请问一下,你这篇博客的思想是直接对整张图片进行匹配是吧。这样的话当人脸角度发生变化,或者拍照的环境变化时,结果应该会受到比较大的影响,是不是去分别比对眼睛,鼻子,嘴巴,结果会更加可信呢?

    1. 你好,文中的方法其实完全是统计意义上的分类,只有样本在外观上很接近时才能保证有很高的精度。你所说的方法属于基于特征的分类,方法也有很多,效果一般也更好。比较经典的从五官位置辅助分类的方法比如香港中文大学的DeepID(2014年):http://mmlab.ie.cuhk.edu.hk/pdf/YiSun_CVPR14.pdf,他们所采用的五官定位算法:http://www.cv-foundation.org/openaccess/content_cvpr_2013/papers/Sun_Deep_Convolutional_Network_2013_CVPR_paper.pdf。

      当然,人脸姿态变化比较大时,五官的变化也很大甚至发生遮挡。最近很流行的思路是用GAN(及其变种)来训练生成网络,顺便也能获得一个高精度的分类器,比如上礼拜中科院的TP-GAN:https://arxiv.org/pdf/1704.04086.pdf

    1. 做过,用的是opencv官方文档里常见的Haar分类器。只不过目前OpenCV Android应用开发的代码更新工作还没有进行到人脸识别的章节。

  2. 你好,我用的vs2010+opencv2.4.9,编译时出现

    fatal error C1083: 无法打开包括文件:“armadillo”: No such file or directory

    想问一下什么原因

  3. 我没有用Intel MKL,而是用LAPACK为Armadillo提供矩阵分解计算。编译通过了,运行时提示warning:svd():decomposition failed,error:Mat::cols(): indices out of  bounds or incorrectly used .测试时arma::svd_econ(U, d, V, R.t());返回0。请问这是什么原因。

      1. 我设置断点单步调试在78行时arma::svd_econ(U, d, V, R.t());svd_econ()返回值已经为0了,

        运行到89行:V_h = V.cols(0, component_num-1);程序出错

        后跳转到了debug.hpp中的这一段代码template<typename T1>
        arma_cold
        arma_noinline
        static
        void
        arma_stop_logic_error(const T1& x)
          {
          #if defined(ARMA_PRINT_ERRORS)
            {
            get_stream_err1() << "\nerror: " << x << std::endl;
            }
          #endif
          
          throw std::logic_error( std::string(x) );
          }

          1. 我设置7,不行

            错误信息:

            Unhandled exception at 0x000007FEFD82A49D in test1.exe: Microsoft C++ exception: std::logic_error at memory location 0x00000000002EDC18.

            这有可能是环境配置的问题吗

          2. 而且也正确读入了?是不是前面哪个步骤出了问题,用断点看一下每一行输出矩阵的尺寸是否和公式符合。实在不行,把你的代码发到我的邮箱johnhany@163.com,我来看看

          3. 我单步调试发现应该是arma::svd_econ(U, d, V, R.t());失败,导致U,V被重置,然后出错。但我不知道为什么会分解 失败,代码已经发到你的邮箱,请把我看一下

    1. 我在我电脑上装上intel的MKL,然后添加到vs2013中,就可以成功运行了。我想之前的错误应该是Armadillo自带的库分解能力问题造成的。

      非常感谢之前的热情帮助

  4. 您好,我正在进行android平台上的人脸识别开发,我用JNI的方法,需要把识别后的图像(Mat类型)存入jintArray型数组,再把jintArray传回java层显示,你有什么好办法啊

    1. 奇异值分解在一般的矩阵理论书籍中都会有介绍,可以看成是对非方阵进行的特征值分解。通俗地讲,矩阵乘以自身的转置,这个乘积矩阵的特征值开根号就是原矩阵的奇异值。

发表评论

电子邮件地址不会被公开。 必填项已用*标注