Lecture 9: Cache-Efficient Algorithms II

• Download English-US transcript (PDF)

The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high-quality educational resources for free. To make a donation or view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.

JOHN DONG: All right, so I'm sure everybody's kind of curious about Project 2.2 Beta, so here's the preliminary performance results. There's still a bit more work to do on finalizing these numbers, but here's how they look. But before I show you guys that, I'd like to yell at you guys for a little bit. So this is a timeline of the submission deadline and when people submitted things. So from this zone to this zone is an hour, and I see like 50% of the commits made during that period.

So a little bit of background, the submission checking system automatically clones a repository up to the deadline and not a second after. So if you look at from this to this, some people took quite a jump back. So just as a warning, please try to submit things on time. 11:59 means 11:59, and to drive that point further, this is an example of a commit that I saw in somebody's repository, who's name I blanked out. Obviously, it's seven seconds past the deadline, so the automatic repository cloner didn't grab it. And the previous commit to that was like 10 days ago when Reid pushed out pentominoes grades. So don't do that either.

A little break down on how often people commit. About a quarter of the class only made one commit to their repository. Half of you guys did three to 10 commits, which seems to be an up. Under 21+, there was somebody who did 100 some commits, which was pretty impressive. Yeah, very smart dude. Committing often is a good idea, so that you don't run into a situation like before. And next time definitely, there is going to be next to zero tolerance for people who don't commit things on time for the deadlines.

Now, the numbers that people want. So for rotate, just as an interesting data point, this is the 512 by 512 case. It seems like not everybody remembered to carry over their optimizations for the 512 case over to the final, otherwise you would expect the numbers to be a little bit more similar and not off by a factor of eight. And for rotate overall, that's the distribution. The speed up factor is normalized to some constant that gives everybody a reasonable number. And there were a lot of groups that had code that didn't build. Yes?

AUDIENCE: Does the speedup only include the rotate dot 64.

JOHN DONG: Yes. Performance was only tested on the rotate dot 64, which is what we said in the handout also. So there were a lot of groups that didn't build, which really surprised me, and I think it's probably because half of you guys pushed things after the deadline, and I presume those things contained important commits towards making your code actually work. But in this case, there was no header files, there was no cross-testing. All we did was run your make file on your code, and we replaced your test bed dot c, and your k timing, so I'm not quite sure why people have code that doesn't build.

For sort, this was the maximum size input that we allowed you guys to run. One group did really well. So all of these are correct. The correctness test is built in, and I replaced that with a clean copy that contains a couple additional checks, by the way. So I tried the other extreme, a relatively small array.

AUDIENCE: Is the top the same person?

JOHN DONG: I'm not sure whether or not the top is the same person. But the distribution wise seems like people didn't quite remember to optimize for the smallest case. And this is the current overall speedup factors for sort averaging. We did about 10 to 15 cases for rotate and sort. And that's the overall speedup distribution.

SAMAN AMARASINGHE: OK, so now you're done with individual project, so you already did the last project individually, and then we are moving into, again, a group project. So the first thing we have, setting up automated system for you to say who your group members are. So we will send you information, and with that, what you have to do is run a script saying who your group members are, both group members have to do it, and then we will basically clear that account in there.

That said, a lot of you didn't know, in the first project, how to work and what's the right mode of operations with the group. OK if we gave you to write 100,000 lines of code, it makes sense to say, OK, I'm going to divide the problem into half, one person do one, the other person do the other half. The reason for doing the group is to try to get you to do pair programming, because talking to a lot of you, getting a lot of feedback, it looks like most of you spent a huge amount of time debugging. And since you're only writing a little amount of code, it makes a lot more sense to sit with your partner next to the screen, one person types, other person looks over, and then you have a much faster way of getting through the debugging process. So the next one, don't try to divide the problem by half. Just try to find some time, sit with each other.

Then the other really disturbing thing is there have been a couple of groups that were completely dysfunctional. We get emails saying, OK, my group member didn't talk to me, or they didn't do any work, or they were very condescending. And that's really sad, because from my experience with MIT students, when you guys go to company, you guys probably will be the best programmers there. There's no question about it. I have seen that. To the point that some people might even resent it, to have this best programmer.

But what I have seen is the fact a lot of you cannot work in the group, if you haven't developed that skill, you will not be the most impactful person. I have seen that again and again in my experience with doing a start-up. Our MIT students, their way of impacting is put all night do the entire project by themselves. Doable when you're doing a small change to a large project, but if you want to have a big change, you can't do that. You have to work with the group, figure out how to impact, how to communicate. This is more important learning than, say, trying to figure out how you can optimize something.

So being individual contributors and able to do amazing things is important, but the fact that you can't work with a group is going to really make the impact you can make much less. So please, please learn how to work with your group members. Some of them might not be as good as you are, and that will be probably true in real life, too, but that doesn't mean you can be condescending towards them, make them feel inferior. That doesn't cut it. You have to learn how to work with these people. So part of your learning is working with others, and that's a large part of your learning. Don't consider that to be this external thing, even though you might think you can do a better job.

Just work with the other person, especially in pair programming, because where both of you are sticking together, it's much easier. Because there are four eyes on your project, they might see something you don't see, and see whether you can work together like that. And I don't want to hear any more stories saying, look, my partner was too dumb, or my partner didn't show up, he couldn't deal with that. Those are not great excuses. And so we wonder. Some of them, we will pay attention to it, because it might be one person's unilateral actions that lead to that. Still, please, try to figure out how we can work with your partners. And hope you have a good partner experience. Use pair programming. Use a lot of good debugging techniques, and the next project will be fine.

CHARLES LEISERSON: Great. We're going to talk more about caches. Whoo-whoo! OK. So for those who weren't here last time, we talked about the ideal cache model. As you recall, it has a two-level hierarchy and a cache size M bytes, and a cache line length of B bytes. It's fully associative and is optimal omniscient replacement strategy. However, we also learned that LRU was a good substitute, and that any of the asymptotic results that you can get with optimal, you could also get with LRU. The two performance measures that we talked about were the work, which deals with what the processor ends up doing, and the cache misses, which is the transfers between cache and main memory. You only have to count one direction, because what goes in basically goes out, so more or less, it's the same number.

