はじめに
D-Waveがコロナウィルスのための研究のために実機を解放してくれたので少しでも役に立てたらとおもう。今後研究を進めていくと大部分の人が、ハイパーパラメータ調整に苦労するとおもうのでこの記事を公開しました。
アニーリングマシンを使う上で、ハイパーパラメータ調整を自動化するためにOptunaをつかった方法を記述していく。
またはじめに、このソースの７割くらいは、２０１８年度IPA未踏ターゲット事業（アニーリング）に開発したものでご支援いただいたことに感謝しています。
D-Wave
https://cloud.dwavesys.com/leap/signup/covid19/
Optuna
https://preferred.jp/ja/projects/optuna/
IPA 未踏ターゲット
https://www.ipa.go.jp/jinzai/target/index.html
ハイパーパラメータ調整の概要
アニーリングマシンで計算を行う場合、マシン性能を最大限引き出すためには、ハイパーパラメータのチューニングが必要不可欠である。本プロジェクトでは、ハイパーパラメータ・チューニングのためのツール開発も行った。図1に開発したツールのハイパーパラメータ調整の概要を示す。アニーリング マシンで計算後出た解について、
１、 Feasible な解が得られているか検証
２、 Feasibleな解が得られる確率をサンプルから推定
３、 Feasibleな解が得られない確率が小さくなるように、PyQUBO のプレイスフォルダーを利用し、Optuna によりペナルティ項のパラメータを調整
以上３つのステップでハイパーパラメータの調整を行っている。また開発したツールにより、３のペナルティ項のパラメータ調整以外のマシン特有のハイパーパラメータの調整も可能である。
図1: ハイパーパラメータ調整の概要
例としてNクイーン問題で行う
Nクイーン問題の詳細は、昔かいた記事の
「Nクイーン問題をPyQUBOのSAで解いてみる」(https://qiita.com/morimorijap/items/5dce94f7c40b04f34165)
を参照してほしい。
PyQUBOのPlaceholderを利用して、 斜めの部分のパラメータ以下では、gammaについて、Optunaで自動的にパラメータを調整します。
このときに、Feasibleな解が得られない確率が小さくなるように判定をしています。
まずは、必要なものをimportする。
n-queen.ipynb
%matplotlib inline
from pyqubo import Array, Placeholder, solve_qubo, Constraint, Sum
from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite
from dwave_qbsolv import QBSolv
import pandas as pd
import numpy as np
import optuna
N-Queenのマス数の設定
n-queen.ipynb
# N-Queenのマスの数
N = 8
#pyQUBOで0、１の変数をNxN正方マス分つくる
q = Array.create('q', (N*N), 'BINARY')
q_shape = np.reshape(q,(N,N))
コストをPyQUBOで記述
n-queen.ipynb
#row
row_const = 0.0
for row in range(N):
row_const += (sum(i for i in q_shape[row, :]) -1)**2
#colum
col_const = 0.0
for col in range(N):
col_const += (sum(i for i in q_shape[:, col]) -1)**2
#print(col_const)
#print('---')
#斜め
# ＼の方
diag_const = 0.0
xi = 0.0
xi_1 = 0.0
for k in range(-N+2,N-1):
xi += (sum(i for i in np.diag(q_shape,k=k)))
xi_1 += (sum(i for i in np.diag(q_shape,k=k))) -1
diag_const += xi * xi_1
#斜め
# /の方のために入れ替えて, ＼ をする。
diag_const_f = 0.0
xi = 0.0
xi_1 = 0.0
for k in range(-N+2,N-1):
xi += (sum(i for i in np.diag(np.fliplr(q_shape),k=k)))
xi_1 += (sum(i for i in np.diag(np.fliplr(q_shape),k=k))) -1
diag_const_f += xi * xi_1
# Hyper parameter and Hamiltonian is using Placeholder
alpha = Placeholder("alpha")
beta = Placeholder("beta")
gamma = Placeholder("gamma")
H = alpha * col_const + beta * row_const + gamma *(diag_const + diag_const_f)
D-Waveでの計算部分、実機をつかわなくても、Qbslovで計算をすることも可能。
Qbsolvである程度パラメータのターゲットを決めて、D-Waveの実機にかけるなどもありとおもいます。
n-queen.ipynb
def calc_dwave(param_a,param_b,param_c,param_num_reads,param_index_label,param_qv):
# 制約項のペナルティウェイト
# alpha = param_a
# beta = param_b
# gamma = param_c
# num_reads = param_num_reads
# index_label= param_index_label #True, False
# pv = param_qv #True, False
#パラメータ
feed_dict = {'alpha': param_a, 'beta': param_b, 'gamma': param_c}
#QUBOのコンパイル
model = H.compile()
qubo, offset = model.to_qubo(feed_dict=feed_dict, index_label= param_index_label)
# Solve QUBO model by D-Wave
# ---------------------------
endpoint = "https://cloud.dwavesys.com/sapi"
token = "CV19-" #<-- your token
solver_name = "DW_2000Q_5"
sampler = EmbeddingComposite(DWaveSampler(endpoint=endpoint,token=token,solver=solver_name))
if param_qv == True:
response = QBSolv().sample_qubo(qubo, num_reads = param_num_reads)
else:
response = sampler.sample_qubo(qubo, num_reads = param_num_reads)
return response
n-queenのチェックプログラム
n-queen.ipynb
def solution_check(res_pd_p):
# n-queenのチェックプログラム
res_pd_p = res_pd_p
# 行数を確認
for kk in range(len(res_pd_p.index)):
# 行のスライス 0:1, q[0]..q[N*N] 切り取り
res_pd_NN = res_pd_p.iloc[kk:kk+1, 0:N*N]
# q[0]..q[N*N]に並べ替え
res_pd_NNsort = res_pd_NN.loc[:, list(map(lambda x: "q[{}]".format(x), range(0,N*N)))]
#pandas から arrayへ
res_pd_NNsortarray = res_pd_NNsort.values
res_pd_NN_shape = np.reshape(res_pd_NNsortarray ,(N,N))
#print(res_pd_NN_shape)
#ここから、判定 もし 満たしてなければ、 param_x_ng へフラグを立てる
#行- の足し算 の判定
for i in np.sum(res_pd_NN_shape, axis=1):
if i >= 2:
res_pd_p.at[kk,'param_a_ng'] = 1
else:
pass
#列| の足し算
for i in np.sum(res_pd_NN_shape, axis=0):
if i >= 2:
res_pd_p.at[kk,'param_b_ng'] = 1
else:
pass
# 斜め/
for k in range(-N+2,N-1):
i = np.sum(np.diag(res_pd_NN_shape,k=k))
if i >= 2:
res_pd_p.at[kk,'param_c_ng'] = 1
else:
pass
# 斜め\
for k in range(-N+2,N-1):
i = np.sum(np.diag(np.fliplr(res_pd_NN_shape),k=k))
if i >= 2:
res_pd_p.at[kk,'param_c_ng'] = 1
else:
pass
#or で param a, b, cにNGあったらparam_ngにまとめる
res_pd_p['param_ng'] = ((res_pd_p['param_a_ng'] == 1) | (res_pd_p['param_b_ng'] == 1) | (res_pd_p ['param_c_ng'] == 1))
#ture, false に *1すると 0,1になる。
res_pd_p['occurrences_prob'] = res_pd_p['num_occurrences'] * res_pd_p['param_ng'] * 1
#全部でどのくらいfissible non-fisibileがあるかカウント num_occurrences
##res_pd_p["num_occurrences"].sum(axis=0)
#失敗率
o_prob = res_pd_p["occurrences_prob"].sum(axis=0)/res_pd_p["num_occurrences"].sum(axis=0)
return o_prob
計算を実際に投げる部分
n-queen.ipynb
def calc_main(xc):
alpha = 1.0
beta = 1.0
gamma = xc
#
# response = calc_dwave(alpha, beta, gamma ,10,False,False)
# (alpha, beta, gamma ,計算回数=10,pyqubo index label=False, D-wave or QBSolv use ex: True = D-Wave, False=QBSolve)
#
response = calc_dwave(alpha, beta, gamma ,10,False,True)
#結果のPandas化
res_pd = response.to_pandas_dataframe()
#行、列、斜めで、満たさない場合がある場合にフラグと立てる場所の項目を作成
res_pd_p = res_pd.assign(param_a_ng = 0, param_b_ng = 0, param_c_ng = 0)
score_m = solution_check(res_pd_p)
return score_m
Optunaによるスコアリングの部分
n-queen.ipynb
# 失敗率でoptunaでパラメータチューニングしてみる。
def objective(trial):
x = trial.suggest_uniform('x', 0,1)
score=calc_main(x)
print('x: %1.3f, score: %1.3f' % (x, score))
return score
OptunaでPyQUBOのパラメータを調整しながらやる部分
n-queen.ipynb
study = optuna.create_study()
study.optimize(objective, n_trials=10)
print("params_{}".format(study.best_params))
print("value_{}".format(study.best_value))
これによって出た結果は、
[I 2020-04-02 14:13:09,287] Finished trial#0 resulted in value: 0.014925373134328358. Current best value is 0.014925373134328358 with parameters: {'x': 0.3716248996269532}.
x: 0.372, score: 0.015
[I 2020-04-02 14:13:11,342] Finished trial#1 resulted in value: 0.918918918918919. Current best value is 0.014925373134328358 with parameters: {'x': 0.3716248996269532}.
x: 0.056, score: 0.919
[I 2020-04-02 14:13:13,497] Finished trial#2 resulted in value: 1.0. Current best value is 0.014925373134328358 with parameters: {'x': 0.3716248996269532}.
x: 0.022, score: 1.000
[I 2020-04-02 14:13:15,488] Finished trial#3 resulted in value: 0.0. Current best value is 0.0 with parameters: {'x': 0.2101892532802433}.
x: 0.210, score: 0.000
[I 2020-04-02 14:13:17,501] Finished trial#4 resulted in value: 0.0. Current best value is 0.0 with parameters: {'x': 0.2101892532802433}.
x: 0.661, score: 0.000
[I 2020-04-02 14:13:19,568] Finished trial#5 resulted in value: 0.0. Current best value is 0.0 with parameters: {'x': 0.2101892532802433}.
x: 0.325, score: 0.000
[I 2020-04-02 14:13:21,845] Finished trial#6 resulted in value: 0.0. Current best value is 0.0 with parameters: {'x': 0.2101892532802433}.
x: 0.277, score: 0.000
[I 2020-04-02 14:13:23,875] Finished trial#7 resulted in value: 0.0. Current best value is 0.0 with parameters: {'x': 0.2101892532802433}.
x: 0.988, score: 0.000
[I 2020-04-02 14:13:25,931] Finished trial#8 resulted in value: 0.0. Current best value is 0.0 with parameters: {'x': 0.2101892532802433}.
x: 0.185, score: 0.000
[I 2020-04-02 14:13:28,039] Finished trial#9 resulted in value: 0.9629629629629629. Current best value is 0.0 with parameters: {'x': 0.2101892532802433}.
x: 0.032, score: 0.963
params_{'x': 0.2101892532802433}
value_0.0
[-]
のようになります。
ここでは、斜めの項のパラメータを調整するだけでしたが、複数のパラメータ調整、ハードウェアのパラメータをつかった調整なども可能とおもいます。
最後に、GitHubにファイルを公開しておきます。