GAMES001课程笔记04-奇异值分解与主成分分析

这个系列是GAMES001-图形学中的数学(GAMES 001: Mathematics in Computer Graphics)的同步课程笔记。课程旨在总结归纳图形学学习过程中重要的数学概念、理论和方法，并将以索引的形式在每一章节中穿插讲解该数学专题在图形学中的应用。本课程既可以作为GAMES系列其它课程学习的基础或「手册」，也可以作为站在不一样的视角复习图形学知识的平台。本节主要介绍奇异值分解与主成分分析的相关内容。

奇异值分解

矩阵有很多不同的理解方式。在前面的课程中我们提到过矩阵可以理解为一种线性变换，除此之外矩阵还可以表示一张图片或是点云这种几何数据。因此我们可以从不同的视角来理解矩阵变换的各种意义。

奇异值分解(singular value decomposition, SVD)是一种对矩阵进行分析和处理的重要计算工具，它在数学和整个计算机科学领域都有着大量的应用。

实对称矩阵

我们从实对称矩阵出发。根据线性代数相关理论，实对称矩阵具有如下重要性质：

实对称矩阵的特征值都是实数。
实对称矩阵不同特征值的特征向量相互正交。

在此基础上我们可以得到实对称矩阵的特征值分解：

\[S = Q \Lambda Q^T = \begin{pmatrix} q_0 & \dots & q_{n-1} \end{pmatrix} \begin{pmatrix} \lambda_0 & & \\ & \ddots & \\ & & \lambda_{n-1} \end{pmatrix} \begin{pmatrix} q_0^T \\ \vdots \\ q_{n-1}^T \end{pmatrix}\]

其中，特征值\(\lambda_i\)有正有负，且非零特征值的个数等于矩阵的秩\(r = \text{rank} (S)\)。

如果\(S\)的所有特征值都是非负的，则称它是一个半正定(semi-positive definite, SPD)矩阵。对于半正定矩阵，我们可以定义实对角矩阵\(\Sigma\)：

\[\Sigma = \begin{pmatrix} \sqrt{\lambda_0} & & \\ & \ddots & \\ & & \sqrt{\lambda_{n-1}} \end{pmatrix}\]

显然它满足\(\Sigma^2 = \Lambda\)。这样就可以定义实对称矩阵的平方根\(P\)：

\[S = Q \Lambda Q^T = Q \Sigma^2 Q^T = (Q \Sigma Q^T)(Q \Sigma Q^T) = P^2\] \[P = S^{\frac{1}{2}}\]

此外，容易验证实对称矩阵的迹等于其特征值之和：

\[\text{Tr} (S) = \text{Tr} (\Lambda) = \sum_i \lambda_i\]

奇异值分解

对于一个更一般的矩阵\(A \in \mathbb{R}^{m \times n}\)，其秩为\(r\)，我们定义一个新的矩阵\(S\)：

\[S = A^T A\]

显然\(S\)是一个实对称半正定矩阵。对\(S\)进行特征值分解，可以得到\(r\)个特征向量\(\{ v_i \}\)构成一组正交基：

\[S v_i = \lambda_i v_i\]

两边同时乘以特征向量\(v_j\)可以得到：

\[v_i^T S v_j = v_i^T A^T A v_j = \lambda_j v_i^T v_j = \lambda_j \delta_{ij}\]

接下来定义向量\(u_i = \frac{1}{\sqrt{\lambda_i}} A v_i\)，容易验证\(\{ u_i \}\)构成一组正交基：

\[\begin{aligned} u_i^T u_j &= \bigg( \frac{1}{\sqrt{\lambda_i}} A v_i \bigg)^T \bigg( \frac{1}{\sqrt{\lambda_j}} A v_j \bigg) \\ &= \frac{1}{\sqrt{\lambda_i \lambda_j}} v_i^T A^T A v_j \\ &= \delta_{ij} \end{aligned}\]

通过Schmidt正交化可以分别对\(\{ v_i \}\)和\(\{ u_i \}\)进行扩充，从而得到\(\mathbb{R}^n\)和\(\mathbb{R}^m\)上的正交基。其中扩充的基满足：

