主成分分析指南

最新推荐文章于 2022-09-28 21:33:13 发布

ybdesire

最新推荐文章于 2022-09-28 21:33:13 发布

阅读量1.8k

点赞数 3

分类专栏： Machine Learning 文章标签：数据分析数学 PCA 机器学习主成分分析

Machine Learning 专栏收录该内容

110 篇文章 17 订阅

订阅专栏

本篇博文根据笔者的理解进行翻译，并不完全按照原作逐字翻译。英文原文见这里

摘要

主成分分析（PCA）是现代数据分析的中流砥柱，但大部分人并不理解它，只是将PCA作为一个黑盒来使用。

本文的主要目的，就是要揭秘PCA背后的魔法。本文能让你直观的理解PCA的作用和原理，与此同时，本文也对PCA背后的数学原理进行了解释。这篇文章并没有觉得用非正式的证明来解释PCA是羞耻的，也没有刻意避免数学公式。作者希望，通过直观解释与数学解释，不同层次的读者都能对PCA有更好的理解，并知道什么时候、为什么与如何运用这种技术。

1. 概述

在应用线性代数领域，PCA已经被认为是一种最好的算法之一。PCA也在分析中被广泛使用：从神经科学到计算机图形学。因为它是一种简单、无参的，能从令人费解的数据集中，提取出有相关性的信息。只需要用最少的数学处理，PCA就能将复杂数据集维度降低。通过降维后的数据，就能看到隐藏在复杂数据集中有价值的信息。

这篇指南文章的目的，有两点：提供PCA的直观理解，并严肃讨论PCA的原理。我们会从一个简单的例子开始，先提供一个对PCA算法目标的直观解释。然后，我们会通过严格的数学描述，与线性代数理论，来彻底解决这个问题。

我们也将看到为什么PCA与奇异值分解（SVD）的关系如此密切。理解了这些，我们就能知道在真实世界中如何正确应用PCA。我们会讨论PCA技术背后的假设，也会告诉大家克服这些局限的方法。

以“指南”的眼光来看，本文的讨论和解释都是非正式的。本文的目的在于“教育”。本文的某些严格数学证明请参考附录。对本文来说，这些严格的数学证明并不是都要读者去理解，它们只是为想深入PCA数学原理的读者准备的。

本文仅假定读者有线性代数基础，并不需要其它深入的数学知识。读者有任何建议、勘误或其它想法，请随时与我联系。

2. 动机：一个最基础的例子

举个例子，假设我们是实验者，我们想通过测量不同的参数（频谱，电压，速度），来了解某种现象。不幸的是，我们有了一堆测量数据，但却没法来认识这个现象。因为数据中有噪声，不明确，甚至有数据冗余。这个例子，并不是一个无价值的问题，这是科学实验中都会遇到的基本障碍。这样的例子在很多复杂系统中都会遇到：图像科学，气象学，海洋学。直接测量得到的数据难以被直接使用，甚至具有很大的迷惑性。而数据背后的潜在事实而又非常简单。

这里举一个最基础的例子，以图1为例。假设我们正在研究物理学中弹簧振子的运动。这套系统由质量为m的球体和无质量、无摩擦系数的弹簧组成。球体已经偏离了平衡状态（弹簧被拉伸）。因为弹簧处于理想情况下，所以球体会沿着x轴以某个频率振动。

这里写图片描述

这是物理学中的一个标准问题。沿着x轴的振动能被表示为时间的函数。换句话说，其动力学方程能被表示为单个变量x的函数。

然而，作为实验者的我们并不知道这个事实。我们并不知道哪些轴（变量）对测量结果有影响。所以我们决定，球体移动过程中，在三个轴上的位移都要测量。所以我们在系统中加上了3个摄像机，它们以200Hz的帧频拍摄视频，每个摄像机能记录球体在三轴系统中两个轴上的投影（见图1，一定要理解这句话）。不幸的是，由于我们的疏忽，我们并不知道哪个摄像头对应x、y、z轴，所以，我们以摄像机方向作为参考系。记住摄像机与各个平面的夹角也不是90度。我们用这三个摄像机记录了2分钟球体的运动。那么，问题来了，我们如何用摄像机记录的数据，得到一个简单的关于x的函数？

