Lecture 15: Iterative Methods and Preconditioners | Video Lectures | Mathematical Methods for Engineers II | Mathematics

Download this transcript - PDF (English - US)

NARRATOR: The following content is provided by MIT OpenCourseWare under a Creative Commons license. Additional information about our license and MIT OpenCourseWare in general is available at ocw.mit.edu.

PROFESSOR: So I'm thinking of large sparse matrices. So the point about a -- the kind of matrices that we've been writing down, maybe two dimensional, maybe three dimensional difference matrices. So it might be five or seven non 0s on a typical row and that means that you don't want to invert such a matrix. These are big matrices. The order could easily be a million. What you can do fast is multiply a by a vector, because multiplying a sparse matrix by a vector, you only have maybe five non 0s times the n -- the vector of length n, so say five n operations -- so that's fast. I'm thinking that we are at a size where elimination, even with a good ordering, is too expensive in time or in storage. So I'm going from direct elimination methods to iterative methods. Iterative means that I never actually get to the perfect answer x, but I get close and what I want to do is get close quickly.

So a lot of thought has gone into creating iterative methods. I'll write down one key word here. Part of the decision is to choose a good preconditioner. I'll call that matrix p. So it's a matrix that's supposed to be -- that we have some freedom to choose and it very much speeds up, so the decisions are, what preconditioner to choose? A preconditioner that's somehow like the matrix a, but hopefully a lot simpler. There are a bunch of iterative methods and that's what today's lecture is about.

Multigrid is an idea that just keeps developing. It really is producing answers to very large problems very fast, so that will be the following lectures and then come these methods -- you may not be familiar with that word Krylov. Let me mention that in this family, the best known method it is called conjugate gradience. The conjugate gradient, so I'll just write that down. Conjugate gradient, so that will come next week. That's a method that's a giant success for symmetric positive definite problem. So I don't assume in general that a is symmetric or positive definite. Often it is if it comes from a finite difference or finite element approximation. If it is, then conjugate gradience is highly recommended.

I'll start with the general picture. Let me put the general picture of what an iteration looks like. So an iteration computes the new x. Let me call it x-new or xk plus 1 from the known x. So we start with x 0, any x 0. In principle, it's not too important for 0 to be close in linear methods. You could start even at 0. In nonlinear methods, Newton methods, you do want to be close and a lot of attention goes into a good start. Here the start isn't too important; it's the steps that you taking. Let me -- so I guess I'm going to -- I'm trying to solve ax equal b. I'm going to split that equation into -- I'm just going to write the same equation this way. I could write it as x equals i minus ax plus b. That would be the same equation, right? Because x and x cancel. ax comes over here and I have ax equal b. I've split up the equation and now actually, let me bring -- so that's without preconditioner. You don't see a matrix p. That's just an ordinary iteration, no preconditioner.

The preconditioner will come. I'll have p and it'll go here too. Then of course, I have to change that to a p. So that's still the same equation, right? px is on both sides, cancels -- and I have ax equal b. So that's the same equation, but it suggests the iteration that I want to analyze. On the right hand side is the known approximation. On the left side is the next approximation, xk plus 1. I multiply by p, so I -- the idea is p should be close to a. Of course, if it was exactly a -- if p equalled a, then this wouldn't be here and I would just be solving in one step ax equal b. So p equal a is one extreme. I mean, it's totally, perfectly preconditioned, but the point is that if p equals a, we don't want to solve that system. We're trying to escape from solving that system, but we're willing to solve a related system that is simpler. Why simpler?

