Faster Set Intersection with SIMD instructions
by Reducing Branch Mispredictions
Hiroshi Inoue†‡, Moriyoshi Ohara†, and Kenjiro Taura‡†IBM Research – Tokyo, ‡University of Tokyo{inouehrs, ohara}@jp.ibm.com, tau@eidos.ic.i.u-tokyo.ac.jp
rich history of research and many algorithms have been proposedto accelerate the operation. Many of these existing algorithms
focus on reducing the number of comparisons to improve theperformance. For example, special data structures, such as hashbased structures [2, 3] or hierarchical data structures [4], can beused to boost performance by skipping redundant comparisons,especially when the sizes of the two input sets are significantly
different. Unfortunately, such algorithms typically require
preprocessing of the input data to represent it in a special form or
to create additional data structures. Also, when the sizes of thetwo inputs are similar, it is much harder to achieve performance
improvements because there are far fewer obviously redundantcomparisons. As a result, simple merge-based implementations of
set intersections are still widely used in the real world, such as thestd::set_intersection implementation of STL in C++. Our goal is
to improve the performance of set intersection even when the two
input sets are of comparable size without preprocessing.
This paper describes our new algorithm to improve theperformance of the set intersection. Unlike most of the existing
advanced techniques, we focus on improving the execution
efficiency of the set intersection on the microarchitectures oftoday’s processors by reducing the branch mispredictions insteadof reducing the number of comparisons. Today’s processors areequipped with a branch prediction unit and speculatively execute
one direction (taken or not-taken) of a conditional branch tomaximize the utilization of processor resources. If the predicteddirection of the branch does not match the actual outcome of the
branch (branch misprediction), the hardware typically wastes
more than ten CPU cycles because it needs to flush speculatively
executed instructions and restart the execution from the fetch of
the next instruction for the actual branch direction.
In our algorithm, we extend the naive merge-based set
intersection algorithm by reading multiple elements, instead ofjust one element, from each of the two input arrays and compare
all of the pairs of elements from the two arrays to find any
elements with the same value. Figure 1 compares a naive
algorithm and our algorithm using 2 as the number of elements
read from each array at one time (block size). Surprisingly, thissimple improvement gave significant performance gains on bothXeon and POWER7+.
Our key insight for the performance improvement is that usinglarger block size yields a smaller number of hard-to-predict
conditional branches in the set intersection. The comparison to
select an input array for the next block, shown in bold in Figure 1,
will be taken in arbitrary order, and therefore it is very hard forbranch prediction hardware to predict the branches. However,most of the conditional branches in the all-pairs comparisons are
not-taken and they do not cause frequent branch mispredictions in
the most frequent cases of many applications, since the number of
output elements is generally much smaller than the number of
input elements for many applications. For example, Ding and
König [2] reported that 94% of the queries in a shopping portalsite select less than 10% of the input elements as the output.
Therefore using the larger block size can reduce the number of

ABSTRACTSet intersection is one of the most important operations for many
applications such as Web search engines or database management
systems. This paper describes our new algorithm to efficiently
find set intersections with sorted arrays on modern processors
with SIMD instructions and high branch misprediction penalties.
Our algorithm efficiently exploits SIMD instructions and can
drastically reduce branch mispredictions. Our algorithm extends a
merge-based algorithm by reading multiple elements, instead ofjust one element, from each of two input arrays and compares allof the pairs of elements from the two arrays to find the elementswith the same values. The key insight for our improvement is that
we can reduce the number of costly hard-to-predict conditional
branches by advancing a pointer by more than one element at atime. Although this algorithm increases the total number ofcomparisons, we can execute these comparisons more efficiently
using the SIMD instructions and gain the benefits of the reduced
branch misprediction overhead. Our algorithm is suitable toreplace existing standard library functions, such as
std::set_intersection in C++, thus accelerating many applications,because the algorithm is simple and requires no preprocessing to
generate additional data structures. We implemented ouralgorithm on Xeon and POWER7+. The experimental results
show our algorithm outperforms the std::set_intersection
implementation delivered with gcc by up to 5.2x using SIMD
instructions and by up to 2.1x even without using SIMD
instructions for 32-bit and 64-bit integer datasets. Our SIMD
algorithm also outperformed an existing algorithm that canleverage SIMD instructions.

1. INTRODUCTION
Set intersection, which selects common elements from two input
sets, is one of the most important operations in many applications,
including multi-word queries in Web search engines and join
operations in database management systems. For example, in Web
search engines the set intersection is heavily used for multi-word
queries to find documents containing two or more keywords byintersecting the sorted lists of matching document IDs from theindividual query words [1]. In such systems, the performance ofthe sorted set intersection often dominates the overall
performance. Due to its importance, sorted set intersection has aThis work is licensed under the Creative Commons AttributionNonCommercial-NoDerivs 3.0 Unported License. To view a copy of thislicense, visit http://creativecommons.org/licenses/by-nc-nd/3.0/.
Obtainpermission prior to any use beyond those covered by the license. Contact
copyright holder by emailing info@vldb.org. Articles from this volume wereinvited to present their results at the 41st International Conference on Very
Large Data Bases, August 31st – September 4th 2015, Kohala Coast, Hawaii.
Proceedings of the VLDB Endowment, Vol. 8, No. 3
Copyright 2014 VLDB Endowment 2150-8097/14/11

293

hard-to-predict conditional branches in exchange for a larger
number of total comparisons and easy-to-predict conditional
branches. When the block size is S, one all-pairs comparison
requires up to S2 conditional branches, but the number of hard-topredict conditional branches is reduced to only one per S input
elements.
The larger number of comparisons from these all-pairs
comparisons is an obvious drawback of our algorithm. However
we can effectively eliminate many of these additional comparisonsby using SIMD instructions. By combining our algorithm with
SIMD instructions, we can reduce the number of branchmispredictions without increasing the total number of executed
instructions by executing multiple comparisons in parallel.
Our algorithm roughly doubled the performance for set
intersection for 32-bit and 64-bit integer datasets even without
using SIMD instructions compared to the std::set_intersection
implementation delivered with gcc. The use of SIMD instructions
further doubled the performance on both processors. We use
SIMD instructions to filter out redundant scalar comparisons byusing only a part of each element instead of finding matchingpairs directly with SIMD comparisons. This approach increases
the data parallelism within each SIMD instruction and leads tohigher performance. It also allows us to use SIMD instructions if
the data type is not natively supported by the SIMD instruction set,e.g. 64-bit integers on POWER7+.
Our new algorithm seeks to accelerate the set intersection
operation when (1) the number of output elements is muchsmaller than the number of input elements and (2) the sizes of thetwo input sets are not significantly different. For datasets that do
not satisfy these assumptions, other algorithms such as binarysearch-based algorithms can outperform our algorithm. This is
why we devised a practical technique to adaptively select the bestalgorithm based on the ratio of the number of output elements
over the number of input elements (selectivity) and the size ratioof two input sets.
The rest of the paper is organized as follows. Section 2
discusses related work. Section 3 describes our technique toreduce branch mispredictions. Section 4 gives a summary of ourresults. Finally, Section 5 summarizes our work.

(a) Schematic and pseudocode for naive merge-based algorithm
input array A

2

5

6

8

3

5

9

1
input array B

1

2
1
2
3
4
5
6
7
8
9
10
11
12

in each iteration
1. to check the equalityof two elements
2. to advance a pointer by 1

while (Apos < Aend && Bpos < Bend) {if (A[Apos] == B[Bpos]) {C[Cpos++] = A[Apos];Apos++; Bpos++;}
else if (A[Apos] > B[Bpos]) {Bpos++;}
else {
Apos++;}}(b) Schematic and pseudocode for our algorithm (block size = 2)
input array A

2

input array B

1

5

6

3

5

8

1

9
2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

