Faster Run-Length Compressed Suffix Arrays

The key idea behind compressed suffix arrays (CSAs) is to store $\mathrm{\Psi }[0..n-1]$ compactly while supporting certain searches on it quickly, where $\mathrm{\Psi }[0..n-1]$ is the permutation of $\{0,\mathrm{\dots },n-1\}$ such that $\mathrm{\Psi }[i]$ is the position of SA entry $(\mathrm{SA}[i]+1)modn$ in $\mathrm{SA}[0..n-1]$ or, equivalently, the position in $L$ of $F[i]$. (This means $\mathrm{\Psi }$ is the inverse of the $\mathrm{LF}$ mapping used in FM-indexes.) By the definition, $\mathrm{\Psi }$ consists of at most $\sigma$ increasing intervals – one for each distinct character that occurs in the text, corresponding to the interval of suffixes starting with that character – and if we can support fast binary searches on these intervals then we can support fast pattern matching.

For example, consider the text

for which $\mathrm{SA}$, $\mathrm{\Psi }$, $F$ and $L$ are shown on the left in Figure 1. If we know $\mathrm{SA}[\mathrm{22..28}]$ is the SA interval for CG (in the green rectangle) and we want the SA interval for GCG, then we can search in the increasing interval

for G (in the red rectangle, with $\mathrm{\Psi }$ values between 22 and 28 shown as orange arrows and the others shown as black arrows) for the successor $\mathrm{\Psi }[41]=23$ of 22 and the predecessor $\mathrm{\Psi }[43]=28$ of 28. We thus learn that the SA interval for GCG is $\mathrm{SA}[\mathrm{41..43}]$ (in the blue rectangle). Knowing this, we can continue backward stepping.

Run-length compressed suffix array (RLCSA) were introduced in [24] for indexing highly repetitive collections. In this section we present an alternative, but functionally equivalent, description of RLCSAs which is more suitable for describing our improvements.

In this paper we view a RLCSA as a data structure storing ${\mathrm{\Psi }}^{\prime }[0..r-1]$ compactly while supporting certain searches on it quickly. By the definition of ${\mathrm{\Psi }}^{\prime }$, it still consists of at most $\sigma$ increasing intervals – one for each distinct character that occurs in $T$, corresponding to the interval of suffixes starting with that character – and if we can still support fast binary searches on these intervals then we can still support fast pattern matching.

For example, consider

again, for which ${\mathrm{\Psi }}^{\prime }$, $F^{\prime }$ and $L^{\prime }$ are shown on the right in Figure 1. If we know the SA interval $\mathrm{SA}[\mathrm{22..28}]$ for CG starts at offset 0 in the $L$ run of character $L^{\prime }[12]$ and ends at offset 1 in the $L$ run of character $L^{\prime }[15]$ (in the green rectangle) and we want the SA interval for GCG, then we can search in the increasing interval

for G (in the red rectangle, with ${\mathrm{\Psi }}^{\prime }$ values between 12 and 15 shown as orange arrows and the others shown as black arrows) for the successor ${\mathrm{\Psi }}^{\prime }[28]=13$ of 12 and the predecessor ${\mathrm{\Psi }}^{\prime }[29]=15$ of 15.

Because $L^{\prime }$ and $F^{\prime }$ do not have the predecessor-successor relationship of $L$ and $F$, we cannot deduce that the SA interval for GCG starts in the $L$ run of character $L^{\prime }[28]$ and ends in the $L$ run of character $L^{\prime }[29]$ (and, in fact, in this example it does not). Instead, we store two $n$-bit SD-bitvectors [21], $B_L$ and $B_F$, with $r$ copies of $1$ each. The 1s in $B_L$ mark the starting positions of runs in $L$ and the 1s in $B_F$ mark the positions in $F$ of the marked characters in $L$. In our example

