Wavelet Tree, Part II: Text Indexing

The CSA represents the text string $𝒯$ by encoding the neighbor function $\mathrm{\Psi }$ [12, 28]. The sequence of $\mathrm{\Psi }$ values suffices to recreate the original text $𝒯$, using some supporting data structures. For example, if we know the symbol that appears at the beginning of the $i$th suffix in lexicographic order, we can easily determine the neighboring symbol in the text (at position $\mathrm{𝑆𝐴}[i]+1$ in the text) by computing $j=\mathrm{\Psi }(i)$. For each $\omega \in \mathrm{\Sigma }$, if we keep the count of the number of symbols in $𝒯$ smaller than or equal to $\omega$ in lexicographic order, then the desired symbol is simply the first symbol in lexicographic order whose count is $\ge j$. We can then continue taking steps forward from one neighbor to the next to recover additional symbols.

To encode $\mathrm{\Psi }$, we conceptually partition $\mathrm{\Psi }$ into context lists according to each symbol $\omega \in \mathrm{\Sigma }$. The context $\omega$ of $\mathrm{\Psi }(i)$ is defined as the symbol $𝒯[\mathrm{𝑆𝐴}[i]]$ that begins the suffix starting in position $\mathrm{𝑆𝐴}[i]$. We denote a $j$-context as a length-$j$ string. Suppose we have two suffixes at indices $i$ and $i^{\prime }$ (where ) in the suffix array that both start with the same symbol (context) $\omega$. Then if we remove the first symbol $\omega$ from each suffix, the two suffixes that remain must maintain the same relative ordering in the suffix array:

We show the neighbor function $\mathrm{\Psi }$ for text $𝒯$ = tcaaaatatatgcaacatatagtattagattgtat# (see Table 1). The suffixes starting with the symbol $𝚊$ span indices $[1,15]$ of the suffix array, and their $\mathrm{\Psi }$ values form an increasing sequence $⟨2$, 4, 5 ,11, 18, 19, 22, 23, 25, 27, 28, 29, 32, 34, $35⟩$, which we call a’s context list. For a given $k$ (to be specified later), we further partition each $\omega$’s context list into sublists $⟨\omega ,𝗑⟩$ according to $k$-context $𝗑\in {\mathrm{\Sigma }}^k$. Here $k$-context $𝗑$ of $\mathrm{\Psi }$ denotes the $k$-symbol prefix of the suffix starting at position $\mathrm{𝑆𝐴}[\mathrm{\Psi }(i)]$ and its preceding symbol is $\mathrm{𝑆𝐴}[i]$. For the text $𝒯$ shown in Table 1, the symbol $\mathrm{𝑆𝐴}[0]=\text{\#}$ precedes the $3$-context $\mathrm{𝚝𝚌𝚊}$ of the suffix starting at $\mathrm{𝑆𝐴}[\mathrm{\Psi }(0)]=\mathrm{𝑆𝐴}[31]$. As another example, the 2-context at associated with indices $[8,15]$ further decomposes a’s context list into the sublist $⟨23,25,27,28,29,32,34,35⟩$. The entries in each sublist also form an increasing sequence.

The next property is key for efficient coding of the neighbor function $\mathrm{\Psi }$:

For coding purposes, this property allows us to normalize each entry $\mathrm{\Psi }(i)$ in each sublist $⟨\omega ,𝗑⟩$ from 1 to $\mathrm{𝑡𝑜𝑡𝑎𝑙}(𝗑)$, where $\mathrm{𝑡𝑜𝑡𝑎𝑙}(𝗑)$ is the number of suffixes whose first $k+1$ symbols are of the form $y𝗑$ for any $y\in \mathrm{\Sigma }$. For example in Table 1, for $𝗑=𝚝$, the union of the sublists $⟨\omega ,𝚡⟩$ for all $\omega$ is the contiguous range $[23,35]$ of 13 values.

In Section 4 we will see how the ability of Wavelet Trees to encode individual sublists with 0th-order entropy results in an cumulative higher-order compression for the full text $𝒯$.

The FM-index involves accessing the Burrows-Wheeler transform (BWT [1]), which appears as the $L$ column in the example in Table 1. Adjacent entries in the $L$ column tend to share similar statistics regarding which symbols follow because the symbols that follow them in the text are adjacent in lexicographic order. As an example, $L[8,15]$ is the substring of $L$ formed by all symbols of $𝒯$ that precede the $2$-context $\mathrm{𝚊𝚝}$, as you can observe in Table 1.