OK, so I'd like to start today by talking about some very basic algorithms that you have seen in your algorithms and data structures class, but which may look new when we start taking caches into account. So the first one here is the problem of merging two sorted arrays. So as you recall, you can basically do this in linear time. The way that the algorithm works is that it looks at the first element of the two arrays to be sorted and whichever is smaller, it puts in the output. And then it advances the pointer to the next element, and then whichever is smaller, it puts in the output. And in every step, it's doing just a constant amount of work, there are n items, so by the time this process is done, we basically spent time proportional to the number of items in the output list. So this should be fairly familiar. So the time to emerge n elements is order n.

Now, the reason merging is useful is because you can use it in a sorting algorithm, a merge sorting algorithm. The way merge sort works is it essentially does divide and conquer on the array. It divides the array into two pieces, and it divides those each into two pieces, and those into two, until it gets down to something of unit size. And then what it does is it merges the two pairs of arrays. So for example here, the 19 and 3 got merged together to be 3 and 19. The 12 and 46 were already in order, but you still had to do work to get it there and so forth. So it puts everything in order in pairs, and then for each of those, it puts it together into fours, and for each of those, it puts it together into the final list. Now of course, the way it does this is not in the order I showed you. It actually goes down and does a walk of this tree. But conceptually, you can see that it essentially comes down to merging pairs, merging quadruples, emerging octuples, and so forth, all the way until the program is done.

So to calculate the work of Merge Sort, this is something you've seen before because it's exactly what you do in your algorithms class. You get a recurrence that says that the work, in this case is-- well, if you have only one element, it's a constant amount of work, and otherwise, I solve two problems of half the size doing order n work, which is the time to merge the two elements. So classic divide and conquer. And I'm sure you're familiar with what the solution to this recurrence is. What's the solution to this recurrence? n log n. I want to, though, step through it, just to get everybody warmed up to the way that I want to solve recurrences so that you are in a position, when we do the caching analysis-- we have a common framework for understanding how the caching analysis will work.

So we're going to solve this recurrence, and if the base case is constant, we usually omit it. It's assumed. So we start out with W of n, and what we do is we replace it by the right hand side, where we put the constant term on the top and then the two children. So here I've gotten rid of the theta, because conceptually when I'm done, I can put a big theta around the whole thing, around the whole tree, and it just makes the math a little bit easier and a little bit clearer to see what's going on. So then I take each of those and I split those, and this time I've got n over 4. Correct? I checked for that one this time. It's funny because it was still there, actually, just a few minutes before class as I was going through. And we keep doing that until we get down to something of size one, until the recurrence bottoms out.

So when you look at a recursion tree of this nature, the first thing that you typically want to do is take a look at what the height of the tree is. In this case, we're taking a problem of size n, and we're halving it at every step. And so the number of times we have to halve the argument-- which turns out to be also equal to the work, but that's just coincidence-- is log n times. So the height is log base 2 of n. Now what we do typically is we add things up across the rows, across the levels. So on the top level, we have n. On the next level, we have n. The next level, hey, n. To add up the bottom, just to make sure, we have to count up how many leaves there are, and the number of leaves, since this is a binary tree, is just 2 to the height. So it's 2 to the log n, which is n. So then I add across all the leaves, and I get the order 1 at the base times n, which is order n.

And so now I'm in a position to add up the total work, which is basically log n levels of n is total of order n log n. So hopefully, this is all review. Hopefully this is all review. You haven't seen this before. It's really neat, isn't it? But you've missed something along the way.

So now with caching. So the first thing to observe is that merge subroutine, the number of cache misses that it has is order n over B. So as you're going through, these arrays are laid out continuously in memory. The number of misses-- you're just going through the data once-- is order n data, all going through contiguously. And so every time you bring in data, you get full spatial locality, there are n elements, it costs n over B. So is that plain? Hopefully that part's plain, because you bring in things of each axis, you get the same factor of B.

So now merge sort-- and this is, once again, the hard part is coming up with a recurrence, and then the other hard part is solving it. So there's two hard parts to recurrences. OK, so the merge sort algorithm solves two problems of size n over 2, and then does a merge. So the second line here is pretty straightforward. I take the cache misses that I need to do, and I take a merge. I may have a few other accesses in there, but they're going to be dominated by the merge. OK, so it's still going to be theta n over B.

Now the hard part, generally, of dealing with cache analysis is how to do the base case, because the base case is more complicated than when you just do running time, you get that to run down to a base case of constant size. Here, you don't get to run down to base case of constant size. So what it says here is that we're going to run down until I have a sorting problem that fits in cache, n is going to be less than some constant times n, for some sufficiently small constant c less than 1. When in finally fits in cache, how many cache misses does it take me to sort it? Well, I only need to have the cold misses to bring that array into memory, and so that's just proportional to n over B, because for all the rest of the levels of merging, you're inside the cache. That make sense? So that's where we get this recurrence. This is always a tricky thing to figure out how to write that recurrence for a given thing.

Then, as I say, the other tricky thing is how do you solve it? But we're going to solve it essentially the same way as we did before. I'm not going to go through all the steps here, except to just elaborate. So what we're doing is we're taking, and we have n over B at the top, and then we divide it into two problems, and for each of those-- whoops, there's a c there that doesn't belong. It should just be n over 2B on both, and then n over 4B, and so forth. And we keep going down until we get to our base case.

Now in our base case, what I claim is that when I hit this base case, it's going to be the case that n is, in fact, a constant factor times n, so that n over B is almost the same as m over B. And the reason is because before I hit the base case, when I was at size twice n, and that was bigger than m. So if twice n is bigger than m-- than my constant times m-- but n is smaller than m, then it's the case that n and m are essentially the same size to within a constant factor, within a factor of two, in fact. And so therefore, here I can say that it's order m over B.

