基于Krylov子空间方法的MINRES算法解对称不定线性方程组
字数 4658 2025-12-23 23:29:20

基于Krylov子空间方法的MINRES算法解对称不定线性方程组

题目描述

我们考虑求解大型稀疏对称线性方程组：

\[Ax = b \]

其中 \(A \in \mathbb{R}^{n \times n}\) 是一个对称矩阵，但不一定正定（即 \(A\) 可能是不定的，其特征值既有正数也有负数），\(b \in \mathbb{R}^{n}\) 是已知右端向量，\(x \in \mathbb{R}^{n}\) 是待求解的未知向量。矩阵 \(A\) 是大型且稀疏的，因此直接法（如LU分解）因内存和计算成本过高而不适用，我们需要迭代法。

核心问题：当 \(A\) 对称但不正定时，经典的共轭梯度法（CG）无法使用，因为CG要求矩阵对称正定（SPD）。我们需要一种适用于对称不定系统的迭代算法。MINRES（最小残差法）正是为此设计的Krylov子空间方法。本题目将详细讲解MINRES算法的原理、构造和求解步骤。

解题过程详解

1. 问题背景与核心思想

  • 为什么CG不能用？
    CG算法在每步迭代中需要保证搜索方向 \(A\)-共轭（即 \(p_i^T A p_j = 0, i \neq j\)）。当 \(A\) 不定时，\(A\)-内积可能为零或负，导致算法中出现零除或数值不稳定，甚至不收敛。

  • MINRES的核心思想：
    在Krylov子空间

\[ \mathcal{K}_m(A, r_0) = \text{span}\{ r_0, A r_0, A^2 r_0, \dots, A^{m-1} r_0 \} \]

中寻找近似解 \(x_m\)，使得残差范数 \(\|r_m\|_2 = \|b - A x_m\|_2\) 最小。
由于只使用2-范数，不涉及 \(A\)-内积，因此不要求 \(A\) 正定，只需对称。

  • 优势：
    对对称不定系统，MINRES在理论上保证残差范数单调不增，且最多 \(n\) 步迭代得到精确解（不考虑舍入误差）。

2. 算法构造的关键步骤

MINRES的实现通常基于Lanczos过程，将 \(A\) 约化为三对角形式。

步骤1：Lanczos三对角化
给定初始向量 \(r_0 = b - A x_0\)（通常 \(x_0 = 0\)），令 \(v_1 = r_0 / \beta_1\)，其中 \(\beta_1 = \|r_0\|_2\)。
Lanczos过程生成正交向量 \(v_1, v_2, \dots, v_m\) 满足：

\[A V_m = V_m T_m + \beta_{m+1} v_{m+1} e_m^T \]

其中 \(V_m = [v_1, \dots, v_m]\) 的列正交，\(T_m\) 是实对称三对角矩阵：

\[T_m = \begin{bmatrix} \alpha_1 & \beta_2 & & \\ \beta_2 & \alpha_2 & \beta_3 & \\ & \beta_3 & \alpha_3 & \ddots \\ & & \ddots & \ddots & \beta_m \\ & & & \beta_m & \alpha_m \end{bmatrix} \]

递推公式：

\[\beta_{j+1} v_{j+1} = A v_j - \alpha_j v_j - \beta_j v_{j-1}, \quad \alpha_j = v_j^T A v_j \]

这里约定 \(v_0 = 0, \beta_1 v_1 = r_0\)。

步骤2：将原问题投影到Krylov子空间
设近似解 \(x_m = x_0 + V_m y_m\)，其中 \(y_m \in \mathbb{R}^m\)。
残差：

\[r_m = b - A x_m = r_0 - A V_m y_m = \beta_1 v_1 - (V_m T_m + \beta_{m+1} v_{m+1} e_m^T) y_m \]

利用 \(V_m^T V_m = I\) 和 \(V_m^T v_{m+1} = 0\)，得：

\[r_m = V_m (\beta_1 e_1 - T_m y_m) - \beta_{m+1} (e_m^T y_m) v_{m+1} \]

由于 \(V_m\) 的列与 \(v_{m+1}\) 正交，残差2-范数：

