FM-index
Run-Length FM-Index
はじめに
本記事は データ構造とアルゴリズム Advent Calendar 2019 13日の金曜日の記事です.
前回(12 日目):Immutableなデータ構造について(@takoshi さん)
次回(14 日目):動的木上の最小シュタイナー木をtoptreeで解く(@niuez_さん)
FM-index は検索対象のテキストから予め索引を構築しておくことで,テキストに含まれる任意のパターン文字列の個数を数えるクエリを,テキストのサイズに依らず高速に実行できるデータ構造です.加えて,suffix array や inverse suffix array の一部を追加で保持しておくことで,パターンの位置の列挙 やテキストの復元 といったクエリを高速に実行することができます(自己索引).
主要な応用としてゲノム解析(例:HISAT2)などが挙げられます.身近なところでは,arXiv をコーパスとした高速な英文コロケーション検索を提供する Hyper Collocation でも用いられています(解説).
FM-Index に関しては,高速文字列解析の世界 や 簡潔データ構造 といった第一人者による書籍に紹介があり,またインターネット上にも id:echizen_tm さんのブログ を初めとする素晴らしい記事があります.なので日本語で読める情報には事欠かきません.しかし,FM-Index に連長圧縮を導入して圧縮効率の改善を図った run-length FM-index については,知る限りでは日本語の解説記事がありませんでした.本記事ではこのデータ構造について紹介します.
本記事で紹介するデータ構造は以下を基にしています:
Veli Mäkinen and Gonzalo Navarro. 2005. Succinct suffix arrays based on run-length encoding. Nordic J. of Computing 12, 1 (March 2005), 40-66.
高速文字列本に倣い,本記事における文字列および配列の添字は 0-based であるものとし,やを以下のように定義します:
:列のスライス中のの出現回数を返す
:列において番目に出現するの位置を返す
FM-Index
FM-index の詳細については,最初に紹介した文献,もしくは拙著のメモをご覧ください.
FM-index は以下2つのデータ構造から成ります.
: テキストを Burrows-Wheeler 変換した文字列について,以下のクエリをサポートするもの:
: のスライス中に出現する文字の個数を数える
: を得る
: テキストに現れる各アルファベットについて,テキスト中のよりも小さい文字の出現回数を与える
テキストのBurrows-Wheeler 変換は,の suffix arrayをもとに
とすれば求めることができますが,を生成せずに直接計算する手法も提案されています (Okanohara & Sadakane, 2009).
また,FM-index 中には直接保存しませんが,今後の議論のために配列をとして定義します.Suffix array の定義から,はを文字の大きさについて(安定)ソートしたものと一致します.下図はについてとを示したものです(太字がの suffix).
さて,の格納に wavelet tree もしくはその亜種,にアルファベットサイズ分の配列を用いることで,テキスト中のパターン文字列の出現回数を数えるクエリを,パターンのサイズとアルファベットのサイズについて時間で計算することができます.アルゴリズムを以下に示します:
let mut s = 0;
let mut e = n;
for i in (0..n).rev() {
let c = P[i];
s = rank(c, L, s) + C[c];
e = rank(c, L, e) + C[c];
if s >= e {
break;
}
}
return e - s
繰り返しごとに suffix の prefix としてパターンが現れる範囲が狭まっていくようなイメージです.
の計算に現れる式は, のとき上の LF-mappingと等しくなります.LF-mapping は におけるある文字がのどの位置にあるかを与える関数です.例えば上図においては及びにあるので,です.またその逆関数を本記事では inverse LF-mapping と呼びます.Inverse LF-mapping も FM-index によってとすることで求めることができます.
Run-Length FM-Index
概要
FM-index を構成するは,使用するデータ構造の種類によりますが,おおよそテキストの長さに比例する容量を必要とします.しかし Burrows-Wheeler 変換された文字列は比較的圧縮しやすいことが知られており,特にテキスト中に同じパターンが頻出しているとには同じ文字が連続して出現しやすくなります.そこで,を連長圧縮 (run-length encoding) することで索引サイズを削減しつつ,クエリのオーバーヘッドを定数倍に抑えたものが run-length FM-index です.
構成
中の連続する同一文字の列を run と呼びます.例えばのとき,はの5つの run からなります.
Run-length FM-index では,の代わりにの各 run の文字 (run head) ,および各 run の長さ (run length) を順番に並べたものを,それぞれととして保存します.例えば run head の列は の例だとになります.でもやといったクエリが必要になるため,それらをサポートする構造を用いて保存します.論文では各文字ごとに出現位置を記録したビットベクトルを用いていますが,FM-index におけると同様に wavelet tree 等を用いる方法についても言及されています.
一方, run length を並べたものは同じ例だと になります.これをそのまま格納する代わりに,各要素について1つの と個のをマップして連結した列をとします.例えばの例だとになります.ここで,の長さはの長さに等しく,に現れるの数は中の run の数に等しいのがポイントです.は簡潔ビットベクトルを用いて格納します.
さらに, run length を並べた配列を,対応する run head の文字の大小関係をキーとして安定ソートしたものを考えます.同様の例だとになります.これにと同様のマッピングを施し,として簡潔ビットベクトルを用いて格納します.例だとになります.
最後に,に現れる各アルファベットについて,中のより小さい文字の出現回数を数えたものをとします.例ではです.
FM-index で登場したの代わりにを保存したものが run-length FM-index です.
LF-Mapping
FM-index においての計算に用いた式 を,run-length FM-index のために変形します.
まず, がにおけるある run の先頭を指している状況を考えます(run の文字は任意).このとき,任意の文字について以下の式が成り立ちます:
右辺の部分式はそれぞれ以下に対応します:
: ある run (下図では) のにおける開始位置(run head の位置)
: その run の run head のにおける位置
: においてその run より前に出現する,でできた run の数
ここで左辺について,はにおいてある run より前に現れるの数,または(または)に現れるより小さい文字の出現回数を表しています.よってそれらの和は,その run の直前にあるの,における位置 + 1 を示しています(下図の☆部分).
一方において,☆で示した位置の手前にある1の個数を数えると,より小さい文字に対応する範囲には個,文字に対応する範囲には個(= ある run より前に出現するでできた run の数)の 1 が出現することがわかります.したがって,それらの和でを取ることで同じ☆の位置を得ることができます.
この結果を利用して,任意のとからが計算できます.の値による場合分けになります.
のとき:
のとき:
まず,を判定するためにの値を知る必要があります.ここでが成り立つので,の要素が取得できれば判定が可能です.
のときは,最初の状況(がある run の先頭位置)と同じに式になっています. において位置を通る run の先頭位置をとすると,, が成り立つので,
となります.
またのときは,最初の式が,をまたぐの run の前にあるのにおける位置 + 1 ,すなわちを通る run の run head (下図において上の位置にある)のにおける位置()を指しています.
を通る run の run head の位置はなので,を上式に足してやれば の LF-mapping になります.
Inverse LF-Mapping
inverse LF-mapping もから求めることができます.これはクエリの実行には必要ありませんが,自己索引として索引からテキストを正順に復元する際に便利です.以下は元論文の内容ではなくこの記事の著者による一つのアイディアです.もっとうまい方法があるかもしれません.
における LF-mapping と異なり,のケースを考慮する必要がないのでの場合のみを考えます.
任意のについて,
ただし.
においてをまたぐ run を考えます.上図において, ,またとにおける run の先頭位置をそれぞれとします.図のように,の中には 1 が個,うちより小さい文字に対応する 1 が個存在するので,その run のにおける位置は になります(上式の第一項).あとは run の先頭からの offset,すなわちを足せば終わりです.
の値は,FM-index では であるようなを二分探索することで求めることができました.Run-length FM-index においては,の単調性との関係から,を満たすようなを二分探索することで求めることができます.
実装と評価
Rust で FM-index と run-length FM-index を実装しました.
やの保存には wavelet matrix (Claude & Navarro, 2012),簡潔ビットベクトルは rsdic でも実装されている Navarro and Providel (2012) の手法を用いています.
ここまで空間計算量に関する議論を避けていました.結果が用いる下部構造に依存するのと,この記事の筆者がその手の計算に不慣れなためです.が連長圧縮された一方で,新しいデータが追加されているため常に圧縮が上がるとは限らなさそうです.そこで実際に索引を構築してサイズを比較してみます.
かなり artificial な例ですが,以下のような確率 またはで次の状態に遷移して状態ラベルの文字を吐く遷移系によって出力されたテキストを用います.を 1 に近づけることで同じ文字列が出現しやすくなります.
テキストを 10,000,000 文字とし, を変化させたときの FM-index と run-length FM-index のサイズは以下のようになりました:
p 索引サイズ (FM-index) (byte) 索引サイズ (run-length FM-index) (byte) run の数 / テキストの長さ
0.5 2182592 3007432 0.3359
0.9 1473872 1585144 0.1303
0.99 1010064 827208 0.0167
0.999 946088 650976 0.0017
Run-length FM-index が有利になるためには,繰り返して出現する文字列がそれなりに多い必要があるようです.
また,の実行速度を測るベンチマークは https://github.com/ajalab/fm-index/blob/master/benches/count.rs にあります.一度のイテレーションで 256通りのパターンを検索しているので,結果を 256 で割ったものをおおよその平均とみて良いでしょう.
筆者の環境では,run-length FM-index は FM-index に比べて 3 ~ 4 倍程度のオーバーヘッドがある様子が見られます.それでも1,000,000 文字のテキスト中に現れるパターンの総数,パターン1文字あたり数 μs 程度のオーダーで数えることができています.ただベンチマークのアルファベットサイズは2文字なので,現実的にはもっと時間がかかります(計算量はアルファベットサイズに対してに比例する).
発展
Run-length FM-index 中のやは, run の数が少ない理想的な状況下だとスパースな配列になるので,それに特化したビットベクトルを用いることでさらなる改善ができるかもしれません.また Gagie et al (2019) は,run-length FM-index の改良の一つとして,の表現に predecessor search に特化した構造を用いています.
さらに,著者らはパターンの位置を列挙する クエリに必要な索引サイズをも run の数 についての空間計算量で抑える手法を同論文で提案しています(普通は suffix array を適当にサンプルして保存するので,時間計算量を抑えるためにはの領域が必要).実装は https://github.com/nicolaprezza/r-index に公開されています.本当はこちらを advent calendar で紹介する予定でしたが時間の都合上叶いませんでした.無念.
まとめ
FM-index に連長圧縮を組み合わせた run-length FM-index を紹介しました.パターンの繰り返しが非常に多いテキストについて全文検索を行いたいときには採用を検討してみるといいかもしれません.
記事でおかしい点を発見された際は twitter 等でご連絡いただければ幸いです.PR も大歓迎です.
修正履歴
参考文献
Veli Mäkinen and Gonzalo Navarro. 2005. Succinct suffix arrays based on run-length encoding. Nordic J. of Computing 12, 1 (March 2005), 40-66.
Daisuke Okanohara and Kunihiko Sadakane. 2009. A Linear-Time Burrows-Wheeler Transform Using Induced Sorting. In Proceedings of the 16th International Symposium on String Processing and Information Retrieval (SPIRE '09), 90-101.
Francisco Claude and Gonzalo Navarro. 2012. The wavelet matrix. In Proceedings of the 19th international conference on String Processing and Information Retrieval (SPIRE'12), 167-179.
Travis Gagie, Gonzalo Navarro, and Nicola Prezza. 2018. Optimal-time text indexing in BWT-runs bounded space. In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA '18), 1459-1477.