最長増加部分列・最長共通部分列

ある数列や文字列 $A=\{a_1,a_2,...,a_n\}$ がある時、部分列(Subsequence)とは、$A$ からいくつかの要素を抽出した列。飛び飛びでもよいが、順番は元から変わってはいけない。

$A$の部分列 $B=\{a_i,a_j,...,a_k\}$、ただし $1 \le i \lt j \lt ... \lt k \le n$

元の数列$\{5,3,7,4,1,4,6\}$
部分列の例$\{5,7,4\},\{1,6\}$

部分列にまつわる問題では、「最長増加部分列」や「最長共通部分列」を求める問題が出題される。これらは動的計画法により高速に求めることができる。また長さだけなら、なお高速に求める方法がある。

最長増加部分列の長さ

数列の部分列のうち、隣接する2要素を見ると常に右の方が大きいものを増加部分列という。同じ値を許すかは定義によるが、ここでは許さないとする。

そのうち最長のものを最長増加部分列(Longest Increase Subsequence)という。

$A$の増加部分列 $B=\{a_i,a_j,...,a_k\}$、ただし $1 \le i \lt j \lt ... \lt k \le n$ かつ $a_i \lt a_j \lt ... \lt a_k$

元の数列$\{5,3,7,4,1,4,6\}$
増加部分列の例$\{5,7\},\{3,6\}$
最長増加部分列$\{3,4,6\},\{1,4,6\}$

この長さを求めましょう、という問題。(このアルゴリズムでは、最長増加部分列自体は求まらない)

  1. $A$ の要素を1つずつ順番に処理する。現在処理中の要素を $a_j$ とする。
  2. $L[i]$ を、「$a_1$~$a_j$ の要素から長さが $i+1$ の増加部分列を作ったとき、末尾要素が取り得る最小値」とし、更新していく。
  3. 最後まで処理した後の $L$ の値が埋まっている所の長さが、最長増加部分列長となる。

from bisect import bisect_left

def lis(A: list):
    L = [A[0]]
    for a in A[1:]:
        if a > L[-1]:
            # Lの末尾よりaが大きければ増加部分列を延長できる
            L.append(a)
        else:
            # そうでなければ、「aより小さい最大要素の次」をaにする
            # 該当位置は、二分探索で特定できる
            L[bisect_left(L, a)] = a
    return len(L)

# -----------------------------------------
print(lis([3, 4, 5, 1, 6, 3, 2, 7])) # => 5

最長共通部分列の長さ

2つの文字列A,Bがある時、「Aの部分列でもあり、Bの部分列でもある」文字列を共通文字列という。

そのうち最長のものを最長共通部分列(Longest Common Subsequence)という。

Aabcbdabc
Bbdcaba
LCSbcba,bdabなど

この長さを求めましょう、という問題。(このアルゴリズムでは、最長共通部分列自体は求まらない)

  1. Bの要素を1つずつ順番に処理する。現在処理中の要素を$b_k$とする。
  2. $L[i]$を、「“$a_1$~$a_j$” と “$b_1$~$b_k$” から長さが $i+1$ の共通部分列が作れるとき、$j$ が取り得る最小値」とし、更新していく。
  3. Bを最後まで処理した後の $L$ の値が埋まっている所の長さが、最長共通部分列長となる。

def lcs(a: str, b: str):
    L = []
    for bk in b:
        bgn_idx = 0  # 検索開始位置
        for i, cur_idx in enumerate(L):
            # ※1
            chr_idx = a.find(bk, bgn_idx) + 1
            if not chr_idx:
                break
            L[i] = min(cur_idx, chr_idx)
            bgn_idx = cur_idx
        else:
            # ※2
            chr_idx = a.find(bk, bgn_idx) + 1
            if chr_idx:
                L.append(chr_idx)
    return len(L)

