第6回: アニーリング基礎5.補助スピン活用基礎1（ナップザック問題）とTTSの活用

-前回、容量$C$のカバンに、夫々値段$p_i$と大きさ$v_i$が定義される$N$個のお菓子を、出来る限り値段が高くなるように、最大量つめる場合、どのような組合せにするか調べる問題を、アニーリングを用いて実施したが、今回は更に補助スピンを使用した計算、並びにTTSを用いた計算精度比較まで行う。
前回：第2回: アニーリング基礎2.ナップザック問題、Optunaによるハイパーパラメタ最適化
Keyword：補助スピン、TTS、ナップザック問題

1. 基本準備/設定

1.1 量子アニーリング用モジュールのインストール

まず、初めに、以下のモジュールをインストールする。

• dwave-ocean-sdk：Dwaveアニーリングマシンのモジュール
• openjij：アニーリングシュミレーター
• pyqubo：Qubo行列作成サポートモジュール
• optuna：ハイパーパラメタの最適化

pip install dwave-ocean-sdk
pip install openjij
pip install pyqubo
pip install optuna

1.2 Dwaveの初期設定

続いて、Dwaveのtoken、enpoint(アクセスURL)、使用するマシンの指定まで実施。併せて、Jijシュミレーターも設定する。

###基本設定
token = "DEV-73ac79954ef21XXXXXXXXXXXXXXXXXXXXXX"
endpoint = "https://cloud.dwavesys.com/sapi/"
solver = "Advantage2_prototype1.1"
###Dwaveの設定
from dwave.system import DWaveSampler, EmbeddingComposite
dw_sampler = DWaveSampler(solver=solver, token=token, endpoint = endpoint)
###Jijシュミレーターの設定
from openjij import SQASampler
sampler_jij = SQASampler()

2. 定式化

2.1 振り返り

前回の「第2回: アニーリング基礎2.ナップザック問題、Optunaによるハイパーパラメタ最適化」では、以下のような問題を解いた。
***「夫々の値段が$p_i$、重さが$w_i$のお菓子$N=8$個を、容量$C=300$の鞄に詰めるとき、お菓子があふれない範囲で値段を最大化する組合せを見つけたい」***夫々の数値は以下の通り。

N, C = 8, 300
p = [80,110,100,90,60,75,40,25]
w = [100,160,120,100,80,70,40,30]

スピン$x_i$は$i$番目のお菓子を入れる($x_i =1$) or 入れない($x_i =0$)と定義して、以下のようにハミルトニアンを定式化した。

H(x) = E(x) + \lambda_1C_1(x) \\
H(x) = -\sum_{i=0}^Nx_i p_i + \lambda_1(C-\sum_{i=0}^Nx_i w_i)^2 \\

この時、コスト関数の第2項 $(C-\sum_{ }x_i w_i)^2$は$C=300$に近づく程エネルギーが小さくなることを意味するが、$\sum_{ }x_i w_i=290$でも$\sum_{ }x_i w_i=310$でも同じ、$10×10=100$のエネルギーを示すことになり、本来罰則を与えるべきではない、$\sum_{ }x_i w_i=290$にも罰則を与えてしまっている為、対策が必要。
具体的に、$0≦\sum_{ }x_i w_i≦300$の範囲では、エネルギーが出来る限り小さくないようなコスト関数に仕立てる必要がある。これにより、$\sum_{ }x_i w_i=310$の時、エネルギーは100になるが、$\sum_{ }x_i w_i=290$の時、エネルギーが略0になるとしたい。
このように制約式を不等号で設定するには、第2のスピンである補助スピン$y_k$を使用することになる。詳細は以下で議論するが、$y_k$を使用することで、Qubo行列は$y_k$の数だけ増大する点に注意。補助スピン$y_k$は出来る限り少なくする必要がある。

2.2 補助スピンyの使い方

まず補助スピン$y_k$は大きく2種類に分類される。

