【问题标题】:BicGStab yields unexpected breakdown flagBicGStab 产生意外故障标志
【发布时间】:2020-03-25 17:46:35
【问题描述】:

我需要解决一系列稀疏线性系统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 的非法输入或故障的文档,因此我是我的问题。

【问题讨论】:

    标签: python-3.x numpy scipy sparse-matrix linear-algebra


    【解决方案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)
    

    (数组([1. , 0. , 0.1 , 0. , 0.1 , 0. , 0.03, 0.07]), 0)

    如果您有兴趣调查 BiCGStab 失败的原因,而 GMRES 在该系统的解决方案中成功,我会将您缩小的问题邀请至Computational Science SE

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 2019-10-27
      • 2023-03-09
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多