Suffix Array を作る - SA-IS の実装

Suffix Array は今若者の間で人気のデータ構造です． マイ suffix array を実装することで，オシャレ度がアップしてモテ系になり，女子力も上がると言われています． その中でも今特に，手軽でクールな SA-IS (アルファベットサイズ固定の下で線形時間で省メモリで suffix array が作れる今最強のアルゴリズム) の実装がブームです． 僕もブームに便乗して，実装してみました．

ところで，SA-IS は流行っているので，日本語でもすでに様々なところで記事が書かれています (日付順)．

• SAIS(Suffix Array - Induced Sorting) - EchizenBlog-Zwei
• SA-IS: SuffixArray線形構築 - sileの日記
• SA-IS - (iwi) ｛ 反省します - TopCoder部
• 接尾辞配列(Suffix Array)の効率的な実装手法について - xxxxxeeeeeの日記

仕掛けが気になる人はまずここら辺を眺めて，え，そこどうなってんの，ってなったら下の論文を読むのが良いです．

• G. Nong, S. Zhang and W. H. Chan, Two Efficient Algorithms for Linear Time Suffix Array Construction, IEEE Transactions on Computers, To Appear

SA-IS の元論文は 2009 年ですが，上の論文 (まだドラフト，PDF が G. Nong さんのページにある) に C のコードがついてるので，それを見ながらがわかりやすいです (というかソース読まないとわかんない)．

実装と実験

この論文についている実装ですが，'\0' を番兵に使うのでバイナリとかだと使えません． ちゃんとした実装は sais が良さそうで，上の G. Nong さんのページからも refer されています． esaxx - C++ enhanced suffix array template library - Google Project Hosting でもこれが使われています．

僕が書いた実装は

beam2d/sara - GitHub

にあります． 本体はヘッダー sara.h 一個です． 使い方は簡単で，文字列 str に対して sara::make(str, ch_max) で suffix array が返ってきます．

上の論文にある C のコード (100 行弱) よりは長い (160 行強) ですが，番兵なしで，入力と出力がテンプレート化されてます (sais もそうですが)． 文字の型は組み込み整数型ならなんでもいいですが，値の範囲は 0 から指定した値 ch_max までに収まっている必要があります．

検証と評価用のサンプルが sara_main.cc (sara) です． 今回僕が唐突に SA-IS を実装し始めたのはレポートのためなんですが (冒頭の件？ 何の話ですか？)，その中では The Canterbury Corpus にあるデータを使って検証と実行時間の計測を行い，std::lexicographical_compare と std::sort を使うだけの単純な suffix array 構築と比較しました． 詳細は省きますが，普通の文書だと 4 倍くらいの差しかなくて，割とクイックソート速いのねという感じです． 一方，似たパターンが連続するバイナリデータだと愚直な文字列比較ソートは激遅で，例えば第一世代 Core i5 2.67GHz で， ptt5 (FAX, 501KB) は std::sort だと 28.1 秒， sara::make だと 0.040 秒でした． 繰り返しを多く含む人工データではさらに差が広がります．

適当な解説

アルゴリズムをひと通り説明します． が，正直がんばってわかりやすく説明するげんきが残っていないので，適当です！！！ わからなくなったら上のブログとか論文とかついてるソースコードとかを見るのが良いです．

各 suffix に対して，右隣の suffix より (辞書式順序で) 大きいか小さいかで L-type と S-type というラベルを振ります． L/S を並べると例えば mmiissiissiippii$ という文字列 ($ は終端文字で最小の文字で必ず S-type とする) なら LLSSLLSSLLSSLLLLS となります． 最終的な suffix array はまず先頭文字についてソートされているので，先頭文字ごとのバケツを作ります． suffix 同士の先頭文字が等しい場合，L-type の方が S-type より必ず先に来ることに注意します． つまり各バケツはさらに前半・後半で L-type / S-type にわかれます．

SA-IS ではまず L-type だけをソートします． 先頭文字と type が等しい場合，右隣の suffix の大小がそのままもとの suffix の大小になることを利用して，「まだ見ていない L-type の右隣の suffix の中で最小のものを含むソートされたリスト」を持っておき，最小のものを取り出し，左隣の suffix が L-type なら対応するバケツに入れて，リストに加える，という操作をリストが空になるまで繰り返します． このリストは，後述する初期値を除いて構築中の suffix array そのものが使えます (まだ埋まっていない部分は飛ばしながら見ます)． 次に S-type では逆向きに，今できたソートされた L-type suffix を大きい方から順に見て，左隣の suffix が S-type ならバケツの後端から順に加える，という操作を繰り返します．