# ---------------------------------------
print(lcs('>abcbdabc', 'bdcaba'))  # => 4

  • 以下の解説にあたっての定義
    • 文字列Bの先頭 $k$ 文字を抽出した文字列を、$B'_k$ とする
    • $AT(S)$ を、「Aの部分文字列$S$ の末尾文字の、Aにおけるindex」とする
      • A:abcabc、S:bb だった場合、$AT(S)=4$
  • ※1:
    • 処理中の文字を $b_k$ とする
    • ※1の時点でbgn_idxに入っているのは、$L[i-1]$、つまり「$AT(A$と$B'_{k-1}$ の長さ $i$ の共通部分列$)+1$」
    • ここで「$A$と$B'_{k-1}$ の長さ $i$ の共通部分列」を $T$ とする(複数ある場合は、$AT(T)$ が最小のもの)
    • ループ内でおこなっているのは、以下の確認
      • $T$ の末尾に $b_k$ を付け足して長さ $i+1$ の共通部分列を作れるか?
        • つまり、$AT(T)$ 以降に、文字 $b_k$ が存在するか?
      • 作れたら、$AT(T+b_k)$ の値は?(これ+1がchr_idx
      • それは、既に発見済みの他の共通部分列$T'$ より短いか?($AT(T+b_k)<AT(T')$ なら更新)
    • 以上の繰り返しで、$L$ は更新される
  • ※2:
    • ループ途中でbreakされないで最後まで進むと、else節に入る
    • ここで、$A$と$B'_{k-1}$ の最長共通部分列を $T_{max}$、その長さを $|T_{max}|$とする
    • ※2の時点でbgn_idxに入っているのは、「$AT(T_{max})+1$」
    • Aにおいて、bgn_idxより右に文字 $b_k$ が存在する ⇔ 最長共通部分列“$T_{max}+b_k$“が存在し、暫定最長共通部分列長を$|T_{max}|+1$に延長できる
    • 延長できる場合、$L$ に追加する

具体的な最長共通部分列

$(|S|+1) \times (|T|+1)$ の配列上でDPを行い、途中の計算結果を保存しておくことで、最長共通部分列と共に、その具体例を構築することが出来る。

DPテーブルの構築

$S_{1 ... i}$ と $T_{1 ... j}$ のLCSの長さを $L_{i,j}$ とする。

$L_{0,j}=0,L_{i,0}=0$ とする。

$S_{i+1}=T_{j+1}$ の時

LCS長は1文字拡張できる。$L_{i+1,j+1}=L_{i,j}+1$

それ以外の時

LCS長は、片方の末尾を1文字削った2通りの長い方となる。 $L_{i+1,j+1}=max(L_{i+1,j},L_{i,j+1})$

これで、$i=1 .. |S|$、$j=1 .. |T|$ の2重ループでDPテーブルを小さい方から埋めていくことで、$L_{|S|,|T|}$ が最終的なLCS長となる。

具体例の構築

具体的な文字列は、DPテーブルから逆向きに復元する。$(i,j)=(|S|,|T|)$ とする。

$L_{i,j}=L_{i,j-1}$ または $L_{i,j}=L_{i-1,j}$ の時

同じである数字の方に$(i,j)$を移動する。両方同じ場合は、具体例が1つ構築できればいい場合は適当に1つ選んで移動する。

それ以外の時

$S_i=T_j$ であるはずなので、それをLCSに加え、$(i-1,j-1)$ に移動する。

      0 1 2 3 4 5 6 7
        M Z J A W X U
0   | 0 0 0 0 0 0 0 0
1 X |⓪ 0 0 0 0 0 1 1
2 M | 0❶① 1 1 1 1 1
3 J | 0 1 1❷ 2 2 2 2
4 Y | 0 1 1② 2 2 2 2
5 A | 0 1 1 2❸③③ 3
6 U | 0 1 1 2 3 3 3❹  丸付きの数字を辿っていき、
7 Z | 0 1 2 2 3 3 3④  黒丸の数字を拾っていくと、LCS(の1つ)は MJAU と復元できる
programming_algorithm/dynamic_programming/longest_common_subsequence.txt · 最終更新: 2019/01/09 by ikatakos
CC Attribution 4.0 International
Driven by DokuWiki Recent changes RSS feed Valid CSS Valid XHTML 1.0