LZ77分解

高速文字列解析の世界を読んでいたら,LZ77分解は接尾辞配列を使ってO(n)で求められると書いてあったので調べた.

LZ77

LZ77は繰り返し出現した文字列を過去に出現した文字列のコピーとして表現することで圧縮を実現する.
文字列のコピーは(出現位置, 長さ)のペアで表現する.
ただし,はじめてでてきた文字は(出現位置,長さ)ではなく(文字, 0)とする.

例としてacaaacatatという文字列のLZ77分解を考える.

0 1 2 3 4 5 6 7 8 9
a c a a a c a t a t
  1. 0番目の文字aは初めてでてきたので('a', 0)
  2. 1番目の文字cは初めてでてきたので('c', 0)
  3. 2番目の文字aは,0番目から1文字一致するので(0, 1)
  4. 3, 4番目の文字aaは,2番目から2文字一致するので(2, 2)
  5. ...

このように計算していくと最終的に('a', 0), ('c', 0), (0, 1), (2, 2), (1, 2), ('t', 0), (6, 2) を求めることができる.

接尾辞配列を用いたLZ77分解

PSV(previous smaller values)とNSV(next smaller values)という配列を使うことでLZ77分解を線形時間で実行することができる.
まず  {PSV_{lex}} {NSV_{lex}} {PSV_{text}} {NSV_{text}}を定義する.

  •  PSV_{lex} [i] = min \{ j \in [1..i - 1 ] \mid  SA [j] \lt SA [i] \}
  •  NSV_{lex} [i] = min \{ j \in [i + 1..n] \mid  SA [j] \lt SA [i]  \}
  •  PSV_{text} [i]  = SA [PSV_{lex} [ISA [i ] ] ]
  •  NSV_{text} [i]  = SA [NSV_{lex} [ISA [i ] ] ]

要素が存在しない場合は-1とする.

acaaacatatの場合は下のようになる.

index SA psv_text nsv_text psv_lex nsv_lex 接尾辞
0 2 -1 -1 -1 2 aaacatat
1 3 0 -1 0 2 aacatat
2 0 -1 0 -1 -1 acaaacatat
3 4 2 0 2 6 acatat
4 8 0 1 3 5 at
5 6 1 -1 3 6 atat
6 1 4 1 2 -1 caaacatat
7 5 5 -1 6 -1 catat
8 9 4 6 7 9 t
9 7 5 7 7 -1 tat

例として {PSV_{text}[6] }を定義に従って求めてみる.
1. まず{ISA[6]}を求める.これは{SA[i] = 6}であるiなので5となる.
2.  {PSV_{lex}[5]}は3であり,また{SA[3]}は4である.
3. よって {PSV_{text}[6]}は4となる.

もう少し直感的にいうと, {PSV_{text}[6] }{ISA[6]}(表のindex5の位置)からSAの値を上にみていったときにでてくる最初の{SA[5]}より小さい値である.
{SA[5] = 6}の位置からSAの値を上にたどっていくと8, 4, 0, 3, 2とあるが,一番最初にでてくる6より小さい値は4なので, {PSV_{text}[6]}は4となる.
そのため {PSV_{text}[6]}は4となる.

結局 {PSV_{text}}がなんなのかというと,文字列の各位置より前にでてきた中で一番近い接尾辞の位置ということになる.
現在acaaacatatの5番目まであるacaaacまで処理済みとする.
この次の目標は6番目以降であるatatと最も近い文字列を探すことである.
接尾辞配列は接尾辞をソートしたものなので,atatである表のindx5から上(または下)を順番にみていけば,atatに近い順に接尾辞をみていくことになる.
ただし,LZ77分解では未処理である6番目以降の文字を参照することはできないので,6より小さいSAを探している.
このようにすることで現在処理している文字列に最も近い過去の文字列を求めることができる.

文字列sの各factorを求めるアルゴリズムは以下のようになる.( {NSV} {PSV}はそれぞれ {NSV_{text}} {PSV_{text}}の値である)
i番目のfactorの開始位置はpsv[i]かnsv[i]のどちらかであるので,lcpが大きい方を採用するということをしている.
LCP(最長共通接頭辞)は手抜き.

# s[i:]とs[j:]の最長共通接頭辞
def lcp(i: int, j: int, s: str):
    if i < 0 or j < 0:
        return 0
    for i, (a, b) in enumerate(zip(s[i:], s[j:])):
        if a != b:
            return i
    return i + 1

