PAC、SVD以及它们的联系之充分数学知识推导

PCA和SVD之间的关系

概述

主成分分析（Principal Component Analysis,简称 PCA）是一种使用广泛的技术，适用于各种应用，例如降维、数据压缩、特征提取和可视化。基本思想是数据集原始的维度可能有许多是线性相关维度的。将数据集投影到被称为主成分（Principal Component，简称PC）的维度上。各主成分是线性无关的，同时保留数据中大部分的可变性（variability）。

奇异值分解(Singular Value Decomposition，简称SVD)是一种经常用于计算数据集主成分的计算方法。采用奇异值分解(SVD)进行主成分分析是一种高效且数值鲁棒的方法。此外，它们之间的亲密关系可以让我们直观地了解我们对PCA实际做了什么，并帮助我们获得对这种技术的额外见解。

在这篇文章中，我将明确描述SVD和PCA之间的数学关系，并强调这样做的一些好处。如果您过去使用过这些技术，但不确定它们如何在内部工作，那么本文将为您提供帮助。最后，您可以理解使用PCA和SVD的动机，并且对如何有效地使用它们有更好的直觉。

数据矩阵化

在我们讨论如何实际执行 PCA 之前，将数据以矩阵形式编码将很有用。由于 PCA 和 SVD 是线性技术，这将使我们能够更轻松地使用线性变换来处理数据。

用例子来引入是最简单的，假设提供了 $n=1000$ 个矩形盒子的宽度 $w$ ,高度 $h$ 和长度 $l$ 。我们将第 $i$ 个矩形盒子编码成 一个 $(w_i,h_i,l_i)$ 元组（tuple），被称为一个样本（sample）。每个样本是 $d=3$ 维的向量（vector），因为有3个分量来描述它。向量 $x_i$ 及 $x_i^T$ 如下所示：
$$x_i=\left(\begin{array}{c}{w}_i\\ h_i\\ l_i\end{array}\right)\text{and}{x}_i^T=\left(w_i,h_i,l_i\right).$$
为了将数据打包到单个对象中，我们只需将样本作为数据矩阵的行进行堆叠，
$$X=\left(\begin{array}{ccc}{w}_1& h_1& l_1\\ w_2& h_2& l_2\\ & ⋮& \\ w_{1000}& h_{1000}& l_{1000}\end{array}\right).$$
这个矩阵的一列对应数据的一个维度，也被称为样本的一个特征。

在一般情况下，我们使用的是数据集有 $n$ 个样本， $d$-维。

我们不使用像$h$、$w$和$l$这样的字母来表示不同的维度，而是简单地枚举每个向量的分量，$x_i^T=(x_{i1},\dots ,x_{id})$ 。这里的向量 $x_i^T$ 表示一个样本， 而它的一个分量表示在一个维度上的取值，比如 $x_{i1}$ 表示第 $i$ 个样本的第 $1$ 个维度的取值。在将这些向量放入数据矩阵之前，为了方便之后的操作，我们实际上要对数据集做中心化（zero-centered），减去来自所有样本数据的平均值
$$\mu =\frac{1}{n}\sum _{i=1}^n{x}_i=(\frac{1}{n}\sum _{i=1}^n{x}_{i1},\dots ,\frac{1}{n}\sum _{i=1}^n{x}_{id}{)}^T$$
得
$$X=\left(\begin{array}{ccccc}& & x_1^T-{\mu }^T& & \\ & & x_2^T-{\mu }^T& & \\ & & ⋮& & \\ & & x_n^T-{\mu }^T& & \end{array}\right).$$
图1给出了$n=1000$， $d=2$的可视化说明。将数据放入矩阵特别方便，因为它允许我们将样本协方差用矩阵表示为
$$S=\frac{1}{n-1}\sum _{i=1}^n(x_i-\mu )(x_i-\mu {)}^T=\frac{1}{n-1}{X}^{T}X.$$
除以$n-1$是典型的使用样本均值而不是真实总体均值的方法。
协方差矩阵将在我们理解主成分和奇异值分解的过程中占据中心地位。

什么是协方差

$X$的第 $j$ 列是上文中数据集中心化后数据集的第 $j$ 维度。
因此，$\frac{1}{n-1}{X}^{T}X$的第 $jk$ 项通过：$X$的第 $j$ 列 和$X$的第 $k$ 列做內积得到的，前者可表示为$x_{\bullet ,j}$，后者可表示为$x_{\bullet ,k}$。即，
$$\frac{1}{n-1}{x}_{\bullet ,j}\cdot x_{\bullet ,k}=\frac{1}{n-1}{x}_{\bullet ,j}^T{x}_{\bullet ,k}=\frac{1}{n-1}\sum _{i=1}^n{x}_{ij}{x}_{ik}.$$
当$k=j$时，这是第 $k$ 个维度的方差。

