原理介紹

PCA

這篇教程講的非常詳細。

以圖片為例，將每張圖片（假設寬高為 m n）轉化為m*n的一維向量，然後減去每張圖片畫素的均值（使得每張圖片均值為零），然後資料集中的所有圖片一起合成一個大的矩陣X。

計算圖片矩陣的方差-協方差矩陣，我們的優化目標就是將這個矩陣對角化，使得欄位（代表一張圖片）兩兩之間的協方差為零，也就是相關性最小，欄位與自身的方差儘可能大。如下：
D=1mYYT=1m(PX)(PX)T=1mPXXTPT=P(1mXXT)PT=PCPT

Y=PX，P為一組基，Y為在P上對映後的資料，D為對映後資料的協方差矩陣，也就是我們的優化目標矩陣，X為資料矩陣。C=

數學上矩陣對角化只需計算出C的特徵值和特徵向量，並取前k大的特徵值對應的特徵向量，也就是投影之後方差最大的那些基，即可得到P矩陣，k就是降維後的維數。

2D PCA

相較於傳統PCA，2D PCA不用將圖片轉換為1D向量，極大的減少了資料的維度。基本思路和PCA相同，也是將資料投影到某組基上，然後使得投影之後資料的協方差矩陣成為對角陣。投影后資料的協方差矩陣為：

Sx=E[(Y−EY)(Y−EY)T]=E[AX−E(AX)][AX−E(AX)]T=E[(A−EA)X][(A−EA)

Y為投影之後的資料，Y=AX，A是原資料，X是一組正交基，Sx也就是協方差矩陣。上述公式變形之後之後可得（利用tr(AB)=tr(BA)，只要保證跡不變即可，這裡還有點沒搞懂為什麼）：

Sx=XT[E(A−EA)T(A−EA)]X

然後再定義矩陣Gt為：

Gt=E(A−EA)T(A−EA)
也就是原資料的協方差矩陣。

類似於傳統PCA，我們優化目標就是將使矩陣S_x成為對角陣，並找出對應的特徵向量和特徵值，也就是將G_t對角化，取前K大的特徵向量。

2D-2DPCA

雙向2D PCA的改進之處在於，普通2D PCA中的變換隻提取了資料矩陣行內的特徵，對行進行了變換，而雙向PCA則加上了對列進行的變換。

Sx=E[(B−EB)(B−EB)T]=E[ZTA−E(ZTX)][ZTA−E(ZTA)]T=ZTE[(A−EA)(A−EA)T]Z

列變換對應的Gt矩陣為：

Gt=E[(A−EA)(A−EA)T]

對列進行2D PCA也就是對上述矩陣提取特徵向量，得到一組基Z，然後對原始資料進行投影。

為了將列和行的變換合併，最後進行如下變換：

C=ZTAX

X為行方向上2D PCA得到的變換矩陣，Z為列方向上2D PCA得到的變換矩陣，A為原資料。要還原圖片資料也很簡單，已知C矩陣：

A=ZCXT

因為特徵向量相互正交，Z，X均為正交陣，正交陣的逆矩陣就是其轉置。

python 實現

用python實現了一下演算法：

# a implementation of 2D^2 PCA algorithm

import numpy as np
from PIL import Image

def PCA2D_2D(samples, row_top, col_top):
    '''samples are 2d matrices'''
    size = samples[0].shape
    # m*n matrix
    mean = np.zeros(size)

    for s in samples:
        mean = mean + s

    # get the mean of all samples
    mean /= float(len(samples))

    # n*n matrix
    cov_row = np.zeros((size[1],size[1]))
    for s in samples:
        diff = s - mean;
        cov_row = cov_row + np.dot(diff.T, diff)
    cov_row /= float(len(samples))
    row_eval, row_evec = np.linalg.eig(cov_row)
    # select the top t evals
    sorted_index = np.argsort(row_eval)
    # using slice operation to reverse
    X = row_evec[:,sorted_index[:-row_top-1 : -1]]

    # m*m matrix
    cov_col = np.zeros((size[0], size[0]))
    for s in samples:
        diff = s - mean;
        cov_col += np.dot(diff,diff.T)
    cov_col /= float(len(samples))
    col_eval, col_evec = np.linalg.eig(cov_col)
    sorted_index = np.argsort(col_eval)
    Z = col_evec[:,sorted_index[:-col_top-1 : -1]]

    return X, Z


samples = []
for i in range(1,6):
    im = Image.open('image/'+str(i)+'.png')
    im_data  = np.empty((im.size[1], im.size[0]))
    for j in range(im.size[1]):
        for k in range(im.size[0]):
            R = im.getpixel((k, j))
            im_data[j,k] = R/255.0
    samples.append(im_data)

X, Z = PCA2D_2D(samples, 90, 90)

res = np.dot(Z.T, np.dot(samples[0], X))
res = np.dot(Z, np.dot(res, X.T))

row_im = Image.new('L', (res.shape[1], res.shape[0]))
y=res.reshape(1, res.shape[0]*res.shape[1])

row_im.putdata([int(t*255) for t in y[0].tolist()])
row_im.save('X.png')

在實現過程中遇到了一個小坑，特此記錄一下。PIL圖片中的size是圖片的（寬，高）二元組，而numpy中的array的shape則是矩陣的（行數，列數）二元組，索引的時候也是按照這種二元組來索引的。

在從圖片讀取畫素然後轉存到array中的時候很容易混淆這兩個概念，實際上矩陣的行數相當於圖片中的高，列數相當於寬，二者的概念剛好顛倒，很容易出錯。

用下面5張圖片做了實驗：

下面分別是X，Z矩陣分別取前5，10，20，40個特徵向量時，圖片還原的結果：

可以看出用前40個特徵向量時，已經能夠還原出很細緻的原圖了，原圖解析度為96*118，就是說，經過雙向2D PCA後，用40*40的特徵就能很好地表達96*118的資料，極大的降低了維度。