奇异值分解 SVD

我们在上一章中研究了矩阵的特征值以及各种求解算法，看到特征值和特征向量反映了矩阵本身的一些特性，但它在数学形式上存在一些限制。比如只能用于方阵，可能存在复特征值的情况。我们在工程中，为了估计某些系统状态，往往会采集大量的样本，构建超定方程。此时，方阵就很难满足我们的分析需求了。大多数时候，我们并不关心具体的特征值是多少，更关心的是它们绝对值的大小，及其对应的特征空间。

奇异值是更一般形式的数学工具，它不仅适用于任何矩阵，而且能够清晰表达矩阵的特征空间。它可以说是现代很多工程问题的数学基石，在数据分析、状态估计、统计学习等诸多方面都发挥着关键作用。

1. SVD 的定义和一些重要性质

任何一个 $m \times n$ 的实数矩阵 $\boldsymbol{A}$ 都可以分解成 $\boldsymbol{A} = \boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T$ 的形式。其中 $\boldsymbol{U}$ 和 $\boldsymbol{V}$ 分别是 $m \times m$ 和 $n \times n$ 的正交矩阵。$\boldsymbol{\Sigma}$ 则是一个对角线元素为 $\sigma_i$ 的 $m \times n$ 的对角矩阵，$\boldsymbol{\Sigma}$ 只有对角线上的 $p$ 个元素有值 $p = \text{min}(m,n)$ 其余全为 0。即， $$ \begin{equation}\label{func_svd} \boldsymbol{U}^T \boldsymbol{A}\boldsymbol{V} = \boldsymbol{\Sigma} = \text{diag}\left(\sigma_1, \cdots, \sigma_p\right), \qquad p = \text{min}(m, n), \qquad \sigma_1 ≥ \sigma_2 ≥ \cdots ≥ \sigma_p ≥ 0 \end{equation} $$ 上式中，$\sigma_i$ 被称为矩阵 $\boldsymbol{A}$ 的奇异值(singular value)。构成矩阵 $\boldsymbol{U}$ 的列向量 $\boldsymbol{u}_i$ 是矩阵 $\boldsymbol{A}$ 的左奇异向量(left singular vector)，矩阵 $\boldsymbol{V}$ 的各列向量 $\boldsymbol{v}_i$ 则是 $\boldsymbol{A}$ 的右奇异向量(right singular vector)。这就是著名的奇异值分解 Singular Value Decomposition，在现代的科学计算中，它几乎无处不在。

根据 SVD 分解的定义，可以直接写出奇异值与左右奇异向量的关系 $\boldsymbol{A}\boldsymbol{v}_i = \sigma_i \boldsymbol{u}_i, \boldsymbol{A}^T\boldsymbol{u}_i = \sigma_i \boldsymbol{v}_i$。我们也可以推出，$\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{v}_i = \sigma_i^2 \boldsymbol{v}_i$, $\boldsymbol{A}\boldsymbol{A}^T\boldsymbol{u}_i = \sigma_i^2 \boldsymbol{u}_i$。这个推论实际反映了奇异值与特征值之间的关系：矩阵 $\boldsymbol{A}$ 的奇异值及对应的左右奇异向量，分别是矩阵 $\boldsymbol{A}\boldsymbol{A}^T$ 和 $\boldsymbol{A}^T\boldsymbol{A}$ 的特征值的平方根和特征向量。此外，SVD 分解还可以写成 $\boldsymbol{A} = \sum_{i=1}^p \sigma_i \boldsymbol{u}_i \boldsymbol{v}_i^T$ 的形式。显然矩阵 $\boldsymbol{A}\boldsymbol{A}^T$ 和 $\boldsymbol{A}^T\boldsymbol{A}$ 是对称的，并且它们的特征值 $\sigma_i^2$ 一定是个实数。根据Schur 分解定理，这两个矩阵一定是可以对角化的。

SVD 可以准确的表达一个矩阵的秩，及其象空间(range)与核空间(nullspace)的正交基。这一特性使得 SVD 占据了现代科学计算中的 C 位。如果矩阵 $\boldsymbol{A}$ 有 $r$ 个奇异值，那么它的秩就是 $r$。并且前 $r$ 个左奇异向量所张成的空间就是它的象空间，后 $n - r$ 个右奇异向量张成了它的核空间，既:

$$ \begin{cases} \text{rank}(\boldsymbol{A}) & = & r \\ \text{range}(\boldsymbol{A})& = & \text{span} \left\{ \boldsymbol{u}_1, \cdots, \boldsymbol{u}_r \right\} \\ \text{null}(\boldsymbol{A}) & = & \text{span} \left\{ \boldsymbol{v}_{r+1}, \cdots, \boldsymbol{v}_{n} \right\} \end{cases} $$

