Lecture 20: Fast Poisson Solver | 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: I thought I would, in this last lecture before the break, speak about one specific topic. It's often referred to as a Fast Poisson Solver, so what does Poisson mean? Poisson means Laplace's equation. So, this is the five-point Laplacian, which could be some other discrete Laplace matrix, but let's take the one we know. So we're in two dimensions, and you use the word Poission's name, instead of Laplace's name, when there's a non-zero right-hand side. So otherwise, it's a Laplace Solver, but here Poisson. OK. So, remember that right-hand side comes from maybe a right-hand side f of xy in the differential equations, but it also comes from non-zero boundary conditions, because non-zero boundary conditions, when the five-points hit a boundary, that known value was moved over to the right-hand side and becomes part of f. So f comes from the non-zero boundary conditions, as well as the non-zero right-hand side. OK.

Important problems, but special on a square or a rectangle or a cube or a box, so we're speaking about special geometry. Today as we've been doing, I'll take the case on a square, and you will see that the whole idea would not fly, if we were on an ellipse or something like that. But a lot of problems on rectangular domains do appear in applications, or we could use, I hope you always think now about possible preconditioners. Any time you have something fast, it's a candidate to be a preconditioner for a real problem. The real problem might not be on a square, or the real problem might not have the constant co-efficient that we have in the Laplacian. In that case, you're not solving the exact problem, but one that could be reasonably close to it. OK. So we've discussed the solution to this problem, and actually I have a little more to say about the movie that is now on the website, because I'm quite excited about that movie. I'll come to the movie in a second.

Let me just say what today's lecture would be. The key idea here is that the eigenvalues and eigenvectors of this giant matrix of order n squared by n squared, sizes n squared by n squared. The eigenvalues and the eigenvectors are known. First of all, they're known. The eigenvectors are nice discrete sine functions; they're sines, because I'm assuming they come to zero at the boundary; so that's why I have a sine, rather than a cosine. First they're known, and second we can work with them very quickly, using the FFT. So the point is that it's quite exceptional; in fact, I don't think I know any comparable example, in which a linear system is solved by using the eigenvalues and eigenvectors. Usually eigenvalues and eigenvectors, they pay off for differential equations that are growing in time. Then it is worth computing them, because then you can just multiply by either the lamda t, and you know what's happening later. Eigenvectors, eigenvalues have other purposes, but very, very rarely are they used to solve a linear system. I mean, it's usually far more work. And of course, it would be incredibly more work if we had to find the eigenvalues and eigenvectors, but for this problem we know them.

And it also would be incredibly more work if the matrix of eigenvectors, the basis matrix, the key matrix that I'll denote by s -- so the eigenvectors go in a matrix s; the eigenvalues go in a matrix capital lambda; so we know lambda. It's going to be a simple matrix, just diagonal, got those n numbers and n squared numbers, I guess, because we're of size n squared. But these eigenvectors we know. So we know them, and we can compute quickly with them, using the FFT, or using, you might want me to say, fast sine transform, discrete sine transform, rather than Fourier transform, which we think of as doing the complex exponential. So this is a small special fast Fourier world. It's a special fast Fourier world, in which the FFT and the related sine and cosine give us a quick answer, faster than elimination, because you all know, it's n log n, if n is the number of Fourier components, and that's hard to beat.

Then I'll mention also, for this same problem, there is another way in which it can be simplified. It's not a Fourier way, just a direct combining neighboring equations, so that'll be number two. OK, but mostly the lecture is about number one. The best example I know in which you would use eigenvectors/eigenvalues to solve an ordinary linear system, and I'll say in a word just how you do it, and then what these eigenvectors are. OK.

Pause. Time out to say something about the movie. So that movie that's now on the website does sparse elimination, and the example it takes is k2d, and you can make it 8 squared by 8 squared, 10 squared by 10 squared, 20 squared by 20 squared, because it shows the order that the nodes are eliminated and the graphs of non-zero's in the matrix. It's a bit slow. If you do 20 by 20, you have to go away for lunch, well maybe not lunch, but at least coffee before it's done, but when it's done, it counts the number of -- it shows the non-zero's that are in the elimination in the factor l from elimination, and it counts the number of non-zero's, and I was just talking to Persson about also getting it to count the number of actual [? flops ?].

