摘要：上一個筆記主要是講了的原理，并給出了二維圖像降一維的示例代碼。當我使用這種方法實現時，程序運行出現錯誤，發現是對負數開平方根產生了錯誤，也就是說對協方差矩陣求得的特征值中包含了負數。而能夠用于任意乘矩陣的分解，故適用范圍更廣。

上一個筆記主要是講了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函數
我盡量使用與書上相同的變量命名，方便對比：

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
人臉圖像數據庫下載