Preconditioner Implementation with matrix-free methods (sparse iterative solvers)

Question

How can I define preconditioners (SPILU, SPAI, etc.) for sparse iterative methods (TFQMR, GMRES, CGS, etc.) for the matrix-free left-hand side? I defined $Ax=b$ using matrix-free $A$ (with LinearOperator and matvec). Therefore, I do not have a matrix $A$ created and kept in memory. For example, in this case, how can I build SPILU preconditioners? I investigated all the tutorials and examples that a preconditioner was created using a matrix LHS, they used matrix $A$, not linear operator.

I created the left-hand side with matrix-free methods as follows:

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)
                       )

Sparse iterative solver can be written as follows:

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).

People use this approach to build preconditioner:

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

However, my left-hand side is not a matrix and it is LinearOperator, linop (matrix-free multiplications), so it is not a square matrix. I think there must be a way to implement a preconditioner in this case, but I couldn’t find how.

Answer 1

Loosely, you can think of a preconditioner as an approximation to the matrix's inverse. (More specifically, it needs to improve various properties such as condition number and eigenvalue distribution, depending on the specific iterative method. But considering approximations to the matrix's inverse is a reasonable starting step.) For example, ILU is basically an approximation of LU factorization in which accuracy is sacrificed for the sake of speed.

So, to find a good preconditioner you need to know something about your matrix. This is why the common preconditioner suggestions require the matrix elements: generic algorithms don't have any other way to get information about the matrix. Ultimately, you'll have to study your matrix and the underlying problem to find a good preconditioner.

To help show you how to get started, I can offer a few comments. Your matvec function indicates that your matrix is decomposed as $L = BAB^T$, and the LinearOperator's size indicates that $A$ and $B$ are square. This tells us that $L^{-1} = B^{-T}A^{-1}B^{-1}$, and so we can approximate $L^{-1}$ by approximating $B^{-T}$, $A^{-1}$, and $B^{-1}$. Since it looks like WWWa and Ap are explicitly stored, you could just do an ILU or SPAI on each of them, then compose a matrix-free preconditioner similar to your matrix-free operator. But you might be able to do even better with a little analysis and experimentation. (ILU and SPAI also have parameters to experiment with.)

• Does $A$ or $B$ have a nice structure that allows for efficient direct solves (i.e., triangular, banded, etc.)?
• Is $A$ or $B$ orthogonal? Then, its inverse is just its transpose.
• Does your underlying physics/chemistry/etc. tell you anything about $A$ or $B$?
  • Can you manipulate the underlying physics/etc. to get a simple approximation? Something like using Fourier space, or replacing relativistic physics with Newtonian physics.
  • Are there separate domains or regions in the problem? Domain decomposition methods might be helpful.
• Search the literature for ideas, even if no one has suggested something for your exact problem, you might get ideas from other problems.