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

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

原作者： John Hany

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

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

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

---

PCA人脸识别

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

        1.假设有400幅尺寸为100*100的图像，构成10000*400的矩阵；

        2.计算均值，令；

        3.根据定义，计算协方差矩阵；

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

        5.矩阵可对数据降维：，Y是h行400列的矩阵，也就是将数据从10000维降为h维。

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

---

利用QR分解减小计算量

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

        1.已知，其中为d*d，H为d*n，d代表原始数据的维数，n代表样本数，d远大于n；

        2.对H作QR分解，，其中Q为d*t，R为t*n，；

        3.则，对作奇异值分解，其中U为n*t，V为t*t，；

        4.于是，其中；

        5.由于，所以QV可将对角化，QV为的特征向量矩阵，为的特征值矩阵；

        6.选取D前h个最大对角元所对应于V中的h个列，构成t*h的矩阵，则降维矩阵；

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

        进一步，进行人脸识别的过程如下：

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

        2.对X计算，将Y（也称特征脸）按类别计算均值，得到c个长度为h的列向量；

        3.对于未知类别的新样本x，计算，y的长度为h；

        4.计算距离，取距离最小的i作为x的类标号。

---

距离度量d

        1.欧式距离

        2.曼哈顿距离

        3.马氏距离

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

---

C++实现

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

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

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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219

