5.6 奇异值分解

奇异值分解 (SVD) 具有十分广泛的应用背景, 因此, 如何更好更快地计算一个给定矩阵的 SVD 是科学与工程计算领域中的一个热门研究课题, 吸引了众多专家进行这方面的研究, 也涌现出了许多奇妙的方法. 本章主要介绍计算 SVD 的常用算法.

对任意矩阵 $A \in \mathbb{R}^{m \times n}$ , 其奇异值与对称矩阵 $A^{\mathrm{T}}A, AA^{\mathrm{T}}$ 和 $\left[ \begin{array}{cc}0 & A^{\mathrm{T}} \\ A & 0 \end{array} \right]$ 的特征值是是密切相关的, 故理论上计算对称特征值的算法都可以用于计算奇异值. 但在实际计算中, 我们需要对算法进行优化, 使得计算过程更加高效, 计算结果更加准确.

当前的SVD算主要可分为两类[33]：一类是基于二对角化的算法，包括Golub-Kahan SVD算法(即QR算法)、DC算法、dqds算法、对分法、MRRR算法等；另一类是Jacobi算法，包括双边Jacobi算法和单边Jacobi算法.

关于SVD的相关参考资料

J. Dongarra, et al., The Singular Value Decomposition: Anatomy of Optimizing an Algorithm for Extreme Scale, 2018. [33]

5.6.1 二对角化

基于二对角化的SVD算法通常分为以下三个步骤：

1. 将 $A$ 二对角化: $B = U_{1}^{\top} A V_{1}$ , 其中 $B$ 为上二对角矩阵, $U_{1}, V_{1}$ 为正交阵;

2. 计算 $B$ 的SVD: $B = U_{2}\Sigma V_{2}^{\mathsf{T}}$ , 其中 $\Sigma$ 为对角阵, $U_{2}, V_{2}$ 为正交阵;

3. 合并得到 $A$ 的SVD: $A = U_{1}BV_{1}^{\mathsf{T}} = (U_{1}U_{2})\Sigma (V_{1}V_{2})^{\mathsf{T}}$ .

我们知道, 对称矩阵可以通过一系列 Householder 变换转化为对称三对角矩阵. 对于一般矩阵 $A \in \mathbb{R}^{m \times n}$ , 我们也可以通过 Householder 变换, 将其转化为二对角矩阵, 即计算正交矩阵 $U_{1}$ 和 $V_{1}$ 使得

其中 $B$ 是一个实 (上) 二对角矩阵. 这个过程就称为二对角化

需要注意的是, 与对称矩阵的对称三对角化不同, $A$ 与 $B$ 是不相似的.

设 $A\in \mathbb{R}^{m\times n}$ ，二对角化过程大致如下：

(1) 首先构造一个 Householder 矩阵 $H_{1} \in \mathbb{R}^{m \times m}$ , 使得 $H_{1}A$ 的第一列除第一个元素外, 其它分量都为零, 即

(2) 再构造一个 Householder 矩阵 $\tilde{H}_1 \in \mathbb{R}^{(n-1) \times (n-1)}$ , 把 $H_1A$ 的第一行的第 3 至第 $n$ 个元素化为零, 即

(3) 重复上面的过程, 直到把 $A$ 最终化为二对角矩阵.

有了分解 (5.6) 以后, 我们可得

即 $V_{1}^{\mathsf{T}}A^{\mathsf{T}}AV_{1} = B^{\mathsf{T}}B.$ 由于 $B^{\mathsf{T}}B$ 是对称三对角的, 所以这就相当于将 $A^{\mathsf{T}}A$ 三对角化.

整个二对角化过程的运算量约为 $4mn^{2} + 4m^{2}n - 4n^{3} / 3$ , 若不需要计算 $U_{1}$ 和 $V_{1}$ , 则运算量约为 $4mn^{2} - 4n^{3} / 3$ .

二对角矩阵的奇异值分解

设 $B \in \mathbb{R}^{n \times n}$ 是一个二对角矩阵

我们假定二对角矩阵的元素都非负, 即 $a_{k}, b_{k} \geq 0$ , 否则可以通过在左右分别乘以一个正交矩阵来实现 (见习题 5.8). 下面三种方法均可将计算 $B$ 的 SVD 转化成计算对称三对角矩阵的特征分解:

