Burrows-Wheeler Transform — Missing Marker Corruption

bzip2 regularly outcompresses gzip by 10-30% on text data. The reason is the Burrows-Wheeler Transform — a reversible rearrangement that clusters similar contexts together, making subsequent compression dramatically more effective. Michael Burrows and David Wheeler published it in 1994 as a DEC Systems Research Center technical report, and it became the core of bzip2.

The BWT also has a second life in bioinformatics. The FM-index — a compressed data structure built on BWT — allows searching a 3-billion-base genome for any short read in milliseconds while using only 1-2 GB of RAM. The entire field of short-read genome alignment (BWA, Bowtie2, HISAT2 — the tools that process most DNA sequencing data) is built on Burrows-Wheeler.

Why Burrows-Wheeler Transform Is Not a Compression Algorithm

The Burrows-Wheeler Transform (BWT) is a reversible permutation that rearranges a string into runs of repeated characters — it does not compress data. The core mechanic: generate all cyclic rotations of the input, sort them lexicographically, then take the last column of the sorted matrix. That last column, plus the index of the original string in the sorted list, is enough to reconstruct the input exactly. The magic is that sorting groups identical characters together, creating long runs that downstream compressors (move-to-front, Huffman, arithmetic coding) exploit efficiently.

BWT runs in O(n log n) time for sorting, but practical implementations use suffix arrays to achieve O(n) with a good constant. The transform itself adds no compression — it only reorganizes. The real value is that after BWT, the entropy of the string drops significantly because runs of identical bytes become common. For example, on English text, BWT can increase the number of zero-byte runs by 10-100x, making subsequent entropy coding vastly more effective.

Use BWT when you need high compression ratios on natural language, DNA sequences, or any data with local structure but no global patterns. It is the foundation of bzip2, which still outperforms gzip on many workloads. BWT is not suitable for small payloads (under a few hundred bytes) because the transform overhead dominates. It also requires the entire input in memory — streaming is possible but complex. In production, BWT-based compressors are a go-to for archival and distribution where decompression speed matters more than compression speed.

Forward BWT

1. Form all rotations of the string (append a unique end-of-string marker $)
2. Sort all rotations lexicographically
3. Output the last column (L) of the sorted matrix and the row index of the original string

def bwt(s: str) -> tuple[str, int]:
    """Compute BWT of s. Returns (transformed, original_row_index)."""
    s = s + '$'  # append end marker (lexicographically smallest)
    n = len(s)
    rotations = sorted(s[i:] + s[:i] for i in range(n))
    last_col = ''.join(r[-1] for r in rotations)
    original_idx = rotations.index(s)
    return last_col, original_idx

result, idx = bwt('banana')
print(f'BWT(banana) = {result!r}, original row = {idx}')
result2, idx2 = bwt('mississippi')
print(f'BWT(mississippi) = {result2!r}, original row = {idx2}')

Why BWT Helps Compression

Notice BWT(banana) = 'annb$aa' — the 'a's and 'n's are grouped. This happens because rotations starting with 'an' are sorted together, and those all have the character preceding 'an' in their last column.

In English text, the character before 'an' is usually 'b' (ban), 'c' (can), 'm' (man), 'p' (pan), etc. — similar characters cluster. bzip2 achieves 10-30% better ratio than gzip on English text because of this clustering.

The clustering effect is strongest when the alphabet is small and patterns repeat often — which is exactly the case for natural language and many structured data formats (e.g., JSON, XML).

Inverse BWT — The LF Mapping

Recovery uses the key insight: sorting the last column L gives the first column F. Since each row is a rotation, (F[i], L[i]) are consecutive characters in the original string. The LF mapping links each occurrence of a character in L to its position in F, effectively walking backwards through the original string.

Algorithm: 1. Build vector of positions: for each character, list its occurrence indices in L. 2. Build first column F = sorted(L). 3. For each character in F, maintain a pointer to the next occurrence in L. 4. Start at the row given by original_idx, and repeatedly output F[pos], then update pos = LF[pos]. 5. Stop after n steps (all characters output), then remove the '$' marker.

