scipy.sparse.linalg.spilu(A).solve() 函数使用什么算法？答案

【问题标题】：What algorithm does the scipy.sparse.linalg.spilu(A).solve() function use?scipy.sparse.linalg.spilu(A).solve() 函数使用什么算法？
【发布时间】：2020-04-20 15:22:31
【问题描述】：

考虑以下代码：

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），这似乎很愚蠢。

【问题讨论】：

你真的很深入，但如果你想追求这个 - 你可以看看源代码。您使用的模块调用： return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr, ilu=True, options=_options) gstrf 是 superlu 包中的一种方法（它不在 python 中) - 您可以尝试在他们的文档中查看他们如何计算 ILU。

标签： python scipy

【解决方案1】：

Docs说：

此函数使用 SuperLU 库。

Code经过：

return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
                      csc_construct_func=csc_construct_func,
                      ilu=True, options=_options)

which is 只是一个 SuperLU 包装器：

static char gstrf_doc[] = "gstrf(A, ...)\n\
\n\
performs a factorization of the sparse matrix A=*(N,nnz,nzvals,rowind,colptr) and \n\
returns a factored_lu object.\n\
\n\
arguments\n\
---------\n\
\n\
Matrix to be factorized is represented as N,nnz,nzvals,rowind,colptr\n\
  as separate arguments.  This is compressed sparse column representation.\n\
\n\
N         number of rows and columns \n\
nnz       number of non-zero elements\n\
nzvals    non-zero values \n\
rowind    row-index for this column (same size as nzvals)\n\
colptr    index into rowind for first non-zero value in this column\n\
          size is (N+1).  Last value should be nnz. \n\
\n\
additional keyword arguments:\n\
-----------------------------\n\
options             specifies additional options for SuperLU\n\
                    (same keys and values as in superlu_options_t C structure,\n\
                    and additionally 'Relax' and 'PanelSize')\n\
\n\
ilu                 whether to perform an incomplete LU decomposition\n\
                    (default: false)\n\
";

看看最后一个参数 -> ilu！

splu 看起来像：

return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
                      csc_construct_func=csc_construct_func,
                      ilu=False, options=_options)

表示将整个完整与不完整的逻辑传递给 SuperLU。

现在让我们看看SuperLU's manual：

2.7 不完全 LU 分解 (ILU) 预条件子

从 SuperLU 4.0 版开始，我们提供 ILU 例程用作预处理器用于迭代求解器。我们的 ILU 方法可以被认为是 ILUTP 方法的变体最初由 Saad [31] 提出，它结合了双重丢弃策略和数值旋转（“T”代表阈值，“P”代表枢轴）。

参考文献[31]为：

萨阿德，优素福。 “ILUT：双阈值不完全 LU 分解。”数值线性代数与应用 1.4 (1994): 387-402.

这对我来说似乎是一种直接方法 (Eigen does use it too)，但我认为您将能够查找您需要知道的内容。

【讨论】：

【解决方案2】：

这种求解方法近似求解Ax=b通过使用A的不完全分解。这就是为什么您看到它在那些迭代方法中用作预处理器。求解方法是使用不完整因子的经典稀疏前/后替换算法，但由于它们不完整，因此不太可能给出准确的结果。

【讨论】：