The Burrows–Wheeler Transform and Backward Search in the FM-Index

The Burrows–Wheeler transform is a permutation of a string’s characters that makes the string easier to compress and, unexpectedly, easier to search. It is invertible with no side information beyond a single index, it does not itself reduce length, and its output is not obviously related to its input. Both properties come from one operation: sorting all cyclic rotations of the text.

The transform is the basis of bzip2’s compression stage and of the FM-index used in read aligners and full-text search. The two uses read the same structure for different purposes — compression uses the clustering of equal characters, search uses the fact that the sorted rotations are the sorted suffixes of the text.

Sorting rotations, reading the last column

Take a text T of length n and append a sentinel character $ that is smaller than every character in T and occurs nowhere else. The sentinel makes all n+1 rotations distinct and gives a unique smallest rotation. Sort the n+1 rotations lexicographically into a matrix, and output the last column L. The first column F is the sorted characters of the text and is recoverable by sorting L.

The sentinel is not a convenience. Without it, a text such as abab has rotations abab, baba, abab, baba — duplicates, so the sort is ambiguous and the transform is not invertible. Any injective termination works, including an out-of-band length, but a sentinel character is the form that keeps the arithmetic uniform.

The reason the output clusters is visible in the matrix. Rows are sorted by rotation, so rows beginning with the same character are adjacent, and within that group rows are sorted by what follows. A row’s last character is the character preceding its first character in the original text. So L lists the predecessors of the first characters in sorted order: wherever the text has many occurrences of a character in similar contexts, the predecessors of those contexts are also similar, and L contains runs of the same symbol. Natural language produces long runs of spaces before common words, and DNA produces runs of the same base before common k-mers. The runs are what the subsequent move-to-front and run-length stages of a compressor exploit. The transform is a permutation, so the character multiset is unchanged and L has the same zero-order entropy as T; what changes is the local order, which is what the later stages turn into compression.

The LF mapping inverts the transform

The inverse uses a correspondence between L and F that is a bijection. In row i, F[i] is the first character of the i-th rotation and L[i] is the character before it. The rotation starting at position j+1 (mod n+1) has first character T[j+1] and last character T[j]; so the last character of row i equals the first character of the row whose rotation starts one position earlier in the text.

Sorting the occurrences of each character in L by their order of appearance, and matching them to the occurrences of the same character in F in sorted order, gives the mapping. Formally,

LF(i) = C[L[i]] + Occ(L[i], i)

where C[c] is the number of characters smaller than c in the text and Occ(c, i) is the number of occurrences of c in L[0..i−1]. LF maps a row to the row one rotation-step earlier. Applying it repeatedly from the row containing the sentinel walks backwards through the text and emits the characters in reverse order, which is the whole inverse transform.

def bwt(s):
    s = s + '\0'                     # sentinel smaller than any byte
    n = len(s)
    rotations = sorted(range(n), key=lambda i: s[i:] + s[:i])
    return ''.join(s[(i - 1) % n] for i in rotations), rotations

def inverse_bwt(last):
    n = len(last)
    counts, start = {}, {}
    for ch in sorted(last):
        counts[ch] = counts.get(ch, 0) + 1
    total = 0
    for ch in sorted(counts):
        start[ch] = total
        total += counts[ch]
    seen, lf = {}, [0] * n
    for i, ch in enumerate(last):
        lf[i] = start[ch] + seen.get(ch, 0)
        seen[ch] = seen.get(ch, 0) + 1
    out, row = [], 0                     # row 0 holds the rotation starting with the sentinel
    for _ in range(n - 1):
        out.append(last[row])            # emit before stepping: the first step yields T[n-1]
        row = lf[row]
    return ''.join(reversed(out))

Emitting before stepping matters. Row 0’s rotation is $T, so its last character is the final character of the text; stepping with LF first and emitting afterwards shifts the whole output by one position and silently drops a character. The loop runs n−1 times because the sentinel itself is not part of the text.

The construction of the transform by sorting rotations is O(n² log n) if each comparison materialises a rotation, which is why real implementations compute it from a suffix array instead: sorting the rotations of T$ is the same as sorting the suffixes of T$, and any suffix array construction gives the order in O(n log n) or O(n). The inverse is linear in n with the counting pass, and the counting pass is also the LF table, computed once.

The same mapping as a search index

The FM-index observes that F is the first column of the sorted rotation matrix, and the sorted rotations are the sorted suffixes of the text. A pattern P occurs in T exactly where a suffix of T begins with P, and those suffixes form a contiguous range of rows in the matrix. The FM-index finds that range by processing P from right to left.