2. Golub-Kahan-Reinsch 算法

Golub-Kahan-Reinsch 算法，简称 GKR 算法，是最早的实用 SVD 算法。它十分稳定并且效率也很高，是经典的基于二对角化的算法。它的大体思想就是，先对矩阵 $\boldsymbol{A}$ 二对角化 $\boldsymbol{A} = \boldsymbol{U}_1\boldsymbol{B}\boldsymbol{V}_1^T$，再对矩阵 $\boldsymbol{B}$ 通过隐式的 QR 迭代进行 SVD 分解 $\boldsymbol{B} = \boldsymbol{U}_2 \boldsymbol{\Sigma} \boldsymbol{V}_2^T$。整个过程与QR算法先初始化一个 Hessenberg 矩阵，再进行 QR 迭代，求得矩阵的所有特征值类似。

2.1. 二对角化直接求解

与 QR 迭代时构造 Hessenberg 矩阵不同，这里对矩阵 $\boldsymbol{A}$ 二对角化并不是一个相似变换的过程。由于矩阵 $\boldsymbol{A}$ 的奇异值实际上是对称矩阵 $\boldsymbol{A}^T\boldsymbol{A}$ 和 $\boldsymbol{A}\boldsymbol{A}^T$ 特征值的平方根。所以我们只要保证变换后的矩阵 $\boldsymbol{B}^T\boldsymbol{B}$ 和 $\boldsymbol{B}\boldsymbol{B}^T$ 分别与$\boldsymbol{A}^T\boldsymbol{A}$ 和 $\boldsymbol{A}\boldsymbol{A}^T$ 相似。二对角化 $\boldsymbol{A} = \boldsymbol{U}_1\boldsymbol{B}\boldsymbol{V}_1^T$ 中的 $\boldsymbol{U}_1$ 和 $\boldsymbol{V}_1$ 都是正交矩阵，有:

$$ \boldsymbol{A}^T\boldsymbol{A} = (\boldsymbol{U}_1\boldsymbol{B}\boldsymbol{V}_1^T)^T\boldsymbol{U}_1\boldsymbol{B}\boldsymbol{V}_1^T = \boldsymbol{V}_1\boldsymbol{B}^T (\boldsymbol{U}_1^T \boldsymbol{U}_1) \boldsymbol{B}\boldsymbol{V}_1^T = \boldsymbol{V}_1\boldsymbol{B}^T\boldsymbol{B}\boldsymbol{V}_1^T $$ $$ \boldsymbol{A}\boldsymbol{A}^T = \boldsymbol{U}_1\boldsymbol{B}\boldsymbol{V}_1^T(\boldsymbol{U}_1\boldsymbol{B}\boldsymbol{V}_1^T)^T = \boldsymbol{U}_1\boldsymbol{B}(\boldsymbol{V}_1^T\boldsymbol{V}_1)\boldsymbol{B}^T\boldsymbol{U}_1^T = \boldsymbol{U}_1\boldsymbol{B}\boldsymbol{B}^T\boldsymbol{U}_1^T $$

任何矩阵都可以进行 SVD 分解，根据矩阵的形状不同，二对角化的方向也略有差异。对于一个 $m \times n$ 的矩阵，如果 $m >= n$ 那么我们就对其上二对角化，即，只有对角线和上次对角线有值其余全为 0。若 $m < n$ 就下二对角化，只保留对角线和下次对角线上的元素。如下面所示。

$$ \left. \begin{bmatrix} * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ \end{bmatrix} \xrightarrow{上二对角化} \begin{bmatrix} * & * & & \\ & * & * & \\ & & * & * \\ & & & * \\ & & & \\ & & & \\ \end{bmatrix} \qquad \right| \qquad \begin{bmatrix} * & * & * & * & * & * \\ * & * & * & * & * & * \\ * & * & * & * & * & * \\ * & * & * & * & * & * \\ \end{bmatrix} \xrightarrow{下二对角化} \begin{bmatrix} * & & & & & \\ * & * & & & & \\ & * & * & & & \\ & & * & * & & \\ \end{bmatrix} $$

