【发布时间】：2016-07-07 20:47:38
【问题描述】：
注意：正如 Warren Weckesser 所说，我在最初发布的代码中犯了一个愚蠢的错误。更正后，一些求解器给出正确答案，但其他求解器给出 NaN 或错误答案。我也忘了在输出中包含运行时警告；他们现在在那里。我已经相应地修改了问题。我也许可以使用有效的求解器，但如果我了解其他人失败的原因，我会更开心。
我正在尝试使用 scipy.sparse.linalg 中的一个或多个求解器来求解稀疏线性方程组。在系统小到可以直接求解的测试用例中，一些稀疏求解器给出的答案不正确，如下例所示：
import numpy as np
import scipy.sparse as ss
A = np.matrix([[ 0., 0., 0., 0., 0., 1., -1., -0., -0., -0., -0.],
[ 0., 0., 0., 0., 0., 2., -0., -1., -0., -0., -0.],
[ 0., 0., 0., 0., 0., 2., -0., -0., -1., -0., -0.],
[ 0., 0., 0., 0., 0., 2., -0., -0., -0., -1., -0.],
[ 0., 0., 0., 0., 0., 1., -0., -0., -0., -0., -1.],
[ 1., 2., 2., 2., 1., 0., -0., -0., -0., -0., -0.],
[-1., 0., 0., 0., 0., 0., -1., -0., -0., -0., -0.],
[ 0., -1., 0., 0., 0., 0., -0., -1., -0., -0., -0.],
[ 0., 0., -1., 0., 0., 0., -0., -0., -1., -0., -0.],
[ 0., 0., 0., -1., 0., 0., -0., -0., -0., -1., -0.],
[ 0., 0., 0., 0., -1., 0., -0., -0., -0., -0., -1.]])
b = np.matrix([0.,0.,0.,0.,0.,1.,0.,0.,0.,0.,0.]).T
As = ss.coo_matrix(A)
# The linear system Ax = b has a solution:
x1 = np.linalg.solve(A,b)
print("Solution to Ax = b:",x1)
print("Ax - b = ",A*x1-b)
print("Info and maximum error in solutions found by various other methods: ")
x2,info = ss.linalg.bicg(As,b)
print("bicg:",info,np.max(np.abs(x2-x1.ravel())))
x2,info = ss.linalg.bicgstab(As,b)
print("bicgstab:",info,np.max(np.abs(x2-x1.ravel())))
x2,info = ss.linalg.cgs(As,b)
print("cgs:",info,np.max(np.abs(x2-x1.ravel())))
x2,info = ss.linalg.gmres(As,b)
print("gmres:",info,np.max(np.abs(x2-x1.ravel())))
x2,info = ss.linalg.lgmres(As,b)
print("lgmres:",info,np.max(np.abs(x2-x1.ravel())))
x2,info = ss.linalg.minres(As,b)
print("minres:",info,np.max(np.abs(x2-x1.ravel())))
x2,info = ss.linalg.qmr(As,b)
print("qmr:",info,np.max(np.abs(x2-x1.ravel())))
当我运行它时，我得到以下输出：
Solution to Ax = b: [[ 0.07142857]
[ 0.14285714]
[ 0.14285714]
[ 0.14285714]
[ 0.07142857]
[-0.07142857]
[-0.07142857]
[-0.14285714]
[-0.14285714]
[-0.14285714]
[-0.07142857]]
Ax - b = [[ 0.00000000e+00]
[ 0.00000000e+00]
[ 2.77555756e-17]
[ 0.00000000e+00]
[ 1.38777878e-17]
[ 0.00000000e+00]
[ 2.77555756e-17]
[ -2.77555756e-17]
[ -5.55111512e-17]
[ 2.77555756e-17]
[ 0.00000000e+00]]
Info and maximum error in solutions found by various other methods:
bicg: 1 nan
bicgstab: 1 nan
cgs: 1 nan
gmres: 0 5.55111512313e-17
lgmres: 0 1.38777878078e-16
minres: 0 0.142857142857
qmr: -11 0.142857142857
/Users/ebunn/anaconda/lib/python3.5/site-packages/scipy/sparse/linalg/isolve/iterative.py:197: RuntimeWarning: invalid value encountered in multiply
work[slice2] *= sclr2
/Users/ebunn/anaconda/lib/python3.5/site-packages/scipy/sparse/linalg/isolve/iterative.py:318: RuntimeWarning: invalid value encountered in multiply
work[slice2] *= sclr2
/Users/ebunn/anaconda/lib/python3.5/site-packages/scipy/sparse/linalg/isolve/minres.py:244: RuntimeWarning: divide by zero encountered in double_scalars
Acond = gmax/gmin
x2 的每次计算都应该是与 x1 相同的线性系统的解，因此最后七行中的所有错误都应该为零（或至少很小）。
gmres 和 lgmres 有效，但其他无效。在大多数情况下，信息正确指示失败，但 minres 指示成功 (info=0)，同时返回全为零的解决方案（不正确）。
以下是一些可能相关的附加信息：
- 矩阵 A 是对称矩阵，但不是正定矩阵。
- 矩阵 A 条件良好。
- 如果我将稀疏表示 As 替换为所有稀疏求解器中的原始矩阵 A，则结果不会改变 - 也就是说，如果我说 ss.linalg.bicg(A,b)而不是 ss.linalg.bicg(As,b) 等
- 实际上包含 bicgstab 是不公平的，因为文档说这仅适用于正定矩阵。我将它包括在内是希望它可以工作，因为这种方法原则上适用于不定矩阵。
- 我正在运行 Python 3.5.1，scipy 版本 0.17.0。
当然，对于这个系统来说，没关系，因为直接解决它是微不足道的。对于稀疏性至关重要的大型问题，这是一个热身问题。
gmres 和/或 lgmres 可能会满足我的需求，但我仍然想了解其他人出了什么问题，部分是为了我自己的理解，但也可以让我有更广泛的方法来从中选择。特别是，这两种有效的方法都没有利用 A 的对称性，如果有一个可以利用的方法可能会很好。
【问题讨论】：
-
小心！因为
A是np.matrix，所以x1也是np.matrix，形状为(11, 1)。另一方面，x2是一个形状为 (11,) 的np.ndarray，因此当您计算x2 - x1时，应用广播并且结果为 (11, 11)。您可以通过多种方式避免这种情况；例如，使用A = np.array(...)（但随后将A*x1替换为A.dot(x1)），或者在进行减法x2 - x1之前将x1展平。 -
很高兴代码可以复制粘贴进行测试。但是您的错误表明您实际上没有查看任何
x2结果。您似乎马上就开始与x1进行比较。 -
谢谢！那是我的愚蠢。正如您可能已经猜到的那样，我对 Python 编程还很陌生。我应该知道的。
标签： python numpy scipy sparse-matrix