Burrows-Wheeler transform (BWT)

Deep dive

How it works

In: a block of bytes (100 KB to tens of MB) · Out: the last column L of the sorted rotations plus the row index of the original string · State: none between blocks · Cost: forward O(n) with SA-IS and about 5n bytes of memory; inverse O(n) with one dependent random access per byte

1. Append a unique sentinel $ (smaller than all bytes) or record the original row instead.
2. Sort the suffixes of s$ (equivalent to sorting rotations); this is the suffix array.
3. Output L[i] = the character before suffix SA[i] ($ for SA[i] = 0); remember the row of $.
4. Inverse (LF mapping): C[c] = number of characters < c; LF(i) = C[L[i]] + (number of L[i] in L[0..i]). Starting at the row beginning with $, repeatedly emit L[i] and jump to LF(i); this yields s backwards.
5. bzip2 back end: move-to-front (recent bytes → small numbers), run-length code zero runs, Huffman with up to 6 tables switched every 50 symbols.

s = banana$      sorted rotations     L
                 $banana              a
                 a$banan              n
                 ana$ban              n
                 anana$b              b
                 banana$              $   <- primary row 4
                 na$bana              a
                 nana$ba              a
BWT = "annb$aa"  (runs: nn, aa)

Code

/// Forward BWT: (last column without the sentinel, row of the sentinel).
/// The suffix sort is O(n^2 log n) worst case; use SA-IS / libdivsufsort for real blocks.
pub fn bwt(s: &[u8]) -> (Vec<u8>, usize) {
    let n = s.len();
    let mut sa: Vec<usize> = (0..n).collect();
    sa.sort_unstable_by(|&a, &b| s[a..].cmp(&s[b..])); // same order as suffixes of s$
    let (mut out, mut primary) = (Vec::with_capacity(n), 0);
    if n > 0 { out.push(s[n - 1]); } // row 0 is the suffix "$", preceded by the last byte
    for (row, &p) in sa.iter().enumerate() {
        if p == 0 { primary = row + 1; } else { out.push(s[p - 1]); }
    }
    (out, primary)
}

pub fn inverse_bwt(l: &[u8], primary: usize) -> Vec<u8> {
    let n = l.len();
    let at = |r: usize| if r < primary { l[r] } else { l[r - 1] }; // row r != primary
    let mut count = [0usize; 256];
    for &c in l { count[c as usize] += 1; }
    let (mut first, mut sum) = ([0usize; 256], 1); // first row starting with c; row 0 is "$"
    for c in 0..256 { first[c] = sum; sum += count[c]; }
    let (mut lf, mut seen) = (vec![0usize; n + 1], [0usize; 256]);
    for r in (0..=n).filter(|&r| r != primary) {
        let c = at(r) as usize;
        lf[r] = first[c] + seen[c]; // i-th c in L is the i-th c in F
        seen[c] += 1;
    }
    let (mut out, mut r) = (vec![0u8; n], 0);
    for k in (0..n).rev() { out[k] = at(r); r = lf[r]; } // walk backwards through s
    out
}

pub fn mtf_encode(data: &[u8]) -> Vec<u8> {
    let mut order: Vec<u8> = (0..=255).collect();
    data.iter().map(|&b| {
        let i = order.iter().position(|&x| x == b).unwrap();
        order[..=i].rotate_right(1); // move b to the front
        i as u8
    }).collect()
}

pub fn mtf_decode(codes: &[u8]) -> Vec<u8> {
    let mut order: Vec<u8> = (0..=255).collect();
    codes.iter().map(|&i| {
        let b = order[i as usize];
        order[..=i as usize].rotate_right(1);
        b
    }).collect()
}

Complexity and constants

• Forward: suffix sorting dominates. libdivsufsort / SA-IS run in O(n) or O(n log n) with about 5n bytes of memory (n bytes input + 4n suffix array) for blocks < 4 GB.
• Inverse: 4-8 bytes per input byte for the LF array and n+1 dependent loads. In L2 (≤ 1 MB blocks) that is a few ns/byte; for 100 MB blocks each step is a DRAM miss (~80-100 ns) unless several blocks or streams are interleaved.
• MTF output on BWT text is dominated by 0s and 1s (often 50-60 % zeros), which is what makes the run-length + entropy stages effective.

