【圧縮文字列索引連載part3】Burrows-Wheeler TransformとFM-index

はじめに

本記事は圧縮文字列索引についての連載第3弾である。

前半は、文字列圧縮のために提案された概念であるBurrows-Wheeler Transform(BWT)*1について、そしてそれを用いた文字列検索アルゴリズムについて説明する。

後半にはBWTを利用した文字列索引の一種であるFM-indexを紹介する。FM-indexの系譜をたどれば、最後には2022年現在トップクラスの性能を誇る索引であるr-index*2へと行き着く。その源流を理解しよう。

はじめに
BWT
FM-index -後方探索を使った文字列索引-
おわりに
補足

BWT

定義

まずBurrows-Wheeler Transform(BWT)を定義しよう。

文字列 $T[1,n$ ]に対するBWTとは、 $T$ を並び替えたシンボル列 $L[1,n$ ]で、以下の条件を満たすものである。

$\displaystyle L[i] = \left\{ \begin{array}{} T[SA[i]-1] & (SA[i]\neq 1) \\ \$ & (SA[i]=1) \end{array} \right.$

定義からだとイメージが掴めないので、例のごとく $T=abracadabra\$$ に対してBWTを構築してみよう。

要は、接尾辞を辞書順にソートして(ここまでは接尾辞配列と同じ。part2参照)、それら接尾辞の直前に出現したシンボルを並べてできるシンボル列がBWTである。なお、該当する辞書順の接尾辞 $T[SA[i], n]$ が文字列 $T$ そのものであった場合は、例外的に $L[i]=\$$ とする。

豆知識

ちなみに、なぜBWTを表す記号として $L$ が使われるかというと、文字列の回転を辞書順に並べてできる行列の、最後の列("L"ast column)がBWTに一致するからである。というより、当初はそもそもそういう定義で提案された。先述の接尾辞配列を用いた機械的な定義は、のちに整理されたものである。

BWTの異なる定義。Tの回転(任意の位置で二分し、順序を入れ替えたもの)を辞書順ソートしたとき、最後の列LがBWTに一致する。最初の列Fは、全てのシンボルを辞書順に並べたものとなっている。

当初はブロックソートという文字列圧縮手法のために提案されたBWTだが、現在ではその枠を超えた応用が考えられている。では次節以降で何に役立つのかを述べていこう。

BWTを使ってできること：後方探索

結論から述べると、BWTはパターンを検索するために使える。それだけならば接尾辞配列でも可能だが、それに加えて接尾辞配列よりも圧縮しやすいことが、BWTの注目に値する特性である。直感的には、辞書順で隣接する接尾辞たちは似た文脈を持っており、それらの直前に出現するシンボルも同じものになりやすい。つまり、冗長な文字列のBWTでは、同じシンボルが何個も続きやすい傾向にある。

$count(P[1,m])$ の計算に入ろう。part2にて言及した、パターンと接尾辞配列上の区間が対応するという性質を思い出してほしい。

この区間をうまいこと求めるのにBWTは使える。このとき、接尾辞配列の場合とは異なり、二分探索を必要とせず、元の文字列 $T$ を持っている必要もない。その代わり、配列 $C[1,\sigma]$ と操作 $\mathrm{rank}_c(L,i)$ を必要とする。これらが何かはすぐに説明する。

あるパターン $Q$ と、接尾辞配列上の区間 $[s,e]$ が対応していると仮定しよう。このとき、 $Q$ の左に任意のシンボルを付け加えた新しいパターン $cQ$ に対応する区間を、

$\displaystyle \left[C[c]+\mathrm{rank}_c(L,s-1)+1, C[c]+\mathrm{rank}_c(L,e)\right]$

という式で計算することができるのだ。この式における見慣れないものたちは

$C[c]=c$ よりも小さいシンボルの出現回数
$\mathrm{rank}_c(L,i) =~$ BWTの先頭 $i$ 個のシンボル $L[1,i$ ]中の $c$ の個数

を意味する。パターンを1文字左に拡張する様子から、この計算に基づく区間の更新操作をleft-extensionと呼ぶことにする。

図を用いて、これが正しい理由を説明しよう。

$cQ$ に対応する接尾辞配列上の区間を求めたい。まず辞書式順序の性質より、シンボル $c$ から始まる接尾辞は、それよりも小さいシンボルから始まる接尾辞の後に並ぶ。この分のずれを $C[c]$ で求める。次に求める必要があるのは、シンボル $c$ から始まりつつもパターン $cQ$ より辞書順が小さい接尾辞の数である。ここで着目するのは、接尾辞 $cA$ と $cB$ の大小は、先頭の同一シンボルを除いた $A$ と $B$ の大小と一致する、という事実である。この事実より、「シンボル $c$ から始まりつつもパターン $cQ$ より辞書順が小さい接尾辞の数」は、「パターン $Q$ よりも辞書順が小さい接尾辞で、直前のシンボルが $c$ であるものの数」と言い換えることが可能になる。パターン $Q$ よりも辞書順が小さい接尾辞は接尾辞配列上で $[1,s-1]$ にあり、このうち直前のシンボルが $c$ であるものとはすなわち、BWTの値が $c$ であるものに他ならない。ゆえに $\mathrm{rank}_c(L,s-1)$ が求める値になるわけだ。区間の右端についても同様で、「辞書順がパターン $Q$ 以下で、直前のシンボルが $c$ である接尾辞の数」を数えればよく、これらをまとめて先述の区間を得る。

なお、「 $C[1,\sigma]$ や $\mathrm{rank}_c(L,i)$ をどうやって計算するのか？」「サイズはどのぐらいになるのか？」という疑問は当然出るが、これは索引の実装に依存する問題なので、ここでは置いておいて後半のFM-indexで説明する。さしあたってはこういった操作がコンパクトかつ高速に計算できるという前提で議論を追ってほしい。

left-extensionができれば、あとは話が早い。

空のパターン $\epsilon$ を表す区間 $[1,n]$ から開始する。
パターンの末尾から1文字ずつleft-extensionを適用していく。
- 途中で区間の大きさが0になれば(区間の左端が右端よりも大きくなることで表される)、出現数0なので $count(P)=0$ で終了する。
$P[1]$ まで拡張したのち、区間の大きさ $e-s+1$ が $count(P)$ の答えになる。

後方探索(backward search)。慌てず挫けずleft-extension。

接尾辞配列ではパターンの先頭から1文字ずつ比較していくのと対照的に、 BWTを用いたパターン検索では、パターンの末尾から1文字ずつ左へマッチングしていくことが分かるだろう。このような特性から、このアルゴリズムは後方探索(backward search)と呼ばれる。

後方探索の時間計算量

後方探索は二分探索を必要としない点で優れている。 left-extensionを計算する時間を $t_{LF}$ とおいたとき、後方探索によるcountの時間計算量は $O(mt_{LF})$ となる。 $t_{LF}$ は実装依存で変わってくるが、基本的には $\log n$ よりも速いことが期待される。

$t_{LF}$ との命名は「LF-mapping」に由来するものだ。以下の関数 $LF$ を考えてほしい。left-extensionと同じ計算を、シンボル $L[i]$ と辞書順 $i$ に対し適用するものだ。

$\displaystyle LF(i) = C[L[i]]+\mathrm{rank}_{L[i]}(L,i)$

これは接尾辞配列上の位置 $i$ を、新たな位置 $LF(i)$ にマッピングする関数になっている。新たな位置 $LF(i)$ は、 $i$ 番目の接尾辞にその直前のシンボル $T[SA[i]-1]$ を足した接尾辞を指し示している。これが、接尾辞の回転行列における最後の列 $L$ と最初の列 $F$ を対応付けるものであることからLF-mappingと呼ばれる。下図に例を示しておこう。

LF-mappingの持つ性質を式にすると以下のようになる。

$\displaystyle SA[LF(i)] = \left\{ \begin{array}{} SA[i] - 1 & (SA[i]\neq 1) \\ n & (SA[i]=1) \end{array} \right.$

接尾辞配列上の位置で見るとかなりざっくばらんとした並び替えだが、文字列位置でみると「一個左にずらす」操作をしている。このように、文字列位置と接尾辞配列上の位置(＝辞書順)という二つの値は、全く異なる挙動をすることに注意されたい。

後方探索について説明が終わったところで、これを利用した初めての索引であるFM-indexを紹介していこう。

FM-index -後方探索を使った文字列索引-

概要

FM-index*3は、2005年に提案された、countとlocateを実現する圧縮文字列索引だ。名前は著者のイニシャルから取ったっぽい(そもそも元論文では名前を付けていないので他者による後付けだ)。

countはここまでで説明してきた後方探索により計算する。 $C[c]$ と $\mathrm{rank}_c(L,i)$ の計算方法については説明を後回しにしてきたが、まず $C[1,\sigma]$ は全ての値を配列に保存すればよい。アルファベットサイズ $\sigma$ は一般的には定数(文字列長 $n$ がいくら大きくなろうと変わらない)と考えられるので、これで十分にコンパクトだ。 $\mathrm{rank}_c(L,i)$ については、 $L$ そのものは文字列と同じサイズになってしまうことに注意する必要がある。これをどうにか圧縮しつつ、その上で $\mathrm{rank}$ という操作を高速に行いたい。この二つの要求に対し、元論文ではBW_RLXというアルゴリズムで応えている。なお、現在では「ウェーブレット木」*4という簡潔データ構造を用いた方が性能が高いことが分かっているので、こちらを採用して性能を評価する。ウェーブレット木はあくまで道具として利用するので、アルゴリズムの詳細が気になった場合は別途調べられたい。既存の構造をコンポーネントとして新たな構造を作っていく連鎖は簡潔データ構造あるある。

BWT $L$ に対してウェーブレット木を構築すると、サイズ $O(nH_k(T))+o(n\log\sigma)$ ビットで、 $L[i]$ へのアクセスおよび $\mathrm{rank}$ の計算が $O(\log\sigma/\log\log n)$ 時間*5になる。ここで $H_k(T)$ は文字列 $T$ の $k$ 次経験エントロピーだ。簡単に言うと、 $T$ の冗長さ・圧縮しやすさを示す値で、 $T$ が冗長であるほど小さくなる。経験エントロピーについても別途エントリを書く予定だ。

次にlocateなのだが、こちらにも工夫が要る。もし仮に、我々が接尾辞配列本体をそのまま持っていたとすれば、後方探索により接尾辞配列上の区間 $[s,e$ ]を求めたのち、 $SA[s,e$ ]を列挙すれば済む話だ。ただ、この話は「接尾辞配列が大きすぎる」というところから始まっている。そのまま持つわけにはいかないのだ。そこで、接尾辞配列の値を全て保存するのではなく、一定間隔ごとにサンプリングして保存するという考えを採用する。サンプリング間隔を $s$ とおくと、 $\lceil n/s\rceil$ 個の値を保存することになる。

なぜサンプリングした値だけで成立するのかというと、LF-mappingが使えるからだ。LFは文字列位置を一個ずらす操作であることを思い出してほしい。LFを何度も計算していれば、いつかはサンプリングされた接尾辞配列の値に行き着く。つまり、 $SA[p$ ]が求めたい値だったとき、 $p$ に対して繰り返しLFを適用し、 $SA[LF^k(p)]$ がサンプリングされた値だった時点で、その値に $k$ を足した $SA[LF^k(p)]+k$ が求める $SA[p]$ の値になるわけだ。このようにして $SA[p]$ を求める過程を $p\in[s,e]$ の全てで行うので、パターンの出現回数を $occ$ としたとき、LFを合計で $O(s\cdot occ)$ 回行うことになる。

時間計算量

$t_{LF}=O(\log\sigma/\log\log n)$ である。

よって、countには $O(m\log\sigma/\log\log n)$ 時間かかる。一方、locateには $O( (m+s\cdot occ)\log\sigma/\log\log n)$ 時間かかる。

空間計算量

$SA$ のサンプルも加えて $O(nH_k(T)+\frac{n}{s}\log n)+o(n\log\sigma)$ ビット必要になる。

おわりに

FM-indexがどうやってcountとlocateを実現しているか詳細に説明してきた。なお、この実装では二点注意するべきところがある。

経験エントロピーという値は頭打ちになる。つまり、巨大で冗長な文字列に対して構築するFM-indexは、およそ $n$ に比例して大きいサイズになってしまう。もっと圧縮したい。
locateの時間計算量を保証するためには、 $SA$ の一定間隔サンプルを保存する必要がある。これも $n$ に比例して大きくなってしまう。一定間隔でのサンプリングを使わないアイディアが欲しい。

次回は、どのようにこれらの問題を解決するかを説明していきたい。その議論は、最新の索引であるr-indexにつながっていく。

補足

ちなみに、FM-indexにちょっと構造を加えるだけで元文字列 $T$ へのランダムアクセスも可能になる。もちろん、 $T$ そのものは持っている必要はないので捨ててしまって構わない。

まず、逆接尾辞配列 $ISA[1,n]$ を考える。これは文字通り接尾辞配列の逆で $SA[ISA[i]]=i$ を満たすものだ。意味合いとしては、「文字列位置を与えられると、そこから始まる接尾辞の辞書順を返す配列」となる。 $SA$ が $\{1,2,\cdots,n\}$ の順列であることから、 $\{1,2,\cdots,n\}$ から $\{1,2,\cdots,n\}$ への関数として見ると全単射であるため、 $ISA$ は必ず存在する。この $ISA$ についても一定間隔のサンプルを保存する。

$T[i]$ にアクセスしたいときは、まず $ISA$ サンプルのうち $i$ の直後にあるものを見つけて $ISA[i+d]$ とする。それに $(d-1)$ 回LFを適用した値 $LF^{d-1}(ISA[i+d])$ が、位置 $i+1$ から始まる接尾辞の辞書順となっている。したがってそこのBWTの値 $L[LF^{d-1}(ISA[i+d])]$ が求めたい $T[i]$ だ。おしまい。

時間計算量は、LFを $O(s)$ 回適用することから $O(s\log\sigma/\log\log n)$ 時間である。なお、部分文字列 $T[i,j]$ に一気にアクセスするときには、いちいち全ての文字で $O(s)$ 回LFを計算する必要がないので $O( (j-i+s)\log\sigma/\log\log n)$ 時間となる。

空間計算量は、これまで述べてきたFM-indexに追加の $ISA$ サンプル $\frac{n}{s}\log n$ ビットを加えたものとなる。

*1:Burrows, Michael, and David J. Wheeler. "A block-sorting lossless data compression algorithm." (1994).

*2:Gagie, Travis, Gonzalo Navarro, and Nicola Prezza. "Fully functional suffix trees and optimal text searching in BWT-runs bounded space." Journal of the ACM (JACM) 67.1 (2020): 1-54.

*3:Ferragina, Paolo, and Giovanni Manzini. "Indexing compressed text." Journal of the ACM (JACM) 52.4 (2005): 552-581.

*4:Grossi, R. "Higher order entropy analysis of compressed suffix arrays." DIMACS Workshop on Data Compression in Networks and Applications, 2003. 2003.

*5:厳密に言うと、これはMultiary wavelet tree(多分岐ウェーブレット木)を使った場合の最新の成果のはず

u-ar’s blog

研究や技術について。もろもろ。