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 hope in time for you to consider this
material also in thinking about projects because there
are lots of experiments to do with these.
I did a MATLAB experiment that I'll mention.
And I guess it was too small because it didn't make
multigrid look as good as it is.
So can I remember the steps of multigrid and then you'll
remember why we take those steps and then make a sort of
basic analysis of how do you decide whether this multiple
step process is successful?
How much does it reduce the error?
So you remember the steps.
Oh, thanks.
So the first step is really a little bit outside multigrid.
It uses an ordinary iteration like Jacobi, a few times and
so you remember that the matrix that governs this
ordinary step, the matrix that's
crucial for this is the--
so what multiplies the error at every step.
So we have initial error and if we take a smoothing step
with a preconditioner, p, then the error is multiplied by i
minus being p inverse A. So that's the matrix.
That's sort of that one step, one iteration matrix for
ordinary iterations where A is of course, we're solving
A sub u equal b.
That's our problem.
And matrix p is close to A. Of course if we took p equal to A
then the error would be reduced to zero right away
because we would be solving the system one shot.
But now, we're left with some error and this--
steps 2, 3, 4 reduce it further.
And well, our hope is to see that this multigrid operation
reduces the low frequency part because this one will act
successfully on the high frequency part, the
oscillating error, but for smooth error the eigenvalues
here are too close to 1 and it's too slow.
OK, so I was just going to track down-- well, first have
to remember, what are the pieces here?
So there's a matrix R, restriction matrix.
A matrix A sub Sh, which is original matrix A or A sub h.
I could call that.
But the original one is on fine mesh.
So I need a coarse mesh version.
So that's what I solve on the coarse mesh, the smaller
problem and then interpolate back with a matrix I back to
the fine mesh.
OK, so can I just remember also what those matrices--
we have all the pieces here.
So if I remember what R and I are--
OK, so let me write down all those pieces.
So A--
let me tell you the MATLAB experiment I ran.
So A is our two's and minus one's matrix.
Our second difference matrix in 1D to keep it simple.
And it should be divided by h squared.
Well, I guess we took the mesh where the coarse mesh just had
two interior points.
The value at the boundaries was zero and the fine mesh has
5 [INAUDIBLE].
So this will be 5 by 5.
Then we needed the restriction and interpolation matrices.
So this was R, the restriction.
So the restriction, which we're going to use right away
will take us from the fine to the coarse.
So that matrix will be--
it'll multiply things with 5 components and produce an
output with only 2 components.
Can always find either of those mesh points and their
numbers were 1, 2, 1, 0, 0, and 0, 0, 1, 2, 1 and needed
to divide by 1/4 because those rows added up to 4.
The idea being that if we restricted all ones we should
get all one and we need that 1/4 factor to do it.
And the interpolation went the other way and it was
essentially the transpose.
That's the nice connection to have and
again, we need a factor--
and maybe the factor is 1/2 there.
Good.
We've now got all the pieces except A sub 2h and I think
just at the very last minute last time I said what A
sub 2h should be.
And the good choice for A sub 2h--
so that's on the coarse mesh.
So that's going to be small matrix.
So the good choice is I, which puts it up on the fine mesh, A
sub h, which is our operator on the fine mesh and then R
that restricts back.
So R, A, I, is the right construction, it's sort of the
natural construction for the coarse mesh operator.
And notice, I don't know if you appreciate having this,
this is the same as R, A, R transposed with some factor so
that it'll come out symmetric, it'll come out positive
definite, and actually, it'll come out in this example-- it
will come out exactly the second difference operator on
the coarse mesh.
So this will turn out to be 2, minus 1, minus 1, 2.
The 2 by 2 operator and since we're on the coarse mesh the
step size is 2.
I won't do that multiplication, but MATLAB
confirmed it.
So that's where we were.
So I've sort of done now a recap of multigrid and I'm
ready to tackle a little understanding now of what
happens to the error.
So at this stage I end up with some error.
U sub h is not the exact answer.
I have an error e-- eh, which is u minus uh.
That's the error at the end.
Now the question is, what happens to that error?
Can we see how and why it's reduced as we do a v-cycle, as
we do a 2, 3, 4 cycle?
So since everything's linear that error is multiplied by
some matrix that we just have to track those steps, see what
that matrix is, and look at its eigenvalues.
So that's my goal here, is to find the matrix, call it s
that takes a step, a multigrid step and what does
it do to the error?
OK, so can I just follow that step?
Here's the key matrix.
Maybe I'll write it underneath.
So what happens to the error?
First of all, let me write the steps down. e sub h then is u
minus u sub h.
That's entering multigrid.
That's my error.
OK, now I compute this residual, r sub h.
I want to say that it's Ae sub h.
Can I check that that's true?
Yes.
A times this is A times u-- that's the b minus
Ah times u sub h.
That's the second term.
So that's OK.
So that's what happened there.
Then it got restricted.
I'll keep a record of what happens to e sub h.
So it's multiplied by A to get the residual.
We work with the residual and with b because b is
the input we know.
Then at this step two, that residual is restricted by R.
Now we've got r sub 2h.
When I multiplied by R I got r sub 2h.
OK, next step.
To get to E sub 2h.
So what do I do there?
I multiply by the inverse.
I solve this system, so of course, I don't actually find
an inverse matrix, but that's the formula.
To get this I multiply the inverse matrix by that, so
I've got an inverse matrix here.
So I write it as A sub 2h inverse.
This was A sub h.
Now when I've multiplied by A sub 2h inverse I'm this far
and the final step, multiplies that by I. So there you go.
It looks pretty awkward, but that's the matrix, s, that
produces E sub h.
So I'm thinking that my matrix is going to be called s.
This matrix here is going to be called s.
So let me write down what s is.
All I want to do is remind us what A sub 2h is by this
choice there, it's RAI.
And then comes R and then comes A. A is now A sub h.
Well, we've got a matrix.
And for this example I can compute what it is.
Can we say anything about that matrix from its formula?
Looks a little messy.
But of course, it's built out of these pieces.
It's built out of the given matrix A and our R and our I.
And there is a neat property hiding there.
Which you can see if you square that matrix.
So look at s squared.
If we look at s squared, so this is the matrix to
understand.
And I don't claim to understand it as well as I
would like, but--
so this is the matrix for these three steps.
We would have a full multigrid v-cycle just before I lose the
track on that.
A full multigrid v-cycle would do M a few times, say twice.
Two smoothers, then it would do a v-cycle
and then smooth again.
Well, I should've said the smooth again would be the one
on the left.
This is the original, so there's two smoothers followed
by a multigrid cycle, followed by a couple of smoothers.
That would be our total overall matrix, but I want to
focus for right now on this matrix, s.
And I guess, I hope that anybody who does a project or
a homework on this topic, you're going to create I, R,
and A and compute s.
And it should show what s squared is.
So could I compute, could I look at s squared?
You'll see why.
So I just write s down twice.
RAI inverse RA.
So that's once, now repeat the whole thing.
I, RAI inverse RA.
I have multiplied s by itself.
And what do you see?
Well it's nice.
If you look at this you see that we've got this
combination RAI sitting in the middle because when I did it
twice the RA came from there and then an I was at the
start, so this is--
so what's left?
S. This matrix S is equal to its own square.
S squared equals S. You know, we're
able to see this without--
this would apply for this example, but for all examples.
Now I'm ready to think about eigenvalues a little bit.
What does that tell me about the eigenvalues of S?
So this is a small class of matrices with this property
that s squared equals s and maybe I say something about
those here.
So if s squared equals s, what does that tell us about
eigenvalues of that matrix?
It's a square matrix.
It doesn't happen to be symmetric.
Very often we're looking at symmetric matrices, but this
is a little lopsided.
We have an A off on the right.
I might be able to make it symmetric by bringing part of
that A off to the left, but this is the main point.
What are the possible eigenvalues of that matrix, S?
Well, zero, 1, 1.
The only eigenvalues are zero and 1.
Why is that?
Because if Se equals lambda e, so I'll put that word if in.
If Se equals lambda e I could multiply by S again so I would
get S squared e equals lambda Se and that's
lambda squared e.
I haven't done anything fancy here, all I'm saying is that
any matrix, if it has an eigenvalue of lambda and an
eigenvector e then S squared has that same eigenvector and
it's eigenvalue is lambda squared.
But now S squared is S. So this is the same as Se, which
is lambda e.
You see that I now have lambda equal lambda squared.
Sorry, I probably made that take more time that it had to.
If S squared equals S then lambda squared equals lambda.
So I have lambda squared equal lambda and the only two
possibilities--
that's a quadratic equation that's got two roots and those
are the roots.
If lambda squared equals lambda then lambda could be
zero or it could be 1.
So what does that tell me?
That tells me that the eigenvectors that have
eigenvalue zero, what happens to those in multigrid?
They're killed.
An eigenvector with an eigenvalue of zero means that
Se is zero.
So whatever error we started with one cycle through
multigrid will remove it.
These are in the null space of S. They're vectors that S
completely removes and these are vectors where Se equals e.
The eigenvalue is 1 and multigrid does nothing.
Well, one part of the message then is that I could use
multigrid--
well, I guess I hadn't thought about this, but one message is
don't repeat 2, 3, 4.
That would be a large waste of time.
If you just repeated--
if you just did another v-cycle and didn't do any
smoothing or anything in between, then the error that
came out of the first v-cycle would survive right through.
Nothing good would happen in the second cycle that didn't
already happened in the first.
Now and another message is of course, I could do multigrid
forever and I'm not going to get convergance.
Actually, that's sort of interesting because we're
totally accustomed to working with matrices M, whose
eigenvalues are spread out, but below 1.
This has its eigenvalues not spread out at all, zero or 1.
That's it.
I'd better say quickly that the key matrix here is I minus
S. Just the way that it was I minus p inverse A. So the
multigrid matrix-- yeah, maybe I could've messed up.
The matrix is really I minus S, so let me be sure I see
where that happens.
I guess I got down here to figure out what E sub h is and
then I didn't do anything with it.
I forgot to write in the key step that at the end of step 4
when I have a piece of the error that I add that to Uh.
So add Uh plus this estimated error.
See if I knew E sub h, if I knew small eh-- the exact
error-- then when I added I would get U, the answer.
Of course, I don't know that.
What I know is capital E. So that's a cap E, which will
agree with e on some of these, but not others.
So all right, I guess I'm seeing the problem.
I may have headed in reverse because it's I
minus S that's crucial.
So these are errors.
Now what happens to them?
Are these the good ones?
I think they're the good ones, is that right?
They're the good ones because the errors in this stage,
multigrid gets them exactly, so that when I add this in
it's right on target.
It's accounted for.
So these are ones that multigrid fixes because of
course, everybody recognizes that I minus-- let me just say
it, though.
That I minus S times those errors is zero.
So the eigenvectors with eigenvalue 1 are the ones that
are perfectly captured--
E sub h has got those components of small e sub h
exactly right when I add them on; perfection.
But then there were other components which it
hasn't got at all.
And when I add those on, no change, no improvement.
By the way, so these eigenvalues are repeated.
Like here, I've got 5 by 5 matrix.
How many zeroes and how many ones?
You could say it's a projection matrix.
S projects on to these eigenvectors.
I minus S projects onto these eigenvectors and how many
zeros and how many ones for this if A is 5 by 5?
Well, what is A?
This A is 2 by 2.
So if I look at my S, look at my nifty formula for S, do you
see all this here?
Let's just get the sizes.
A is 5 by 5.
R is 2 by 5.
It shrinks us down to the coarse mesh, which is only
where this is only 2 by 2 and this takes us back up.
So the total result, 5 by 2, 2 by 2, 2 by 5, 5
by 5 is 5 by 5.
But what's the rank of S?
Small prize for telling me the rank.
I'm interested in this count.
How many vectors are killed?
How many go through unchanged when I multiply by S?
So what's your guess for the rank of S?
2, exactly.
So S will have 2 of these guys.
2.
This has multiplicity 2 in the 5 by 5 case and this one has
multiplicity 3.
So multigrid is only fixing 2-- getting 2
vectors exactly right.
And not getting these right.
Well, I mean, 2 is half of 5.
Well, in fact 2 is exactly half of 5 because if I went up
it would be 3 and 7, 4 and 9.
The point is that the number of coarse grid points is about
1/2 the number of fine grid points.
So 2 is about 1/2 of 1, 2 3, 4, 5.
Of course, how could I expect better then getting 2 perfect?
I can't expect more because the problem I ultimately solve
here is on the 2 by 2 size.
So when I actually do something useful, I'm only
doing it on the coarse grid and that's got dimension 2 and
so there will be 2 eigenvectors that get
perfectly fixed, but the others don't.
I do have to think which are the eigenvectors?
What eigenvectors are there?
I guess what I hope is, what I believe is that these
eigenvectors, which multigrid fixes for me, those E's--
well, the E's have 5 components.
Our matrix S is 5 by 5.
I guess I'm thinking, I'm expecting that the
eigenvectors that these 2 are somehow basically on the
coarse mesh.
And that these 3 that don't get fixed are somehow
basically orthogonal to the coarse mesh I would say.
That these will be low frequency vectors and these
will be high frequency vectors.
That's what I expect and that's what somebody--
I mean, there are whole books on multigrid, which where the
main step in the analysis is to see that those eigenvectors
are as I described, sort of coarse mash vectors that are
fixed and fine mesh vectors that don't get fixed by
multigrid and those are the ones that a
smoother better deal with.
I mentioned that just before class I did a MATLAB
experiment and I actually did it on these matrices with a
Jacobi smoother.
And I broke down the numbers here somewhere.
and actually I was a little disappointed.
But I came away saying, OK, 5 by 5 just wasn't big enough.
So I'll put down what I've got.
What did I get for the eigenvectors?
eigenvalues of S I got 1 and zero.
That was a good thing.
Now did I write them down?
I sure hope I did.
Yep.
So the question is, what about when M is in there?
So the exercises were find the eigenvalues--
well, it'd be interesting first of all, to know the
eigenvalues of M. This is the Jacobi
smoother, weighted Jacobi.
Here's one we're really interested in.
We're interested in doing a smoother, doing a multigrid
cycle, and doing a smoother.
I'm hoping that those eigenvalues are below this.
So this is one smoother and this is smoothers with
multigrid and this would be 2
smoothers and 2 post smoothers.
Repeating Jacobi twice and M cubed would be interesting.
Here is what I got for these things.
What do I know about the eigenvalues of M?
Do you remember that we actually worked those out?
For ordinary Jacobi they're cosines and for weighted
Jacobi that weight omega, which I chose to be
2/3 comes into it.
Anyway, so in a 5 by 5 case-- well, it's the biggest one
that we all really want to know.
So this is maximum eigenvalue.
Let me just write down the maximum.
So the maximum eigenvalue of M was 0.91.
So that says that ordinary Jacobi on a 5 by 5 problem
reduces the error by that factor at every step.
And of course that factor is totally satisfactory.
If we could always do that we wouldn't need to do multigrid.
But that 0.91, if I go to 7 by 7 or 9 by 9 it's going to
quickly move up toward 1.
That's an indication that my little experiment here is a
bit small because that's not as near 1 as a true
[? vector. ?]
OK, what was the result after this?
I'm sorry to say that it was 0.88.
Disappointing.
So I did this one.
I did 2 smoothers.
Well, of course that's going to give me
0.91 squared in there.
But anyway, I got 0.64.
Which of course is very good.
I'll say very good, but actually for that size I had
hoped for better because my impression of multigrid from
reading about it is that--
I mean that a good multigrid cycle, but this has an
eigenvalue of about 0.1, is the goal.
So how is that achieved?
Well, I don't want to achieve it by just doing a whole lot
more smoothing.
Actually it's achieved by doing a whole lot more
multigrid instead of doing just a little v-cycle this
would be the same for a good, deep v-cycle that went down
maybe 4 and maybe a w-cycle; do it twice.
Then of course, this is for a big problem, but really making
multigrid work like a w-cycle.
Or even there's something called full multigrid.
Should I draw the picture that corresponds
to w for full multigrid?
So full multigrid you start on a coarse mesh.
Go up, back up, up, back, back, maybe twice, up, up, up
to this level, back, back, back, and up again.
So that's called full multigrid.
That gives the best numbers.
Anyway, we get numbers like this for when we really do
more and it's very much worth it.
I mean, it's completely a success.
I would be interested to see a little--
well, I guess what I'm going to do now is, for this model
problem that we worked on I could find these
eigenvectors, e.
And analyze this model problem.
OK, so how to find these eigenvectors.
Again, for these particular matrices.
Well, let me start with the eigenvectors of A itself.
Do we know the eigenvectors of A?
I guess if I had one important matrix, whose eigenvectors--
and nontrivial.
I can't look at the 5 by 5 there and guess at sight, say
what the eigenvectors are.
So everybody sees it's 5 by 5 with 4 minus ones below, 5
two's on the diagonal, 4 minus ones above the diagonal, and
the eigenvectors of those special matrices have a
special form, which allows us to do everything here.
So what can I do here?
I can do the computation of what S is.
Can I report on that, which will be in the notes.
Let me tell you what S is.
If I take that as a R, that as I, that as A sub 2h I get S to
be this matrix.
Column of zeroes, 1/2, 1, 1/2.
Column of zeroes.
1/2, 1, 1/2 and column of zeroes.
So it's rather remarkable in this particular case.
Now what can we recognize from this?
Well can we see eigenvalues of S explicitly from the matrix?
Well, we certainly see that it has rank 2.
It has 2 independent columns so its
rank is 2 as we predicted.
So it will have 3 zero eigenvalues and we can see
what the eigenvectors are actually.
What are the eigenvectors if I go up above here?
Can I do that?
I want to go up and say something about the
eigenvectors with eigenvalue zero.
What's the null space of that matrix?
If I multiple S by vectors I want to get zero vectors.
The eigenvector's e for lambda equals zero.
If I multiply by anything with a first component, with these
three components, you see that if I multiply Se is zero there
for these e's.
In a way that's exactly what I hoped for.
That these vectors, e are the vectors that live on the fine
grid, but with no component on the coarse grid.
And they're the ones which have Se equal zero and
multigrid doesn't help.
Multigrid doesn't help.
I minus S times those vectors is those vectors, no
improvement.
So they're eigenvectors with eigenvalue zero for S, but
with 1 for I minus S. See my general thought is those would
be the high frequency eigenvectors,
high frequency vectors.
They're sort of high frequency, but they're not low
frequency anyway.
Because low frequency these would be a smooth vector.
OK, so what are the eigenvectors with the other
two eigenvectors--
the ones that multigrid fixes?
OK, so I'm looking for the other two eigenvectors here
and how do I find them?
Let me see if I can even do it for this matrix.
I'm looking for the eigenvectors with
eigenvalue 1 now.
So I should subtract the identity from this matrix.
When I subtract the identity the ones become zeroes and the
zeroes become minus ones.
That's S minus I now, trying to keep a running equation
that's correct.
That's S minus I. And I'm looking for the null
space of that now.
Am I going to find it?
I don't know.
Do see immediately what-- that has a
two-dimensional null space.
There are 2 vectors and can you see what they are?
What do you think?
Yeah, we certainly could.
I don't know whether--
anybody see a vector that's in the null space of that matrix?
Let me just start and try to get one.
So I'm looking at that matrix, looking for a vector in its
null space.
And that will be an error vector that
multigrid totally fixes.
Say if I start with a 1 up in that first component then the
first row will force me to have a 2 there.
Is that right?
And then I fixed the first row, the second row is all
zeroes already.
Then I have 0.5 of that, which is 1.
Oh, well let me start from the bottom here.
Well, it's probably symmetric.
What do you think?
Maybe something like 1, 2, and what's that middle number?
Maybe 2.
I think that that is in the null space of S minus I
because if I multiply by the first row I get minus 1, plus
1 and the third row gives me 1, minus 2, plus 1 and yes--
so that's one of them anyway.
And the point about it is it's pretty smooth.
That's what I anticipate vectors that are smooth,
that's the type of vectors we're going to
find multigrid fixing.
So I've got one more eigenvector to find.
Maybe I'll leave that as a puzzle, somewhere there's
another eigenvector.
Guess I'm not sure where it is, but I think
it should be smooth.
Somewhere there is one.
Can I now go away from this numerical example to give you
the answer for this problem, but for any size?
Oh, it's because we know the eigenvectors of A. I started
to talk about the eigenvectors of A and I didn't finish.
If I want to get formulas for these I need to know the
eigenvectors of A. So A, this matrix of
twos and minus ones---
the eigenvectors of that matrix are
discreet sine functions.
So an eigenvector would be for example, sine of pi on 6, sine
of 2 pi on 6 up to sine of 5 pi on 6.
It's worth knowing that those are eigenvectors.
And other eigenvectors would double the frequency, triple
the frequency, four times the frequency and five times.
The eigenvectors of this matrix are what I will call,
discrete sines.
Discrete sine vectors.
And why?
Because this matrix builts in the boundary condition since
there should be a minus 1 there, but it's chopped off.
That minus 1 if it was there would
multiply the zeroes component.
But the zeroes component is removed, so we have to be sure
that the zeroes component would have been zero anyway.
And that's the point.
You see if I follow this pattern the zero component
that I'm missing would be sine of zero.
The component I'm missing up here would be sine of zero pi
over 6, which is zero.
And the component that's missing down at this boundary
would be sine 6 pi over 6, which is sine
pi, which is zero.
That choice of vectors is the perfect Fourier vector.
I expect Fourier when I see constant coefficients there
with the right boundary conditions zero at both ends.
Anyway, those are the eigenvectors.
The k'th eigenvector will have sine k pi, sine k2 pi, down to
sine k5 pi, all over 6.
And I'll just, in the remaining moment write down
the result of following through those vectors and
finding finally, the eigenvectors of S itself.
So this is like multigrid in a nutshell if I can find it.
Multigrid in a nutshell.
All right, better be on this final board.
I minus S times this sine vector , the k'th sine vector
that I wrote over there and I minus S times the n minus k'th
sine vector.
So these will be low frequencies for small k and
these will be high frequencies.
And it turns out that these are 1/2, 1 minus the cosine of
k pi over 6 or n plus 1 times yk plus--
sk, sorry.
Sk plus S sub n minus k.
And this acting on the high frequencies is 1/2 of 1 plus
cosine k pi over 6.
Sk plus S sub n minus 1.
What am I seeing there?
Let's just accept the algebra that went into this and see,
does that tell me what the eigenvectors of I minus S are
because that's my goal.
Well, it does.
Look, if I add these two equations, let me look first
just at this equation.
So this says that if I start with a low frequency, k small,
low value of k that it's nearly removed.
If k is very small the cosine of k pi over 6 is
somewhere near 1.
Of course, 6 was a little feeble, so let me take n;
bigger matrices.
Then the low frequencies are the ones where this
cosine is near 1.
This number is near zero and they're practically removed.
Where the high frequencies, not much happes because the
cosine is near 1, this number in near 2.
I divide by 2 and it's near 1.
So the overall picture is right.
Now to get to pin it down, if I add that two equations then
I get this combination Sk plus S minus k.
Oh in fact, I guess I see it.
What happens if I add the two equations?
I get I minus S times Sk plus S sub n minus k and what do I
have on the right-hand side?
Well, 1.
These add to 1.
So I get Sk plus S sub n minus k.
Those are the eigenvectors of I minus S with eigenvalue 1.
They're the ones I found here.
They're the eigenvalues of S with eigenvalue zero.
These are the Sk's plus S sub n minus k.
And they live completely on the mesh points that are in
the fine, but not the coarse mesh.
And then if I play around one more step I can find the
eigenvectors for the other eigenvalue.
The ones that are basically low frequency and where we
have a numerical example.
OK, I've run out of time before I can do that last bit
of algebra.
That will be in the notes, which we'll get onto the
website, Monday at the latest.
So I hope you have a sort of picture of what multigrid does
and the experiments that you might do to check it out.
OK, I'll see Monday and we'll take a little more time to
discuss projects because we'll have material on
the web to work with.
OK, thanks.