主成分分析总结

原创码字不易，如有引用或转载，请注明出处！！！

1. 主成分分析简介

在进行数据分析任务时，我们总会遇到下面一些问题：1. 维数灾难，随着数据维数的增加，机器学习算法处理数据所需要的内存和时间会呈指数增长；2. 数据冗余，很多数据没有明显作用，这些没有明显作用的数据掩盖了重要特征； 3.变量间的相关性，变量之间是互相影响的，当其中的一个变量发生改变时，其他的也会发生相应的变化；4. 数据噪音。噪音是所有数据集中都不可避免的问题，在数据的采集、存储、压缩和传输过程中都有可能产生噪音。对数据进行去噪有助于提高机器学习模型的训练效率和模型的预测能力；5. 数据可视化，由于人脑只能处理三维空间图形，我们是无法直观观察数据在高维空间中是如何分布的，因此如何寻找合适的角度取观察数据是可视化任务面临的一个问题。

针对数据分析中的这些问题，有一种叫主成分分析（Principal Component Analysis，以下简称PCA）的方法可以为上述所有问题提供一种解决思路。PCA有着简洁优美的理论基础，活跃在机器学习的各个领域，是一种非常重要的数据处理工具。

目前，学术界普遍认为PCA是由英国统计学家Karl Pearson和美国统计学家Harold Hotelling提出。Pearson的研究领域十分广泛，是一位百科全书式的学者，被认为是现代统计学的奠基人之一，统计学中的很多公式都用Pearson的名字命名。Hotelling是美国统计学家、经济学家，主要贡献在计量经济学、多元分析和统计推断领域。Hotelling是多元分析领域的领军人物，1970年被评为美国国家科学院院士。

Pearson在1901年发表了一篇文章¹，文章中解决了如何用直线或平面拟合高维空间中的数据样本，这个问题与数据降成一维或二维等价。Pearson的文章从点到直线的距离最小化出发求解最优直线，并且给出了一个数值算法，这是最早的PCA数值算法。

1933年，Hotelling在《British Journal of Educational Psychology》上发表了一篇文章²，与Pearson的几何背景不同，Hotelling从多元正态分布的角度看待主成分的几何含义，从方差最大化的角度研究了PCA，通过Lagrange数乘法求出最优解是特征值和特征向量，除了没有使用矩阵的概念进行运算，文章的方法与现代教材中的PCA非常接近。

PCA在信息损失最小的情况下可以达到消除多重共线性和降维的目的，但是通过PCA降维有四个问题：第一个问题是PCA在降维时没有考虑数据的概率分布，如果是多种分布混合的数据PCA降维效果不好；第二个问题是PCA提取的主成分是原始变量的线性组合，其实际意义很难直观理解，这对模型的解释造成困难；第三个问题是PCA对稀疏的强噪音、异常值或缺失值敏感；第四个问题是PCA只能进行线性降维，对大规模非线性数据降维效果不好。针对上面这些问题，学者们沿着不同的思路对PCA进行改进。

第一个问题，Michael E. Tipping和Christopher M. Bishop于1999年³提出了概率主成分分析（Probabilistic Principal Component Analysis，简称PPCA），可以处理缺失数据和混合模型，PPCA模型中使用了EM算法确定相关参数。第二个问题，统计学家发展的因子分析（Factor Analysis）可以解决隐含变量的直观解释，因子分析也是一种降维方法，从数据中找出影响变量的隐含的因素，通过因子旋转对这些因素进行解释，这些因素可能无法被直接观测，可以根据经验或逻辑分析对这些因素进行验证，因子分析在多元统计分析的教材中都有详细介绍。第三个问题，学者们提出了稳健主成分分析（Robust Principal Component Analysis，简称RPCA），RPCA通过将矩阵分解为一个低秩矩阵加一个稀疏的噪音矩阵进行降维，这种方法要求噪音是稀疏的，对噪音强度没有要求，RPCA在图像降噪和缺失值处理上有很好的效果，Emmanuel J. Candes等人的文章⁴详细介绍了RPCA的方法原理。第四个问题，Bernhard Scholkopf等人在1997年⁵提出了核主成分分析（Kernel Principal Component Analysis，简称KPCA），很多情况下，数据在空间中的分布是非线性的，KPCA首先通过非线性映射将数据映射到高维特征空间中，通过核函数求解变换后的特征值和特征向量。

