基于Krylov子空间的LSQR算法求解线性最小二乘问题
题目描述
给定一个超定线性方程组 \(A \mathbf{x} = \mathbf{b}\)，其中 \(A \in \mathbb{R}^{m \times n} \ (m \geq n)\)，\(\mathbf{b} \in \mathbb{R}^{m}\)。我们希望通过最小二乘法求解 \(\min_{\mathbf{x}} \| A\mathbf{x} - \mathbf{b} \|_2\)。当矩阵 \(A\) 规模非常大且稀疏时，直接使用正规方程 \(A^T A \mathbf{x} = A^T \mathbf{b}\) 可能因计算 \(A^T A\) 导致数值条件数平方增长而变得不稳定。LSQR算法是一种基于Krylov子空间的迭代方法，它通过隐式的Lanczos双对角化过程，在不显式计算 \(A^T A\) 的情况下，稳定地求解最小二乘问题。
解题过程
1. 问题转化
最小二乘问题等价于求解正规方程：
\[ A^T A \mathbf{x} = A^T \mathbf{b} \]
但如果直接求解，\(A^T A\) 的条件数是 \(\kappa(A)^2\)，可能引发数值不稳定。LSQR的核心思想是：通过构造Krylov子空间 \(\mathcal{K}_k(A^T A, A^T \mathbf{b})\) 来迭代逼近解，但避免显式构造 \(A^T A\)。
2. Golub-Kahan双对角化过程
LSQR的基础是对矩阵 \(A\) 进行双对角化。其步骤如下：
- 初始化：
\[ \beta_1 = \|\mathbf{b}\|_2, \quad \mathbf{u}_1 = \mathbf{b} / \beta_1 \]
\[ \alpha_1 = \|A^T \mathbf{u}_1\|_2, \quad \mathbf{v}_1 = A^T \mathbf{u}_1 / \alpha_1 \]
令 \(\mathbf{w}_1 = \mathbf{v}_1\)，并设初始解估计 \(\mathbf{x}_0 = \mathbf{0}\)。
- 迭代双对角化（对 \(k = 1, 2, \dots\)）：
- 前向乘法：
\[ \hat{\mathbf{u}}_{k+1} = A \mathbf{v}_k - \alpha_k \mathbf{u}_k \]
\[ \beta_{k+1} = \|\hat{\mathbf{u}}_{k+1}\|_2, \quad \mathbf{u}_{k+1} = \hat{\mathbf{u}}_{k+1} / \beta_{k+1} \]
- 后向乘法：
\[ \hat{\mathbf{v}}_{k+1} = A^T \mathbf{u}_{k+1} - \beta_{k+1} \mathbf{v}_k \]
\[ \alpha_{k+1} = \|\hat{\mathbf{v}}_{k+1}\|_2, \quad \mathbf{v}_{k+1} = \hat{\mathbf{v}}_{k+1} / \alpha_{k+1} \]
这一过程可写为矩阵形式：
\[ A V_k = U_{k+1} B_k, \quad A^T U_{k+1} = V_k B_k^T + \alpha_{k+1} \mathbf{v}_{k+1} \mathbf{e}_{k+1}^T \]
其中 \(B_k\) 是 \((k+1) \times k\) 的双对角矩阵：
\[ B_k = \begin{pmatrix} \alpha_1 & & \\ \beta_2 & \alpha_2 & \\ & \beta_3 & \ddots & \\ & & \ddots & \alpha_k \\ & & & \beta_{k+1} \end{pmatrix} \]
且 \(U_{k+1} = [\mathbf{u}_1, \dots, \mathbf{u}_{k+1}]\)、\(V_k = [\mathbf{v}_1, \dots, \mathbf{v}_k]\) 的列向量相互正交。
3. 子空间投影与约化问题
在第 \(k\) 步，我们在子空间 \(\mathcal{K}_k(A^T A, A^T \mathbf{b}) = \mathrm{span}\{V_k\}\) 中寻找近似解 \(\mathbf{x}_k = V_k \mathbf{y}_k\)。代入最小二乘目标：
\[ \min_{\mathbf{x} \in \mathrm{span}\{V_k\}} \|A \mathbf{x} - \mathbf{b}\|_2 = \min_{\mathbf{y}_k} \|A V_k \mathbf{y}_k - \mathbf{b}\|_2 \]
利用 \(A V_k = U_{k+1} B_k\) 和 \(\mathbf{b} = \beta_1 \mathbf{u}_1 = U_{k+1} (\beta_1 \mathbf{e}_1)\)（其中 \(\mathbf{e}_1 = (1,0,\dots,0)^T \in \mathbb{R}^{k+1}\)），得：
\[ \|A V_k \mathbf{y}_k - \mathbf{b}\|_2 = \|U_{k+1} (B_k \mathbf{y}_k - \beta_1 \mathbf{e}_1)\|_2 = \|B_k \mathbf{y}_k - \beta_1 \mathbf{e}_1\|_2 \]
于是问题转化为求解小规模最小二乘问题：
\[ \min_{\mathbf{y}_k} \| B_k \mathbf{y}_k - \beta_1 \mathbf{e}_1 \|_2 \]
4. 通过QR分解迭代更新解
为了避免保存所有 \(V_k\)，LSQR利用递推关系更新解。对双对角矩阵 \(B_k\) 进行QR分解：
\[ Q_k B_k = \begin{pmatrix} R_k \\ \mathbf{0}^T \end{pmatrix}, \quad Q_k (\beta_1 \mathbf{e}_1) = \begin{pmatrix} \mathbf{f}_k \\ \phi_{k+1} \end{pmatrix} \]
其中 \(R_k\) 是 \(k \times k\) 上三角矩阵，\(\mathbf{f}_k \in \mathbb{R}^k\)，\(\phi_{k+1}\) 是标量。则最小二乘解为 \(\mathbf{y}_k = R_k^{-1} \mathbf{f}_k\)，且残差范数 \(\rho_k = |\phi_{k+1}|\)。
- 更新解：设 \(P_k = V_k R_k^{-1}\)（其列可递推计算），则 \(\mathbf{x}_k = P_k \mathbf{f}_k\)。
- 递推公式（推导略去细节）：
\[ \mathbf{x}_k = \mathbf{x}_{k-1} + \phi_k \mathbf{p}_k \]
其中 \(\mathbf{p}_k\) 根据 \(\mathbf{v}_k\) 和前一个 \(\mathbf{p}_{k-1}\) 更新得到。
实际算法中，我们维护标量 \(\rho_k, \theta_{k+1}, \phi_k, \bar{\phi}_k\) 等，并通过类似Lanczos系数的递推避免存储所有基向量。
5. 算法步骤摘要
- 初始化：
\[ \beta_1 = \|\mathbf{b}\|, \ \mathbf{u}_1 = \mathbf{b}/\beta_1, \ \alpha_1 = \|A^T \mathbf{u}_1\|, \ \mathbf{v}_1 = A^T \mathbf{u}_1/\alpha_1 \]
\[ \mathbf{w}_1 = \mathbf{v}_1, \ \mathbf{x}_0 = \mathbf{0}, \ \bar{\phi}_1 = \beta_1, \ \bar{\rho}_1 = \alpha_1 \]
- 迭代（对 \(k = 1, 2, \dots\) 直到收敛）：
- 双对角化：
\[ \mathbf{u}_{k+1} = (A \mathbf{v}_k - \alpha_k \mathbf{u}_k)/\beta_{k+1} \]
\[ \mathbf{v}_{k+1} = (A^T \mathbf{u}_{k+1} - \beta_{k+1} \mathbf{v}_k)/\alpha_{k+1} \]
- 构造并旋转Givens平面：
\[ \rho_k = \sqrt{\bar{\rho}_k^2 + \beta_{k+1}^2}, \ c_k = \bar{\rho}_k / \rho_k, \ s_k = \beta_{k+1} / \rho_k \]
\[ \theta_{k+1} = s_k \alpha_{k+1}, \ \bar{\rho}_{k+1} = -c_k \alpha_{k+1}, \ \phi_k = c_k \bar{\phi}_k, \ \bar{\phi}_{k+1} = s_k \bar{\phi}_k \]
- 更新解与搜索方向：
\[ \mathbf{p}_k = (\mathbf{v}_k - \theta_k \mathbf{p}_{k-1}) / \rho_k \quad (\text{其中 } \theta_1 = 0) \]
\[ \mathbf{x}_k = \mathbf{x}_{k-1} + \phi_k \mathbf{p}_k \]
- 检查收敛：若 \(|\bar{\phi}_{k+1}|\) 小于容忍度，则停止。
6. 收敛性与特点
- LSQR本质上是共轭梯度法（CG）在正规方程上的应用，但通过双对角化避免平方条件数。
- 每次迭代主要计算两次矩阵-向量乘法（\(A \mathbf{v}_k\) 和 \(A^T \mathbf{u}_{k+1}\)），适合大规模稀疏矩阵。
- 可自然估计残差范数 \(\|\mathbf{b} - A \mathbf{x}_k\| = |\bar{\phi}_{k+1}|\)，便于设置停止准则。
- 对于离散不适定问题，可通过早期停止实现正则化效果。
总结：LSQR通过隐式双对角化在Krylov子空间中构建近似解，避免了正规方程的病态问题，同时以较低的存储和计算成本迭代求解大规模稀疏最小二乘问题，是科学与工程计算中常用的稳定算法。