最長共通部分列問題 (Longest Common Subsequence)

部分列 (Subsequence) は系列のいくつかの要素を取り出してできた系列のことです。二つの系列の共通の部分列を共通部分列 (Common Subsecuence)と言います。共通部分列のうち、もっとも長いものを最長共通部分列 (Longest Common Subsequence, LCS) と言います。

X = <A, B, C, B, D, A, B>
Y = <B, D, C, A, B, A>

という二つの系列から得られる LCS は

<B, C, B, A>

で、その長さは 4 です。長さ 2 の<B, D> の長さ 3 の <A, B, A> なども共通部分列ですが、最長ではないのでこれらは LCS ではありません。また、LCS は最長であれば位置はどこでも良いので、この場合 <B, D, A, B> も LCS です。 LCS は動的計画法 (Dynamic Programmming, DP) を利用すると効率良く求めることができることが知られています。LCS を解説したウェブ上の文書では Common Subsequence 解説 (www.deqnotes.net) が分かりやすかったです。

アルゴリズムイントロダクションの第2巻第15章は DP の解説ですが、その中で LCS の話があります。はじめに読んだときはあまり理解できなかったのですが、編集距離を求める問題によく似ていることやプログラミングコンテストでよく LCS を基礎にした問題が出るというような話を聞いて興味が沸いたので、改めて整理してみました。アルゴリズムイントロダクション第15章の要約は、id:nitoyon さんが作成された神資料 (slideshare) がありますので、興味のある方はこちらもどうぞ。

LCS のキーポイント

LCS の問題では <B, C, B, A> のように具体的な最長共通部分系列を求める場合と、最長共通部分系列の長さを求める場合があります。他の DP の問題でも良くあるように、後者がわかれば前者はその解から構成できることが分かっているので、まずは最長部分系列長 (LCS長) を求めることから考えます。

LCS 長を DP で求めるにあたって鍵になるのは、二つの系列 X, Y があったとき i 番目までと j番目までの LCS 長、つまり X_i と Y_i のLCS 長が、より小さな系列の LCS 長から求められるのを上手く利用する点です。X_i と Y_j の LCS 長 LCS(i, j) はより小さな部分問題

LCS(i -1, j -1)
LCS(i - 1, j)
LCS(i, j - 1)

の3つが分かっていれば求められるという関係 (ある問題の最適解 LCS(i, j) は、その部分問題の最適解である LCS(i -1, j -1), LCS(i - 1, j), LCS(i - 1, j) から導かれるという部分問題最適性) から再帰的な漸化式が得られるので、あとは DP のパターンに当てはめて解けば良いのです。

LCS の部分問題最適性は以下のように例を考えるとわかります。

(i) 一番後ろの要素が同じである場合

例として

X_i = <..., A>
Y_j = <..., A>

だったとします。

そしてX, Y の最後の要素 (i番目, j 番目) の手前、つまり ... で表されている系列の LCS 長 LCS(i -1, j - 1) が既に分かっているとしましょう。上記の例では最後の要素は X, Y ともに "A" で同一です。この場合明らかに LCS 長 LCS(i, j) は LCS(i - 1, j - 1) +1 です。

これは例えば

X_i = <A, B, C, B, D, A>
Y_j = <B, D, C, A, B, A>

のようなケースです。LCS(i -1, j -1) は最後の要素 "A" を除いた X_{i - 1} = <A, B, C, B, D> と Y_{j - 1} = <B, D, C, A, B> の LCS <B, C, B> の長さなので 3 です。この LCS(i -1, j - 1) = 3 が分かっているときに X_i と Y_j の最後の要素が等しかったわけですから、当然 LCS(i, j) は 3 + 1 = 4 になるでしょう。

この例では X と Y の要素数が等しくなっていますが、そうでなくても同じことが言えます。

X = <A, B, C, A>
Y = <A, B, A>