And now the question is, how many levels did I have to go down cutting things in half before I got to something of size m over B? Well, the way that I usually think about this is-- you can do it by taking the difference as I did before. The height of the whole tree is going to be log base 2 of n, and the height of this is basically log of-- what's the size of n? It's going to be log of the size of n when this occurs. Well, n at that point is something like cm. So it's basically log n minus log cm, which is basically log of n over cm. How about some questions. Yeah, question.

AUDIENCE: --the reason for why you just substituted on the left, [INAUDIBLE] over B, but on the right [INAUDIBLE].

CHARLES LEISERSON: Here?

AUDIENCE: No.

CHARLES LEISERSON: Or you mean here?

AUDIENCE: [INAUDIBLE].

CHARLES LEISERSON: On the right side.

AUDIENCE: On the right-most leaf, it's n over B. Is that all of the leaves added up because [INAUDIBLE]?

CHARLES LEISERSON: No, no, no, this is going to be all of the leaves added up here. This is the stack I have on the right hand side. So we'll get there. So the point is, the number of leaves is 2 to this-- so that's just n over cm-- times m over B. Well, n over cm times m over B, the m's cancel, and I get essentially n over b with whatever that constant is here. And so now I have n over b across every level, and then when I add those up, I have to go log of n over cm levels, which is the same as log of n over m. So I have n over B times log of n over m. Yeah, question?

AUDIENCE: Initial assumption that c is some sufficiently small number, so 1 over c would be a rather large factor.

CHARLES LEISERSON: It could potentially be a large factor, but it's a constant. In other words, it can't vary with n. So in fact, typically for something like merge sort, the constant is going to be very close to-- for most things, the constant there is typically only a few. Because the question is, how many other things do you need in order to make sure it fits in? In this case, you have n. Well, you have both the input and the output, so here, it's going to be, you have to fit both the input and the output into cache in order not to have the cache it. So it's basically going to be like a factor of 2 for merge sort. For the matrix multiplication, it was like a factor of three. So generally a fairly small number. Question?

AUDIENCE: I guess that makes sense. But I [INAUDIBLE]. So on the order of the size of the leaves, can you assume that n is more than cm, so you can substitute--

CHARLES LEISERSON: Yeah, because basically, when it hits this condition--

AUDIENCE: Right. I understand that. Then you're also saying, why isn't there more or less just one B, because there's n over cm, n is the same as cm. That's [INAUDIBLE]. At the bottom level, there should be--

CHARLES LEISERSON: Oh, did I do something wrong here? Number of leaves is--

AUDIENCE: [INAUDIBLE].

CHARLES LEISERSON: Right, right. Sorry. This is the n at the top. You always have to be careful here. So this is the n in this case, so this is some n, little n. So this is not the same n. This is the n that we had at the top. This is the notion of recurrences, like the n keeps recurring, but you know which ones-- so it can be confusing, because if you're-- so yeah. So this is not the n that's here. This is the n that started out at the top. So we're analyzing it in terms of the n. Some people write these things where they would write this in terms of k, and then analyze it for n, and for some people, that can be helpful to do, to disambiguate the two things. I always find it wastes a variable, and you know, those variables are hard to come by. There's only a finite number of them.

OK, so are we good for this? So here, we ended up with n over b, log of n over m cache misses. So how does that compare? Let's just do a little thinking about this. Here's the recurrence, and I solved it out to this. So let's just look to see what this means. If I have a really big n, much bigger than the size of my cache, then I'm going to do a factor of B log n less misses than work, because the work is n over B log n. So if n compared to m-- let's say n was m squared or something-- then n over m would still be n, if n was as big as m squared. So if n was as big as m squared, this log of n over m would still be log n. And so I would basically have n over B log n for a factor of B log n less misses than work. If they're about the same-- did I get this right? If they're about the same size, if n is approximately m, maybe it's just a little bit bigger, then the log here disappears completely, and so I basically just have n over B misses.

AUDIENCE: [INAUDIBLE].

CHARLES LEISERSON: Yeah, but if n is like--

AUDIENCE: [INAUDIBLE].

CHARLES LEISERSON: Yeah. In fact, for this, you have to be careful as you get to the base cases. Technically, for some of this, I should be saying 1 plus log of n over m, and in some of the things I do later, I will put in the ones. But if you're looking at it asymptotically and n gets big, you don't have to worry about those cases. That just handles the case whether you're looking at n getting large or whether you're trying to handle a formula for all n, even if and n is small. Question?

AUDIENCE: [INAUDIBLE]?

CHARLES LEISERSON: The work was n log n, yes. The work was n log n. So here we basically have n over B log n, so I'm saving a factor of B in the case where it's about the same. Did I get this right? I'm just looking at this, and now I'm trying to reverse engineer what my argument is. So we're looking at n log n versus n over B log n over m.

AUDIENCE: [INAUDIBLE PHRASE]. So that you get a factor of B less misses, because you would be getting n over B like log of n, that's the only way you're getting a factor of B less misses. So I don't understand how you're saying, for n more or less equal to m. You would want something more like, for n--

CHARLES LEISERSON: Well, if n and m are about the same size, the number of cache misses is just n over B. And the number of cache misses is n over B, and the work is n log n. So I've saved a factor of B times log n, OK? What did I say?

AUDIENCE: [INAUDIBLE].

CHARLES LEISERSON: B log m? No, I was saying that's for the case when n is much bigger than m. So let's take a look at the case-- let me just do it on the board here. Let's suppose that n is like m squared, just as an example, big number. So I'm going to look at, essentially, n over B times log of n over m. So log of n over m, so what is n over m is about m, which is about square root of n, right? So this basically ends up being approximately n over B log of square root of n, which is the same as log n, to within a constant factor. I'm going to leave out the constant factors here.

