見出し画像

QAP.11:「DQM:Discrete Quadratic Models」で等式/不等式制約【量子コンピュータ/アニーリング@Python/D-Wave】

【はじめに】

※再掲
DQM」では「制約条件に相当する数式」として「linearパターン」と「quadraticパターン」がある。

・linear:
 各変数を「足し算する形」で制約条件を表現する式の場合。

・quadratic:
 変数同士を「掛け算する形」で制約条件を表現する式の場合。

▲DQMの数式イメージと制約条件に相当する数式の構成

前回は、「制約条件となる数式」を考えておき、それを「quadraticパターン」側に1つ1つ数式として追加していた。

「制約条件」は、「等式制約」の他にも「不等式制約」「論理制約」…など様々あり、『与えられた条件をどのように「数式として表現するのか」は難しい』

【DQMの等式・不等式制約条件設定】

これに対し、「DQM」では
・等式制約の設定用 :add_linear_equality_constraint()
不等式制約の設定用:add_linear_inequality_constraint()
用意している。

今回は「add_linear_equality_constraint()」を使って簡単な例題を解いてみる。

【例題】:合計6になる数の組み合わせ(複数回の数字選択OK)

問題:
(1,2,3)の中から数字を3回取り出して合計6になるような組み合わせを探そう。※同じ数字を複数回ピックアップしてもOK

【考え方】

以下の「等式制約」をみたす「x_0, x_1, x_2」の組み合わせを探せばよい

x_0 + x_1 + x_2 = 6
※ただし「x_0, x_1, x_2」はそれぞれ「1, 2, 3」のいずれかの値をとる

【1】ライブラリのインストール

#ライブラリのインストール
!pip install dwave-ocean-sdk

(※再掲:ネット上を検索すると「dwave-system」というよく似たライブラリも出てくるので困惑するかもしれないが、「dwave-ocean-sdk」の「pip」時に一緒に入るので気にしなくてよい。)

【2】トークン設定


※再掲:
トークンを使った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')]

【3】DQMの生成と変数用意

( i )まずは「DQMオブジェクト」を作成する。

import dimod


# DQMオブジェクトの作成
dqm = dimod.DiscreteQuadraticModel()

※DQMオブジェクトの使い方詳細は以下参照


( ii )次に変数を用意する。

DQMを使う場合は「add_variable()」で「変数名(ラベル)」と「必要なバイナリ変数の数」を設定する。

※「add_variable()」

■変数の用意

# 定義する変数につける係数(今回は全部0)
linear_coef =[0, 0, 0]

# 変数定義(1つ目)
x_0 = dqm.add_variable(len(linear_coef), label="x_0")
#dqm.set_linear(x_0, linear_coef) # 今回は係数0なので不要
print(dqm.get_linear(x_0))

# 変数定義(2つ目)
x_1 = dqm.add_variable(len(linear_coef), label="x_1")
#dqm.set_linear(x_1, linear_coef) # 今回は係数0なので不要
print(dqm.get_linear(x_1))

# 変数定義(3つ目)
x_2 = dqm.add_variable(len(linear_coef), label="x_2")
#dqm.set_linear(x_2, linear_coef) # 今回は係数0なので不要
print(dqm.get_linear(x_2))

【ここまでの実行結果】
array([0., 0., 0.])
array([0., 0., 0.])
array([0., 0., 0.])

※注意点
「各変数は1,2,3のいずれかの離散値をとる」ので、生成する変数に対して係数を設定したくなるが、今回の問題では、目的関数相当のものはいらない。
最初に生成する変数は係数を0とするのを基本としていい。

というのも、「生成した変数」も「エネルギー式(数式)に組み込まれる」ので、係数を与えてしまうと目的関数が設定されたような感じなってしまう。

ようは計算結果がわかりやすくなるように生成するDQM変数の係数は0にしている、ということ。

【4】DQMに制約条件を仕込む

DQMの制約条件は基本的に「制約条件に対応する数式(ペナルティ関数)」を考えることになる。

※再掲

ここで「add_linear_equality_constraint()」を使うことで「等式制約に対応する数式」がいい感じに作られる。

簡単にまとめると、

「変数と係数のペア」のリストで渡すと、全部足し算した形をつくる
②さらに「定数C」を加算して「=0」となるような等式制約式をつくる

というもの。
例えば今回の場合は、

$$
(x_0[0] + 2x_0[1] + 3x_0[2] )+(x_1[0] + 2x_1[1] + 3x_1[2] )+(x_2[0] + 2x_2[1] + 3x_2[2] )- 6= 0
$$

という式を作らせることになる。

