双共轭梯度法解线性方程组(BiCG)

                     

贡献者: addis

  • 本文处于草稿阶段。
预备知识 多元函数的极值线性方程组解的结构厄米矩阵、自伴矩阵

1. 共轭梯度法

  1要求解对称正定矩阵(symmetric and positively defined, SPD) $ \boldsymbol{\mathbf{A}} $,的线性方程组

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

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

\begin{equation} f( \boldsymbol{\mathbf{x}} ) = \frac{1}{2} \boldsymbol{\mathbf{x}} ^\dagger \boldsymbol{\mathbf{A}} \boldsymbol{\mathbf{x}} ~, \end{equation}
同样可以证明 $f$ 取最小值时式 1 成立。注意此时梯度不太好定义,那么我们可以分别令
\begin{equation} \frac{\partial f}{\partial \operatorname{Re} [x_i]} = 0~, \qquad \frac{\partial f}{\partial \operatorname{Im} [x_i]} = 0~ \end{equation}
即可证明。

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

2. 双共轭梯度法

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

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

\begin{equation} (\tilde { \boldsymbol{\mathbf{A}} }^{-1} \boldsymbol{\mathbf{A}} ) \boldsymbol{\mathbf{x}} = \tilde { \boldsymbol{\mathbf{A}} }^{-1} \boldsymbol{\mathbf{b}} ~. \end{equation}
其中 $\tilde { \boldsymbol{\mathbf{A}} }$ 和 $ \boldsymbol{\mathbf{A}} $ 接近,但更易于求解。这样就有 $\tilde { \boldsymbol{\mathbf{A}} }^{-1} \boldsymbol{\mathbf{A}} \approx \boldsymbol{\mathbf{I}} $。这里 $\tilde { \boldsymbol{\mathbf{A}} }$ 通常称为 preconditioner。该方法称为 preconditioned biconjugate gradient method (PBCG)。如果你找不到更好的 preconditioner,通常可以用 $ \boldsymbol{\mathbf{A}} $ 的对角线充当。若选择不使用 preconditioner,也可以直接令 $\tilde { \boldsymbol{\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.

致读者: 小时百科一直以来坚持所有内容免费无广告,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 20 元,我们一周就能脱离亏损, 并在接下来的一年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利