2. RELATED WORK
Sorted set intersection has a rich history of research and manyalgorithms have been proposed to improve its performance.
Previous studies, such as [5], showed that simple algorithms
based on a linear merge (e.g. Figure 1(a)) performed best when
the sizes of the two input sets were similar. Another advantage ofa merge-based algorithm is that it works well with input datasetsthat are compressed with various algorithms, such as delta
encoding, due to its simplicity. Our algorithm extends the mergebased algorithm and improves its efficiency on today’s processors.
This means our algorithm inherits the advantages of the mergebased algorithms and performs well on input sets of similar size.
When one input set is much larger than another, the mergebased algorithm is less efficient and many techniques have been
proposed to improve its performance. For example, algorithms
based on binary search [5-7] reduce the numbers of comparisonsand memory accesses by picking values from the smaller set, andefficiently searching for the matching value in the larger set.
Similarly, hash-based techniques [2, 3] or techniques using
hierarchical data representations [3] improve the performance byreducing the number of comparisons. However, most of thesetechniques are effective only when the sizes of the two input sets
are significantly different (as by an order of magnitude or more).
Our algorithm focuses on improving the performance of the input

in each iteration:
1. to find any matching pairsin blocks of S elements,
here S=2(all‐pairs comparison)
2. to advance a pointer by S

while (1) {Adat0 = A[Apos]; Adat1 = A[Apos + 1];
Bdat0 = B[Bpos]; Bdat1 = B[Bpos + 1];
if (Adat0 == Bdat0) {C[Cpos++] = Adat0;}
else if (Adat0 == Bdat1) {C[Cpos++] = Adat0;
goto advanceB;
all‐pairs}
comparison
else if (Adat1 == Bdat0) {(using S2
C[Cpos++] = Adat1;
goto advanceA;
comparisons)}if (Adat1 == Bdat1) {C[Cpos++] = Adat1;
goto advanceAB;}
else if (Adat1 > Bdat1) goto advanceB;
else goto advanceA;
advanceA:
Apos+=2;
if (Apos >= Aend) { break; } else { continue; }advanceB:
Bpos+=2;
if (Bpos >= Bend) { break; } else { continue; }advanceAB:
Apos+=2; Bpos+=2;
if (Apos >= Aend || Bpos >= Bend) { break; }}// fall back to naive algorithm for remaining elements

Figure 1. Overview of set intersection without and with ourtechnique. Red bold (line 6 in (a) and line 19 in (b)) shows ahard-to-predict conditional branch used to select the pointer toadvance.
sets when the sets are roughly the same size. We can easily takeadvantages of both our algorithm and these existing techniques byselecting a suitable algorithm based on the size of the input sets
before executing the operations as we describe in Section 4.
Recently, Ding and König [2] proposed a hash-based
technique that is effective even for input sets of similar size. They
partitioned the input set into small partitions and encoded thevalues in each partition into a machine-word size bitmap in the
preprocessing phase to efficiently calculate the intersection using
bit-wise AND instructions at run time. Though both theiralgorithm and ours are effective for arrays of similar sizes, anadvantage of our algorithm is that we do not need preprocessing
for additional data structures before an intersection operation.

294

Table 1. Summary of set intersection algorithms
the sizes of the two
input sets (Na, Nb)

Similar size

Significantly
different(more than 10x)

without SIMD instructions

With SIMD instructions

Merge-based algorithms (e.g. STL): is simple and widelyused, but suffer from branch mispredictions.
If-conversion: eliminates the branch mispredictions in tradefor longer path length. is hard to SIMDize.
Ding and König [2]: requires preprocessing.
Our block-based algorithm: is an extention to the mergebased technique, reduces the branch mispredictions
Binary-search-based techniques (e.g. galloping [10]): do
not require preprocessing
Techniques using additional data structures (skip list,
hash etc): require preprocessing

Schlegel [8]: is mainly targetting 8-bit or 16-bit integers
supported by STTNI instructions
Lemire [9] (V1 algorithm): can operate on 32-bit integers
supported by using SIMD comparisons
Our block-based algorithm with SIMD: yields larger gain
by increased data parallelism. can support data types, evennot natively supported by the SIMD instruction set.

Schlegel et al. [8] exploited a special instruction called STTNI(STring and Text processing New Instruction) included in the
SSE 4.2 instruction set of recent Intel processors for sorted set
intersections. They compared multiple values read from each
input array by using the special instruction to execute an all-pairs
comparison in one step. They showed up to 5.3x and 4.8x
accelerations for 8-bit and 16-bit data, respectively, on a Core i7
processor. Because the STTNI instruction does not support datatypes larger than 16 bits, for 32-bit integer data, they intersect theupper 16 bits and lower 16 bits separately. This techniqueachieved good speedups only when the value domain was limited,
so that a sufficient number of elements shared the same value intheir upper 16 bits, which limits the real-world use of thisapproach. Our algorithm does not have this limitation on the value
domain even when we use the SIMD instructions because we use
SIMD instructions as a filter to reduce scalar comparisons insteadof finding matching pairs directly with SIMD comparisons. Also,
our non-SIMD algorithm does not use any unique instructions,which makes it more portable among processors.
Lemire et al. [9] accelerated decompression and set
intersection used in index search systems by SIMD instructions.
For the set intersection, they used three algorithms designed for
SIMD instructions; so called V1 algorithms, V3 algorithm andgalloping [10] with SIMD. They selected the algorithm based on
the ratio of sizes of the two input arrays. When the sizes of the
input sets are of similar size, the V1 algorithm performed best
among the three algorithms. When the lengths of the two input
arrays are significantly different, SIMD galloping was the best.
Our algorithm focuses on the case when the sizes of the two
inputs are similar. As showed in Section 4, our SIMD algorithm
achieved about 2x better performance than their V1 SIMD
algorithm. Also our SIMD algorithm can be used even if the data
type is not natively supported by the SIMD instruction set, suchas 64-bit integers on POWER7+, while their technique requires
special handling for such case. When the gap between the two
input sizes becomes large (by more than an order of magnitude),
we switch to their SIMD galloping algorithm.
Schlegel et al. [8] also showed that replacing the branch
instruction in the merge-based set intersection with the predicatedinstructions of the Intel processor improved the performance overthe branch-based implementation, though this was not the main
focus of their work. This optimization, which replaces control
flow by data flow, is called if-conversion [11], and is aconventional technique to reduce branch misprediction overhead.
Figure 2 shows an example of set intersection implementation
without using a conditional branch to advance a pointer. As
shown in the performance comparisons later, our techniqueyielded much better performance than a branchless
implementation by the if-conversion approach on both of thetested processors. Although both techniques reduce the branchmisprediction overhead, our algorithm achieved better

1
2
3
4
5
6
7
8
9
10
11
12

Lemire [9] (SIMD galloping): is an extention to the gallop
algorithm [10] with SIMD instructions

while (Apos < Aend && Bpos < Bend) {Adata = A[Apos];
Bdata = B[Bpos];
if (Adata == Bdata) { // easy-to-predict branch
C[Cpos++] = Adata;
Apos++; Bpos++;}
else { // advance pointers without conditional branches
Apos += (Adata < Bdata);
Bpos += (Bdata < Adata);}}

Figure 2. Pseudocode for set intersection without hard-topredict conditional branches (the branchless algorithm).
performance even without SIMD instructions by having a shorter
path length. Table 1 highlights these existing and our new
algorithms.
We can improve the performance of the set intersection by
reducing the branch mispredictions. Branch misprediction
overhead is recognized as a major performance constraint in someworkloads. For example, Zhou et al. [12] reported that theperformance of many database operations, such as B+ tree searchand nested loop joins for unsorted arrays, can be improved by
reducing the branch misprediction overhead using SIMD
instructions. Also, some sorting algorithms [13-15] improved theperformance in sorting random input by reducing the branchmisprediction overhead using the predicated instructions or SIMD
instructions.
Some set intersection algorithms exploit other characteristicsof the modern hardware architectures beyond their branchperformance, such as the large amount of cache memory andmultiple cores [16-18]. For example, Tsirogiannis et al. [16]improve the performance of a search-based intersection algorithm
by using a small index called a micro-index, which can fit into theprocessor’s cache memory.

3. OUR ALGORITHM
In this section, we first describe our block-based set intersection
algorithm that reduces branch mispredictions with two sorted andunique-element arrays without using SIMD instructions. Then weshow how we exploit SIMD instructions to execute multiple
comparisons in parallel to further improve the performance.
Although our algorithm is for intersecting two input sets, we canintersect multiple sets by using this algorithm as a building block,repeatedly executing our algorithm for the two smallest among allof the input data sets, which is a frequently used technique.
3.1

Key Observation and Assumptions for Input Data

