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)


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.

| cite | improve this answer | |
  • 1
    $\begingroup$ In particular, you don't want to explicitly form $(U^{-1}L^{-1})$ because this is likely to be a dense matrix! $\endgroup$ – Brian Borchers Nov 6 '19 at 16:34
  • $\begingroup$ Good clarification, thank you. $\endgroup$ – rchilton1980 Nov 6 '19 at 20:29

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)
| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.