The interval $B_F[\mathrm{41..43}]$ in $B_F$ starting immediately before the bit with offset 0 in the block whose starting position is marked with the 29th copy of 1 and ending immediately before the bit with offset 1 in the block whose starting position is marked with the 30th copy of 1, is shown in blue. (We are interested in the blocks marked with the 29th and 30th copies of 1 because we count from 0 in the $j$ column in Figure 1, so those blocks correspond to ${\mathrm{\Psi }}^{\prime }[28]$ and ${\mathrm{\Psi }}^{\prime }[29]$.) We can find this interval with 2 ${\mathrm{select}}_1$ queries on $B_F$, which take constant time.

The corresponding interval $B_L[\mathrm{41..43}]$ in $B_L$ is also shown in blue, starting immediately before the bit with offset 3 in the block whose starting position is marked with the 22nd copy of 1 and ending immediately before the bit with offset 1 in the block whose starting position is marked with the 24th copy of 1. We can find the 2 indices 22 and 24 with 2 ${\mathrm{rank}}_1$ queries on $B_L$, which take $O(\log \log n)$ time. This means the SA interval for GCG is $\mathrm{SA}[\mathrm{41..43}]$ and it starts at offset 3 in the $L$ run of character $L^{\prime }[21]$ and ends at offset 1 in the $L$ run of character $L^{\prime }[23]$. Knowing this we can continue backward stepping.

The RLCSA in Sirén’s PhD thesis [24] for a text $T[0..n-1]$ with $r$ BWT runs takes $O\left(\text{<span class="ltx\_rule" style="width:0.0pt;height:8.6pt;background:black;display:inline-block;"></span>}r\log (n/r)+r\log \sigma +\sigma \log n\right)$ bits. Given a character $a$ and the SA interval for $P$, it can find the SA interval for $aP$ in $O(\log n)$ time.

Suppose we are given an increasing list ${\mathrm{\ell }}_1,\mathrm{\dots },{\mathrm{\ell }}_k$ of $k$ integers in the range $[0..n-1]$. To encode them with interpolative coding [19], we first write ${\mathrm{\ell }}_{\lceil k/2\rceil }$ using $\lfloor \lg (n-1)\rfloor +1$ bits (except that we write 0 using 1 bit). All the numbers ${\mathrm{\ell }}_1,\mathrm{\dots },{\mathrm{\ell }}_{\lceil k/2\rceil -1}$ are in the range $[0..{\mathrm{\ell }}_{\lceil k/2\rceil }-1]$, so we can encode them recursively. All the numbers ${\mathrm{\ell }}_{\lceil k/2\rceil +1},\mathrm{\dots },{\mathrm{\ell }}_k$ are in the range $[{\mathrm{\ell }}_{\lceil k/2\rceil }+1..n-1]$, so we can encode them recursively as ${\mathrm{\ell }}_{\lceil k/2\rceil +1}-{\mathrm{\ell }}_{\lceil k/2\rceil }-1,\mathrm{\dots },{\mathrm{\ell }}_k-{\mathrm{\ell }}_{\lceil k/2\rceil }-1$. Each encoding has $O(\log n)$ bits, so we can read them in $O(1)$ time. If we imagine the list stored as keys in a balanced binary search tree then we encode the keys according to a pre-order traversal: when we reach each key ${\mathrm{\ell }}_i$, we know ${\mathrm{\ell }}_i$ lies between the numbers shown to the left and right of ${\mathrm{\ell }}_i$ and we encode ${\mathrm{\ell }}_i$ using the maximum number of bits we would need for any key in that range.

