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

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.

• ILU preconditioners typically require at least the sparsity pattern of $A$ or some power of $A$. Is this something that is easy to access for your problem? Commented Jul 24 at 5:32
• There is a vast literature on building preconditioners for matrix-free methods. What have you already looked at? Commented Jul 24 at 16:47
• The point, of course, of using matrix-free methods is that you want to avoid the memory cost and the data transfer of matrices. Using a preconditioner that requires you to store matrix elements (e.g., for an ILU decomposition) defeats that purpose. If you want to store an ILU, you might of course just store the original matrix as well. Commented Jul 24 at 16:48
• All of the factorization-based preconditioners are ill-suited to matrix-free implementations for the reasons Wolfgang mentioned. Some "easy" preconditioners like Jacobi can be implemented with only access to some of the matrix elements, but going completely matrix-free for everything requires a lot of problem-specific information and there are no black box options available typically. Commented Jul 24 at 19:14

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$$?