第3章

鞍点探索を改良する

(やさしい版) Pearls of Functional Algorithm Design(関数プログラミングによるアルゴリズム設計の真珠)

どんな問題?

関数 f(x, y) と数 z が与えられたとき、f(x, y) = z となるすべての組 (x, y) を見つけるのがこの章の問題です。

ただし前提として、f は自然数の組から自然数への関数で、各引数について「厳密に増える」ものとします(xy を大きくすると値も必ず大きくなる)。それ以上のことは何も仮定できません。

ポイント f の計算コストは非常に大きいかもしれないので、目標は「f を評価する回数をできるだけ少なくする」ことです。

この章では、先生と 4 人の生徒(アン、ジャック、メアリー、テオ)の対話形式で、探索アルゴリズムを段階的に改良していきます。

ステップ1:素朴な探索(ジャックの案)

まずは全部試す

f は増加関数なので、f(x, y) = z なら必ず xz かつ yz。だから (z+1) × (z+1) の範囲を全部見ればよいです。

invert f z = [(x, y) | x ← [0 .. z], y ← [0 .. z], f (x, y) == z]
Dart typedef Fn = int Function(int x, int y); List<(int, int)> invertNaive(Fn f, int z) { final results = <(int, int)>[]; for (var x = 0; x <= z; x++) { for (var y = 0; y <= z; y++) { if (f(x, y) == z) results.add((x, y)); } } return results; }
問題点 f の評価が (z+1)2もかかる。f のコストが大きいと現実的でない。

対角線より下だけを見る(テオの改良)

f(x, y) ≥ x + y が成り立つので、探索を対角線以下に限定できます:

invert f z = [(x, y) | x ← [0 .. z], y ← [0 .. z − x], f (x, y) == z]
Dart List<(int, int)> invertDiagonal(Fn f, int z) { final results = <(int, int)>[]; for (var x = 0; x <= z; x++) { for (var y = 0; y <= z - x; y++) { if (f(x, y) == z) results.add((x, y)); } } return results; }

これで評価回数はおよそ半分に。さらに上限を z − f(0, 0) 系に置き換えれば、z < f(0, 0) の場合は即終了できます。

ステップ2:鞍点探索(アンの案)

左上の角から始める

正方形の左上の角 (0, z) から始めるのがミソです。探索中は、常に左上 (u, v) と右下 (z, 0) を対角とする長方形の中を探します。

(0, z) ─────────── (z, z)
   │                │
   │  (u, v) ┌──┐   │
   │         │  │   │
   │         └──┘   │
(0, 0) ─────────── (z, 0)

3 通りの場合分け

角 (u, v) で f(u, v) を計算して、z と比べます:

invert f z = find (0, z) f z find (u, v) f z | u > z ∨ v < 0 = [] | z' < z = find (u+1, v) f z | z' == z = (u, v) : find (u+1, v−1) f z | z' > z = find (u, v−1) f z where z' = f (u, v)
Dart List<(int, int)> invertSaddleback(Fn f, int z) { final results = <(int, int)>[]; var u = 0, v = z; while (u <= z && v >= 0) { final zp = f(u, v); if (zp < z) { u++; // 列 u の残りは全て z 未満なので捨てる } else if (zp > z) { v--; // 行 v の残りは全て z より大きいので捨てる } else { results.add((u, v)); u++; v--; } } return results; }
評価回数 最悪でも 2z+1 回、最良で z+1 回。二次から線形に大幅改善!

二分探索で範囲を絞る(テオの追加)

そもそも (z+1)×(z+1) の正方形は大きすぎるので、まず境界 m, n を二分探索で求めます:

m = bsearch (λy → f (0, y)) (−1, z + 1) z n = bsearch (λx → f (x, 0)) (−1, z + 1) z
Dart // 架空値 f(0, -1) = 0, f(-1, 0) = 0 で f を拡張してから int fExt(Fn f, int x, int y) => (x < 0 || y < 0) ? 0 : f(x, y); int mOf(Fn f, int z) => bsearch((y) => fExt(f, 0, y), -1, z + 1, z); int nOf(Fn f, int z) => bsearch((x) => fExt(f, x, 0), -1, z + 1, z);

そして (m+1)×(n+1) の長方形だけを探索。二分探索の実装はこう:

bsearch g (a, b) z | a+1 == b = a | g m ≤ z = bsearch g (m, b) z | otherwise = bsearch g (a, m) z where m = (a + b) div 2
Dart // g(m) ≤ z < g(m+1) となる m を返す。呼び出し側は g(a) ≤ z < g(b) を保証。 int bsearch(int Function(int) g, int a, int b, int z) { while (a + 1 < b) { final m = (a + b) ~/ 2; if (g(m) <= z) { a = m; } else { b = m; } } return a; }

この版では、f の評価は最悪 2 log z + m + n 回。たとえば f(x, y) = 2x + 3y のように m, nz よりずっと小さければ、O(log z) で済みます。

豆知識 この戦略は 鞍点探索(saddleback search)と呼ばれ、David Gries が名付けました。f の三次元プロットが左下最小・右上最大・両翼のような「鞍」の形に見えることから来ています。

ステップ3:分割統治法(メアリーの提案)

中央の要素で 2 つに分ける

「二分探索の二次元版」を考えたらどうか?というのがメアリーのアイデアです。長方形 (u, v)〜(r, s) について、中央の (p, q) を調べます:

問題点 残るのは L 字形で、単一の長方形ではない。関数型プログラマは再帰で対応できるので、これは大きな問題ではない。

計算量を比べてみる(ジャック vs テオ)

T(m, n) を m × n の長方形を探索する評価回数とします。

テオの案(3 つの長方形に分ける):

T(m, n) = 1 + T(⌈m/2⌉, ⌊n/2⌋) + T(⌈m/2⌉, ⌈n/2⌉) + T(⌊m/2⌋, ⌈n/2⌉)

解くと T(m, n) ≤ m1.59 log(2n/m)。

ジャックの案(2 つの長方形に分ける、水平カット):

T(m, n) = 1 + T(⌊m/2⌋, ⌈n/2⌉) + T(⌈m/2⌉, n)

解くと T(m, n) ≤ m log(2n/√m)。mn なら水平、mn なら垂直カット。どちらも鞍点探索より速い。

ステップ4:下界と最適解

下界を求める(アンの分析)

m × n の長方形で可能な答え(解のパターン)の数を A(m, n) とすると、これは階段状経路の数と一致し、二項係数を使って表せます:

A(m, n) = Σk=0..m C(m, k) · C(n, k) = C(m + n, n)
導出のヒント ヴァンデルモンドの畳み込み(Vandermonde's convolution)を使うと C(m+n, n) にまとまります。

各比較には 3 通りの結果があるので、三分木の高さは log3 A(m, n) 以上。対数を取ると:

log A(m, n) = Ω(m log(1 + n/m) + n log(1 + m/n))

mn のときは Ω(m log(n/m)) が下界。ジャックの解はまだこの下界を達成していません。

最適な解(メアリーの最終案)

メアリーは 中央の行に二分探索を使う戦略を提案します。mn(列の方が多い)と仮定して、中央の行 q = (v+s) div 2 で二分探索し、f(p, q) ≤ z < f(p+1, q) となる p を求める:

find (u, v) (r, s) f z | u > r ∨ v < s = [] | v−s ≤ r−u = rfind (bsearch (λx → f (x, q)) (u−1, r+1) z) | otherwise = cfind (bsearch (λy → f (p, y)) (s−1, v+1) z) where p = (u+r) div 2 q = (v+s) div 2 rfind p = (if f (p, q) == z then (p, q) : find (u, v) (p−1, q+1) f z else find (u, v) (p, q+1) f z) ++ find (p+1, q−1) (r, s) f z cfind q = find (u, v) (p−1, q+1) f z ++ (if f (p, q) == z then (p, q) : find (p+1, q−1) (r, s) f z else find (p+1, q) (r, s) f z)
Dart // 長方形 (u, v)〜(r, s) を再帰的に探索。results に副作用で追加していく。 void find(Fn f, int z, int u, int v, int r, int s, List<(int, int)> results) { if (u > r || v < s) return; if (v - s <= r - u) { // 列の方が多い → 中央の行 q で二分探索 final q = (v + s) ~/ 2; final p = bsearch((x) => f(x, q), u - 1, r + 1, z); if (f(p, q) == z) { results.add((p, q)); find(f, z, u, v, p - 1, q + 1, results); } else { find(f, z, u, v, p, q + 1, results); } find(f, z, p + 1, q - 1, r, s, results); } else { // 行の方が多い → 中央の列 p で二分探索 final p = (u + r) ~/ 2; final q = bsearch((y) => f(p, y), s - 1, v + 1, z); find(f, z, u, v, p - 1, q + 1, results); if (f(p, q) == z) { results.add((p, q)); find(f, z, p + 1, q - 1, r, s, results); } else { find(f, z, p + 1, q, r, s, results); } } }

全体の invert はこう:

invert f z = find (0, m) (n, 0) f z where m = bsearch (λy → f (0, y)) (−1, z+1) z n = bsearch (λx → f (x, 0)) (−1, z+1) z
Dart List<(int, int)> invertMary(Fn f, int z) { final m = mOf(f, z); final n = nOf(f, z); final results = <(int, int)>[]; find(f, z, 0, m, n, 0, results); return results; }

計算量の証明

最悪ケースの漸化式は:

T(m, n) = log n + 2 T(m/2, n/2)

これを解くと T(m, n) = O(m log(n/m))。アンの下界と一致し、漸近的に最適です!

実測結果

5 つの関数について invert fi 5000 を計算した結果:

評価回数(図 3.2)

アルゴリズムf0f1f2f3f4
Anne(鞍点)75015011666850689989
Theo(+二分探索)25373817491575025
Mary(分割統治)12142445181134

絶対実行時間(秒、図 3.3)

アルゴリズムf0f1f2f3f4
Anne(鞍点)0.420.400.170.150.54
Theo(+二分探索)0.060.010.050.010.15
Mary(分割統治)0.010.010.020.020.01
結論 メアリーの分割統治法が圧倒的に速い。特に f0 では、鞍点探索の 7501 回に対してわずか 121 回と、60 倍以上の高速化。
Dart // 動作確認: f0(x, y) = 2^y * (2x + 1) - 1 void main() { int f0(int x, int y) => (1 << y) * (2 * x + 1) - 1; print(invertSaddleback(f0, 15)); // (0, 4), (1, 3), (3, 2), (7, 1), (15, 0) print(invertMary(f0, 15)); // 同じ 5 要素(順序は異なる場合あり) }

まとめ

著者の裏話 この真珠は元々、著者が Oxford の入試面接で使っていた問題。多くの受験生が「二分探索を使いたい」と言うのを、著者は「鞍点探索が金字塔だ」と思って軌道修正させていた。ところが、あとになって受験生たちの方が正しかったと気づいたそうです。常識を疑うことの大切さを示すエピソードです。