双共轭梯度法解线性方程组（BiCG）

                     

贡献者： addis

1. 共轭梯度法

　　¹要求解对称正定矩阵（symmetric and positively defined, SPD） $\mathbf{A}$，的线性方程组

只需要令 容易证明 $\mathbf{\nabla }f=\mathbf{A}\mathbf{x}-\mathbf{b}$。注意 $f(x)$ 是一个凹二次函数，所以取最小值当且仅当梯度为零。这样，解方程组的问题就转化为求函数极小值问题。我们可以用梯度法来求最小值，即从出发点 ${\mathbf{x}}_0$ 开始，在梯度方向搜索函数最小值的位置 ${\mathbf{x}}_1$，再在其梯度方向搜索最小值的位置 ${\mathbf{x}}_2$……

　　 同理，若 $\mathbf{A}$ 是一个厄米矩阵且正定，把式 2 变为

同样可以证明 $f$ 取最小值时式 1 成立。注意此时梯度不太好定义，那么我们可以分别令 即可证明。

　　 要找现成的程序，Eigen 提供 Eigen::ConjugateGradient，适用于对称或者厄米矩阵。密矩阵和稀疏矩阵都可以。

2. 双共轭梯度法

　　 梯度法可以拓展为双共轭梯度法，以适用于任意方矩阵的线性方程组。该方法的优势在于用户只需要向解算器提供矩阵乘矢量的函数，而不需要提供矩阵本身。这样，矩阵可以具有任意的数据结构，例如各种稀疏矩阵。

　　 当矩阵 $\mathbf{A}$ 接近于单位矩阵时，该方法收敛更快，因此，我们可以选择不直接求解式 1 而是求解

其中 $\tilde{\mathbf{A}}$ 和 $\mathbf{A}$ 接近，但更易于求解。这样就有 ${\tilde{\mathbf{A}}}^{-1}\mathbf{A}\approx \mathbf{I}$。这里 $\tilde{\mathbf{A}}$ 通常称为 preconditioner。该方法称为 preconditioned biconjugate gradient method (PBCG)。如果你找不到更好的 preconditioner，通常可以用 $\mathbf{A}$ 的对角线充当。若选择不使用 preconditioner，也可以直接令 $\tilde{\mathbf{A}}$ 为单位矩阵。

　　 PBCG 的 C++ 代码见 [1] 的 linbcg.h。另外 C++ 的 Eigen 库也提供 Eigen::BiCGSTAB 算法。Matlab 也有 x = bicgstab(A,b) 函数（支持复数）。其中 STAB 意思是 stablized，单纯的 BiCG 据说并不稳定，见 Wikipedia 相关页面（Matlab 的 bicgstab 函数和这里给出的伪代码非常相似）。

---

1. ^ 本文参考 [1]。

---

[1] ^ W. H. Press, et al. Numerical Recipes 3rd edition.