\[\sqrt{\lambda_i} u_j = A v_j = 0\]

这样两组基的关系写成矩阵形式有：

\[A V = U \Sigma\]

其中，\(V\)和\(U\)都是正交矩阵。

上面的推导即为矩阵的奇异值分解，其更常见的形式为：

\[A = U \Sigma V^T\]

其中，矩阵\(A\)为一个任意m×n维矩阵(不需要是方阵，也不需要满秩)；\(U\)和\(V\)分别是\(\mathbb{R}^{m \times m}\)和\(\mathbb{R}^{n \times n}\)的正交矩阵；\(\Sigma\)是\(\mathbb{R}^{m \times n}\)的「对角」矩阵，其对角元均是非负数称为矩阵\(A\)的奇异值。

矩阵的奇异值分解是线性代数中的重要概念。当矩阵的奇异值分解已知时，我们可以快速地求解矩阵相关的各种问题。

矩阵的逆与SVD

对于一个可逆矩阵\(A \in \mathbb{R}^{n \times n}\)，我们可以通过SVD来直接获得它的逆：

\[A^{-1} = V \Sigma^{-1} U^T\]

列满秩

对于非方阵的情况，我们可以推广逆阵的概念从而得到矩阵的伪逆。假设矩阵\(A\)是列满秩矩阵而且它的行数大于列数\(m \gt n\)，我们的目标是求解矩阵方程：

\[A x = b\]

由于方程数大于未知量的个数，它的解可能是不存在的。此时我们可以转而寻找它在最小二乘意义下的解，这样构造最小二乘的优化目标：

\[\min_x \Vert Ax - b \Vert^2\]

利用SVD可以对优化目标进行变形：

\[\begin{aligned} \Vert Ax - b \Vert^2 &= \Vert U \Sigma V^T x - b \Vert^2 \\ &= \Vert U (\Sigma V^T x - U^T b) \Vert^2 \\ &= \Vert \Sigma V^T x - U^T b \Vert^2 \end{aligned}\]

然后通过换元\(x' = V^T x\)，\(b' = U^T b\)得到新的优化目标：

\[\Vert \Sigma V^T x - U^T b \Vert^2 = \Vert \Sigma x' - b' \Vert^2\]

由于\(\Sigma\)是一个包含\(n\)个非零奇异值的对角矩阵，\(\Vert \Sigma x' - b' \Vert^2\)取最小值时需要满足以下条件：

\[\begin{cases} \begin{aligned} \sigma_0 x_0' &= b'_0 \\ \vdots \\ \sigma_{n-1} x_{n-1}' &= b'_{n-1} \end{aligned} \end{cases}\]

也即

\[x' = \hat{\Sigma}^{-1} b'\]

其中\(\hat{\Sigma}^{-1}\)是一个n×m维对角矩阵，对角元为相应奇异值的倒数。最后解出未知量\(x\)的表达式为：

\[x = V \hat{\Sigma}^{-1} U^T b\]

其中的系数矩阵\(V \hat{\Sigma}^{-1} U^T\)就称为矩阵\(A\)的伪逆(pseudo inverse)。

除了上面的方法，我们还可以直接对最小二乘的优化目标进行求导来得到伪逆。此时可以得到矩阵方程：

\[A^T A x = A^T b\]

假设\(A^T A\)是满秩的，我们可以得到\(x\)的表达式：

\[x = (A^T A)^{-1} A^T b\]

其中的系数矩阵\((A^T A)^{-1} A^T\)称为矩阵\(A\)的左逆(left inverse)。容易验证当矩阵\(A\)是列满秩的时候，左逆和伪逆是等价的。

行满秩

另一种可能出现的情况是矩阵\(A\)是行满秩，且\(m \lt n\)。此时方程组的个数小于未知量的个数，方程组有无穷多的解。在这种情况下我们可以构造一个新的最小二乘优化目标：

\[\begin{aligned} \min_x &\Vert x \Vert^2 \\ \text{s.t.} \ &Ax = b \end{aligned}\]

