QuadRank: Engineering a High Throughput Rank

A primary application of Rank queries is in the FM-index [13], a succinct data structure that can efficiently locate all occurrences of a pattern in a text and is used in tools such as BWA-MEM [32] and Bowtie [30, 29]. Whereas most of the literature on rank structures assumes a binary alphabet ($\sigma =2$), in this case the DNA alphabet has size $\sigma =4$. Indeed, BWA-MEM implements its own rank structure over a 2-bit alphabet²²2https://github.com/lh3/bwa/blob/master/bwt.c, and this paper started as an attempt to speed this up.

For alphabets of arbitrary size, wavelet trees [21] or the wavelet matrix [8] can be used for succinct rank queries. They both need ${\lg }_2\sigma$ queries to a binary rank structure. Recently, quad wavelet trees [6] have been introduced, following earlier theoretical [14] and practical [3] results on multi-ary wavelet trees. Quad wavelet trees use rank over a quad vector as a building block, and thus need only ${\log }_4\sigma$ rank queries, leading to 2$\times$ to 3$\times$ speedups over binary wavelet trees.

In the past, increasing CPU frequencies led to faster code, but nowadays, improvements are mostly in increasing parallelism. Furthermore, total compute of a CPU increases faster than the available memory bandwidth, resulting in the need for communication-avoiding algorithms that minimize their memory bandwidth [11]. Indeed, in bioinformatics applications, one often has many independent queries (DNA sequences) that need to be processed (searched in an FM-index) as fast as possible. In particular, this allows using all cores/threads of the CPU as well as processing queries in batches inside each thread, to hide the memory latency.

Current benchmarks usually measure the throughput of answering rank queries in a for loop on a single thread, but this does not take into account the possibility for batching, nor does it capture the effects of running many threads in parallel. As we will see, in a high-throughput setting, many existing methods become bottlenecked by the total memory bandwidth of the CPU. We specifically design our data structures to use the memory bandwidth maximally efficient.

We develop two data structures, BiRank and QuadRank, that support high-throughput rank queries over texts over alphabets of size 2 and 4.

Both of them integrate a number of existing techniques (see next section), and are not designed to support select queries, since these are not needed for the motivating FM-index application, thus allowing for more optimizations. Specifically, BiRank has 3.28% overhead and integrates (1) inlining of counts into the bitvector [31], which reduces cache misses, (2) pairing with mask-lookup [18], halving the number of popcounts, and (3) an additional third layer [47] that is modified to be only half its usual size.

QuadRank extends the ideas of BiRank, but has roughly 4$\times$ larger space overhead (14.4%) since it stores metadata for each symbol. It combines the cache locality of the implementation in BWA-MEM [32] with the low overhead of quad vectors [5] and a transposed bit layout for faster queries [1, 18]. QuadRank is optimized for returning ranks for all 4 symbols at once by using AVX2 instructions, which is useful for approximate pattern matching in an FM-index.

Both data structures need only a single cache line from RAM to answer queries as long as the input is less than roughly 16 GiB. The main novelty is in combining all these existing ideas and applying them to the size-4 alphabet, while also adding support for prefetching of cache lines to enable much more efficient batch-processing of queries. As a side-effect, we also added prefetching to some other libraries, yielding up to 2$\times$ speedups.

For both data structures, we implement a number of variants that have different space-time trade-offs. When used in a for loop, BiRank is up to 1.5$\times$ faster than the next-fastest rust implementation of equal size, with the speedup being larger when using many threads. Prefetching memory improves the throughput of many libraries by around 1.5$\times$, and improves BiRank by 2$\times$. In this setting, all methods are bottlenecked by the memory throughput, and BiRank is 2$\times$ faster than all others because it only needs to read 1 instead of 2 cache lines from RAM. Similarly, QuadRank is at least 1.5$\times$ faster than the next-fastest Rust library, QWT [6], and 2$\times$ faster after adding prefetch instructions, again being bottlenecked by the RAM throughput.

Inspired by genedex [12], we develop QuadFm, a toy-implementation of a count-only FM-index that uses batching, prefetching, and multithreading. At 14.4% overhead (2.29 bits/bp), our implementation is over 1.5$\times$ faster using QuadRank compared to using QWT’s quad vector, and at 100% space overhead, QuadFm is 4$\times$ faster than genedex, a state-of-the-art FM-index implementation.

For later reference, we summarize our terminology. The raw data is split into superblocks (the second level, L2) that are further split into blocks (L1). Block is used for both the raw bits themselves, as well as the cache line containing them. For each superblock, an L2 offset is stored, representing the number of 1-bits before the superblock. For each block, an L1 delta is stored, typically representing the number of 1-bits preceding it inside the superblock. Conceptually, these levels form a summary tree, named as such by Kurpicz et al. [27]. We follow their notation, and number the levels of the tree bottom-up, breaking tradition with e.g. the presentation in [47] and [26]. The root has all superblocks (L2) as children, each superblock has its contained blocks (L1) as children, and each block has the contained bits as children.

