基于Krylov子空间的LSMR算法求解线性最小二乘问题
字数 3218 2025-12-15 03:55:18

基于Krylov子空间的LSMR算法求解线性最小二乘问题

问题描述
LSMR算法是一种迭代算法，用于求解线性最小二乘问题 \(\min_x \|Ax - b\|_2\)，其中 \(A \in \mathbb{R}^{m \times n}\)（通常 \(m \geq n\)），\(b \in \mathbb{R}^{m}\)。它基于Krylov子空间方法，通过构造Krylov子空间 \(\mathcal{K}_k(A^\top A, A^\top b)\) 来逐步逼近最小二乘解。LSMR与经典的LSQR算法类似，但通过隐式构造残差范数最小化问题，往往在数值稳定性和收敛性上更具优势，尤其适用于病态问题。

解题过程循序渐进讲解

步骤1：问题转化与核心思想

线性最小二乘问题等价于求解正规方程 \(A^\top A x = A^\top b\)。
直接求解正规方程可能面临条件数平方放大的问题（\(\kappa(A^\top A) = \kappa(A)^2\)），导致数值不稳定。LSMR的核心思想是：

  • 利用Krylov子空间方法，避免显式计算 \(A^\top A\)。
  • 通过Golub-Kahan双对角化过程将 \(A\) 转化为双对角矩阵，从而将原问题转化为一个小规模的子问题。
  • 在子问题中隐式最小化残差范数 \(\|A x_k - b\|_2\) 和正规方程残差范数 \(\|A^\top (A x_k - b)\|_2\)。

步骤2：Golub-Kahan双对角化过程

双对角化过程是LSMR的基础，它通过迭代构造两个正交矩阵 \(U_k\) 和 \(V_k\)，将 \(A\) 转换为双对角矩阵 \(B_k\)：

  1. 初始化：
    设 \(\beta_1 u_1 = b\)（其中 \(\beta_1 = \|b\|_2\), \(u_1\) 为单位向量）。
    设 \(\alpha_1 v_1 = A^\top u_1\)（其中 \(\alpha_1 = \|A^\top u_1\|_2\), \(v_1\) 为单位向量）。

  2. 迭代步骤（对于 \(i=1,2,\dots\)）：

    • 计算 \(\beta_{i+1} u_{i+1} = A v_i - \alpha_i u_i\)（确保 \(u_{i+1}\) 与之前的 \(u_j\) 正交）。
    • 计算 \(\alpha_{i+1} v_{i+1} = A^\top u_{i+1} - \beta_{i+1} v_i\)（确保 \(v_{i+1}\) 与之前的 \(v_j\) 正交）。

    经过 \(k\) 步后，得到关系：

\[ A V_k = U_{k+1} B_k, \quad A^\top U_{k+1} = V_k B_k^\top + \alpha_{k+1} v_{k+1} e_{k+1}^\top \]

其中 \(B_k\) 是 \((k+1) \times k\) 的双对角矩阵：

\[ B_k = \begin{bmatrix} \alpha_1 & & & \\ \beta_2 & \alpha_2 & & \\ & \beta_3 & \ddots & \\ & & \ddots & \alpha_k \\ & & & \beta_{k+1} \end{bmatrix}. \]

步骤3：构造子问题

将解 \(x\) 表示为 \(x = V_k y_k\)（其中 \(y_k \in \mathbb{R}^k\)）。
原始最小二乘问题的残差可写为：

\[\|A x - b\|_2 = \|A V_k y_k - b\|_2 = \|U_{k+1} B_k y_k - \beta_1 u_1\|_2. \]

由于 \(U_{k+1}\) 是列正交矩阵，最小化 \(\|A x - b\|_2\) 等价于最小化 \(\|B_k y_k - \beta_1 e_1\|_2\)，其中 \(e_1 = (1,0,\dots,0)^\top \in \mathbb{R}^{k+1}\)。
这样，原问题被转化为一个小规模的子问题：

\[\min_{y_k} \|B_k y_k - \beta_1 e_1\|_2. \]

步骤4：LSMR的隐式求解策略

LSMR并不直接求解上述子问题，而是通过隐式方式逐步更新近似解 \(x_k\)。具体过程：

  1. 利用QR分解处理子问题：
    对双对角矩阵 \(B_k\) 进行QR分解：\(B_k = Q_k R_k\)，其中 \(R_k\) 是上三角矩阵。
    子问题变为 \(\min \|R_k y_k - Q_k^\top (\beta_1 e_1)\|_2\)。

  2. 递推更新：
    每一步迭代中，通过Givens旋转更新 \(R_k\) 和右端项，避免存储整个矩阵。
    设当前近似解 \(x_k = V_k y_k\)，正规方程残差 \(r_k = A^\top (A x_k - b)\)。
    LSMR的关键在于：通过双对角化过程，可以递推计算 \(\|r_k\|_2\) 并确保其单调递减。

  3. 终止条件：
    当 \(\|r_k\|_2\) 小于给定阈值，或达到最大迭代次数时停止。
    由于 \(r_k\) 与 \(A^\top (A x_k - b)\) 成比例，这保证了近似解满足正规方程的残差足够小。

