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, it mentions
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 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 multi grid, and I
computed the wrong thing.
You remember that I computed m. m was the jacobian matrix.
So m was the jacobian matrix, the iteration matrix with
jacobian which is i minus the 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 jacobian now.
And this is what we would get for ordinary simple iteration.
And it's not satisfactory.
Remember it had an Eigen value lamda max.
The spectra radius was about point 9.
Then I reported about the Eigen values m s, The multi
grid matrix sm and I was a little disappointed in it's
Eigen values.
The reason is, I should have been doing i minus sn.
When I do that, I'm not disappointed anymore.
What are the Eigen values of this?
You remember s is the matrix that tells us how much error
we're capturing. s is that matrix that came from the
error we captured at the end of multi grid, 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 multi grid does it smoother.
Steps two, three, four leave me i minus s.
And step five did another smoother.
And this is just with one jacobian step, which you would
probably do three.
Anyway the Eigen values of this were, well --
this has two zero Eigen values, and three 1's.
So this matrix has Eigen values 0, 0, 1,1,1.
Because, why's that?
Well, you remember the Eigen values had to be zero 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 Eigen value.
I was astonished to discover 1/9.
Of course, zero, zero 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 multi grid, it's not showing
it true colors.
But actually it does.
This is a kind Eigen value we expect.
Like 1/10, or a multi grid cycle.
So much better if we just did 3 m's for example, I would get
them down the point 9 q, which is more than point 7, where by
doing multi grid instead, we're down to point 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 3, no class.
So it'll be Wednesday that I'll see you again.
But Mr. [? Cho ?] will be available his 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 to 1:30, 1:00.
Otherwise taking 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 Eigen values.
Well, Eigen vectors.
Those will be in the notes too.
I guess I found one Eigen vectors, 1,2,22,1.
And glancing at the matrix, I didn't spot the other one, but
it has one oscillation.
It's 1,2,0, minus 2 minus 1.
Anyway, this is the main thing.
A much better result.
So that's multi grid, which we now can use as our iteration.
Or we can use it as a preconditioner.
Many of the recommended methods use multi grid as a
preconditioner before something even more powerful.
Well, I don't know about even more powerful.
Multi grid people would say not possible.
But some situations we may want to go to a different
method but we want a preconditioner,
and multi grid excellent.
By multi grid, 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 au equals b.
There's a Krylov matrix that created, exactly as you see
here, this is the js the one with j columns.
b, Ab, 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 span by is the same thing is saying all linear
combinations 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 I 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 ab.
After three steps, Jacobi produces some combination of
b, ab, a squared b.
In other words, all of those iterative methods are picking
their j'th approximation in this space, kj.
I should put a little j on it.
So 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 xj or
uj 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.
Lets not just use as simple iteration that built 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 xj --
since I think this section we use x, I'm going
to change that x.
They'll all the iterative.
They'll all started with xzero as zero 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 will be the conjugate gradient method.
CG.
The conjugate gradient method.
What determines xj?
It it chooses xj in kj, the space that we always look for
our approximation in.
Let me not forget to say these vectors, b, ab, 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
sparce matrixes.
And mostly, and especially, sometimes, especially
symmetric ones.
So just let me put in here, before I even finish that
sentence, 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 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 min res.
Maybe I'll just mention one more on our eventual list.
Actually that tells us probably what choice it makes.
Min res chooses the xj to minimize the residual.
Minimum r.
Remember r is b minus ax.
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
fantastic method.
Superior, just quicker than min res for a nice reason so
it chooses xj in kj and I think the rule it uses is that
I think xj should be orthogonal I'll go and check
it in my notes to the residual rj.
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 orthoganal to the whole space kj.
It turns out that we can choose xj.
When I make a choice of xj its in kj.
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, kj.
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 sparce 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 to be good.
Not in the normal.
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 orthoganalize the Krylov
basis, which is b, ab, and so on, producing orthonormal
basis q1, q2, up to qj.
It's just beautiful.
So arnoldi's taking like a preliminary step Krylov has to
make to get something that numerically
reasonable to work with.
Then if it's fast, 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 ab, in the direction of ab, which won't be
orthoganal 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 it's 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 it's projections on to the ones that are already set,
subtract those components off, you're left with a piece that
and you find it's 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 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 and
it's 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.
hj minus 1j.
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 plus 1j 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 the Legendre polynomials,
Chebyshev polynomials, all those classical things.
The reason is actually the same as it is here.
That some thing 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 symmetric and then you'll see its short
recurrence, and then we want to see what?
OK.
So I work 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, its diagonal.
So of course, to find a inverse b is not different.
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 k4, 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 r4 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 product methods, after four
steps, I will have the exact answer.
Why?
Because Krylov methods take the best x4 out of the space
spanned by those four columns.
Well, the whole space, r4 is spanned by those columns and
taking the best one must be the answer.
So x4 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 method.
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 part way, 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 ten 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 its 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
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?
Vandermon, yeah.
Vandermon.
So I could call it v for vandermonde.
And vandermonde matrixes are not very well conditioned.
So now a little time out to say, because it's so
important, how do you judge the good or poor conditioning
of a matrix?
Those vectors are linealally independent.
The determinant of v is not zero.
A matrix is not singular, but it's too close to singular.
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, --
I don't want to call them v because -- well I guess I
could call them v, they come out of Vandermon.
v1, v2, v3, v4.
That's our matrix.
Essentially I want its condition number.
And to find it's condition number, it's not symmetric.
So I can't just take the Eigen value of that.
You might say, look at Lamda max and Lamda min.
The condition number is associated with
max divided by min.
But when the matrix isn't symmetric, just taking its
Eigen values directly is not cool.
It's not reliable.
I could have a matrix that's badly conditioned, but all
it's Eigen values were 1.
I mean a matrix was 1s on the diagonal, and zillions up
above the diagonal would be have Eigen values of 1, but it
would be badly conditioned.
So the right way to take it is v transpose v. Look at v1
transpose v2 transpose, 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 Eigen values.
And those Eigen values, the Eigen values 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
Eigen values of v. It's the singular values you can trust.
And the way to find singular values is form v transpose v.
That give you a symmetric matrix.
Its Eigen values, so the i'th Eigen value would be the i'th
singular value squared.
So the condition number, is signma max over cignma 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 disasterious.
Completely disastrous actually.
10 by 10 vandermonde matrix, with 1, 2 3, 4, up to 10 as
the points 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 independent of the vs. 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 qs, if we've arnoldi and
got a good basis, it's one.
What's the gram matrix?
If this is q transpose q, with the qs the orthonormal basis
in the columns.
Then q transpose q is the identity.
So a grand matrix.
So q transpose q would be the identity, and it's condition
number is the best possible one.
Lamda max is lamda min is one.
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 10 arnoldi steps on the board.
The notes give 10 comments.
I'm quite proud my MATLAB is unreliable.
Actually you'll find a few errors like after end, I
accidentally put a semicolon [INAUDIBLE]
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 q2 is.
Why do we have a short recurrence?
This is a key point here.
And in this example, if I work out the hs, I'll discover sure
enough h13 will be zero. h14 will be zero.
Here's the key equation.
Here's arnoldi in matrix language.
Let me see if I can remember arnoldi in matrix language.
Yeah.
So arnoldi in matrix language is going to be this.
It's going to aq equals qh.
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 qh is a little bit like gram schmidt.
Do you remember gram schmidt is described by q times r.
q again is orthonormal.
So its 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 recurrence.
It's the 3 term recurrence relation is what I'm seeing
here in matrix language, because there were 3 non zeros
in the columns of h.
All right.
I was going to write out, I 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 q1.
q1 is the first column vector in the first basis vector.
And what do I have here when I take qs times hs?
Do you do matrix multiplication
a column at a time?
You should.
OK.
So this says take 5 halves of the first column.
And this says 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 aq2, the
next would an aq3.
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, as 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?
Its q transpose.
Anytime we see q, that's my letter. o is for
an orthogonal matrix.
So this is q transposed, aq.
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 matrixes.
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 transfer this
whole thing, what will happen?
To transpose things their transposes come in the
opposite order.
So q, it's transpose comes first. a, it's 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 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 its 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 kylov
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 then 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.