The overhead of a data structure is the increase in space consumption relative to the size of the input data. We use bp (base pair) as the unit for 2-bit encoded DNA characters, and occasionally use the Rust syntax u64 for a 64-bit variable and u64x4 for a 256-bit SIMD register containing 4 64-bit words. A symbol is an element of the alphabet $\mathrm{\Sigma }$, whereas a character is an element of a string.

As a baseline, Jacobson [22] stores the bitvector, and additionally two levels of blocks alongside this. Blocks consist of $\lfloor \log (n)/2\rfloor$ bits, and $\lfloor \log n\rfloor$ blocks form a superblock. Level L2 of the tree then contains a $\lceil \log n\rceil$ bit offset for each superblock, counting the number of set bits preceding it. Level L1 stores for each block a $\lceil \log \log n\rceil$ bit delta counting the number of 1-bits preceding it inside its superblock. The number of 1-bits in a (prefix of) a block is obtained via a lookup in a precomputed table of size $2^{(\log n)/2}=\sqrt{n}$.

González et al. [17] observe that the classic method above has 66.85% overhead in practice for $n=2^{30}$. They replace the large lookup table by a smaller table of per-byte popcounts. (Meanwhile, CPUs natively support 64-bit popcount instructions.) They use 256-bit superblocks with a 32-bit offset, containing 8 32-bit blocks, each with their own 8-bit delta. Alternatively, they introduce a single-level tree storing a 32-bit L2 offset after every e.g. $4\cdot 32$ bits and omitting L1. This requires popcounting more words, but has the benefit of improved cache locality compared to a two-level tree.

Rank9 [43] has 25% overhead and interleaves the L2 and L1 levels of the classic tree. Each block is 64 bits, and 8 blocks form a 512-bit superblock, exactly matching a cache line. For each superblock, the interleaved tree stores a 64-bit integer with the offset of the superblock, and 7 9-bit deltas (for all but the first block) in an additional 64-bit word. This needs two cache misses per query (for the L2 array and bits), and is very fast in practice. Specifically it only needs to popcount a single 64-bit word, which is done using broadword programming (also known as SWAR, SIMD Within A Register).

Poppy [47] is optimized for space and has only 3.125% overhead. First, it makes the observation that performance is largely determined by the number of cache misses. Thus, it uses larger blocks of 512 bits. It then re-uses Rank9’s interleaved index with two modifications. Each superblocks contains 4 blocks, and for each superblock it stores a 32-bit offset (L2) followed by 3 10-bit popcounts of the first 3 blocks. Queries then require a prefix-sum over these counts. To handle 64-bit outputs, it stores an additional layer (L3) of the tree, with a full 64 bit offset after every $2^{32}$ input bits.

BWA-MEM [32] implements a 100% overhead rank data structure on $\sigma =4$ DNA. It interleaves L2 offsets with the data, and requires only a single cache miss per query. In each cache line, it stores 4 64-bit offsets (one for each DNA character), followed by 256 bits encoding 128 bp.

The succinct data structure library (SDSL) [15] implements Rank9 and introduces rank_support_v5, which has 6.25% overhead. It uses superblocks of 2048 bits. For each, it interleaves a 64-bit offset (L2) and 5 11-bit deltas (packed into 64 bits) to all but the first of 6 blocks covering $6\cdot 64$ bit. rank_support_il interleaves 64-bit offsets with 512-bit blocks.

EPR-dictionaries [38] work for arbitrary alphabet. For $\sigma =4$, they use 64-bit (32 bp) blocks and have 42% overhead, and interleave an independent 2-level rank structure for each character. Compared to earlier work, space is saved by storing a packed representation of the text instead of $\sigma$ (1-hot) encoded bitvectors that each indicate which text positions contain each symbol $c\in \mathrm{\Sigma }$.

Pibiri and Kanda [37] diverge from the classic approach and introduce a rank and select structure based on highly tuned B-trees that have 3.6% overhead. Each rank query traverses roughly ${\log }_{16}n$ levels of the tree, with the middle levels packing 16 32-bit values in a cache line. Due to efficient caching of the top levels of the tree, performance is similar to poppy, although not as fast as rank9.

The AWFM-index and its Rust implementation AWRY [1] builds an FM-index on a size $\sigma =6$ alphabet of 4 DNA characters as well as a sentinel and ambiguity symbol. It uses blocks of 256 3-bit characters, preceded by 5 64-bit offsets that are padded to 512 bits. Each block is encoded using a similar strided or transposed layout: instead of concatenating the 3 bits of each character, it stores 3 256-bit vectors containing bit 0, bit 1, and bit 2 of each character. This allows for more efficient popcounting. The FM-index processes queries in batches of size 4, and prefetches memory needed for the next rank operation as soon as possible.