For example, if $n=66$, $k=13$ and the list is $6,9,14,15,16,23,24,28,29,30,42,46,63$ then, as illustrated in Figure 2, we start by encoding ${\mathrm{\ell }}_7=24$ using $\lfloor \lg 65\rfloor +1=7$ bits as 0011000. We then encode ${\mathrm{\ell }}_3=14$ using $\lfloor \lg 23\rfloor +1=5$ bits as 01110. We then encode ${\mathrm{\ell }}_1,{\mathrm{\ell }}_2,{\mathrm{\ell }}_5,{\mathrm{\ell }}_4,{\mathrm{\ell }}_6=6,9,16,15,23$ as $\mathrm{𝟶𝟷𝟷𝟶},\mathrm{𝟶𝟷𝟶},\mathrm{𝟶𝟶𝟶𝟷},\mathrm{𝟶},\mathrm{𝟷𝟷𝟶}$, and ${\mathrm{\ell }}_{10},{\mathrm{\ell }}_8,{\mathrm{\ell }}_9,{\mathrm{\ell }}_{12},{\mathrm{\ell }}_{11},{\mathrm{\ell }}_{13}$ as $\mathrm{𝟶𝟶𝟶𝟷𝟶𝟷},\mathrm{𝟶𝟷𝟷},\mathrm{𝟶},\mathrm{𝟶𝟶𝟷𝟷𝟷𝟷},\mathrm{𝟷𝟶𝟷𝟷},\mathrm{𝟷𝟶𝟶𝟶𝟶}$. When we reach 46, say, in a pre-order traversal of the tree in Figure 2, we know it lies between 31 and 65, so we encode it using $\lfloor \lg (65-31)\rfloor +1=6$ bits as ${(46-31)}_2=\mathrm{𝟶𝟶𝟷𝟷𝟷𝟷}$.

The binary search tree has height $\lfloor \lg k\rfloor$ and the bottom level contains at most $k$ keys. By Jensen’s Inequality, we encode those keys using $O(k\log (n/k)+k)$ bits. Similarly, there at most $k/2^h$ keys at height $h$ and we encode those keys using $O\left(\frac{k}{2^h}\log \frac{n}{k/2^h}+\frac{k}{2^h}\right)=O\left(\frac{k}{2^h}\log (n/k)+\frac{k(h+1)}{2^h}\right)$ bits. Since

we use $O(k\log (n/k)+k)$ bits in total.

In this paper we want to perform binary search on the list – the reader may have noticed that our example is $\mathrm{\Psi }[\mathrm{36..48}]$ from Figure 1 – so we want fast access to the encodings of the numbers in it in the order we check them in a binary search. We can store the encodings according to an in-order traversal instead of a pre-order traversal, and store an uncompressed bitvector with as many bits as there are in the concatenation of the encodings and 1s marking where the encodings start. Since the bitvector delimits the encodings, however, we can delete the leading 0s from each encoding before concatenating them and building the bitvector. The in-order encodings for our example are shown below the tree in Figure 2, with the leading 0s removed, and the bitvector is shown below them. Since the bitvector uses at most as many bits as the encodings, we still use $O(k\log (n/k)+k)$ bits in total and – although random access still takes $O(\log k)$ time – we can perform binary search in $O(\log k)$ total time. This scheme is similar to Teuhola’s [25] and Claude, Nicholson and Seco’s [5].

To find the successor of 22 in the list, we start at the root knowing $n=66$ and $k=13$ and perform ${\mathrm{select}}_1(7)$ and ${\mathrm{select}}_1(8)$ queries on the bitvector to find the starting and ending positions of the encoding 0011000 of ${\mathrm{\ell }}_7=24$ in the range $[\mathrm{0..65}]$. Since , we then perform ${\mathrm{select}}_1(3)$ and ${\mathrm{select}}_1(4)$ queries to find the starting and ending positions of the encoding 01110 of ${\mathrm{\ell }}_3=14$ in the range $[\mathrm{0..23}]$. Since $22>14$, we then perform ${\mathrm{select}}_1(5)$ and ${\mathrm{select}}_1(6)$ queries to find the starting and ending positions of the encoding 0001 of ${\mathrm{\ell }}_5=16$ in the range $[\mathrm{15..23}]$. Since $22>16$, we then perform ${\mathrm{select}}_1(6)$ and ${\mathrm{select}}_1(7)$ queries to find the starting and ending positions of the encoding 110 of ${\mathrm{\ell }}_6=23$ in the range $[\mathrm{17..23}]$. Since , we know the successor of 22 in $L$ is 23. We can find the predecessor of 28 in $O(\log k)$ time symmetrically.

