2-アニーリングアルゴリズムの評価

アニーリングアルゴリズムはヒューリスティクスの一つなので毎回最適解を出せるとは限りません。近似解を出すためのアルゴリズムです。 また、確率的アルゴリズムなので毎回解が異なります。
そのためアルゴリズムの評価をする際は様々な平均値を持って解を評価します。よく用いられる手法が
  • 成功確率 : Success probability

  • TTS : Time to solution

  • 残留エネルギー : Resudial energy

の3つです。特にTTSは計算時間を表すので評価としてよく用いられます。残留エネルギーは最適解にどれくらい近づけたかを平均的に示す値です。

Time to solution

アニーリングアルゴリズムはどの計算時間でも何かしらの解を出すことができます。しかし、いくら計算が速くても、誤った解が求まるようだと、意味がありません。

そこで、必要な確率で最適解がでるための計算時間を指標とします。

“1-Introduction”でも行なったようにアニーリングアルゴリズムは複数回行って解を探すため、計算時間を評価するには複数回行うことも考慮に入れなければなりません。

複数回アニーリングを行うことも考慮に入れたのがTime to solution(TTS)です。TTSは簡単に導くことができます。

ここで1回のアニーリング時間を\(\tau\)とし、1回の\(\tau\)時間のアニーリングで最適解を得る確率を\(p_s(\tau)\)としましょう。この\(p_s(\tau)\)が、アルゴリズムの評価に用いられる成功確率となります。 すると\(R\)\(\tau\)時間のアニーリングを行って1回でも最適解を得る確率\(p_R\)

\[p_R = 1-(1-p_s(\tau))^R\]

となります。第二項は\(R\)回やっても最適解を得られない(間違える)確率なのでそれを1から引けば\(R\)回やって最適解を得られる確率になります。 この式を\(R\)について解きましょう。すると

\[R = \frac{\ln(1-p_R)}{\ln(1-p_s(\tau))}\]

となります。これに1回の計算時間\(\tau\)をかけるとトータルの計算時間となるのでそれがTime to solution (TTS)です。

\[TTS(\tau, p_R) = \tau R = \tau \frac{\ln(1-p_R)}{\ln(1-p_s(\tau))}\]

\(TTS(\tau, p_R)\)が示すのは計算時間が\(\tau\)の一回のアニーリングで確率\(p_s(\tau)\)で最適解が得られる時、確率\(p_R\)で最適解が得られるまでの繰り返し回数も考慮した総計算時間になります。

実際の計算の評価では、\(p_R\)を定数で与えます。\(p_R=0.99\)などにすることが多いかと思います。そして様々なアニーリング時間\(\tau\)\(p_s(\tau)\)を計算して、それらを使って\(TTS(\tau, p_R)\)を計算します。

OpenJij で実験を行ってTTSを測る

では以下では実際にOpenJijを使ってTTSを図ってみましょう。
1次元反強磁性モデルという問題を考えます。 1次元反強磁性イジングモデルというのは以下のエネルギーです。
\[H(\{\sigma_i\}) = \sum_{i=0}^{N-1} J_{i, i+1}\sigma_i \sigma_{i+1} + \sum_{i=0}^{N-1} h_i \sigma_i\]

ここで\(J_{ij} \in [0.1, 1.0]\)とし、\(h_0 = -1\)で他の縦磁場は0とします。\(J_{ij} > 0\)(反強磁性)なので各スピンは違う向きに向いた方がエネルギーが下がります。なので最適解となる:math:{sigma_i}\(1, -1, 1,\cdots\) のように値が交互になります。また\(h_0=-1\)としているので、最初のスピンは\(\sigma_0 =1\)になります。よって最適解は \(1, -1, 1, \cdots\)です。

つまり\(TTS\)\(1, -1, 1, \cdots\)を得るためにかかる総計算時間となります。
以下では上記のIsing モデルをといて、1回の計算時間を延ばすとTTSがどう変化するかを見てみましょう。
[ ]:
#!pip install -U cmake
#!pip install openjij
!pip show openjij  # openjij == 0.0.8 を仮定します
[1]:
import random
import numpy as np
import matplotlib.pyplot as plt
import openjij as oj
[2]:
# 反強磁性1次元イジングモデル を作る
N = 30
h = {0: -10}
J = {(i, i+1): 1 for i in range(N-1)}

TTSを計算してみる

TTSを計算してみます。openjijのsample_isingまたはsample_qubo から返ってくる Response クラスは info というメンバ変数を持っています。
infoにはsamplerごとに異なった情報が辞書型で格納されています。ほとんどのSamplerには ['execution_time']という名前で各アルゴリズムの1回の実行時間 (μs) が格納されています。

SASampler の場合は Simulated annealing の1回あたりの計算時間を格納しています。

[3]:
# 最適解
correct_state = [(-1)**i for i in range(N)]

