主成分分析 (Principal components analysis，PCA)

主成分分析（PCA）是一种线性降维算法，它将数据的变化由高维压缩到更少的维度（主成分）。

设想存在这样一个情景：存在p个不同的菌，每个菌有n个基因，每个基因都有一个表达量，现在想要知道：

1、根据基因的表达量，哪些菌是更相近的？

2、不同的菌是否呈现轨迹变化？

这对应了PCA两种常用的功能：数据的质量控制和结构探索。

简单了解：

在绪论的情景中，我们当然可以画一个三维图片，x轴和y轴代表不同的菌和不同的基因，z轴代表这些基因的表达量。但是它的可读性、直观性和美观性都会非常差。所以我们使用PCA去处理这种问题。

在最常见的场景，一个PCA的结果往往被可视化为一个散点图。其横坐标和纵坐标分别一般被命名为PC1和PC2。

以绪论中的情景为例，在p×n个变量中，PC1的含义是数据中最大变异模式对应的线性方向，代表最主要的结构性差异来源；PC2代表在排除 PC1 后剩余变化中最显著的正交模式，代表第二层结构性差异。

如果还是难以理解，我们可以将绪论的问题简化为p个菌中2个基因的表达量。这样对于每个菌都有两个基因和其表达量，我们可以将每个菌两个表达量变为两个坐标，即x坐标和y坐标，这样每个菌就变为了一个二维图像内的点(x,y)。PCA的核心就是**在这个二维平面中寻找一条线，使得所有点投影到该方向上的方差最大；第二条线是在与第一条线垂直的前提下，使投影方差最大的方向。**这两条方向即分别对应 PC1 和 PC2。

简言之，PC1代表最主要的差异，PC2代表次主要的差异。在绪论的例子中，可视化的散点图共有p个点，每个点PC1 或 PC2 的数值，本质上是： **这个样本（菌）在某种变化趋势上的得分。**数值越大，说明越符合该趋势；数值越小（或为负），说明越偏向相反趋势。

PC1和PC2后往往带有一个百分比数字，如PC1（45%）、PC2（30%），这代表可解释性，即差异最大的方向解释了45%的数据，次大的方向解释了30%的数据，那么理论上PC1+PC2的数据（即从散点图上的数据）理论上可以重现45%+35%=80%的原始数据。

深度了解：

下列内容涉及高等代数知识，但实操方法部分仅需要本科线性代数知识

以绪论的内容为例，给定任意一个p个样本、n个观测的数据对象，可以构成如下样本观测矩阵x

$$ \mathbf{x} =\begin{bmatrix} x_{11}&x_{12} &\dots &x_{1n} \\ x_{21} & x_{22}& \dots&x_{2n} \\ \vdots &\vdots & \ddots &\vdots \\ x_{p1} &x_{p2} & \dots &x_{pn}\end{bmatrix} $$

1、基本样本统计量

先让我们了解一些基本量

定义1 样本均值（μ）、样本方差（σ2）：给定某个样本x的n个观测量

$$ x_{1},x_{2},\dots ,x_{n} $$

分别称

$$ \mu = mean(\mathbf{x} ) = \frac{1}{n} \sum_{i=1}^{n}x_{i}= \frac{1}{n}\mathbf{1} ^{\top} \mathbf{x}\\ \sigma ^{2}=var(x)=\frac{1}{n}\sum_{i=1}^{n} (x_{i}-\mu)^{2}=\frac{1}{n}(\mathbf{x} -\mu \mathbf{1} )^{\top }(\mathbf{x} -\mu \mathbf{1} )\\ $$

为随机变量的样本均值（简称均值）和样本方差（简称方差），其中

$$ \mathbf{x} =\left [ x_{1} \quad x_{2}\quad\cdots \quad x_{n}\right ] ^{\top},\mathbf{1} =\left [ 1 \quad 1\quad\cdots \quad 1\right ] ^{\top} $$

定义2 样本协方差（cov）：给定某个样本x的n个观测量以及另一个样本y的n个观测量

$$ x_{1},x_{2},\cdots ,x_{n} \\y_{1},y_{2},\cdots ,y_{n} $$

称

