第3章
鞍点探索を改良する
(やさしい版)
Pearls of Functional Algorithm Design(関数プログラミングによるアルゴリズム設計の真珠)
どんな問題?
関数 f(x, y) と数 z が与えられたとき、f(x, y) = z となるすべての組 (x, y) を見つけるのがこの章の問題です。
ただし前提として、f は自然数の組から自然数への関数で、各引数について「厳密に増える」ものとします(x や y を大きくすると値も必ず大きくなる)。それ以上のことは何も仮定できません。
ポイント
f の計算コストは非常に大きいかもしれないので、目標は「f を評価する回数をできるだけ少なくする」ことです。
この章では、先生と 4 人の生徒(アン、ジャック、メアリー、テオ)の対話形式で、探索アルゴリズムを段階的に改良していきます。
ステップ1:素朴な探索(ジャックの案)
まずは全部試す
f は増加関数なので、f(x, y) = z なら必ず x ≤ z かつ y ≤ z。だから (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 と比べます:
- f(u, v) < z → 列 u の下側はすべて z より小さいので捨てる(右へ進む)
- f(u, v) > z → 行 v の右側はすべて z より大きいので捨てる(下へ進む)
- 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, n が z よりずっと小さければ、O(log z) で済みます。
豆知識
この戦略は 鞍点探索(saddleback search)と呼ばれ、David Gries が名付けました。f の三次元プロットが左下最小・右上最大・両翼のような「鞍」の形に見えることから来ています。
ステップ3:分割統治法(メアリーの提案)
中央の要素で 2 つに分ける
「二分探索の二次元版」を考えたらどうか?というのがメアリーのアイデアです。長方形 (u, v)〜(r, s) について、中央の (p, q) を調べます:
- f(p, q) < z → 左下の長方形 A を捨てる
- f(p, q) > z → 右上の長方形 B を捨てる
- f(p, q) = z → 両方捨てる
問題点
残るのは 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)。m ≤ n なら水平、m ≥ n なら垂直カット。どちらも鞍点探索より速い。
ステップ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))
m ≤ n のときは Ω(m log(n/m)) が下界。ジャックの解はまだこの下界を達成していません。
最適な解(メアリーの最終案)
メアリーは 中央の行に二分探索を使う戦略を提案します。m ≤ n(列の方が多い)と仮定して、中央の行 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)
| アルゴリズム | f0 | f1 | f2 | f3 | f4 |
| Anne(鞍点) | 7501 | 5011 | 6668 | 5068 | 9989 |
| Theo(+二分探索) | 2537 | 38 | 1749 | 157 | 5025 |
| Mary(分割統治) | 121 | 42 | 445 | 181 | 134 |
絶対実行時間(秒、図 3.3)
| アルゴリズム | f0 | f1 | f2 | f3 | f4 |
| Anne(鞍点) | 0.42 | 0.40 | 0.17 | 0.15 | 0.54 |
| Theo(+二分探索) | 0.06 | 0.01 | 0.05 | 0.01 | 0.15 |
| Mary(分割統治) | 0.01 | 0.01 | 0.02 | 0.02 | 0.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 要素(順序は異なる場合あり)
}
まとめ
- 鞍点探索は 25 年以上「最適」と信じられてきたが、実は漸近的に最適ではない。
- 各行の中央で二分探索を行う 分割統治法が漸近的に最適な O(m log(n/m)) を達成する。
- 関数型プログラミングは「再帰+リスト連結」で分割統治法を自然に書ける。配列とループしかないと書きにくい。
著者の裏話
この真珠は元々、著者が Oxford の入試面接で使っていた問題。多くの受験生が「二分探索を使いたい」と言うのを、著者は「鞍点探索が金字塔だ」と思って軌道修正させていた。ところが、あとになって受験生たちの方が正しかったと気づいたそうです。常識を疑うことの大切さを示すエピソードです。