如何在 scipy 中实现 ILU 预处理器？

Question

对于scipy.sparse.linalg 中的迭代求解器，例如bicg、gmres 等，可以选择为矩阵A 添加前置条件。但是，documentation 不太清楚我应该作为预处理器给出什么。如果我使用ilu = sp.sparse.linalg.spilu(A)，ilu 不是任何矩阵，而是包含许多事物的对象。

有人问过类似的问题here 用于 Python 2.7，但我不适合我（Python 3.7，scipy 版本 1.1.0）

所以我的问题是如何为那些迭代算法合并不完整的 LU 预处理器？

Answer 1

作为前置条件，bicg 或 gmres 接受

• 稀疏矩阵
• 密集矩阵
• 线性算子

在您的情况下，预条件子来自分解，因此它必须作为线性运算符传递。

因此，您可能希望从通过spilu 获得的 ILU 分解中显式定义一个线性运算符。类似的东西：

sA_iLU = sparse.linalg.spilu(sA)
M = sparse.linalg.LinearOperator((nrows,ncols), sA_iLU.solve)

这里，sA 是 CSC 格式的稀疏矩阵，M 现在将是您将提供给迭代求解器的预处理器线性运算符。

---

一个完整的例子based on the question you mentioned:

import numpy as np
from scipy import sparse
from scipy.sparse import linalg

A = np.array([[ 0.4445,  0.4444, -0.2222],
              [ 0.4444,  0.4445, -0.2222],
              [-0.2222, -0.2222,  0.1112]])
sA = sparse.csc_matrix(A)

b = np.array([[ 0.6667], 
              [ 0.6667], 
              [-0.3332]])

sA_iLU = sparse.linalg.spilu(sA)
M = sparse.linalg.LinearOperator((3,3), sA_iLU.solve)

x = sparse.linalg.gmres(A,b,M=M)

print(x)

注意事项：

• 我实际上以密集矩阵为例，而在您的情况下，从具有代表性的稀疏矩阵开始会更有意义。
• 线性运算符M 的大小是硬编码的。
• ILU 无论如何都没有配置，而是使用默认值。
• 这与 in the comments to that aforementioned answer 的建议差不多，但是，我必须进行一些小调整以使其与 Python3 兼容。