(1) 构造增广矩阵 $T = \left[ \begin{array}{cc}0 & B^{\mathsf{T}}\\ B & 0 \end{array} \right]$ , 取置换矩阵 $P = [e_1,e_{n + 1},e_2,e_{n + 2},\ldots ,e_n,e_{2n}]$ , 则 $T_{ps} = P^{\mathsf{T}}TP$ 是对称三对角矩阵, 且主对角线元素全为 0, 次对角线元素为 $a_1,b_1,a_2,b_2,\ldots ,a_{n - 1}$ , $b_{n - 1},a_{n}$ . 设 $(\lambda_i,x_i)$ 是 $T_{ps}$ 的一个特征对, 则

其中 $\sigma_{i}$ 为 $B$ 一个奇异值, $u_{i}$ 和 $v_{i}$ 分别为对应的左和右奇异向量.

(2) 令 $T_{BB^{\intercal}} = BB^{\intercal}$ , 则

$T_{BB^{\intercal}}$ 的特征值为 $B$ 的奇异值的平方, 且 $T_{BB^{\intercal}}$ 的特征向量为 $B$ 的左奇异向量.

(3) 令 $T_{B^{\intercal}B} = B^{\intercal}B$ , 则

$T_{B^{\intercal}B}$ 的特征值为 $B$ 的奇异值的平方, 且 $T_{B^{\intercal}B}$ 的特征向量为 $B$ 的右奇异向量.

Golub & Kahan (1965) 采用前面两种做法. Reinsch (1970) 采用第三种做法, 其巧妙之处是,只需对 $B$ 做 Givens 变换即可, 类似隐式 QR 迭代.

理论上, 我们可以直接使用 QR 迭代、分而治之法或带反迭代的对分法, 计算对称三对角矩阵 $T_{ps}, T_{BB^{\intercal}}$ 和 $T_{B^{\intercal}B}$ 的特征值和特征向量. 但一般来说, 这种做法并不是最佳的, 原因如下:

(1) 对 $T_{ps}$ 做 QR 迭代并不划算, 因为 $QR$ 迭代计算所有的特征值和特征向量, 而事实上只要计算正的特征值即可;
(2) 直接构成 $T_{BB^{\top}}$ 或 $T_{B^{\top}B}$ 是数值不稳定的. 事实上, 这样做可能会使得 $B$ 的小奇异值的精度丢失一半.

下面是一些计算奇异值分解的比较实用的算法

1. Golub-Kahan SVD 算法: 由 Golub 和 Kahan [53] 于 1965 年提出, 并在 1970 年进行了完善 [54], 是一种十分稳定且高效的计算 SVD 的算法, 也是最早的实用 SVD 算法. 主要思想是将带位移的对称 QR 迭代算法隐式地用到 $B^{\mathsf{T}} B$ 上, 在该算法中, 并不需要显示地把 $B^{\mathsf{T}} B$ 计算出来. 该算法也通常就称为 SVD 算法, 是一个基本且实用的算法, 目前仍然是计算小规模矩阵奇异值分解的常用算法. 关于这个算法的详细描述, 也可以参见 [57, 152].

2. dqds 算法: 由 Fernando 和 Parlett [44] 于 1994 年提出, 是计算二对角矩阵所有奇异值的最快算法, 而且能达到很高的相对精度, 包括奇异值很小的情形. 该算法主要基于对 $B^{\mathsf{T}} B$ 的 Cholesky 迭代, 可以看作是 LR 迭代算法的改进. 由于 LR 迭代算法在一定条件下与对称 QR 算法是等价的, 因此该算法也可以看作是 QR 迭代的变形. 该算法是 LAPACK 中计算二对角矩阵奇异值 (不计算奇异向量) 的标准算法 [33].

3. 分而治之法: 该算法是计算较大规模矩阵的所有奇异值和奇异向量的最快算法, 但不能保证小奇异值的相对精度, 即 $\sigma_{i}$ 的相对精度为 $\mathcal{O}(\varepsilon)\sigma_{1}$ , 而不是 $\mathcal{O}(\varepsilon)\sigma_{i}$ .

4. 对分法和反迭代: 主要用于计算某个区间内的奇异值及对应的奇异向量, 能保证较高的相对精度.

5. Jacobi 迭代: 可隐式地对 $AA^{\mathsf{T}}$ 或 $A^{\mathsf{T}}A$ 实施对称 Jacobi 迭代, 能保证较高的相对精度. 最近, Z. Drmac 和 K. Veselic [34, 35] 改进了最初的 Jacobi 算法, 使其变成一个速度快、精度高的实用算法.