步骤5：算法伪代码

  1. 初始化：
    \(\beta_1 = \|b\|_2\), \(u_1 = b/\beta_1\)
    \(\alpha_1 = \|A^\top u_1\|_2\), \(v_1 = A^\top u_1 / \alpha_1\)
    \(w_1 = v_1\), \(x_0 = 0\), \(\bar{\rho}_1 = \beta_1\), \(\bar{\phi}_1 = \alpha_1 \beta_1\)

  2. 对于 \(k=1,2,\dots\) 直到收敛：

    • 双对角化：
      \(\beta_{k+1} u_{k+1} = A v_k - \alpha_k u_k\)
      \(\alpha_{k+1} v_{k+1} = A^\top u_{k+1} - \beta_{k+1} v_k\)
    • 构造并更新旋转矩阵，计算标量 \(\rho_k, \theta_{k+1}, \phi_k\) 等。
    • 更新解：\(x_k = x_{k-1} + (\phi_k / \rho_k) w_k\)
    • 更新辅助向量：\(w_{k+1} = v_{k+1} - (\theta_{k+1} / \rho_k) w_k\)
    • 检查 \(\|A^\top (A x_k - b)\|_2\) 是否小于阈值。

步骤6：算法特点与优势

  • 数值稳定性：隐式残差最小化避免了正规方程的条件数放大问题。
  • 适用性：尤其适合大型稀疏矩阵、病态问题，或需要提前终止迭代的场景。
  • 与LSQR的关系：LSMR可视为LSQR的改进版本，通过更精细的残差控制，通常能更快降低正规方程残差。
  • 计算开销：每步迭代需两次矩阵-向量乘（\(A v_k\) 和 \(A^\top u_{k+1}\)），内存占用为 \(O(m+n)\)。

总结
LSMR算法通过Golub-Kahan双对角化将最小二乘问题投影到Krylov子空间，并利用隐式QR分解递推求解，在保证数值稳定性的同时高效逼近最优解。该算法在稀疏矩阵、病态问题及大规模计算中具有重要应用价值。

