上一個筆記主要是講了PCA的原理，並給出了二維影像降一維的示例程式碼。但還遺留了以下幾個問題：
在計算協方差和特徵向量的方法上，書上使用的是一種被作者稱為
compact trick的技巧，以及奇異值分解（SVD），這些都是什麼東西呢？如何把PCA運用在多張圖片上？
所以，我們需要進一步的瞭解，同時，為示例對多張圖片進行PCA，我選了一個跟書相似但更有趣的例子來做——人臉識別。
人臉識別與特徵臉
一個特徵臉（Eigenface，也叫標準臉）其實就是從一組人臉影像應用PCA獲得的主成分特徵向量之一，下面我們能驗證，每個這樣的特徵向量變換為二維影像後看起來有點像人臉，所以才被稱為特徵臉，計算特徵臉是進行人臉識別的首要步驟，其計算過程其實就是PCA。
特徵臉計算步驟
準備一組（假設10張）具有相同解析度（假設:100 × 100）的人臉影像，把每張影像打平（numpy.flatten）成一個向量，即所有畫素按行串聯起來, 每個影像被當作是10000維空間的一個點。再把所有打平的影像儲存在 10000 × 10 矩陣X中，矩陣的每一列就是一張圖片，每個維為一行，共10000維
對X的每一維（行）求均值，得到一個10000 × 1的向量，稱為
平均臉（因為把它變換為二維影像看起來像人臉），然後將X減去平均臉，即零均值化計算X的協方差矩陣
C=XX^(X^表示X的轉置)計算C的特徵值和特徵向量，這組特徵向量就是一組特徵臉。在實際應用中，我們只需要保留最主要的一部分特徵向量做為特徵臉即可。
有了上個筆記的基礎知識後，上面計算過程不難理解。但在實現程式碼之前，我們先來看看上面提到的計算X的協方差矩陣C=XX^引發的一個問題。
協方差計算技巧
對於上面舉例的矩陣X，它有10000行（維），它的協方差矩陣將達到 10000 × 10000，有10000個特徵向量，這個計算量是很大的，消耗記憶體大，我嘗試過不可行。
這個數學問題終究還是被數學解決了，解決的方法可以見維基百科特徵臉內容描述。簡言之，就是把本來計算XX^的協方差矩陣（設為C）變換為計算X^X的協方差矩陣（C'），後者的結果是 10 × 10（10為樣本數量），很快就可以算出來。當然，通過這個變換分別算出來的特徵向量不是等價的，也需要變換一下：設E為從C算出來的特徵向量矩陣，E'為從C'算出來的特徵向量矩陣，則E = XE'，最後再把E歸一化。
這個技巧就是書上PCA示例使用的，被稱為compact trick的方法。
但要看明白書上的示例程式碼，還要搞清一點：
對原始影像資料集矩陣的組織方式，我們用行表示維度，列表示樣本，而書上和《Guide to face recognition with Python》（見底部參考連結）使用的是行表示樣本，列表示維度。就是因為這兩種組織方式的不同，導致了PCA演算法的程式碼看起來有些不同。這一點很容易讓人困惑，所以寫到這裡，我應該特別的強調一下。
我之所以在上個筆記，包括上面對特徵臉的計算步驟描述，都認定以行表示維度，列表示樣本的方式，是為了與數學原理的詳解保持一致（注：下面的程式碼示例還是使用這種方式），當我們明白了整個原理之後，我們便知道使用這兩種矩陣表達方式都可以，兩者實現的PCA程式碼差別也很小，下面會講到。
準備人臉影像樣本
網上有不少用於研究的人臉資料庫可以下載，我在參考連結給出了常被使用的一個。下載解壓後在目錄orl_faces下包含命名為s1,..,s40共40個資料夾，每個資料夾對應一個人，其中儲存10張臉照，所有臉照都是92 × 112的灰度圖，我把部分照片貼出來：
接下來，我們按照特徵臉計算步驟中的第1點所述，把這400張影像組成矩陣（影像組織不分先後），程式碼：
def getimpaths(datapath):
paths = []
for dir in os.listdir(datapath):
try:
for filename in os.listdir(os.path.join(datapath, dir)):
paths.append(os.path.join(datapath, dir, filename))
except:
pass
return paths
impaths = getimpaths('./orl_faces')
m,n = np.array(Image.open(impaths[0])).shape[0:2] #圖片的解析度，下面會用到
X = np.mat([ np.array(Image.open(impath)).flatten() for impath in impaths ]).T
print 'X.shape=',X.shape #X.shape= (10304, 400)我們把每個影像都打平成行向量，把所有影像從上到下逐行排列，最後轉置一下，便得到一個10304 × 400 的矩陣，其中10304 = 92 × 112
實現PCA函式求解特徵臉
PCA函式
我儘量使用與書上相同的變數命名，方便對比：
def pca(X):
dim, num_data = X.shape #dim: 維數， num_data: 樣本數
mean_X = X.mean(axis=1) #求出平均臉，axis=1表示計算每行（維）均值，結果為列向量
X = X - mean_X #零均值化
M = np.dot(X.T, X) #使用compact trick技巧，計算協方差
e,EV = np.linalg.eigh(M) #求出特徵值和特徵向量
print 'e=',e.shape,e
print 'EV=',EV.shape,EV
tmp = np.dot(X, EV).T #因上面使用了compact trick，所以需要變換
print 'tmp=',tmp.shape,tmp
V = tmp[::-1] #將tmp倒序，特徵值大的對應的特徵向量排前面，方便我們挑選前N個作為主成分
print 'V=',V.shape,V
for i in range(EV.shape[1]):
V[:,i] /= np.linalg.norm ( EV[:,i]) #因上面使用了compact trick，所以需要將V歸一化
return V,EV,mean_X執行PCA並畫圖
對上面得到的X矩陣呼叫pca函式，並畫出平均臉和部分特徵臉：
V,EV,immean = pca(X)
#畫圖
plt.gray()
plt.subplot(2,4,1) #2行4列表格，第一格顯示`平均臉`
plt.imshow(immean.reshape(m,n))
#以下選前面7個特徵臉，按順序分別顯示到其餘7個格子
for i in range(7):
plt.subplot(2,4,i+2)
plt.imshow(V[i].reshape(m,n))
plt.show()顯示效果圖如下：
希望不會被這些特徵臉嚇到：）
這些所謂的臉事實上是特徵向量，只不過維數與原始影像一致，因此可以被變換成影像顯示出來，不同的特徵臉代表了與均值影像差別的不同方向。
當然，我們求特徵臉，並不是為了顯示他們，而是保留部分特徵臉來獲得大多數臉的近似組合。因此，人臉便可通過一系列向量而不是原始數字影像進行儲存，節省了很多儲存空間，也便於後續的識別計算。
與書上pca的實現對比
上面我給出的pca函式程式碼，是按照我們一路學習PCA的思路寫出來的，雖然跟書上pca實現很接近，但還是有幾個點值得分析：
如上提到，我們對X矩陣的組織是以行為維、列為樣本的方式，即一個列對應一張打平的影像，而書上的例子是以行為樣本、列為維的方式，每一行對應一張打平的影像，而且參考連結裡的例子也都是以後者進行組織的，但沒關係，我們只需要對上面的程式碼作一點修改即可：
def pca_book(X):
num_data, dim = X.shape #注意：這裡行為樣本數，列為維
mean_X = X.mean(axis=0) #注意：axis=0表示計算每列（維）均值，結果為行向量
X = X - mean_X
#M = np.dot(X.T, X) #把X轉置後代入，得到
M = np.dot(X, X.T) #跟書上一樣
e,EV = np.linalg.eigh(M) #求出特徵值和特徵向量
#tmp = np.dot(X, EV).T #把X轉置後代入，得到
tmp = np.dot(X.T, EV).T #跟書上一樣
V = tmp[::-1]
#以下是對V歸一化處理，先省略，下面講所以我們看到，其實演算法的本質是一樣的，只不過要注意的地方是維數和樣本數反過來了，另外，對X的運算換成X的轉置即可。類推的，如果X使用我們的上面的組織方式，呼叫pca_book函式的程式碼為V,EV,immean = pca_book(X.T)
歸一化演算法不同。因為使用書上的方法，在對特徵值求平方根(np.sqrt(e))的時候會產生一個錯誤（負數不能開平方根），所以我這裡使用的歸一化方法是從《Guide to face recognition with Python》抄來的。書上的pca演算法多了一個判斷分支，當
dim <= num_data即維數少於樣本數的時候直接使用SVD（奇異值分解）演算法，顯然在一般的人臉識別的例子裡，不會被用到，因為單張92 × 112影像打平後維數為10304，而樣本數為400，遠遠低於維數。
歸一化
原先以為歸一化是一種比較簡單的運算，一瞭解才發現原來是一種不簡單的思想，在機器學習中常被使用，看回上面的程式碼：
for i in range(EV.shape[1]):
V[:,i] /= np.linalg.norm ( EV[:,i])首先讀者得自行了解範數(norm)的概念, 範數（norm）還分L0、L1、L2等好幾種，而函式np.linalg.norm就是用來求矩陣或向量的各種範數，預設就是求L2範數，具體可查閱《linalg.norm API說明》
上面程式碼的作用就是對V每一列歸一化到單位L2範數。
而書上使用的歸一化方法是：
S = sqrt(e)[::-1] #計算e的平方根再對結果倒序排列
for i in range(V.shape[1]):
V[:,i] /= S我在網上找到了關於compat trick後如何對求得的向量歸一化的數學推薦，截圖如下：
這跟左奇異值有關，屬於SVD中的內容，有興趣的話自行研究。
當我使用這種方法實現時，程式執行出現錯誤：FloatingPointError: invalid value encountered in sqrt，發現是對負數開平方根產生了錯誤，也就是說對協方差矩陣求得的特徵值中包含了負數。那麼，如果要使用這種歸一化方法的話，是否只要排除掉負的特徵值及其對應的特徵向量就可以了？
SVD（奇異值分解）
我們的程式碼示例的PCA方法使用的是特徵分解，線性代數中，特徵分解（Eigendecomposition），又稱譜分解（Spectral decomposition）,是將矩陣分解為由其特徵值和特徵向量表示的矩陣之積的方法，但需要注意只有對可對角化矩陣才可以施以特徵分解。
而SVD（singular value decomposition）能夠用於任意m乘n矩陣的分解，故SVD適用範圍更廣。
但是，如果矩陣維數很大，如我們之前所舉例10000維的時候，計算SVD也是很慢的，所以我們看到書上的例子增加了一個分支判斷，當維度 < 樣本數的時候，才使用SVD，否則使用compat trick方法的PCA。
回顧一下上篇筆記舉的PCA應用例子：把一張二維影像變換成一維，維數為2，對於這個例子，直接使用SVD是比較合適的，這樣PCA函式將變得很簡潔。
利用特徵臉進行人臉識別
這裡不會詳細地討論如何實現一個好的人臉識別演算法，而是為了示例PCA的運用，所以只是簡單的介紹一下。
上面我們已經得出了400張人臉的特徵臉（特徵向量），首先第一個問題，我們得選擇多少個特徵向量作為主成分？
特徵值貢獻率
假設我們選擇k個特徵向量，其對應的特徵值之和與所有特徵值之和的比值，就是這k個特徵值的貢獻率。所以主成分的選擇問題就轉化為選擇k個特徵向量，使得特徵值的貢獻率大於等於某個值（如90%）即可。我們把選定的k個特徵向量組成的矩陣設為W。
識別步驟
第一步：將每個人臉樣本影像減去平均影像後，投影到主成分上
W = EV[:,k] #假設k已經根據特徵值貢獻率算出來了
projections = [] #存放每個人臉樣本的投影
for xi in X.T: #X為我們之前組織的所有人臉樣本的 10000 × 400矩陣
projections.append(np.dot(W.T, xi - mean_X)) #mean_X為之前我們求得的平均影像第二步：設要識別的影像為D，將D也投影到主成分上得到Q，然後計算Q與各個樣本人臉投影的歐幾里得距離，得出最小的歐幾里得距離
def euclidean_distance(p, q): #求歐幾里得距離
p = np.array(p).flatten()
q = np.array(q).flatten()
return np.sqrt(np.sum(np.power((p-q) ,2)))
minDist = np.finfo('float').max
Q = np.dot(W.T, D - mean_X)
for i in xrange (len(projections)):
dist = euclidean_distance( projections[i], Q)
if dist < minDist :
minDist = dist如果要識別的影像是樣本影像之一，那麼求得的最小的歐幾里得距離對應的樣本影像與要識別的影像是同一個人。若要識別的影像非樣本之一或根本不是人臉，我們就需要有一個閥值與求得的最小的歐幾里得距離作比較，若在閥值之外，則可以判斷要識別的影像不在樣本中。關於如何設定閥值，本文不再討論。
小結
人臉識別的方法有很多種，基於特徵臉的識別只是其中一種。但要實現一個可用的人臉識別，還有很多問題要考慮。另外PCA本身對某些特定情況的原始資料集也存在一些缺點。
至此，關於PCA，將不再深入探討。
在PCA的學習過程中，深感一種技術應用的背後，必有驚豔的數學原理支撐，體會了一把數學之美:)，但因本人數學水平有限（後悔大學時沒好好學），對以上理解必存在錯漏和不詳之處，所以也是很歡迎你的批評指正。
參考文件
維基百科：特徵臉
compute pca with this useful trick
Guide to face recognition with Python
FACE RECOGNITION HOMEPAGE
人臉影像資料庫下載