15.4 最長共通部分列
生物学応用には異なる2体(またはそれ以上)の生物のDNAの比較がよく出現する。DNAの一本鎖は塩基と呼ばれる分子列から構成される。塩基はA(アデニン)、C(シトシン)、G(グアニン)、T(チミン)。
DNAの一本鎖を有限集合{A,C,G,T}上の文字列として表現する。
たとえば、ある生物のDNAはS_1=ACCCGTCGAGであり、別のある生物のDNAはS_2=ACGGTCGAAAAであるかもしれない。
2本のDNAの一本鎖を比較する目的の1つは
2体の生物の近接度を量る尺度の1つとして2本のDNAの一本鎖間の類似度を決定することである。
類似度は異なる多くの方法で定義可能であり、事実、異なる多くの方法で定義されてきた。
たとえば、あるDNA一本鎖が別のDNA一本鎖の部分文字列であるとき、これら2本のDNA一本鎖は互いに類似であるということができる。
別の可能性として、一方を他方に変換するために必要な変更の回数が少ないときにこれら2本の一本鎖は類似であるということもできる。
さらに別の方法として、S_1とS_2に共通して出現するある第3の一本鎖S_3を見つけることで一本鎖S_1とS_2の類似度を量ることもできる。S_3に現れる塩基はS_1とS_2のどちらの中でも同じ順序で現れる必要があるが、殻ずしも引き続いて現れる必要はない。長い一本鎖S_3が発見できれば類似度が高い。
ex.
S_1 = ACG'T'ACGT,
S_2 = ACG'A'ACGT
=> S_3 = ACGACGT
最長共通部分列問題
最後に説明した類似度の概念を最長共通部分列問題として定式化する。ある列の部分列はその列から0以上の個数の要素を取り去ったものである。正確には
列 X = [x_1, x_2, ..., x_m] が与えられたとき、ある列 Z = [z_1, z_2, ..., z_k] がXの部分列であるとは、真に増加する X の添字の列 [i_1, i_2, ..., i_k] が存在して、すべての j = 1, 2, ..., k に対して x_{i_j} = z_j を満たすときをいう。たとえば
Z = [B, C, D, B] は添字の列 [2, 3, 5, 7] によって対応付けられる X = [A, B, C, B, D, A, B] の部分列である。
2つの列 X と Y が与えられているとする。ある列 Z が X と Y の両方の部分列であるとき、 Z を X と Y の共通部分列という。
X = [A, B, C, B, D, A, B],
Y = [B, D, C, A, B, A]
とすると、列 [B, C, A] は X と Y の共通部分列である。しかし、最長共通部分列(Longest Common Subsequence)ではない。なぜなら長さ4の共通部分列 [B, C, B, A] が存在するからである。長さが5以上の共通部分列が存在しないから、列 [B, C, B, A] は X と Y の LCS の1つである。列 [B, D, A, B] もまた LCS の1つである。
最長共通部分列問題は与えられた2つの列 X = [x_1, x_2, ..., x_m] と Y = [y_1, y_2, ..., y_n] の最長共通部分列を求める問題である。本節では動的計画法を用いるとLCS問題が効率よく解けることを示す。
第1段階: 最長共通部分列の特徴づけ
LCS問題を総当りで解くにはXのすべての部分列を生成し、Yの部分列であるかをチェックし、最長のものを記録しておけばよい。
Xの各部分列はXの添字集合 {1, 2, ..., m} のある部分集合に対応するので、2^mの部分列が存在し、指数時間かかる。現実的ではない。
LCS問題は次の定理で示す部分構造最適性を持つ。そして、2つの入力列の「接頭語」の組が部分問題の自然な族であることをすぐに説明する。
正確に定義する。任意の列 X = [x_1, x_2, ..., x_m] とする。各 i = 0, 1, ..., m に対して、 X の i 番目の接頭語を X_i = [x_1, x_2, ..., x_i] と定義する。たとえば、
X = [A, B, C, B, D, A, B] とすると、
X_4 = [A, B, C, B]
X_1 = [A]
X_0 = []
LCSの部分構造最適性 - 定理 15.1
X = [x_1, x_2, ..., x_m]
Y = [y_1, y_2, ..., y_n]
を列、
Z = [z_1, z_2, ..., z_k]
を X と Y の任意の LCS とする。
- x_m = y_n ならば z_k = x_m = y_n であり、 Z_{k-1} は X_{m-1} と Y_{n-1} の LCS である
- x_m \neq y_n のとき、 z_k \neq x_m ならば Z は X_{m-1} と Y の LCS である
- x_m \neq y_n のとき、 z_k \neq y_n ならば Z は X と Y_{n-1} の LCS である
証明
(1) z_k \neq x_m を仮定すると、 Z に x_m = y_n を付け加えて長さ k + 1 の X と Y の共通部分列が構成できる。これは Z が LCS であることに矛盾する。よって、 z_k = x_m = y_n である。
接頭語 Z_{k-1} は X_{m-1} と Y_{n-1} の長さが (k-1) の共通部分列である。これが LCS であることを示したい。長さが k 以上の X_{m-1} と Y_{n-1} の共通部分列 W が存在すると仮定して矛盾を導く。このとき、 W に x_m = y_n を付加すると長さが k+1 以上の X と Y の共通部分列が構成できるが、 Z が長さ k で LCS であることに矛盾する。
(2) z_k \neq x_m ならば Z は X_{m-1} と Y の共通部分列である。長さが k+1 以上の X_{m-1} と Y の共通部分列 W が仮定すると、W は X_m と Y の共通部分列でもあるから、 Z が X と Y の LCS であるという仮定に矛盾する。よって、Z は X_{m-1} と Y の LCS である。
(3) (2) と同様
以上から2つの列の LCS が、その一部として接頭語の LCS を含むことが分かる。したがって、部分構造最適性を持つ。
第2段階: 再帰的な解
定理15.1から列 X と列 Y のLCSを求めるには1個または2個の部分問題を検討する必要がある。
x_m = y_n ならば X_{m-1} と Y_{n-1} の LCS を。
x_m \new y_n ならば X_{m-1} と Y の LCS と X と Y_{n-1} の LCS を。その2つのLCSの長い方がXとYのLCS。
部分問題の最適解がXとYのLCSの中に必ず出現することがわかった。
LCS問題が部分重複性を持つことを確かめる。XとYのLCSを発見するには、XとY_{n-1}のLCSおよびX_{m-1}とYのLCSを発見する必要性が生じるかもしれない。そして、これら2つの部分問題はX_{m-1}とY_{n-1}のLCSを発見するという「部分部分問題」を共有する。
同様にして他の多くの部分問題も部分部分問題を共有している。
列X_iとY_jのLCSのの長さをc[i,j]と定義する。i=0またはj=0ならばLCSの長さ0は明らか。LCS問題の部分構造最適性から漸化式
c[i, j] = \left\{
\begin{array}{ll}
0 & (i = 0 または j = 0)\\
c[i-1, j-1] + 1 & (i, j > 0 かつ x_i = y_i)\\
\max(c[i, j-1], c[i-1, j]) & (i, j >0 かつ x_i \neq y_j)
\end{array}
\right.
を得る。
// 2019/12/26
第3段階: LCSの長さの計算
2つの列のLCSの長さを指数時間で計算する再帰的アルゴリズムは簡単に記述できる。
しかし、異なる部分問題が \Theta (mn) 個しか存在しないことに注目すると、動的計画法を用いて解をボトムアップ的に計算できる。
LCS-LENGTHは2つの列
X = [x_1, x_2, ..., x_m]
Y = [y_1, y_2, ..., y_n]
を入力として受け取る。
表c[0..m, 0..n]
に c[i, j] を保存するが、その計算は行優先で行う。
すなわち、cの第一行を左から右に埋めて、その次に第二行を左から右に・・・と
手続きは最適解を容易に構成するために
表b[1..m, 1..n]
も同時に管理する。
直感的に言うと
b[i, j] は c[i, j] を計算するときに選択した部分問題の最適解に対応する要素を指す。
手続きは表bとcを返す。この時 c[m, n] にはXとYのLCSの長さが格納されている。
LCS-LENGTH(X, Y)
m = X.length
n = Y.length
b[1..m, 1..n] と c[0..m, 0..n] を新しい表とする
for i = 1 to m
c[i, 0] = 0
for j = 0 to n
c[0, j] = 0
for i = 1 to m
for j = 1 to n
if x_i == y_j
c[i, j] = c[i-1, j-1] + 1
b[i, j] = "↖(左上)"
elseif c[i-1, j] >= c[i, j-1]
c[i, j] = c[i-1, j]
b[i, j] = "↑(上)"
else
c[i, j] = c[i, j-1]
b[i, j] = "←(左)"
return c and b
表の各要素の計算に \Theta (1) しかかからないので、この手続きの実行時間は \Theta (mn) である。
X = [A, B, C, B, D, A, B]
Y = [B, D, C, A, B, A]
| y_j | B | D | C | A | B | A | |
|---|---|---|---|---|---|---|---|
| x_i | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| A | 0 | ↑,0 | ↑,0 | ↑,0 | ↖,1 | ←,1 | ↖,1 |
| B | 0 | ↖,1 | ←,1 | ←,1 | ↑,1 | ↖,2 | ←,2 |
| C | 0 | ↑,1 | ↑,1 | ↖,2 | ←,2 | ↑,2 | ↑,2 |
| B | 0 | ↖,1 | ↑,1 | ↑,2 | ↑,2 | ↖,3 | ←,3 |
| D | 0 | ↑,1 | ↖,2 | ↑,2 | ↑,2 | ↑,3 | ↑,3 |
| A | 0 | ↑,1 | ↑,2 | ↑,2 | ↖,3 | ↑,3 | ↖,4 |
| B | 0 | ↖,1 | ↑,2 | ↑,2 | ↑,3 | ↖,4 | ↑,4 |
第4段階: LCSの構成
LCS-LENGTHが出力する表bを用いると、2つの列XとYのLCSを高速に構成できる。b[m, n]から開始し、矢印に従って表を辿るだけでよい。
要素b[i, j]で"↖"にぶつかれば、それはx_i = y_j がLCS-LENGTHが発見したLCSの要素であることを意味している。この方法ではLCSの要素に逆順で出会うことになる。つぎの再帰的手続きはXとYのLCSを正しい順序でprintする。初期呼び出しはPRINT-LCS(b, X, X.length, Y.length)である。
PRINT-LCS(b, X, i, j)
if i == 0 or j == 0
return
if b[i, j] == "↖"
PRINT-LCS(b, X, i-1, j-1)
print x_i
elseif b[i, j] == "↑"
PRINT-LCS(b, X, i-1, j)
else PRINT-LCS(b, X, i, j-1)
各再帰呼び出しの中でiとjのうちの少なくとも一方は減るので、この手続の実行時間はO(m+n)である。
コードの改良
アルゴリズムを開発し終わった後になって、その時間や領域量を改善できることに気づくことがある。ある変更はコードを単純化し定数倍の改良をもたらすが、効率の漸近的な改良に至らない。しかし、時間と領域で漸近的にも十分な改良となる変更もある。
たとえば、LCSアルゴリズムから表bをすべて削除できる。c[i, j]の各要素は表cの他の3つの要素 c[i-1, j-1], c[i-1, j], c[i, j-1]だけに依存する。そこでc[i,j]の値が与えられると表bを調べなくてもこれら3つの値の中からc[i, j]の計算に使われた値をO(1)時間で決定できる
(漸化式15.9の条件を見ればよい)。
したがって、PRINT-LCSと似た手続きを用いてO(m+n)の時間でLCSを再構成できる。
この方法によって領域を \Theta (mn) だけ節役えきるが、表cに Theta (mn) 領域が必要なので、LCSの計算に必要な補助領域は漸近的は変わらない。
しかし、どの時点でも表cの2つの行、計算中の行と直前の行だけが必要であることに気がつけば、LCS-LENGTHの漸近的領域計算量を削減できる。この改善はLCSの長さだけが求められている時に有効である。LCSを再構成する必要がある場合にはこの小さい表は最適解を構成してきた道をO(m+n)時間で辿れるだけの情報を蓄えていない。
練習問題
15.4-1
[1,0,0,1,0,1,0,1]と[0,1,0,1,1,0,1,1,0]のLCSを決定せよ
| y_j | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 1 | 0 | |
|---|---|---|---|---|---|---|---|---|---|---|
| x_i | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 1 | 0 | ↑,0 | ↖,1 | ←,1 | ↖,1 | ↖,1 | ←,1 | ↖,1 | ↖,1 | ←,1 |
| 0 | 0 | ↖,1 | ←,1 | ↖,2 | ←,2 | ←,2 | ↖,2 | ←,2 | ←,2 | ↖,2 |
| 0 | 0 | ↖,1 | ←,1 | ↖,2 | ←,2 | ←,2 | ↖,3 | ←,3 | ←,3 | ↖,3 |
| 1 | 0 | ↑,1 | ↖,2 | ←,2 | ↖,3 | ↖,3 | ←,3 | ↖,4 | ↖,4 | ←,4 |
| 0 | 0 | ↖,1 | ↑,2 | ↖,3 | ←,3 | ←,3 | ↖,4 | ←,4 | ←,4 | ↖,5 |
| 1 | 0 | ↑,1 | ↖,2 | ↑,3 | ↖,4 | ↖,4 | ←,4 | ↖,5 | ↖,5 | ←,5 |
| 0 | 0 | ↖,1 | ↑,2 | ↖,3 | ↑,4 | ←,4 | ↖,5 | ←,5 | ←,5 | ↖,6 |
| 1 | 0 | ↑,1 | ↖,2 | ↑,3 | ↖,4 | ↖,5 | ←,5 | ↖,6 | ↖,6 | ←,6 |
[0,1,0,1,0,1]
15.4-2
完成された表cおよび列X,Yから表bを用いずにO(m+n)時間でXとYのLCSを構成する擬似コードを書け
PRINT-LCS(c, X, Y, i, j)
if c[i, j] == 0
return
if X[i] == Y[i]
PRINT-LCS(c, X, Y, i-1, j-1)
PRINT X[i]
else if c[i-1, j] > c[i, j-1]
PRINT-LCS(c, X, Y, i-1, j)
else PRINT-LCS(c, X, Y, i, j-1)
15.4-3
O(mn)時間で走る履歴管理版のLCS-LENGTH
MEMOIZED-LCS-LENGTH(X, Y, i, j)
if c[i, j] >= 0
return c[i, j]
if i == 0 or j == 0
return c[i, j] = 0
return c[i, j] = x[i] == y[j] ? MEMOIZED-LCS-LENGTH(X, Y, i-1, j-1) + 1 : max(MEMOIZED-LCS-LENGTH(X, Y, i-1, j), MEMOIZED-LCS-LENGTH(X, Y, i, j-1))
15.4-4
表cの2 * min(m, n)個の要素に加えてO(1)の領域を用いてLCSの長さを計算する方法を示せ
短い方をmとする(違うならば転置する)。
前記アルゴリズムはc[i, j]の計算にc[i-1, j-1]とc[i-1, j]とc[i, j-1]しか必要としないので、c[i, j]をc[i % 2, j]に保存し、計算時の参照も同様にすればよい。
min(m, n)個の要素に加えてO(1)の領域を用いて同じことを実行する方法を示せ
同じく短い方をmとする(違うならば転置する)。
前記アルゴリズムはc[i, j]の計算にc[i-1, j-1]とc[i-1, j]とc[i, j-1]しか必要としないので、c[i-1, j-1]からc[i, j-2]までc[i, j]をc[0, j]に保存し、c[i, j-1]とc[i, j]だけx, yに保存して適宜ローリングバッファ的に使っていけばよい。
15.4-5
長さがnの自然数の列から最長増加部分列を発見するO(n^2)時間アルゴリズムを示せ
1
sort(O(nlogn))して、元の列とsortした列のLCS(O(n*n)) => O(n^2)
2
c[i]を先頭からi番目までの最長部分増加列の長さとする。
l[i]を先頭からi番目までの最長部分増加列の最後の数字とする。
c[n] = \left\{
\begin{array}{ll}
0 & (n = 0)\\
\max(c[i] + l[i] < l[n] ? 1 : 0) & (n > i >= 0) ただし、c[i]が同じになる場合はl[i]が小さくなる方を優先する
\end{array}
\right.
| l_i | c_i | LCS | |
|---|---|---|---|
| 8 | 8 | 1 | 8 |
| 1 | 1 | 1 | 1 |
| 3 | 3 | 2 | 1,3 |
| 5 | 5 | 3 | 1,3,5 |
| 4 | 4 | 3 | 1,3,4 |
| 9 | 9 | 4 | 1,3,4,9 |
| 7 | 7 | 4 | 1,3,4,7 |
| 2 | 7 | 4 | 1,3,4,7 |
| 6 | 6 | 4 | 1,3,4,6 |
| 10 | 10 | 5 | 1,3,4,6,10 |
15.4-6
長さがnの自然数の列から最長の単調増加部分列を発見するO(n log n)時間アルゴリズムを示せ