第13章 Burrows–Wheeler 変換（やさしい版） — 関数プログラミングによるアルゴリズム設計の真珠

Burrows–Wheeler 変換（略して BWT）とは、リストの並び方を工夫して、同じ文字が近くに集まるように並べ替える手法です。データ圧縮の「前処理」として使うのが本来の目的です。

なぜ「近くに集める」と得？同じ文字が連続していれば、たとえば「a が 5 回連続」を「a×5」と書くだけで済みます（ランレングス符号化）。BWT を通した後の文字列に、こういう簡単な圧縮を掛け、さらに Huffman 符号や算術符号など高度な圧縮を掛けると、非常によく縮む、というわけです。

ただソートするだけではダメ？

「同じ文字を集めたいなら、単に文字順にソートすれば？」と思うかもしれません。でも、それでは元の文字列に戻す方法がなくなってしまいます。圧縮なのに復元できないのでは意味がありません。

本章では BWT を定義し、なぜ戻せるのかの本質的な理由を突き止め、そこから逆変換のアルゴリズムを導きます。

BWT の定義

やっていること

例として文字列 yokohama を考えてみましょう。左が回転一覧、右がそれを辞書順にソートしたものです。

0	y o k o h a m a
1	o k o h a m a y
2	k o h a m a y o
3	o h a m a y o k
4	h a m a y o k o
5	a m a y o k o h
6	m a y o k o h a
7	a y o k o h a m

0	a m a y o k o h
1	a y o k o h a m
2	h a m a y o k o
3	k o h a m a y o
4	m a y o k o h a
5	o h a m a y o k
6	o k o h a m a y
7	y o k o h a m a

図 13.1 回転（左）とソート後の回転（右）

Haskell での定義

transform :: Ord a ⇒ [a] → ([a], Int) transform xs = (map last xss, position xs xss) where xss = sort (rots xs)

Dart // 素朴版 transform: 回転を全部作ってソートし、最終列と元の行番号を返す。 (List<T>, int) transformNaive<T extends Comparable<Object?>>(List<T> xs) { final xss = rots(xs).toList() ..sort((a, b) => _lexCompare(a, b)); final lastCol = [for (final row in xss) row.last]; final k = position(xs, xss); return (lastCol, k); } int _lexCompare<T extends Comparable<Object?>>(List<T> a, List<T> b) { final n = a.length < b.length ? a.length : b.length; for (var i = 0; i < n; i++) { final c = a[i].compareTo(b[i]); if (c != 0) return c; } return a.length.compareTo(b.length); }

行の位置を返す関数と、回転を作る関数はこう定義します。

position xs xss = length (takeWhile (≠ xs) xss) rots :: [a] → [[a]] rots xs = take (length xs) (iterate lrot xs) where lrot (x : xs) = xs ++ [x]

Dart // position: xs と一致しない行を数え上げ、xs 自身の位置を返す。 int position<T>(List<T> xs, List<List<T>> xss) { for (var i = 0; i < xss.length; i++) { if (_listEq(xss[i], xs)) return i; } return xss.length; } bool _listEq<T>(List<T> a, List<T> b) { if (a.length != b.length) return false; for (var i = 0; i < a.length; i++) { if (a[i] != b[i]) return false; } return true; } // lrot: 先頭要素を末尾へ回す（左回転）。 List<T> lrot<T>(List<T> xs) => [...xs.skip(1), xs.first]; // rots: xs の全ての左回転を長さ順に返す。 Iterable<List<T>> rots<T>(List<T> xs) sync* { var cur = List<T>.of(xs); for (var i = 0; i < xs.length; i++) { yield List<T>.of(cur); cur = lrot(cur); } }

lrot は 1 回の左回転（先頭の文字を末尾へ）です。この transform は仕様通りには動きますが、効率は良くありません。ひとまず効率の話は後回しにします。²

なぜ圧縮に効くのか（直感）英語の文には “the”, “this”, “that” など h の前に t が来る単語や、“which”, “where”, “when” など h の前に w が来る単語がたくさん出ます。回転を辞書順に並べると、行の先頭が h で始まる行が固まります。そしてそれらの行の末尾（＝行の先頭文字の直前にあった文字）は、t や w に偏るのです。だから最終列を取ると、t や w が近くに集まって出てきます。

逆変換の仕様