# termsの形式は(variable, case, bias)のtuple型
# biasが係数に該当する
answer_discrete = [1 , 2 , 3]
terms =[
        (x_0,0, answer_discrete[0]),
        (x_0,1, answer_discrete[1]),
        (x_0,2, answer_discrete[2]),
        (x_1,0, answer_discrete[0]),
        (x_1,1, answer_discrete[1]),
        (x_1,2, answer_discrete[2]),
        (x_2,0, answer_discrete[0]),
        (x_2,1, answer_discrete[1]),
        (x_2,2, answer_discrete[2]),
]



# 等式制約を設定する
dqm.add_linear_equality_constraint(terms,lagrange_multiplier=1.0, constant=-6)

イメージ図

これで等式制約条件が設定されている。内部的には数式が展開されており、linear部分、quadratic部分に加算されている。

■現在の数式のlinear部分の確認

# 数式のlinear部分の出力確認
print(dqm.get_linear(x_0))
print(dqm.get_linear(x_1))
print(dqm.get_linear(x_2))

【ここまでの実行結果】
[-11. -20. -27.]
[-11. -20. -27.]
[-11. -20. -27.]

▲もともと[0,0,0]だったものに加算されている。

■現在の数式のquadratic部分の確認

# 数式のquadratic部分の出力確認
print( dqm.get_quadratic(x_0, x_1) )
print( dqm.get_quadratic(x_0, x_2) )
print( dqm.get_quadratic(x_1, x_2) )

【ここまでの実行結果】
{(0, 0): 2.0,
(0, 1): 4.0,
(0, 2): 6.0,
(1, 0): 4.0,
(1, 1): 8.0,
(1, 2): 12.0,
(2, 0): 6.0,
(2, 1): 12.0,
(2, 2): 18.0}
… (略) …

あとはSampler経由で量子アニーリング計算をすればよい。

【5】Sampler経由で量子アニーリング計算する

※再掲:Samplerオブジェクト
「D-Wave(Ocean SDK)」では、「解きたい問題に対するパラメータを設定したModel」を「Samplerオブジェクト」に渡す。これは大きく分けて5つある。

「DQM」に対応する「Sampler」は「LeapHybridDQMSampler」である。

LeapHybridDQMSampler

これを使って量子アニーリング計算をする。

from dwave.system import LeapHybridDQMSampler

# トークンを仕込んだSamplerにmodelを投げ込んで計算する
sampler = LeapHybridDQMSampler(token=MY_DWAVE_TOKEN)
sampleset = sampler.sample_dqm(dqm, label='my_total_six')

【6】結果を確認する

「Sampler」経由で「D-Wave上のSolver」が計算した結果は「SampleSetオブジェクト」として返ってくる。

■SampleSetオブジェクトについて

※今回は変数が少なく単純な組み合わせだが、変数が多いのに制約条件がゆるいと膨大な量になってしまうのでDQMでは特に注意しよう。

■求められた組み合わせの数を確認

print( len(list(sampleset.data())) )

【実行結果例】
19

■結果をざっと確認(list化して整形)

print( list(sampleset.data()) )

【実行結果例】
[Sample(sample={'x_0': 0, 'x_1': 2, 'x_2': 1}, energy=0.0, num_occurrences=1), Sample(sample={'x_0': 2, 'x_1': 0, 'x_2': 1}, energy=0.0, num_occurrences=1), Sample(sample={'x_0': 1, 'x_1': 2, 'x_2': 0}, energy=0.0, num_occurrences=1),
… …
Sample(sample={'x_0': 1, 'x_1': 0, 'x_2': 2}, energy=0.0, num_occurrences=1), Sample(sample={'x_0': 1, 'x_1': 0, 'x_2': 2}, energy=0.0, num_occurrences=1), Sample(sample={'x_0': 0, 'x_1': 1, 'x_2': 2}, energy=0.0, num_occurrences=1)]

▲見た目上同じ組み合わせでも、内部的につかっている変数の値に違いがある

エネルギー値が最小のものを取り出し、numpy.unique()で重複削除してみる

import numpy as np

data_arr = list(sampleset.lowest().record.sample)
unique_arr = np.unique(data_arr, axis=0)

print(unique_arr)

【実行結果例】
array([
[0, 1, 2],
[0, 2, 1],
[1, 0, 2],
[1, 1, 1],
[1, 2, 0],
[2, 0, 1],
[2, 1, 0]], dtype=int32)

▲各変数で「1」をたてる位置(インデックス番号)を取得できている

各変数で1を立てるインデックスがわかったので、実際の数字で出力確認してみる。

■インデックスから実際の数字取り出して出力確認する

