摘要：從而將高維圖像識別問題轉化為低維特征向量的識別問題，大大降低了計算復雜度，同時也減少了冗余信息所造成的識別誤差。向量上的大量藍色圓點白色邊緣表示二維數據在其上的投影。

原書對于PCA的講解只有一小節，一筆帶過的感覺，但我發現PCA是一個很重要的基礎知識點，在機器機視覺、人臉識別以及一些高級圖像處理技術時都被經常用到，所以本人自行對PCA進行了更深入的學習。

PCA（Principal Component Analysis，主成分分析或主元分析）是一種算法，PCA的結果是用盡可能少的特征數據來表達最多的原始圖像的本質結構特征。即使它會丟失一部分原始圖像的特征表達，但它仍然是很有用的處理技巧，也很常用，特別在計算機視覺和人臉識別方面。

假設我們有一個二維數據集，它在平面上的分布如下圖：

如果我們想要用一個一維向量來表達此數據集，就會丟失一部分此數據集的信息，但我們的目標是讓求得的這個一維向量可以盡可能多地保留這個數據集的特征信息，那么這個求解過程就是PCA。

通過PCA我們可以找到若干個1維向量，如圖：

直觀上看出，向量u1是數據集變化的主方向，而u2是次方向，u1比u2保留了更多的數據集的結構特征，所以我們選擇u1作為主成分，并把原數據集投影到u1上就可以得出對原數據集的一維近似重構：

以上只是一種直觀的示例，把二維數據降為用1維來表示，當然，PCA通常是應用在高維數據集上。

假設我們有10張100 × 100像素的灰度人臉圖，我們目標是要計算這10張圖的主成分來作為人臉特征，這樣就可以基于這個‘特征臉’進行人臉匹配和識別。但即使一個100 × 100像素的灰度圖像維度就達到10,000維，10張圖像的線性表示可以達到100,000維，如此高維的數據帶來幾個問題：

對高維數據集進行分析處理的計算量是巨大的，消耗資源太大，時間太長

高維數據包含了大量冗余和噪聲數據，會降低圖像識別率

所以通常對于高維數據集，首先需要對其進行降維運算，以低維向量表達原數據集最多最主要的結構特征。從而將高維圖像識別問題轉化為低維特征向量的識別問題，大大降低了計算復雜度，同時也減少了冗余信息所造成的識別誤差。PCA其實就是最常用的一種降維算法。

PAC也可用于高維數據壓縮、高維數據可視化（轉二維或三維后就可以畫圖顯示）等方面，也是其它很多圖像處理算法的預處理步驟。

關于PCA，網上一搜還是不少的，但我仔細看了幾篇文章之后，發現這些文章講的跟書上講的有些地方不一致，甚至連計算時的公式都不一樣，這讓我產生了很多困惑。所以我想最重要的還是要理解PCA的數學原理，數學原理才是根，掌握了數學原理，再來寫代碼。恰好找到一篇文章專門講PCA數學原理，作者的數學功底和邏輯表達能力非常棒，讓我很容易看明白。另外，我也找到一篇老外寫的文章（見底部參考文章），這兩篇文章對PCA的計算描述是一致的，所以我決定在這兩篇文章的基礎上，結合書上的示例代碼進行學習和驗證。

本文不是要要把PCA的數學原理及推導寫出來，而是通過理解PCA的數學原理，總結PCA的計算步驟，因為計算步驟是代碼實現的必備知識。

PCA的計算過程涉及到幾個很重要的數學知識點：

零均值化

矩陣的轉置及乘法

協方差與協方差矩陣

特征值及特征向量

現在來看PCA的計算步驟：
1）將原始數據按列組成d行n列矩陣X
重要說明：d對應的就是數據的字段（或叫變量、特征、維，下稱’維‘），而n表示n條記錄（或叫樣本、觀察值，下稱’樣本‘），即每1列對應1個樣本，之所以行和列這樣安排，是為了與數學公式保持一致，很多文章對這一點都沒有明確的說明，導致計算的方式各有不同，讓人產生不必要的困惑
2）將X的每個維（行）進行零均值化，即將行的每個元素減去這一行的均值
3）求出X的協方差矩陣C，即 X 乘 X的轉置
4）求出C所有的特征值及對應的特征向量
5）將特征向量按對應特征值大小從上到下按行排列成矩陣，取前k行組成矩陣E
6）Y=EX即為降維到k維后的數據

下面用一個例子來驗證一下這個計算過程。

書中的例子是用PCA計算特征臉（人臉識別中的一步），它應用在多張圖片上面。為直觀起見，我用另外一個例子——對單張黑白圖像進行PCA，相當于把二維圖像降為一維。

把圖像轉成二維矩陣
這張圖像是原書封面：

下面代碼是將圖中‘怪魚’部分截取出來，并轉成黑白圖像顯示：

from PIL import Image
pim = Image.open("cover.png").crop((110,360,460,675)).convert("1")
pim.show()