これをやるためには，L-type をソートするときに「まだ見ていない L-type の右隣の suffix の中で最小のものを含むソートされたリスト」を常に持っている必要があります． これは，L の右隣にいる S (LMS: Left Most S-type という) 全部のソートされたリストをはじめに作り，上の方法で最小のものを取り出しては左隣の L-type を入れる，という操作を繰り返せばできます．

理由ですが，L-type が並んでいるところ (...LLLLLS... みたいな． この S が LMS) では右の suffix ほど小さくなります． 特にこれらの L-type の直後に来る LMS の suffix はこの L-type suffix のどれよりも小さいです． 同じように，S-type が並んでいるところでは右に行くほど大きくなります． つまり LMS は suffix の大小という意味で谷になっています． なので谷を全部持っておけば，必ず最小のものが含まれるし，取り出した最小の suffix の一個前 (山を後ろに一歩登ったところ) を加えていけば「まだ見ていない L-type の右隣の suffix の中で最小のものを含むリスト」が常に保たれることになります． S-type のソートでも同じように，今度は山のてっぺんから左に下っていく感じで進めていけばよく，これは単純にソート済みの L-type suffix とすでにソートされた S-type suffix を後ろから順に見ていけばよいです．

そのためにはじめにソートされた LMS-suffix のリストを作るのですが，ここが SA-IS 最大のポイントになります． 各 LMS 文字に対して，そこから次の LMS 文字までの substring を LMS-substring と呼びます (右端の LMS 文字を含みます)． LMS-substring は各文字と L/S-type とのペアに関する辞書式順序で比較します． LMS-substring をソートして，各 LMS-substring に対して LMS-substring 同士の大小と一致するような整数を振り，それを出現順に並べた文字列 S1 を考えます． すると，各 LMS-suffix の順序と対応する S1 の suffix の順序が一致します． なので S1 の suffix array を作れば LMS-suffix のソートされたリストが得られます． S1 の suffix array は，重複する文字がない場合にはただ文字をソートすればよいです (これは実は LMS-substring を整数に変換するところですでにやっています)． S1 に重複する文字がある場合には SA-IS を再帰的に適用します． LMS 文字は全体の半分以下しかないので，問題サイズは半分以下になっており，再帰させても線形時間で終了します．

まだ続きます． 上でサラっと LMS-substring のソートと書いた部分です． ここでも四段落上で書いたのとほとんど同じやり方が使えます． 今度は LMS-substring ごとの suffix を考え，それを全部集めたものをソートします． L/S-type はもとのと同じものが使えます． この場合，最初の「左隣が L-type な suffix のうち最小のものを含むソート済みリスト」の初期値として，LMS 文字一文字からなる suffix をバケツに入れたものが使えます． あとは上と同じように L-type，S-type の順でソートすれば，全 LMS-substring の suffix のソートが完了します． LMS-substring に対応する部分だけを取り出せば，LMS-substring のソートが得られます．

まとめると，

• 各 suffix に対して L/S-type のラベルを振る (ここは文字列を後ろから前へ一回走査するだけでできます)．
• LMS-substring をソートする． これはまず LMS 文字一文字からなる LMS-substring の suffix をバケツに入れ，そこから上記の方法で L-type，S-type の順にソートする．
• ソートされた LMS-substring に整数値を振る． 等しい LMS-substring 同士は同じ値になるように，小さい方から順に 0, 1, 2, ... と振る． これをもとの文字列での LMS の出現順と同じ順序で並べて S1 を作る．
• 整数値が被らなければ (= 重複する LMS-substring がなければ) すでに作った LMS-substring のソートが LMS-suffix のソートを与える． そうでなければ S1 に対して SA-IS を再帰的に適用して S1 の suffix array を作り，対応する LMS-suffix のソートを得る．
• LMS-suffix のソートされたものを対応する各バケツの後端に入れる． これを初期状態として，上記の方法で L-type，S-type の順にソートする．

となります．

アルゴリズムの概要はこんな感じです． 何言ってるのかわからない部分もあると思いますが，上のリンクとか論文とかを参考にしてください． 本当は山と谷の図をたくさん書いて，具体的な文字列で手でアルゴリズムを回してみるとわかりやすいです． 論文は，ソースと本文を行ったり来たりしてればわかると思います． やってることはシンプルで，ソースも短くて (100 行弱)，コメントもたくさん入っています．