/************************************************************************/ /*                   QR-PCA-FaceRec  by John Hany                       */ /* */ /* A face recognition algorithm using QR based PCA.               */ /* */ /* Released under MIT license. */ /* */ /* Contact me at johnhany@163.com */ /* */ /* Welcome to my blog http://johnhany.net/, if you can read Chinese:) */ /* */ /************************************************************************/   #include <opencv2/opencv.hpp> #include <armadillo> #include <iostream>   using namespace std ;   int component_num = 7 ;   string orl_path = "G:\\Datasets\\orl_faces" ;   enum distance_type { ECULIDEAN = 0 , MANHATTAN , MAHALANOBIS } ; //double distance_criterion[3] = { 10.0, 30.0, 3.0}; double distance_criterion [ 3 ] = { 1000.0 , 1000.0 , 1000.0 } ;   bool compDistance ( pair < int , double > a , pair < int , double > b ) ; double calcuDistance ( const arma :: vec vec1 , const arma :: vec vec2 , distance_type dis_type ) ; double calcuDistance ( const arma :: vec vec1 , const arma :: vec vec2 , const arma :: mat cov2 , distance_type dis_type ) ;   int main ( int argc , const char * argv [ ] ) { int class_num = 40 ; int sample_num = 10 ;   int img_cols = 92 ; int img_rows = 112 ; cv :: Size sample_size ( img_cols , img_rows ) ;   arma :: mat mat_sample ( img_rows* img_cols , sample_num* class_num ) ;   //Load samples in one matrix `mat_sample`.   for ( int class_idx = 0 ; class_idx < class_num ; class_idx ++ ) { for ( int sample_idx = 0 ; sample_idx < sample_num ; sample_idx ++ ) {   string filename = orl_path + "\\s" + to_string ( class_idx + 1 ) + "\\" + to_string ( sample_idx + 1 ) + ".pgm" ; cv :: Mat img_frame = cv :: imread ( filename , CV_LOAD_IMAGE_GRAYSCALE ) ; cv :: Mat img_sample ; cv :: resize ( img_frame , img_sample , sample_size ) ;   for ( int i = 0 ; i < img_rows ; i ++ ) { uchar* pframe = img_sample . ptr < uchar > ( i ) ; for ( int j = 0 ; j < img_cols ; j ++ ) { mat_sample ( i* img_cols + j , class_idx* sample_num + sample_idx ) = ( double ) pframe [ j ] / 255.0 ; } } } } // cout << mat_sample.n_rows << endl << mat_sample.n_cols << endl << mat_sample(img_rows*img_cols/2, 0) << endl;   //Calculate PCA transform matrix `mat_pca`.   arma :: mat H = mat_sample ; arma :: mat mean_x = arma :: mean ( mat_sample , 1 ) ;   for ( int j = 0 ; j < class_num * sample_num ; j ++ ) { H . col ( j ) -= mean_x . col ( 0 ) ; } H /= 1.0 / sqrt ( sample_num - 1 ) ;   arma :: mat Q , R ; arma :: qr_econ ( Q , R , H ) ;   arma :: mat U , V ; arma :: vec d ; arma :: svd_econ ( U , d , V , R . t ( ) ) ;   // cout << "d" << endl << d << endl;   // arma::rowvec vec_eigen = d.head(component_num).t(); // cout << "vec_eigen" << endl << vec_eigen << endl;   arma :: mat V_h ( V . n_rows , component_num ) ; if ( component_num == 1 ) { V_h = V . col ( 0 ) ; } else { V_h = V . cols ( 0 , component_num - 1 ) ; }   arma :: mat mat_pca = Q * V_h ;   //Calculate eigenfaces `mat_eigen_vec`.   arma :: mat mat_eigen = mat_pca . t ( ) * mat_sample ; // cout << "mat_eigen" << endl << mat_eigen << endl; arma :: mat mat_eigen_vec ( component_num , class_num , arma :: fill :: zeros ) ; vector < arma :: mat > mat_cov_list ;   for ( int class_idx = 0 ; class_idx < class_num ; class_idx ++ ) {   arma :: vec eigen_sum ( component_num , arma :: fill :: zeros ) ; for ( int sample_idx = 0 ; sample_idx < sample_num ; sample_idx ++ ) { eigen_sum += mat_eigen . col ( class_idx* sample_num + sample_idx ) ; } eigen_sum /= ( double ) sample_num ; mat_eigen_vec . col ( class_idx ) = eigen_sum ;   mat_cov_list . push_back ( arma :: cov ( ( mat_eigen . cols ( class_idx* sample_num , class_idx* sample_num + sample_num - 1 ) ) . t ( ) ) ) ;   // cout << mat_cov_list[class_idx] << endl;   }   // cout << "mat_eigen_vec" << endl << mat_eigen_vec << endl;   /* cout << "dis within class" << endl; for(int class_idx = 0; class_idx < class_num; class_idx++) { for(int sample_idx = 0; sample_idx < sample_num; sample_idx++) { double dis = calcuDistance(mat_eigen.col(class_idx*sample_num+sample_idx), mat_eigen_vec.col(class_idx), mat_cov_list[class_idx], distance_type::MAHALANOBIS); cout << dis << " "; } cout << endl; }   cout << "dis between classes" << endl; for(int class_idx = 0; class_idx < class_num; class_idx++) { for(int sample_idx = 0; sample_idx < class_num; sample_idx++) { double dis = calcuDistance(mat_eigen.col(sample_idx*sample_num), mat_eigen_vec.col(class_idx), mat_cov_list[class_idx], distance_type::MAHALANOBIS); cout << dis << " "; } cout << endl; } */   //Classify new sample.   int correct_count = 0 ;   double max_dis = 0.0 ;   for ( int class_idx = 0 ; class_idx < class_num ; class_idx ++ ) { for ( int sample_idx = 0 ; sample_idx < sample_num ; sample_idx ++ ) { arma :: mat mat_new_sample ( img_rows* img_cols , 1 ) ;   string filename = orl_path + "\\s" + to_string ( class_idx + 1 ) + "\\" + to_string ( sample_idx + 1 ) + ".pgm" ; cv :: Mat img_new_frame = cv :: imread ( filename , CV_LOAD_IMAGE_GRAYSCALE ) ; cv :: Mat img_new_sample ; cv :: resize ( img_new_frame , img_new_sample , sample_size ) ;   for ( int i = 0 ; i < img_rows ; i ++ ) { uchar* pframe = img_new_sample . ptr < uchar > ( i ) ; for ( int j = 0 ; j < img_cols ; j ++ ) { mat_new_sample ( i* img_cols + j , 0 ) = ( double ) pframe [ j ] / 255.0 ; } }   arma :: mat mat_new_eigen = mat_pca . t ( ) * mat_new_sample ;   vector < pair < int , double >> dis_list ; for ( int new_class_idx = 0 ; new_class_idx < class_num ; new_class_idx ++ ) { double dis = calcuDistance ( mat_new_eigen . col ( 0 ) , mat_eigen_vec . col ( new_class_idx ) , mat_cov_list [ new_class_idx ] , distance_type :: MAHALANOBIS ) ; dis_list . push_back ( make_pair ( new_class_idx , dis ) ) ; } sort ( dis_list . begin ( ) , dis_list . end ( ) , compDistance ) ;   if ( dis_list [ 0 ] . first == class_idx && dis_list [ 0 ] . second <= distance_criterion [ distance_type :: MAHALANOBIS ] ) { correct_count ++ ; }   if ( dis_list . back ( ) . second > max_dis ) { max_dis = dis_list . back ( ) . second ; } } }   cout << "Maximum distance: " << max_dis << endl ;   double correct_ratio = ( double ) correct_count / ( class_num * sample_num ) ; cout << "Correctness ratio: " << correct_ratio * 100.0 << "%" << endl ;   cin . get ( ) ;   return 0 ; }   bool compDistance ( pair < int , double > a , pair < int , double > b ) { return ( a . second < b . second ) ; }   double calcuDistance ( const arma :: vec vec1 , const arma :: vec vec2 , distance_type dis_type ) {   if ( dis_type == ECULIDEAN ) { return arma :: norm ( vec1 - vec2 , 2 ) ; } else if ( dis_type == MANHATTAN ) { return arma :: norm ( vec1 - vec2 , 1 ) ; } else if ( dis_type == MAHALANOBIS ) { arma :: mat tmp = ( vec1 - vec2 ) . t ( ) * ( vec1 - vec2 ) ; return sqrt ( tmp ( 0 , 0 ) ) ; }   return - 1.0 ; }   double calcuDistance ( const arma :: vec vec1 , const arma :: vec vec2 , const arma :: mat cov2 , distance_type dis_type ) {   if ( dis_type == ECULIDEAN ) { return arma :: norm ( vec1 - vec2 , 2 ) ; } else if ( dis_type == MANHATTAN ) { return arma :: norm ( vec1 - vec2 , 1 ) ; } else if ( dis_type == MAHALANOBIS ) { arma :: mat tmp = ( vec1 - vec2 ) . t ( ) * cov2 . i ( ) * ( vec1 - vec2 ) ; return sqrt ( tmp ( 0 , 0 ) ) ; }   return - 1.0 ; }

---

分类测试

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

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

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

        即使将h设为最大的400（样本数），其分类正确率也只能达到99%。

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

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

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

        在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.