Lecture 24: Advanced Topics (cont.)

Related Resources

Assignments
Exams

• Download English-US transcript (PDF)

-- week of 6.046. Woohoo!

The topic of this final week, among our advanced topics, is cache oblivious algorithms. This is a particularly fun area, one dear to my heart because I've done a lot of research in this area. This is an area co-founded by Professor Leiserson. So, in fact, the first context in which I met Professor Leiserson was him giving a talk about cache oblivious algorithms at WADS '99 in Vancouver I think. Yeah, that has to be an odd year. So, I learned about cache oblivious algorithms then, started working in the area, and it's been a fun place to play.

But this topic in some sense was also developed in the context of this class. I think there was one semester, probably also '98-'99 where all of the problem sets were about cache oblivious algorithms. And they were, in particular, working out the research ideas at the same time. So, it must have been fun semester. We consider doing that this semester, but we kept it to the simple.

We know a lot more about cache oblivious algorithms by now as you might expect. Right, I think that's all the setting. I mean, it was kind of developed also with a bunch of MIT students in particular, M.Eng. student, Harold Prokop.

It was his M.Eng. thesis.

There is all the citations I will give for now.

I haven't posted yet, but there are some lecture notes that are already on my webpage.

But I will link to them from the course website that gives all the references for all the results I'll be talking about.

They've all been done in the last five years or so, in particular, starting in '99 when the first paper was published. But I won't give the specific citations in lecture. And, this topic is related to the topic of last week, multithreaded algorithms, although at a somewhat high level.

And then it's also dealing with parallelism in modern machines.

And we've had throughout all of these last two lectures, we've had this very simple model of a computer where we have random access. You can access memory at a cost of one. You can read and write a word of memory. There is some details on how big a word can be and whatnot. It's pretty basic, simple, flat model. And, at the multithreaded algorithm is the idea that, well, maybe you have multiple threads of computation running at once, but you still have this very flat memory. Everyone can access anything in memory at a constant cost. We're going to change that model now. And we are going to realize that a real machine, the memory of a real machine is some hierarchy. You have some CPU, you have some cache, probably on the same chip, level 1 cache, you have some level 2 cache, if you're lucky, maybe you have some level 3 cache, before you get to main memory.

And then, you probably have a really big disk and probably there's even some cache out here, but I won't even think about that. So, the point is, you have lots of different levels of memory and what's changing here is that things that are very close to the CPU are very fast to access. Usually level 1 cache you can access in one clock cycle or a few.

And then, things get slower and slower.

Memory still costs like 70 ns or so to access a chunk out of.

And that's a long time. 70 ns is, of course, a very long time. So, as we go out here, we get slower. But we also get bigger.

I mean, if we could put everything at level 1 cache, the problem would be solved. But what would be a flat memory. Accessing everything in here, we assumed takes the same amount of time.

But usually, we can't afford, it's not even possible to put everything in level 1 cache.

I mean, there's a reason why there is a memory hierarchy.

Does anyone have a suggestion on what that reason might be?

It's like one of these limits in life.

Yeah? Fast memory is expensive.

That's the practical limitations indeed, that you could try to build more and more at level 1 cache and maybe you could try to, well, yeah.

Expenses is a good reason, and practically, that's why they may be the sizes are what they are.

But suppose really fast memory were really cheap.

There is a physical limitation of what's going on, yeah? The speed of light.

Yeah, that's a bit of a problem, right?

No matter how much, let's suppose you can only fit so many bits in an atom. You can only fit so many bits in a particular amount of space. If you want more bits, and you need more space, and the more space you have, the longer it's going to take for a round-trip.

So, if you assume your CPU is like this point in space, so it's relatively small and it has to get the data in, the bigger the data, the farther it has to be away.

But, you can have these cores around the CPU that are, we don't usually live in 3-D, and chips were usually in 2-D, but never mind. You can have the sphere that's closer to the CPU that's a lot faster to access.

And as you get further away it costs more.

And that's essentially what this model is representing, although it's a bit approximated from the intrinsic physics and geometry and whatnot.

But that's the idea. The latency, the round-trip time to get some of this memory has to be big.

In general, the costs to access memory is made up of two things.

There's the latency, the round-trip time, which in particular is limited by the speed of light.

And, plus the round-trip time, you also have to get the data out. And depending on how much data you want, it could take longer. OK, so there's something.

There could be, get this right, let's say, the amount of data divided by the bandwidth.

OK, the bandwidth is at what rate can you get the data out?

And if you look at the bandwidth of these various levels of memory, it's all pretty much the same.

If you have a well-designed computer the bandwidths should all be the same. OK, as you can still get data off disc really, really fast, usually at about the speed of your bus, and that the bus gets the CPU hopefully as fast as everything else.

So, even though they're slower, they're really only slower in terms of latency. And so, this part is maybe reasonable. The bandwidth looks pretty much the same universally. It's the latency that's going up. So, if the latency is going up but we still get to divide by the same amount of bandwidth, what should we do to make the access cost at all these levels about the same? This is fixed.

Let's say this is increasing, but this is still staying big.

What could we do to balance this formula?

Change the amounts. As the latency goes up, if we increase the amount, then the amortized cost to access one element will go down. So, this is amortization in a very simple sense. So, this was to access a whole block, let's say, and this amount was the size of the block. So, the amortized cost, then, to access one element is going to be the latency divided by the size of the block, the amount plus one over the bandwidth. OK, so this is what you should implicitly be thinking in your head.

So, I'm just dividing here by the amounts because the amount is how many elements you get in one access, let's suppose.

OK, so we get this formula for the amortized cost.

The one over bandwidth is going to be good no matter what level we are on, I claim. There's no real fundamental limitation there except it might be expensive.

And the latency week at the amortized out by the amounts, so whatever the latency is, at the latency gets bigger out here, we just get more and more stuff and then we make these two terms equal, let's say. That would be a good way to balance things. So what particular, disc has a really high latency. Not only is there speed of light issues here, but there's actually the speed of the head moving on the tracks of the disk.

That takes a long time. There's a physical motion.

Everything else here doesn't usually have physical motion.

