第16章 主成分分析
本文是李航老師的《統計學習方法》一書的代碼複現。作者：黃海廣
備注：代碼都可以在github中下載下傳。我将陸續将代碼釋出在公衆号“機器學習初學者”，可以在這個專輯線上閱讀。
1.假設
為
維随機變量，其均值為
，協方差矩陣為
。
考慮由
維随機變量
到
維随機變量
的線性變換
其中
。
如果該線性變換滿足以下條件，則稱之為總體主成分：
（1）
；
（2）
;
（3）變量
是
的所有線性變換中方差最大的；
是與
不相關的
的所有線性變換中方差最大的；一般地，
是與都不相關的
的所有線性變換中方差最大的；這時分别稱
為
的第一主成分、第二主成分、…、第
主成分。
2.假設
是
維随機變量，其協方差矩陣是
，
的特征值分别是
，特征值對應的機關特征向量分别是
，則
的第2主成分可以寫作
并且，
的第
主成分的方差是協方差矩陣
的第
個特征值，即
3.主成分有以下性質：
主成分
的協方差矩陣是對角矩陣
主成分
的方差之和等于随機變量
的方差之和
其中
是
的方差，即協方差矩陣
的對角線元素。
主成分
與變量
的相關系數
稱為因子負荷量（factor loading），它表示第
個主成分
與變量
的相關關系，即
對
的貢獻程度。
4.樣本主成分分析就是基于樣本協方差矩陣的主成分分析。
給定樣本矩陣
其中
是
的第
個獨立觀測樣本，
，
。
的樣本協方差矩陣
給定樣本資料矩陣
，考慮向量
到
的線性變換
這裡
如果該線性變換滿足以下條件，則稱之為樣本主成分。樣本第一主成分
是在
條件下，使得
的樣本方差
最大的
的線性變換；
樣本第二主成分
是在
和
與
的樣本協方差
條件下，使得
的樣本方差
最大的
的線性變換；
一般地，樣本第
主成分
是在
和
與
的樣本協方差
條件下，使得
的樣本方差
最大的
的線性變換。
5.主成分分析方法主要有兩種，可以通過相關矩陣的特征值分解或樣本矩陣的奇異值分解進行。
（1）相關矩陣的特征值分解算法。針對
樣本矩陣
，求樣本相關矩陣
再求樣本相關矩陣的
個特征值和對應的機關特征向量，構造正交矩陣
的每一列對應一個主成分，得到
樣本主成分矩陣
（2）矩陣
的奇異值分解算法。針對
樣本矩陣
對矩陣
進行截斷奇異值分解，保留
個奇異值、奇異向量，得到
的每一列對應一個主成分，得到
樣本主成分矩陣
本章代碼直接使用Coursera機器學習課程的第六個程式設計練習。
PCA（principal components analysis）即主成分分析技術旨在利用降維的思想，把多名額轉化為少數幾個綜合名額。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sb
from scipy.io import loadmat data = loadmat('data/ex7data1.mat')
# data X = data['X']
fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(X[:, 0], X[:, 1])
plt.show() 
PCA的算法相當簡單。在確定資料被歸一化之後，輸出僅僅是原始資料的協方差矩陣的奇異值分解。
def pca(X):
# normalize the features
X = (X - X.mean()) / X.std()
# compute the covariance matrix
X = np.matrix(X)
cov = (X.T * X) / X.shape[0]
# perform SVD
U, S, V = np.linalg.svd(cov)
return U, S, V U, S, V = pca(X)
U, S, V (matrix([[-0.79241747, -0.60997914],
[-0.60997914, 0.79241747]]),
array([1.43584536, 0.56415464]),
matrix([[-0.79241747, -0.60997914],
[-0.60997914, 0.79241747]])) 現在我們有主成分（矩陣U），我們可以用這些來将原始資料投影到一個較低維的空間中。對于這個任務，我們将實作一個計算投影并且僅選擇頂部K個分量的函數，有效地減少了維數。
def project_data(X, U, k):
U_reduced = U[:,:k]
return np.dot(X, U_reduced) Z = project_data(X, U, 1)
Z matrix([[-4.74689738],
[-7.15889408],
[-4.79563345],
[-4.45754509],
[-4.80263579],
[-7.04081342],
[-4.97025076],
[-8.75934561],
[-6.2232703 ],
[-7.04497331],
[-6.91702866],
[-6.79543508],
[-6.3438312 ],
[-6.99891495],
[-4.54558119],
[-8.31574426],
[-7.16920841],
[-5.08083842],
[-8.54077427],
[-6.94102769],
[-8.5978815 ],
[-5.76620067],
[-8.2020797 ],
[-6.23890078],
[-4.37943868],
[-5.56947441],
[-7.53865023],
[-7.70645413],
[-5.17158343],
[-6.19268884],
[-6.24385246],
[-8.02715303],
[-4.81235176],
[-7.07993347],
[-5.45953289],
[-7.60014707],
[-4.39612191],
[-7.82288033],
[-3.40498213],
[-6.54290343],
[-7.17879573],
[-5.22572421],
[-4.83081168],
[-7.23907851],
[-4.36164051],
[-6.44590096],
[-2.69118076],
[-4.61386195],
[-5.88236227],
[-7.76732508]]) 我們也可以通過反向轉換步驟來恢複原始資料。
def recover_data(Z, U, k):
U_reduced = U[:,:k]
return np.dot(Z, U_reduced.T) X_recovered = recover_data(Z, U, 1)
X_recovered matrix([[3.76152442, 2.89550838],
[5.67283275, 4.36677606],
[3.80014373, 2.92523637],
[3.53223661, 2.71900952],
[3.80569251, 2.92950765],
[5.57926356, 4.29474931],
[3.93851354, 3.03174929],
[6.94105849, 5.3430181 ],
[4.93142811, 3.79606507],
[5.58255993, 4.29728676],
[5.48117436, 4.21924319],
[5.38482148, 4.14507365],
[5.02696267, 3.8696047 ],
[5.54606249, 4.26919213],
[3.60199795, 2.77270971],
[6.58954104, 5.07243054],
[5.681006 , 4.37306758],
[4.02614513, 3.09920545],
[6.76785875, 5.20969415],
[5.50019161, 4.2338821 ],
[6.81311151, 5.24452836],
[4.56923815, 3.51726213],
[6.49947125, 5.00309752],
[4.94381398, 3.80559934],
[3.47034372, 2.67136624],
[4.41334883, 3.39726321],
[5.97375815, 4.59841938],
[6.10672889, 4.70077626],
[4.09805306, 3.15455801],
[4.90719483, 3.77741101],
[4.94773778, 3.80861976],
[6.36085631, 4.8963959 ],
[3.81339161, 2.93543419],
[5.61026298, 4.31861173],
[4.32622924, 3.33020118],
[6.02248932, 4.63593118],
[3.48356381, 2.68154267],
[6.19898705, 4.77179382],
[2.69816733, 2.07696807],
[5.18471099, 3.99103461],
[5.68860316, 4.37891565],
[4.14095516, 3.18758276],
[3.82801958, 2.94669436],
[5.73637229, 4.41568689],
[3.45624014, 2.66050973],
[5.10784454, 3.93186513],
[2.13253865, 1.64156413],
[3.65610482, 2.81435955],
[4.66128664, 3.58811828],
[6.1549641 , 4.73790627]]) fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(list(X_recovered[:, 0]), list(X_recovered[:, 1]))
plt.show() 
下載下傳位址
https://github.com/fengdu78/lihang-code
參考資料：
[1] 《統計學習方法》: https://baike.baidu.com/item/統計學習方法/10430179
[2] 黃海廣: https://github.com/fengdu78
[3] github: https://github.com/fengdu78/lihang-code
[4] wzyonggege: https://github.com/wzyonggege/statistical-learning-method
[5] WenDesi: https://github.com/WenDesi/lihang_book_algorithm