Here are some possible preconditioners. p equals the diagonal. Just take the diagonal of a and that choice is named after Jacobi. That's Jacobi's method. Actually, it's already a good idea, because the effect of that is somehow to properly scaled the equations. When I had the identity here, I had no idea whether the identity and a are in the right scaling. a may be divided by delta x squared or something. It could be totally different units, but p -- if I take the diagonal of a, at least they're comparable. Then you see what -- we'll study Jacobi's method. So that's a very simple choice. Takes no more work than the identity, really. Second choice would be -- so that would be p for Jacobi. A second choice that you would think of pretty fast is the diagonal and the lower triangular part. So it's the -- I'll say the triangular -- I'm going to say lower triangular, triangular part of a, including the diagonal. This is associated with Gauss's great name and Seidel. So that's the Gauss-Seidel choice. So why triangular? We'll analyze the effect of this choice. If p is triangular, then we can solve -- I'm never going to write out an inverse matrix. That's understood, but I can compute the new x from the old x 1 just by back substitution. I find the first component of x, k plus 1, then the second, then the third. It's just because the first equation -- if p is triangular, the first equation will only have one term, one new term on the left side and everything else will be over there. The second equation will have a diagonal and something that uses the number I just found for the first component and onward. So in other words, triangular matrices are good.

Then people thought about combinations of these. A big lot of activity went in something called overrelaxation, where you did a -- you had a weighting of these and it turned out to be if you adjusted the weight, you could speed up the method, sometimes called SOR-- successive overrelaxation. So when I was in graduate school, that was a very -- that was sort of the beginning of progress beyond Jacobi and Gauss-Seidel. Then after that came a very natural -- ilu is the next one that I want -- the i stands for incomplete, incomplete lu. So what are l and u? That's the letters I always use as a signal for the lower triangular and upper triangular factors of a, the ones that exact elimination would produce. But we don't want to do exact elimination. Why? Because non 0s fill in. Non 0s appear in those spaces where a itself is 0. So the factors are much -- have many more non 0s than the original a. They're not as sparse.

So the idea of incomplete lu is -- and I'll come back to it -- the idea of incomplete lu -- you might say, should have occurred to people earlier -- only keep the non 0s that are non 0s in the original a. Don't allow fill in or allow a certain amount of fill in. You could have a tolerance and you would include a non 0 if it was big enough, but the many, many small non 0s that just fill in and make the elimination long, you throw them out. So p is -- in this case, p is l, is an approximate l times an approximate u. So it's close to a, but because we've refused some of the fill in, it's not exactly a. You have a tolerance there where if you set the tolerance high, then these are exact, but you've got all that fill in. If you set the tolerance to 0, you've got the other [? instance. ?]

So that would give you an idea of preconditions, of reasonable choices. Now I'm going to think about -- first, how do -- how to decide whether the xs approach the correct answer? I need an equation for the error and because my equations are linear, I'm just going to subtract the upper equation from the lower one and I'll call the -- so the error ek will be x minus the true answer. This is x exact. Minus x, my case approximation. So when I subtract that from that, I have ek plus 1 and here I have p minus a. I subtract that from that. I have ek and when I subtract that from that, 0. So very simple error equation. so this is an equation, this would be the error equation. You see the b has disappeared, because I'm only looking at differences now. Let me write that over on this board, even slightly differently.

Now I will multiply through by p inverse. That really shows it just so clearly. So the new error is, if I multiply by p inverse, p inverse times p is i and I have p inverse aek. So what does that say? That says that at every iteration, I multiply the error by that matrix. Of course, if p equals a, if I pick the perfect preconditioner, then p inverse a is the identity. This is the 0 matrix and the error in the first step is 0. I'm solving exactly. I don't plan to do that. I plan to have a p that's close to a, so in some way, this should be close to i. This whole matrix should be close to 0, close in some sense. Let me give a name for that matrix. Let 's call that matrix m. So that's the iteration matrix, which we're never going to display. I mean, it would be madness to display. We don't want to compute p inverse explicitly or display it, but to understand it -- to understand convergence or not covergence, it's m that controls everything.

So this is what I would call pure iteration or maybe another word that's used instead of pure is stationary iteration. What does stationary mean? It means that you do the same thing all the time. At every step, you just use the same equation, where these methods will be smarter. They'll adapt. They'll use more information. They'll converge faster, but this is the simple one -- and this one plays a part in those.

