基于Krylov子空间的COCR算法解复对称线性方程组
题目描述
我们考虑求解一个大规模稀疏的复对称线性方程组：
\[A\mathbf{x} = \mathbf{b} \]
其中 \(A \in \mathbb{C}^{n \times n}\) 是一个复对称矩阵（即 \(A = A^T\)，注意不是 Hermitian 矩阵，因为不要求共轭转置相等），且通常 \(n\) 很大，\(A\) 稀疏。\(b \in \mathbb{C}^n\) 是已知向量，我们需要求解 \(x \in \mathbb{C}^n\)。
这类问题出现在许多工程与科学计算领域，例如频域波动方程离散、电磁场计算、量子力学等。因为矩阵是复对称的，传统的共轭梯度法（CG）不适用（CG 要求矩阵是 Hermitian 正定），而 GMRES 等通用方法又可能效率不高。COCR（Conjugate Orthogonal Conjugate Residual）算法正是针对复对称矩阵设计的一种 Krylov 子空间方法，它通过短递归实现，内存开销小，且具有残差正交性。
解题过程循序渐进讲解
我们一步步推导 COCR 算法的原理和迭代格式。
第一步：问题与 Krylov 子空间框架
对于复对称矩阵 \(A\)，我们希望构造迭代解序列 \(\{x_k\}\)，使得第 \(k\) 步的解 \(x_k\) 位于仿射空间：
\[x_k \in x_0 + \mathcal{K}_k(A, r_0) \]
其中初始残差 \(r_0 = b - A x_0\)，Krylov 子空间：
\[\mathcal{K}_k(A, r_0) = \mathrm{span}\{r_0, A r_0, \dots, A^{k-1} r_0\}. \]
COCR 要求残差向量 \(r_k = b - A x_k\) 满足 Galerkin 条件（也称为正交条件）：
\[r_k \perp \mathcal{K}_k(A, r_0) \]
但注意，这里的“垂直”是在标准复内积下吗？不是的，因为 \(A\) 不是 Hermitian，我们不能直接用标准内积空间的正交性。实际上，COCR 算法要求残差满足：
\[r_k \perp A \mathcal{K}_k(A, r_0) \]
这等价于：
\[r_k \perp \mathcal{K}_k(A, A r_0) \]
这称为正交残差条件（与 MINRES 对实对称矩阵的条件类似），但这里针对复对称矩阵设计。
第二步：Lanczos 双正交化与短递归
由于 \(A\) 是复对称，我们考虑 Lanczos 双正交化过程。定义两个 Krylov 子空间：
\[\mathcal{K}_k(A, r_0) \quad \text{和} \quad \mathcal{K}_k(A, A r_0) \]
我们需要构造基向量 \(\{p_0, p_1, \dots, p_{k-1}\}\) 和 \(\{s_0, s_1, \dots, s_{k-1}\}\)，使得它们满足双正交性：
\[(p_i, A p_j) = 0, \quad i \ne j \]
实际上，COCR 采用的是一种特殊的短递归三项递推，类似于对称 Lanczos 过程，但由于矩阵对称不 Hermitian，递归中会出现两个序列：
\[A P_k = S_k T_k \quad \text{或类似形式} \]
其中 \(P_k = [p_0, \dots, p_{k-1}]\), \(S_k = [s_0, \dots, s_{k-1}]\) 是双正交基，\(T_k\) 是三对角矩阵。
推导细节（简化）：
设初始搜索方向 \(p_0 = r_0\)，定义 \(s_0 = A p_0\)。迭代中维持：
\[r_k \perp \mathrm{span}\{s_0, \dots, s_{k-1}\} \]
并且搜索方向满足 \(A P_k = S_k\) 的某种三对角化关系。通过短递归得到：
\[p_{k} = r_{k} + \beta_k p_{k-1} \]
\[ s_{k} = A p_{k} = A r_{k} + \beta_k s_{k-1} \]
然后 \(x_{k+1} = x_k + \alpha_k p_k\)，\(r_{k+1} = r_k - \alpha_k A p_k\)。
第三步：系数确定（正交条件）
利用正交条件 \(r_{k+1} \perp s_j\) 对 \(j=0,\dots,k\)，可得系数公式。
由 \(r_{k+1} = r_k - \alpha_k s_k\)，要求 \((r_{k+1}, s_k) = 0\) 得：
\[\alpha_k = \frac{(r_k, s_k)}{(s_k, s_k)} \]
为了保证 \(r_{k+1} \perp s_{k-1}\)，需要适当选择 \(\beta_k\)。由 \(s_k = A r_k + \beta_k s_{k-1}\) 和 \(r_{k+1} = r_k - \alpha_k s_k\)，代入正交条件 \((r_{k+1}, s_{k-1}) = 0\) 可推导出：
\[\beta_k = -\frac{(A r_k, s_{k-1})}{(s_{k-1}, s_{k-1})} \]
但为了短递归，通常会利用 \(A\) 对称性简化表达式。
经推导（Sogabe 和 Zhang 在 2007 年提出的 COCR 算法），得到更简洁的系数：
\[\alpha_k = \frac{(r_k, r_k)}{(A p_k, p_k)} \]
\[ \beta_k = \frac{(r_{k+1}, r_{k+1})}{(r_k, r_k)} \]
注意，这要求内积是标准复内积 \((u,v) = u^T v\)（不取共轭！），因为 \(A\) 是对称非 Hermitian，所以 \((A p_k, p_k) = p_k^T A p_k\) 是实数吗？不一定，但算法仍然有效，因为正交条件是在这个内积下定义的。
实际上，标准 COCR 算法采用的内积是 \((u,v) = u^T v\)（不共轭），因此称为“复内积”但实为双线性型。
在复对称条件下，有 \((A u, v) = (u, A v)\) 对此内积成立。
第四步：COCR 算法迭代格式
给定初始解 \(x_0\)，计算 \(r_0 = b - A x_0\)，设 \(p_0 = r_0\)。
对 \(k=0,1,2,\dots\) 直到收敛，执行：
- 计算 \(A p_k\)。
- 计算 \(\alpha_k = \frac{(r_k, r_k)}{(A p_k, p_k)}\)，其中内积 \((x,y) = x^T y\)。
- 更新解：\(x_{k+1} = x_k + \alpha_k p_k\)。
- 更新残差：\(r_{k+1} = r_k - \alpha_k A p_k\)。
- 计算 \(\beta_k = \frac{(r_{k+1}, r_{k+1})}{(r_k, r_k)}\)。
- 更新搜索方向：\(p_{k+1} = r_{k+1} + \beta_k p_k\)。
注意：残差范数 \((r_k, r_k)\) 在复数下不一定是实数正数，但算法中直接作复数运算，迭代中会保持正交性。
第五步：算法性质
- 短递归：只需存储几个向量，类似 CG 法。
- 有限步终止：理论上在 \(n\) 步内得到精确解（忽略舍入误差）。
- 残差正交性：满足 \(r_k \perp A \mathcal{K}_k(A, r_0)\)，即 \((r_k, A p_j) = 0\) 对 \(j < k\)。
- 适用于复对称矩阵，不要求正定。
- 可能发生“严重中断”（serious breakdown），当 \((A p_k, p_k) = 0\) 时算法失效，但概率低，有修补方法。
第六步：与 COCG 的关系
COCG（Conjugate Orthogonal CG）算法也是解复对称线性方程组的著名方法，但 COCG 要求矩阵是复对称正定（即 \(x^T A x\) 的实部正定），而 COCR 是 COCG 的变体，更类似 MINRES 的思想，不要求正定性，但需要复对称性。
第七步：预处理
为了加速收敛，可采用预处理技术。设预处理子 \(M \approx A\)，且 \(M\) 容易求逆，将原方程转换为：
\[M^{-1} A x = M^{-1} b \]
但要注意，预处理后矩阵 \(M^{-1} A\) 一般不再是复对称的。为了保持对称性，可采用分裂预处理形式，比如：
\[M = L L^T \quad \text{（不完全复对称分解）} \]
然后解方程：
\[L^{-1} A L^{-T} \tilde{x} = L^{-1} b, \quad x = L^{-T} \tilde{x} \]
此时 \(L^{-1} A L^{-T}\) 仍为复对称。在算法中应用预处理时，需在迭代中替换内积和矩阵向量乘。
预处理 COCR 算法称为 PCOCR，形式与 COCR 相似，但内积和向量需替换为预处理空间中的相应量。
总结
COCR 算法是一种针对复对称线性方程组设计的 Krylov 子空间迭代法，具有短递归、低存储、残差正交等特点，是 CG 方法在复对称情形的自然推广。通过合理选择内积和正交条件，保证了迭代的稳定性和有限步收敛性。在实际应用中，常结合预处理技术加速收敛。