基于Krylov子空间的CGS（共轭梯度平方）算法解非对称线性方程组
字数 3881 2025-12-19 14:19:47

基于Krylov子空间的CGS（共轭梯度平方）算法解非对称线性方程组

题目描述
给定一个大型稀疏非对称线性方程组 \(A \mathbf{x} = \mathbf{b}\)，其中 \(A \in \mathbb{R}^{n \times n}\) 是非对称矩阵，\(\mathbf{b} \in \mathbb{R}^n\) 是已知向量。我们希望求解未知向量 \(\mathbf{x} \in \mathbb{R}^n\)。由于矩阵 \(A\) 非对称，无法直接使用共轭梯度法（CG）。CGS（共轭梯度平方）算法是Krylov子空间方法的一种，它通过巧妙的迭代格式，避免了非对称矩阵带来的正交性限制，适用于求解非对称线性方程组。

解题过程循序渐进讲解

步骤1：问题转化与算法思路
CGS算法的核心思想是在Krylov子空间 \(\mathcal{K}_m(A, \mathbf{r}_0) = \text{span}\{\mathbf{r}_0, A\mathbf{r}_0, A^2\mathbf{r}_0, \dots, A^{m-1}\mathbf{r}_0\}\) 中寻找近似解，其中 \(\mathbf{r}_0 = \mathbf{b} - A \mathbf{x}_0\) 是初始残差向量，\(\mathbf{x}_0\) 是初始猜测解。CGS通过多项式加速技术，构造迭代解，使得残差满足 \(\mathbf{r}_m = P_m(A)^2 \mathbf{r}_0\)，其中 \(P_m\) 是某个多项式。这可以显著加速收敛，尤其当矩阵 \(A\) 的特征值分布有利时。

步骤2：算法推导的基础——双共轭梯度法（BiCG）回顾
CGS的推导基于双共轭梯度法（BiCG）。BiCG在构造向量序列时，需要同时处理原方程组和伴随方程组 \(A^T \mathbf{\tilde{x}} = \mathbf{\tilde{b}}\)，并引入“影子残差”向量 \(\mathbf{\tilde{r}}_m\) 来强制双正交性。BiCG的残差满足：

\[\mathbf{r}_m = P_m(A) \mathbf{r}_0, \quad \mathbf{\tilde{r}}_m = P_m(A^T) \mathbf{\tilde{r}}_0 \]

其中 \(P_m\) 是m次多项式。CGS的关键观察是：如果我们令影子残差初始向量 \(\mathbf{\tilde{r}}_0 = \mathbf{r}_0\)，那么从BiCG的迭代关系中可以导出 \(\mathbf{r}_m = P_m(A)^2 \mathbf{r}_0\)，从而避免显式计算伴随矩阵 \(A^T\) 的乘积，提高效率。

步骤3：CGS算法的迭代公式构造
CGS的迭代过程定义以下向量序列：

  • 近似解 \(\mathbf{x}_m\)
  • 残差 \(\mathbf{r}_m = \mathbf{b} - A \mathbf{x}_m\)
  • 辅助向量 \(\mathbf{p}_m\) 和 \(\mathbf{q}_m\)，用于更新。

具体迭代公式（从初始猜测 \(\mathbf{x}_0\) 开始）：

  1. 初始化：\(\mathbf{r}_0 = \mathbf{b} - A \mathbf{x}_0\)，设置 \(\mathbf{\tilde{r}}_0 = \mathbf{r}_0\)（影子残差初始值）。
  2. 令 \(\mathbf{p}_0 = \mathbf{r}_0\)，\(\mathbf{q}_0 = \mathbf{0}\)。
  3. 对于 \(m = 0, 1, 2, \dots\) 直到收敛，执行：
    a. 计算内积：\(\alpha_m = \frac{(\mathbf{r}_m, \mathbf{\tilde{r}}_0)}{(A \mathbf{p}_m, \mathbf{\tilde{r}}_0)}\)
    b. 计算辅助向量：\(\mathbf{q}_{m+1} = \mathbf{r}_m - \alpha_m A \mathbf{p}_m\)
    c. 更新解：\(\mathbf{x}_{m+1} = \mathbf{x}_m + \alpha_m (\mathbf{p}_m + \mathbf{q}_{m+1})\)
    d. 计算新残差：\(\mathbf{r}_{m+1} = \mathbf{r}_m - \alpha_m A (\mathbf{p}_m + \mathbf{q}_{m+1})\)
    e. 计算内积：\(\beta_m = \frac{(\mathbf{r}_{m+1}, \mathbf{\tilde{r}}_0)}{(\mathbf{r}_m, \mathbf{\tilde{r}}_0)}\)
    f. 更新方向向量：\(\mathbf{p}_{m+1} = \mathbf{r}_{m+1} + \beta_m \mathbf{q}_{m+1} + \beta_m^2 \mathbf{p}_m\)

