# Efficient implementation of preconditioners for iterative solvers

I am struggling a bit with the concept of preconditioners for iterative solvers and how to implement them efficiently. The literature mostly provides methods to create a preconditioner matrix $M$ (incomplete LU, incomplete Cholesky etc.). But the conjugate gradient methods need $M^{-1}$ during the iterations. As inverting a matrix is not a so good idea, I do not know how to get $M^{-1}$. I looked at some libraries that use iterative solvers and was able to see that the equation where $M^{-1}$ is needed, is actually solved by forward and backward substitution. Is this really pratical? When I think that I need to create the preconditioner and do forward and backward substitution, I just can use direct methods? Doing it that way I would save the time of all the iterations.

I also have read somewhere (unfortunately cannot find it anymore) that the LDL Cholesky decomposition can be easily inverted. The same literature also provides information that it might be better to solve $Mx=b$ instead of $Ax=b$ as $M$ should be very similar to $A$. I tried a couple of matrices but the results of $Mx=b$ are way different than $Ax=b$.

What do you think? What options are available for solving with preconditioners?

• Finding a good preconditioner $M$ means balancing two things: $M^{-1}A$ should be close to the identity and $M^{-1}$ should be cheap to apply. The extremes are $M=A$ and $M=I$, both of which are obviously not the best choice; the middle bear ("just right!") will very much depend on the structure of $A$ -- there's no "one size fits all"! (And your literature doesn't sound correct -- using a preconditioner means solving $M^{-1}Ax = M^{-1}b$, not $Mx=b$, so I'm not surprised it doesn't work.) Aug 2, 2018 at 22:58
• in addition, you do not need the explicit inverse $M^{-1}$. you just need to calculate its action on a vector, i.e. $M^{-1}z$ for some vector $z$. this is equivalent to solving the system $My = z$ (for $y$), which should be "cheap" to do. a good reference is the book of Saad - see chapters 9 and 10. Aug 3, 2018 at 5:36
• Thanks for the help and the link to the book. I read a few chapters in the book but I think I am still missing something. Assume I want to solve $My = z$ where M is the incomplete LU decomposition. As far as I was able to find an answer in the book, it is stated that the best way would be forward and backward substitution. Is this correct? If not, what am I missing? Aug 3, 2018 at 7:51
• Assume you have $M = LU$ where $L$ lower triangular and $U$ upper triangular, and you want to solve $My = z$. You can instead solve the problem $Lq=z$ for the right hand side of $Uy = q$, and the best way to do this is by forward substitution. Then you need to solve $Uy=q$ for $y$, and the best way to do this is by back substitution. Aug 3, 2018 at 13:33
• Just for my understanding. How can an iterative solver be faster when it needs to form the preconditioner and solve it with above mentioned methods compared to a direct solver which needs to factorize the system and solve this one? Aug 3, 2018 at 21:45

I don't particularly care for the notation $M^{-1}$ precisely because of the confusion you find yourself in. I (and others) simply call the preconditioner $P$.

The point, however, is that for iterative solvers, you generally never need a matrix or its inverse explicitly. All you need to be able to do is apply it to a vector. So, let's say you want to solve $$A x = b$$ and you've computed something like an incomplete LU factorization (ILU) of $A$, i.e., you have two triangular matrices $\tilde L,\tilde U$ so that $$A \approx \tilde L\tilde U$$ It is important to remember that computing an incomplete LU decomposition is substantially cheaper than computing the full LU decomposition because in general the LU factors of $A$ are dense.

Then call $P=(\tilde L\tilde U)^{-1}=\tilde U^{-1}\tilde L^{-1}$. In iterative solvers, you then have to "apply" this preconditioner, i.e., you will be given a vector $z$ and asked to return $y=Pz=\tilde U^{-1}\tilde L^{-1}z$. You may think that for this, you will need to know the inverses of $\tilde U,\tilde L$, but that's not true. You can compute it by doing the following steps:

• Compute $a=\tilde L^{-1}z$ by solving $\tilde La = z$, which is easy because $\tilde L$ is triangular and so this only requires (forward) substitution.
• Compute $y=Pz=\tilde U^{-1}\tilde L^{-1}z=\tilde U^{-1}a$ by solving $\tilde U y = a$, which is again easy because $\tilde U$ is also triangular and so only requires (backward) substitution.

So once you have $\tilde L, \tilde U$, computing $y=Pz$ is not actually that difficult. Now, whether you call the matrix $P=(\tilde L\tilde U)^{-1}$ and say that $P$ approximates $A^{-1}$, or whether you call it $M^{-1}=(\tilde L\tilde U)^{-1}$ and say that $M$ approximates $A$ does not make a difference -- that's just notation. In either case, we never need to actually know the entries or $P$ or $M$ or their inverses.

• Thank you very much for this detailed answers. Now it is clearer for me:) Aug 4, 2018 at 16:56
• My attempt at summarizing the answer: you can view the preconditioners as products of triangular matrices (except for the case of diagonal matrices). Then applying these preconditioners is a matter of solving triangular linear systems, which are straightforward and efficient for sparse matrices (assuming your code uses caching smartly). Sep 13, 2021 at 13:14
• @Breno That is not generally true. Preconditioners like Jacobi, SSOR, Gauss-Seidel can be seen as products of triangular matrices, but this is no longer true for things such as multigrid. Sep 13, 2021 at 19:42
• @WolfgangBangerth Thank you Wolfgang. I have just finished seeing yours lectures on this, on YouTube. Specifically on block (or physics-based) preconditioning. I would like to shamelessly request an answer to a question I made on this lecture, please. :) Sep 13, 2021 at 19:51