基于Krylov子空间的CGS算法解非对称线性方程组
字数 2859 2025-11-07 12:33:00

基于Krylov子空间的CGS算法解非对称线性方程组

我将为你讲解共轭梯度平方法（Conjugate Gradient Squared, CGS），这是一种用于求解非对称线性方程组的Krylov子空间方法。我们将从一个具体问题开始。

问题描述
考虑一个非对称的线性方程组 \(A\vec{x} = \vec{b}\)，其中：

  • \(A\) 是一个 \(n \times n\) 的非对称矩阵（即 \(A \neq A^T\)）。
  • \(\vec{b}\) 是一个已知的 \(n \times 1\) 向量。
  • \(\vec{x}\) 是我们需要求解的 \(n \times 1\) 未知向量。

我们的目标是找到一个向量 \(\vec{x}\)，使得方程成立。当矩阵 \(A\) 规模很大且稀疏时，直接解法（如LU分解）计算成本高昂，迭代法如CGS则更为高效。

解题过程

CGS算法是双共轭梯度法（BiCG）的一种变体，它通过平方化的思想来避免BiCG中与 \(A^T\) 相关的运算，从而在某些情况下能更快地收敛。其核心思想是在Krylov子空间 \(\mathcal{K}_m(A; \vec{r}_0) = \text{span}\{\vec{r}_0, A\vec{r}_0, A^2\vec{r}_0, \dots, A^{m-1}\vec{r}_0\}\) 中寻找近似解，其中 \(\vec{r}_0\) 是初始残差。

第一步：算法初始化
算法从一个初始猜测解 \(\vec{x}_0\) 开始。通常，如果没有任何先验知识，可以设 \(\vec{x}_0 = \vec{0}\)。

  1. 计算初始残差：\(\vec{r}_0 = \vec{b} - A\vec{x}_0\)。
  2. 选择一个辅助向量 \(\vec{\tilde{r}}_0\)，它必须与 \(\vec{r}_0\) 非正交（即内积不为零）。最简单且常见的选择是令 \(\vec{\tilde{r}}_0 = \vec{r}_0\)。
  3. 初始化一些中间向量：令 \(\vec{p}_0 = \vec{r}_0\)，\(\vec{q}_0 = \vec{0}\)。
  4. 设置初始参数：\(\rho_0 = \alpha_0 = \omega_0 = 1\)。
  5. 开始迭代，对于 \(k = 0, 1, 2, \dots\)，直到残差 \(\vec{r}_k\) 的范数足够小为止。

第二步：迭代步骤（第 k 步）
CGS的迭代过程涉及多个向量的更新，我们一步步来看。

  1. 计算迭代参数 \(\rho_k\) 和 \(\beta_k\)：

    • \(\rho_k = (\vec{\tilde{r}}_0, \vec{r}_k)\)（计算辅助残差与当前残差的内积）。
    • 如果 \(\rho_k\) 接近于0，算法可能会失败（称为“breakdown”）。
    • \(\beta_k = (\rho_k / \rho_{k-1}) \times (\alpha_{k-1} / \omega_{k-1})\)。这里用到了上一步的参数。

  2. 更新搜索方向向量 \(\vec{u}_k\) 和 \(\vec{p}_k\)：

    • \(\vec{u}_k = \vec{r}_k + \beta_k \times \omega_{k-1} \times (\vec{u}_{k-1} - \omega_{k-1} A \vec{p}_{k-1})\)（对于k=0，相关历史向量视为0）。
    • \(\vec{p}_k = \vec{u}_k + \beta_k \times (\vec{p}_{k-1} - \omega_{k-1} A \vec{p}_{k-1})\)。
    • 这些更新公式看似复杂，但其目的是构造一个在某种意义下“共轭”的搜索方向，其推导源于对BiCG算法的多项式平方化。

  3. 进行矩阵向量乘法并计算步长 \(\alpha_k\)：

    • 计算 \(\vec{v}_k = A \vec{p}_k\)。
    • 计算步长：\(\alpha_k = \rho_k / (\vec{\tilde{r}}_0, \vec{v}_k)\)。分母是辅助残差与 \(\vec{v}_k\) 的内积。

  4. 更新解和残差（第一部分）：

    • \(\vec{q}_k = \vec{u}_k - \alpha_k \vec{v}_k\)。
    • 这是一个中间步骤，用于更精确地更新解。

  5. 进行第二次矩阵向量乘法并计算权重 \(\omega_k\)：

    • 计算 \(\vec{w}_k = A \vec{q}_k\)。
    • 计算权重：\(\omega_k = (\vec{w}_k, \vec{r}_k - \alpha_k \vec{w}_k) / (\vec{w}_k, \vec{w}_k)\)。这个 \(\omega_k\) 可以看作是一个局部最速下降步长，用于加速收敛。

  6. 最终更新解和残差：

    • 更新解：\(\vec{x}_{k+1} = \vec{x}_k + \alpha_k \vec{p}_k + \omega_k \vec{q}_k\)。
    • 更新残差：\(\vec{r}_{k+1} = \vec{r}_k - \alpha_k A \vec{p}_k - \omega_k A \vec{q}_k = \vec{q}_k - \omega_k \vec{w}_k\)。
    • 检查 \(\| \vec{r}_{k+1} \|\) 是否小于预设的容忍度（tolerance）。如果满足，则 \(\vec{x}_{k+1}\) 即为所求近似解；否则，令 \(k = k+1\)，返回第二步继续迭代。