# TTS を計算するのに必要なp^R
pR = 0.99

# Samplerの引数のstep_num というパラメータに渡すリスト(step_num_list)
# step_num はアニーリング中のパラメータ(温度, 横磁場)を下げていくときの分割数
# なので増やせば増やすほどゆっくりアニーリングすることになってアニーリング時間が伸びる
step_num_list = list(range(10, 51, 10))  # [10, 20, 30, 40, 50]

TTS_list = []  # 各計算時間に対するTTSを格納しておくリスト
tau_list = []  # 計算時間を格納しておくリスト

#  計算の過程で、成功確率が求まるので、ついでにリストにしておきましょう
ps_list = []  # 成功確率を格納しておくリスト

iteration = 100  # 確率を計算するために1回のアニーリングを行う回数

for step_num in step_num_list:
    sampler = oj.SASampler(num_sweeps=step_num, iteration=iteration)
    response = sampler.sample_ising(h, J)

    # 返ってきた解であっている状態の数を数えて最適解を得た確率を計算する
    tau = response.info['execution_time']
    ps = sum([1 if state == correct_state else 0 for state in response.states])/iteration

    # ps=0だとTTSが無限大になってしまうのでそこは回避
    if ps == 0:
        continue

    # TTSを計算する
    TTS_list.append(np.log(1-pR)/np.log(1-ps)*tau)
    tau_list.append(tau)

    # 成功確率を計算する
    ps_list.append(ps)
[4]:
fig, (axL, axR) = plt.subplots(ncols=2, figsize=(15,3))
plt.subplots_adjust(wspace=0.4)

fontsize = 10
axL.plot(tau_list, TTS_list)
axL.set_xlabel('annealing time', fontsize=fontsize)
axL.set_ylabel('TTS', fontsize=fontsize)

axR.plot(tau_list, ps_list)
axR.set_xlabel('annealing time', fontsize=fontsize)
axR.set_ylabel('Success probability', fontsize=fontsize)
[4]:
Text(0, 0.5, 'Success probability')
../_images/ja_2-Evaluation_12_1.svg

TTSと成功確率を計算することができました。 横軸は、1回のアニーリングにかかる計算時間となっています。 アニーリング時間が増加すると、TTSも増加する傾向にあります。 しかし、必要な成功確率が補償できた時点で、アニーリングをやめるなどできるので、これらの指標は有効に使えそうですね!

OpenJijにはデフォルトでTTS, 残留エネルギー, 成功確率を評価してくれる関数が付いています。
openjij.solver_benchmark です。
次は、このベンチマーク関数を用いて、これらを求めてみたいと思います。

solver_benchmark 関数

solver_benchmark 関数では TTSと残留エネルギー、成功確率を計算してくれます。

引数は以下です

  • solver: function Response クラスを返す関数です。time と iteration という引数を持つ必要があります。 time は計算時間を制御するパラメータにしてください。SASamplerの場合は step_num に相当します。 iteration は TTSや残留エネルギーなどを計算するときに必要なサンプリング回数を指定します。 また 関数の戻り値のResponse.info には、time_nameという引数で指定する文字列をキーワードで持ち、time_nameに紐づく値は1回あたりの計算時間が格納されている必要があります。

  • time_list: list solver の time 引数に入れる値のリストです。

  • solutions : list(list: state) 基底状態となる状態のリストです。 縮退している場合は[state1, state2]のようにして複数入れます。

  • args : dict solver にオプションで必要な場合につけます。デフォルトでは args = {} です。

  • p_r: 0 < float <= 1 TTSを計算するために必要な値です。上記のTTSの説明のp_R に相当します。

  • ref_energy: float 参照エネルギー。次のmeasure_with_energyと合わせて使います。デフォルトでは ref_energy = 0 です。

  • measure_with_energy :bool False : スピンの状態が基底状態と一致していれば、成功とカウントします。 True : エネルギーが先ほどのref_energy以下だと、成功とカウントします。つまり、基底状態がわからない場合に使えます。 デフォルトでは、Falseです。

  • time_name: str Response.info に入っている実行時間と紐づくキーを指定します。デフォルトでは"execution_time"です。

返り値はbenchmark計算の結果が以下のように辞書型で格納されています。

  • time: 実行時間のリスト。

  • success_prob: 成功確率のリスト。

  • tts: ttsのリスト。

  • residual_energy: 残留エネルギーのリスト。

  • info: (dict) ベンチマーク関数のパラメータ情報。

[ ]:
# 最適解
correct_state = [(-1)**i for i in range(N)]

# ステップ数とアニーリングの反復数を与えます
step_num_list = list(range(10, 51, 10))  # [10, 20, 30, 40, 50]
iteration = 100

# benchmark 関数で TTS 残留エネルギー 成功確率を計算
result = oj.solver_benchmark(
                      solver=lambda time, **args: oj.SASampler(step_num=time, iteration=iteration).sample_ising(h,J),
                      time_list=step_num_list, solutions=[correct_state], p_r=0.99
            )
