A common problem within databases and search engines is to compute the intersection between two sorted array. Typically one array is much smaller than the other one.

The conventional strategy is the “galloping intersection”. In effect, you go through the values in the small arrays and then do a binary search in the large array. A binary search is a simple but effective algorithm to search through a sorted array. Given a target, you compare with with the midpoint value. If your target is smaller than the midpoint value, you search in the first half of the array, otherwise you search in the second half. You can recurse through the array in this manner, cutting the search space in half each time. Thus the search time is logarithmic.

If the small array has M elements and the large array has N elements, then the complexity of a galloping search is O(M log N). In fact, you can be more precise: you never need more than M * log N + M comparisons.

Can you do better? You might.

Let me describe an improved strategy which I call “shotgun intersection”. It has been in production use for quite some time, through the CRoaring library, a C/C++ implementation of Roaring Bitmaps.

The idea is that galloping search implies multiple binary searches in sequence through basically the same array. Doing them consecutively might not be best. A binary search, when the large array is not in cache, is memory-bound: it waits for the memory subsystem to deliver the data. So you are constantly waiting. What if you tried to do something else while you wait. What about starting right away on the next binary search?

That is how a shotgun search works. You take, say, the first four values from the small array. You load the midpoint value from the large array, then you compare all of your four values against this midpoint value. If the target value is larger, you set a corresponding index so that the next search will hit the second half of the array. And so forth. In effect, shotgun search does many binary searches at once.

I make my Java code available, if you want a full implementation.

Does it help? It does. Sometimes it helps a lot. Let us intersect an array made of 32 integers with an array made of 100 million sorted integers. I use a cannonlake processor with Java 8.

1-way | 1.3 microseconds |

4-way | 0.9 microseconds |

**Credit**: Shotgun intersections are based on an idea and an initial implementation by Nathan Kurz. I’d like to thank Travis Downs for inspiring discussions.

Daniel Lemire, "Faster intersections between sorted arrays with shotgun," in *Daniel Lemire's blog*, January 16, 2019.

Why not do the galloping algorithm but alternating from the small array head and tail (1, N, 2, N-1, 3, …) and cache the upper and lower indices of the last binary search on the large array? Not a big saving, approximately M*log 2 fewer ops but possibly comparable.

I’d love to see a benchmark of your idea. Did you try to code it up?

I tried something along the suggested line recently. For each element in the small array I kept an upper and lower limit which I updated depending on the values I encountered while probing. The small array was in memory, the large array on disk. (Since disk probes are very expensive, I figured I could waste some time on updating the limits.) It worked quite nicely in the sense that it saved a lot of disk probes in typical cases.

However, the running time was less than ideal. It took me a while to figure out what was happening. Because the elements of the small array now all had different lower and upper limits, the midpoints would differ ever so slightly all the time, thereby frustrating the disk cache.

This led me to discard the idea of adjusting the limits and instead add a LRU cache to the binary search. Since a large part of the “search tree” is repeated for each lookup, this resulted in a handsome speedup!

It has been fascinating to observe the trend from “let’s design the fastest algorithm” to “let’s design an algorithm that ||-es better”.

If the two arrays are sorted, can’t the same be accomplished using a ‘merge’ operation with complexity O(M) faster? (no comparison with N-array needs to be made twice)

To clear up my last comment, it seems your benchmark code is doing a ‘merge’ in the sense of line 115:

`idx_l = index1;`

and then line 84 says:`int base1 = idx_l;`

. My proposed optimization would be to have 115 use`index4`

to really get ahead in the long array.In any case the O(M*log N + M) remains the complexity.

It depends on whether you expect the arrays to be strictly sorted or merely sorted (possibly with repeated values).

Still can’t see how repeats could make

`idx_l = index4`

fail, since`index4`

points to`target4`

if found and next target value will be at or after`index4`

. If`target4`

not found`index4`

will be just one cell after lower-than-`target4`

