焼きなまし法（simulated annealing）の理論について

概要

この記事では、焼きなまし法を支える理論について数学的な側面から解説することを目指します。
焼きなまし法をマルコフ過程としてモデリングすることにより、アルゴリズムが大域的最適解に収束するための条件について調べよう、と言うのが全体的な流れです。
連続と離散のギャップにより、ここで紹介する最適解への収束条件は実装されたアルゴリズムの場合には成立しないことを最初に注意しておきます。

焼きなまし法の背景には統計物理の知識も関係してきますが、本稿では触れないこととします。興味のある方は参考文献に目を通されてください。

こちらの記事で全体的な解説をしているので参考にどうぞ。

導入

Simulated annealingとは、物性物理学における「焼きなまし法」からのアナロジーとして生み出された近似解法を指します。 焼きなまし法では、高温の物体を平衡状態までゆっくり冷却することによって、元々の状態よりもさらに内部エネルギーが低い状態を得る可能性があります。これを最適化に応用することで生まれたアルゴリズムがSimulated annealingです。

Simulated annealingは、各状態からその近傍への確率的な遷移を繰り返すアルゴリズムと見なすことができ、それはマルコフ連鎖によってうまくモデリングすることができます。

まずはマルコフ連鎖についての簡単な説明を行い、simulated annealingの理論、実装の説明に移っていこうと思います。

記号

• ​$\mathcal{R}$​：状態空間
• ​$\mathcal{R}_{\rm opt}$：最適値を与える状態全体の集合

マルコフ連鎖（Markov chain）

​$\mathcal{S}$を可算集合として、​​$\mathcal{S}$上に値を取る確率過程​$(X_n,P)\quad(n=0,1,\dots)$が次を満たすとき、$(X_n,P)$​をマルコフ連鎖といいます。

簡単に表現すれば、マルコフ過程は「次の行動が現在の状態のみによって決まるような確率過程」のことです。

ではこのマルコフ過程を用いて、焼きなまし法をモデル化していきましょう。

​$P_{ij}(k-1,k)$を$k-1$​回目の試行を行なった後の状態$i$​から$k$​回目の試行によって状態$j$​に遷移する確率とします。 このとき、​$a_i(k)$を$k$​回目の試行後に状態​$i$である確率とすると、

$$
a_i(k)=\sum_{l}a_l(k-1)P_{li}(k-1,k)
$$

のように表現することができます。

さらに、​$X(k)$を$k$​回目の試行が終わった後の状態を表す確率変数とすると、

のように書くことができます。

条件付き確率​$P_{ij}(k-1,k)$が試行回数​に依存しない場合をhomogeneousなマルコフ連鎖、依存する場合をinhomogeneousなマルコフ連鎖といいます。

Homogeneousマルコフ連鎖が既約とは、「$\forall i,j,k$​に対し、ある有限の整数​$n$が存在して、$(P)^n>0$​を満たす」ことを指します。

またHomogeneousマルコフ連鎖が非周期的とは「​$\forall i$について、$(P_{ii})^n>0$​を満たす全ての$n$​の最大公約数が1になる」ことをいいます。

ここで紹介した概念を元に、最適解へ収束するための条件が導かれていきます。

Simulated annealingの場合

上記をsimulated annealingに適用すると以下のようになります。 状態空間を$\mathcal{R}$​、制御パラメーターを$c$​とします。

この時​$P_{ij}(k-1, k)$は、​$i,j\in\mathcal{R}$に対して、状態$i$​から状態$j$​に遷移する確率を表し、​$X(k)$は​$k$回の遷移の後に観測される状態を表します。 このことから、​$P_{ij}(k-1, k)$を遷移確率、​各成分が遷移確率であるような行列$P(k-1, k)$​を遷移行列と呼びます。

焼きなまし法とのアナロジーから、遷移確率は制御パラメーター$c$​に依存し、​$c$が定数の時、遷移行列$P$​に対応するマルコフ連鎖はhomogeneousになります。この場合、遷移確率は以下のように書き下すことができます。