如果我们是足够聪明的实验者，我们知道这个先验知识（其动力学方程能被表示为单个变量x的函数），那我们只需要测量球体在x轴的运动。但真实世界中，我们并不知道这样的先验知识。我们常常是不知道测量什么参数才能最好的反应系统的动力学方程。所以，我们会测量更多的参数，而其中很多参数是不需要的。

除此之外，在真实世界中，我们还要和噪声做斗争。在刚才那个例子中，这意味着我们还要处理空气阻力，摄像机的成像误差，以及非理想弹簧的质量和摩擦。这些噪声都会对我们的测量数据产生干扰。这个例子就是实验者们每天都要面对的挑战。在深入抽象概念之前，我们都会以这个例子来讲解。希望在文章的最后，我们能很好的理解如何用PCA来科学的提取出x。

3. 框架：基变换

目标：PCA计算最有意义的基底，并用它来重新表示误差和混乱的数据集。希望新的基底能过滤掉噪声，并将内在的数据释放出来。在弹簧振子的例子中，PCA的终极目标是：“动力学方程仅与x轴有关”。换句话说，PCA的目标是，给出结论：沿着x轴的单位基底向量，是最重要的。能得到这个事实，就能帮助实验者们分辨出哪些数据是重要的，哪些是冗余的。

3.1 原始基底

为了得到我们目标的精确定义，我们也需要对数据集做一个更精确的定义。对每一个实验，实验者们会记录一系列的数据（例如电压、位置）。待测量数据量的个数，就是数据集的维度。在弹簧振子的例子中，数据集有12000 个6维向量组成，这其中每个摄像机记录了球体的两维投影位移。

通常来说，每个数据样本都是m维的向量，m是待测量的数据量。这等效于，每个时刻，样本是一个基于正交基生成的m维向量空间中的一个向量。在这个向量空间中，所有测量得到的向量都是由单位长度为1的正交基生成的。最简单的基底B就是单位矩阵I。

这里写图片描述

单位阵中每个行向量都是都是m维的基底。

总结一下，在弹簧振子例子中，每个时刻，摄像机A记录了球体的位置 $(x_{A}(t), y_{A}(t))$ 。三个摄像机的数据可以用六维列向量来表示。

这里写图片描述

其中，每个摄像机记录了两个位置，所以整个向量X是单位基底B乘以不同的系数。

3.2 基底变换

有了基底这个更严格的定义，我们现在可以更精确的说明PCA的：能不能找到另一组基底，来更好的表示我们的数据集？而这组新基底是原始基底的线性组合。

细心的读者可能会注意到，上一段话中有一个显眼的关键词“线性”。实际上，PCA有一个严格的假设：（新基底由老基底）线性（组合而成）。线性假设在两个方面简化了问题：（1）限制了潜在基底的集合（有限集合）（2）暗含了数据集连续性的假设。线性假设只是一个微妙的点，但我们已经含蓄的说明，例子中测量数据集能表征系统的动力学方程，这就做了线性假设（因为我们已经知道，球体运动的动力学方程能被表示为单个变量x的函数）。换句话说，我们已经知道了先验知识。

有了这个假设，PCA就有了局限，它仅能将数据表征为基底向量的线性组合。用Y表示X基于P的线性变换，X是原始测量数据，Y是数据的重新表示。

这里写图片描述

并做如下定义：

$p_{i}$ 是矩阵 $P$ 的行向量
$x_{i}$ 是矩阵 $X$ 的列向量
$y_{i}$ 是矩阵 $Y$ 的列向量

等式（1）表示了基底变换，因此还能从中看出更多含义：

(1). 矩阵P能将X变换为Y
(2). 在几何上，P只是对X做了旋转、拉伸
(3). P的行向量，就是X列向量的新基底

接下来的解释并没有那么直观，但也能够通过PX的点表示法来说明：

这里写图片描述

Y的每一列可以表示如下：

这里写图片描述

