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 Poisson'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 x, y in the differential
equation, 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 coefficient
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 --
size is 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 e to the lambda*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 -- 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-zeros 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-zeros that are in the elimination,
in the factor L from elimination,
and it counts the number of non-zeros,
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 growing with.
I don't know whether the number of non-zeros -- 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.
So K*U 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 y_1, second eigenvector,
n-th 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 coefficients.
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 c_i by the eigenvalue lambda_i.
OK, so it's by eigenvector, so I'm assuming that K*y_i is
lambda_i*y_i 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 scalar;
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 c_1 over lambda_1 y_1 and c_2 over lambda_2 y_2 up to c_n
over lambda_n y_n, 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 lambda_i,
I'm getting the output here is U,
when the input was F. And do you see that that
is the correct answer?
All I have to do is check that K*U equals F.
So check that this answer from step three is the right answer.
OK, so I multiply by K. When I multiply y_1 by K,
a factor lambda_1 appears, the eigenvalue; it cancels that;
well that's y divided, so it would cancel,
and I have c_1*y_1.
When I multiply this by K, K*y_2 is lambda_2*y_2;
cancel the lambda_2's, and you're left with c_2*y_2,
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 y_1, y_2, y_n in its 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_1*y_1,
K times y_2 gives me lambda_2*y_2,
K*y_n is lambda_n*y_n, and if I look to see what this is,
this is the same as y_1 to y_n multiplied by lambda.
If I just multiply on the right by lambda,
it will take lambda_1 times the first column,
lambda_2 times the second, lambda_n 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 K*S 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 S*c.
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 operations,
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 F.
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*lambda 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 sine 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's and minus 1's and minus 1's.
Its eigenvectors are discrete sines,
if I multiply that by sine k*h -- well,
let me just take the first one, sine h, sine 2h,
sine n minus 1 h, that will turn out to be an eigenvector.
So this is K*y, K*y_1, 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 equals sign?
Some eigenvalue times the same vector, sine h, sine 2h,
down to sine N*h, 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 sine 2h.
Well, the 2 sine 2h give us a 2, and then the minus sine h
and the minus sine 3h combine into sine 2h times,
I think, it's a cosine of h or something,
it's that eigenvector that's near zero,
because the cosine of h is near 1.
Does that look familiar?
That a 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 h*k 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 five-point two-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 print 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, are -- so they 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 k*h, I guess, l*h.
Those are the first components, I guess I have to tell you what
all the components are: k -- the seventh component in the x
direction, there'd be a factor 7 -- so k*m*h sine l*n*h.
This is the (m, n) component of y_(k, l).
It's just what we had in 1D.
In 1D there was no l, the components were just sine
k*m*h.
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 eigenvectors.
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-2) plus 2*U_(i-1) minus U_i,
that would be F_(i-1); that would be equation number i
minus 1, right?
With a minus 1, 2, minus 1, centered there.
And then the next one will be a minus U_(i-1) plus 2*U_i --
I better move this guy over further -- minus U_(i+1),
that will be F_i, right?
That's equation number i, and now I just want to look
at equation number -- the next equation,
minus U_i plus 2*U_(i+1) minus U_(i+2) is F_(i+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 U_(i-1)'s and the
U_(i+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 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-2) --
these cancel, that was the point -- plus 4 minus a couple,
so that's two u_i's, minus this guy,
u_(i-2) is that sum F_(i-1), two F_i's, and F_(i+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 unknowns.
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 bandwidth 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 U_i 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 --
U_(i-2) and 2K, right?
Oh no, what is the middle -- what's on the --
so multiplying U_i, U_(i-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 the K2D one?
Times U_(i-1) minus I*U_i 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, they're difference is
in the vertical direction, and this K*U_i the second
difference is along the row.
OK, so same equation at the next row.
So the next row is minus I*U_(i-1), K plus 2 U_i --
because that's now the diagonal block --
minus I*U_(i+1) is F_(i+1).
And it'll just take me a second to write this one.
This is K plus 2I U_(i+1), minus I*U_(i+2) is F_(i+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, U_(i+2),
and I have F_i and K plus 2I F_(i+1) and F_(i+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 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 U_i.
Yeah, that was minus 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 bandwidth is growing. k was tri-diagonal,
but when we square it, it will have five diagonals,
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 bandwidth.
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 its 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 a 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.