注意：上述公式中，\((\cdot, \cdot)\) 表示向量内积。CGS只需计算矩阵 \(A\) 与向量的乘积，无需 \(A^T\)，计算量比BiCG更小。

步骤4：算法执行的详细步骤
下面我们逐步拆解CGS的一次迭代（以第 \(m\) 步为例）：

  1. 输入：矩阵 \(A\)，右端项 \(\mathbf{b}\)，初始解 \(\mathbf{x}_0\)，最大迭代次数 \(M\)，容差 \(\epsilon\)。
  2. 初始化：
    • 计算初始残差 \(\mathbf{r}_0 = \mathbf{b} - A \mathbf{x}_0\)。
    • 设置影子残差 \(\mathbf{\tilde{r}}_0 = \mathbf{r}_0\)（通常直接赋值，不实际解伴随系统）。
    • 令 \(\mathbf{p}_0 = \mathbf{r}_0\)，\(\mathbf{q}_0 = \mathbf{0}\)。
    • 计算初始残差范数 \(\rho_0 = \|\mathbf{r}_0\|_2\)，若 \(\rho_0 < \epsilon\) 则输出 \(\mathbf{x}_0\) 并终止。
  3. 迭代循环（对于 \(m = 0, 1, \dots, M-1\)）：
    a. 计算矩阵向量积 \(\mathbf{v} = A \mathbf{p}_m\)。
    b. 计算内积 \(\sigma_m = (\mathbf{v}, \mathbf{\tilde{r}}_0)\)。若 \(\sigma_m = 0\)，算法可能中断（通常重启或报错）。
    c. 计算步长 \(\alpha_m = (\mathbf{r}_m, \mathbf{\tilde{r}}_0) / \sigma_m\)。
    d. 计算辅助向量 \(\mathbf{q}_{m+1} = \mathbf{r}_m - \alpha_m \mathbf{v}\)。
    e. 计算矩阵向量积 \(\mathbf{u} = A (\mathbf{p}_m + \mathbf{q}_{m+1})\)。
    f. 更新解 \(\mathbf{x}_{m+1} = \mathbf{x}_m + \alpha_m (\mathbf{p}_m + \mathbf{q}_{m+1})\)。
    g. 更新残差 \(\mathbf{r}_{m+1} = \mathbf{r}_m - \alpha_m \mathbf{u}\)。
    h. 计算残差范数，若 \(\|\mathbf{r}_{m+1}\|_2 < \epsilon\)，输出 \(\mathbf{x}_{m+1}\) 并终止。
    i. 计算内积 \(\beta_m = (\mathbf{r}_{m+1}, \mathbf{\tilde{r}}_0) / (\mathbf{r}_m, \mathbf{\tilde{r}}_0)\)。
    j. 更新方向向量 \(\mathbf{p}_{m+1} = \mathbf{r}_{m+1} + \beta_m \mathbf{q}_{m+1} + \beta_m^2 \mathbf{p}_m\)。
  4. 输出：最终近似解 \(\mathbf{x}_M\)，或未达到容差时的警告。

步骤5：CGS算法的特点与注意事项

  • 优点：CGS完全避免使用 \(A^T\)，适合 \(A^T\) 难以获取或计算代价高的问题；每次迭代仅需两次矩阵向量乘积，计算量较小。
  • 缺点：残差范数可能振荡剧烈，收敛不稳定；可能发生“溢出”或除零错误（当内积接近零时）。
  • 数值稳定性改进：实际中常使用更稳定的变体如BiCGSTAB（双共轭梯度稳定法），它在CGS基础上引入局部最小化，平滑收敛过程。

总结
CGS算法是一种基于Krylov子空间的迭代法，通过平方多项式加速，高效求解非对称线性方程组。其核心在于利用BiCG的思想但不需 \(A^T\)，迭代公式紧凑。尽管存在数值不稳定的风险，但在许多实际问题中，结合适当的预处理技术，CGS能有效求解大规模稀疏非对称系统。

