PCA——主成分分析

　　PCA（Principal Components Analysis）主成分分析是一個簡單的機器學習算法，利用正交變換把由線性相關變量表示的觀測數據轉換為由少量線性無關比變量表示的數據，實現降維的同時盡量減少精度的損失，線性無關的變量稱為主成分。大致流程如下：

　　首先對給定數據集（數據是向量）進行規范化，使得數據集的平均值為0，方差為1（規范化是為了使數據散布在原點附近，而不是遠離原點的某塊區域，便於后面的計算）。之后對每個數據進行正交變換，把數據投影到幾個少量的相互正交的方向（這些方向構成了數據空間的一個子空間）上。數據在每個方向上都有對應的坐標，而用這些方向和對應的坐標（坐標×方向的累加）就能近似表示原來高維的數據。因此這些方向的求解就是PCA的關鍵，這些方向有如下性質：

　　如果再由這些坐標通過對應的方向映射回原來的數據，精度損失是同等方向數量的方向集合（或者叫同維度的子空間吧）中最小的，而數據集在各個方向上的坐標的方差之和是同等方向數量的方向集合中最大的，也正對應着方差越大表示這個方向上保存數據的信息量越大（方差越小所有數據越集中在一個點上信息量當然越小）。數據集在這些方向的上的坐標的方差，從大到小排序，就把這每一個方向稱為第一主成分、第二主成分……

主成分推導

　　接下來推導什么樣的方向是主成分，即什么樣的方向集合能保存原數據集更多的信息。下面把數據映射為各個方向上的坐標稱為編碼（即降維），再反映射回來稱為解碼。

　　定義矩陣$X\in R^{m\times n}$為數據集，由$m$個$n$維行向量數據$\{x_1,x_2,...,x_m\}$組成。首先求得數據均值與方差（以下是按元素進行的運算）：

$\displaystyle\overline{x}=\sum\limits_{i=1}^{m}\frac{x_i}{m}$

$\displaystyle s=\frac{1}{m-1} \sum\limits_{i=1}^{m}\left(x_i-\overline{x}\right)^2 $

　　然后規范化樣本：

$\displaystyle x_i^\ast = \frac{x_i-\overline{x}}{\sqrt{s}}, i=1,2,\ldots,m$

　　獲得規范化后的數據集$X^*$。隨后PCA將把$X^*$編碼成矩陣$Y^*\in R^{m\times k},k<n$，通過右乘列向量單位正交的方向矩陣$D\in R^{n\times k}$得到編碼后的矩陣$Y^*$：

$Y^* = X^*\cdot D$

　　而再解碼回來是再乘上$D^T$：

$\hat{X^*} = Y^*\cdot D^{T}$

　　因為矩陣$D$並不可逆，所以編碼是不可逆的，也就是說解碼並不能完全恢復回$X^*$。為了最小化精度損失，容易想到，優化$D$來最小化$\hat{X^*}$與$X^*$之差的Frobenius范數：

$ \begin{aligned} &D = \displaystyle\text{arg}\min\limits_D{||\hat{X^*} - X^*||^2_F}\\ &D = \text{arg}\min\limits_D{|| X^*DD^T - X^*||^2_F}\\ &D = \text{arg}\min\limits_D{\text{Tr}[(X^*DD^T - X^*)^T(X^*DD^T - X^*)]}\\ &D = \text{arg}\min\limits_D{\text{Tr}[DD^TX^{*T}X^*DD^T-DD^TX^{*T}X^*-X^{*T}X^*DD^T+X^{*T}X^*]}\\ \end{aligned} $

　　去掉與$D$不相關的$X^{*T}X^*$項，再由跡的性質，矩陣乘法循環右移跡的值不變：

 $ \begin{aligned} &D = \text{arg}\min\limits_D{\text{Tr}[X^*DD^TDD^TX^{*T}-X^*DD^TX^{*T}-D^TX^{*T}X^*D]}\\ &D = \text{arg}\min\limits_D{\text{Tr}[X^*DD^TX^{*T}-X^*DD^TX^{*T}-D^TX^{*T}X^*D]}\\ &D = \text{arg}\max\limits_D{\text{Tr}[D^TX^{*T}X^*D]}\\ \end{aligned} $

　　容易發現當$D$的列向量取$X^{*T}X^*$的前$k$大特征值對應的特征向量時，有這個跡的值最大，等於前$k$大特征值之和，且此時壓縮的精度損失最少。

性質

　　經過上面的推導，我們可以發現，$X^{*T}X^*$實際上就是$X$行向量每個元素之間的協方差矩陣$\Sigma_x$（沒有除以$m-1$）。定義$\lambda_1,\lambda_2,...,\lambda_n$為$\Sigma_x$從大到小排列的特征值。我們可以推出下面兩個性質。

　　性質一、編碼得到的矩陣$Y^*\in R^{m\times k}$的行向量各個元素對應的協方差矩陣$\displaystyle\Sigma_y = \text{diag}(\lambda_1,\lambda_2,...,\lambda_k)$ 。