It's just electric. So, this is really, really slow and latency, so when you read something out of disk, you might as well read a lot of data from disc, like a megabyte or so. It's probably even old these days. Maybe you read multiple megabytes when you read anything from disk if you want these to be matched. OK, there's a bit of a problem with doing that. Any suggestions what the problem would be? So, you have this algorithm.

And, whenever it reads something off of desk, it reads an entire megabyte of stuff around the element it asked for. So the amortized cost to access is going to be reasonable, but that's actually sort of assuming something. Yeah?

Right. I'm assuming I'm ever going to use the rest of that data. If I'm going to read 10 MB around the one element that asked for, I access A bracket I, and I get 10 million items from A around I, it would be kind of good if the algorithm actually used that data for something.

It seems reasonable. So, this would be spatial locality. So, we want, I mean the goal of this world in cache oblivious algorithms and cache efficient algorithms in general is you want algorithms that perform well when this is happening.

So, this is the idea of blocking.

And we want the algorithm to use all or at least most of the elements in a block, a consecutive chunk of memory.

So, this is spatial locality.

Ideally, we'd use all of them right then.

But I mean, depending on your algorithm, that's a little bit tricky. There is another issue, though. So, you read in your thing into, read your 10 MB into main memory, let's say, and your memory, let's say, is at least, these days you should have a 4 GB memory or something.

So, you could read and actually a lot of different blocks into main memory. What you'd like is that you can use those blocks for as long as possible.

Maybe you don't even use them. If you have a linear time algorithm, you're probably only going to visit each element a constant number of times. So, this is enough.

But if your algorithm is more than linear time, you're going to be accessing elements more than once.

So, it would be a good idea not only to use all the elements of the blocks, but use them as many times as you can before you have to throw the block out. That's temporal locality.

So ideally, you even reuse blocks as much as possible.

So, I mean, we have all these caches.

So, I didn't write this word. Just in case I don't know how to spell it, it's not the money. We should use those caches for something. I mean, the fact that they store more than one block, each cache can store several blocks. How many?

Well, we'll get to that in a second.

OK, so this is the general motivation, but at this point the model is still pretty damn ugly.

If you wanted to design an algorithm that runs well on this kind of machine directly, it's possible but pretty difficult, and essentially never done, let's say, even though this is what real machines look like.

At least in theory, and pretty much in practice, the main thing to think about is two levels at a time.

So, this is a simplification where we can say a lot more about algorithms, a simplification over this model. So, in this model, each of these levels has different block sizes, and a different total sizes, it's a mess to deal with and design algorithms for. If you just think about two levels, it's relatively easy. So, we have our CPU which we assume has a constant number of registers only.

So, you know, once it has a couple of data items, you can add them and whatnot.

Then we have this really fast pipe.

So, I draw it thick to some cache.

So this is cache. And, we have a relatively narrow pipe to some really big other storage, which I will call main memory. So, I mean, that's the general picture. Now, this could represent any two of these levels. It could be between L3 cache and make memory. That's maybe, what? The naming corresponds to best.

Or cache could in fact be main memory, what we consider the RAM of the machine, and what's called a memory over there to be the disk. It's whatever you care about.

And usually, if you have a program, that's what usually we assume everything fits in main memory.

Then you care about the caching behavior.

So you probably look between these two levels.

That's probably what matters the most inner program because the cost differential here is really big relative to the cost differential here. If your data doesn't even fit it main memory, and you have to go to disk, then you really care about this level because the cost differential here is huge. It's like six orders of magnitude, let's say. So, in practice you may think of just two memory levels that are the most relevant.

OK, now I'm going to define some parameters.

I'm going to call them cache and make memory just for clarity because I like to think of main memory just the way it used to be. And now all we have to worry about is this extra thing called cache.

It has some bounded size, and there's a block size.

The block size is B. and a number of blocks is M over B. So, the total size of the cache is M. OK, main memory is also blocked into blocks of size B. And we assume that it has essentially infinite size. We don't care about its size in this picture. It's whatever is big enough to hold the size of your algorithm, or data structure, or whatever. OK, so that's the general model. And for strange, historical reasons, which I don't want to get into, these things are called capital M and capital B.

Even though M sounds a lot like memory, it's really for cache, and don't ask. OK, this is to preserve history. OK, now what do we do with this model? It seems nice, but now what do we measure about it?

What I'm going to assume is that the cache is really fast.

So the CPU can access cache essentially instantaneously.

I still have to pay for the computation that the CPU is doing, but I'm assuming cache is close enough that I don't care.

And that may memory is so big that it has to be far away, and therefore, this pipe is a problem.

I mean, what I should really draw is that pipe is still thick, but is really long. So, the latency is high.

The bandwidth is still high. OK, and all transfers here happened as blocks. So, when you don't have something, so the idea is CPU asks for A of I, as for something in memory, if it's in the cache, it gets it. That's free.

Otherwise, it has to grab the entire block containing that element from main memory, brings it into cache, maybe kicks somebody out if the cache was full, and then the CPU can use that data and keep going.

Until it accesses something else that's not in cache, then it has to grab it from main memory.

When you kick something out, you're actually writing back to memory. That's the model.

So, we suppose the accesses to cache are free.

But we can still think about the running time of the algorithm. I'm not going to change the definition of running time. This would be the computation time, or the work if you want to use multithreaded lingo, computation time. OK, so we still have time, and T of N will still mean what it did before.

This is just an extra level of refinement of understanding of what's going on. Essentially, measuring the parallelism that we can exploit out of the memory system, that when you access something you actually get B items. So, this is the old stuff.

Now, what I want to do is count memory transfers.

These are transfers of blocks, so I should say block memory transfers between the two levels, so, between the cache and main memory. So, memory transfers are either reading reads or writes. Maybe I should say that.

These are number of block reads and writes from and to the main memory. OK, so I'm going to introduce some notation. This is new notation, so we'll see how it works out. MT of N I want to represent the number of memory transfers instead of just normal time of the problem of size N. Really, this is a function that depends not only on N but also on these parameters, B and M, in our model. So, this is what it should be, MT_B,M(N), but that's obviously pretty messy, so I'm going to stick to MT of N.