​$G_{ij}(c)$を生成確率（状態​から状態​が発生する確率）、​$A_{ij}(c)$を受理確率（状態​から状態​に実際に遷移する確率）と呼び、これらを成分に持つ行列​$G(c),A(c)$をそれぞれ生成行列、受理行列と呼びます。​​$G_{ij}(c), ​A_{ij}(c)$はともに条件付き確率です。

​$P_{ij}(c)$の定義から、遷移行列​は確率行列になっています。

simulated annealingでは、アルゴリズム中で制御パラメータの値を小さくしていきますが、その方法によって二つに分類することができます。

1. Homogeneous アルゴリズムは複数のhomogeneousなマルコフ連鎖で構成され、それぞれのマルコフ連鎖はある固定された制御パラメータによって生成されます。制御パラメーターは、次のマルコフ連鎖に移行するタイミングで更新します。
2. inhomogeneous アルゴリズムが一つのinhomogeneousなマルコフ連鎖で表現されます。制御パラメーターは各遷移ごとに更新します。

疑似コードで書くと以下のようになります。

# homogeneous algorithm
while T > TFINAL:
    for _ in range(MAXSTEP):# 同一のマルコフ連鎖内での遷移
        j # 現在の状態をiとして、次の状態jを選択する
        swap(i, j, T)# 状態iから状態jに遷移
    T *= C
                
# inhomogeneous algorithm
while T > TFINAL:
    j # 現在の状態をiとして、次の状態jを選択する
    swap(i, j, T)# 状態iから状態jに遷移
    T *= C

最適性を保証するための条件

注意：全ての実装は離散的になるため、ここで述べる最適性の議論は実装上あまり役に立つものではない。（実用上最適性は保証されない）

Homogeneous algorithm

homogeneousなアルゴリズムにおける遷移確率は以下のようにかけます。

ポイントは、

