C++: 編集距離を求めるアルゴリズム

編集距離(edit distance)とは二つの文字列がどの程度異なっているかを示す数値であり、レーベンシュタイン距離(Levenshtein distance)を指すことが多い。文字の挿入、削除、置換それぞれを一つの操作として必要な操作の最小数を求めるものだ。例えば、kittenとsittingの編集距離を求める場合、下記のように3回の操作でkittenをsittingに変更できるので編集距離は3となる。

1. sitten (k を s に置換） 2. sittin (e を i に置換） 3. sitting (g を挿入)

そこで今回は編集距離を求める複数のアルゴリズムについてC++で実装してみた。

動的計画法

編集距離を求めるもっとも一般的なアルゴリズムは、動的計画法(dynamic programming)だろう。計算時間はO(mn)であり、手軽だ。C++で書いたコードを下に示す。尚、コード中に出てくる定数SIZEは扱う文字列にあわせて適当に決めておく。動的にメモリを確保することもできるが処理が遅くなるので、あらかじめ決め打ちしておく方が効率がよい。

int edit_distance_dp(const string& str1, const string& str2) { static int d[SIZE][SIZE]; for (int i = 0; i < str1.size() + 1; i++) d[i][0] = i; for (int i = 0; i < str2.size() + 1; i++) d[0][i] = i; for (int i = 1; i < str1.size() + 1; i++) for (int j = 1; j < str2.size() + 1; j++) d[i][j] = min(min(d[i-1][j], d[i][j-1]) + 1, d[i-1][j-1] + (str1[i-1] == str2[j-1] ? 0 : 1)); return d[str1.size()][str2.size()]; }

O(ND)アルゴリズム

追記(2009/7/8): コメントで指摘を受けたので、O(ND)/O(NP)アルゴリズムについて調べたところ、挿入と削除のみを使った操作の最小数(論文中ではshortest edit "script")を求めるアルゴリズムだった。つまり、置換は削除・挿入の2手順になる。上記の動的計画法の三項演算子の部分を、(str1[i-1] == str2[j-1] ? 0 : 1)から(str1[i-1] == str2[j-1] ? 0 : 2)と変更すれば同じ操作数を求めることができる。

次に、O(ND)アルゴリズムを示す。Nは2つの要素数(m, n)の合計、Dは編集距離になる。つまり、編集距離が少なければ少ないほど高速に計算されるわけだ。原著論文は、E. W. Myers. An O(ND) difference algorithm and its variations. Algorithmica, 1, 251-266 (1986)で、日本語による解説では、後述のO(NP)アルゴリズムとともに、文書比較アルゴリズムが参考になる。以下にコードを示す。

int edit_distance_ond(const string& str1, const string& str2) { static int V[SIZE]; int x, y; int offset = str1.size(); V[offset + 1] = 0; for (int D = 0; D <= str1.size() + str2.size(); D++) { for (int k = -D; k <= D; k += 2) { if (k == -D || k != D && V[k-1+offset] < V[k+1+offset]) x = V[k+1+offset]; else x = V[k-1+offset] + 1; y = x - k; while (x < str1.size() && y < str2.size() && str1[x] == str2[y]) { x++; y++; } V[k+offset] = x; if (x >= str1.size() && y >= str2.size()) return D; } } return -1; }

O(NP)アルゴリズム

O(ND)をさらに改良したアルゴリズムが、O(NP)アルゴリズムとなる。Pは、m, n (m >= n)をそれぞれ要素数とした場合、P=(D-(m-n))/2で表される。経験的にはO(ND)の半分ほどになるらしい。原著論文は、S. Wu, U. Manber, and G. Myers. An O(NP) sequence comparison algorithm. Information Processing Letter. 35(6) 317-332 (1990)だ。以下にコードを示す。

int snake(int k, int y, const string& str1, const string& str2) { int x = y - k; while (x < str1.size() && y < str2.size() && str1[x] == str2[y]) { x++; y++; } return y; } int edit_distance_onp(const string& str1, const string& str2) { // required: s1->size() <= s2->size() const string* const s1 = str1.size() > str2.size() ? &str2 : &str1; const string* const s2 = str1.size() > str2.size() ? &str1 : &str2; static int fp[SIZE]; int x, y, k, p; int offset = s1->size() + 1; int delta = s2->size() - s1->size(); for (int i = 0; i < SIZE; i++) fp[i] = -1; for (p = 0; fp[delta + offset] != s2->size(); p++) { for(k = -p; k < delta; k++) fp[k + offset] = snake(k, max(fp[k-1+offset] + 1, fp[k+1+offset]), *s1, *s2); for(k = delta + p; k > delta; k--) fp[k + offset] = snake(k, max(fp[k-1+offset] + 1, fp[k+1+offset]), *s1, *s2); fp[delta + offset] = snake(delta, max(fp[delta-1+offset] + 1, fp[delta+1+offset]), *s1, *s2); } return delta + (p - 1) * 2; }

ビットパラレル法

最後にビットパラレル(bit-parallel)法を示す。これは動的計画法を応用しており、ビット演算を利用することで並列処理が可能となり、O(m/w n)で計算ができる。ここでwはワード長を示す。つまり、要素数mがワード長wと同じかそれ以下ならO(n)で計算ができることになり、非常に高速に処理ができる。原著論文は、G. Myers. A fast bit-vector algorithm for approximate string matching based on dynamic progamming. JACM, 46(3) 395-415, (1999)である。以下にコードを示す。

template<typename T> int edit_distance_bit(const string& str1, const string& str2) { char s1[sizeof(T) * 8] = { ' ' }; char s2[sizeof(T) * 8] = { ' ' }; char *p1, *p2; // required: str1.size() >= str2.size() if (str1.size() >= str2.size()) { p1 = s1; p2 = s2; } else { p1 = s2; p2 = s1; } copy(str1.begin(), str1.end(), p1 + 1); copy(str2.begin(), str2.end(), p2 + 1); int m = strlen(s1); int n = strlen(s2); const T ONE = 1; T Peq[256] = { 0 }; T Pv, Eq, Xv, Xh, Ph, Mh; T Mv = 0; int Score = m; for (int i = 0; i < m; i++) { Peq[s1[i]] |= ONE << i; Pv |= (ONE << i); } for (int j = 0; j < n; j++) { Eq = Peq[s2[j]]; Xv = Eq | Mv; Xh = (((Eq & Pv) + Pv) ^ Pv) | Eq; Ph = Mv | ~(Xh | Pv); Mh = Pv & Xh; if (Ph & (ONE << (m - 1))) Score++; else if (Mh & (ONE << (m - 1))) Score--; Ph <<= ONE; Pv = (Mh << ONE) | ~(Xv | Ph); Mv = Ph & Xv; } return Score; }

ついでに、C++のbitsetを利用して、何ビットでも演算できるようにしてみた。bitsetではほとんどの場合、通常のビット演算が可能なので、コードは上記のコードとあまり変わらない。ただ、足し算の演算が用意されていないので、演算子オーバーロードでoperator+()を用意した。もっとも、bitsetの利用はビット演算の高速処理が失われるのであまり意味はないかもしれない。コードは以下の通り。

template<size_t N> const bitset<N> operator+(const bitset<N>& lhs, const bitset<N>& rhs) { bitset<N> a(lhs), b(rhs), ret(lhs ^ rhs); for (b = (a & b) << 1, a = ret; b.any(); b = (a & b) << 1, a = ret) ret ^= b; return ret; } template<size_t N> int edit_distance_bitset(const string& str1, const string& str2) { char s1[N] = { ' ' }; char s2[N] = { ' ' }; char *p1, *p2; // required: str1.size() >= str2.size() if (str1.size() >= str2.size()) { p1 = s1; p2 = s2; } else { p1 = s2; p2 = s1; } copy(str1.begin(), str1.end(), p1 + 1); copy(str2.begin(), str2.end(), p2 + 1); int m = strlen(s1); int n = strlen(s2); const bitset<N> ONE(1); bitset<N> Peq[256]; bitset<N> Pv, Mv, Eq, Xv, Xh, Ph, Mh; int Score = m; for (int i = 0; i < m; i++) { Peq[s1[i]] |= ONE << i; Pv |= (ONE << i); } for (int j = 0; j < n; j++) { Eq = Peq[s2[j]]; Xv = Eq | Mv; Xh = (((Eq & Pv) + Pv) ^ Pv) | Eq; Ph = Mv | ~(Xh | Pv); Mh = Pv & Xh; if ((Ph & (ONE << (m - 1))).any()) Score++; else if ((Mh & (ONE << (m - 1))).any()) Score--; Ph <<= 1; Pv = (Mh << 1) | ~(Xv | Ph); Mv = Ph & Xv; } return Score; }

速度比較

さて、これらのアルゴリズムがどれほどの速度で、どれだけの差があるのか、実際に計測してみた。環境は、Windows XP、Intel Core2 6600@2.4GHz、RAM 2GBで、"cl edist_test.cpp /EHsc /O2"としてコンパイルした。結果は以下の通りで、ビット長以下の文字列長ではやはりビットパラレル法が最も速い。今回の結果では動的計画法の20分の1以下の時間で計算が完了した。次いでO(NP)となり、O(ND)、動的計画法の順で遅くなった。遅いと云えばbitsetを利用したビットパラレル法も遅かったがそれでも動的計画法よりは速かった。また、文字列長やその組み合わせによって速度は大きく変わるので、今回の結果はそれらの一例として考えてもらいたい。因みに、文字列はDNAの塩基配列(acgt)となっている。

文字列1: agtcaaaagtcagtcagtcagtcagtcacagtcagaaggcatccaaccga 文字列2: ccgttagtcagaaacagtcagtcagtcagtcagtccagtcttaggcccgga 動的計画法: 10万回の繰り返しで2.859秒 O(ND)アルゴリズム: 10万回の繰り返しで0.500秒 O(NP)アルゴリズム: 10万回の繰り返しで0.359秒 ビットパラレル法(unsigned long long : 64ビット): 10万回の繰り返しで0.125秒 ビットパラレル法(bitset : 60ビット): 10万回の繰り返しで1.281秒

参考までに、測定に用いたコードを示しておく。

void run(int (*func)(const string&, const string&), const string& arg1, const string& arg2, int num, const string& msg) { clock_t start, finish; start = clock(); for (int i = 0; i < num - 1; i++) (*func)(arg1, arg2); cout << msg << " : " << (*func)(arg1, arg2) << endl; finish = clock(); cout << "Time: " << (double)(finish - start) / CLOCKS_PER_SEC << "s (" << num << " times)" << endl; cout << endl; } int main() { string str1 = "agtcaaaagtcagtcagtcagtcagtcacagtcagaaggcatccaaccga"; string str2 = "ccgttagtcagaaacagtcagtcagtcagtcagtccagtcttaggcccgga"; cout << str1 << endl; cout << str2 << endl; run(edit_distance_dp, str1, str2, 100000, "dp "); run(edit_distance_ond, str1, str2, 100000, "ond "); run(edit_distance_onp, str1, str2, 100000, "onp "); run(edit_distance_bit<unsigned long long>, str1, str2, 100000, "bit "); run(edit_distance_bitset<60>, str1, str2, 100000, "bitset"); return 0; }

ところで、アルゴリズムによって編集距離が微妙に違ってくるが、これは仕方がないことなのかなぁ。
上でも書いたが、一操作を「挿入、削除、置換」とするか(動的計画法、ビットパラレル法)、「挿入、削除」にするか(O(ND)/O(NP)アルゴリズム)の違いだった。いくつかのバグと一緒に修正しておいた。