基于Krylov子空间的LSMR算法求解线性最小二乘问题 问题描述 LSMR算法是一种迭代算法，用于求解线性最小二乘问题 \(\min_ x \|Ax - b\|_ 2\)，其中 \(A \in \mathbb{R}^{m \times n}\)（通常 \(m \geq n\)），\(b \in \mathbb{R}^{m}\)。它基于Krylov子空间方法，通过构造Krylov子空间 \(\mathcal{K}_ k(A^\top A, A^\top b)\) 来逐步逼近最小二乘解。LSMR与经典的LSQR算法类似，但通过隐式构造残差范数最小化问题，往往在数值稳定性和收敛性上更具优势，尤其适用于病态问题。 解题过程循序渐进讲解 步骤1：问题转化与核心思想 线性最小二乘问题等价于求解正规方程 \(A^\top A x = A^\top b\)。 直接求解正规方程可能面临条件数平方放大的问题（\(\kappa(A^\top A) = \kappa(A)^2\)），导致数值不稳定。LSMR的核心思想是： 利用Krylov子空间方法，避免显式计算 \(A^\top A\)。 通过Golub-Kahan双对角化过程将 \(A\) 转化为双对角矩阵，从而将原问题转化为一个小规模的子问题。 在子问题中隐式最小化残差范数 \(\|A x_ k - b\|_ 2\) 和正规方程残差范数 \(\|A^\top (A x_ k - b)\|_ 2\)。 步骤2：Golub-Kahan双对角化过程 双对角化过程是LSMR的基础，它通过迭代构造两个正交矩阵 \(U_ k\) 和 \(V_ k\)，将 \(A\) 转换为双对角矩阵 \(B_ k\)： 初始化 ： 设 \(\beta_ 1 u_ 1 = b\)（其中 \(\beta_ 1 = \|b\|_ 2\), \(u_ 1\) 为单位向量）。 设 \(\alpha_ 1 v_ 1 = A^\top u_ 1\)（其中 \(\alpha_ 1 = \|A^\top u_ 1\|_ 2\), \(v_ 1\) 为单位向量）。 迭代步骤 （对于 \(i=1,2,\dots\)）： 计算 \(\beta_ {i+1} u_ {i+1} = A v_ i - \alpha_ i u_ i\)（确保 \(u_ {i+1}\) 与之前的 \(u_ j\) 正交）。 计算 \(\alpha_ {i+1} v_ {i+1} = A^\top u_ {i+1} - \beta_ {i+1} v_ i\)（确保 \(v_ {i+1}\) 与之前的 \(v_ j\) 正交）。 经过 \(k\) 步后，得到关系： \[ A V_ k = U_ {k+1} B_ k, \quad A^\top U_ {k+1} = V_ k B_ k^\top + \alpha_ {k+1} v_ {k+1} e_ {k+1}^\top \] 其中 \(B_ k\) 是 \((k+1) \times k\) 的双对角矩阵： \[ B_ k = \begin{bmatrix} \alpha_ 1 & & & \\ \beta_ 2 & \alpha_ 2 & & \\ & \beta_ 3 & \ddots & \\ & & \ddots & \alpha_ k \\ & & & \beta_ {k+1} \end{bmatrix}. \] 步骤3：构造子问题 将解 \(x\) 表示为 \(x = V_ k y_ k\)（其中 \(y_ k \in \mathbb{R}^k\)）。 原始最小二乘问题的残差可写为： \[ \|A x - b\|_ 2 = \|A V_ k y_ k - b\| 2 = \|U {k+1} B_ k y_ k - \beta_ 1 u_ 1\| 2. \] 由于 \(U {k+1}\) 是列正交矩阵，最小化 \(\|A x - b\|_ 2\) 等价于最小化 \(\|B_ k y_ k - \beta_ 1 e_ 1\| 2\)，其中 \(e_ 1 = (1,0,\dots,0)^\top \in \mathbb{R}^{k+1}\)。 这样，原问题被转化为一个小规模的子问题： \[ \min {y_ k} \|B_ k y_ k - \beta_ 1 e_ 1\|_ 2. \] 步骤4：LSMR的隐式求解策略 LSMR并不直接求解上述子问题，而是通过隐式方式逐步更新近似解 \(x_ k\)。具体过程： 利用QR分解处理子问题 ： 对双对角矩阵 \(B_ k\) 进行QR分解：\(B_ k = Q_ k R_ k\)，其中 \(R_ k\) 是上三角矩阵。 子问题变为 \(\min \|R_ k y_ k - Q_ k^\top (\beta_ 1 e_ 1)\|_ 2\)。 递推更新 ： 每一步迭代中，通过Givens旋转更新 \(R_ k\) 和右端项，避免存储整个矩阵。 设当前近似解 \(x_ k = V_ k y_ k\)，正规方程残差 \(r_ k = A^\top (A x_ k - b)\)。 LSMR的关键在于：通过双对角化过程，可以递推计算 \(\|r_ k\|_ 2\) 并确保其单调递减。 终止条件 ： 当 \(\|r_ k\|_ 2\) 小于给定阈值，或达到最大迭代次数时停止。 由于 \(r_ k\) 与 \(A^\top (A x_ k - b)\) 成比例，这保证了近似解满足正规方程的残差足够小。 步骤5：算法伪代码 初始化： \(\beta_ 1 = \|b\|_ 2\), \(u_ 1 = b/\beta_ 1\) \(\alpha_ 1 = \|A^\top u_ 1\|_ 2\), \(v_ 1 = A^\top u_ 1 / \alpha_ 1\) \(w_ 1 = v_ 1\), \(x_ 0 = 0\), \(\bar{\rho}_ 1 = \beta_ 1\), \(\bar{\phi}_ 1 = \alpha_ 1 \beta_ 1\) 对于 \(k=1,2,\dots\) 直到收敛： 双对角化： \(\beta_ {k+1} u_ {k+1} = A v_ k - \alpha_ k u_ k\) \(\alpha_ {k+1} v_ {k+1} = A^\top u_ {k+1} - \beta_ {k+1} v_ k\) 构造并更新旋转矩阵，计算标量 \(\rho_ k, \theta_ {k+1}, \phi_ k\) 等。 更新解：\(x_ k = x_ {k-1} + (\phi_ k / \rho_ k) w_ k\) 更新辅助向量：\(w_ {k+1} = v_ {k+1} - (\theta_ {k+1} / \rho_ k) w_ k\) 检查 \(\|A^\top (A x_ k - b)\|_ 2\) 是否小于阈值。 步骤6：算法特点与优势 数值稳定性 ：隐式残差最小化避免了正规方程的条件数放大问题。 适用性 ：尤其适合大型稀疏矩阵、病态问题，或需要提前终止迭代的场景。 与LSQR的关系 ：LSMR可视为LSQR的改进版本，通过更精细的残差控制，通常能更快降低正规方程残差。 计算开销 ：每步迭代需两次矩阵-向量乘（\(A v_ k\) 和 \(A^\top u_ {k+1}\)），内存占用为 \(O(m+n)\)。 总结 LSMR算法通过Golub-Kahan双对角化将最小二乘问题投影到Krylov子空间，并利用隐式QR分解递推求解，在保证数值稳定性的同时高效逼近最优解。该算法在稀疏矩阵、病态问题及大规模计算中具有重要应用价值。