上式的意义在于寻找满足约束\(Ax = b\)条件下模长最小的一组未知量。

通过SVD分解，可以得到新的优化目标的解同样是伪逆：

\[x = V \hat{\Sigma}^{-1} U^T b\]

当然，我们同样可以通过优化的方法来求解这个欠定方程组。此时可以得到矩阵\(A\)的右逆(right inverse)：

\[A^T (A^T A)^{-1}\]

简单总结一下：

当矩阵\(A\)是列满秩时，其左逆\((A^T A)^{-1} A^T\)存在，可以用来求解超定方程。
当矩阵\(A\)是行满秩时，其右逆\(A^T (A^T A)^{-1}\)存在，可以用来求解欠定方程。
任何情况下，矩阵的伪逆\(V \hat{\Sigma}^{-1} U^T\)都一定存在。

几何变换

从几何变换的角度来看，SVD分解说明一个线性变换可以分解为两次旋转加上一次拉伸缩放。

除了SVD外，任意实方阵\(A\)还可以通过极分解(polar decomposition)来得到一个正交阵和半正定矩阵：

\[A = R S\]

极分解可以通过SVD来显式构造得到：

\[R = U V^T\] \[S = V \Sigma V^T\]

形状匹配

SVD在图形学中的一个经典应用是形状匹配(shape matching)，它的目标是寻找一个刚体变换使得形状\(q\)能够变形到形状\(p\)上。这个问题等价于如下优化问题：

\[R, t = \arg \min_i \sum_i \Vert p_i - (R q_i + t) \Vert^2\]

这个优化问题可以利用SVD来进行求解。首先假设旋转\(R\)是已知的从而只考虑位移\(t\)，通过对优化目标求偏导可以得到：

\[\begin{aligned} \frac{\partial}{\partial} \sum_i \Vert p_i - (R q_i + t) \Vert^2 = \sum_i R q_i + t - p_i = 0 \end{aligned}\] \[t = \frac{1}{n} \bigg( \sum_i p_i - R \sum_i q_i \bigg)\]

取质点集合变形前后的质心分别为\(q_c = \frac{1}{n} q_i\)以及\(p_c = \frac{1}{n} p_i\)，则有：

\[t = p_c - R q_c\]

接下来考虑求解旋转\(R\)。把位移的表达式带入优化目标中可以得到：

\[\sum_i \Vert p_i - (R q_i + t) \Vert^2 = \sum_i \Vert (p_i - p_c) - R (q_i - q_c) \Vert^2\]

记中心化后的坐标分别为\(p_i' = p_i - p_c\)以及\(q_i' = q_i - q_c\)，这样优化目标可以转换为：

\[R = \arg \min \sum_i \Vert p_i' - R q_i' \Vert^2\]

将上式展开并且只考虑关于旋转\(R\)的部分可以得到新的最优化问题：

\[\min \sum_i \Vert p_i' - R q_i' \Vert^2 \Rightarrow \max \sum_i p_i'^T R q_i'\]

利用矩阵迹的性质继续进行变形：

\[\sum_i p_i'^T R q_i' = \text{Tr} \bigg( \sum_i p_i'^T R q_i' \bigg) = \text{Tr} \bigg( R \sum_i q_i' p_i'^T \bigg)\]

记求和号内的矩阵为\(H = \sum_i q_i' p_i'^T\)，得到最终的优化目标：

\[\max_R \text{Tr} (R H)\]

此时旋转矩阵\(R\)可由\(H\)矩阵的SVD给出：

\[R = V U^T\]

总结一下，形状匹配算法的流程如下：

计算两组质点的质心并分别进行中心化\(q_i' = q_i - q_c\)、\(p_i' = p_i - p_c\)。
构造矩阵\(H = \sum_i q_i' p_i'^T\)，并对其进行分解\(H = U \Sigma V^T\)。
计算最优旋转\(R = V U^T\)。
计算相对位移\(t = p_c - R q_c\)。