Well, why am I interested? I'm interested, because I don't know what power of n those numbers are going there. I don't know whether the number of non-zero's, and I did 10, 20, 30, and it looked not too far from the power capital n cubed, but the notes say for nested dissection, which would be another ordering n squared log n. So, and of course, what power we get depends on what algorithm we use for ordering, so nested dissection is one ordering, which hopefully we could put into another movie. The movie now has exact minimum degree, real MMD, which takes as the next node the one with absolutely the minimum degree, not just close to it. Anyway have a look at that movie, and if you have any interest, see how the count increases as n increases, and also you could change, slightly adapt, the algorithm that creates the ordering, creates the permutation and see what happens there. There's a lot to do with that, and it's a pretty fundamental problem. Actually, we're talking there about what MATLAB's backslash operation will do for this equation. So MATLAB'S backslash operation will use elimination; it won't use fast Poisson solvers, but now let me come to the fast Poisson solver. OK.

So I guess the main point is I have to say what are the eigenvalues and eigenvectors, and how do they get used. So let me say, how do they get used. So how can I use eigenvalues and eigenvectors to solve a problem like that. Let me just call it k rather than k2d, sorry. ku equals f, OK, by eigenvectors. OK, there are three steps. Step one. Expand the right-hand side as a combination of the eigenvectors. OK. So expand f, this right-hand side vector, as some combination of the first eigenvector, maybe I'm going to call it y1. Second eigenvector, nth eigenvector, OK, good. That means what do I mean, I mean you have to find those c's, that's the job there. Find the co-efficients. So that's a job, a numerical -- it's a linear system to solve and we'll see what it amounts to. OK, but suppose the right-hand side is a combination of the eigenvectors, how can you use that? Well, step two is the real quick step. Divide each ci by the eigenvector, eigenvalue, lambda. OK, so it's my eigenvector, so I'm assuming that kyi is lambda iyi and that we know these guys. So this is known. And now the question, I'm just saying, how do we assume known?

So my question now is how do we use it? OK, step one, the idea is going to be write everything in terms of eigenvectors. Work with the eigenvectors, because if I've got eigenvectors, the step is scaler; I just divide these numbers by those numbers, and then I've got the answer. And then construct, the correct answer will be u, will be c1 over lambda 1y1 and c2 over lambda 2y2 up to cn over lambda nyn, a combination of those same eigenvectors with the same coefficients, just divided by lambda. But this is, of course, another numerical job; this is like adding up a Fourier series; this is like finding the Fourier coefficients, this is like reconstructing the input, only because I've divided by lamba i, I'm getting the output here is u, when the input was x. And do you see that that is the correct answer? All I have to do is check that ku equals f, so check that this answer from step three is the right answer. OK, so I multiplied by k. When I multiply y1 by k, a factor lambda 1 appears, the eigenvalue; it cancels that; well that's y divided, so it would cancel, and I have c1y1. When I multiply this by k, ky2 is lambda 2y2; cancel the lambda 2s, and you're left with c2y2, and so on. So, when I multiplied by k, I got f. That's it.

So that's the whole idea, written out, but now what actual computations go into steps one and three? Step two is pretty simple. Well, actually this is a good way to look it. I want to write that same algorithm in matrix language. OK, so in matrix form we have the matrix of eigenvectors, and that's what I'm calling s. And it's got the eigenvectors y1 y2, yn, and it's columns, and the eigenvalue matrix, we need a name for that too, and that we decided to call lambda, and that's got the eigenvalues on its diagonal. So this is 18.06 linear algebra. The matrix of eigenvectors if I multiply k by s, then I'm multiplying k by all its little eigenvectors, and k times this gives me lambda 1y1, k times y2 gives me lambda 2y2, kyn is lambda nyn, and if I look to see what this is, this is the same as y1 to yn multiplied by lambda. If I just multiply [UNINTELLIGIBLE] on the right by lambda, it will take lambda 1 times the first column, lamda 2 times the second lambda n times the third, which is, times the last, which is what we want, so it's s lambda. This is all n eigenvalues and eigenvectors in one matrix equation, that's all that is. It's just ks equal s lambda i just says, this for all i at once, all i at the same time. OK.

