Levenshtein distance

diffの動作原理を知る〜どのようにして差分を導き出すのか を読んでからずいぶん時間を経たけれども Diff のアルゴリズムを理解できた*1
理解するためにいろいろ検索したが Diff algorithm の解説が自分にはいちばん分かりやすかった*2。このエントリがなければまだ理解していないと思う。id:Constellation 氏に感謝したい。


Diff のアルゴリズムはレーベンシュタイン距離を求める部分がメインになるので今回はそこだけに注目する。編集グラフ上でレーベンシュタイン距離を求める際に通ったパスをたどれば Diff がわかる。
レーベンシュタイン距離は編集距離 (Edit distance) や SES (Shortest edit script) など呼ばれ方はいろいろあるみたいだ。
前回の LCS が共通部分を求めるのに対して、レーベンシュタイン距離はどれだけ変更が必要かを求める。
計算方法は2通りあるようだ

  • 追加、削除、置換をそれぞれ 1 と数える
  • 置換は追加+削除で表現できるので 2 で数える

Levenshtein distance - Rosetta Code は前者だが学んだアルゴリズムに合わせて後者で数える。

Dynamic programming O(NM)

// O(NM)
def levenshtein(a, b) {
  int N = a.size()
  int M = b.size()
  def d = new int[N+1][M+1]
  for (i in 0..N) d[i][0] = i
  for (j in 0..M) d[0][j] = j
  for (i in 0..<N)
    for (j in 0..<M)
      d[i+1][j+1] = a[i] == b[j] ?
        d[i][j] :
        Math.min(d[i+1][j]+1, d[i][j+1]+1)
  d[N][M]
}

assert 5 == levenshtein("abcabba", "cbabac")

//        a  b  c  a  b  b  a
//   [*0,*1,*2, 3, 4, 5, 6, 7]
// c [ 1, 2, 3,*2,*3, 4, 5, 6]
// b [ 2, 3, 2, 3, 4,*3, 4, 5]
// a [ 3, 2, 3, 4, 3,*4, 5, 4]
// b [ 4, 3, 2, 3, 4, 3,*4, 5]
// a [ 5, 4, 3, 4, 3, 4, 5,*4]
// c [ 6, 5, 4, 3, 4, 5, 6,*5]

DP によるものは LCS のときとほとんど変わらない。
d の表はその文字までの編集距離で (N,M) の値が最終的な編集距離になる。
* は距離を求める際に通ったパス。

Myers のアルゴリズム O(ND)

DP の処理では最終的な編集距離が 5 なのに 6 や 7 の部分も計算していた。
その部分の計算を省略するために貪欲法で編集距離を0から増やしながら (N,M) に辿り着くまで調べる。

def levenshtein1(a, b) {
  int N = a.size()
  int M = b.size()
  int MAX = M+N
  // v[k] は k で到達可能な x の位置
  def v = new int[MAX*2+1]
  // 編集距離 d は 0 (一致) から M+N (共通部分が存在しない) の間に存在する
  for (d in 0..MAX) {
    // x-y = k とし、編集距離 d のときの k 線上で到達可能な位置を調べる
    // k と d の偶奇は一致するので step は 2
    for (k in (-d..d).step(2)) {
      // k == -d の場合、k-1 は範囲外になる
      // k == d  の場合、k+1 は範囲外になる
      int x = k == -d || (k != d && v[k-1] < v[k+1]) ?
                v[k+1] :  // k+1 から y 方向へ移動して k へ
                v[k-1]+1  // k-1 から x 方向へ移動して k へ
      int y = x - k
      // a[x] == b[y] なら斜めに移動する
      while (x < N && y < M && a[x] == b[y]) (x,y) = [x+1,y+1]
      // 到達点を更新する
      v[k] = x
      if (x >= N && y >= M) return d
    }
  }
  // Bug
  assert false
}

assert 5 == levenshtein1("abcabba", "cbabac")


DP で求めた表の部分に限定する

def levenshtein2(a, b) {
  int N = a.size()
  int M = b.size()
  // a も b も空文字の場合は編集距離は 0
  if (N == 0 && M == 0) return 0
  // -M..N の範囲に限定
  def v = new int[M+N+1]
  for (d in 0..M+N) {
    for (k in (-d..d).step(2)) {
      if (k < -M || N < k) continue
      int x = k == Math.max(-d,-M) ? v[k+1] :
              k == Math.min( d, N) ? v[k-1]+1 :
                                     Math.max(v[k+1], v[k-1]+1)
      int y = x - k
      while (x < N && y < M && a[x] == b[y]) (x,y) = [x+1,y+1]
      v[k] = x
      // (N,M) に到達したときの編集距離を返す
      if (x == N && y == M) return d
    }
  }
  assert false
}

assert 5 == levenshtein2("abcabba", "cbabac")

Wu のアルゴリズム O(NP)

Myers のアルゴリズムでは 0 から調べていたが a と b の長さが違えば編集距離はその差 D 以上になるので D から探索する。

def levenshtein3(a, b) {
  int M = a.size()
  int N = b.size()
  // N >= M でなければ a と b を入れ替える
  if (N < M) return levenshtein3(b, a)
  assert N >= M
  int D = N-M
  // どちらかが空文字なら編集距離は D
  if (N == 0 || M == 0) return D
  // k = y-x とする
  // k 線上で a[x] == b[y] の部分を斜めに移動する
  def snake = { k, int y ->
    int x = y-k
    while (x < M && y < N && a[x] == b[y]) (x,y) = [x+1,y+1]
    return y
  }
  // p は k = D から離れるときのコスト
  // fp[k] はコスト p のときに k 線上で到達可能な y の位置
  // p と p-1 について必要だが更新順序を D < k と D > k で変えることによって1次元配列で足りる
  def fp = new int[M+N+3]
  Arrays.fill(fp, -1)
  int p = -1
  while (fp[D] != N) {
    p = p+1
    // -p <= k < D の場合、x 方向への移動はコストがかかる
    // max(f[p][k-1]+1, f[p-1][k+1])
    for (int k = -p; k < D; k++)
      fp[k] = snake(k, Math.max(fp[k-1]+1, fp[k+1]))
    // D < k <= D+p の場合、y 方向への移動はコストがかかる
    // max(f[p-1][k-1]+1, f[p][k+1])
    for (int k = D+p; k > D; k--)
      fp[k] = snake(k, Math.max(fp[k-1]+1, fp[k+1]))
    // k = D 場合、移動のコストはかからない
    // max(f[p][k-1]+1, f[p][k+1])
    fp[D] = snake(D, Math.max(fp[D-1]+1, fp[D+1]))
  }
  // 文字列長の差 D に 2p を加えたものが編集距離 (p は離れる分だけなので戻る分が必要)
  return D+2*p
}

assert 5 == levenshtein3("abcabba", "cbabac")


コードは論文のものを Groovy に書き換えただけなので読み返したときに思い出すためのコメントがメイン。
配列の部分は offset を指定した方がいいと思うが本質的な問題ではないので楽した。

groovy:000> Arrays.toString(new int[7].with { for (i in -3..3) it[i]=i; it})
===> [0, 1, 2, 3, -3, -2, -1]

他にもアルゴリズムがあるみたいだ。

*1:読むのに必要な時間を軽くオーバーした

*2:参考文献には論文へのリンクもある