当前位置：首页 > news >正文

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

news 来源：原创 2025/8/4 3:01:43

主成分分析

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{equation*} A= \begin{pmatrix} \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_n \end{pmatrix} \end{equation*}$
若：

$\operatorname{Cov}(\mathbf{Y})$ 是一个对角矩阵，即 $\operatorname{Cov}(\mathbf{Y}_i,\mathbf{Y}_j)=0,\;i\ne j$ ；
$\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{equation*} \operatorname{Var}(\mathbf{Y}_i)=\operatorname{Cov}(\mathbf{Y})_{(i,i)}=\operatorname{Cov}(A\mathbf{X})_{(i,i)}=(A\Sigma A^T)_{i,i}=\alpha_i\Sigma\alpha_i^T \end{equation*}$
若 $\operatorname{Var}(\mathbf{Y}_i)>0$ ，取矩阵 $B = k A$ ， $\mathbf{Z}=kA\mathbf{X}$ ，则：
$\begin{equation*} \operatorname{Var}(\mathbf{Y}_i)=(k\alpha_i)\Sigma(k\alpha_i)^T=k^2\alpha_i\Sigma\alpha_i \end{equation*}$
改变 $k$ 的值，即可对 $\mathbf{Y}_i,\;i=1,2,\dots,n$ 的方差进行任意的放缩。
因此，我们需要对 $A$ 进行相应的限制，在这里我们人为地选择要求 $A$ 是一个正交矩阵，也就是让 $\alpha_i\alpha_i^T=1$ 。

总体主成分完整定义

$AA^T=I$ ；
$\operatorname{Cov}(\mathbf{Y})$ 是一个对角矩阵，即 $\operatorname{Cov}(\mathbf{Y}_i,\mathbf{Y}_j)=0,\;i\ne j$ ；
$\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$ ， $\varphi_1,\dots,\varphi_n$ 为对应的标准正交化特征向量，则 $\mathbf{X}$ 的第 $i$ 个主成分以及其方差为：
$\begin{equation*} \mathbf{Y}_i=\varphi_i^T\mathbf{X},\;\operatorname{Var}(\mathbf{Y}_i)=\varphi_i^T\Sigma\varphi_i=\lambda_i \end{equation*}$
考虑到：
$\begin{equation*} \operatorname{Var}(\mathbf{Y}_i)=\alpha_i\Sigma\alpha_i^T,\quad \operatorname{Cov}(\mathbf{Y}_i,\mathbf{Y}_j)=\alpha_i\Sigma\alpha_j^T \end{equation*}$
求解主成分的过程即为求解：
$\begin{gather*} \alpha_i=\arg\max\alpha_i\Sigma\alpha_i^T \\ \operatorname{s.t.} \begin{cases} ||\alpha_i||=1,\;&i=1,2,\dots,n\\ \alpha_i\Sigma\alpha_j=0,\;&j<i \end{cases} \end{gather*}$
在这里插入图片描述
于是上述结论成立。

因子载荷的定义

将第 $i$ 个主成分 $\mathbf{Y}_i$ 与变量 $\mathbf{X}_j$ 的相关系数 $\rho(\mathbf{Y}_i,\mathbf{X}_j)$ 称为因子载荷。可推得：
$\begin{equation*} \rho(\mathbf{Y}_i,\mathbf{X}_j)=\frac{\sqrt{\lambda_i}\alpha_{ij}}{\sqrt{\sigma_{jj}}},\;i,j=1,2,\dots,n \end{equation*}$
由相关系数的定义：
$\begin{align} \rho(\mathbf{Y}_i,\mathbf{X}_j) &=\frac{\operatorname{Cov}(\mathbf{Y}_i,\mathbf{X}_j)}{\sqrt{\operatorname{Var}(\mathbf{Y}_i)\operatorname{Var}(\mathbf{X}_j)}}=\frac{\operatorname{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{align}$

总体主成分的性质

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

$\operatorname{Cov}(\mathbf{Y})=\operatorname{diag}\{\lambda_1,\dots,\lambda_n\}$ ；
$\mathbf{Y}$ 的方差之和等于 $\mathbf{X}$ 的方差之和，即 $\sum\limits_{i=1}^{n}\lambda_i=\sum\limits_{i=1}^{n}\sigma_{ii}$ ；
第 $i$ 个主成分与原变量的因子负荷量满足：
$\begin{equation*} \sum_{j=1}^{n}\sigma_{jj}\rho^2(\mathbf{Y}_i,\mathbf{X}_j)=\lambda_i \end{equation*}$
原变量的第 $j$ 个分量与所有主成分的因子负荷量满足：
$\begin{equation*} \sum_{i=1}^{n}\rho^2(\mathbf{Y}_i,\mathbf{X}_j)=1 \end{equation*}$
证明：
1. 由PCA求解过程直接可得。
2. 显然：
  $\begin{align*} \sum_{i=1}^{n}\operatorname{Var}(\mathbf{Y}_i) &=\operatorname{tr}[\operatorname{Cov}(\mathbf{Y})]=\operatorname{tr}[\operatorname{Cov}(A\mathbf{X})]=\operatorname{tr}(A\Sigma A^T) \\ &=\operatorname{tr}(\Sigma A^TA)=\operatorname{tr}(\Sigma)=\sum_{i=1}^{n}\operatorname{Var}(\mathbf{X}_i) \end{align*}$
3. 显然：
  $\begin{equation*} \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{equation*}$
4. 因为 $A$ 是正交矩阵，所以 $A$ 可逆，于是 $\mathbf{X}$ 可以表示为 $\mathbf{Y}_1,\dots,\mathbf{n}$ 的线性组合，所以二者的复相关系数为 $1$ 。由复相关系数性质可直接得出结论。

贡献率定义

称第 $i$ 个主成分 $\mathbf{Y}_i$ 的方差与所有主成分方差之和为 $\mathbf{Y}_i$ 的方差贡献率，记为 $\eta_i$ ，即：
$\begin{equation*} \eta_i=\frac{\lambda_i}{\sum\limits_{j=1}^{n}\lambda_j} \end{equation*}$
将：
$\begin{equation*} \frac{\sum\limits_{i=1}^{k}\lambda_i}{\sum\limits_{i=1}^{n}\lambda_i} \end{equation*}$
称为主成分 $\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{equation*} R^2=\sum_{i=1}^{k}\frac{\lambda_i\alpha_{ij}^2}{\sigma_{ii}} \end{equation*}$

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

样本主成分分析

假设对 $n$ 维随机变量 $\mathbf{X}$ 进行 $m$ 次独立观测，得到 $m$ 个 $n$ 维样本 $x_1,\dots,x_m$ 。在样本主成分分析中，我们使用样本来估计 $\mathbf{X}$ 的协方差矩阵，即：
$\begin{equation*} 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{equation*}$
其中：
$\begin{equation*} \hat{\mathbf{X}}_i=\frac{1}{m}\sum_{j=1}^{m}x_{ji},\;i=1,2,\dots,n \end{equation*}$
其余步骤与总体主成分分析一致。

注意事项

多重共线性问题

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

算法流程

在这里插入图片描述

主成分分析的应用

主成分分类

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

主成分回归

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

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")