利用PCA進行資料降維

Babyface Killer發表於2020-11-10
PCA

PCA原理


在介紹PCA之前首先要熟悉一下數學推導過程。

特徵多項式：

設A為一個方陣，則該方陣的特徵多項式就為該方陣減去\lambda倍的單位矩陣後構成的矩陣的行列式。而該多項式的所有解即為\lambda的值，也就是該方陣的特徵值。

解得特徵值之後如何求得特徵向量：

找到特徵值後，根據上式定義我們可推出Ax-\lambda x=0,即(A-\lambda )x=0，該式中x即為特徵向量，A與\lambda已知，解出上式即可求出方陣A的特徵向量，且每個特徵值\lambda對應一個特徵向量。

特徵分解：

對於任一方陣A，其可根據上式被分解，P為由方陣的特徵向量構成的方陣，D為一對角矩陣，其中當方陣A的特徵向量是n維的基向量時對角矩陣D的值為方陣A的特徵值。

注：特徵分解和奇異值分解的區別

特徵分解公式：

奇異值分解公式：

  1. SVD對任意矩陣都適用，而特徵分解只適用於方陣
  2. 特徵分解中的P矩陣不一定正交，而SVD中的U和V矩陣一定正交
  3. SVD中的U和V代表不同含義因此大部分情況下兩者不互為對方的逆矩陣

實踐中PCA應用的關鍵步驟：

  1. 標準化，對於資料的每一個特徵減去其對應均值併除以其對應標準差
  2. 計算協方差矩陣，C=A^{T}\cdot A
  3. 對協方差矩陣做特徵分解得到其特徵值和特徵分解（在使用大部分機器學習庫時這一步的輸出為按照特徵值由高到低排序後的結果）
  4. 根據選取的特徵向量做投影，即把標準化後的資料投影到選取的特徵向量上，Projection=A\cdot P
  5. 在下面的例子中我對投影后的資料進行了重建，即把投影后的資料轉化為和原始資料相同的維度：Reconstruction=Projection\cdot P^{T} ，這一步在大部分實踐中是不需要的，這一步的目的是為了比較使用PCA前後資料的變化

PCA的應用

介紹完了PCA的原理下面我們來看看PCA在實際應用中的效果。我使用的資料集是不同國家每天的新冠確診人數,在進行PCA之前我已經對資料集做了標準化處理。

我這裡寫了一個自定義函式，輸入為每個國家對應的字串，輸出為對應國家每日確診病例變化曲線（標準化後），依次使用不同數量的PC重建資料集的影像（max 10），使用不同數量PC重建資料與原始資料的殘差影像（max 10），使用不同數量PC的RMSE曲線（max10），RMSE降低到1，0.1，0.01分別需要的PC數量。（這裡PC為principal component，也就是特徵向量）

from sklearn.metrics import mean_squared_error
import numpy as np

def pca(country):
    #找出輸入國家對應的index
    country_index=np.where(cases_standardized.index==country)
    fig,axes=plt.subplots(4,1,figsize=[10,20])
    #提取出該國家對應的資料
    country_data=cases_standardized.iloc[country_index[0][0],:]
    axr=axes.ravel()
    #在第一張圖中畫出對應國家的每日確診病例曲線
    axr[0].plot(country_data)
    axr[0].set_title('Standardized Time Series')
    labels=[x for x in range(0,265,30)]
    axr[0].set_xticks(labels)
    axr[0].set_xticklabels(dates[labels],rotation=30)
    n,m=data_standardized.shape
    #計算協方差矩陣
    C = np.dot(data_standardized.T, data_standardized) / (n-1) 
    #使用numpy.linalg.eigh()計算協方差矩陣的特徵值和特徵向量
    eigenValues, eigenVectors = np.linalg.eigh(C) 
    #提取排序後的index
    args = (-eigenValues).argsort()
    #把特徵值和特徵向量按照從高到低的資料排序
    eigenValues = eigenValues[args]
    eigenVectors = eigenVectors[:, args]
    #提取該國家原始資料
    true_value=data_standardized[country_index]
    RMSE=[]
    #使用不同數量的特徵向量進行十次投影和重建
    for i in range(10):
        W=eigenVectors[:,:i+1]
        projX=np.dot(data_standardized,W)
        ReconX = np.dot(projX, W.T)
        result=ReconX[country_index][0]
        residual_error=true_value-result
        #計算RMSE
        RMSE.append(np.sqrt(mean_squared_error(true_value[0],result)))
        #在第二張圖中繪製重建後的資料影像
        axr[1].plot(ReconX[country_index][0],label='CumPC{}'.format(i+1))
        #在第三張圖中繪製殘差影像
        axr[2].plot(residual_error[0],label='CumPC{}'.format(i+1))
    axr[1].set_title('Cumulative Reconstruction')
    labels=[x for x in range(0,265,30)]
    axr[1].set_xticks(labels)
    axr[1].set_xticklabels(dates[labels],rotation=30)
    axr[1].legend()
    axr[2].set_title('Residual Error')
    axr[2].set_xticks(labels)
    axr[2].set_xticklabels(dates[labels],rotation=30)
    axr[2].legend()
    #在第四張圖中繪製RMSE曲線
    axr[3].plot(RMSE)
    axr[3].set_title('RMSE')
    axr[3].set_xlabel('Number of component')
    #找出使RMSE低於1需要的最少的特徵向量數
    for i in range(10):
        if RMSE[i] < 1:
            print('Number of PCs required for RMSE < 1: {}'.format(i+1))
            break
    #找出使RMSE低於0.1需要的最少的特徵向量數
    for i in range(10):
        if RMSE[i] < 0.1:
            print('Number of PCs required for RMSE < 0.1: {}'.format(i+1))
            break
    #找出使RMSE低於0.01需要的最少的特徵向量數
    for i in range(10):
        if RMSE[i] < 0.01:
            print('Number of PCs required for RMSE < 0.01: {}'.format(i+1))
            break
    #如果十個特徵向量不能滿足RMSE低於0.01的條件則繼續使用更多的特徵向量進行投影和重建直到找到能使RMSE低於0.01的特徵向量數
    if RMSE[9] >= 0.01:
        for i in range(10,265):
            W=eigenVectors[:,:i+1]
            projX=np.dot(data_standardized,W)
            ReconX = np.dot(projX, W.T)
            result=ReconX[country_index][0]
            RMSE=np.sqrt(mean_squared_error(true_value[0],result))
            if RMSE < 0.01:
                print('Number of PCs required for RMSE < 0.01: {}'.format(i+1))
                break

我們把‘China’輸入這個函式來看一下輸出。

可以看到只使用一個PC時對於原資料的還原不是很理想，但是在使用兩個PC時重建後的資料已經和原始資料非常相似了，而使用5個PC就可以使原資料與重建資料之間的RMSE低於0.01。

通過以上的例子我們可以發現，使用PCA可以把一個非常高維的資料集降為低維，在上面的例子中原資料集中包含了189個國家265天的資料，而我們只使用5個主成分就可以幾乎完美地擬合原資料集。但PCA最大的缺陷就在於把原資料集投影到主成分後資料的可解釋性就變得極低，因此在使用PCA前就需要考慮特徵的可解釋性對於之後的建模和分析到底是不是必要的。

本文中公式推導部分來自於：Mathematics for Machine Learning https://github.com/mml-book/mml-book.github.io

相關文章