Then I want to compare that with n log n. So I get a factor of B less misses. So the first one, yes, OK. So I get a factor of B less misses, you're right. Then I get a factor of B less misses. So I think I've got these switched. So this is the case I'm doing is for n much bigger than m. So let's do the other case. I think I've got the two things switched. I'll fix it in the notes.

If n and m are approximately the same, then the log is a constant, right? So this ends up being approximately n over B. And now when I take a look at the difference between the number of things, I get B log n. So I've got the two things mixed. Yeah?

AUDIENCE: As n approaches m, then the log would approach zero, but were you talking about how it technically should be--

CHARLES LEISERSON: 1 plus n, yes.

AUDIENCE: So technically, that approaches one when the log approaches zero.

CHARLES LEISERSON: Yeah.

AUDIENCE: These things are really hard for me, because they are really arbitrary. And then you're like, oh yeah, you can just put a 1 on top of there. And for example, I always miss those, because I usually try to do the math as rigorously as I can, and those ones generally do not appear, and you're like, oh, sure whatever. So how am I supposed to know that the log is actually not going to be zero, and I'm going to be like, yeah, you're not going to do any caches.

CHARLES LEISERSON: Because generally, what we're doing is we're looking at how things scale, so we're generally looking at n being big, in which case it doesn't matter. These things only matter if n's-- for example, notice here that if n goes less than m, we're in real trouble, right? Because now the log is negative. Wait, what does that mean? Well, the answer is the analysis was assuming that n was sufficiently large compared with m.

AUDIENCE: Why can't you just be like, oh, when n is for less than one, you can assume, well, n is 2n. In that case, you get log of two, which is still something or other.

CHARLES LEISERSON: Yeah, exactly. So what happens in these things is if you get right on the cusp of fitting in memory, then these analyses like, well, what exactly is the answer, is dicey. But if you assume that it doesn't fit in, what's going to happen? Or that does fit in, what is going to happen? And then the analysis right on the edge is somewhere between there. Good. So I switched these. I said this the other way around. That's funny. I went through this, and then in my notes, I had them switched, and I said, oh my gosh, I did this wrong. And I've just gone through it, and it turns out I was right in my notes.

Now, one of the things, if you look at what's going on-- let's just go back to this picture here. What's going on here is each one of the passes that we're doing to do a merge is basically taking n over B misses to do a binary merge. We're going through all the data to merge just two things, and traversing all the data. So you can imagine, what would happen if I did, say, a four-way merge? With a four-way merge, I could actually merge four things with only a little bit more than n over B misses. In fact, that's what we're going to analyze in general. So the idea is that we can improve our cache efficiency by doing multi-way merging.

So the idea here is, let's merge R, which is, let's say, less than n subarrays with a tournament. So here are R subarrays, and here's they're each, let's say, is of size n over R. And what we're going to do is merge them with a tournament, so this is a tournament where we say, who's the winner of these two, who's the winner of these two, et cetera. And then whoever wins at the top here, we take them and put them in the output, and then we repeat the tournament.

Now let's just look what happens. It takes order R work to produce the forced output. So we got R things here. To playoff this tournament, there are R nodes here. They each have to do a constant amount of comparing before I end up with a single value to put in the output. So it costs me R to get this thing warmed up. But once I find the winner, and I remove the winner whatever chain he might have come along, how quickly can I repopulate the tournament with the next guy? The next guy only has to play the tournament on the path that was there. All the other matches, we know who won. So the second guy only cost me log R to produce the next guy. And the next guy is log R, and so once we get going, to do an R-way merge only costs me log R work per element.

So the total work in merging is R, to get started, plus n log R. Well, R is less than n, so that's just n log R total to do the merging. That's the work. Now, let's take a look at what happens if I now do merge sort with R-way merges. So what I do is if I have only one element, then it's going to cost me order one time to merge it, because there's nothing to do, just put it in the output. Otherwise, I've got R problems of size n over R that I'm going to merge, and my merge takes n log R time to do the merge. So if I look at the recursion tree, I have n log R here starting here, then I branch R ways, and then I have n over R log R to do the R-way branching at the next level, n over R squared log R at the next level, et cetera.

AUDIENCE: You said that the cost of processing is R.

CHARLES LEISERSON: --is n log R.

AUDIENCE: But is--

CHARLES LEISERSON: Upfront, there's an order R cost, but the order R cost is dominated by the n log R, so we don't have to count that separately. We just have to worry about this guy. So as I go through here, I basically end up having a tree which is only log base R of n tall, because I'm dividing things into R pieces each time, rather than into two pieces. So I only go log base R steps till I get to the base case. But I'm doing an R-way merge, so the number of leaves is still n. But now when I add across here, I get n times log R, and I go across here, I get n times log R, because I got R copies of the same thing. Now I've got R squared times n over R squared log R, and so forth. And so at every level, I have n log R. So I have n log R times the number of levels here, which is log base R of n, plus the order n work at the bottom, which we can ignore because it's going to be dominated.

And so what you notice here is, what's log base R of n? That's just log n over log R. So the log Rs cancel, and I get n log n plus n, which is just n log n. So after all that work, we still do the same amount of work, whether I do binary merging or R-way merging, the work is the same. But there's kind of a big difference when it comes to caching. So it's the same work as binary merge sort.

So let's take a look at the caching. So let's assume that my tournament fits in the cache. So I want to make sure that R is less than m over B for some constant R. So when I do constant way, when I consider the R-way merging of contiguous arrays of total size n, the entire tournament plus 1 block from each array can fit in cache. So the tournament is never going to be responsible for generating cache misses, because I'm going to leave the tournament in memory. So if I'm the optimal algorithm, it's going to say, let's just leave the tournament in memory and bring in all the other things as we do the operation. Question?

AUDIENCE: [INAUDIBLE]

CHARLES LEISERSON: Those circles that I had, the tree.

AUDIENCE: Is that a cumulative list of the elements that you've merged in already?

CHARLES LEISERSON: I'm sorry. Is the--

AUDIENCE: Is it a cumulative list of the length arrays that you've merged already?