PastaFlat [26, 25] has the same 3.125% space overhead as Poppy, but improves query time by 8% by avoiding Poppy’s need to take a prefix sum over L1 counts. Pasta doubles the metadata for each superblock to 128 bits, covering 8 512-bit blocks of 4096 bits in total. It stores a 44-bit offset (L2) followed by 7 12-bit deltas (L1) from the start of the superblock to each block. A second structure, PastaWide (3.198% overhead) uses 16-bit values for L1, which allows faster select queries using SIMD instructions. Here, each superblock covers 128 blocks and stores a 64-bit L2 value, this time not interleaved with the L1 values, and the L3 level is dropped.

Quad wavelet trees internally use quad vectors [5, 6], which have a layout very similar to PastaFlat. Super blocks cover eight 512 bp blocks and stores 128 bits of data for each of the 4 symbols. This takes 4$\times$ more space per character, but since the text doubles in space as well, the overhead only doubles to 6.25%. Alternatively, 256 bp (512 bit) blocks can be used to reduce the number of cache misses, using 12.5% overhead.

SPIDER [31] has only 3.3% overhead and reduces the number of cache misses from 2 to (nearly) 1 by interleaving L1 with the bitvector itself (like BWA-MEM), instead of interleaving L1 with L2: each cache line stores a 16-bit L1 delta, and 496 input bits. L2 superblocks store a 64-bit offset for each 128 blocks, taking only 0.1% extra space and thus likely fitting in a cache.

Pairing (pfBV) [18] is an idea that halves the memory overhead again, to 1.6%. Compared to PastaWide, instead of storing 16-bit (L1) deltas to the start of each 512-bit block, here we store 16-bit deltas to the middle of each pair of 512-bit blocks. Then, the second block can add a prefix-popcount to this as usual, while the first block can subtract a suffix-popcount instead. Similarly, the 64-bit L2 offset is to the middle of a twice-as-large superblock. This is similar to the alternate counters idea for the FM-index [7], where, for alphabet size 4, each block stores half the offsets. A small complication with this design is that conditionally popcounting a prefix or suffix of bits is slightly slower. Instead, Gottlieb and Reinert [18] introduce a lookup table that stores a precomputed mask for each position. Lastly, for $\sigma =4$, this paper uses the transposed layout of AWFM, but calls it scattered instead.

We aim to minimize the number of cache misses on large (many GB) inputs, to enable efficient usage in high-throughput settings where the memory bandwidth is the bottleneck. A single cache miss is inevitable, and so we must avoid any further cache misses. This means that any additional data should fit in L3 cache, which is not the case for non-interleaved layouts. For example, pairing has 1.6% overhead, which would only support 1 GiB of input with a 16 MiB L3 cache.

Thus, like SPIDER [31], BiRank inlines a 16-bit L1 delta $b_j$ into each block/cache line $j$ (lines 7-8 of Listing 1), so that each of $\lceil n/B\rceil$ blocks covers $B:=512-16=496$ bits. The $\lceil n/S\rceil$ superblocks cover $S:=128\cdot B$ bits each, and a second much smaller array stores a 32-bit L2 offset $s_i$ for each superblock $i$.

Poppy [47] uses an additional 64-bit third level L3, so that 32-bit L2 values are sufficient. Even though this L3 layer is already very small, we remove it completely. Rather than directly encoding $\mathrm{𝗋𝖺𝗇𝗄}(i\cdot S)$, the rank at the start of a superblock, we store

in the 32-bit L2 value. The remainder $\mathrm{𝗋𝖺𝗇𝗄}(i\cdot S)mod{2}^{11}$ will be added to the 16-bit $b_j$ delta for each block in the superblock. This configuration supports inputs up to $2^{43}$ bits, or 1 TiB, since implies .

Each superblock must contain at most $2^{16}$ bits, so that the 16-bit $b_j$ can represent their deltas. Thus, we could fit $\lfloor 2^{16}/B\rfloor =132$ blocks inside each superblock, but we round this down for computational efficiency: $S:=128\cdot B$.

We also implemented a variant of the pairing of superblocks technique of [18], which doubles the superblock size $S$ and halves the cache usage, see Appendix B. We decided not to use it though: while the reduced cache usage could be beneficial, in practice, the gains are inconsistent and small at best.

The overhead of the block deltas is $16/B=3.226\%$, whereas the superblocks have an overhead of $32/S=0.05\%$. Thus, the total overhead is 3.28%, and the superblock array fits in a 16 MiB L3 cache for inputs up to 32 GiB.

