從QR分解到PCA，再到人臉識別

時間  2021-07-13

  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。

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 [email protected] */
/* */
/* 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。

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.


轉載http://johnhany.net/2016/05/from-qr-decomposition-to-pca-to-face-recognition/

联系我们 沪ICP备13005482号-10