CHARLES LEISERSON: No, no, no. You haven't merged them. Let's just go back and make sure we understand the algorithm. The algorithm says that what we do is we compare the head of this pair and we produce a single value here, for which whatever is-- because these are already sorted to do the merge. These are already sorted. So I just have the minimum of these two here, and the minimum of these two here, and the minimum of all four of them here. So it's repeated.

When we get to the top, we have now the minimum of all of these guys, and that's the guy that's the minimum overall, we put him in the output array. And now we walk back down the path that he came from, and what we do is we say, oh, this guy had a-- let's walk down the path, let's say we get to this guy. Let's advance the pointer in here and bring out another element. And now we play off the tournament here, play off the guy here, and he advances, and he advances, whatever. And now some other path may be the minimum.

But it only took me log n work. I'm only keeping copies of the element, if you will, or the results of the comparisons along this in the tree here. And that tree, we're saying, fits in the cache, plus 1 block, the first block. Whatever cache block fits in each of these arrays, no matter how much we've gone down, one of those is fitting in memory. So then what happens here is the entire tournament plus one block from each memory can fit in cache, and so therefore the number of cache misses that I'm going to have when I do the merge it's just essentially the time to take faults on that one cache block whenever I exceed it in each array, plus the one for the output that's similar.

And so the total number of cache misses that I'm going to have is going to be n over B, because I'm just striding straight through memory, and that tournament, I don't have to worry about, because it's sitting in cache. And there's enough sitting in cache that all the other stuff, I can just keep one block from each of them in memory and still expect to get it. In fact, you need the tall cache assumption to assume that they all fit in memory. So therefore, the R-way merge sort is then, if it's sufficiently small, once again, we have the case that it fits in memory so I only have the cold misses to get there, n over B, if n is less than cm. And otherwise, it's R copies of the number of cache misses for n over R, plus n over B. Because this is what it took us here to do the merge. We get only n over B faults when we merge, as long as the tournament fits in cache. If the tournament doesn't fit in cache, it's a more complicated analysis.

AUDIENCE: --n over B, that's cold misses. You're getting the stuff--

CHARLES LEISERSON: Yeah, basically, it's the cold misses on the data, yes, basically. Good. So now, let's do the recursion tree for this. So we basically have n over B that we're going to pay at every level, dividing by R, et cetera, down to the point where things fit in cache. And by the time it fits in cache, it's going to be m over B, because n will be approximately m, just as we had before when we were doing the binary case. As soon as the subarray completely fits in memory, I don't have to when I'm doing the sorting.

So this is now analyzing not the merging, this is analyzing the sorting now. This is the sorting, not the merging. So we get down to m over B, and I've gone now log base R, not log base 2 as we did before, but log base R n over cm. The number of leaves is n over cm, and so when I multiply this out, I get the same n over B here, and I've got n over B at every level here. So where's the win? The win is that I have only log base R of n over cm, rather than log base 2 levels in the tree, because the amount that every level cost me was the same, asymptotically. So when I add it up, I get n over B log base R of n over m, instead of n over B log base 2 of n over m.

So how do we tune R? Well if we just look at this formula here, if I want to tune R, what should I do to R to make this be as small as possible? Make it as big as possible. But I had to assume that R was less than some constant times m so that it fits in cache. So that's, in fact, what I do, is I say R is m over B. So it fits in cache, we have at least one block for each thing that we're merging. And then when we do the analysis now, I take log base m over B here, and that's compared to the binary one, which was log base 2, which is a factor of-- because this is just log 2 over log base-- of log of m over B savings in cache misses.

Now, is that a significant number? Let's take a look. So if your l one cache is 32 kilobytes, and we have cache lines of 64 bytes, that is basically the difference in the exponents, 9x savings. For l two cache, we get about a 12x savings. For l three, we get about a 17x savings. Now of course, there are some other constants going on in here, so you can't be absolutely sure that it's exactly these numbers, but it's going to be proportional to these numbers. So that's pretty good savings to do multi-way merging. So generally when you merge, don't merge pairs. Not a very good way of doing it if you want to take good advantage of cache. May give you some ideas for how to improve some sorts that you might have looked at.

Now it turns out that there's a cache oblivious sorting algorithm, where you don't actually have-- that was a cache aware algorithm that knew the size of the cache, and we tune R to get there. There is an algorithm called funnelsort, which is based on recursively sorting and n to the 1/3 groups of n to the 2/3 items. And then you merge the sorted groups with a merging process called an n to the 1/3 funnel.

So this is more for fun, although the sorting algorithm, in my experience, from what others have told me about implementing it and so forth, is probably about 30% slower than the best hand-tuned algorithm. Whereas with matrix multiplication, the cache oblivious algorithms are as good as any cache aware algorithm as a practical matter, here, they're off by about 20% or 30%. So interesting research topic is build one of these things and make it really efficient so that it can compete with real sorts.

So the k funnel merges k cubed items in k sorted lists, incurring this many cache funnels. Here, I did put in the one, for people who are concerned about the ones. And so then, you get this recurrence for the cache misses, because you solve n to the 1/3 problems of size n to the 2/3 recursively, plus this amount for merging. And that ends up giving you this bound, which turns out to be asymptotically optimal. And the way it works is there's basically, a k funnel is constructed recursively.

And the idea is that what we have is, we have recursive k funnels, so this is going to be a merging process, that is going to produce k cubed items by having k to the 3/2 buffers that each are taking the square root of k guys and producing k guys out. So each of these guys is going to produce k, and the square root of k of them for a total of k to the 3/2, but each of these is going to be length square root of k, so we end up with k cubed. and And they basically feed each other, and then they get merged with their own k thing, and each of these then recursively is constructed the same way.

And the basic idea is that you keep filling the buffers, I think I say this here, so that all these buffers end up being in contiguous storage. And the idea is, rather than going and just getting one element out as you do in a typical tournament, as long as you're going to go merge, let's merge a lot of stuff and put it into our buffer so we don't have to go back here again. So you sort of batch your merging in local regions, and that ends up using the cache efficiently in the local regions. Enough of sorting. Let's go on to physics.