\[\|r_m\|_2^2 = \|\beta_1 e_1 - T_m y_m\|_2^2 + |\beta_{m+1} (e_m^T y_m)|^2 \]

MINRES通过选择 \(y_m\) 来最小化 \(\|\beta_1 e_1 - T_m y_m\|_2\)，从而最小化 \(\|r_m\|_2\)。

步骤3：通过QR分解最小化残差
对 \(T_m\) 进行QR分解：\(Q_m T_m = R_m\)，其中 \(R_m\) 是上三角矩阵，\(Q_m\) 是正交矩阵。由于 \(T_m\) 是三对角的，其QR分解可通过Givens旋转高效递推计算。
最小化问题转化为：

\[\min_{y_m} \|\beta_1 e_1 - T_m y_m\|_2 = \min_{y_m} \|Q_m (\beta_1 e_1 - T_m y_m)\|_2 = \min_{y_m} \|g_m - R_m y_m\|_2 \]

其中 \(g_m = Q_m (\beta_1 e_1) = ( \gamma_1, \gamma_2, \dots, \gamma_m, \gamma_{m+1})^T\)。
由于 \(R_m\) 是上三角矩阵，最小二乘解可通过回代得到。更重要的是，在迭代中，我们可以通过更新 \(x_m = x_0 + V_m y_m\) 而不显式存储 \(V_m\)。

步骤4：递推更新解向量
实际实现中，我们通过递推Givens旋转来更新解。设第 \(j\) 步的Givens旋转作用于 \(T_m\) 的 \(2\times2\) 子块，将次对角线消为零。这些旋转可同时作用于右端向量 \(g_m\)，从而得到 \(y_m\) 的递推关系。
解更新公式为：

\[x_j = x_{j-1} + c_j \tau_j w_j \]

其中 \(w_j\) 是递推生成的向量，\(c_j, \tau_j\) 是旋转产生的系数。具体细节涉及较多记号，但核心是：每步迭代只需常数次内积和向量更新，存储开销为 \(O(n)\)。

3. MINRES算法伪代码（简化版）

输入: 对称矩阵 \(A\)，右端项 \(b\)，初始猜测 \(x_0\)（通常为零），容差 \(\epsilon\)，最大迭代次数 \(m_{\max}\)
输出: 近似解 \(x\)

  1. 初始化：

\[ r_0 = b - A x_0, \; \beta_1 = \|r_0\|_2, \; v_1 = r_0 / \beta_1, \; w_0 = w_1 = 0, \; x_0 \text{给定} \]

\[ \eta = \beta_1, \; c_0 = c_1 = 1, \; s_0 = s_1 = 0 \]

  1. For \(j = 1, 2, \dots, m_{\max}\) do：
    • Lanczos步骤：

\[ \alpha_j = v_j^T A v_j \]

\[ \hat{v}_{j+1} = A v_j - \alpha_j v_j - \beta_j v_{j-1} \quad (\text{约定 } v_0=0, \beta_1 v_1 = r_0) \]

\[ \beta_{j+1} = \|\hat{v}_{j+1}\|_2 \]

\[ v_{j+1} = \hat{v}_{j+1} / \beta_{j+1} \]

  • 应用前一步的旋转到新列（更新三对角矩阵的分解）：

\[ \rho_{j-1} = s_{j-1} \alpha_j + c_{j-1} \beta_{j+1} \]

\[ \delta_j = c_{j-1} \alpha_j - s_{j-1} \beta_{j+1} \]

\[ \gamma_j = \sqrt{ \delta_j^2 + \beta_{j+1}^2 } \]

  • 计算当前Givens旋转：

\[ c_j = \delta_j / \gamma_j, \quad s_j = \beta_{j+1} / \gamma_j \]

  • 更新右端向量变换：

\[ \eta_j = -s_j \eta_{j-1}, \quad \eta_{j-1} = c_j \eta_{j-1} \]

  • 更新解向量：

\[ w_j = (v_j - \rho_{j-2} w_{j-2} - \rho_{j-1} w_{j-1}) / \gamma_j \quad (\text{约定 } \rho_0=\rho_{-1}=0) \]

