基于双对角化的Lanczos方法在求解最小二乘问题中的核心作用与步骤
问题描述
对于大型稀疏的最小二乘问题：
\[\min_{x \in \mathbb{R}^n} \|Ax - b\|_2, \]
其中 \(A \in \mathbb{R}^{m \times n}\)（\(m \geq n\)）是一个大型稀疏矩阵，传统的直接法（如正规方程 \(A^T A x = A^T b\) 或QR分解）会破坏稀疏性且计算成本高。基于Krylov子空间的迭代法（如LSQR算法）能有效利用矩阵的稀疏结构，其核心是将 \(A\) 通过双对角化（Bidiagonalization） 转化为上双对角矩阵，从而在低维子空间中逐步逼近解。这一过程本质上是Lanczos方法在最小二乘问题中的变体，我们详细讲解双对角化过程及其在LSQR算法中的关键作用。
1. 核心思想：双对角化过程
双对角化的目标是构造正交矩阵 \(U_k \in \mathbb{R}^{m \times (k+1)}\) 和 \(V_k \in \mathbb{R}^{n \times k}\)，以及上双对角矩阵 \(B_k \in \mathbb{R}^{(k+1) \times k}\)，使得：
\[A V_k = U_{k+1} B_k, \]
其中 \(U_{k+1} = [u_1, u_2, \dots, u_{k+1}]\) 和 \(V_k = [v_1, v_2, \dots, v_k]\) 的列向量分别构成标准正交基，\(B_k\) 的形式为：
\[B_k = \begin{bmatrix} \alpha_1 & \beta_1 & & & \\ & \alpha_2 & \beta_2 & & \\ & & \ddots & \ddots & \\ & & & \alpha_k & \beta_k \\ & & & & \alpha_{k+1} \end{bmatrix}. \]
这样，原始问题被投影到低维Krylov子空间 \(\mathcal{K}_k(A^T A, A^T b)\)，且 \(B_k\) 的稀疏结构使得子问题易于求解。
2. 双对角化的递推构造（Golub-Kahan过程）
我们从初始向量 \(b\) 开始，令 \(u_1 = b / \beta_1\)（其中 \(\beta_1 = \|b\|_2\)），然后交替使用 \(A\) 和 \(A^T\) 生成正交基。
步骤1：初始化
- 计算 \(\beta_1 = \|b\|_2\)，设置 \(u_1 = b / \beta_1\)。
- 计算 \(v_1 = A^T u_1 / \alpha_1\)，其中 \(\alpha_1 = \|A^T u_1\|_2\)。
此时有：
\[A v_1 = \alpha_1 u_1 + \beta_2 u_2, \]
其中 \(\beta_2 u_2 = A v_1 - \alpha_1 u_1\)，且 \(\beta_2 = \|A v_1 - \alpha_1 u_1\|_2\)，\(u_2\) 是单位向量。
步骤2：递推第 \(j\) 步（\(j \geq 1\)）
假设已得到 \(u_j\) 和 \(v_j\)：
- 计算中间向量：
\[ \tilde{v}_{j+1} = A^T u_j - \beta_j v_{j-1} \quad (\text{约定 } v_0 = 0)。 \]
实际上，由于正交性，更精确的公式为：
\[ \tilde{v}_{j+1} = A^T u_j - \alpha_j v_j. \]
这里需要纠正：标准递推关系是：
- 已知 \(u_j\) 和 \(v_j\)，先计算：
\[ \tilde{v}_{j+1} = A^T u_j - \beta_j v_j, \]
其中 $ \beta_j $ 是上一步中与 $ u_j $ 相关的系数（见下文完整流程）。
正确递推流程（Golub-Kahan双对角化）：
- 输入：矩阵 \(A\)，向量 \(b\)。
- 初始化：
\[ \beta_1 = \|b\|_2, \quad u_1 = b / \beta_1, \quad \alpha_1 v_1 = A^T u_1, \quad \alpha_1 = \|A^T u_1\|_2. \]
- 对于 \(j = 1, 2, \dots, k\)：
- 计算 \(v_j\) 和 \(\alpha_j\)（若 \(j=1\) 已初始化，则跳过）：
\[ \alpha_j = \|A^T u_j - \beta_j v_{j-1}\|_2, \quad v_j = (A^T u_j - \beta_j v_{j-1}) / \alpha_j. \]
- 计算 \(u_{j+1}\) 和 \(\beta_{j+1}\)：
\[ w = A v_j - \alpha_j u_j, \quad \beta_{j+1} = \|w\|_2, \quad u_{j+1} = w / \beta_{j+1}. \]
- 此时得到 \(B_k\) 的第 \(j\) 列：对角线元 \(\alpha_j\)，下次对角线元 \(\beta_{j+1}\)。
经过 \(k\) 步，我们有：
\[A V_k = U_{k+1} B_k, \quad U_{k+1}^T U_{k+1} = I_{k+1}, \quad V_k^T V_k = I_k. \]
3. 利用双对角化求解最小二乘问题
将原问题投影到子空间：设解近似为 \(x_k = V_k y_k\)，其中 \(y_k \in \mathbb{R}^k\)。则残差：
\[r_k = b - A x_k = \beta_1 u_1 - A V_k y_k = U_{k+1} (\beta_1 e_1 - B_k y_k), \]
其中 \(e_1 = [1, 0, \dots, 0]^T \in \mathbb{R}^{k+1}\)。由于 \(U_{k+1}\) 列正交，最小化 \(\|r_k\|_2\) 等价于最小化：
\[\min_{y_k \in \mathbb{R}^k} \|\beta_1 e_1 - B_k y_k\|_2. \]
这是一个小型且结构简单的上双对角最小二乘问题，可用QR分解（通过Givens旋转）在 \(O(k)\) 时间内高效求解。
逐步更新解的过程（LSQR算法核心）：
- 每进行一步双对角化，得到新的 \(\alpha_j, \beta_{j+1}\)。
- 用Givens旋转更新 \(B_k\) 的QR分解（不显式存储 \(B_k\)）：
- 维护一个上三角矩阵 \(R_k\) 和变换后的右端项。
- 从QR分解直接得到 \(y_k\) 并更新近似解 \(x_k = V_k y_k\)。
- 监控残差范数 \(\|\beta_1 e_1 - B_k y_k\|_2\)，当满足容忍度时终止迭代。
4. 算法伪代码（LSQR简化版）
- 输入：\(A, b\)，初始猜测 \(x_0 = 0\)，容忍度 \(\epsilon\)，最大迭代步 \(k_{\max}\)。
- 初始化：
- \(\beta_1 = \|b\|_2, \quad u_1 = b / \beta_1\)
- \(\alpha_1 = \|A^T u_1\|_2, \quad v_1 = A^T u_1 / \alpha_1\)
- \(w_1 = v_1\)（用于更新解）
- \(\phi_1 = \beta_1, \quad \rho_1 = \alpha_1\)
- 对于 \(j = 1, 2, \dots, k_{\max}\)：
- 双对角化：
\[ \begin{aligned} &\tilde{u}_{j+1} = A v_j - \alpha_j u_j, \\ &\beta_{j+1} = \|\tilde{u}_{j+1}\|_2, \quad u_{j+1} = \tilde{u}_{j+1} / \beta_{j+1}, \\ &\tilde{v}_{j+1} = A^T u_{j+1} - \beta_{j+1} v_j, \\ &\alpha_{j+1} = \|\tilde{v}_{j+1}\|_2, \quad v_{j+1} = \tilde{v}_{j+1} / \alpha_{j+1}. \end{aligned} \]
- 用Givens旋转更新QR分解：
\[ \begin{aligned} &\rho_j = \sqrt{(\rho_j)^2 + \beta_{j+1}^2}, \quad c_j = \rho_j / \tilde{\rho}_j, \quad s_j = \beta_{j+1} / \tilde{\rho}_j, \\ &\theta_{j+1} = s_j \alpha_{j+1}, \quad \rho_{j+1} = -c_j \alpha_{j+1}, \\ &\phi_j = c_j \phi_j, \quad \phi_{j+1} = s_j \phi_j. \end{aligned} \]
- 更新解：
\[ x_j = x_{j-1} + (\phi_j / \rho_j) w_j, \quad w_{j+1} = v_{j+1} - (\theta_{j+1} / \rho_j) w_j. \]
- 若 \(|\phi_{j+1}| < \epsilon\)，终止并输出 \(x_j\)。
5. 关键性质与优势
- 存储效率：只需存储最近的 \(u_j, v_j, w_j\)，避免大型稠密矩阵。
- 数值稳定性：递推中显式正交化（或通过Lanczos过程隐式保持正交性）。
- 适用于稀疏矩阵：仅需矩阵-向量乘积 \(A v\) 和 \(A^T u\)。
- 自动正则化：双对角化过程可视为对 \(A\) 的逐步正则近似，特别适合病态问题。
6. 总结
基于双对角化的Lanczos方法将最小二乘问题投影到逐步扩张的Krylov子空间，并通过上双对角结构快速求解子问题。LSQR算法是这一思想的经典实现，广泛应用于大型稀疏最小二乘问题（如CT成像、地质反演）。其核心在于交替作用 \(A\) 和 \(A^T\) 构建双对角基，每一步仅需少量内积和向量更新，兼具高效性与稳定性。