Cache Oblivious Algorithms Piyush Kumar? Departmentof ComputerScience State Universityof New Yorkat StonyBrook StonyBrook, NY11790, USA [email protected] http://www.compgeom.com/co-chap/y 1 Introduction The cache oblivious model is a simple and elegant model to design algorithms thatperformwellinhierarchicalmemorymodelsubiquitousoncurrentsystems. This model was (cid:12)rst formulated in [20] and has since been a topic of intense research.Analyzing and designing algorithmsand data structures in this model involvesnotonlyanasymptoticanalysisofthenumberofstepsexecutedinterms of the input size, but also the movement of data optimally among the di(cid:11)erent levels of the memory hierarchy. This chapter is aimed as an introduction to the \ideal-cache" model of [20] and techniques used to design cache oblivious algorithms. The chapter also presents some experimental insights and results. A dream machine would be fast and would never run out of memory. Since an in(cid:12)nite sized memory wasnever built, one hasto settle for varioustrade-o(cid:11)s in speed, size and cost. In both the past and the present, hardware suppliers seem to have agreed on the fact that these parameters are well optimized by buildingwhatiscalledamemoryhierarchy(SeeFigure1,Chapter??).Memory hierarchiesoptimize the three factorsmentioned aboveby being cheapto build, trying to be as fast as the fastest memory present in the hierarchy and being almost ascheap asthe slowestlevel of memory. The hierarchyinherently makes use of the assumption that the access pattern of the memory has locality in it and can be exploited to speed up the accesses. The locality in memory access is often categorized into two di(cid:11)erent types, code reusing recently accessed locations (temporal) and code referencing data items that are close to recently accessed data items (spatial) [22]. Caches use both temporal and spatial locality to improve speed. Surprisingly many things canbecategorizedascaches,forexample,registers,L1,L2,TLB,Memory,Disk, Tape etc (Chapter ??). The whole memory hierarchycan be viewed as levels of caches,eachtransferringdatatoitsadjacentlevelsinatomicunitscalledblocks. When data that is needed by a process is in the cache, a cache hit occurs. A cache miss occurs when data can not be supplied. Cache misses can be very y Anupdatedversion of the chapter can befound atthe webpage ? Part of this work was done while the author was visiting MPI-Saarbru(cid:127)cken. The authorispartiallysupportedbyNSF(CCR-9732220,CCR-0098172) andbyagrant from SandiaNational Labs. 2 Piyush costly in terms of speed. Cache misses can be reduced by designing algorithms that use locality of memory access. TheRandom Access Model (RAM) inwhich wedoanalysisof algorithmsto- daydoesnottakeintoaccountdi(cid:11)erencesinspeedsofrandomaccessofmemory dependinguponthelocalityofaccess[16].Althoughthereexistmodelswhichcan dealwithmulti-levelmemoryhierarchies,theyarequitecomplicatedtouse[1,3, 6,5,2,23,28,29].Itseemsthereisatrade-o(cid:11)betweentheaccuracyofthe model and the ease of use. Most algorithmic work has been done in the RAM model which models a ’(cid:13)at’ memory with uniform access times. The external memory model (with which the reader is assumed to be familiar, See Chapter ??) is a two level memory model, in the context of memory and disk. The ’ideal’ cache oblivious model is a step towards simplifying the analysis of algorithms in light of the fact that localaccessesin memory arecheaperthan non-localonesin the whole memory hierarchy(and not just two-levelsof memory).It helps takeinto account the whole memory hierarchyand the speed di(cid:11)erences hidden therein. In the external memory model, algorithm design is focussed on a particular levelofmemory(usuallythedisk)whichisthebottleneckfortherunningtimeof the algorithminpractice.Recall(From Chapter??),thatintheexternalmem- orymodel,processingworksalmostasinthe RAMmodel,exceptthatthereare only M words of internal memory that can be accessed quickly. The remaining memorycanonlybeaccessedusingI/OsthatmoveB contiguouswordsbetween external and internal memory. The I/O complexity of an algorithm amounts to counting the number of I/Os needed. The cache oblivious model was proposed in [20] and since then has been used in more than 30 papers already. It is becoming popular among researchers in external memory algorithms, parallel algorithms, data structures and other related(cid:12)elds.Thismodelwasbornoutofthenecessitytocapturethehierarchical natureofmemoryorganization.(ForinstancetheIntelItaniumhas7levelsinits memory hierarchy, See Chapter ??). Although there have been other attempts to capture this hierarchical information the cache oblivious model seems to be oneofthemostsimpleandelegantones.Thecacheobliviousmodelisatwolevel model (like the [4] model that has been used in the other chapters so far) but with the assumption that the parameters M,B are unknown to the algorithm (See Figure 1). It can work e(cid:14)ciently on most machines with multi-level cache hierarchies. Note that in the present volume, all external memory algorithms that have been dealt with yet, need the programmer/userto specify M,B. This chapter is intended as an introduction to the design and analysis of cache oblivious algorithms, both in theory and practice. Chapter Outline:We introducethe cacheobliviousmodel insection 2. In section3weelaboratesomecommonlyuseddesigntoolsthatareusedtodesign cache oblivious algorithms. In section 4 we choose matrix transposition as an example to learn the practical issues in cache oblivious algorithm design. We studythecacheobliviousanalysisofStrassen’salgorithminsection 5.Section6 discusses a method to speed up searching in balanced binary search trees both in theory and practice. In section 7, a theoretically optimal, randomized cache Cache Oblivious Algorithms 3 oblivious sorting algorithm along with the running times of an implementation is presented. In section 8 we enumerate some practicalities not caught by the model.Section9presentssomeofthebestknownboundsofothercacheoblivious algorithms. Finally we conclude the chapter by presenting some related open problems in section 10. 2 The Model The ideal cache oblivious memory model is a two level memory model. We will assume that the faster level has M size and the slower level always transfers B words of data together to the faster level. These two levels could represent the memory and the disk, memory and the cache, or any two consecutive levels of the memoryhierarchy(SeeFigure1).Inthischapter,M andB canbe assumed to be the sizes of any two consecutive levels of the memory hierarchy subject to some assumptions about them (For instance the inclusion propertywhich we will see soon). We will assume that the processor can access the faster level of memorywhichhassizeM.Iftheprocessorreferencessomethingfromthesecond levelofmemory,anI/OfaultoccursandB wordsarefetchedintothefasterlevel ofthe memory.Wewillrefertoablockasthe minimumunitthatcanbe present orabsentfromalevelinthetwolevelmemoryhierarchy.WewilluseB todenote the size of a block as in the external memory model. If the faster level of the memory is full (i.e. M is full), a block gets evicted to make space. The ideal cache oblivious memory model enables us to reason about a two level memory model like the external memory model but prove results about a multilevel memory model. Compared with the external memory model it seems surprisingthatwithoutanymemoryspeci(cid:12)cparametrization,orinotherwords, without specifying the parameters M;B, an algorithm can be e(cid:14)cient for the wholememoryhierarchy,neverthelessitspossible.Themodelisbuiltuponsome basic assumptions which we enumerate next. Assumptions:The following four assumptions are key to the model. Optimal replacement The replacement policy refers to the policy chosen to replaceablockwhenacachemissoccursandthecacheisfull.Inmosthard- ware, this is implemented as FIFO, LRU or Random. The model assumes thatthecachelinechosenforreplacementistheonethatisaccessedfurthest in the future. The strategy is called optimal o(cid:11)-line replacement strategy. 2 levels of memory Therearecertainassumptionsinthemodelregardingthe twolevelsofmemorychosen.Theyshouldfollowtheinclusionpropertywhich says that data cannot be present at level i unless it is present at level i+1. In most systems, the inclusion property holds. Another assumption is that the sizeofleveliof the memoryhierarchyisstrictlysmallerthan leveli+1. Full associativity When ablockof dataisfetchedfromthe slowerlevelof the memory, it can reside in any part of the faster level. Automatic replacement When a block is to be brought in the faster level of thememory,itisautomaticallydonebytheOS/hardwareandthealgorithm 4 Piyush Faster CPU Smaller Registers On Chip Cache On Board Cache Main Memory Disk Tapes Slower Bigger Two consecutive levels of the hierarchy Cache Level 1 Misses Level 2 CPU Memory Disk Size =M Size =B Cache Lines Block Transfer Fig.1.The ideal cache oblivious model designer does not have to care about it while designing the algorithm. Note that we could access single blocks for reading and writing in the external memory model, which is not allowed in the cache oblivious model. Wewillnowexamineeachofthe assumptionsindividually.Firstweconsider the optimal replacementpolicy. The most commonlyused replacementpolicy is LRU (least recently used). In [20] the following lemma, whose proof is omitted here, was proved using a result of [30]. Lemma 1 ([20]). An algorithm that causes Q (n;M;B) cache misses on a (cid:3) problem of size n using a (M;B)-ideal cache incurs Q(n;M;B) 2Q (n;M;B) (cid:20) (cid:3) 2 cache misses on a (M;B) cache that uses LRU, FIFO replacement. This is only true for algorithms which follow a regularity condition. An algorithm whose cache complexity satis(cid:12)es the condition Q(n;M;B) (cid:20) O(Q(n;2M;B)) is called regular (All algorithms presented in this chapter are regular).Intuitively,algorithmsthatslowdownbyaconstantfactorwhenmem- ory (M) is reduced to half, are called regular. It immediately follows from the abovelemmathatifanalgorithmwhosenumberofcachemissessatis(cid:12)esthereg- ularity condition does Q(n;M;B) cache misses with optimal replacement then thisalgorithmwouldmake(cid:2)(Q(n;M;B)) cachemissesonacachewith LRUor FIFO replacement. The automatic replacement and full associativity assumption can be imple- mentedinsoftwarebyusingLRUimplementationbasedonhashing.Itwasshown in [20] that a fully associative LRU replacement policy can be implemented in O(1) expected time using O(M) recordsof size O(B) in ordinarymemory. Note B Cache Oblivious Algorithms 5 that, the above description about the cache oblivious model also proves that anyoptimalcacheobliviousalgorithmcanalsobeoptimallyimplementedinthe external memory model. Wenowturnourattentiontomulti-levelidealcaches.Weassumethatallthe levels of this cache hierarchy follow the inclusion property and are managed by anoptimalreplacementstrategy.Thusoneachlevel,anoptimalcacheoblivious algorithm will incur an asymptotically optimal number of cache misses. From Lemma1,thisbecomestrueforcachehierarchiesmaintainedbyLRUandFIFO replacement strategies. Apart from not knowing the values of M;B explicitly, some cache oblivious algorithms(forexampleoptimalsortingalgorithms)requireatallcacheassump- tion. The tall cache assumption states that M = (cid:10)(B2) which is usually true in practice. Its notable that regular optimal cache oblivious algorithms are also optimalinSUMH[5]andHMM[1]models.Recently,compilersupportforcache oblivious type algorithms have also been looked into [34,35]. 3 Algorithm design tools In cache oblivious algorithm design some algorithm design techniques are used ubiquitously. One of them is a scan of an array which is laid out in contiguous memory. Irrespective of B, a scan takes at most 1+ N I/Os. The argument dBe is trivial and very similar to the external memory scan algorithm Chapter ??. The di(cid:11)erence is that in the cache oblivious setting the bu(cid:11)er of size B is not explicitly maintained in memory. In the assumptions of the model, B is the size of the data that is always fetched from level 2 memory to level 1 memory. The scan does not touch the level 2 memory until its ready to evict the last loaded bu(cid:11)er of size B already in level 1. Hence, the total number of times the scan algorithm will force the CPU to bring bu(cid:11)ers from the level 2 memory to level 1 memory is upper bounded by 1+ N . dBe Exercise 1 Convince yourself that indeed a scan of an array will at most load O(N) bu(cid:11)ers from the level 2. Re(cid:12)ne your analysis to 1+ N bu(cid:11)er loads. B dBe 3.1 Main Tool: Divide and conquer Chances are that the reader is already an expert in divide and conquer algo- rithms. This paradigm of algorithm design is used heavily in both parallel and externalmemoryalgorithms.Itshouldnotcomeasasurprisetothe readerthat cache oblivious algorithms make heavy use of this paradigm and lot of seem- ingly simple algorithms that were based on this paradigm are already cache oblivious! For instance, Strassen’s matrix multiplication, quicksort, mergesort, closestpair[16],convexhulls[7],medianselection[16]areallalgorithmsthatare cache oblivious, though not all of them are optimal in this model. This means that they might be cache oblivious but can be modi(cid:12)ed to make fewer cache misses thantheydoin the currentform.Inthe sectiononmatrixmultiplication 6 Piyush wewillseethatStrassen’smatrixmultiplicationalgorithmisalreadyoptimalin the cache oblivious sense. Whydoesdivideandconquerhelpingeneralforcacheobliviousalgorithms? Divide and conquer algorithms split the instance of the problem to be solved into several sub problems such that each of the sub problems can be solved independentally.Sincethealgorithmrecursesonthesubproblems,atsomepoint of time, the sub problems (cid:12)t inside M and subsequent recursion, (cid:12)ts the sub problems into B. For instance lets analyze the average number of cache misses in a randomized quicksort algorithm. This algorithm is quite cache friendly if not cache optimal and is described in Section 8.3 of [16]. Lemma 2. RandomizedversionofquicksortincursanexpectednumberofO(N log (N)) B 2 B cache misses. Proof. Choosing a random pivot makes at most one cache miss. Splitting the inputsetintotwooutputsets,suchthattheelementsofonearealllessthanthe pivot and the other greater than or equal to the pivot makesat most O(1+ N) B cache misses by Exercise 1. Assoonasthesizeoftherecursion(cid:12)tsintoB,therearenomorecachemisses for that subproblem (Q(B) = O(1)). Hence the average cache complexity of a randomizedquicksortalgorithmcanbegivenbythe followingrecurrence(which is very similar to the average case analysis presented in [16]). 1 N Q(N)= (Q(i)+Q(N i))+(1+ ) (1) N 2 (cid:0) dBe 3 i=1X::NB(cid:0)1 4 5 which solves to O(N log N). B 2 B Exercise 2 Solve the recurrence 1. Hint 1 First simplify the recurrence to 2 N Q(N)= Q(i)+(cid:2)(1+ ) (2) N 2 B 3 i=1X::NB(cid:0)1 4 5 and then you can take help from [21]. Exercise 3 Prove that mergesort also does O(N log N) cache misses. B 2 B As we said earlier, the number of cache misses randomized quicksort makes is not optimal. The sortinglowerbound in the cacheoblivious model is also the same as the external memory model (See Chapter ??). We will see later in the chapter a sorting algorithm that does sorting in O(N log N) cache misses. B M B B We now will analyze another cache oblivious algorithm that is also optimal, the median (cid:12)ndingasgivenin[16](Section10.3).Therecurrencefor(cid:12)nding the median of N entities is given by: Cache Oblivious Algorithms 7 N 7N N Q(N)=Q( )+Q( +6)+O(1+ ) 5 10 B Solving which we get, Q(N) = O(1+ N). Note that again Q(B) = O(1) so B therecursionstopsaddingcosttocachemisseswhentheinstanceoftheproblem becomes smaller than B. 4 Matrix transposition Matrix transposition is a fundamental operation in linear algebra and in fast fourier transforms and has applications in numerical snalysis, image processing and graphics. The simplest hackfor transposing a N N square matrix in C++ (cid:2) could be: for (i = 0; i < N; i++) for (j = i+1; j < N; j++) swap(A[i][j], A[j][i]) In C++ matrices are stored in \row-major" storage, i.e. the rightmost di- mension varies the fastest ( Chapter ??). In the above case, the number of cache misses the code could do is O(N2). The optimal cache oblivious matrix transposition makes O(1+ N2) cache misses. Before we go into the divide and B conquerbasedalgorithmformatrixtranspositionthatiscacheoblivious,letssee an experimental results. (See Figure 2). The (cid:12)gure implements a blocked cache obliviousimplementation,westopthe recursionwhenthe problem sizebecomes less than a certain block size and then use the simple for loop implementation insidethe block.Notethatusingdi(cid:11)erentblocksizeshaslittlee(cid:11)ectontherun- ning time. Thisexperiment wasdone onWindows NT running on1Ghz/512Mb RAM notebook. The code was compiled with g++ on cygwin. Here is the C/C++ code for cache oblivious matrix transposition. The fol- lowing code takes as input a submatrix given by (x;y) (x+delx;y+dely) in (cid:0) the inputmatrixI andtransposesittotheoutputmatrixO.ElementType1can be any element type , for instance long. 1 In all ourexperiments, ElementType was set to long. 8 Piyush void transpose(int x, int delx, int y, int dely, ElementType I[N][P], ElementType O[P][N]){ // Base Case of recursion // Should be tailored for specific machines if one wants // this code to perform better. if((delx == 1) && (dely == 1)) { O[y][x] = I[x][y]; return; } // Divide the transposition into two sub transposition // problems, depending upon which side of the matrix is // bigger. if(delx >= dely){ int xmid = delx / 2; transpose(x,xmid,y,dely,I,O); transpose(x+xmid,delx-xmid,y,dely,I,O); return; } // Similarly cut from ymid into two subproblems ... } Thecodeworksbydivideandconquer,dividingthebiggersideofthematrix in the middle and recursing. Exercise 4 ProvethatonanN N inputmatrix,theabovecodecausesatmost O(1+ N2) cache misses. (cid:2) B Hint 2 Let the input be a matrix of N P size. There are three cases: (cid:2) Case I max N;P (cid:11)B In this case, f g(cid:20) NP Q(N;P) +O(1) (cid:20) B Case II N (cid:11)B <P In this case, (cid:20) O(1+N) if (cid:11)B P (cid:11)B Q(N;P) 2 (cid:20) (cid:20) (3) (cid:20)(2Q(N;P=2)+O(1) N (cid:11)B <P (cid:20) Case III P (cid:11)B <N Analogous to Case III. (cid:20) Case IV min N;P (cid:11)B f g(cid:21) O(N +P + NP) if (cid:11)B N;P (cid:11)B B 2 (cid:20) (cid:20) Q(N;P) 2Q(N;P=2)+O(1) P N (4) 8 (cid:20) (cid:21) ><2Q(N=2;P)+O(1) N P (cid:21) >: Cache Oblivious Algorithms 9 Fig.2. The graph compares a simple for loop implementation with a blocked cache oblivious implementation. The above recurrence solves to Q(N;P)=O(1+ NP). B Thereisasimplerwaytovisualizetheabovemess.Oncetherecursionmakes the matrix small enough such that max(N;P) (cid:11)B (cid:12)pM (here (cid:12) is a (cid:20) (cid:20) suitableconstant),orsuchthatthesubmatrix(ortheblock)weneedtotranspose (cid:12)ts in memory, the number of I/O faults is equal to the scan of the elements in thesubmatrix.Apackingargumentofthesenotsosmallsubmatrices(blocks)in the large input matrix shows that we do not do too many I/O faults compared to a linear scan of all the elements. Remark:Figure2showsthee(cid:11)ectofusingblockedcacheobliviousalgorithm for matrix transposition. Note that in this case, the simple for loop algorithm is almost always outperformed. This comparison is not really fair. The cache oblivious algorithm gets to use blocking whereas the naive for loop moves one elementatatime.Acarefulimplementationofablockedversionofthesimplefor loopmightbeatthe blockedcacheoblivioustranspositionalgorithmin practice. (See the timings of Algorithm 2 and 5 in [13]). The same remark also applies to matrix multiplication. (Figure 3) 5 Matrix multiplication Matrix multiplication is one of the most studied computational problems: We aregiventwomatricesof m n andn pand wewanttocompute the product (cid:2) (cid:2) 10 Piyush matrixof m psize.Inthissectionwewillusen=m=N althoughthe results (cid:2) can easily be extended to the case when they are not equal. Thus, we are given twoN N matricesx=(x ),y =(y ),andwewishtocomputetheirproduct i;j i;j (cid:2) z, i.e. there are N2 outputs where the (i;j)’th output is N z = x y i;j i;k k;j (cid:1) k=1 X In1969’StrassensurprisedtheworldbyshowinganupperboundofO(Nlog27)[31] using a divide and conquer algorithm. This bound was later improved and the best upper bound today is O(N2:376) [15]. We will use Strassen’s algorithm as given in [16] for the cache oblivious analysis. The algorithm breaks the three matrices x;y;z into four submatrices of size N N, rewriting the equation z =xy as : 2 (cid:2) 2 r s a c e g = t u b d (cid:2) f h (cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19) Nowwehave8subproblems(matrixmultiplications)of size N N andfour 2 (cid:2) 2 additions (which will lead us to a O(N3) algorithm). This can be done more cleverly by only solving 7 subproblems and doing more additions and subtrac- tions but still keeping the number of additions and subtraction a constant (See [16] for more details). The recurrence for the internal work now becomes : N T(N)=7T( )+(cid:2)(N2) 2 which solves to T(N) = O(Nlg7) = O(N2:81). The recurrence for the cache faults is: O(1+N + N2) if N2 (cid:11)M Q(N) B (cid:20) (5) (cid:20)(7Q(N)+O(1+ N2) otherwise: 2 B whichsolvestoQ(N) O(N+N2 + Nlg7 ).Wealsoimplemented ablocked (cid:20) B BpM cacheobliviousmatrixtranspositionroutinethatworksverysimilartothetrans- positionroutineintheprevioussection.Itbreaksthelargestofthethreedimen- sionsofthe matricestobemultiplied(dividesitby2)andrecursestill theblock sizeisreached.NotethatthisalgorithmisnotoptimalanddoesO(1+N2+ N3 ) B BpM cache misses. We also present here experiments of a simpler divide and conquer matrix multiplication algorithm. The worst case of the implementation is O(N +N2+ N3 ). The experimental results are shown in Figure 3. We compare a blocked BpM cacheobliviousmatrixmultiplicationimplementationwithsimpleforloopbased matrix multiplication. This experiment was done on an dual processor Intel Itanium with 2Gb RAM. Only one processor was being used. Note that on this platform, the blocked implementation consistently outperforms the simple loop