当前位置：首页 > news >正文

矩阵奇异值分解（SVD）中Golub–Kahan 双对角化 + 对双对角矩阵的隐式QR详解

news 2025/10/14 8:08:04

总体思路（两阶段）：

正交降阶（reduction）：用一系列正交变换（通常是 Householder 反射）把任意矩阵 $A∈Rm×nA\in\mathbb R^{m\times n}$ 变换成上（或下）双对角（bidiagonal）矩阵 $B$ ，同时累积这些正交变换得到 $U_0,V_0$ ，使得

$U_0 \, B \, V_0^\top,$

其中 $V0∈Rn×nU_0\in\mathbb R^{m\times m},\ V_0\in\mathbb R^{n\times n}$ 是正交矩阵（通常只显式保存产生的 Householder 向量以便以后生成奇异向量）。
对双对角矩阵做迭代（compute SVD of B）：对上步得到的双对角矩阵 $B$ （通常为上双对角：主对角 $αi\alpha_i$ ，上超对角 $βi\beta_i$ ）做迭代求 SVD。常用的数值方法是把问题转为对称三对角矩阵 $T=B⊤BT=B^\top B$ 求特征值（因为 $spec⁡(B⊤B)=σi(B)2\operatorname{spec}(B^\top B)=\sigma_i(B)^2$ ），并使用隐式位移（implicit shift）QR 或 bulge-chasing with Givens 在不形成 $T$ （或显式避免平方）情况下稳定求解奇异值与奇异向量。LAPACK 中的方案（历史上）是先用 bidiagonalize，再用 bdsqr / dbdsqr 或 dgesvd / dgesdd 的细化方法求解。

下面把两部分分别深入展开并给出具体代数变换与伪代码。

A. Householder 双对角化（Golub–Kahan reduction）：如何把 $A$ 变成上双对角 $B$

目标：通过若干左右正交变换把 $A$ 变为上双对角（upper bidiagonal）矩阵 $B$ ：

$\begin{bmatrix} \alpha_1 & \beta_1 & 0 & \cdots\\[4pt] 0 & \alpha_2 & \beta_2 & \cdots\\[4pt] 0 & 0 & \alpha_3 & \ddots\\ \vdots & & \ddots & \ddots \end{bmatrix}.$

基本步骤（第 $k$ 步， $k=1,…,min⁡(m,n)k=1,\dots,\min(m,n)$ ）

用左 Householder 把当前矩阵的第 $k$ 列在第 $k$ 行以下的元素全部零掉（保留位置 $(k, k)$ ）。设当前矩阵为 $A^{(k)}$ ，取列向量

$A^{(k)}_{k:m,\ k}\in\mathbb R^{m-k+1}.$

构造 Householder 向量 $v$ （例如 $\pm \|x\| e_1$ ）并得到 $Hk=I−2vv⊤v⊤vH_k = I - 2\frac{v v^\top}{v^\top v}$ ，将其扩展为 $ \tilde H_k = \begin{bmatrix}I_{k-1} & 0\ 0 & H_k\end{bmatrix}$。然后

$A^{(k+\tfrac12)} = \tilde H_k A^{(k)}.$

这一步把第 $k$ 列中第 $k + 1$ 到 $m$ 行置零；主对角位置 $(k, k)$ 变成某个数，记作 $αk\alpha_k$ .
用右 Householder（作用在列方向）把第 $k$ 行在第 $k + 1$ 列以后的元素全部清零，保留上三角位置 $(k, k + 1)$ （即生成 $βk\beta_k$ ）。取行向量（转置后作为列向量）

$KaTeX parse error: Double superscript at position 34: …2)}_{k,\ k+1:n}^̲\top \in \mathb…$

构造右向量 $w$ 做 Householder $Gk=I−2ww⊤w⊤wG_k = I - 2\frac{w w^\top}{w^\top w}$ （内部维度为 $n - k$ ），扩展为 $G~k=[Ik00Gk]\tilde G_k = \begin{bmatrix} I_k & 0 \\ 0 & G_k\end{bmatrix}$ ，然后

$A^{(k+1)} = A^{(k+\tfrac12)} \tilde G_k.$

这一步把第 $k$ 行在 $k+2,…,nk+2,\dots,n$ 的元素置零，只留下主对角（ $αk\alpha_k$ ）和超对角（ $βk\beta_k$ ）的位置。

循环直到 $k=min⁡(m,n)k=\min(m,n)$ 。最终 $B=A^{(1)}$ 的形式即为上双对角。

代数要点