尽管PCA是一百多年前的提出的统计方法，但是其并没有因为新理论的出现而沉寂，是数据分析学科中一颗愈久弥新的珍珠。

2. 主成分分析的基本思想

我们用一个简单的例子来说明主成分分析的基本思想。有一组二维向量：
$$(-2,-4),(-1,-2),(1,2),(2,4),(3,6)$$

从图像上容易观察出这5个数据在同一条直线 $y=2x$ 上，因此，虽然这5个点是二维数据，但如果我们用向量所在的直线 $y=2x$ 作为坐标轴，只需要一个坐标就可以表示这5个二维数据，同时不损失原数据的任何信息，换句话说就是如果知道坐标轴的方向和点在坐标轴上的位置，我们可以完全用勾股定理重构出原始数据的二维坐标。从线性代数的角度看，这组二维向量其实是在一个一维线性子空间中，找出这组向量所在的直线方向就是数据降维。

上面的例子是一个最简单的理想化情况，如果这组向量并不是在一条直线之上，而是在直线附近出现，是否可以降维？假设这组数据为：
$$(-2,-3.2),(-1,-2.6),(1,1.3),(2,4.9),(3,5.5)$$

观察图像发现，虽然这些点不在任何相同直线上，但是它们都在直线 $y=2x$ 附近，如果我们用 $y=2x$ 作为降维后的子空间，用这些点在直线上的投影作为刻度，在损失了点到直线的距离这部分信息后，仍然可以用一组在同一直线上的新数据近似代替原数据，因为有部分数据损失，新数据无法重构原数据，但是可以作为原数据的近似替换，这样我们就实现了二维数据的降维。这就是Person在论文中使用的降维思想。

但是，上面的过程仍然有几点问题需要明确：一是这些点附近有无数条直线可以作为投影方向，这些方向都可以用来降维，我们应该用什么方法评价这些投影方向的好坏？ $y=2x$ 是否是最好的那一个方向？二是如果有了明确的评价方法，如何计算出最优的投影方向？以上就是主成分分析的基本思想和要解决的核心问题。对于第一个问题，通过数据重构时的信息最大化和损失最小化来确定最优方向。对于第二个问题，通过线性代数的特征值理论可以求出最优方向。

在降维之后，我们只保留投影的信息，忽略与投影垂直的向量信息，如果我们想通过投影重构原向量，则投影的长度越大保留的信息越多，损失的信息越小，因此所有样本投影的长度的平方可以作为衡量投影方向好坏的标准，如下图所示，

由勾股定理，投影长度的平方与损失向量长度的平方和等于原向量长度的平方，所以也可以用损失向量长度的平方和最小化作为衡量标准，两种评价标准其实是同一件事情。有些资料是从方差最大化的角度去求最优投影方向，方差越大的方向包含的信息越多，当方差越小，数据越集中，当方差为0时，在这个方向上所有数据都集中为一个点，不包含任何信息。下图显示了数据在不同方向上的投影：

用这种方法计算出来的主成分有什么含义？主成分可以看成是原变量合成的一些新的特征，如何解释这些新特征需要结合问题的实际背景去分析，大部分情况下我们无法对这些抽象的新特征做出解释。PCA的计算思路可以概括为：1、用一组向量表示投影方向或降维后的空间，2、将样本中心化得到新的样本。3、计算新样本在这个方向或空间中的投影。4、计算投影与原数据的误差平方和。5、求出在什么方向上误差平方和最小。

3. PCA的数学原理

先介绍如何计算向量在 $W$ 上的投影，设 $W^T$ 为 $W$ 的正交补空间，则 $W^T$ 中的向量与 $W$ 中的向量正交。根据线性代数的知识，$m$ 维空间中的向量 $x$ 可以被唯一分解为 $W$ 和 $W^T$ 中的向量之和，$x=u+v$，$u\in W$，$v\in W^T$。向量 $u$ 就是 $x$ 在 $W$ 上的投影。