So if I use these matrices in describing steps one, two, three, I'll see what's happening matrix language. OK, let me just do that. Step one, step one is looking for f as a combination of the columns of s. So step one is just f equals s times c. The vector of coefficients multiplies the columns of s and adds to give f. Then step two, which just divides everything by -- divides by the lambdas. Step two just creates lambda inverse, sc. So I took what I had -- let's see, no, the lambda inverse better be multiplying the c. Well, actually I can do it all in -- well step two, that's the easiest step, I should be able to do it. c is changed to lambda inverse c, so that c becomes lambda inverse c, OK. And then step three uses lambda inverse c to construct u. So step three is the answer u, what do I have here? I've got a combination of these vectors, so they're the columns of s, and what are they multiplied by? They're multiplied by the c's over lambdas, which is what I have here. That's s lambda inverse c. Those are the three steps.

And what's the work involved? Here, the work is solving a linear system with the matrix s. Here, the work is taking a combination of the columns of s, doing a matrix multiplication. Those two steps are usually full-scale matrix operation, and of course, the s , if I just complete this, I'll see that I get the right thing that's s lambda inverse and c is s inverse s. There's the answer. u is s lambda inverse, that's a lambda inverse there: s inverse f. That's the correct answer in matrix language. Right. This is k inverse, that's k inverse. k is s lamda s inverse, and if I take the inverse of that, I get s lambda inverse, s inverse to multiply f. Well, I doubt if you're much impressed by this lower board, because the upper board was the same thing written out. It took some indication of what the separate pieces were, but it's pretty clear.

OK, now the million dollar question is, is it fast? And the answer is, almost certainly no. But for the particular matrix s, which by good fortune s could also stand for sine, this matrix of eigenvectors for this particular problem are sines. These are the discrete sines, so this is the discrete sine transform. That's we're doing, we're doing the discrete sign transform, because those discrete sine vectors are the eigenvectors of k. OK, now let me say what that means. First I'm thinking of k in 1d. So this is my 2 2's and minus 1's and minus 1. Its eigenvectors are discrete sines, if I multiply that by sine kh -- well, let me just take the first one, sine h, sine 2h, sine n minus 1h, that will turn out to be an eigenvector. So this is ky, ky1, the first eigenvector. The eigenvectors are -- let me draw them. The eigenvectors for that second difference matrix k are -- here's the interval, 0 to 1, I chop it up in steps of h, and I plot the sine, which starts at zero and ends at zero, because those are the boundary conditions, and here is sine h, sine 2h, sine 3h, sine 4h, sine 5h, so for the five by five case -- maybe I should just be using n here, or maybe I should even be using capital n, so capital n is 5 in this example. Good.

What's on the other side of that equal sign? Some eigenvalue times the same vector, sine h, sine 2h, down to sine nh, and that eigenvalue -- oh, let me just write lambda 1 for it. We know what it is. The fact that this is an eigenvector is just trig; you know, I multiply minus 1 of that plus 2 of that minus 1 of sine 3h, and I use a little trig identity. So minus that, 2 of that, minus that, turns out that to be a multiple of sign 2h. Well sine 2h give us a 2, and then the minus sign h and the minus sine 3h combine into sign 2h times, I think, it's s cosine of h or something, it's that eigenvector that's near zero. But the cosine of h is near 1. Does that look familiar? That combination of sine h and sine 3h should give us twice sine 2h times some cosine, yep, Elementary trig identity. OK, so those are eigenvectors. That's the first one, the next one would have h times the k'th one would have h times k instead of just h itself, it would take jumps of every k sine, and then a cosine hk would show up there, and this would still be an eigenvector. OK. I'm making a little bit explicit these vectors, but the main point is they're sines, they're discrete sines at equally spaced points.

That's what the real version of the FFT just lives on, and it would also live on discrete cosines if we had different boundary conditions, we could do those, too. So this isn't the only -- these zero boundary conditions are associated with the name of Dirichlet, where zero slopes are associated with the name of Neumann, and both, this one gives sines, Neumann gives cosines, the FFT deals with both. OK. So, that's the fast solution, and it would take n squared -- well I have to go to 2d. sorry, I guess I have a little more to say because I have to get from this one-dimensional second difference to the 5.2-dimensional second difference, and that's what's going to happen over here. I wrote up some stuff about the Kronecker operation, which is the nifty way for these special problems to go from 1d to 2d. You remember the deal, that k2d, our 2d matrix was this Kronecker product of k and i, that gave us second differences in one direction, and then we have to add in the Kronecker product of i and k to get second differences in the other direction. And then we better print, because that matrix is going to be large, I don't want to [UNINTELLIGIBLE] it, yeah. What's the point?

