ざっくり解説: 焼きなまし法はなぜうまくいくのか？

概要

焼きなまし法は遺伝的アルゴリズムなどと並ぶ最適化問題の発見的解法の一つとして有名であり, 競技プログラミングの文脈ではマラソン(厳密解を求めるのが困難とされる一つの最適化問題に長時間取り組み最も最適値に近い解を得るという種目) においては常套手段の一つとして用いられています. 従ってこの文脈では様々な解説記事があるのですが, いわゆるマラソン競技特有の「コツ」や実装方法の説明がほとんどで, 焼きなまし法自体を数理的に深掘りした記事もなくはないのですが, やはり「なぜうまくいくのか」の観点で説明している記事はほとんど見たことがないので, なんとなくこの辺りが気になるという人向けにこの記事で解説していきます (公開時は把握してなかったのですが, 本記事の内容はこちらの記事とほぼ同じ内容になっていますので, 数理的な説明をする記事が全くないというわけではありません). 内容としてはメトロポリス・ヘイスティング法の説明になりますが, 数学的な厳密性は犠牲にして自分の理解で執筆していますので, 厳密な議論まで気になる方は論文を読んでください (もしかしたら他の説明の仕方があるかもしれません).

結論から述べると

• 焼きなまし法は, 遷移確率が時間とともに変動するランダムウォークを実用的な観点から修正したものとみなすことができる.
• このランダムウォークは, 遷移確率の変動が各時刻において十分に微小であれば, そのランダムウォークの収束先は最適解上の分布になることが保証される.
• しかし, 理論的に保証される収束は非常に遅いため, 実用の場面では, 適当なところでランダムウォークを打ち切ったり, より良い解を得るための(理論保証のない)修正を施す. これがマラソンなどの場面でよく見かける焼きなまし法である.

この記事にアクセスする人は焼きなまし法についてはなんとなく知ってるということを仮定します. ただ, 念の為, 初見の人にもアクセスできるよう, 焼きなまし法に関する説明を次の章で与えておきます. 焼きなまし法を知っている人はその次の2章までスキップして構いません.

1. 焼きなまし法とは?

焼きなまし法とは最適化問題に対し, 最適に近いであろう実行可能解を得る手法の一つで, 局所探索と呼ばれるカテゴリーに属します. 局所探索とは, 何かしらの実行可能解を一つ用意し, その解を逐次更新していくというアルゴリズムの総称であり, 焼きなまし法以外には山登り法などが知られています. ここでは以下の記法を用います.

• 実行可能解全体の空間は有限集合とし記号$S$で表します. 
• 実行可能解を$x,y,z\in S$などと表します
• 目的関数値を$f\colon S \to \mathbb{R}_{\ge 0}$で表し, 最小化を目指すこととします.
• 実行可能解$x\in S$の近傍集合を$N(x) \subseteq S$とし, その要素$y\in N(x)$を近傍と呼ぶ.

すなわち, $S$という有限集合と$S$上に非負値をとる関数$f$があり, $f(x)$が最小となる$x\in S$を見つけることが目的となります. 直感的には近傍$y \in N(x)$は$x$と「ちょっとだけ違う実行可能解」です. 実行可能解と目的関数の定義は考える最適化問題に依存しますが, 実行可能解はアルゴリズムの設計者が定義するため, 同じ最適化問題であっても様々な近傍を考えることができます.

山登り法では, 現在の実行可能解$x\in S$に対し, $f(y)$が最小となる近傍$y\in N(x)$を選び, それを$x$に置き換えて同じ処理を繰り返すというアルゴリズムです (または$N(x)$の元$y$を順番に一つずつ見ていき, $f(x)>f(y)$を満たす近傍$y$が見つかったらその瞬間に$x:=y$とするようなものも考えられる. これは, $|N(x)|$の要素数が非常に大きい時は$\min_{y\in N(x)} f(y)$を計算に時間がかかる際に有効な方法である).

このアルゴリズムだと谷底 (全ての近傍の目的関数値が今の解より真に良くならないような解; 局所最適解と呼ぶ) にいきつくとそれ以上動かなくなるのでアルゴリズムは終了します. 上図の例では偶然にも最適解に辿り着きますが, 一般に谷底は複数あり, どの谷底に至るかは一番最初に選んだ解に依存するため, 最適性の保証はありません (ただし空間が連続的で目的関数値が凸性を持つなど, 特定の状況下では最適性が保証できます).

