Lecture 16: General Methods for Sparse Systems | Video Lectures | Mathematical Methods for Engineers II | Mathematics

Download this transcript - PDF (English - US)

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: This morning's Siam News. Well, it's an article on something that's not totally up to date because it's a method that was developed 50 years ago called alternating direction and it was exactly developed-- first you did an iteration in the x directions. Of course, we're tridiagonal along every row, we're one-dimensional. And then the second half-step is an iteration in the y direction. Again, tridiagonal, so that iteration's very fast. You flip-- alternating direction describes it perfectly. And that method is still used. It got developed in the oil reservoir simulation world. That is, Houston. And the article says the evolution of techniques for oil reservoir simulation has continued. And then mentions 3-- so it went from this alternating direction to line over relaxation, LSOR over relaxation and I only spoke briefly about that. What's that word line mean? That means that we're taking hole rows at once. Working with blocks and I'll add some notes about that today. So that was the next step, line SOR. Then to the Newton-Krylov schemes, that's what's coming. Conjugate [? gradients ?] with ILU type preconditions. So again, we're talking about preconditioners and one important choice is this incomplete LU. So notes on Gauss-Seidel and I should maybe add incomplete LU.

By the way, experiments badly needed here too. Because incomplete LU has this tolerance and as you move the tolerance you get closer or far away to the exact L and U, but you get faster or slower. So what's the balance? What's the right tolerance? How do these methods compare-- and overall, how do these iterative methods that we're speaking about compare with these direct methods? So we really are in this new chapter solving large systems, facing a whole lot of possible new experiments. So just to continue the history or the future as described in Siam News, so there was alternating direction, line SOR. Then to Krylov schemes with ILU preconditioners in use today. So this is a big use area and the multiscale, multigrid solvers that are quote, now starting to come of age. So that's a third topic for today's lecture, multigrid. That of course, it's been around for some time, but it involved, especially at the start more effort in programming. So it didn't quickly-- I mean, it gave great results on model problems. On other problems, too, but still it kind of hasn't yet, but is penetrating into production code or the oil reservoir industry or other big industries would use it. So here you see really a pretty fair picture of where we are today with large solvers. These are very important.

Minimum degree-- if the size is not overwhelming direct methods are going to be right up there. And these are the pure iterations, stationary iterations that with this idea, this incomplete LU idea get pretty good. But the formula's always simple there. The new x comes from the old x by the same iterations. Whereas, if we use the ideas there as preconditioners and look to see-- so the point was, and I'll repeat it again, that Gauss-Seidel and Jacobi, the standard iterations were quite satisfactory on high frequencies. They damped those out quite well in the error, but they don't tackle the low frequency, the smooth part of the error, very efficiently. And that's what multigrid does incredibly well, so it's a nice mixture of the two. These are the smoothers, this is the change of scale, change of grid that's coming now and will come. So all this will be my subject for this week and next, prior to spring break and I think you will have enough to go on to do experiments. It could be highly interesting. OK, so that's the Siam News report.

OK, what do I want to say about minimum degree? And more will go on the web. Just a few more notes about what we did on Monday after class. If you take a matrix like that, you take a matrix and how is it described in sparse format? Sparse format you can see-- normally you don't have to see it, it's usually inside the sparse math lab, but if you want to see what's happening you could ask for i, j, and s. Now what are those? What will that produce? i will be the list of nonzeroes, the list of the row numbers of the nonzeroes. So i will be-- maybe I'll write it as a row so I'm transposing, So the row numbers would be row 1, row 2, not row 3. Row, 1, 2, 3, row 2, 3. The column numbers would be-- these are going to be the pairs i,j. So that's column 1 twice, column 2 three times. Column 3 twice. And what's s? s is the actual numbers in position i,j is the number 2. In position 2,1 is the number minus 1. So those numbers minus 1, 2, minus 1 is that column. And then minus 1, 2 is that column. So of course, we've got all the information in the matrix. We know every nonzero position and we know what that entry is. And of course, our matrices are much bigger than this one, but already we can see one useful point. Point being that key word that the column numbers-- this I mean, imagine we have a large matrix, so this j is not very efficiently recorded here because what's j? We're looking at nonzeroes a column at a time, so of course, we'll have a few ones for the nonzeroes in column 1 and then some twos and then some threes and some fours, but all this row is-- the real information in that row is a pointer to-- in other words, I don't have to repeat 2 three times.

