python-3.x - BicGStab 产生意外故障标志

Question

我需要求解级联的稀疏线性系统 Ax=b。第一个系统的解决方案 x 是第二个系统的输入，第二个系统是第三个系统的输入，依此类推。由于数值误差的叠加和其他原因，我必须使用 scipy.sparse.linalg.bicgstab 作为我的线性求解器。然而，对于一个甚至没有病态并且肯定有逆的系统，求解器给我一个标志:“非法输入或故障”。

import numpy as np
from scipy.sparse.linalg import bicgstab, inv
from scipy import sparse

A = np.array(
    [[ -1.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
     [  0.,  -1.,   0.,   0.,   0.,   0.,   0.,   0.],
     [  0.,   0., -10.,   0.,   0.,   0.,   0.,   0.],
     [  0.,   0.,   0., -10.,   0.,   0.,   0.,   0.],
     [  0.,   0.,   3.,   0.,  -3.,   0.,   0.,   0.],
     [  0.,   0.,   0.,   3.,   0.,  -3.,   0.,   0.],
     [  0.,   0.,   0.,   7.,   3.,   0., -10.,   0.],
     [  0.,   0.,   7.,   0.,   0.,   3.,   0., -10.]]
)
A = -sparse.csc_matrix(A)
b = np.array([ 1.,  0., 10.,  0.,  0.,  0.,  0.,  0.])
x, flag = bicgstab(A=A, b=b, maxiter=40, tol=1e-6)
x, flag

>>> (array([1.        , 0.        , 1.        , 0.        , 1.00118012,
        0.        , 0.3004875 , 0.70009946]), -10)

只是为了证明这一点

inv(A).dot(b)

>>> array([1. , 0. , 1. , 0. , 1. , 0. , 0.3, 0.7])

上面的输出正是我所期望的。有谁知道为什么 bicgstab 没有给我想要的输出？我找不到关于 bicgstab 的非法输入或故障的文档，因此我是关于 SO 的问题。

Answer 1

-10 错误代码并不一定意味着您输入错误；在您的情况下，故障很可能发生在迭代求解期间。

通过稍微改变你的 RHS:

b = np.array([ 1.,  0., 0.,  0.,  10.,  0.,  0.,  0.])

即使没有预处理器，scipy.bicgstab 也没有收敛问题:

x, flag = bicgstab(A=A, b=b, maxiter=40, tol=1e-6)
print (x, flag)

(array([1.        , 0.        , 0.        , 0.        , 3.33333333,
   0.        , 1.        , 0.        ]), 0)

矩阵有一个逆矩阵和一个合适的条件数

print(np.linalg.cond(A))

14.616397823169317

不保证很容易获得特定 RHS 的解决方案，尤其是使用迭代求解器或特定迭代求解器。在我看来(没有对矩阵谱及其内核空间进行详分割析)，您的 RHS 恰好位于这样一个“坏区域”。

如果您只是对解决方案感兴趣，我建议切换到 GMRES:

x, flag = gmres(A=A, b=b, maxiter=40, tol=1e-6)

(array([1. , 0. , 0.1 , 0. , 0.1 , 0. , 0.03, 0.07]), 0)

如果您有兴趣调查为什么 BiCGStab 失败，而 GMRES 成功解决了这个系统，我会邀请您将问题缩小到 Computational Science SE .