So many of you are probably studying in your linear algebra class or elsewhere, the heat equation. So, people familiar with heat diffusion? So it's a common one to do, and these were-- I have a former student, Matteo Frigo, who is a brilliant coder on anything cache oblivious. He's got the best code out there. So the 2D heat equation, what we do is let's let u(t, x, y) be the temperature at time t at point (x, y). And now you can go through the physics and come up with an equation that looks like this, which says that basically the partial of u with respect to t is proportional to the sum of the second partials with respect to x and with respect to y.

So basically what that says is, the hotter the difference between two things, the quicker things are going to adjust, the quicker the temperature moves between them. And alpha is the thermal diffusivity, which has-- different materials have different thermal diffusivities. Say that three times fast. So if we do a simulation, we can end up with heats, say, put like this, and after it looks like this. See if we can get this running here. So now, let me see. So I can move my cursor around and make things. You can just sort of see that it simulates. You can see the simulation is actually pretty slow. Now, on my slide, I have a thing here that says-- let's see if this breaks when we do it again. There we go. So we're getting around 100 frames per minute in doing this simulation.

And so how does this simulation work? So let's take a look at that. It's kind of a neat problem. So this is what happened when I did 6.172 for a little while. It basically gave me that after a while, because it just sort of averages things, smears it out. So what's going on? Let's look at it in one dimension, because it's easier to understand than if we take on two dimensions. So assuming that we have, say, a bar which has no differential in this direction, only in this direction. So then we get to drop the partials with respect to y.

So if I take a look at that, what I can do is what's called a finite difference approximation, which you probably-- who's studied finite differences? So a few people. It's OK if you haven't. That's OK if you haven't, I'll teach it to you now. And then you're free to forget it, because that's not the part that I want you to understand, but it is interesting. So what I can do is look at the partial, for example, with respect to t, and just do an approximation the says, well, let me perturb t a little bit-- that's what it means. So I add delta t minus u of t divided by t plus delta t minus t, which gives me delta t. And I can use that as an approximation for this partial.

Then on the right hand side-- well first of all, let me get the first derivative with respect to x. And basically here what I'll do is I'll do an approximation where I take x plus delta x over 2 minus x minus delta x over 2, and once again, the differences in the terms there ends up being delta x. And now I use that to take the next one. So basically, to take this one, I basically take the partial with respect to delta x over 2, minus the partial with x minus delta x over 2, and take the partial of that, do the approximation. And what happens is, if you look at it, when I take a partial here I'm adding delta x over 2 twice, so I end up getting just a delta x here, and then the two things on either side combined give me my original one, 2 times u(t, x), and then another one here, and now the whole thing over delta x squared.

And so what I can do is to reduce this heat equation, which is continuous, to something that we can handle in a computer, which is discrete, by saying OK, let's just do this approximation that says that this term must be equal to that term. And if you've studied the linear algebra that said that there are all kinds of conditions on convergence, and stability, and stuff like that, that are actually quite interesting from a numerical point of view, but we're not going to get into it.

But basically, I've just taken that equation here and said, OK, that's my approximation for this one. And now what do I have here? I've got u of t plus delta t, and u things u of t, and then over here, they're all with t, but now the deltas are in-- whoops, that should have been a delta x there. I don't know how that got there. That should be a delta x there. They're all in spatial over here.

So what I can do is take this, and do an iterative process to compute this. And so the idea is, let me take this and throw this term onto the right hand side, and look at u of t plus delta t as if it's t plus 1. Let me make my delta be one, essentially. Throw the delta t over here times the alpha over delta x1, and then I get basically u of tx is based on u of t of x plus 1 and of x and x minus 1. As I say, there's a typo here. That should be a delta t.

So what that says is that if I look at my one-dimensional process proceeding through time, what I'm doing is updating every point here based on the three points below it, diagonally to the right, and diagonally to the left. So this guy can be updated because of those. These we're not going update, because they're the boundary. So these can be fixed. In a periodic stencil, they may even wrap around like a torus. So basically, I can go through and update all these with whatever that hairy equation is. And this is basically what the code is that I showed you is doing. It just keeps updating everyone based on three until I've gone through a bunch of time, and that's how the system evolves. So any questions about how I got to here?

So we're going to now look at this purely computer sciencey now. We don't have to understand any of those equations. We just have to understand the structure. The structure is that we're updating t plus 1 based on stuff on three points with some function that some physicist oracle gave us out of the blue. And so here is a pretty simple algorithm to do it. I basically have what's called the kernel, which does this updating, basically updating each one based on things. And what I'm going to do for computer science is I don't need to keep all the intermediate values.

And so what I'm going to do is do what's called an even-odd trick. Basically if I have one row, I compute the next row into another array, and then I'll reuse that first array-- it's all been used up-- and go back to the first one. So basically here, I'm just going to update t plus 1 mod 2, and just allocate two arrays of size n, and just do modding all the way up. Is that clear? And other than that, it's basically doing the same thing, and I'm doing a little bit of fancy arithmetic here by passing-- see stuff, where I'm passing the pointer to where I am in the array, so I only have to update it locally within the array. So I don't have to double indexing once I'm in the array, because I'm already indexed into the part of the array that I'm going to use, and then I am doing flipping. So this is just a little bit of cleverness. You might want to study this later.

So what's happening then is I have this double nested loop where I have a time loop on the outside, and a space loop on the inside, and I'm basically going through and using a stencil of this shape, this is called a three point stencil, because you're basically taking three points to update one point. And now if I imagine that this dimension is bigger, n here is bigger than my cache size, what's going to happen? I'm going to take a cache fault here, these are all cold misses, et cetera. But when I get back to the beginning here, if I use LRU, nothing is going to be in memory that I happened to update over here. So I have to go and I take a cache fault on every cache line.