As shown in Figure 1(a), a naive merge-based algorithm for set
intersection reads one value from each of the input arrays (A and
B in the Figure) and compares the two values. If the two values
are equal, the value is copied into the output array (C) and the295pointers for both arrays are advanced to load the next values. If
the two values are not equal, the pointer for the array whose
element is smaller is advanced. This approach requires up to (Na +Nb - 1) iterations to complete the operation. Here Na is the numberof elements in the array A and Nb is the number of elements in thearray B. Each iteration includes one if_equal conditional branch(Line 2 in Figure 1(a)) and one if_greater conditional branch(Line 6 in Figure 1(a)).
Here, the if_greater conditional branches in the set
intersection are hard to predict and incur significant branchmisprediction overhead, while the if_equal conditional branches
rarely cause mispredictions. Hence our algorithm focuses on
reducing the number of costly if_greater conditional branches atthe cost of using more if_equal conditional branches to optimize
for the common case. The if_greater conditional branches cause
frequent mispredictions because they will be taken in arbitraryorder with roughly 50% probability when the sizes of the input
sets are similar. This makes it very hard for branch prediction
hardware to predict the branch directions correctly. Most of theconditional branches in the all-pairs comparisons are not-taken
and they do not cause frequent branch mispredictions in typicalcases for many applications, since it is known that the number of
output elements is much smaller than the number of input
elements in practice [2].
To achieve its speedup, our algorithm assumes that:

(a)

A[0]

A[1]

A[2]

A[3]

A[4]

A[5]

input array A

1

5

6

8

11

14

2

3

5

9

10

12

B[0]

B[1]

B[2]

B[3]

B[4]

B[5]

(b)

A[0]

A[1]

A[2]

A[3]

A[4]

A[5]

input array A

1

5

6

8

11

14

input array B

2

3

5

9

10

12

B[0]

B[1]

B[2]

B[3]

B[4]

B[5]

(c)

A[0]

A[1]

A[2]

A[3]

A[4]

A[5]

input array A

1

5

6

8

11

14

input array B

2

3

5

9

10

12

B[0]

B[1]

B[2]

B[3]

B[4]

B[5]

(d)

A[0]

A[1]

A[2]

A[3]

A[4]

A[5]

input array A

1

5

6

8

11

14

input array B· the number of output elements is much smaller than the numberof input elements, andinput array B· the sizes of input sets are not significantly different (as by anorder of magnitude or more).

3

5

9

10

12

B[1]

B[2]

B[3]

B[4]

B[5]

( {A[0], A[1]}with {B[0], B[1]} )2. advance the pointerfor the arraywith smaller value
3. compare next twovalues from each
input array( {A[0], A[1]}with {B[2], B[3]} )4. find a matching pairin third comparisonand copy the value
into the output array
5. skip the fourthcomparison in the
all‐pair comparisons(A[1] with B[3])and advancethe pointer for A
because B[3] mustbe larger than A[1](A[1] = B[2] < B[3])Figure 3. An example of set intersection of two sorted arrays with
our technique using 2 as the block size.

Our algorithm performs well for datasets that satisfy theseassumptions. Otherwise, we adaptively switch to anotheralgorithm to combine the advantages of our algorithm with thestrengths of such algorithms as the binary-search-basedalgorithms. The first assumption is important to avoid frequent
mispredictions in the if_equal conditional branches, which we
assume are not costly. In Section 4, we describe an adaptive
fallback technique to validate the first assumption at runtime. The
second assumption ensures that the if_greater conditional
branches cause lots of mispredictions. We select the bestalgorithm and a parameter (block size) based on the sizes on thetwo input sets before executing the operations. If the sizes of thetwo input arrays are significantly different, we switch to a binarysearch-based algorithm.
3.2

2
B[0]

1. search matchingpairs in two values
from each input arrayby comparing all pairs

Figure 3 shows a step-by-step example of our algorithm with theblock size of 2. Because the block size is 2, the all-pairs
comparison shown in Figure 3 uses up to 4 if_equal conditional
branches (Lines 4-18 in Figure 1(b)). Then, in Figure 3(b), we
compare the second value from each array (A[1] = 5 and B[1] =3) and advance the pointer for the array B because B[1] is smaller
than A[1]. This step uses only one if_greater conditional branch(Line 19 in Figure 1(b)). Typically, this if_greater conditional
branch is hard to predict and hence it causes frequent branchmispredictions. Our algorithm aims to reduce the number of thehard-to-predict conditional branches executed. If there is nomatching pair found in the all-pairs comparison (the mostfrequent case), then we repeat the steps shown in Figures 3(a) and3(b). If a matching pair is found in an all-pairs comparison(Figure 3(c)), then we copy the value (5 in the figure) into the
output array and skip the following comparisons, which are nolonger possible matches. In the Figure, B[3] (= 9) must be larger
than B[2] and A[1] (= 5) because each array was sorted, andhence we can advance the pointer for the array A without
comparing A[1] and B[3] (Figure 3(d)).
If the number of total elements in an input array is not a
multiple of the block size S, we can just fall back to the naive
approach for the remaining elements. If the number of elements inan array is large enough, this does not measurably affect theoverall performance of the set intersection operation.

Our Basic Approach without SIMD instructions

Our technique extends the naive merge-based algorithm shown in
Figure 1(a) by reading multiple values from each input array andcompares all of their pairs using if_equal conditional branches asshown in Figure 1(b). We call the number of elements comparedat one time the block size (S). We repeat the following steps until
we process all of the elements in the two input arrays (A and B):(1) read S elements from each of two input arrays,(2) compare all possible S2 pairs of elements (e.g. four pairs in
Figure 1(b), where S = 2) by using if_equal conditional
branches to find any matching pairs,(3) if there is one or more matching pairs, copy the value orvalues of the found pairs into the output array (C),(4) compare the last elements of the two arrays used in Step 2,(5) advance the pointer by S elements for the array whose lastelement is smaller in Step 4.

Identifying the best block size: With our scalar (non-SIMD)
algorithm, the total number of if_equal conditional branches
increases for larger block sizes. The number of if_equal
conditional branches, which are typically easy to predict, involvedin one all-pairs comparison is up to S2. To complete the set
intersection for the entire array, we need to execute the all-pairs296

Table 2. Summary of the number of conditional branches without SIMD instructions using the same block size S for both input arrays
Approach

Naive algorithm
Figure 1(a)Our algorithm
Figure 1(b)Number of hard-topredict if_greater
conditional branches(may cause frequent
mispredictions)

Number of if_equal
conditional branches(mispredictions
infrequent)

Total number ofconditional branches

up to Na+Nb - 1

up to Na+Nb - 1

up to 2 (Na+Nb - 1)

up to (Na+Nb)/S - 1
about S times lessup to S·(Na+Nb) - S2
about S times more

up to (S2+1)·((Na+Nb)/S - 1)
about (S2+1)/S times more

Performance
characteristics- shorter path length- larger misprediction overhead- longer path length- smaller misprediction overhead

Na, Nb: the number of elements in the two input arrays A and B. S: the block size.
comparisons up to ((Na + Nb)/S – 1) times, so the total number of
if_equal conditional branches is up to S2·((Na + Nb)/S – 1) = S·(Na+ Nb) - S2. This is almost S times larger than the naive approach.
At the same time, the number of if_greater conditional
branches, which are more costly than if_equal conditional
branches due to their frequent branch mispredictions, decreases as
the block size increases. We advance the pointer by S elements at
a time instead of by just one element as in the naive approach.
Hence, the number of if_greater conditional branches is ((Na +Nb)/S – 1). This is almost S times smaller compared to the naive
approach. We summarize the total numbers of comparisons inTable 2. Because our algorithm is a generalization of the naive
algorithm, which is equivalent to our algorithm with a block sizeof 1, the number of comparisons in Table 1 is the same for bothalgorithms when S = 1.
The key parameter to find the best block size is the penalty of abranch misprediction compared to the number of CPU cycles to
execute a conditional branch without a misprediction. To find thebest block size, we calculate the cost of total branch instructions
including the misprediction overhead. We assume that only the
hard-to-predict conditional branches cause mispredictions fortypical input. The best block size with this assumption is the S
that minimizes this cost function f (S):f (S) = branch_cycles × (S2+1) / S +mpred_penalty × mpred_ratio / S,array may give additional performance gains, especially when thesizes of the two input arrays are very different. When we use
different block sizes Sa and Sb for each input array, the number of
if_greater conditional branches is up to Na/Sa + Nb/Sb, and thenumber of if_equal conditional branches is (Na/Sa + Nb/Sb) ×(Sa·Sb). Hence the cost function becomesf (Sa, Sb) = branch_cycles × (Sa·Sb+1) × (q/Sa+(1-q)/Sb) +mpred_penalty × mpred_ratio × (q/Sa+(1-q)/Sb). (2)Here, q shows how the sizes of two input arrays different, q =Na/(Na+Nb). When Sa = Sb, equation (2) becomes equivalent withequation (1) regardless of q. When the sizes of two input arrays
are significantly different, the misprediction rate at the hard-topredict conditional branch is much difficult to estimate anddepends on the branch predictor implementation. When we
assume a simple predictor, which just predicts the more frequentside of the two branch directions, taken or not-taken, themisprediction rate is min(q, 1-q). With this assumption and themisprediction penalty of 16 cycles, for example, the best block
sizes for two arrays with sizes of Na and 2Na, i.e. q = 1/(1+2) =1/3, are Sa = 2 and Sb = 4, while Sa = Sb = 3 is the best if the two
input arrays have the same size, i.e. q = 1/2, as already discussed.
Our scalar algorithm selects the block sizes before starting theoperation based on the ratio of the sizes of the two input arrays.
We use Sa = 2 and Sb = 4 if the size of the larger array Nb is morethan twice the size of the smaller array Na. Otherwise, we use Sa =Sb = 3. We show how the block size affects the performance inSection 4.

