4
$\begingroup$

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?

$\endgroup$
5
  • 6
    $\begingroup$ 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.) $\endgroup$ Aug 2 '18 at 22:58
  • 4
    $\begingroup$ 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. $\endgroup$
    – GoHokies
    Aug 3 '18 at 5:36
  • $\begingroup$ 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? $\endgroup$
    – vydesaster
    Aug 3 '18 at 7:51
  • $\begingroup$ 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. $\endgroup$
    – amarney
    Aug 3 '18 at 13:33
  • $\begingroup$ 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? $\endgroup$
    – vydesaster
    Aug 3 '18 at 21:45
7
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Thank you very much for this detailed answers. Now it is clearer for me:) $\endgroup$
    – vydesaster
    Aug 4 '18 at 16:56
  • $\begingroup$ 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). $\endgroup$
    – Breno
    Sep 13 at 13:14
  • 1
    $\begingroup$ @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. $\endgroup$ Sep 13 at 19:42
  • 1
    $\begingroup$ @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. :) $\endgroup$
    – Breno
    Sep 13 at 19:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.