And so if I'm going t steps into the future from where I started, I basically have n times t updates, and I save a factor of B, because I get the spatial locality because the u of t minus 1, u of t, and u of t plus 1, are all generally on the same-- are nearby, and all within one cache line. Question?

AUDIENCE: The x's, what are the x's for?

CHARLES LEISERSON: Sorry. I should have put the legend on here. The x's are a miss. So I do a miss when I update these, and then these I don't miss on, because it was brought in when I accessed that. And then I do a miss, and I'll do it-- so basically I do it, then I shift over the stencil by one, and then I won't get a miss. So I'm just looking at the misses on the reads, not misses on the writes. I should have made that clear, too. But the point is, the writes don't help you, because it's all out of memory by the time I get up here. To the second row, if this is longer than my cache size, none of that's there if I'm using LRU.

AUDIENCE: You have also a miss, like you need to get two [INAUDIBLE].

CHARLES LEISERSON: Yeah, but what I'm saying is I'm only looking at the read misses. Yes, there are write misses as well, but basically, I'm only doing the read misses. You can look at the write misses as well. It makes the picture messier. So we've basically have nt over b. However this, let me tell you, is the way that everybody codes it. and And if you have a machine where you have any bandwidth issues to memory, especially for these large problems, this is not a very good way to do it, as it turns out.

So it turns out that what you want to do is, as we've seen, divide and conquer is a really good way to do it. But in this case, when we're doing divide and conquer, we're actually not going to use rectangles, we're going to use trapezoids. And the reason is that a trapezoid has the nice property that-- notice that if I have all these points in memory, then notice that I can compute all the guys that are read on the next level, and then I can compute all the guys that are next on the next level. And so for example, if you imagine that this part here fit within cache, I could actually keep going. I didn't have to stop here, I could keep going right up to a triangle if I wanted to, and compute all the values without having any more cache misses than those needed to bring in, essentially, one row-- two rows, actually, because I'm reusing the rows as I go up.

So what we're going to do is traverse trapezoidal regions of space-time points such that the points are between an upper limit, T1, and a low one, T0, and between an x0 and an x1, where now I have slopes here that are going to be, in general, this is plus 1 minus 1. And in fact, sometimes it will be straight, in which case we'll call it 0. It's really the inverse of the slope, but we'll still call it zero rather than infinity. So it's 1 over the slope. There's a name for that, right? Is that called the run or something? I forget, I don't remember my calculus.

So that's what we're going to do. And we're going to leave the upper and right borders undone and include, so it's going to be a sort of half open trapezoid on the left and bottom, closed on the left and bottom, and open on the top and right. So the width is basically the midpoint here, and the height is the height, because they're always going to have parallel axes here. So here's how are our recursion is going to work. If the height is 1, then we can compute all space-time points in any way we want. I can just go through them if I want, because they're all independent. None depends on anybody else. So that's going to be our base case.

If the width is greater than twice the height, however, then what we're going to do is we're going to cut the trapezoid through the middle of the slope of minus 1. And that will produce two new trapezoids, which we then will recursively compute all the elements of. So I'll start out with a trapezoid. Basically, if it ends up that it's a long and wide one, I'm going to make what's called a space cut, and cut it this way, and then I'm going to recursively do this one and then this one. And notice that I can do that because-- all these guys I can do, but then when I get to the border here, this will already have been done by the time I'm computing these guys.

So the requirement is that I've got to do things according to that map of triples that I showed you before, but I don't have to do them in the same order. I don't have to do the whole bottom row first. In this case, I can compute the whole trapezoid here, and then I can compute this trapezoid here, and then all the values that I'll need will have already been computed over here, that are on the boundary of this trapezoid. The other type of cut I'll do is what happens when a trapezoid gets too tall for me. So if the trapezoid is too tall, then what we'll do is we'll slice it through the middle, but the other way. We call that a time cut. So we won't take it all the way through time, we'll only take it partially through time. Now you can show, and I'm not going to show this in detail, but you can show that if I do this, my trapezoids are always sort of medium sized. I never get long, long skinny ones. If I start with something that's sort of got a good aspect ratio, I maintain a good aspect ratio through the entire code.

So here's the implementation. This is what Matteo Frigo wrote, and I've modified it a little bit. So basically, we pass in it the values that let us identify the trapezoid, t0, t1, x0, and then the slope on the left side, x1 in the slope on the right side, where the dx0 and the dx1s are all either 0, 1, or minus 1. And then what I do is I look at what the height is that my trapezoid is going to operate on. And if the height is 1, well, then I just run through all the elements, and I just compute the kernel-- that program that I showed you before, that kernel-- on all the elements. Nothing really to be done there, just go through and compute them individually with a four loop.

Otherwise, if I've got the situation where the trapezoid is big, then I do this comparison, which I promise you-- you can work out the math if you wish-- which I promise you tells you whether or not it's more than twice the height, as I said before, whether the width is more than twice the height. And if so, I compute the middle point, and then I partition it into two trapezoids, and I recursively call them. And otherwise, I simply cut the time in half, and then I do the bottom half and then the upper half.

So getting all those parameters exactly right takes a little bit of thinking, makes my brain hurt, but Matteo is brilliant at this kind of coding. So let's see how well this does. So I'm not going to do a detailed analysis that I did before, but basically what's going on is at this level, if I'm doing divide and conquering, I'm only doing a constant amount of work managing this stuff. So my caches that I'm taking in the internal part of the tree are all going to be order one.

Now each leaf is going to represent a trapezoid, which is going to be approximately h times w, where h and w are approximately equal, because they're going to be shaped-- This is assuming I start out with a number of iterations to do that is at least as large as the number of points that I have to go on. If I start out with something that's really thin and flat, then it's not going to be the case. But if I start out with something that's deep enough, then I'm going to be able to make progress in an unconventional order into time by moving the time non-uniformly through the space. So each leaf represents a fairly balanced trapezoid.

Each leaf basically is going to-- if you look that the direction of the trapezoid is in time, so this represents the spatial dimension, and if I have something of size w, I can access it with only w over B misses. And when that fits in cache, where w is some constant less than m, so w is order m. So each of these things that's a leaf is only going to occur w over B misses.