${\sum }_{i=1}^n(x_i-\mu )(x_i-\mu {)}^T=X^{T}X$是一个简单的矩阵乘法。注意，等号左边的项是表示为$d\ast 1$和$1\ast d$矩阵的向量乘积的和，产生大小为$d\ast d$的结果。

注:概率统计中讲协方差是针对不同随机变量的。在本文中，数据集的一个维度（一个特征）是一个随机变量。而不同样本在同一个随机变量中会有不同的观测值。

主成分分析

最早介绍 PCA 的论文之一是由 Hotelling 于 1933 年发表的。作者的动机是将一组可能线性相关的变量转换为“一些更基本的自变量集…决定了[原始变量]将采用的值。” 受到心理学影响，他命名它们的首选是“因子” （factors），但考虑到这个术语在数学中已经有了意义，他将简化的变量集称为**“成分”**（components），并将找到它们的技术称为“主成分分析”。这些成分是按顺序选择的，让它们“对[原始变量]方差的贡献尽可能大。”

数学上，主成分分析的目标是找到一个集合，它由 $k\le d$ 个单位向量 $v_i\in {\mathbb{R}}^d（i\in 1,\dots ,k$)构成。每个单位向量被称为主成分。满足，

1. 投影到 $v_i$ 决定的方向上的数据集的方差是最大的
2. $v_i$与$v_1，v_2,\dots ,v_{i-1}$正交。

现在，向量 $x\in {\mathbb{R}}^d$在由 $v_i$ 所决定的直线上的投影被表示为内积$v_i^{T}x$ (单位向量模长为1)。这意味着数据集的投影到第一个主成分$v_1$上的方差可以写成
$$\frac{1}{n-1}\sum _{i=1}^n(v_1^T{x}_i-v_1^T\mu {)}^2=v_1^{T}S{v}_1.$$
为了找到$v_1$，我们必须最大化这个量，满足额外单位向量约束，即$\Vert v_1\Vert =1$。

为了解决该优化问题，使用拉格朗日乘子法，得
$$S{v}_1={\lambda }_1{v}_1,\text{\hspace{1em}(1)}$$
这就意味着$v_1$是协方差矩阵$S$的特征向量。事实上，由于$\Vert v_1\Vert =v_1^T{v}_1=1$，在式（1）等号两边同时左乘$v_1^T$ ，我们也得出相应的特征值正好等于数据集沿$v_1$的方差，即，
$$v_1^{T}S{v}_1={\lambda }_1.$$
您可以继续这个过程，将数据集投射到新的方向$v_2$，同时强制附加的约束$v_1\perp v_2$，然后再强制$v_3\perp v_1、v_2$时将数据投射到$v_3$

最终的结果是$X$的前$k$主成分完全对应于协方差矩阵$S$的特征向量按其特征值大小排序。

此外，特征值完全等于数据集投影到相应特征向量方向上计算出来的方差。

您可能已经注意到，这个结果表明$S$ 存在一组完整的标准正交特征向量。实际上，**因为$S$是一个实对称矩阵，这意味着$S=S^T$，事实上谱定理（THE SPECTRAL THEOREM）
说明了这一点。这是一个非常重要的结果，我们将在后面的文章中使用它，所以让我们对它进行一点扩展。

注：

• 实对称矩阵的不同特征值所对应的特征向量是正交的。
• 实对称矩阵的特征值都是实数，特征向量都是实向量。
• n阶实对称矩阵必可对角化。

考虑 $k=d<n$ 的情况。如果我们的目的是通过一个低维子空间来理解$X$，那么使用$k=d$个主成分可能看起来是一个奇怪的选择。但是这样做可以让我们构造一个 $d\times d$ 的矩阵$V$，它的列是$S$的特征向量，因此它可以对角化 $S$，即:
$$S=V\Lambda V^T=\sum _{i=1}^d{\lambda }_i{v}_i{v}_i^T=\sum _{i=1}^r{\lambda }_i{v}_i{v}_i^T,$$
其中$\Lambda =\text{diag}({\lambda }_1,\dots ,{\lambda }_d)$ 而且 $r=\text{rank}(X)$。注意，上式成立是由于$S$是实对称矩阵，附录A有详细说明。由附录B，可知${\lambda }_{r+1}={\lambda }_{r+2}=\dots {\lambda }_d=0$。