So what's the story on convergence now? What about that matrix m? If I take a starting vector, e 0 -- my initial error could be anything -- I multiply by m to get e 1. Then I multiply m by m again to get e 2. Everu step, I multiply by m. So after k steps, I've multiplied k times by m. I want to know -- so the key question is, does that go to 0? So the question, does it go to 0 and if it does, how fast? The two parts, does it go to 0 and how fast are pretty much controlled by the same property of m. So what's the answer? When does -- if I take -- I don't know anything about e 0, so what property of m is going to decide whether its powers get small, go to 0, so that the error goes to 0? That's convergence. Maybe wall up or maybe stay away from 0. What controls it?

One word is the answer. The eigenvalues. It's the eigenvalues of m, and in particular, it's the biggest. If I have an eigenvalue of m, as far as I know, e-not could be the eigenvector and then every step would multiply by that eigenvalue lambda. So the condition is all eigenvalues of m have to be smaller than 1 and in magnitude of course, because they could be negative, they could be complex. So that's the total condition, that's exactly the requirement on m. How do we say how fast they go to 0? How fast -- if all the eigenvalues are below 1, then the slowest convergence is going to come for the eigenvalue that's the closest to 1. So really it will be -- and this is called the spectral radius and it's usually denoted by a Greek row of m, which is the maximum of these guys. It's the max of the eigenvalues and that has to be below 1. So it's this number, that number, because it's that number which really says what I'm multiplying by it every step. Most likely, e-not, my initial unknown error has some component in the direction of this biggest eigenvalue, in the direction of its eigenvector. Then that component will multiply at every step by lambda, but that one will be the biggest guy, rho. So it's the maximum eigenvalue.

That's really what it comes to. What is the maximum eigenvalue? Let me take Jacobi. Can we figure out the maximum eigenvalue for Jacobi? Of course, if I don't know anything about the matrix, I certainly don't want to compute eigenvalues. I could run Jacobi's method. So again, Jacobi's method, I'm taking p to be the diagonal part, but let me make that explicit and let me take the matrix that we know the best. I'll call it k. So k or a -- this is the a -- is my friend with 2s on the diagonal minus 1s above. I apologize -- I don't really apologize, but I'll pretend to apologize for bringing this matrix in so often. Do we remember its eigenvalues? We computed them last semester, but that's another world. I'll say what they are. The eigenvalues of a -- this is an m now -- the eigenvalues of this matrix a -- I see 2s on the diagonal. So that's a 2. That accounts for the diagonal part, 2 times the identity. I have a minus and then the part that comes from these 1s, which are forward and back.

What would von Neumann say? What would von Neumann do with that matrix? Of course, he would tested on exponentials and he would get 2 minus e to the ik minus e to the minus ik. He would put the two exponentials together into cosines and that's the answer. It's 2 minus 2 cosines of whatever angel's theta. Those are the eigenvalues. The j'th eigenvalue is -- the cosine is at some angle theta j. It's worth since -- maybe just to take again a minute on this matrix. So the eigenvalues are between 0 and 4. The eigenvalues never go negative, right? Because the cosine can't be bigger than 1. Actually, for this matrix the cosine doesn't reach 1. So the smallest eigenvalue is 2 minus 2 times that cosine close to 1. It's too close for comfort, but it is below 1. The eigenvalues are positive and they're between 0 and 4. At the other extreme, theta is near pi. The cosine is near minus 1 and we have something near 4. So the eigenvalues of that -- if you look at that matrix, you should see. Special matrix eigenvalues in this interval 0 to 4. The key point is, how close to 0?

Now what about m -- what about p first? What's p? For Jacobi, so p for Jacobi is the diagonal part, which is exactly 2 i. That's a very simple diagonal part, very simple p. It produces the matrix m. You remember m is the identity matrix minus p inverse a. That's beautiful. So p is just 2 i. So I'm just dividing a by 2, which puts a 1 there and then subtracting from the identity. So it has 0s on the diagonal. That's what we expect. It will be 0s on the diagonal and off the diagonal, what do I have? Minus -- with that minus, p is 2 so it's 1/2. It's 1/2 off the diagonal and otherwise, all 0s. That's the nice matrix that we get for m in one dimension.