在这里, 我们介绍 Golub-Kahan SVD 算法, dqds 算法和 Jacobi 迭代.

5.6.2 Golub-Kahan SVD算法

该算法主要思想是将带位移的对称QR迭代算法隐式地用到 $B^{\mathsf{T}}B$ 上，而无需将 $B^{\mathsf{T}}B$ 显式地计算出来. Golub-Kahan SVD算法有时也简称SVD算法,其基本框架是：

• 将矩阵 $A$ 二对角化, 得到上二对角矩阵 $B$ ;

• 用隐式 QR 迭代计算 $B^{\mathsf{T}} B$ 的特征值分解, 即

• 计算 $BQ$ 的列主元QR分解, 即

其中 $P$ 是置换矩阵, $U$ 是正交矩阵, $R$ 是上三角矩阵

由 (5.8) 可知

因此 $BQ$ 是列正交矩阵 (但不是单位列正交). 再由 (5.9) 可知 $R = U^{\mathsf{T}}(BQ)P$ 也是列正交矩阵. 又 $R$ 是上三角矩阵, 所以 $R$ 必定是对角矩阵. 令 $V = QP$ , 则由 (5.9) 可知

这就是二对角矩阵 $B$ 的奇异值分解

算法的具体实现可参见[57,152].

5.6.3 dqds算法

我们首先介绍针对实对称正定矩阵的LR算法，该算法思想与QR迭代算法类似，但提出时间更早.

算法5.9.带位移的LR算法

1: Let $T_0$ be a given real symmetric positive definite matrix
2: set $i = 0$
3: while not converge do
4: choose a shift $\tau_i^2$ satisfying $\tau_i^2 < \min \{\lambda (T_i)\}$
5: compute $B_i$ such that $T_i - \tau_i^2 I = B_i^\top B_i$ . % Cholesky factorization
6: $T_{i + 1} = B_iB_i^\top +\tau_i^2 I$
7: $i = i + 1$
8: end while

LR迭代算法在形式上与QR迭代算法非常类似.事实上，对于不带位移的LR迭代算法，我们可以证明，两步LR迭代等价于一步QR迭代.

引理5.13 设 $\tilde{T}$ 是不带位移的LR算法迭代两步后生成的矩阵, $\hat{T}$ 是不带位移的QR算法迭代一步后生成的矩阵, 则 $\tilde{T} = \hat{T}$ .

(1) LR 算法中要求 $T_{0}$ 对称正定, 但并不一定是三对角矩阵;

(2) 由该引理可知, QR 算法与 LR 算法有相同的收敛性.

下面我们介绍 dqds (differential quotient difference with shifts) 算法. 该算法是针对三对角的对称正定矩阵 $B^{\mathsf{T}}B$ , 其中 $B$ 是二对角矩阵. 在数学上, dqds 算法与 LR 算法是等价的, 但在该算法中, 我们是直接通过 $B_{i}$ 来计算 $B_{i+1}$ , 从而避免计算中间矩阵 $T_{i+1}$ , 这样也就尽可能地避免了由于计算 $B_{i}B_{i}^{\mathsf{T}}$ 而可能带来的数值不稳定性.

下面推导如何从 $B_{i}$ 直接计算 $B_{i + 1}$ . 设

为了书写方便, 我们记 $b_{0} = b_{n} = \tilde{b}_{0} = \tilde{b}_{n} = 0$ . 由 LR 算法 5.9 可知

比较等式两边矩阵的对角线和上对角线元素，可得

记 $\delta = \tau_{i + 1}^{2} - \tau_{i}^{2},p_{k} = a_{k}^{2},q_{k} = b_{k}^{2},\tilde{p}_{k} = \tilde{a}_{k}^{2},\tilde{q}_{k} = \tilde{b}_{k}^{2}$ ，则可得下面的qds算法.

算法5.10.qds算法的单步 $(B_{i}\to B_{i + 1})$

1: $\delta = \tau_{i + 1}^{2} - \tau_{i}^{2}$
2: for $k = 1$ to $n - 1$ do
3: $\tilde{p}_k = p_k + q_k - \tilde{q}_{k - 1} - \delta$
4: $\tilde{q}_k = q_k\cdot (p_{k + 1} / \tilde{p}_k)$
5: end for
6: $\tilde{p}_n = p_n - \tilde{q}_{n - 1} - \delta$