Figure 1: 从多元高斯分布中抽取的1000个样本的原始和不线性相关的观点。

奇异值分解

奇异值分解是一种矩阵分解方法，用于线性代数的许多数值应用，例如 PCA。这种技术增强了我们对主成分是什么的理解，并提供了一个强大的计算框架，让我们能够为更多数据集准确计算它们。

让我们先回顾一下针对任意一个$n\ast d$矩阵$ A $的 SVD。
SVD的动机是基于这样一个事实，当 $A$ 被看作一个线性变换时，$A$ 将单位球 ${\mathbb{S}}^d\subset {\mathbb{R}}^d$ 映射到 ${\mathbb{R}}^n$ 中的(超)椭圆。让我们考虑一个$n=d=2$的例子，以便更容易地理解这一事实

图2：在2维欧式空间对单位球体进行[1 2; 0 1]的线性变换

从这个图中，我们可以提取出一些有用的定义，这些定义适用于任意维度，但需要注意以下几点:

• 椭圆的半轴$A{\mathbb{S}}^d在{\mathbb{R}}^n$中的长度${\sigma }_i$是**$A$的奇异值**。

• 沿椭圆半轴的单位向量$u_i$称为$A$的“左”奇异向量。

• （很重要）单位向量$v_i$ 满足
  $$A{v}_i={\sigma }_i{u}_i\text{\hspace{1em}(2)}$$

  $v_i$称为$A$的“右”奇异向量。

注：“左”奇异向量在等号右边，“右”奇异向量在等号左边，请注意。

证明奇异值和奇异向量的存在性和唯一性有点超出了本文的范围，但是证明是相对简单的，并且可以在任何关于数值线性代数的书中找到。然而，图2中的二维图像有望提供足够的提示，使结果变得直观。
一般来说，存在 $r=\text{rank}(A)$个奇异值和相同数量的左奇异向量。准确地说，$r=\text{rank}(A)$是非零奇异值的数目。

将向量$v_i$和$u_i$分别叠加到矩阵$\hat{V}$和$\hat{U}$的列中，可以将$A{v}_i={\sigma }_i{u}_i$写成
$$A\hat{V}=\hat{U}\hat{\Sigma },$$
其中$\hat{\Sigma }=\text{diag}({\sigma }_1,\dots ,{\sigma }_r)$.

个人理解，矩阵方程的关键在于形式是矩阵，其本质还是内部的线性方程组，只不过将内部的线性方程组组成成矩阵方程的形式。然后可以对其进行一些矩阵分析的操作。

通过在$\hat{\Sigma }$中填充0并在$\hat{U}$和$\hat{V}$中添加任意正交列，我们得到了一个更方便的分解式子
$$AV=U\Sigma .$$
因为$V$是幺正（unitary）的——它有单位长度，正交列，满足$V^{-1}=V^T$，那么乘以$V^T$，我们得到$A$的奇异值分解:
$$A=U\Sigma V^T.$$
图2是这个结果的视觉总结。

图3：矩阵A的奇异值分解

注意这个结果最终和上面的椭圆图是一致的，因为三个变换的效果是：用$V^T$切换到新的基，用$\Sigma$拉伸轴得到（超）椭圆，用$U$旋转和反射（不影响椭圆的形状）。

SVD和PCA的关系

因为任何矩阵都有奇异值分解，让我们取$ A = X$ 并记
$$X=U\Sigma V^T.$$
到目前为止，我们认为$A$是一个线性变换，但是没有什么可以阻止我们在任何数据矩阵上使用SVD。**注：**后者是以数据样本集的视角而不是线性变换的视角。

事实上，请注意，从分解我们有
$$X^{T}X=(U\Sigma V^T{)}^T(U\Sigma V^T)=V{\Sigma }^T{U}^{T}U\Sigma V^T=V({\Sigma }^T\Sigma )V^T,$$
这意味着$X^{T}X$和${\Sigma }^T\Sigma$是相似的。

相似矩阵具有相同的特征值，因此协方差矩阵$S=\frac{1}{n-1}{X}^{T}X$的特征值${\lambda }_i$与矩阵$X$的奇异值${\sigma }_i$有关， 通过
$${\sigma }_i^2=(n-1){\lambda }_i,\text{\hspace{1em}(3)}$$
i ∈ { 1 , … , r } i \in \{1, \ldots, r\} i∈{1,…,r}, $r=\text{rank}(A)$.