So really if I put inside here a compressed version of j, a pointer would be a short factor that just has a 1. It says look in position 1 for the start of column 1. Then a 3, see I didn't need that because that just continued column 1. That just repeated the j I already had, but this tells me that in the third position I start down column 2. And this would tell me that in position number 4, 5, 6, I guess, I start down column 3. You see the point of that 6. That 6 picks out-- yes-- 1, 2, 3, 4, 5, the sixth row number is the beginning of information on column 3. And then it's conventional to have a pointer, 8, to say finish. So that's pointing to empty space. So in other words, j got compressed to 1, 3, 6, 8 and in a large problem it would get seriously compressed. OK, so that's the form in which the code keeps the matrix and does the reoderings. Let me just mention that the opposite of this would create the matrix out of these i, j, s. What command would that be? You could create the matrix a out of the command sparse, would be good. Sparse of the inputs now would be i, j, s. If I input i, j, and s then that lab creates a matrix. And sometimes I may want to include the m and n, the shape of the matrix as further parameters, but actually here I wouldn't have to. So this is the opposite command from this one that we're interested in. You can imagine.

Suppose by using minimum degree or some other decision, I've eliminated up to a certain point and I want to put the remaining columns in a different order. I can do that just by playing with the pointers. So it's a very efficient structure for sparse matrices. And then comes the question, OK, what's the good order? I was surprised to realize how open a problem that still is, even from the expert who's developing the key code. So approximate minimum degree of course by that word you see that it allows freedom. There are also decisions to be made when minimum degree is a tie between several nodes as it commonly is. And it's rather nice to get a math lab movie that shows the order in which nodes are eliminated, edges are removed from the graph, from the mesh. You'll see that. All right, I'm ready to go ahead now to some comment on Gauss-Seidel. Maybe I'll put that here. In fact, maybe I'll just take the same matrix. So now I want to remember, what is the Gauss-Seidel method for ax equal b. What's the iteration now? And again, a will be the same: 2, minus 1, zero minus 1, 2, minus 1. My favorite, OK.

So all these iterations you remember are splitting of the matrix. Some of the matrix goes on the left-hand side. That's the preconditioner. So it's p, the preconditioner that's sort of close to a in some sense, but easy to work with, multiplies the new thing and on the other side is p minus a-- the rest of the matrix-- that multiplies the old plus right-hand side b. Let's just remember again, that if we converge to the point where this is the same as this then we have the right answer. We have the right answer because if that's the same as that, px is the same as px, this comes over on the other side, ax equals b. So when it converges it converges to the right answer, but you remember that the key equation was that the new error is the old error multiplied by this matrix i minus p inverse a. That's the iteration matrix you could say. I just get it when I multiply both sides by p inverse. And so the problem is choose p so that this is easy to do and at the same time, this has small Eigen values or as small as you can get.

So the Gauss-Seidel is a particular choice, which takes the lower-- p is the lower triangular part here. So p is, so I just thought I'd better write down explicitly Gauss-Seidel takes that with zeroes there as multiplying-- that's p. And what's p minus a? It's the rest and it's been moved to the other side so the rest is going to be the strictly upper triangular part. And because it's moved over to the other side it'll have a plus sign. It'll be 1 and 1. x sub k plus b. I just write it out so that you see, totally explicitly how being triangular makes the solution step immediate. Because the first equation will tell us the first component, right? The first equation, because it's triangular there's only an entry here. And by the way, all these codes, including Tim Davis' minimum degree codes, their very first step is reorder the equations if necessary to guarantee that the main diagonal has no zeroes. We want to know that in advance. So that you just do, let's assume of course, for us it happens automatically. OK, so that's not zero. And that first equation gives us the first component then we know the first component so we use it here and the second equation tells us what the second component is. We use that second component in the third equation to find the third component. So there's no loss of speed compared to diagonal and actually, it's faster because the storage we're changing-- we're using the new first component to find the second and the new second component to find the third, so we can overwrite x sub k by x sub k plus 1. Let's see. I guess if I had space and had prepared I would figure out what this matrix is. Maybe you could do that. Figure out what-- here's p. It's invertable. P inverse is going to be quite simple. Actually, it would be extremely simple. You'll be able to see its Eigen values immediately. We could stop to do it, but I think if you do it it's more valuable. So we would find that its Eigen values were below 1 and that they were the squares of the Jacobi Eigen values, so the method is twice as fast as Jacobi. But it has the same good feature of damping the high frequencies and the same bad feature of not very much damping the low, smooth part of the error.

