稀疏迭代求解器无矩阵方法预处理器

如何为无矩阵左侧的稀疏迭代方法（TFQMR、GMRES、CGS 等）定义预处理器（SPILU、SPAI 等）？我使用无矩阵 A（使用 LinearOperator 和 matvec）定义了 Ax=b。因此，我没有创建矩阵 A 并将其保存在内存中。例如，在这种情况下，我如何构建 SPILU 预处理器？我在所有教程和示例中看到预处理器是使用矩阵 LHS 创建的。

我使用无矩阵方法创建了左侧，如下所示：

def matvec(x, A, B):  
    First = B.dot(x)  
    Second = A.dot(First)  
    result = B.T.dot(Second)  
    return result

linop = LinearOperator((WWWa.shape[0], WWWa.shape[0]), 
                        matvec=lambda x: matvec(x, WWWa, Ap), 
                        rmatvec=lambda x: rmatvec(x, WWWa, Ap)
                       ) 

稀疏迭代求解器可以编写如下：|| |但是，当我想添加预处理器时，我无法弄清楚如何创建预处理器，因为我没有 LHS 的直接矩阵（而且我也不想获得矩阵；我有很大的稀疏矩阵）具有数十亿个元素的矩阵，我只想使用矩阵向量乘法）。

l_finest = A.dot(g) 
x_gmres, info_gmres = gmres(linop, l_finest.flatten(), callback=callback, tol=tolerance, atol=0,callback_type='x')

However, when I wanted to add the preconditioner, I could not figure out how to create the preconditioner since I don't have a direct matrix for the LHS (and I don't want to get a matrix either; I have large sparse matrices with billions of elements and I just want to use matrix-vector multiplication).

是对的，大多数预处理器（如 SPILU 和 SPAI）都需要对矩阵进行显式访问才能计算不完全分解或逼近。在不构建完整矩阵的情况下，可以在无矩阵设置中使用的一些预处理器方法包括：

1. 基于算子的预处理器：

• 多重网格（MG）： MG 方法特别适合于偏微分方程（PDE）产生的问题，它可以作为无矩阵预处理器工作。需要实现平滑算子（例如，加权 Jacobi 或 Gauss-Seidel）和粗网格层次结构，这些操作可以使用的矩阵向量乘积执行。
• 代数多重网格（AMG）： 与需要几何信息的 MG 不同，AMG 方法直接从矩阵构建粗网格层次结构。虽然经典的 AMG 需要矩阵，但也存在无矩阵 AMG 变体，它们仅依靠矩阵向量乘积，例如， 无矩阵代数多重网格 .
• 基于低秩近似的预处理器： 这些方法使用低秩矩阵近似原始算子，并且可以以无矩阵方式应用。示例包括基于随机抽样的方法或使用 Krylov 子空间方法（例如 Lanczos 或 Arnoldi 算法）。

2. 预处理矩阵向量乘积：

• 固定点迭代： 使用像 Jacobi 或 Gauss-Seidel 这样的简单迭代方法作为预处理器。每次预处理步骤都需要几次迭代，但这可以在没有显式矩阵的情况下完成。
• 近似逆预处理器： 这些方法尝试使用矩阵向量乘积构建 A 的逆的近似。示例包括因子化稀疏近似逆（FSAI）和基于 Krylov 子空间的方法。

实现示例（使用 Python 中的 PyAMG）：

import pyamg
import scipy.sparse.linalg as spla

# ... 定义的 matvec 函数和 LinearOperator ...

# 使用 AMG 作为预处理器
strength = 'symmetric'  # 选择 AMG 强度
smoother = 'jacobi'   # 选择 AMG 平滑器
presmoother = (smoother, {'omega': 1.0})
postsmoother = (smoother, {'omega': 1.0})
# 创建 AMG 预处理器
aml = pyamg.smoothed_aggregation_solver(linop, strength=strength, 
                                        presmoother=presmoother, 
                                        postsmoother=postsmoother)

# 创建预处理的线性算子
precond = spla.LinearOperator(linop.shape, aml.solve)

# 使用预处理器求解系统
x_gmres, info_gmres = spla.gmres(linop, l_finest.flatten(), M=precond,
                                callback=callback, tol=tolerance, atol=0, 
                                callback_type='x')

此示例演示了如何使用 PyAMG 库作为 GMRES 求解器的预处理器创建一个无矩阵代数多重网格（AMG）预处理器。

选择最佳预处理器取决于的特定问题和可用资源。通常值得尝试不同的方法以查看哪种方法最有效。