The close relation between the CSA and BWT follows from the fact that the neighbor function $\mathrm{\Psi }$ is the inverse of the so-called $\mathrm{𝐿𝐹}$ function of the BWT (which maps the index of a symbol in the $L$ column to the index where it appears when rotated forward in the $F$ column). In other words, $\mathrm{\Psi }(i)=j$ iff $\mathrm{𝐿𝐹}(j)=i$, as can be seen from Table 1.

Building $\mathrm{𝐿𝐹}$ is straightforward for symbols that occur only once, as it is the case of # in our running example of Table 1. But computing $\mathrm{𝐿𝐹}$ efficiently is no longer trivial regarding recurring symbols. Nonetheless, it can be solved in optimal $𝒪(n)$ time thanks to a clever algorithm that uses an auxiliary vector $C$ of size $\sigma$ and is shown in Algorithm 1. For the sake of description, we assume that array $C$ is indexed by a symbol rather than by an integer. The first for-loop in Algorithm 1 computes, for each symbol $c$, the number $n_c$ of its occurrences in $L$, and thus it sets $C[c]=n_c$. Then, the second for-loop turns these symbol-wise occurrences into a cumulative sum, namely . We notice that $C[c]$ gives the first position in $F$ where the symbol $c$ occurs. Therefore, before the last for-loop starts, $C[c]$ is the landing position in $F$ of the first $c$ in $L$ (we thus know the $LF$-mapping for the first occurrence of every alphabet symbol). Finally, the last for-loop scans the column $L$, and whenever it encounters the symbol $L[i]=c$, it sets $\mathrm{𝐿𝐹}[i]=C[c]$ and then increments $C[c]$. So the algorithm maintains the invariant that , after that $k$ occurrences of $c$ in $L$ have been processed. The time complexity of such computation is $𝒪(n)$.

Like the neighbor function $\mathrm{\Psi }$ in the CSA, the BWT transform suffices to recreate the original text $𝒯$. To do so, it uses a backward scan supported by the $\mathrm{𝐿𝐹}$-mapping in $𝒪(n)$ time and space. The algorithm starts from the last symbol of $𝒯$, which occurs at $L[0]$; and then it proceeds by moving one symbol at a time to the left in $𝒯$, deploying the properties of the $\mathrm{𝐿𝐹}$-mapping: if we know the symbol in the $F$ column at index $i$, then the preceding symbol in the text is the symbol in the $L$ column on the same row $i$. Given this observation, the backward reconstruction of $𝒯$ first maps the current symbol occurring in $L$, say $L[i]$, to its corresponding copy in $F$, hence $j=LF[i]$; and then takes the symbol $L[j]$ that is ensured by the properties of the cyclic rotations of the rows of the BWT to precede $F[j]=L[i]$ in $𝒯$ (see above). This double step, which returns the algorithmic focus on $L$, allows to move one symbol leftward in $𝒯$. Hence, repeating this for $n$ steps, we can reconstruct the original input string $𝒯$ in $𝒪(n)$ time. As far as the construction of the BWT is concerned, we mention that it can be done by building the suffix array data structure via a plethora of time and space-efficient algorithms that now abound in the literature (see e.g. [17, 16, 22]) or directly via several elegant approaches (see e.g. [14, 25, 18, 23]).

But how useful is it to permute $𝒯$ via the BW-Transform? If we read the first column $F$ of Table 1, we notice a sequence of possibly long runs of equal symbols. It is clear that by compressing that sequence with a simple RLE-compressor (namely, the one that substitutes each run of equal symbols $x^h$ with a pair $⟨x,h⟩$) we would get a high compression ratio. But $F$ does not allow us to return to the original string $𝒯$ because all strings consisting of 15 a, 3 c, 4 g, and 13 t symbols would obtain the same $F$. It can be proved [1] that the only column allowing to return to $𝒯$, for all possible $𝒯$, is properly $L$. That column could be highly compressible because it satisfies the so-called locally homogeneous property, namely, that substrings of $L$ consist of a few distinct symbols. The reason comes from the fact that if we take a string (context) of length $k$ (as in the Definition 3 in PART I of $k$th-order entropy), say $𝗑\in {\mathrm{\Sigma }}^k$, the $n_{𝗑}$ symbols ${𝒯}_{𝗑}$ immediately preceding each occurrence of $𝗑$ in $𝒯$ occur contiguously in $L$. The nice issue here is that many real sources (they are called Markovian) do exist that generate data sequences, other than texts, that can be turned to be locally homogeneous via the Burrows-Wheeler Transform and thus can be highly compressed by simpler compressors.

