The following content it's 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: Several pieces of good news.
Since Brett is back now, all of chapters five and six
of my notes, that is to say, all that we've covered in class,
plus a little more, are on the web page now.
So, I was cautious at the beginning of the semester
and didn't put that up, and then forgot that it wasn't there.
So now you have something to work with,
particularly on any experiments, for example,
minimum degree ordering.
Oh and also, there's a movie now.
Tim Davis from Florida and [? Perils ?] [? Person ?],
my friend here, created a little movie to show the order that
elimination occurs.
And I had a lot of email with Tim Davis.
On that one-page handout, for 2D Laplace's equation,
it mentions N cubed for the nested dissection, where
you cut the region in half, and half, and half, by separators.
I believe, and I think experiments will show,
that minimum degree is even better.
But the analysis is extremely weak, apparently,
on minimum degree.
So anyway, one possible project either for now
or for the final project would be some experiments
to see what the actual convergence is like.
And maybe a little understanding of minimum degree.
So we'll talk more about that.
Anyway, there's a movie and a lot more material
there on the website.
It'll get updated and corrected as I do the edits.
Second bit of good news.
I made a mistake.
That's not normally good news, but this time, I
reported a little calculation for multigrid,
and I computed the wrong thing.
You remember that I computed M. M was the Jacobi matrix.
So M was the Jacobi matrix, the iteration matrix with Jacobi,
which is I minus D inverse A This is
for the 1D second differences.
In fact, this was the 5 by 5 matrix
that we've been talking about.
And now this is multiplied by D inverse which is a half,
but then with a weighting factor, that becomes a third.
So that's weighted Jacobi now.
And this is what we would get for ordinary, simple iteration.
And it's not satisfactory.
Remember it had an eigenvalue lambda_max, the spectra radius,
was about 0.9.
Then I reported about the eigenvalues of M*S,
the multigrid matrix S, M, and I was a little disappointed
in its eigenvalues.
The reason is, I should have been doing I
minus S M. When I do that, I'm not disappointed anymore.
What are the eigenvalues of this?
You remember, S is the matrix that tells us
how much error we're capturing.
S is the matrix that came from-- the error we captured
at the end of multigrid, was S times
the error that we entered with.
And so the error that's left, the remaining error,
the new error, is the rest.
It's the part we don't get, and that's when I forgot.
That's I minus S e.
So that's why I should have been using I minus S.
So then, the step one of multigrid does a smoother.
Steps two, three, four leave me I minus S.
And step five did another smoother.
And this is just with one Jacobi step, which you would probably
do three.
Anyway, the eigenvalues of this were, well --
this has two zero eigenvalues, and three 1's.
So this matrix has eigenvalues 0, 0, 1, 1, 1.
Because, why's that?
Well, we get two 0's-- you remember the eigenvalues had
to be 0 or 1, because S squared equaled S,
and I minus S squared equals I minus S. So.
We remember that.
Now, the question is, what are these?
And the answer is 1/9 is a triple eigenvalue.
I was astonished to discover 1/9-- of course, 0, 0 survived.
The rank is 3.
And this is now looking much better.
What I reported last time was some more like 7/8,
or something, and I thought oh, multigrid,
it's not showing its true colors.
But actually it does.
This is a kind eigenvalue we expect, like 1/10,
for a multigrid cycle.
So much better.
If we just did three M's, for example,
I would get down to 0.9 cubed, which is more than 0.7,
where by doing multigrid instead, we're down to 0.1.
OK.
So that's the good news.
And, of course, if you experiment a little,
you'll find different numbers for different sizes
and really see what's happening.
And, in fact, you could do 2D.
I mentioned all this partly for now,
if you're wanting to do it now, or partly for eventually,
later.
Oh, and thinking about projects, just one word to repeat,
that the Monday after spring break,
I think that's April the third, no class.
So it'll be Wednesday that I'll see you again.
But Mr. [? Cho ?] will be available.
He has office hours, so that's April 3,
if you wanted to discuss issues with your project
with Mr. [? Cho ?], that would be at class time in his office
2-130, at one o'clock.
Otherwise, take an extra spring break,
and I'll see you Wednesday.
OK.
So that's various bits of good news.
Oh.
One other thing that I guess I didn't finish in that lecture,
was to identify the two eigenvalues.
Well, eigenvectors.
Those will be in the notes too.
I guess I found one eigenvector, [1, 2, 2, 2, 1].
And glancing at the matrix, I didn't spot the other one,
but it has one oscillation.
It's [1, 2, 0, -2, -1].
Anyway, this is the main thing.
A much better result. So that's multigrid,
which we now can use as our iteration.
Or we can use it as a preconditioner.
Many of the recommended methods use
multigrid as a preconditioner before something even more
powerful.
Well, I don't know about even more powerful.
Multigrid people would say not possible.
But some situations, we may want to go to a different method
but we want a preconditioner, and multigrid excellent.
By multigrid, I include the smoother.
So it'd be that.
That would be a possible preconditioner
for something else.
OK.
So where are we now?
We're into the new section starting this moment.
And what's the section about?
It's about things associated with this name, Krylov.
So, I'm going to use the letter K for the things that are
associated-- and I'm always solving A*u equals b.
There's a Krylov matrix that's created exactly
as you see here.
This is the j-th-- the one with j columns.
b, A*b, A squared b, up to A to the the j minus 1 b.
And the Krylov space is the combinations of those vectors.
That's what this word, span, means.
I could erase span, and say, "all combinations of".
Maybe that even more familiar.
So "spanned by" is the same thing
is saying "all linear combinations of".
Of those j vectors.
And, of course, that's the same as saying
it's the column space of the matrix,
because the column space is all combinations.
OK.
So, why am I interested?
Why was Krylov interested?
Why is everybody interested in these vectors?
Because actually, that's what an iteration like Jacobi produces.
If I use Jacobi's method-- or Gauss-Seidel, any of those--
after one step, I've got b.
After two steps, there's a multiplication by A in there,
right?
And some combination is taken, depending
on the particular method.
So after two steps, I've got a combination of b and A*b.
After three steps, Jacobi produces some combination of b,
A*b, A squared b.
In other words, all of those iterative methods
are picking their j-th approximation
in this space K_j, actually.
I should put a little j on it, to indicate that that's-- so we
have these spaces are growing.
They grow by one dimension, by one new basis
vector at each iteration.
And the point is, Jacobi makes a particular choice
of x_j or u_j, the approximation after j steps.
But does it make the best choice?
Probably not.
In fact, pretty definitely not.
And so the idea is let's make the best choice.
Let's not just use a simple iteration
that builds in a choice that might not be so terrific.
Let's make the best choice.
There are a whole lot of Krylov methods, which all choose x_j--
since I think this section will use x,
I'm going to change that to x.
They'll all be iterative.
They'll all start with x_0 as 0, say.
And then the first guess will be b, maybe.
Or maybe not.
Maybe some multiple of b.
Actually, a good Krylov method will
take the best multiple of b as the first guess,
and not necessarily 1 times b.
And onwards.
So we have this word, best, coming up.
What's the best vector in that space?
And there are different methods depending
on what I mean by best.
Oh, let me tell you the name of the method
that I'm going to concentrate on first and most.
Will be the conjugate gradient method.
CG.
The conjugate gradient method.
What determines x_j?
It chooses x_j in K_j, the space that we always
look for our approximation in.
Let me not forget to say: these vectors, b, A*b, A squared b,
and so on, they're easy to compute,
because each one is just a matrix multiplication from
the previous one.
And the matrix is-- we're assuming we're
working with sparse matrices.
And mostly, and especially, sometimes,
especially symmetric ones.
So just let me put in here, before I even
finish that sentence, that CG, this
is for A transpose equal A symmetric.
It's only for those.
And positive definite so symmetric, positive definite.
So it's a limited class of problems, but a highly, highly
important class.
And you may say what do we do if the matrix is
symmetric, but indefinite?
Well, that comes next.
Or what if the matrix is not to symmetric at all.
Well, if you're brave, you might try conjugate gradients anyway.
But if you're cautious, then you would
use one of the other methods like MINRES.
Maybe I'll just mention one more on our eventual list.
Actually that tells us probably what choice it makes.
MINRES chooses the x_j to minimize the residual.
Minimum r.
Remember r is b minus A*x.
The residual r will always denote
the error in the equation.
And so it's minimum residual.
OK.
So that's a natural choice, but you'll
see that this is a fantastic method.
Superior, just quicker than MINRES, for a nice reason.
So it chooses x_j in K_j, and I think the rule it uses
is that I think x_j should be orthogonal-- I'll go and check
it in my notes-- to the residual r_j.
Let me just be sure-- if I'm going to write that down,
I better be sure I said it correctly.
I actually, I didn't say it correctly.
r_j is orthogonal to the whole space K_j.
It turns out that we can choose x_j--
when I make a choice of x_j, it's in K_j.
Now when I compute r, there's a multiplication by A
to get the residual, so that moves us up to the next space,
and I'm going to make a choice where
this is orthogonal to the whole space K_j.
You have to do conjugate gradient.
And I can't start on that for a good reason.
I have to start with something called Arnoldi.
And and what's Arnoldi about?
Well.
Let me come back to these vectors.
Quick to compute.
But what's the other property?
Everything in applied numerical mathematics,
you're choosing basis vectors, and you're
looking for two properties.
And I guess this is like a general rule, whenever
you meet a whole new problem, in the end,
you're going to construct some basis vectors
if you're going to compute.
And what properties are you after?
You're after speed.
So they have to be quick to compute and work with.
Well these are.
Multiplying by A is sparse matrix multiplication,
you can't beat that.
But the other property that you need
is some decent independence.
And not just barely independent.
You want the condition number of the basis,
somehow-- if I use that word-- to be good.
Not enormous.
And best of all, you would like the basis vectors
to be orthonormal.
That's a condition number of one.
That's the best basis you can hope for.
And the point is, that this construction
produces a lousy basis from the point of view of condition.
So Arnoldi's job is orthogonalize the Krylov basis,
which is b, A*b, and so on, producing orthonormal basis
q_1, q_2, up to q_j.
It's just beautiful.
So Arnoldi's taking like a preliminary step
that Krylov has to make to get something that numerically
reasonable to work with.
Then if it's fast, and Arnoldi doesn't-- we have to do
a multiplication by A, and you'll see one here
in the middle of Arnoldi-- and then, of course,
if you look at those ten lines of MATLAB,
and we could go through them a little carefully,
but let's just take a first look.
The first step is like Gram-Schmidt, right?
It's sort of a Gram-Schmidt idea.
You take that first vector b, you accept that direction,
and the only step remaining is to normalize it.
Divide by its length.
Then you go on to the next.
Then you have some trial vector t,
which would in this case be A*b, in the direction of A*b,
which won't be orthogonal to the original b, right?
So this won't be orthogonal to that one.
So, of course, you've got to compute
an inner product between them.
And subtract off the right multiple
of the previous one from the current t
to get the improved t.
And then you're going to normalize it again.
So you compute its length and divide by the length.
You see that overall pattern in Arnoldi?
It's the familiar idea of Gram-Schmidt of,
take a new vector, find its projections
onto the ones that are already set,
subtract those components off, you're left with a piece,
and you find its length and normalize
that to be a unit factor.
It's very familiar.
Its exactly the type of algorithm
that you would write down immediately.
And it just involved the one multiplication by A
so that it's not going to cost us
a lot to make the vectors orthonormal.
Well.
Wait a minute.
Is it going to be expensive or not?
That's like the key question.
It's certainly going to produce a good result.
Producing orthonormal vectors is going to be a good thing to do.
But is it expensive or not?
That depends.
It's certainly not expensive if the new trial vector
t has only components in one or two
of the already settled directions.
In other words, if I only have to subtract off
a couple of earlier components, then I'm golden.
And that's the case when A is symmetric.
So that will be the key point here.
And I just make it here.
If A is symmetric, A transpose equals A, then I only need,
in this subtracting off, where I'm headed for the new q,
I only need to subtract off h_(j, j), multiplying the q_j,
the q that was just set, and the one before.
h_(j-1, j).
I'll call that a short recurrence.
It's short because it only has two terms.
And then there'll be a new h_(j+1, j),
which is just the length.
So there'll be one of these magic things in so many parts
of mathematical analysis, a three-term recurrence relation
will hold-- in fact, I guess the reason
we see three-term recurrence relations in all-- Legendre
polynomials, Chebyshev polynomials,
all those classical things-- the reason
is actually the same as it is here.
That something in the background is symmetric.
I want to put a few numbers on the board
so you see what is typical input and output to Arnoldi.
And you'll see it's symmetric and then you'll
see it's short recurrence, and then we want to see what?
OK.
So I worked out a typical Arnoldi example.
Not that one.
Here.
OK.
I think this is a good example to look at.
So matrix A is not only symmetric, it's diagonal.
So of course, to find A inverse b is not difficult.
But we're going to do it through conjugate gradients.
So that means that we don't figure out A inverse,
which, of course, we easily could.
We just use A, very sparse.
Here's our right-hand side.
Here's our Krylov matrix.
It has b in its first column.
It has A times b in the second column.
Multiply by A again to get the third column.
And A one more time to get the fourth column.
So that's K_4, right?
And for a 4 by 4 problem, that's the end.
So actually here I'm carrying Krylov to the end of it.
There's no more to do.
Because once I've got this K_j, the combinations of those,
are all of R^4, all of four-dimensional space.
Oh, yeah.
That raises an important point about how
we're improving on iterations.
If I was using Jacobi, or Gauss-Seidel, or one
of those others, then-- well, I won't
say for this particular A, but usually-- after four steps,
I am by no means finished.
I've taken four steps, I've gotten closer,
but I haven't got the exact answer.
But now, in this the world of Krylov methods,
after four steps, I will have the exact answer.
Why?
Because Krylov methods take the best x_4 out of the space
spanned by those four columns.
Well, the whole space, R^4 is spanned by those columns
and taking the best one must be the answer.
So x_4 will actually be A inverse
b, the answer we're looking for.
So these methods are on the one hand, direct methods.
They finish after four steps, finish.
You've got the answer.
Now.
That's actually how they were born.
As direct methods that gave the answer
from an interesting iteration and they nearly
died for that reason.
They were born and died, or were on death's door,
as direct methods.
Because as direct methods, there are better ones.
Then some years later, people went back to it, back
to conjugate gradients, and noticed
that thought of as just going partway,
gave very successful answers.
So we're going to think of Krylov-- we're planning
to stop before j equals n.
In this picture, j equals n equal 4 here.
But, we plan to stop for j much below n.
We're thinking of applications where n is 10 the fifth,
and we're going to stop at 10 squared, 100 steps, maybe.
Or ten steps.
So that meant a total reconsideration,
reanalysis of conjugate gradients,
and it is now back to a big success.
Well it's a rather unusual thing that a method that people
have put aside gets picked up again and becomes a favorite,
as conjugate gradients have.
First of all, I was saying that the basis, those basis vectors
are not very independent.
That's not a good basis.
Of course, it's a detractive basis and maybe,
do you know the name for a matrix
that has this particular form?
The columns are constant, then, that's the key column, there,
and then it's first powers, second powers, third powers.
Do you remember whose name is associated with that matrix?
It comes up in interpolation, and it starts with a V?
Anybody know?
Vandermonde, yeah.
Vandermonde.
So I could call it V for Vandermonde.
And Vandermonde matrices are not very well conditioned.
So now a little timeout to say, because it's so important,
how do you judge the good or poor conditioning of a matrix?
Those vectors are linearly independent.
The determinant of V is not zero.
The matrix is not singular, but it's too close to singular,
and how do you test the-- suppose you
have a basis, as we have here.
And you want to know is it good or not?
So, of course, you always put your bases vectors,
let me call them-- well,I don't want to call them v because --
well I guess I could call them v,
they come out of Vandermonde. v_1, v_2, v_3, v_4.
That's our matrix.
Essentially I want its condition number.
And to find its condition number,
I-- it's not symmetric, so I can't just
take the eigenvalues of that.
You might say, look at lambda_max and lambda_min.
The condition number is associated
with max divided by min.
But when the matrix isn't symmetric,
just taking its eigenvalues directly is not cool.
It's not reliable.
I could have a matrix that's badly conditioned,
but all its eigenvalues were 1.
I mean a matrix with 1's on the diagonal
and zillions up above the diagonal
would have eigenvalues of 1, but it would be badly conditioned.
So the right way to take it is V transpose V.
Look at v_1 transpose, v_2 transpose...
As always, if a matrix isn't symmetric,
if the matrix V is not symmetric,
good idea to form V transpose V. That is symmetric.
It does have positive eigenvalues.
And those eigenvalues, the eigenvalues of V transpose V,
are the singular values, or rather the singular values
squared, of V. So I guess I'm saying,
you can't trust the eigenvalues of V.
It's the singular values you can trust.
And the way to find singular values
is form V transpose V, that gives you
a symmetric matrix, its eigenvalues--
so the i-th eigenvalue would be the i-th singular
value squared.
So the condition number, is sigma_max over sigma_min.
And, well, it's not enormous for this 4 by 4 matrix,
but if I go up to 10 by 10 or 100 by 100.
100 by 100 would just totally wipe me out.
10 by 10 would already be disastrous.
Completely disastrous actually.
10 by 10 Vandermonde matrix, with 1, 2 3, 4, up to 10
as the points, would be-- well, it would have entries--
if that ended in a 10 and I had 10 rows,
would that be something like 10 to the ninth power.
I mean the dynamics scale would be terrible.
So finally, just to-- V transpose V
is called the Gram matrix.
So that guy Gram is coming back in as giving the good measure.
So the point was then, that this measures the dependence
or independence of the v's.
That ratio.
The bigger that ratio is, the more dependency they are.
What's the ratio if they're orthonormal?
What's the ratio of the q's?
What's the condition number of the q's?
If we've run Arnoldi and got a good basis, it's one.
What's the Gram matrix?
If this is Q transpose Q, with the q's, the orthonormal basis
in the columns, then Q transpose Q is the identity.
So it's a Gram matrix.
So Q transpose Q would be the identity,
and its condition number is the best possible one.
Lambda_max is 1, lambda_min is 1.
The condition number is one.
OK.
So that's just some comments about why Arnoldi gets brought
in to fix the situation.
OK.
So I'll just leave those ten Arnoldi steps on the board.
The notes give ten comments.
I'm quite proud of-- my MATLAB is unreliable.
Actually you'll find a few errors
like after end, I accidentally put
a semicolon, which was absurd.
But the comments I'm pleased with,
and then the numerical example that
runs through one cycle of this, with these numbers,
to see what q_2 is, is, I hope, useful.
Why do we have a short recurrence?
This is a key point here.
And in this example, if I work out the h's, I'll discover sure
enough, h_(1, 3) will be 0. h_(1, 4) will be 0.
Here's the key equation.
Here's Arnoldi in matrix language.
Let me see if I can remember Arnoldi in matrix language.
So, Arnoldi is taking the matrix-- yeah.
So Arnoldi in matrix language is going to be this.
It's going to A*Q equals Q*H. I can't write out all of Q.
So that's the big equation.
Its a very important equation that we now
have all the pieces for.
So A is our original matrix that we were given.
Symmetric let's say.
Q is our basis out of Arnoldi and H is the multipliers
that gave that basis.
So this Q*H is a little bit like Gram-Schmidt.
Do you remember, Gram-Schmidt is described by Q times R. Q,
again, is orthonormal.
So it's an orthogonal matrix.
In Gram-Schmidt R is upper triangular.
Here it's not.
Here it's Hessenberg.
So H stands for Hessenberg.
I'll write down what the actual H
is for these for these numbers.
I won't write down the Q. I'll just write down
what H turned out to be for those numbers
if I did it right.
Five halves.
Oh, interesting.
The lengths all turned out to be five halves.
And this turned out to be root 5 on 2.
This turned out to be the square root of 4 over 5,
and this turned out to be the square root of 9 over 20.
And the point is from here, this one is just below the diagonal.
And it will show up as symmetric.
Root 5 over 2, root 4/5 and root 9/20.
OK.
So, what am I seeing from that particular H which
somehow can't be an accident?
It must be that it's built in.
It's the fact that H is symmetric and tridiagonal.
And what does that tridiagonal tell us?
It tells us that we have short recurrences.
It's the three-term recurrence relation,
is what I'm seeing here in matrix language,
because there were three non-zeros in the columns of H.
All right.
I was going to write out-- I was going
to try to understand this one by looking at the first column.
What if I take the first column of both sides?
That'll be A times Q_1, right?
Q_1 is the first column vector in the-- first basis vector.
And what do I have here when I take q's times h's?
Do you do matrix multiplication a column at a time?
You should.
OK.
So this says take 5/2 of the first column.
And this says take that factor times the second column.
And I could track back and see, yes that's
what Arnoldi has produced.
And then the second one, the next one, would be an A*q_2.
the next would an A*q_3.
Well, look, here it is.
I want to show that H is symmetric when A is.
Could you do that?
We know what the property of Q is here.
We know that Q is, because Arnoldi made it that way,
has Q transpose Q equal I. And I can, in one step,
show that H is symmetric if A is symmetric.
How do you do that?
I guess we need a formula for H, so I just
multiply by Q inverse.
So that's good.
And even better is to recognize what Q inverse is.
So what is Q inverse here?
It's Q transpose.
Anytime we see Q, that's my letter, always,
for an orthogonal matrix.
So this is Q transpose A Q. And now what?
The argument's finished.
We're here.
If A is symmetric, what can you say about that combination?
If A is a symmetric matrix?
It is symmetric.
Right?
That's how you get the symmetric matrices.
You start with one.
You multiply on one side by a matrix, and on the other side
by its transpose.
The thing has to be symmetric, because if I
transpose this whole thing, what will happen?
To transpose things their transposes
come in the opposite order.
So Q, its transpose comes first.
A, its transpose comes in the middle, but what's that?
The transpose of A is A. We're assuming A to be symmetric.
And then the transpose of Q transpose is?
Is Q.
So we got this back again when we transpose,
so it's symmetric.
So H is symmetric.
So the conclusion is, H equals H transposed.
And then we know immediately that it's tridiagonal,
because every H from Arnoldi is Hessenberg.
We know that these zeros are here.
The Arnoldi cycle ended produced h i j in column j
but it's stopped one below the diagonal.
So we know these are zero, but now, if we know
the matrix is symmetric, then we know these are zero.
And, of course, it works out that way.
So the conclusion is that we can orthogonalize the Krylov
basis, quickly, easily, and work with that basis,
either explicitly by computing it
or by implicitly by keeping things orthogonal
and that's what conjugate gradients will do.
So next time, I'm going to make an effort
to describe the conjugate gradients method.
I'll pick the highlights of it.
It has fantastic properties, and to verify those properties
in full detail is often more confusing than not.
If you see today's lecture, you're
seeing the important points: the role of symmetry
of A and the Arnoldi algorithm.
OK. so that's our first lecture on the Krylov ideas.
And next time, we'll probably complete that topic,
and it will be on the web.
So I'll see you Wednesday for the end of Krylov.
Thanks.
Good.