效果如圖：

之所以截取這部分的圖片，是因為我們大概能猜到這幅圖像降到一維后，其一維表示的向量應該跟怪魚的方向大概一致。

使用黑白圖像是因為黑點才是我們關心的數據，因為是這些黑點描繪了圖像，每個黑點有唯一確定的行和列位置，對應平面上的（x,y）坐標，于是我們就可以得到此圖像的 2乘n 矩陣表示：第一行表示x維，第二行表示y維，每一列表示一個點。參考代碼：

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

im = np.array(Image.open("cover.png").crop((110,360,460,675)).resize((256,230)).convert("L"))
n,m = im.shape[0:2]
points = []
for i in range(n):
    for j in range(m):
        if im[i,j] < 128.0:  #把小于128的灰度值當作黑點取出來
            points.append([float(j), float(n) - float(i)]) #坐標轉換一下
    
im_X = np.mat(points).T; #轉置之后，行表示維度（x和y），每列表示一個點（樣本）
print "im_X=",im_X,"shape=",im_X.shape

現在，我們按上面說明的計算步驟來實現PCA：

def pca(X, k=1): #降為k維
  d,n = X.shape
  mean_X = np.mean(X, axis=1) #axis為0表示計算每列的均值，為1表示計算每行均值
  print "mean_X=",mean_X
  X = X - mean_X
  #計算不同維度間的協方差，而不是樣本間的協方差，方法1：
  #C = np.cov(X, rowvar=1) #計算協方差，rowvar為0則X的行表示樣本，列表示特征/維度
  #方法2：
  C = np.dot(X, X.T)
  e,EV = np.linalg.eig(np.mat(C)) #求協方差的特征值和特征向量
  print "C=",C
  print "e=",e
  print "EV=",EV
  e_idx = np.argsort(-e)[:k] #獲取前k個最大的特征值對應的下標（注：這里使用對負e排序的技巧，反而讓原本最大的排在前面）
  EV_main = EV[:,e_idx]   #獲取特征值（下標）對應的特征向量，作為主成分
  print "e_idx=",e_idx,"EV_main=",EV_main      
  low_X = np.dot(EV_main.T, X)    #這就是我們要的原始數據集在主成分上的投影結果
  return low_X, EV_main, mean_X

OK,現在我們調用此PCA函數，并把原圖像和投影到一維向量后的結果也描繪出來：

low_X, EV_main, mean_X = pca(im_X)
print "low_X=",low_X
print "EV_main=",EV_main
recon_X = np.dot(EV_main, low_X) + mean_X  #把投影結果重構為二維表示，以便可以畫出來直觀的看到
print "recon_X.shape=",recon_X.shape

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(im_X[0].A[0], im_X[1].A[0],s=1,alpha=0.5)
ax.scatter(recon_X[0].A[0], recon_X[1].A[0],marker="o",s=100,c="blue",edgecolors="white")
plt.show()

畫散點圖函數pyplot.scatter說明：

matplotlib.pyplot.scatter(x, y, ...)
x: 數組，樣本點在X軸上的坐標集
y: 數組，樣本點在Y軸上的坐標集
s: 表示畫出來的點的縮放大小
c: 表示畫出來的點（小圓圈）的內部顏色
edgecolors: 表示小圓圈的邊緣顏色

運行以上代碼打印：

im_X= [[  23.   24.   25. ...,  215.  216.  217.]
 [ 230.  230.  230. ...,    5.    5.    5.]] shape= (2, 19124)
mean_X= [[ 133.8574566 ]
 [ 123.75941226]]
C= [[ 2951.65745054 -1202.25277635]
 [-1202.25277635  3142.71830026]]
e= [ 1841.14567037  4253.23008043]
EV= [[-0.73457806  0.67852419]
 [-0.67852419 -0.73457806]]
e_idx= [1] EV_main= [[ 0.67852419]
 [-0.73457806]]
low_X= [[-153.26147057 -152.58294638 -151.90442219 ...,  142.29523704
   142.97376123  143.65228541]]
EV_main= [[ 0.67852419]
 [-0.73457806]]
recon_X.shape= (2, 19124)

并顯示如下圖像：

從圖中看出，向量的方向跟位置與我們目測的比較一致。向量上的大量藍色圓點（白色邊緣）表示二維數據在其上的投影。

以上實現的PCA算法，跟我參考的兩篇文章所講的原理一致。但跟書中的PCA計算方法有一定的不同，但也因為使用的例子不一樣，對原始數據集的定義不一樣導致的，由于避免文章過長，這些放在后面再講。
至此，通過對PCA的計算過程的學習，了解了一些線性代數知識、numpy和pyplot模塊的一些接口的用法，接下來我打算做點更有興趣的事情——就是使用PCA來實現人臉識別。

參考鏈接：
PCA的數學原理
A Tutorial on Principal Component Analysis
NumPy函數索引