PCA——主成分分析数学原理及代码

主成分分析

PCA的目的是：对数据进行一个线性变换，在最大程度保留原始信息的前提下去除数据中彼此相关的信息。反映在变量上就是说，对所有的变量进行一个线性变换，使得变换后得到的变量彼此之间不相关，并且是所有可能的线性变换中方差最大的一些变量（我们认为方差体现了信息量的大小）。

总体主成分分析

设$\mathbf{X}$是一个$n$维随机向量，其均值向量为$\mu ={\mu }_1,\dots ,{\mu }_n$、协方差矩阵为$\Sigma =({\sigma }_{ij}),i,j=1,2,\dots ,n$。对$\mathbf{X}$进行一个线性变换$\mathcal{T}$得到一个$n$维随机向量$\mathbf{Y}=({\mathbf{Y}}_1,\dots ,{\mathbf{Y}}_n)$，$\mathcal{T}$的矩阵为：
$$\begin{array}{c}A=\left(\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ ⋮\\ {\alpha }_n\end{array}\right)\end{array}$$
若：

1. $\mathrm{Cov}(\mathbf{Y})$是一个对角矩阵，即$\mathrm{Cov}({\mathbf{Y}}_i,{\mathbf{Y}}_j)=0,i\ne j$；
2. ${\mathbf{Y}}_1$是所有对$\mathbf{X}$进行线性变换后得到的随机变量中方差最大的随机变量，${\mathbf{Y}}_2$是与${\mathbf{Y}}_1$不相关的所有对$\mathbf{X}$进行线性变换后得到的随机变量中方差第二大的随机变量，以此类推。

则分别称${\mathbf{Y}}_1,{\mathbf{Y}}_2,\dots ,{\mathbf{Y}}_n$是第一、第二、……、第$n$主成分。
这一定义是否足够？

若不对$\mathcal{T}$的矩阵$A$作出相应的限制，对$\mathbf{X}$进行线性变换后得到的${\mathbf{Y}}_i,i=1,2,\dots ,n$的方差可以任意大。
$$\begin{array}{c}\mathrm{Var}({\mathbf{Y}}_i)=\mathrm{Cov}(\mathbf{Y}{)}_{(i,i)}=\mathrm{Cov}(A\mathbf{X}{)}_{(i,i)}=(A\Sigma A^T{)}_{i,i}={\alpha }_i\Sigma {\alpha }_i^T\end{array}$$
若$\mathrm{Var}({\mathbf{Y}}_i)>0$，取矩阵$B=kA$，$\mathbf{Z}=kA\mathbf{X}$，则：
$$\begin{array}{c}\mathrm{Var}({\mathbf{Y}}_i)=(k{\alpha }_i)\Sigma (k{\alpha }_i{)}^T=k^2{\alpha }_i\Sigma {\alpha }_i\end{array}$$
改变$k$的值，即可对${\mathbf{Y}}_i,i=1,2,\dots ,n$的方差进行任意的放缩。
因此，我们需要对$A$进行相应的限制，在这里我们人为地选择要求$A$是一个正交矩阵，也就是让${\alpha }_i{\alpha }_i^T=1$。

总体主成分完整定义

设$\mathbf{X}$是一个$n$维随机向量，其均值向量为$\mu ={\mu }_1,\dots ,{\mu }_n$、协方差矩阵为$\Sigma =({\sigma }_{ij}),i,j=1,2,\dots ,n$。对$\mathbf{X}$进行一个线性变换$\mathcal{T}$得到一个$n$维随机向量$\mathbf{Y}=({\mathbf{Y}}_1,\dots ,{\mathbf{Y}}_n)$，$\mathcal{T}$的矩阵为：
$$\begin{array}{c}A=\left(\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ ⋮\\ {\alpha }_n\end{array}\right)\end{array}$$
若：

1. $A{A}^T=I$；
2. $\mathrm{Cov}(\mathbf{Y})$是一个对角矩阵，即$\mathrm{Cov}({\mathbf{Y}}_i,{\mathbf{Y}}_j)=0,i\ne j$；
3. ${\mathbf{Y}}_1$是所有对$\mathbf{X}$进行线性变换后得到的随机变量中方差最大的随机变量，${\mathbf{Y}}_2$是与${\mathbf{Y}}_1$不相关的所有对$\mathbf{X}$进行线性变换后得到的随机变量中方差第二大的随机变量，以此类推。

