主成分分析(PCA)
假设xxx为mmm 维随机变量,其均值为μ\muμ,协方差矩阵为Σ\SigmaΣ
考虑由mmm维随机变量xxx到mmm维随机变量yyy的线性变换
yi=αiTx=∑k=1mαkixk,i=1,2,⋯,my _ { i } = \alpha _ { i } ^ { T } x = \sum _ { k = 1 } ^ { m } \alpha _ { k i } x _ { k } , \quad i = 1,2 , \cdots , myi=αiTx=k=1∑mαkixk,i=1,2,⋯,m
其中αiT=(α1i,α2i,⋯,αmi)\alpha _ { i } ^ { T } = ( \alpha _ { 1 i } , \alpha _ { 2 i } , \cdots , \alpha _ { m i } )αiT=(α1i,α2i,⋯,αmi)
如果该线性变换满足以下条件,则称之为总体主成分:
αiTαi=1,i=1,2,⋯,m\alpha _ { i } ^ { T } \alpha _ { i } = 1 , i = 1,2 , \cdots , mαiTαi=1,i=1,2,⋯,m;
cov(yi,yj)=0(i≠j)\operatorname { cov } ( y _ { i } , y _ { j } ) = 0 ( i \neq j )cov(yi,yj)=0(i=j);
变量y1y_1y1是xxx的所有线性变换中方差最大的
y2y_2y2是与y1y_1y1不相关的xxx的所有线性变换中方差最大的
一般地,yiy_iyi是与y1,y2,⋯,yi−1,(i=1,2,⋯,m)y _ { 1 } , y _ { 2 } , \cdots , y _ { i - 1 } , ( i = 1,2 , \cdots , m )y1,y2,⋯,yi−1,(i=1,2,⋯,m)都不相关的xxx的所有线性变换中方差最大的
这时分别称y1,y2,⋯,ymy _ { 1 } , y _ { 2 } , \cdots , y _ { m }y1,y2,⋯,ym为xxx的第一主成分、第二主成分、…、第mmm主成分
假设xxx是mmm维随机变量
其协方差矩阵是Σ\SigmaΣ
Σ\SigmaΣ的特征值分别是λ1≥λ2≥⋯≥λm≥0\lambda _ { 1 } \geq\lambda _ { 2 } \geq \cdots \geq \lambda _ { m } \geq 0λ1≥λ2≥⋯≥λm≥0
特征值对应的单位特征向量分别是α1,α2,⋯,αm\alpha _ { 1 } , \alpha _ { 2 } , \cdots , \alpha _ { m }α1,α2,⋯,αm,则xxx的第2主成分可以写作
yi=αiTx=∑k=1mαkixk,i=1,2,⋯,my _ { i } = \alpha _ { i } ^ { T } x = \sum _ { k = 1 } ^ { m } \alpha _ { k i } x _ { k } , \quad i = 1,2 , \cdots , myi=αiTx=k=1∑mαkixk,i=1,2,⋯,m
并且,xxx的第iii主成分的方差是协方差矩阵Σ\SigmaΣ的第iii个特征值,即var(yi)=αiTΣαi=λi\operatorname { var } ( y _ { i } ) = \alpha _ { i } ^ { T } \Sigma \alpha _ { i } = \lambda _ { i }var(yi)=αiTΣαi=λi
主成分有以下性质:
主成分yyy的协方差矩阵是对角矩阵cov(y)=Λ=diag(λ1,λ2,⋯,λm)\operatorname { cov } ( y ) = \Lambda = \operatorname { diag } ( \lambda _ { 1 } , \lambda _ { 2 } , \cdots , \lambda _ { m } )cov(y)=Λ=diag(λ1,λ2,⋯,λm)
主成分yyy的方差之和等于随机变量xxx的方差之和
∑i=1mλi=∑i=1mσii\sum _ { i = 1 } ^ { m } \lambda _ { i } = \sum _ { i = 1 } ^ { m } \sigma _ { i i }i=1∑mλi=i=1∑mσii
其中σii\sigma _ { i i }σii是x2x_2x2的方差,即协方差矩阵Σ\SigmaΣ的对角线元素
主成分yky_kyk与变量x2x_2x2的相关系数ρ(yk,xi)\rho ( y _ { k } , x _ { i } )ρ(yk,xi)称为因子负荷量(factor loading)
它表示第kkk个主成分yky_kyk与变量xxx的相关关系,即yky_kyk对xxx的贡献程度
ρ(yk,xi)=λkαikσii,k,i=1,2,⋯,m\rho ( y _ { k } , x _ { i } ) = \frac { \sqrt { \lambda _ { k } } \alpha _ { i k } } { \sqrt { \sigma _ { i i } } } , \quad k , i = 1,2 , \cdots , mρ(yk,xi)=σiiλkαik,k,i=1,2,⋯,m
样本主成分分析就是基于样本协方差矩阵的主成分分析
给定样本矩阵
X=[x1x2⋯xn]=[x11x12⋯x1nx21x22⋯x2n⋮⋮⋮xm1xm2⋯xmn]X = \left[ \begin{array} { l l l l } { x _ { 1 } } & { x _ { 2 } } & { \cdots } & { x _ { n } } \end{array} \right] = \left[ \begin{array} { c c c c } { x _ { 11 } } & { x _ { 12 } } & { \cdots } & { x _ { 1 n } } \\\\ { x _ { 21 } } & { x _ { 22 } } & { \cdots } & { x _ { 2 n } } \\\\ { \vdots } & { \vdots } & { } & { \vdots } \\\\ { x _ { m 1 } } & { x _ { m 2 } } & { \cdots } & { x _ { m n } } \end{array} \right]X=[x1x2⋯xn]=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡x11x21⋮xm1x12x22⋮xm2⋯⋯⋯x1nx2n⋮xmn⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
其中xj=(x1j,x2j,⋯,xmj)Tx _ { j } = ( x _ { 1 j } , x _ { 2 j } , \cdots , x _ { m j } ) ^ { T }xj=(x1j,x2j,⋯,xmj)T是xxx的第jjj个独立观测样本,j=1,2,…,nj=1,2,…,nj=1,2,…,n。
XXX的样本协方差矩阵
S=[sij]m×m,sij=1n−1∑k=1n(xik−x‾i)(xjk−x‾j)i=1,2,⋯,m,j=1,2,⋯,m\left. \begin{array} { c } { S = [ s _ { i j } ] _ { m \times m } , \quad s _ { i j } = \frac { 1 } { n - 1 } \sum _ { k = 1 } ^ { n } ( x _ { i k } - \overline { x } _ { i } ) ( x _ { j k } - \overline { x } _ { j } ) } \\\\ { i = 1,2 , \cdots , m , \quad j = 1,2 , \cdots , m } \end{array} \right.S=[sij]m×m,sij=n−11∑k=1n(xik−xi)(xjk−xj)i=1,2,⋯,m,j=1,2,⋯,m
给定样本数据矩阵XXX,考虑向量xxx到yyy的线性变换y=ATxy = A ^ { T } xy=ATx
这里
A=[a1a2⋯am]=[a11a12⋯a1ma21a22⋯a2m⋮⋮⋮am1am2⋯amm]A = \left[ \begin{array} { l l l l } { a _ { 1 } } & { a _ { 2 } } & { \cdots } & { a _ { m } } \end{array} \right] = \left[ \begin{array} { c c c c } { a _ { 11 } } & { a _ { 12 } } & { \cdots } & { a _ { 1 m } } \\\\ { a _ { 21 } } & { a _ { 22 } } & { \cdots } & { a _ { 2 m } } \\\\ { \vdots } & { \vdots } & { } & { \vdots } \\\\ { a _ { m 1 } } & { a _ { m 2 } } & { \cdots } & { a _ { m m } } \end{array} \right]A=[a1a2⋯am]=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡a11a21⋮am1a12a22⋮am2⋯⋯⋯a1ma2m⋮amm⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
如果该线性变换满足以下条件,则称之为样本主成分。样本第一主成分y1=a1Txy _ { 1 } = a _ { 1 } ^ { T } xy1=a1Tx是在a1Ta1=1a _ { 1 } ^ { T } a _ { 1 } = 1a1Ta1=1条件下
使得a1Txj(j=1,2,⋯,n)a _ { 1 } ^ { T } x _ { j } ( j = 1,2 , \cdots , n )a1Txj(j=1,2,⋯,n)的样本方差a1TSa1a _ { 1 } ^ { T } S a _ { 1 }a1TSa1最大的xxx的线性变换
样本第二主成分y2=a2Txy _ { 2 } = a _ { 2 } ^ { T } xy2=a2Tx
是在a2Ta2=1a _ { 2 } ^ { T } a _ { 2 } = 1a2Ta2=1和a2Txja _ { 2 } ^ { T } x _ { j }a2Txj与a1Txj(j=1,2,⋯,n)a _ { 1 } ^ { T } x _ { j } ( j = 1,2 , \cdots , n )a1Txj(j=1,2,⋯,n)的样本协方差
a1TSa2=0a _ { 1 } ^ { T } S a _ { 2 } = 0a1TSa2=0条件下
使得a2Txj(j=1,2,⋯,n)a _ { 2 } ^ { T } x _ { j } ( j = 1,2 , \cdots , n )a2Txj(j=1,2,⋯,n)的样本方差a2TSa2a _ { 2 } ^ { T } S a _ { 2 }a2TSa2最大的xxx的线性变换
一般地,样本第iii主成分yi=aiTxy _ { i } = a _ { i } ^ { T } xyi=aiTx是在aiTai=1a _ { i } ^ { T } a _ { i } = 1aiTai=1和aiTxja _ { i } ^ { T } x _ { j }aiTxj与akTxj(k
akTSai=0a _ { k } ^ { T } S a _ { i } = 0akTSai=0条件下
使得aiTxj(j=1,2,⋯,n)a _ { i } ^ { T } x _ { j } ( j = 1,2 , \cdots , n )aiTxj(j=1,2,⋯,n)的样本方差akTSaia _ { k } ^ { T } S a _ { i }akTSai最大的xxx的线性变换
主成分分析方法主要有两种,可以通过相关矩阵的特征值分解或样本矩阵的奇异值分解进行
相关矩阵的特征值分解算法
针对m×nm \times nm×n样本矩阵XXX,求样本相关矩阵
R=1n−1XXTR = \frac { 1 } { n - 1 } X X ^ { T }R=n−11XXT
再求样本相关矩阵的kkk个特征值和对应的单位特征向量,构造正交矩阵
V=(v1,v2,⋯,vk)V = ( v _ { 1 } , v _ { 2 } , \cdots , v _ { k } )V=(v1,v2,⋯,vk)
VVV的每一列对应一个主成分,得到k×nk \times nk×n样本主成分矩阵
Y=VTXY = V ^ { T } XY=VTX
矩阵XXX的奇异值分解算法
针对m×nm \times nm×n样本矩阵XXX
X′=1n−1XTX ^ { \prime } = \frac { 1 } { \sqrt { n - 1 } } X ^ { T }X′=n−11XT
对矩阵X′X ^ { \prime }X′进行截断奇异值分解,保留kkk个奇异值、奇异向量,得到
X′=USVTX ^ { \prime } = U S V ^ { T }X′=USVT
VVV的每一列对应一个主成分,得到k×nk \times nk×n样本主成分矩阵YYY
Y=VTXY = V ^ { T } XY=VTX
PCA(principal components analysis)即主成分分析技术旨在利用降维的思想,把多指标转化为少数几个综合指标。
PCA的算法相当简单。 在确保数据被归一化之后,输出仅仅是原始数据的协方差矩阵的奇异值分解。
现在我们有主成分(矩阵U),我们可以用这些来将原始数据投影到一个较低维的空间中
对于这个任务,我们将实现一个计算投影并且仅选择顶部K个分量的函数,有效地减少了维数
我们也可以通过反向转换步骤来恢复原始数据。
第一主成分的投影轴基本上是数据集中的对角线
当我们将数据减少到一个维度时,我们失去了该对角线周围的变化,所以在我们的再现中,一切都沿着该对角线
代码实现
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sb
from scipy.io import loadmatX = data['X']
def pca(X):# normalize the featuresX = (X - X.mean()) / X.std()# compute the covariance matrixX = np.matrix(X)cov = (X.T * X) / X.shape[0]# perform SVDU, S, V = np.linalg.svd(cov)return U, S, V
def project_data(X, U, k):U_reduced = U[:,:k]return np.dot(X, U_reduced)
def recover_data(Z, U, k):U_reduced = U[:,:k]return np.dot(Z, U_reduced.T)
U, S, V = pca(X)
Z = project_data(X, U, 1)
X_recovered = recover_data(Z, U, 1)
fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(list(X_recovered[:, 0]), list(X_recovered[:, 1]))
plt.show()