But this will always, because mainly I care about the growth in terms of N. well, I care about the growth in terms of all things, but the only thing I could change is N. So, most of the time I only think about, like when we are writing recurrences, only N is changing. I can't recurse on the block size. I can't recurse on the size of cache. Those are given to me.

They're fixed. OK, so we'll be changing N mainly. But B and M always matter here.

They're not constants. They're parameters of the model. OK, easy enough.

This is something called the disk access model, if you like DAM models, or the external memory model, or the cache aware model. Maybe I should mention that; this is the cache aware. In general, you have some algorithm that runs on this kind of model, machine model.

That's a cache aware algorithm. OK, we're not too interested in cache aware algorithms. We've seen one, B trees. B trees are cache aware data structure. You assume that there is some block size, B, underlying.

Maybe you didn't see exactly this model.

In particular, it didn't really matter how big the cache was because you just wanted to know.

When I read B items, I can use all of them as much as possible and figure out where I fit among those B items, and that gives me log base B of N memory transfers instead of log N, which would be, if you just threw your favorite balanced binary search tree. So, log base B of N is definitely better than log base 2 of N.

B trees are a cache aware algorithm.

OK, what we would like to do today and next lecture is get cache oblivious algorithms. So, there's essentially only one difference between cache aware algorithms and cache oblivious algorithms. In cache oblivious algorithms, the algorithm doesn't know what B and M are.

So this is a bit of a subtle point, but very cool idea.

You assume that this is the model of the machine, and you care about the number of memory transfers between this cache of size M with blocking B, and main memory with blocking B. But you don't actually know what the model is. You don't know the other parameters of the model. It looks like this, but you don't know the width. You don't know the height.

Why not? So, the analysis knows what B and M are. We are going to write some algorithms which look just like boring old algorithms that we've seen throughout the lecture. That's one of the nice things about this model. Every algorithm we have seen is a cache oblivious algorithm, all right, because we didn't even know the word cache in this class until today.

So, we already have lots of algorithms to choose from.

The thing is, some of them will perform well in this model, and some of them won't.

So, we would like to design algorithms that just like our old algorithms that happened to perform well in this context, no matter what B and M are. So, another way this is the same algorithm should work well for all values of B and M if you have a good cache oblivious algorithm.

OK, there are a few consequences to this assumption.

In a cache aware algorithm, you can explicitly say, OK, I'm blocking my memory into chunks of size B.

Here they are. I was going to store these B elements here, these B elements here, because you know B, you can do that.

You can say, well, OK, now I want to read these B items into my cache, and then write out these ones over here. You can explicitly maintain your cache. With cache oblivious algorithms, you can't because you don't know what it is.

So, it's got to be all implicit.

And this is pretty much how caches work anyway except for disk. So, it's a pretty reasonable model. In particular, when you access an element that's not in cache, you automatically fetch the block containing that element.

And you pay one memory transfer for that if it wasn't already there. Another bit of a catch here is, what if your cache is full? Then you've got to kick some block out of your cache. And then, so we need some model of which block gets kicked out because we can't control that.

We have no knowledge of what the blocks are in our algorithm.

So what we're going to assume in this model is the ideal thing, that when you fetch a new block, if your cache is full, you evict a block that will be used farthest in the future.

Sorry, the furthest. Farthest is distance.

Furthest is time. Furthest in the future.

OK, this would be the best possible thing to do.

It's a little bit hard to do in practice because you don't know the future generally, unless you're omniscient.

So, this is a bit of an idealized model.

But it's pretty reasonable in the sense that if you've read the reading handout number 20, this paper by Sleator and Tarjan, they introduce the idea of competitive algorithms.

So, we only talked about a small portion of that paper that moved to front heuristic for storing a list.

But, it also proves that there are strategies, and maybe you heard this in recitation.

Some people covered it; some didn't, that these are called paging strategies.

So, you want to maintain some cache of pages or blocks, and you pay whenever you have to access a block that's not in your cache. The best thing to do is to always kick out the block that will be used farthest in the future because that way you'll use all the blocks that are in there. This turns out to be the offline optimal strategy if you knew the future.

But, there are algorithms that are essentially constant competitive against this strategy.

I don't want to get into details because they're not exactly constant competitive. But they are sufficiently constant competitive for the purposes of this lecture that we can assume this, not have to worry about it.

Most of the time, we don't even really use this assumption. But there it is.

That's the cache oblivious model.

It makes things cleaner to think about just anything that should be done, will be done.

And you can simulate that with least recently used or whatever good heuristic that you want to that's competitive against the optimal. OK, that's pretty much the cache oblivious algorithm. Once you have the two level model, you just assume you don't know B and M.

You have this automatic request in writing, and whatnot.

A little bit more to say, I guess, it may be obvious at this point, but I've been drawing everything as tables.

So, it's not really clear what the linear order is.

Linear order is just the reading order.

So, although we don't explicitly say it most of the time, a typical model is that memory is a linear array.

Everything that you ever store in your program is written in this linear array. If you've ever programmed in Assembly or whatever, that's the model.

You have the address space, and any number between here and here, that's where you can actually, this is physical memory. This is all you can write to.

So, it starts at zero and goes out to, let's call it infinity over here. And, if you allocate some array, maybe it occupies some space in the middle.

Who knows? OK, we usually don't think about that much. What I care about now is that memory itself is blocked in this view.

So, however your stuff is stored in memory, it's blocked into clusters of length B.

So, if this is, let me call it one and be a little bit nicer. This is B.

This is position B plus one. This is 2B, and 2B plus one, and so on. These are the indexes into memory. This is how the blocking happens. If you access something here, you get that chunk from U, round it down to the previous multiple of B, round it up to the next multiple of B. That's what you always get.

OK, so if you think about some array that's maybe allocated here, OK, you have to keep in mind that that array may not be perfectly aligned with the blocks.

But more or less it will be so we don't care too much.

But that's a bit of a subtlety there.

OK, so that's pretty much the model.

So every algorithm we've seen, except B trees, is a cache oblivious algorithm. And our question is, now, we know how everything runs in terms of running time.

