-容量$C$のカバンに、夫々値段$p_i$と大きさ$v_i$が定義される$N$個のお菓子を、出来る限り値段が高くなるように、最大量つめる場合、どのような組合せにするか調べる問題を、アニーリングを用いて実施する。
-更に、今回はハイパーパラメタをOptunaを用いて探索する。
前回: 第1回: アニーリング基礎1. 荷物の分担
Keyword ：ナップザック問題、Optuna、ハイパーパラメーターの最適化
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)、使用するマシンの指定を行う。
token = "DEV-73ac79954ef2135ada52dXXXXXXXXXXXXXXXXXX"
endpoint = "https://cloud.dwavesys.com/sapi/"
solver = "Advantage_system4.1"
更に、Dwaveマシンを指定するDwaveSamplerと、Qbitへの割り当て行うEmbeddingCompositeをインポート。続いて、DWaveSamplerを用いてDwaveマシンを指定する。尚、Dwaveの動かしかた概要は、Qiitaの記事を参考/ https://qiita.com/y_daito/items/aad74517927e393564ee
from dwave.system import DWaveSampler, EmbeddingComposite
dw_sampler = DWaveSampler(solver=solver, token=token, endpoint = endpoint)
1.3 各種モジュールのインポート
openjijからSQASampler、pyquboからArray、その他としてnumpy、pandas、matplotlib.pyplotをインポートする。
from openjij import SQASampler
from pyqubo import Array, Placeholder
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
これで下準備は完了。
2. 定式化
まず始めに、お菓子を入れる入れないを$x_i$が0であれば入れない、1であれば入れると定義する。
今回の最適化は、値段の最大化である。単純に$\vec{x}$と$\vec{p}$の合計が最小化されればよい。従い、エネルギ―式を以下とする。
E(x) =-\sum_{i=0}^Nx_i p_i \\
続いて、コスト関数を定義する。
選択した重量が容量$C$から異なれば、ハミルトニアンが増加するように設定する。
C_1(x) =(C-\sum_{i=0}^Nx_i w_i)^2 \\
C_1(x) =C^2-2C\sum_{i=0}^Nx_i w_i+\sum_{i=0}^N\sum_{j=0}^Nx_i w_ix_j w_j
最後に、エネルギー式とコスト関数を足し合わせて、ハミルトニアンが完成。
H(x) = E(x) + α_1C_1(x) \\
H(x) = -\sum_{i=0}^Nx_i p_i + α_1(C-\sum_{i=0}^Nx_i w_i)^2 \\
H(x) = -\sum_{i=0}^Nx_i p_i + α_1(C^2-2C\sum_{i=0}^Nx_i w_i+\sum_{i=0}^N\sum_{j=0}^Nx_i w_ix_jw_j) \\
定式化は以上。続いて実装する。
3. 実装
3.1 Qubo行列の作成からDwaveでの計算
今回のお菓子の数$N$を8個、カバンの容量$C$を300として、夫々の値段$p_i$、夫々の重さ$w_i$を以下の通り設定する。
N = 8
C = 300
p = [80,110,100,90,60,75,40,25]
w = [100,160,120,100,80,70,40,30]
まずPyqubo用いて$\vec{x}$を作成する。
from pyqubo import Array
x = Array.create(name="x", shape =(N), vartype="BINARY")
また、今回はPyquboを用いてコスト関数の影響度合い$a_1$を設定する。
from pyqubo import Placeholder
a1 = Placeholder("a1")
続いて、Pyqubo用いて下式のハミルトニアンを表現する。
H(x) = -\sum_{i=0}^Nx_i p_i + α_1(C-\sum_{i=0}^Nx_i w_i)^2 \\
#第1項
E = - np.dot(x,p)
#第2項
C1 = a1*(C - np.dot(x,w))**2
#ハミルトニアン
H = E + C1
続いて、Qubo行列を取り出す。
model = H.compile()
Qubo, offset = model.to_qubo(feed_dict={"a1":5})#feed_dictでa1の数値設定を実施
DwaveにQubo行列を投げる。
sampler_dw = EmbeddingComposite(dw_sampler)
solution = sampler_dw.sample_qubo(Qubo, num_reads=100)
print(solution.record[:5])
<Output>
[([1, 0, 0, 1, 0, 1, 0, 1], -450270., 2, 0.)
([0, 0, 1, 0, 1, 1, 0, 1], -450260., 3, 0.)
([0, 1, 0, 0, 0, 1, 1, 1], -450250., 2, 0.)
([1, 1, 0, 0, 0, 0, 1, 0], -450230., 1, 0.)
([0, 0, 1, 0, 1, 1, 1, 0], -449775., 4, 0.)]
結果の可視化を関数化する。
def view(solution, n=5):
for i in range(n):
answer = solution.record[i][0]
price, weight = np.dot(answer,p), np.dot(answer,w)
print(f"{i+1}番目の結果: 重量合計は{weight}。値段合計は{price}")
view(solution,10)
<Output>
1番目の結果: 重量合計は300。値段合計は270
2番目の結果: 重量合計は300。値段合計は260
3番目の結果: 重量合計は300。値段合計は250
4番目の結果: 重量合計は300。値段合計は230
5番目の結果: 重量合計は310。値段合計は275
6番目の結果: 重量合計は290。値段合計は265
7番目の結果: 重量合計は290。値段合計は265
8番目の結果: 重量合計は290。値段合計は255
9番目の結果: 重量合計は310。値段合計は255
10番目の結果: 重量合計は290。値段合計は255
以上。
3.2 Optunaを用いたハイパーパラメタの最適化
ハイパーパラメタの設定とは、何だかの目的値を最大化する最適化の一種。つまり、目的値を設定しないといけない。
そこで、今回は得られた結果の内、値段合計が最大であり、重量合計が300以下となっている組合せ数を計算数(num_readsの数)で割ったものを正答率と定義して、この正答率を最大化する$a1$をOptunaで探索する。
はじめに、以下の通り、正答率を算出する関数を定義する。
def rate(solution):
### Dwaveの計算結果を重量合計と値段合計のリストに返還する
result = []
_temp = []
for i in range(len(solution.record)):
_w, _p = np.dot(solution.record[i][0],w), np.dot(solution.record[i][0],p)
### 制約条件を満たさない組合せ結果は[0, 0]とする
if _w > C:
_temp = [0,0]
else:
_temp = [_w, _p]
result.append(_temp)
### 得られた結果を転置して、値段の最大値を取得する。
price_list = np.array(result).T[1]
pmax = price_list.max()
### pmaxを含む結果数をカウントする
count_true = np.count_nonzero(price_list == pmax)
###正答率
rate = count_true / len (solution.record) * 100
return(rate)
続いて、アニーリングによる計算を実行する関数を定義する。関数ではDwave以外に、Jijシュミレーターも選択できるようにする。
def a_cal(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を用いてハイパーパラメタの最適化を実施する。手順は以下。
- ハイパーパラメタが取り得る範囲を指定し、ハイパーパラメタを変数とする。
- 変数としたハイパーパラメタを使用し、ハミルトニアンをコンパイルしてQubo行列を取り出す。
- 取り出したQubo行列を作成したa_cal関数を用いてアニーリング計算を実施する。
- 計算結果を作成したrate関数に渡して正答率を計算する。
- 上記の一連の流れを関数として定義する。
- 上記で定義した関数をoptunaに渡して最適化する。
Optunaに就いては、Qiita: Optuna入門を参考にした。
コーディングは次の通り。
###ライブラリの呼び出し
import optuna
###関数の定義
def objective(trial):
###ハイパーパラメタの設定
_a1 = trial.suggest_uniform('_a1', 1, 1000)
###コンパイルとQuboの取り出し
model = H.compile()
Qubo, offset = model.to_qubo(feed_dict={"a1":_a1})
###a_calでアニーリング計算実施し、それをrateで正答率にする。
scoure = rate(a_cal("Jij"))*-1
return scoure
### Optunaでの最適化
study = optuna.create_study()
study.optimize(objective, n_trials=30)
### Optunaでの最適解を表示
best_a1 = [i for i in study.best_params.values()][0]#Optunaで得られたハイパーパラメタは辞書型で整理されるため、辞書から数値を取りだしている。
best_value = study.best_value
print(f"最適なa1:{best_a1}、その時の正答率は:{best_value}")
<Output>
[I 2022-11-20 01:05:50,320] A new study created in memory with name: no-name-c04904XXXXXX/ipykernel_launcher.py:8: FutureWarning: suggest_uniform has been deprecated in v3.0.0. This feature will be removed in v6.0.0. See https://github.com/optuna/optuna/releases/tag/v3.0.0. Use :func:`~optuna.trial.Trial.suggest_float` instead.
[I 2022-11-20 01:05:51,253] Trial 0 finished with value: -5.0 and parameters: {'_a1': 149.12878823094317}. Best is trial 0 with value: -5.0.
[I 2022-11-20 01:05:52,133] Trial 1 finished with value: -6.0 and parameters: {'_a1': 978.200186229481}. Best is trial 1 with value: -6.0.
[I 2022-11-20 01:05:52,579] Trial 2 finished with value: -4.0 and parameters: {'_a1': 853.832892239228}. Best is trial 1 with value: -6.0.
[I 2022-11-20 01:05:52,963] Trial 3 finished with value: -6.0 and parameters: {'_a1': 611.5789182373358}. Best is trial 1 with value: -6.0.
[I 2022-11-20 01:05:53,358] Trial 4 finished with value: -3.0 and parameters: {'_a1': 947.7254202512732}. Best is trial 1 with value: -6.0.
[I 2022-11-20 01:05:53,780] Trial 5 finished with value: -1.0 and parameters: {'_a1': 66.21224854613251}. Best is trial 1 with value: -6.0.
[I 2022-11-20 01:05:54,186] Trial 6 finished with value: -3.0 and parameters: {'_a1': 227.86168781214587}. Best is trial 1 with value: -6.0.
[I 2022-11-20 01:05:54,579] Trial 7 finished with value: -3.0 and parameters: {'_a1': 19.03203831820299}. Best is trial 1 with value: -6.0.
[I 2022-11-20 01:05:54,960] Trial 8 finished with value: -4.0 and parameters: {'_a1': 53.749831319506065}. Best is trial 1 with value: -6.0.
[I 2022-11-20 01:05:55,355] Trial 9 finished with value: -5.0 and parameters: {'_a1': 246.13097126050434}. Best is trial 1 with value: -6.0.
[I 2022-11-20 01:05:55,735] Trial 10 finished with value: -3.0 and parameters: {'_a1': 695.1231998157043}. Best is trial 1 with value: -6.0.
[I 2022-11-20 01:05:56,124] Trial 11 finished with value: -7.000000000000001 and parameters: {'_a1': 510.2109774419712}. Best is trial 11 with value: -7.000000000000001.
[I 2022-11-20 01:05:56,530] Trial 12 finished with value: -2.0 and parameters: {'_a1': 420.5500673231111}. Best is trial 11 with value: -7.000000000000001.
[I 2022-11-20 01:05:56,930] Trial 13 finished with value: -10.0 and parameters: {'_a1': 455.067830482922}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:05:57,334] Trial 14 finished with value: -5.0 and parameters: {'_a1': 425.0029813430113}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:05:57,711] Trial 15 finished with value: -5.0 and parameters: {'_a1': 575.743786226125}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:05:58,099] Trial 16 finished with value: -2.0 and parameters: {'_a1': 369.1900256115749}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:05:58,493] Trial 17 finished with value: -3.0 and parameters: {'_a1': 731.5965806629351}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:05:58,880] Trial 18 finished with value: -1.0 and parameters: {'_a1': 536.589333972217}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:05:59,291] Trial 19 finished with value: -3.0 and parameters: {'_a1': 334.77138022203997}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:05:59,668] Trial 20 finished with value: -4.0 and parameters: {'_a1': 719.7084086904817}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:06:00,052] Trial 21 finished with value: -1.0 and parameters: {'_a1': 990.4407637607262}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:06:00,455] Trial 22 finished with value: -4.0 and parameters: {'_a1': 837.3201048742774}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:06:00,836] Trial 23 finished with value: -4.0 and parameters: {'_a1': 498.3347919942139}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:06:01,219] Trial 24 finished with value: -2.0 and parameters: {'_a1': 634.3885305472178}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:06:01,617] Trial 25 finished with value: -2.0 and parameters: {'_a1': 858.3007458023199}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:06:02,003] Trial 26 finished with value: -3.0 and parameters: {'_a1': 303.96005052654306}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:06:02,428] Trial 27 finished with value: -3.0 and parameters: {'_a1': 469.59039901120616}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:06:02,807] Trial 28 finished with value: -4.0 and parameters: {'_a1': 774.1313095739922}. Best is trial 13 with value: -10.0.
[I 2022-11-20 01:06:03,201] Trial 29 finished with value: -1.0 and parameters: {'_a1': 133.1334195978052}. Best is trial 13 with value: -10.0.
最適なa1:455.067830482922、その時の正答率は:-10.0
得られた最適化なパラメタを用いてハミルトニアンを再度コンパイルし、Dwaveを用いて計算してみる。
model = H.compile()
Qubo, offset = model.to_qubo(feed_dict={"a1":best_a1})
solution = a_cal("Dwave")
print(f"rate:{rate(solution)}")
view(solution,10)
<Output>
rate:3.389830508474576
1番目の結果: 重量合計は300。値段合計は260
2番目の結果: 重量合計は300。値段合計は250
3番目の結果: 重量合計は300。値段合計は240
4番目の結果: 重量合計は300。値段合計は230
5番目の結果: 重量合計は310。値段合計は285
6番目の結果: 重量合計は310。値段合計は275
7番目の結果: 重量合計は290。値段合計は265
8番目の結果: 重量合計は290。値段合計は265
9番目の結果: 重量合計は290。値段合計は255
10番目の結果: 重量合計は290。値段合計は255
改めて、ハイパーパラメタを10などと設定して、結果を見比べてみる。
model = H.compile()
Qubo, offset = model.to_qubo(feed_dict={"a1":5})
solution = a_cal("Dwave")
print(f"rate:{rate(solution)}")
view(solution,10)
<Output>
rate:1.7543859649122806
1番目の結果: 重量合計は300。値段合計は270
2番目の結果: 重量合計は300。値段合計は250
3番目の結果: 重量合計は300。値段合計は240
4番目の結果: 重量合計は300。値段合計は230
5番目の結果: 重量合計は310。値段合計は285
6番目の結果: 重量合計は290。値段合計は265
7番目の結果: 重量合計は290。値段合計は265
8番目の結果: 重量合計は290。値段合計は255
9番目の結果: 重量合計は290。値段合計は255
10番目の結果: 重量合計は290。値段合計は255
複数回繰り返してみても、あまり正答率は変わらない。ただし、ハイパーパラメタを理屈でもって説明するには丁度良さそうである。
おわりに
今回のナップザック問題のコスト関数は、$C_1(x) =(C-\sum_{i=0}^Nx_i w_i)^2$となっている。できる限り$\sum_{i=0}^Nx_i w_i$が$C$に近づくわけだが、コストの観点からは$C+10$になっても$C-10$になっても、コスト関数は2次式なので、同じコストになる。
例えば、容量$C=300$として、条件を満たす$\sum_{i=0}^Nx_i w_i=290$も、条件を満たしていない$\sum_{i=0}^Nx_i w_i=310$も同じコストという扱いになる。従い、$\sum_{i=0}^Nx_i w_i≦300$場合は、$C_1(x)=0$とするような条件分岐を追加する必要がある。
アニーリングでは補助スピン$y_i$を用いる方法がある。次回検討する。
他方、今回Optunaを使用した際、任意に設定した正答率を用いて、これを最大化するハイパーパラメタを探索した。一般的には、アニーリングで正解を見つけ出すまでに要する時間TTSを用いて探索することが多い。これも次回検討する。
今後の課題 : 補助スピン$y_i$の活用, TTSを用いた探索
積み残し : ハミルトニアンが2次式以上になる場合の対応