Now we're ready for multigrid, which is intended to solve it. So multigrid is going to be-- the point is by changing grid, by changing from the fine grid to a coarse grid the smooth oscillation doesn't look so smooth. On the coarse grid, it's effective frequency is effectively doubled and it's moving over in the direction where the smoother can attack it. We'll see that happen.

So now, really this begins the lecture on multigrid, which is going to take certainly, today and Friday. There's a very good book by Briggs and others. He wrote the first edition of the book and coauthored the second one. And it's beautifully short and clear and simple and this presentation will follow the same exposition that the book does. So the key idea is to see what-- how to go from a fine grid problem, so this is on the fine grid with step size a. So this A sub h is our problem. A sub h is A. A sub h, b sub h, we're looking for u sub h on the fine grid. Fine grid means lots of mesh points, lots of unknowns, big problem. OK, so the idea is going to be somehow to-- let me just start with two grids. Well, you see what you do on the fine grid at the start. You iterate. You use Gauss-Seidel or Jacobi, whatever. Maybe three times. Maybe three iterations, but don't continue to a thousand iterations, it's a waste of time. Then here is the real-- this is the multigrid part. The multigrid part is going-- you get an answer after a few iterations, after you've smoothed it. You compute this residual, which is the amount you're wrong in the equation, the difference between the right side and the left side. I can put h is here too, to emphasize. This is residuals being computed with what we have, which is on the fine grid.

Now here's the change of grid. Number two says, restrict to a coarse grid. Take that rh, that residual which is probably a bit smooth and move it to the coarse grid. So the coarse grid we identify by the fact that it's mesh width is twice as big, 2h. And so we're going to need some restriction, so the input to this multigrid is to decide on a restriction matrix, how shall we take values that are given on a fine grid and restrict them to a coarse grid. Of course, one restriction is-- and it's not that bad-- is just to take the values of r that are on the coarse grid at every other grid point, use those. Or we could take into account the neighbors. You see the question. So this is coming now in this restriction. So the restriction is, restriction r, this is fine to course, and we have this question of use neighbors or not? That will answer it. Well, the decision on what the matrix r is. That's the key.

OK, so we choose an r, we restrict, and we have now a half-size problem. Well, half-size in 1D, quarter-size in 2D, eighth-size in 3D because the mesh width has doubled in every coordinate. Now I've put solve. I better put quotes around that. Well, solve the problem on the coarse mesh. So the idea is that that's much faster Well, you might say, what problem? Well, we have to create somehow the matrix on the coarse mesh. That's not necessarily given to us. We start the problem with A, I mean, A sub h, the matrix on the fine. Can I just mention that this is probably, this v-cycle-- that letter v is supposed to suggest going down to the coarse mesh and back up to the fine mesh. And the standard notation in all the multigrid books is a capital V. And why do I use a small v? This is my educational contribution. Capital V is appropriate when you have several meshes. You go to 2h, you go to 4h, you go to 8h, back to 4h, back to 2h, back to h. So a deeper multigrid. So it would look more complicated. I would repeat this idea instead of solve here. Instead of solve at the 2h level, which of course, I'm not going to do exactly anyway. But if I was in a big V-cycle I would iterate a few times to smooth here, just as I did here. I would iterate using weighted Jacobi or something and then go down to the 4h mesh. So restrict to the 4h mesh. So you can see that I would stay in this little loop to get to the bottom, 8h and then I would start back up. So I think it's a good idea to use a small leter v to tell us rather than saying two grid capital V-cycle. I'm just going to use a small v to signal right away that it's two grids. Fine, coarse, fine.