This is the same as constructing the full rotation table column by column, but the LF method is O(n) in time after O(n) preprocessing, versus O(n^2) for the naive table approach shown next.

def ibwt_lf(last: str, orig_idx: int) -> str:
    """Reconstruct original string from BWT using LF mapping O(n)."""
    n = len(last)
    # Build occurrence index for each character in L
    from collections import defaultdict
    occ = defaultdict(list)
    for i, ch in enumerate(last):
        occ[ch].append(i)
    # Build first column F and mapping LF
    first = sorted(last)
    # For each character in F, we need a pointer into its occurrence list
    # We'll build an array 'next_in_L' of length n
    next_in_L = [0] * n
    # Temporary dict: for each char, current pointer into occ[char]
    ptr = dict()
    for ch in occ:
        ptr[ch] = 0
    for i, ch in enumerate(first):
        # get the position in L of the ptr[ch]-th occurrence of ch
        pos = occ[ch][ptr[ch]]
        ptr[ch] += 1
        next_in_L[i] = pos
    # Now reconstruct: start at original_idx in F, follow next_in_L
    result = []
    pos = orig_idx
    for _ in range(n):
        result.append(first[pos])
        pos = next_in_L[pos]
    # The result includes the marker; remove it
    return ''.join(result).rstrip('$')

print(ibwt_lf('annb$aa', 3))        # banana
print(ibwt_lf('ipssm$pissii', 2))  # mississippi

BWT Pipeline — bzip2

bzip2 uses BWT as the first stage of a multi-step pipeline: 1. BWT: Rearrange for character clustering 2. Move-to-Front (MTF): Encode each character as its index in a list of recently seen characters — this turns clusters of the same character into runs of 0s. 3. RLE: Compress runs of zeros (from MTF) using run-length encoding. 4. Huffman coding: Entropy-code the result with multiple Huffman tables.

This pipeline achieves better compression than gzip (DEFLATE) on typical text at the cost of more memory and slower speed.

Bzip2 splits input into blocks (default 900 KB) and processes each independently. This allows streaming and limits memory per block, but reduces compression ratio slightly compared to processing the whole file.

# Simulate the bzip2 pipeline steps on the BWT output
bwt_result = 'annb$aa'

# Step 1: Move-to-Front transform
def mtf(s: str) -> list[int]:
    alphabet = list(set(s))
    alphabet.sort()  # deterministic order
    result = []
    for ch in s:
        idx = alphabet.index(ch)
        result.append(idx)
        # move to front
        alphabet.pop(idx)
        alphabet.insert(0, ch)
    return result

mtf_output = mtf(bwt_result)
print('MTF output:', mtf_output)  # many zeros in clustered regions

# Step 2: Simple run-length encode zeros
# (bzip2 uses a more sophisticated variant)
# Step 3: Huffman coding (not shown)
# The high frequency of zeros makes Huffman very efficient.

BWT in Bioinformatics — The FM-Index

The FM-index (Ferragina-Manzini index) uses the BWT of a reference genome to enable fast substring search with small memory footprint. Instead of storing the full genome, we store the BWT plus a few auxiliary data structures (suffix array samples, occurrence counts).

Key insight: The LF mapping used in inverse BWT also allows counting occurrences of a query string in the original text without decompressing. By walking the query backwards through the BWT, you can determine the range of rows where the query is a prefix.

This is the foundation of modern short-read aligners: BWA, Bowtie2, HISAT2. They can align 100 million reads to the human genome in a few hours using < 4 GB RAM.

Performance and Memory Trade-offs

BWT computation and inversion have different costs depending on the algorithm: - Naive rotation sort: O(n log n) time, O(n^2) memory — only for small strings. - Suffix array construction: O(n) to O(n log n) time, O(n) memory — production choice. Libraries like divsufsort (used by bzip2) are widely available. - Block-based processing: bzip2 processes 900 KB blocks independently, limiting memory to ~ 9 MB per block.

Memory usage during decompression is about 4x the block size due to the need for multiple stages.

Alternative: For better compression ratio, use xz (LZMA) which achieves 20-30% better than bzip2 but is even slower. For speed, zstd outperforms both gzip and bzip2 on most data.

