Okay. Good morning, everybody. I'm not sure if there's anybody new here, but this is scientific computing for Dphil students,
the second of 12 lectures in this term and then 12 more in the next term.
Let me remind you of a few things.
One is that there's an assignment due next Tuesday, which certainly requires you to make sure you've got MATLAB working.
Once you've got MATLAB working, it's not a very large assignment.
And you turn them in in the boxes across the atrium here.
160 or 161. Let me also remind you that we have a course Web page.
So there is the Web page, and I hope you'll get in the habit of looking at it in particular.
It's worthwhile just going and after a lecture, if you like, either usually not posted long in advance,
but just go to the lecture and spend a few minutes looking over the fairly terse notes that should be there.
I hope it works not working. Well, I think it works eventually.
My laptop has a funny viewer. Let's make sure it's there.
Well, it's there somewhere. So you'll find that the lecture notes are fairly terse.
They're typically just five or six pages of outline like this,
but it's certainly worth spending some minutes looking that over to make sure you're in touch with things.
Now let's see. Today there are three handouts, three new handouts.
The one with the blue is just two more copies of what we handed out last time.
But there's three new pages. Let's take a moment and look at the one called Journals and Books.
So this is just a list of some journals and books that are well known to numerical analysts like me.
You remember, one of my goals in this course is to immerse you a little bit in a field that isn't your own.
So you get a sense of how we think and what our points of reference are.
Well, these are some of our points of reference. So I've listed a dozen of the journals that are well known in numerical analysis
and a couple of dozen books that are relevant to this term of the course.
Most of this stuff is online. In fact, this is a copy of a HTML file that I'll probably put up on the website.
So most of these things are freely available in one form or another.
For example, let me just mention one or two things on this list of books that are particularly interesting.
There are the two books called Templates for Eigenvalue Problems and for Systems of Equations,
and those are kind of a fun way to see those problems through the software window of a kind.
The book by Bjork is a standard book on least squares problems with many interesting things in it.
The book by Davis is a important book on direct methods for sparse, big linear systems of equations.
The book by Government Van Loan is sort of the Bible of numerical linear algebra.
The book by Know Sedol and write numerical optimisation is absolutely the standard text worldwide for numerical optimisation.
And then you'll notice books by Saad and Vandervoort. Those are two of the leaders in Iterative Methods,
and I think you'll find both of their books online for linear equations and also eigenvalue problems.
So Saad and Vandervoort are two very big names in this area.
Okay. Let's take a look at the quiz. I hope you've all done the quiz.
Problem one What is the inverse of Spurs? What's the answer? In full.
That's right. So let's see here. The quiz. The answer to one is C.
What's the inverse of N? It's enough.
That's right. So the answer to two is D, which inequality is invalid?
What's that? The second one. That's. That's right. It's the second one that should be a backslash or a backslash.
Then it would in principle be valid, but it's a forward slash, which is something else entirely.
So that's the invalid one. And finally, did anyone compute the output of program in number four?
You all did, I'm sure. But let me tell you, it's a two by two matrix containing the numbers, two in all four positions.
Okay.
So today we're talking about conjugate gradients, which is one of the greatest hits of, well, really computer science computing, numerical analysis.
It's really one of the great algorithms of all time whose impact is absolutely huge conjugate gradients.
So we call this section 1.4. It's an idea that had a very clear beginning at the beginning of the 1950s.
There's a article getting on in age now on the history of conjugate gradients that's kind of interesting.
And that's by Gallivan O'Leary. In 1989, I guess it is Gottlieb and O'Leary in Siam Review in 1989.
Anyway, the, the starting point of the algorithm was two people hesston us and Stiefel, who independently had more or less the same idea.
In fact, there were three key people. Heston, he Stiefel and Lantos were all doing very closely related things at the beginning of the 1950s,
and in particular, Heston, even Stiefel communicated in time that they ended up writing the main paper together.
So this is the Great Paper Yesterdays and Stiefel in 1952.
And certainly by any standard, that's one of the ten most important numerical analysis papers ever written.
Let's see, there's a handout which shows you a bit of the flavour of that first paper.
So methods of conjugate gradients for solving linear systems.
I like to give you some flavour of the past. It's a remarkable paper, partly because the two of them collaborated on it.
It's amazing how much is in this paper.
People keep rediscovering parts of what's in this paper because of course, not so many people read the actual paper anymore.
It's a remarkable paper, but it also has an oddity,
which is that it can't quite decide whether conjugate gradients is an iterative method or a direct method.
And that's a flavour we'll be talking about a lot in the next few weeks.
And today in particular, a direct method is an algorithm where you do a certain amount of computation and then you're done.
You've solved the problem. An iterative method is one where as you compute, you get closer and closer to the solution.
And once you've got six or 16 digits, then you stop.
Very different style. So a direct method, you sort of think of it as going like this, an iterative method you think of as going like this.
The paradox of the conjugate gradient algorithm is that in principle it's both.
If you do enough work, you solve a matrix problem X equals B, but also as you go, you get closer and closer to the solution.
So has Farnese and Stiefel were very proud of inventing an algorithm that could do both of these things.
And as a result, they wrote a paper which didn't focus on the important one, which is the iterative side.
So if you just look at this first page that I've copied, the second word in the abstract is iterative.
So obviously they know. It has an iterative flavour.
On the other hand, if you look at the second column and you look at item A,
it's pointing out that like Gaussian elimination, the method gives the exact solution after N steps.
So they're hedging their bets. It's iterative, it's direct. And it was funny what happened historically.
They wrote this great paper. People didn't focus very much on the iterative side of conjugate gradients.
They focussed on it as a direct algorithm. And the fact is it's not a good direct algorithm.
It's nowhere near as reliable as Gaussian elimination.
So it sort of lost the battle in the fifties and sort of got put aside.
Everybody knew about it, but people didn't use it much.
And it really took 20 years for conjugate gradients to be kind of rediscovered as an iterative algorithm.
I'd like to draw you a plot to indicate why that happened.
It wasn't just that people were being dumb. It has to do with how things evolve with time.
So this is time, but it's also the dimension of a problem.
Because of course, as computers get more powerful, you can solve bigger problems.
So this is the 1950s, this is 2015.
And of course, n gets bigger with time.
As we discussed last lecture, this access is the amount of work involved in solving an end by end system of equations.
So where solving X equals B and it's end by end.
Now if you do this with the obvious method of Gaussian elimination.
Which is too familiar for me to want to talk about it.
The amount of work is NN cubed, basically. So that's a fine algorithm.
And back in 1950, let's say, and might have been ten in 2015 and might be ten to the eighth or something.
Now conjugate gradients depends on the matrix, but loosely speaking, it's got a picture that looks more like this.
Let me put it in quotes. I'll say C, G, and I'll put in quotes O of n squared.
That's not a general figure. There's there's no sense in which it is generally of complexity and squared,
but maybe that's a typical amount of work involved when things are going pretty well with conjugate gradients.
On the other hand, the constant is works. If you have a small matrix of five, ten or 20 or 30, you would never use conjugate gradients.
So it really didn't make sense in those early days for people to use conjugate gradients.
Somewhere in the early 1970s, and let's to be precise, say 1971,
these two curves crossed over and people started to notice that for big problems, at least conjugate gradient was a very powerful idea.
So it was a remarkable change of perspective there.
The person who made that observation was John Reid, who is still around.
He's at the Rutherford Appleton Laboratory in ten miles south of here.
He wrote a paper which pointed out how powerful conjugate gradient was.
And the funny thing is all the foundations are already there and has been is in freefall in 1952 but didn't get focussed on.
It's a very curious bit of history. In fact, one reason I'm so interested in this is that you'll often hear people say that numerical computation is,
let's say, ten to the 12th times faster than it was way back when.
And then they'll say six of those orders of magnitude come from computers and the other six come from algorithms.
And in some sense that's true. The algorithms are much faster and the computers are much faster.
But still, there's an asymmetry because the computers drive the algorithms, you know.
Yes, it's true that we didn't focus on conjugate gradients back here, but that's not really because we were stupid.
It's because it didn't make sense to focus on it back here. So in some sense, as the machines get faster, that affects what the right algorithms are.
And in the end, you get your tremendous improvements, you get a curving of a plot like that because of the interaction of machines and algorithms.
Anyway, conjugate gradient is a huge thing. I checked on the web on if you Google search it, you get 573,000 hits today.
The third hit is a very popular one because it's a paper with a nice title.
It's a paper by Jonathan Schuchat called an introduction to the Conjugate Gradient Method Without the Agonising Pain.
Chuck. Chuck, by the way, is a well known guy at Berkeley who won the Wilkinson Software Prize for his package called Triangle.
Okay. So what I want to do today is show you the conjugate gradient method and derive a
certain part of the theory that explains why it's so spectacularly fast sometimes.
So it applies to symmetric, positive, definite matrices.
So let a which is n. By n. B.
Well, we always say speedy, which means symmetric, positive, definite.
Symmetric means if you transpose it, you get the same matrix.
So symmetric means a equals a transpose and positive definite.
Well, you can define that in various ways. One way to define it is to say that all the eigenvalues are positive.
A real symmetric matrix always has real eigenvalues, but they might be zero or negative.
A real symmetric, positive, definite one has positive eigenvalues.
Another way to say it is that x transpose a x is always greater than zero.
If x is a non-zero vector. So X transpose x is a very standard way that we would write a thing like this.
You should think of that as a picture of a row vector.
Times a. Matrix times a. Column vector.
A lot of the notations we use are equivalent to though different from what physicists in quantum mechanics use.
So physicists like to talk about bras and cats, which have to do with sort of left and right vectors,
and they would sort of split this operator into a left half and the right half.
This is the way numerical analyst like to talk about such operations.
Anyway. In applications, symmetric positive, definite problems come up all the time,
essentially because her mission differential equations come up all the time.
Self Adjoint Differential Equations. When you have a self a joint differential equation, then almost certainly you can find a way to democratise it,
which turns it into a symmetric, positive, definite matrix.
So we want to solve x equals B, so let B be an end vector.
And of course, we're interested in the case where end is pretty big. So we want X Star is what I'll call it, such that a x star equals B.
That solution exists. That is unique because of the assumption that the matrix is positive.
Definite. That implies that it's non singular. All the eigenvalues are positive, so it's a non singular matrix.
There's a unique solution. So X Star is our unique solution.
And the way conjugate gradients works is it starts from an initial vector, an initial gets x nought and it iteratively gets closer to x star.
So c g finds x star iteratively.
Starting from some vector which I'll call x not.
Now. It could be anything you like, but we're going to take it to be zero.
It doesn't have to be zero, except it simplifies the formulas.
So without. Well, with only slight loss of generality, let's say that x nought is the zero vector, the N by one zero vector.
And the way conjugate gradient does this it finds x one in the next two and x three and these iterates x
little n belong to very interesting sub spaces so they so x and belongs to the so-called cry of space.
Which is a subspace of capital and dimensional vectors.
And the way we write it is we'll write K or K seven is the space spanned by the
right hand side vector B and A times the right hand side vector dot dot dot up to.
A race to the end, minus one power times the right hand side vector.
So that idea of a krylov subspace is absolutely at the heart, not just of conjugate gradients,
but of all sorts of iterative algorithms for solving systems of equations and for finding eigenvalues and singular values, all sorts of things.
So let me say a little more about that. This notation means the set of vectors spanned by these vectors, linear combinations.
So if I just have one of them, the first kind of subspace, that's a one dimensional space.
If I have two of them, that's probably a two dimensional space and so on.
So this isn't generically an end dimensional space which is contained in the set of all end vectors.
And as little then increases, you get a bigger and bigger subspace of your giant space.
So conjugate gradients is searching in bigger and bigger spaces for a good approximation to the solution.
Now let me show you some of its properties.
So first of all, we define the so-called a norm of a vector to be the size of the vector measured in a special way.
We say that the a norm of a vector x is equal to the square root of x transpose x.
So that's a number and it's a positive number if X is non-zero.
Because by definition, A is symmetric, positive, definite, and therefore this thing is a positive number so we can square root it.
That's going to be the measure of our error. And let's define the error at each step like this.
I'll say that the error the error at step n I'll define to be sub n.
And that will be the exact solution. X star minus the gas explosion and the residual, as we call it, is what you get if you apply A to that.
So it's our seven hour residual equals A times seven.
Now eight times X star is B minus eight times x n.
So throughout numerical analysis, you always have these two ideas of an error and a residual and residuals, a kind of a dual of the error.
It's what you get when you apply the error applied to the error, the problem you're trying to solve,
the operator you're working with as a rule in principle, you would like a small error, but the thing you can measure is very often the residual.
So you minimise the residual as a means of minimising the error all over numerical analysis.
Okay, so let me move over here and write down an amazing property that I'm not going to prove of conjugate gradients.
So I'll call this theorem one. And the amazing property is that CG conjugate gradients magically find the very
best vector in this space in the sense of minimising the norm of the error.
All right. So at the end of the span, yes.
So the set of all linear combinations of these vectors.
Other questions. I didn't come to this extent. How do I come up with it?
I haven't told you yet. Yeah, I'm going to tell you. It's an amazingly slick and obscure.
It's not at all obvious when you see the formulas, how it works. But I'm about to tell you, of course.
But first, let me tell you the theorem, which is that among all vectors in this kind of space, I guess I should have given it a subscript.
So k seven is the end dimensional cry of space.
The one that conjugate gradients is going to give us when I tell you the formulas is the unique one that minimises the norm of the error.
87. Measured in the a norm. That's the theorem.
So this optimality property which has two needs and Stiefel of course, knew about.
So at every step conjugate gradients find you the best approximation in a certain measure in a space.
And those spaces are growing. They're getting bigger and bigger,
which means that at some rate we must be approaching the correct solution because after n steps will have an end capital n dimensional subset of this.
So we will have it all. So once we if we really take capital and steps, we've got to solve the problem.
Nobody would ever take capital n steps in practice anymore. So indeed,
we have immediately from this theorem the beautiful corollary that the initial error
measured in the a norm is greater than or equal to the first error in the norm,
and so on. So the errors keep going down monotonically, at least in the norm.
Very beautiful stuff. I think they also go down in the two norm, but that's less obvious, if I remember.
But this is the fundamental property. So I'm going to now give you the formulas.
And as I do that, let me let me tell you the context in which I like to present them.
The way most people think of conjugate gradients is as a kind of an optimisation search.
You're in this giant 10,000 dimensional space.
Let's say you're searching for a particular vector in this space,
and you do that by starting with one vector and then moving along in some direction to improve the gas and then moving in another direction.
So you're searching around in this huge dimensional space and conjugate gradients can be described very compactly in terms of how it does that search.
So let me be specific. So I'll call it the CG iteration and let's in particular relate it to the idea
of minimising something which we always call optimisation we want to minimise.
A quantity, and that's equivalent to finding a place where a certain gradient is zero.
Right. You imagine a surface, you're trying to minimise it. That's that's amounts to making some gradient equal to zero.
Let's be specific. Let's define F of any vector X to be one half times x, transpose a x minus x, transpose b.
That's a scalar. We already saw this scalar that has the structure of rho vector times, matrix times,
column vector x, transpose b has the structure of rho vector times column vector.
So this is a scalar and you can compute that the gradient of that scalar function of x is equal to A, x minus B.
How do you do that? I won't go through that. But because X appears twice, you get two terms, you add them up, they're the same.
That's so the one half gets doubled and turns into a one. And then the gradient of this term is B, so you have this 10,000 dimensional space,
you have the scalar function of this vector in that space, the gradient of that function is a vector, and this is that vector.
Now we want to be equal to B, which means we want the gradient to be zero, right?
So our goal is to make that zero.
In other words, to make the gradient zero, in other words, to minimise F of f's f of x.
So F takes a minimum uniquely.
At the point ex star.
I'm not very good at drawing pictures in a thousand or 10,000 dimensions, but in two dimensions anyway.
Of course, you imagine you've got some kind of a parabola here.
Well, that's one dimension, I guess you'd say. And we're trying to find.
That Poindexter. At which the parabola is minimised.
More generally, it's a parabolic surface in capital and dimensions, so the zero gradient gives us the minimum.
So what conjugate gradient does is to cleverly search to minimise that function.
And now I write down the formulas. So here they are.
And we're going to treat them as magic.
I'm not going to tell you where they come from, but I will tell you how to interpret each line of the algorithm.
So so here's the algorithm. We start with an initial guess, which, as I say, we're going to take to be the zero vector.
Corresponding to that is a residual and the residual will be it's B minus A times the initial vector, but the initial vector zero.
So it's B, so the initial residual is B.
And then we're also going to work with a search direction vector, which we'll take to be the initial residual.
So in our thousand dimensional space, we're going to be searching in all sorts of directions.
And this is a vector that will define those search directions.
Now, here's how the algorithm goes for little n equals 1 to 3, and so on until we're happy.
Will define Alpha seven to be the quotient.
Of R and minus one transpose R and minus one divided by P and minus one transpose a pound minus one.
So that's a scalar divided by a scalar.
So Alpha is a scalar. I'm going to interpret each of these expressions in a moment.
Xand Well, of course you know what the interpretation of that is. It's the next guess, and that's X and minus one, plus alpha, minus one.
And then our end, the next residual is equal to R and minus one, minus alpha n, a p and minus one.
And then there's another scalar called beta or beta if you're English,
which is the quotient of our n transpose r n divided by R and minus one transpose R and minus one.
And then the final line is we update the vector p by the formula p.
N equals r. N plus beta MP and minus one.
And so that's just magic formula.
It's not at all obvious why they work.
But let me at least tell you what they're doing. So you can get the optimisation interpretation here.
You already know that we're in a big dimensional space and we're trying to get better and better approximations to the point X star.
Well, alpha being a scalar is the step length, which means it's how far we're going along a search direction to improve our current guess.
Xin is of course the approximate solution at step n.
Our seven is of course, the residual corresponding to that.
Well, it's obvious from this formula how we can interpret better.
Subban. It's the improvement. This is the square of the residual at the last step, and this is the square of the residual at the current step.
So this measure is, in some sense the improvement in the step we're now taking.
A scalar. And finally, this is a search direction.
So we're in this high dimensional space and at every step we take a step from the current position to our new position along this direction.
So this is a search direction.
And maybe the most interesting formula to look at in this magic algorithm is this one where you can see that at every step to get the new estimate X,
then we take the old one and then we take a step in the direction point minus one of length alpha.
And so it's crazy and magical, but this has the property that Theorem One says it has of giving you an optimal solution in a life space.
So I want to show you that a little bit numerically, and then I want to say a bit about what that implies about convergence rates.
So we're now at the MATLAB handout. There's a sheet of paper there with CG and M three on it.
Now of course,
real software for conjugative gradients would be much fancier than this with all sorts of bells and whistles and convergence criteria and so on.
But this little code, SIG, is designed to show you that iteration in its simplest form, and you can see it has exactly the structure.
I just wrote down an iteration involving five or six update formulas at each step.
We should say a word about where the work is involved in this iteration.
Most of this is just manipulating vectors.
The only part of the CG operation that is more than just a vector is the multiplication of the matrix A by the step P and minus one.
So that's a matrix vector product. And it's the same here and here.
So you only do it once at each step. So every step of conjugate gradients involves a matrix vector product.
If that were a dense matrix, that would be and squared work. If it's a sparse matrix, it's less than that.
Now this particular example I'm going to run or do the first half of M three convergence.
So let's take A to B, a diagonal matrix of dimension 100 with entries going from 100 to 199 along the diagonal.
So it's a trivial matrix. Mathematically, it's just diagonal 101 up to 199.
So the size of a is 100 by 100. And let's take a trivial right hand side to be the vector.
What's happening here. The vector of length, 100 of ones.
And now let's say X equals C G of A, comma B.
So I'm just running this code. Now what you see is these are the residuals at each step and you see they're beautifully converging.
Now, what does that mean? We start with an error of one, and after 20 steps or so, we've got about ten digits improvement.
It's 100 by 100 matrix. And yet it only took us ten or 20 steps to get this very nice improvement.
So the amount of work we've done is small.
Of course it's a trivial problem. We could have done it on paper, but you see the idea.
And if I look at the solution, of course, it's accurate to ten digits or so.
Let's look at the first four entries. The diagonal matrix was 100, 100 and 102.
So the first four entries are one over 101 over 101, one over 102, and so on to ten digits.
So there we just used conjugate gradient to solve a rather trivial 100 by 100 problem.
Let's do another one here. So let's say A equals the diagonal matrix consisting of the vectors from let's see, I'm going to put in six zeros there.
One, two, three, four, five, 61999999.
Hope that's right. So I've now got a million by million matrix represented sparsely.
Of course, again, this is a problem we can do on paper because it's just a diagonal matrix.
But it illustrates the idea. Let's take a right hand side to be the trivial vector of length.
2 million. And now let's again do X equals C, g of a B.
So this is a million by millions system of equations. And you see it chugs along with essentially the same performance.
So this gives you a hint of how powerful conjugate gradients can be when things are working.
So we get solved a million by million system of equations. Admittedly, it was diagonal.
Okay, so I want to say a word about convergence.
The convergence of conjugate gradients is an example where some really beautiful mathematics has some very concrete consequences.
And the part of mathematics we're talking about is approximation, theory, approximation of polynomials, of certain things.
So let's call this section 1.5. This is not well in this course.
Everything we do is mathematics in some sense, but mostly it's not really a mathematics course very much.
Maybe the next 10 minutes are a little more mathematical than some.
So it turns out that just from the optimality property of Theorem One,
you can get everything that I haven't proved that theorem, but from that everything follows.
And let me remind you that the initial error is x da minus x, not which we've taken to be zero.
So that's x star. Similarly, the error would be x star minus x, and so that would be E, nought minus x n.
And I want to say something about how small the end is.
So. Using our theorem one optimality.
We well we we note that X seven, as I've told you, is a member of this cry loft space.
So it's in the space B A, B up to a, n minus one B but B is eight times X star, eight times E not so.
Another way to write that is that it's a times E, not a squared times e not dot dot dot up to A to the n e not because eight times e not is B.
So these are just two different ways to write the same little n dimensional space.
So examine our guests at the step. At the end step is in this space.
In other words, it's a linear combination of these vectors. And now there's really only one central idea to this analysis.
A linear combination of these vectors can be interpreted as a polynomial of a applied to inward.
In fact, the set of all linear combinations of these vectors is precisely well,
it's the set of all polynomials of degree, little n applied to E, not with a with no constant term in that polynomial.
So in conjugate gradients we're implicitly dealing with polynomial of a times initial error.
And that's why the analysis turns out so beautifully. So this tells us that e sub ne7 is e not minus x ten and x and is in this space.
So e seven is equal to a polynomial of degree.
N in a times e sub not where.
P seven is a polynomial. Of degree and or at most end because of this linear combination of these terms.
But we know one other thing about that polynomial, which is that the degree zero term, a constant term just comes from there.
So it's one. So it's a polynomial of degree.
N With the property that P of zero is one.
So I've now told you really the crucial things about conjugate gradients, the rest is just working out the consequences.
The crucial things are. Two. One is this theorem that at each step we're optimising the norm in a certain space.
And the other is this bit of algebra that being in that space amounts to being equal to a polynomial and a times the initial error.
All the rest is fiddling with polynomials.
So that's really a beautiful, elegant property, I think.
So let me restate theorem one.
I'll call it Theorem one prime. It's the same theorem just restated E and is equal to some polynomial and a times e not.
And it's the one that's optimal. So that is to say where that point of a in art is minimal.
And we have to say, of course, it's minimal among polynomials p of degree and with P of not equal to one.
So that's what conjugate gradient is magically doing at every step.
Now there's one other piece I have to bring in in order to really show you how it's an approximation problem that determines the convergence rate.
And that other piece is the eigenvalues of The Matrix. So I hope you all know about eigenvalues.
I'm just going to use them without really talking about them much.
But in order to get quantitative about what this means, we need to look at the eigenvalues of A and the eigenvectors of A.
Now as a symmetric matrix.
A symmetric matrix always has real eigenvalues and orthogonal eigenvectors, so the eigenvectors provide an orthogonal basis of the space our capital.
N So that's what we do now.
We say Let a have eigenvalues. I always write that e w it's German.
I convert for eigenvalue and eigenvectors e v So let a have eigenvalues,
which by the way we know are going to be positive Lambda one through Lambda Capital, F and eigenvectors, which are orthogonal.
In fact, let's make them all the normal. They won through Van.
Now the thing to do is expand the vectors we care about in this basis of all the normal eigenvectors.
So we then expand and we say that e not.
Well, any vector can be expanded in this basis. So we're going to write it not as the sum from one to capital n of some number a sub k times v.
It's a basis you can always expand in the eigenvector basis.
And now we look at the polynomial of a applied to a knot and we find that n is equal to the polynomial of a times e not.
Now the whole definition of an eigenvector is that a times vk equals Lambda K VK.
That Generalises a polynomial and a applied to VK is that polynomial of lambda k multiplied by v k.
So when I do this I get the sum from one to capital and of a k times.
PM of Lambda K applied to VK.
So now. We've made the connection between a polynomial and a particular basis.
We want this to be small, right? Well, how could this thing be small if all of these are small?
So that's basically where the approximation theory comes in.
Conjugate gradients will converge quickly if there exist polynomials that go to zero quickly at the eigenvalues of the matrix.
And thanks to the miracle of Theorem one, it's always finding the perfect polynomial.
It's always optimising. So if there exist good polynomials, the algorithm will find them.
So that's our conclusion.
If there exist polynomials that are very small at the eigenvalues.
Then C g will converge quickly.
And this is where you get the kind of picture that is often drawn in write ups of conjugate gradient.
People like to draw a picture of a polynomial, let's say P of lambda.
It always has to go through one at the origin.
But we want it to be small at the eigenvalues. So A is a big matrix with a thousand or a million eigenvalues and they're all positive.
So they're on the real axis somewhere, and there's just a lot of them. We?
The question we ask in considering the convergence of conjugate gradients is does
there exist a polynomial that passes through this point and yet is very small?
And all of these eigenvalues. And of course, I've been writing these things down quickly, but I've sort of told you the whole logic of the argument,
if that polynomial exists and thanks to Theorem one prime or one,
we know that there will be a small error because conjugate gradients will find that polynomial.
Okay, so that's what it all comes down to. And I just want to show you a corollary and one more computation and then we'd better stop.
So what this all leads to, I'll call it theorem two.
Is that the ratio of the error at the 10th step to the initial error.
When you do conjugate gradients is bounded by the minimum overall polynomials with the usual
restriction times the maximum overall eigenvalues of the matrix times point of lambda.
Okay. So as always, pain of zero has to be one.
Now there's more to say, but I don't really have time to say it. Let me show you a final demo.
So let's turn to the second half of the MATLAB page.
Just to give you a sense of how when the eigenvalues of a matrix behave, you can get beautiful behaviour.
So suppose I say a equals ran down of 500.
That's now a 500 by 500 random matrix.
It's not symmetric, positive, definite. But if I say A equals A transpose times A, then it is symmetric, positive, definite.
So I've got a 500 by 500 symmetric, positive, definite matrix.
If I make the usual right hand side. What we'll find here is that the eigenvalues of A are not so well behaved.
So if I say eigenvalues of A and I look at the smallest eigenvalue,
it's very small and the largest eigenvalue is au of one or au of a thousand, actually.
So that's bad. That's the eigenvalues are badly distributed for this matrix.
If I try to run conjugate gradients.
It eventually gets there, but it doesn't do well.
It's very slow convergence. On the other hand, suppose I shift the matrix to make it a better behaved matrix.
Suppose it happens to have eigenvalues that are bigger by 100.
So I've now taken my badly behaved matrix and moved it over by 100 times the identity.
If I define l mean to be the minimal eigenvalue it's now about a hundred and l max is 100 more than what it was.
We now have a good set of eigenvalues which are well separated from zero.
So there exist polynomials that are very small on the eigenvalues, and now the convergence is completely different.
If I say C, g. So let's put more on you know about that in MATLAB.
If you say more on it, just fill the screen. So let's do that again.
So you see, I've got this 500 by 500, but it could be 5000 by 5000 matrix.
And we're converging at a reasonable rate after 20 or 30 steps.
We have a few digits and after 20 or 30 more steps, we have six digits.
So this is working just fine. Okay, I guess it's time.
We'd better stop.