[ ]:
fig, (axL,axC,axR) = plt.subplots(ncols=3, figsize=(15,3))
plt.subplots_adjust(wspace=0.4)

fontsize = 10
axL.plot(result['time'], result['tts'])
axL.set_xlabel('annealing time', fontsize=fontsize)
axL.set_ylabel('TTS', fontsize=fontsize)

axC.plot(result['time'], result['residual_energy'])
axC.set_xlabel('annealing time', fontsize=fontsize)
axC.set_ylabel('Residual energy', fontsize=fontsize)

axR.plot(result['time'], result['success_prob'])
axR.set_xlabel('annealing time', fontsize=fontsize)
axR.set_ylabel('Success probability', fontsize=fontsize)
Text(0, 0.5, 'Success probability')
../_images/ja_2-Evaluation_15_1.png
このようにしてベンチマークをとることができます。
先ほどと同様、反強磁性1次元イジングモデルを用いていますので、TTSと成功確率は、同様に単調増加のグラフが出力されたのではないでしょうか? (ヒューリスティックスな解法なので、厳密な計算結果は毎回変わります。)
残留エネルギーについてもアニーリング時間を増加すれば、いつかは収束しそうですね。

solver に入れる関数は Responseクラスを返して、.info に計算時間を格納しておけばどのような関数でもよいです。適当なやつを作ってみましょう。

次は、[1,1,1,…]と[1,-1,1,…]と[-1,1,-1,…]の3つのスピン状態から、ランダムに一つの状態を返す関数です。最適解を[1,-1,1,…]として、アニーリングアルゴリズムのベンチマークを取ってみましょう。

[ ]:
import time
def anti_ferro_solver(time_param, iteration, h, J):
    """
    すべて 1 と [1,-1,1,...] と [-1,1,-1,...] の3つの状態からランダムに選ぶ
    """

    # 入力された h と J から添字の集合をつくる
    indices = set(h.keys())
    indices = list(indices | set([key for keys in J.keys() for key in keys]))

    ones_state = list(np.ones(len(indices), dtype=int))        # all 1

    minus_plus_state = np.ones(len(indices), dtype=int)
    minus_plus_state[::2] *= -1                                         # -1, 1, -1, 1, ...
    plus_minus_state = -1 * minus_plus_state                   # 1, -1, 1, -1

    start = time.time()
    _states = [ones_state, list(minus_plus_state), list(plus_minus_state)]
    state_record = [_states[np.random.randint(3)] for _ in range(iteration)]   # 3つの状態からランダムにひとつ選ぶ
    exec_time = (time.time()-start) * 10**6 * time_param                                # 適当に計算時間を伸ばしておきます

    energies = [sum(state) for state in state_record]                                     # エネルギーの計算はてきとうです


    # Responseクラスに状態とエネルギーを格納します
    response = oj.Response(indices=indices, var_type='SPIN')
    response.update_ising_states_energies(state_record, energies)
    response.info['execution_time'] = exec_time

    return response
[ ]:
# 最適解
correct_state = [(-1)**i for i in range(N)]  # 1, -1, 1, -1

# ステップ数とアニーリングの反復数を与えます
time_param_list = list(range(10, 51, 10))  # [10, 20, 30, 40, 50]
iteration = 2000

# benchmark 関数で TTS 残留エネルギー 成功確率を計算。
result = oj.solver_benchmark(
                      solver= lambda time_param, **args: anti_ferro_solver(time_param, iteration, h, J),
                      time_list=time_param_list, solutions=[correct_state], p_r=0.99
              )
[ ]:
fig, (axL,axC,axR) = plt.subplots(ncols=3, figsize=(15,3))
plt.subplots_adjust(wspace=0.4)

fontsize = 10
axL.plot(result['time'], result['tts'])
axL.set_xlabel('annealing time', fontsize=fontsize)
axL.set_ylabel('TTS', fontsize=fontsize)

axC.plot(result['time'], result['residual_energy'])
axC.set_xlabel('annealing time', fontsize=fontsize)
axC.set_ylabel('Residual energy', fontsize=fontsize)

axR.plot(result['time'], result['success_prob'])
axR.set_xlabel('annealing time', fontsize=fontsize)
axR.set_ylabel('Success probability', fontsize=fontsize)
Text(0, 0.5, 'Success probability')
../_images/ja_2-Evaluation_19_1.png
今、最適解は[1, -1, 1, -1, …] です。
3つからランダムにひとつ1つ正解の状態を選んでいるので成功確率はだいたい1/3くらい。
上図右の成功確率もだいたいそれくらいになっています。

このようにうまくsolver関数を作っておけばOpenJijのソルバーに限らず TTSと残留エネルギーと成功確率を計算してくれます。