Correct use of scipy's sparse.linalg.spilu

Question

I'm attempting to use scipy's spilu routine as a preconditioner and I'm finding bad performance for my application (solving a global linear system arising from a DG discretization of an time-dependent ADR PDE).

Before I write tests to tune the options and compute the effect of the preconditioner on the eigenvalues of my linear system, I'd like to make sure I'm applying the preconditioning matrix in the scipy intended manner, since the documentation for spilu is pretty spartan...

elif solver == 'spILUPCG': P = sp.sparse.linalg.spilu(A.tocsc()) P = P.L * P.U x = sp.sparse.linalg.cg(A, b, M=P)

Answer 1

If $\mathbf L$ and $\mathbf U$ give an approximate factorization of $\mathbf A$, you wouldn't want to use $\mathbf P = \mathbf L\cdot \mathbf U$ as a preconditioner (that's approximately $\mathbf A$), you'd want to use $\mathbf P = \mathbf U^{-1}\cdot \mathbf L^{-1}$ (that's approximately $\mathbf A^{-1}$). Even then, you would not want to multiply them out explicitly into $\mathbf P$, it's probably dense. Instead think of $\mathbf P$ as an operation/operator that cascades two solve steps (first solve by $\mathbf L$, then solve by $\mathbf U$).

I would expect there to be some way to inject that operation as a closure/callback to precondition cg(). Sadly I'm not a scipy expert, so I don't know the exact way to do it. But somewhere scipy should have an abstraction for applying the action of a linear operator, and I would expect both their matrix objects and their preconditioner objects to fulfill that abstraction. I suppose I'd try just passing your $\mathbf P$ object itself.

Answer 2

This question has an example of how to create the preconditioner M with a scipy sparse matrix A of shape NxN

from scipy.sparse.linalg import LinearOperator, spilu
ilu = spilu(A)
Mx = lambda x: ilu.solve(x)
M = LinearOperator((N, N), Mx)