基于Krylov子空间的LSMR算法求解线性最小二乘问题
我将为你详细讲解这个算法，包括其背景、基本原理、推导过程和实现细节。
1. 问题背景与问题描述
线性最小二乘问题是线性代数中的经典问题：
给定矩阵 \(A \in \mathbb{R}^{m \times n}\)（通常 \(m \geq n\)）和向量 \(b \in \mathbb{R}^m\)，我们希望找到 \(x \in \mathbb{R}^n\) 使得残差范数最小：
\[\min_{x} \|Ax - b\|_2 \]
传统解法的问题：
- 正规方程法：求解 \(A^T A x = A^T b\)，但当 \(A\) 条件数大时，\(A^T A\) 的条件数是 \(A\) 条件数的平方，数值不稳定
- QR分解：数值稳定但计算成本高，对大规模稀疏矩阵不高效
- 共轭梯度法（CG）：只能用于对称正定系统
LSMR的诞生：
LSMR（Least Squares Minimal Residual）是Christopher C. Paige和Michael A. Saunders在2011年提出的算法，专门为求解最小二乘问题而设计，是著名LSQR算法的改进版本。
2. 算法核心思想
LSMR的核心思想是：在Krylov子空间中寻找最小残差的解
Krylov子空间定义：
\[\mathcal{K}_k(A^T A, A^T b) = \text{span}\{A^T b, (A^T A)A^T b, \dots, (A^T A)^{k-1}A^T b\} \]
但LSMR实际上工作在更优的空间上。
3. 算法推导步骤
步骤1：双对角化过程（Golub-Kahan双对角化）
这是算法的关键预处理步骤：
初始化：
\[\beta_1 u_1 = b, \quad \alpha_1 v_1 = A^T u_1 \]
其中 \(\beta_1 = \|b\|_2\)，\(\alpha_1 = \|A^T u_1\|_2\)，且 \(u_1\) 和 \(v_1\) 是单位向量。
迭代过程（对 \(i = 1, 2, \dots\)）：
- 正向迭代：
\[\beta_{i+1} u_{i+1} = A v_i - \alpha_i u_i \]
其中 \(\beta_{i+1} = \|A v_i - \alpha_i u_i\|_2\)
- 反向迭代：
\[\alpha_{i+1} v_{i+1} = A^T u_{i+1} - \beta_{i+1} v_i \]
其中 \(\alpha_{i+1} = \|A^T u_{i+1} - \beta_{i+1} v_i\|_2\)
结果：
经过 \(k\) 步迭代后，我们得到：
\[U_{k+1} = [u_1, u_2, \dots, u_{k+1}] \in \mathbb{R}^{m \times (k+1)}, \quad U_{k+1}^T U_{k+1} = I \]
\[ V_k = [v_1, v_2, \dots, v_k] \in \mathbb{R}^{n \times k}, \quad V_k^T V_k = I \]
\[ B_k = \begin{bmatrix} \alpha_1 \\ \beta_2 & \alpha_2 \\ & \beta_3 & \alpha_3 \\ & & \ddots & \ddots \\ & & & \beta_k & \alpha_k \\ & & & & \beta_{k+1} \end{bmatrix} \in \mathbb{R}^{(k+1) \times k} \]
满足：
\[A V_k = U_{k+1} B_k \]
\[ A^T U_{k+1} = V_k B_k^T + \alpha_{k+1} v_{k+1} e_{k+1}^T \]
步骤2：问题转化
原始问题转化为：
\[\min_{x \in \mathcal{K}_k(A^T A, A^T b)} \|Ax - b\|_2 \]
由于 \(x = V_k y\) 对于某个 \(y \in \mathbb{R}^k\)，我们有：
\[\|Ax - b\|_2 = \|A V_k y - b\|_2 = \|U_{k+1} B_k y - \beta_1 e_1\|_2 = \|B_k y - \beta_1 e_1\|_2 \]
因为 \(U_{k+1}\) 是列正交的。
关键点：原 \(m \times n\) 问题转化为 \((k+1) \times k\) 的双对角最小二乘问题，维度大大降低。
步骤3：LSMR的优化准则
与LSQR不同，LSMR最小化的是：
\[\|A^T (Ax - b)\|_2 \]
等价于正规方程的残差范数。
推导：
\[A^T (A V_k y - b) = A^T (U_{k+1} B_k y - \beta_1 u_1) \]
\[ = (V_k B_k^T + \alpha_{k+1} v_{k+1} e_{k+1}^T) B_k y - \beta_1 A^T u_1 \]
\[ = V_k B_k^T B_k y + \alpha_{k+1} v_{k+1} e_{k+1}^T B_k y - \beta_1 \alpha_1 v_1 \]
定义 \(T_k = B_k^T B_k\)，这是一个对称三对角矩阵：
\[T_k = \begin{bmatrix} \alpha_1^2 + \beta_2^2 & \alpha_2 \beta_2 \\ \alpha_2 \beta_2 & \alpha_2^2 + \beta_3^2 & \alpha_3 \beta_3 \\ & \alpha_3 \beta_3 & \alpha_3^2 + \beta_4^2 & \ddots \\ & & \ddots & \ddots & \alpha_k \beta_k \\ & & & \alpha_k \beta_k & \alpha_k^2 + \beta_{k+1}^2 \end{bmatrix} \]
步骤4：递推求解
LSMR通过QR分解和Givens旋转高效求解：
定义矩阵：
\[R_k = \begin{bmatrix} \rho_1 & \theta_2 \\ & \rho_2 & \theta_3 \\ & & \ddots & \ddots \\ & & & \rho_{k-1} & \theta_k \\ & & & & \rho_k \end{bmatrix}, \quad \bar{R}_k = \begin{bmatrix} R_k \\ \phi_{k+1} e_k^T \end{bmatrix} \]
递推关系：
- 初始化：
\[\rho_1 = \alpha_1, \quad \bar{\phi}_1 = \beta_1, \quad \phi_1 = \beta_1 \]
- 对于 \(i = 1, 2, \dots\):
a. 应用Givens旋转消去 \(\beta_{i+1}\):
\[ \begin{bmatrix} c_i & s_i \\ -s_i & c_i \end{bmatrix} \begin{bmatrix} \rho_i \\ \beta_{i+1} \end{bmatrix} = \begin{bmatrix} \hat{\rho}_i \\ 0 \end{bmatrix} \]
其中 \(c_i^2 + s_i^2 = 1\)
b. 更新：
\[ \theta_{i+1} = s_i \alpha_{i+1}, \quad \rho_{i+1} = -c_i \alpha_{i+1} \]
c. 更新残差：
\[ \begin{bmatrix} \bar{\phi}_i \\ \phi_{i+1} \end{bmatrix} = \begin{bmatrix} c_i & s_i \\ -s_i & c_i \end{bmatrix} \begin{bmatrix} \phi_i \\ 0 \end{bmatrix} \]
5. 完整LSMR算法
算法伪代码：
输入: A, b, 最大迭代次数k_max，容差tol
输出: 近似解x
1. 初始化:
β₁ = ||b||₂, u₁ = b/β₁
α₁ = ||Aᵀu₁||₂, v₁ = Aᵀu₁/α₁
w₁ = v₁
x = 0
ρ̄₁ = α₁, ϕ̄₁ = β₁
2. for i = 1 to k_max do
// 双对角化
βᵢ₊₁uᵢ₊₁ = Avᵢ - αᵢuᵢ
αᵢ₊₁vᵢ₊₁ = Aᵀuᵢ₊₁ - βᵢ₊₁vᵢ
// 更新QR分解
ρᵢ = sqrt(ρ̄ᵢ² + βᵢ₊₁²)
cᵢ = ρ̄ᵢ/ρᵢ, sᵢ = βᵢ₊₁/ρᵢ
θᵢ₊₁ = sᵢαᵢ₊₁
ρ̄ᵢ₊₁ = -cᵢαᵢ₊₁
ϕᵢ = cᵢϕ̄ᵢ
ϕ̄ᵢ₊₁ = sᵢϕ̄ᵢ
// 更新解
ζᵢ = ϕᵢ/ρᵢ
x = x + ζᵢwᵢ
// 更新搜索方向
wᵢ₊₁ = vᵢ₊₁ - (θᵢ₊₁/ρᵢ)wᵢ
// 检查收敛
if |ϕ̄ᵢ₊₁|/β₁ < tol then
break
end if
// 重新正交化（可选，提高数值稳定性）
for j = 1 to i do
wᵢ₊₁ = wᵢ₊₁ - (wⱼᵀwᵢ₊₁)wⱼ
end for
wᵢ₊₁ = wᵢ₊₁/||wᵢ₊₁||₂
3. end for
6. LSMR的优势与特性
- 数值稳定性：相比LSQR，LSMR能更准确地求解病态问题
- 单调收敛：残差范数 \(\|A^T(Ax_k - b)\|_2\) 单调递减
- 无需额外存储：只需存储最近几个向量
- 自适应正则化：自动处理秩亏问题
7. 收敛性分析
重要性质：
- LSMR等价于MINRES方法应用于正规方程 \(A^T A x = A^T b\)
- 残差满足：
\[ \|A^T(Ax_k - b)\|_2 = |\bar{\phi}_{k+1}| \]
- 误差估计：
\[ \|x_k - x_*\|_2 \leq \frac{|\bar{\phi}_{k+1}|}{\sigma_{min}(A)} \]
其中 \(\sigma_{min}(A)\) 是A的最小奇异值
8. 实际应用中的考虑
预处理技术：
为了加速收敛，常用预处理：
\[\min_{x} \|AM^{-1}y - b\|_2, \quad x = M^{-1}y \]
其中M是预处理子，通常取为对角矩阵或不完全Cholesky分解。
停止准则：
- 相对残差：\(\|A^T(Ax - b)\|_2 / \|A^T b\|_2 < \text{tol}\)
- 迭代次数限制
- 残差不再显著下降
内存优化：
- 使用矩阵-向量乘积而非显式存储A
- 只存储最近的几个向量而非全部基向量
- 实现短递归（3项递推关系）
9. 与相关算法的比较
| 特性 | LSQR | LSMR | CGNR |
|---|---|---|---|
| 最小化目标 | |Ax-b|₂ | |Aᵀ(Ax-b)|₂ | |Aᵀ(Ax-b)|₂ |
| 稳定性 | 好 | 更好 | 较差 |
| 存储需求 | O(m+n) | O(m+n) | O(m+n) |
| 适合问题 | 良态最小二乘 | 病态最小二乘 | 对称正定正规方程 |
10. 应用场景
- 大规模稀疏最小二乘：如CT重建、图像处理
- 正则化问题：Tikhonov正则化的求解
- 矩阵补全：协同过滤推荐系统
- 压缩感知：稀疏信号恢复
- 大地测量：GPS定位等大规模平差问题
LSMR算法将大规模最小二乘问题转化为一系列小型三对角子问题，通过Krylov子空间方法高效求解，兼具数值稳定性和计算效率，是大规模稀疏最小二乘问题的首选算法之一。