一方, 焼きなまし法では, 現在の実行可能解$x\in S$に対し, 近傍$y\in N(x)$を順番に見ていき,

1. 目的関数値が改善, すなわち$f(x)>f(y)$を満たすとき, $x:=y$として次の反復へ移行
2. そうでない場合, $f(x),f(y)\in\mathbb{R}_{\ge 0}$, そしてこれまでの反復回数$t$のみから定まるある確率$p$で$x:=y$として次の反復へ移行

というアルゴリズムです. 例えば$p=0$ならば山登り法に一致します.

焼きなまし法では確率的に登ることがある.

ここで確率$p$は

$$ p = \exp\left(-\frac{f(y) - f(x)}{T}\right) \tag{1} $$

とします. ここで$T=T(t)>0$は温度と呼ばれるパラメータであり, これまでの反復回数$t$の関数であり, $t\to\infty$としたとき$T\to 0$となります. ここで, $f(y)\ge f(x)$を仮定しているため常に$p\in[0,1]$です.

確率$p$については以下の観察ができます:

• 解の改悪幅$f(y)-f(x)$が大きいほど$p$が小さい, すなわち$x\to y$の遷移が発生しにくい.
• 改悪幅を正の値で固定したとき, 温度$T$が$0$に近いほど$p$は小さくなります. つまり, $T$が大きいときは改悪($f(x)<f(y)$なる$y\in N(x)$への更新)が発生しやすくなるが, 反復を繰り返し$t\to\infty$としていくと$T\to 0$に近づくため, より山登り法に近い振る舞いをするようになる.

温度$T(t)$の定め方は「非常にゆっくり$T\to 0$に収束する関数」であることが理想ですが, 実際に焼きなまし法を走らせる際は手っ取り早く良い解を得たい場合がほとんどなので, あらかじめ反復回数を決定しておき, 最終的に$T\approx 0$になるように線形で下がっていくようなものを考えることもあります. ただし理論的な保証があるのは「非常にゆっくり」収束するという条件が必要なので, この時点で理論保証から外れてしまうことになります.

2. メトロポリス-ヘイスティング法

突然ですが, 有限集合$S$上の関数

$$ f\colon S \to \mathbb{R}_{\ge 0} $$

が定まっているとき, $f(x)$を最小にする$x\in S$を求めるにはどうすれば良いでしょうか? 唐突ですが, 温度と呼ばれるパラメータ$T>0$に対して$S$上の確率分布$\pi\in[0,1]^S$を以下で定めます:

$$ \pi (x) = \frac{\exp(- f(x)/T)}{Z}. \tag{2}$$

ここで$Z = \sum_{x\in S}\exp(-f(x)/T)$は$\pi$を分布に正規化する値 (様々な文脈で分配関数と呼ばれる値) です. また, この分布$\pi$を(エネルギー$f$に対する)ボルツマン分布と呼ぶらしいです.

分布$\pi$に従って生成される確率変数$X \sim \pi$を考えると, 以下の観察が芽生えます:

• 温度$T$が非常に大きい値, 極端に言えば$T=\infty$ のとき, $f(x)$の影響は無視されるため, $\pi$は$S$上で一様ランダムとなります.
• 温度$T$が非常に小さい値, すなわち$T\approx 0$ のとき, $\pi$は最適解$S^*=\{x\in S\colon f(x) = \min_{y\in S}f(y)\}$上の一様分布となります.

従って, 理想的には$T\approx 0$において$\pi$に従う確率変数を生成できれば最適化問題を解いたことになります!! ところが, 分配関数$Z$が簡単に計算できそうな見た目をしてないので, $\pi$に従う確率変数をサンプリングするということ自体が非自明なタスクになっています. メトロポリス-ヘイスティング法とは, 目標となる分布$\mu$があったとき, ランダムウォークを用いてこの分布に非常に近い分布に従う確率変数を生成する手法の一つです.

具体的には, $S$上の何らかのランダムウォークの遷移確率行列$P \in [0,1]^{S\times S}$が与えられているとします. ここで$P(x,y)$は$x$から$y$に遷移する確率を表しており, 例えば焼きなまし法の例でいうと, 一様ランダムな近傍に遷移する, というランダムウォークだと思えばOKです. 一般にランダムウォークはその遷移がとある条件 (例えば遷移の強連結性など) を満たせば, 定常分布と呼ばれるある分布に収束するという事実が知られています. 今考えているランダムウォーク$P$もその条件を満たしているとします. ただし$P$の定常分布は目標としている分布$\pi$とは異なっていることもあります.