OK, so we have some solution, not necessarily exact and actually, there isn't that much point in getting it exact because what we want is to move toward the right answer. And notice that I'm looking here at the error rather than looking at A 2h, u 2h-- the actual solution. What I'm computing with is the correction term. So this is the correction term that I find. So this is the residual over here. this is the A 2h that I still have to choose and then solving that will give me a correction, but that correction is on coarse mesh. I've only had that correction E 2h defined on the 2h mesh. So now comes the other part of multigrid. Get back to the fine mesh, climb back up. I take that E 2h and apply an interpolation. I need to create an interpolation. This is going to be coarse to fine. And I have to decide how to do that. Once I've decided it that gives me a correction at the fine mesh level, which is where I'm really working, so I make the correction and I have a better answer. And probably, I iterate a few times-- that would be called a post smoother on that correction. And maybe that's where I stop or maybe I repeat little v-cycles, but that's not brilliant. If you're going to do multiple, if you're going to repeat multigrid it's much more efficient to move to more and more coarse meshes because the calculations on these are way faster than on the fine mesh.

You see, rather than-- I much prefer to go way down and back up then a lot of v-cycles. I mean, that would be a small w-cycle. And not brilliant. Much better to use a big V-cycle or a big W-cycle. There's a place for W-cycles but they're capital W-cycles because you want to get down where it is very fast, very inexpensive and in fact, so efficient that multigrid achieves this holy grail of giving you an answer with whatever accuracy you want-- giving you an answer in o of n operations. n being the matrix size. Captial N squared in our example. So that's the fantastic result from multigrid. So I'll come back to this, but just say if the size of A is N, which is capital N in 1D, captial N squared, capital N cubed, then the multigrid works in o of in flops, Floating Point Operation. Not even log n, which we think of for the-- n log n we think of for the FFP, for the Fast Fourier Transform. But here's it's actually order of n. That really is a goal worth achieving.

So we know what the smoothers might be. It's the i and the r that are new. And then the analysis of convergence, why does it converge? But to use the method, and so by the way, projects using multigrid are totally welcome also. So let me take the interpolation i. So what's that? Coarse to fine. Let me imagine I'm in 1D, so there's the x direction. Here's one end of the interval, here's the other end of the interval and suppose my boundary conditions are zero. So this is the coarse mesh and let me put in here, make it the fine mesh. So at this point I have an answer on the coarse mesh, which satisfies the boundary condition and let's say it's there, there, and there. OK, and I need values on the fine mesh at the other points. Of course, I'm going to use these values when I interpolate. I mean, that word interpolation means implies, keep what you have on the given points and I'm just going to do linear interpolation. So halfway there, halfway there, and halfway there. So ignoring the boundary conditions, the zeroes, it's going to be 5 by 2 I think, this matrix i in this example will be 5 by 2 because it has 2 inputs, these and it has 5 outputs. And what's the matrix? Well, for their second and fourth outputs it uses exactly-- this multiplies x or I should say u, u sub 1 and u sub 2. Follow the notes again. The notes use v sub 1 and v sub 2 for the coarse mesh values to avoid h 2h, extra subscripts. So this height is v sub 1. This height is v sub 2 and I'm going to take-- this'll be the new u sub 1, u sub 3, and u sub 5 And u sub 3 I'm going to save. I better put i and make this into a true equation, i sub v. This is going to give me the values u sub 1, u sub 2, u sub 3, u sub 4 and u sub 5. So u sub 2 is just v sub 1. u sub 4 is just v sub 2. No problem. And what do I do with the other ones?

Let's see. I guess u sub 1 is halfway between v1 and 0. OK, so it's 1/2. u sub 1 is just 1/2 v1 It doesn't involve v2. u3, which is sitting here is 1/2 of these so it's 1/2 of this and 1/2 of that where this was nothing. Then u sub 4 was the same as v sub 2, just take v sub 2. And what's u sub 5? u sub 5 is 1/2 of v sub 2 and 1/2 of zero. So it's a simple matrix. The interpolation matrix. It's just reasonable to use the letter I and we all recognize here it's a rectangular matrix, not the identity. I doesn't stand for identity here, it stands for interpolation. OK, so I hope you can focus on the matrix. And if it was much, much bigger, but still in 1D it wouldn't look very different. Each column would have 1/2, 1, 1/2. That's a typical part of the interpolation matrix. It uses the value v sub 1, saves it at the center point, uses 1/2 of it at the previous point, and 1/2 of it at the right-hand point. And similarly, you see it's sort of a typical rectangular matrix in signal processing as well. I have to say something about what happens in 2D because our problems are in 2D. Maybe I just draw a typical 2D mesh. These are the coarse values, this is the coarse mesh and now of course, we'll save those values. I will have a 1, so it'll save those values when it goes to the fine mesh, the fine mesh being 1/2. So I just have to decide what to do there. So question, what value shall I take? Oh, I suppose I have to decide there, too. These are all new points and what shall I do. I'm just going to stay with linear interpolation, which is quite fast and satisfactory. This value will be the average of those two. This value will be the average of those and what will this one be in the center? It's the average of these four, of course. And that's what we'll get. Actually, that Kronecker product business, a Kronecker product of a 1D matrix like that with itself would be a 2D interpolation that would do exactly this. We can write it out more fully, but it's not fun to write out really fully because if this is 5 by 2 then I 2D would be 25 by 4. So we're right away, even on this tiny, tiny problem, we're right away up to 25 internal mesh points on the fine mesh. That's the matrix that I want you to see and the natural choice for R is to transpose.