The point is that if I know the eigenvectors of k, then I can find the -- if I know the 1d eigenvectors, I can find the 2d eigenvectors, and you don't have to know Kronecker products to do that. All you have to do is just make a sensible guess, so the eigenvectors in 2d, have a double index, k and l, and their components are sines in one direction times sines in the other direction. So what are those sines, there's a kh, I guess, lh. Those are the first components, I guess I have to tell you what all the components are: k, the seventh component in the x'th direction, there'd be a factor 7, so kmh sine lnh. This is the mn component of ykl. It's just what we had in 1d, in 1d there was no l, the components were just sine kmh. Now we've got two, one in the x direction. These are the analogs of the -- the continuous case would be sine, k pi x, times sine l pi y. Those are the eigenvectors as eigenfunctions, functions of x and y. And the point is that, once again, with these beautiful matrices, I can sample these at the equally spaced points, and I get discrete sines that the FFT is ready to go with, OK. I'm giving this much detail partly because the continuous case, of course, our operators d second by the dx squared and d second by dy squared, and the whole idea of separating variables, of looking for solutions, u -- here is the eigenvalue problem, the continuous eigenvalue problem. The Laplacian of u, maybe I do a minus, the Laplacian of u equal lambda u, and that's a partial differential equation, usually it's not easy to solve, but if I'm on a square, and I have zero boundary conditions, then I've solved it, by separation of variables, a function of x times a function of y. And that function of x times function of y is exactly what Kronecker product is doing for matrices, yep. I thought maybe this is good to know when the problem is easy, and as I say, the possibility of using the easy case on a square for preconditioning a not so easy case is attractive. All right, so that's what I wanted to say about number one, that's the main suggestion, and again, the point was just to take these three simple steps, provided we know and like the eigenvector. Here we know them and we like them very much, because they're those discrete sines.

OK, now just to finish comes, what's up with odd even reduction. I'll use the same 1d problem first. It works great in -- that's not good english, but it works very well in 1d odd even reduction, you'll see it, you'll see, oh boy, simple idea. But of course, don't forget that ordinary elimination is a breeze with a tri-diagonal matrix, so nothing I could do could be faster than that. But let's see what you can do. I just want to write down the -- OK, keep your eye on this matrix, so I'm going to write out the equations. So it'll be minus u, i minus 2 plus 2ui minus 1 minus ui, that would be fi minus 1; that would be equation number -- i minus 1, right with a minus one, two minus 1, centered there; and then the next one will be a minus ui minus 1 plus 2ui -- I better move this guy over further -- minus ui plus 1, that will be fi, right? That's equation number i, and now I just want to look at equation number, the next equation, minus ui plus 2ui plus 1 minus ui plus 2 is fi plus 1.

So I've written down three consecutive equations, three consecutive rows from my simple matrix: a row, the next row, and the next row. So the right-hand sides are coming in order, and the diagonals are there, and if I look at that, do I get any idea? Well, there is an idea here. I'd like to remove these guys, the ui minus 1's and the ui plus 1's, so that's where this word odd even reduction is coming in. I'm going to reduce this system by keeping only every second unknown and eliminating those. How to do that? Well, you can see how to eliminate, if I'm multiply this equation by 2. If I multiply that middle equation by 2, that becomes a 4, this becomes a minus 2, this becomes a 2, and now what shall I do? Add. If I add the equations together, I get minus u, i minus 2, these cancel, that was the point plus 4 minus a couple, so that's 2 ui's minus this guy, ui minus 2 is that sum fi minus 1, 2 fi's and fi plus 1. Well, sorry it's squeezed down here, but this is the point.

The main point is look at this. What do we have? We've got a typical equation, but now we've removed half the unknown. The problem is now half-sized. We've only got the even-numbered unknowns at a small cost in updating the right-hand side, and the problem's cut in half. So that's this odd even reduction, and it cuts a problem in half, and everybody knows right away what to do next. The one mantra of computer science, "Do it again." So that is the same problem on the even, we do it again, so I should really call it cyclic reduction, we just cycle with this reduction, and in the end, we have a 2 by 2 problem. That seems pretty smart, pretty successful move, and I guess if we do an operation count, well, I haven't thought that through. What does the operation count look like? It must be pretty quick, right? Well, undoubtedly we're solving this linear system in o of n steps. Well, no big surprise to be able to deal with that matrix in o of n steps, because elimination would take o of n steps, it's size n, but it's bandwith is 1, so 2n or something steps would do it, and maybe, I don't know how many steps we have here. I guess, when we cut it in half that required us to do that much. It's order n, it's order n.