注：$\Sigma$是对角矩阵。$X$是半正定，对称矩阵。对称矩阵的相应奇异值等于相应特征值的绝对值。

再注：${\lambda }_i$ 是协方差矩阵$S$ 的特征值， ${\sigma }_i$ 是 数据矩阵$X$的奇异值，它们不是针对同一对象的。

为了充分联系SVD和PCA，我们也必须描述主成分和奇异向量之间的对应关系。

根据式（2），可以先取右奇异向量，再得到左奇异向量。式（2）的等号右边理解为左边的值。

对于右奇异向量，取
$${\hat{V}}^T=\left(\begin{array}{ccccc}& & v_1^T& & \\ & & ⋮& & \\ & & v_r^T& & \end{array}\right)\text{\hspace{1em}(4)}$$
其中$v_i$是$X$的主成分。

对于我们取的左奇异向量，在式（2）等号两边同时除以${\sigma }_i$，又根据式（3）可得
$$u_i=\frac{1}{\sqrt{(n-1){\lambda }_i}}X{v}_i.\text{\hspace{1em}(5)}$$
在证明这些选择是正确的之前，我们先来验证一下它们至少在直觉上是有意义的。

从图2中可以看出，右奇异向量 { v i ∣ \{v_i\big| {vi​∣∣​ , i ∈ { 1 , … , r } } i \in \{1, \ldots, r\}\} i∈{1,…,r}}是一个张成（span）数据集行空间的$d$维向量的标准正交集合。

我们还看到$X$的列空间是由左奇异向量$u_i$张成的。

数据矩阵的列空间是数据集不同维度的汇总，因此我们选择$u_i$作为$X$在第$i$个主成分方向上的(缩放)投影是有意义的。

列空间： 因为$X$的秩为$r$ ，所以前 $r$ 个左奇异向量 { u 1 , u 2 , . . . , u r } \{ u_{1},u_{2},...,u_{r}\} {u1​,u2​,...,ur​} （$U$的列向量），为$A$的列空间提供一组标准正交基。

行空间： 类似地，所以前 $r$ 个右奇异向量 { v 1 , v 2 , . . . , v r } \{ v_{1},v_{2},...,v_{r}\} {v1​,v2​,...,vr​} （$V$的列向量），为$A$的行空间提供一组标准正交基。

零空间： 剩下的右奇异向量 { v r + 1 , v r + 2 , . . . , v d } \{ v_{r+1},v_{r+2},...,v_{d}\} {vr+1​,vr+2​,...,vd​} （$V$的剩下的列向量）,为A的零空间提供了一组标准正交基。

上面的证明见Statistical Modeling and Analysis of Neural Data (NEU 560)

这些观察结果令人鼓舞，但并不精确，所以让我们实际证明我们的***选择***在代数上是正确的。

结果需要一个引理。

引理： 设$A=U\Sigma V^T$是实矩阵$A$的奇异值分解。则,
$$A=\sum _{i=1}^r{\sigma }_i{u}_i{v}_i^T.\text{\hspace{1em}(6)}$$
注：这里求和的上界是$r=rank(A)$。引理可在附录A的参考文献中找到。但实际上根据附录B的定理，联系A的奇异值分解展开矩阵乘法即可立即得到。

让我们把这个引理应用到$X$上。 我们得到

第一个等号由式（6）得到，第二个等号由式子（4），（5）得到

第四个等号由 $I=V^{T}V={\sum }_{i=1}^r{v}_i{v}_i^T$ 得到。

我们由此建立了SVD和PCA之间的联系，关键在于式（4），（5）

便利

让我们探讨一下这种对应关系的一些结果。
根据 $U$ 的定义，我们已经知道它将数据投影并稍微缩放到主成分上。
然而，将这个结果推导出为能够写出 $X=U\Sigma V^T$ 的一般结果更令人满意，所以让我们注意这一点。
将$X$的SVD乘以$V$我们得到
$$XV=U\Sigma V^{T}V=U\Sigma .$$
从等号左边我们可以看到这个矩阵的第 $i$ 列是数据集在第 $i$个主成分上的投影
$$z_i=X{v}_i=((x_1-\mu {)}^T{v}_i,\dots ,(x_n-\mu {)}^T{v}_i{)}^T$$
从等号右边我们看到这对应于 $z_i=u_i{\sigma }_i$，这是我们试图恢复的。

我们可以通过考虑$U$的提取另一个有用的部分。

回想一下，每个奇异值对应于椭圆的一条半轴的长度。椭圆描述$X$的列空间。