メトロポリス-ヘイスティングス法のイメージ

メトロポリス-ヘイスティング法とは, 手元にあるランダムウォーク$P$に適当な重さの自己ループの追加などにより修正し, 目標分布$\pi$を定常分布とするような別のランダムウォーク$P'$をうまく定めるという手法です.

より詳細を述べると, 修正後のランダムウォーク$P'$は以下のような形をしています: ランダムウォーク$P$により, $x\to y$という遷移が発生したとき, 確率$a(x,y)$でその遷移を受理し実際に遷移し, 確率$1-a(x,y)$でその遷移を拒否し$x$にとどまる.

簡単のため$P(x,x)=0$とすると, $P'\in[0,1]^{S\times S}$の成分は

$$
P'(x,y) = \begin{cases}
P(x,y)\cdot a(x,y) & & \text{if }x\ne y,\\
1-\sum_{y} P(x,y)\cdot a(x,y) & & \text{otherwise}
\end{cases}
$$

となります. 肝心の受理確率$a(x,y)$の定め方ですが, $P'$が定常分布$\pi$を持つように, そして$a(x,y)$の値ができるだけ大きくなるように (自己ループの確率が少ない方が遷移する回数が多くなるため) 設計します. 詳細は省きますが釣り合い条件と呼ばれる条件から, $P'$は次の等式

$$ {}^{\forall} x,y\in S,\quad\pi(x) P'(x,y) = \pi(y) P'(y,x) $$

を満たしていなければなりません. これに$P'$の成分を代入すると

$$ \pi(x) a(x,y) P(x,y) = \pi(y) a(y,x) P(y,x) $$

となります. 

計算の簡単のため, $P$が対称行列だとしましょう. すると上の等式は$\pi(x)a(x,y) = \pi(y)a(y,x)$となります. ここで, できるだけ$a(\cdot,\cdot)$の値を大きくとると, 片方は$1$に設定できます. 従って

\begin{align*}
&a(x,y) = \min\left\{ \frac{\pi(y)}{\pi(x)},\,1\right\},\\
&a(y,x) = \min\left\{ \frac{\pi(x)}{\pi(y)},\,1\right\}
\end{align*}

を得ます. このようにして定まる受理確率$a(x,y)$は$\pi(x)$と$\pi(y)$の比のみに依存するという素晴らしい性質を持っています. 実際, $\pi$が(2)のような形になっているとき, その比は分配関数$Z$に依存せず,

$$ \frac{\pi(y)}{\pi(x)} = \exp\left( - \frac{f(y)-f(x)}{T} \right) $$

で簡単に計算できます. この値はまさに(1)で考えていた値と一致します. つまり, 焼きなまし方とは, メトロポリス-ヘイスティングス法の受理確率を模した確率的な遷移を許容する局所探索といえます.

なお, 実際の焼きなまし法では温度パラメータ$T$は遷移回数とともに変動します. 直感的にはこの変動は十分ゆっくりであれば局所的には$T$が固定された値とみなせるので解析できそうですが, この直感を厳密な形にするにはそれなりの苦労が必要のようです (例えばこの記事. この辺は私もまだわかっていない)

なお, メトロポリス-ヘイスティングス法で得られるランダムウォーク$P'$は, 収束先の定常分布が$\pi$であるということのみが保証されるため, その収束の速さについては何も保証がありません. しかし, 適当な仮定の下では何かしらの収束時間のバウンドを与えることは可能となります.

3. まとめ

冒頭にも述べましたが, 要するに

1. 最適化問題を解くというタスクをある目標分布に従うサンプリングに置き換える
2. このサンプリングタスクをメトロポリス-ヘイスティング法で解く
3. プロセス1,2の構造をできるだけ保ちつつ, 実用的に手っ取り早く良い解を得るための工夫を施して実装したものが焼きなまし法である

ということになります (もしもうちょっと深い理論の話があったら教えていただけると幸いです). 

なお, 最適化問題は$f$が凸性を持つときは「局所最適解は最適解の一つである」が成り立つため非常に解きやすくなるのですが, 同様に分布$\pi$は対数凹性を持つときはサンプリングしやすいという事実があるようです. 実際, $f$が凸のとき, 対応するボルツマン分布$\pi$は$\log \pi(x) = -\frac{f(x)}{T}$より, 対数凹性を持ちます.