QAP.07：BQMに制約条件を仕込む(pyqubo版)【量子コンピュータ/アニーリング@Python/D-Wave】
【はじめに】
※再掲
「D-Wave（Ocean SDK）」では、計算させたい問題の種類にあわせて、いい感じに取り扱ってくれる「Model（計算モデルオブジェクト）」を用意している。大きく分けて「３つ」ある。
今回は「制約条件」の中でも「BQMに対しての制約条件の仕込み方(pyqubo版)」について述べる。（※）
※2022年2月現在、数式の通りに制約条件を記述して追加していく「CQM」が提供されているため、わざわざ「BQM」を使って制約条件を仕込む、ということはあまりないかもしれない。
あくまで「制約条件の仕込み方の一つ」程度でおさえておこう。
【BQMに制約条件を仕込む(pyqubo版)】
結論から言うと、、、
「D-Wave」の「BQM」の場合は「Fixstars Amplify」と同じように
『目的関数:ｆの後ろにペナルティ関数:ｇ（＋重み）』をつけた数式（エネルギー式）』を作成すればよい。
■イメージ図（再掲）
以上を踏まえて、次の例題をやってみる（「Fixstars Amplify」と同じ例題）
【例題】x + 2y - 3zに制約条件ある場合（pyqubo版）
次の「目的関数」の値をできるだけ小さくする「(x, y, z)の組み合わせ」を求めよう。(※x, y, zは「0または1」のどちらかの値をとる)
・目的関数：x + 2y - 3z
・制約条件：x + y + z = 2 （※３つのうち２つの変数が１になる）
【１】ライブラリのインストール
#ライブラリのインストール
!pip install dwave-ocean-sdk（※再掲：ネット上を検索すると「dwave-system」というよく似たライブラリも出てくるので困惑するかもしれないが、「dwave-ocean-sdk」の「pip」時に一緒に入るので気にしなくてよい。）
【２】トークン設定
※再掲：
トークンを使ったD-Waveへのアクセスについては以下の通り
上記ドキュメント記載の通り、本来はアクセストークン漏洩防止など、セキュリティの観点から「コンフィグファイル」や「環境変数」に仕込む。
【注意】
今回は簡単なプログラムをちょっと試してみるだけなので、「Samplerオブジェクトに仕込む」というやり方で手抜きをしている。
リリースアプリ、実運用ではSamplerへのトークン埋め込みはしないこと。
■トークン設定
# Sampler埋め込み用トークン（本番では行わないこと）
# 各自登録して取得したアクセストークンを指定する
MY_DWAVE_TOKEN = "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"トークンを使った接続確認として、使用可能なSolver情報を出力してみる。
from dwave.cloud import Client
# トークンを使った接続確認
# 使用可能なSolver情報を出力する
client = Client.from_config(token=MY_DWAVE_TOKEN)
print(client.get_solvers())【実行結果例】
[BQMSolver(id='hybrid_binary_quadratic_model_version2'), DQMSolver(id='hybrid_discrete_quadratic_model_version1'), CQMSolver(id='hybrid_constrained_quadratic_model_version1'), StructuredSolver(id='DW_2000Q_6'), StructuredSolver(id='Advantage_system4.1')]
【３】変数を用意する（pyqubo）
※再掲
「Fixstars Amplify」で「QUBO変数/イジング変数」と呼ばれていたものは、「D-Wave」では次のように呼ばれている。
Spin：「-1」 または「1」をとる （※イジング変数相当）
Binary：「0」または「1」をとる （※QUBO変数相当）
今回の例題は「Binaryオブジェクト」を使って変数を生成する。
また、Binaryオブジェクトには２種類ある。
・pyqubo版
・dimod版
今回は「pyqubo版のBinaryオブジェクト」を使用して変数を生成し「BQMオブジェクト」を作っていく。
※参考：第３回ではdimod版を使用した
■pyqubo.Binaryオブジェクトで変数生成
import pyqubo
# pyqubo.Binaryオブジェクトの生成
q0 = pyqubo.Binary('q0')
q1 = pyqubo.Binary('q1')
q2 = pyqubo.Binary('q2')【４】目的関数:ｆの作成
数式の表記通り記述して目的関数を作成する
■目的関数:ｆの作成
# 目的関数の作成
f = q0 + 2*q1 -3*q2
print(f)【ここまでの実行結果例】
((Binary('q0') + (2.000000 * Binary('q1'))) + (-1.000000 * (3.000000 * Binary('q2'))))
【５】ペナルティ関数:ｇの作成、Constraintオブジェクト
「pyqubo」では「制約条件を管理するオブジェクト」として「Constraintオブジェクト」を用意している
これを使って「制約条件」に対する「ペナルティ関数」を作る。
今回の「制約条件：x + y + z = 2」に対する「ペナルティ関数(数式)」は「(x + y + z - 2)^2」となる。
■ペナルティ関数:ｇを作成する
# 「制約条件」に対応する「ペナルティ関数(数式)」をセットする
g = pyqubo.Constraint((q0+q1+q2-2)**2, "const1")
print(g)【ここまでの実行結果例】
Constraint(((((Binary('q0') + Binary('q1')) + Binary('q2')) + -2.000000) * (((Binary('q0') + Binary('q1')) + Binary('q2')) + -2.000000)), 'const1')
※Note
「BQM」では「制約条件を表現する数式（ペナルティ関数）」を知らないとプログラムを作ることができない。
「Constraintオブジェクト」にはあくまで、「制約条件に対応する数式（ペナルティ関数）」を設定することに注意する。
これに対し、制約条件式の表記をそのまま記述して設定できるのが「CQM」である。
【６】重みを決めて数式（エネルギー式）を完成させる
次に「ペナルティ関数ｇ」に対する「重み」を決めて「最終的な数式（エネルギー式）」を完成させる。
今回は「重み：1.0」にしてみる。
■重みを決めて数式を完成させる
# 最終的な数式、重み:1
h = f + 1.0*g
print(h)【ここまでの実行結果例】
(((Binary('q0') + (2.000000 * Binary('q1'))) + (-1.000000 * (3.000000 * Binary('q2')))) + Constraint(((((Binary('q0') + Binary('q1')) + Binary('q2')) + -2.000000) * (((Binary('q0') + Binary('q1')) + Binary('q2')) + -2.000000)), 'const1'))
▲最終的な数式（エネルギー式）が作られているのもわかる。
【７】pyquboで作成した数式をコンパイルしてBQM化する
「pyqubo」で作成した数式は「compile()」を使ってコンパイルする必要がある。
これにより「QUBO（BQMオブジェクト）」作成のための情報をもつ「Modelオブジェクト」を生成する。
これを「D-Wave」で取り扱えるように「to_bqm()」を使って「BQMオブジェクト」にする
■pyquboで作成した数式をコンパイルし、BQMオブジェクトを生成する
model = h.compile() # コンパイル
bqm = model.to_bqm() # bqmを取得
print(bqm)【ここまでの実行結果例】
BinaryQuadraticModel({'q0': -2.0, 'q1': -1.0, 'q2': -6.0}, {('q1', 'q0'): 2.0, ('q2', 'q0'): 2.0, ('q2', 'q1'): 2.0}, 4.0, 'BINARY')
【８】Sampler経由で量子アニーリング計算する
※再掲：Samplerオブジェクト
「D-Wave（Ocean SDK）」では、「解きたい問題に対するパラメータを設定したModel」を「Samplerオブジェクト」に渡す。これは大きく分けて５つある。
ここでは「Sampler」として「DWaveSampler」を使用してみる。
（※繰り返しになるが、今の段階では細かいことは気にせず「DWaveSampler」を使う時は「EmbeddingComposite」と組み合わせて使う、ということをごま塩程度に覚えておいてくれればいい）
from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite
# DWaveSampler＋EmbeddingCompositeでsamplerを生成する
sampler = EmbeddingComposite(DWaveSampler(token=MY_DWAVE_TOKEN))
# num_reads：アニーリングサイクルは500（QPUが試行する回数、問題を解く回数のこと）
response = sampler.sample(bqm,num_reads=1000)■EmbeddingComposite
■DWaveSampler
【９】計算結果を確認する
「Sampler」経由で「D-Wave上のSolver」が計算した結果は「SampleSetオブジェクト」として返ってくる。
■SampleSetオブジェクトについて
for sample, energy, num_occurrences, chain_break_fraction in list(response.data()):
print(sample, "Energy:", energy, "Occurrences:", num_occurrences)【実行結果例】
{'q0': 0, 'q1': 0, 'q2': 1} Energy: -2.0 Occurrences: 604
{'q0': 1, 'q1': 0, 'q2': 1} Energy: -2.0 Occurrences: 396
▲「制約条件」を満たしているかどうかは別として、「最小値（エネルギー）」となる組み合わせを出力している。（最小値となる組み合わせとしては正解）
【１０】制約条件を満たしているかチェックする
制約条件を満たしているかどうかは、次のように確認できる
コンパイルした「pyquboのModelオブジェクト」に出力された解を渡す
※「decode_sample()」をコールして「DecodedSampleオブジェクト」を取得する「DecodedSampleオブジェクト」の「constraints()」でチェックする
■Model.decode_sample()
■DecodedSample
■pyquboで制約条件の充足具合をチェックする
# 制約条件の充足状況確認
for sample, _ , _ , _ in list(response.data()):
dec = model.decode_sample(sample, vartype='BINARY') # dict型で確認したいデータを渡せばよい
print(dec.constraints(), sample)【実行結果】
{'const1': (False, 1.0)} {'q0': 0, 'q1': 0, 'q2': 1}
{'const1': (True, 0.0)} {'q0': 1, 'q1': 0, 'q2': 1}
▲制約条件：(x + y + z = 2)を満たしていない組み合わせは「False」となる。
このような感じで「D-WaveのBQM」を使った「量子アニーリング計算」について、
「pyqubo」で変数を定義して「目的関数」や「ペナルティ関数」を作成すると、「pyquboの機能」で「制約条件の満たし具合をチェック」することができる。
【１１】全体コード
最後に全体コードを示しておく
【例題】
・x + 2y - 3zを最小にするx, y, zを探す。
ただし
(i)「x + y + z = 2（制約条件）」を満たすこと。
【実行環境】
・D-Wave Leap ＋ Google Colab
■事前準備：ライブラリのインストール
#ライブラリのインストール
!pip install dwave-ocean-sdk■ここからが全体コード
import pyqubo
from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite
from dwave.cloud import Client
# Sampler埋め込み用トークン（本番では行わないこと）
# 各自登録して取得したアクセストークンを指定する
MY_DWAVE_TOKEN = "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
# トークンを使った接続確認
# 使用可能なSolver情報を出力する
client = Client.from_config(token=MY_DWAVE_TOKEN)
print(client.get_solvers())
# pyqubo.Binaryオブジェクトの生成
q0 = pyqubo.Binary('q0')
q1 = pyqubo.Binary('q1')
q2 = pyqubo.Binary('q2')
# 目的関数の作成
f = q0 + 2*q1 -3*q2
print(f)
# 「制約条件」に対応する「ペナルティ関数(数式)」をセットする
g = pyqubo.Constraint((q0+q1+q2-2)**2, "const1")
print(g)
# 最終的な数式、重み:1
h = f + 1.0*g
print(h)
# pyquboで作成した数式をコンパイルしてbqm取得
model = h.compile() # コンパイル
bqm = model.to_bqm() # bqmを取得
print(bqm)
# Sampler経由で量子アニーリング計算する
# DWaveSampler＋EmbeddingCompositeでsamplerを生成する
sampler = EmbeddingComposite(DWaveSampler(token=MY_DWAVE_TOKEN))
# num_reads：アニーリングサイクルは500（QPUが試行する回数、問題を解く回数のこと）
response = sampler.sample(bqm,num_reads=1000)
# 結果の確認
for sample, energy, num_occurrences, chain_break_fraction in list(response.data()):
print(sample, "Energy:", energy, "Occurrences:", num_occurrences)
# 制約条件の満たし具合を確認
for sample, _ , _ , _ in list(response.data()):
dec = model.decode_sample(sample, vartype='BINARY') # dict型で確認したいデータを渡せばよい
print(dec.constraints(), sample)
【実行結果例】
{'q0': 0, 'q1': 0, 'q2': 1} Energy: -2.0 Occurrences: 604
{'q0': 1, 'q1': 0, 'q2': 1} Energy: -2.0 Occurrences: 396
{'const1': (False, 1.0)} {'q0': 0, 'q1': 0, 'q2': 1}
{'const1': (True, 0.0)} {'q0': 1, 'q1': 0, 'q2': 1}
▲重み調整をわざと失敗させ、ペナルティ関数があまり影響しないパターンにした。そのため、制約条件を違反した変数の組み合わせも観測されている。（今回の実行結果例では604回観測された）
この後は実際のところ、「重みを調整して再計算する」「制約条件の満たし具合でフィルタリング」、、等のようなことをして最終的な答えどれにするか決める。
もっと応援したいなと思っていただけた場合、よろしければサポートをおねがいします。いただいたサポートは活動費に使わせていただきます。