This basic principle is deployed by the design of the well-known $\mathrm{𝖻𝗓𝗂𝗉}$ tool (now¹¹1See https://en.wikipedia.org/wiki/Bzip2 and https://github.com/kspalaiologos/bzip3 arrived to version 3), whose algorithmic pipeline we briefly sketch: Let us for a moment focus on two simple algorithms, namely, the Move-To-Front (shortly, MTF) and the Run-Length Encoding (shortly RLE, mentioned above). The former maps symbols into integers, and the latter maps runs of equal symbols into pairs. We observe that RLE is a compressor indeed, because the output sequence may be reduced in length in the presence of long runs of equal symbols; while MTF is a transform that can be turned into a compressor by encoding the integers via proper variable-length encoders. In general, the compression performance of those algorithms is very poor; while BWT turns them magically to be very effective if combined together: first apply MTF on $L$, then input the result on RLE, and finally encode its output with a statistical coder like Huffman or arithmetic coding. That’s it²²2For a scholarly discussion of the BWT and the FM-index, the reader is directed to the book [2]..

Manzini [21] showed that this, and other pipelines based upon the BWT, can represent a text string $𝒯$ of $n$ symbols using space proportional to $n$ times its $k$th-order empirical entropy ${\mathscr{H}}_k(𝒯)$, with guarantees stronger than the ones ensured by Lempel-Ziv-based compressors (such as gzip). Later analysis using Wavelet Trees reduced the constant factor in front of the $n{\mathscr{H}}_k(𝒯)$ term to 1 [11, 8]. Manzini’s result generated increased interest in the BWT transform, whose impact on data compression, and not just compressed text indexing, was highlighted by [5] as a sort of compression booster.

In the previous section, we decomposed the string being compressed (such as the neighbour sequence $\mathrm{\Psi }$ of the CSA or the $L$ column in the FM Index) into contexts. Each individual context ${𝒯}_{𝗑}$ was then encoded to its 0th-order entropy, and by Definition 3 in PART I of $k$th-order entropy, the sum of the resulting individual $n_{𝗑}{\mathscr{H}}_0(T_{𝗑})$ terms resulted in a total of $n{\mathscr{H}}_k(𝒯)$ bits (the boosting effect). In some cases, the results hold for all $k$ in the allowed range for some constant .

Foschini et al. [9] showed how order $k$th-order entropy can be achieved by using a single Wavelet Tree without explicit partitioning into $k$-contexts. Intuitively, there are two reasons:

• $\mathrm{■}$

  One reason is that the individual Wavelet Trees that would be constructed for each context can be found within the single Wavelet Tree constructed on the entire BWT string $L$, as illustrated in Figure 1.

• $\mathrm{■}$

  The other reason follows from using a method like run-length encoding (RLE) or gap encoding (Gap) (in conjunction with a prefix coder such as the $\gamma$ or $\delta$ code to encode the lengths) in order to encode the bit arrays of the nodes of the single Wavelet Tree. On the plus side, there is thus no need for data structures to keep track of the statistics of many individual contexts. On the negative side, the coding method RLE or Gap does not know the $k$-context boundaries. However, Foschini et al. [9] showed that RLE can adapt to the statistics for each $k$-context it encounters, and the positives make up for the negatives. By superimposing hypothetical $k$-context boundaries on the string for purposes of analysis, the resulting RLE implicitly encodes each $k$-context in roughly twice $0$th-order space, and by Definition 3 in PART I of $k$th-order entropy, the resulting overall encoding achieves $2n{\mathscr{H}}_k({𝒯}_k)$ bits leading space term. In effect, it achieves implicit boosting.

In implicit boosting, there is no need to specify $k$ in advance; the compression results hold for all $k$ in the allowed range ,  . Foschini et al. [9] used RLE encoding with $\gamma$ encoding, but other coding methods could be used, possibly switching dynamically among several methods (along with the bits that indicate the choices made) [3, 13, 15].

Gog et al. [10] compress the BWT by using the RRR method [26] to encode the Wavelet Trees of fixed-size blocks of the BWT rather than using partitions based upon higher-order contexts; they showed that the resulting compression achieves the $n{\mathscr{H}}_k(𝒯)$ leading space term. Mäkinnen and Navarro [19, 20] showed the surprising result that RRR [26] applied to a single Wavelet Tree of the entire BWT sequence or $\mathrm{\Psi }$ sequence (i.e., without any partitioning) also achieves the $n{\mathscr{H}}_k({𝒯}_k)$ bits leading space term, for a similar reason as in [9] described above; careful analysis in [19, 20] showed that the sublocks formed by [26] within each Wavelet Tree node implicitly encode each $k$th-order context in 0th-order entropy space plus lower-order terms. Grossi, Vitter, and Xu [13] present interesting experiments on the practical performance of these approaches.