qds算法中的每个循环仅需5个浮点运算，所以运算量较少

为了提高算法的精确性, 我们引入一个辅助变量 $d_k \triangleq p_k - \tilde{q}_{k-1} - \delta$ , 则

于是就可得到 dqds 算法:

算法5.11.dqds算法的单步 $(B_{i}\to B_{i + 1})$

1: $\delta = \tau_{i + 1}^{2} - \tau_{i}^{2}$
2: $d_{1} = p_{1} - \delta$
3: for $k = 1$ to $n - 1$ do
4: $\tilde{p}_k = d_k + q_k$
5: $t = p_{k + 1} / \tilde{p}_k$
6: $\tilde{q}_k = q_k\cdot t$
7: $d_{k + 1} = d_k\cdot t - \delta$
8: end for
9: $\tilde{p}_n = d_n$

dqds算法的运算量与dqs差不多，但更精确

关于 dqds 算法中位移的选取, 以及如何判断收敛性, 可以参见 [44].

$\mathcal{A}$ dqds算法是求解二对角矩阵所有奇异值的主要方法之一，2000年就被加入到LAPACK中[2,102].关于dqds的最新改进可以参见[90].

5.6.4 Jacobi算法

本节讨论对矩阵 $M = A^{\mathrm{T}}A$ 实施隐式的Jacobi算法来计算 $A$ 的奇异值.

我们知道, Jaboci 算法的每一步就是对矩阵作 Jacobi 旋转, 即 $A^{\mathsf{T}}A \to J^{\mathsf{T}}A^{\mathsf{T}}AJ$ , 其中 $J$ 的选取将某两个非对角元化为 0. 在实际计算中, 我们只需计算 $AJ$ , 故该算法称为单边 Jacobi 旋转.

算法5.12.单边Jacobi旋转的单步

$\%$ 对 $M = A^{\mathsf{T}}A$ 作Jacobi旋转, 将 $M(i,j), M(j,i)$ 化为0

1: Compute $m_{ii} = (A^{\mathsf{T}}A)_{ii}$ , $m_{ij} = (A^{\mathsf{T}}A)_{ij}$ , $m_{jj} = (A^{\mathsf{T}}A)_{jj}$ .
2: if $m_{ij}$ is not small enough then
3: $\tau = (m_{ii} - m_{jj}) / (2\cdot m_{ij})$
4: $t = \mathrm{sign}(\tau) / (|\tau| + \sqrt{1 + \tau^2})$
5: $c = 1 / \sqrt{1 + t^2}$
6: $s = c\cdot t$
7: $A = AG(i,j,\theta)$ $\% G(i,j,\theta)$ 为Givens变换
8: if eigenvectors are desired then
9: $J = J\cdot G(i,j,\theta)$

10: end if
11: end if

在上面算法的基础上, 我们可以给出完整的单边 Jacobi 算法.

算法5.13.单边Jacobi:计算 $A = U\Sigma V^{\mathsf{T}}$

1: while $A^{\mathsf{T}}A$ is not diagonal enough do
2: for $i = 1$ to $n - 1$ do
3: for $j = i + 1$ to $n$ do
4: 调用单边 Jacobi 旋转
5: end for
6: end for
7: end while
8: compute $\sigma_{i} = \| A(:,i)\|_{2}, i = 1,2,\ldots n$
9: $U = [u_{1}, \dots, u_{n}]$ with $u_{i} = A(:, i) / \sigma_{i}$
10: $V = J$

Jacobi算法的特点：

• 不需要双对角化,这样可以避免双对角化引入的误差;

• 可达到相对较高的计算精度；

• 速度较慢. (目前已有快速的改进算法, 参见 [34, 35])

定理5.14 设 $A = DX \in \mathbb{R}^{n \times n}$ , 其中 $D$ 为非奇异对角阵, $X$ 非奇异. 设 $\hat{A}$ 是按浮点运算单边 Jacobi 旋转 $m$ 次后所得到的矩阵. 若 $A$ 和 $\hat{A}$ 的奇异值分别为 $\sigma_{1} \geq \sigma_{2} \geq \ldots \geq \sigma_{n}$ 和 $\hat{\sigma}_{1} \geq \hat{\sigma}_{2} \geq \ldots \geq \hat{\sigma}_{n}$ , 则

故 $X$ 的条件数越小, 计算矩阵 $A$ 的奇异值时相对误差越小.