基于Krylov子空间的CGNR算法解最小二乘问题
字数 2977 2025-12-13 16:08:19

基于Krylov子空间的CGNR算法解最小二乘问题

题目描述

考虑线性最小二乘问题：给定矩阵 \(A \in \mathbb{R}^{m \times n}\)（通常 \(m \ge n\)）和向量 \(b \in \mathbb{R}^{m}\)，寻找 \(x \in \mathbb{R}^{n}\) 使得残差范数 \(\| b - A x \|_2\) 最小化。当 \(A\) 是大型稀疏矩阵时，直接法（如正规方程法或SVD）可能计算代价过高。本题目要求解释基于Krylov子空间的共轭梯度法应用于正规方程（Conjugate Gradient on the Normal Equations, CGNR） 的算法原理、推导过程和迭代步骤，并分析其收敛性和适用场景。

解题过程

第一步：问题转化为正规方程

最小二乘问题 \(\min_x \| b - A x \|_2^2\) 等价于求解正规方程：

\[A^T A x = A^T b \]

这是因为目标函数 \(f(x) = \| b - A x \|_2^2\) 的梯度为 \(\nabla f(x) = 2A^T (A x - b)\)，令梯度为零即得正规方程。注意矩阵 \(C = A^T A\) 是对称半正定的，当 \(A\) 列满秩（即 \(\text{rank}(A) = n\)）时，\(C\) 对称正定。

第二步：共轭梯度法（CG）的适用性

共轭梯度法（CG）是求解对称正定线性方程组 \(C x = d\)（其中 \(d = A^T b\)）的经典迭代算法。然而，直接对 \(C = A^T A\) 应用CG存在两个潜在问题：

  • \(C\) 的条件数为 \(\kappa(A^T A) = \kappa(A)^2\)，当 \(A\) 病态时，\(C\) 的条件数平方放大，可能导致CG收敛缓慢。
  • 计算 \(C\) 的矩阵-向量乘积需要两次矩阵-向量乘：一次是 \(A p\) 一次是 \(A^T (A p)\)，每次迭代代价稍高。

尽管如此，CGNR 仍然是一种常用算法，因为它避免了直接存储 \(C\)（对稀疏 \(A\) 节省内存），且当 \(A\) 相对良态时表现良好。

第三步：CGNR 算法推导

设初始猜测解为 \(x_0\)，对应初始残差 \(r_0 = A^T (b - A x_0) = A^T (b_0)\)，其中 \(b_0 = b - A x_0\)。定义初始搜索方向 \(p_0 = r_0\)。

对于第 \(k\) 次迭代（\(k = 0, 1, 2, \dots\)），标准CG步骤应用到方程 \((A^T A) x = A^T b\) 如下：

  1. 计算步长：

\[\alpha_k = \frac{r_k^T r_k}{p_k^T (A^T A) p_k} = \frac{\| r_k \|_2^2}{\| A p_k \|_2^2} \]

分母中 \(p_k^T (A^T A) p_k = (A p_k)^T (A p_k) = \| A p_k \|_2^2\)。

  1. 更新解：

\[x_{k+1} = x_k + \alpha_k p_k \]

  1. 更新残差：

\[r_{k+1} = r_k - \alpha_k (A^T A) p_k \]

注意 \(r_{k+1} = A^T (b - A x_{k+1})\)，这是正规方程的残差。

  1. 计算系数：

\[\beta_k = \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k} \]

  1. 更新搜索方向：

\[p_{k+1} = r_{k+1} + \beta_k p_k \]

迭代直到残差范数 \(\| r_k \|_2\) 或 \(\| b - A x_k \|_2\) 小于给定容差。

第四步：CGNR 的实用迭代格式

