整数変数に対するSimulated Annealingについて①

この記事はJij Inc. Advent Calendar 2024の20日目の記事です．

概要

この記事では整数変数で構成される制約なし最適化問題に対する擬似焼きなまし法(Simulated Annealing: SA)について解説します。SAはバイナリ変数で構成される二次制約無し二値最適化問題(QUBO)に対して説明されることが多い^[1]ですが、アルゴリズム自体はとても単純なため、バイナリ変数を整数変数に拡張してもやることは明確です。ですが、ある程度SAを高速に動作させるためには変数を更新した際のエネルギー差分を保持している必要があり、整数変数の場合はこれが少し複雑になります。本記事ではこの部分を重点的に解説します。とは言っても書いていたら長くなったので、今回はアルゴリズムの簡単な紹介にとどめて、次回の記事で詳細に実装方法を説明することにします。

二次制約無し多値最適化問題

対象とする最適化問題としては以下のような、目的関数が整数変数の１次および２次式で構成されたものを考えます。

$$E(z)=∑i,jJi,jzizj+∑ihizi,Li≤zi≤Ui E(\bm{z})=\sum_{i, j}J_{i,j}z_{i}z_{j} + \sum_{i} h_i z_{i},\quad L_i \leq z_{i} \leq U_i$$

ここで、$z={z1,z2,…,zN}\bm{z}=\{z_1,z_2,\ldots,z_N\}$は整数変数、$NN$は変数の個数、$LiL_i$と$UiU_i$は$ii$番目の整数変数$ziz_i$の取りうる値の最小値と最大値を表します。$ziz_i$は整数変数なので、$zi∈{Li,Li−1,…,Ui}z_i\in\{L_i,L_i - 1,\ldots,U_i\}$ということになります。
以下では、$z\bm{z}$を状態と呼び^[2]、$E(z)E(\bm{z})$を状態$z\bm{z}$に対応するエネルギーと呼ぶことにします。特に$E(z)E(\bm{z})$が最小となる状態のことを基底状態と呼びます。

擬似焼きなまし法

擬似焼きなまし法(Simulated Annealing: SA)とは、目的関数^[3]などと呼ばれる関数の最小値とそれを与える変数の組を求めることを目的とした確率的なアルゴリズムです。もちろんこの記事で対象とする目的関数は上で説明した$E(z)E(\bm{z})$です^[4]。

アルゴリズムの概略

はじめにアルゴリズムの概略について説明しましょう。まず、$z\bm{z}$の初期値$z(0)\bm{z}^{(0)}$を用意します。これは通常ランダムに変数の値を決めれば良いです。そして、温度と呼ばれる非負実数$TT$を用意してこれの初期値を十分大きな値に設定しておきます^[5]。
これで初期条件が設定されました。以降は状態$z(0)\bm{z}^{(0)}$をある規則に従って更新していきます。$kk$ステップ目の状態を$z(k)\bm{z}^{(k)}$とし、温度を$TT$としましょう^[6]。
$z(k)\bm{z}^{(k)}$と何らかの意味で「近い」状態を生成し、これを$z′(k)\bm{z'}^{(k)}$とします。この状態は$k+1k+1$ステップの状態$z(k+1)\bm{z}^{(k+1)}$の候補状態です。
そして$z(k)\bm{z}^{(k)}$をある確率$pp$で$z′(k)\bm{z'}^{(k)}$に更新します。この確率は各状態に対応するエネルギーと温度$TT$のみに依存して定まります^[7]。つまり$p=p(E(zk),E(z′(k)),T)p=p(E(\bm{z}_k),E(\bm{z'}^{(k)}), T)$ということです。
この手続きを繰り返して状態を更新していきます。また、温度$TT$は状態を更新して行く中でゆっくりと0に近づけていきます。どのような頻度でどうのように0に近づけていくかは後ほど説明します。こうして得られた最終的な状態は$E(z)E(\bm{z})$が最小となる状態になっていることが期待されます。

候補状態の選び方