Now we want to measure the number of memory transfers, MT of N. I want to mention one other fact or theorem. I'll put it in brackets because I don't want to state it precisely.

But if you have an algorithm that is efficient on two levels, so in other words, what we're looking at, if we just think about the two level world and your algorithm is cache oblivious, then it is efficient on any number of levels in your memory hierarchy, say, L levels. So, I don't want to define what efficient means. But the intuition is, if your machine really looks like this and you have a cache oblivious algorithm, you can apply the cache oblivious analysis for all B and M.

So you can analyze the number of memory transfers here, here, here, here, and here.

And if you have a good cache oblivious algorithm, the performances at all those levels has to be good.

And therefore, the whole performance is good.

Good here means asymptotically optimal up to constant factors, something like that. OK, so I don't want to prove that, and you can read the cache oblivious papers.

That's a nice fact about cache oblivious algorithms.

If you have a cache aware algorithm that tunes to a particular value of B, and a particular value of M, you're not going to have that problem.

So, this is one nice feature of cache obliviousness.

Another nice feature is when you are coding the algorithm, you don't have to put in B and M.

So, that simplifies things a bit.

So, let's do some algorithms. Enough about models.

OK, we're going to start out with some really simple things just to get warmed up on the analysis side.

The most basic thing you can do that's good in a cache oblivious world is scanning. So, scanning is just visiting the items in an array in order. So, visit A_1 up to A_N in order. For some notion of visit, this is presumably some constant time operation.

For example, suppose you want to compute the aggregate of the array. You want to sum all the elements in the array. So, you have one extra variable using, but you can store that in a register or whatever, so that's one simple example. Sum the array.

OK, so here's the picture. We have our memory.

Each of these cells represents one item, one element, log N bits, one word, whatever.

Our array is somewhere in here. Maybe it's there.

And we go from here to here to here to here.

OK, and so on. So, what does this cost?

What is the number of memory transfers?

We know that this is a linear time algorithm.

It takes order N time. What does it cost in terms of memory transfers? N over B, pretty much.

We like to say it's order N over B plus two or one in the big O. This is a bit of worry.

I mean, N could be smaller than B.

We really want to think about all the cases, especially because usually you're not doing this on something of size N. You're doing it on something of size k, where we don't really know much about k.

But in general, it's N over B plus one because we always need at least one memory transfer to look at something, unless N is zero. And in particular, it's plus two if you care about the constants.

If I don't write the big O, then it would be plus two at most because you could essentially waste the first block and that everything is fine for awhile.

And then, if you're unlucky, you essentially waste the last blocked. There is just one element in that block, and you're not getting much out of it.

Everything in the middle, though, every block between the first and last block has to be full.

So, you're using all of those elements.

So out of the N elements, you only have N over B blocks because the block has B elements.

OK, that's pretty trivial. Let me do something slightly more interesting, which is two scans at once.

OK, here we are not assuming anything about M.

we're not assuming anything about the size of the cache, just that I can hold a single block.

The last block that we visited has to be there.

OK, you can also do a constant number of parallel scans.

This is not really parallel in the sense of multithreaded, bit simulated parallelism. I mean, if you have a constant number, do one, do the other, do the other, come back, come back, come back, all right, visit them in turn round robin, whatever. For example, here's a cute piece of code. If you want to reverse an array, OK, then you can do it. This is a good puzzle.

You can do it by essentially two scans where you repeatedly swapped the first and last element.

So I was swapping A_i with N minus i plus one, and just restart at one. So, here's your array.

Suppose this is actually my array.

I swap these two guys, and I saw these two guys, and so on. That will reverse my array, and this should work hopefully the middle as well if it's odd.

It should not do anything. And you can view this as two scans. There is one scan that's coming in this way. There's also a reverse scan, ooh, some more sophisticated, coming back this way.

Of course, reverse scan has the same analysis.

And as long as your cache is big enough to store at least two blocks, which is a pretty reasonable assumption, so let's write it. Assuming the number of blocks in the cache, which is M over B, is at least two in this algorithm, the number of memory transfers is still order N over B plus one.

OK, the constant goes up maybe, but in this case it probably doesn't. But who cares.

OK, as long as you're doing a constant number of scans, and some constant number of arrays, it happens to be one of them's reversed, whatever, it will take, we call this linear time. It's linear in the number of blocks in your input. OK, great.

So now you can reverse an array: exciting.

Let's try another simple algorithm on another board.

Let's try binary search. So just like last week, we're going back to our basics here.

Scanning we didn't even talk about in this class.

Binary search is something we talked about a little bit.

It was a simple divide and conquer algorithm.

I hope you all remember it. And if we look at an array, and I'm not going to draw the cells here because I want to imagine a really big array, binary search, but suppose it always goes to left.

It starts by visiting this element in the middle.

Then ago so the quarter marked. Then it goes to the one eighth mark. OK, this is one hypothetical execution of a binary search. OK, and eventually it finds the element it's looking for. It finds where it fits at least. So x is over here.

So, we know that it takes log N time.

How many memory transfers of the take?

Now, I blocked this array into chunks of size B, blocks of size B. How many blocks do I touch?

This one's a little bit more subtle.

It depends on the relative sizes of N and B, yeah. Log base B of N would be a good guess. We would like it to be, let's say, hope, is that it's log base B of N because we know that B trees can search in what's essentially a sorted list of N items in log base B of N time.

That turns out to be optimal in the cache oblivious model or in the two level model you've got to pay log base B of N.

I won't prove that here. The same reason you need log N comparisons to do a binary search in the normal model.

Alas, it is possible to get log base B of N even without knowing B. But, binary search does not do it. Log of N over B, yes. So the number of memory transfers on N items is log of N over B also known as, let's say, plus one, also known as log N minus log B. OK, whereas log base B of N is log N divided by log B, OK, clearly this is much better than subtracting. So, this would be good, but this is bad. Most of the time, this is log N, which is no better, I mean, you're not using blocks at all essentially.

The idea is, out here, I mean, there's some little, tiny block that contains this thing. I mean, tiny depends on how big B is. But, each of these items will be in a different block until you get essentially within one block worth of x. When you get within one block worth of x, there's only like a constant number of blocks that matter, and so all of these accesses are indeed within the same block. But, how many are there?

