We describe an external memory sux array construction algorithm based on constructing sux arrays for blocks of text and merging them into the full sux array. The basic idea goes back over 20 years and there has been a couple of later improvements, but we describe several further improvements that make the algorithm much faster. In particular, we reduce the I/O volume of the algorithm by a factor O(log n). Our experiments show that the algorithm is the fastest sux array construction algorithm when the size of the text is within a factor of about five from the size of the RAM in either direction, which is a common situation in practice.
The sux array [
For internal memory, there exists an essentially optimal sux array construction
algorithm (SACA) that runs in linear time using little extra space in addition to
what is needed for the input text and the output sux array [
In this paper, we focus on an algorithm, which we call SAscan, that lies between internal memory SACAs and eSAIS in scalability. SAscan is a true external memory algorithm in the sense that it can handle texts that do not fit in internal memory, but its time complexity ⌦( n2/M ) makes it hopelessly slow when the text is much larger than the RAM. However, when the text is too large for an internal memory SACA, i.e., larger than about one fifth of the RAM size, but not too much bigger than the RAM, SAscan is probably the fastest SACA in practice. SAscan is also lightweight in the sense that it uses less than half of the disk space of eSAIS, and can be implemented to use little more than what is needed for the text and the sux array.
The basic approach of SAscan was developed already in the early days of sux
arrays by Gonnet, Baeza-Yates and Snider [
The early version of SAscan depended heavily on the text not having long repeats
and thus had a poor worst case time complexity. This problem was solved by
Crauser and Ferragina [
In this paper, we describe several further improvements to the SAscan algorithm.
The first improvement is a new merging technique that reduces the I/O complexity
of SAscan to O n2/(M B log n) + n/B and provides a substantial speedup in
practice too. Another target of improvement is the rank data structure that plays
a key role in the algorithm. In theory, we observe that the time complexity can
be reduced to O (n2/M ) log(2 + (log / log log n)) by plugging in the rank data
structure of Belazzougui and Navarro [
We show experimentally that our practical improvements reduce the running time of SAscan by more than a factor of three. We also show that the algorithm is faster than eSAIS when the text size is less than about six times the available RAM, at which point the disk space usage of eSAIS is already well over 150 times the available RAM. 2
Let X = X[0..m) be a string over an integer alphabet [0.. ). For i = 0, . . . , m 1 we write X[i..m) to denote the sux of X of length m i, that is X[i..m) = X[i]X[i + 1] . . . X[m 1]. Similarly, we write X[i..j) to denote the substring X[i]X[i + 1] . . . X[j 1] of length j i. If i = j, the substring X[i..j) is the empty string, also denoted by ".
The sux array SAX of X contains the starting positions of the non-empty
suxes of X in the lexicographical order, i.e., it is an array SAX[0..m) which contains
a permutation of the integers [0..m) such that X[SAX[0]..m) < X[SAX[
The partial sux array SAX:Y is the lexicographical ordering of the suxes of
XY with a starting position in X, i.e., it is an array SAX:Y[0..m) that contains a
permutation of the integers [0..m) such that X[SAX:Y[0]..m)Y < X[SAX:Y[
Let a string T[0..n) be the text. It is divided into blocks of size (at most) m, where m is chosen so that all the in-memory data structures of O(m log n) bits fit in the RAM. The blocks are processed starting from the end of the text. Assume that so far we have processed Y = T[i..n) and constructed the sux array SAY. Next we will construct the partial sux array SAX:Y for the block X = T[i m..i) and merge it with SAY to form SAXY.
The suxes in SAX:Y and SAY are in the same relative order as in SAXY and we just need to know how to merge them. For this purpose, we compute the gap array gapX:Y[0..m], where gapX:Y[i] is the number of suxes of Y that are lexicographically between the suxes SAX:Y[i 1] and SAX:Y[i] of XY. Formally, for i 2 [1..m), gapX:Y[0] = |{j 2 [0..|Y|) : Y[j..|Y|) < X[SAX:Y[0]..m)Y}| gapX:Y[i] = |{j 2 [0..|Y|) : X[SAX:Y[i
1]..m)Y < Y[j..|Y|) < X[SAX:Y[i]..m)Y}| gapX:Y[m] = |{j 2 [0..|Y|) : X[SAX:Y[m 1]..m)Y < Y[j..|Y|)}| .
The construction of gapX:Y scans Y and is the computational bottleneck of the algorithm.
The merging of SAX:Y and SAY is trivial with the help of gapX:Y but involves a lot of I/O for reading SAY and writing SAXY. The total I/O volume is O n2/m in units of O(log n)-bit words. However, we can reduce the I/O by delaying the merging. We write SAX:Y and gapX:Y to disk and then proceed to process the next block. Once all partial sux arrays and gap arrays have been computed, we perform one multiway merging of the partial sux arrays with the help of the gap arrays. The I/O volume for merging is reduced to O(n) words.
Suppose that during the construction of SAX:Y or gapX:Y we need to compare two suxes of SAXY, at least one of which begins in X. In the worst case, the suxes could have a very long common prefix, much longer than m, making it impossible to perform the comparison without a lot of I/O — unless we have extra information about the order of the suxes. In our case, that extra information is provided by a bitvector gtY, which tells whether each sux of Y is lexicographically smaller or larger than Y itself. Formally, for all i 2 [0..|Y|), gtY[i] =
The output bitvector gtXY is needed as input for the next block. The two other arrays SAX:Y and gapX:Y are stored on disk until all blocks have been processed. The final phase of the algorithm takes all the partial sux arrays and gap arrays as input and produces the full sux array SAT. 4
Details and Analysis
The first stage in processing a block X is constructing the partial sux array
SAX:Y. In the full paper, we show how to construct a string Z such that SAZ =
SAX:Y. We can then construct the sux array using any standard SACA; In the
implementation we use Yuta Mori’s divsufsort [
The partial Burrows–Wheeler transform [
BWTX:Y[i] =
⇢ X[SAX:Y[i]
$
if SAX:Y[i] > 0
if SAX:Y[i] = 0
where $ is a special symbol that does not appear in the text. For a character
c and an integer i 2 [0..m], the answer to the rank query rankBWTX:Y (c, i) is the
number of occurrences of c in BWTX:Y[0..i). Rank queries can be answered in
O(log(2 + (log / log log n))) time using a linear space data structure [
For a string S, let sufrankX:Y(S) be the number of suxes of XY starting in X that are lexicographically smaller than S. Let C[0.. ) be an array, where C[c] is the number of positions i 2 [0..m) such that X[i] < c. In the full paper, we prove the following lemma. Lemma 4.1 Let k = sufrankX:Y(S) for a string S. For any symbol c, sufrankX:Y(cS) = C[c] + rankBWTX:Y (c, k) +
Note that when S = Y[j..|Y|), we can replace the comparison Y < S with gtY[j] = 1. Thus, given sufrankX:Y(Y[j..|Y|)), we can easily compute sufrankX:Y(Y[j 1..|Y|)) using the lemma, and we only need to access Y[j 1] in Y and gtY[j] in gtY. Hence, we can compute sufrankX:Y(Y[j..|Y|)) for j = |Y| 1, . . . , 0 with a single sequential pass over Y and gtY. This is all that is needed to compute gapX:Y and gtXY, which are the output of the second stage of processing X.
The final phase of the algorithm is merging the partial sux arrays into the full sux array SAT. For k 2 [0..dn/me), let Xk = T[km..(k + 1)m), Yk = T[(k + 1)m..n), SAk = SAXk:Yk and gapk = gapXk:Yk . The algorithm shown below moves suxes from the input sux arrays to the output sux array in ascending lexicographical order. 1: for k = 0 to dn/me 2: for i = 0 to n 1 do 3: k 0 4: while gapk[ik] > 0 do 5: gapk[ik] gapk[ik] 6: k k + 1 7: SAT[i] SAk[ik] + km 8: ik ik + 1 1 do ik 1
0 The correctness of the algorithm is based on the following invariants maintained by the algorithm: (1) ik is the number of suxes already moved from SAk to SAT, and (2) gapk[ik] is the number of suxes remaining in SAk+1, SAk+2, . . . , SAdn/me 1 that are smaller than SAk[ik].
Theorem 4.2 SAscan can be implemented to construct the sux array of a
text of length n over an alphabet of size in O⇣ nM2 log ⇣2 + lolgolgog n ⌘⌘ time and
O MB log n + Bn log MB Bn ⌘ I/Os in the standard external memory model (see [
In the full paper, we describe an implementation that needs 5.2m bytes of RAM and at most 11.5n bytes of disk space, and could be implemented to use just 6.5n bytes of disk space. 5
Experimental Results We performed experiments on a machine with a 3.16GHz Intel Core 2 Duo CPU with 6144KiB L2 cache running Linux (Ubuntu 12.04, 64bit, kernel 3.2). All 4 6 Input size GiB ● ● ● 8 eSAIS bwtdisk SAscan ● ● 10 8 6 4 2 0 ● ● 2 ● ● ●
● 4 6 Input size GiB ● ● eSAIS bwtdisk SAscan ● ● 8 enwik hg programs were compiled using g++ version 4.6.4 with -O3 -DNDEBUG options. To reduce the time to run the experiments, we artificially restricted the RAM size to 2GiB using the Linux boot option mem, and the algorithms were allowed to use at most 1.5GiB. We used two big test files: a concatenation of three di↵erent Human genomes1,2,3 (hg) and a prefix of the English Wikipedia dump4 (enwik).
The first experiment measures the scalability of our new algorithm. We computed the sux array for varying length prefixes of each testfile using our algorithm and compared to eSAIS – currently the fastest algorithm for building SA in external memory. In addition, we also show the runtimes of bwtdisk5, which is essentially the Ferragina–Gagie–Manzini version of SAscan though it constructs the BWT instead of the sux array. The sux array version of bwtdisk would be slower as it needs more I/O during merging. The results are given in Figure 1. SAscan is faster than eSAIS up to input sizes of about 9GiB, which is about six times the size of the RAM available to the algorithms. It is this ratio between the input size and the RAM size that primarily determines the relative speeds of SAscan and eSAIS. Note that the main limitation to the scalability of eSAIS is the disk space requirement, which is about 170 times the available RAM at the crossing point. The runtime of bwtdisk is always at least 3.5 times the runtime of SAscan, showing the dramatic e↵ect of our improvements. In the second experiment, we take a closer look at the e↵ect of our improvements on the runtime. More precisely, we show a detailed runtime breakdown after turning on individual improvements one by one: the fast merging of partial sux arrays, the space-ecient representation of the gap array (which reduces the RAM usage from 8m to 5.2m bytes), and the optimized rank data structure. The results are presented in Figure 2. Each of the improvements produces a significant speedup. The combined e↵ect is more than a factor of three. 1http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/ 2ftp://public.genomics.org.cn/BGI/yanhuang/fa/ 3ftp.ncbi.nlm.nih.gov/genbank/genomes/Eukaryotes/vertebrates_mammals/Homo_ sapiens/ 4http://dumps.wikimedia.org/enwiki/ 5http://people.unipmn.it/manzini/bwtdisk/ enwik (4GiB) hg (4GiB) 4
B s i M3 e iTm 2 6 5 1 0 merge gap other 6 5 4 3 2 1 0 merge gap other e n i l e s a b e g r e m t s a f