逆変換 untransform は、変換したものを元に戻す関数と定めます。

untransform · transform = id

問題は、どうやって計算するかです。素朴なアイデアはこうです。

そこで、次を満たすような関数 recreate を作れたと仮定します。

recreate · map last · sort · rots = sort · rots

すると untransform (ys, k) = (recreate ys) !! k とすればよいのです。しかし、最終列 1 本だけから本当に行列全体を復元できるのでしょうか？——できます。しかもその導出は、プログラム計算の見事な演習になります。

復元の計算（Recreational calculation）

列ごとに復元していく作戦

入力の長さを n とします。recreate は n × n の行列を復元するわけですが、これを列ごとに 1 列ずつ足していく方針で組み立てます。最初の j 列だけを取り出す関数を用意します。

takeCols :: Int → [[a]] → [[a]] takeCols j = map (take j)

Dart // takeCols: 行列の各行から先頭 j 要素を取り出す。 List<List<T>> takeCols<T>(int j, List<List<T>> xss) => [for (final row in xss) row.take(j).toList()];

takeCols n は行列全体を取り出す（＝恒等関数）ので、目標を「recreate j で最初の j 列を復元する」に置き換えます。仕様は次のとおり。

recreate j · map last · sort · rots = takeCols j · sort · rots(13.1)

基底ケース（j = 0）

recreate 0 = map (const [ ])

Dart // recreate の基底: ys（最終列）と同じ本数の空行を返す。 List<List<T>> recreate0<T>(List<T> ys) => [for (var i = 0; i < ys.length; i++) <T>[]];

n 個の空リストが返ります。takeCols 0 も同じく n 個の空リストを返すので、この段階では OK です。

帰納ケースに必要な 3 つの道具

rrot xs = [last xs] ++ init xs hdsort = sortBy cmp where cmp (x : xs) (y : ys) = compare x y consCol (xs, xss) = zipWith (:) xs xss

Dart // rrot: 末尾要素を先頭へ回す（右回転）。 List<T> rrot<T>(List<T> xs) => [xs.last, ...xs.take(xs.length - 1)]; // hdsort: 各行の先頭要素だけで安定ソートする（同点は元の順序を保つ）。 List<List<T>> hdsort<T extends Comparable<Object?>>(List<List<T>> xss) { final indexed = [for (var i = 0; i < xss.length; i++) (i, xss[i])]; indexed.sort((a, b) { final c = a.$2.first.compareTo(b.$2.first); return c != 0 ? c : a.$1.compareTo(b.$1); // 同点は元順（安定化） }); return [for (final e in indexed) e.$2]; } // consCol: 列 xs を行列 xss の先頭列としてくっつける。 List<List<T>> consCol<T>((List<T>, List<List<T>>) pair) { final xs = pair.$1, xss = pair.$2; return [for (var i = 0; i < xs.length; i++) [xs[i], ...xss[i]]]; }

map rrot · rots = rrot · rots(13.2)

意味「全部の回転を、それぞれ 1 個ずつ右に回す」ことは、「全部の回転リスト自体を 1 個右に回す」ことと同じ、ということです。回転のリストは巡回的に閉じているので、こういう入れ替えができます。

証明は次のように rots の定義を展開して、rrot · lrot = id を使い、また rrot xs = lrotⁿ⁻¹ xs を使って書き換えていくだけです。

map rrot (rots xs) = {rots の定義} map rrot [xs, lrot xs, lrot² xs, …, lrotⁿ⁻¹ xs] = {map の定義、および rrot · lrot = id を用いて} [rrot xs, xs, lrot xs, …, lrotⁿ⁻² xs] = {rrot の定義} rrot [xs, lrot xs, …, lrotⁿ⁻² xs, rrot xs] = {rrot xs = lrotⁿ⁻¹ xs より} rrot [xs, lrot xs, …, lrotⁿ⁻² xs, lrotⁿ⁻¹ xs] = {rots の定義} rrot (rots xs)

役に立つ 3 つの恒等式

takeCols (j+1) · map rrot = consCol · fork (map last, takeCols j)(13.3)

ここで fork (f, g) x = (f x, g x) です。左辺は「各行を右回転してから最初の j+1 列を取る」、右辺は「最終列と最初の j 列を取り、最終列を先頭にくっつける」。同じ結果になります。