(1)Here, branch_cycles is the number of cycles to execute a
conditional branch and mpred_penalty is the penalty of a branchmisprediction (in terms of cycles). The mpred_ratio is the branchmisprediction ratio for the hard-to-predict conditional branches.
We assume the mpred_ratio is 50% when the sizes of the two
input set are comparable because the hard-to-predict conditional
branches are taken in arbitrary order and hence no branchpredictor can predict them correctly. When the misprediction
penalty is more than twice the cost of a successfully-predictedconditional branch, our technique improves the performance overthe naive algorithm by using a block size of 2. When the relative
cost of the misprediction is between 10 and 22, as is true for manyof today’s processors, the best block size is 3. The branchmisprediction penalties for POWER7+ and Xeon were both about
16 cycles as measured with a micro-benchmark and the cost of abranch instruction is expected to be 1 cycle. We predict that theblock size of 3 yields the best performance and the block size of 4
is a close second best. We empirically confirmed this estimate in
Section 4. Because the best block size also depends on the input
data and not just the processor, we used an adaptive control
technique with a runtime check to detect pathological
performance cases for our algorithm.

3.3

Exploiting SIMD Instructions

Our algorithm reduces the branch misprediction overhead butwith an increased number of easy-to-predict conditional branches,as discussed in Section 3.2. To further improve the performance,
we employ SIMD instructions to reduce the number of
instructions by executing the all-pairs comparisons within eachblock in parallel. Unlike the previous SIMD-based set intersection
algorithms, we use SIMD instructions to filter out unnecessary
scalar comparisons by comparing only a part of each element.
This approach allows us to use SIMD instructions if the data type
is not natively supported by the SIMD instruction set. For
example, we can use processors without 64-bit integer
comparisons in their SIMD instructions to intersect 64-bit integer
arrays, e.g. 64-bit integers on POWER7+. Also it increases thedata parallelism within one SIMD instruction by using only a partof the elements.
For our SIMD algorithm, we used a multiple of four as theblock size S (or Sa and Sb when using different block sizes for twoarrays) so we could fully exploit the SIMD instructions of theprocessors, which can compare up to 16 or more data pairs of 1byte elements in parallel with one instruction. Our SIMD
algorithm selects Sa = Sb = 4, if the size of the larger array Nb isless than twice the size of the smaller array Na. Otherwise, we use

Using different block sizes for each input array: Up to now, wehave been assuming that we use the same block size S for both
input arrays. However, using different block size for each input

297

A[i]

A[i+1]

A[i+2]

A[i+3]

input array A
vector register A(128 bit)
vector register B(128 bit)compare eachbyte pair

input array B
B[j]B[j+1]A[i]A[i+1]

B[j+2]

input array A
vector register A(128 bit)
vector register B(128 bit)compare eachword pair

input array B
B[j]B[j+1]

B[j+3]

Byte‐wise check‐ check four elementsin parallel‐ gather one byte fromeach element as
illustrated and compare
each byte pair‐ store result ina vector register as
bitmask

Word‐wise check‐ check two elementsin parallel‐ gather one word fromeach element as
illustrated and compare
each word pair‐ store result ina vector register as
bitmask

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

while (1) {// do byte-wise check: step (1)-(3) of the hierarchical checkcompare_result = bytewise_check1(&A[Apos], &B[Bpos]) &bytewise_check2(&A[Apos], &B[Bpos]);
if (!is_all_bit_zero(compare_result)) { // step (4)// found a potential matching value// do word-wise check and scalar check: step (5) - (7)...
}else if (A[Apos+3] > B[Bpos+3]) goto advanceB;
else goto advanceA;
advanceA:
Apos+=4;
if (Apos >= Aend) { break; } else { continue; }advanceB:
Bpos+=4;
if (Bpos >= Bend) { break; } else { continue; }advanceAB:
Apos+=4; Bpos+=4;
if (Apos >= Aend || Bpos >= Bend) { break; }}Figure 5. Pseudocode of our SIMD algorithm for block size of 4x4.
(1) Do a byte-wise check for A[i .. i+3] and B[j .. j+3] using theleast significant byte,(2) Do a byte-wise check for A[i .. i+3] and B[j .. j+3] using thesecond-least significant byte,(3) Do a bit-wise AND operation for the results of Steps 1 and 2,(4) If every bit is zero in Step 3, then skip further checks because
there is no matching pair (most frequent case),(5) Do a word-wise check for A[i .. i+1] and B[j .. j+1] using thethird to sixth bytes,(6) Do a scalar check for Step-5 matches, and(7) Repeat Steps 5 and 6 for {A[i .. i+1] and B[j+2 .. j+3]},{A[i+2 .. i+3] and B[j .. j+1]}, and {A[i+2 .. i+3] and B[j+2 ..
j+3]}.

Figure 4. Overview of byte-wise check and word-wise check.
Sa = 4 and Sb = 8. As discussed in Section 3.2, the block size of 3
is best for our scalar (non-SIMD) algorithm. However we canreduce the number of comparisons by using SIMD instructions,which justifies using a larger block size than used in the scalar
algorithm.
The vector sizes of the SIMD instruction sets of today’s
processors, such as SSE/AVX of Xeon or VSX of POWER7+, arelimited to 128 bits or 256 bits. This means we can execute only
two or four comparisons of 64-bit elements in parallel. This
parallelism is insufficient for an all-pairs comparison of large
blocks in one step. To execute the all-pairs comparisons for larger
blocks efficiently by increasing the parallelism available in one
SIMD instruction, we use a parallel SIMD comparison, which
compares only a part of each element, to filter out all of the values
with no outputs before executing the all-pairs comparison using
scalar comparisons. Unlike the previous SIMD approaches [8, 9],
we did not fully replace the scalar comparisons with parallel
SIMD comparisons. Because the number of matching pairs are
typically much smaller than the number of input elements inpractice, our filtering technique is effective to avoid most of thescalar comparisons and hence achieves higher overall
performance. Zhou et al. [12] also used a similar idea ofcomparing only a part of each element to increase the data
parallelism with one SIMD instruction for the nested loop join for
unsorted arrays.
We use two different types of checks to find a matching pair in
the all-pairs comparison hierarchically. Figure 4 shows an
overview of our byte-wise check and word-wise check using
SIMD parallel compare instructions. In this example, we assume
64-bit integers as the data type and an SIMD instruction set with
128-bit registers. These checks execute only a partial comparisonof each element. It means that if the check does not find any
matching byte or word pair, there cannot be any matching element
pairs (no false negatives). However, if the check finds a matching
byte or word pair, the matching pair may still be a false positive.
To reduce the number of false-positive matches, we hierarchicallydo two different types of checks. If the data type of each element
is a 64-bit integer and the block size S is 4, our hierarchicalfiltering uses these steps:

Alternatively, we can replace Steps 5 to 7 with a countleading-zero instruction to identify the location of the matchingpairs found in Step 3. When we use a 32-bit integer data type, we
use the first to fourth bytes in Step 5. Figure 5 shows thepseudocode for the set intersection algorithm with our hierarchicalfiltering with SIMD.
On Xeon, we can use the STTNI (pcmpestrm) instruction,
which is unique to the Xeon processor, to execute the all-paircomparison efficiently. This instruction can execute the all-paircomparisons of eight 2-byte characters in one vector register
against eight characters in another vector register with only one
instruction. Thus we can implement Steps 1 to 3 with a block sizeof 8 by using the STTNI instruction very efficiently. Unlike
Schlegel’s algorithm [8], our algorithm uses the STTNI
instruction to filter out redundant comparisons and thus we do notneed to limit the data types to the 8- or 16-bit integer supportedby this instruction. When we use the STTNI in our SIMD
algorithm, we use the popcnt instruction to identify the position ofmatching pair efficiently because the processor does not support
the count-leading-zero instruction. We can get the position of theleast significant non-zero bit in the result of STTNI, x, by
popcnt((~x) & (x-1)).
When the SIMD instruction set supports a wider vector, suchas a 256-bit vector in AVX, one way to exploit the wider vector isdoing multiple byte-wise checks at in once step. For example, we
could do Steps 1 and 2 with just one parallel comparison using
256-bit vector registers, with 16 pairs of 2-byte elements.
By using our hierarchical filtering with SIMD instructions, we
can avoid increasing the number of instructions and still gain thebenefits of the reduced branch mispredictions.

298

using STTNI instruction unique to Xeon

1,600

higher is faster

performance (input entries / usec)

1,800
1,400
1,200
1,000
800
600
400
200
0

Xeon
(for 32-bit integers)
STL
Our scalar algorithm (block size = 3)
Our SIMD algorithm (block size = 4)
Our algorithm w/ STTNI (block size = 8)POWER7+(for 32-bit integers)
Naive (block size = 1)
Our scalar algorithm (block size = 4)
Branchless algorithm
Schlegel's algorithm w/ STTNI [8]POWER7+(for 64-bit integers)
Our scalar algorithm (block size = 2)
Our scalar algorithm (block size = 5)V1 SIMD algorithm (Lemire et al. [9])Figure 6. Performance for set intersection of 32-bit and 64-bit random integer arrays of 256k elements on Xeon and POWER7+. The
selectivity was set to 0 (as the best case). The error bars show 95% confidence intervals.
results shown as naive in the figures are the performance of thecode shown in Figure 1(a). As already discussed, our algorithm isequivalent to the naive algorithm when the block size is 1. We
also evaluated and compared the performances of the existingalgorithms including the widely used std::set_intersection librarymethod in the STL delivered with gcc, the branchless algorithm
shown in Figure 2, a galloping algorithm [10] (as a popularbinary-search-based algorithm), the two SIMD algorithms by
Lemire et al. [9] (which are the V1 SIMD algorithm and an SIMD
galloping algorithm), and Schlegel’s algorithm that uses the
STTNI instruction of Xeon [8]. We picked these algorithms forthe comparisons because, like our algorithm, they need no
preprocessing. Among these evaluated algorithms, the two
galloping algorithms based on binary searches are tailored forpaired arrays of very different sizes. The other algorithms,
including ours, are merge-based algorithms, which are known to
work better when the sizes of the two inputs are similar.

4. EXPERIMENTAL RESULTS
We implemented and evaluated our algorithm on Intel Xeon and
IBM POWER7+ processors with and without using SIMD
instructions. On Xeon, we also evaluated our SIMD algorithm
implemented using the Xeon-only STTNI instruction. We
implemented the program in C++ using SSE intrinsics on Xeon
and Altivec intrinsics on POWER7+, but the algorithm is thesame for both platforms. The POWER7+ system used for ourevaluation was equipped with a 4.1-GHz POWER7+ processor.
Redhat Enterprise Linux 6.4 was running on the system. We
compiled all of the programs using gcc-4.8.3 included in the IBM
Advance Toolkit with the –O3 option. We also evaluated theperformance of our algorithm on a system equipped with two 2.9GHz Xeon E5-2690 (SandyBridge-EP), also with Redhat
Enterprise Linux 6.4 as the OS, but the compiler on this system
was gcc-4.8.2 (still with the –O3 option) on this Xeon system. We
disabled the dynamic frequency scaling on both systems for more
reproducible results.
In the evaluation, we used both artificial and more realistic
datasets. With the artificial datasets generated using a random
number generator, we assessed the characteristics of our algorithm
for three key parameters, (1) the ratio of the number of output
elements over the input elements (selectivity), (2) the difference inthe sizes of the two input arrays and (3) the total sizes of the input
arrays. We define the selectivity as the number of output elements
over the number of elements in the smaller input array. To create
input datasets with a specified selectivity, we first generate a longenough array of (unsorted) random numbers without duplicates.
We then trim two input arrays from this long array with thespecified number of elements included in both arrays. Each arrayis then sorted and the pair is used as an input for experiments. We
executed the measurements 16 times using different seeds for the
random number generator and averaged the results. For therealistic datasets, we used arrays generated from Wikipedia data.
We generated a list of document IDs for the articles that include a
specified word. Then we executed the set intersection operations
for up to eight arrays to emulate the set intersection operation forthe multi-word queries in a query serving system. We also
averaged the results from 16 measurements for the real-world data.
We show the performance of our block-based algorithm withand without SIMD instructions and with various block sizes. The

4.1

Performance Improvements from Our Algorithm

Figure 6 compares the performance of the set intersection
algorithms for two datasets of 256k integers based on a random
number generator. The selectivity was set to zero. We used 32-bit
integer as the data type for both Xeon and POWER7+ and also
tested 64-bit integers for POWER7+.
The results show that our block-based algorithm improved theperformance over the naive merge-based algorithms (STL andnaive) on both platforms even without using the SIMD
instructions. When comparing how the block size affected theperformance of our scalar algorithm on these two platforms, thebest performance was when the block size was set to 3. On both
platforms, the block sizes of 3 and 4 gave almost comparable
performance. Our prediction based on the simple model discussed
in Section 3.2, which predicts the block size of 3 is the best and 4
is a close second best, seems reasonably accurate for bothprocessors, although they have totally different instruction sets
and implementations. The performance gains over the widely used
STL were 2.1x on Xeon and 1.8x on POWER7+ with the block
size of 3. Compared to our algorithm, the branchless algorithm
did not yield large performance gains over STL, although it
caused a smaller number of branch mispredictions than ouralgorithm (as shown later).

299

When we used the SIMD instructions, there were additional
performance improvements of about 2.5x over our scalar
algorithm on both platforms, where the total improvement was
4.8x to 5.2x better than STL for 32-bit integers and 4.2x better for
64-bit integers. Although the V1 SIMD algorithm [9] alsoachieved performance improvements over the STL using SIMD
instructions, the performance of our SIMD algorithm was 2.9x
and 3.0x better than V1 algorithm on Xeon and POWER7+
respectively. Also, the V1 algorithm cannot support 64-bit integer
on POWER7+ because POWER7+ does not have SIMD
comparisons for 64-bit integers, while our SIMD algorithm
achieved good performance improvements even for 64-bit integers
on POWER7+. This is because the V1 algorithm uses the SIMD
comparison for the entire elements to find the matching pairs, but
our algorithm uses the SIMD comparisons to filter out
unnecessary scalar comparisons by comparing only a part of eachelement. Schlegel's algorithm [8] did not achieve good
performance for the artificial datasets generated by the random
number generator. As the authors noted, Schlegel's algorithm isefficient only when the value domain is limited, so that asufficient number of elements share matching values in their
upper 16 bits, and this is not true for our artificial datasets.

0.6

Shorter is better

branch mispredictions per input element

0.7

0.5
0.4
0.3
0.2
0.1
0.0
Xeon

POWER7+

Cycles Per instruction (CPI)
2.0

Microarchitectural Statistics