在实现中，我们通常避免显式计算 \(A^T A\)，而是通过中间变量简化运算。具体步骤如下：

  1. 初始化：

    • 输入 \(A, b\)，初始解 \(x_0\)（常取零向量）。
    • 计算初始残差 \(s_0 = b - A x_0\) 和 \(r_0 = A^T s_0\)。
    • 设置初始搜索方向 \(p_0 = r_0\)。
    • 设定容差 \(\epsilon\) 和最大迭代次数。

  2. 对 \(k = 0, 1, 2, \dots\) 直到收敛：
    a. 计算矩阵-向量乘积 \(q_k = A p_k\)。
    b. 计算步长 \(\alpha_k = \frac{\| r_k \|_2^2}{\| q_k \|_2^2}\)。
    c. 更新解：\(x_{k+1} = x_k + \alpha_k p_k\)。
    d. 更新残差向量：\(s_{k+1} = s_k - \alpha_k q_k\)（这是原始残差 \(b - A x_{k+1}\)）。
    e. 计算正规方程残差：\(r_{k+1} = A^T s_{k+1}\)。
    f. 计算系数 \(\beta_k = \frac{\| r_{k+1} \|_2^2}{\| r_k \|_2^2}\)。
    g. 更新搜索方向：\(p_{k+1} = r_{k+1} + \beta_k p_k\)。
    h. 检查收敛条件，例如 \(\| r_{k+1} \|_2 < \epsilon\) 或 \(\| s_{k+1} \|_2 < \epsilon\)。

第五步：收敛性与适用性分析

  • 收敛速度：CGNR 的收敛速率依赖于 \(A^T A\) 的特征值分布。由于 \(\kappa(A^T A) = \kappa(A)^2\)，当 \(A\) 条件数大时，收敛可能很慢。预处理技术（如不完全QR预处理）常用于加速。
  • 适用场景：
    • \(A\) 是大型稀疏矩阵，且 \(A^T A\) 无法显式存储。
    • 问题中等规模，且 \(A\) 相对良态。
    • 当 \(A\) 是列满秩时，CGNR 可得到唯一解；若秩亏损，则收敛到最小范数最小二乘解（假设初始猜测为零）。
  • 对比其他方法：
    • 与LSQR算法相比，CGNR 是LSQR的一种变体（数学等价于LSQR在精确算术下），但LSQR在数值上更稳定，能更好地处理病态问题。
    • 与直接法（如QR分解）相比，CGNR 适合大规模问题，内存需求低。

总结

CGNR 将最小二乘问题转化为对称正定线性方程组，利用共轭梯度法迭代求解，通过矩阵-向量乘积避免构造 \(A^T A\)。它在处理大型稀疏最小二乘问题时是一种有效工具，但需注意条件数平方放大可能影响收敛速度，可通过预处理技术改善。理解CGNR 的推导和实现细节，有助于掌握Krylov子空间方法在最小二乘问题中的灵活应用。