设 $U=(u_1,u_2,\cdots ,u_k)$ 为 $W$ 一组基构成的矩阵，投影 $u$ 可以表示为 $U$ 的线性组合，设$u={\sum }_{i=1}^k{c}_i{u}_i=Uy$
$y=(c_1,c_2,\cdots ,c_k{)}^T$，只要求出 $y$ 就能求出投影 $u$。$v=x-u=x-Uy\in W^T$，所以 $v$ 与 $W$ 中的向量正交，$U^{T}v=U^T(x-Uy)=(u_1,u_2,\cdots ,u_k{)}^{T}v=(u_1^{T}v,u_2^{T}v,\cdots ,u_k^{T}v{)}^T=0$。$U^{T}x=U^{T}Uy$，$y=(U^{T}U{)}^{-1}{U}^{T}x$，投影$u=Uy=U(U^{T}U{)}^{-1}{U}^{T}x$。当 $U$ 为规范正交基时$U^{T}U$为单位矩阵，$u=U{U}^{T}x$。下面计算中假设 $U$ 为规范正交基。

有了上面的概念后可以开始PCA的推导。首先将PCA算法转换成对应的数学问题：给定一组 $m$ 维样本$Z=\left(Z_1,Z_2,\cdots ,Z_n\right)$，将 $Z$ 中心化得到样本 $X$，要在 $m$ 维线性空间中找出一个线性子空间 $W$，使 $X$ 在 $W$ 上的投影与 $X$ 的误差平方和要尽可能小，同时 $W$ 的维数要尽可能小，要找出 $W$ 的一组基 $u_1,u_2,\cdots ,u_k$。我们按照上面提到的5个步骤寻找最优的投影空间。

将样本 $Z$ 中心化，$X=\left(X_1,X_2,\cdots ,X_n\right)=\left(Z_1-\bar{Z},Z_2-\bar{Z},\cdots ,Z_n-\bar{Z}\right)$，$X$ 满足$\bar{X}=0$。再求 $X$ 中的个体在 $W$ 上的投影，设个体 $X_i$ 的投影为$x_i$，$x_i=U{U}^T{X}_i$，所有投影的平均值$\bar{x}=0$。设 $Err$ 为投影与 $X$ 的误差平方和，则有
$$Err=\sum _{i=1}^n{\left(X_i-x_i\right)}^T\left(X_i-x_i\right)=\sum _{i=1}^n{X}_i^T{X}_i-\sum _{i=1}^n{X}_i^{T}U{U}^T{X}_i$$
上式前一项是常数，因此求Err的最小值等价于求 ${\sum }_{i=1}^n{X}_i^{T}U{U}^T{X}_i$ 的最大值。因为
$$X_i^{T}U{U}^T{X}_i=X_i^{T}U{\left(X_i^{T}U\right)}^T=\sum _{j=1}^k{u}_j^T{X}_i{X}_i^T{u}_j$$
所以，
$$\sum _{i=1}^n{X}_i^{T}U{U}^T{X}_i=\sum _{i=1}^n\sum _{j=1}^k{u}_j^T{X}_i{X}_i^T{u}_j=\sum _{j=1}^k\sum _{i=1}^n{u}_j^T{X}_i{X}_i^T{u}_j=\sum _{j=1}^k{u}_j^T\left(\sum _{i=1}^n{X}_i{X}_i^T\right)u_j$$
样本 $X$ 的协方差矩阵为 $\Sigma =\frac{1}{n-1}{\sum }_{i=1}^n{X}_i{X}_i^T$，所以求 ${\sum }_{i=1}^n{X}_i^{T}U{U}^T{X}_i$ 的最大值等价于求 ${\sum }_{j=1}^k{u}_j^T\Sigma u_j$ 的最大值。又因为 $tr\left(U^T\Sigma U\right)={\sum }_{j=1}^k{u}_j^T\Sigma u_j$，所以问题转化为求矩阵 $U$ 使 $tr\left(U^T\Sigma U\right)$取最大值。