takeCols (j+1) · hdsort = hdsort · takeCols (j+1)(13.4)

第 1 列で並べ替えるのと、正の数の列を取るのは、順番を入れ替えてもよい、ということです。

hdsort · map rrot · sort · rots = sort · rots(13.5)

大事な恒等式 (13.5) 「ソート済み回転行列に対して、最終列を先頭に移して第 1 列で安定ソートし直す」と、元の行列に戻る、ということです。ただしこれは hdsort が安定ソート（同点の場合、元の順序を保つ）である場合にのみ成立します。ここが BWT の逆変換が成り立つ根本的な理由です。

この (13.5) の背景には基数ソートの性質があります。n × n 行列に対して、

sort = (hdsort · map rrot)ⁿ

すなわち「最終列を先頭に回して第 1 列で安定ソートする」を n 回繰り返せば全体のソートと同じ結果になります。これは基数ソートそのものです。

recreate の構成

recreate (j+1) · map last · sr = {仕様 (13.1)} takeCols (j+1) · sr = {性質 (13.5)} takeCols (j+1) · hdsort · map rrot · sr = {性質 (13.4)} hdsort · takeCols (j+1) · map rrot · sr = {性質 (13.3)} hdsort · consCol · fork (map last, takeCols j) · sr = {仕様 (13.1)} hdsort · consCol · fork (map last, recreate j · map last) · sr = {fork (f · h, g · h) = fork (f, g) · h より} hdsort · consCol · fork (id, recreate j) · map last · sr

recreate (j+1) = hdsort · consCol · fork (id, recreate j)

Dart // recreate（素朴版）: 最終列 ys から j 列ぶんの部分行列を再構築する。 List<List<T>> recreateNaive<T extends Comparable<Object?>>(int j, List<T> ys) { if (j == 0) return recreate0(ys); final prev = recreateNaive(j - 1, ys); // recreate (j−1) ys return hdsort(consCol((ys, prev))); // 先頭に ys を挿し込み再ソート }

より高速なアルゴリズム

今のままでは遅すぎる

とはいえ、行列全体をまるごと復元してから 1 行を選ぶのはあまりに無駄です。

目標は untransform を Θ(n log n) ステップで計算することです。

hdsort は 1 回で十分

まず、単集合（要素 1 個のリスト）だけの行列に対しては、先頭ソートは普通のソートと同じです。

hdsort · map wrap = map wrap · sort

次に、各段のソートは毎回同じ置換になることに気付きます。ソートに対応する置換 p を最初に 1 回だけ求めれば、それを毎回使い回せます。

sort ys = apply p ys p = map snd (sort (zip ys [0 .. n−1])) apply p xs = [xs !! (p !! i) | i ← [0 .. n−1]]

Dart // 「ys を安定ソートする置換 p」を 1 回だけ求めておく。 List<int> permOf<T extends Comparable<Object?>>(List<T> ys) { final n = ys.length; final ix = [for (var i = 0; i < n; i++) i]; ix.sort((a, b) { final c = ys[a].compareTo(ys[b]); return c != 0 ? c : a.compareTo(b); // 安定性のため同点は元順 }); return ix; } // apply p xs: p の指示どおりに xs を並べ替える。 List<T> apply<T>(List<int> p, List<T> xs) => [for (var i = 0; i < p.length; i++) xs[p[i]]];

apply p · consCol = consCol · pair (apply p)(13.6)

recreate 0 = map (const [ ]) recreate (j+1) = consCol · fork (apply p, apply p · recreate j)

Dart // hdsort の代わりに、事前計算した置換 p を毎回 apply するだけで済む。 // apply p は「行列そのもの（＝行のリスト）」を並べ替える。 List<List<T>> recreateApply<T extends Comparable<Object?>>( int j, List<T> ys, List<int> p) { if (j == 0) return recreate0(ys); final prev = recreateApply(j - 1, ys, p); return consCol((apply(p, ys), apply(p, prev))); }

これで hdsort の繰り返しが apply p の繰り返しに置き換わりました。

再帰を「解く」

次に、この再帰を閉じた形の式に直します。次が解になります。

recreate j = tp · take j · tail · iterate (apply p)(13.7)