• 独立型: $y_k$は$x_i$と完全に独立。つまり、$y_k$と$x_i$間で関係式は成り立たない。独立型はコスト関数を不等号で示す場合に使用する。
• 関連型: 例えばハミトルニアンが$x_i×x_n×x_p$のように3次式になった場合、計算不可になるので、$y_k$=$x_n×x_p$などと定義して、$x_i×y_k$とすることで、無理矢理2次式にする場合。通常、$x$と$y_k$の関係性を担保する条件をコスト関数に加えることになり、問題が一挙に複雑化する。

今回は範囲指定する問題なので、よりシンプルな独立型で補助スピンを使用する。尚、より複雑な関連型は別途検討する。

まず、下準備として、$W = \sum_{ }x_i w_i$が取り得る値を列挙する。通常$n$が膨大な、本来アニーリングが対象とする問題では、取り得る値を列挙することは不可能であるが、本問題は$n=8$であり、組合せは$2^8=256$通りのみなので、無理矢理計算してみる。

###intertoolsで[0,1]*8の組み合わせを抽出する。
import itertools
l01 = [0,1]
combination = list(itertools.product(l01,l01,l01,l01,l01,l01,l01,l01))

###続いて重さの組み合わせリスト作成してソートする
import pandas as pd
import numpy as np
weight_list =np.array([np.dot(w,_comb) for _comb in combination])
weight_list =np.sort(weight_list)
#print(weight_list)

###続いて重複を除外する
weight_list_unique = np.unique(weight_list)
weight_list_unique

###続いて容量C=300以下のみを抽出する
under_C = list(np.delete(weight_list_unique,np.where(weight_list_unique>C)))#300以上を削除した
print("データ数: "+str(len(under_C)), under_C)

<Output>
データ数: 26 [0, 30, 40, 70, 80, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250, 260, 270, 280, 290, 300]

コスト関数の範囲は$0≦\sum_{ }x_i w_i≦300$としたいが、理解しやすくする為に、$250≦\sum_{ }x_i w_i≦300$の範囲とする場合のコスト関数の設定を考える。
まず、取り得る値リストを確認すると、上記範囲内に収まる重さは$w_L=[250, 260, 270, 280, 290, 300]$の6パータンとなる。$250≦\sum_{ }x_i w_i≦300$だからと言って、$[250, 251, 252…]$と考える必要はない。
$\sum_{ }x_i w_i=260$の時、コスト関数が0になるような設計をするわけだが、まず、$K=6$個の補助スピン$[y_1,y_2,y_3,y_4,y_5,y_6]$を設定する。勿論、$y_k$はスピンなので、取り得る値は0 or 1。この6個のスピンを用いて以下を表現したい。

$y_k$の状態	表現したいこと
$y_1$のみ「1」	250
$y_2$のみ「1」	260
$y_3$のみ「1」	270
$y_4$のみ「1」	280
$y_5$のみ「1」	290
$y_6$のみ「1」	300

これは以下の式で表現できる。

Target = \sum_{k=1}^Ky_k w_L (ただし、y_kは1つのみ「1」をとる)\\

例えば、260を表現したい時、以下のようになる。

H(x) = [0,1,0,0,0,0]×[250, 260, 270, 280, 290, 300] = 260

ここで、上の式の「ただし、$y_k$は1つのみ「1」をとる」であるが、これは制約条件によって表現することが出来る。

C_{sub1}(y) = \lambda_{sub1}(\sum_{k=1}^Ky_k-1)^2\\

一方で、コスト関数は、以下の従来のコスト関数のCをTargetに置き換え、さらに「ただし、$y_k$は1つのみ「1」をとる」を表現する$C_{sub1}(y)$を追加すればよい。

C_1 = \lambda_1(C-\sum_{i=0}^Nx_i w_i)^2 \\
C_1 = \lambda_1(\sum_{k=1}^Ky_k w_L-\sum_{i=0}^Nx_i w_i)^2+\lambda_{sub1}(\sum_{k=1}^Ky_k-1)^2\\