实际上上述形状匹配算法即为点云匹配中经典的ICP算法。

图像压缩

SVD的一个重要应用是对图像进行压缩。我们可以把一幅图像看作是一个大型稠密矩阵，根据SVD的形式把它写成一系列秩为1的矩阵之和：

\[A = U \Sigma V^T = \sum_i \sigma_i u_i v_i^T\]

对于通常的图像，其奇异值的大小往往只几种在前若干阶上。随着奇异值数量的增加，奇异值的大小一般会显著减少。因此，我们可以只考虑那些奇异值比较大的特征向量来重建原始的图像。这样可以显著减少存储所需要的空间。

实际上前面叙述的对矩阵进行低秩重建的思想在如今火热的大模型领域也有着重要的应用。

SVD求导

除此之外，SVD实际上还是可导的运算。在很多应用中涉及到对SVD进行求导。

然而遗憾的是，SVD的求导一般是比较难处理的。

一般来说，SVD的复杂度还是比较高的。好在目前绝大多数的数值计算库中都有比较好的实现，在实际工程中一般不会需要从头进行编程。

主成分分析

主成分分析(principal component analysis, PCA)是SVD的一个经典应用。PCA的一大作用是对高维的数据进行降维，使得我们可以在相对低维的空间来分析它们。

从几何的角度来看，PCA的目标是寻找到数据分布的主方向\(v \in \mathbb{R}^n\)，使得我们把数据投影到主方向上后有着最大的分布。对于样本\(x_i\)，它在直线\(x = m + t v\)上的投影为：

\[x_i' = (x_i - m) \cdot v\]

其中\(m = \frac{1}{n} \sum_i x_i\)为数据的中心，也是样本数据的零阶近似。此时，数据点在垂直直线方向上的距离为：

\[d_i^2 = \Vert x_i - m \Vert^2 - \Vert x_i' \Vert^2\]

最小化\(d_i\)等价于最大化\(\Vert x_i' \Vert^2\)。这样，我们可以定义一个最优化问题：

\[\max v^T \bigg( \sum_i (x_i - m) (x_i - m)^T \bigg) v\]

其中，\(S = \sum_i (x_i - m) (x_i - m)^T\)即为数据分布的方差矩阵。

显然，方差矩阵\(S\)是一个半正定矩阵，其全部特征值非负。对它进行特征值分解有：

其中，\(Q^T Q = I\)，\(\lambda_0 \geq \lambda_1 \geq \cdots \geq 0\)为\(S\)的全部特征值。把\(S\)的特征值分解代入优化目标中，可以得到：

\[v^T S v = v^T Q \Lambda Q^T v = (Q^T v)^T \Lambda (Q^T v)\]

显然，当\(v = q_0\)时上式取到最大值\(\lambda_0\)。

找到主方向并完成投影后，我们可以对剩下的分量继续寻找主方向。实际上这些分量的主方向对应着方差矩阵\(S\)的第二个特征向量，因此基于特征值分解我们可以得到原始样本数据的全部主方向。

总结一下，\(\mathbb{R}^n\)中样本数据组成的数据集其前\(m\)个主方向为方差矩阵\(S\)的前\(m\)个特征向量。这些向量相互正交，张成了\(m\)维的子空间。将数据投影到这个子空间中，就构成了原始数据的低维表达。

应用

我们可以利用PCA来计算点云的法向。此时对于点云中的每个点我们需要提取它和它附近的一些样本点组成数据集，然后在这个数据集上计算三个主方向就能够得到该点的法向。

在人脸识别领域，PCA的一个经典案例是eigenface。我们可以对人脸图像数据使用PCA进行降维，然后提取特征值大的那些特征向量组成人脸数据的特征脸。

在数值模拟中，我们也可以利用PCA来对模型进行降阶从而降低仿真的计算量。

最后需要说明的是PCA本质上是一种线性降维方法，对于非线性的数据PCA可能不会得到特别好的效果。

要处理高度非线性的数据可以考虑kernel PCA或是自编码器这样的方法。

Reference

Lecture 05: 奇异值分解与主成分分析