Dart // 「apply p を繰り返し掛けた列」を先頭 j 本取り、転置すれば復元行列になる。 List<List<T>> recreateClosed<T extends Comparable<Object?>>( int j, List<T> ys, List<int> p) { final cols = <List<T>>[]; var cur = ys; for (var i = 0; i < j; i++) { cur = apply(p, cur); // tail · iterate なので 1 段ぶん先から取る cols.add(cur); } // 転置: cols は j 本の列 → 行列（各行の長さ j）に組み直す。 return [ for (var r = 0; r < ys.length; r++) [for (var c = 0; c < j; c++) cols[c][r]] ]; }

ここで tp は行列を転置する Haskell の transpose の略記です。tp について次の 3 つの性質を使います。

tp · take 0 = map (const [ ])(13.8) tp · take (j+1) = consCol · fork (head, tp · take j · tail)(13.9) apply p · tp = tp · map (apply p)(13.10)

(13.10) を平易に言い換えると、「行に置換を掛けたければ、転置して列に掛けて戻せばよい」ということです。

j = 0 では (13.8) から直ちに従います。帰納ステップでは、(13.10) と map と iterate の基本法則、そして (13.9) を順に使って書き換えていくと、次が得られます。

recreate (j+1) = tp · take (j+1) · tail · iterate (apply p)

最後の詰め

untransform (ys, k) = (recreate n ys) !! k でしたから、k 番目を取る操作と (13.7) の右辺を組み合わせて

(!!k) · recreate n = {(13.7)} (!!k) · tp · take n · tail · iterate (apply p) = {(!!k) · tp = map (!!k) より} map (!!k) · take n · tail · iterate (apply p) = {map と take, tail は入れ替えられる} take n · tail · map (!!k) · iterate (apply p)

ここで iterate の便利な法則があります。f x ⊕ y = x ⊕ g y が成り立つとき、

map (⊕y) (iterate f x) = map (x⊕) (iterate g y)

この法則と (apply p ys) !! k = ys !! (p !! k) を組み合わせると、

map (!!k) (iterate (apply p ys)) = map (ys!!) (iterate (p!!) k)

take (length ys) (tail (map (ys!!) (iterate (p!!) k)))

Dart // リストのままの untransform: 添字 k から p[k], p[p[k]], … を辿り // その各段で ys[…] を読み取れば元の列が復元される（tail で 1 段飛ばす）。 List<T> untransformList<T extends Comparable<Object?>>(List<T> ys, int k) { final n = ys.length; final p = permOf(ys); final out = <T>[]; var j = k; for (var i = 0; i < n; i++) { j = p[j]; // tail 相当（0 段目は捨てる） out.add(ys[j]); } return out; }

なぜ配列を使うのかこの式のままだと、リスト索引付け (!!) がリストの長さに比例する時間かかるため線形時間になりません。そこで Haskell の Data.Array の配列に置き換え、定数時間の (!) を使います。

untransform (ys, k) = take n (tail (map (ya!) (iterate (pa!) k))) where n = length ys ya = listArray (0, n−1) ys pa = listArray (0, n−1) (map snd (sort (zip ys [0..])))

Dart // 配列版 untransform: Dart の List は既に O(1) 添字なので、 // あとは p を 1 回計算しておけば全体 Θ(n log n) で終わる。 List<T> untransform<T extends Comparable<Object?>>((List<T>, int) input) { final ys = input.$1; final k = input.$2; final n = ys.length; final pa = permOf(ys); // Θ(n log n) final out = List<T>.filled(n, ys.first); var j = k; for (var i = 0; i < n; i++) { j = pa[j]; out[i] = ys[j]; } return out; }

pa の計算は Θ(n log n) ステップ、残りは Θ(n) ステップで済みます。全体で Θ(n log n) です。

transform の再考

アイデア：回転のソートを「接尾辞のソート」に置き換える

元の transform も改良できます。カギは、回転をソートするのを、接尾辞をソートすることに置き換えられるという点です。

tag xs = xs ++ [eof] と置きます。eof（end of file）は xs に絶対現れない番兵記号です（文字列ならヌル文字、自然数のリストなら −1 など）。

なぜこれで OK? eof は他のどの要素とも異なるので、tag xs の先頭 n 個の接尾辞をソートする置換は、xs の n 個の回転をソートする置換とちょうど一致します。長さの異なる接尾辞どうしを比べても、番兵で必ず順序が決まるからです。

sort (rots xs) = apply p (rots xs) p = map snd (sort (zip (tails (tag xs)) [0 .. n−1]))