定式は上記で完了だが、これが意味する所を考えてみる。$\sum_{ }x_i w_i=260$の時、$\sum_{ }y_k w_L$を動かしてコスト関数$C_1$がどのように変化するか確認してみる。

### wLの設定
weight = 260
wL=[250,260,270,280,290,300]
_lambda1 = 1

###夫々の場合のC1を計算
for _k,_wL in enumerate(wL):
  _c1 = _lambda1*(_wL-weight)**2
  print(f"y{_k}のみ1の時、C1は{_c1}となる")

<Output>
y0のみ1の時、C1は100となる
y1のみ1の時、C1は0となる
y2のみ1の時、C1は100となる
y3のみ1の時、C1は400となる
y4のみ1の時、C1は900となる
y5のみ1の時、C1は1600となる

計算結果を確認すると、C1が6パターン出てくるので頭が混乱するが、アニーリングはエネルギーが最小になるパターンが抽出されるので、上から2点目の、$y_1$のみ「1」の時の$C_1=0$が採用されることになると理解できる。
一方で、今回$w_L$は250～300の値をとるが、$\sum_{ }x_i w_i=240$の時、240に最も近い$w_L=250$が採用されて、$C_1=100$となる。一方で、$\sum_{ }x_i w_i=310$の時、310に最も近い$w_L=300$が採用されて、$C_1=100$となり、先程と等価になる。従い、$0≦\sum_{ }x_i w_i≦300$の範囲に設定したいが、その場合、補助スピン$y_k$を26個準備する必要があるが、計算量が増大するので避けたい。
一方で、エネルギー式自体は指定の範囲内で計算される為に、あまりシビアに範囲設定をしなくても良い。実際に、第2回: アニーリング基礎2.ナップザック問題、Optunaによるハイパーパラメタ最適化」で説いた時は$C=300$と範囲設定しなかったが、ある程度の結果は得られた。
そこで、計算数を減らしつつ、範囲を広げる為にはやはり$w_L=[260,270,280,290,300]$などと設定すれば良い。
最後に、今回は$n=8$であったので、全ての場合を計算した上で$w_L$を設定したが、実問題ではある程度狙いを定めて色々と数値を試しながら範囲を設定していくことになる。

2.3 定式化

それでは、上で整理した結果を踏まえて改めて定式化行うと以下となる。

H(x) = E(x) + \lambda_1C_1(x) \\
H(x) = -\sum_{i=0}^Nx_i p_i+
\lambda_1(\sum_{k=1}^Ky_k w_L-\sum_{i=0}^Nx_i w_i)^2+
\lambda_{sub1}(\sum_{k=1}^Ky_k-1)^2 \\
ここで、\\
w_{L} = [220,230,240,250,260,270,280,290,300], K=9 とする。

3. 実装

3.1 Qubo行列の作成からDwaveでの計算

再掲になるが問題並びにスピンを設定する。

### 問題の設定
N, C, K = 8, 300, 9
p = [80,110,100,90,60,75,40,25]
w = [100,160,120,100,80,70,40,30]
wl = [220,230,240,250,260,270,280,290,300]

### スピンの設定
from pyqubo import Array, Placeholder
x = Array.create(name="x", shape =N, vartype="BINARY") 
y = Array.create(name="y", shape =K, vartype="BINARY")

### lamdaの設定
lambda_1 = Placeholder("lambda_1")
lambda_s1 = Placeholder("lambda_s1")

今回は、補助スピンを使用する場合、しない場合の双方の定式化を行う。

H_{w.s}(x) = -\sum_{i=0}^Nx_i p_i+
\lambda_1(\sum_{k=1}^Ky_k w_L-\sum_{i=0}^Nx_i w_i)^2+
\lambda_{s1}(\sum_{k=1}^Ky_k-1)^2 \\
H_{w.o.s}(x) = -\sum_{i=0}^Nx_i p_i
+ \lambda_1(C-\sum_{i=0}^Nx_i w_i)^2 \\

