この記事はJij Inc. Advent Calendar 2024の20日目の記事です．
概要
この記事では整数変数で構成される制約なし最適化問題に対する擬似焼きなまし法(Simulated Annealing: SA)について解説します。SAはバイナリ変数で構成される二次制約無し二値最適化問題(QUBO)に対して説明されることが多いですが、アルゴリズム自体はとても単純なため、バイナリ変数を整数変数に拡張してもやることは明確です。ですが、ある程度SAを高速に動作させるためには変数を更新した際のエネルギー差分を保持している必要があり、整数変数の場合はこれが少し複雑になります。本記事ではこの部分を重点的に解説します。とは言っても書いていたら長くなったので、今回はアルゴリズムの簡単な紹介にとどめて、次回の記事で詳細に実装方法を説明することにします。
二次制約無し多値最適化問題
対象とする最適化問題としては以下のような、目的関数が整数変数の１次および２次式で構成されたものを考えます。
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
ここで、\bm{z}=\{z_1,z_2,\ldots,z_N\}は整数変数、Nは変数の個数、L_iとU_iはi番目の整数変数z_iの取りうる値の最小値と最大値を表します。z_iは整数変数なので、z_i\in\{L_i,L_i - 1,\ldots,U_i\}ということになります。
以下では、\bm{z}を状態と呼び、E(\bm{z})を状態\bm{z}に対応するエネルギーと呼ぶことにします。特にE(\bm{z})が最小となる状態のことを基底状態と呼びます。
擬似焼きなまし法
擬似焼きなまし法(Simulated Annealing: SA)とは、目的関数などと呼ばれる関数の最小値とそれを与える変数の組を求めることを目的とした確率的なアルゴリズムです。もちろんこの記事で対象とする目的関数は上で説明したE(\bm{z})です。
アルゴリズムの概略
はじめにアルゴリズムの概略について説明しましょう。まず、\bm{z}の初期値\bm{z}^{(0)}を用意します。これは通常ランダムに変数の値を決めれば良いです。そして、温度と呼ばれる非負実数Tを用意してこれの初期値を十分大きな値に設定しておきます。
これで初期条件が設定されました。以降は状態\bm{z}^{(0)}をある規則に従って更新していきます。kステップ目の状態を\bm{z}^{(k)}とし、温度をTとしましょう。
\bm{z}^{(k)}と何らかの意味で「近い」状態を生成し、これを\bm{z'}^{(k)}とします。この状態はk+1ステップの状態\bm{z}^{(k+1)}の候補状態です。
そして\bm{z}^{(k)}をある確率pで\bm{z'}^{(k)}に更新します。この確率は各状態に対応するエネルギーと温度Tのみに依存して定まります。つまりp=p(E(\bm{z}_k),E(\bm{z'}^{(k)}), T)ということです。
この手続きを繰り返して状態を更新していきます。また、温度Tは状態を更新して行く中でゆっくりと0に近づけていきます。どのような頻度でどうのように0に近づけていくかは後ほど説明します。こうして得られた最終的な状態はE(\bm{z})が最小となる状態になっていることが期待されます。
候補状態の選び方
仮に現在の状態が\bm{z}^{(k)}=(z^{(k)}_1,z^{(k)}_2,\ldots,z^{(k)}_N)だとしましょう。この状態の遷移先\bm{z'}^{(k)}はどのようして決めればいいのでしょうか。\bm{z'}^{(k)}と\bm{z}^{(k)}は何らかの意味で「近い」必要があります。様々な方法が考えられますが、ここでは最も単純な、\bm{z}^{(k)}のうちある一つの変数の値のみを変更することにしましょう。変数の値が一つだけ違うという意味で「近い」状態です。
この更新すべき変数は\bm{z}^{(k)}=(z^{(k)}_1,z^{(k)}_2,\ldots,z^{(k)}_N)のN個の変数の中からランダムに選べば良さそうです。
その変更する変数をz^{(k)}_mとしましょう。ですが、今の場合z^{(k)}_m\in\{L_m,L_m+1,\ldots,U_m\}なのでU_m - L_m > 1である限り、z^{(k)}_mの値の変更先が複数あることになります。したがって、z^{(k)}_mをいくつにするかまで決めないといけません。これについては次で説明します。
遷移確率の計算法
状態を\bm{z}^{(k)}から\bm{z'}^{(k)}に更新する遷移確率p=p(E(\bm{z}^{(k)}), E(\bm{z'}^{(k)}), T)をどのように計算するかについてはある程度の任意性があります。ここでは代表的な以下の３つの方法を考えます。
メトロポリス法
最もメジャーと言っても過言ではないほどよく用いられる更新方法です。
まず、候補状態\bm{z'}^{(k)}を決めるにあたって、z^{(k)}_mをどの値に更新するかですが、これは現在の値を除いた取りうる値の中からランダムに決めます。例えばL_m=-1、U_m=2であった場合z^{(k)}_m\in\{-1,0,1,2\}ですが、現在の状態がz^{(k)}_m=0であった場合は\{-1,1,2\}の中からランダムに決めるということです。
これで\bm{z'}^{(k)}が決まりました。実際にこの状態に更新するかは以下の確率で決定します。
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}
熱浴法
こちらもとても有名な状態更新法の一つです。候補状態\bm{z'}^{(k)}を決めるにあたって、z^{(k)}_mをどの値に更新するかですが、熱浴法ではz^{(k)}_m\in\{L_m,L_m+1,\ldots,U_m\}の取りうる値全てを考慮してどの値に更新するかを確率的に決めます。
これらU_m - L_m + 1個ある取りうる値に対して、その値が小さい順に\mu=1,2,\ldots,U_m - L_m + 1とラベル付けしましょう。そして、\bm{z}^{(k)}において、m番目の変数z^{(k)}_mを\muでラベル付されたある値に変更した状態を\bm{z}^{(k)}_{\mu}と書きます。これが先程説明した候補状態\bm{z'}^{(k)}に対応します。なお、\muに値によっては現在の値\bm{z}^{(k)}と等しくなることに注意してください。
この前提のもとで、状態\bm{z}^{(k)}が\bm{z}^{(k)}_{\mu}に更新される確率は以下のように計算します。
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年に提案された手法で、状態が更新されない確率である棄却率が平均的に最小となるような更新方法となっています。詳細は論文に書いてるのでそちらを読んでいただいた方がいいのですが、ここではこれまでの記法でどのように確率pが表現されるか見ておきます。
熱浴法の説明の際と同様に、更新する変数z^{(k)}_m\in\{L_m,L_m+1,\ldots,U_m\}の取りうる値に対して、その値が小さい順に\mu=1,2,\ldots,U_m - L_m + 1とラベル付けします。そして、\bm{z}^{(k)}において、m番目の変数z^{(k)}_mを\muでラベル付されたある値に変更した状態を\bm{z}^{(k)}_{\mu}と書きます。さらに、\bm{z}^{(k)}_{\mu}が現在の値\bm{z}^{(k)}と等しくなる\muを\lambdaと書きます。\bm{z}^{(k)}_{\lambda}=\bm{z}^{(k)}ということです。
この前提のもとで、状態\bm{z}^{(k)}が\bm{z}^{(k)}_{\mu}に更新される確率は以下となります。
p_{\mu} = \max\left[0, \min(\Delta_{\lambda,\mu}, w_{\lambda} + w_{\mu} - \Delta_{\lambda,\mu}, w_{\lambda}, w_{\mu})\right]
ここで、
\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番目の状態\mu=1が最も重みw_{\mu}が大きくないといけません。また、この更新確率はU_m - L_m = 1、つまり各変数が２状態しか取れない場合はメトロポリス法と同じ更新確率になります。
温度Tの下げ方
最後に温度Tの下げ方を説明します。ここでは最もよく使われているであろう幾何冷却について説明します。まず、sweepという単位を導入します。これは変数の数N回だけ変数の更新を試みることを1-sweepと定義します。そしてn-sweep目の温度をT_nと書きます。
幾何冷却とは以下の式に従って温度を下げることを指します。
毎回定数aをかけるだけです。このaはどのように決まるかというと、はじめに初期温度T_{\text{max}}と終温度T_{\text{min}}を決めておきます。そしてアルゴリズム全体でのsweep数を\text{num\_sweeps}とすればこれらの条件から
a = \left(\frac{T_{\text{min}}}{T_{\text{max}}}\right)^{\frac{1}{\text{num\_sweeps} - 1}}
と決まります。これら初期温度、終温度、sweep数をどう決めるかですが、これについては次回、実装して実験する際に考えることにします。
アルゴリズムのまとめ
ここまでで一応SAというアルゴリズムの概略は説明しました。最後にまとめておきます。
- 初期状態\bm{z}^{(0)}と初期状態T_{\text{min}}を決める。
- 次を\text{num\_sweeps}回行う。
2-1. 候補状態を生成し、その状態に確率pで更新する。このプロセスを変数の数Nだけ繰り返す。
2-2. 温度を下げる。
- 最終的に得られた状態を最適化問題の答えとする。
まとめと次回予告
この記事では整数変数で構成された最適化問題に対する疑似焼きなまし法(Simulated annealing: SA)のざっくりとした解説を行いました。
一応ここまでの知識で実際にSAを行うことができるはずですが、このアルゴリズムの中で最も重たい処理が遷移確率pの計算です。というのもここで説明した方法ではpを計算するために状態\bm{z}に対するエネルギーE(\bm{z})を計算する必要があります。この計算をそのまま行うと時間がかかってしまうので工夫を考える必要があります。具体的には少し式変形すると遷移確率pは状態を遷移した際のエネルギー差分\Delta E=E(\bm{z'}) - E(\bm{z})で表現することができます。そこで、このエネルギー差分をあらかじめ計算しておけば遷移確率pの計算を高速に行うことができます。ただし、このあらかじめ計算したエネルギー差分は、実際に状態が更新された際に改めて計算する必要があります。次回はここら辺の実装上の注意点を解説しようと思います。
追記: 続きはこちら。
募集してます！！！
Jijでは各ポジションを絶賛採用中です！
ご興味ございましたらカジュアル面談からでもぜひ！ご連絡お待ちしております。
数理最適化エンジニア（採用情報）
量子コンピューティングソフトウェアエンジニア（採用情報）
Rustエンジニア[アルゴリズムエンジニア]（採用情報）
研究チームPMO（採用情報）
オープンポジション等その他の採用情報はこちら
Discussion