主成分分析简介

已注销
230 阅读4分钟

本文主要介绍总体及样本的主成分的概念，如何可视化，以及一些最基本的性质。

1 总体的主成分

考虑x∼(μ,Σ)x\sim (\mu,\Sigma)，其中Σ\Sigma可分解为Σ=ΓΛΓ′\Sigma=\Gamma\Lambda\Gamma'，Γ=(η1,…,ηd)\Gamma=(\eta_1,\ldots,\eta_d)，各列为单位特征向量，Λ\Lambda为特征值降序排列的对角矩阵，协方差矩阵的秩rank(Σ)=r\text{rank}(\Sigma)=r，对于k=1,…,rk=1,\ldots,r，定义如下概念：

  • 第kk个principal component score（主成分得分）为wk=ηk′(x−μ)w_k=\eta_k' (x-\mu)；
  • 第kk个principal component vector（主成分向量）为w(k)=(w1,…,wk)′=Γk′(x−μ)w^{(k)}=(w_1,\ldots,w_k)'=\Gamma_k' (x-\mu)；
  • 第kk个principal component projection (vector)（主成分映射（向量））为pk=ηkηk′(x−μ)=wkηkp_k=\eta_k\eta_k'(x-\mu)=w_k\eta_k。

我们来仔细看以上几个概念。ηk\eta_k也叫载荷（loading），score xkx_k其实就是xx在方向ηk\eta_k上的贡献，也就是将xx投影到ηk\eta_k，而将前kk个score放到一起，就组成了x(k)x^{(k)}。pkp_k是一个dd维向量，它的方向就是ηk\eta_k的方向，长度/欧式范数为∣xk∣\vert x_k\vert。

2 样本的主成分

假设我们有nn个样本，将它们排成d×nd\times n的矩阵X=(x1,…,xn)X=(x_1,\ldots,x_n)，记样本均值xˉ=1nXℓn\bar x=\dfrac{1}{n}X\ell_n，样本协方差矩阵S=1n−1(X−xˉℓn′)(X−xˉℓn′)′S=\dfrac{1}{n-1}(X-\bar x \ell_n')(X-\bar x \ell_n')'，并且rank(S)=r≤d\text{rank}(S)=r\leq d。

我们可以对样本协方差矩阵进行谱分解S=Γ^Λ^Γ^′S=\hat\Gamma \hat\Lambda \hat\Gamma'。与在总体中的定义一样，可以在样本中定义如下概念：

  • 第kk个principal component score为wk=η^k′(X−xˉℓn′)w_k=\hat\eta_k' (X-\bar x\ell_n')，这是一个1×n1\times n的行向量；
  • principal component data为W(k)=(w1′,…,wk′)′=Γ^k′(X−xˉℓn′)W^{(k)}=(w_1',\ldots,w_k')'=\hat\Gamma_k' (X-\bar x\ell_n')，它是k×nk\times n的矩阵；
  • 第kk个principal component projection为Pk=η^kη^k′(X−xˉℓn′)=η^kwkP_k=\hat\eta_k\hat\eta_k'(X-\bar x\ell_n')=\hat\eta_k w_k，它是d×nd\times n的矩阵。

3 主成分的可视化

本节考察如何将PC的贡献、PC scores W(k)W^{(k)}、PC projection PkP_k进行可视化。

3.1 Scree plot

对于总体数据，rank(Σ)=r\text{rank}(\Sigma)=r，我们定义第kk个PC，对总体协方差的贡献为

λk∑j=1rλj=λktr(Σ)\dfrac{\lambda_k}{\sum_{j=1}^{r}\lambda_j}=\dfrac{\lambda_k}{\text{tr}(\Sigma)}

再定义前κ\kappa个PC对总体协方差的累积贡献为

∑k=1κλk∑j=1rλj=∑k=1κλktr(Σ)\dfrac{\sum_{k=1}^{\kappa}\lambda_k}{\sum_{j=1}^{r}\lambda_j}=\dfrac{\sum_{k=1}^{\kappa}\lambda_k}{\text{tr}(\Sigma)}

对于不同kk画出前kk个PC的贡献和累积贡献，就得到了scree plot，如下图：

image.png

Scree plot中，有时可以找到elbow，传统方法是找出能比较好地代表数据的κ\kappa个点，就是elbow在的地方。有些文献中也叫knee或kink。在上图的例子中，elbow就是在第33个特征值处出现。当然，在很多情况下，也可能没有elbow。

3.2 PC score plot

如果我们画出各PC的score，就叫PC score plot，但由于我们的视觉只能接受到三维的事物，因此对于超过三维，需要进行处理。比如对于W(4)W^{(4)}，我们可以成对地在二维平面上画出各PC的score，如下图：

image.png

也可以直接画出三维图，如下图：

image.png

3.3 Projection plot

接下来看如何将PkP_k可视化。PkP_k为d×nd\times n的矩阵，它的第ii列就是第ii个样本在η^k\hat\eta_k方向的贡献。对于每个kk，都可以用PkP_k画出平行坐标图，这样就很方便地将PkP_k可视化了。

