好的，作为“线性代数算法”领域的无所不知大神，我将为你讲解一个尚未在列表中出现的重要算法。
题目：基于Krylov子空间的MINRES算法解对称不定线性方程组
1. 问题描述与背景
我们要求解一个大规模、稀疏的线性方程组：
\[A\mathbf{x} = \mathbf{b} \]
其中 \(A \in \mathbb{R}^{n \times n}\) 是一个对称矩阵，但它不一定正定（因此可能含有零特征值或负特征值），矩阵规模 \(n\) 非常大（例如 \(n > 10^5\) ），并且矩阵 \(A\) 是稀疏的（即大部分元素为零）。向量 \(\mathbf{b} \in \mathbb{R}^{n}\) 是已知的右端项。
为什么这是一个挑战？
- 对称不定：对于对称正定矩阵，我们有非常高效且稳定的共轭梯度法（CG）。但CG算法要求矩阵正定，当矩阵不定时，CG算法在迭代过程中可能会遇到“分母为零”的数值问题而崩溃。
- 大规模稀疏：直接求解法（如LU分解）因内存和时间消耗过大而不适用，必须使用迭代法。我们需要一个能利用矩阵稀疏性、内存友好且能处理对称不定性的迭代算法。
MINRES（最小残差法）正是为解决此类问题而设计的Krylov子空间方法。
2. 核心思想：在Krylov子空间中寻找最优解
MINRES的核心思想与GMRES类似，但专门针对对称矩阵进行了优化。其基本步骤如下：
- 构建Krylov子空间：
从初始猜测解 \(\mathbf{x}_0\)（通常取零向量）出发，计算初始残差 \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\)。
算法在第 \(m\) 步迭代时，构造一个维度为 \(m\) 的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{x}_m \in \mathbf{x}_0 + \mathcal{K}_m(A, \mathbf{r}_0) $，使得残差 $ \mathbf{r}_m = \mathbf{b} - A\mathbf{x}_m $ 的**2-范数最小**，即求解：
\[ \min_{\mathbf{x} \in \mathbf{x}_0 + \mathcal{K}_m(A, \mathbf{r}_0)} \| \mathbf{b} - A\mathbf{x} \|_2 \]
- 利用对称性简化计算：
对于对称矩阵 \(A\)，我们可以使用Lanczos过程（而不是用于非对称矩阵的Arnoldi过程）来为Krylov子空间生成一组标准正交基 \(\{ \mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_m \}\)，其中 \(\mathbf{v}_1 = \mathbf{r}_0 / \| \mathbf{r}_0 \|_2\)。
Lanczos过程会产生一个三对角矩阵 \(T_m\)：
\[ A V_m = V_m T_m + \beta_{m+1} \mathbf{v}_{m+1} \mathbf{e}_m^T \]
其中 $ V_m = [\mathbf{v}_1, \dots, \mathbf{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} \]
$ \alpha_i, \beta_i $ 是Lanczos迭代中产生的标量。
3. 算法推导与步骤详解
现在，我们将最小化问题转化为一个更易求解的小规模问题。
第一步：问题转化
令近似解 \(\mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m\)，其中 \(\mathbf{y}_m \in \mathbb{R}^m\) 是系数向量。
则残差为：
\[\mathbf{r}_m = \mathbf{b} - A(\mathbf{x}_0 + V_m \mathbf{y}_m) = \mathbf{r}_0 - A V_m \mathbf{y}_m \]
代入Lanczos关系式 \(A V_m = V_m T_m + \beta_{m+1} \mathbf{v}_{m+1} \mathbf{e}_m^T\) 和 \(\mathbf{r}_0 = \| \mathbf{r}_0 \|_2 \mathbf{v}_1\)：
\[\mathbf{r}_m = \| \mathbf{r}_0 \|_2 \mathbf{v}_1 - (V_m T_m + \beta_{m+1} \mathbf{v}_{m+1} \mathbf{e}_m^T) \mathbf{y}_m = V_m (\| \mathbf{r}_0 \|_2 \mathbf{e}_1 - T_m \mathbf{y}_m) - \beta_{m+1} (\mathbf{e}_m^T \mathbf{y}_m) \mathbf{v}_{m+1} \]
由于 \(V_m\) 的列相互正交且与 \(\mathbf{v}_{m+1}\) 正交，残差的范数平方为：
\[\| \mathbf{r}_m \|_2^2 = \| V_m (\| \mathbf{r}_0 \|_2 \mathbf{e}_1 - T_m \mathbf{y}_m) \|_2^2 + \| \beta_{m+1} (\mathbf{e}_m^T \mathbf{y}_m) \mathbf{v}_{m+1} \|_2^2 \]
因为 \(V_m\) 是正交矩阵，其作用不改变2-范数（即 \(\|V_m \mathbf{z}\|_2 = \|\mathbf{z}\|_2\)），所以：
\[\| \mathbf{r}_m \|_2^2 = \| \| \mathbf{r}_0 \|_2 \mathbf{e}_1 - T_m \mathbf{y}_m \|_2^2 + |\beta_{m+1} (\mathbf{e}_m^T \mathbf{y}_m)|^2 \]
我们的最小化问题等价于：
\[\min_{\mathbf{y}_m \in \mathbb{R}^m} \| \| \mathbf{r}_0 \|_2 \mathbf{e}_1 - T_m \mathbf{y}_m \|_2 \]
这是一个规模为 \(m\) 的最小二乘问题，系数矩阵是三对角矩阵 \(T_m\)。
注意，第二项 \(|\beta_{m+1} (\mathbf{e}_m^T \mathbf{y}_m)|^2\) 在优化 \(\mathbf{y}_m\) 时无法控制，但MINRES通过上述转化，巧妙地将大规模问题降维为一个仅与 \(T_m\) 相关的小规模问题。
第二步：利用QR分解高效求解
为了稳定高效地求解这个小规模最小二乘问题，MINRES使用了对三对角矩阵 \(T_m\) 进行QR分解的策略，并且采用递推更新，避免每次迭代都重新进行分解。
- 对 \(T_m\) 进行QR分解： \(T_m = Q_m R_m\)，其中 \(Q_m\) 是正交矩阵（一系列Givens旋转的乘积）， \(R_m\) 是上三角矩阵（由于 \(T_m\) 是三对角的， \(R_m\) 是带状上三角矩阵，只有两条非零对角线上方有元素）。
- 最小二乘问题变为：
\[ \min \| \| \mathbf{r}_0 \|_2 \mathbf{e}_1 - Q_m R_m \mathbf{y}_m \|_2 = \min \| Q_m^T (\| \mathbf{r}_0 \|_2 \mathbf{e}_1) - R_m \mathbf{y}_m \|_2 \]
令 $ \mathbf{g}_m = Q_m^T (\| \mathbf{r}_0 \|_2 \mathbf{e}_1) $。
- 由于 \(R_m\) 是上三角矩阵，最优解 \(\mathbf{y}_m\) 可以通过回代法轻松求得。更重要的是，当从第 \(m\) 步迭代到第 \(m+1\) 步时， \(T_{m+1}\) 只是在 \(T_m\) 的右下角添加了一行一列。MINRES算法利用此特性，通过更新已有的QR分解（即追加一个Givens旋转来消去新产生的非零次对角元）来高效计算 \(\mathbf{y}_{m+1}\)，而无需重新分解整个 \(T_{m+1}\)。
第三步：递推更新近似解
得到 \(\mathbf{y}_m\) 后，近似解为 \(\mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m\)。同样，我们也可以避免存储所有基向量 \(V_m\)。通过分析更新公式，可以发现 \(\mathbf{x}_m\) 可以通过一个简单的三项递推关系从前两个迭代解 \(\mathbf{x}_{m-1}\) 和 \(\mathbf{x}_{m-2}\) 计算出来，只需存储少量向量。这使得MINRES的内存开销是常数级的 \(O(n)\)，非常适合大规模问题。
4. 算法流程总结
下面是MINRES算法的简化伪代码流程：
-
初始化：
- 输入：矩阵 \(A\)，向量 \(\mathbf{b}\)，初始猜测 \(\mathbf{x}_0\)，容差 \(\epsilon\)，最大迭代次数 \(m_{\text{max}}\)。
- 计算 \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\)， \(\beta_1 = \| \mathbf{r}_0 \|_2\)。
- 若 \(\beta_1 < \epsilon\)，则 \(\mathbf{x}_0\) 即为解，算法终止。
- 令 \(\mathbf{v}_1 = \mathbf{r}_0 / \beta_1\)。
-
Lanczos过程与QR更新：
- 对于 \(j = 1, 2, \dots, m_{\text{max}}\)：
a. Lanczos步骤：计算 \(\mathbf{w}_j = A \mathbf{v}_j\)，进行正交化得到 \(\alpha_j = \mathbf{v}_j^T \mathbf{w}_j\)， \(\mathbf{w}_j := \mathbf{w}_j - \alpha_j \mathbf{v}_j - \beta_j \mathbf{v}_{j-1}\)（令 \(\beta_1 \mathbf{v}_0 = \mathbf{0}\)），然后 \(\beta_{j+1} = \| \mathbf{w}_j \|_2\)， \(\mathbf{v}_{j+1} = \mathbf{w}_j / \beta_{j+1}\)。
b. QR更新：利用新计算的 \(\alpha_j, \beta_{j+1}\) 以及前一步的旋转矩阵，通过追加1-2个Givens旋转来更新对 \(T_j\) 的QR分解，并同时更新变换后的右端项向量 \(\mathbf{g}_j\)。
c. 检查收敛：计算当前残差范数估计（可从 \(\mathbf{g}_j\) 的最后一个分量直接得到）。若小于 \(\epsilon\)，则退出循环。
d. 更新解：利用递推公式由前两步的解向量 \(\mathbf{x}_{j-1}\)、 \(\mathbf{x}_{j-2}\) 和QR分解得到的系数计算 \(\mathbf{x}_j\)。
- 对于 \(j = 1, 2, \dots, m_{\text{max}}\)：
-
输出：迭代终止时的近似解 \(\mathbf{x}_m\)。
5. 关键特性与注意事项
- 适用性：严格适用于对称矩阵（包括不定矩阵）。对于非对称矩阵，应使用GMRES。
- 最优性：在每一步迭代中，它在当前Krylov子空间内最小化残差2-范数。
- 短递归：由于对称性和递推求解策略，算法只需要存储最近几个向量，内存占用固定且小。
- 收敛性：对于对称矩阵，其收敛速度与矩阵特征值的分布有关。即使矩阵不定（有负特征值），只要特征值不聚集在原点附近，通常也能很好收敛。
- 预处理：为了加速收敛，通常会结合预处理技术，例如使用一个对称正定的预处理矩阵 \(M\)，将原系统转化为等价的预处理系统 \(M^{-1}A\mathbf{x} = M^{-1}\mathbf{b}\)（注意保持对称性，通常要求 \(M\) 是 \(A\) 的一个好的近似）。预处理后的MINRES算法需要稍作调整以适应预处理Lanczos过程。
总结：MINRES算法巧妙地结合了对称矩阵的Lanczos过程、QR分解的递推更新和三项递推的解更新公式，为求解大规模稀疏对称不定线性方程组提供了一个内存高效、数值稳定且理论上最优（在Krylov子空间内）的迭代解法。