まず、should beの数値を確認する為、補助スピンを使用しない場合を解いてみる。尚、$\lambda_1$は「第2回: アニーリング基礎2.ナップザック問題、Optunaによるハイパーパラメタ最適化」にて455.068が最適数であったことから、これを流用する。

### 補助スピン無
##定式化
#第1項
Hwos = 0
for i in range(N):
  Hwos = Hwos - x[i]*p[i]
#第2項
_x = 0
for i in range(N):
  _x = _x + x[i]*w[i]
Hwos = Hwos + lambda_1*(C-_x)**2 
##Qubo行列抽出
model_wos = Hwos.compile()
Qubo_wos, offset = model_wos.to_qubo(feed_dict={"lambda_1":455.068})
##Dwaveでの計算
sampler_dw = EmbeddingComposite(dw_sampler)
solution_wos = sampler_dw.sample_qubo(Qubo_wos, num_reads=30)

続いて、結果の評価と可視化関数を定義する。

def view(_solution,n=1000):
  #Dwaveの結果の数だけ繰り返す
  _result = []
  for _i in range(len(_solution.record)):
    #補助スピンの有無で場合分けし、条件満たす場合だけ[重量、値段]を返す
    _temp = [0,0,0]
    _price, _weight = np.dot(_solution.record[_i][0][0:N],p), np.dot(_solution.record[_i][0][0:N],w)
    _temp[1],_temp[2] = _weight, _price
    if len(_solution.record[0][0]) == N:#補助スピンを使用しない
      if _weight <=C:
        _temp[0] = 1
    else:
      if sum(_solution.record[_i][0][N+1:N+K])==1:
        if _weight <=C:
          _temp[0] = 1
    #Output
    _result.append(_temp)
  return _result[:n]
print(view(solution_wos,10))

<Output>
[[1, 300, 240], [1, 300, 230], [0, 310, 275], [1, 290, 255], [1, 290, 255], [0, 310, 235], [0, 320, 270], [0, 320, 250], [1, 280, 210], [0, 330, 305]]

260前後が最適値のようである。続いて補助スピン有の場合を検証する。

### 補助スピン有
##定式化
#第1項
Hws = 0
for i in range(N):
  Hws = Hws -x[i]*p[i]

#第2項
_x, _y = 0, 0
for k in range(K):
  _y = _y + y[k]*wl[k]
for i in range(N):
  _x = _x + x[i]*w[i]
Hws = Hws + lambda_1*(_y-_x)**2

#第3項
_y = 0
for k in range(K):
  _y = _y+y[k]
Hws = Hws + lambda_s1*(_y-1)**2

##Qubo行列抽出
model_ws = Hws.compile()
Qubo_ws, offset = model_ws.to_qubo(feed_dict={"lambda_1":5,"lambda_s1":10000})

##Dwaveでの計算
sampler_dw = EmbeddingComposite(dw_sampler)
solution_ws = sampler_dw.sample_qubo(Qubo_ws, num_reads=30)

##結果確認
print(view(solution_ws,10))

<Output>
[[0, 280, 230], [0, 420, 350], [1, 250, 215], [1, 290, 265], [1, 210, 205], [0, 540, 440], [0, 480, 400], [1, 200, 150], [0, 430, 385], [0, 430, 385]]

計算結果を見比べてみると補助スピンを考慮したことでかえって結果が悪くなっている模様。ただし、ハイパーパラメタを最適化していない。従い、Optunaによる最適化を実行する。

3.2 Optunaによるハイパーパラメタの最適化

まず、アニーリングの計算を行う関数を定義。

def a_cal(_Qubo, machine = "Dwave", case =100):
### Dwaveかどうかの選択
  if machine == "Dwave":
    sampler_dw = EmbeddingComposite(dw_sampler)
    _solution = sampler_dw.sample_qubo(_Qubo, num_reads=case)
    
### Jijかどうかの選択
  elif machine =="Jij":
    sampler_jij = SQASampler()
    _solution = sampler_jij.sample_qubo(_Qubo,num_reads=case)
  return(_solution)