由于二对角化之后的矩阵 $\boldsymbol{B}$ 只有对角线和次对角线有值，所以 $\boldsymbol{B}^T\boldsymbol{B}$ 是一个 $n \times n$ 的三对角对称方阵 $\begin{bmatrix} * & * & & \\ * & * & * & \\ & * & * & * \\ & & * & * \\ \end{bmatrix}$。与 QR 分解类似，二对角化可以通过不断构造 Householder 矩阵来实现。不同之处是需要构造右乘的 Householder 矩阵，用于消除每行中对角或者次对角右侧的元素。假设 $m >= n$，首先构造矩阵 $\boldsymbol{H}_1 \in \mathbb{R}^{m\times m}$，消除第一列第一个行以下的元素。在构造矩阵 $\begin{bmatrix} 1 & \\ & \tilde{\boldsymbol{H}}_1 \end{bmatrix} \in \mathbb{R}^{n \times n}$，消除第一行第二列右侧的元素。

$$ \boldsymbol{A} = \begin{bmatrix} * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ \end{bmatrix} \longrightarrow \boldsymbol{H}_1 \boldsymbol{A} = \begin{bmatrix} * & * & * & * \\ & * & * & * \\ & * & * & * \\ & * & * & * \\ & * & * & * \\ & * & * & * \\ \end{bmatrix} \longrightarrow \boldsymbol{H}_1 \boldsymbol{A}\begin{bmatrix} 1 & \\ & \tilde{\boldsymbol{H}}_1 \end{bmatrix} = \begin{bmatrix} * & * & & \\ & * & * & * \\ & * & * & * \\ & * & * & * \\ & * & * & * \\ & * & * & * \\ \end{bmatrix} $$

如此重复下去就可以完成上二对角化。当 $m < n$ 时，我们可以先构造右乘的 Householder 矩阵消除右侧元素，再构造左乘的矩阵消除下方元素，实现下二对角化。完整的实现可以参考例程Bidiagonal.hpp。接下来我们需要求解矩阵 $\boldsymbol{B}^T\boldsymbol{B}$ 的所有特征值和特征向量，这可以通过我们之前研究的隐式QR迭代算法完成。由于矩阵 $\boldsymbol{B}^T\boldsymbol{B}$ 是对称的，所以它是可对角化的，对其进行 QR 迭代最后一定收敛到特征值分解的形式， $$ \begin{equation}\label{btb} \boldsymbol{B}^T\boldsymbol{B} = \boldsymbol{Q}\boldsymbol{\Lambda}\boldsymbol{Q}^T, \qquad \boldsymbol{\Lambda} = \text{diag}(\sigma_1^2, \sigma_2^2, \cdots, \sigma_n^2) \end{equation} $$ 其中，$\boldsymbol{\Lambda}$ 是个对角矩阵，$\boldsymbol{Q}$ 是正交矩阵。根据上式，有 $(\boldsymbol{B}\boldsymbol{Q})^T(\boldsymbol{B}\boldsymbol{Q}) = \boldsymbol{\Lambda}$。所以矩阵 $\boldsymbol{B}\boldsymbol{Q}$ 是列正交的(不是单位正交)。如果我们对矩阵 $\boldsymbol{B}\boldsymbol{Q}$ 进行 QR 分解，则有 $\boldsymbol{B}\boldsymbol{Q} = \boldsymbol{U}\boldsymbol{R}$，其中$\boldsymbol{U}$是个正交矩阵，$\boldsymbol{R}$是上三角矩阵。因为 $\boldsymbol{B}\boldsymbol{Q}$ 列正交，$\boldsymbol{U}$ 正交，所以 $\boldsymbol{R} = \boldsymbol{U}^T\boldsymbol{B}\boldsymbol{Q}$ 也是列正交的，而$\boldsymbol{R}$ 还是一个上三角阵，所以它一定是对角的。所以 $$ \boldsymbol{U}^T\boldsymbol{B}\boldsymbol{Q} = \boldsymbol{R} $$ 就是矩阵$\boldsymbol{B}$ 的 SVD 分解。矩阵 $\boldsymbol{R}$ 中对角线上的元素就是要求的奇异值。将上式带入二对角化，有: $$ \boldsymbol{U}^T\boldsymbol{U}_1^T\boldsymbol{A}\boldsymbol{V}_1\boldsymbol{Q} = \boldsymbol{R} $$ 那么 $\boldsymbol{U}_1\boldsymbol{U}$ 就是 $\boldsymbol{A}$ 的左奇异向量矩阵，$\boldsymbol{V}_1\boldsymbol{Q}$ 则是其右奇异向量矩阵。详细的求解过程参见SVD_Naive.hpp

2.3. 隐式迭代