For more insight into the improvements from our algorithm withand without SIMD instructions, Figure 7 displays somemicroarchitectural statistics of each algorithm for the artificial
datasets in the 32-bit integer arrays as measured by the hardware
performance monitors of the processors. We studied the branchmisprediction rate (the number of branch mispredictions per input
element), the CPI (cycles per instruction), and the path length (thenumber of instructions executed per input element).
We begin with the microarchitectural statistics of ouralgorithm when not using the SIMD instructions. When
comparing the statistics for our scalar algorithm and the naive
algorithm, which is equivalent to our algorithm with a block sizeof 1, the branch mispredictions are reduced as intended by using
the larger block sizes. The reduction in the branch mispredictions
directly affected the overall CPI, which was improved when weused the larger block sizes. The improvements in CPI were
especially significant when we enlarged the block size from 1 to 2
and from 2 to 3. By using our scalar algorithm with the block sizeof 3, the branch mispredictions were reduced by more than 75%
compared to the naive algorithm on both platforms, which was
higher than the predicted reduction of 66%.
In contrast, the path lengths increased steadily with theincreasing block sizes. As a result of the reduced CPI and theincreased path length, our best performance without SIMD
instructions was with the block sizes of 3 and 4. When the block
size increased beyond 4, the benefits of reduced branchmispredictions were not significant enough to compensate for theincreased path length. This supports our belief that the best block
size might be larger on processors with larger branchmisprediction overhead. Since most of today’s high performance
processors use pipelined execution units and typically have large
branch misprediction overhead, we expect that our algorithm
would be generally effective for most of the modern processors,
not just the two tested processors.
The branchless algorithm showed the smallest number of the
branch mispredictions by totally replacing the hard-to-predict
conditional branches with arithmetic operations. However, as
shown in Figure 7, the path length of the branchless algorithm
was larger than our algorithm. Due to this long path length, thebranchless algorithm did not outperform our scalar algorithm in

Shorter is better

1.8cycles per instruction (CPI)

1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
Xeon

POWER7+

instructions executed per input element

Path length (Instructions executed per input element)
18.0
16.0
Shorter is better

4.2

Branch mispredictions per input element

14.0
12.0
10.0
8.0
6.0
4.0
2.0
0.0
Xeon
Naive (scalar, block size = 1)
Ours (scalar, block size = 3)
Ours (scalar, block size = 5)
Branchless (scalar)
Ours (STTNI, block size = 8)POWER7+
Ours (scalar, block size = 2)
Ours (scalar, block size = 4)
Ours (SIMD, block size = 4)
V1 SIMD algorithm [9]Figure 7. Branch misprediction rate, CPI, and path length.
spite of its small number of branch mispredictions. Our algorithm
achieved comparable or even better CPI than the branchlessalgorithm even with the larger numbers of branch mispredictions.
This is due to our algorithm’s higher instruction-level parallelism,
since all of the comparisons in the all-pair comparisons can bedone in parallel on the hyperscalar processors.
For our algorithm with the SIMD instructions, we observed
significant improvements in the path lengths. Because the parallelcomparisons of the SIMD instructions make the all-pairs
comparisons of the costly scalar comparisons unnecessary in mostcases, this greatly reduced the number of instructions executed,
even with the large block size of 4. When we use STTNI on Xeon,
our algorithm achieved further reductions in the path lengths. Dueto the shorter path lengths, the SIMD instructions showed huge
boosts in the overall performance.

300

Scalar algorithms on Xeon

Scalar algorithms on POWER7+
1600

STL

higher is fasterOur scalar algorithm (block size = 2 x 4)
branchless algorithm
galloping algorithm

800
600
400
200

1400
1200

STL
Our scalar algorithm (block size = 3 x 3)
Our scalar algorithm (block size = 2 x 4)
branchless algorithm
galloping algorithm

1000
800
600
400
200

0

0
4k/4k

4k/8k

4k/16k

4k/32k

4k/64k4k/128k 4k/256k 4k/512k

4k/1M

4K/4K

4K/8K

input array sizes (number of elements)

SIMD algorithms on Xeon

STL
4000

Our SIMD algorithm (block size = 4 x 4)

higher is faster

Our SIMD algorithm w/ STTNI (block size = 8 x 8)
V1 SIMD algorithm [9]SIMD galloping algorithm [9]

2000
1500
1000

performance (input entries / usec)

performance (input entries / usec)

4K/1M

SIMD algorithms on POWER7+Our SIMD algorithm (block size = 4 x 8)

2500

4K/64K 4K/128K 4K/256K 4K/512K

4500
STL

3000

4K/32K

input array sizes (number of elements)

4000
3500

4K/16K

500

3500
3000

Our SIMD algorithm (block size = 4 x 4)
Our SIMD algorithm (block size = 4 x 8)

higher is faster1000

Our scalar algorithm (block size = 3 x 3)

performance (input entries / usec)

performance (input entries / usec)

1200

higher is faster

1400

V1 SIMD algorithm [9]SIMD galloping algorithm [9]

2500
2000
1500
1000
500
0

0
4k/4k

4k/8k

4k/16k

4k/32k

4k/64k4k/128k 4k/256k 4k/512k

4K/4K

4k/1M

4K/8K

4K/16K

4K/32K

4K/64K 4K/128K 4K/256K 4K/512K

4K/1M

input array sizes (number of elements)

input array sizes (number of elements)

Figure 8. Performance of scalar and SIMD algorithms for intersecting 32-bit integer arrays on Xeon and POWER7+ when the sizes of thetwo input arrays are different. The selectivity was set to 0.
4.3

than STL. The branchless algorithm does not incur the branchmisprediction overhead and hence its performance was notaffected by the size ratio of the two arrays. As shown in many
previous studies, when the ratio of the input sizes exceeds anorder of magnitude, binary-search-based algorithms, such as thegalloping algorithm in the figure, outperform the merge-basedalgorithms, including our algorithm.
For our SIMD algorithms, the block size of 4x8 yielded better
performance than the block size of 4x4 (or the block size of 8x8
with STTNI on Xeon) when the sizes of the two arrays are
significantly different, while the block size of 4x4 gave the best
performance when the two arrays are of the same size (theleftmost point in the figure). When the ratio of the input arraysizes is very large, the V1 SIMD algorithm and the SIMD
galloping algorithm [9] had better performances than our SIMD
algorithm with the block size of 4x4 or 4x8. The V1 algorithm is amerge-based algorithm and is very similar to our algorithm with ablock size of 1x8 implemented with SIMD instructions. Although
this block size gave better performance for two arrays with very
different sizes than 4x4 or 4x8, the binary-search-based galloping
algorithm implemented with SIMD outperformed any mergebased algorithms we tested with such inputs.
Based on these observations, we combined our block-basedalgorithm with the galloping algorithm, so as to improve theperformance for datasets with very different sizes. We select thebest algorithm based on the ratio of the sizes of the two input
arrays. When using SIMD, we use the SIMD galloping algorithm
when the ratio of input sizes is larger than 32. Otherwise we use
our new SIMD algorithm. We use a block size setting of 4x8
when the ratio of input sizes is larger than 2.0, and a block size
setting of 4x4 when the difference is smaller than this threshold.

Performance For Two Arrays of Various Sizes

In this section, we show how the differences in the sizes of thetwo input arrays and the total sizes of the input arrays affect theperformance of each algorithm.
Figure 8 compares the performances of scalar and SIMD
algorithms for 32-bit integer arrays with changing ratios between
the sizes of the two input arrays. When comparing our scalar
algorithm with different block sizes, it worked best with the block
size of Sa = Sb = 3 (we denote this block size as 3x3) when thesizes of two input arrays are the same (the leftmost point in the
figure), while Sa = 2 and Sb = 4 (block size 2x4) worked better
than the block size of 3x3 when the larger of the two input arrays
was at least twice as large as the smaller array, as predicted in
Section 3.2. For two input arrays with very different sizes, thenumbers of branch mispredictions with merge-based algorithms,
STL and ours, became much smaller than for the two arrays of thesame size. When the sizes of the two input arrays are different, thehard-to-predict conditional branches to select which array’s
pointer to advance, e.g. the branches shown in bold in Figure 1,
become relatively easier to predict because the frequency of one
branch direction (taken or not-taken) becomes much higher thanthe other direction on average. This means there were fewer
opportunities to improve the performance with our scalar
algorithm. As a result, the absolute performances became higher
for these algorithms and also the benefits of the reduced branchmispredictions with our algorithm became smaller with the larger
gaps between the sizes of the two arrays. However, even for thelargest differences between the sizes of the two arrays, our scalar
algorithm with the block size of 2x4 achieved higher performance301

on Xeon

1500

Fit in L3 (20 MB)

1000

500

1400

1K/1K

4K/4K