If we apply interpolative coding with fast binary search to the increasing interval of $\mathrm{\Psi }$ for a character $a$ in a text $T$ of length $n$, then we use $O(n_a\log (n/n_a)+n_a)$ bits and can support binary search in $O(\log n_a)$ time, where $n_a$ is the frequency of $a$ in $T$. If we do this for all the characters then we use $O(n(H_0(T)+1))$ bits, where $H$ is the 0th-order empirical entropy of $T$. If we encode the increasing interval of ${\mathrm{\Psi }}^{\prime }$ for $a$ with interpolative coding, then we use $O(r\log (r/r_a)+r_a)$ bits and can support binary search in $O(\log r_a)$ time, where $r_a$ is the number of runs of copies of $a$ in the BWT of $T$ (and, equivalently, in $L$). If we do this for all the characters then we use $O(r(H_0(L^{\prime })+1))$ bits, where $L^{\prime }$ is again the sequence of $r$ characters in the runs of the run-length encoding of $T$. To be able to find the increasing interval for $a$ in ${\mathrm{\Psi }}^{\prime }$, we store an $r$-bit uncompressed bitvector with 1s marking where the intervals start.

To use Theorem 4 in an RLCSA, we store

• $\mathrm{■}$

  an uncompressed bitvector marking which distinct characters occur in $T$, in $O(\sigma )$ bits;

• $\mathrm{■}$

  the SD-vectors $B_F$ and $B_L$ in $O(r\log (n/r))$ bits;

• $\mathrm{■}$

  an uncompressed bitvector with 1s marking where the intervals for the characters start in ${\mathrm{\Psi }}^{\prime }$, in $O(r)$ bits;

• $\mathrm{■}$

  ${\mathrm{\Psi }}^{\prime }$ in $O(r\log \sigma )$ bits.

If we are given $a$ and the SA interval for $P$ then we can find the SA interval for $aP$ by

• $\mathrm{■}$

  using a rank query on the first uncompressed bitvector to find $a$’s rank among the distinct characters that occur in $T$, in $O(1)$ time;

• $\mathrm{■}$

  using rank queries on $B_L$ to find the runs in $L$ overlapping the SA interval for $P$, in $O(\log \log n)$ time;

• $\mathrm{■}$

  using select queries on the second uncompressed bitvector to find the interval for $a$ in ${\mathrm{\Psi }}^{\prime }$, in $O(1)$ time;

• $\mathrm{■}$

  using binary search in the interval for $a$ in ${\mathrm{\Psi }}^{\prime }$ to find the successor and predecessor of the ranks of the first and last runs in $L$ overlapping the SA interval for $P$, in $O(\log r_a)$ time;

• $\mathrm{■}$

  using select queries on $B_F$ and arithmetic to find the SA interval for $aP$ in $O(1)$ time.

We store $O(r\log (n/r)+r\log \sigma +\sigma )$ bits in total and find the SA interval for $aP$ in $O(\log r_a+\log \log n)$ total time. Notice that the $O(\log \log n)$ term in the query time comes only from the rank query on $B_L$.

Nishimoto and Tabei [20] showed how we can split the runs in $L$ such that no block in $B_F$ overlaps more than a constant number of blocks in $B_L$ without increasing the number of runs by more than a constant factor, and then store $\mathrm{LF}$ in $O(r\log n)$ bits and evaluate it in constant time. Brown, Gagie and Rossi [4] slightly generalized their key theorem:

If $L[i]=L[i-1]$ then $\mathrm{LF}(i)=\mathrm{LF}(i-1)+1$, so