I have to say, of course, as you know, we wouldn't be doing any of this for this one dimensional matrix a. It's tridiagonal. We can't ask for more than that. Direct elimination would be top speed, but the 2D case is what we understand that it has a big bandwidth because it's two dimensional, three dimensional even bigger -- and the eigenvalues will come directly from the eigenvalues of this 1D case. So if we understand the 1D case, we'll knock out the two and three dimensional easily. So I'll stay with this one dimensional case, where the matrix m -- and I'll even maybe shouldn't write down exactly what the iteration does -- but if I'm looking now at convergence, do I have or don't I have convergence and how fast? What are the eigenvalues and what's the biggest one? The eigenvalues of a are this and p is just -- p inverse is just 1/2 times the identity. So what are the eigenvalues? Again, I mean, I know the eigenvalues of a, so I divide by 2. Let me write down what I'm going to get now for the eigenvalues of m, remembering this connection. So the eigenvalues of m are 1 minus 1/2 times the eigenvalues of a. Better write that down. 1 minus 1/2 times the eigenvalues of a. Of course, it's very special that we have this terrific example matrix where we know the eigenvalues and really can see what's happening and we know the matrix.

So what happens? I take half of that -- that makes it a 1. When I subtract it from that one, it's gone. 1/2 of that, that 2 is cancelled -- I just get cos theta. Theta sub whatever that angle was. It's a cosine and it's less than 1. So convergence is going to happen and convergence is going to be slow or fast according as this is very near 1 or not -- and unfortunately, it is pretty near 1, because the first, the largest one -- I could tell you what the angles are. That's -- to complete this information, I should have put in what -- and von Neumann got them right. That's j -- is there maybe a pi over n plus 1? I think so. Those are the angles that come in, just the equally spaced angles in the finite case. So we have n eigenvalue. The biggest one is going to be when j is 1. So the biggest one, rho, the maximum of these guys, of these lambdas will be the cosine of when j is 1 pi over n plus 1. That's the one that's slowing us down. Notice the eigenvector that's slowing us down, because it's the eigenvector that goes with that eigenvalue that's the biggest eigenvalue that's going to hang around the longest, decay the slowest. That eigenvector is very smooth. It's the low frequency. It's frequency 1. So this is j equal to 1. This is the low frequency, j equal 1, low frequency. The eigenvector is just samples of the sine, of a discrete sine. One lowest frequency that the grid can sustain.

Why do I emphasize which eigenvector it is? Because when you know which eigenvector is slowing you down -- and here it's the eigenvector you would think would be the simplest one to handle -- this is the easiest eigenvector, higher frequency eigenvector, have eigenvalues that start dropping. Cosine of 2 pi. Cosine 3 pi over n plus 1 is going to be further away from 1. Now I have to admit that the very highest frequency, when j is capital n, so when j is capital n, I have to admit it happens to come back with a minus sign, but the magnitude is still just as bad. The cosine of n pi over n plus 1 is just the negative of that one. So actually we have two eigenvalues that are both hitting the spectral radius and both of them are the worst, the slowest converging eigenvectors. This eigenvector looks very smooth. The eigenvector for that high frequency one is the fastest oscillating. If I graph the eigenvectors, eigenvalues, I will. I'll graph the eigenvalues. You'll see that it's the first and the last that are nearest in magnitude to 1.

I was just going to say, can I anticipate multigrid for a minute? So multigrid does -- is going to try to get these guys to go faster. Now this one can be handled, you'll see, by just a simple adjustment of Jacobi. I just wait Jacobi a little bit. Instead of 2 i, I take 3 i or 4 i. 3 i would be very good for the preconditioner. That will have the effect on this high frequency one that'll be well down now, below 1. This will still be the slowpoke in converging. That's where multigrid says, that's low frequency. That's too low. That's not converging quickly. Take a different grid on which that same thing looks higher frequency. We'll see that tomorrow.

