Blog
最近オープンウォーターダイバーのライセンスを取りました。徳永です。
今日はBurrows Wheeler Transform（BW変換もしくはBWT）の逆変換において用いられるLF mappingを説明します。
BWTはデータ圧縮の前処理などに使われるテクニックです。Burrows Wheeler Transformはとても簡単でわかりやすい（高速な実装は複雑ですが……）のですが、逆変換で用いられるLF mappingは、実装は簡単なものの、なぜそれでよいのかは少しわかりにくいところがあります。また、私はこれまで、LF mappingがなぜあれでうまくいくのか、わかりやすい説明を日本語でも英語でも見た記憶がありません。そこで今回はLF mappingを中心に説明します。なお余談ですが、BTWのMichael Burrowsは現在はGoogle勤務で、ChubbyやBigTableなどのソフトウェアに関わっているそうです。
Burrows Wheeler Transformとは
まずはBWT自体について簡単に説明しておきましょう。BWTとは文字列に対する一種の変換です。逆変換が存在し、変換後のデータに逆変換の操作を施すと、変換前のデータを得ることができます。いろいろ良い性質があるのですがそれについてはまた後ほど書きます。
今回の説明では、サンプル文字列としてabracaの6文字を使います。これは[1]で使われていたものです。
後段での説明の都合上、abracaの末尾に仮想的な文字$を追加します。これでabraca$になります。このabraca$を一文字ずつ左に循環シフトさせたものを縦に並べてみます。
abraca$ braca$a raca$ab aca$abr ca$abra a$abrac $abraca
次に、それぞれの行を辞書式順序で比較しながらソートします。ただし、$はどの文字と比較しても小さいものとします。結果は以下のようになります。
$abraca a$abrac abraca$ aca$abr braca$a ca$abra raca$ab
一番左の列をF、一番右の列をLと呼ぶことにします。ソート後の各行の一番最後の文字を取ってくるとac$raabとなりますが、これが元の文字列abraca$のBWTの結果となります。各業の一番最後の文字を取ってくるという事は、Lの部分を取ってくるという事ですね。
BWTには
- 圧縮しやすい
- 元のテキストが復元できる
という2つの性質があります。圧縮しやすいのは、各行が辞書式順序でソートされており、BWTの各文字というのは辞書式順序でソートされる直前の文字なので、自然言語の性質を考えると圧縮されやすいのはまぁなんとなくそうだろうね、という感じです。自然言語は完全にランダムではないので、BWTをかけると、変換後のデータでは同じ文字が大量に並びやすくなります。
元のテキストの復元については、BWT後の文字列からさっきの行列が復元できれば後は自明であると言えます。行列の復元法については以下で説明します。
まず、BWT後の文字列をソートすると最初の列が復元できます。行列の各行は辞書式順序でソートされているから、一列目は辞書式順序でソートされていなければならないわけです。一番最後の列（L）はBWT変換後の文字列そのものですから、結局、一番最初の列と最後の列だけはわかっている、という状態になります。その状態の行列を以下に示します。?の部分はまだどんな文字になるのかわからない状態です。
$?????a a?????c a?????$ a?????r b?????a c?????a r?????b
ここで、最初に行列は元の文字列を循環シフトさせて作ったことを思い出しつつ、右に1つ循環シフトしてみます。
a$????? ca????? $a????? ra????? ab????? ac????? br?????
この行列は循環シフトを一つ進めただけなので、元の行列がソートされていない状態になってしまっただけだという事がわかります。（ここは直観的ではないので本当はもう少し丁寧に説明すべきなのですが、今回は省略します。Don’t think, feel!!）このため、これをソートすると、2列目が復元できます。以下はソート後の行列です。
$a????? a$????? ab????? ac????? br????? ca????? ra?????
辞書式順序でソートされていることから、最後の列はBWT後の文字列で置き換えられます。（2つ上の行列と比較してみて下さい。）最後の列をBTW後の文字列で埋めた行列は以下のようになります。
$a????a a$????c ab????$ ac????r br????a ca????a ra????b
以下同様に、循環シフトとソートを繰り返し、3列目、4列目、5列目…と行列を復元していきます。
原理的にはこのようにして復元していけるわけですが、スペースが文字列長の2乗分必要になりますし、ソートのコストもばかになりません。そこで、もっと高速に復元を行う手法としてLF mappingが出てきます。
LF mappingとは
LF mappingは、L[i]の（元の文字列での）1文字前の文字はL上でどこにあるかを返す、という操作です。以下の説明を読めばわかりますが、これは結局L上のある文字がF上でどこに出現するかを探してその添字を返す操作になっています。LとFを行き来する様からLF mappingと言います。（Lは一番右の列、Fは一番左の列のことでしたね。）
LF mappingのためには、C[ch]とocc(ch, beg, end)という2つの操作が必要になります。Cは配列で、C[ch]は、chという文字よりも小さな文字の出現回数を表すものとします。例としてC[b]の値を計算してみましょう。今回はbより小さな文字は$とaの2種類なので、C[b] = 1 + 3 = 4となります。occはchという文字がbegとendの間の範囲で何回出現するかで、これは簡潔データ構造の世界ではrankと呼ばれる操作で実現できます。occの効率的な実装はWavelet TreeやWavelet Matrixが提案されています。ここではなにか効率的な実装がすでにあるものと仮定します。
LF mappingはCとoccがあれば後は非常に簡単で、以下のように書けます。ただし、LはBWT後の文字列です。
LF(i) = C[L[i]] + occ(L[i], 0, i)
なぜこれで求まるのかは後で詳しく説明します。
LF mappingを用いた逆変換
LF mappingは、L[i]の（元の文字列での）1文字前はL上ではどこにあるかを返す、という操作ですから、一番最後の文字がどこにあるかを記録してあれば、あとはLF mappingを順繰りに適用していくことで、後ろから文字列を復元していくことができます。
LF mappingで使われるCとoccは定数時間で実行できるので、行列を2列目、3列目……と復元していくよりはだいぶ高速に復元することができます。結局、逆変換は元の文字列の長さに対し線形時間で実現できます。
特に難しいことはないのでコードは示しませんが、普通にfor文かなにかでループを回してLF mappingを実行していくだけです。
なぜこれでうまくいくのか
LF mappingではひとつ前の文字のL上におけるインデックス（添字）が求まるわけですが、なぜさっきの式で求められるのでしょうか。先ほどの行列を使って説明します。
$abraca a$abrac abraca$ aca$abr braca$a ca$abra raca$ab
まず、簡単な例として、L[1] = cの一つ前の文字がどこにあるかを求めてみましょう。上の行列において、一列目（Fの列）にはcはひとつしかなく、Fの行は辞書式順序でソートされているので、cが初めて出現する位置はC、つまり$とaとbの3つの文字の出現回数の和と同じになります。（+1しないでいいのは配列の添字が0から始まるからですね。）Cは5ですので、ca$abraの行が見つかります。なので、LF(1) = 5となります。
問題は、aのように複数回現れる文字です。複数の候補の中から適切な行を見つける必要があります。aは3回出現しますので、L上の3つのaとF上の3つのaの間の対応関係の情報をどこかに保存しておかなければならないように見えます。が、実はかんたんな方法で対応が取れるので、いちいち情報を保存しておく必要はありません。
ここからがこの記事で一番大事なポイントなので、目を皿にしてよく理解してください。
aの場合で説明を続けますが、F側で1番目に出てくるaとL側で1番目に出てくるa、F側で2番目に出てくるaとL側で2番目に出てくるa、… というふうに対応関係があります。例えば、L側がaである$abracaという行はF側がaであるa$abracと対応します。実際、$abracaを1つ右に循環シフトするとa$abracになりますね。つまり、L側の3つのaをa1, a2, a3として区別することにすると、F側での3つのaもa1, a2, a3の順に出現するのです。なぜこのように対応が取れるのでしょうか、その不思議をもう少し詳しく調べてみましょう。
F側でaが出る行は、当たり前ですが辞書式順序でソートされています。一方、L側でaが出てくる行は、循環シフトで一つ右に動かすと、F側でaが出てくるようになります。L側でaが出てくる3つの行を取り出して右に循環シフトしたものが下の行列です。
a$abrac abraca$ aca$abr
ソートはしていませんが、辞書式順序でソート済みの結果と同じ状態になっています。これは偶然ではなく、1列目は全部aで同じなのでソート順序に影響を与えず、2列目以降はもともとシフト前に辞書式順序でソートされていたわけですから、ソート済みなわけです。（ちょっとトートロジー的ですが。）
という訳で、L側でaが出てくる行は循環シフトで右に動かすとF側でaが出てくるようになり、しかもソート済みなわけなので、F側でi番目に出てくるaとL側でi番目に出てくるaは対応しているのです。という訳で、L側で何個目に出てくるaなのかがわかれば、それを足してやればF側での対応するaがどこにあるのかが求まります。つまり、まずC[a]で1つ目のaの出現位置を求め、さらに何個目のaなのかをocc(a, 0, i)で求めてたしてやれば、F上でどこにあるのかがわかるわけです。
上記の議論はaだけではなくすべての文字について適用できることは、一度理解できれば簡単にわかります。が、理解するのは最初は難しいと思いますので、何度か読んでみてください。
参考文献
2,3本の論文がここに書かれる予定でしたが、なくしてしまって見つからないので、見つかったら書きます。