我们能看到， $y_{i}$ 中的每个系数，都是P中行与X中列相乘所得。换句话说， $y_{i}$ 中的第j个系数，是在P的第j行上的投影。这就是最正式的含义： $y_{i}$ 是X在基底上的投影。所以，矩阵P的行向量，就是用于重新表示X列向量的新基底。

3.3 遗留问题

通过将问题限制在线性范围内，从而将问题简化为寻找合适的基底变换。在变换中，行向量 $\{p_{1}, ..., p_{m}\}$ 将变为X的主成分。现在，问题来了。

什么是最好的重新表示X的最好方式？
基底P的最好选择是什么？

这些问题呢，通过问我们自己，我们想让Y表示什么特征，接下来就能得到解决。最终，除了线性假设，还需要其它假设，才能得到合适的结果。如何选择其它假设，是下一部分的主题。

4 差异与目标

现在，最重要的问题来了：什么是数据的“最优表示”？本章会对这个问题提供一个直观的解答，同时给出其它的假设。在本章的结尾，我们也会给出一个数学上的解释“畸变数据”的目标。

在线性系统中，畸变表示两种情况：噪声，冗余。让我们分别处理这两种情况。

4.1 噪声

无论使用什么分析技术，数据中的噪声都必须要足够低，否则是没法从数据中提取出信息的。关于噪声，现实中并没有绝对的噪声衡量方法（区分出每一种噪声）。但可以测量到所有噪声的总和。噪声的通用计算方法是SNR（信噪比）。

这里写图片描述

如果数据的SNR很高（远大于1），说明数据很精确。反之，若数据的SNR很低，说明数据中含有的噪声很高。

假设我们在图二中，画出了来自摄像机A采集的数据。

这里写图片描述

每一个摄像机记录的球体移动轨迹，理论上都是一条直线。因此，不在直线上的数据，就是受到了噪声的干扰。信号与噪声的偏差，都在图中能看到。所以，信噪比SNR也可以根据图中椭圆的形状来判断，形状越瘦，说明信噪比越大；若形状是一个圆形，说明信噪比为1，甚至情况更糟。
总结一下，我们必须要保证我们的测量设备是可靠的；从量化分析上来说，要保证有很高的SNR。

4.2 冗余量

冗余量是一个更棘手的问题。这个问题在弹簧振子例子中是很明显的。在这个例子中，多个传感器共同记录了同一个动力系统的运动信息。假设有两个任意的测量参数r1和r2，见图3.

这里写图片描述

图三(a)中，是两个参数无冗余的情况。换句话说，r1和r2不相关。当你画出两个参量如（ $x_{A}$ ，湿度），就是这种情况。
图三(c)中，两个参量高度相关。这种情况可能有以下可能造成：

两个参量其实是同一个数据，只是单位不同
两个摄像机离的太近

图三(c)中，记录了r1和r2的线性关系，这比单纯记录一个参量是更有意义的。它即能更简洁的表示数据，又能降低传感器的数量。实际上，这就是维度降低背后的基本思想。

4.3 协方差矩阵

信噪比SNR是通过计算方差得到的。一种等效，但更简洁的衡量两组向量之间冗余量的方式，是进行方差计算。

考虑两组测量系统同时记录的数据（做了零均值处理后的）。

A = {a 1, a 2, a 3, . . ., a n}

$A=\{a_{1},a_{2},a_{3},...,a_{n}\}$

B = {b 1, b 2, b 3, . . ., b n}

$B=\{b_{1},b_{2},b_{3},...,b_{n}\}$

A和B的方差用如下公式定义

δ 2 A = < a i a i > i

$\delta_{A}^{2} = <a_{i}a_{i}>_{i}$

δ 2 B = < b i b i > i

$\delta_{B}^{2} = <b_{i}b_{i}>_{i}$

其中，期望是n个变量的平均数。A与B的协方差用如下公式计算

A 与 B 的 协 方 差 = δ 2 A B = < a i b i > i

$A与B的协方差 = \delta_{AB}^{2} = <a_{i}b_{i}>_{i}$

关于协方差，有两点要注意的：