Start with the range [0, n) covering all rows, and for each character c of P taken from last to first, narrow the range to rows whose first character is c and whose remaining prefix matches the already processed suffix of P. The narrowing is a rank query on L:

lo' = C[c] + Occ(c, lo)
hi' = C[c] + Occ(c, hi)

If lo’ equals hi’, the pattern does not occur. If the range survives all m characters, its width is the number of occurrences, and the rows in it are the starting positions in sorted order. Each step is two rank queries on L, so the search is O(m) rank operations with no dependence on the text length except through the rank structure’s own query cost.

The correctness of that recurrence is the LF mapping read in the other direction. Rows in [lo, hi) share the prefix P[i+1..m−1]. Applying LF to them gives the rows whose rotations are one position earlier, which share the prefix P[i..m−1], provided the character prepended is c. The rank query counts exactly the rows in the range whose last character is c, and C[c] offsets them into the block of rows starting with c.

What a rank query costs

Occ(c, i) is a rank query on a string over an alphabet of size σ. Three implementations are common.

A flat occurrence table stores, for every block of b positions and every character, the count at the block boundary; the counts within a block are recovered by scanning at most b characters of L. The boundary table costs (n/b)·σ entries. With b = 64, σ = 256, and 4-byte counters that is 16n bytes, sixteen bytes per text character, which is why practical indexes replace the table with a wavelet tree when the alphabet is large, or with a bitvector per character and rank support when the alphabet is small. The query is then a couple of memory accesses plus that scan.

A wavelet tree over L answers rank in O(log σ) with n log σ bits plus rank directories, and it answers range quantile and range reporting queries that a flat table cannot. Its space is closer to the information-theoretic minimum, at the cost of log σ steps per query and pointer-chasing through the levels.

A compressed rank structure — the original FM-index formulation — stores L in blocks compressed by the same BWT machinery that produced it, with sampled rank values at block boundaries. Query cost becomes O(log n) per rank in the simple form, or constant with more sampling, and total space can fall below the text size for repetitive input. This is the variant that makes the index attractive for genomes, where the text is large and highly repetitive.

Structure	Space	Rank query	Extra queries
flat table, block 64	16n bytes	O(1) + scan ≤ 64	none
wavelet tree	n log σ bits + overhead	O(log σ)	quantile, reporting
compressed blocks	below n bytes for repetitive text	O(log n) or O(1)	none

The flat table’s space is why FM-indexes in practice tune the block size against the counter width — a larger block shrinks the table but lengthens the in-block scan — and why the table is often paired with a sampled suffix array rather than used alone.

Locating rather than counting

Counting occurrences needs only the range. Reporting their positions needs the text positions themselves, which the matrix does not store. The standard solution samples every k-th suffix position in the suffix array and stores the sampled positions alongside. To locate a row, apply LF repeatedly until the row corresponds to a sampled position; the number of steps is at most k, and each step is one rank query. Locating one occurrence therefore costs O(k) rank queries, and locating all of them costs O(occurrences · k). The sample rate k is the space-time knob: a small k means a large sample table and fast locating, a large k the reverse. Sampling every 32 positions with 4-byte entries costs n/8 bytes, which is the usual compromise.

A second option stores the text itself and walks LF until the position is known by other means, which is only useful when the text is present anyway.

Limitations that shape real use

The transform is defined on a block, and the block must be sorted as a whole. For a text larger than memory, the suffix array construction is the bottleneck, not the transform. bzip2 restricts blocks to a size chosen from 100 kB to 900 kB by the compression level, so that the sort fits comfortably and the compressor’s memory bound is predictable; the format’s block size is encoded in the header, and the largest block is under a megabyte.

Backward search requires the pattern to be processed right to left, which is fine for a batch query and awkward for a streaming one. It also requires the rank structure over L, which must be rebuilt when the text changes; the FM-index is static in the same way the suffix array is. Appending text requires either rebuilding or a dynamic variant with a logarithmic penalty per update.

Finally, the sentinel must be a symbol outside the alphabet. With byte text, a zero byte is only safe if the input cannot contain one; with a general byte stream, an out-of-band length or an escape convention is required. Getting this wrong produces a transform that is invertible on the test data and silently ambiguous on inputs containing the sentinel value, a failure that shows up as corrupted output rather than as a crash.

The transform’s compression role and its search role are separable. bzip2 uses it purely for the clustering, then applies move-to-front, run-length encoding, and Huffman coding, and never builds a rank index. Search systems build the rank index and never compress. The structure that serves both is the sorted rotation matrix; what differs is which column is read and what is asked of it. For the dictionary-based alternative that trades the global sort for local matches, see LZ77 and LZ78 parsing.