So if I stay with Jacobi, what have I learned? I've learned that this is rho. For this matrix, the spectral radius is that number and let's just see how big it is. How big is -- that's a cosine of a small number. So the cosine of a small number, theta near -- theta near 0, the cosine of theta is 1 minus -- so it's below 1, of course -- 1/2 theta squared. Theta is pi over n plus 1 squared. All these simple estimates really tell you what happens when you do Jacobi iteration. You'll see it on the output from -- I mean, you can -- I hope will -- code up Jacobi's method, taking p to be a multiple of the identity, so very fast and you'll see the convergence rate you can graph.

Here's an interesting point and I'll make it again. What should you graph in seeing -- to display convergence and rate of convergence? Let me draw a possible graph here. Let me draw two possible graphs. I could graph -- this is k as I take more steps. I could graph the error, the size of the error. It's a vector, so I take its length. I don't know where I start -- at e-not. This is k equals 0. This is the initial guess. If I use -- I'm going to draw a graph, which you'll draw much better by actually getting MATLAB to do it. Before I draw it, let me mention the other graph we could draw. This tells us the error in -- this is the error in the solution. If I think that's the natural thing to draw, it is, but we could also draw the error in the equation, the residual error. So I'll just put up here what I mean by that. If I put it in -- so ax equals b exactly. axk is close. The residual is -- let me call it r, as everybody does -- is the difference ax minus axk. It's how close we are to b, so it's b minus axk. The true ax gets b exactly right. The wrong ax gets b nearly right. So you see, because we have this linear problem, it's a times the error. x minus xk is ek. So this is the residual and I could grapj that. That's how close my answer comes to getting b right, how close the equation is. Where this is how close the xk is. So that would be the other thing I could graph. rk is the size of the residual, which is aek.

So now I have to remember what I think these graphs look like, but of course, it's a little absurd for me to draw them when -- I think what we see is that the residual drops pretty fast. The residual drops pretty fast and keeps going. Let's say Jacobi's method or in a minute or maybe the next time, another iterative method -- but the error in the solution starts down fast and what's happening here is the high frequency components are getting killed, but those low frequency components are not disappearing, because we figured out that they give the spectral radius the largest eigenvalue. So somehow it slopes off and it's -- the price of such a simple method is that you don't get great convergence.

Of course, if I put this on a log-log plot, the slope would be exactly connected to the spectral radius. That's the rate of decay. It's the -- ek is approximately -- the size of ek is approximately that worst eigenvalue to the kth power. Do you see how this can happen? You see, this can be quite small, the residual quite small. I think it's very important to see the two different quantities. The residual can be quite small, but then from the residual, if I wanted to get back to the error, I'd have to multiply by a inverse and a inverse is not small. Our matrix a is pretty -- has eigenvalues close to 0. So a inverse is pretty big. a inverse times aek -- this is small but a inverse picks up those smallest eigenvalues of a -- picks up those low frequencies and it's bigger, slower to go to 0 as we see in this picture. So it's this picture that is our problem. That shows where Jacobi is not good enough.

So what more to say about Jacobi? For this matrix k, we know its eigenvalues. The eigenvalues of m, we know exactly. The matrix m we know exactly. You see in that matrix -- now if I said, look at that matrix and tell me something about its eigenvalues, what would you say? You would say, it's symmetric, so the eigenvalues are real and you would say that every eigenvalue is smaller than 1. Why? Because every eigenvalue is -- if you like to remember this business of silly circles, every eigenvalue is in a circle centered at this number - in other words, centered at 0 and its radius is the sum of those, which is 1. So we know that every eigenvalue is in the unit circle and actually, it comes inside the unit circle, but not very much -- only by 1 over n squared. What do I learn from this 1 over n squared? How many iterations do I have to take to reduce the error by, let's say, 10 or e or whatever? If the eigenvalues are 1 minus a constant over n squared, I'll have to take the n squared power. So if I have a matrix of order 100, the number of iterations I'm going to need is 100 squared, 10,000 steps. Every step is fast, especially with the Jacobi choice. Here's my iteration -- and p is just a multiple of the identity here. So every step is certainly fast. Every step involves a matrix multiplication and that's why I began the lecture by saying, matrix multiplications are fast. a times xk -- that's the work and it's not much. It's maybe five non 0s and a -- so five n. That's fast, but if I have to do it 10,000 times just to reduce the error by some constant factor that's not good. So that's why people think Jacobi gives us the idea. It does knock off the high frequencies quite quickly, but it leaves those low frequencies.