$\delta_{AB}^{2}=0$ ，说明A与B不相关
$\delta_{AB}^{2}=\delta_{A}^{2}$ ，说明A=B

所以，我们现在可以用矩阵表示协方差：

这里写图片描述

最终，我们能从两个向量的计算中，得到一个数值结果。我们把行向量x1, x2重命名为a, b，然后增加其他行向量x3, …, xm。于是我们得到一个新的m行n列矩阵X。

这里写图片描述

矩阵X的解释如下：X中的行代表测量某种类型（xi）得到的数据；列代表每一次实验得到的测量数据（3.1节中的向量X）。因此我们得到协方差矩阵 $S_{X}$ 的定义如下。

这里写图片描述

$XX^{T}$ 就是计算 $S_{X}$ 第(i,j)位置的值。 $S_{X}$ 第(i,j)个值，就是第i次测量值的向量，与第j次测量值的向量乘积。
* $S_{X}$ 第(i,j)个值，等效于将xi，xj带入式（3）。
* $S_{X}$ 是一个m行m列的矩阵。
* $S_{X}$ 的对角线元素，就是测量值之间的方差。
* $S_{X}$ 的非对角线元素，就是测量值之间的协方差。

计算 $S_{X}$ ，就等效于计算所有测量值之间的相关性。任意两次测量值，若它们的协方差太大，就等效于图3的情况(c)；协方差为0，就等效于图3的情况(a)。

$S_{X}$ 是很特殊的。协方差矩阵描述了各测量值之间的所有关系。假设我们可以操纵 $S_{X}$ ，那我们就能定义一个经过操纵后的协方差矩阵 $S_{Y}$ 。那么问题来了，什么样的特征是我们想在 $S_{Y}$ 中优化的呢？

如果我们的目标是降低冗余量，我们希望每个变量都与其他变量相关性越低越好。更精确一点说，为了去除冗余量，我们希望每两个变量间的协方差都为0。经过优化后的协方差矩阵 $S_{Y}$ 长什么样呢？很明显， $S_{Y}$ 的非对角线元素，就是测量值之间的协方差，都为0，这就表示冗余量都被去掉了。

学界有很多得到 $S_{Y}$ 的方法，PCA只是其中最容易的方法，这也是PCA得到广泛应用的原因。

PCA假设所有的偏置向量 $\{p_{1},p_{2},...,p_{m}\}$ 都是正交的，即 $p_{i}p_{j}=\delta_{ij}$ 。在线性代数中，PCA假设P是正交矩阵。

其次，PCA假设方差最大位置的向量，是最重要的向量，即主成分。为什么说这些假设使得问题得到简化呢？

想象一下PCA是如何工作的。通过方差假设，PCA首先选取了m维空间中的归一化备选向量。再者，它还找到了所有的备选向量选择方式，即方差最大。然而，由于正交话条件，它限制了找主成分的范围，只能在上一步方差假设中选择到的备选向量中选取。

原则上，这个简单算法是可行的。然而，它工作的原因是有了明智的正交性假设。也就是说，这个假设的真正好处，是它使得我们的解决方案适应线性代数。这里存在高效、明确的代数解法。

注意我们也从第二个假设（方差最大位置为主成分）中获益，我们得到了一个评价数据“重要性”的方法。也就是说，与每个方向 $p_{i}$ 相关联的方差，决定了主成分的选取。我们因此可以对偏置向量根据相应的方差值大小进行排序。

我们现在暂停讨论达到这个数学目标所做的所有假设的影响。

4.5 假设条件极其限制的总结

本节带大家理解，什么时候PCA会表现不佳。理解了这点，你才会对PCA有进深入的认识。最后一节也提供了最近对于PCA扩展的讨论。

（1）线性
线性将问题变为偏置的改变。有几个领域的研究，已经探讨了如何用非线性先验知识做PCA计算，从而扩展算法，这被称为核PCA。

（2）
均值和方差是统计学的概念。已经有足够的统计数据表明，均值和方差足够描述一个概率分布。能够完全用方差描述的，并满足0均值的概率分布，就是高斯分布。