但是直接通过转置相乘得到矩阵 $\boldsymbol{B}^T\boldsymbol{B}$ 在数值上是不稳定的，因为它相当于一个平方运算，当奇异值接近 0 时会损失计算精度。GKR 算法直接在矩阵$\boldsymbol{B}$左右两侧不断地进行 Givens 变换完成隐式迭代，得益于隐式 Q 定理，该算法不需要直接计算$\boldsymbol{B}^T\boldsymbol{B}$。

$$ \text{上二对角阵} \begin{bmatrix} d_1 & & & \\ f_1 & d_2 & & \\ & f_2 & d_3 & \\ & & f_3 & d_4 \\ \end{bmatrix} \begin{bmatrix} d_1 & f_1 & & \\ & d_2 & f_2 & \\ & & d_3 & f_3 \\ & & & d_4 \\ \end{bmatrix} = \begin{bmatrix} d_1 d_1 & d_1 f_1 & & \\ f_1 d_1 & f_1^2+d_2^2 & d_2 f_2 & \\ & f_2 d_2 & f_2^2+d_3^2 & d_3 f_3 \\ & & f_3 d_3 & f_3^2+d_4^2 \\ \end{bmatrix} $$

假设 $\boldsymbol{A}$ 是一个 $m \times n$ 的矩阵，$m >= n$，我们对它进行的是上二对角化。矩阵 $\boldsymbol{B}$ 和 $\boldsymbol{B}^T\boldsymbol{B}$ 的形式如上式所示。我们可以轻松写出矩阵 $\boldsymbol{B}^T\boldsymbol{B}$ 的右下角$2\times 2$ 的子阵 $\boldsymbol{T}$

$$ \boldsymbol{T} = \begin{bmatrix} d_{n-1}^2 + f_{n-2}^2 & d_{n-1} f_{n_1} \\ f_{n-1} d_{n-1} & d_{n}^2 + f_{n-1}^2 \end{bmatrix} = \begin{bmatrix} t_{11} & t_{12} \\ t_{21} & t_{22} \end{bmatrix} $$

对于带位移的 QR 迭代而言，我们可以选择右下角的元素 $t_{22}$作为 Rayleigh 偏移，或者直接根据一元二次方程求解公式写出 $\boldsymbol{T}$ 的两个特征值，选择更接近 $t_{22}$ 的那个作为 Wilkinson 偏移。Golub 和 Reinsch 认为 Wilkinson 偏移更好。选择了偏移量 $\mu$ 之后，就可以根据 $\boldsymbol{B}^T\boldsymbol{B}$ 的第一列构造第一个 Givens 变换，$G(1,2,\theta)$，满足如下关系:

$$ \underbrace{\begin{bmatrix} \cos {\theta} & \sin {\theta} \\ -\sin {\theta} & \cos {\theta} \end{bmatrix}}_{G(1,2,\theta)} \begin{bmatrix} d_1^2 - \mu \\ f_1 d_1 \end{bmatrix} = \begin{bmatrix} * \\ 0 \end{bmatrix} $$

显然不需要完整计算 $\boldsymbol{B}^T\boldsymbol{B}$ 我们就可以写出偏移量 $\mu$ 和 $G(1,2,\theta)$。根据隐式 Q 定理，在QR分解过程中一旦我们确定了矩阵Q的第一列，Q矩阵就确定了。当我们把 $G(1,2\theta)$ 右乘到矩阵 $\boldsymbol{B}$ 上时，就会打破它的二对角矩阵结构。GKR 算法接下来要做的就是用 Givens 变换，整理矩阵$\boldsymbol{B}$，恢复它的二对角结构。比如，我们对一个 $4 \times 4$ 的上二对角阵右乘 $G(1,2,\theta)$，就会得到如下的结构:

$$ \begin{bmatrix} * & * & & \\ & * & * & \\ & & * & * \\ & & & * \\ \end{bmatrix} \xrightarrow{\boldsymbol{B} \leftarrow \boldsymbol{B}G(1,2,\theta)} \begin{bmatrix} * & * & & \\ + & * & * & \\ & & * & * \\ & & & * \\ \end{bmatrix} $$

接下来的动作与隐式QR迭代算法类似，通过一系列左乘和右乘 Givens 变换不断地将新增的“气泡”从矩阵中赶出去。