则分别称${\mathbf{Y}}_1,{\mathbf{Y}}_2,\dots ,{\mathbf{Y}}_n$是第一、第二、……、第$n$主成分。

主成分求解定理

设$\mathbf{X}$是一个$n$维随机向量，$\Sigma$是其协方差矩阵，$\Sigma$的特征值\footnote{若特征多项式有重根，则标准正交化特征向量组不唯一，主成分也不唯一。}从大到小记作${\lambda }_1,\dots ,{\lambda }_n$，${\phi }_1,\dots ,{\phi }_n$为对应的标准正交化特征向量，则$\mathbf{X}$的第$i$个主成分以及其方差为：
$$\begin{array}{c}{\mathbf{Y}}_i={\phi }_i^T\mathbf{X},\mathrm{Var}({\mathbf{Y}}_i)={\phi }_i^T\Sigma {\phi }_i={\lambda }_i\end{array}$$
考虑到：
$$\begin{array}{c}\mathrm{Var}({\mathbf{Y}}_i)={\alpha }_i\Sigma {\alpha }_i^T,\mathrm{Cov}({\mathbf{Y}}_i,{\mathbf{Y}}_j)={\alpha }_i\Sigma {\alpha }_j^T\end{array}$$
求解主成分的过程即为求解：
$$\begin{array}{c}{\alpha }_i=\arg \max {\alpha }_i\Sigma {\alpha }_i^T\\ s.t.\{\begin{array}{ll}||{\alpha }_i||=1,& i=1,2,\dots ,n\\ {\alpha }_i\Sigma {\alpha }_j=0,& j<i\end{array}\end{array}$$

于是上述结论成立。

因子载荷的定义

将第$i$个主成分${\mathbf{Y}}_i$与变量${\mathbf{X}}_j$的相关系数$\rho ({\mathbf{Y}}_i,{\mathbf{X}}_j)$称为因子载荷。可推得：
$$\begin{array}{c}\rho ({\mathbf{Y}}_i,{\mathbf{X}}_j)=\frac{\sqrt{{\lambda }_i}{\alpha }_{ij}}{\sqrt{{\sigma }_{jj}}},i,j=1,2,\dots ,n\end{array}$$
由相关系数的定义：
$$\begin{array}{rl}\rho ({\mathbf{Y}}_i,{\mathbf{X}}_j)& =\frac{\mathrm{Cov}({\mathbf{Y}}_i,{\mathbf{X}}_j)}{\sqrt{\mathrm{Var}({\mathbf{Y}}_i)\mathrm{Var}({\mathbf{X}}_j)}}=\frac{\mathrm{Cov}({\alpha }_i\mathbf{X},e_j^T\mathbf{X})}{\sqrt{{\lambda }_i{\sigma }_{jj}}}\\ & =\frac{{\alpha }_i\Sigma e_j}{\sqrt{{\lambda }_i{\sigma }_{jj}}}=\frac{e_j^T\Sigma {\alpha }_i}{\sqrt{{\lambda }_i{\sigma }_{jj}}}=\frac{e_j^T{\lambda }_i{\alpha }_i}{\sqrt{{\lambda }_i{\sigma }_{jj}}}=\frac{\sqrt{{\lambda }_i}{\alpha }_{ij}}{\sqrt{{\sigma }_{jj}}}\end{array}$$

总体主成分的性质

总体主成分具有如下性质：

1. Cov ⁡ ( Y ) = diag ⁡ { λ 1 , … , λ n } \operatorname{Cov}(\mathbf{Y})=\operatorname{diag}\{\lambda_1,\dots,\lambda_n\} Cov(Y)=diag{λ1​,…,λn​}；
2. $\mathbf{Y}$的方差之和等于$\mathbf{X}$的方差之和，即$\sum _{i=1}^n{\lambda }_i=\sum _{i=1}^n{\sigma }_{ii}$；
3. 第$i$个主成分与原变量的因子负荷量满足：
   $$\begin{array}{c}\sum _{j=1}^n{\sigma }_{jj}{\rho }^2({\mathbf{Y}}_i,{\mathbf{X}}_j)={\lambda }_i\end{array}$$
4. 原变量的第$j$个分量与所有主成分的因子负荷量满足：
   $$\begin{array}{c}\sum _{i=1}^n{\rho }^2({\mathbf{Y}}_i,{\mathbf{X}}_j)=1\end{array}$$
   证明：
   1. 由PCA求解过程直接可得。
   2. 显然：
      $$\begin{array}{rl}\sum _{i=1}^n\mathrm{Var}({\mathbf{Y}}_i)& =\mathrm{tr}[\mathrm{Cov}(\mathbf{Y})]=\mathrm{tr}[\mathrm{Cov}(A\mathbf{X})]=\mathrm{tr}(A\Sigma A^T)\\ & =\mathrm{tr}(\Sigma A^{T}A)=\mathrm{tr}(\Sigma )=\sum _{i=1}^n\mathrm{Var}({\mathbf{X}}_i)\end{array}$$
   3. 显然：
      $$\begin{array}{c}\sum _{j=1}^n{\sigma }_{jj}{\rho }^2({\mathbf{Y}}_i,{\mathbf{X}}_j)=\sum _{j=1}^n{\lambda }_i{\alpha }_{ij}^2={\lambda }_i{\alpha }_i{\alpha }_i^T={\lambda }_i\end{array}$$
   4. 因为$A$是正交矩阵，所以$A$可逆，于是$\mathbf{X}$可以表示为${\mathbf{Y}}_1,\dots ,\mathbf{n}$的线性组合，所以二者的复相关系数为$1$。由复相关系数性质可直接得出结论。

贡献率定义

称第$i$个主成分${\mathbf{Y}}_i$的方差与所有主成分方差之和为${\mathbf{Y}}_i$的方差贡献率，记为${\eta }_i$，即：
$$\begin{array}{c}{\eta }_i=\frac{{\lambda }_i}{\sum _{j=1}^n{\lambda }_j}\end{array}$$
将：
$$\begin{array}{c}\frac{\sum _{i=1}^k{\lambda }_i}{\sum _{i=1}^n{\lambda }_i}\end{array}$$
称为主成分${\mathbf{Y}}_1,\dots .{\mathbf{Y}}_k$的累计方差贡献率。称主成分${\mathbf{Y}}_1,\dots .{\mathbf{Y}}_k$与变量${\mathbf{X}}_j$之间的复相关系数的平方$R^2$为${\mathbf{Y}}_1,\dots .{\mathbf{Y}}_k$对${\mathbf{X}}_j$的贡献率，其计算公式为：
$$\begin{array}{c}{R}^2=\sum _{i=1}^k\frac{{\lambda }_i{\alpha }_{ij}^2}{{\sigma }_{ii}}\end{array}$$

由前述，我们一般通过选择主成分的个数来实现对数据的降维，即选择主成分的个数使它们的累计方差贡献率达到一定比例（一般为$85\%$）。

样本主成分分析

假设对$n$维随机变量$\mathbf{X}$进行$m$次独立观测，得到$m$个$n$维样本$x_1,\dots ,x_m$。在样本主成分分析中，我们使用样本来估计$\mathbf{X}$的协方差矩阵，即：
$$\begin{array}{c}S=(s_{ij}),s_{ij}=\frac{1}{m-1}\sum _{k=1}^m(x_{ki}-{\hat{\mathbf{X}}}_i)(x_{kj}-{\hat{\mathbf{X}}}_j),i,j=1,2,\dots ,n\end{array}$$
其中：
$$\begin{array}{c}{\hat{\mathbf{X}}}_i=\frac{1}{m}\sum _{j=1}^m{x}_{ji},i=1,2,\dots ,n\end{array}$$
其余步骤与总体主成分分析一致。

注意事项

多重共线性问题

当原始变量出现多重共线性时，PCA的效果会受到影响，这是因为重复的信息在方差占比中重复进行了计算。我们可以通过计算协方差矩阵的最小特征值来判断是否出现多重共线性的情况。若最小特征值趋于$0$，则需要对纳入研究的变量进行考察与筛选。

相关矩阵导出主成分

上面我们都是对协方差矩阵的特征值分解进行计算，但在现实中，我们可能会对数据进行标准化处理来消除量纲带来的影响，注意到标准化后数据的协方差矩阵即为相关矩阵，此时将相关矩阵作对应的特征值分解即可。但需要注意：标准化后各变量方差相等均为$1$，损失了部分信息，所以会使得标准化后的各变量在对主成分构成中的作用趋于相等。因此，取值范围在同量级的数据建议使用协方差矩阵直接求解主成分，若变量之间数量级差异较大，再使用相关矩阵求解主成分。

算法流程

主成分分析的应用

主成分分类

可以实现对变量之间的分类。将变量进行主成分分析，得到第一主成分与第二主成分，然后画出各变量与两个主成分载荷的二维平面图，即将各变量画在以两个主成分为轴的平面上，变量的两个坐标是主成分在变量上的载荷。可以认为相似变量会聚在平面图中聚在一起。

主成分回归

若数据存在高度的多重共线性，对数据进行主成分筛选，然后用主成分去进行回归。

R语言代码

> # 以鸢尾花数据集作为示例
> data <- iris[1:4]
> # 变量间量纲差异不大，使用协方差矩阵进行分析
> Sigma <- cov(data)
> # 计算特征值，因为是示例就不做处理了
> eigen(Sigma)$values
[1] 4.22824171 0.24267075 0.07820950 0.02383509
> # center和scale.控制是否对原始数据进行标准化
> # 注意scale.，后面有一个点，不是scale
> x <- prcomp(data, center=FALSE, scale.=FALSE)
> # 呈现结果
> summary(x)
Importance of components:PC1     PC2     PC3     PC4
Standard deviation     7.8613 1.45504 0.28353 0.15441
Proportion of Variance 0.9653 0.03307 0.00126 0.00037
Cumulative Proportion  0.9653 0.99837 0.99963 1.00000
> # 输出因子载荷矩阵，显示出的矩阵其实就是A^T
> # 可以自己尝试证明一下，利用协方差的性质以及Y之间的线性无关性
> x$rotationPC1        PC2         PC3        PC4
Sepal.Length -0.7511082  0.2841749  0.50215472  0.3208143
Sepal.Width  -0.3800862  0.5467445 -0.67524332 -0.3172561
Petal.Length -0.5130089 -0.7086646 -0.05916621 -0.4807451
Petal.Width  -0.1679075 -0.3436708 -0.53701625  0.7518717
> # 如果输入x，则输出训练样本经过主成分分析后的结果，即Y值
> # 如果输入一组新的数据newdata，则输出newdata经过主成分分析后的结果
> # 这里因为数据太多就用head函数控制只显示前10行
> head(predict(x), 10)PC1      PC2          PC3          PC4[1,] -5.912747 2.302033  0.007401536  0.003087706[2,] -5.572482 1.971826  0.244592251  0.097552888[3,] -5.446977 2.095206  0.015029262  0.018013331[4,] -5.436459 1.870382  0.020504880 -0.078491501[5,] -5.875645 2.328290 -0.110338269 -0.060719326[6,] -6.477598 2.324650 -0.237202487 -0.021419633[7,] -5.515975 2.070904 -0.229853120 -0.050406649[8,] -5.850929 2.148075  0.018793774 -0.045342619[9,] -5.158920 1.775064  0.061039220 -0.031128633
[10,] -5.645001 1.990001  0.224852923 -0.057434390
> # 画碎石图，即样本主成分与方差贡献率之间的曲线图或条形图
> screeplot(x, type = "barplot")
> screeplot(x, type = "lines")

主成分分析在R语言中还有另一个函数叫做princomp，但官方更推荐使用prcomp。主要是因为求解算法问题，prcomp使用svd分解进行求解，后面我会对这里进行一些补充。

Python代码

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris# 加载数据
iris = load_iris()
X = iris.data
y = iris.target
feature_names = iris.feature_names# 对数据进行标准化（可选）
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)# 初始化 PCA，保留所有主成分
pca = PCA(n_components=4)
X_pca = pca.fit_transform(X_scaled)print("各主成分解释的方差比例：")
print(pca.explained_variance_ratio_)
print("累计解释方差：", np.cumsum(pca.explained_variance_ratio_))# 每个主成分上的特征权重（特征向量）
loadings = pd.DataFrame(pca.components_.T,columns=[f'PC{i+1}' for i in range(pca.n_components_)],index=feature_names)print(loadings)

参考文献

1. 薛毅，统计建模与R软件
2. 何晓群，多元统计分析——基于R语言
3. 李航，统计学方法
4. Richard A.Johnson，Applied Multivariate Statistical Analysis

后续会补充的内容

1. 具体数值求解算法
2. R语言两个函数的具体差异
3. 为什么特征值接近于0就意味着多重共线性
4. R2的性质，即为什么R2等于所有自变量与因变量相关系数平方的加和