Why the Sentinel Matters — Don’t Skip the $

Every BWT implementation needs a unique end-of-string marker. Usually '$'. It’s not a suggestion — it’s a requirement for reversibility. Without it, the LF mapping breaks and you’ll reconstruct garbage.

Here’s why: The sentinel guarantees that the last character of the original string is uniquely identified during the inverse transform. It forces the cyclic rotations to sort in a way that the LF-mapping property holds. No sentinel means you cannot distinguish between rotations that produce identical last columns — a fatal ambiguity.

In production, pick a character that is lexicographically smallest in your alphabet. For DNA sequences, use '\0' or some non-nucleotide character. For text, '$' is common. The sentinel adds exactly one byte of overhead per block — negligible for the correctness it buys.

The trap: reusing a character that already exists in the string. That breaks the transform. Always validate your input doesn’t contain the sentinel before encoding.

// io.thecodeforge — dsa tutorial

public class SentinelBWT {
    public static String encode(String input, char sentinel) {
        if (input.indexOf(sentinel) != -1) {
            throw new IllegalArgumentException(
                "Sentinel " + sentinel + " already exists in input");
        }
        String text = input + sentinel;
        int n = text.length();
        String[] rotations = new String[n];
        for (int i = 0; i < n; i++) {
            rotations[i] = text.substring(i) + text.substring(0, i);
        }
        java.util.Arrays.sort(rotations);
        StringBuilder lastColumn = new StringBuilder();
        for (String rot : rotations) {
            lastColumn.append(rot.charAt(n - 1));
        }
        return lastColumn.toString();
    }

    public static void main(String[] args) {
        String result = encode("banana", '$');
        System.out.println("BWT: " + result);
    }
}

LF-Mapping on the Ground — Reverse Without a Matrix

The inverse BWT seems magical until you trace the LF-mapping manually. It lets you walk backwards through the original string without storing the full rotation matrix. Here’s the mechanics:

LF-mapping says: for row i in the sorted rotations, the next row in the original sequence is given by the index of the character at the last column in the first column. More concretely: given last column L and first column F (which is just sorted L), rank[i] = count of character L[i] seen so far in L. Then LF[i] = C[L[i]] + rank[i], where C[c] is the number of characters in L that are lexicographically smaller than c.

Let’s reverse the BWT “annb$aa” from our banana example. F column (sorted L): $ a a a b n n. L: a n n b $ a a. Compute C: (before a: 1 char '$', before b: 4 chars a+$, before n: 5 chars a+b+$). Build ranks: a:0, n:0, n:1, b:0, $:0, a:1, a:2. LF[i] = C[L[i]] + rank[i]. Starting from row of '$' (index 4), follow LF links backwards to reconstruct "banana$".

This is O(n) in time and O(σ) in space for the C array. That’s why bzip2 can decompress gigabytes on a phone.

// io.thecodeforge — dsa tutorial

public class LFInverse {
    public static String inverseBWT(String lastColumn) {
        int n = lastColumn.length();
        char[] firstColumn = lastColumn.toCharArray();
        java.util.Arrays.sort(firstColumn);

        // Build C array
        int[] C = new int[256];
        for (char c : firstColumn) C[c]++;
        int cumulative = 0;
        for (int i = 0; i < 256; i++) {
            int count = C[i];
            C[i] = cumulative;
            cumulative += count;
        }

        // Build rank array
        int[] rank = new int[n];
        int[] seen = new int[256];
        for (int i = 0; i < n; i++) {
            char c = lastColumn.charAt(i);
            rank[i] = seen[c]++;
        }

        // Find sentinel row
        int row = 0;
        for (int i = 0; i < n; i++) {
            if (lastColumn.charAt(i) == '$') {
                row = i;
                break;
            }
        }

        // Reconstruct
        StringBuilder sb = new StringBuilder();
        for (int i = 0; i < n - 1; i++) {
            row = C[lastColumn.charAt(row)] + rank[row];
            sb.append(lastColumn.charAt(row));
        }
        return sb.reverse().toString();
    }

    public static void main(String[] args) {
        String original = inverseBWT("annb$aa");
        System.out.println("Original: " + original);
    }
}