2
$\begingroup$

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)

$\endgroup$
0

2 Answers 2

6
$\begingroup$

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.

$\endgroup$
1
  • 2
    $\begingroup$ In particular, you don't want to explicitly form $(U^{-1}L^{-1})$ because this is likely to be a dense matrix! $\endgroup$ Nov 6, 2019 at 16:34
5
$\begingroup$

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)
$\endgroup$
1
  • 1
    $\begingroup$ I believe the Mx lambda definition is unnecessary. You can just do M = LinearOperator(A.shape, ilu.solve) $\endgroup$
    – AlexQueue
    May 12 at 17:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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