第kk个PC projection表示的其实就是用第kk个score wkw_k对向量η^k\hat\eta_k进行伸缩变换，因此也可以看一下这些score的分布，这就是密度估计（density estimates）。

下图是一个d=30d=30的示例，左图为某个PC的projection plot，右图为非参数密度估计（可用Matlab工具curvdatSM）：

image.png

4 PC的性质

4.1 结构性质

本节考虑PC的一些性质，先看XX和它的PC的相关结构。这里采用第1节的设定。

首先来看w(k)=Γk′(x−μ)w^{(k)}=\Gamma_k' (x-\mu)，它的期望必满足E(w(k))=0\text{E}(w^{(k)})=0，这说明PC vectors都是中心化了的。而它的方差

Var(w(k))=Γk′ΣΓk=Γk′ΓΛΓ′Γk=Ik×dΛId×k=Λk\begin{aligned} \text{Var}(w^{(k)})=&\Gamma_k' \Sigma\Gamma_k\\ =&\Gamma_k' \Gamma\Lambda\Gamma'\Gamma_k\\ =&I_{k\times d}\Lambda I_{d\times k}\\ =&\Lambda_k \end{aligned}

即PC scores互不相关。当然，如果xx是正态分布，那么PC scores也相互独立，否则需要用其他方法做出相互独立。

而w(k)w^{(k)}是由wk=ηk′(x−μ)w_k=\eta_k' (x-\mu)作为元素组成的列向量，因此Var(wk)=λk\text{Var}(w_k)=\lambda_k，Cov(wk,wl)=ηk′Σηl=λkδkl\text{Cov}(w_k, w_l)=\eta_k'\Sigma\eta_l=\lambda_k\delta_{kl}，这里δkl\delta_{kl}是Kronecker delta function（δkk=1\delta_{kk}=1，对于k≠lk\neq l有δkl=0\delta_{kl}=0）。由于Λ\Lambda为按特征值降序排列的对角矩阵，因此必有

Var(w1)≥Var(w2)≥⋯≥Var(wk)\text{Var}(w_1)\ge \text{Var}(w_2)\ge \cdots\ge \text{Var}(w_k)

对于Σ\Sigma的特征值，有 ∑j=1dλj=∑j=1dσj2\sum_{j=1}^{d}\lambda_j=\sum_{j=1}^{d}\sigma^2_j 其中σj2\sigma_j^2为Σ\Sigma的第jj个对角线元素。这是因为Σ\Sigma与Λ\Lambda为相似矩阵，所以它们的迹一定相同。

如果用第2节中有关样本的设定，同样由于S=Γ^Λ^Γ^′S=\hat\Gamma \hat\Lambda \hat\Gamma'，SS与Λ^\hat\Lambda为相似矩阵，因此必有

∑j=1dλ^j=∑j=1dsj2\sum_{j=1}^{d} \hat\lambda_j = \sum_{j=1}^{d} s_j^2

4.2 最优化性质

我们从另一个角度出发，来看主成分分析。

对于第1节中设定的总体xx，假设Σ\Sigma满秩，我们能不能找到一个dd维的单位向量uu，使得u′xu'x的方差最大？

由于Σ\Sigma满秩，那么它的dd个特征向量可构成正交基，不妨设u=c1η1+⋯+cdηdu=c_1 \eta_1+\cdots+c_d \eta_d，而u′u=1u'u=1就等价于∑ici2=1\sum_i c_i^2=1，则方差

Var(u′x)=u′Σu=∑j,kcjckηj′Σηk=∑j,kcjckλkηj′ηk=∑jcj2λj≤λ1∑jcj2=λ1\begin{aligned} \text{Var}(u'x) =& u'\Sigma u\\ =& \sum_{j,k} c_j c_k \eta_j'\Sigma\eta_k\\ =& \sum_{j,k} c_j c_k \lambda_k \eta_j'\eta_k\\ =& \sum_j c_j^2\lambda_j\\ \leq & \lambda_1 \sum_j c_j^2\\ =& \lambda_1 \end{aligned}

假设各特征值互不相同，那么上式只有当c12=1c_1^2=1且cj=0c_j=0（j≠1j\neq 1）时等号成立，此时u=±η1u=\pm \eta_1。

接下来，我们寻找下一个u2u_2，要求u2′xu_2'x与η1′x\eta_1' x不相关，并能使u′xu'x的方差最大。同样可设u2=c1η1+⋯+cdηdu_2=c_1 \eta_1+\cdots+c_d \eta_d，u2′xu_2'x与η1′x\eta_1' x不相关等价于u2′η1=0u_2'\eta_1=0，此时又等价于c1=0c_1=0，重复上述过程，可以解出u2′xu_2'x的最大方差为λ2\lambda_2，此时u2=±η2u_2=\pm \eta_2。

不断重复上述过程，可以找出dd个uu。

如果rank(Σ)=r≤d\text{rank}(\Sigma)=r\leq d，那么可找出rr个特征值。如果特征值有重复，那么找出的uu可能有多个解。

avatar
文章
阅读
粉丝