value in long list, and next target in short list must surely be higher by non-strict monotonicity (probably should be coding and running this myself).You might be right.

I think you are right. I have updated the code.

It might be a bit faster. New timing?

Deliberately, I am giving rough numbers with one or two significant digits. My timings are not going to be affected by micro-optimizations.

I think the old complexity definition largely doesn’t work nowadays. For posting processing we have a bunch of “same-complexity” methods to do unions and intersections, but they have very different run-times easily orders of magnitude different.

I agree. Even on the JVM, the same or worse complexity can win when some exposed feature of the underlying hardware can be utilised.

Indeed.

Sometimes the problem is simply that the wrong thing is being counted. E.g., most searching algorithms give the complexity in the number of “comparisons”.

However, for searching in large arrays or many small arrays (such that the aggregate is large), what dominates is the memory access time, so you’d be better of counting something like “unique memory locations accessed”, with the idea that once you’ve accessed a location once, it’s basically free to access it again (it’s cached). Of course, even that’s not exactly right: since, for example, caching isn’t forever and also adjacent locations (within a cache line) will also come along for the ride, so you maybe you should count “unique cache lines accessed” – but if you go far enough down that rabbit hole you end up building a complete machine model.

This, but also the effect of parallelization that is hard to take into account many times. Things that run in parallel are often for run by using little or no resources.

Looks like interleaved execution but done by hand. I don’t think that interleaving all operation one-by-one gives any benefit. Only the memory access in parallel results in a speedup.

There are recent publications using C++ coroutines for interleaved execution which provide a nicer abstraction and leave the algorithm mostly unchanged.

http://www.vldb.org/pvldb/vol11/p230-psaropoulos.pdf

There was also a recent CppCon talk about the topic:

https://www.youtube.com/watch?v=j9tlJAqMV7U

Thanks, the cppcon talk is well worth watching.

I understand your comment better now. Yes, this is even a lower level of interleaving than the approaches presented there, which break the work into resumable “tasks” – either by hand using an array of state objects, or using coroutines TS and then “schedule” these tasks round-robin using a more-or-less generic scheduler.

Here, the algorithm is simply modified to process within the main loop, more than one element at a time. It’s a less generic approach (in particular, you cannot easily vary the number of interleaved elements), but in principle it could be since the barrier between the scheduler and the tasks is removed. For example, either task based approach will probably never keep state values in registers even if there are few enough tasks that it would be feasible, but in this type hard-coded batched loop approach it is possible.

Now the relative results here are much worse: a speedup of ~1.4x vs more than 3x for co-routines, but this have more to do with Java vs C++ than any inherent deficiency.

There is still a fundamental difference: the number of computations between memory accesses.

In your shotgun approach you have all memory accesses right after each other. The next operations require the loaded values. There is not enough computation to load the first value from main memory, resulting in memory stalls. By grouping all operations before the memory access together and only issue a prefetch, we increase the number of computations before we use the loaded value. If we have N searches, (N-1) other searches execute all the operations leading to the next prefetch. We hide the memory latency with useful computation resulting in reduced CPU stalls.

In this example of binary search (that’s the same example as the linked CppCon talk), there is essentially no computation. Maybe something like 1 or 2 cycles of work (a shift, compare and branch) pretty much, compared to the 100s of cycles needed for the memory access. So the compared to the memory access time, the time needed for computation rounds down to 0.

So the goal here is not to overlap memory access with computation, but memory access with more memory accesses. That’s how the coroutine method gets its speedup.

In either Daniels code or the coroutine case, ignoring the complication of prefetching vs loading, the pattern is the same: a bunch of loads, then a bunch of “computation” (not actually a bunch at all). So having “all memory accesses right after each other” is fine (and coroutine does the same thing, at least if the arrays are the same size) – the only requirement is that the adjacent loads are