続いて、Optunaを用いてハイパーパラメタの最適化を実施する。補助スピン無の場合では、既に最適化出来ているので補助スピン有の場合の最適化なパラメタを算出する。

###ライブラリの呼び出し
import optuna

###関数の定義
def objective(trial):

  ###ハイパーパラメタの設定
  _lambda1, _lambda2 = trial.suggest_uniform('lambda1', 0.001, 10), trial.suggest_uniform('lambda2', 100, 100000)

  ###コンパイルとQuboの取り出し
  _model = Hws.compile()
  _Qubo, offset = _model.to_qubo(feed_dict={"lambda_1":_lambda1,"lambda_s1":_lambda2})

  ###a_calでアニーリング計算実施し、それをrateで正答率にする。
  _scoure = 0
  _solution = a_cal(_Qubo, "Dwave")
  for _i in range(15):
    _temp = view(_solution)[_i]
    _scoure = _scoure + _temp[0]*_temp[1]*_temp[2]
    _scoure = -1*_scoure
  return _scoure

### Optunaでの最適化
study = optuna.create_study()
study.optimize(objective, n_trials=50)

### Optunaでの最適解を表示
best_lambda_1, best_lambda_s1= [i for i in study.best_params.values()][0], [i for i in study.best_params.values()][1]
print(f"lambda1: {best_lambda_1}, lambda2: {best_lambda_s1}")

<Output>
[I 2022-12-01 09:55:21,048] A new study created in memory with name: no-name-db81f3d6-e2a7>:8: FutureWarning: suggest_uniform has been deprecated in v3.0.0. This feature will be removed in v6.0.0. See . Use :func:`~optuna.trial.Trial.suggest_float` instead.
  _lambda1, _lambda2 = trial.suggest_uniform('lambda1', 0.001, 10), trial.suggest_uniform('lambda2', 100, 100000)
[I 2022-12-01 09:55:22,881] Trial 0 finished with value: -71050.0 and parameters: {'lambda1': 7.931853736284838, 'lambda2': 66037.35469172239}. Best is trial 0 with value: -71050.0.
[I 2022-12-01 09:55:23,612] Trial 1 finished with value: -103300.0 and parameters: {'lambda1': 0.8377312197086793, 'lambda2': 59787.12989846234}. Best is trial 1 with value: -103300.0.
[I 2022-12-01 09:55:24,971] Trial 2 finished with value: 133900.0 and parameters: {'lambda1': 6.1098598745932176, 'lambda2': 67683.42029376556}. Best is trial 1 with value: -103300.0.
[I 2022-12-01 09:55:26,249] Trial 3 finished with value: 0.0 and parameters: {'lambda1': 7.520078359350096, 'lambda2': 70730.92471344008}. Best is trial 1 with value: -103300.0.
[I 2022-12-01 09:55:27,080] Trial 4 finished with value: 0.0 and parameters: {'lambda1': 9.790003763225371, 'lambda2': 31893.623278203253}. Best is trial 1 with value: -103300.0.
...
[I 2022-12-01 09:56:04,905] Trial 49 finished with value: 27500.0 and parameters: {'lambda1': 1.5117423748222834, 'lambda2': 77285.09710571305}. Best is trial 21 with value: -180250.0.
lambda1: 1.1086266828212719, lambda2: 98964.87401621634

得られた最適化なパラメタを用いて再度計算行いそれぞれを見比べてみる。

##補助スピン有
model_ws = Hws.compile()
Qubo_ws, offset = model_ws.to_qubo(feed_dict={"lambda_1":best_lambda_1,"lambda_s1":best_lambda_s1})
sampler_dw = EmbeddingComposite(dw_sampler)
solution_ws = sampler_dw.sample_qubo(Qubo_ws, num_reads=30)

##補助スピン無
model_wos = Hwos.compile()
Qubo_wos, offset = model_wos.to_qubo(feed_dict={"lambda_1":455.068})
sampler_dw = EmbeddingComposite(dw_sampler)
solution_wos = sampler_dw.sample_qubo(Qubo_wos, num_reads=30)