for data in unique_arr:
  print(f"{answer_discrete[ data[0] ]}, {answer_discrete[ data[1] ]}, {answer_discrete[ data[2] ]},")

【実行結果例】
1, 2, 3,
1, 3, 2,
2, 1, 3,
2, 2, 2,
2, 3, 1,
3, 1, 2,
3, 2, 1,

▲合計6になるパターンが取り出せている。

今回はたまたま7つの組み合わせが出てきたが、この結果より少ない組み合わせの数しか求められない時もある。

答えが間違っている、というわけではない。
量子アニーリング計算では複数回の試行で観測された組み合わせをあげているだけである。

※関連

【7】全体コード

最後に全体コードを示しておく。

問題:
(1,2,3)の中から数字を3回取り出して合計6になるような組み合わせを探そう。※同じ数字を複数回ピックアップしてもOK

【実行環境】
・D-Wave Leap + Google Colab

■事前準備:ライブラリのインストール

#ライブラリのインストール
!pip install dwave-ocean-sdk

■ここからが全体コード

from dwave.cloud import Client
from dwave.system import LeapHybridDQMSampler

import dimod
import numpy as np


# Sampler埋め込み用トークン(本番では行わないこと)
# 各自登録して取得したアクセストークンを指定する
MY_DWAVE_TOKEN = "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"


# トークンを使った接続確認
# 使用可能なSolver情報を出力する
client = Client.from_config(token=MY_DWAVE_TOKEN)
print(client.get_solvers())


# DQMオブジェクトの作成
dqm = dimod.DiscreteQuadraticModel()


# 定義する変数につける係数(今回は全部0)
linear_coef =[0, 0, 0]

# 変数定義(1つ目)
x_0 = dqm.add_variable(len(linear_coef), label="x_0")
#dqm.set_linear(x_0, linear_coef) # 今回は係数0なので不要
print(dqm.get_linear(x_0))

# 変数定義(2つ目)
x_1 = dqm.add_variable(len(linear_coef), label="x_1")
#dqm.set_linear(x_1, linear_coef) # 今回は係数0なので不要
print(dqm.get_linear(x_1))

# 変数定義(3つ目)
x_2 = dqm.add_variable(len(linear_coef), label="x_2")
#dqm.set_linear(x_2, linear_coef) # 今回は係数0なので不要
print(dqm.get_linear(x_2))


# termsの形式は(variable, case, bias)のtuple型
# biasが係数に該当する
answer_discrete = [1 , 2 , 3]
terms =[
        (x_0,0, answer_discrete[0]),
        (x_0,1, answer_discrete[1]),
        (x_0,2, answer_discrete[2]),
        (x_1,0, answer_discrete[0]),
        (x_1,1, answer_discrete[1]),
        (x_1,2, answer_discrete[2]),
        (x_2,0, answer_discrete[0]),
        (x_2,1, answer_discrete[1]),
        (x_2,2, answer_discrete[2]),
]



# 等式制約を設定する
dqm.add_linear_equality_constraint(terms,lagrange_multiplier=1.0, constant=-6)


# 数式のlinear部分の出力確認
print(dqm.get_linear(x_0))
print(dqm.get_linear(x_1))
print(dqm.get_linear(x_2))


# 数式のquadratic部分の出力確認
print( dqm.get_quadratic(x_0, x_1) )
print( dqm.get_quadratic(x_0, x_2) )
print( dqm.get_quadratic(x_1, x_2) )


# トークンを仕込んだSamplerにmodelを投げ込んで計算する
sampler = LeapHybridDQMSampler(token=MY_DWAVE_TOKEN)
sampleset = sampler.sample_dqm(dqm, label='my_total_six')


# 結果の確認
print( len(list(sampleset.data())) )
print( list(sampleset.data()) )

# 重複を削除して結果を整理する
data_arr = list(sampleset.lowest().record.sample)
unique_arr = np.unique(data_arr, axis=0)

print(unique_arr)

for data in unique_arr:
  print(f"{answer_discrete[ data[0] ]}, {answer_discrete[ data[1] ]}, {answer_discrete[ data[2] ]},")

【実行結果例】※抜粋
[[0 1 2]
[0 2 1]
[1 0 2]
[1 1 1]
[1 2 0]
[2 0 1]]

1, 2, 3,
1, 3, 2,
2, 1, 3,
2, 2, 2,
2, 3, 1,
3, 1, 2,

▲改めて実行した今回は6つの組み合わせが出てきた。




もっと応援したいなと思っていただけた場合、よろしければサポートをおねがいします。いただいたサポートは活動費に使わせていただきます。