第三步：算法总结与特点

  • 优势：CGS完全避免了与 \(A^T\) 的运算，这在 \(A^T\) 难以显式表示或应用时是一个巨大优点。相比于BiCG，它有时能达到“平方”的收敛效果，即收敛所需的迭代次数更少。
  • 劣势：由于是平方化过程，CGS的收敛行为可能比BiCG更不规则，残差范数的震荡可能更剧烈。它也可能比BiCG更容易发生“breakdown”或接近breakdown的情况（即 \(\rho_k\) 或分母内积接近0）。
  • 关键操作：每次迭代需要进行两次矩阵向量乘法（计算 \(\vec{v}_k\) 和 \(\vec{w}_k\)），这是主要的计算开销。

通过以上步骤，CGS算法在Krylov子空间中系统地构建近似解，最终收敛到原非对称线性方程组的解。

基于Krylov子空间的CGS算法解非对称线性方程组 我将为你讲解共轭梯度平方法（Conjugate Gradient Squared, CGS），这是一种用于求解非对称线性方程组的Krylov子空间方法。我们将从一个具体问题开始。 问题描述 考虑一个非对称的线性方程组 \( A\vec{x} = \vec{b} \)，其中： \( A \) 是一个 \( n \times n \) 的非对称矩阵（即 \( A \neq A^T \)）。 \( \vec{b} \) 是一个已知的 \( n \times 1 \) 向量。 \( \vec{x} \) 是我们需要求解的 \( n \times 1 \) 未知向量。 我们的目标是找到一个向量 \( \vec{x} \)，使得方程成立。当矩阵 \( A \) 规模很大且稀疏时，直接解法（如LU分解）计算成本高昂，迭代法如CGS则更为高效。 解题过程 CGS算法是双共轭梯度法（BiCG）的一种变体，它通过平方化的思想来避免BiCG中与 \( A^T \) 相关的运算，从而在某些情况下能更快地收敛。其核心思想是在Krylov子空间 \( \mathcal{K}_ m(A; \vec{r}_ 0) = \text{span}\{\vec{r}_ 0, A\vec{r}_ 0, A^2\vec{r}_ 0, \dots, A^{m-1}\vec{r}_ 0\} \) 中寻找近似解，其中 \( \vec{r}_ 0 \) 是初始残差。 第一步：算法初始化 算法从一个初始猜测解 \( \vec{x}_ 0 \) 开始。通常，如果没有任何先验知识，可以设 \( \vec{x}_ 0 = \vec{0} \)。 计算初始残差：\( \vec{r}_ 0 = \vec{b} - A\vec{x}_ 0 \)。 选择一个辅助向量 \( \vec{\tilde{r}}_ 0 \)，它必须与 \( \vec{r}_ 0 \) 非正交（即内积不为零）。最简单且常见的选择是令 \( \vec{\tilde{r}}_ 0 = \vec{r}_ 0 \)。 初始化一些中间向量：令 \( \vec{p}_ 0 = \vec{r}_ 0 \)，\( \vec{q}_ 0 = \vec{0} \)。 设置初始参数：\( \rho_ 0 = \alpha_ 0 = \omega_ 0 = 1 \)。 开始迭代，对于 \( k = 0, 1, 2, \dots \)，直到残差 \( \vec{r}_ k \) 的范数足够小为止。 第二步：迭代步骤（第 k 步） CGS的迭代过程涉及多个向量的更新，我们一步步来看。 计算迭代参数 \( \rho_ k \) 和 \( \beta_ k \) ： \( \rho_ k = (\vec{\tilde{r}}_ 0, \vec{r}_ k) \)（计算辅助残差与当前残差的内积）。 如果 \( \rho_ k \) 接近于0，算法可能会失败（称为“breakdown”）。 \( \beta_ k = (\rho_ k / \rho_ {k-1}) \times (\alpha_ {k-1} / \omega_ {k-1}) \)。这里用到了上一步的参数。 更新搜索方向向量 \( \vec{u}_ k \) 和 \( \vec{p}_ k \) ： \( \vec{u} k = \vec{r} k + \beta_ k \times \omega {k-1} \times (\vec{u} {k-1} - \omega_ {k-1} A \vec{p}_ {k-1}) \)（对于k=0，相关历史向量视为0）。 \( \vec{p} k = \vec{u} k + \beta_ k \times (\vec{p} {k-1} - \omega {k-1} A \vec{p}_ {k-1}) \)。 这些更新公式看似复杂，但其目的是构造一个在某种意义下“共轭”的搜索方向，其推导源于对BiCG算法的多项式平方化。 进行矩阵向量乘法并计算步长 \( \alpha_ k \) ： 计算 \( \vec{v}_ k = A \vec{p}_ k \)。 计算步长：\( \alpha_ k = \rho_ k / (\vec{\tilde{r}}_ 0, \vec{v}_ k) \)。分母是辅助残差与 \( \vec{v}_ k \) 的内积。 更新解和残差（第一部分） ： \( \vec{q}_ k = \vec{u}_ k - \alpha_ k \vec{v}_ k \)。 这是一个中间步骤，用于更精确地更新解。 进行第二次矩阵向量乘法并计算权重 \( \omega_ k \) ： 计算 \( \vec{w}_ k = A \vec{q}_ k \)。 计算权重：\( \omega_ k = (\vec{w}_ k, \vec{r}_ k - \alpha_ k \vec{w}_ k) / (\vec{w}_ k, \vec{w}_ k) \)。这个 \( \omega_ k \) 可以看作是一个局部最速下降步长，用于加速收敛。 最终更新解和残差 ： 更新解：\( \vec{x}_ {k+1} = \vec{x}_ k + \alpha_ k \vec{p}_ k + \omega_ k \vec{q}_ k \)。 更新残差：\( \vec{r}_ {k+1} = \vec{r}_ k - \alpha_ k A \vec{p}_ k - \omega_ k A \vec{q}_ k = \vec{q}_ k - \omega_ k \vec{w}_ k \)。 检查 \( \| \vec{r} {k+1} \| \) 是否小于预设的容忍度（tolerance）。如果满足，则 \( \vec{x} {k+1} \) 即为所求近似解；否则，令 \( k = k+1 \)，返回第二步继续迭代。 第三步：算法总结与特点 优势 ：CGS完全避免了与 \( A^T \) 的运算，这在 \( A^T \) 难以显式表示或应用时是一个巨大优点。相比于BiCG，它有时能达到“平方”的收敛效果，即收敛所需的迭代次数更少。 劣势 ：由于是平方化过程，CGS的收敛行为可能比BiCG更不规则，残差范数的震荡可能更剧烈。它也可能比BiCG更容易发生“breakdown”或接近breakdown的情况（即 \( \rho_ k \) 或分母内积接近0）。 关键操作 ：每次迭代需要进行两次矩阵向量乘法（计算 \( \vec{v}_ k \) 和 \( \vec{w}_ k \)），这是主要的计算开销。 通过以上步骤，CGS算法在Krylov子空间中系统地构建近似解，最终收敛到原非对称线性方程组的解。