Dart // 接尾辞ソートによる置換 p の計算。 // eof は入力に絶対現れない番兵で、他のどの要素よりも「大きい」値にする // （そうすると短い接尾辞が同点比較で後ろに回り、回転をソートしたのと同順序になる）。 List<int> suffixPerm(List<int> xs, int eof) { final n = xs.length; final tagged = [...xs, eof]; final ix = [for (var i = 0; i < n; i++) i]; // tag xs の先頭 n 個の接尾辞を辞書順にソートして得た元位置の並びが p。 ix.sort((a, b) { var i = a, j = b; while (true) { final va = tagged[i], vb = tagged[j]; if (va != vb) return va.compareTo(vb); i++; j++; } }); return ix; }

と書けます（tails は接尾辞を長さの降順に返す標準関数）。あとは同様の書き換えで、

map last · sort · rots = apply p · rrot

transform xs = ([xa ! (pa ! i) | i ← [0 · · n−1]], k) where n = length ys k = length (takeWhile (≠ 0) ps) xa = listArray (0, n−1) (rrot xs) pa = listArray (0, n−1) ps ps = map snd (sort (zip (tails (tag xs)) [0 .. n−1]))

Dart // 効率版 transform: 接尾辞ソートによる置換 ps を使って // (rrot xs)[ps[i]] を並べたものが最終列。k は ps 内の 0 の位置。 (List<int>, int) transformFast(List<int> xs, {int eof = 0x7fffffff}) { final n = xs.length; final ps = suffixPerm(xs, eof); final xa = rrot(xs); // 右回転しておくと最終列は xa[ps[i]] で取れる final lastCol = [for (var i = 0; i < n; i++) xa[ps[i]]]; var k = 0; while (k < n && ps[k] != 0) { k++; } return (lastCol, k); }

ボトルネックは ps の計算で、それ以外は Θ(n) です。前の pearl で見たように、長さ n のリストの接尾辞のソートは Θ(n log n) でできます。さらに、有限アルファベット上のリストなら接尾辞木を作ることで線形時間まで下げられます（Gusfield (1997) を参照）。

Dart // 動作確認: 'yokohama' を int リストとして BWT にかけ、逆変換で復元する。 void main() { final xs = 'yokohama'.codeUnits; // 効率版で順変換 → 逆変換 final (ys, k) = transformFast(xs); print(String.fromCharCodes(ys)); // hmooakya print(k); // 7 final restored = untransform((ys, k)); print(String.fromCharCodes(restored)); // yokohama }

最後に

BWT は Burrows and Wheeler (1994) の技術報告で発表されましたが、実際のアルゴリズムは 1983 年に Wheeler によって発見されていました。Nelson (1996) は BWT を世に紹介した記事で、当時の商用圧縮プログラムを凌駕できることを示しました。今では BWT は高性能圧縮ユーティリティ bzip2（www.bzip.org）に組み込まれています。基数ソートについては Gibbons (1999) が、ポイントレス計算による導出を扱っています。

まとめ

BWT は「回転をソートして最終列と行番号を返す」変換で、元に戻せる並べ替えを行う。
逆変換の本質は基数ソート性 (13.5)：最終列を先頭に回して安定ソートすると元の並びに戻る。
素朴には Ω(n² log n) だが、置換 p の 1 回計算と再帰の「閉じた形」への変形、配列の導入によって Θ(n log n) にできる。
順変換も接尾辞ソートで書き直せ、接尾辞木を使えば線形時間になる。

参考文献

Burrows, M. and Wheeler, D. J. (1994). A block-sorting lossless data compression algorithm. Research report 124, Digital Systems Research Center, Palo Alto, USA.

Gibbons, J. (1999). A pointless derivation of radix sort. Journal of Functional Programming 9 (3) 339–46.

Gusfield, D. (1997). Algorithms on Strings, Trees and Sequences. Cambridge University Press, Cambridge, UK.

Nelson, M. (1996). Data compression with the Burrows–Wheeler transform. Dr. Dobb’s Journal, September.

1 算術符号化については、後続の 2 つの pearl（Pearl 24 と Pearl 25）で扱われる。

2 ここでは transform [ ] = ([ ], 0) となる事実は無視し、transform の引数としては空でないリストのみを扱う。

どんな話？