\[ x_j = x_{j-1} + (\eta_{j-1} / \gamma_j) w_j \]

  • 检查收敛：如果 \(|\eta_j| < \epsilon \cdot \beta_1\)，则终止并输出 \(x_j\)。
  1. 输出 \(x_j\) 作为近似解。

4. 算法特性与注意事项

  • 对称性要求：MINRES严格需要矩阵 \(A\) 对称。如果不对称，残差最小化方法应为GMRES。
  • 存储开销：每步只需存储几个 \(n\) 维向量（\(v_j, v_{j-1}, w_j, w_{j-1}, x_j\)），内存友好。
  • 收敛性：残差范数 \(\|r_j\|_2\) 单调不增。收敛速度依赖于 \(A\) 的特征值分布。若特征值聚集在正负两个区间，且远离零，则收敛较快；若有特征值接近零，则可能收敛缓慢。
  • 与CG的关系：当 \(A\) 对称正定时，MINRES的残差范数曲线与CG的残差范数曲线相同，但MINRES的计算量稍大。
  • 预处理：为加速收敛，常用预处理。对称预处理矩阵 \(M\) 应使 \(M^{-1}A\) 的特征值更聚集。注意预处理后系统应为对称形式（如 \(M^{-1}A\) 对称）。

5. 总结

MINRES是求解大型稀疏对称不定线性方程组的经典Krylov子空间方法。它通过Lanczos过程构建正交基，并利用Givens旋转递推地求解子空间上的最小二乘问题，从而保证残差2-范数最小。算法具有短递推格式，内存开销小，适合大规模稀疏问题。理解MINRES的关键在于掌握Lanczos过程、Givens旋转的递推应用，以及如何从投影子问题中更新近似解。