每一步左 / 右 Householder 都是正交（ $Hk⊤Hk=IH_k^\top H_k=I$ ），所以总体上

$\bigl(\tilde H_1\tilde H_2\cdots\bigr)^\top \; B \; \bigl(\tilde G_1\tilde G_2\cdots\bigr) .$

你可以把 $U0=∏H~i⊤U_0 = \prod \tilde H_i^\top$ 与 $V0=∏G~iV_0 = \prod \tilde G_i$ 累积起来以便最终得到奇异向量（若需要完整的 $U, V$ ）。
实现注意：不需要显式形成 $U_0, V_0$ 的完整矩阵（代价高），只存 Householder 向量，并在需要奇异向量时回溯应用 Householder。

伪代码（简洁版）

Input: A ∈ R^{m×n}
for k = 1..min(m,n):# left reflect to zero below diagonal in col kx = A[k:m, k]v = householder_vector(x)           # v produces H = I - 2 vv^T/(v^T v)A[k:m, k:n] = (I - 2 vv^T/(v^T v)) * A[k:m, k:n]store v in array for U accumulationα_k = A[k,k]if k < n:# right reflect to zero beyond superdiagonal in row ky = A[k, k+1:n]^Tw = householder_vector(y)A[1:m, k+1:n] = A[1:m, k+1:n] * (I - 2 ww^T/(w^T w))store w for V accumulationβ_k = A[k, k+1]
end
Result: B is upper bidiagonal with diagonals α_i and superdiagonals β_i

B. 从双对角到 SVD：为何把问题转为 $T=B⊤BT=B^\top B$ （对称三对角）以及隐式 QR 的思想

已得上双对角 $B$ （假设 $m≥nm\ge n$ ， $B$ 为 $m×nm\times n$ 上双对角）。要得到 $B$ 的奇异值与奇异向量，可以通过 $T=B⊤BT=B^\top B$ （ $n×nn\times n$ 对称三对角）来做特征分解，因为

$B^\top B = T \quad\Rightarrow\quad \text{eig}(T)=\sigma_i(B)^2.$

$T=B⊤BT=B^\top B$ 的显式三对角结构（下标从 1 开始）

设 $B$ 的主对角为 $α1,…,αn\alpha_1,\dots,\alpha_n$ ，超对角为 $β1,…,βn−1\beta_1,\dots,\beta_{n-1}$ （注意某些索引处 $β0\beta_0$ 视作 0）。则

$B^\top B = \begin{bmatrix}\alpha_1^2 & \alpha_1\beta_1 & & \\\alpha_1\beta_1 & \beta_1^2+\alpha_2^2 & \alpha_2\beta_2 & \\& \ddots & \ddots & \ddots \\& & \alpha_{n-1}\beta_{n-1} & \beta_{n-1}^2+\alpha_n^2 \end{bmatrix}.$

（第 $i$ 个对角项一般是 $αi2+βi−12\alpha_i^2+\beta_{i-1}^2$ ，边界处 $β0=βn=0\beta_0=\beta_n=0$ 处理即可。）

为什么不直接对 $T$ 做普通 QR（而要做隐式、在 $B$ 上 chase）？

在数值实现上，直接构造 $T=B⊤BT=B^\top B$ 会平方条件数并放大舍入误差（ $B$ 的小奇异值的数值不稳定）。更稳定的做法是在不显式构造 $T$ 的前提下，利用 $T=B⊤BT=B^\top B$ 的结构，在 $B$ 上以 Givens 旋转“追踪（bulge-chasing）”来实现对 $T$ 的隐式 QR 步。这样既保留了矩阵的稀疏（带）结构，也能用数值上稳定的 Givens 旋转来实现相似变换。

隐式位移 QR 的核心思想（对称三对角）

隐式 QR（with shift $μ\mu$ ）对对称三对角矩阵 $T$ 做一次迭代相当于：

$T-\mu I = Q R\quad\Rightarrow\quad T' = RQ + \mu I = Q^\top T Q,$

于是特征值不变。对三对角矩阵，这个 QR 步可通过一系列 Givens 旋转（消除第一列下元素，然后“bulge chase”）高效完成。把这些旋转映射回 $B$ 上，会对应于在 $B$ 上轮流右乘列旋转（对 $V$ 的更新）和左乘行旋转（对 $U$ 的更新），整个过程把三对角隐式地做了相似变换，而在 $B$ 上只改变对角 $α\alpha$ 与 $β\beta$ 并且保持其“几乎双对角”的结构（中间会有短时’bulge’，会被追赶下去，最终恢复双对角并完成一次带位移的 QR 步）。