$$ \begin{bmatrix} * & * & & \\ + & * & * & \\ & & * & * \\ & & & * \\ \end{bmatrix} \xrightarrow{\boldsymbol{B} \leftarrow G(1,2,\theta)\boldsymbol{B}} \begin{bmatrix} * & * & + & \\ 0 & * & * & \\ & & * & * \\ & & & * \\ \end{bmatrix} \xrightarrow{\boldsymbol{B} \leftarrow \boldsymbol{B}G(2,3,\theta)} \begin{bmatrix} * & * & 0 & \\ & * & * & \\ & + & * & * \\ & & & * \\ \end{bmatrix} \xrightarrow{\boldsymbol{B} \leftarrow G(2,3,\theta)\boldsymbol{B}} \begin{bmatrix} * & * & & \\ & * & * & + \\ & 0 & * & * \\ & & & * \\ \end{bmatrix} \xrightarrow{\boldsymbol{B} \leftarrow \boldsymbol{B}G(3,4,\theta)} \begin{bmatrix} * & * & & \\ & * & * & 0 \\ & & * & * \\ & & + & * \\ \end{bmatrix} $$

2.3. 对角分块

在迭代的过程中，一定会遇到次对角线上的元素收敛到 0 的情况。此时 $\boldsymbol{B}$ 具有如下的结构。求解矩阵 $\boldsymbol{B}$ 的奇异值，就可以拆分为求解两个规模更小的矩阵 $\boldsymbol{B}_{11}, \boldsymbol{B}_{22}$ 的奇异值。

$$ \boldsymbol{B} = \begin{bmatrix} \boldsymbol{B}_{11} & \\ & \boldsymbol{B}_{22} \end{bmatrix} $$

有时，如果矩阵$\boldsymbol{B}$不满秩，我们也会遇到对角线上的元素近似为 0 的情况。为了完整体现消解过程，我们以一个 $5 \times 5$ 的矩阵为例，假设其中(3,3)的元素为 0。此时，我们需要连续左乘一系列的 Givens 变换，在行空间上将第三行完全消解为 0。

$$ \begin{bmatrix} * & * & & & \\ & * & * & & \\ & & 0 & * & \\ & & & * & * \\ & & & & * \\ \end{bmatrix} \xrightarrow{\boldsymbol{B} \leftarrow G(4,3,\theta)\boldsymbol{B}} \begin{bmatrix} * & * & & & \\ & * & * & & \\ & & & 0 & + \\ & & & * & * \\ & & & & * \\ \end{bmatrix} \xrightarrow{\boldsymbol{B} \leftarrow G(5,3,\theta)\boldsymbol{B}} \begin{bmatrix} * & * & & & \\ & * & * & & \\ & & & & 0 \\ & & & * & * \\ & & & & * \\ \end{bmatrix} \longrightarrow \left[ \begin{array}{c|c} \begin{matrix} * & * & \\ & * & * \\ & & \\ \end{matrix} & \begin{matrix} & \\ & \\ \end{matrix} \\ \hline \begin{matrix} \end{matrix} & \begin{matrix} * & * \\ & * \\ \end{matrix} \end{array} \right] $$

仍然按照次对角线的元素是否为 0，对矩阵$\boldsymbol{B}$进行分割，就会得到 $\boldsymbol{B}_{11} = \begin{bmatrix} * & * & \\ & * & * \\ & & \end{bmatrix}$， $\boldsymbol{B}_{22} = \begin{bmatrix} * & * \\ & * \end{bmatrix}$。这里的 $\boldsymbol{B}_{11}$ 还需要再处理一下，对它在列空间上将第三列消解为 0。这样 $\boldsymbol{B}_{11}$ 实质上就是一个 $2 \times 2$ 的二对角矩阵。

$$ \begin{bmatrix} * & * & \\ & * & * \\ & & \end{bmatrix} \xrightarrow{\boldsymbol{B} \leftarrow \boldsymbol{B}G(2,3,\theta)} \begin{bmatrix} * & * & + \\ & * & 0 \\ & & 0 \\ \end{bmatrix} \xrightarrow{\boldsymbol{B} \leftarrow \boldsymbol{B}G(1,3,\theta)} \begin{bmatrix} * & * & 0 \\ & * & \\ & & \\ \end{bmatrix} $$

当把矩阵$\boldsymbol{B}$完全分割成 $1\times 1$ 的对角块，就完成了SVD的迭代。这些对角块中的元素就是 $\boldsymbol{B}$ 的奇异值。将迭代过程中左乘和右乘的 Givens 变换收集起来，就得到了奇异值分解 $\boldsymbol{U}_2^T\boldsymbol{B}\boldsymbol{V}_2 = \boldsymbol{\Sigma}$。

对于 $m < n$ 的情形，数学上只需要转置一下皆可以直接应用上述的讨论了。工程实现上，考虑到存储空间和额外的拷贝动作，我在 SVD_GKR.hpp 中，根据矩阵的形状区别对待 $m ≥ n$ 和 $m < n$ 两种情况。