因为 $\Sigma$ 为对称矩阵，由线性代数知识 $\Sigma$ 可以分解为 $\Sigma =Q\Lambda Q^T$，$\Lambda$ 为对角矩阵，对角线上元素为 $\Sigma$ 的特征值，$Q$ 为正交矩阵，$Q$ 的列向量为特征向量。令$\Lambda =diag\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_m\right)$，$C=Q^{T}U=\left(c_1,c_2,\cdots ,c_k\right)$，$c_i={\left(c_{1i},c_{2i},\cdots ,c_{mi}\right)}^T$
$$tr\left(U^T\Sigma U\right)=tr\left(U^{T}Q\Lambda Q^{T}U\right)=tr\left(C^T\Lambda C\right)=\sum _{i=1}^m\sum _{j=1}^k{\lambda }_i{c}_{ij}^2$$
因为 $C^{T}C=U^{T}Q{Q}^{T}U=I_k$，所以 $c_i^T{c}_i=1$，${\sum }_{i=1}^m{c}_{ij}^2=1$
$$\sum _{j=1}^k\sum _{i=1}^m{c}_{ij}^2=k=\sum _{i=1}^m\sum _{j=1}^k{c}_{ij}^2$$
$C$ 是一个 $m\times k$ 矩阵，$C$ 的列向量组是规范正交向量组，可以将 $C$ 的列向量组扩充成 $m$ 维的规范正交向量组，设扩充后的列向量组构成的向量为 $D$，$D$ 的前 $k$ 列是 $C$ 的 $k$ 个列向量，$D$ 是一个正交矩阵，其行向量和列向量都是规范正交向量组，$D$ 的行向量模长为1。因此 $C$ 的第 $i$ 行 $\left(c_{i1},c_{i2},\cdots ,c_{ik}\right)$ 和 $D$ 的第 $i$ 行$$\left(d_{i1},d_{i2},\cdots ,d_{im}\right)=\left(c_{i1},c_{i2},\cdots ,c_{ik},d_{i\left(k+1\right)},\cdots ,d_{im}\right)$$满足如下关系：
$$\sum _{j=1}^k{c}_{ij}^2\le \sum _{j=1}^m{d}_{ij}^2=1$$
令 ${\delta }_i={\sum }_{j=1}^k{c}_{ij}^2$，$tr\left(U^T\Sigma U\right)={\sum }_{i=1}^m{\sum }_{j=1}^k{\lambda }_i{c}_{ij}^2={\sum }_{i=1}^m{\lambda }_i{\delta }_i$，求 $tr\left(U^T\Sigma U\right)$ 最大值问题可以转换成下面的优化问题：
$$\begin{array}{r}min\sum _{i=1}^m{\lambda }_i{\delta }_i\\ s.t\sum _{i=1}^m{\delta }_i=k\\ 0\le {\delta }_i\le 1\end{array}$$
设 ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_m$。可以证明上面的优化问题最大值为${\lambda }_1+{\lambda }_2+\cdots +{\lambda }_k$。证明如下：
$$\begin{array}{rl}& {\lambda }_1+{\lambda }_2+\cdots +{\lambda }_k-{\lambda }_1{\delta }_1-{\lambda }_2{\delta }_2-\cdots {\lambda }_m{\delta }_m\\ =& {\lambda }_1\left(1-{\delta }_1\right)+{\lambda }_2\left(1-{\delta }_2\right)+\cdots {\lambda }_k\left(1-{\delta }_k\right)-{\lambda }_{k+1}{\delta }_{k+1}-\cdots -{\lambda }_m{\delta }_m\\ \ge & {\lambda }_k\left(1-{\delta }_1\right)+{\lambda }_k\left(1-{\delta }_2\right)+\cdots {\lambda }_k\left(1-{\delta }_k\right)-{\lambda }_{k+1}{\delta }_{k+1}-\cdots -{\lambda }_m{\delta }_m\\ =& {\lambda }_k\left(k-{\delta }_1-{\delta }_2-\cdots -{\delta }_k\right)-{\lambda }_{k+1}{\delta }_{k+1}-\cdots -{\lambda }_m{\delta }_m\\ =& {\lambda }_k\left({\delta }_{k+1}+{\delta }_{k+2}+\cdots +{\delta }_m\right)-{\lambda }_{k+1}{\delta }_{k+1}-\cdots -{\lambda }_m{\delta }_m\\ =& \left({\lambda }_k-{\lambda }_{k+1}\right){\delta }_{k+1}+\left({\lambda }_k-{\lambda }_{k+2}\right){\delta }_{k+2}+\cdots +\left({\lambda }_k-{\lambda }_m\right){\delta }_m\\ \ge & 0\end{array}$$
当 ${\delta }_i$ 满足
$${\delta }_i=\{\begin{array}{ll}1,& 1\le i\le k\\ 0,& k+1\le i\le m\end{array}$$
时取最大值。当 $U$ 为 $Q$ 的前 $k$ 列构成的矩阵时，
$$C=Q^{T}U=\left(\genfrac{}{}{0ex}{}{I_k}{O_{\left(n-k\right)\times k}}\right)$$
此时 ${\delta }_i$ 满足上面条件，优化问题取最大值。所以最优的线性子空间是 $U$ 的列向量构成的线性子空间，也就是协方差矩阵 $\Sigma$ 的前 $k$ 个最大的特征值对应的特征向量。
设 $Q=\left(q_1,q_2,\cdots ,q_m\right)$，$q_1$ 是 ${\lambda }_1$ 对应的特征向量，所有样本在 $q_1$ 上的投影被称为第一主成分 $F_1$，$F_1=q_1^{T}Z$。同理， ${\lambda }_2$ 对应的特征向量 $q_2$，所有样本在 ${\lambda }_2$ 上的投影称为第二主成分 $F_2$，$F_2=q_2^{T}Z$，依次类推，可以获得 $m$ 维空间的 $m$ 个主成分，新的主成分为
$$F=\left(\begin{array}{c}{F}_1\\ F_2\\ ⋮\\ F_m\end{array}\right)=Q^{T}Z$$
不同主成分之间正交，因此各主成分之间不存在多重共线性的问题。
有了主成分之后需要确定降维后的维数 $k$ 。从主成分的推导过程可以看出 $\Sigma$ 的特征根大小衡量了包含信息的多少，因此，我们先将特征根从大到小排列 ${\lambda }_1\ge {\lambda }_2\ge \cdots {\lambda }_m\ge 0$，表示将不同主成分的按包含信息的多少从高到低排列，定义
$$\frac{{\lambda }_i}{{\sum }_{i=1}^n{\lambda }_i}$$
为第 $i$ 个主成分的贡献率，定义
$$\frac{{\sum }_{i=1}^k{\lambda }_i}{{\sum }_{i=1}^m{\lambda }_i}$$
为前 $k$ 个主成分的贡献率。降维的过程就是确定选择多少主成分可以尽可能的保存原数据的信息，常用的确定 $k$ 的方法如下：

1. 根据具体问题确定一个阈值，当累计贡献率高于这个阈值时选择最小的维数。
2. 如果数据经过了标准化，可以将特征根大于1作为提取主成分的标准。这是因为当数据经过标准化后，协方差矩阵与相关系数矩阵相同，此时的特征根之和等于 $n$，因此提取的主成分的贡献率应大于所有主成分的均值。
3. 从大到小排列将主成分的贡献率连成一条折线，曲线拐点位置即为的主成分位置，将前几个主成分提取出来。

确定 $k$ 主成分之后，将最大的 $k$ 个特征值对应的特征向量取出，令$Q_k=\left(q_1,q_2,\cdots ,q_k\right)$，最终降维后的结果为$$F=Q_k^{T}Z$$
我们将主成分分析的步骤归纳如下：

1. 将样本中心化，$X=\left(Z_1-\bar{Z},Z_2-\bar{Z},\cdots ,Z_n-\bar{Z}\right)$.
2. 计算样本的协方差矩阵$\Sigma =\frac{1}{m-1}X{X}^T$.
3. 计算协方差矩阵 $\Sigma$ 的特征值和特征向量，将特征值按照从大到小的顺序进行排列${\lambda }_1\ge {\lambda }_2\ge \cdots {\lambda }_m\ge 0$，计算贡献度
   $$\frac{{\lambda }_i}{{\sum }_{i=1}^m{\lambda }_i}$$
4. 确定主成分的个数 $k$，将前 $k$个特征值对于的特征向量取出，$Q_k=\left(q_1,q_2,\cdots ,q_k\right)$
   降维后的 $k$ 个主成分为 $F=Q_k^{T}Z$.

4. 算法实现

主成分分析的步骤并不复杂，涉及到的矩阵运算可用Numpy库中的函数实现。我们通过前面提到的鸢尾花（iris）数据集来验证主成分分析的降维过程。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn import datasets

iris_data = datasets.load_iris()
iris = pd.DataFrame(iris_data.data, columns=iris_data.feature_names)
iris['target'] = iris_data.target
iris.head()

输出

	sepal length (cm)	sepal width (cm)	petal length (cm)	petal width (cm)	target
0	5.1	3.5	1.4	0.2	0
1	4.9	3.0	1.4	0.2	0
2	4.7	3.2	1.3	0.2	0
3	4.6	3.1	1.5	0.2	0
4	5.0	3.6	1.4	0.2	0

接下来尝试主成分分析降维效果：

1. 数据中心化，将每条样本作为列向量构成矩阵 $Z$，每一列数据减去样本平均值 $X=Z-\bar{Z}$.

   # 原数据一行是一条样本，首先将其转换成列向量
   Z = iris_data.data.T
   # 将Z中心化
   X = Z - Z.mean(axis=1).reshape((4, 1))
   

2. 计算协方差矩阵的特征值和特征向量，协方差矩阵 $\Sigma =\frac{1}{n-1}X{X}^T$，特征值和特征向量用numpy.linalg.eig()求出，输入参数是矩阵，输出结果是一个元组，第一个元素是特性值构成的向量，第二个元素是矩阵，矩阵的列向量是特征向量。

   # 协方差矩阵
   Sigma = X.dot(X.T)/(X.shape[1]-1)
   # 将协方差矩阵分解
   r, Q = np.linalg.eig(Sigma)
   

3. 计算累计贡献度，选择主成分。

   r.cumsum()/r.sum()
   

   输出结果:

   array([0.92461872,  0.97768521,  0.99478782,  1.  ])
   

   通过计算发现2个主成分可以达到0.977以上的累计贡献度，3个主成分可以达到0.994以上的主成分，与原数据之间的差别非常小.
4. 计算降维后的坐标，选定 $k$ 个主成分后，计算降维后的数据，设 $Q_k$ 表示特征向量的前 $k$ 列，$Y=Q_k^{T}Z$，矩阵的每一列表示样本在新的方向上的坐标。对于鸢尾花数据，我们选定前两个主成分，则新的主成分 $Q_2^{T}Z$。画出新坐标下的散点图，

   # 计算主成分
   Q2 = Q[:, :2]
   F = Q2.T.dot(Z)
   # 画散点图
   plt.scatter(F[0,:], F[1,:], c=iris['target'], alpha=0.6)
   

   输出结果如下，
   
   从上图可以看到在新坐标下，三类花可以用清晰的划分为3类.
   如果选择3个主成分，画出散点图，结果如下，

   # 计算主成分
   Q3 = Q[:, :3]
   F = Q3.T.dot(Z)
   # 画散点图
   figure = plt.figure(figsize=[5, 5])
   ax = figure.add_subplot(111, projection='3d')
   ax.scatter(F[0,:], F[1,:],F[2,:], c=iris['target'], alpha=0.6)
   

   输出图片
   

5. 小结

主成分分析法是常用的数据处理工具，在使用主成分分析降维时，有一些需要注意的细节。在进行数据分析时，会遇到变量之间数值绝对差距过大的情况，在处理数据时，需要对数据进行标准化。

主成分分析法的一个局限是处理数值型的矩阵，有很多分类问题中，类别是字符型变量，虽然可以用整型代替，但是对这些变量分解或线性变换是没有意义的。

主成分分析的另一个局限是只能处理线性问题，对于非线性数据的降维效果不好，此时可以用核主成分分析。

原创码字不易，如有引用或转载，请注明出处！！！

参考文献

---

1. Pearson K . On lines and planes of closest fit to systems of points in space. ↩︎

2. Hotellings H . Analysis of a complex of statistical variables into principal components. ↩︎

3. Tipping M E , Bishop C M . Probabilistic Principal Component Analysis. ↩︎

4. Candes E J , Xiaodong L I , Yl M A , et al. Robust Principal Component Analysis? ↩︎

5. Scholkopf B , Smola A . Kernel Principal Component Analysis. ↩︎