C. 在 $B$ 上做一次隐式位移 QR 的「操作级」说明（bulge-chasing；右旋／左旋交替）

下面给出一份常见实现（即 Golub–Kahan/Reinsch 风格）的操作级步骤与关键更新公式（对实现者友好）。

数据结构约定

有数组 $⁣n]\alpha[1\!:\!n]$ 表示主对角（ $n$ 个）； $⁣n−1]\beta[1\!:\!n-1]$ 表示超对角；
我们隐式维护 $B$ 的上双对角结构；迭代过程中 $βi\beta_i$ 会逐步变小以实现对角化（ $βi→0\beta_i\to 0$ 表示问题分裂 / deflation）。

一次带位移的循环（高层伪代码）

while not converged:# (1) 处理可拆分点：若某 β_j 足够小 -> problem splits; 处理每个子块 separately# (2) 选择位移 μ（通常使用底部 1x1 或 2x2 子块的 Wilkinson shift）μ = choose_shift(alpha, beta)  # (3) 用 μ 在 B 上做隐式 QR 步（在 B 上做 Givens bulge-chase）# 初始化：构造用于开始 chase 的初始右旋 R_1compute initial (c,s) from (alpha[1]^2 - μ, alpha[1]*beta[1])apply right rotation R_1 to columns 1&2 of B  # 更新 alpha[1], beta[1], alpha[2] 并产生 bulgefor j = 1 .. n-1:# 消除由右旋产生的 bulge：用左旋 L_j 作用于行 j+1 和 j+2compute (cL,sL) to annihilate bulge; apply left rotation -> update alpha[j+1], beta[j], etc.# 接着右旋 R_{j+1} 作用于列 j+1 & j+2 以继续 chase bulgecompute (cR,sR) from current (alpha[j+1], beta[j+1], ...) ; apply right rotation -> update entries,产生下一个 bulgeendfor# 经过一次完整 chase，B 被转换到新的上双对角形式，相应地 alpha/beta 被更新
endwhile

关键代数：右旋与左旋对 $α,β\alpha,\beta$ 的更新（局部公式）

（下面只列出初始几步的清晰公式；一般步骤为这些公式的循环推广。实现时建议用 hypot / 稳定求旋转分量。）

右旋（作用于列 j 与列 j+1）：设一次右旋参数 $(c, s)$ （二维正交矩阵 $[cs−sc]\begin{bmatrix} c & s\\ -s & c\end{bmatrix}$ ），作用只影响列 $j$ 与 $j + 1$ ：

初始（旋转前）在相关行列位置的非零元素（局部）看作：
- $αj\alpha_j$ （出现在行 $j$ , 列 $j$ ）
- $βj\beta_j$ （出现在行 $j$ , 列 $j + 1$ ）
- $αj+1\alpha_{j+1}$ （出现在行 $j + 1$ , 列 $j + 1$ ）
右旋后局部更新（局部代数）：

$\begin{aligned} \alpha_j' & = c\,\alpha_j + s\,\beta_j,\\ \beta_j' & = -s\,\alpha_j + c\,\beta_j,\\ \text{bulge } g & = s\,\alpha_{j+1},\\ \alpha_{j+1}' & = c\,\alpha_{j+1}. \end{aligned}$

右旋会在位置 $(j + 1, j)$ 生成一个“bulge”值 $g$ （即在下三角出现一个短暂非零），破坏了双对角形式；接下来用左旋来消除该 bulge。

