Burrows-Wheeler変換の線形時間アルゴリズム

研究紹介です。今夏のSPIRE 2009という学会で

"A Linear-Time Burrows-Wheeler Transform using Induced Sorting", D. Okanohara, K. Sadakane, SPIRE 2009 pdf(draft)

というのを発表します。これは与えられた文字列に対し接尾辞配列を経ないでBurrows-Wheeler変換を直接行うというもので、アルファベットサイズによらず入力長に対して線形時間で行えます。基本的なアイディアは昨年のInduced Sortingによる接尾辞配列の線形時間構築アルゴリズム（いわゆるSAIS）を接尾辞配列を使わないでシミュレートするものです。pushとpop操作だけからなり、そのまま外部記憶上での構築とかにも対応できるようになっています。

Burrows-Wheeler変換（BWT, Block Sortingとも呼ばれるen-wikipedia）は、94年（発見自体は80年代）にBurrows, Wheelerらによって提案された文字列の可逆変換であり、データ圧縮、全文検索、データマイニングなど広い分野で使われる変換手法です。なにか大規模な文字列処理をするなら、とりあえずBWTをかけとけば大抵楽に処理できます。

--
今回の論文の位置づけを示すために、文字列処理用のデータ構造の歴史的背景とかを書いてみます。

まず、接尾辞木（Suffix Tree, ja-wikipedia）が1973年にWeinerによって提唱されます。接尾辞木は次のように作られます。まず長さnの入力文字列T[1, n]に対し、i文字目からj文字目までの部分文字列をT[i, j]で表すとします。この時、Si = T[i, n]をTのi番目の接尾辞と呼びます。入力Tには全部でn個の接尾辞があるのですが、この接尾辞それぞれをキーとしたtrieを作ります。次に、子が一つしかない節点を取り除き、更に文字列を枝上でそのまま持つ代わりに、葉に何番目の接尾辞に対応するかの情報を格納すると共に、各節点には深さ情報のみを記録するとします。するとT中に出現する全ての部分文字列はO(n^2)個あるのですが、節点の数は高々n個で、さらに葉の数もn個しかないようなコンパクトなデータ構造ができます。これが接尾辞木です。接尾辞木を利用することで非常に多くの様々な文字列処理を高速に行えることが知られています
（suffix treeを利用した様々なアルゴリズムについては有名なGusfield本にいろいろ書いてありますAlgorithms on Strings, Trees and Sequences: Computer Science and Computational Biology）

こうしたデータ構造は、大量のデータ（30億文字かならるヒトゲノムとか）を扱うことが主な目的なので、高速に構築できることが重要になります。接尾辞木も提案されてから様々な手法が提案されてきました。様々な改良が重ねられた後に、95年 Ukkonenによってアルファベットサイズが定数の場合の現実的な線形時間アルゴリズムが、97年にFarachによって任意のアルファベットサイズに対する線形時間アルゴリズムが提案されています。

接尾辞木は非常に強力なツールなのですが大きな問題がありました。それは接尾辞木を保存するのに必要な作業領域量、つまりメモリを非常にたくさん必要とすることです。入力サイズに対し（ポインタのサイズが定数の場合）線形のサイズとはいえ効率的な実装でも入力長に対し数十倍のメモリが必要とされてました。

こうした問題を解決したのが91年にManberらによって提案された接尾辞配列（Suffix Array)です。これは接尾辞木の葉の部分のみを使ったもので作業領域量を入力長*ポインタサイズ（普通は4バイト）と大幅に減らすことに成功しました。接尾辞配列を正確にいうと、入力Tに対し、接尾辞配列SA[1, n]をT[SA[i], n] < T[SA[i+1], n]が全ての1 <= i < n-1に対し成り立つような配列と定義されます（文字列間の<は辞書式順序とする）。SAを使えば、接尾辞が辞書式順序に並んでいるので、ある文字列qを検索したい場合にはT[SA[l], n] <= q <= T[SA[r], n]が成り立つようなl, rを二分探索で求めることでqの出現場所を全て発見することができ、SA[l, r]がqの出現場所であり、r-l+1がqの出現回数となります。任意の文字列に対しO(log n)時間で検索できるデータ構造が約5n byte（入力文字列n + 接尾辞配列4nバイト）で保持できるということになります。

大幅に作業領域量を減らすことに成功した接尾辞配列でしたが、構築時間が大きくなるという問題がありました。接尾辞木を最初に作り葉の部分だけを回収すれば線形時間で構築できるのですが、結局作業領域量が大きくかかってしまいます。接尾辞木を経ないで接尾辞配列を直接構築するには、接尾辞を辞書式順序にソートする必要があり、最悪計算量がO(n^2 log n)かかるという問題がありました（入力にT = aaaaaa...a$というものを考えた場合、毎回の比較にO(n), 比較回数にO(n log n)かかるため）。この構築時間は97年に定兼先生によってO(n log n)に改善され、2003年には接尾辞配列を直接線形時間で構築するアルゴリズムが3つ同時に提案されます。その後もさらに高速な構築アルゴリズムが提案されていき、現在の構築時間では1秒あたり4～8MB程度構築できるようになっています（ベンチマーク@white page）