基于Krylov子空间的CGNR算法解最小二乘问题 题目描述 考虑线性最小二乘问题：给定矩阵 \( A \in \mathbb{R}^{m \times n} \)（通常 \( m \ge n \)）和向量 \( b \in \mathbb{R}^{m} \)，寻找 \( x \in \mathbb{R}^{n} \) 使得残差范数 \( \| b - A x \|_ 2 \) 最小化。当 \( A \) 是大型稀疏矩阵时，直接法（如正规方程法或SVD）可能计算代价过高。本题目要求解释 基于Krylov子空间的共轭梯度法应用于正规方程（Conjugate Gradient on the Normal Equations, CGNR） 的算法原理、推导过程和迭代步骤，并分析其收敛性和适用场景。 解题过程 第一步：问题转化为正规方程 最小二乘问题 \( \min_ x \| b - A x \|_ 2^2 \) 等价于求解正规方程： \[ A^T A x = A^T b \] 这是因为目标函数 \( f(x) = \| b - A x \|_ 2^2 \) 的梯度为 \( \nabla f(x) = 2A^T (A x - b) \)，令梯度为零即得正规方程。注意矩阵 \( C = A^T A \) 是对称半正定的，当 \( A \) 列满秩（即 \( \text{rank}(A) = n \)）时，\( C \) 对称正定。 第二步：共轭梯度法（CG）的适用性 共轭梯度法（CG）是求解对称正定线性方程组 \( C x = d \)（其中 \( d = A^T b \)）的经典迭代算法。然而，直接对 \( C = A^T A \) 应用CG存在两个潜在问题： \( C \) 的条件数为 \( \kappa(A^T A) = \kappa(A)^2 \)，当 \( A \) 病态时，\( C \) 的条件数平方放大，可能导致CG收敛缓慢。 计算 \( C \) 的矩阵-向量乘积需要两次矩阵-向量乘：一次是 \( A p \) 一次是 \( A^T (A p) \)，每次迭代代价稍高。 尽管如此，CGNR 仍然是一种常用算法，因为它避免了直接存储 \( C \)（对稀疏 \( A \) 节省内存），且当 \( A \) 相对良态时表现良好。 第三步：CGNR 算法推导 设初始猜测解为 \( x_ 0 \)，对应初始残差 \( r_ 0 = A^T (b - A x_ 0) = A^T (b_ 0) \)，其中 \( b_ 0 = b - A x_ 0 \)。定义初始搜索方向 \( p_ 0 = r_ 0 \)。 对于第 \( k \) 次迭代（\( k = 0, 1, 2, \dots \)），标准CG步骤应用到方程 \( (A^T A) x = A^T b \) 如下： 计算步长： \[ \alpha_ k = \frac{r_ k^T r_ k}{p_ k^T (A^T A) p_ k} = \frac{\| r_ k \|_ 2^2}{\| A p_ k \|_ 2^2} \] 分母中 \( p_ k^T (A^T A) p_ k = (A p_ k)^T (A p_ k) = \| A p_ k \|_ 2^2 \)。 更新解： \[ x_ {k+1} = x_ k + \alpha_ k p_ k \] 更新残差： \[ r_ {k+1} = r_ k - \alpha_ k (A^T A) p_ k \] 注意 \( r_ {k+1} = A^T (b - A x_ {k+1}) \)，这是正规方程的残差。 计算系数： \[ \beta_ k = \frac{r_ {k+1}^T r_ {k+1}}{r_ k^T r_ k} \] 更新搜索方向： \[ p_ {k+1} = r_ {k+1} + \beta_ k p_ k \] 迭代直到残差范数 \( \| r_ k \|_ 2 \) 或 \( \| b - A x_ k \|_ 2 \) 小于给定容差。 第四步：CGNR 的实用迭代格式 在实现中，我们通常避免显式计算 \( A^T A \)，而是通过中间变量简化运算。具体步骤如下： 初始化： 输入 \( A, b \)，初始解 \( x_ 0 \)（常取零向量）。 计算初始残差 \( s_ 0 = b - A x_ 0 \) 和 \( r_ 0 = A^T s_ 0 \)。 设置初始搜索方向 \( p_ 0 = r_ 0 \)。 设定容差 \( \epsilon \) 和最大迭代次数。 对 \( k = 0, 1, 2, \dots \) 直到收敛： a. 计算矩阵-向量乘积 \( q_ k = A p_ k \)。 b. 计算步长 \( \alpha_ k = \frac{\| r_ k \| 2^2}{\| q_ k \| 2^2} \)。 c. 更新解：\( x {k+1} = x_ k + \alpha_ k p_ k \)。 d. 更新残差向量：\( s {k+1} = s_ k - \alpha_ k q_ k \)（这是原始残差 \( b - A x_ {k+1} \)）。 e. 计算正规方程残差：\( r_ {k+1} = A^T s_ {k+1} \)。 f. 计算系数 \( \beta_ k = \frac{\| r_ {k+1} \| 2^2}{\| r_ k \| 2^2} \)。 g. 更新搜索方向：\( p {k+1} = r {k+1} + \beta_ k p_ k \)。 h. 检查收敛条件，例如 \( \| r_ {k+1} \| 2 < \epsilon \) 或 \( \| s {k+1} \|_ 2 < \epsilon \)。 第五步：收敛性与适用性分析 收敛速度 ：CGNR 的收敛速率依赖于 \( A^T A \) 的特征值分布。由于 \( \kappa(A^T A) = \kappa(A)^2 \)，当 \( A \) 条件数大时，收敛可能很慢。预处理技术（如不完全QR预处理）常用于加速。 适用场景 ： \( A \) 是大型稀疏矩阵，且 \( A^T A \) 无法显式存储。 问题中等规模，且 \( A \) 相对良态。 当 \( A \) 是列满秩时，CGNR 可得到唯一解；若秩亏损，则收敛到最小范数最小二乘解（假设初始猜测为零）。 对比其他方法 ： 与LSQR算法相比，CGNR 是LSQR的一种变体（数学等价于LSQR在精确算术下），但LSQR在数值上更稳定，能更好地处理病态问题。 与直接法（如QR分解）相比，CGNR 适合大规模问题，内存需求低。 总结 CGNR 将最小二乘问题转化为对称正定线性方程组，利用共轭梯度法迭代求解，通过矩阵-向量乘积避免构造 \( A^T A \)。它在处理大型稀疏最小二乘问题时是一种有效工具，但需注意条件数平方放大可能影响收敛速度，可通过预处理技术改善。理解CGNR 的推导和实现细节，有助于掌握Krylov子空间方法在最小二乘问题中的灵活应用。