So one possible choice-- and it will, the answer will be yes, we do use [? Naver. ?] So R will be the transpose of I. Actually, I'll need a factor here. I think it's 1/2. Let me just see what it is. Let me see why we need a little multiple there to make things right. So there's my matrix I, so I'm just going to transpose it and let me factor 1/2 out of that so you see these matrices without fractions inside. That would be 1, 2, 1. So I'm going to transpose that. 1, 2, 1, zero, zero. Zero, zero, 1, 2, 1. So I factored 1/2 out of I when I did that. And then I believe that I need another 1/2 here. I think I need 1/4. And what do I mean by needing 1/4? Of course, the main point is that this is the transpose of this. And you'll see that that preserves the symmetry that's just sort of the natural thing to pick. Not the only possibility, but it's very good to connect the restriction with the interpolation.

Of course, you can't take the restriction to be the inverse of the interpolation. Why not? I mean that sort of mental idea is that the restriction goes in one direction and the interpolation in the other direction. So why not just let one be the inverse of the other? Well of course, that matrix doesn't have an inverse. I mean it's rank is only 2. When we do a restriction we're going to lose information, the interpolation can't put it all back because we threw it away when we restrict it to the coarse mesh.

Now why 1/4? I think it's just if we had all constants, like all ones as our mesh values, we would want to get all ones for the restriction. So if I multiply this matrix by the vector of all ones-- so assuming you remember where R comes in by just copying this equation, you remember that when I go from the h mesh, it's fine to coarse. So R times ones on the fine mesh, all ones on the fine mesh gives ones on the coarse mesh. That's why I wanted that-- you see, if I multiply by all ones then that multiplication gives me 1 plus 2 plus 1 equal 4 and I need to divide by 4 to get back to the 1. You quickly discover that that just arises because of a change of scale. It's just a factor of 1/2. In 2D that 1/2 there would change to 1/4. In 3D it'd be 1/8. It's just the scaling to preserve constants.

And now the remaining question is, what's A 2h. That's the other matrix here that we're not given. We're given the A sub h, our problem. We choose an R and we choose an I and if we're smart, ones that transpose of the other. And then the question is what's A 2h? And there's a beautiful answer, so let me just say what it is. A 2h, a coarse mesh matrix is the fine mesh matrix, but I need to do first-- if this is going to apply to the coarse mesh guys, the v's I have to do an interpolation. You'll see it here. This is 5 by 5 in my problem. So I need to first do the interpolation to get from a vector blank 2 up to a vector blank 5. Then this and then the restriction. That's terrific. That's A 2h. So this is A 2h and we'd better do an example. I guess it's going to be next time. So that'll be the first thing for next time, would be to see OK, for the standard second difference A sub h, for these piecewise linear interpolation and restriction, what comes out as A 2h? Do we get the standard second difference on the coarse mesh, which you really hope we do and we do. This will be the same difference in our example on the coarse mesh. So that this middle step of multigrid is exactly what we would have expected. What we woud've had if we set up the problem originally on the coarse mesh. But now we're going to get an answer that we take back to the fine mesh.

OK, so I'll do that next time. Some homeworks here that I held for an extra two days and now I'm returning and some thoughts to go up on to the website about possible projects, but actually, I expressed a lot of them at the beginning of the lecture. OK, thanks. Good.