Simulated Annealing(焼きなまし法)の実装

あくまで研究などとは分野外なので素人ですが、現在は趣味でイジングモデルに基づいた、臨界点で複雑性を最大化するような手法を自作でコーディングしています。その過程で量子アニーリング手法などの原点となっているSimulated annealingを実装するのも悪くないと思い、ウェブサイトなどを参考にして実装した備忘録になります。

Simulated Annealing(焼きなまし法)とは

最適化問題を解くための手法で、ヒューリスティックに解を求めて行くのは良くある基本的な勾配法と同じです。しかし、温度パラメータTを導入していることが大きな特徴になっています。エネルギー変化を-Tで割ったものをexpで括った値である確率p次第では、もしも変更後のコストの値が前に算出したコストより悪かったとしても更新する方策をとります。一見損をしている手法に思えますが、このような手法を取ることによって、更新の際に局所解に陥ることを防いでくれる役割を持っています。 $p=\exp(-\frac{\Delta E}{T})$

この手法の重要なパラメータTですが、通常は繰り返し更新すればするほど、減少するように設定されています。これは、最初の段階はコストが悪くても更新する確率pが高いですが、更新するたびに徐々にpの値が小さくなることを意味しています。

実装

実装の際に使用したコスト関数は、以下のサイトに倣ってHimmelblau's functionを使うことにしました。これは、最適化アルゴリズムのパフォーマンスを評価するための関数であり、最小値の位置は分析的に見つけることが可能になっています。

orako-column.com

今回は、上のサイトなどを参考にして3Dに拡張してみました。

import random
import numpy as np
from matplotlib import cm
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt

def cost_function(x, y=0):
    """Himmelblau's function"""
    return ((x**2) + y -11)**2 + (x + (y**2) - 7)**2

INITIAL_STATE = -5 # Initial position
T = 10000 # Temperature
Tmin = 0.0001 # Minimum temperature
STEP = 0.1 # How change variable
COOL = 0.999 # Cooling variable

prev_state = np.array([INITIAL_STATE,INITIAL_STATE])
prev_energy = cost_function(prev_state, 0)
state = np.zeros(2)
while T > Tmin:
    a=0
    #d_state = np.random.randn(2)
    d_state = np.random.beta(1, 1, 2)
    for i in d_state:
        if i > 0.5:
            d_state[a] = STEP
        else:
            d_state[a] = -STEP
        a+=1
    a=0
    state = prev_state + d_state
    energy = cost_function(state) # Calculate cost
    p = np.exp(-np.abs(energy - prev_energy) / T) # Calculate p
    # Update in a case of probability p or smaller energy
    if (energy[0] < prev_energy[0] or random.random() < p[0]):
        if (energy[1] < prev_energy[1] or random.random() < p[1]):
            prev_state = state
            prev_energy = energy
    T = T * COOL


# plot
fig = plt.figure(figsize=(8,8))

ax = fig.add_subplot(111, projection='3d')
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y) 
Z = cost_function(X, Y)
plt.title("Himmelblau's function")
plt.xlabel("x")
plt.ylabel("y")
ax.plot_surface(X, Y, Z, cmap=cm.coolwarm, linewidth=0, rstride=1, cstride=1, alpha=.5)
cset = ax.contour(X, Y, Z, zdir='z', offset=-5, cmap=cm.coolwarm)
#cset = ax.contour(X, Y, Z, zdir='x', offset=-5, cmap=cm.coolwarm)
#cset = ax.contour(X, Y, Z, zdir='y', offset=5, cmap=cm.coolwarm)

# Initial end points in 3D space (x,y,z)
mi = (INITIAL_STATE, INITIAL_STATE, cost_function(INITIAL_STATE, INITIAL_STATE))
print(mi)
me = (prev_state[0], prev_state[1], cost_function(prev_state[0], prev_state[1]))
print(me)

#plot points.
res = np.array([mi,me])
print(res)
ax.plot(res[:,0], res[:,1], res[:,2], marker="o", ls="", c=cm.coolwarm(1.))
ax.view_init(azim=-45, elev=15)
plt.show()

結果

f:id:YuKR:20191209090952p:plain

最初の配列が初期位置とコスト関数の値、次の行の配列が終点の位置とコスト関数の値で、それらを合わせた行列が最後の出力結果になっています。 3Dプロットを見ていただくとお分りいただけると思うのですが、終点が最小値の位置に近い値に落ち着いています。STEPやCOOLの値や、何回か試行しても色々値が変わるので、興味がある人は弄ってパフォーマンスを確認して頂ければと思います。

まとめ

今回はSimulated Annealingの実装をちょっとやってみました。イジングモデルを元にした時系列解析手法も作っているので、機会があればそちらも公開できたらと思っています。

たけのこブログ

凡人が頑張って背伸びするブログ

Simulated Annealing(焼きなまし法)の実装

Simulated Annealing(焼きなまし法)とは

実装

結果

まとめ