仮に現在の状態が$z(k)=(z1(k),z2(k),…,zN(k))\bm{z}^{(k)}=(z^{(k)}_1,z^{(k)}_2,\ldots,z^{(k)}_N)$だとしましょう。この状態の遷移先$z′(k)\bm{z'}^{(k)}$はどのようして決めればいいのでしょうか。$z′(k)\bm{z'}^{(k)}$と$z(k)\bm{z}^{(k)}$は何らかの意味で「近い」必要があります^[8]。様々な方法が考えられますが、ここでは最も単純な、$z(k)\bm{z}^{(k)}$のうちある一つの変数の値のみを変更することにしましょう。変数の値が一つだけ違うという意味で「近い」状態です。
この更新すべき変数は$z(k)=(z1(k),z2(k),…,zN(k))\bm{z}^{(k)}=(z^{(k)}_1,z^{(k)}_2,\ldots,z^{(k)}_N)$の$NN$個の変数の中からランダムに選べば良さそうです^[9]。
その変更する変数を$zm(k)z^{(k)}_m$としましょう。ですが、今の場合$zm(k)∈{Lm,Lm+1,…,Um}z^{(k)}_m\in\{L_m,L_m+1,\ldots,U_m\}$なので$Um−Lm>1U_m - L_m > 1$である限り、$zm(k)z^{(k)}_m$の値の変更先が複数あることになります。したがって、$zm(k)z^{(k)}_m$をいくつにするかまで決めないといけません。これについては次で説明します。

遷移確率の計算法

状態を$z(k)\bm{z}^{(k)}$から$z′(k)\bm{z'}^{(k)}$に更新する遷移確率$p=p(E(z(k)),E(z′(k)),T)p=p(E(\bm{z}^{(k)}), E(\bm{z'}^{(k)}), T)$をどのように計算するかについてはある程度の任意性があります。ここでは代表的な以下の３つの方法を考えます^[10]。

メトロポリス法