　　性質二、如果不進行壓縮，即編碼矩陣形如$D\in R^{n\times n}$，從而$Y^*\in R^{m\times n}$，那么$Y^*$的行向量各個元素方差之和等於$X^*$（$X$也一樣）的行向量各個元素方差之和，也就是：

$\displaystyle\sum\limits_{i=1}^n \Sigma_{y_{ii}}= \sum\limits_{i=1}^n \Sigma_{x_{ii}} = \sum\limits_{i=1}^n (X^{*T}X^*)_{ii}$

　　首先證明性質一。

　　由之前的定義，$X^*$是$X$經過行標准化后的矩陣，所以對$X^*$的任意列進行求和有

$\displaystyle\sum\limits_{i=1}^mX^*_{ij}=0,\;\;j=1,2,...,n$

　　而對於$Y^*=X^*D$，對$Y^*$的任意列進行求和有

$ \begin{aligned} \sum\limits_{i=1}^m Y^*_{ij} &= \sum\limits_{i=1}^mX^*_{i:}D_{:j} = \sum\limits_{i=1}^m\sum\limits_{k=1}^nX^*_{ik}D_{kj}\\ &= \sum\limits_{k=1}^nD_{kj}\sum\limits_{i=1}^mX^*_{ik}= 0,\;\;j=1,2,...,k \end{aligned} $

　　每列求和都為0，說明，$Y^*$行向量每個元素的均值都為0，無需再計算行向量的均值進行標准化。因此，可以直接計算$Y^*$的協方差矩陣：

$ \begin{aligned} \Sigma_y = Y^{*T}Y^* = D^TX^{*T}X^*D=\text{diag}(\lambda_1,\lambda_2,...,\lambda_k) \end{aligned} $

　　有性質一以后，性質二就容易了：

　　如果不進行壓縮，由性質一可得，$Y^*$行向量各元素方差之和就是$X^{*T}X^*$的特征值之和，矩陣特征值之和又等於矩陣的跡，而$X^{*T}X^*$的跡就是$X^*$行向量各個元素的方差之和。所以得到性質二。 

　　另外，映射矩陣$D$的列向量在人臉識別中就是所謂的特征臉（PCA+SVM實現人臉識別）。

手動代碼實現

　　實驗代碼：

import matplotlib.pyplot as plt
import pylab                #原本在jupyter里才能顯示的圖片，可以用窗口顯示
import numpy  as np 

def PCA_img(img,n):
    t = np.mean(img,0)                                          #以行向量為壓縮單元，求行平均
    s = np.std(img,0)                                           #求行標准差
    img=(img-t)/s                                               #每行減去平均值進行規范化
    v,vectors = np.linalg.eig(np.dot(img.T,img)/(len(img)-1))   #生成規范化矩陣后，轉置乘自身再除以（行數-1），即為行向量各個元素的協方差矩陣，求協方差矩陣的特征值與特征向量
    indv = v.argsort()                                          #將特征值排序，從小到大
    vectors = vectors[:,indv]                                   #將特征向量按特征值大小排序
    vectors = vectors[:,-n:]                                    #選擇前n大特征值對應的特征向量，作為映射矩陣。在人臉識別中就是特征臉
    img = np.dot(np.dot(img,vectors),vectors.T)                 #規范化矩陣先乘映射矩陣再乘映射矩陣的轉置，也就是先壓縮再解壓縮
    img = np.multiply(img,s)+t                                  #把減掉的均值加回去變為原矩陣
    return img
        
img=plt.imread("aaa.jpg") #讀取圖片  

img.flags["WRITEABLE"] = True
img = img.astype(float) 
n = 50
img[:,:,0] = PCA_img(img[:,:,0],n)
img[:,:,1] = PCA_img(img[:,:,1],n)
img[:,:,2] = PCA_img(img[:,:,2],n)
fig = plt.figure()           
ax = fig.add_subplot(111)    
ax.imshow(img/255)  
print(img)           
pylab.show()              

　　原圖是1200×1920的圖：

　　壓縮為1200×50后再解壓回去：

sklearn包PCA

　　sklearn中計算PCA比直接使用numpy矩陣計算快得多，這是因為它里面實現了比較高效的隨機算法，雖然每次算出來的值不一樣，但是相差不大。代碼如下：

import sklearn.decomposition as dc 

pca = dc.PCA(n_components = 10)       #規定PCA要取前k大特征值對應的特征向量的個數，即要映射到幾維，PCA默認壓縮行，比如100*100壓縮到100*10
pca.fit(A)                            #將要用於壓縮的矩陣傳入，用來“學習”壓縮矩陣
D = pca.components_                   #獲取“學習”到的編碼矩陣（壓縮矩陣）D。原矩陣A經過AD^T會得到壓縮后的矩陣（這里的D是上面推算的D^T）
C = pca.transform(B)                  #傳入新的矩陣B，用“學習”到的編碼矩陣進行壓縮，即BD^T。在人臉特征提取中，同一個人的不同人臉照片，用同一個壓縮矩陣壓縮后的坐標向量通常相似，因此可以用於人臉識別的降維