Efficiency

• Cost: forward: suffix sorting dominates, O(n) with SA-IS or libdivsufsort at about 5n bytes of memory (n input + 4n suffix array) for blocks under 4 GB; bzip2 compresses at order of 10-20 MB/s. Inverse: 4-8 bytes per input byte for the LF array and n+1 dependent loads, a few ns per byte when the block fits L2 (1 MB or less) and a DRAM miss (about 80-100 ns) per byte on 100 MB blocks. Bound: memory latency on decode, compute on encode.
• Levers:
  • linear-time suffix sorting (SA-IS, libsais) — O(n) instead of O(n² log n) on repetitive input — every encoder
  • interleave 2-4 inverse chains split at stored checkpoint rows — overlapping loads hide DRAM latency — large blocks
  • bigger blocks (libbsc: tens of MB) — much better ratio than bzip2's 900 KB — when memory allows
  • adaptive context model instead of MTF + Huffman — better ratio (bsc) — ratio-critical archives
  • sampled rank tables and SA entries (FM-index) — pattern count in |P| steps independent of text size — search
• Naive → standard → state of the art:
  • sort full rotations · O(n² log n) on runs such as aaaa…, minutes for MB blocks
  • bzip2: 900 KB blocks, MTF + RLE + Huffman · order of 10-20 MB/s compress, ratio between gzip and xz
  • SA-IS/libsais forward, interleaved inverse, context-model back end (libbsc) · O(n), tens of MB blocks

Implementation notes

• Block size trades ratio against memory and cache behaviour: bzip2 uses 900 KB at -9; libbsc goes to tens of MB with much better ratios.
• Speed up the inverse by decoding 2-4 independent chains at once: split the walk at known rows (store a few checkpoint rows) so loads overlap.
• Better back ends than MTF: modern BWT compressors (bsc) code the BWT with adaptive context models / arithmetic coding and skip MTF.
• Run-length before BWT (bzip2 RLE1) exists to cap the cost of sorting long runs with the old comparison sort; with SA-IS it is unnecessary.
• FM-index: store L plus sampled rank tables (occurrences of c up to every 64th row) and sampled SA entries. Counting a pattern P takes |P| backward-search steps of 2 rank queries each, independent of text size.

Pitfalls

• Naive rotation sort (sort with full rotation comparison) is O(n² log n) on repetitive input such as aaaa… or zero pages: minutes for MB blocks.
• Forgetting that the sentinel must be strictly smallest and unique; using byte 0 as a sentinel breaks on binary data.
• Off-by-one between "row of $" and "row of the original string" conventions (bzip2 stores the original-string row); mix them and the output is a rotation of the input.
• The inverse's random access pattern makes large blocks slow to decompress even when compression was fine.

Variants and alternatives

In production

Archives: bzip2 (and pbzip2/lbzip2), libbsc, Linux kernel images compressed with bzip2 (legacy). Genomics: BWA and Bowtie/Bowtie2 DNA read aligners (FM-index over the reference genome; Li and Durbin 2009, Langmead et al. 2009). Text and data search: SDSL and other succinct data structure libraries (FM-index).

Implementations

Sources

• M. Burrows, D. J. Wheeler, A Block-sorting Lossless Data Compression Algorithm, DEC SRC Research Report 124, 1994.
• P. Ferragina, G. Manzini, Opportunistic Data Structures with Applications, FOCS, 2000. paper
• G. Nong, S. Zhang, W. H. Chan, Linear Suffix Array Construction by Almost Pure Induced-Sorting, DCC, 2009. paper
• H. Li, R. Durbin, Fast and accurate short read alignment with Burrows-Wheeler transform, Bioinformatics, 2009. paper · code
• B. Langmead, C. Trapnell, M. Pop, S. L. Salzberg, Ultrafast and memory-efficient alignment of short DNA sequences to the human genome, Genome Biology, 2009. paper · code