为了满足这个假设，xi的概率分布必须是高斯分布。不满足高斯分布，则会使这个假设无效。另一方面，这个假设也保证了SNR与协方差矩阵能完全表征噪声和冗余。

（3）方差越大，重要性越大
这个假设也包含了这样一个信息：数据的信噪比很高。方差越大的主成分，重要性越大。反差小的，就是噪声。

（4）主成分是正交的
这个假设提供了一种直观的简化，让PCA与线性代数解耦技术融合在一起。这些技术会在下面两节重点介绍。

我们已经讨论了PCA的方方面面，剩下的就是怎么用线性代数求解了。下面给出两个解法，第一个比较直接，第二个涉及到理解一个重要的解耦技术。

5. 求解PCA：协方差的特征向量

我们先用线性代数推导PCA的第一种数值解法。这种解法是基于特征向量分解的重要属性。这里再一次说明，数据集是X，一个m行n列的矩阵，m是测量数据类别的个数，n是数据测量实验的次数。目标总结如下：

找到正交矩阵 $P$ ， $Y=PX$ ，其中 $S_{Y}\equiv \frac{1}{n-1}YY^{T}$ 是对角阵。则 $P$ 的行就是X的主成分。

我们以选择 $P$ 的方式，来重写 $S_{Y}$ 。

这里写图片描述

注意我们在此定义了一个新矩阵A，根据相关定理，A是对称阵。

对称阵A可以由正交阵将其特征向量对角化得到。

这里写图片描述

其中D是对角阵，E由A的特征向量组成。

矩阵A有r（r<=m）个正交特征向量，其中r是矩阵的秩。当A退化，或选取的r小于m时，A的维度低于m。我们可以通过选择(m-r)个附加的正交向量填满矩阵E，来弥补“保持约束的正交性”这种情况。而这些额外的向量，不影响最终解决方案，因为方差与这些向量的关联是0。

下面，好把戏要来了。P中的每一行 $p_{i}$ ，是矩阵 $XX^{T}$ 的特征向量。所以 $P=E^{T}$ 。由式（6）可知 $A=P^{T}DP$ ，因此

这里写图片描述

很明显，P的选择，使得 $S_{Y}$ 对角化。这就是PCA的目标，我们可以用P和 $S_{Y}$ 来总结PCA。

X的主成分是 $XX^{T}$ 的特征向量，或者是P的行；
$S_{Y}$ 的第i个值，是X沿着 $p_{i}$ 的方差；

在实践中，计算数据集X的PCA，需要：（1）减去每种测量值的平均值；（2）计算 $XX^{T}$ 的特征向量。这个解决方案的演示，在附录B的matlab代码中。

6 更通用的解法：SVD（奇异值分解）

本节含有很多数学公式，即便你跳过本节，也不会有太多损失。它只是为了提供完整性而存在。我们在这里提供PCA的另一种数学解法，PCA与SVD是很接近的。事实上，它们两者是如此密切，以至于经常互相使用。我们会看到，SVD是理解基底变换的一种更通用的方法。

我们从快速导出奇异值分解开始。在下一节中，我们解释奇异值分解，最后一节中再降它们的结果与PCA关联到一起。

6.1 SVD（奇异值分解）

令X是任意一个 $n \times m$ 的矩阵。 $XX^{T}$ 是一个对称方阵（ $n \times n$ ），其秩为r。让我们定义一些接下来会用到的概念：

$\{\widehat{v_{1}}, \widehat{v_{2}}, ..., \widehat{v_{r}}\}$ 是正交集，它们是由对称阵 $XX^{T}$ 的m+1个具有特征值xxx的特征向量所对应的正交集。

$(XX^{T})\widehat{v_{i}} = \lambda_{i}\widehat{v_{i}}$

\sigma {i} = \sqrt{\lambda {i}}，是正实数，也就是所谓奇异值。
$\{\widehat{u_{1}}, \widehat{u_{2}}, ..., \widehat{u_{r}}\}$ 是 $n \times 1$ 向量的正交集，定义为： $\widehat{u_{i}}=\frac{1}{\sigma _{i}}X\widehat{v_{i}}$ 。

