2
注意：正如Warren Weckesser指出的那樣，我在最初發布的代碼中犯了一個愚蠢的錯誤。如果糾正，一些解決者給出正確的答案，但其他人給出NaN或不正確的答案。我也忘記在輸出中包含運行時警告;他們現在在那裏。我相應地修改了這個問題。我可能會使用解算器，但如果我明白其他人失敗的原因，我會更加高興。爲什麼scipy的稀疏求解器給出了不正確的答案？
我想解決使用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相同的線性系統的解，所以最後7行中的所有錯誤都應該爲零（或至少很小）。
gmres和lgmres工作，但其他人沒有。在大多數情況下，info正確指示失敗，但minres指示成功（info = 0），同時返回全零（不正確）的解決方案。
下面是一些潛在相關的其他信息：
- 矩陣A是對稱的，但不是正定的。
- 矩陣A是相當有條件的。
- 如果我用所有稀疏求解器中的原始矩陣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.dot（x1）'替換'A * x1'），或者在執行減法之前將'x1'平坦化'x2 - x1' 。 –
代碼可以複製粘貼進行測試，這很好。但是你的錯誤表明你沒有看到任何'x2'結果。你似乎馬上跳到與'x1'比較。 – hpaulj
謝謝！那真是愚蠢。正如你可能已經猜到的那樣，我對Python編程還很陌生。我應該明白這一點。 –