ハクビシンにもわかる全文検索

高速な全文検索アルゴリズムであるFM-indexについて解説する。理解しがたい点や間違っている点があれば是非コメントで指摘してほしい。

概要

FM-indexはリニアな文字列に対して検索をするアルゴリズムで、主に簡潔データ構造とBWT(およびLF mapping)という二つのアイデアから成り立っている。BWTはBurrows-Wheeler変換のことで、文字列を特殊な並び順に変換するという可逆関数である。BWTされた文字列を簡潔データ構造固有の操作をすることで、クエリ文字列の長さに比例した短い時間で文字列を探し出すのがFM-indexだ。

簡潔データ構造

簡潔データ構造に関してはFM-indexで必要となる二つの関数だけ説明して、詳細は次の機会に譲るとする。さて、二つの関数はともに文字列のある位置より前の部分に含まれている文字の数を数え上げるというものでrank()とrankLessThan()だ。

rank関数は二つの引数を渡す。位置と文字コードだ。指定した位置より手前の、指定した文字コードの数を返す関数だ。例えば "abracadabra" という文字列に対して rank(5, "a") とすると、"abrac" に含まれる "a" の数なので 2 が帰る。

rankLessThan関数はrank関数同様に二つの引数、位置と文字コードを渡す。指定した位置より手前にある、指定した文字コードよりも小さい文字コードの数を返す。例えば "abracadabra" に対してrankLessThan(5, "c") とすると、"a" "b" "a" がヒットするので 3 が帰る。

なんと検索を実現するのに必要な操作はこの文字を数え上げる二つの関数だけでできてしまうのだ。その秘密は次に説明するBWTの特殊な構造にある。

BWT

BWT, あるいはブロックソートと呼ばれる技法は、文字列を名前の通りソートして並び替える技法だ。まずは名前にも出てくるブロックを作り上げる。ブロックは文字列を一文字ずつシフトしてゲームの横スクロールみたく操作して作る。ここでは "abracadabra" という文字列を元に説明するが、末尾にマーカーとして "$" という文字を加える。このマーカーはほかの文字よりも小さいコードであることが重要だ。

abracadabra$
bracadabra$a
racadabra$ab
acadabra$abr
cadabra$abra
adabra$abrac
dabra$abraca
abra$abracad
bra$abracada
ra$abracadab
a$abracadabr
$abracadabra

マーカーを合わせて12文字 x 12行のブロックができあがった。次はこれを辞書順にソートする。

$abracadabra
a$abracadabr
abra$abracad
abracadabra$
acadabra$abr
adabra$abrac
bra$abracada
bracadabra$a
cadabra$abra
dabra$abraca
ra$abracadab
racadabra$ab

このブロックの最後の列を順番に並べると、"ard$rcaaaabb" となるが、これが目的の BWT された文字列となる。aが連続して並んでいる通り、BWT された文字列は仕組み上同じ文字が連続しやすい。そのため、bzip2 などの圧縮ソフトではBWTを前段の処理として使われる。今回はその特性に関しては関係が無いので置いておく。

圧縮の目的とは違うもう一つの性質が全文検索においては重要になってくる。つまり、ブロックの前方は辞書順にソートされているということだ。

BWT逆変換

BWTは可逆変換である。ではソートしてしまった文字列がどうやれば元に戻るのかだ。BWTされた文字列はブロックのうち末尾の列を抜き出したものだがこれを"L"と呼ぶ。これに対してブロックの先頭の文字は、BWT文字列をソートすれば得ることができる。これを"F"と呼ぶ。たとえば、abracadabra$の行はソート済みブロックの上から四つ目にあたるが、これは、Lが元文字列の末尾"$"であり、Fが元文字列の先頭"a"を指す。

F L
---
$ a
a r
a d
a $
a r
a c
b a
b a
c a
d a
r b
r b

また、ブロックの文字列はシフトしながら作られているためFの文字は、Lの文字の一つ先の文字である。"a$" "ra" "da" "$a" "ra" ca" "ab" "ab "ac" "ad" "br" "br" をいい感じにつなげれば元の文字列が完成する。

ではどうつなげればいい感じになるのか。少なくとも同じ文字が存在しない文字同士のつなげ方は簡単だ。一行目のFの$と、四行目のLの$は同じモノを指している。問題となるのが、同じ文字がいくつも存在している、例えば"a"のような文字の場合どうすれば良いか。答えを言うと、上から数えていくつ目かで対応できる。例えば、一番上の"a"同士がつながっている。一行目のLの"a"と二行目のFの"a"だ。七行目のLの"a"と三行目のFの"a"もつながっている。

なぜか？ブロックの先頭は綺麗にソートされているが、ほかの列は別の法則に従って並んでいるわけではない。F列あるいはその次以後の列がソートされたのに合わせてシフトされた文字列が並んでいるのである。一番上のF列の"a"と一番上以外のL列の"a"がつながったりなどはしない。そんなことをすればねじ曲がってしまうからだ。ここらへんは僕も明確な理解ができているわけではないので、きっちりした説明ができないのがもどかしいところだ。

F L
---
$ a   a1 -> $
a r   r1 -> a1
a d   d  -> a2
a $   $  -> a3
a r   r2 -> a4
a c   c  -> a5
b a   a2 -> b1
b a   a3 -> b2
c a   a4 -> c
d a   a5 -> d
r b   b1 -> r1
r b   b2 -> r2

LF mapping

