I think. Everybody, let's get going. I've made a few more copies of a sign up sheet.
Please keep them circulating. So anybody who didn't sign last lecture, please put your name down.
Now, a comment about that.
I normally put up a list on the web of people in the course by department because it's interesting to see who's in the course.
If for any reason you don't want your name on that list, send me an email and I will commit your name from the list.
I think, by the way, the filming here probably doesn't include your faces, so you're fairly anonymous anyway.
I don't think we're yet available online, but the hope is that before too long you'll be able to look at lectures online if you want to.
This thing that I've put up here has nothing to do with this course,
but it's just something I've been interested in recently thought might interest you.
So the centrepiece is this graph. I've been involved for many years with writing papers and serving as a referee and an editor and so on.
And I got interested in the phenomenon that papers are getting longer.
In particular, my organisation is Siam, the Society for Industrial and Applied Mathematics,
and I took a look at three of the main journals from Siam, Siam Journal of Applied Mathematics and Numerical Analysis and Scientific Computing.
And this plot shows as a function of time the average length of a paper in these Siam journals.
It's just astonishing how clean the data is. So basically when my career was getting going, a paper was 12 pages long.
And as your career is getting going, a paper is 24 pages long.
It's a doubling of the length of papers.
There are many things going on here, and indeed it's fairly controversial what the essential explanations are.
But it all happened unintentionally. Nobody planned this.
It's just one of these cultural shifts. I think it's happened in other fields too, to some extent.
And actually, if any of you have views on this matter, please send me email.
I'm very interested in the subject and I plan to write some kind of summary of that subject.
I know, for example, it's also happened in economics, in pure mathematics.
The papers are three times as long and they've been growing in proportion.
Okay. So we're talking about conjugate gradients.
That's what the last lecture was about. And today we're going to continue on that subject.
It's such an important subject and talk in particular about Preconditioning, but we didn't really do things properly last time.
It all went rather quickly.
So I want to bring you back to the last lecture and emphasise this picture which captures the essence of the behaviour of conjugate gradients.
Now you remember that conjugate gradient is an iterative algorithm which
constructs a succession of approximations x sub j to an exact solution x star.
And those approximations lie in these so-called cri loft spaces which are getting bigger and bigger.
So at each step you're in a space of dimension one greater, and as a result, your approximations can't get worse.
They have to get better because conjugate gradients magically find the optimal approximation in the kirilloff space.
So as time goes by, it gets better and better. And the question is how fast?
Now we showed last time you can relate that to a polynomial and it's a polynomial in the matrix.
A we're solving x equals B, a is a very large matrix and the polynomial of A is involved.
But that in turn we related to the polynomial applied to the eigenvalues of a and then these
are the crucial pictures the red lines are supposed to suggest the eigenvalues of your matrix.
And here's a case where the eigenvalues, they're all positive, but they come close to zero.
And here's the case the eigenvalues are all positive and they don't come so close to zero.
Now remember that the key fact is how small can a polynomial p be at the eigenvalues
subject to the constraint that it has to take the value one at the origin?
So conjugate gradients is implicitly finding an optimal polynomial in some sense that small on the spectrum but takes the value one at the origin.
Now it's intuitively plausible that if the spectrum comes very close to the origin, that's hard.
So your polynomial is relatively big, whereas if the spectrum is well separated from the origin,
then it's easy and your polynomial is relatively small.
So this is the case where things converge quickly. All of that can be made very precise.
And in particular, the standard thing one says is based on comparing the maximum and the minimum eigenvalue.
So if that's the maximum eigenvalue and that's the minimum eigenvalue.
A is a symmetric positive, definite matrix.
So its eigenvalues are real and positive. We define the condition number of a to be the ratio.
So Kappa of A is the ratio of Lambda Max to Lambda Min.
And that's a measure of how difficult this approximation problem is.
If that ratio is high, then you're in this situation with slow convergence and if the ratio is small, you're in this situation with fast convergence.
Now let me write down the theorem we didn't quite get to in the last lecture on that subject.
So the theorem goes like this. It's really a corollary of the theorem I did write down.
And it says this. If you look at the a norm of the error at the end step of conjugate gradients and divided by the norm of the initial error,
then that ratio is bounded by twice.
Square root of Kappa minus one over square root of Kappa, plus one to the N.
So what's going on there? Is that. We have geometric convergence.
That is to say, as you take more and more steps with each step, you get at least a constant factor improvement in the error.
The error goes down by at least a constant at each step, and the constant depends on the condition number.
If cap is large, this is close to one, so that's not so good.
If Kappa is small, it's well below one and that's very good.
And incidentally, if Kappa is large, then this is approximately two times E.
To the minus two NW over the square root of Kappa.
So you get at worst exponential convergence e to the minus constant n as you take more and more steps.
And the timescale of that convergence is the square root of Kappa.
So loosely speaking, it takes square root of Kappa iterations of conjugate gradient to converge.
More precisely, it takes square root of Kappa iterations to get a constant factor improvement.
And then, depending how far you want to go, you have to act accordingly.
It's a convenient fact that E to the minus two is approximately 0.1.
It's about 0.14. So thanks to the two here, what that means is that.
Basically when you take square root of Kappa steps of conjugate gradient, you get one digit improvement or nearly one digit improvement in your error.
So that's an amazing, beautiful fact about conjugate gradients and it's a starting point of all sorts of things.
And let me illustrate that now coming back to what we did in the last lecture.
So there's no new handout here. Let's see.
So we'll go back to MATLAB and.
There was a code called M three, M three, c g convergence.
And let me run that again. Okay. So just to remind you, it looked like this.
We set a equal to ram down of 500 of 500 by 500 matrix matrix with random normally distributed entries.
And then we took a equal a transpose A, which gives us a 500 by 500 symmetric positive definite matrix.
And then we set B equal to one 500 comma one to give us a right hand side.
And now the point is that the smallest eigenvalue of this matrix, AA, is very close to the orange, and we're in this situation.
Let's see how small it is. If equals I give a and lambda min equals min of a very small.
If Lambda max equals max of E. Rather big.
So the condition number is the ratio of Lambda Max divided by Lambda Min.
And that's a bad number, right? So ten to the fifth indicates we're going to take a thousand steps even to get one digit of improvement.
The square root is about 1000. If I type in MATLAB, kind of a I get the same number.
MATLAB computes the condition number using whatever methods it uses.
So that's a bad situation. The square root of Kappa.
Is 656. So that means that if we take 656 steps, we can expect to reduce the error by about a factor of .14.
That's awful, right? And sure enough, if you run C, G of a B using the code I showed you before, it eventually converges.
But that's actually because the matrix is only 500 by 500.
It's not really converging due to this effect, but due to the finite dimensionality of the matrix.
The initial phase is awful. So they're the first steps.
Not much is happening. So more steps. The more steps you can see, the error is going down very, very slowly.
That's an example of an IL condition matrix, for which conjugate gradients is essentially useless.
Now let me do again what I did last time. I then shifted the matrix.
I said A equals a plus a 100 times the identity of size 500.
So now, of course, the eigenvalues shift accordingly, and the minimum eigenvalue is now 100 bigger than it used to be.
The maximum value is also 100 bigger. So now we have a much better ratio.
The ratio of the biggest to the smallest eigenvalues is only 20.
The square root of 20 is 4.5.
So now we're in this situation, and if we take 4.5 steps of conjugate gradients, we can expect an improvement by a factor of 0.14.
Let's do that. So I'll say CG.
Now you can see one screen full. We're already getting a couple of digits.
Let's take another screen full. So after two screen falls, far fewer than 500 steps, I've already got six or eight digits of convergence.
That's kind of get gradients when things are going well. Okay.
I meant to say all that last lecture, and I didn't quite get around to it.
Now. Let's digress for a moment and play with the handout you did get today, which is called M4 DOT Indices.
This is just sort of exploring MATLAB,
just to remind you of some of the stuff you can do in MATLAB and in other languages too, if they're at a high level.
The point being to encourage you to think at the matrix and vector level rather than the individual entry level.
Of course, ultimately computers are mostly working with individual entries,
but you shouldn't be thinking at that level and you shouldn't be writing triply nested for loops ever.
So let's just play with that handout. If I say one coal and six.
I get a vector of length six. If I say one coal and two coal and six.
It's based by truth. When MATLAB was invented in the late seventies, there existed one or two things that could have similar structure like APL maybe,
but this way of user interaction was really non-standard.
It seems so elementary to us now, but it was a big transition.
A few more examples. So if I say six colon minus one one, I get a vector of decreasing numbers.
Let's say a equal 1 to 16. Let's say a of end minus three colon end.
Oops. A minus three. So end always refers to the last entry in a vector.
Let's suppose I say ants of 1 to 3 ants always refers to the most recent thing that you've computed and not given a name to.
Suppose I say a of two colon five. Colon 16.
Well, two. 712. Suppose I say a equals reshape a little A for for little AA is.
A vector, but if you reshape it, you can turn it into a matrix.
So I've taken that data and strung it out column wise into a four by four matrix.
Suppose I say a of two 516.
Now with a capital AA, well, I get the same numbers as before,
which means that in MATLAB you can index a two dimensional matrix by a single index if you wish to do that.
Suppose I say rep matt. Repeat matrix ance six.
One. Well.
And since this most recent thing we've computed this row vector and rep Matt will take that that matrix or vector and panel it in a six by one array.
So we get that. Now let's take a and change it a little bit.
I'll say A of two colon. Five colon 16 equals.
97, 98, 99.
So there you see we've modified some entries of The Matrix or I suppose I say A of two 516 equals zero.
So the point here is that if you're setting a bunch of entries to a single number, you don't have to make that a vector.
On the right hand side. You can just write the number. Suppose I say a of colon two.
Then you get the second column of the matrix. Similarly a of two comma colon would give you the first row.
Or suppose I say a of colon two equals the empty vector.
Then you see A has now been changed. It has removed the second column.
Suppose I say find a double equals zero?
A logical equals. So it finds those indices where A takes the value zero and notice by default it's doing that in the single index mode.
So the second entry is a zero and also the eighth entry in the matrix is a zero.
Suppose I say a ants equals 15?
Well, it will take these two indices we've just found and change the value of a at those points to 15.
Suppose I say find a equals 15.
There are three of them now. On the other hand, suppose I want to indices, I could say I j equals find a equals 15.
And then you see it's giving us the row and column indices of the three entries of a that are 15.
Or suppose I say a equals 15 equals,
so a logical condition that will give me a matrix of zeros and ones with a one in each position where a equals 15.
Suppose I say a greater than 15? That's another logical condition.
So again, I'll get a matrix of zeros and ones showing me where A is bigger than 50.
Suppose I say find of a sorry find of mod of a two equals zero that will find the entries of a.
That r zero mod to. In other words, even. So there they are.
Here's a again. So those are the indices of the even entries of a.
Suppose I say a of ants equals 100 times a events.
We've now taken the the even entries of the matrix and multiplied them by 100.
Suppose I say who's round of five?
There we have some numbers. Let me make it format short. So there we have a five by five array of random numbers between zero and one.
Suppose I say a of five colon for colon 21 equals one.
Colon five. Squared. Now what I've done is I've taken the cross diagonal entries and I've changed them from random numbers into 149 1625.
Now suppose I say a of a not equal round of a equals 999.
Round of a of course rounds a number to an integer. So the places where A is not equal to its rounded value are the places where it's not an integer.
So you see, I guess turned all the random numbers into 999.
This is the way in high level languages you can do a lot of stuff and it gets convenient always to try.
As a matter of principle, to do things the slick way. Okay.
That was the end of the digression. I know there's a lot of variety in this room.
Some of you have used MATLAB for years, others not so long.
Okay. So back to our main story, Preconditioning.
The lesson here is that if you have a matrix problem whose spectrum is well separated from the origin, then conjugate gradients will converge quickly.
So that suggests given a matrix problem X equals B, we should add 100 times the identity to it.
The trouble is, you can't do that. If X equals B, you can't change AA and still have an equivalent problem.
There's no way automatically to turn that picture into that picture.
If there were, it would be a very different world.
So adding a multiple of the identity is not a successful idea in general, but multiplying is preconditioning is a multiplicative idea.
Which sort of came along in the 1970s and their variations on this.
And the one I'll talk about is what you could call left preconditioning.
Suppose you have a x equals b is the problem you want to solve.
Well, that's of course equivalent to m inverse.
A X equals m inverse B for any matrix m that's non singular.
And the whole idea of preconditioning is defined matrices m.
Which have two properties. One is that you can effectively compute the inverse or at least solve matrix problems involving them reasonably fast.
And the other is that this improves the picture to a better condition matrix problem.
The surprising thing is that in almost all applications there are ways to do this.
So there's a huge business of finding good preconditions, and that's more or less what I'm talking about today.
So let me say more. Now, of course, if A is symmetric, positive, definite, which it had better be for conjugate gradients, then.
If we take. M and factor it.
In that form. That's called a collective factorisation.
So if we have a symmetric positive, definite pre conditioner and we factor it into a product of one matrix,
transpose times that matrix, then another way to write the equation turns out to be c inverse.
A C minus transpose.
C transpose x equals C inverse.
B, I'm not going to go through that and indeed, I'm not going to use this formulation.
But the point is that in practice, symmetry can often be preserved.
So if you have a symmetric problem and you want to precondition it, there are ways to do that, to preserve the symmetry and.
That's just mechanics. I'm not going to talk about it. Symmetry can be.
Preserved when you precondition. But I'm just going to talk about the fundamental idea using this basic mathematical idea.
So we want to improve the conditioning of a matrix by pre multiplication by something.
Let's take a look just for fun at the paper. That's among the handouts by my shrink and then divorced.
So this is truly one of the greatest hits of numerical analysis.
There's a lot of controversy about who invented Preconditioning.
And without a doubt, these were not the first people to have ideas in this direction.
On the other hand, they were the ones who put it on the map. This paper really did change computational science slowly.
It was an agonising story for the authors. I think they wrote the paper in like 72 or something.
Had a lot of trouble with referees. There was one particularly famous referee whose name I will not utter who tried to squelch the paper.
Eventually it got published years later, and slowly it took off.
So starting in the late seventies,
everybody was talking about pre-condition conjugate gradients and I'll say a word in a moment about their preconditions.
Which is just one of many methods, but a particularly good one called incomplete factorisation.
So Vandervoort in particular has been one of the great figures in iterative linear algebra for many years.
Henk Van divorced from the University of Utrecht in the Netherlands.
So you've now seen two papers that changed the world.
One was destinies and Stiefel on gradients, and then this on preconditioning conjugate gradients 25 years later.
Okay. So what I want to do is illustrate a bit of the effect, a precondition or can have.
So we're now going to look at the codes called PRC and M5.
So let me remind you that MATLAB has built into it a precondition conjugate gradient code called Pcgg.
And then in addition, I've written a few elementary ones which are not software,
they're just illustrations of how simple the preconditioning idea can be.
So you look at my codes in order to get sort of the essential mathematical idea, but you would never use them as software.
So c g is un precondition and Prague is precondition with a precondition or m just has written there.
And let's play with that. So let's say.
Suppose a is the diagonal matrix from 1 to 100.
And let's suppose B is the usual trivial one right hand side.
And let's try CG of a. So that's an example of not so good convergence.
It's going but not so great. Let's take AI to be the identity of size 100.
And now let's try various preconditions. So first of all, suppose I try precondition conjugate gradients and my m matrix is the identity.
So I precondition by replacing that by this where m is the identity, of course it will have no effect.
So that's one extreme you get exactly what have I done?
Oops. The same numbers as before. Here's a plot of the convergence.
Eventually it gets there, but it's not especially great.
Let's go to the other extreme. So suppose I say p.r c g of a b.
Hey. So now I'm going to precondition this problem by setting em equal to A so it becomes x equals adverse B, it will converge in one step.
Oops. What am I doing wrong here? Okay, it converge in one step.
The. The X axis should really say 0 to 1 instead of 1 to 2.
Basically in one step it converts. So that's the perfect precondition when M equals A itself.
But of course, the whole point is that we assume A is hard to invert.
So although it's a perfect precondition, it's a unusable precondition.
The whole point is to find a compromise sum matrix m which can be inverted quickly, which improves the conditioning.
So in this simple illustration, let's just try as the square root of the diagonal.
So let's say. Let's say d equals diag of square root of 1 to 100.
Now our matrix happened to be diagonal. So if I took the equals the diagonal of the matrix, it would be a itself.
That's why I've put in the square root just to mess things up a little bit.
So here's an example of a matrix that captures some of the flavour of a but not all of it, and it's a pretty good precondition.
So if I say. PR c.
G of a. B. D.
This is more typical of the kind of behaviour you strive for in using conjugate gradients.
You see, we're having good smooth convergence at a good geometric rate.
In 30 steps we've got ten digits of accuracy. This matrix was only 100 by 100, so that's not too exciting.
But it might have been 10,000 by 10,000. Let's do a bigger problem.
Oh, I'm sorry. I forgot that I was going to do the same thing using the built in MATLAB precondition conjugate gradient.
So this is now the built in MATLAB code. If I say PC G, a, b, and then I give it a tolerance and I tell it how many steps it can take.
Then eventually it converges in an automated way.
So let's now do a bigger problem. I'll say an equal 30,000 and I'll allocate some memory sparse allocation of memory Ballack of end to end for end.
So we've just allocated we've told MATLAB that we're going to be working with an end by end sparse matrix and give it for and.
Double precision memory locations and I'm going to set I to be an index set going from end 30,000 in steps of minus eight up to N squared.
And I'm going to set a of I equal to one. So you see I'm now doing the kind of fancy stuff with indices that I showed you a moment ago.
Fancy fooling around. We won't know what that looks like until we draw a picture.
But first let me say a equals a plus a transpose plus diag of sparse of one to n.
So what I'm doing here is just fooling around and constructing an example matrix.
It's 30,000 by 30,000. Unless you're very quick, you probably can't visualise its sparsity pattern on the diagonal.
It has one to n and on the off diagonal it has the number one appearing in various positions.
So it's a sparse matrix and it's symmetric. Here's what it looks like.
To me that looks like some kind of a cake that people have decorated with frosting or something.
So on the diagonal it goes from 1 to 30000 and all the off diagonal blue dots are the numbers one.
So a typical sparse matrix with and it's not at all clear what kind of structure it's got, whether it's a usable structure or what.
Now let's solve a problem with that.
So I'll say B equals ones of 30,000 comma one and I'll say c, g of a B, so we're now trying unfree preconditioned conjugate gradients.
The first screen fall immediately tells you that it's not doing very well.
We've got about a factor four in one screen full. It's converging slowly, so without a precondition or conjugate gradients is not satisfactory.
Now let's precondition and here's what I'll do. I'm going to take a diagonal precondition.
And this is a very standard, important idea. In fact, this is almost the basic first precondition.
Or you should always try for every problem. A is some complicated matrix with non zeros all over the place.
That constructs a diagonal matrix with the diagonal entries of a and surprisingly often that's a pretty good precondition.
So let's try it. Let's say X equals prc g of a B with the diagonal precondition.
So look at that. It's beautiful behaviour for this matrix.
In seven steps we got 12 digits of accuracy so that's better than usual.
If you achieve this you are in very good shape, spectacularly effective preconditioning just by using the diagonal of the matrix.
Even this idea wasn't really around until the 1970s.
There are precursors, of course, but nobody was really taking advantage of this kind of thing until the seventies.
Of course one can also do it using the built in software. So suppose I say pc gab1e.
Or am I. If you press the wrong thing, it can be annoying.
One E, minus ten. 100 D So using the built in matlab pc g code you see in five steps it's happy.
So that looks like a great example of a matrix where iterative methods are spectacularly good.
We could even do a tick tock on that. Suppose I said talk.
Tick. So you see in less so about a 10th of a second, we get this 12 digit accuracy.
Let's try for contrast doing a backslash b.
So this is Mount Labs, the standard built in X equals B solver backslash.
And exactly what that does is an interesting story that we haven't talked about.
So. You might expect it to take a staggeringly long amount of time.
Remember the first lecture we did a 4000 by 4000 matrix and it took several seconds.
This is a 30,000 by 30,000 matrix. So who knows?
It may take minutes or not. Let's put a tick and talk around it to see what it takes.
Well, look at that. That's also spectacularly fast.
This, it turns out, is using a completely different method. It's not iterative.
It's using a method that cleverly reorders the rows in the columns of the matrix to bring those ones closer to the diagonal.
So this is a surprise. When I tried it, it surprised me.
I thought I had cooked up an example in which you really needed an iterative method.
But this experiment shows that for this example, you don't need an iterative method.
And in general, you should never assume a problem needs an iterative method.
It's always worth trying these more elementary so-called direct methods, which in MATLAB are encoded with backslash.
A rule of thumb is that often problems that come from a two dimensional spatial domain don't need iterative methods.
Often, backslash type algorithms work fine. If you come from a three dimensional spatial domain more often, you do need iterative methods.
Let's compare the accuracy. Suppose I say max of abs of x minus x two.
So I'm now comparing the result of preconditioned conjugate gradients to the result of backslash,
which is a direct non-intuitive method and I find perfectly good agreement.
Okay. So what can one do in a course like this with preconditions?
I can tell you they're incredibly important. It's kind of embarrassing how important they are.
A computational science is not in principle about linear algebra, and yet in practice it often is.
And linear algebra is not in principle about preconditioning. And yet in practice it often is.
It's just remarkable how much of computational science often comes down to having a good precondition.
And what I want to do in the rest of this lecture is just sort of talk about
14 classes of preconditions that have proved important in the last 40 years.
Just to give you a sense of some of the breadth of this. Now, maybe I can get away with even shutting that down.
If I shut down the computer, will that cause trouble? Yes.
I won't get away with that then. So. Okay, we'll stick with this screen.
So this is what I'm going to call Section 1.7 Whistle Stop Survey.
Of preconditions. And before we do our survey, let me mention that Oxford, as it happens,
has an expert in the field, and that's Andy Watson in the Numerical Analysis Group upstairs.
So he is in the field of numerical linear algebra and in particular, he's an expert in Preconditioning.
So if you are trying to precondition your problem, this is really one of the world's most knowledgeable people right here in this building.
Okay. Oh, another remark to make is. For reading, you might try.
Chapter 40 of Profession in Beau is more or less what I'm just about to mention to you, minus a few things.
That's a discussion of a number of preconditions. Okay, so our whistle stop tour, we have about a minute per precondition.
The first one dimension is diagonal scaling.
And there's another word that's often used for that.
It's called Jacoby because the great mathematician Jacoby used ideas like this.
He was amazingly modern in the middle of the 19th century. He was doing things which still have numerical application today.
So this is absolutely the first thing to try for almost any problem.
It's often good. At the other extreme in complexity is the idea that my shrink and that divorced introduced namely incomplete factorisation.
Now, let me just say a word about what that means. So we haven't talked about Gaussian elimination because it's too familiar, basically.
But the standard method by which people solve X equals B is to factor A into a product of a lower triangular and an upper triangular matrix.
That's Gaussian elimination. Now, if you can do that, that's great.
The trouble is that even if A is sparse, the factors L and you are typically not so sparse.
So you destroy a lot of the potential of your problem if you can factor it.
Moreover, if you aren't so sparse, a lot of work is going to be involved.
So the idea that my shrink and Vandervoort had was to do a factor ization in which you force a lot of the entries to be zero,
even though mathematically they shouldn't be you instead of a equals L2, which would be too expensive and not sparse enough.
They do some kind of an A equals an approximate L and an approximate you wear and it's not equal.
It's approximately equal where L and U tilde are forced to have a sparsity pattern.
That turns out to be a very powerful technique. It's much more complicated.
Of course, you have to decide exactly what your strategy will be, and the analysis of it also turns out to be not at all trivial.
So it took many years for people to be comfortable with exactly why it worked as well as it does.
But it really is a powerful technique. And in MATLAB.
There are some built in capabilities for this, in particular.
C l u i n c which is a non symmetric incomplete factorisation code and then ch0lincchol
stands for LSC and that refers to the symmetric case of incomplete factorisation.
The third one I'd like to mention is the what I call course grid approximation.
So in computational science, very often we have a big matrix problem because there's some grid that we've constructed in two or three dimensions.
And of course we want the grid to be as fine as possible if it's 100 by 100 by 100,
and there's just a scalar variable, that's a million unknowns already.
If there are three variables at each point, that's 3 million unknowns.
If it's 200 by 200 by 200, that's 24 million unknowns.
So you very quickly get big matrix problems. However, as you might imagine.
If you have a very coarse grid that might give you some kind of crude approximation to the solution on the grid, you ultimately really care about it.
So just to draw the simplest picture, the four by four problem.
Might be some kind of approximation to the HIPAA problem and so on.
That might be some kind of an approximation to the 16 by 16 problem.
So this might be a precondition for that.
At each step, we might somehow solve a problem involving this relatively easy matrix.
And if we do things right, use that approximation as our precondition.
That's a great idea in many cases. And indeed, that's the idea which eventually led to what's called multi grid methods.
So if you do this just once. Precondition a fine grid by a coarser grid.
Well, that's a precondition. But if you recurs on that idea going down a level to coarser and coarser grids, that's called multi grid iteration.
And that's another one of the greatest hits ideas of the last few decades multi grid
iterations for solving the linear algebra problems that arise in computational science.
One more thing to say about that. If you think about what's involved, it's intuitively clear that this is has a hope of getting,
as it were, the long range structure, the the big wavelengths.
Obviously, you're not going to resolve the small wavelengths on the course grid.
So this preconditioning idea is one of preconditioning the problem by getting the large scale structure somehow.
Now the next idea is almost a dual of that. I could call it local approximation.
So any scientific problem with a spatial aspect to it?
Probably has some physics that's essentially local and some other things going on that are reasonably global.
The boundaries may be affecting you in fundamental ways,
but also just whatever locally is diffusing or convecting or whatever will of course, be important.
And another idea for Preconditioning is that if you can do everything sort of locally,
correct, ignoring the long range that may correspond to a close to diagonal matrix.
So this may give you a precondition or which is close to diagonal and therefore implementable.
So in some sense, if you omit long range interactions, that may be an idea for a precondition.
Let me mention then number as number five, what I'll call block precondition is.
And more or less the same idea. Domain decomposition.
This comes up in various context. Physically, you can imagine a system, let's say you have an aeroplane, it's got a left wing,
it's got a right wing, it's got the body, it's got the engines and so on.
Maybe you can solve the body relatively easy, but then coupling it in with the others leads to complexity.
So the idea here is to use as your precondition for these various sub systems of the complicated system.
Another way to say that, though, is you've got a giant matrix and maybe the first few columns correspond to the body of the aeroplane,
and the next few columns sort of correspond to the left wing and the next few columns correspond to the right wing and so on and so on.
You can imagine they're all sorts of connections between these different systems.
But if we approximate by pretending that the body is only connected with the body and the left wing is only connected with the left wing,
and the right wing is only connected with the right wing and so on. Then we're getting some kind of a block diagonal structure.
We're omitting crucial parts of the physics that are out here, but hopefully of a lower importance.
So that's another generalisation, if you like, of the idea of a diagonal precondition.
Are these all of these ideas are just unbelievably big in practice.
For example, there's a domain decomposition conference that has happened every two years for the last 25 years.
You know, little bits of this have whole research communities and big scientific impact.
The next one I want to mention is low order discontinuation.
Now what I mean by that, if you take a differential equation and approximate it by finite differences,
if you just look for second order accuracy, you probably get a reasonably sparse matrix.
But bad accuracy. If you look for fourth or sixth order accuracy, you get a worse matrix, but better accuracy.
So a natural idea is to formulate your dissertation at a reasonably high order to get accuracy on,
not to find a grid, but then use the sparser low order matrix as your precondition.
And if you take that idea to the limit, there are things called spectral methods which are sort of of infinite order.
And some people would approximate a spectral method by a finite element method or something.
So you can have an infinite order method preconditioned by a second order method, for example.
Okay. Next thing, let's mention the idea of taking a problem and a constant coefficient approximation to it.
And I'll put in the same category or a symmetric approximation.
So in many cases, if your problem had constant coefficients or were symmetric, there might be special fast ways to solve it.
For example, a the Poisson equation with constant coefficients.
There are things called fast Poisson solvers which take advantage of the fast Fourier transform to get spectacular speed.
If only my problem had constant coefficients, one would say I could solve it very quickly.
Well, whenever you're in an if only situation, then maybe you've got a good precondition.
So use the constant coefficient thing to precondition.
Or if your problem is dominated by some symmetric physics that can be encoded in a symmetric matrix, and yet it also has some non symmetric physics.
Maybe the symmetric part can be solved more quickly and can be used as a precondition.
The next one is splitting. Related to that.
So I guess the term we like these days would be multi physics.
Suppose you have a problem, as people so often do, which has more than one bit of physics going on.
And the classic example would be in fluid mechanics, you have convection and you also have diffusion viscosity.
So these are different bits of physics and they turn into different bits of mathematics.
They're discrete times by different linear algebra problems that all combine into one.
It may be that if you only had one of the problems or the other, you could solve it quickly.
Maybe if you only had the diffusion, you could use a fast Poisson solver.
But the convection makes that difficult. Well, that's a very standard way of devising a precondition.
You've set up the big problem. Hive off part of it that you can solve as your precondition.
Next idea is dimensional splitting, and related to that is the idea of Add II, which stands for alternating direction implicit methods.
Here, as you can guess from the name. You have a problem in 2D or 3D.
If only it had one dimension, then it would be a diagonal matrix.
It's the multi dimensionality that makes it non diagonal.
So the preconditioning idea is to only look at the interactions in X, for example, as your precondition to solve the X part of the problem.
Exactly. That's not the whole problem. But maybe it's a good precondition.
Another idea is to take one step as your precondition of a classical iteration.
Now, by a classical iteration, I mean the things that go by the names of Gauss Seidel.
Or successive over relaxation or symmetric successive over relaxation.
These are things that were big in the 1950s before people realised how powerful conjugate gradients was.
They're much less important now. That's an L. It's much rarer these days for people to solve a problem just by using successive over relaxation.
However, these things make great preconditions for higher order problems,
and this has to do again with separating the local effects from the global effects.
These can sort of take your problem and globalise it in a in a good way.
11. Is the idea of a periodic approximation to a non periodic problem.
This was invented by Gil Strang at MIT.
He was actually visiting Oxford this term. I'm curious who have seen lectures by Gil Strang on the web?
Okay. One or two? Yeah. He has these linear algebra and other things that have been very popular at MIT.
If you Google him, you'll find that. So it's trying.
Years ago proposed this idea that take a non periodic problem like a toeplitz matrix which has constant long diagonals.
Imagine the Associated Periodic version of the problem, which in that example would be a circular matrix.
It may be you have special methods for the periodic problem, like a fast Fourier transform, so you can solve that one very quickly.
It turns out in some cases it's a spectacularly good precondition.
We're a little short on time, but let me quickly at least name the final three.
One idea is to use an unstable method. We haven't talked about stability yet, but some numerical methods are fast but unstable.
They amplify rounding errors. The classic example would be Gaussian elimination without pivoting.
If you don't interchange rows, you get an unstable method of solving X equals B.
Probably not good by itself, but it might make a great precondition.
And Bob Skeel is a famous name in that area. Another idea is polynomial preconditions.
And the idea here. Is to approximate a inverse depth,
to approximate a inverse somehow by a polynomial of a and that has the advantage of being a beautifully parallelism idea.
So it's an idea that it's somehow it's never the best precondition or but it's always the most paralysing parallelism of one.
So the idea keeps having importance.
And actually one of the guys in the Numerical Analysis Group is working with things along these lines to solve eigenvalue problems.
Jarrett aren't working with use of on that. And finally, let me mention the idea of sparse approximate inverses.
Well, the idea is if we do the inverse. Exactly, that would be great. At least if it were sparse enough.
In this business, people solved least squares problems to explicitly construct a matrix, which is a pretty good inverse.
Not exact, but close to an inverse and is also quite sparse.
So least squares technology turns into making a precondition.
So I think this list probably covers more than half of the big ideas out there.
It's just a huge, huge subject. And that's pretty much all we're going to say about it except for the next assignment.
Okay. Thanks.