基于Krylov子空间方法的MINRES算法解对称不定线性方程组 题目描述 我们考虑求解大型稀疏对称线性方程组： \[ Ax = b \] 其中 \(A \in \mathbb{R}^{n \times n}\) 是一个 对称 矩阵，但 不一定正定 （即 \(A\) 可能是不定的，其特征值既有正数也有负数），\(b \in \mathbb{R}^{n}\) 是已知右端向量，\(x \in \mathbb{R}^{n}\) 是待求解的未知向量。矩阵 \(A\) 是大型且稀疏的，因此直接法（如LU分解）因内存和计算成本过高而不适用，我们需要迭代法。 核心问题 ：当 \(A\) 对称但不正定时，经典的共轭梯度法（CG）无法使用，因为CG要求矩阵对称正定（SPD）。我们需要一种适用于对称不定系统的迭代算法。MINRES（最小残差法）正是为此设计的Krylov子空间方法。本题目将详细讲解MINRES算法的原理、构造和求解步骤。 解题过程详解 1. 问题背景与核心思想 为什么CG不能用？ CG算法在每步迭代中需要保证搜索方向 \(A\)-共轭（即 \(p_ i^T A p_ j = 0, i \neq j\)）。当 \(A\) 不定时，\(A\)-内积可能为零或负，导致算法中出现零除或数值不稳定，甚至不收敛。 MINRES的核心思想 ： 在Krylov子空间 \[ \mathcal{K}_ m(A, r_ 0) = \text{span}\{ r_ 0, A r_ 0, A^2 r_ 0, \dots, A^{m-1} r_ 0 \} \] 中寻找近似解 \(x_ m\)，使得 残差范数 \(\|r_ m\|_ 2 = \|b - A x_ m\|_ 2\) 最小 。 由于只使用2-范数，不涉及 \(A\)-内积，因此不要求 \(A\) 正定，只需对称。 优势 ： 对对称不定系统，MINRES在理论上保证残差范数单调不增，且最多 \(n\) 步迭代得到精确解（不考虑舍入误差）。 2. 算法构造的关键步骤 MINRES的实现通常基于 Lanczos过程 ，将 \(A\) 约化为三对角形式。 步骤1：Lanczos三对角化 给定初始向量 \(r_ 0 = b - A x_ 0\)（通常 \(x_ 0 = 0\)），令 \(v_ 1 = r_ 0 / \beta_ 1\)，其中 \(\beta_ 1 = \|r_ 0\| 2\)。 Lanczos过程生成正交向量 \(v_ 1, v_ 2, \dots, v_ m\) 满足： \[ A V_ m = V_ m T_ m + \beta {m+1} v_ {m+1} e_ m^T \] 其中 \(V_ m = [ v_ 1, \dots, v_ m]\) 的列正交，\(T_ m\) 是实对称三对角矩阵： \[ T_ m = \begin{bmatrix} \alpha_ 1 & \beta_ 2 & & \\ \beta_ 2 & \alpha_ 2 & \beta_ 3 & \\ & \beta_ 3 & \alpha_ 3 & \ddots \\ & & \ddots & \ddots & \beta_ m \\ & & & \beta_ m & \alpha_ m \end{bmatrix} \] 递推公式： \[ \beta_ {j+1} v_ {j+1} = A v_ j - \alpha_ j v_ j - \beta_ j v_ {j-1}, \quad \alpha_ j = v_ j^T A v_ j \] 这里约定 \(v_ 0 = 0, \beta_ 1 v_ 1 = r_ 0\)。 步骤2：将原问题投影到Krylov子空间 设近似解 \(x_ m = x_ 0 + V_ m y_ m\)，其中 \(y_ m \in \mathbb{R}^m\)。 残差： \[ r_ m = b - A x_ m = r_ 0 - A V_ m y_ m = \beta_ 1 v_ 1 - (V_ m T_ m + \beta_ {m+1} v_ {m+1} e_ m^T) y_ m \] 利用 \(V_ m^T V_ m = I\) 和 \(V_ m^T v_ {m+1} = 0\)，得： \[ r_ m = V_ m (\beta_ 1 e_ 1 - T_ m y_ m) - \beta_ {m+1} (e_ m^T y_ m) v_ {m+1} \] 由于 \(V_ m\) 的列与 \(v_ {m+1}\) 正交，残差2-范数： \[ \|r_ m\|_ 2^2 = \|\beta_ 1 e_ 1 - T_ m y_ m\| 2^2 + |\beta {m+1} (e_ m^T y_ m)|^2 \] MINRES通过选择 \(y_ m\) 来最小化 \(\|\beta_ 1 e_ 1 - T_ m y_ m\|_ 2\)，从而最小化 \(\|r_ m\|_ 2\)。 步骤3：通过QR分解最小化残差 对 \(T_ m\) 进行QR分解：\(Q_ m T_ m = R_ m\)，其中 \(R_ m\) 是上三角矩阵，\(Q_ m\) 是正交矩阵。由于 \(T_ m\) 是三对角的，其QR分解可通过Givens旋转高效递推计算。 最小化问题转化为： \[ \min_ {y_ m} \|\beta_ 1 e_ 1 - T_ m y_ m\| 2 = \min {y_ m} \|Q_ m (\beta_ 1 e_ 1 - T_ m y_ m)\| 2 = \min {y_ m} \|g_ m - R_ m y_ m\| 2 \] 其中 \(g_ m = Q_ m (\beta_ 1 e_ 1) = ( \gamma_ 1, \gamma_ 2, \dots, \gamma_ m, \gamma {m+1})^T\)。 由于 \(R_ m\) 是上三角矩阵，最小二乘解可通过回代得到。更重要的是，在迭代中，我们可以通过更新 \(x_ m = x_ 0 + V_ m y_ m\) 而不显式存储 \(V_ m\)。 步骤4：递推更新解向量 实际实现中，我们通过递推Givens旋转来更新解。设第 \(j\) 步的Givens旋转作用于 \(T_ m\) 的 \(2\times2\) 子块，将次对角线消为零。这些旋转可同时作用于右端向量 \(g_ m\)，从而得到 \(y_ m\) 的递推关系。 解更新公式为： \[ x_ j = x_ {j-1} + c_ j \tau_ j w_ j \] 其中 \(w_ j\) 是递推生成的向量，\(c_ j, \tau_ j\) 是旋转产生的系数。具体细节涉及较多记号，但核心是：每步迭代只需常数次内积和向量更新，存储开销为 \(O(n)\)。 3. MINRES算法伪代码（简化版） 输入 : 对称矩阵 \(A\)，右端项 \(b\)，初始猜测 \(x_ 0\)（通常为零），容差 \(\epsilon\)，最大迭代次数 \(m_ {\max}\) 输出 : 近似解 \(x\) 初始化： \[ r_ 0 = b - A x_ 0, \; \beta_ 1 = \|r_ 0\|_ 2, \; v_ 1 = r_ 0 / \beta_ 1, \; w_ 0 = w_ 1 = 0, \; x_ 0 \text{给定} \] \[ \eta = \beta_ 1, \; c_ 0 = c_ 1 = 1, \; s_ 0 = s_ 1 = 0 \] For \(j = 1, 2, \dots, m_ {\max}\) do： Lanczos步骤 ： \[ \alpha_ j = v_ j^T A v_ j \] \[ \hat{v} {j+1} = A v_ j - \alpha_ j v_ j - \beta_ j v {j-1} \quad (\text{约定 } v_ 0=0, \beta_ 1 v_ 1 = r_ 0) \] \[ \beta_ {j+1} = \|\hat{v} {j+1}\| 2 \] \[ v {j+1} = \hat{v} {j+1} / \beta_ {j+1} \] 应用前一步的旋转到新列 （更新三对角矩阵的分解）： \[ \rho_ {j-1} = s_ {j-1} \alpha_ j + c_ {j-1} \beta_ {j+1} \] \[ \delta_ j = c_ {j-1} \alpha_ j - s_ {j-1} \beta_ {j+1} \] \[ \gamma_ j = \sqrt{ \delta_ j^2 + \beta_ {j+1}^2 } \] 计算当前Givens旋转 ： \[ c_ j = \delta_ j / \gamma_ j, \quad s_ j = \beta_ {j+1} / \gamma_ j \] 更新右端向量变换 ： \[ \eta_ j = -s_ j \eta_ {j-1}, \quad \eta_ {j-1} = c_ j \eta_ {j-1} \] 更新解向量 ： \[ w_ j = (v_ j - \rho_ {j-2} w_ {j-2} - \rho_ {j-1} w_ {j-1}) / \gamma_ j \quad (\text{约定 } \rho_ 0=\rho_ {-1}=0) \] \[ x_ j = x_ {j-1} + (\eta_ {j-1} / \gamma_ j) w_ j \] 检查收敛 ：如果 \(|\eta_ j| < \epsilon \cdot \beta_ 1\)，则终止并输出 \(x_ j\)。 输出 \(x_ j\) 作为近似解。 4. 算法特性与注意事项 对称性要求 ：MINRES严格需要矩阵 \(A\) 对称。如果不对称，残差最小化方法应为GMRES。 存储开销 ：每步只需存储几个 \(n\) 维向量（\(v_ j, v_ {j-1}, w_ j, w_ {j-1}, x_ j\)），内存友好。 收敛性 ：残差范数 \(\|r_ j\|_ 2\) 单调不增。收敛速度依赖于 \(A\) 的特征值分布。若特征值聚集在正负两个区间，且远离零，则收敛较快；若有特征值接近零，则可能收敛缓慢。 与CG的关系 ：当 \(A\) 对称正定时，MINRES的残差范数曲线与CG的残差范数曲线相同，但MINRES的计算量稍大。 预处理 ：为加速收敛，常用预处理。对称预处理矩阵 \(M\) 应使 \(M^{-1}A\) 的特征值更聚集。注意预处理后系统应为对称形式（如 \(M^{-1}A\) 对称）。 5. 总结 MINRES是求解大型稀疏对称不定线性方程组的经典Krylov子空间方法。它通过Lanczos过程构建正交基，并利用Givens旋转递推地求解子空间上的最小二乘问题，从而保证残差2-范数最小。算法具有短递推格式，内存开销小，适合大规模稀疏问题。理解MINRES的关键在于掌握Lanczos过程、Givens旋转的递推应用，以及如何从投影子问题中更新近似解。