という系列だったとすると、最後の要素 "A" を除いたとの LCS 長は 2 です。除いていた最後の要素は二つの系列で等しいので明らかに i - 1, j - 1 から LCS は 1 増えます。すなわち LCS(i, j) = LCS(i - 1, j - 1) + 1 です。

(ii) 一番後ろの要素同士が異なっていた場合

例を変えて

X_i = <..., A>
Y_j = <..., B>

だったとします。最後の要素は X, Y で異なるので先ほどのように既知の LCS(i - 1, j - 1) から ... というわけにはいきませんが、ほかに LCS(i - 1, j) と LCS(j - 1, i) も既知だったとすると、LCS 長 LCS(i, j) は LCS(i - 1, j) もしくは LCS(i, j - 1) のいずれかで大きい方になります。自分はここを理解するのに手間取ったのですが、以下のような例を考えたところ理解できました。

例えば

X_i = <A, B, A, B, D>
Y_j = <B, A, B>

のようなケースです。LCS は <B, A, B> です。このとき X か Y の末尾の要素を削ってみます。X の末尾を削ったときは LCS には影響がないのに対し、Y の末尾を削ると LCS は <B, A> になってしまいます。この場合 LCS(i, j) = LCS(i - 1, j) が成り立ちます。つまり LCS(i - 1, j) > LCS(i, j - 1) のときには LCS(i, j) = LCS(i - 1, j) です。

この例で X, Y は適当に選んでいるので、逆のケース、X の末尾を削ったら LCS に影響があるけど Y の末尾を削っても影響はない ... というのは同じように考えることができます。(X を削っても Y を削っても LCS に影響が出てしまうということはあり得ません。それは要素が等しい場合すなわち (i) のケースです) 結局 LCS(i, j) = max{LCS(i - 1, j), LCS(i, j - 1)} です。

ほかの例として

X_i = <A, B, A, B, D>
Y_j = <B, A, B, C>

という場合では X と Y の末尾の要素をいずれかを削っても LCS であるには影響がありません。この場合は LCS(i, j) = LCS(i - 1, j) でも LCS(i, j -1) でもどちらでも同じです。

LCS(i, j) はひとつ前の要素の LCS との関係で決まる

こうして結局 LCS(i, j) は

(i) のケースの場合
- LCS(i - 1, j - 1) + 1
(ii) のケースの場合
- LCS(i, j - 1) か LCS(i - 1, j) の大きい方

なので、

LCS(i - 1, j - 1)
LCS(i, j - 1)
LCS(i - 1, j)

の 3 つが分かっていれば求まることが分かります。

LCS(i, j) の i, j は任意ですから、LCS(i - 1, j - 1) や LCS(i, j -1), LCS(i - 1, j) もまたより小さな LCS(i', j') から求まります、つまり LCS 長を求める問題は部分問題の最適解を再帰的に求めていくことで計算できるわけです。なお LCS(0, 0), LCS(i, 0), LCS(0, j) はいずれかの系列の要素が空なのでいずれもLCS長が 0 です。

まとめると X_i = < x₁, x₂, ... x_i >, Y_j = < y₁, y₂, ... y_j > のとき

i = 0 または j = 0 のとき: LCS(i, j) = 0
i, j > 0 かつ x_i = y_j のとき: LCS(i, j) = LCS(i - 1, j - 1) + 1
i, j > 0 かつ x_i ≠ y_j のとき: LCS(i, j) = max( LCS(i - 1, j), LCS(i, j - 1) )

です。これは再帰的な漸化式です。

LCS 長を求めるには i 行 j 列の2次元表を用意してその要素を各々の LCS長で埋めて行けば良いのですが、LCS(i, j) は二つの系列における一つ前要素との関係だけで求まるので

LCS(1, 1) を LCS(0, 0), LCS(0, 1), LCS(1, 0) から求める
LCS(1, 2) を LCS(0, 1), LCS(0, 2), LCS(1, 1) から求める
LCS(1, 3) を LCS(0, 2), LCS(0, 3), LCS(1, 2) から求める
...
LCS(1, j) を LCS(0, j -1), LCS(0, j), LCS(1, j - 1) から求める

