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.

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.

Frequently Asked Questions