has cardinality $r$. If $\mathrm{LF}(i)=\mathrm{LF}(i-1)+1$ then, since $\mathrm{\Psi }$ and $\mathrm{LF}$ are inverse permutations, $\mathrm{\Psi }[j]=\mathrm{\Psi }[j-1]+1$ where $j=\mathrm{LF}(i)$. Therefore,

also has cardinality $r$ and applying Theorem 6 with $d=2$ to $\mathrm{\Psi }$ splits the runs in the BWT such that no block in $B_F$ overlaps more than 3 blocks in $B_L$, without increasing the number of runs by more than a factor of $3/2$. In our example, the number of runs increases by only 1, from 40 to 41, as shown below with the split block – corresponding to the first run of 6 copies of T in $L$ – in red:

Suppose we apply Theorem 6 with $d=2$ to $\mathrm{\Psi }$ and then store, for , the index of the block in $B_L$ containing $\mathrm{LF}(i_b)$ and $\mathrm{LF}(i_b)$’s offset in that block, where $i_b$ is the starting position of block $b$ in $B_L$. Nishimoto and Tabei called this the move table for $\mathrm{LF}$ (see also [4, 26]) and it takes a total of $O(r\log n)$ bits. If we know $B_L[j]$ is in block $b$ in $B_L$ with offset $j-i_b$ then, since the block in $B_F$ to which $\mathrm{LF}$ maps block $b$ in $B_L$ now overlaps at most the block containing $B_L[\mathrm{LF}(i_b)]$ and the next 2 blocks in $B_L$, we can find the index of the block in $B_L$ containing $B_L[\mathrm{LF}(j)]=B_L[\mathrm{LF}(i_b)]+j-i_b$ and $B_L[\mathrm{LF}(j)]$’s offset in that block with at most 2 constant-time select queries on $B_L$. We could use at most 2 constant-time lookups instead if we have the starting positions of the blocks in $B_L$ stored explicitly in another $O(r\log n)$ bits.

Recall that the $O(\log \log n)$ term in the query-time bound in Corollary 5 comes only from the use of rank queries on an SD-vector. Since rank and select queries can be combined to support predecessor queries and select queries on sparse bitvectors can easily be supported in constant time and space polynomial in the number of 1s, rank queries on compact sparse bitvectors inherit lower bounds from predecessor queries [3] – so they cannot be implemented in constant time. Therefore, to get rid of that $O(\log \log n)$ term, we must somehow avoid rank queries.

We could replace the rank queries with a move table, but that would result in an $O(r\log n)$ term in our space bound. Instead, we introduce an uncompressed $2r$-bit bitvector $B_{FL}$ indicating how the starting positions of the blocks in $F$ and $L$ are interleaved. Specifically, we scan $B_F$ and $B_L$ simultaneously – assuming we have already applied Theorem 6 to them so that no block in $F$ overlaps more than 3 blocks in $L$ (so $r$ is a constant factor larger than it was before the application of the theorem) – and

• $\mathrm{■}$

  if we see 0s in both bitvectors in position $i$ then we write nothing;

• $\mathrm{■}$

  if we see a 1 in $B_F$ and a 0 in $B_L$ then we write a 1 (indicating that a block starts in $F$);

• $\mathrm{■}$

  if we see a 0 in $B_F$ and a 1 in $B_L$ then we write a 0 (indicating that a block starts in $L$);

• $\mathrm{■}$

  if we see 1s in both bitvectors then we write a 0 and then a 1 (indicating that blocks start in both $L$ and $F$).

This way, $B_{FL}.{\mathrm{select}}_1(j)$ tells us which at most 3 blocks in $L$ – those corresponding to the 0 preceding the $j$th copy of 1 in $B_{FL}$ and possibly to the next 2 copies of 0 – could overlap block $j$ in $F$ (counting from 1). We can then find the starting positions of those blocks in $L$ using at most 3 select queries on $B_L$.

