A bloated FM-index reducing the number of cache misses during the search

Document text

A bloated FM-index reducing the number of
cache misses during the search
Szymon Grabowskiyand Aleksander Cis lakz
yLodz University of Technology, Institute of Applied Computer Science,
Al. Politechniki 11, 90{924  L od z, Poland, [email protected]
zWarsaw University of Technology, Faculty of Mathematics and Information Science,
ul. Koszykowa 75, 00{662 Warsaw, Poland, [email protected]
Abstract. The FM-index is a well-known compressed full-text index,
based on the Burrows{Wheeler transform (BWT). During a pattern
search, the BWT sequence is accessed at \random" locations, which is
cache-unfriendly. In this paper, we are interested in speeding up the FM-
index by working on q-grams rather than individual characters, at the
cost of using more space. The rst presented variant is related to an
inverted index on q-grams, yet the occurrence lists in our solution are in
the sorted sux order rather than text order in a traditional inverted
index. This variant obtains O(m=jCLj+ lognlogm) cache misses in the
worst case, where nandmare the text and pattern lengths, respectively,
andjCLjis the CPU cache line size, in symbols (typically 64 in mod-
ern hardware). This index is often several times faster than the fastest
known FM-indexes (especially for long patterns), yet the space require-
ments are enormous, O(nlog2n) bits in theory and about 80 n{95nbytes
in practice. For this reason, we dub our approach FM-bloated . The second
presented variant requires O(nlogn) bits of space.
1 Introduction
The FM-index [8, 9] invented by Ferragina and Manzini is a well-known com-
pressed data structure based on the Burrows{Wheeler transform (BWT) which
can answer full-text queries. Given a range of text suxes starting a sequence S,
it allows to nd quickly all text suxes starting with cS, for any alphabet symbol
c, using a mechanism called LF-mapping. An inherent property of LF-mapping
is however non-local access to the BWT sequence, which typically results in

(m) CPU cache misses during the search for a pattern of length m, even for a
small alphabet. (Handling a large alphabet of size gives an extra time factor,
e.g., login popular wavelet tree based implementations; see [1,14,15] for more
details and known tradeos.)
The problem of cache misses during the FM-index backward search was iden-
tied as the main performance limiter by Chac on et al. [2], who proposed to
perform the LF-mapping with several symbols at a time (in practice, at most 3,
for the 4-symbol DNA alphabet for which the scheme was only described). This
solution allowed, for example, to improve the search speed by a factor of 1.5 for
the price of doubling the index space.arXiv:1512.01996v1  [cs.DS]  7 Dec 2015
In this paper we address the problem of cache misses during a pattern search
(count query) in a way related to the Chac on et al. solution. We also work on
q-grams, yet the algorithmic details are dierent. We propose two variants, one
withO(m=jCLj+lognlogm) misses in the worst case, but also O(nlog2n) bits of
space, which makes it hardly practical, the other with O(m=jCLj+(m=q0) logn)
misses on average and O(nlogn) bits of space, where q01 is some parameter
andjCLjis the cache line size (nowadays, typically 64 bytes). While in practice
characters over an alphabet of size are \embedded" in chunks of blog
bits (it often holds that b= 8, i.e. the chunks are bytes), it is convenient to
assume that b=(log), and we use this assumption henceforth. Note that the
well-known String B-tree data structure [7] achieves O((m=B ) logBn) I/Os (plus
(occ=B ) I/Os for reporting the matches), which is optimal in the I/O comparison
model. The parameter B, the block size (expressed in the number of items stored
in it), corresponds to our cache line size.
We use a standard notation throughout the paper, namely S[i:::j ] for a
sequenceSand a pair of indexes iandjdenotes a concatenation of symbols
S[i];S[i+ 1];:::;S [j]. In particular, S[i:::j ] is an empty string when i>j . The
index is built for the text T[1:::n], allowing to search for the pattern P[1:::m ].
All logarithms are in base 2.
A preliminary version of the rst presented variant was introduced in the
Master thesis of the second author [4].
2 An index with superlinear space
In accordance with the title we present here a truly bloated index (later referred
to as FM-bloated), namely requiring O(nlog2n) bits of space. Like in any FM-
index, we sort all suxes of T, but instead of taking the single symbols standing
before the suxes in the sorted order (which form the BWT of T), we take all
suchq-grams for all relevant values qbeing the power of 2. Namely, for each
suxT[i+ 1:::n] we take the q-grams:T[i] (1-gram), T[i1:::i] (2-gram),
T[i3:::i] (4-gram), ..., T[i2blogic+ 1:::i]. LetQdenote a collection of all
suchq-grams for all i. For each distinct item xfromQwe create a list Lxof its
occurrences in the sorted sux order (later we use the term sux array order
or simply SA order ). This resembles an inverted index on q-grams, yet the main
dierence is that the elements of the lists are arranged in SA order rather than
text order in a traditional inverted index.
For a given pattern P[1:::m ], we start the LF-mapping with its longest sux
of length being a power of 2, namely, P[m2blogmc+ 1:::m ]. The following
backward steps deal with the remaining prex of Pin a similar way. Note that
the number of LF-mapping steps is the number of 1s in the binary representation
ofm, which isO(logm). Figure 1 illustrates this approach.
In our representation, each LF-mapping step translates to performing two
predecessor queries on a list Lx. A na ve solution is a binary search with O(logn)
worst-case time (or even a linear search, which may turn out to be faster for very
short lists), yet the predecessor query can essentially be handled in O(log logn)
time [16] using an y-fast trie [19]. In terms of cache misses we can however
resort to the the optimal search in the I/O comparison model [18], obtaining
O(logjCLjlog=logn(n)) cache misses.
We also test a variant in which the values of qare Fibonacci numbers instead
of powers of 2, with the starting values of 1 and 2. For example, if m= 20, the
pattern is decomposed to substrings of length 13, 5, and 2.
Bloated1-Count-Occs( n,C,L,P,m)
1i m
2q0 2blogic
3 sp 1;ep n
4while ((spep)and (i1)do
5qGram P[iq0+ 1:::i]
6sp C[qGram ] +Occ(LqGram;sp1) + 1
7ep C[qGram ] +Occ(LqGram;ep)
8i iq0
9q0 2blogic
10 if(ep<sp)then return \not found" else return \found ( epsp+ 1) occs"
Fig. 1. Counting the number of occurrences of pattern PinTwith the
superlinear-space FM-bloated variant. The function Occ(
) is dened as
Occ(Lx;pos) = arg max
i2f1;:::;jLxjg(Lx[i]pos). To simplify the notation, we assume
that log 0 = 0 (line 9).
3 An index with linear space
In this variant we use at most two q-grams preceding any T[i+ 1:::n] sux,
hence reducing the overall space to O(nlogn) bits.
The solution makes use of so-called minimizers, proposed in 2004 by Roberts
et al. [17] and seemingly rst overlooked in the bioinformatics (or string match-
ing) community, only to be revived in the last years (e.g., [3, 6, 11, 13, 20]). A
minimizer for a sequence Sof lengthris the lexicographically smallest of its all
(rp+ 1)p-grams (or p-mers, in the term commonly used in bioinformatics);
usually it is assumed that pr. For a simple example, note that two DNA
sequencing reads with a large overlap are likely to share the same minimizer, so
they can be clustered together. That is, the smallest p-mer may be the identier
of the bucket into which the read is then dispatched.
In our variant, we pass a sliding window of length qoverTand calculate
the lexicographically smallest substring of length pin each window (i.e., its
minimizer). Ties are resolved in favor of the leftmost of the smallest substrings.
The positions of minimizers are start positions of the q-grams into which the
text is partitioned (note it resembles the recently proposed SamSAMi index,
Bloated2-Count-Occs( n,C,L,P,m,q,p)
1 nd 1i1<:::<i m0msuch that
eachP[ij] is a start position of a ( q;p)-minimizer in P
2 sp 1;ep n
3forj mdowntoim0do /* pattern sux */
4sp C[j] +Occ(LP[j];sp1) + 1
5ep C[j] +Occ(LP[j];ep)
6forj m01downto 1do
7qGram P[ij:::ij+11]
8sp C[qGram ] +Occ(LqGram;sp1) + 1
9ep C[qGram ] +Occ(LqGram;ep)
10 forj i11downto 1do /* pattern prex */
11sp C[j] +Occ(LP[j];sp1) + 1
12ep C[j] +Occ(LP[j];ep)
13 if(ep<sp)then return \not found" else return \found ( epsp+ 1) occs"
Fig. 2. Counting the number of occurrences of pattern PinTwith the linear-
space FM-bloated variant. The function Occ(
) is dened as Occ(Lx;pos) =
arg max
i2f1;:::;jLxjg(Lx[i]pos).
a sampled sux array on minimizers [12]). The values of qandp,pq, are
construction-time parameters.
Let us assume that 1 i1< i2< ::: < i n0np+ 1,n0n, are the
start positions of the obtained q-grams inT. We add to QallnunigramsT[i],
1in, and for each position ij, 2j <n0, also theq-gramT[ij1:::ij1]
assuming its length is greater than 1 (as all unigrams from Tare already added).
For a given pattern P[1:::m ] we pass a sliding window over it nding the
start positions of all its minimizers, which partition the pattern into several
(variable-length) q-grams. IfPoccurs inT, all thoseq-grams must belong to
Q. The following LF-mapping, together with the discussion on representation
of theLxlists, is identical as in the previous variant. The LF-mapping on the
pattern boundaries works on (zero or more) unigrams. Figure 2 illustrates this
approach.
Assume that the average distance between the start positions of the found
minimizers in T(or inP) is about ( qp+ 1)=2; let us denote it with q0. This
gives usO(m=q0) LF-mappings on average, if q0m. As we don't know mat
the index construction time, we cannot choose a proper q0(or ratherqandp)
to guarantee o(m) LF-mappings on average. This can be helped with increasing
the index space by a factor of (log logn). To this end, we look for minimizers
overT(andP, in the query time) several times; rst for some constant q1(and
naturally a constant p1,p1< q1), then for q2=q2
1, then forq3=q2
2, and so
on, as long as the current window size is not greater than n. Assuming that
the average distances q0
ibetween the start positions of the minimizers also grow
approximately quadratically, there exists q0
jsuch thatq0
jpm=2 andq0
jm=4.
This givesO(m11=2) =O(m1=2) LF-mappings per pattern on average.
4 Implementation and experimental results
As regards the implementation, the index uses a hash table in order to store (for
each selected q-gram) the value of the count table and the occurrence positions
which are used by the Occ(
) function. Rather than having bitvectors of length n
(e.g., as the wavelet tree components), as is usually the case for the FM-index, we
explicitly store only the occupied positions (as 32-bit integers), since otherwise
the space requirements would be prohibitive ( O(n2logn) bits). The search itself
is either a linear search if the list size is, roughly speaking, not much greater than
one cache line, or otherwise a binary search (BS). Moreover, in order to speed
up this search, we store an additional quick access (QA) list. This list points to
values corresponding to consecutive percentiles on the list | we rst perform a
linear search on a short QA list, and then jump to the specic position on the
positions list. Optimal values such as the length of the QA list and the binary
search threshold (i.e. the minimum list length which is required to use the BS
instead of a linear search) were determined empirically. The buckets, which store
theq-grams and pointers to their respective lists, have a contiguous layout, as
illustrated in Figure 3.
The experimental results were obtained on a machine equipped with the
Intel i5-3230M processor running at 2.6 GHz and 8 GB of DDR3 memory, and
the C++ source code was compiled (as a 32-bit version) with clang v. 3.6.2-1
and run on the Ubuntu 15.04 OS. In the following paragraphs, we present the
results for the superlinear variant.
As regards the hash function, xxhash ( https://code.google.com/p/xxhash/ )
was used. We have considered two pattern splitting schemes, namely one based
on the powers of 2, and the other based on the sequence of Fibonacci numbers.
The results for both cases dier signicantly due to dierent relations between
the sizes of stored q-grams and specic patterns, still, we can see the charac-
teristic spikes, with the lowest search times reported for pattern lengths that
matched the length of one of the stored q-grams (i.e. when it was equal to the
power of 2 or to one of the Fibonacci numbers). The times are given per charac-
ter and they represent the average times calculated for 106queries which were
extracted from the input text.
In order to evaluate both schemes, we have calculated the total sum of all
search times for the pattern lengths from 10 to 128 (both inclusive). The former
turned out to produce a relatively smaller index (around 80 n) with the total
time being equal to 283.18 ms, and the latter produced a relatively bigger index
(around 95n) with the lower total time of 225.91 ms; consult Figure 4 for details.
In Figure 5, we can see a comparison with other FM-index-based struc-
tures. We used the implementations from the sdsl library ( https://github.
com/simongog/sdsl-lite ) and the implementations of FM-dummy structures
by Grabowski et al. [10] ( https://github.com/mranisz/fmdummy/releases/
Fig. 3. Layout of the hash table used by the FM-bloated index for full-text
indexing. We can see the insertion of a q-gram \ab". The words are shown for
illustrative purposes, and in the actual implementation only pointers to the orig-
inal text are stored. The buckets contain q-gram sizes (stored as 8-bit characters)
as well as pointers to their respective lists. In this case, the value of the count
tableC[\ab"] is equal to 2 (rst value on the list), and there are 5 positions
(the number of positions is stored as the second value on the list) which will
be used for the calculation of Occ(
), namely 4, 5, 11, 13, and 77 (all arbitrary
values). Moreover, we can see the quick access (QA) list of length 2 whose be-
ginning and end are indicated by the dashed lines (the length of each QA list
is xed | in reality, QA lists are useful only for much longer positions lists,
e.g., ones that do not t into a single cache line). Let us note that this lay-
out is based directly on the layout that we have devised for another quick-access
text-based data structure called a split index [5]. Adapted from Wikimedia Com-
mons (author: Jorge Stol; available at http://en.wikipedia.org/wiki/File:
Hash_table_3_1_1_0_1_0_0_SP.svg ; CC A-SA 3.0).
tag/v1.0.0 ). As regards the space usage, the FM-bloated structure (just as the
name suggests) is roughly two orders of magnitude bigger than other indexes
(the index size for other methods ranged from approximately 0 :6nto 4:25n).
Acknowledgement
The work was supported by the Polish National Science Centre upon decision
DEC-2013/09/B/ST6/03117 (the rst author).
References
1. D. Belazzougui and G. Navarro. Alphabet-independent compressed text indexing.
ACM Transactions on Algorithms , 10(4):article 23, 2014.
2. A. Chac on, J. C. Moure, A. Espinosa, and P. Hern andez. n-step FM-index for faster
pattern matching. In Proc. International Conference on Computational Science
ICCS , pages 70{79. Elsevier, 2013.
3. R. Chikhi, A. Limasset, S. Jackman, J. T. Simpson, and P. Medvedev. On the
representation of de Bruijn graphs. Journal of Computational Biology , 22(5):336{
352, 2015.
4. A. Cis lak. Full-text and keyword indexes for string searching. Master's thesis, TU
M unchen, 2015. CoRR abs/1508.06610.
5. A. Cis lak and Sz. Grabowski. A practical index for approximate dictionary match-
ing with few mismatches. CoRR , abs/1501.04948, 2015.
6. S. Deorowicz, M. Kokot, Sz. Grabowski, and A. Debudaj-Grabysz. KMC 2: fast
and resource-frugal k-mer counting. Bioinformatics , 31(10):1569{1576, 2015.
7. P. Ferragina and R. Grossi. The String B-tree: A new data structure for string
search in external memory and its applications. J. ACM , 46(2):236{280, 1999.
8. P. Ferragina and G. Manzini. Opportunistic data structures with applications. In
Proc. FOCS , pages 390{398. IEEE, 2000.
9. P. Ferragina and G. Manzini. Indexing compressed text. Journal of the ACM ,
52(4):552{581, 2005.
10. S. Grabowski, M. Raniszewski, and S. Deorowicz. FM-index for dummies. CoRR ,
abs/1506.04896, 2015.
11. Sz. Grabowski, S. Deorowicz, and  L. Roguski. Disk-based compression of data
from genome sequencing. Bioinformatics , 31(9):1389{1395, 2015.
12. Sz. Grabowski and M. Raniszewski. Sampling the sux array with minimizers.
Accepted to SPIRE, 2015.
13. Y. Li, P. Kamousi, F. Han, S. Yang, X. Yan, and S. Suri. Memory ecient minimum
substring partitioning. In Proc. VLDB , pages 169{180. VLDB Endowment, 2013.
14. G. Navarro. Wavelet trees for all. Journal of Discrete Algorithms , 25:2{20, 2014.
15. G. Navarro and V. M akinen. Compressed full-text indexes. ACM Computing
Surveys , 39(1):2, 2007.
16. M. P atra scu and M. Thorup. Predecessor search. In M. Kao, editor, Encyclopedia
of Algorithms . Springer, 2008.
17. M. Roberts, W. Hayes, B. R. Hunt, S. M. Mount, and J. A. Yorke. Reducing storage
requirements for biological sequence comparison. Bioinformatics , 20(18):3363{
3369, 2004.
18. J. S. Vitter. Algorithms and data structures for external memory. Foundations
and Trends R
in Theoretical Computer Science , 2(4):305{474, 2008.
19. D. E. Willard. Log-logarithmic worst-case range queries are possible in space (n).
Information Processing Letters , 17(2):81{84, 1983.
20. D. E. Wood and S. L. Salzberg. Kraken: ultrafast metagenomic sequence classi-
cation using exact alignments. Genome Biology , 15(3):R46, 2014.
0 20 40 60 80 100 120 140
Pattern length050100150200250300Query time per char (ns)
0 20 40 60 80 100 120 140
Pattern length050100150200250Query time per char (ns)Fig. 4. Query time per character vs pattern size (from 10 to 128, inclusive) for
the superlinear FM-bloated index for powers of 2 (upper diagram) and Fibonacci
numbers (lower diagram) for the English text of size 25 MB. A hash table with
quick access lists was used, which allowed for a slight speedup with respect to a
hash table without said lists at the cost of a negligible increase in index size.
20 40 60 80 100 120
Pattern length101102103104Query time per char (ns)
FM-bloated superlinear
FM-index Huffman WT
CSA compressed bit vector
FM-dummy 2 256c
FM-dummy WT4 512Fig. 5. Query time per character vs pattern size (16, 24, 64, 80, and 120) for
dierent methods for the English text of size 25 MB. The FM-bloated superlinear
version with powers of 2 was used. Note the logarithmic y-scale.