$$ cov(\mathbf{x},\mathbf{y})=\frac{1}{n} \sum_{i=1}^{n} (x_{i}-\mu_{x})(y_{i}-\mu_{y})=\frac{1}{n}(\mathbf{x} - \mu_{x}\mathbf{1})^{\top}(\mathbf{y} - \mu_{y}\mathbf{1}) $$

为样本x和y的样本协方差（简称协方差），其中x和1定义如上，y定义与x相同，这里不做赘述。

显而易见的，对p个样本，应该有p×p个协方差。

定义3 样本均值向量（μ）和样本协方差矩阵（Σ） ：给定p维样本向量X的n个观测值（即p×n个样本），称

$$ \boldsymbol{\mu }=mean(\mathbf{x })=\begin{bmatrix} mean(\mathbf{x_{1} })\\mean(\mathbf{x_{2} }) \\ \vdots\\mean(\mathbf{x_{p} })\end{bmatrix}=\begin{bmatrix} \mu_{1}\\ \mu_{2}\\\vdots \\\mu_{p}\end{bmatrix}=\frac{1}{n}\mathbf{x1} $$

为样本向量X的样本矩阵x的样本均值向量（简称均值向量或均值)。其中

$$ \mu _{i}=mean(\mathbf{\mathit{X_{i}} } )=\frac{1}{n} \sum_{k=1}^{n} x_{ik}=\frac{1}{n}\mathbf{1}^{\top} \mathbf{x} _{i},i=1,2,\cdots ,p $$

为样本Xi的均值。

根据上面协方差的定义，对p个样本，计两两之间的协方差为

$$ \sigma_{ij}= cov(\mathbf{\mathit{X} }_{i}, \mathbf{\mathit{X} }_{j} )=\frac{1}{k} \sum_{i=1}^{k} (\mathbf{\mathit{X} }_{i}-\mu_{i})(\mathbf{\mathit{X} }_{j}-\mu_{j}),i,j=1,2,\cdots ,p $$

则

$$ \Sigma =cov(\mathbf{x} )=\begin{bmatrix} \sigma _{11} &\sigma _{12} &\dots &\sigma _{1p}\\\sigma _{21} &\sigma _{22} &\dots &\sigma _{2p}\\ \vdots &\vdots &\ddots &\vdots \\\sigma _{p1} &\sigma _{p2} &\dots &\sigma _{pp}\end{bmatrix}\\=\begin{bmatrix} cov(X_{1},X_{1}) &cov(X_{1},X_{2}) &\dots &cov(X_{1},X_{p})\\cov(X_{2},X_{1}) &cov(X_{2},X_{2}) &\dots &cov(X_{2},X_{p})\\ \vdots &\vdots &\ddots &\vdots \\cov(X_{p},X_{1}) &cov(X_{p},X_{2}) &\dots &cov(X_{p},X_{p})\end{bmatrix}\\=\frac{1}{n} (\mathbf{x}-\boldsymbol{\mu1}^{\top } )(\mathbf{x}-\boldsymbol{\mu1}^{\top } )^{\top}=\frac{1}{n} (\mathbf{xx}^{\top}-\boldsymbol{\mu\mu}^{\top}) $$

为样本向量X的样本观测矩阵x的样本协方差矩阵（简称x的协方差矩阵）

2、定义任意的方向方差

由前面的内容可以得到，给定的p个样本、n个观测的样本观测矩阵x，PCA的核心目标是寻找一个能够最大程度保留数据方差信息的投影方向。

为了简化运算，我们通过中心化处理将坐标原点平移至样本均值，使数据的均值向量为零。则对于任意方向向量u，观测数据在任意方向u上的投影方差为

$$ var(\mathbf{u}^{\top} \mathbf{x})=\frac{1}{n} \mathbf{u}^{\top} \mathbf{x}(\mathbf{u}^{\top} \mathbf{x})=\mathbf{u}^{\top}(\frac{1}{n}\mathbf{x}\mathbf{x}^{\top} )\mathbf{u}=\mathbf{u}^{\top}\Sigma\mathbf{u} $$

这个公式是如此的优美。它表明，数据在任意方向的方差都可以通过协方差矩阵Σ以及方向向量u以双线性的形式表达。