LCS(2, 1) を LCS(1, 0), LCS(1, 1), LCS(2, 0) から求める
LCS(2, 2) を LCS(1, 1), LCS(1, 2), LCS(2, 1) から求める
...
LCS(i, j) を LCS(i - 1, j - 1), LCS(i - 1, j), LCS(i, j -1) から求める

という風に小さな要素から順に計算していくことで表を埋めることができます。i または j が 0 の値は既知で、そこを起点に i, j の昇順に解を求めていくと、常に既知の値だけで計算ができるところがポイントです。これは DP そのものです。

このような再帰アルゴリズムにおいてある問題から生じる部分問題の総数が少ない (多項式程度) 場合は部分問題重複性が高いと言えます。LCS はその部分問題が非常に限定的ですから、そういう意味でも DP が有効であると言えます。

LCS 問題を Perl で実装

LCS 長(の表)を DP で求めるコードを Perl で書いてみました。

#!/usr/bin/env perl
# Longest common subsequence
use strict;
use warnings;
use Perl6::Say;
use List::Util qw/max/;

my $s1 = shift;
my $s2 = shift or die "usage: $0 <string1> <string2>";

my $a = [ unpack('C*', $s1) ];
my $b = [ unpack('C*', $s2) ];

my $lcs = [];
for my $i (0..@$a) {
    $lcs->[$i]->[0] = 0;
}

for my $j (0..@$b) {
    $lcs->[0]->[$j] = 0;
}

for my $i (1..@$a) {
    for my $j (1..@$b) {
        if ($a->[$i - 1] == $b->[$j - 1]) {
            $lcs->[$i]->[$j] = $lcs->[$i - 1]->[$j - 1] + 1;
        } else {
            $lcs->[$i]->[$j] = max(
                $lcs->[$i]->[$j - 1],
                $lcs->[$i - 1]->[$j],
            );
        }
    }
}

show_lcs_table($a, $b, $lcs);

sub show_lcs_table {
    my ($a, $b, $lcs) = @_;
    say sprintf "    %s", join ' ', map { chr } @$b;

    for (my $i = 0; $i < @$lcs; $i++) {
        printf
            "%s %s\n",
            $i > 0 ? chr $a->[$i - 1] : ' ',
            join ' ', @{$lcs->[$i]};
    }
}

コマンドライン引数から二つの文字列 (系列) を受け取り、その LCS 長を計算した二次元表を出力します。結果は以下のようになります。

% perl lcs.pl ABCBDAB BDCABA
    B D C A B A
  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

ABCBDAB と BDCABA の LCS 長は一番右下の要素に格納されますので、4 ということが分かります。

コードの核になるのは二重ループの中の

for my $i (1..@$a) {
    for my $j (1..@$b) {
        if ($a->[$i - 1] == $b->[$j - 1]) {
            $lcs->[$i]->[$j] = $lcs->[$i - 1]->[$j - 1] + 1;
        } else {
            $lcs->[$i]->[$j] = max(
                $lcs->[$i]->[$j - 1],
                $lcs->[$i - 1]->[$j],
            );
        }
    }
}

の箇所です。二つの系列の最後の要素を比較 (系列は 0 が空を表すので、0 .. n になりますが配列は 0 .. n - 1 なので、比較は $a->[$i] == $b->[$j] ではなく $a->[$i - 1] == $b->[$j - 1] になる点に注意) して同じであれば LCS(i - 1, j - 1) + 1, そうでなければ max を取って LCS(i, j) を求めています。先に見た条件そのままです。

なお、LCS(i - 1, j - 1) が LCS(i - 1, j), LCS(i, j - 1) より大きくなることはあり得ませんので

for my $i (1..@$a) {
    for my $j (1..@$b) {
        my $match = ($a->[$i - 1] == $b->[$j - 1]) ? 1 : 0;

        $lcs->[$i]->[$j] = max(
            $lcs->[$i - 1]->[$j - 1] + $match,
            $lcs->[$i]->[$j - 1],
            $lcs->[$i - 1]->[$j],
        );
    }
}

