PCA原理及代碼實現

PCA(Principle Component Analysis)主成分分析是廣泛使用的降維算法，由PCA的名字就可以知道，PCA的主要目標是把數據維度降下來，使得減少數據冗余，降低數據處理帶來的計算資源消耗。

1 PCA原理

PCA的基本思想是將數據的最主要成分提取出來代替原始數據，也就是將\(n\)維特征映射到，由\(k\)維正交特征組成的特征空間就是主成分，這里使用的降維方法就是投影。問題是怎樣抽取數據的主要成分，如何衡量投影后保存的信息呢？PCA算法使用方差來度量信息量，為了確保降維后的低維度數據盡可能多的保留原始數據的有效信息，需要使降維后的數據盡可能的分散，從方差角度理解就是保留最大的方差。那么如何得到包含最大差異性的主成分呢？實際上，計算數據矩陣的協方差矩陣，得到協方差矩陣的特征值和特征向量，然后選擇特征值最大的\(k\)個特征對應的特征向量組成的矩陣，就將原始數據矩陣投影到了新的\(k\)維特征空間，實現了數據特征的降維。下面介紹方差和協方差的計算過程，關於特征值和特征向量的計算查看SVD原理。
樣本均值：$$\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_{i}$$
樣本方差：$$S^{2}=\frac{1}{n-1}\sum_{i=1}^{n}\left ( x_{i}-\bar{x} \right )^{2}$$
樣本協方差：$$Cov\left ( X,Y \right )=E\left [ \left ( X-E(X) \right ) \left ( Y-E(Y) \right )\right ]=\frac{1}{n-1}\sum_{i=1}^{n}\left ( x_{i} -\bar{x}\right )\left ( y_{i}-\bar{y} \right )$$

2 PCA算法

輸入： \(n\)維數據集\(D=\left \{ x^{(1)},x^{(2)},\cdots ,x^{(m)} \right \}\)，降維到\(k\)維
輸出： 降維后的數據集\({D}'\)
1）對所有的樣本數據去中心化：\(x^{(i)}=x^{(i)}-\frac{1}{m}\sum_{j=1}^{m}x^{(j)}\)
2）計算數據集的協方差矩陣\(XX^{T}\)
3）分解協方差矩陣\(XX^{T}\)得到特征值和特征向量
4）取出最大的\(k\)個特征值對象的特征向量\(\left ( w_{1},w_{2},\cdots ,w_{k} \right )\)，將所有特征向量標准化得到特征向量矩陣\(W\)
5）對數據集中的每一個樣本\(x^{(i)}\)，轉換為新的樣本\(z^{(i)}=W^{T}x^{(i)}\)
6）得到輸出樣本集\({D}'=\left ( z^{(1)},z^{(2)},\cdots ,z^{(m)} \right )\)

3 PCA代碼實現

PCA降維

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat

#2D-->1D
mat = loadmat('D:/Python/Andrew-NG-Meachine-Learning/machine-learning-ex7/ex7/ex7data1.mat')
X = mat['X']
print(X.shape)   #(50, 2)
plt.scatter(X[:,0], X[:,1], facecolors='none', edgecolors='b')

#X 均值歸一化
def featureNormal(X):
    means = X.mean(axis=0)
    stds = X.std(axis=0, ddof=1)
    X_norm = (X - means)/stds
    return X_norm, means, stds


#PCA
def pca(X):
    sigma = (X.T@X)/len(X)
    U, S, V = np.linalg.svd(sigma)
    return U, S, V


X_norm, means, stds = featureNormal(X)
U, S, V = pca(X_norm)


print(U[:,0])


plt.figure(figsize=(7, 5))
plt.scatter(X[:,0], X[:,1], facecolors='none', edgecolors='b')
plt.plot([means[0], means[0]+1.5*S[0]*U[0,0]],
         [means[1], means[1]+1.5*S[0]*U[0,1]],
         c='r', linewidth=3, label='First Principal Component')


plt.plot([means[0], means[0]+1.5*S[1]*U[1,0]],
         [means[1], means[1]+1.5*S[1]*U[1,1]],
         c='g', linewidth=3, label='Second Principal Component')  
plt.grid()
plt.axis("equal")
plt.legend()      


#Dimensionality Reduction with PCA


def projectData(X, U, K):
    Z = X @ U[:,:K]
    return Z


Z = projectData(X_norm, U, 1)
Z[0]
#print(Z[0]) #[ 1.48127391]


#Reconstructing an approximation of the data 重建數據
def recoverData(Z, U, K):
    X_rec = Z @ U[:,:K].T
    return X_rec
           
X_rec = recoverData(Z, U, 1)
X_rec[0]  
#print(X_rec[0])     #[-1.04741883 -1.04741883]        


#Visualizing the projections
plt.figure(figsize=(7,5))
plt.axis("equal")
plot = plt.scatter(X_norm[:,0], X_norm[:,1], s=30, facecolors='none',
                   edgecolors='b',label='Original Data Points')


plot = plt.scatter(X_rec[:,0], X_rec[:,1], s=30, facecolors='none',
                   edgecolors='r',label='PCA Reduced Data Points')


plt.title("Example Dataset: Reduced Dimension Points Shown", fontsize=14)
plt.xlabel('x1 [Feature Normalized]',fontsize=14)
plt.ylabel('x2 [Feature Normalized]', fontsize=14)
plt.grid(True)


for x in range(X_norm.shape[0]):
    plt.plot([X_norm[x,0], X_rec[x,0]],[X_norm[x,1], X_rec[x,1]], 'k--')
    #輸入第一項全是X坐標 第二項全是y坐標
plt.legend()

可視化 PCA將數據從2D減少到1D：

PCA應用 Face Image Dataset 人臉識別圖像上運行PCA 實踐中使用PCA減少維度

大數據集實現PCA

import numpy as np
import pandas as pd
from scipy.io import loadmat
import matplotlib.pyplot as plt

mat = loadmat('D:/Python/ex7faces.mat')
X = mat['X']
print(X.shape)  #(5000, 1024)


def displayData(X, row, col):
    fig, axs = plt.subplots(row, col, figsize=(8,8))
    for r in range(row):
        for c in range(col):
            axs[r][c].imshow(X[r*col + c].reshape(32,32).T, cmap = 'Greys_r')
            axs[r][c].set_xticks([])
            axs[r][c].set_yticks([])
            
displayData(X, 10, 10)

def featureNormalize(X):
    means = X.mean(axis=0)
    stds = X.std(axis=0, ddof=1)
    X_norm = (X - means) / stds
    return X_norm, means, stds


def pca(X):
    sigma = (X.T @ X) / len(X)
    U, S, V = np.linalg.svd(sigma)
    
    return U, S, V    

X_norm, means, stds = featureNormalize(X)
U, S, V = pca(X_norm)    
#print(U.shape, S.shape)  #(1024, 1024) (1024,)

displayData(U[:,:36].T, 6, 6)

#Dimensionality Reduction

def projectData(X, U, K):
    Z = X @ U[:,:K]
    
    return Z

z = projectData(X_norm, U, K=36)
def recoverData(Z, U, K):
    X_rec = Z @ U[:,:K].T
    
    return X_rec
X_rec = recoverData(z, U, K=36)
displayData(X_rec, 10, 10)

參考： 吳恩達機器學習