双共轭梯度稳定法

双共轭梯度稳定法（biconjugate gradient stabilized method，BiCGSTAB）是求解一般非对称线性系统的短递推 Krylov 方法。它在双共轭方向更新后增加一个局部最小残差步，以减轻残差的剧烈振荡。它不需要保存整个 Krylov 正交基，也不需要显式乘以转置矩阵，但存在 breakdown 和非单调收敛的可能。

两个连续的残差修正

取 $A=\left(\begin{array}{cc}4& 1\\ 2& 3\end{array}\right)$、$b=(1,2{)}^T$，从 $x_0=0$ 开始。选固定影子残差 $\tilde{r}=r_0=b$，首方向 $p=r_0$。计算 $v=Ap=(6,8{)}^T$，并取 $$\alpha =\frac{{\tilde{r}}^T{r}_0}{{\tilde{r}}^{T}v}=\frac{5}{22},s=r_0-\alpha v=(-4/11,2/11{)}^T.$$ 这一步使 ${\tilde{r}}^{T}s=0$。接着 $t=As=(-14/11,-2/11{)}^T$，沿 $s$ 再作修正，取 $$\omega =\frac{t^{T}s}{t^{T}t}=\frac{13}{50}.$$ 完整一步为 $x_1=\alpha p+\omega s=(73/550,138/275{)}^T$，残差为 $r_1=s-\omega t=(-9/275,63/275{)}^T$。精确解是 $(1/10,3/5{)}^T$；第一步尚未到达它，第二步在精确算术下可完成。算法不是每轮简单重复最速下降，而是保留一组双共轭递推参数。

短递推的各个量

一般第 $k$ 轮先算 ${\rho }_k={\tilde{r}}^T{r}_k$。第一轮令 $p_0=r_0$，以后使用 $${\beta }_k=\frac{{\rho }_k}{{\rho }_{k-1}}\frac{{\alpha }_{k-1}}{{\omega }_{k-1}},p_k=r_k+{\beta }_k(p_{k-1}-{\omega }_{k-1}{v}_{k-1}).$$ 然后依次计算 $$v_k=A{p}_k,{\alpha }_k=\frac{{\rho }_k}{{\tilde{r}}^T{v}_k},s_k=r_k-{\alpha }_k{v}_k.$$ 若 $s_k=0$，当前半步已经得到精确解，不再除以后续的零范数。否则计算 $$t_k=A{s}_k,{\omega }_k=\frac{t_k^T{s}_k}{t_k^T{t}_k},x_{k+1}=x_k+{\alpha }_k{p}_k+{\omega }_k{s}_k,r_{k+1}=s_k-{\omega }_k{t}_k.$$ 这里的 $r_k$ 按递推表示残差，浮点运行中应定期直接计算 $b-A{x}_k$，防止递推残差与真实残差偏离。SciPy 的 BiCGSTAB 接口列明了一般系统、预条件与不收敛状态。

在残差平面中沿 t=As 的方向减去 ωt，使新残差与 t 垂直。它最小化这一条直线上的残差，不能保证相对上一完整迭代的范数每次都下降。

稳定化步骤的最小残差证明

固定非零 $t=As$，新残差为 $s-\omega t$。其平方范数展开为 $$\psi (\omega )=s^{T}s-2\omega t^{T}s+{\omega }^2{t}^{T}t.$$ 由于 $t^{T}t>0$，这是严格凸的二次函数。导数为零给出 $\omega =t^{T}s/(t^{T}t)$，最小值为 $$\Vert r_{new}{\Vert }_2^2=\Vert s{\Vert }_2^2-\frac{(t^{T}s{)}^2}{\Vert t{\Vert }_2^2}\le \Vert s{\Vert }_2^2.$$ 同时 $t^T{r}_{new}=0$。主例中这个正交关系可直接用分数核验。它是“稳定化”局部作用的精确含义；前一半步可能增大残差，所以不能把它推广成整条残差历史必定单调。