と書くこともできます。

計算量は m x n の二次元表を埋めるだけ、その表の各セル LCS(m, n) は部分問題3つの max だけで求まるので定数とみなせるため結局 Θ(mn) です。

LCS の構成

LCS 長ではなく LCS を求めるには以下の関数を実装します。引数にふたつの系列、先に求めた LCS 長の表、求めたい LCS の i と j を渡します。

sub print_lcs {
    my ($a, $b, $lcs, $i, $j) = @_;
    if ($i == 0 or $j == 0) {
        return;
    }

    if ($a->[$i - 1] == $b->[$j - 1]) {
        print_lcs($a, $b, $lcs, $i - 1, $j - 1);
        print chr $a->[$i - 1];
    } else {
        if ($lcs->[$i - 1]->[$j] >= $lcs->[$i]->[$j - 1]) {
            print_lcs($a, $b, $lcs, $i - 1, $j);
        } else {
            print_lcs($a, $b, $lcs, $i, $j - 1);
        }
    }
}

LCS(i, j) は例によって一つ前の LCS 長にのみ依存しているので、LCS(i, j) が分かれば $a->[$i - 1] や $b->[$j - 1] の要素が LCS の要素であるか否か、また次にどこの要素を見ればよいかが分かる、というのを利用して LCS(i, j) から逆に遡って行って LCS を求めます。この辺りは http://en.wikipedia.org/wiki/Longest_common_subsequence_problem を見た方が早いでしょう。

print_lcs() の結果は以下のようになります。

% perl lcs_02.pl ABCBDAB BDCABA
BCBA

めでたしめでたし。

アルゴリズムイントロダクションはどの章も最初は難しいというかとっつきにくいなと感じるのですが、後からよくよく読んでみると、これ以上にないくらい簡潔かつ分かりやすく書いてるようにも感じるという不思議な本です。本を読む本 (講談社学術文庫) によれば「後から読み返すたび新たな発見がある書籍こそが名著」だそうですから、こういう本こそが良著であるとも言えそうです。

追記: Python で実装

この手の二次元配列を使うようなアルゴリズム的なコードは Perl ですと記号やブロックが多くなりがちなため、Python でも書いてみました。

#!/usr/bin/pyton
import sys

def print_lcs (a, b, lcs, i, j):
    if (i == 0 or j == 0):
        return
    if (a[i-1] == b[j-1]):
        print_lcs(a, b, lcs, i - 1, j - 1)
        print a[i - 1],
    else:
        if lcs[i-1][j] >= lcs[i][j-1]:
            print_lcs(a, b, lcs, i - 1, j)
        else:
            print_lcs(a, b, lcs, i, j - 1)

a = sys.argv[1]
b = sys.argv[2]

lcs = [ [0 for j in range(len(b) + 1)] for i in range(len(a)+1) ]

for i in xrange(1, len(a) + 1):
    for j in xrange(1, len(b) + 1):
        if a[i - 1] == b[j - 1]:
            x = 1
        else:
            x = 0
        lcs[i][j] = max(lcs[i-1][j-1] + x, lcs[i-1][j], lcs[i][j-1])

print_lcs(a, b, lcs, len(a), len(b))

% python lcs.py XMJYAUZ MZJAWXU
M J A U

慣れないもので、おかしなところがあるかも。

参考文献

Charles E. Leiserson, Ronald L. Rivest, Clifford Stein, Thomas H. Cormen (共著) 浅野哲夫, 岩野和正, 梅尾博司, 山下雅史, 和田幸一 (共訳) 『アルゴリズムイントロダクション改訂2版第2巻アルゴリズムの設計と解析手法』, 近代科学者, 2007
Kenji Inoue『Common Subsequence 解説』,http://www.deqnotes.net/acmicpc/1458/ , 2006