前述のやり方には一つ問題がある。ソートなり、上から文字を探すなりが重い操作であるという点だ。そこで、BWTのソートされたブロックという性質を利用して高速にLとFをたどる方法を紹介しよう。それがFM-indexでとても重要となる概念であるLF mappingだ。

前述の逆変換とは逆に、後ろからたどるやり方となる。一行目のF($)からL(a)を求めるのは単なる配列アクセスなので高速だ。元文字列の最後尾は、この"a"だ。次にこのL(a)に対応するのは先ほどの考え方と同じく、一つ目の"a"なので、

二行目の"a$abracadabr"となる。F(a)に対してL(r)で、最後尾から二文字目は"r"だ。同様にたどっていく。

さて、F -> L は配列アクセスという定数時間で実現できるとして、L -> F をどうやって実現するかだが、ソートされたF列に対してある二つの性質が見えてくるのではないだろうか。

一つ目。"b"を例に取ると、一番上のbの行番号は、"$"と"a"の数に等しいという点だ。つまりある文字群の一番上の行は、その文字よりも小さい文字の数と一致するということ。

二つ目は、L列における上から数えていくつ目かという問題は、BWT文字列においてその文字より手前にある同じ文字の数に等しいということだ。例えば、8行目の"a"の上には2つの"a"があるが、それは言い換えるとBWT文字列の8文字目よりも手前にある"a"の数なのだ。

iを現在注目している位置、cを現在見ている文字、全体の文字列長をlenとすると、以下の計算式で L -> F が求まる。

rank(i, c) + rankLessThan(len, c)

説明は省いたが、ある簡潔データ構造では、rankやrankLessThanは一定の短い時間で求まるため、L -> Fも高速に算出することができる。実際には、rankLessThanは、文字列全体での数値なので、事前に求めることもできる。

LF mappingを使った検索

さてLF mappingはBWTの逆変換だけに使われるアイデアではない。ソートされたブロックとLF mappingによって、文字列検索も可能なのだ。それが今回の話のキモになるFM-indexである。このアルゴリズムは backward とも言われ、クエリ文字列の後ろから順にたどるものだ。LF mappingまでを理解できた読者であればもう一歩だ。

まずは"a"一文字の検索を行おう。

考えるべきこととしてFが"a"である行は並んでいる。その範囲を見つけることが1文字の検索と同義となる。"a"の始まりの行は、"a"よりも小さい文字の総数に等しい。"a"の終わりの次の行は、"a"よりも小さい文字の総数とBWT文字列に出てくる"a"の数に等しい。以下のコードとなる。

begin = rank(0, "a") + rankLessThan(len, "a")
end = rank(len, "a") + rankLessThan(len, "a")

注意しなければならないのは、実際に文字列検索として使うためには、この「範囲」から、元文字列の位置に変換する必要があるということだ。LF mappingを繰り返して先頭にたどり着くまでの操作回数(=位置)を算出するという方法があるが、当然これは文字列長が長くなればなるほど遅くなってしまう。手としては、一定間隔で位置情報を配列に持たせることだ。位置情報を持った行にたどり着いたら配列を参照することで計算のショートカットができる。

$abracadabra
a$abracadabr < begin
abra$abracad
abracadabra$
acadabra$abr
adabra$abrac
bra$abracada < end
bracadabra$a
cadabra$abra
dabra$abraca
ra$abracadab
racadabra$ab

複数文字の検索

複数文字の検索もアルゴリズム的にはそう変わらない。

"ab"で検索するのを例に取ると、まず全体に対して"b"だけで前述の方法で検索するとbegin, endが求まる。次にbegin, endの範囲に対して"a"で検索を行う。

begin = rank(begin, "a") + rankLessThan(len, "a")
end = rank(end, "a") + rankLessThan(len, "a")

さて、なぜこれで二文字の検索ができるのか。二つの性質を思い出してほしい。まずはabを探すというのはabから始まる行を探すということ。もう一つはFの一つ前の文字がLであるということだ。つまり、"ab"をさがすためには"b"で始まり"a"で終わる行を探せば良いということだ。

一文字目の"b"を探した段階でbegin = 6, end = 8となっている。6行目よりも以前にL列にaが登場するのは1つ、8行目よりも以前は3つとなる。"a"よりも小さい文字である"$"は1つなので、begin = 2, end = 4 となる。

$abracadabra
a$abracadabr
abra$abracad < 二文字目 begin
abracadabra$
acadabra$abr < 二文字目 end
adabra$abrac
bra$abracada < 一文字目 begin
bracadabra$a
cadabra$abra < 一文字目 end
dabra$abraca
ra$abracadab
racadabra$ab

後はこれを繰り返すだけで複数文字の検索は可能である。また、begin >= end となった時点で、対象の文字列は存在しないことが確定する。

最後に

• BWTの生成は実際にはもっと高速なやり方で作られる。興味のある人は Induced Sort を調べると良いだろう
• 今回説明を省略した簡潔データ構造は、完備辞書とwavelet matrixの二つだ
• 僕が作成中の erukiti/cerebrums というアプリで FM-Index を CoffeeScript で実装したものを https://github.com/erukiti/cerebrums/blob/master/src/fm_index.coffee として公開している

ふとしたことから知ったFM-indexを実現するために、完備辞書・wavelet matrix・BWT・FM-index と順に実装してきて、アルゴリズムとはかくも楽しいモノだったのかと感動したので、是非ともその感動を分かち合いたく、この記事を書いた。