1
$\begingroup$

I was reading these slides about preconditioners. I believe I grasp the idea of how they work but there is something that is still not making sense.

If we have the system $Ax=b$ and use a preconditioner we end up with a system $M^{-1}Ax = M^{-1}b$.

Now say I decide to use an $LU$ factorization, it means that I need to compute $M = LU$ and to use it as preconditioner I need to invert it.

Here is where I am a little bit missing, does this mean that we are actually computing the inverse for $M$? I guess there is some property of lower/upper matrices that makes the inverse easy but I can't seem to find it. If we are actually computing the inverse and no property of lower/upper matrices makes the task easy, why do we bother to compute the inverse of the preconditioner instead of computing the inverse (or approximate the inverse) of $A$?

$\endgroup$
  • 2
    $\begingroup$ the key observations are (1) $M^{-1} A$ is (much) better conditioned than $A$, so the iterative method will converge much faster than in the un-preconditioned case (where we just "invert" $A$); and (2) $M$ should be cheap (i.e., computationally "easy") to invert, i.e., the $M = LU$ is cheap to compute and invert. The simplest example is (block-)Jacobi preconditioning, where $M$ is (block-)diagonal thus trivial to invert. As to how $M$ is applied, that is shown in slides 5-6 of the notes you linked to. $\endgroup$ – GoHokies May 12 '17 at 12:23
  • $\begingroup$ Just to make sure you are clear on this, you never actually compute an inverse of $M$. Once you have $M=LU$, you simply use the $L$ and $U$ factors to solve for the relevant vectors. $\endgroup$ – Bill Greene May 12 '17 at 12:57
  • $\begingroup$ @BillGreene I think that is what messes with my head, so if we do not compute the inverse how do we actually apply it when solving the system? Does it depend on the iterative solver? Any resource you can point me to where I it explains how the LU (without invert) is used to solve the problem? $\endgroup$ – BRabbit27 May 12 '17 at 13:05
  • $\begingroup$ @GoHokies I just see $L^{-1}$ and $U^{-1}$ I can't see how using the non-inverted versions come into the solver :/ $\endgroup$ – BRabbit27 May 12 '17 at 13:07
  • 1
    $\begingroup$ @BRabbit27 whenever you see $L^{-1}$ applied to a vector $y$ read it as "solve $Lx = y$ (for $x$)". A linear solve involving the matrix $L$ is "cheap" because $L$ is lower triangular (hence the total cost is quadratic in the number of rows/columns of $L$). Same goes for $U$ (upper triangular). $\endgroup$ – GoHokies May 12 '17 at 13:19
4
$\begingroup$

OK - so your original equation is $$Ax = b$$

Say you've come up with a good preconditioner for $A$, call this $M$. Also, say you have pre-computed an LU-decomposition for this $M$, i.e.

$$M = L_m U_m$$

where I've used the subscript $m$ to indicate that it is the LU decomposition of $M$ that we're talking about here (and not that of $A$). Then, assuming left preconditioning, you will be solving the preconditioned linear system

$$ M^{-1} A x = M^{-1}b,$$

or, equivalently

$$ (L_m U_m)^{-1} A x = (L_m U_m)^{-1} b,$$

or

$$ U_m^{-1} L_m^{-1} A x = U_m^{-1} L_m^{-1} b.$$

Thus the preconditioned system matrix is now $U_m^{-1} L_m^{-1}A$, which is (if you have chosen a good preconditioner) much better conditioned than $A$.

Now, you can

(a) (as a purely academic exercise,) solve the system with a non-preconditioned Krylov method (say GMRES). Because the system matrix is well conditioned, you can expect quick convergence.

The Krylov algorithm will need to multiply this matrix to an arbitrary (known) vector, call that $r$. In other words, what the Krylov procedure needs from you is a routine that, given $r$, computes $$z = U_m^{-1} L_m^{-1} A r$$

This is equivalent to solving the linear system $L_m U_m z = A r$ (for $z$). This system is easy (=cheap) to solve because the matrices $L_m$ and $U_m$ are triangular.

or

