In the lectures. So we'll see how that goes. This is scientific computing for Dphil students.
I'm the the professor of numerical analysis. So there's a numerical analysis group here, and I'm the head of that.
Of course, I'm curious to know who you are. So could we begin?
Raise your hand if you're a dphil student in mathematics.
Physics. Chemistry. That's less than usual.
That's interesting. Engineering, science. Computer science.
Plant science. What have I forgotten?
Which science? Animal science. Okay, good.
There are some mafia students here, I think. Who are you right now?
But who else have I overlooked? Okay.
Well, anyway, I hope you have fun and I hope you meet some of the other people in the room, you know?
One of the problems at Oxford is outside our colleges. We're not so good at meeting people from other department as well.
This is a chance because all of you are interested in numerical computing.
So let's see. Let's take a moment and look over just a few operational details.
So there's a pack of seven handouts. I hope you all got them. They're scattered about a few didn't.
I hope you've done the quiz at the end. We'll go over the answers to the quiz.
If you turn to the blue sheet, this contains all sorts of information about how the course runs.
There's also an outline of the lectures of the course on another sheet.
So as I say, I'm Nick Trevathan and there are two Tas for the course.
One of them is Mikhail's Levinsky in the back row there, and the other is Adrian Martinelli in the back row there.
So Mikel and Adrian know everything about scientific computing.
And by all means, send them email if you'd like to set up an appointment, or if you have a question they can answer by email.
And of course, you're also welcome to contact me by email. We're all in this building and not that hard to find.
Mikhail and Adrian. We'll be marking the problem, sheets of which there are four.
So each Tuesday of an even numbered week, we hand out an assignment and it's due the next Tuesday, the odd numbered week.
The first assignment doesn't count for very much while you're getting up to speed and the others count for more.
These are mostly computational work which you do in MATLAB, and I'm assuming you have access to MATLAB.
If you don't, you'll see there's a note here about how to get it from I.T services.
Anyway, we have these assignments and the course marks are based on those and they have to be turned in on schedule.
They can't be late because we hand out solutions at the lecture where the assignments are due.
So we do not accept late assignments. Now, I think handwritten on your blue sheet you'll see a mention of boxes 160 and 161.
Those refer to some hand in boxes that are across this atrium out there.
So that's where you turn in the assignments in those boxes across the atrium.
Michaela, Adrian, is there anything you need to tell us about the mechanics?
Okay. Right. Uh oh.
Well, there's a web page so you can see the the URL is here, and basically all the lecture notes will be put up on the Web page.
They're fairly terse lecture notes.
In principle, you could learn this material by looking at the lecture notes, but in practice, no human being can really do that.
You want to be here at the lectures? Things move quickly.
You want to be here at the lectures. There are, of course, 12 lectures this term and then another 12 in the next term.
Any questions about mechanical things? Okay.
I don't know details of whether your department counts this course for this or that.
That's up to you and your department. Okay.
So the way the course is structured, it's in this term, we're doing some linear algebra,
hopefully from a point of view that isn't entirely standard to you and then some optimisation.
But of course, I'd like to begin with some big picture about scientific computing.
So if you'll allow me, I'll write down the structure the first part of the course.
Part one I call sparse matrices and iterative methods.
And let me begin immediately, not on that subject, but with a view of the field.
To give you an idea of how I see this fitting in to scientific computing more generally.
So we'll call that 1.1 view of the field.
So the first question is what is the field scientific computing, numerical analysis?
I regard those as synonyms. Not everybody shares that view.
But to me, numerical analysis and scientific computing are both the discipline concerned with
algorithms of a numerical sort for mathematical problems of a continuous sort.
And the difference in the phrases maybe is what they connote.
Numerical analysis, of course, sounds like the more theoretical and scientific computing, the more applied end.
But I don't really make a distinction between those two. Now,
what I'd like to mention is that I'm trying to give you a very quick review
of a lot of material in these lectures and fairly classical in orientation.
And in particular, my view of the classical structure of this field is that there are three big parts of scientific computing.
One is linear algebra. One is differential equation.
So let's call that Odie and PD, by which we mean ordinary and partial differential equations.
And then the other is optimisation.
Now, of course, not everything fits this picture, but these are the three great areas I think, of classic scientific computing.
And we're beginning on the left here with linear algebra. Let me mention three particularly big parts of linear algebra.
Three big problems. One is solving systems of linear equations.
Another is least squares problems where your data fitting if you like.
Or another way to describe it is solving rectangular rather than square system of equations.
There's plenty of room scattered about, so please sit down.
And the third would be Eigenvalues Eigenvectors and singular values the SVD.
Once again, of course, those are not the only three parts of this field, but they are the three biggest, very important.
And there's an important distinction to make in all of this numerical linear algebra,
which is that in all of these three bubbles, there are sort of two main styles of work.
One is dense matrices, by which we mean matrices more or less all of whose entries are non-zero.
And you could call this classical numerical linear algebra, which means 1950s and 1960s, 1970s, that kind of era.
And then the other point of view would be the sparse numerical linear algebra.
And another aspect of that would be iterative methods.
And now we're getting into more modern aspects, if you like, parallel computing, large scale problems tend to be sparse.
These of course, are themes throughout scientific computing, and we'll certainly see them in linear algebra.
There's I'm you'll find I'm very interested in history, but to me, history gives you a perspective on a field that's interesting.
And there's a fascinating historical aspect of all this.
If you look at mathematics over the last few centuries, linear algebra was not important until computers were invented.
Of course, linear algebra was a part of mathematics, but it was never a hot field.
It was never a very big field. And then computers came along.
And since then, linear algebra has just grown and grown in its importance.
So, for example, every undergraduate degree on Earth nowadays probably requires linear algebra that might not have been true before the war.
Or shall we say. I'd like to say a word about why that happens.
And the word I'd like to say is schematic.
So let's imagine that we have a big problem of computational science, which, let's say is a nonlinear partial differential equation,
probably in three dimensions or something might be the navier-stokes equations and fluid mechanics.
So a big problem. How do we solve it computationally?
Well, whenever anything is nonlinear, you probably linear ising, which means locally you pretend it's linear and then somehow you have to iterate.
Take a few steps in order to converge to the nonlinear solution.
So given a nonlinear PD, the standard algorithms we all know about turn it into linear problems of a succession of linear problems, typically.
Now what about the PD? So the non linearity turns into linearity through linear ization, but the PD is discrete sized one way or another.
That's the standard method to solve a PD. You have a grid or finite elements or some sort of expansion, some sort of democratisation.
And that turns a problem of analysis into a problem of algebra.
So it really is a standard picture that Nonlinear PD e's turn into problems of linear algebra because of these two processes.
Von Neumann in the old days, you know, the late forties, early fifties, that kind of thing.
He began to see this, but I think those guys just had no idea how big linear algebra would grow.
They had no idea how powerful computers would grow.
Actually, von Neumann is another interesting point of history, so I think he predicted fairly accurately how good weather prediction would be in 2015.
He knew that we'd be good at predicting weather, but that the reason he got it right is that he made two errors that cancel.
It turns out it's about a million times harder than he would have thought,
and our machines are surely a million times more powerful than he would have predicted.
So in the end, he made the right prediction. Okay.
Just about this philosophy. Let me emphasise that I'm well aware that you in the room have different particular interests.
And chances are, most of you don't regard yourself as interested in linear algebra per se.
But let me just assure you that this is foundational material that you can't know too well.
Anyone doing computational science really needs to understand these things.
Okay. So that's it for that beginning. Now I think I'll switch over here and start talking about speed of solving linear algebra problems.
So let's ask the question how fast can we solve A X equals B?
So X equals B is the standard notation for a square system of linear equations.
A is a matrix an end by and matrix x and B are column vectors.
We're given B we want to find x, so I suppose I should write that down.
A is a capital N by capital n matrix.
I think I'll be inconsistent about capital versus lowercase b is a vector,
so we could call it an end by one vector just to emphasise that it's a column vector and X is an N by one vector of unknowns.
So we're given A and B and we want to find X.
And the question is how big an end is feasible.
It's just amazing to read some of the old stuff from the early days of computing when even end of ten or 20 seemed pretty impressive.
I have a quote here from 1951 from the great numerical analyst Wilkinson.
So he's describing the pilot ace computer, which is the one that Turing was in charge of setting up at the National Physical Laboratory.
And Wilkinson says by 1949, the major components of the pilot ace were complete and undergoing trials.
During 1951, a program for solving simultaneous linear algebraic equations was used for the first time on the 26th of June.
1951 was a landmark in the history of the machine.
For on that day, it first rivalled alternative computing methods by yielding by 3 p.m.
The solution of a set of 17 equations submitted the very same morning.
So in 1951, it takes 6 hours to do a 17 by 17 system of equations.
It's quicker now. Let me give some idea very schematic here of big values of N so following from
Wilkinson then let's say around 1950 and equals 20 would have been a big value.
1965 and equals 200 would have been a big value.
1980. And equals 2000 would have been a big value.
1995 I think the pattern is pretty clear and equals 20,000 would have been a big value.
2010 and equals 200,000.
Those are big matrices, reasonably close to the limit of what people can do.
Bigger than ordinary people ever do. No normal person in 1950 was solving 20 by 20 systems, and no normal person in 2010.
2010 was solving 200,000 by 200,000. But you can see it's a dramatic increase.
In fact, how big an increase is it? If you take these numbers? It's about a factor of ten to the fourth, right?
Every 15 years, roughly speaking, we're getting an increase in end of a factor of ten.
Now, of course, computers have speeded up a lot more than that. So this is how much PN has increased?
How much have computers increased? Well, more like ten to the 12th.
So. In this business, we talk about flops, floating point operations and flops and mega flops and giga flops and so on.
And roughly speaking, in this era, computers increased from flops to giga flops or to teraflops, I guess you'd say.
So you being young, probably no more of these prefixes than I do.
I get stuck after around exa. But they keep going, don't they?
As yet, they all sound like characters in Star Wars to me.
But it keeps going flops.
It means one floating point operation per second. TeraFlops means 1010 of the 12 floating point operations.
So let's say it would be flops and then killer flops.
I've never heard that expression. And then mega flops and then giga flops and then teraflops and petaflops and then X of flops.
And what, ten to the 21? Somebody remind me. We're at Terra heading for exa.
Well, actually, yacht. What is it?
Yacht, huh? That sounded right. It's something like that. Yeah, it's just around the corner.
This is about a factor of ten to the 12th. I'm being very rough here with my numbers, but.
Computers got 10 to 12 times faster and matrices got ten to the four times bigger.
And there's a reason for that. The algorithms, the standard classical algorithms take care of and cubed flops to solve A and B problem.
So we see it right there. Error during recording, Steve. Should I worry about that?
Oh. Okay. Good. That was easy.
So I think it's marvellous here how you can see in history something about algorithms.
If you didn't know then cubed you could maybe guess from looking at the data.
If you look at the top 500 list now, you can Google that top 500.
This is the list that's been maintained for 15 or 20 years.
Of the 500 fastest computers on Earth.
Most of them are so parallel that they bear little resemblance to the type of machines that we use.
That wasn't true 15 years ago, but now it is true.
The whole top 500 thing has become a little like Arnold Schwarzenegger in his muscle era, a little bit divorced from reality.
But it is incredible how things have speeded up. Let's see, I made a note here.
The fastest machine, the number one, which is a Chinese thing, is way into the of flops.
I think it's sort of ten to the 19th. Even the 500 fastest machine is in the hundreds of teraflops.
So if a teraflop is ten to the 12th out, ten to the 14th would put you, you know, not even the 500th fastest machine on earth.
It's just amazing. Let me make two comments about that.
There exist o of end to the 2.376 algorithms.
So although the standard algorithms are n cubed, in principle you can do better.
However, these have no practical value. They only exist theoretically are never used because the constants are far too awful.
So as far as anybody knows, at the moment, no practical algorithms are better than an cubed or maybe end of the 2.8.
And the other thing I want to mention is that this classical point of view is all about operation counts.
And that's considered these days a very classical point of view, because communication as opposed to floating point operation, also matters.
And. You see this on your own machine, but you certainly see it on these big,
massively parallel machines where they might have hundreds of thousands of cores.
And most of the concern is how to move the information from one core to another.
So as time evolves, this kind of end cube discussion is becoming a little further removed from high performance computing.
But I think now it still does tell us a lot about what most of us do most of the time.
Okay. So I like to play on the computer. And of course, in this course we play in MATLAB.
So let me. To give you an idea, just explore how fast things can go.
Now, the way I do this is I put together sheets telling you what commands I'm going to
type so that you can make notes if you like or go back and do it in your own time.
And these also end up online. And of course, part of the purpose here is also to teach you a bit of MATLAB.
I hope you've all used MATLAB, but I'm sure there are some people here who haven't used it before.
So for example, let me say a equals rand of five.
That produces a five by five matrix of random numbers, each one between zero and one.
And let me say B equals ones of five comma one that produces a right hand side vector B and let me shrink the font a little.
This command, by the way, is not a standard MATLAB command, so you can't do that.
You can do it, but much more slowly. And now let's solve the system of equations.
So if I say X equals a backslash, B, an algorithm is invoked and you don't even know what it is at this point necessarily.
But that's the MATLAB syntax for solving a system of equations.
Mathematically, it's a inverse times B, although the algorithm doesn't actually construct a inverse.
So if I do that, it solves the system. And there's the solution.
By the way, I'm in format, short mode, but if I type format long, then you could see all the digits, more or less that are stored.
Let's see how good it is. Suppose I say eight times X minus B.
So mathematically that should be zero because x equals B on the computer it's not zero, it's ten to the minus 15th times various things.
Because standard precision on most computers is.
Ten to the -16, more or less. If I type who's in MATLAB, it tells me what variables I've got.
And if you look at those variables, you'll see that there's a five by five matrix and a B and an X.
The number of bytes is always eight times the number of numbers.
So a five. Entry vector is 40 bytes because a standard word, a standard number is two words and a word is four bytes.
So 64 bits for each number. Let's make it a little bigger.
So suppose I now say a equals round end, which means normal random numbers of size 500 by 500.
And now I could say B equals ones of 501.
And I could say X equals a backslash b. No problem.
Oh, there is a problem. No problem.
If I now say who's you can see that the 500 by 500 matrix A each and there's 250,000 entries and each one takes eight bytes.
So we get 200,000 sorry, 2 million bytes, which when I was your age of course would have been a big matrix, but now completely trivial on any laptop.
What if I say A, X minus B, I get a long vector.
It looks like big numbers, but it's not really because at the top of this vector, it said ten to the -15.
In fact, if I look at the end entry of that long vector, you can see it's, well, ten to the -14.
So 500 by 500 matrix, we lost another digit to random effects, presumably.
So let's look at the speed. Let's try.
I'll say N equals 1000 and I'll say A equals around N of n.
So a thousand by thousand random matrix and I'll make a right hand side and then I'll say tick a backslash, B torque.
So that's how long it takes on my laptop to solve a 1000 by thousand system.
Let's try 2000. So on my laptop.
According to the end cubed, it should be eight times longer.
It's not quite there. We're not in the asymptotic regime. Timings are always complicated, but it's certainly more than four times longer.
Let's try 4000. So that should be eight times longer than half a second.
So the prediction is about 4 seconds. Works.
Okay, so the Mncube is certainly visible in these experiments.
By the way, on my MATLAB, on my laptop, you can also do 8000 and that works.
So this says a lot about storage.
Just a few years ago, an 8000 by 8000 matrix would not have fit on your laptop so conveniently.
But now that's fine. Let's look at the data more systematically.
So following the handout here, I'm going to try bank format, which means dollars and cents and I'll say four.
I equals 1 to 12.
Let's let RN be true to the AI and let's let a b round and then and let's say take a backslash ones of n1t equals talk and display an anti.
So here you see, we're doing an experiment in which we keep doubling and see how long it takes in seconds,
and most of them aren't visible in bank format.
But if if I displayed more digits, they'd not be very informative anyway because it's hard to time quick operations.
Throughout this course, I want to urge you always to interact with problems like this.
You always explore, do experiments all the time.
It's so easy in modern computing environments.
Let's do the same thing again, except make a plot. So I'll say again, for I equals 1 to 12 and equals two to the I.
And now I'll say tick a backslash ones of and comma one and I'll say T equals talk and I'll say, and then of I equals n and t t of I equals t.
And I've made a mistake, Bones. That's a new one.
Okay. I when I was getting ready for this lecture, at one point I had to take the word zeros and I left off the Z and it said Eros Unknown Command.
And I thought, Why did I type arrows? That confused me. Okay.
So now it's accumulating the data in a couple of vectors and we can now plot them.
So for example, suppose, I suppose I said plot and then against t t well then we get a plot of an n cubed
relationship of some sort number dimension of the matrix against the time.
That's kind of hard to see too. Well, let's improve it a bit. So I could say, for example, an empty like this and then.
S-H gee for show graph and I've just docked to the graph.
So now we can see at least the dots are visible as well as the lines connecting them.
But maybe a log log scale would be more enlightening.
So let's say log log of NN against t.
So there you can see kind of a typical experiment where the small sizes give you meaningless data.
Of course, it's not really true that a two by two matrix poses a real difficulty,
but when you try to time things, you always find funny stuff at the low end.
But at the high end, things are looking nice and clean.
There does seem to be some kind of an end cubed relationship. Let's see how close to end cubed it really does look.
So let's say hold on. So then when I plot the next curve, it doesn't erase the old curve.
And I'll say log log of and then one E minus ten.
Times and then. Oops, what have I done? One E minus ten times and then cubed, and I'll do it as a dashed red line.
So there you see, we've plotted a multiple of cubed superimposed on the plot.
It's no sense going down to ten to the minus ten. So let's say y lim of one E minus five, 1e2.
We're changing the Y limits of the plot there.
So now you can see more cleanly that it is reasonably close to an end cubed relationship with the slope getting a little bigger as we go up.
If we kept going, other things would happen due to caching and memory problems.
And just to have fun a little bit with MATLAB graphics, of course you can do things like put a label on,
I could say X label dimension and I could say Y label time in seconds and I could make it bigger.
I could say font size 20. And then I could say title.
The speed of MATLAB backslash font size 24 and I could say g text that's a MATLAB
command for pointing for putting something a bit of text on a plot interactively.
So let's say g text and cubed. Font size 24 in the colour red.
So when I issue that command, it then waits for me to click somewhere and I'll click right here and it puts that bit of text on the plot.
You might have heard that I was actually the first customer who ever bought MATLAB back in 1985.
So I've been doing this for a long time.
It's strange when I talk to people at the company because my closeness to their product is much greater than theirs.
Of course there are a few old timers, but as a rule, the people that work for MathWorks,
you know, have known MATLAB for ten years and I've known it for 35 years.
And it's it's very strange.
Okay. So let's move on and say a word about sparse matrices.
And maybe we should begin by thinking, why do we want end to be so large?
Where in the world would you ever get a matrix of size? 20,000?
You could get it from data of some sort. But traditionally, more often, you get it by disk sizing, a PDA.
That's where you get big matrices, you discrete high something. And the more space dimensions your discrete rising, the bigger the dimensions grow.
So let's do an example. Yeah, a two dimensional example.
Why do we want GN to be so large? Well, let's look at a small, small grid in 2D.
So here's a three by three grid.
I will imagine there are nodes like that.
And let's suppose they're numbered. One, two, three, four, five, six, seven, eight, nine.
And they're linked somehow because we're doing some kind of a finite difference approximation to a differential equation.
So if we're approximating the heat equation, there might be one value at each of those nodes.
If we're doing the navier-stokes equations, we might have a pressure and a couple of velocities or whatever you can if you're
modelling whether you might have a dozen different variables at each node.
So already with this trivial course grid in only two dimensions, we have a nine by nine matrix.
And the reason for numbering them is when you try to make a matrix, it depends on the numbering.
So the kind of matrix that would be associated with that grid, with those linkages would be a nine by nine matrix.
And. It would have a what we call a block three by three structure.
So it would be a nine by nine matrix, which naturally fell into three by three blocks,
each of which was a well into a three by three array of blocks, each of which is a three by three matrix.
So, for example, one is connected to two and four and two itself, let's say.
So there would be a nonzero entry in the second and the fourth column two is connected to one, three and five.
So the second row would have a non-zero n three in one, three and five.
And also itself three is connected to two and six.
So the third row in addition to itself would be connected to two and six.
So you keep going through this and you find that the block three by three matrix can be described as having diagonal blocks,
which are tri diagonal matrices, all with the same structure.
And then off diagonal blocks, which are diagonal matrices.
And these are all just acts. All right. These are zero.
So the matrix associated with that connectivity pattern is a nine by nine matrix.
It's a block tri diagonal. The blocks on the diagonal are themselves tri diagonal.
And the seven super diagonal blocks are diagonal.
Well, you can imagine this kind of thing scales up in all sorts of ways, in two dimensions as we make it,
and the dimension in each direction bigger than we have a quadratic effect here.
In three dimensions it would be cubic. So you very easily get big grids.
If you have 100 by 100 by 100, that's a million variables right there.
So it's to get a million by million. Matrix is nothing in computational science very easily done.
So we're always trying to deal with that problem, and there are two ways to deal with it.
So how do we beat the o of n cubed bottleneck?
In practice there are two ways, and one of them is to use a better machine or let's say, parallel computing.
And that idea has been around for a very long time.
It's not new, but it's getting more and more powerful with every passing year.
These days we have GPUs which make it easy for many people to use many, many cores.
So this is a very powerful strategy, just getting more important with the years, the other approach.
It is to take advantage of the sparsity. So special algorithms.
Exploiting Sparsity. Because often, though not always, the matrices you want to solve have mostly zeros in them.
Sparse means mostly zero. Enough zeros that you can take advantage of it.
As you can imagine, that's a highly developed field. There are some beautiful algorithms and we'll talk about a few of them.
And before I explore in MATLAB a little bit, I wanted to pass around a couple of books.
So this is a fantastic book by Trevathan and Bough on Numerical Linear Algebra.
And David Beau, my co-author, was a graduate student of mine when I was at Cornell,
and he then went on to Microsoft and then Google have done very well.
This is a beautiful little book by Toby Driscoll about learning MATLAB.
You don't need a book to learn MATLAB. Everything is online, of course, but it is very nice if you want to look at it.
So this is very appealing. Okay.
I wanted to go back and explore some sparse matrices before doing that.
Let me draw your attention to a couple of the handouts in the packet there.
So. The first one is the paper by Gilbert, Mueller and Schreiber from 1992.
This is a beautiful paper by three remarkable people.
Mueller is the inventor of MATLAB. Gilbert and Schreiber are leading figures in linear algebra and algorithms and graph algorithms and so on.
To me, this this was the moment at which MATLAB made the transition from being an educational tool to being a serious tool for scientific computing.
And it's gradually grown in practical importance with the decades they took
algorithms for sparse matrices that had been invented somewhat academically,
if you like, and then found a way elegantly to put them into MATLAB and wrote this beautiful paper about it.
It has a personal resonance in that, as you can see, it's dedicated to Gene Ghalib on the occasion of his 60th birthday.
Ghalib was a great, great figure in numerical linear algebra, and I remember going to that 60th birthday conference,
and now 23 years later, I've just had my 60th birthday conference so I can when I look at this, I feel the generations moving along.
Anyway, one of the themes of this course is to remind you that the tools we use have an intellectual background.
And this is an example of how a big part of MATLAB is based on some very major research over the years, carried forth by some very major people.
On the back of that handout is a pointer to unpack, which is the software system that MATLAB calls these days for sparse linear algebra.
It wasn't there when Gilbert and Mueller and Schreiber were writing their paper.
But now this is one of the standard fast packages for sparse linear algebra.
It's due to Timothy Davis, it says here, University of Florida.
He's now at Texas A&M, I think. Anyway, AfPak is something to know about if you care about the details of algorithms and in fact,
in MATLAB instead of backslash, you can type lind solve and you get more direct access to einfach.
If you want to tune the knobs a little bit. There's also a hand out of the MATLAB guide.
This is a book by Nick Higham and does Higham beautiful, thorough introduction to many aspects of MATLAB.
And what I've done here is to copy on two pages illegally, I guess the entire chapter on sparse matrices.
So it's actually worth reading this.
I would recommend you find a few minutes and read these eight pages copied on to to to get an idea of how MATLAB stores and handles sparse matrices.
Well, let's play around a little bit. Let's look at a sparse matrix in MATLAB.
So I'll say format short and I'll say a equals 1 to 0 for.
Now, if we look at a it's obvious from the display what it is.
And if I say who's no surprises there, it's just a two by two matrix for entry.
So 32 bytes. Now let's take another matrix called B, which is the specifications of A.
So it's the same matrix mathematically, but stored in a sparse format.
And you can see that there are three lines of data because I only had three non-zero entries.
It's only telling you the non-zero entries.
If I say who's you can see B actually is more wasteful of storage per entry.
It takes more bites. But of course, if you have a big matrix that's mostly sparse, it will be a big win to store it in the sparse fashion.
So let's do one other thing. Well, let's go up to dimension one.
Let's say A equals. I have 100 and let's say B equals sparse identity SPI of 100.
And let's do A who's. So now you see that A has 10,000 entries.
So 80,000 bytes, whereas B 10,000 entries, but only 100 of them are zero.
So for whatever reason, it has only 2408 bytes.
So there is sparse storage, which enables you to hold much bigger mattresses in memory.
But of course, also algorithms will take advantage of the zero structure to go faster.
Let's play around a little bit to give you a sense of MATLAB syntax.
So if I say Spy of A, I get a picture of the nonzero entries of a.
Sorry. I'm having trouble with my mouse here. Okay.
If I say Spy of Bee, it's the same because mathematically they're the same matrix.
Notice at the bottom NZ does not stand for a rugby team.
It stands for the number of non-zero-sum. Now suppose I say be of 40 coal and 6010, 20, 40 equals three.
I've just taken a bunch of entries in the Matrix B and changed their value from 0 to 3.
And of course, the data structure has been what has done whatever is necessary to accommodate that.
If I do a spy plot, you can see where the non zeros have appeared.
Or suppose I say spy of B plus B transpose b prime means B transpose.
Then we get a symmetric structure. Or suppose I say spy of B plus B transpose squared.
So before pressing return, I can ask you to think about what that might look like,
what will happen when I square that matrix to the non-zero structure?
Well, that's what happens. Now let's solve a system of equations or two.
So just to give you a sense of how MATLAB sets up some random matrices, let me say a equals ESP ran the SIM of 7.25.
Now what that means is. That was the command sp round sim of 7 to 5.
It means make a random matrix of seven by seven, sparsely stored and approximately 0.25 of the entries should be non-zero.
So SPI of a is a seven by seven matrix and about a quarter of the entries are non zeros.
If I said full are they then you could actually see this seven by seven matrix.
Now let's make a big one. I'll say equals ESP round sim of 10,000.
And only a fifth of a percent of them should be non-zero.
So now if I say Spy of A, it looks dense, doesn't it?
But that's because the dots are too big. It's not really dense.
Most of those entries are zero. We could zoom in by a factor of two.
Why? Zoom in another factor of two. Every time we zoom in, we begin to see.
So here we have a sparse random matrix where 99.8% of the entries are zero.
Let's add a multiple of the identity to that matrix to make it a better behaved matrix.
So I'm going to add ten times the sparse identity of size 10,000,
and I'm going to make a right hand side ones of 10,000, comma one, and I'm going to solve a x equals a backslash p.
Now earlier we solved a 4000 by 4000 matrix.
Right. And that took 4 seconds. This is a 10,000 by 10,000 matrix.
Sure enough. It might finish eventually, but it's not going to finish quickly.
However we can solve it using a sparse algorithm. When I say backslash, it isn't using an algorithm that exploits this type of sparsity.
So it takes a long time, but I'm going to use a iterative algorithm called preconditioned conjugate gradients.
I'll say x equals PC G of A, comma B, and we'll be playing with this.
Salman. Preconditioned to conjugate gradients you can see in with its default parameters.
It doesn't succeed, but it will succeed if we.
Increase the number of allowed iterate.
So I'll say x equals pre-condition conjugate gradients of a b, I'll give it a tolerance of ten to the -12 and I'll allow it to take 200 steps.
What did I do wrong? All right. Lowercase, right?
Yeah. So you can see how quickly that happened.
It didn't take 20 or 50 seconds. It converged in a fraction of a second.
Let's do it again, say less than a second. It converged and presumably got the right answer.
Let's see if I say norm of A times X minus B,
norm means I'm taking the product of the 10,000 by 10,000 matrix times that 10,000 vector subtracting off an and another 10,000 vector.
Let's look at the infinity norm. The biggest entry. So the biggest entry in the error is about ten to the minus 12th.
If you're curious to learn more. Of course, there are many ways to learn things these days.
One of them, if you're just doing it within MATLAB, is to put on the more features.
So it only fills one screen at a time.
And then if I say Help Spa fun for spa matrices, you can get a sense of some of the commands that are available for dealing with sparse matrices,
random, random, symmetric diagonal matrices, identity and so on and so on.
Here's another example. If I say a equals gallery of Watson.
Well, Andy Watson is a faculty member here in the Numerical Analysis Group.
He's Oxford's main expert on numerical linear algebra.
And MATLAB has a matrix named after him. I won't go into the details, but here's what it looks like.
You can see it's nice structure there. Let's square it.
Let's quit. So there's a lot of stuff you can find in gallery.
Okay, so we're at the end of this lecture, except that I want to hand out the solutions to the quiz.
I'm guessing most of you knew about five of these acronyms.
Who knows? Here's what they stand for.
And let me just emphasise that the first assignment is due on Tuesday at 9:00 and it's in the box at 160 or 161 across the atrium from here.
Okay. Thank you very much.