I was going to mention weighed Jacobi. Weighted Jacobi is to take instead of p -- and put in a weight. So weighted Jacobi -- let me put this here and we'll come back to it. So weighted Jacobi simply takes p to be the diagonal over a factor omega. For example, omega equals -- 2/3 is a good choice. Let me draw the -- you can imagine that if I -- this is 2 times the identity, so it's just 2 over omega i -- I'm just choosing a different constant. I'm just choosing a different constant in p and has a simple affect on n. I'll draw the pictures. The eigenvalues of m are these cosines. So I'm going to draw those cosines, the eigenvalues of m. So they start -- I'll just draw the picture of cos theta if I can. What does a graph of cos theta look like? From 0 to pi, so I'm going to see -- this is theta equals 0. This is theta equal pi. If I just draw cos theta -- please tell me if this isn't it. Is it something like that? Now where are these theta 1. the bad one. Theta 2, theta 3, theta 4, theta five, equally bad down here.

So there's cos, there's the biggest one. That's rho and you see it's pretty near 1. Here I'm going to hit -- at this frequency, I'm going to hit minus rho. It'll be, again, near 1. So the spectral radius is that, that's rho. The cosine of that smallest angle, theta over n plus 1. What happens when I bring in these weights? It turns out that now the eigenvalues of m with the weights will be 1 -- it turns out 1 minus omega plus omega cos theta. You'll see it in the notes, no problem. So it's simply influenced by omega. If omega is 2/3 -- if omega is 2/3, then I have 1 minus omega as 1/3 plus 2/3 cos theta. Instead of graphing cos theta, which I've done -- let me draw -- let me complete this graph for cos theta. The eigenvalues unweighted with weight omega equal 1. Now it's this thing I draw. So what happens? What's going on here? For a small theta, cosine is near 1. I'm again near 1. In fact, I'm probably even little worse, probably up here. I think it goes down like that. So it starts at near 1, but it ends when theta is pi. This is 1/3 minus 2/3. It ends at minus 1/3. That's a lot smarter. This one ended at minus 1, that one ended at near minus 1 but this one will end with omega -- this is the omega equal 2/3 of weighted 1. All you've done is fix the high frequency. When I say, that was a good thing to do, the high frequencies are no longer a problem. The low frequencies still are. So the only way we'll get out of that problem is moving to multigrid.

Can I, in the remaining minute, report on Gauss-Seidel and then I'll come back to it? But quick report on Gauss-Seidel. So what's the difference in Gauss-Seidel? What does Gauss-Seidel do? It uses the latest information. As soon as you compute an x, you use it. Jacobi compueted all -- had to key all these xs until it found all the new xs, but the Gauss-Seidel idea is, as soon as you have a better x, put it in the system. So the storage is cut in half, because you're not saving a big old vector while you compute the new one. I'll write down Gauss-Seidel again. So it's more up to date and the effect is that the Jacobi eigenvalues get squared. So you're squaring numbers that are less than 1. So Gauss-Seidel, this gives the Jacobi eigenvalues squared. The result is that if I draw the same curve for Gauss-Seidel, I'll be below, right? Essentially a Gauss-Seidel step is worth two Jacobi steps. Eigenvalues get squared as they would do in two Jacobi steps. If I do Jacobi twice, I square its eigenvalues. So it seems like absolutely a smart move. It's not brilliant, because this becomes cosine squared. I lose the 1/2, but I'm still very, very close to 1. It still takes order n squared iterations to get anywhere. so Jacobi is, you could say, replaced by Gauss-Seidel as being one step better, but it's not good enough and that's why the subject kept going. This gives rho as 1 minus order of 1 over n first power, where these were 1 minus 1 over n squared. So this is definitely further from 1 and then this can be way further from 1.

So you can see some of the numerical experiments and analysis that's coming. I'll see you Wednesday to finish this and move into multigrid, which uses these guys.