(b) (what you'd do in practice) run preconditioned CG on the preconditioned system. CG will apply the inverse preconditioner to the residual vector $r := b - Ax$. Again, this just means it solves the linear system $L_m U_m z = r$ (for $z$).

So, to summarize: never invert matrices explicitly (you can't do that in practice for large matrices anyway), instead solve the corresponding linear system (which is mathematically equivalent to applying the inverse matrix operator to a known vector).

$\endgroup$
  • $\begingroup$ Actually the preconditioner routine should just compute $U_m^{-1} L_m^{-1} r$ (without the $A$). The $A$ comes in implicitly because CG will apply this routine to a residual. $\endgroup$ – cfh May 16 '17 at 18:16
  • $\begingroup$ But we have no reason to believe that $M^{-1} A$ is symmetric, so CG won't work if you apply it directly to that. $\endgroup$ – cfh May 16 '17 at 18:25
  • 1
    $\begingroup$ @cfh thanks for spotting this. i've edited the answer. $\endgroup$ – GoHokies May 16 '17 at 18:43
  • $\begingroup$ The fact that the preconditioner is an LU factorization (and not Cholesky) suggests to me that $A$ is not symmetric positive definite. If that were the case, then PCG would not be applicable. $\endgroup$ – Richard Zhang May 16 '17 at 22:38
4
$\begingroup$

Based on your comments, it looks like your main difficulty is how having an LU factorisation for a given matrix $M$ makes it easy to find a vector $\mathbf{y}$ such that $M\mathbf{y}=\mathbf{b}$ for a given $\mathbf{b}$. Hopefully you're comfortable that if we have an easy method to do that, then this bit of the problem is solved, since we can treat $M^{-1}$ as a black box, where we're given an input and the method gives back an output.

Let's suppose for a second that we (somehow) know how to find a $\mathbf{z}$ such that $L\mathbf{z} = \mathbf{b}$ for a given $\mathbf{b}$. Then providing we can find a $\mathbf{y}$ which gives $U\mathbf{y} = \mathbf{z}$, then $$LU\mathbf{y} = L\mathbf{z} = \mathbf{b}, $$ and we're done. But how to find $\mathbf{y}$? Since $U$ is upper triangular, then writing things as simultaneous equations, the last equation is just $$ U_{n,n}y_n = z_n $$ This is easy to solve in one go as $y_n = \frac{z_n}{U_{n,n}}$. But now the equation one back from the last can be written as $$ U_{n-1,n-1}y_{n-1} = - U_{n-1,n}y_n + z_{n-1},$$ which again is a one liner to solve, given that we've just found $y_n$ (this method is called backwards substitution if you want to look it up). Since getting the i-th component means doing $n-i$ multiplications, and one division, and there are $n$ lines, we're doing about $n^2$ operations, which is better than general matrix inversion. Even better, if $U$ actually has fewer entries than a general upper triangular matrix, then we have even less maths to do.

So now we need to go back to finding $\mathbf{z}$. But the first equation in $L\mathbf{z}=\mathbf{b}$ is just $$ L_{1,1} z_1 = b_1, $$ so we can play exactly the same game we played before, except using forward substitution.

$\endgroup$
  • $\begingroup$ Uhm ... I think this is not what I'm missing. So, suppose we have the following equation to solve $Ax = b$ where we are solving for $x$. I know it is possible to find an $LU$ factorization for $A$, i.e. $A = LU$ so that we can transform the original problem to $LUx = b$. That I understand and know how to solve. Is this what you tried to describe? $\endgroup$ – BRabbit27 May 13 '17 at 12:40
  • $\begingroup$ Now, when using an iterative solver to solve for $x$ in the equation $Ax = b$ we can use a preconditioner such that $M^{-1}Ax = M^{-1}b$. In the slides I saw it is possible to use the $LU$ factorization as preconditioner, i.e. $M = LU$, so the iterative solver will solve the following $(LU)^{-1}Ax = (LU)^{-1}b$ but in the commets above someone said we do not actually compute the inverse, therefore, how to use it in an iterative solver? $\endgroup$ – BRabbit27 May 13 '17 at 12:43
  • $\begingroup$ @BRabbit27 We do not compute the inverse. We play smart—we define how $M^{-1}$ acts on a vector. For example, if one uses (I)LU decomposition as a preconditioner, one defines the action of $M^{-1}$ via forw and back substitutions. $\endgroup$ – 56th May 13 '17 at 21:15
  • 1
    $\begingroup$ @BRabbit27 just to restate the same thing in different language, inverse can either mean "find the matrix which gives $M^{-1}M=I$" or "the vector function where you tell me the $\mathbf{b}$ and I'll tell you an $\mathbf{x}$ such that $M\mathbf{x}=\mathbf{b}$". When people are saying we don't do the inverse, they mean we don't do the first one. When you see $M^{-1}$ turning up in your iterative method, you do do the second one. $\endgroup$ – origimbo May 13 '17 at 22:58

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.