Okay. Good morning, everybody. So assignment two is due now.
I hope most of you manage to solve the big million by million problem.
I found that the with a diagonal precondition artist of 20 seconds on my lot on my desktop, which kind of amazed me.
The eigenvalues took longer, a couple of minutes. But anyway, the Million Buy Million worked for me.
Is there anyone who had to go to a smaller dimension to make it work on your computer?
Yeah. Okay. Sorry about that. Next year you could do 2 million.
Okay, so today we're going to talk about floating point arithmetic, which is an amazingly fundamental subject and actually pretty interesting.
It's, in my view, not quite the heart of numerical analysis, but it's one of the foundations.
Well, the foundation. And it's pretty interesting. Let me mention three books, which I'll pass around.
And they're also listed on the website, very different books, each in their own way.
Fantastic. So this very small book by Michael Overton is called Numerical Computing with Tripoli Floating Point Arithmetic.
And it's just an incredibly clear, straightforward way to find out all sorts of things.
You had no idea about how our computers work with numbers, so pass that around very much.
At the other end of the complexity.
Thank to is this thick book by Nick Higham called Accuracy and Stability of Numerical Algorithms.
And this is the book about careful analysis of all sorts of fundamental numerical algorithms.
And in particular, what effect do the rounding errors have? How does floating point arithmetic affect the accuracy of the final result?
Hi is at Manchester and he's the world leader on that subject.
And then finally, this is another thick book called Handbook of Floating Point Arithmetic by Mueller et al.
This is a group in Lyon, in France, which is a remarkable, very interesting group.
They work with Intel and other people designing the details of the elementary functions, special functions that make computers so accurate.
So if you're interested in the details of how sine of X or something might be computed, this is a great starting point.
Okay. So to talk about floating point arithmetic, I want to begin by thinking in terms of physics and one of the handouts,
the one with a colour picture on it, discrete or continuous, conveys a bit of this flavour.
So I want us to think about the physical world that we live in, like this room which contains a lot of solid matter and a lot of gas in the middle.
And of course, you and I know that this is all about molecules that bounce around and do their thing.
That's a very old idea in some sense. But it's really in the last couple of hundred years that it became definite.
I guess the atom was pretty much discovered in Manchester by John Dalton.
And statistical mechanics came around in the 19th century, Maxwell and Boltzmann and others.
So these are absolutely fundamental ideas in physics.
Everything is made of molecules or atoms.
The point of view I want you to take is that molecules are, as it were, the implementation of the classic laws of continuous physics.
Just as floating point arithmetic is the implementation of the laws of real mathematics on a computer.
It's a very interesting analogy.
In both cases, you have something discrete that is necessary at the hardware level in order to achieve an effect which feels effectively continuous.
Now, let me give you a couple of examples of how physics implements the laws of physics.
So, for example, suppose you have a box filled with gas.
So that's a box filled with some air in it, if you like.
And what happens if you have the volume of the box and suppose you keep the temperature constant so you make it the same temperature?
Well, you know that the pressure doubles.
In fact, that was discovered in Oxford, I suppose, in the 1660s.
So the gas laws tell you that if you have the volume at constant temperature, the pressure doubles.
Now that's a continuous statement, but there's an implementation underlying that because there are a lot of molecules in that box.
And the point is, if you have the volume, then the molecules hit the boundary twice as often.
So they're twice as dense. So you have molecules hitting the boundary twice as often.
That's the implementation. That's how this law works, because twice as many molecules per unit time per unit area hit the wall.
So the amount of momentum transferred is doubled, the pressure doubles basically.
So that's implementation of one of the gas laws. Let's take another example.
Suppose you have the same box. But now let's double the temperature instead of having the volume.
Now, you know, another gas law says that again, if you do that, the pressure doubles.
Why does it do that? Well, again, there's an implementation at the hardware level that makes that happen.
The implementation now is actually a little more interesting. This one was kind of trivial, but this is more interesting.
Again, you have these. Molecules in there.
But if you double the temperature, each one has twice as much energy, which means the square root of two times as much momentum.
So the velocities increase by the square root of two.
Because the energy of each particle is doubling. But at the same time, if the velocities are increasing by the square root of two,
so that means that the square root of two times more momentum is imparted to the boundary with each impact.
But there are also the square root of two times as many impacts.
So more square root of two times more impact squared or two times more momentum per impact.
That means, again, that the pressure doubles.
So this kind of thing was among the great discoveries of the 19th century, figuring out how our continuous world is based on discrete hardware.
Now we all know about Avogadro's number. Let me remind you of some numbers.
So we all know ten to the 23, and those of us who are chemists know exactly what that means,
and everybody else probably forgets exactly what it means. It's the number of things in a mole.
And then the question is, what's a mole?
And if you're a typical mathematician like me, you have to think a little to remind yourself of what a mole is.
Anyway, ten to the 23, roughly speaking, is the number of molecules in a handful of matter.
And that applies. It's precisely a mole, but a mole is something on a human scale.
And let me tell you what that translates into under normal conditions of pressure and density in the atmosphere.
Avogadro's number. Under normal conditions turns into about.
Three times ten to the 25th.
Molecules in a cubic metre of gas.
So in a cubic metre of gas in this room, that's how many molecules there are.
And that's called Low. Schmitz Constant. I didn't know that number.
But anyway, there it is. So there are some real scientists in the room.
I'm curious, has anyone heard of Low? Schmitz considers that standard. Okay.
Anyway, of course, it follows from Avogadro's number with under standard conditions.
Now that's a volume.
So in a certain volume of air, you have that many molecules, which means along a length you have the cube root of that many molecules.
So the cube root of three times ten to the 25th is approximately three times ten to the eighth.
So that's how many molecules there are per metre in a linear sense in air.
Typically. And solids are about a thousand times denser than gases, which means linearly.
They're about ten times denser. So three times ten to the eighth would be the number of molecules in a metre of air.
Three times ten to the ninth, roughly, is the number of molecules in a metre of solid.
So let's average three times ten to the eighth and three times ten to the ninth.
And what we find is you get about ten to the ninth molecules or atoms.
I'm being very vague here, of course, per metre in in gas or solid.
That's the order of magnitude we're talking about under normal conditions.
So if this podium is a metre long, roughly speaking from this end to this end is ten to the ninth atoms,
or if my glass of water is a 10th of a metre long, it would have ten to the eighth atoms in a row.
So that's the resolution of continuum mechanics.
To get from the microscopic to the human scale is about a scaling of ten to the ninth.
Now computers. Approximate the real line by a discrete thing, and the resolution is much finer.
It's ten to the 16th, so computers, depending exactly how you count it,
are six or seven orders of magnitude finer resolution of the real line than the human scale.
Physical world is a resolution of, you know, by atoms and molecules.
I find that interesting. So you shouldn't think of computers as very discrete.
You should think of the real world is very discrete and computers are practically a continuum.
Or another way to put it is it's something like ten to the seventh metres from here to Australia.
So if you took went to Australia through the earth and counted the molecules as you went, then you'd get about ten to the 16th,
which would be comparable to the number of floating point numbers between one and two on your computer.
So here to Australia is like 0 to 1 on a computer.
Let me then say a bit about how computers represent numbers.
So this is a beautiful aspect of computing.
Since about the 1980s, in the olden days, computers were all different.
Every manufacturer would make a machine and it would do its own computer arithmetic.
And they were in various bases. So pre, let's say pre 1980s.
You. Most computers would have been probably in basis two or 16, but you would also find bases eight, ten and even three.
There's a famous example of a Soviet computer that was made in base three as an experiment.
People like to laugh at that as showing how silly the Soviets were, but actually it was probably a pretty interesting experiment.
But anyway, certainly different powers of two and also Base ten have appeared over the years in various computers.
In the 1980s, it was two things were agreed upon.
One was that some of the arithmetic involved here was a mess.
Another was that the mere diversity between all these different platforms caused great confusion.
So for both of those reasons, people decided to make a standard of floating point arithmetic, and that was in practice came in around 1985.
The Tripoli Floating Point Arithmetic Standard.
And there was an individual who was the key in making that happen.
William Tahan, who is a professor at Berkeley. He twisted arms all around the world and persuaded the world to adopt this standard.
And ever since then, most computers have more or less adhered to the floating point standard from my Tripoli.
It's not perfect. There are always some deviants and who knows what's happening in the future for all sorts of reasons.
But so far, for the last 30 years, there's been a tremendous amount of consistency, completely different from the previous 30 years.
So I want to tell you how I Tripoli represents our standard double precision number on your computer.
Some of you know this, but probably most of you have a sort of vague idea of how it works.
So. Let's draw one.
And to. So that's the number one and that's the number two on your computer.
And that represented exactly on your computer. And between one and two, the spacing is exactly uniform.
And it's two to the -52.
Exactly. So those are the numbers on your computer between one and two.
There are two to the 52 of them or two to the 52 plus one if you count both endpoints.
And they're exactly distributed by two to a distance of two to the -52.
And then the beautiful thing is that that pattern is repeated many, many times on a larger and smaller scale.
So if that's three and that's four.
What you find is that between two and four you again have an exactly equal spacing,
but now that spacing would be two to the -51 and between four and eight it would be equal spacing two to the -50.
And then between one and a half, it's equal spacing.
Two to the minus. 53. So that is the system of floating point numbers.
Of course, there are also the negatives of those. And most of the time we act as if this goes on forever.
So we imagine that between a quarter and a half it's the same.
And between an eighth and a quarter it's the same, and so on and so on. I did have an item down to zero.
Of course that's not really true. It has to stop somewhere at both ends.
Similarly, at the big end, 2244288 216.
They all have the same number of numbers between them and it goes a long way, but it doesn't go forever.
So let me tell you a little more about what it does do. A floating point number.
And I'm always referring here to what the double precision standard.
Consist of two words. Which is eight bytes, which is 64 bits.
And the conventional way that we talk about how those bits are used is that we say one bit is used for the sign of the number plus or -11 bits.
I used for the exponent. Which is interpreted as a binary exponent.
And those 11 bits include the sign on the exponent. And then 52 bits.
For the fraction. Which again, is interpreted as a binary base to.
Fraction. So this 52 matches, this 52.
These bits are what tell you which of those you're in. The 11 bits are what?
Tell you whether you're between one and two or two and four or four and eight and so on.
And then the one bit for the sign, of course, tells you if you're positive or negative and if you want to learn this stuff over ten.
This book is absolutely beautiful for that. In a moment, we're going to do experiments on the computer.
And you'll see, of course, that if if there are only 11 bits for the exponent and one of them is a bit,
that means the maximum exponent is two to the 10th. So the this keeps going up to a scale of two to the 10th.
So it's two to the 1024 basically.
So this scale invariant persists up to two to the 1024, which is about ten to the 308.
And then it similarly persists at the low end down to about two to the -1 to 4, which is ten to the -308.
So these are two numbers. Two, it's very important to be aware of the distinction between 16 and 308.
16 is the number of digits.
A number is represented. Two has nothing to do with the scale, it's just the relative precision.
Whereas 308, roughly speaking, is the maximum exponent.
So you can have a number of size ten to the zero, ten to the 110 to the 200.
No problem. Whatever size it is, it will always be represented to about 16 digits of relative precision.
Why 16? Because ten to the 16th is approximately two to the 52.
So always bear in mind that your computer represents numbers to about 16 digits or exactly 52 bits.
There are wonderful examples of computers that failed to do this.
Cohan's bete noire was always the Cray machines, which would do certain subtractions with and lose all the information in effect.
And then there was another example in my career was the great Pentium bug scandal of 1994.
I think it was. None of you were born, hardly. But in 1994, Intel produced with great fanfare, its Pentium chip.
And there was a bug which took a long time to discover because the bug did not make the arithmetic terrible.
It merely meant that instead of 16 digits, it had five digits of accuracy in certain cases, usually 16, but occasionally five.
And eventually a number theorist who was doing very delicate number theory computations discovered that it was inaccurate,
thus proving that number theory is useful. Intel lost a huge amount of money on this.
They had to, of course, do all sorts of things to recall and get their reputation back.
But it changed the world in all sorts of ways.
The whole field of hardware verification, which is an important area of computer science, was more or less created by this Pentium bug in 1994.
So when things go wrong to from other points of view, that can have positive effects.
Now, the way that chips are designed is much more rigorous than before 1994.
Okay. So I want to show you on the computer. Some stuff from the handout called M12.
I'm sorry. I'm not quite ready for that. My apologies. I want to say a little more on the board.
What I've been telling you about is what are floating point numbers.
And that is an answer to the question of how do we represent real numbers on a computer.
What I've forgotten to tell you before turning to the demo is, is how do we compute with real numbers?
So what we've just talked about is representing.
Numbers. Now, let's talk about computing with these representations.
Now, the problem is that if you take one floating point number, let's say,
and divide it by another floating point number, the answer will not be a floating point number, typically.
So you have to do something. And of course, what we do is around at every step when you do something on your computer,
if necessary, the result is rounded to a floating point number.
And this is where the I Tripoli standard is very firm and has improved the world greatly.
Here's the principle. This is the principle of essentially all computer arithmetic.
Since the 1980s. Let's suppose that X and Y are floating point numbers.
And let's suppose we compute one of the basic operations, and people generally regard those as plus minus times a divide.
So suppose we compute X plus Y or x minus Y or x times y or x divided by Y.
The principle is that on your computer, the computed result.
Is the mathematically exact result. Correctly rounded to your floating point system.
Now we have to say a word about what that means.
Of course, rounding means you take a number like one third and you turn it into point three, three, three, three, actually binary, not decimal.
But that's the idea. Now, more precisely.
Well, the idea is that you take a number and you find the nearest floating point number.
Now, what if there's a tie? How do you break the tie? There are very precise rules about this.
And moreover, in my Tripoli arithmetic, there are different rounding modes you can specify.
You can round up, you can round down, you can do various things. So all of that is spelled out very fully.
So correctly rounded is a slightly adjustable parameter, but a fully prescribed adjustable parameter.
And according to this principle, there is no ambiguity whatsoever about the result.
When you do these fundamental operations on the computer, your computer is required to produce exactly the number prescribed by this principle.
And this is pretty much true in most hardware implementations of floating point arithmetic these days.
There might be a few arguments at the end. At the end, you know, down ten to the -308 and things like that.
But more or less this is true. There's a corollary of this principle.
So supposing this principle is upheld. It follows that each of these operations gives a result we can talk about mathematically.
So I'll put it this way. I'll say each operation yields. Computed.
X operation y.
Equals the exact value of X operation.
Y. Times one plus some number epsilon whose absolute value.
Is less than what we call machine epsilon. Which is two to the -53 in the standard I Tripoli.
And that's about 1.1 times ten to the -60.
Now, let me say a word about what I mean here.
I'm assuming the default rounding mode which is round to the nearest with ties broken in a systematic fashion.
And that's where two to the -53 comes from. If the gap between two numbers is two to the -52, then the nearest one is two to the -53 away.
So whenever you round a real number to a floating point number, you make a relative error of at most two to the -53.
So that's what's that being said here. Every computation involves a relative error of scale, at most two to the -53.
Now I'm ignoring in all of this discussion what happens if you go over the top or down below the bottom?
I'm assuming no overflow, no under flow.
And most of the time things don't overflow or under flow.
So no, none of our numbers get bigger than ten to the 308 or smaller than ten to the -308 if they do,
of course, that's prescribed also what's supposed to happen. So this is a very firm principle.
It's really specifying exactly what every operation has to deliver.
It is possible to implement it.
And this is successfully done and it yields this very fruitful corollary on which you can do all kinds of mathematical reasoning.
So just to emphasise the point is that. Computers don't make sort of approximately a rounding error at each step of roughly this size.
They rigorously precisely make a rounding error of strictly no greater than exactly this size.
Now, let's go to the experiments. So.
I'm going to take us through the first half of this M12 handout.
The second half uses the symbolic toolbox, and my laptop doesn't have the symbolic toolbox, so I won't be able to do that.
But if you have access to that, I give it a try because symbolic stuff is always fun.
But of course, it's not what computational science normally uses.
I want to do some experiments which. Are probably things you've never tried.
And what I hope to reassure you with is how deterministic and understandable floating point arithmetic really is.
So let's begin by computing one over 4096, which of course is two to the -12.
You get a number printed out that is just a digital display.
It's not the exact number. The exact number is 64 bits in the split up in the way I've described.
So there's a digital approximation to an exact number.
Now you know that one over 4096 is going to be represented exactly on the computer
because it's precisely a power of two and it's nowhere near the under flow limit.
So whatever that is. It will be exact.
Suppose I multiply it by 4096?
You know what the answer has to be? It has to be exactly one, because we know we've got exactly one over four nine, six.
And our principal says. The we should take the mathematically correct answer, which is one and round it.
So that will be exactly one. If you're not sure, it's exactly one you could.
Well, in MATLAB, you can tell that's one. If it weren't exactly, I would say 1.000.
But another way to look at it is to look at the answer minus one.
And you see that's represented as zero. It's not ten to the -60.
Now, by contrast, suppose I say one over 49.
And then multiply that by 49. Now a priori, you don't know what this answer will be.
Of course you can guess that it won't be exactly one because I'm showing it to you.
But I priori it's a pretty safe bet that one over 49 will be an approximate number.
It'll be a floating point number that is not mathematically equal to 1/49 when you then multiply that by 49.
Do you get something that round to one? Well, I'm not quick enough to do that computation.
But as it happens, you don't. You can see there it's not exactly one.
And if I compare that answer to one, I see there's been a slight error, which, by the way, is probably two to the -53.
Notice how exact that is? It really is exactly two to the -53.
The error we've just made, 49 is the first number for which this is true.
So amazingly, all of this stuff, one seventh works out, one ninth, 1/11, 1/13.
They're all fine until 1/49. Okay.
Let's do a few more examples of rounding. Suppose I say the square root of a half.
Now that's of course not. Exactly represent a ball in the floating point system.
This was proved by the Greeks more than 2000 years ago.
So if I take it square and subtract the half, you don't know whether that will be exactly zero or not.
You know that the answer squared mathematically will not be exactly two, but oh, sorry a half.
But whether it rounds to a half, you don't know. Now as it happens, it doesn't.
So there again, we've got two to the -53. Exactly.
On the other hand, if I say a square root of a quarter, well, you know, that will be exact, because the quarter is represented in the binary system.
The square root of it will be calculated. Exactly. And then rounded to floating point.
Well, the exact square root is, of course, one half, which is already a floating point number.
So we know that if I square that and subtract a quarter, I'll get truly zero.
Exactly zero. Let's do a few more examples.
Suppose I say sine of zero. Now, this stuff I told you, concerns the four fundamental operations plus minus times and divide.
I haven't told you anything about what computers are obliged to do with the sign or the cosine or the Riemann Zeta function or anything like that.
And I don't know the precise status of that.
I don't actually remember whether I prescribed anything.
I think not. But most computers try to do their best to give you the correctly rounded sign function, at least for a broad range of arguments.
Not necessarily for all arguments, but certainly for numbers like zero that are of one, loosely speaking.
Computers will give you the correctly rounded sine. Obviously sine of zero is going to give you the correctly rounded number zero.
That's no surprise. Suppose I say pie?
Well, of course, that will be the correctly rounded number pie in your computer.
If I say sign of pie. Well, we have to think about what's happening here.
Pi is not exactly pi. It's the floating point approximation.
And then sine of that number introduces a new question. So you simply can't tell whether that's going to be exactly zero or not.
We have no reason to be sure it will be exactly zero. And in fact, it's not exactly zero.
Notice it's not even exactly two to the -53, because we've now composed two things an approximation to pi, followed by an approximation to the sun.
Suppose I say sign of a thousand times pi. Well, we should be able to predict what will happen a thousand times pi.
Well, let's see. Pi will have an hour at the end of the minor 16th level, so a thousand times pie will be off by about ten to the -13.
When we take the sign of that, we will probably get the correctly rounded sign of some number that's about ten to the -13 from a thousand PI.
So we'll probably get something about ten to the -13.
And there it is, ten to the -13. So similarly, if I said the sign of ten to the 20th times pi.
Well, ten of the 20th times pie is a number of size, ten of the 20 of from which we've lost a well represented in 16 digit precision.
So of course there will be no relevant accuracy left.
So this is a random number. Basically, there's your random number.
Now there's a lot of philosophical conversation of whether it's important on a computer for that to be exactly the right random number.
Maybe times PI is confusing matters. But suppose I said sign of ten to the 20.
Another random number one point of view is that.
There's no point in getting that number to any accuracy at all.
We're working in 16 digit arithmetic. Of course, we can't meaningfully represent the sign of ten to the 20.
We simply don't know enough. We lose all the digits when we take ten of the 20.
On the other hand, there are those who say that. It's important for other reasons that that should be essentially exactly the same of ten to the 20.
So we can test this. What happens if you take the sign of ten to the 20 squared plus the cosine of ten to the 20 squared?
So if these are essentially random numbers, then this result will be essentially a random number.
On the other hand, if they're if they have more structure than that, then the sum will be approximately one.
It is approximately one. So that tells us that the people designing the computer arithmetic were really, you know, maniacs.
It wasn't enough for them. In some sense, they really do try to give you the correctly rounded sign of ten to the 20th.
Even though you know the next floating point,
number ten to the 20th plus machine epsilon basically will be a totally different there's no continuity there at all.
They're random numbers, but they're the correct random numbers. Okay.
Let's do just a couple more things at the high and low ends.
These are aspects of Tripoli arithmetic that got a lot of attention in the eighties and are still a little controversial today.
One over zero is infinity, and that is a well-defined concept in floating point arithmetic.
One over infinity is zero, and it's just the usual zero.
So if I add two thirds to it, for example, I get two thirds.
If i00 divided by zero, I say not a number.
This is an introduction of I Tripoli arithmetic. Nine is not a number.
It obeys a very precisely prescribed arithmetic.
It's very convenient to have nans when these guys introduced Tripoli arithmetic.
Their hope was that nans would become an absolutely central part of how people compute.
You can do some really neat things with them. It hasn't turned out completely that way.
They're useful around the edges, but not as central, I think, as expected.
And sometimes manufacturers try to cut corners to speed things up in how they implement things like M&Ms and infinities.
Let's do another example. Suppose I say three times that you get na na na means we don't know.
It's a meaningless expression. If you multiply it by three, we still don't know.
You don't get three nine, you just get nine. So it's impossible to destroy a nine.
If I say zero times nine, it's still nine because you don't know, you know, was nine infinity?
Was it finite? You don't know. So you can't kill a man. It's just there forever to to the 1023 is a big number to to the 1024 is infinity.
So there you can see in this very precise sense that around ten to the 308 is the ceiling on floating point arithmetic, two to the -1 to 4.
You might expect to be zero, but it's not actually.
And that has to do with what are called D normalised floating point numbers down at the low end.
They allow the 52 bits to have some zeros up front.
And that's another controversial aspect of my Tripoli arithmetic.
It has some slick properties, but it adds to the complexity of implementations.
So in fact, if you go down to two to the -1 to 7 four, it's still non-zero.
But if you go to two to the -1 to 7 five, it's zero.
So floating point arithmetic represents numbers between ten to the minus three, eight and ten to the plus three, eight to full 16 digit precision.
And then there's this weird zone between ten to the minus three, three,
24 and ten to the -308, where the precision is fading out gradually rather than abruptly.
Okay. I think that's all I'm going to show you on this demo.
Most of the time, numerical people like me don't look at the details of floating point arithmetic,
but it's tremendously reassuring to know that they're all precisely prescribed in a deterministic fashion.
Now for the final few minutes, I want to give an absurdly speedy discussion of the consequences of rounding errors.
And let me mention three key terms stability, numerical stability, conditioning.
And backward error analysis.
Now, these are the three concepts that numerical analysts have found are at the heart of analysing the effects of rounding errors.
So if you think physically again, most of the time the fact that the continuum is made up of a lot of molecules doesn't matter.
But of course there are special cases where the quantum world really does come up to the macroscopic world.
So Superfluidity or Bose-Einstein condensation or something are cases where the microscopic and the macroscopic begin to blend together.
Analogously, when you do numerical computing in floating point, there are cases where the floating point stuff kills you.
In fact, it's more common than superfluidity. It's certainly something you have to be aware of.
You have to make sure that your algorithms don't have the property of amplifying the rounding errors in a exponential fashion, roughly speaking.
So if. I'm being very loose here.
If the rounding errors. Are amplified enough to cause trouble.
That's called numerical instability. And of course, that idea has been made much more precise with all sorts of things said about it.
And here's a slightly more precise way of thinking of that.
We say that a stable algorithm. Is one that gives the exact answer.
The exactly correct answer to a slightly perturbed problem.
So that's a very fundamental idea of numerical analysis.
Most of the time, it's not so much appropriate to ask.
Can we give nearly the right answer? What's more appropriate is to ask, can we give the right answer to nearly the posed problem?
Think, for example, of the sign of ten to the 20th. On a 16 digit computer.
A good algorithm is one that will give you the exact sign of some number that's within a relative error of ten to the -16 of ten to the 20th,
which means a random number. A stable algorithm is going to give you a random number when you ask for the sign of ten to the 20th.
That's okay. It wouldn't be reasonable to expect a computer to really do better than that.
So this is a the right compromise between what we would like to happen and what's reasonable to expect on a computer.
So slightly perturbed always means down at the level of rounding errors ten to the minus 16th kind of perturbation.
And all across numerical analysis, people make stable algorithms to do all sorts of things.
We then combine that idea of stability with the idea of conditioning.
So rounding errors amplified to cause trouble that we call instability.
If you have great sensitivity of answers to data.
That's called ill conditioning. So for example,
the sign of ten to the 20th is a very ill conditioned problem because a relative change in the
number ten to the 20 by one part in ten of the 16 is enough to change the answer completely.
So ten to the 20th is an extraordinarily ill condition problem, and condition number is ten to the 20th or something.
So any old sign algorithm is probably going to be stable.
Nevertheless, you can't expect sign of ten to the 20 necessarily to be correct because that problem is so ill conditioned.
And combining together these two ideas, you can analyse the accuracy of all sorts of problems in numerical analysis.
And the as it were, the fundamental axiom of numerical analysis is that stability plus conditioning implies accuracy.
Plus good conditioning. Well, conditioning. So that's a five minute summary of 50 years of very detailed, rigorous academic work.
All of all the time we ask, do we have a stable algorithm that doesn't amplifiers too much?
And is it a well conditioned problem where it's reasonable to expect that we should be able to compute a good answer?
Now, I just want to show you these things a little bit with these demos of the last three codes.
Let me just emphasise that some names connected with this.
Some of these ideas really go back to Turing, who was actually an even greater genius than they tell you.
A key person was Wilkinson, who really put them on the map in our field.
And he, incidentally, was Turing's assistant working on the pilot ace computer.
And then let me mention this this person and the book that we passed around by Nick Higham,
who is currently a leading figure in understanding this kind of analysis.
Okay. So codes 13, 14 and 15, just aim to illustrate some of these effects.
So for example, suppose I say X equals -40 and I'm interested in each of the x.
Now, of course, we don't know what algorithm our computer uses, but you can be sure it's it's a good one.
This is one of the most important of all functions. That's probably the correctly rounded result.
It's probably exactly e to the -40, correctly rounded.
I'm not certain of that, but I suspect it is. Let's look at an example of an unstable algorithm.
And the the this is a classic example. If you add up the Taylor series for E to the X and set X equal to -40,
that's a wonderfully unstable algorithm because the Taylor series has the huge terms before.
They've got small, huge terms of alternating sine powers of -40 divided by factorial.
Let's do it. Suppose I say index equals 0 to 150 and I'll say Taylor equals X to the I.
Divided by factorial of I. So I've just made the Taylor coefficients.
The first hundred and 50 Taylor coefficients, but that's plenty. For each of the x at x equals -40.
So in principle, if I sum them up, I should get about four times ten to the minus 18th.
But in fact, I get 1.8. So that has no relation whatsoever to the correct answer.
And the reason is, I've tried to sum up a bunch of things that are huge numbers of oscillating sine rounding errors.
This is rigorously derivable from our assumptions about floating point arithmetic.
Nothing has gone wrong. It's just an unstable algorithm.
Now, luckily, there do exist stable algorithms.
For example, if I said the same thing except with minus x.
Suppose I now add up the Taylor series?
I'll get eight of the 40, won't I? And that's a stable algorithm.
You can see it's beautifully computed. Eight of the 40.
So here's a stable algorithm for computing E to the -41 over E to the 40.
So you can see there we've illustrated the built in E to the -40 compared with the stable calculation of it, but also an unstable calculation of it.
The next example illustrates well and ill conditioned matrix problems.
So this is 14. Suppose I say a equals a random matrix of size 12 by 12.
There it is. That's an example of a well conditioned problem.
Probably. Of course it's random. I could have been unlucky, but I'm confident.
Let's say B equals the sum of the columns.
So I've just added up the columns, which means that we know mathematically what the answer to a backslash B should be.
Remember, if you add up the columns,
that means you're taking a linear combination of the columns with coefficients one therefore the solution of x equals B four.
That B will be the vector of all ones. And if I compute it, I find their very accurate format long.
Right. So we lost a digit maybe. But here's an example of a well-behaved problem.
Of course, backslash is a stable algorithm because MATLAB is a good code, a good language.
But moreover, this experiment shows that the problem was not too well conditioned.
Now, by contrast, let's look at the so called Hilbert matrix.
So I'll say A equals Hilbert of 12. It's an example of a high level condition matrix.
You can learn more about it by typing help. And let's again take a right hand side corresponding to the solution vector of all ones.
So if I now solve a x equals a backslash b again, mathematically I should get all ones.
But you see, I get numbers that are about one but not accurate.
So there we have a terribly ill conditioned matrix. It's a stable algorithm applied to a highly ill conditioned problem.
And in fact, the condition number of this matrix is bad.
What have I done? So the condition number of AA is something like ten to the 16th.
So it's amplifying small perturbations by six orders of magnitude.
Finally, I'll spend one more minute on 15, which is illustrating one of Wilkinson's great discoveries.
Suppose I say our equals 0.01 times.
I am going quickly because we're at the end. Oops, not too quickly.
Plus nine times the upper triangular ones of seven.
I'm just copying here. That turns out to be an ill conditioned upper triangular matrix.
If I say I don't know why I've done this anyway, let's do exactly what's on the sheet.
I'm now going to make a QR factorisation.
Of. A random matrix in order to get a random or orthogonal matrix.
Sorry this is going so fast and I'm now going to take a equals Q are.
So let me tell you what I've just done. I've just made a matrix whose q r factorisation involves an il condition r.
The details don't matter. Now I want to show you something astonishing that Wilkinson discovered.
You can forget the details of the matrix. It's just a seven by seven matrix.
Let's compute its QR factorisation. We've used whatever algorithm is built into MATLAB.
I can assure you it's a stable algorithm. Based on householder.
Now there's R and there's R two. If you look at them, you'll see they're totally different.
So our stable algorithm has given us a garbage.
Are the difference between the computed hour and the correct hour is huge.
Now, here's Q and here's Q2. If you look especially on the bottom right, you'll see those are totally different.
So we have a garbage. Q So the norm of Q minus Q two is big.
So we just executed a QR factorisation on the computer.
The Q we've gotten is no good and the R we've gotten is no good.
They're both very far from correct. Let's compare now.
The difference between A and Q.
Times are we find these garbage factors are actually fine.
Their product is ten to the -14. So this is all going by so fast.
At the end of a lecture, you probably don't appreciate how incredible it is.
But this was sort of the great discovery of Wilkinson's career that when you do calculations in linear algebra,
even though the individual answers may not be correct, they're correlated in such a way as to give you the finally correct thing.
In other words, Q2 and R2 are the exactly correct factors of a slightly perturbed matrix.
Okay, let's stop. Thank you.