3、求解方向向量

对于刚才的公式可以发现，不对向量u的长度加以限制的话，var将随着u长度的增加而增加，并且上限趋于无穷大。因此，为了消除向量u的长度的影响，需要将其限制为单位向量，这便有了如下主成分分析的优化模型

$$ \left\{\begin{matrix} \underset{u}{\max}\mathbf{u}^{\top}\mathbf{\Sigma u} \\ \mathrm{s.t.} \ \mathbf{u}^{\top}\mathbf{u}=1 \end{matrix}\right. $$

显然，该模型为等式约束下的目标函数的极值问题，可以用拉格朗日乘数法构建如下目标函数

$$ \mathcal{L}(\mathbf{u},\lambda ) =\frac{1}{2} \mathbf{u}^{\top}\mathbf{\Sigma u}+\frac{\lambda }{2}(1-\mathbf{u}^{\top}\mathbf{u}) $$

两侧对u求导，得到

$$ \frac{\partial \mathcal{L}(\mathbf{u},\lambda ) }{\partial \mathbf{u}} =\mathbf{\Sigma u} -\lambda \mathbf{u} $$

由于L(u,λ)在极值点处的偏导数必定为零向量，因此可得

$$ \mathbf{\Sigma u} =\lambda \mathbf{u} $$

神奇的事情出现了，上述公式突然变成了协方差矩阵Σ的特征值与特征向量问题，我们可以很简单的求得u和λ。

4、进行主成分分析

在得到数据协方差矩阵Σ的所有特征值λ和特征向量u之后，即可对原始数据进行主成分分析。

$$ \mathbf{Y}= \mathbf{U}^{\top}\mathbf{x}\\ \mathbf{U}=\begin{bmatrix} \mathbf{u}_{1} & \mathbf{u}_{2} &\cdots &\mathbf{u}_{p}\end{bmatrix}; \mathbf{Y}=\begin{bmatrix} \mathbf{Y}_{1}\\\mathbf{Y}_{2} \\ \cdots\\\mathbf{Y}_{p}\end{bmatrix} $$

其中U为特征向量u构成的特征向量矩阵,Yi (i=1,2,…,p)称作x的第i个主成分。

假设协方差矩阵的p个特征值满足λ1 ≥ λ2 ≥ … ≥ λp，基于协方差矩阵Σ与特征值λ和特征向量u的关系可得

$$ \mathbf{u}_{i}^{\top} \boldsymbol{\Sigma}\mathbf{u}_{i}=\lambda _{i}\mathbf{u}_{i}^{\top}\mathbf{u}_{i}=\lambda _{i} $$

即协方差矩阵的第i个特征值λi正是数据在第i个特征向量方向上的方差。因此，可以定义数据前k个主成分的贡献率为

$$ \eta (k)=\frac{\lambda _{1}+\lambda _{2}+\cdots +\lambda _{k}}{\lambda _{1}+\lambda _{2}+\cdots +\lambda _{p}} $$

在应用中，可以根据给定的贡献率来确定需要保留的主成分的个数。

实操

接下来，我们通过一个简单的例子，来熟悉主成分分析的主要流程。

对以下一组观测数据进行PCA分析：

$$ \mathbf{x} =\begin{bmatrix} 1 & 2 &3 \\ 1& 1&3 \end{bmatrix} $$

显然，该数据可以视为二维平面上的三个点（即 p = 2、n = 3），坐标为（1 , 1）、（2 , 1）、（3 , 3）

步骤 1 观测矩阵中心化：

首先来计算均值向量μ：

$$ \boldsymbol{\mu}= \frac{1}{3} \mathbf{X1} =\frac{1}{3}\begin{bmatrix} 1 & 2&3 \\ 1& 1 &3 \end{bmatrix}\begin{bmatrix} 1\\ 1\\ 1 \end{bmatrix}=\begin{bmatrix} 2\\ \frac{5}{3} \end{bmatrix} $$

对数据中心化：

$$ \mathbf{X} =\mathbf{X}-\boldsymbol{\mu1}^{\top }\begin{bmatrix} -1 & 0 & 1\\ -\frac{2}{3} & -\frac{2}{3} & \frac{4}{3} \end{bmatrix} $$

步骤2 计算协方差矩阵

计算协方差矩阵为：

$$ \boldsymbol{\Sigma}=\frac{1}{3}\mathbf{XX}^{\top } =\frac{1}{9}\begin{bmatrix} 6 & 6\\ 6&8 \end{bmatrix} $$

步骤 3 计算特征矩阵和特征值

对协方差矩阵进行特征分解

$$ \boldsymbol{\Sigma}=\mathbf{U} \boldsymbol{\Lambda }\mathbf{U}^{\top }\\\boldsymbol{\Lambda }=\begin{bmatrix} \lambda _{1}&0 \\ 0 &\lambda _{2} \end{bmatrix}=\begin{bmatrix} 1.4536 & 0\\ 0&0.1019 \end{bmatrix},\mathbf{U}=\begin{bmatrix} 0.6464 & -0.7630\\ 0.7630&0.6464 \end{bmatrix} $$

步骤 4 主成分变换

使用特征向量对数据进行主成分变换：

$$ \mathbf{Y}= \mathbf{U}^{\top }\mathbf{x}=\begin{bmatrix} -1.1551 &-0.5087 & 1.6637\\ 0.3321 & -0.4309 &0.0988 \end{bmatrix} $$

其中矩阵的第一行为Y1，及第一主成分；第二行为Y2，及第二主成分。

由步骤 3可得，两个特征值分别为λ1 = 1.4536，λ2 = 0.1019，故第一主成分的贡献率为：

$$ \eta _{1}=\frac{\lambda _{1}}{\lambda _{1}+\lambda _{2}} =\frac{1.4536}{1.4536+0.1019} =0.9345 $$

即第一主成分可以解释93.45%的内容，数据主要分布在第一主成分，只需要一个维度就可以描述这组数据，即 (-1.1551 , -0.5087 , 1.6637)。维度由二维降到了一维。

代码运行

下面给出一段python代码用于计算一组数据的PCA

import numpy as np

def main():
    # 1. 输入数据
    p = int(input("输入变量个数 p: "))
    n = int(input("输入观测个数 n: "))
    print(f"\n输入 {p} 行数据，每行 {n} 个数字，用空格分隔：")    
    X = []
    for i in range(p):
        row = list(map(float, input(f"第 {i+1} 行: ").split()))
        if len(row) != n:
            raise ValueError("输入数据长度与 n 不一致")
        X.append(row)   
    X = np.array(X)    
    print("\n观测矩阵 X:")
    print(X)  
    
    # 2. 计算样本均值向量 μ
    mu = np.mean(X, axis=1, keepdims=True)   
    print("\n样本均值向量 μ:")
    print(mu)
    
    # 3. 中心化
    X_centered = X - mu    
    print("\n中心化后的矩阵:")
    print(X_centered)
    
    # 4. 计算协方差矩阵
    Sigma = (1/n) * X_centered @ X_centered.T   
    print("\n协方差矩阵 Σ:")
    print(Sigma)
    
    # 5. 求解特征值和特征向量
    eigenvalues, eigenvectors = np.linalg.eig(Sigma)
    print("\n特征值 λ:")
    print(eigenvalues)   
    print("\n特征向量 U (每一列是一个特征向量):")
    print(eigenvectors)
    
    # 6. 按特征值从大到小排序
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]
    print("\n排序后的特征值:")
    print(eigenvalues)
    print("\n排序后的特征向量矩阵 U:")
    print(eigenvectors)
    
    # 7. 主成分
    Y = eigenvectors.T @ X_centered
    print("\n主成分矩阵 Y:")
    print(Y)
    
    # 8. 贡献率
    total_variance = np.sum(eigenvalues)
    contribution_ratio = eigenvalues / total_variance   
    print("\n各主成分贡献率:")
    print(contribution_ratio)
    cumulative = np.cumsum(contribution_ratio)
    print("\n累计贡献率:")
    print(cumulative)
    print("\n===== 计算完成 =====")
    
if __name__ == "__main__":
    main()

简单了解：#

深度了解：#

1、基本样本统计量#

2、定义任意的方向方差#

3、求解方向向量#

4、进行主成分分析#

实操#

代码运行#