前提知識
文字列関連の定義
文字列 $S$ の subsequence とは、文字列 $S$ から $0$ 個以上 $|S|$ 個以下の文字を取り除いて得られる文字列です。文字列 $S$ と $T$ の LCS (longest common subsequence) とは、$S$ の subsequence でもあり $T$ の subsequence でもある文字列のうち、最も長いものです（複数ある場合はすべてが該当します）。
文字列 $S$ と $T$ の編集距離とは、文字列 $S$ に対して「ある位置の文字の削除」「ある位置の文字の置換」「ある文字の任意位置への追加」を繰り返して $T$ を作るときの操作回数の最小回数です。
文字列 $S$ から $T$ への SES (shortest edit script) とは、文字列 $S$ に対して「ある位置の文字の削除」「ある文字の任意位置への追加」を繰り返して $T$ を作るときの操作列 (edit script) のうち、長さが最小のものです（複数ある場合はすべてが該当します）。編集距離の場合と違い、置換はできないことに注意してください。
アルゴリズム関連
多くのアルゴリズムでは、入力の長さ（や入力の値）で計算量を見積もりますが、別のパラメータで見積もることもあります*1。 特に出力の長さや個数に計算量が依存するアルゴリズムを output-sensitive algorithm と呼ぶことが多いです。列挙系の分野の人が好む傾向にありそうです。
（出力は入力に依るので）出力の長さも入力に関するパラメータのような気はしますが、アルゴリズムを実行した後に得られるパラメータであり、またその分野では出力長を気にすることが多いので、名前が特についているのかなという気がしています。もっと総称して “(implicit) parameter”-sensitive algorithm とか呼びたい文脈がなくもない気もします。
簡単な事実
SES の長さは編集距離の 2 倍で上から抑えられます（編集距離の操作における置換を、削除と追加に変えるとよい）。
LCS は SES と対応する問題であることがわかります。 LCS に使われなかった文字たちを削除・追加する操作列が SES となります。 $S$ と $T$ の LCS の長さを $\lambda(S, T)$、SES の長さを $\delta(S, T)$ と書くと、$|S| + |T| = 2\cdot \lambda(S, T) + \delta(S, T)$ が成り立ちます。 また、$0\le \delta(S, T)\le |S|+|T|$ です。
お詫び
タイトルで「編集距離に対して二乗時間」と言った部分は $\Theta(\delta(S, T)^2)$ time のことでした。紙面の都合でそう書きましたが、定数倍しか違わないことは上記の通りなので、嘘ではないでしょう。
準備
Edit graph
文字列 $S$ から $T$ への edit graph $G = (V, E)$ を定義します。
\[ \gdef\hrz{\text{horizontal}} \gdef\vrt{\text{vertical}} \gdef\diag{\text{diagonal}} \]
$$ \begin{aligned} V &= \{(i, j)\in\Z^2 \mid 0\le i\le |S|, 0\le j\le |T|\}; \\ E &= E_{\hrz} \cup E_{\vrt} \cup E_{\diag}. \end{aligned} $$
ただし、$E$ の各構成要素は次の通りです。
$$ \begin{aligned} E_{\hrz} &= \{( (i, j), (i, j+1)) \in \Z^2\times \Z^2 \mid 0\le i\le |S|, 0\le j\lt |T|\}; \\ E_{\vrt} &= \{( (i, j), (i+1, j)) \in \Z^2\times \Z^2 \mid 0\le i\lt |S|, 0\le j\le |T|\}; \\ E_{\diag} &= \{( (i, j), (i+1, j+1)) \in \Z^2\times \Z^2 \mid 0\le i\lt |S|, 0\le j\lt |T|, S_i = T_j\}. \end{aligned} $$
また、重み関数 $w: E\to E$ を $w(e) = 0$ ($e \in E_{\diag}$), $w(e) = 1$ ($e \notin E_{\diag}$) で定義します。
さらに、頂点 $(i, j)\in V$ のうち $j-i = k$ ($-|S|\le k\le |T|$) を満たす部分集合を $k$-diagonal と呼びます。
$S = \text{“\texttt{coco{}a}”}$, $T = \text{“\texttt{concave}”}$ における edit graph は以下の通りです。いくつかの頂点にはラベルをつけています。また赤の破線は各 $k$-diagonal を表しています。
さて、LCS・SES は edit graph 上での $(0, 0)$ から $(|S|, |T|)$ への最短経路に対応します。 たとえば愚直に 0/1-BFS などを行うことで、$O(|S|\cdot |T|)$ 時間で LCS・SES を求めることができます*2。 たぶん $\Theta( (|S|+|T|)\cdot \delta(S, T))$ 時間とかになると思います。
↑ ここあまり考えてません。参考文献ひとつめの論文の figure 2 にあるアルゴリズムと似たような感じになるのかな〜くらいの気持ちで書きました。
$D$-path
edit graph において $(0, 0)$ を始点とするパスのうち、辺の重みの総和が $D$ のものを $D$-path と呼ぶことにします。
「共通でない文字を選んだ回数が $D$ であるようなパス」「ここまでの edit script の長さが $D$ であるようなパス」などに相当するものです。
$(D, k)$-furthest path
終端 $(i, j)$ が $k$-diagonal にある $D$-path のうち、$i$ が最大のものを $(D, k)$-furthest path と呼ぶことにします。
ここで、$(D, |T|-|S|)$-furthest path の終端が $(|S|, |T|)$ にあるような最小の $D$ が $\delta(S, T)$ となることがわかります。
添字の記法
$\gdef\suffix#1#2{#1[#2\dots]}$
文字列 $S$ の先頭 $i$ 文字を取り除いて得られる文字列を $\suffix{S}{i}$ と書くことにします。
Longest common prefix
$\gdef\lcp{\operatorname{lcp}}$
文字列 $S$, $T$ に共通する接頭辞のうち最長のものの長さを $\lcp(S, T)$ と書きます。 たとえば、$S = \text{“\texttt{echo}”}$, $T = \text{“\texttt{chat}”}$ のとき、$\lcp(\suffix{S}{1}, \suffix{T}{0}) = 2$ です。
考察
まず、$(D, k)$-furthest path の構造を調べます。
$(0, 0)$-furthest path は、$p = \lcp(S, T)$ 回だけ斜めの辺を辿ったパスになります。終端は $(p, p)$ です。
$(D, k)$-furthest path ($D \gt 0$) は、以下の 2 本の furthest な方を選んで得られます。
- $(D-1, k-1)$-furthest path から右向きの辺を 1 回辿って（その行き先を $(i, j)$ とする）、そこから $\lcp(\suffix{S}{i}, \suffix{T}{j})$ 回だけ斜めの辺を辿ったパス
- $(D-1, k+1)$-furthest path から下向きの辺を 1 回辿って（その行き先を $(i, j)$ とする）、そこから $\lcp(\suffix{S}{i}, \suffix{T}{j})$ 回だけ斜めの辺を辿ったパス
証明は帰納的に行えばよいです。
アルゴリズム
方針
そこで、$(D, k)$-furthest path の終端に関する次のような DP を構成できます。
$\gdef\dp{\mathrm{dp}}$
- $\dp[0][0] = \lcp(S, T)$
- for $d$ in $(0, 1, \dots, |S|+|T|)$ do
- for $k$ in $(-d, -d+1, \dots, d-1, d)$ do
- $i \gets \dp[d][k]$
- $j \gets i + k$
- $\dp[d+1][k-1] \xleftarrow{\max} \lcp(\suffix{S}{i + 1}, \suffix{T}{j})$
- $\dp[d+1][k+1] \xleftarrow{\max} \lcp(\suffix{S}{i}, \suffix{T}{j + 1})$
- for $k$ in $(-d, -d+1, \dots, d-1, d)$ do
$S = \text{“\texttt{coco{}a}”}$, $T = \text{“\texttt{concave}”}$ における DP は以下の通りです。 頂点の上の文字は、赤が $S$ のみに含まれる文字（SES における削除）、紫が共通の文字（LCS として選んだ文字）、緑が $T$ のみに含まれる文字（SES における追加）を表します。頂点の下の値は $D$ であり、DP 配列の値ではありません。少しトリッキーに感じるかもしれません。
高速化
ここで、$\lcp(\suffix{S}{i}, \suffix{T}{j})$ の形のクエリに高速に答えたいなとなります。 このクエリは、文字列 $S \concat \bullet \concat T$ に対する LCP array に対する区間最小値クエリに帰着できます。ここで $\concat$ は文字列の連結、$\bullet$ は $S$ にも $T$ にも属さない文字（その 1 文字からなる文字列）を意味します。
LCP array は suffix array を求めていれば線形時間で求められますが、suffix array は SA-IS というアルゴリズムを用いて線形時間で求められます。 また、区間最小値クエリも、線形時間の前処理をすることで定数時間で処理できます。
SA-IS を使う関係で、入力アルファベットはサイズ $n^{O(1)}$ の integer alphabet を仮定します。
アルファベットについて
有限の文字集合 $\Sigma$ を固定し、入力される文字列は $\Sigma$ の要素を任意個つなげたものとします。 この $\Sigma$ を入力アルファベットと呼びます。
アルゴリズムによっては $\Sigma$ が持つ性質を仮定したりします。
$x, y \in \Sigma$ について $x = y$ や $x \neq y$ は定義されているとするのが普通*3だと思いますが、$x \lt y$ や $x \le y$ などは定義されるとは限りません*4。
$x \le y$ が定義されているのを仮定するアルゴリズムとしては「文字列リストの重複を除くために、辞書順ソートして連続する要素を比較する」といったものがありそうです。順序がなければソートができないわけですね。
順序が定義されていないアルファベットとしては、たとえば $\Sigma = \{\dag, \ddag, \ast\}$ や $\Sigma = \{-1, 0, 1, \sqrt{-1}, -\sqrt{1}\}$ などをイメージすればよいでしょうか。 そうした集合でも適当に順序を定義してやればうまくいくことも多そうです。
あるいは、integer alphabet というものもあります。$\Sigma = \{0, 1, \dots, \sigma-1\}$ のような非負整数の集合です。順序が定義されているほか、配列の添字として利用したりする前提のときはこれを仮定しそうです。
いわゆる座圧を行うことで、順序が定義されていれば $O(|\Sigma|\log(|\Sigma|))$ 時間で integer alphabet に帰着できそうです。 帰着に log がつくのがボトルネックにならないケースでは、integer alphabet を仮定しなくてもよさそうという感じです。
また、$k$ のループ範囲は parity を考慮することで $(-d, -d+2, \dots, d-2, d)$ でよいとか、DP 配列を一本で済ませる（これも parity に基づく）などの高速化ができます。
もちろん、$(|S|, |T|)$ に辿りついたらそれ以降の $D$ を調べる必要はないので、そこで打ち切ることができます。 よって（$\Theta(1)$ 時間である）LCP のクエリなどが $\Theta(\delta(S, T)^2)$ 回行われることになります。 前処理は $\Theta(|S|+|T|)$ 時間なので、全体で $\Theta(|S|+|T|+\delta(S, T)^2)$ 時間となります。
復元
$(D, k)$-furthest path を作るときに（下方向ではなく）右方向の辺を用いた $(D, k)$ のリストを保持しておきます*5。右方向の辺を使うたびにリストの末尾に加えていく感じです。
DP が完了したとき、（DP のループ順から）$(D, k)$ の辞書順に昇順ソートされた状態のリストが得られることになります。 復元は「いま見ている $(D, k)$ がこのリストに含まれるか？」を調べながら行うわけですが、このクエリは $(D, k)$ の降順になります。 そこで、このリストを後ろから見て、尺取り法のような要領で log をつけずに行うことができます。
よって、復元まで含めても上記と同じ計算量オーダーで可能です。
実測
普通の DP 提出 #37630508 と、今回紹介した方法 提出 #37630463 を比べます。
$\delta(S, T)$ の最大値は $|S|+|T|$ であり、$|S|\cdot |T| \le (|S|+|T|)^2$ を考えると、最大ケースだけ比較するのはうれしくないです。 そもそも parameter-sensitive algorithm ですので。
そこで、テストケースを調べて得たパラメータと実行時間のリストを載せます。
| ケース名 | $|S|$ | $|T|$ | $\delta(S, T)$ | #37630463 |
#37630508 |
|---|---|---|---|---|---|
0_00 |
5 | 4 | 3 | 11 ms | 9 ms |
0_01 |
5 | 2 | 3 | 5 ms | 6 ms |
0_02 |
1 | 1 | 2 | 6 ms | 5 ms |
0_03 |
12 | 11 | 9 | 7 ms | 7 ms |
1_00 |
1 | 1 | 0 | 7 ms | 5 ms |
1_01 |
1 | 1 | 2 | 8 ms | 6 ms |
1_02 |
3000 | 3000 | 0 | 9 ms | 110 ms |
1_03 |
3000 | 3000 | 6000 | 353 ms | 108 ms |
1_04 |
3000 | 3000 | 5998 | 354 ms | 109 ms |
1_05 |
3000 | 3000 | 5998 | 353 ms | 110 ms |
1_06 |
2985 | 3000 | 15 | 8 ms | 109 ms |
1_07 |
2994 | 2966 | 1150 | 39 ms | 118 ms |
1_08 |
2975 | 2963 | 1726 | 72 ms | 120 ms |
1_09 |
2994 | 2973 | 2083 | 90 ms | 122 ms |
1_10 |
2963 | 2944 | 2579 | 120 ms | 117 ms |
1_11 |
2992 | 2949 | 3045 | 158 ms | 114 ms |
1_12 |
2976 | 2889 | 3373 | 181 ms | 110 ms |
1_13 |
2976 | 2994 | 4046 | 242 ms | 113 ms |
テストケースがおあつらえ向きでびっくりしたのですが、1_06 から 1_13 で徐々に $\delta(S, T)$ が増えているのに応じて #37630463 の方の実行時間も増えているのがわかりますね（二乗っぽいかはちょっと微妙ですが）。1_02 と 1_03 の差も顕著ですね。
入力文字列によっては SA-IS などの処理の重さも変わったりするので、ある程度はブレたりするのかもしれませんね。
おきもち
こうしたアルゴリズムは、たとえば $|S|, |T| \approx 10^6$ かつ $\delta(S, T) \le 10^4$ のような制約があるときにうれしそうです（そんなことある？）
また、「解きたい問題が LCS の特殊ケースに帰着できるが、SES が短くなることがわかっている」のような場合にもうれしそうです（そんなことある？）
でも、転倒数が小さく済むことを利用して高速にソートする のようなものはあるし、ないとは言えないのではないでしょうか。
おためし
$|S| = |T| = 10^6$ かつ $\delta(S, T) = 6\cdot 10^3$ になるような入力を作り、コードテストで複数回試したところ、ブレはあったものの 8393–10034 ms 程度で終わりました。そこそこいけますね。$\delta(S, T) = 5\cdot 10^3$ にすると 5206 ms 程度でした。
下界との関連
アルファベットが $=$, $\neq$ の比較しかできない場合の LCS は $\Omega(n^2)$ の下界が、$<$, $>$ の比較もできる場合の LCS は $\Omega(n\log(n))$ が知られているようです。比較ソートの下界を示すとき同様に決定木の高さを用いて証明するようです。
また、SETH の下では、$\varepsilon\gt 0$ に対して $O(n^{2-\varepsilon})$ 時間で編集距離を求めるアルゴリズムが存在しないことが示せるらしいですが、このときに仮定しているアルファベットがなんなのかはいまいち知りません。
ともかく、今回は最悪 $\Theta(n^2)$ 時間になってしまうため、これらを violate していません。比較どころか integer algorithm を仮定しているのに $\Omega(n\log(n))$ 時間には程遠いですね。もっとよい下界が知られていたりするのでしょうか。
参考文献
- Myers, Eugene W. “An $O(ND)$ difference algorithm and its variations.” Algorithmica 1, no. 1 (1986): 251-266.
- Ullman, J. D., A. V. Aho, and D. S. Hirschberg. “Bounds on the complexity of the longest common subsequence problem.” Journal of the ACM (JACM) 23, no. 1 (1976): 1–12.
- $\Omega(n^2)$ の下界の話が載っている
- Hirschberg, Daniel S. “An information theoretic lower bound for the longest common subsequence problem.” Rice University ECE Technical Report 7705 (1977).
- $\Omega(n\log(n))$ の下界の話が載っている
ひとつめの Myers の文献が今回の話です。当時は SA-IS などがなかったためか、「suffix tree などの重いデータ構造を使えば $O(N\log(N)+D^2)$ 時間のアルゴリズムが作れるが practical ではない」のような旨が書かれていました。今回の感じだと、だいぶ practical になっているのではないでしょうか。
ふたつめの Ullman の文献内の [9] によれば、最悪 $O(n^2\log(n))$ 時間になるが、諸々の応用で出てくる多くのケースでは $O(n\log(n))$ 時間で動作する、実装も簡単な LCS のアルゴリズムが存在するらしいです。調べようとしたら Unpublished memorandum と書かれていて参ってしまいました。
おわり
にゃんこ。
