# 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.) – Christian Clason Aug 2 '18 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. – GoHokies Aug 3 '18 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? – vydesaster Aug 3 '18 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. – amarney Aug 3 '18 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? – vydesaster Aug 3 '18 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.