请考虑以下代码:
import numpy as np
import scipy.sparse.linalg
# Setup
A = scipy.sparse.csc_matrix([[1, -3], [-1, 4]])
b = np.array([1, 0])
spilu = scipy.sparse.linalg.spilu(A) # Find ILU decomposition
x = spilu.solve(b) # Use an iterative method to solve Ax = b ?
这导致x
是Ax = b
的解。
我的理解是scipy.sparse.linalg.spilu(A)
计算的是A的不完全LU分解。
我的问题是:spilu.solve(b)
使用什么算法来求解Ax = b
?
我希望它使用一些迭代方法,如共轭梯度法,因为这似乎是使用不完全LU分解的正常方式。然而,我一直无法找到支持或不同意这一点的文档。此外,我感到困惑,因为我看到一些人将LinearOperator
和scipy.sparse.linalg.cg
/scipy.sparse.linalg.cg
与scipy.sparse.linalg.spilu(A)
结合使用。如果我的假设是正确的(for example),那么这看起来是愚蠢的。
3条答案
按热度按时间wecizke31#
医生说:
此函数使用SuperLU库。
代码通过:
它只是一个SuperLU Package 器:
"看最后一个参数-〉ilu!"
splu看起来像:
其指示整个完整的逻辑与不完整的逻辑被传递到SuperLU。
现在,我们来看看SuperLU's manual:
2.7不完全LU分解(ILU)预条件子
从SuperLU版本4.0开始,我们提供了ILU例程,用作迭代求解器的预处理器。我们的ILU方法可以被认为是最初由Saad [31]提出的ILUTP方法的变体,它结合了双重丢弃策略和数值主元(“T”代表阈值,“P”代表主元)。
参考文献[31]是:
萨阿德,优素福。“ILUT:双阈值不完全LU分解〉,《数值线性代数与应用》1.4(1994):387-402.
在我看来这是一个DIRECT方法(Eigen does use it too),但我认为您将能够查找您需要了解的内容。
lqfhib0f2#
此求解方法 * 通过使用
A
的不完全因式分解近似求解 *Ax=b
。这就是为什么您会看到它在那些迭代方法中用作预处理器的原因。此求解方法是使用不完全因式的经典稀疏前向/后向替换算法,但由于它们是不完全的,因此不太可能给予准确的结果。70gysomp3#
萨莎补充道:
_superlu.gstrf
返回的对象的solve
函数可以在scipy/sparse/linalg/_dsolve/_superluobject. c中找到。主要部分是对gstrs
的调用(请注意,它是gstr**s**
而不是gstr**f**
):查看documentation of SuperLU:
所以Python代码基本上是调用
*gstrf
来计算LU分解,然后把它给*gstrs
来计算一个系统的解。看一下
dgstrs
的源代码,看起来它只是使用了LU分解,不管它是否是不完整的(也没有参数允许指定),所以它只是用不完整的分解给予你一个近似的解。