接尾辞配列(Suffix Array)の効率的な実装手法について
Ge Nong, Sen Zhang and Wai Hong Chan: “Linear Suffix Array Construction by Almost Pure Induced-Sorting,” Data Compression Conference (2009)で発表された手法の論文を読むメモ
ぶっちゃけSAIS(Suffix Array - Induced Sorting) - EchizenBlog-Zweiを読んだほうが早い気がするからちょっと読んでくる。
さらにアリ本のiwiwi氏がSA-IS - (iwi) ｛ 反省します - TopCoder部を書いているので更に読んでくる。
手法の全体像
接尾辞配列において最も計算量がかかるのは文字列の比較である。文字列の比較は通常の数字の比較と異なって文字列の長さに比例したオーダの計算量がかかってしまう。つまり、通常の配列のソーティングにO(N log N)かかるなら接尾辞配列の構築にはO(N^2 log N)かかってしまうということである。ここで必要になるのは再帰的なアプローチを用いて一度比較した部分を2度と比較しない工夫が必要になる。つまり最適なアルゴリズムは何らかの形で再帰的な形を取るということである。
このプログラムの全体像は以下のようになる。意味不明ですね。
def SA-IS(S, SA):
# S is the input string;
# SA is the output suffix array of S;
t: array [0..n − 1] of boolean;
S1 : array [0..n1 − 1] of integer;
P1 : array [0..n1 − 1] of integer;
B: array [0.. Σ(S) − 1] of integer;
Scan S once to classify all the characters as L- or S-type into t;
Scan t once to find all the LMS-substrings in S into P1;
Induced sort all the LMS-substrings using P1 and B;
Name each LMS-substring in S by its bucket index to get a new shortened string S1;
if Each character in S1 is unique
then
Directly compute SA1 from S1;
else
SA-IS(S1, SA1); # where recursive call happens
Induce SA from SA1;
return
L-type, S-typeの定義
まず、文字列全体をS[n]とし、終端S[n-1]には任意のアルファベットより辞書順で小さい番兵$が入っているとする。ここで接尾辞S[i:n-1]をsuf(S,i)と表現する。ここで任意の接尾辞に対してsuf(S,i)
mississippi$ LSLLSLLSLLLS abracadabra$ SSLSLSLSSLLS
のようになる。このとき、S[i]がSなのはS[i]
two stage法
ここで先にL, Sの概念だけで理解できるtwo-stage法を説明する。http://homepage3.nifty.com/DO/two-stage.htmを参考にした。
まず、性質「S[i]==S[j]でありt[i]==S, t[j]==Lの時suf(S,i)>suf(S,j)が成り立つ」を示す。まず、$が1回しか出現しない以上S[i]=$とはなりえない。この時、S[i]==S[j]よりsuf(S,i)とsuf(S,j)の大小はsuf(S,i+1)とsuf(S,j+1)の大小に等しい。またS[i],S[j]のタイプよりsuf(S,i)
- 構築済みのSuffix Arrayの上から順に1つsuffix suf(S,i)を持ってくる。
- そのsuffixの1つ前suf(S,i-1)がLタイプであればSuffix ArrayのS[i-1]に該当するbucketの一番上を埋める。
これを繰り返すだけである。なぜこれが成り立つかというと、S[i]がLタイプである時suf(S,i)>suf(S,i+1)が言える。つまりLタイプsuffixは計算済みのsuffixの1個長いやつを調べれば必ずあるってわけですね。次にi=!jに対してS[i]=S[j]をみたすようなLタイプのS[i], S[j]に対してsuf(S,i)
Left Most S-type
次にLMS character, LMS-substringを定義しよう。LMSはLeft Most S-typeの略であり、LMS characterとはS[i], i=1,2,...,n-1に対してS[i]がSタイプでありかつ左隣S[i-1]がLタイプであるものを表す。なんでS[0]は対象外なんでしょうね。さらにLMS-substringを次のように定義する。
- 番兵$自身はLMS-substringである
- S[i:j]でS[i], S[j]がLMS characterであり、S[i+1]からS[j-1]の間にLMS characterが存在せず、更にi
具体例で考えると、LMS characterは下でアスタリスクで示した文字。
mississippi$ LSLLSLLSLLLS * * * *
LMS-substringは
- issi
- issi
- ippi$
- $
である。
ここでLMS-substringを文字列の単位と見なすことで再帰的にアルゴリズムを構築できる（なんでやねん(´・ω・`)）ので、LMS-substringの順序比較を考える。2つのLMS-substringを比較するにはそれぞれの構成文字の対を左から右へ順に見ていく。つまりa[0]とb[0], a[1]とb[1]…のように。
- 文字自体がa
- 文字自体はa=bならばそのタイプを比較する。S>Lとする。
定義から考えて当然だが、この比較順序では順位が一致することはあり得る。ここでこの順位を用いてLMS-substringをソートするが、このソートはLMS-substringのポインタの配列P1を用いてO(n)で行うことができる。また、唐突な。P1はLMS-substringを元の文字列に登場する順に格納したポインタの配列である。ここでLMS-substringをbucketソートする。（頭1文字みてどうのこうのなのかは不明）その後P1の中身をソートした時のbucketのindexで名前付けしてS1を得る。ここでS1の長さはSの半分以下である。なぜなら左端の文字がLMSでなく、S中のどの隣接する2文字もともにLMSではありえないからである。さらに、S[n-1], $がLMS-substringである以上S1の最後も$を一意に指すインデックスになる。
次にCoverageという補題を示す。S1[i]=S1[j]ならばP1[i+1]-P1[i]=P1[j+1]-P1[j]である。証明は以下の通り。まずS1の最後の文字は$に対応するインデックスなのでP1[i+1],P1[j+1]は存在する。この時P1はS中のLMSを順に保持しているので、S1[i]=S1[j]ということは、S[P1[i]:P1[i+1]]=S[P1[j]:P1[j+1]]かつt[P1[i]:P1[i+1]]=t[P1[j]:P1[j+1]]が成り立つ。これの長さが一致しないといけないのでP[i+1]-P[i]=P[j+1]=P[j]が成り立つ。
次に順序の保存を示す。S1の任意の2つの接尾辞suf(S1,i)とsuf(S1,j)の順序はSの接尾辞suf(S,P1[i])とsuf(S,P1[j])の順序と一致する。証明は以下に示す。
- S1[i] != S1[j]が成り立つ場合を考える。この時元のSにおけるLMS-substringであるS[P1[i]:P[i+1]], S[P1[j]:P1[j+1]]が辞書順もしくはタイプで異なっている。辞書順で異なっているときは明らかにsuf(S,P1[i])とsuf(S,P1[j])の順序に一致する。またタイプが異なっている場合を考える。ここでS1[i]
- S1[i] == S1[j]が成り立つ場合を考える。このときはsuf(S1,i), suf(S1,j)の大小はsuf(S1,i+1), suf(S1,j+1)の大小に一致し、これをS1[i+k] != S2[j+k]が成り立つまで再帰的に進めるとsuf(S1,i), suf(S1,j)の大小はsuf(S1,i+k), suf(S1,j+k)の大小に一致する。ここで、Lemma: Coverageの証明にあるようにS[P1[i]:P1[i+k]]=S[P1[j]:P1[j+k]]が言える。またsuf(S1,i+k), suf(S1,j+k)の大小は先の場合分けよりsuf(S,P1[i+k]), suf(S,P1[j+k])に一致するのでこれに共通の接頭辞をつけてもsuffixの大小は変化せず、結局suf(S,P1[i]), suf(S,P1[j])の大小に一致する。
これだけ掛けて何を証明したかったかというと、SのLMS characterで始まるsuffixの大小関係とS1のsuffixの大小関係は等しいということである。つまり、半分の計算量で計算できるわけである。ではS1のSuffix Array SA1からSのSuffix Array SAを構築する方法を見る。
LMS-substringの接尾辞配列SA1から本来の接尾辞配列SAを作成する方法: Induced sort
当然のことながら、SAの中で同じ文字から始まる接尾辞は連続して位置している。ここでtwo-stage法同様に同じ文字から始まるSAの部分配列をbucketと呼ぶ。さらに、同じ文字を共有するsuffix中でLタイプのものはSタイプよりも順位が大きいので、これらは塊となって存在していて、早い方、左がLで右がSになる。つまり1つのbucketの中にL-bucketとS-bucketの2つのsub bucketが存在することになる。
この時S1の接尾辞配列の要素SA1[i]をSAのbucketに入れるということはP1[SA1[i]]からなる接尾辞suf(S,P1[SA1[i]])をSAの適当なバケツに登録するということである。この定義を利用して、SA1からSAを作るアルゴリズムは次のように表される。
- SA1の各要素をSAのSバケツの最後に入れていく。;
- SAの全要素SA[i]に対して前から順番に、S[SA[i]-1]がLタイプならSA[i]-1を適合するLバケツの先頭に入れていく。これを繰り返す。two-stage法と同様の議論？
- さっきと逆の事をSタイプに対しても行う。SA[i]を後ろから順に、条件に合えばSバケツの末尾に。
どうみてもO(n)なのでこれの正当性を示す。まず3ステップ目はtwo-stage法まんまである。次に2ステップ目も同様に言える。Lの右にはLかLMSしかいない＝LMSさえあれば任意のL-suffixは取り出せるからかな？
Induced sortのLMS-substringへの応用
で、この手法がまんまLMS-substringの比較へ応用できるらしい。違いはさっきのアルゴリズムのステップ1を変更するだけである。
- 各Sバケツの最後を見つける; Sの全LMS-suffixを頭1文字を見てSバケツに放り込む。
ここで今後の議論のために任意の0<=i
定理:この改造したInduced sortで2文字以上のLMS-prefixをソートできる
証明を示す。まず初期的に＝ステップ1で長さ1のSタイプのLMS-prefixはSAのバケツに入っている。（よくわからない）
ステップ2でL-typeのLMS-prefixをソートできることを示す。まず最初のL-type LMS-prefixは適切に配置される。次にk番目までのL-type LMS-prefixは適切にソートされている時にk+1番目のL-type LMS-prefix pre(S,i)をバケツに入れた時に既にそれより前により大きなL-type LMS-prefix pre(S,j)が存在したとする。この時同じバケツであることからS[i]=S[j]でpre(S,j+1)>pre(S,i+1)が言えて、pre(S,i)を入れる前に既に順番が間違っていたということになる。これは仮定に反するので、順番は揃う。
ステップ3の証明は2の時と同様である。注意すべきなのは、ステップ1でバケツに入れたsizeが1のLMS-prefixは上書きされる。番兵$はステップ1で入ってそのままなので大丈夫。
この定理から
- 任意のLMS-substringは2文字以上のLMS-prefixか番兵$なのでLMS-prefixと番兵がソートできることはLMS-substringがソートできることに等しい
- 任意のS-substringはLMS-prefixのprefixか番兵なので（ｒｙ
が言える。
具体例
"mmiissiissiippii$"という文字列に対してSAを構築することを考える。
0. L/Sの配列tを構成してLMSを探す。
Index: 00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
S : m m i i s s i i s s i i p p i i $
t : L L S S L L S S L L S S L L L L S
LMS : * * * *
1. 各文字に対してバケツを用意、放り込む
Bucket: $ i m p s
SA : {16} {-1 -1 -1 -1 -1 10 06 02} {-1 -1} {-1 -1} {-1 -1 -1 -1}
これが長さ1のLMS-prefixのソートに相当する
2. Lバケツ埋めを開始: Vでバケツの先頭を、今読んでるSAを@で示す
V@ V V V V
SA : {16} {-1 -1 -1 -1 -1 10 06 02} {-1 -1} {-1 -1} {-1 -1 -1 -1}
SA[0]-1は15でL-typeなのでiバケツの先頭へ放り込む
V @ V V V V
SA : {16} {15 -1 -1 -1 -1 10 06 02} {-1 -1} {-1 -1} {-1 -1 -1 -1}
V @ V V V V
SA : {16} {15 14 -1 -1 -1 10 06 02} {-1 -1} {-1 -1} {-1 -1 -1 -1}
V V @ V V V
SA : {16} {15 14 -1 -1 -1 10 06 02} {-1 -1} {13 -1} {-1 -1 -1 -1}
V V @ V V V
SA : {16} {15 14 -1 -1 -1 10 06 02} {-1 -1} {13 -1} {09 -1 -1 -1}
V V @ V V V
SA : {16} {15 14 -1 -1 -1 10 06 02} {-1 -1} {13 -1} {09 05 -1 -1}
V V @ V V V
SA : {16} {15 14 -1 -1 -1 10 06 02} {01 -1} {13 -1} {09 05 -1 -1}
V V @ V V
SA : {16} {15 14 -1 -1 -1 10 06 02} {01 00} {13 -1} {09 05 -1 -1}
V V V @ V
SA : {16} {15 14 -1 -1 -1 10 06 02} {01 00} {13 12} {09 05 -1 -1}
V V V @ V
SA : {16} {15 14 -1 -1 -1 10 06 02} {01 00} {13 12} {09 05 08 -1}
V V V V
SA : {16} {15 14 -1 -1 -1 10 06 02} {01 00} {13 12} {09 05 08 04}
3. Sバケツ埋めを開始: Vでバケツの末尾を、今読んでるSAを@で
V V V V @V
SA : {16} {15 14 -1 -1 -1 10 06 02} {01 00} {13 12} {09 05 08 04}
V V V V @ V
SA : {16} {15 14 -1 -1 -1 10 06 03} {01 00} {13 12} {09 05 08 04}
V V V @V V
SA : {16} {15 14 -1 -1 -1 10 07 03} {01 00} {13 12} {09 05 08 04}
V V @ V V V
SA : {16} {15 14 -1 -1 -1 11 07 03} {01 00} {13 12} {09 05 08 04}
V V @ V V V
SA : {16} {15 14 -1 -1 02 11 07 03} {01 00} {13 12} {09 05 08 04}
V V @ V V V
SA : {16} {15 14 -1 06 02 11 07 03} {01 00} {13 12} {09 05 08 04}
V V V V V
SA : {16} {15 14 10 06 02 11 07 03} {01 00} {13 12} {09 05 08 04}
4. この結果として次が得られる（わからん）
S1: 2 2 1 0
*1:smallとlargeの頭文字か
