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

1 总体的主成分

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

  • \(k\)principal component score主成分得分)为\(w_k=\eta_k' (x-\mu)\)
  • \(k\)principal component vector主成分向量)为\(w^{(k)}=(w_1,\ldots,w_k)'=\Gamma_k' (x-\mu)\)
  • \(k\)principal component projection (vector)主成分映射(向量))为\(p_k=\eta_k\eta_k'(x-\mu)=w_k\eta_k\)

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

2 样本的主成分

假设我们有\(n\)个样本,将它们排成\(d\times n\)的矩阵\(X=(x_1,\ldots,x_n)\),记样本均值\(\bar x=\dfrac{1}{n}X\ell_n\),样本协方差矩阵\(S=\dfrac{1}{n-1}(X-\bar x \ell_n')(X-\bar x \ell_n')'\),并且\(\text{rank}(S)=r\leq d\)

我们可以对样本协方差矩阵进行谱分解\(S=\hat\Gamma \hat\Lambda \hat\Gamma'\)。与在总体中的定义一样,可以在样本中定义如下概念:

  • \(k\)principal component score\(w_k=\hat\eta_k' (X-\bar x\ell_n')\),这是一个\(1\times n\)的行向量;
  • principal component data\(W^{(k)}=(w_1',\ldots,w_k')'=\hat\Gamma_k' (X-\bar x\ell_n')\),它是\(k\times n\)的矩阵;
  • \(k\)principal component projection\(P_k=\hat\eta_k\hat\eta_k'(X-\bar x\ell_n')=\hat\eta_k w_k\),它是\(d\times n\)的矩阵。

3 主成分的可视化

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

3.1 Scree plot

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

\[\dfrac{\lambda_k}{\sum_{j=1}^{r}\lambda_j}=\dfrac{\lambda_k}{\text{tr}(\Sigma)} \]

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

\[\dfrac{\sum_{k=1}^{\kappa}\lambda_k}{\sum_{j=1}^{r}\lambda_j}=\dfrac{\sum_{k=1}^{\kappa}\lambda_k}{\text{tr}(\Sigma)} \]

对于不同\(k\)画出前\(k\)个PC的贡献和累积贡献,就得到了scree plot,如下图:

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

3.2 PC score plot

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

也可以直接画出三维图,如下图:

3.3 Projection plot

接下来看如何将\(P_k\)可视化。\(P_k\)\(d\times n\)的矩阵,它的第\(i\)列就是第\(i\)个样本在\(\hat\eta_k\)方向的贡献。对于每个\(k\),都可以用\(P_k\)画出平行坐标图,这样就很方便地将\(P_k\)可视化了。

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

下图是一个\(d=30\)的示例,左图为某个PC的projection plot,右图为非参数密度估计(可用Matlab工具curvdatSM):

4 PC的性质

4.1 结构性质

本节考虑PC的一些性质,先看\(X\)和它的PC的相关结构。这里采用第1节的设定。

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

\[\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互不相关。当然,如果\(x\)是正态分布,那么PC scores也相互独立,否则需要用其他方法做出相互独立。

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

\[\text{Var}(w_1)\ge \text{Var}(w_2)\ge \cdots\ge \text{Var}(w_k) \]

对于\(\Sigma\)的特征值,有

\[\sum_{j=1}^{d}\lambda_j=\sum_{j=1}^{d}\sigma^2_j \]

其中\(\sigma_j^2\)\(\Sigma\)的第\(j\)个对角线元素。这是因为\(\Sigma\)\(\Lambda\)为相似矩阵,所以它们的迹一定相同。

如果用第2节中有关样本的设定,同样由于\(S=\hat\Gamma \hat\Lambda \hat\Gamma'\)\(S\)\(\hat\Lambda\)为相似矩阵,因此必有

\[\sum_{j=1}^{d} \hat\lambda_j = \sum_{j=1}^{d} s_j^2 \]

4.2 最优化性质

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

对于第1节中设定的总体\(x\),假设\(\Sigma\)满秩,我们能不能找到一个\(d\)维的单位向量\(u\),使得\(u'x\)的方差最大?

由于\(\Sigma\)满秩,那么它的\(d\)个特征向量可构成正交基,不妨设\(u=c_1 \eta_1+\cdots+c_d \eta_d\),而\(u'u=1\)就等价于\(\sum_i c_i^2=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} \]

假设各特征值互不相同,那么上式只有当\(c_1^2=1\)\(c_j=0\)\(j\neq 1\))时等号成立,此时\(u=\pm \eta_1\)

接下来,我们寻找下一个\(u_2\),要求\(u_2'x\)\(\eta_1' x\)不相关,并能使\(u'x\)的方差最大。同样可设\(u_2=c_1 \eta_1+\cdots+c_d \eta_d\)\(u_2'x\)\(\eta_1' x\)不相关等价于\(u_2'\eta_1=0\),此时又等价于\(c_1=0\),重复上述过程,可以解出\(u_2'x\)的最大方差为\(\lambda_2\),此时\(u_2=\pm \eta_2\)

不断重复上述过程,可以找出\(d\)\(u\)

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