independentso they can all execute without waiting for any other one.Within a single search, the loads are all dependent – when considering the whole group of searches, there are many independent loads (i.e., each new search involves loads independent of all the other ones), but the CPU can’t “see” them in the middle of a single search because all in can see is a stream of dependent loads all part of the current search. Its internal resources will often fill up before it speculates into the next search (and even then you are only at 2 parallel searches: you want many more).

So interleaving is all about reorganizing the stream of accesses so any random window of say 50-200 instructions contains a lot of independent accesses. The code Daniel presents here and coroutines both accomplish that.

Looks like you are right about the computation. There is just not enough in a binary search. I did some quick&dirty measurements with the following code:

https://github.com/tetzank/interleaved-algorithms

Group prefetching (simpler version of what coroutines do) and plain interleaving (called shotgun in the article) have pretty much the same performance on my machine. I have not thoroughly profiled it and only briefly checked the assembly (gcc creates branchless code (cmov) for both). I also threw in an AVX2 version. On Haswell it’s slower as the gather instruction is not performant enough. Might be better on Skylake.

Just a testbed to get some numbers into the discussion. And, eventually, to try it on other memory latency bound algorithms.

I tested your SIMD code and it is good on Skylake, the 2-way version being at least twice as fast as the fastest “plain4” interleaved variant.

The 2-way SIMD version was nearly twice as fast as the 1-way SIMD version which was itself slightly faster than the fastest plain variants.

Well I think Daniel’s result here shows that it does give a benefit, right?

In particular, the parallel memory access is made

possible(or at least “more possible”) by interleaving since there are more independent memory operations visible in the out-of-order window of the processor. Typical binary search is kind of a worst case for memory parallelism since each access in the search depends on the prior one, either via a control dependence (with a probability of a wrong prediction and hence a pointless probe) or a data dependency (even worse since now all accesses are serialized).So yeah, “only the memory access in parallel results in a speedup”, but that’s like saying “only the sugar and chocolate in this brownie results in deliciousness” or something to that effect.

I will check out the work you linked – but my initial question is always “how much does the abstraction cost” – even a few instructions limits how fine-grained you can do the switching. So I am very interested in this stuff but still in “see it when I believe it” mode.

The problem is that with interleaving each operation you carry a lot of state around. If every value still fits in the registers, it’s probably faster than other approaches. But if it doesn’t, one loads and stores too often from and to the stack. Sure stack is in L1, but still. I would assume that it limits the number of searches you can interleave too much. On Intel CPUs, we can have up to 10 memory loads in flight. Getting close to 10 interleaved searches would be ideal.

Yes there is definitely a limit to the interleaving before code quality drops and you start spilling registers, and this limit depends in the number of registers available (obviously) and also the code in question. Sometimes you can come up with a level transformation that reduces the number of registers used.

In this example, I think you only need one additional register per interleaved search, plus a few other registers (whose number doesn’t depend on the number of searches) as temporaries and to track the current search span – so you could in principle get to your desired limit of 10 on an 16-reigster architecture.

As you point out, spilling isn’t necessarily all that bad if you fighting main memory latencies as your primary limiting factor, since spills to L1 take only a small fraction of the time of a miss to main memory (but spills also tend to clog up the out of order window with more instructions, which may also limit the achieved memory parallelism).

Coroutines don’t have any magic way around spills though, do they? They also have to dealing with saving and restoring registers, and I don’t see them having access to any tricks that the register allocator in a compiler would not (and on the contrary, the “single moment in time” switch implies additional restrictions which result in more spills.

But that’s always the case. Optimal code has to target the particular properties of the machine, and there are limits to how much you can hope for the compiler intuiting your goal and appropriately matching HW!

For example back in the day (early 90s) I wrote a DFT to handle AAC audio. Most practical books presented radix-2 code, and for x86 that was probably the optimal solution given the paucity of x87 registers. But for PPC I could just fit a radix-8 DFT into the 32 FP registers.

This seems to me analogous — optimal structure (ie optimal “interleaving”) for x86 may differ from optimal structure for “generic” ARM which will in turn differ from optimal structure for Apple ARM.