接尾辞配列と大体同時期にBurrows Wheeler変換というものが提案されました。これは入力T[1, n]をB[i] := T[SA[i]-1]（SA[i]=0の時はB[i]=T[n])で定義される文字列B[1, n]に変換するというものです。
この変換は一見すると直感的ではありませんが、非常に重要な性質を持っていることが分かっています。代表的なものを三つあげると、

(1)可逆変換である

BW変換後のB[1, n]のみから入力文字列T[1, n]に戻すことがO(n)時間で可能です。これはLF-mappingと呼ばれるものを利用します。B[i]=cに対応する接尾辞はrank(B, i, c) + C(c)番目にあるというものです、但しrank(B, i, c)はB[1, i]中のcの出現回数、C(c)はB中のcより小さい文字の出現回数です（つまりT中と言い換えても同じ）。（実際の逆変換のコードはC++20行弱程度で例えばpdfの20枚目あたりにあります）

(2)元の入力が圧縮しやすい場合、BW変換後は同じ文字が連続して出現しやすい

Bは定義より、各文字をそれに後続する文字列をキーとしてソートした結果となっています。よって似た文脈の直前に同じ文字が出現しやすいのであれば、Bでは連続して同じ文字が出現することになります。例えば英語では"the"が頻出するので"he .."の文脈の直前には"t"が出現しやすくB中でも"t"が連続して出現するようになります。

この直感は様々な形で解析されています。まず圧縮しやすいというのは(k-1)文字を見たら次の文字は非常に予測しやすいということでいえて、これはHk（k次経験エントロピー）という数値で測ることができます。そして、BW変換後は文脈情報を全く無視して先頭移動法などローカルな情報だけを使って符号化してもその後のサイズがHkに漸近することが知られています（"A simpler analysis of Burrows-Wheeler-based compression, H. Kaplan., and et. al.pdf）

(3)（圧縮）全文検索に利用できる

元々データ圧縮用に作られたBW変換は、2001年にFerraginaらによるFM-indexと呼ばれるデータ構造によって全文検索に利用できるということが分かりました。BW変換の定義には接尾辞配列が使われているので接尾辞の情報が入っているのですが実際にこれを、LF-mappingの性質を使うことで全文検索ができることが示されました。
このFM-indexではLF-mapping、つまりはrank(B, i, c)を求める必要があり、FM-indexが最初に提案された時には、この部分の効率的かつ高速な計算の実現が困難でした。それが03年にGrossiらによって提案されたwavelet treeを使うことで高速に行うことができ、現在ではFM-indexはテキスト長に依存せずクエリ長に比例する時間で全文検索が行えます。しかも必要な作業領域量はBと同じだけのスペースなので約nバイトで済みます。

詳細を追うことはやめておきますが、最長共通接頭辞配列や、木を表すテーブルと組み合わせることで接尾辞配列（拡張接尾辞配列、Enhanced Suffix Arrays）やBWTは接尾辞木と同じように複雑な文字列検索を行えることがわかっています。

ここまでの流れを追うと、文字列処理用のデータ構造は、効率的に操作できる能力を変えずにどんどんコンパクトになってきていました。

- 接尾辞木　1973年 O(n log n)、元文書の約20倍でヒトゲノムの場合は約50GB
- 接尾辞配列　1991年 n log n 、元文書の約5倍でヒトゲノムの場合は約13GB
- BW変換 1994年(2001年）、nHk、元文書の等倍以下でヒトゲノムの場合は約750MB

また、それらのデータ構造を直接高速に線形時間で構築するようなアルゴリズムも提案されてきました。今回の発表はこの流れに沿うものになってます。

線形時間の構築
- 接尾辞木 Ukkonen 1995年、Farach 1997年
- 接尾辞配列 Aluruなど 2003年
- BW変換 今回の 2009年

さて、今後もこの流れは続くかですが、コンパクトになっていく流れは、このままでは難しいといえます。すでにBW変換は元の文書を圧縮したサイズに近いサイズまできています。全文検索において、任意の部分文字列に対し正しい結果が得られるということは、元の文書を復元できるということを意味しているので、データ圧縮とやっていることは同じです。ですので、元の文書が持っている情報を捨てない限りはこれ以上小さくするのはかなり難しいと思われます（もしできたらデータ圧縮においても革命）。結果の正確さを100%ではなく、何らかの形で保障しながらサイズを劇的に減らす方向については十分ありうると思います。

それに対し、構築時間については、理論上は線形時間よりは速くするのはちょっと大変そうですが、実際には改善が続けられています。今回提案した手法も実装に落すところでは更に工夫が必要です（その実験/論文はそのうち）。実装以外にも例えば接尾辞配列において並列に計算するような方法や、検索時間を犠牲にして構築時間を減らすような方法もでています。秒間数MBを処理できるとは数年前から比べたら各段の進歩ですが、今流行りのマイクロブログのリアルタイム検索などを実現しようと考えたり、もうすぐ実用化されると噂される次々世代シーケンサは秒間1GB超という速さでゲノム配列を読めるということで、まだまだ改良しなければならなさそうです。