Well, just log B because you're only spending log B within a, if you're within an interval of size k, you're only going to spend log k steps in it. So, you're saving log B in here, but overall you're paying log N, so you only get log N minus log B plus some constant. OK, so this is bad news for binary search. So, not all of the algorithms we've seen are going to work well in this model.

We need a lot more thinking before we can solve what is essentially the binary search problem, finding an element in a sorted list, in log base B of N without knowing B.

OK, we know we could use B trees.

If you knew B, great, that works, and that's optimal. But without knowing B, it's a little bit harder. And this gets us into the world of divide and conquer. Also like last week, and like the first few weeks of this class, divide and conquer is your friend. And, it turns out divide and conquer is not the only tool, but it's a really useful tool in designing cache oblivious algorithms.

And, let me say why.

So, we'll see a bunch of divide and conquer based algorithms, cache oblivious. And, the intuition is that we can take all the favorite algorithms we have, obviously it doesn't always work.

Binary search was a divide and conquer algorithm.

It's not so great. But, in general, the idea is that your algorithm can just do the normal divide and conquer thing, right?

You divide your problem into subproblems of smaller size repeatedly, all the way down to problems of constant size, OK, just like before. But, if you're recursively dividing your problem into smaller things, at some point you can think about it and say, well, wait, I mean, the algorithm divides all the way, but we can think about the point at which the problem fits in a block or fits in cache. OK, and that's the analysis.

OK, we'll think about the time when your problem is small enough that we can analyze it in some other way.

So, usually, we analyze it recursively.

We get a recurrence. What we're changing, essentially, is the base case.

So, in the base case, we don't want to go down to a constant size. That's too far.

I'll show you some examples. We want to consider the point in recursion at which either the problem fits in cache, so it has size less than or equal to M, or it fits in order one blocks. That's another natural time to do it. Order one blocks would be even better than fitting in cache. So, this means a size order B.

OK, this will change the base case of the recurrence, and it will turn out to give us good answers instead of bad ones. So, let's do a simple example.

Our good friend order statistics, in particular, for finding medians. So, I hope you all know this by heart. Remember the worst case linear time, median finding algorithm by Bloom et al.

I'll write this fast. We partition our array.

It turns out, this is a good algorithm as it is. We partition our array conceptually into N over five, five tuples into little groups of five. This may not have been exactly how I wrote it last time. I didn't check.

But, it's the same algorithm. You compute the median of each five tuple. Then you recursively compute the median of the medians of these medians.

Then, you partition around x. So, that gave us some element that was roughly in the middle. It was within the middle half, I think. Partition around x, and then we show that you could always recurse on just one of the sides.

OK, this was our good old friend for computing, order statistics, or medians, or whatnot.

OK, so how much time does this, well, we know how much time this takes. It should be linear time.

But how many memory transfers does this take?

Well, conceptually partitioning that, I can do, in zero. Maybe I have to compute N over five, no big deal here. We're not thinking about computation. I have to find the median of each tuple. So, here it matters how my array is laid out. But, what I'm going to do is take my array, take the first five elements, and then the next five elements and so on.

Those will be my five tuples. So, I can implement this just by scanning, and then computing the median on those five elements, which I stored in the five registers on my CPU.

I'll assume that there are enough registers for that.

And, I compute the median, write it out to some array out here. So, it's going to be one element. So, the median of here goes into there. The median of these guys goes into there, and so on. So, I'm scanning in here, and in parallel, I'm scanning an output in here.

So, it's two parallel scans. So, that takes linear time.

So, this takes order N over B plus one memory transfers.

OK, then we have recursively compute the median of the medians. This step used to be T of N over five. Now it's MT of N over five, OK, with the same values of B and M.

Then we partition around x. Partitioning is also like three parallel scans if you work it out.

So, this is also going to take linear memory transfers, N over B plus one. And then, we recurse on one of the sides, and this is the fun part of the analysis which I won't repeat here. But, we get MT of, like, three quarters N. I think originally it was seven tenths, so we simplified to three quarters, which is hopefully bigger than seven tenths.

Yeah, it is. OK, so this is the new analysis. Now we get a recurrence.

So, let's do that.

So, the analysis is we get this MT of N is MT of N over five plus MT of three quarters N plus, this is just as before.

Before we had linear work here. And now, we have what we call linear number of memory transfers, linear number of blocks. OK, I'll sort of ignore this plus one. It's not too critical.

So, this is our recurrence. Now, it depends what our base case is. And, usually we would use a base case of constant size. So, let's see what happens if we use a base case of constant size just so that it's clear why this base case is so important. OK, this describes a recurrence as one of these hairy recurrences.

And, I don't want to use substitution.

I just want the intuition of why this is going to solve to something rather big. OK, and for me, the best intuition always comes from recursion trees.

If you don't know the solution to recurrence and you need a good guess, use recursion trees. And today, I will only give you good guesses. I don't want to prove anything with substitution because I want to get to the bigger ideas.

So, this is even messy from a recursion tree point of view because you have these unbalanced sizes where you start at the root with some of size N over B.

Then you split it into something size one fifth N over B, and something of size three quarters N over B, which is annoying because now this subtree will be a lot bigger than this one, or this one will terminate faster. So, it's pretty unbalanced.

But, summing per level doesn't really tell you a lot at this point. But let's just look at the bottom level. Look at all the leaves in this recursion tree. So, that's the base cases.

How many base cases are there? This is an interesting question. We've never thought about it in the context of this recurrence. It gives a somewhat surprising answer. It was surprising to me the first time I worked it out. So, how many leaves does this recursion tree have? Well, we can write a recurrence. The number of leaves in a problem of size N, it's going to be the number of leaves in this problem plus the number of leaves in this problem plus zero. So, that's another recurrence.

We'll call this L of N. OK, now the base case is really relevant. It determines the solution to this recurrence. And let's, again, assume that in a problem of size one, we have one leaf.

That's our only base case. Well, it turns out, and here you need to guess, I think.

This is not particularly obvious.

Any of the TA's have guesses of the form of this solution?

Or anybody, not just TA's. But this is open to everyone.