Now, the total space number of points I have to go after is n times t. N is going to be the full dimension this way, t is the height that way. And so since each leaf has hw points, I have nt over hw leaves. And the number of internal nodes is just the leaves minus 1, so they can't contribute substantially, because there's only order one misses I'm taking here, whereas I've got something on the order of w over B misses for this. So therefore, now I can do my math. The number of cache misses I'm going to take is-- well, how many leaves do I have? nt over hw. And what does each one cost us? W over B.

And so now, when I multiply that out, well, hw is about m squared, and w over B is about m over B. And so I get nt over MB as being the total number of savings. so whereas the original algorithm only got nt over B, we've got this factor of a memory cache size in there showing us that we have far fewer cache misses. So the cache misses end up not being an issue for this. Any questions about that?

So I want to show you a simulation of this three point stencil and comparing the two things. So this is going to be the looping version, where the red dots are where the cache misses are, and this is going to be the trapezoidal one. And basically, I have an n of 95 and a t of 87, and what I'm going to do is assume a fully associative LRU cache that fits four points on a cache line, where the cache size is 32, two to the fifth as opposed to two to the 15th, it's really little. If I get a cache hit, I'm going to call it one cycle. If I get a cache miss, I'm going to call it 10 cycles. We're going to race them. So on the left is the current world champion, the looping algorithm. And on the right is the cache oblivious trapezoid algorithm. So let's go.

So you can see that it's basically, it's made a space cut there, but it's made a time cut across the top there. It said, this is too tall, so let me cut it this way. And that guy's, meanwhile, taking all those-- you can see how many cache misses he's taking. Let's speed him up. That's one way you can do it, is make it think.

So let's see what happens if I have a cache of size eight. So here we go. I think I'm just doing the same thing. I know I can show you the cuts. Can I show you the cuts? I know. I think it's because I'm not-- OK, let's try it. There we go. Now I'm showing the cuts as they go on. Let's do that again. We'll go fast and do it again. So those are the cuts that it's making to begin with as it's doing the divide and conquer. So I think this is the same size cache. So now I think I'm doing a bigger cache. I think I did a bigger cache, but I'm not sure I gave the other guy a bigger cache. Yeah, because it doesn't matter for the guy on the left, right? As long as the cache line is the same length and as long as it's not big enough, he's going to do the same thing. He didn't get to take advantage of the fact that the cache was bigger, because it was still smaller than the array that he's striping out there. Anyway, we can play with these all day.

So if you make the cache lines bigger, then of course it'll go faster, because he'll have fewer misses. He'll get to bring it in. So let's see here. So let's now do it for real. So this is a two-dimensional problem. You can use the same thing to do what end up being three-dimensional trapezoids. In fact, you can generalize this trapezoid method to multiple dimensions. So this is the looping one. So let's start that one out. So it's going about 104 frames a minute. I think by resizing it, the calibration is off.

But in any case, let's switch to the cash oblivious version. Anybody notice something? Slower. Why is that? I gave code exactly as I had up there. No, it's not because it's two dimensions.

AUDIENCE: [INAUDIBLE].

CHARLES LEISERSON: I'm sorry?

[INTERPOSING VOICES]

CHARLES LEISERSON: Yeah, so now it's the trapezoiding at only 86 frames. What do you suppose is going on there?

AUDIENCE: You have [INAUDIBLE].

CHARLES LEISERSON: Yeah. So this is a case where if you look at the code I wrote, I went down to a t, a delta t, of one in my recursion. I recursed all the way down. Let's see what happens if instead of going all the way down, playing the trapezoid game on little tiny trapezoids, suppose I go down only to, say, when t is 10, and then do essentially the row major ones. So I'm basically coarsening the leaves of the thing. So to do that, I do this. So now we go-- ah.

So I have to coarsen in order to overcome the procedure call overhead. It has nothing to do with the cache. It has to the fact that the way that you implement recursion, recursion and function calls have a cost to them. And if what you're going to do is do a little tiny update of those few floating point operations-- let's go back to the looping just to see. The looping is going about 107, 108, and trapezoiding at 136. So unfortunately, you need a voodoo variable, but it's a voodoo variable not to overcome the cache, but rather to deal with what's the overhead in using the processor when you do function calls.

So let's see. How coarse can we make it? Let's try five, a coarsening of five? That's still pretty good. That's still pretty good. How about four? Still doing 131 frames a minute. How about three? Oh, we lost something there. How about two? So at a coarsening of two, I go 138, whereas the looping goes at about the same. I can't do 20. I didn't program that in. I just programmed up to 10. So if I go down to one, however, then you see it's not that efficient. But if I pick any number that's even slightly larger, that gives me just enough that the function call overhead ends up not being a substantial cost of the things.

So let me just wrap up now. So I just have a couple more things. So I'm not going to really talk about these, but there are lots of cash oblivious algorithms that have been discovered in the last 10 or 15 years for doing things like matrix transposition, which is similar to rotating a matrix. You can do that in a cache oblivious fashion. Strassen's algorithm, which does matrix multiplication using fewer than n cubed operations. The FFT can be computed in a cache oblivious fashion. And LUD composition is a popular thing to solve systems.

In addition, there are cache oblivious data structures, and here are just a few of them. There's cache oblivious B-Trees and priority queues, and doing things called ordered-file maintenance. There's a whole raft. There's probably now several hundred papers written on cache oblivious algorithms, so something you should be aware of and understand how it is that you go about designing an algorithm of this nature. Not all of them are straightforward. For example, the FFT one does divide and conquer but not by dividing it into two. It divides it into square root of n pieces of size square root of n each in order to get a good cache efficient algorithm that doesn't have any tuning parameters. But almost all of them, since they're recursive, do have this annoying thing that you have to still coarsen the base case in order to get really good performance if you're not doing a lot of work in the leaves of the recursion. So any questions?