通过除以奇异值(并乘以$\sqrt{n-1}$)，我们得到了一个转换后的数据集，它在某种意义上是球形的:
$$Y=\sqrt{n-1}U{\Sigma }^{-1}.$$

Figure 4: 白化可用于轻松减轻以不同单位给出的坐标的影响。在左图中，您可以看到 X_2与 X_1 的比例不同。中间图描绘了居中和缩放的数据集。在右图中，您可以看到白化的结果。

最后，让我们讨论一下使用SVD的数值优势。

在计算机上实际计算PCA的基本方法是直接执行$X^{T}X$的特征值分解。

结果表明，这样做会引入一些潜在的严重数值问题，而SVD可以避免这些问题。

问题在于，在计算机上实数以有限精度表示。
由于每个数字都存储在有限的内存中，因此连续可表示的数字之间必然存在间隙。
例如，对于双精度浮点数，实数与其最接近的浮点近似值之间的相对误差约为 $\epsilon \approx 10^{-16}$。
考虑到这种限制的算法称为向后稳定算法。
基本思想是，这样的算法为输入值生成正确的输出，该输入值在您实际询问的值的 $\epsilon$ 内。

如果 ${\tilde{\sigma }}_i$ 是计算 $X$ 奇异值的后向稳定算法的输出，而 ${\sigma }_i$ 是真正的奇异值，则可以证明
$$|{\tilde{\sigma }}_i-{\sigma }_i|=O(\epsilon \Vert X\Vert ),$$

$\Vert X\Vert$ 是 $X$ 大小的合理度量。
另一方面，计算 $X^{T}X$ 的特征值 $l_i$ 的后向稳定算法可以证明满足
$$|{\tilde{l}}_i-l_i|=O(\epsilon \Vert X^{T}X\Vert )=O(\epsilon \Vert X{\Vert }^2),$$
当以 $X$ 的奇异值给出时，它变成
$$|{\tilde{\sigma }}_i-{\sigma }_i|=O(\epsilon \frac{\Vert X{\Vert }^2}{{\sigma }_i}).$$
这里的最终结果是，如果您对相对较小的奇异值感兴趣，例如 ${\sigma }_i\ll \Vert X\Vert$，直接使用 $X$ 会比使用 $X^{T}X$ 产生更准确的结果。
幸运的是，SVD 可以让我们做到这一点。

附录

附录A

定理：设Ａ 为ｎ 阶实对称矩阵，则存在一个正交矩阵$P=(p_1,p_2,\dots ,p_n)$使得
$$A=P\Lambda P^T\text{\hspace{1em}(1)}$$
其中$\Lambda$ 为对角元素是$A$ 的特征值${\lambda }_1,{\lambda }_2,\dots {\lambda }_n$所构成的对角阵；$p_1,p_2,\dots p_n$为特征值${\lambda }_1,{\lambda }_2,\dots {\lambda }_n$ 对应的单位特征向量。式（１）还可以写成如下的形式：
$$A=P\Lambda P^T={\lambda }_1{p}_1{p}_1^T+{\lambda }_1{p}_2{p}_2^T+\cdots +{\lambda }_1{p}_n{p}_n^T=\sum _{i=1}^d{\lambda }_i{p}_i{p}_i^T．$$
[1]张丽静,刘白羽,申亚男.实对称矩阵对角化教学的应用案例[J].大学数学,2019,35(02):116-121.

附录B

$n$阶矩阵一定有$n$个特征值(特征值可相同)
矩阵 $A$ 为实对称矩阵 $\Rightarrow$ $A$可相似对角化 $\iff$ $A$有$n$个线性无关的特征向量
矩阵 $A$ 可相似对角化，$\lambda$ 为$A$的$k$重特征值 $\iff$ $\lambda$ 对应有$k$个线性无关的特征向量
若矩阵 $A$ 为$n$阶实对称矩阵，$r=\text{rank}(A)$, 则$\lambda =0$为 $A$的$n-r$重特征值
若矩阵 $A$ 能相似对角化, $r=\text{rank}(A)$ ,则$\lambda =0$为 $A$ 的 $n-r$ 重特征值
若矩阵 $A$ 为$n$阶方阵,$r=\text{rank}(A)$,则$\lambda =0$ 至少为 $A$ 的 $n-r$ 重特征值.

[2]秦川,李小飞.方阵的秩与特征值的关系[J].课程教育研究：学法教法研究,2015(27):120-120.

全文翻译自

HOW ARE PRINCIPAL COMPONENT ANALYSIS AND SINGULAR VALUE DECOMPOSITION RELATED?