从多项式角度看，末步将中间残差再乘以 $I-\omega A$。连续轮次叠加这些一阶最小残差因子，改善原双共轭过程可能出现的振荡。完整收敛还依赖矩阵、影子向量、舍入和预条件，不由单个投影不等式保证。Netlib 给出相应算法和 breakdown 条件。

失败条件与成本

若非零残差却有 ${\rho }_k=0$ 或 ${\tilde{r}}^T{v}_k=0$，双共轭步无法定义；${\omega }_k=0$ 会破坏下一轮递推。极小分母也可能放大误差。这些情况不等于原系统无解，可尝试不同影子向量、适当预条件或换用 GMRES 等方法，并核验真实残差。

每个完整迭代需要两次矩阵向量积、若干内积和向量更新。稀疏系统每轮约为 $O(\mathrm{nnz}(A)+n)$，只保存常数个长度为 $n$ 的向量；与完整 GMRES 的扩展基存储形成区别。并行环境中，内积归约通信也占成本，不能只计算乘法次数。

可运行的 NumPy 演示

函数输入实矩阵、右端、正容限和最大轮数，从零开始，返回解和轮数，不修改输入。每轮验收真实残差；分母为零、非有限结果或轮数耗尽均明确报告失败。生产程序还会按尺度判断近 breakdown，并提供预条件接口。

Python 3

import numpy as np


def bicgstab_demo(A, b, tol=1e-12, maxiter=100):
    A, b = np.array(A, dtype=float, copy=True), np.array(b, dtype=float, copy=True)
    if b.ndim != 1 or not b.size or A.shape != (b.size, b.size):
        raise ValueError("square matrix and matching vector required")
    if not np.isfinite(A).all() or not np.isfinite(b).all() or not np.isfinite(tol) or tol <= 0:
        raise ValueError("finite data and positive tolerance required")
    if not isinstance(maxiter, int) or maxiter < 0:
        raise ValueError("nonnegative integer maxiter required")
    x = np.zeros_like(b); r = b.copy(); shadow = r.copy()
    norm_b = np.linalg.norm(b)
    if not np.isfinite(norm_b):
        raise ArithmeticError("initial residual norm outside floating-point range")
    target = tol*max(1.0, norm_b)
    if not np.isfinite(target):
        raise ArithmeticError("residual threshold outside floating-point range")
    p = np.zeros_like(b); v = p.copy(); rho_old = alpha = omega = 1.0
    for k in range(maxiter+1):
        if np.linalg.norm(b-A@x) <= target:
            return x, k
        if k == maxiter:
            break
        rho = np.dot(shadow, r)
        if rho == 0 or omega == 0:
            raise ArithmeticError("BiCGSTAB breakdown")
        beta = (rho/rho_old)*(alpha/omega)
        p = r+beta*(p-omega*v)
        v = A@p; denominator = np.dot(shadow, v)
        if denominator == 0:
            raise ArithmeticError("zero shadow denominator")
        alpha = rho/denominator; s = r-alpha*v
        candidate = x+alpha*p
        if np.linalg.norm(b-A@candidate) <= target:
            return candidate, k+1
        t = A@s; tt = np.dot(t, t)
        if tt == 0 or not np.isfinite(tt):
            raise ArithmeticError("zero or nonfinite stabilization norm")
        omega = np.dot(t, s)/tt
        x = candidate+omega*s; r = s-omega*t
        if not np.isfinite(x).all() or not np.isfinite(r).all():
            raise ArithmeticError("nonfinite iteration")
        rho_old = rho
    raise RuntimeError("true residual did not reach tolerance")


print(bicgstab_demo([[4, 1], [2, 3]], [1, 2]))  # ([0.1, 0.6], 2), up to rounding

有解变式取 $A=2I$。首方向为 $b$，$\alpha =1/2$，中间残差已经为零，程序应立即返回 $x=b/2$。若仍继续求 $\omega$，会出现无意义的 $0/0$；这个分支是算法正确性的一部分。