def lz_factor(i: int, psv: int, nsv: int, s: str):
    if lcp(i, psv, s) > lcp(i, nsv, s):
        p, l = psv, lcp(i, psv, s)
    else:
        p, l = nsv, lcp(i, nsv, s)

    # special factor
    if l == 0:
        p = s[i]

    factor = (p, l)
    next_idx = i + max(l, 1)
    return factor, next_idx

i = 0
while i < n:
    factor, i = lz_factor(i, psv[i], nsv[i], s)
    print(factor)

PSVとNSVを求める

PSVとNSVを線形時間でもとめるKKP3について述べる(Linear Time Lempel-Ziv Factorization: Simple, Fast, Small1アルゴリズムは下のようになる.

def suffix_array(s: str):
    l = [(s[i:], i) for i in range(len(s))]
    return [i for suffix, i in sorted(l)]

s = "acaaacatat"
n = len(s)
sa = [-1] + suffix_array(s) + [-1]
top = 0
nsv, psv = [None] * n, [None] * n
for i in range(1, n + 2):
    while sa[top] > sa[i]:
        psv[sa[top]] = sa[top - 1]
        nsv[sa[top]] = sa[i]
        top -= 1
    top += 1
    sa[top] = sa[i]

空間計算量を節約するためにSAの配列をスタックとして利用している.
1. SAを上から順番にみていく.このときSAの値が大きくなるようならスタックに積む.
2. SAの値が小さいならスタックからとりだし,値に応じてNSV, PSVを埋めるということを繰り返していく.

アルゴリズムにしたがってPSV,NSVに値を埋めていってみる.
1. スタックの底として-1をいれておく
2. SA[0]=2は-1より大きいのでスタックに積む
3. SA[1]=3は2より大きいのでスタックに積む
4. SA[2]=0は3より小さいので3をとりだす.
このとき,SAが3の位置(表のindex1の箇所)から下にみていったときにでてくる最初の3より小さい値は0であるということがわかるので,PSV[3] = 0となる.
また,SAが3の位置から上にみていったときにでてくる最初の3より小さい値は2であるということがわかるので(スタックは大きければ積んでいるため),NSV[3] = 2となる.
5. SA[2]=0は2より小さいので2をとりだす.
このとき,SAが2の位置から下にみていったときにでてくる最初の2より小さい値は0であるということがわかるので,PSV[2] = 0となる.
また,SAが2の位置から上にみていったときにでてくる最初の2より小さい値は-1であるということがわかるので,PSV[2] = -1となる.
6. ...

これを繰り返すことですべてのPSVとNSVを求めることができる.

まとめ

接尾辞配列からNSVとPSVをO(n)でを求めることができ,NSVとPSVからLZ77分解をO(n)で実行することができる.
よって全体としてO(n)でLZ77分解をすることができた気がする.

コード全体は以下の通り.

def suffix_array(s: str):
    l = [(s[i:], i) for i in range(len(s))]
    return [i for suffix, i in sorted(l)]


# s[i:]とs[j:]の最長共通接頭辞
def lcp(i: int, j: int, s: str):
    if i < 0 or j < 0:
        return 0
    for i, (a, b) in enumerate(zip(s[i:], s[j:])):
        if a != b:
            return i
    return i + 1


def lz_factor(i: int, psv: int, nsv: int, s: str):
    if lcp(i, psv, s) > lcp(i, nsv, s):
        p, l = psv, lcp(i, psv, s)
    else:
        p, l = nsv, lcp(i, nsv, s)

    # special factor
    if l == 0:
        p = s[i]

    factor = (p, l)
    next_idx = i + max(l, 1)
    return factor, next_idx


def main():
    s = "acaaacatat"
    n = len(s)
    sa = [-1] + suffix_array(s) + [-1]
    top = 0
    psv, nsv = [None] * n, [None] * n
    for i in range(1, n + 2):
        while sa[top] > sa[i]:
            psv[sa[top]] = sa[top - 1]
            nsv[sa[top]] = sa[i]
            top -= 1
        top += 1
        sa[top] = sa[i]
    i = 0
    while i < n:
        factor, i = lz_factor(i, psv[i], nsv[i], s)
        print(factor)


if __name__ == '__main__':
    main()

参考

高速文字列解析の世界――データ圧縮・全文検索・テキストマイニング (確率と情報の科学)

高速文字列解析の世界――データ圧縮・全文検索・テキストマイニング (確率と情報の科学)


  1. 論文ではKKP3とKKP2が提案されている.