So the key question is, can we do it 2d? Can we do the same thing in 2d? So I want to follow that plan in two dimensions where now u will be a whole row at a time, so I'm doing block 2d, block 2d. So, can I write down the equations in block 2d for whole rows, so ui is the vector? So now this is the 2d problem, so it'll be minus the identity, where instead of instead of 1, I have to write the identity, ui minus 2 and 2k, right? Oh no, what is the middle -- what's on the -- so multiplying ui, ui minus 1. I wanted to just say the same thing, but I have to write down the -- this is going to be n squared equations, n at a time, a whole row at a time, and what's on the diagonal of k2d? Not 2k. It's 2i, is it 2i plus k? Yeah. Yeah, k plus 2i. Isn't that what we have on the diagonal of k2d1 times ui minus 1 minus iui is equal to some -- that's a whole row at a time. These are all vectors with n components now, right? The minus i, the 2i, and the minus i are the second differences of one of rows; their difference is in the vertical direction, and this kui the second difference is along the row.

OK, so same equation at the next row. So the next row is minus i, ui minus 1, k plus 2ui, ui, because that's now the diagonal block minus iui plus 1 is fi plus 1, and it's just taking a second to write this one. This is k plus 2i, ui plus 1, minus i ui plus 2 is fi plus 2, OK. Exactly the same, but now a row at a time in 2d. So the same idea is going to work, right, what do I do? I want to cancel this, so I multiply that row by k plus 2i. Before I multiplied it by 2, but now I have to multiply it by k plus 2i times the row, the whole row. So when I do that this cancels, so I have minus this, this guy wasn't affected. And this will cancel, the k plus 2i will cancel this one, just as before. This will not be affected, ui plus 2, and I have fi and k plus 2i, fi plus one and fi plus 2. That should have been i minus 1, and this should have been i, and this should have been i plus 1, sorry, mis-labeled label the f's, but no big deal. The point is the left side, and the point is what's in that space. What goes in that space? Well, it's minus i, k plus 2i squared, minus i, so this is k plus 2i squared, which before was so easy, and then the minus i, and the minus i is the minus 2i is multiplying the ui. Yeah, that was [? my last i. ?]

OK, this is my matrix from odd even reduction. It's just natural to try the idea, and the idea works, but not so perfectly, because previously in 1d, that just gave us the answer 2. It was 4 minus 2, but now in 2d, we have a matrix, not surprising, but what we don't like is the fact that the bandwith is growing. k was tri-diagonal, but when we square it, it will have five diagonal, and when we repeat the odd even cycles. When we do it again. Those five diagonals will be nine diagonals, and onwards. So, I'm getting half-sized problems. All the odd-numbered rows in my square are eliminated, this just involves the even-numbered rows, but the matrix is not tri-diagonal anymore, it's growing in bandwith. So, and then you have to keep track of it, so the final conclusion is that this is a pretty good idea, but it's not quite as good as this one. It's not as good as the FFT-based idea.

Also, if you look to see what is most efficient -- see the eigenvectors are still here, so I could do three steps of this and then go to Fourier, and that probably is about right. So, if you really wanted to polish off a fast Poisson solver, you could do maybe two steps or three of odd even cyclic reduction, but then your matrix is getting messy and you switch to the fast Poisson solver, so it's not quite Poisson anymore, because it's has a messier matrix, but it still has the same eigenvectors. As long as we stay with the matrix k, we know it's eigenvectors, and we know that they're sines and that they're quick to work with. OK. Anyway, there you go. That's a fast algorithm for the lecture before the spring break.

So after the break is, first of all discussion of projects. If your project could include page -- and you could maybe email the whole project to Mr. [? Cho ?]. Maybe also, could you email to me a sort of summary page that tells me what you did, so I'll save the summary pages, and I'll have for Mr. [? Cho ?] the complete project with printout and graph, as far as appropriate, and so we'll spend some time on that, and then move to the big topic of the rest of the course, which is solving optimization problems, minimization, maximization in many variables. OK, so have a good spring break and see you a week from Wednesday. Good.