基于Krylov子空间的CGS（共轭梯度平方）算法解非对称线性方程组 题目描述 给定一个大型稀疏非对称线性方程组 \(A \mathbf{x} = \mathbf{b}\)，其中 \(A \in \mathbb{R}^{n \times n}\) 是非对称矩阵，\(\mathbf{b} \in \mathbb{R}^n\) 是已知向量。我们希望求解未知向量 \(\mathbf{x} \in \mathbb{R}^n\)。由于矩阵 \(A\) 非对称，无法直接使用共轭梯度法（CG）。CGS（共轭梯度平方）算法是Krylov子空间方法的一种，它通过巧妙的迭代格式，避免了非对称矩阵带来的正交性限制，适用于求解非对称线性方程组。 解题过程循序渐进讲解 步骤1：问题转化与算法思路 CGS算法的核心思想是在Krylov子空间 \(\mathcal{K}_ m(A, \mathbf{r}_ 0) = \text{span}\{\mathbf{r}_ 0, A\mathbf{r}_ 0, A^2\mathbf{r}_ 0, \dots, A^{m-1}\mathbf{r}_ 0\}\) 中寻找近似解，其中 \(\mathbf{r}_ 0 = \mathbf{b} - A \mathbf{x}_ 0\) 是初始残差向量，\(\mathbf{x}_ 0\) 是初始猜测解。CGS通过多项式加速技术，构造迭代解，使得残差满足 \(\mathbf{r}_ m = P_ m(A)^2 \mathbf{r}_ 0\)，其中 \(P_ m\) 是某个多项式。这可以显著加速收敛，尤其当矩阵 \(A\) 的特征值分布有利时。 步骤2：算法推导的基础——双共轭梯度法（BiCG）回顾 CGS的推导基于双共轭梯度法（BiCG）。BiCG在构造向量序列时，需要同时处理原方程组和伴随方程组 \(A^T \mathbf{\tilde{x}} = \mathbf{\tilde{b}}\)，并引入“影子残差”向量 \(\mathbf{\tilde{r}}_ m\) 来强制双正交性。BiCG的残差满足： \[ \mathbf{r}_ m = P_ m(A) \mathbf{r}_ 0, \quad \mathbf{\tilde{r}}_ m = P_ m(A^T) \mathbf{\tilde{r}}_ 0 \] 其中 \(P_ m\) 是m次多项式。CGS的关键观察是：如果我们令影子残差初始向量 \(\mathbf{\tilde{r}}_ 0 = \mathbf{r}_ 0\)，那么从BiCG的迭代关系中可以导出 \(\mathbf{r}_ m = P_ m(A)^2 \mathbf{r}_ 0\)，从而避免显式计算伴随矩阵 \(A^T\) 的乘积，提高效率。 步骤3：CGS算法的迭代公式构造 CGS的迭代过程定义以下向量序列： 近似解 \(\mathbf{x}_ m\) 残差 \(\mathbf{r}_ m = \mathbf{b} - A \mathbf{x}_ m\) 辅助向量 \(\mathbf{p}_ m\) 和 \(\mathbf{q}_ m\)，用于更新。 具体迭代公式（从初始猜测 \(\mathbf{x}_ 0\) 开始）： 初始化：\(\mathbf{r}_ 0 = \mathbf{b} - A \mathbf{x}_ 0\)，设置 \(\mathbf{\tilde{r}}_ 0 = \mathbf{r}_ 0\)（影子残差初始值）。 令 \(\mathbf{p}_ 0 = \mathbf{r}_ 0\)，\(\mathbf{q}_ 0 = \mathbf{0}\)。 对于 \(m = 0, 1, 2, \dots\) 直到收敛，执行： a. 计算内积：\(\alpha_ m = \frac{(\mathbf{r}_ m, \mathbf{\tilde{r}}_ 0)}{(A \mathbf{p}_ m, \mathbf{\tilde{r}} 0)}\) b. 计算辅助向量：\(\mathbf{q} {m+1} = \mathbf{r}_ m - \alpha_ m A \mathbf{p} m\) c. 更新解：\(\mathbf{x} {m+1} = \mathbf{x} m + \alpha_ m (\mathbf{p} m + \mathbf{q} {m+1})\) d. 计算新残差：\(\mathbf{r} {m+1} = \mathbf{r} m - \alpha_ m A (\mathbf{p} m + \mathbf{q} {m+1})\) e. 计算内积：\(\beta_ m = \frac{(\mathbf{r} {m+1}, \mathbf{\tilde{r}} 0)}{(\mathbf{r} m, \mathbf{\tilde{r}} 0)}\) f. 更新方向向量：\(\mathbf{p} {m+1} = \mathbf{r} {m+1} + \beta_ m \mathbf{q} {m+1} + \beta_ m^2 \mathbf{p}_ m\) 注意：上述公式中，\((\cdot, \cdot)\) 表示向量内积。CGS只需计算矩阵 \(A\) 与向量的乘积，无需 \(A^T\)，计算量比BiCG更小。 步骤4：算法执行的详细步骤 下面我们逐步拆解CGS的一次迭代（以第 \(m\) 步为例）： 输入 ：矩阵 \(A\)，右端项 \(\mathbf{b}\)，初始解 \(\mathbf{x}_ 0\)，最大迭代次数 \(M\)，容差 \(\epsilon\)。 初始化 ： 计算初始残差 \(\mathbf{r}_ 0 = \mathbf{b} - A \mathbf{x}_ 0\)。 设置影子残差 \(\mathbf{\tilde{r}}_ 0 = \mathbf{r}_ 0\)（通常直接赋值，不实际解伴随系统）。 令 \(\mathbf{p}_ 0 = \mathbf{r}_ 0\)，\(\mathbf{q}_ 0 = \mathbf{0}\)。 计算初始残差范数 \(\rho_ 0 = \|\mathbf{r}_ 0\|_ 2\)，若 \(\rho_ 0 < \epsilon\) 则输出 \(\mathbf{x}_ 0\) 并终止。 迭代循环 （对于 \(m = 0, 1, \dots, M-1\)）： a. 计算矩阵向量积 \(\mathbf{v} = A \mathbf{p}_ m\)。 b. 计算内积 \(\sigma_ m = (\mathbf{v}, \mathbf{\tilde{r}}_ 0)\)。若 \(\sigma_ m = 0\)，算法可能中断（通常重启或报错）。 c. 计算步长 \(\alpha_ m = (\mathbf{r} m, \mathbf{\tilde{r}} 0) / \sigma_ m\)。 d. 计算辅助向量 \(\mathbf{q} {m+1} = \mathbf{r} m - \alpha_ m \mathbf{v}\)。 e. 计算矩阵向量积 \(\mathbf{u} = A (\mathbf{p} m + \mathbf{q} {m+1})\)。 f. 更新解 \(\mathbf{x} {m+1} = \mathbf{x} m + \alpha_ m (\mathbf{p} m + \mathbf{q} {m+1})\)。 g. 更新残差 \(\mathbf{r} {m+1} = \mathbf{r} m - \alpha_ m \mathbf{u}\)。 h. 计算残差范数，若 \(\|\mathbf{r} {m+1}\| 2 < \epsilon\)，输出 \(\mathbf{x} {m+1}\) 并终止。 i. 计算内积 \(\beta_ m = (\mathbf{r} {m+1}, \mathbf{\tilde{r}} 0) / (\mathbf{r} m, \mathbf{\tilde{r}} 0)\)。 j. 更新方向向量 \(\mathbf{p} {m+1} = \mathbf{r} {m+1} + \beta_ m \mathbf{q} {m+1} + \beta_ m^2 \mathbf{p}_ m\)。 输出 ：最终近似解 \(\mathbf{x}_ M\)，或未达到容差时的警告。 步骤5：CGS算法的特点与注意事项 优点 ：CGS完全避免使用 \(A^T\)，适合 \(A^T\) 难以获取或计算代价高的问题；每次迭代仅需两次矩阵向量乘积，计算量较小。 缺点 ：残差范数可能振荡剧烈，收敛不稳定；可能发生“溢出”或除零错误（当内积接近零时）。 数值稳定性改进 ：实际中常使用更稳定的变体如BiCGSTAB（双共轭梯度稳定法），它在CGS基础上引入局部最小化，平滑收敛过程。 总结 CGS算法是一种基于Krylov子空间的迭代法，通过平方多项式加速，高效求解非对称线性方程组。其核心在于利用BiCG的思想但不需 \(A^T\)，迭代公式紧凑。尽管存在数值不稳定的风险，但在许多实际问题中，结合适当的预处理技术，CGS能有效求解大规模稀疏非对称系统。