To reduce the amount of work needed for popcounting, we apply a variant of the pairing technique [18]: the 16-bit L1 value is not the delta from the start of the superblock to the start of the current block ($\mathrm{𝗋𝖺𝗇𝗄}(j\cdot B)$), but instead to the middle of the current block ($\mathrm{𝗋𝖺𝗇𝗄}(j\cdot B+240)$, after taking into account a 16-bit padding, line 9):

By construction, these values are indeed bounded by (line 8).

A query for position $q$ first determines the superblock $i_q=\lfloor q/S\rfloor$ and block $j_q=\lfloor q/B\rfloor$. Then, we compute the rank of the middle of the block as $\mathrm{𝗋𝖺𝗇𝗄}(j\cdot B+240)=2^{11}\cdot s_{i_q}+b_{j_q}$. We then make a case distinction on whether $q$ lies left or right of the middle of its block to determine the final value

In code, we popcount up to 256 bits (lines 15 and 18): either a suffix of the first half, or a prefix of the second half. The conditional negation (line 20) is optimized to a branchless cmov instruction.

Instead of a for loop over the 64-bit words in the block and bit-shifting, we prefer a branchless technique that always covers the full 256 bits. Uncounted bits are masked out (line 18) via a 256-bit mask that is precomputed (line 11) for each , again following [18]. These are simply stored as a 16 KiB array [u256; 512], which fits in a typical 32 KiB L1 cache.

In order to facilitate efficient batching algorithms (see Appendix D), we provide a prefetch(q) function that starts loading the two cache lines containing $s_{\lfloor q/S\rfloor }$ and $b_{\lfloor q/B\rfloor }$ that are needed for $\mathrm{𝗋𝖺𝗇𝗄}(q)$. For simplicity and reliability, we prefetch into all levels of the cache hierarchy.

The parallel construction algorithm used for both BiRank and QuadRank builds on the rayon crate. First, we count the number of 1-bits in each superblock, which is trivial to parallelize. Then, we take a prefix sum over these counts in a (non-parallel) linear pass, so that we know the number of 1-bits preceding each superblock. This then allows us to fully construct all superblocks (and their contained blocks) in parallel.

Compared to the layout for binary input, the main difference is that we now store the input data in transposed (or strided) layout [1, 18] (as opposed to packed). Ignoring the inlined $b_{j,c}$ deltas for the moment, the 256 characters in a block are split into 4 groups of 64. Each group of 64 characters is encoded as two 64-bit values, one consisting of the negation of all low bits, and one of the negation all high bits. The 4 16-bit deltas replace the 64 bits corresponding to the 32 first characters (Listing 1, line 8), as shown in the bottom row of Figure 1. The positions matching a symbol are then found by and-ing the two values together (line 18), after possibly negating one or both (lines 13-17). This layout makes more efficient use of popcount instructions, since each counted word now contains up to 64 1 bits, compared to 32 with the packed layout.

Computing the rank for a single character is similar to before (1): we retrieve the superblock offset $s_{i_q}$, multiply it by $2^{13}$, and then add the $b_{j_q,c}$ for the current block and character. Lastly, we add or remove the count for up to 128 characters in either the first or second half of the cache line, processed in two chunks of 64 characters.

To return the rank of all 4 symbols, we essentially do the above method 4 times in parallel in u64x4 256-bit AVX2 SIMD registers. In particular, we use a single SIMD lane for each symbol (Listing 1). To popcount the number of 1-bits in each lane, we use Mula’s algorithm [34, 35]. Essentially, this splits each byte into two 4-bit nibbles and for each does a _mm256_shuffle_epi8 instruction to do a 16-way lookup returning the precomputed number of 1 bits in each nibble. It then adds these two values, resulting in per-byte popcounts, and finally uses the _mm256_sad_epu8 instruction to take a horizontal sum of the 8 bytes in each 64-bit lane. We convert the counts to u32x4 and then conditionally negate them using _mm_sign_epi32(counts, u32x4::splat(pos-96)), which multiplies each lane by the sign of $(qmodB)-96$ (i.e., -1, 0, or 1).

Future work remains in generalizing and optimizing the library for other platforms than AVX2. As the current code was optimized for Skylake (2015), it is likely that more modern platforms (Golden Cove, 2021 or zen 5, 2024) admit different trade-offs. For AVX512, there are dedicated popcount instructions that could be used, while ARM NEON only supports 128-bit instructions and will need further work.

Additional features could be in-place parallel construction and a $\mathrm{𝗉𝗋𝖾𝖿𝗂𝗑}\mathrm{\_}\mathrm{𝗋𝖺𝗇𝗄}(q,c)$ operation that counts the number of occurrences of characters at most $c$ to support the bidirectional FM-index. Lastly, our FM-index could be extended with support for $\mathrm{𝗅𝗈𝖼𝖺𝗍𝖾}$ queries.