SVD 的三步走：双对角化 Givens 收敛 排序

写在前面：万能的 SVD，缺席的算法

SVD 是线性代数的瑞士军刀。

你做主成分分析（PCA），底层是 SVD；你做推荐系统的协同过滤，底层是 SVD；你算伪逆、解最小二乘，底层是 SVD；你做图像压缩、信号去噪、潜在语义分析（LSA），底层还是 SVD。统计软件里凡是涉及"降维""求秩""解超定方程组"的功能，几乎绕不开它。

但翻开任何一本线性代数教材，SVD 这一章的写法几乎一模一样：先证明对任意矩阵 $A$ 都存在分解 $A=U\mathrm{\Sigma }{V}^T$，其中 $U$、$V$ 是正交阵、$\mathrm{\Sigma }$ 是对角阵、对角元叫奇异值，然后画一个示意图，再举一个 2×2 的例子手算一遍，翻篇。

没人告诉你，怎么把一个 1000×1000 的矩阵数值地算成 $U\mathrm{\Sigma }{V}^T$。

这是教科书和生产代码之间最大的一道鸿沟之一。因为真要算出来，你立刻撞上一堆教科书从来不提的问题：

• 能不能直接算 $A^{T}A$ 的特征值？为什么实际代码不这么干？
• $A$ 不是对称矩阵，怎么办？
• 怎么保证算出来的 $U$、$V$ 真的是正交阵？
• 迭代收敛到什么程度才算够？什么时候该认输？
• 奇异值算出来是无序的，怎么保证降序输出？

今天我们用 这份约 300 行的真实 C 源码 svd.c，把这些问题一个一个填平。它是 LINPACK 的 hsvdc 的 C 移植，忠实保留了 1970 年代末期那批数值分析大师（其中就有 MATLAB 之父 Cleve Moler）写下的每一行工程决策。

第一个陷阱：为什么不能直接算 $A^{T}A$ 的特征值

先说数学。奇异值的定义是这样的：$A^{T}A$ 是一个对称半正定矩阵，它的特征值全是非负实数，把这些特征值开方，就是奇异值 ${\sigma }_i=\sqrt{{\lambda }_i(A^{T}A)}$

​，对应的特征向量就是 $V$。然后 $U=AV{\mathrm{\Sigma }}^{-1}$。

看起来这条路最直接——算特征值而已，对称矩阵的特征值算法一抓一大把，为什么要单独发明 SVD？

因为条件数会平方。

回忆条件数的定义：$\kappa (A)={\sigma }_{\max }\mathrm{/}{\sigma }_{\min }$，衡量矩阵在数值上"病态"的程度。一个 $\kappa (A)=10^8$ 的矩阵，在 double 精度（约 15~16 位有效数字）下，做线性求解就会损失 8 位精度，几乎不能用了。

现在你想通过 $A^{T}A$ 来算奇异值。注意 $A^{T}A$ 的特征值是 ${\sigma }_i^2$，于是 $\kappa (A^{T}A)={\sigma }_{\max }^{}\mathrm{/}{\sigma }_{\min }^{}=\kappa (A{)}^2$。原本 $10^8$ 的条件数，平方之后变成 $10^{16}$——直接超出了 double 的精度极限。换句话说，原本还能勉强算的矩阵，一旦你先算 $A^{T}A$ 再做特征分解，最小那几个奇异值就完全淹没在浮点噪声里了。

举个直观的数：若 ${\sigma }_{\min }=10^{-4}$、${\sigma }_{\max }=10^4$，那么 $\kappa (A)=10^8$，仍然在 double 的可处理范围；但 ${\sigma }_{\min }^{}=10^{-8}$、${\sigma }_{\max }^{}=10^8$，$\kappa (A^{T}A)=10^{16}$，那些小特征值会被舍入误差完全吞掉。

所以生产代码从不直接形成 $A^{T}A$。 正确的做法是对 $A$ 本身做正交变换，把奇异值直接"挤"出来，不让它先被平方。这就是 SVD 数值算法的全部出发点。

svd.c 开头的文件注释就写得很明白：

/*============================================================================* *svd.c —— 奇异值分解 X = U·S·Vᵀ 的 C 移植 *对应 Fortran：hsvdc.for * *数学： *第一阶段「双对角化」：交替用左 Householder（清列）与右 Householder *（清行次对角），把任意 n×p 矩阵化成双对角阵 B（仅主、次两条对角）。 *第二阶段「Givens QR 收敛」：对双对角阵做带 Wilkinson 隐式位移的 *QR 迭代（成对 Givens 旋转），把次对角 e 逐步清零，主对角 s 收敛为 *奇异值；同步把变换累积进 U（左奇异向量）与 V（右奇异向量）。 *============================================================================*/

注意"双对角化"和"Givens QR 收敛"两个词，这就是 SVD 的两阶段。但真正读懂源码，你会发现还有隐藏的第三阶段——排序。所以本文标题叫"三步走"。

三阶段算法总览

在读代码之前，先把整体框架画清楚。SVD 数值算法分成三步：

第一阶段：Householder 双对角化。 用一系列 Householder 反射（左右交替），把任意 $n\times p$ 矩阵 $A$ 化成一个只有主对角线和次对角线有非零元素的"双对角阵" $B$：