最もメジャーと言っても過言ではないほどよく用いられる更新方法です。
まず、候補状態$z′(k)\bm{z'}^{(k)}$を決めるにあたって、$zm(k)z^{(k)}_m$をどの値に更新するかですが、これは現在の値を除いた取りうる値の中からランダムに決めます。例えば$Lm=−1L_m=-1$、$Um=2U_m=2$であった場合$zm(k)∈{−1,0,1,2}z^{(k)}_m\in\{-1,0,1,2\}$ですが、現在の状態が$zm(k)=0z^{(k)}_m=0$であった場合は${−1,1,2}\{-1,1,2\}$の中からランダムに決めるということです。
これで$z′(k)\bm{z'}^{(k)}$が決まりました。実際にこの状態に更新するかは以下の確率で決定します。

$$p={exp⁡(−(E(z′k)−E(zk))/T)if E(z′k)>E(zk),1if E(z′k)≤E(zk). p= \begin{cases} \exp(-(E(\bm{z'}^{k}) - E(\bm{z}^{k}))/T) & \text{if } E(\bm{z'}^{k}) > E(\bm{z}^{k}), \\ 1 & \text{if } E(\bm{z'}^{k}) \leq E(\bm{z}^{k}). \end{cases}$$

熱浴法

こちらもとても有名な状態更新法の一つです。候補状態$z′(k)\bm{z'}^{(k)}$を決めるにあたって、$zm(k)z^{(k)}_m$をどの値に更新するかですが、熱浴法では$zm(k)∈{Lm,Lm+1,…,Um}z^{(k)}_m\in\{L_m,L_m+1,\ldots,U_m\}$の取りうる値全てを考慮してどの値に更新するかを確率的に決めます。
これら$Um−Lm+1U_m - L_m + 1$個ある取りうる値に対して、その値が小さい順に$μ=1,2,…,Um−Lm+1\mu=1,2,\ldots,U_m - L_m + 1$とラベル付けしましょう。そして、$z(k)\bm{z}^{(k)}$において、$mm$番目の変数$zm(k)z^{(k)}_m$を$μ\mu$でラベル付されたある値に変更した状態を$zμ(k)\bm{z}^{(k)}_{\mu}$と書きます。これが先程説明した候補状態$z′(k)\bm{z'}^{(k)}$に対応します。なお、$μ\mu$に値によっては現在の値$z(k)\bm{z}^{(k)}$と等しくなることに注意してください。
この前提のもとで、状態$z(k)\bm{z}^{(k)}$が$zμ(k)\bm{z}^{(k)}_{\mu}$に更新される確率は以下のように計算します。

$$pμ=exp⁡(−E(zμ(k))/T)∑μ=1Um−Lm+1exp⁡(−E(zμ(k))/T) p_{\mu}= \frac{\exp(-E(\bm{z}^{(k)}_{\mu})/T)}{\sum^{U_m - L_m + 1}_{\mu=1}\exp(-E(\bm{z}^{(k)}_{\mu})/T)}$$

諏訪藤堂法

2010年に提案された手法で、状態が更新されない確率である棄却率が平均的に最小となるような更新方法となっています^[11]。詳細は論文に書いてるのでそちらを読んでいただいた方がいいのですが、ここではこれまでの記法でどのように確率$pp$が表現されるか見ておきます。
熱浴法の説明の際と同様に、更新する変数$zm(k)∈{Lm,Lm+1,…,Um}z^{(k)}_m\in\{L_m,L_m+1,\ldots,U_m\}$の取りうる値に対して、その値が小さい順に$μ=1,2,…,Um−Lm+1\mu=1,2,\ldots,U_m - L_m + 1$とラベル付けします。そして、$z(k)\bm{z}^{(k)}$において、$mm$番目の変数$zm(k)z^{(k)}_m$を$μ\mu$でラベル付されたある値に変更した状態を$zμ(k)\bm{z}^{(k)}_{\mu}$と書きます。さらに、$zμ(k)\bm{z}^{(k)}_{\mu}$が現在の値$z(k)\bm{z}^{(k)}$と等しくなる$μ\mu$を$λ\lambda$と書きます。$zλ(k)=z(k)\bm{z}^{(k)}_{\lambda}=\bm{z}^{(k)}$ということです。
この前提のもとで、状態$z(k)\bm{z}^{(k)}$が$zμ(k)\bm{z}^{(k)}_{\mu}$に更新される確率は以下となります。

$$pμ=max⁡[0,min⁡(Δλ,μ,wλ+wμ−Δλ,μ,wλ,wμ)] p_{\mu} = \max\left[0, \min(\Delta_{\lambda,\mu}, w_{\lambda} + w_{\mu} - \Delta_{\lambda,\mu}, w_{\lambda}, w_{\mu})\right]$$

ここで、

$$Δλ,μ=Sλ−Sμ−1+w1(1≤λ,μ≤Um−Lm+1),Sμ=∑k=1μwμ,S0=SUm−Lm+1,wμ=exp⁡(−E(zμ(k))/T). \Delta_{\lambda,\mu}=S_{\lambda} - S_{\mu - 1} + w_1\quad (1 \leq \lambda,\mu \leq U_m - L_m + 1),\\ S_{\mu}=\sum^{\mu}_{k=1}w_{\mu},\quad S_0 = S_{U_m - L_m + 1},\\ w_{\mu}=\exp(-E(\bm{z}^{(k)}_{\mu})/T).$$

となります。ただし、1番目の状態$μ=1\mu=1$が最も重み$wμw_{\mu}$が大きくないといけません。また、この更新確率は$Um−Lm=1U_m - L_m = 1$、つまり各変数が２状態しか取れない場合はメトロポリス法と同じ更新確率になります。

温度$TT$の下げ方

最後に温度$TT$の下げ方を説明します。ここでは最もよく使われているであろう幾何冷却について説明します^[12]。まず、sweepという単位を導入します。これは変数の数$NN$回だけ変数の更新を試みることを1-sweepと定義します。そして$nn$-sweep目の温度を$TnT_n$と書きます。
幾何冷却とは以下の式に従って温度を下げることを指します。

$$Tn+1=aTn T_{n+1}=aT_{n}$$

毎回定数$aa$をかけるだけです。この$aa$はどのように決まるかというと、はじめに初期温度$TmaxT_{\text{max}}$と終温度$TminT_{\text{min}}$を決めておきます。そしてアルゴリズム全体でのsweep数を$num_sweeps\text{num\_sweeps}$とすればこれらの条件から

$$a=(TminTmax)1num_sweeps−1 a = \left(\frac{T_{\text{min}}}{T_{\text{max}}}\right)^{\frac{1}{\text{num\_sweeps} - 1}}$$

と決まります。これら初期温度、終温度、sweep数をどう決めるかですが、これについては次回、実装して実験する際に考えることにします。

アルゴリズムのまとめ

ここまでで一応SAというアルゴリズムの概略は説明しました。最後にまとめておきます。

1. 初期状態$z(0)\bm{z}^{(0)}$と初期状態$TminT_{\text{min}}$を決める。
2. 次を$num_sweeps\text{num\_sweeps}$回行う。
   2-1. 候補状態を生成し、その状態に確率$pp$で更新する。このプロセスを変数の数$NN$だけ繰り返す。
   2-2. 温度を下げる。
3. 最終的に得られた状態を最適化問題の答えとする。

まとめと次回予告

この記事では整数変数で構成された最適化問題に対する疑似焼きなまし法(Simulated annealing: SA)のざっくりとした解説を行いました。
一応ここまでの知識で実際にSAを行うことができるはずですが、このアルゴリズムの中で最も重たい処理が遷移確率$pp$の計算です。というのもここで説明した方法では$pp$を計算するために状態$z\bm{z}$に対するエネルギー$E(z)E(\bm{z})$を計算する必要があります。この計算をそのまま行うと時間がかかってしまうので工夫を考える必要があります。具体的には少し式変形すると遷移確率$pp$は状態を遷移した際のエネルギー差分$ΔE=E(z′)−E(z)\Delta E=E(\bm{z'}) - E(\bm{z})$で表現することができます。そこで、このエネルギー差分をあらかじめ計算しておけば遷移確率$pp$の計算を高速に行うことができます。ただし、このあらかじめ計算したエネルギー差分は、実際に状態が更新された際に改めて計算する必要があります。次回はここら辺の実装上の注意点を解説しようと思います。
追記: 続きはこちら。

募集してます！！！

Jijでは各ポジションを絶賛採用中です！
ご興味ございましたらカジュアル面談からでもぜひ！ご連絡お待ちしております。
数理最適化エンジニア（採用情報）
量子コンピューティングソフトウェアエンジニア（採用情報）
Rustエンジニア[アルゴリズムエンジニア]（採用情報）
研究チームPMO（採用情報）
オープンポジション等その他の採用情報はこちら

脚注

1. 弊社が開発しているOpenJijもSAはQUBOを対象としています。

2. 筆者は物理出身のためこのような言葉遣いをしてしまいます。数理最適化など他の分野でこれらの言葉遣いがどこまで一般的か筆者は知りません。

3. 他にはコスト関数、エネルギー関数などと呼ばれることがあります。

4. SAのアルゴリズムは非常に一般的なので目的関数の形に制限はありません。ある変数の組$x\bm{x}$とその変数に対する実数値$E(x)E(\bm{x})$が定まっていれば適用することができます。

5. 十分大きいというのは何に対してかというと２次項や１次項の係数$Ji,jJ_{i,j}$、$hi,jh_{i,j}$に対してということになります。

6. 後で説明するように温度$TT$もアルゴリズムのステップが進行するにつれて更新されていきます。しかしこの更新は変数$z\bm{z}$の更新よりも頻度が低く、更に単純な規則なのでここではアルゴリズムのステップ依存性を持たせず単に温度$TT$と表記しています。

7. 確率$pp$は$E(zk),E(zk+1),TE(\bm{z}_k),E(\bm{z}_{k+1}), T$のみから定まる０以上1以下の関数であれば何でも良い、、、わけではありません。もちろんアルゴリズムを実行することはできますが、へんな確率を選ぶと基底状態が得られる保証がなくなってしまいます。

8. 「近く」なくてもアルゴリズムは実行できますが、適当にやると基底状態への収束性が悪化したりします。

9. 筆者の経験上はランダムに選ぶのではなく、変数の番号順に逐次的に更新する変数を選ぶほうが収束性がいい場合が多いです。

10. 正直なところ添字が多くなって分かりづらくなっている気もするので、そう感じた方はタイトルでweb検索をしてみてください。色々と記事がヒットするはずです。

11. こちらで無料で読むことができます。

12. 他には線形冷却$Tn+1=a+TnT_{n+1}=a + T_{n}$や対数冷却$Tn+1=Tn/log⁡(n)T_{n+1}=T_{n}/\log(n)$などがあります。対数冷却は基底状態、つまり$E(z)E(\bm{z})$を最小にする$z\bm{z}$に漸近的に収束する保証がありますが、あまりに冷却速度が遅いため、実際にはほとんど用いられません。