If Charles were here, I would ask him.

I had to think for a while, and it's not linear, right, because you're somehow decreasing quite a bit.

So, it's smaller than linear, but it's more than a constant.

OK, it's actually more than polylog, so what's your favorite function in the middle? N over log N, that's still too big. Keep going.

You have an oracle here, so you can, N to the k, yeah, close. I mean, k is usually an integer. N to the alpha for some real number between zero and one. Yeah, that's what you meant.

Sorry. It's like the shortest mathematical joke. Let epsilon be less than zero or for a sufficiently large epsilon.

I don't know. So, you've got to use the right letters. So, let's suppose that it's N to the alpha. Then we would get this N over five to the alpha, and we'd get three quarters N to the alpha. When you have a nice recurrence like this, you can just try plugging in a guess and see whether it works, OK, and of course this will work only depending on alpha. So, we should get an equation on alpha here. So, everything has an N to the alpha, in fact, all of these terms.

So, I can divide through my N to the alpha.

That's assuming that it's not zero or something.

That seems reasonable. So, we have one equals one fifth to the alpha plus three quarters to the alpha.

This is something you won't get on a final because I don't know any good way to solve this except with like Maple or Mathematica. If you're smart I'm sure you could compute it in a nicer way, but alpha is about 0.8, it turns out. So, the number of leaves is this sort of in between constant and linear.

Usually polynomial means you have an integer power.

Let's call it a polynomial. Why not?

There's a lot of leaves, is the point, and if we say that each leaf costs a constant number of memory transfers, we're in trouble because then the number of memory transfers has to be at least this.

If it's at least that, that's potentially bigger than N over B, I mean, bigger than in an asymptotic sense. This is little omega of N over B if B is big. If B is at least N to the 0.2 something, OK, or one seventh something.

But if, in particular, B is at least N to the 0.2, then this should be bigger than that.

So, this is a bad analysis because we're not going to get the answer we want, which is N over B.

The best you can do for median is N over B because you have to read all the element, and you should spend linear time. So, we want to get N over B.

This algorithm is N over B plus one.

So, this is why you need a good base case, all right?

So that makes the point. So, the question is, what base case should I use?

So, we have this recurrence What base case should I use? Constant was too small.

We have a couple of choices listed up here.

Any suggestions? B, OK, MT of B is?

The hard part. So, if my problem, if the size of my array fits in a block and I do all this stuff on it, how many memory transfers could that take?

One, or a constant, depending on alignment.

OK, maybe it takes two memory transfers, but constant.

Good. That's clearly a lot better than this base case, MT of one equals order one, clearly stronger. So, hopefully, it gives the right answer, and now indeed it does.

I love this analysis. So, I'm going to wave my hands.

OK, but in particular, what this gives us, if we do the previous analysis, what is the number of leaves?

So, in the leaves, now L of B equals one instead of L of one equals one. So, this stops earlier.

When does it stop? Well, instead of getting N to the order of 0.8, whatever, we get N over B to the power of 0.8 whatever. OK, so it turns out the number of leaves is N over B to the alpha, which is little o of N over B. So, we don't care.

It's tiny. And, if you look at the root cost is N over B in the recursion tree, the leaf cost is little o of N over B, and if you wave your hands, and close your eyes, and squint, the cost should be geometrically decreasing as we go down, I hope, more or less. It's a bit messy because of all the things terminating, but let's say cost is roughly geometric. Don't do this in the final, but you won't have any messy recurrences like this.

So, don't worry. Down the tree, so you'd have to prove this formally, but I claim that the root cost dominates. And, the root cost is N over B.

So, we get N over B. OK, so this is a nice, linear time algorithm for order statistics for cache oblivious.

Great. This may turn you off a little bit, but even though this is like the simplest algorithm, it's also probably the most complicated analysis that we will do. In the future, our algorithms will be more complicated, and the analyses will be relatively simple. And usually, it's that way with cache oblivious algorithms.

So, I'm giving you this sort of as the intuition of why this should be enough. Then you have to prove it.

OK, let's go to another problem where divide and conquer is useful, our good friend, matrix multiplication.

I don't know how many times we've seen this in this class, but in particular we saw it last week with a recursive matrix multiply, multithreaded algorithm.

So, I won't give you the algorithm yet again, but we're going to analyze it in a very different way.

So, we have C and we have A, and actually up to you.

So, I could cover standard matrix multiplication, which is when you do it row by row, and column by column.

And, we could see why that's bad.

And then, we could do the recursive one and see why that's good. Or, we could skip the standard algorithm. So, how many people would like to see why the standard algorithm is bad?

Because it's not totally obvious.

One, two, three, four, five, half?

Wow, that's a lot of votes. Now, how many people want to skip to the chase? No one.

One, OK. And, everyone else is asleep.

So, that's pretty good, 50% awake, not bad.

OK, then, so standard matrix multiplication.

I'll do this fast because it is, I mean, you all know the algorithm, right? To compute this value of C; in A, you take this row, and in B you take this column.

Sorry I did a little bit sloppily.

But this is supposed to be aligned.

Right? So I take all of this stuff, I multiply it with all of the stuff, add them up, the dot product. That gives me this element.

And, let's say I do them in this order row by row.

So for every item in C, I loop over this row and this column, B, multiply them together.

That is an access pattern in memory.

So, exactly how much that costs depends how these matrices are laid out in memory. OK, this is a subtlety we haven't had to worry about before because everything was uniform. I'm going to assume to give the standard algorithm the best chances of being good, I'm going to store C in row major order, A in row major order, and B in column major order.

So, everything is nice and you're scanning.

So then this inner product is a scan.

Cool. Sounds great, doesn't it? It's bad, though.

Assume A is row major, and B is column major.

And C, you could assume is really either way, but if I'm doing it row by row, I'll assume it's row major.

So, this is what I call the layout, the memory layout, of these matrices. OK, it's good for this algorithm, but the algorithm is not good.

So, it won't be that great.

So, how long does this take? How many memory transfers?

We know it takes M^3 time. Not going to try and beat M^3 here. Just going to try and get standard matrix multiplication going faster.