16K/16K 64K/64K256K/256K 1M/1M

STL

1200

Our scalar algorithm (block size = 3 x 3)

1000

0

higher is faster2000

higher is fasterhigh performance at small
sizes by very low branchmisprediction rate

1600

performance (input entries / usec)

performance (input entries / usec)

2500

on POWER7+STL
Our scalar algorithm (block size = 3 x 3)
Our SIMD algorithm (block size = 4 x 4)
Our SIMD algorithm w/ STTNI (block size = 8 x 8)
V1 SIMD algorithm [9]Our SIMD algorithm (block size = 4 x 4)

800

V1 SIMD algorithm [9]

600
400
200

Fit in L3 (10 MB)

0

4M/4M 16M/16M 32M/32M

1K/1K

4K/4K

16K/16K 64K/64K256K/256K 1M/1M

input array sizes (number of elements)

4M/4M 16M/16M 32M/32M

input array sizes (number of elements)

Figure 9. Performance for intersecting 32-bit integer arrays of various sizes. The selectivity was set to 0.

800
700

higher is faster

500
400
300
200

0

STL
Our SIMD algorithm (block size = 4 x 4)
Our SIMD algorithm w/ STTNIV1 SIMD algorithm [9]
0%

10%

20%

30%

40%

50%

Our scalar algorithm (block size = 3 x 3)
Our adaptive SIMD algorithm
Our adaptive SIMD algorithm w/ STTNI( 1497 for Ours adaptive SIMD algorithm
1501 for Our SIMD algorithm at 0%)
our adaptive algorithm
switches to scalar algorithm

600
500
400
300
200
100

STL

Our scalar algorithm (block size = 3 x 3)

Our SIMD algorithm (block size = 4 x 4)

Our adaptive SIMD algorithm

V1 SIMD algorithm [9]

0
60%

70%

80%

90%

0%

100%

10%

20%

selectivity (output size / input size)

higher is faster

500
400
300
200
100

STL

Our scalar algorithm (block size = 2 x 4)

Our SIMD algorithm (block size = 4 x 8)

Our adaptive SIMD algorithm

V1 SIMD algorithm [9]performance (input entries / usec)

performance (input entries / usec)

1200

our adaptive algorithm
switches to scalar algorithm

600

40%

50%

60%

70%

80%

90%

100%

256K/1M input size on POWER7+( 1145 for Our adaptive algorithm fallback1293 for Our SIMD algorithm at 0%)

700

30%

selectivity (output size / input size)

256K/1M input size on Xeon
800

switch to
STL0( 1744 for Ours with adaptive fallback
1707 for Our SIMD algorithm at 0%)
higher is faster100

performance (input entries / usec)

performance (input entries / usec)

600

256K/256K input size on POWER7+( 989 for Our adaptive SIMD algorithm, 1053 for Our SIMD algorithm
1664 for Our adaptive SIMD algorithm w/ STTNI,
1612 for Our SIMD algorithm w/ STTNI at 0%)
our adaptive algorithm
switch toswitches to scalar algorithm
STL

higher is faster

256K/256K input size on Xeon

1000

our adaptive algorithm
switches to scalar algorithm

800
600
400
200

STL

Our scalar algorithm (block size = 2 x 4)

Our SIMD algorithm (block size = 4 x 8)

Our adaptive SIMD algorithm

V1 SIMD algorithm [9]
0

0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

0%

selectivity (output size / input size)

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

selectivity (output size / input size)

Figure 10. Performance of each algorithm for random 32-bit integers with various selectivity on Xeon and POWER7+.
For the scalar algorithm, we used the (non-SIMD) galloping
algorithm if the ratio is larger than 32. Otherwise, we use ourblock-based algorithm. The block size setting is 2x4 if the ratio islarger than 2.0 or otherwise the block size setting is 3x3.
Figure 9 shows how the total size of the two input arrays
affects the performance. We used 32-bit integer arrays and theselectivity was zero. On both platforms, we observed small
performance degradations when the total size exceeded the lastlevel (L3) cache of the processor because of the stall cycles towait for data to be loaded from the main memory. The effects ofthe stall cycles due to the cache misses were not significant,because the memory accesses in the set intersection are almost

sequential and this means the hardware prefetcher of theprocessors worked well to hide that latency by automaticallyreading the next data into the cache. The performance advantagesof our scalar and SIMD algorithms over the other algorithms, theSTL and V1 SIMD algorithms, were unchanged, even with thelargest datasets we tested. On Xeon, the performance of four outof five tested algorithms was significantly improved when the
input size was very small (left side of the figure). This was caused
by very low branch misprediction overhead rather than the
reduced cache miss stall cycles. The Xeon seems to employ abranch prediction mechanism that is very effective only when the
input size is very small. POWER7+ did not exhibit this behavior.

302

4.4 Adaptive Fall Back Based on Selectivity to Avoid
Performance Degradations

select algorithm based on the differencein the sizes of the two input arrays
Start of SIMD
algorithm

Figure 10 shows how the selectivity affected the performance of
our algorithm using 32-bit integers based on a random number
generator with various selectivity values. Our algorithm worked
best when the selectivity was small, which is true for many realworld applications. For example, Ding and König [2] reported
that the selectivity was less than 10% for 94% of the queries andless than 1% for 76% of the queries in the 10,000 most frequentqueries in a shopping portal site. For this frequent situation, ouralgorithm worked well, especially with the SIMD instructions.
However, the performance of our algorithm was worse than
STL when the selectivity was high. To avoid these performancedegradations, we added an adaptive fallback mechanism to ouralgorithm. We start execution with our SIMD algorithm, but witha periodic runtime check of the selectivity that may trigger the
fallback mechanism. We calculate the selectivity after each 1,024
output elements by checking the numbers of input elements
processed to generate these output elements. When the numbersof input elements is larger than the threshold in at least one array,we fall back to another algorithm. This insures the overheadcaused by the runtime check is not significant when there are few
output elements. An adaptive fallback using a runtime check is astandard heuristic technique to avoid worst case performance inmany algorithms. For example, introsort [19] used in the STL’s
std::sort library method uses quicksort with adaptive fallback toheapsort to avoid the O(N2) worst-case performance of quicksort.
From the results shown in Figure 10, for two input arrays withcomparable sizes, we start execution with our SIMD algorithm
using the block size setting of 4x4 (or 8x8 if we use STTNI on
Xeon). We switch to the STL when the selectivity is higher than
65%. When the selectivity is higher than 15% but lower than 65%,
we use our scalar algorithm with the block size setting of 3x3. We
also execute a periodic check in our scalar algorithm that may fallback to the STL algorithm. If one of the input arrays is more than
twice as large as the other, we start execution with our SIMD
algorithm using the block size setting of 4x8 and fall back to ourscalar algorithm with the block size setting of 2x4 if the selectivity
becomes higher than 35%. We do not switch to STL because out
scalar algorithm consistently outperformed STL in Figure 10. We
summarize how we select the algorithm and the block size basedon the size of two input arrays and the selectivity with andwithout using SIMD instructions in Figure 11. We call theseoverall algorithms the adaptive SIMD algorithm and the adaptive
scalar algorithm. Figure 10 shows that our fallback mechanisms
selected the appropriate algorithm for each selectivity.
4.5

use STTNI
if enabled
selectivity
> 65%

< 1:2

1:2~1:32

SIMD algorithm(block size 4x4)
selectivity> 15%

SIMD algorithm(block size 4x8)

> 1:32SIMD galloping[9]

selectivity> 35%

Start of scalar
algorithm< 1:2

1:2~1:32

scalar algorithm(block size 3x3)
selectivity> 65%

STL (naive
merge-based)scalar algorithm(block size 2x4)

> 1:32

galloping [10]a dashed line shows
adaptive fallback witha runtime check of selectivity

our adaptive scalar algorithm
our adaptive SIMD algorithm