For our example, taking $B_F$ and $B_L$ to be as shown just after Theorem 6,

(with the grey numbers only to show positions) we have

In position 38 we see 1s in both $B_F$ and $B_L$, so we write 01 in $B_{FL}$ (in positions 49 and 50, respectively); in positions 39 and 40 we see 0s in both in $B_F$ and $B_L$, so we write nothing; in position 41 we see a 1 in $B_F$ and a 0 in $B_L$, so we write a 1 in $B_{FL}$; in position 42 we see a 0 in $B_F$ and a 1 in $B_L$, so we write a 0 in $B_{FL}$; in position 43 we see 1s in both $B_F$ and $B_L$, so we write 01 in $B_{FL}$; in position 44 we see 0s in both $B_F$ and $B_L$, so we write nothing; and in position 45 we see a 0 in $B_F$ and a 1 in $B_L$, so we write a 0 in $B_{FL}$ (in position 55). Admittedly, when $n=66$ and after applying Theorem 6 $2r=82$, it seems foolish to store a $2r$-bit uncompressed bitvector instead of simply storing $B_L$ uncompressed. This is due to the small size of our example, however; for massive and highly repetitive datasets, $r$ can easily be hundreds of times smaller than $n$.

Suppose we know the SA interval $\mathrm{SA}[\mathrm{41..43}]$ for $aP$ starts at offset 0 in block 28 in $F$ and ends at offset 1 in block 29 in $F$ and we want to find which blocks contain its starting and ending positions in $L$ and the offsets of those positions. In Section 2, we performed 2 rank queries on $B_L$, but now we perform queries $B_{FL}.{\mathrm{select}}_1(29)=51$ and $B_{FL}.{\mathrm{select}}_1(30)=54$ (with arguments 29 and 30 instead of 28 and 29 because we mark with a 1 the starting of the first block in $F$, which we index with 0; the results 51 and 54 are indexed from 0 as well). Since the 29th and 30th copies of 1 are $B_{FL}[51]$ and $B_{FL}[54]$ (shown in red above), they are preceded by the $51-29+1=23$rd and $54-30+1=25$th copies of 0, respectively.

Because we applied Theorem 6, this means the 29th and 30th blocks in $F$ (shown in red in $B_F$ above) overlap the 23rd block in $L$ and possibly the 24th and 25th blocks (shown in blue in $B_L$), and the 25th block and possibly the 26th and 27th blocks (also shown in blue in $B_L$). Notice that, because we split the 34th block in $F$ but the first block in $L$ for Theorem 6, the block numbers we find in $F$ are the same as in Section 2 but the block numbers we find in $L$ will be incremented. Although in general we need 6 select queries on $B_L$, in this case we can use only 5 – $B_L.{\mathrm{select}}_1(23),\mathrm{\dots },B_L.{\mathrm{select}}_1(27)$ – to find where these blocks begin in constant time, and determine which contain the starting and ending positions of the SA interval $\mathrm{SA}[\mathrm{41..43}]$: the 23rd and the 25th, respectively.

In short, we replace a rank query on SD-bitvector $B_L$ by queries on uncompressed bitvector $B_{FL}$ and constant-time select queries on $B_L$. This gives us the following theorem:

Instead of viewing $B_{FL}$ as replacing slow rank queries while using the overall same space, we can also view it (and $B_F$ and $B_L$) as replacing an $O(r\log n)$-bit move table while using the same overall query time. Brown, Gagie and Rossi [4] implemented a similar approach to speeding up $\mathrm{LF}$ computations in an RLFM-index, but only alluded to it briefly in their paper – the path to Bitvector in their Figure 3 – and gave no analysis nor bounds. We conjecture that a similar approach can also be applied to reduce the size of fast move tables for $\varphi$ and ${\varphi }^{-1}$ [13], which return $\mathrm{SA}[i-1]$ and $\mathrm{SA}[i+1]$ when given $\mathrm{SA}[i]$.