So, well, for each item over here I pay N over B to do the scans and get the inner product. So, N over B per item.

So, it's N over B, or we could go with the plus one here, to compute each c_ij. So that would suggest, as an upper bound at least, it's N^3 over B.

OK, and indeed that is the right bound, so theta.

This is memory transfers, not time, obviously.

That is indeed the case because if you look at consecutive, I do this c_ij, then this one, this one, this one, this one, keep incrementing j and keeping I fixed, right?

So, the row that I use stays fixed for a long time.

I get to reuse that if it happens, say that that fits a block maybe, I get to reuse that row several times if that happens to fit in cache. But the column is changing every single time. OK, so every time I moved here and compute the next c_ij, even if a column could fit in cache, I can't fit all the columns in cache.

And the columns that I'm visiting move, you know, they just scan across.

So, I'm scanning this whole matrix every time.

And unless you're entire matrix fits in cache, in which case you could do anything, I don't care, it will take constant time, or you'll take M over B time, enough to read it into the cache, do your stuff, and write it back out. Except in that boring case, you're going to have to pay N^2 over B for every row here because you have to scan the whole collection of columns.

You have to read this entire matrix for every row over here.

So, you really do need N^3 over B for the whole thing.

So, it's usually a theta. So, you might say, well, that's great. It's the size of my problem, the usual running time, divided by B.

And that was the case when we are thinking about linear time, N versus N over B. It's hard to beat N over B when your problem is of size N. But now we have a cubed.

And, this gets back to, we have good spatial locality.

When we read a block, we use the whole thing.

Great. It seems optimal.

But we don't have good temporal locality.

It could be that maybe if we stored the right things, we kept them around, we could them several times because we're using each element like a cubed number of times.

That's not the right way of saying it, but we're reusing the matrices a lot, reusing those items.

If we are doing N^3 work on N^2 things, we're reusing a lot.

So, we want to do better than this, and that's the recursive algorithm, which we've seen. So, we know the algorithm pretty much. I just have to tell you what the layout is. So, we're going to take C, partition of C_1-1, C_1-2, and so on.

So, I have an N by N matrix, and I'm partitioning into N over 2 by N over 2 submatrices, all three of them times whatever. And, I could write this out yet again but I won't. OK, we can recursively compute this thing with eight matrix multiplies, and a bunch of matrix additions. I don't care how many, but a constant number. We see that at least twice now, so I won't show it again. Now, how do I lay out the matrices? Any suggestions how I lay out the matrices? I could lay them out in row major order. I'll call it major order.

But that might be less natural now.

We're not doing anything by rows or by columns.

So, what layout should I use? Yeah?

Quartet major order, maybe quadrant major order unless you're musically inclined, yeah.

Good idea. You've never seen this order before, so it's maybe not so natural.

Somehow I want to cluster it by blocks.

OK, I think that's about all. So, I mean, it's a recursive layout. This was not an easy question.

It's OK. Store matrices or lay out the matrices recursively by block. OK, I'm cheating a little bit.

I'm redefining the problem to say, assume that your matrices are laid out in this way. But, it doesn't really matter.

We can cheat, can't we?

In fact, it doesn't matter. You can turn a matrix into this layout without too much linear work, almost linear work.

Log factors, maybe.

OK, so if I want to store my matrix A as a linear thing, I'm going to recursively defined that layout to be recursively store the upper left corner, then store, let's say, the upper right corner.

It doesn't matter which order I do these.

I should have drawn this wider, then store the lower left corner, and then store the lower right corner recursively.

So, how do you store this? Well, you divide it in four, and lay out the top left, and so on.

OK, this is a recursive definition of how the element should be stored in a linear array.

It's a weird one, but this is a very powerful idea in cache oblivious algorithms.

We'll use this multiple times. OK, so now all we have to do is analyze the number of memory transfers.

How hard could it be? So, we're going to store all the matrices in this order, and we want to compute the number of memory transfers on an N by N matrix.

See, I lapsed and I switched to lowercase n.

I should, throughout this week, be using uppercase N because for historical reasons, any external memory kinds of algorithms, to level algorithms, always talk about capital N.

And, don't ask why. You should see what they define little n to be. OK, so, any suggestions on what the recurrence should be now? All his fancy setup with the recurrence is actually pretty easy.

So, definitely it involves multiplying matrices that are N over 2 by N over 2. So, what goes here?

Eight, thank you. That you should know.

And that the tricky part is what goes here.

OK, what goes here is, now, the fact that I can even write this, this is the matrix additions.

Ignore those for now. Suppose there weren't any.

I just have to recursively multiply.

The fact that this actually is eight times memory transfers of N over 2 relies on this layout. Right, I'm assuming that the arrays that I'm given are given as contiguous intervals and memory. If they aren't, I mean, if they're scattered all over memory, I'm screwed. There's nothing I can do.

So, but by assuming that I have this recursive layout, I know that the recursive multiplies will always deal with three consecutive chunks of memory, one for A, one for B, one for C, OK, no matter what I do.

Because these are stored consecutively, recursively I have that invariant.

And I can keep recursing. And I'm always dealing with three consecutive chunks of memory.

That's why I need this layout is to be able to say this.

OK, Now what does addition cost?

I'll just give you two matrices.

They're stored in some linear order, the same linear order among the three of them. Do I care what the linear order is? How should I add two matrices, get the output?

Yeah?

Right, if each of the three arrays I'm dealing with are stored in the same order, I can just scan in parallel through all three of them and just add corresponding elements, and output it to the third. So, I don't care what the order is, as long as it's consistent and I get N^2 over B.

I'll ignore plus one here. That's just looking at the entire matrix. So, there we go: another recurrence. We've seen this with N^2, and we just got N^3. But, it turns out now we get something cooler if we use the right base case.

So now we get to the base case, ah, the tricky part.

So, any suggestions what base case I should use?

The block size, good suggestion.

So, if we have something of size order B, we know that takes a constant number of memory transfers.

It turns out that's not enough. That won't solve it here.

But good guess. In this case, it's not the right answer. I'll give you some intuition why. We are trying to improve on N^3 over B. If you were just trying to get it divided by B, this is a great base case.