1. アルゴリズムに対応するマルコフ連鎖が定常分布を持ち、
2. その定常分布が最適解集合（もしくはその部分集合）上の一様分布$​​​\pi=(\pi_i)=\left\{\begin{array}{ll}|\mathcal{R}_{opt}|^{-1}&(i\in\mathcal{R}_{opt})\\0&(\text{o.w})\end{array}\right.$​になる

ことです。この節ではそのための条件を求めていきます。

まず、マルコフ連鎖の定常分布​$q$は、各成分が以下の式で定義されるようなベクトルで定義されます。
$$
q_i=\lim_{k\to\infty}P{X(k)=i|X(0)=j}(\forall j\in\mathcal{R})
$$

この定義から、homogeneousなマルコフ連鎖の遷移行列$P$を​、初期状態における確率分布を$a(0)$​として

のようにかけます。このことから、定常分布はsimulated annealingにおいて無限回の遷移を行なったあとの各状態が従う確率分布だと言うことがわかります。

マルコフ連鎖が定常分布を持つための（十分）条件

定理1

後ほど$\forall i,j,c$​に対して​$A_{ij}(c)>0$を仮定するため、アルゴリズムに対応するマルコフ連鎖の既約性を保証するためには​によって定義されるマルコフ連鎖の既約性を保証すれば十分です。

また、非周期性を言うために、以下の事実を使います。

以上を踏まえると、アルゴリズムに対応するマルコフ連鎖が定常分布を持つためには以下の二つを仮定すれば十分。

1. ​$\forall i\in\mathcal{R},\exists p\geq1, \exists l_0,\dots, l_p(l_0=iかつl_p=j)\space\text{s.t.}\space G_{l_k l_{k+1}}(c)>0(k=0,1,\dots,p-1)$
2. $\forall c>0,\exists i_c,j_c\in\mathcal{R}\space\text{s.t}\space A_{i_cj_c}(c)<1$​

本節の締めくくりとして、二つ目の条件が成り立つための条件を考えていきます。

定常分布が最適解（の部分）集合上の一様分布になっているための（十分）条件

まずは、以下の事実に基づく最小限の要求から紹介する。

事実

この事実を用いると、$\lim_{c\to+0}q(c)=\pi$​が成り立つための十分条件は、

となります。以上から、上記条件が成り立つときに定常分布が最適解（の部分）集合上の一様分布に収束することがわかりましたが、この条件から​を陽に得ることは大変難しいです。そこで、​$A(c), G(c)$に更なる条件を与えることで$\phi$を明示的に与えます。

定理2（Folklore）

​証明：略

注意：​$A_{ij}(c)$がコスト関数の値のみに依存していることを仮定している。

Note：

• 定常分布が$\pi$​に収束するためには以下の条件が十分
  5. $\forall i,j\in\mathcal{R},C(i)>C(j)\Rightarrow \lim_{c\to+0}A_{ij}(c)=0$
• 条件3と条件5によって定常分布の収束が保証されている。
• 条件1は以下の条件1’で置き換え可能
  1′. $\forall i\in\mathcal{R},G_{ij}=\left\{\begin{array}{ll}|\mathcal{R}_i|^{-1}&j\in\mathcal{R}_i\\0&\rm elsewhere\end{array}\right.$ ​
  ただし、​$\mathcal{R}_i=\left\{j\in\mathcal{R} | G_{ij}\neq0\right\}$
  この条件に置き換えた場合、$q_i(c)=\frac{|\mathcal{R}_i|A_{i_0i}(c)}{\sum_{j\in\mathcal{R}}|\mathcal{R}_j|A_{i_0j}(c)}$​となる。
• 条件1と条件1’がともに成り立つとき、$|\mathcal{R}_i|$​は​$i\in\mathcal{R}$とは独立な値をとる。

以上の議論から、homogeneousなアルゴリズムが大域的最適解に収束するための十分条件を求めることができました。

2. Inhomogeneous

遷移行列は以下のようにかけます。

この節では、

1. $\lim_{k\to\infty}c_k=0$
2. $c_k\geq c_{k+1},\space k=0,1,2,\dots$

を仮定します。
まずは定常分布が最適解の部分集合上に収束するための十分条件を与え、その後必要十分条件を考えていきます。

十分条件

定義1（弱エルゴード）

定義2（強エルゴード）​

※​$X(k)$が従う分布が収束するということ
注：homogeneousマルコフ連鎖には強弱の区別はない。

定理3

​定理4

受理確率、生成確率に、homogeneousなアルゴリズムに対する条件と同じものを仮定すると、各$k\geq0$​に対して$P(k,k+1)$​を遷移確率、​$q(c_k)$を定常分布に持つマルコフ連鎖が存在します。（​$q(c_k)$は​$P$の固有ベクトル）

さらにいくつか仮定すると、​$c_k\to0$の元で​​$q(c_k)$が$\pi\left(\pi_i=\left\{\begin{array}{ll}|\mathcal{R}_{opt}|^{-1}&i\in\mathcal{R}_{opt}\\0&\text{O.W.}\end{array}\right.\right)$​に収束するようにできます。​

$\pi(k)=q(c_k)$として定理4を適用すると、​$q(c_k)$を定常分布にもつマルコフ連鎖が強エルゴード性を持つことは、

• マルコフ連鎖が弱エルゴードである
• $q(c_k),k=0,1,2,\dots$が​(1)を満たす

ことだと言い換えることができます。 ​
$\pi_i$の定義と強エルゴード性より、
$$
\begin{align}
&\quad\lim_{k\to\infty}P\left(X(k)\in\mathcal{R}_{opt}\right)\\ &=\lim _{k\to\infty}\sum_{j\in\mathcal{R}_{opt}} P\left\{X(k)=j\right\}\\
&=\sum_{j\in\mathcal{R}_{opt}}\lim _{k\to\infty}P\left\{X(k)=j\right\}(有限和だからok)\\
&=1
\end{align}
$$
がわかり、大域的最適解に収束することがわかります。