主成分分析（PCA）學習筆記

　　這兩天學習了吳恩達老師機器學習中的主成分分析法（Principal Component Analysis， PCA)，PCA是一種經常使用的降維方法。這裏對PCA算法作一個小筆記，並利用python完成對應的練習(ps:最近公式有點多，開始沒找到怎麼敲公式，前面幾篇都是截的圖^_^，後面問了度娘，原來是支持latex的)。代碼和數據見github

1、PCA基本思路

　　將數據從原來的座標系轉換到新的座標系，新座標系的選擇由數據自己決定。第一個新座標軸選擇的是原始數據中方差最大的方向，第二個新座標軸選擇的是和第一個座標軸正交且具備最大方差的方向，依次類推，直到找出\(k\)個新座標軸。也就是將原始數據投影到一個低維的座標系中。

2、PCA目標函數

　　以最小化投影偏差爲目標函數，這裏注意與線性迴歸的區別，線性迴歸是最小化垂直距離。下圖左邊圖中藍線是線性迴歸的目標，右圖中藍線是PCA的目標。

3、PCA算法步驟

一、數據預處理，數據維度爲 $ m*n $ 。

　　對 $ X= x^{(1)}, x^{(2)}, ... , x^{(m)} $ 計算平均值和方差，\[ u_j = \frac{1}{m}\sum_{i=1}^{m}x_j^{(i)} \] \[ s_j = \frac{1}{m}\sum_{i=1}^{m}(x_j^{(i)}-u_j)^2 \]
對原始數據進行歸一化處理獲得 \[ x_j = \frac{x_j-u_j}{s_j} \]
其中 $ u_j $ 表示均值， $ s_j $ 表示方差。

二、計算協方差矩陣 $ \Sigma $ 。

\[ \Sigma = \frac{1}{m} \sum_{i=1}^{m}(x^{(i)}.(x^{(i)})^T) \]
　　其中 $ x^{(i)} $ 是 $ n \times 1 $ 維的，則其轉秩是$ 1 \times n $ 維的，因此 $ \Sigma $ 是 $ n \times n $ 維的。

三、對 $ \Sigma $ 進行奇異值分解。

\[ U,S,V = svd(\Sigma) \]
奇異值分解能夠參考博客，博主講的比較清楚。上式中
\[ U = [u^{(1)}, u^{(2)},...,u^{(n)}]\in R^{n*n} \]
從中選取 $ k $ 個主要成分
\[ U_k = [u^{(1)}, u^{(2)},...,u^{(k)}] \]
則 $ x^{(i)} $ 在新座標上的投影能夠表示爲
\[ z^{(i)} = [u^{(1)}, u^{(2)},...,u^{(k)}]^T\cdot x^{(i)}= \begin{bmatrix} u^{(1)}\\ u^{(2)}\\ \vdots \\ u^{(k)}\\ \end{bmatrix}\cdot x^{(i)} \]
其中 $ U_k $ 是 $ k \times n $ 維的，其轉秩爲 $ n \times k $ 維的，$ x^{(i)} $ 是 $ n \times 1 $ 維，因此 $ z^{(i)} $ 是 $ k \times 1 $ 維的。投影 $ Z $ 能夠表示爲
\[ Z = U_k^T \cdot X \]

四、計算重構後特徵 $ x_{approx}^{(i)} $ 。

\[ x^{(i)}_{approx} = U_k \cdot z^{(i)} \]
其中 $ U_k $ 是 $ n \times k $ 維，$ z^{(i)} $ 是 $ k \times 1 $ 維，則 $ x_{approx}^{(i)} $ 是 $ n \times 1 $ 維。

五、根據投影偏差檢驗選擇的 $ k $ 個主要成分是否知足要求。

\[ \frac{\frac{1}{m}\sum_{i=1}^{m}||x^{(i)}-x^{(i)}_{approx}||^2}{\frac{1}{m}\sum_{i=1}^{m}||x^{(i)}||^2} \leq 0.01 \]
若是不知足上式，則從步驟3中從新選擇 $ k $ 個主成分，繼續第4和第5步，直到知足要求爲止。小於等於0.01代表保留了原始數據99%信息，這裏能夠根據需求進行更改。
這裏若是不用奇異值分解後的 $ U $ 矩陣，也能夠根據奇異值矩陣 $ S $ 計算 ，$ S $ 是主對角線上爲從大到小排列的奇異值，其餘元素全爲0的對角矩陣。
\[ 1-\frac{\sum_{i=1}^{k}S_{ii}}{\sum_{i=1}^{n}S_{ii}}\leq 0.01 \]

PCA應用實例

一、二維數據投影到一維，熟悉PCA流程

第一步 引入相關庫並導入數據

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

data1 = loadmat('./ex7/ex7data1.mat')
data1

X = data1['X']
X.shape

原始數據展現

plt.scatter(X[:,0], X[:,1])
plt.show()

第二步 數據預處理

# 定義歸一化函數featureNormalize
def featureNormalizse(x):
    mean = x.mean(axis=0)
    std = x.std(axis=0)
    return (x-mean)/std, mean, std

# test
x_norm, means, stds = featureNormalizse(X)
x_norm[:5]

第三步 計算協方差矩陣

print(x_norm.shape)
sigma = (x_norm.T.dot(x_norm))/x_norm.shape[0]
sigma

第四步 奇異值分解

U,S,V = np.linalg.svd(sigma)
U,S,V

# 整合sigma和svd
def pca(x):
    sigma = (x.T.dot(x))/x.shape[0]
    U,S,V = np.linalg.svd(sigma)
    return U,S,V

可視化主成分

x_norm, means, stds = featureNormalizse(X)
U, S, V = pca(x_norm)

plt.figure(figsize=(8, 6))
plt.scatter(X[:,0], X[:,1], label='sample data')  # 樣本數據點

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()
plt.show()

第五步 計算 $ x_{approx} $

# 根據U_reduce計算x_norm的投影Z
def compute_z(X, U, k):
    Z = X.dot(U[:,:k])
    return Z

# test
Z = compute_z(x_norm, U, 1)
Z[:5]

# 計算x_approx
def compute_x_approx(U, k, Z):
    x_approx = Z.dot(U[:,:k].T)  # 50*1 * 1*2
    return x_approx

# test
x_approx = compute_x_approx(U, 1, Z)
x_approx[:5]

第六步 可視化投影效果

plt.figure(figsize=(8,6))
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_approx[:,0], x_approx[:,1], s=30, facecolors='none', 
                   edgecolors='r',label='PCA Reduced Data Points')  # 畫出通過PCA後構造的估計值

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_approx[x,0]],[x_norm[x,1],x_approx[x,1]],'k--')
    # 輸入第一項全是X座標，第二項都是Y座標
plt.legend()
plt.show()

二、利用PCA減小圖片維度

第一步 導入數據

# 導入數據
face_data = loadmat('./ex7/ex7faces.mat')
face = face_data['X']
face.shape

第二步 展現數據

def showFace(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([])
            
showFace(face, 10, 10)
plt.show()

第三步 降維並重構

face_norm, means, stds = featureNormalizse(face) # 歸一化
U, S, V = pca(face_norm)  # 奇異值分解
Z = compute_z(face_norm, U, 16)  # 計算Z
face_approx = compute_x_approx(U, 16, Z) # 計算降維重構後的數據
face_approx.shape

第四步 先後對比圖

showFace(face, 10, 10)
showFace(face_approx, 10, 10)
plt.show()