Figure 11. Overall scheme of our adaptive algorithm.
words and generated 2-word, 3-word, 6-word, and 8-word queries
by randomly selecting the keywords from the 16 prepared words.
The size of the list for each keyword ranged from 10,733 elements
to 528,974 elements. For each class, we generated 100 queries
and measured the total execution time of these queries. Forintersecting multiple words, we repeatedly picked the two smallest
sets and did set intersection for the two arrays, a technique that isfrequently described in the literature.
Figure 12 shows the relative performance of each algorithm
over the baseline (STL + galloping). On both platforms, our
SIMD algorithm more than doubled the baseline performance.
The V1 SIMD + SIMD galloping algorithm also accelerated theoperation by exploiting SIMD instructions, but its gain was about
60% on both platforms and hence our SIMD algorithm
outperformed V1 SIMD + SIMD galloping by from 24% (3-word
queries on Xeon) to 44% (8-word queries on Xeon). Although V1
SIMD + SIMD galloping and our SIMD algorithm use the same
SIMD galloping algorithm when intersecting two arrays with very
different sizes, our algorithm achieved higher performance forarrays with similar sizes and this mattered for the overall
performance. On Xeon, our SIMD algorithm can achieve evenhigher performance with STTNI. Schelegel’s algorithm also
accelerated the set intersection using the STTNI instruction, while
the algorithm performed much worse than STL for the artificial
dataset generated by the random number generator, as shown in
Figure 6. This is because the value domain for the Wikipedia
dataset was smaller than the artificial datasets and hence moreelements shared the same values in their upper 16 bits. This isimportant for Schelegel’s algorithm because they use STTNI tofind matching pairs in the lower 16 bits within the elements thatshare the same value in the upper 16 bits. However, theperformances of Schelegel’s algorithm were not as good as our
SIMD algorithm or the V1 SIMD algorithm.
Our scalar algorithm improved the performance by about 50%
over the baseline in spite of not using the SIMD instructions. The
performance of STL alone was significantly lower than the otheralgorithms because STL, or merge-based-algorithms in general,
performed poorly when the sizes of two arrays are quite different
and this configuration is known to be important for intersecting
multiple sets.

Performance of our algorithm with realistic datasets

Finally, we evaluated the performance of our algorithms forrealistic datasets generated from a Wikipedia database dump toemulate the set intersection operation in a query serving system.
We compare the performance of set intersections of multiple
arrays to emulate multi-word queries. Here we compare our SIMD
and scalar algorithm against a combination of existing SIMD
algorithms, the V1 SIMD algorithm with SIMD galloping [9]. We
switched between these two algorithms based on the difference inthe two input arrays and we used 1:50 as the selection thresholdbased on their results. We also compared the performance of our
SIMD algorithm with the STTNI instruction against Schlegel’s
algorithm [8], which also exploits the STTNI, combined with
SIMD galloping on Xeon. As a baseline, we also measured a
combination of STL (as a representative merge-based algorithm)
and a scalar galloping algorithm (as a binary-search-basedalgorithm). We prepared a list of document IDs for 16 search303

using STTNI unique to Xeon

2.0
1.5
1.0
0.5
0.0
2-word query

3-word query

Our adaptive SIMD algorithm
V1 SIMD + SIMD galloping [9]Our adaptive SIMD algorithm w/ STTNISTL only

6-word query

8-word query

Our adaptive scalar algorithm
STL + galloping (baseline)
Schlegel's algorithm [8] + SIMD galloping

2.5

on POWER7+

higher is faster

2.5

relative performance over (STL + galloping)on Xeonhigher is faster

relative performance over (STL + galloping)

3.0

2.0
1.5
1.0
0.5
0.0
2-word query

3-word query

Our adaptive SIMD algorithm
V1 SIMD + SIMD galloping [9]STL only

6-word query

8-word query

Our adaptive scalar algorithm
STL + galloping (baseline)

Figure 12. Performance of set intersection algorithms using the datasets generated from Wikipedia database.

5. SUMMARY[8] B. Schlegel, T. Willhalm, and W. Lehner. Fast Sorted-Set
Intersection using SIMD Instructions. In Proceedings of theSecond International Workshop on Accelerating DataManagement Systems Using Modern Processor and Storage
Architectures, 2011.
[9] D. Lemire, L. Boytsov and N. Kurz. SIMD Compression andthe Intersection of Sorted Integers. arXiv:1401.6399, 2014[10] J. L. Bentley and A. C. Yao. An almost optimal algorithm forunbounded searching. Information processing letters, 5(3),pp. 82-87, 1976.
[11] J. R. Allen, K. Kennedy, C. Porterfield, and J. Warren.
Conversion of control dependence to data dependence. In
Proceedings of the Symposium on Principles ofProgramming Languages, 1983.
[12] J. Zhou and K. A. Ross. Implementing database operations
using SIMD instructions. In Proceedings of the ACM
SIGMOD International Conference on Management of Data,
pp. 145–156, 2002.
[13] P. Sanders and S. Winkel. Super Scalar Sample Sort. In
Proceedings of the European Symposium on Algorithms,
Volume 3221 of LNCS, pp. 784–796, 2004.
[14] H. Inoue, T. Moriyama, H. Komatsu, and T. Nakatani. AASort: A New Parallel Sorting Algorithm for Multi-CoreSIMD Processors, In Proceedings of the International
Conference on Parallel Architectures and Compilation
Techniques, pp. 189–198, 2007.
[15] J. Chhugani, A. D. Nguyen, V. W. Lee, W. Macy, M. Hagog,
Y. Chen, A. Baransi, S. Kumar, and P. Dubey. Efficient
implementation of sorting on multi-core SIMD CPU
architecture. In Proceedings. VLDB Endow., 1(2), pp. 1313–
1324, 2008.
[16] D. Tsirogiannis, S. Guha, and N. Koudas. Improving theperformance of list intersection. In Proceedings of VLDB
Endow., 2(1), pp. 838–849, 2009.
[17] S. Tatikonda, F. Junqueira, B. B. Cambazoglu, and V.
Plachouras. On efficient posting list intersection withmulticore processors. In Proceedings of the 32nd
International ACM SIGIR Conference on Research andDevelopment in Information Retrieval, 2009.
[18] N. Ao, F. Zhang, D. Wu, D. S. Stones, G. Wang, X. Liu, J.
Liu, and S. Lin. Efficient parallel lists intersection and index
compression algorithms using graphics processing units. In
Proceedings of VLDB Endow., 4(8), pp. 470–481, 2011.
[19] D. R. Musser. Introspective Sorting and Selection
Algorithms. Software Practice and Experience, 27(8), pp.
983–993, 1997.

This paper described our new highly efficient algorithm for set
intersections on sorted arrays on modern processors. Our
approach drastically reduces the number of branch mispredictions
and efficiently exploits the SIMD instructions. Our algorithm is
not only efficient but also portable, easy to implement, andrequires no preprocessing. Our results show that our simple andportable scalar algorithm improved the performance of theimportant set intersection operation by reducing the branchoverhead. The use of the SIMD instructions in our algorithm gave
additional performance improvements by reducing the path length
significantly for many datasets. We believe our algorithm will bequite effective to improve the performance of the set intersection
operations for many workloads.

6. ACKNOWLEDGEMENTS
We are grateful to the anonymous reviewers for their valuablecomments and suggestions. We thank Toshio Nakatani, Tamiya
Onodera, and Takanori Ueda for their useful feedback on earlier
drafts on this work.

7. REFERENCES[1] L. A. Barroso, J. Dean, and U. Hölzle. Web Search for a
Planet: The Google Cluster Architecture. IEEE Micro 23(2),pp. 22–28. 2003.
[2] B. Ding and A. C. König. Fast set intersection in memory. In
Proceedings of VLDB Endow., 4, pp. 255–266, 2011.
[3] R. Baeza-Yates and A. Salinger. Experimental Analysis of a
Fast Intersection Algorithm for Sorted Sequences. In
Proceedings of the International Conference on String
Processing and Information Retrieval, pp. 13–24, 2005.
[4] R. Baeza-Yates. A Fast Set Intersection Algorithm for Sorted
Sequences. In Proceedings of the Annual Symposium on
Combinatorial Pattern Matching, pp. 400–408, 2004.
[5] E. D. Demaine, A. López-Ortiz, and J. I. Munro. Adaptive
set intersections, unions, and differences. In Proceedings of
the Annual Symposium on Discrete Algorithms, pp. 743–752,
2000.
[6] P. Bille, A Pagh, and R. Pagh. Fast evaluation of unionintersection expressions. In Proceedings of the International
Conference on Algorithms and Computation, pp. 739–750,
2007.
[7] P. Sanders, F. Transier. Intersection in Integer Inverted
Indices. In Proceedings of the Workshop on Algorithm
Engineering and Experiments, pp. 71–83, 2007.

304