But here, we know that just the improvement afforded by the block size is not enough. We have to somehow use the fact that the cache is big. It's M, so however big M is, it's that big. OK, so if we want to get some improvement on this, we've got to have M in the formula somewhere, and there's no M's yet.

So, it's got to involve M. What's that?

MT of M over B? That would work, but MT of M is also OK, I mean, some constant times M, let's say. I want to make this constant small enough so that the entire problem fits in cache.

So, it's like one third. I think it's actually, oh wait, is it the square root of M actually?

Right, this is an N by N matrix.

So, it should be C times the square root of M.

Sorry. So, the square root of M by square root of M matrix has M entries.

If I make C like one third or something, then I can fit all three matrices in memory. Actually, one over square root of three would do, but who cares?

So, for some constant, C, now everything fits in memory. How many memory transfers does it take? One?

It's a bit too small, because I do have to read the problem in. And now, I mean, here was one because there's only one block to read.

Now how many blocks are there to read?

Constants? No.

B? No.

M over B, good. Get it right eventually.

That's the great thing about thinking with an oracle.

You can just keep guessing. M over B because we have cache size M. There are M over B blocks in that cache to read each one, OK?

This is maybe, you forgot what M was because we haven't used it for a long time.

But M is the number of elements in cache.

This is the number of blocks in cache.

OK, some of was saying B, and it's reasonable to assume that M over B is about B. That's like a square cache, but in general, we don't make that assumption.

OK, where are we? We're hopefully done, just about, good, because we have three minutes.

So, that's our base case. I have a square root here; I just forgot it. Now we just have to solve it.

Now, this is an easier recurrence, right?

I don't want to use the master method, because master method is not going to handle these B's and M's, and these crazy base cases. OK, master method would prove N^3. Great.

Master method doesn't really think about these kinds of cases. But with regression trees, if you remember way back to the proof of the master method, just look at the recursion tree as geometric up or down where everything is equal, and then you just add them up, every level. The point is that this is a nice recurrence. All of the sub problems are the same size, and that analysis always works, I say, when everything has the same size, all the children.

So, here's the recursion tree. We have N^2 over B at the top.

We split into eight subproblems where each one, the cost is one half N^2 over B.

I'm not going to write them all.

There they are. You add them up.

How much do you get? Well, there's eight of them.

Eight times a half is two. Four.

[LAUGHTER] Thanks. Four, right?

OK, I'm bad at arithmetic. I probably already said it, but there are three kinds of mathematicians, those who can add, and those who can't.

OK, why am I looking at this? It's obvious.

OK, so we keep going. This looks geometrically increasing. Right?

You just know in your heart that if you work out the first two levels, you can tell whether it's geometrically increasing, decreasing, or they're all equal, or something else.

And then you better think. But I see this as geometrically increasing. It will indeed be like 16 at the next level, I guess.

OK, it should be. So, it's increasing.

That means the leaves matter. So, let's work out the leaves.

And, this is where we use our base case.

So, we have a problem of size square root of M.

And so, yeah, you have a question?

Oh, indeed. I knew there was something.

I knew it was supposed to be two out here.

Thanks. This is why you're here.

It's actually N over two squared over B.

Thanks. I'm substituting N over 2 into this. OK, so this is actually N^2 over 4 B. So, I get two, because there are eight times one over four.

OK, I wasn't that far off then. It's still geometrically increasing, still the case, OK?

But now, it actually doesn't matter.

Whatever the cost is, as long as it's bigger than one, great. Now we look at the leaves.

The leaves are root M by root M.

I substitute root M into this: I get M over B with some constants. Who cares?

So, each leaf is M over B, OK, lots of them.

How many are there? This is the only, deal with recursion trees, counting the number of leaves is always the annoying part. Oh boy, well, we start with an N by N matrix. We stop when we get down to root N by root N matrix. So, that sounds like something.

Oh boy, I'm cheating here. Really?

That many? It sounds plausible.

OK, the claim is, and I'll cheat.

So I'm going to use the oracle here, and we'll figure out why this is the case. N over root N^3 leaves, hey what? I think here, it's hard to see the tree. But it's easy to see in the matrix. Let's enter the matrix.

We have our big matrix. We divided in half.

We recursively divide in half. We recursively divide in half.

You get the idea, OK?

Now, at some point these sectors, let's say one of these sectors, and each of these sectors, fits in cache.

And three of them fit in cache. So, that's when we stop the recursion in the analysis. The algorithm goes all the way.

But in the analysis, let's say we stop at M.

OK, now, how many leaves or problems are there?

Oh man, this is still not obvious.

OK, the number of leaf chunks here is, like, I mean, the number of these things is something like N over root M, right, the number of chunks.

But, it's a little less clear because I have so many of these.

But, all right, so let's just suppose, now, I think of normal, boring, matrix multiplication on chunks of this size. That's essentially what the leaves should tell me. I start with this big problem, I recurse out to all these little, tiny, multiply this by that, OK, this root M by root M chunk. OK, how many operations, how many multiplies do I do on those things?

N^3. But now, N, the size of my matrix in terms of these little sub matrices, is N over root M. So, it should be N over root M^3 subproblems of this size. If you work it out, normally we go down to things of constant size and we get exactly N^3 of them. Now we are stopping at this short point in saying, well, it's however many there are, cubed. OK, this is a bit of hand waving. You could work it out with the recurrence on the number of leaves.

But there it is. So, the total here is N over, let's work it out. N^3 over M to the three halves, that's this number of leaves, times the cost at each leaf, which is M over B. So, some of the N's cancel, and we get N^3 over B root M, which is a root M factor better than N^3 over B. It's actually quite a lot, the square root of the cache size.

That is optimal. The best two level matrix multiplication algorithm is N^3 over B root M memory transfers.

Pretty amazing, and I'm over time.

You can generalize this into all sorts of great things, but the bottom line is this is a great way to do matrix multiplication as a recursion. We'll see more recursion for cache oblivious algorithms on Wednesday.

Free Downloads

Free Streaming

• VideoLectures.net

Video

• iTunes U (MP4 - 178MB)
• Internet Archive (MP4 - 349MB)

Caption

• English-US (SRT)