##結果確認
print(f"補助スピン有:{view(solution_ws,10)}")
print(f"補助スピン無:{view(solution_wos,10)}")

<Output>
補助スピン有:[[0, 330, 295], [0, 300, 260], [0, 310, 275], [0, 270, 235], [0, 220, 200], [0, 320, 290], [0, 190, 175], [0, 190, 175], [0, 400, 350], [0, 220, 200]]
補助スピン無:[[1, 300, 250], [1, 300, 250], [0, 310, 285], [0, 310, 275], [1, 290, 265], [0, 310, 245], [0, 320, 290], [0, 320, 280], [0, 320, 270], [0, 320, 250]]

得られた結果を見比べると、補助スピン無の方が良い結果を示している。補助スピンを使用すると制約条件が増加するのに加えて、勿論$\lambda$の設定が必要になる為、例題のような小規模な問題には適さないのかもしれない。もしくは、Optunaによる$\lambda$の設定を改良する必要がありそうである。
一方で、今回補助スピンを用いたが、Qubo行列がどのような形になっているか不透明なので、確認してみる。

3.3 Qubo行列の形の確認

Pyqubo用いて抜き出すQubo行列は辞書型になっている為、これを行列に変換し、更にSeabornを用いて可視化する関数を定義する。

###Pyquboの辞書型のQubo行列をリストに変換する
import seaborn as sns

def transform_qubo(_n,_k,_target_qubo,_show=False):

  ##辞書型の全座標作成
  _all_dict=[]
  for _i in range(_n+_k):
    if _i >=_n:
      _i = _i - _n
      _le_i = "y"
    else:
      _le_i = "x"
    for _j in range(_n+_k):
      if _j >=_n:
        _j = _j - _n
        _le_j = "y"
      else:
        _le_j = "x"
      _all_dict.append((f"{_le_i}[{_i}]",f"{_le_j}[{_j}]"))

  ##行列型のリスト作成
  _all_list = [ ]
  for _i in range(_n+_k):
    for _j in range(_n+_k):
      _all_list.append([_i,_j])

  ##Quboへ変更
  _qubo = np.zeros((_n+_k)**2).reshape(_n+_k,_n+_k)
  for _k in range((_n+_k)**2):
    try:
      _qubo[_all_list[_k][0],_all_list[_k][1]]=_target_qubo[_all_dict[_k]]
    except:
      _qubo[_all_list[_k][0],_all_list[_k][1]]=0

  #可視化
  if _show == True:
    _qubo_mask = (_qubo==0)
    sns.heatmap(_qubo, square=True,robust=True, linewidths=0.01,linecolor="Blue",cmap="Reds",mask = _qubo_mask)  

  else:
    return _qubo

補助スピン有のQubo行列を可視化する。

transform_qubo(N,K,Qubo_ws,True)

補助スピン無のQubo行列を可視化する。

transform_qubo(N,0,Qubo_wos,True)

補助スピンを使用した場合、計算量は追加したスピンの2乗で増大する。つまり8x8のQubo行列に補助スピンを9個追加すれば、今回のように17x17になることが分かる。
また、今回は独立型の補助スピンを用いた為、スピン$x_i$とスピン$y_k$の間に関係性は確認されない。特に$y_k$は1つだけ「1」とするには強い重みを与える必要があり、スピン$y_k$の箇所に特に大きい数値を示す結果となった。

4. TTSによる評価

これまでの検証では、アニーリングの結果精度確認は、独自にEvaluation関数等を定義して評価してきた。具体的には、制約条件満たさない場合は「0」、満たす場合は最大化/最小化したい対象の計算結果を複数回算出してその合計を評価値としてきた。
一方で正式な評価方法としてTime To Solution/TTSなるものがある。これは、「正解とする数値を個別に設定して、アニーリングにより少なくとも1回、正解とする数値を算出するまでに要する時間」を意味する。
尚、正解とする数値は必ずしも真値である必要はなく(そもそもアニーリングの対象は真値が不明な場合が多い)、制約条件を満たした上での、狙いとする値などを設定できる。上記のナップザック問題であれば、250円以上などと設定すればよい。

