Efficient implementation of preconditioners for iterative solvers

Question

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?

Answer 1

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.