我们通过参考附录A的定理5，得到最终定义。最终定义包括两个新的，意想不到的特性（如下）。

$\widehat{u_{i}}\widehat{u_{j}}=\delta _{ij}$
$\parallel X\widehat{v_{i}} \parallel = \sigma _{i}$

这两个特性都被定理5所证明。我们现在有了所有进行分解需要的工具。奇异值分解的数值化版本（如下）只是定义3的重述。

这里写图片描述

上式这个结果就足以说明很多。X乘以 $XX^{T}$ 的特征向量，等效于标量乘以另一个向量。特征向量 $\{\widehat{v_{1}}, \widehat{v_{2}}, ..., \widehat{v_{r}}\}$ 的集合，以及向量 $\{\widehat{u_{1}}, \widehat{u_{2}}, ..., \widehat{u_{r}}\}$ 的集合，都是r维空间的正交集。

我们先构建一个新的对角阵 $\sum$ 。

这里写图片描述

其中 $\sigma _{\tilde{1}}>\sigma _{\tilde{2}}>...>\sigma _{\tilde{r}}$ 是按奇异值大小排序的。因此我们构建了相应的正交阵V和U。

这里写图片描述

其中，我们还在V和U中填充了（m-r），（n-r）个正交向量，以填满V和U。图4对这个填充过程进行了直观的展示。

这里写图片描述

其中V和U的每一列，都进行了“数值”版本的分解（等式7）。因为V是正交的，所以 $V^{-1}=V^{T}$ ，我们可以将等式变换为：

这里写图片描述

这个分解是很强大的，虽然暂时看不出来它的强大。等式9也说明了，任意矩阵X，都能被分解为一个正交阵，一个对角阵，与另一个正交阵的乘积。理解等式9是理解下一部分的关键。

这里写图片描述

6.2 SVD解释

SVD的最终形式（等式9）很简洁，但不容易理解。让我们重新解释等式7。

X a = k b

$Xa=kb$

其中a与b是列向量，k是一个常数。集合 $\{\widehat{v_{1}}, \widehat{v_{2}}, ..., \widehat{v_{m}}\}$ 相当于a，集合 $\{\widehat{u_{1}}, \widehat{u_{2}}, ..., \widehat{u_{n}}\}$ 相当于b。独特的地方在于， $\{\widehat{v_{1}}, \widehat{v_{2}}, ..., \widehat{v_{m}}\}$ 与 $\{\widehat{u_{1}}, \widehat{u_{2}}, ..., \widehat{u_{n}}\}$ 是向量的正交集，这能相应的被扩展到m维或n维空间中。这些集合能扩展到所有可能的输入(a)和输出(b)。我们能否公式化 $\{\widehat{v_{1}}, \widehat{v_{2}}, ..., \widehat{v_{m}}\}$ 与 $\{\widehat{u_{1}}, \widehat{u_{2}}, ..., \widehat{u_{n}}\}$ 能扩展所有的输入和输出？

对式9进行一些调整，使上面这个模糊的假设更精确。

X = U Σ V T

$X=U\Sigma V^{T}$

U T X = Σ V T

$U^{T}X=\Sigma V^{T}$

U T X = Z

$U^{T}X=Z$

这里我们定义了 $Z=\Sigma V^{T}$ 。注意之前的 $\{\widehat{u_{1}}, \widehat{u_{2}}, ..., \widehat{u_{n}}\}$ 现在是 $U^{T}$ 。与式1进行比较，这里的 $\{\widehat{u_{1}}, \widehat{u_{2}}, ..., \widehat{u_{n}}\}$ 与 $\{\widehat{p_{1}}, \widehat{p_{2}}, ..., \widehat{p_{m}}\}$ 相同。因此， $U^{T}$ 是X变换为Z的基底变换。与之前相同，我们在转换列向量。事实上，正交基底 $U^{T}$ 转换列向量，意味着， $U^{T}$ 是一个基底，能将X的列进行扩展。基于这个扩展，这些列向量变为X的列空间。列空间将什么是任何矩阵的可能输出，进行了公式化。