TTSは以下の式で定義される。

TTS = \frac{\log(1-p_{d})}{\log(1-p_{s})}×t\\
ここで、\\
p_{d}：要求精度(0.99が標準的)、p_{s}：正解割合、t：1回の平均アニーリング時間

実際にTTSを計算してみる。
まず、Dwaveで計算を行い、計算時間、サンプリング回数、最大値を返す関数を定義する。

import time

def get_result(_target_qubo, _num_reads = 100):
  ###埋込前後の時間計測
  _time_0= time.time()
  sampler_dw = EmbeddingComposite(dw_sampler)
  _time_1 = time.time()

  ###Dwaveでの計算前後の時間計測
  _time_2 = time.time()
  _solution = sampler_dw.sample_qubo(_target_qubo, num_reads=_num_reads)
  _time_3 = time.time()

  ###計算時間の取得
  _time_embedding = (_time_1 - _time_0)
  _time_dwave = (_time_3 - _time_2)

  ###計算回数の取得
  _case = len(_solution.record)

  ###計算結果
  _ans = max([_ans[2] for _ans in view(_solution) if _ans[0]==1 and _ans[1]<=C])#何を正解とするかのif文は問題毎に変更必要

  ###Output
  return  [_time_embedding,_time_dwave,_case,_ans]

###テスト
get_result(Qubo_wos)

<Output>
[8.440017700195312e-05, 0.029740333557128906, 39, 270]

続いて、TTSを計算する関数を定義する。尚、TTSの時間はEmbeddingとDwave計算の双方とする。

import math
def tts (_target_qubo,_n,_tar,_num_reads=100,_pd=0.99):
  ###get resultで得られた結果をリストに代入
  _result = np.array([get_result(_target_qubo,_num_reads) for n in range(_n)])

  ###リストの1列、2列目の埋め込み時間、Dwave時間のそれぞれを合計する
  _time = sum(sum(_result)[0:2])/_n

  ###リストの4列目の結果の内、正解に該当する数を取得し、計算回数で割り、正答率を出す。
  _ps = sum(np.where(np.array(_result).T[3]>=_tar,1,0))/_n

  ###リストの3列目のサンプル回数平均値を取得
  _case = sum(sum(_result)[2:3])/_n

  ###TTSの計算
  _tts = math.log(1-_pd)/math.log(1-_ps)
  _tts = _tts*_time

  ###Output
  return [_time,_ps,_case,_tts,_result]

それではTTSを計算してみる。

###補助スピン無の場合
result_wos = tts(Qubo_wos,20,270,1000)
print(f"補助スピン無の場合_ TTS: {result_wos[3]}, Rate: {result_wos[1]}")

###補助スピン有の場合
result_ws = tts(Qubo_ws,20,270,1000)
print(f"補助スピン有の場合_ TTS: {result_ws[3]}, Rate: {result_ws[1]}")

補助スピン無の場合_ TTS: 0.04862314782014121,Rate: 0.95
補助スピン有の場合_ TTS: 0.9433376911659825,Rate: 0.7

これまで補助スピンを使用した場合、計算精度が下がることが確認できていたが、改めてTTSで比較してみるとやはり計算精度が低下することが確認された。

おわりに

今回は独立型の補助スピンを使用した。更に、補助スピンの使用有無での計算精度の違いをTTSを用いて評価した。計算精度の低下が確認されたが、これは$\lambda$の設定がうまく行ってないものと考えられる。Optunaを用いて最適化トライしているが小規模の計算回数では精度上がらず、大規模な繰り返しが必要となる。
一方で、次回は関連型の補助スピンの使用方法を検討する。具体的にはシフト最適化に於いて、3連勤禁止などの条件に活用することが出来る。

今後の課題 : 関連型の補助スピンの活用
積み残し : 巡回セールスマン問題での分割法適用、時間制限付き巡回セールスマン問題