左旋（作用于行 j+1 与 j+2）：左旋用来消除位置 $(j + 1, j)$ 的 bulge。假设要消掉的向量为 $αj+1′]⊤[g,\ \alpha_{j+1}']^\top$ ，选 $c_L,s_L)$ 使得

$\begin{bmatrix} c_L & s_L\\ -s_L & c_L\end{bmatrix} \begin{bmatrix} g \\ \alpha_{j+1}' \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix}.$

左旋作用后会更新（示意）：

把 bulge 清零；
更新 $βj′\beta_j'$ （因为第 j 列与第 j+1 行的乘积受影响）；
产生可能的新 bulge（在下一位置），這個新 bulge 随后被下一个右旋处理。

注：上面右旋 / 左旋更新写法是局部的、示意性的；实现里每次 Givens 旋转都用 hypot(a,b) 稳定求 $c, s$ ，然后就按矩阵乘法在受影响的两个行/列上更新常数（常数个数目很少，复杂度 O(1) per rotation）。整个 chase 对长度 $n$ 的双对角块需要 O(n) 次旋转来推进一个 bulge 到底部。

初始选择位移 μ（Wilkinson shift）

实际算法中常选取底部 1×1 或 2×2 子块的 Wilkinson shift（对对称三对角的经典选择），这是加速收敛且保留数值稳定性的经验/理论做法。

deflation 与结束条件

如果某个 $∣βi∣|\beta_i|$ 足够小（比如与周围 $α\alpha$ 的量级相比低很多或低于机器精度阈值），就认为该位置已对角化（deflated），可以把问题分成左右两块分别处理（这显著加速实际收敛与复杂度）。
结束条件一般是当所有 $∣βi∣<ε|\beta_i| < \varepsilon$ （机器精度相关）时， $αi\alpha_i$ 就是奇异值（对数值细节可再正号化为非负）。

如何得到奇异向量 $U, V$

在 chase 的过程中同时累积所用的右旋会更新 $V$ （列旋转对应 $V$ 左乘/右乘的更新），累积左旋更新 $U$ （行旋转对应 $U$ ）——这正是保持 $V^\top$ 的变换关系的方式。
真实实现中通常不在每次旋转就去更新完整 $U, V$ （那样代价 O(n^3)），而是记录旋转并在需要奇异向量时后向应用这些旋转以构造 $U, V$ 。

D. 小例子（n=3 的“手算”式追踪一次初始 step，便于理解）

设初始上双对角 $B$ 的参数为 $α1,β1,α2,β2,α3\alpha_1,\beta_1,\alpha_2,\beta_2,\alpha_3$ 。取某位移 $μ\mu$ 。计算首个右旋（针对列 1&2）：

先形成向量 $α1β1)⊤\big(\alpha_1^2-\mu,\ \alpha_1\beta_1\big)^\top$ ，用 Givens 得到 $(c, s)$ 使得第二分量为 0（这与在 $T=B⊤BT=B^\top B$ 上对第一列做 QR 等价）。
右旋后更新（按上面公式）：

$\begin{aligned} \alpha_1 &\leftarrow c\alpha_1 + s\beta_1,\\ \beta_1 &\leftarrow -s\alpha_1 + c\beta_1,\\ \text{bulge } g &= s\alpha_2,\\ \alpha_2 &\leftarrow c\alpha_2. \end{aligned}$
对 $g,α2g,\alpha_2$ 做左旋（求 $c_L,s_L$ 使 $[g,α2]⊤[g,\alpha_2]^\top$ 的第二分量被消掉），左旋会影响 $β1\beta_1$ 与 $α2\alpha_2$ 并可能对 $β2\beta_2$ 造成影响（产生下一级 bulge）——如此交替向下追赶直到到达底部。底部 bulge 消失后，一次完整的隐式位移 QR 步完成。

实现时不要把这些更新写成大量矩阵乘法（那样慢）；正确的做法是对于每个 Givens 用 4–6 个标量更新式局部修改 $α,β\alpha,\beta$ （见上面的局部公式），这使每次旋转是 O(1) 操作，每次完整的 chase 是 O(n)。

E. 数值细节、稳定性与实践建议

在构造 Givens 参数时用 hypot(a,b) 来获得 $(c, s)$ 的稳定计算，避免直接做 $a2+b2\sqrt{a^2+b^2}$ 导致上溢/下溢。
在操作中尽量避免显式形成 $T=B⊤BT=B^\top B$ ，通过在 $B$ 上的旋转来实现等价的相似变换（这减少误差倍增）。
在实现中要检测并做 deflation（ $)|\beta_i| < \epsilon(\,|\alpha_i|+|\alpha_{i+1}|\;)$ ），这样把问题切分成更小子问题并节省计算量。
对于求奇异值谱（而非全部奇异向量），可以只运行求值分支（不累积 $U, V$ ），这样更节省内存与时间。
大型 dense SVD 的高效实现一般使用分治法（divide-and-conquer，LAPACK gesdd）或 dqds / bidiag-SVD 的改良版本；隐式 QR（Golub–Kahan）在很多实现里仍用于较小块或精细控制数值稳定性。

小结 / 要点回顾

Golub–Kahan + 隐式 QR：数值上 SVD 实现分为两阶段（Householder bidiagonalization + 对 bidiagonal B 的迭代求 SVD）。后者通过在 $B$ 上做交替的右/左 Givens 旋转（bulge-chasing），隐式等价于对 $T=B⊤BT=B^\top B$ 做对称 QR 且数值更稳定。实现关键在于局部旋转的稳定计算、恰当选择位移（Wilkinson）、及时的 deflation，以及累积旋转（如果需要奇异向量）。