Go to Qiita Advent Calendar 2024 Top
LoginSignup
3
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@opt_d(D D)

第8回: アニーリング基礎7. PyquboのDecode活用

Last updated at Posted at 2022-12-10

これまで色々なアニーリング問題を整理してきたが、PyquboのDecodeクラスを有効活用できていなかったので、今回整理してみる。

前回: 第7回: アニーリング基礎6.補助スピン活用基礎2（シフト最適化問題）
参考: 公式チュートリル
Keyword ：Pyqubo, Decode

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()
###基本的なライブラリのインポート
import numpy as np
import pandas as pd
import seaborn as sns

2. 問題設定、定式化、最適化計算

2.1 問題

複数の制約式があり、エネルギー式が明確な以下の問題を例にとって、PyquboのDecodeクラスを整理する。
image.png
スコア表の値$S$、スコア表の行数$R$、列数$C$を以下のように定義する。
また、選択箇所を示すスピン$x_{rc}$も設定する。$r行c列$を選択する時、$x_{rc}=1$となる。

S = [[10,20,15],[30,40,35],[50,60,55]]
R, C = len(S),len(S)
from pyqubo import Array
x = Array.create(name="x", shape=(R,C), vartype="BINARY")

2.2 定式化

エネルギ-式$E(x)$は選択したスコアの合計値として、制約1、2併せて以下のように定式化できる。

E(x)=-\sum_{r=1}^R\sum_{c=1}^CS_{rc}x_{rc}\\
C_{1}(x) = \sum_{c=1}^C(\sum_{r=1}^Rx_{rc}-1)^2\\
C_{2}(x) = \sum_{r=1}^R(\sum_{c=1}^Cx_{rc}-1)^2\\
H(x)=\lambda_{e1}E(x)+\lambda_{c1}C_{1}(x)+\lambda_{c2}C_{2}(x)

### 第1項の定式化
E = 0
for r in range(R):
  for c in range(C):
    E = E - S[r][c]*x[r][c]

### 第2項の定式化
C1 = 0
for c in range(C):
  _c = 0
  for r in range(R):
    _c = _c + x[r][c]
  C1 = C1 + (_c-1)**2

### 第3項の定式化
C2 = 0
for r in range(R):
  _c = 0
  for c in range(C):
    _c = _c + x[r][c] 
  C2 = C2 + (_c-1)**2

### lambdaの設定
from pyqubo import Placeholder
lambda_e1,lambda_c1,lambda_c2 = Placeholder("lambda_e1"),Placeholder("lambda_c1"),Placeholder("lambda_c2")

PyquboのDecodeクラスは、Dwaveで得られた結果を見やすく整理してくれる。具体的には、.sampleスピン毎の結果、.energyエネルギー値、.constraints()/.subh制約成立有無、を取得することが出来る。ただし、事前準備として、ハミルトニアンを作成する時に、エネルギー式、制約式を区別する必要がある。Pyquboの制約式をラベリングするConstraint関数により区別が可能。コーディングは以下の通り。Constraint(制約式,"制約式名")でラベリングされる。尚、ラベリングされなかったものがエネルギー式として認識される。

### 制約式をラベリングするConstraintクラスのインポート
from pyqubo import Constraint

### ハミルトニアンの作成
H = lambda_e1*E + lambda_c1*Constraint(C1,"C1") + lambda_c2*Constraint(C2,"C2")

2.3 最適化計算

Qubo行列の取出しから、求解まで一挙に進める。

###Qubo行列の取り出し
model = H.compile()
feed_dict = {"lambda_e1":1,"lambda_c1":100,"lambda_c2":100}
Qubo,offset = model.to_qubo(feed_dict=feed_dict)
###Dwaveでの計算
sampler_dw = EmbeddingComposite(dw_sampler)
result = sampler_dw.sample_qubo(Qubo, num_reads=10)

result.record
----------------------------
<Output>
rec.array([([1, 0, 0, 0, 0, 1, 0, 1, 0], -705., 1, 0.),
           ([0, 1, 0, 1, 0, 0, 0, 0, 1], -705., 1, 0.),
           ([0, 0, 1, 0, 1, 0, 1, 0, 0], -705., 6, 0.),
           ([0, 1, 0, 0, 0, 1, 1, 0, 0], -705., 2, 0.)],
          dtype=[('sample', 'i1', (9,)), ('energy', '<f8'), ('num_occurrences', '<i8'), ('chain_break_fraction', '<f8')])

3. PyquboのDecodeの検証

3.1 Decodeの使い方

decode_sampleset関数にDwaveでの計算結果と、$\lambda$を渡すことでdecodeされた結果を得ることができる。

decoded_result = model.decode_sampleset(sampleset=result,feed_dict=feed_dict)
decoded_result
----------------------------
<Output>
[DecodedSolution({x[0][1]:0, x[0][0]:0, x[1][0]:0, x[2][1]:0, x[0][2]:1, x[2][2]:0, x[1][2]:0, x[2][0]:1, x[1][1]:1}, energy=-105.000000),
 DecodedSolution({x[0][1]:1, x[0][0]:0, x[1][0]:0, x[2][1]:0, x[0][2]:0, x[2][2]:0, x[1][2]:1, x[2][0]:1, x[1][1]:0}, energy=-105.000000),
 DecodedSolution({x[0][1]:1, x[0][0]:0, x[1][0]:1, x[2][1]:0, x[0][2]:0, x[2][2]:1, x[1][2]:0, x[2][0]:0, x[1][1]:0}, energy=-105.000000),
 DecodedSolution({x[0][1]:0, x[0][0]:1, x[1][0]:0, x[2][1]:1, x[0][2]:0, x[2][2]:0, x[1][2]:1, x[2][0]:0, x[1][1]:0}, energy=-105.000000)]

試しに1回目の結果を、decode有り無しで比較すると以下の通り。decodeすることにより、スピンの座標と値の双方が示されるとともに、ハミルトニアン合計値ではなく、エネルギー値のみが返される。

print(result.record[0])
print(decoded_result[0])
----------------------------
<Output>
([0, 0, 1, 0, 1, 0, 1, 0, 0], -705., 6, 0.)
DecodedSolution({x[0][1]:0, x[0][0]:0, x[1][0]:0, x[2][1]:0, x[0][2]:1, x[2][2]:0, x[1][2]:0, x[2][0]:1, x[1][1]:1}, energy=-105.000000)

得られた結果から、.sampleスピン毎の結果、.energyエネルギー値、.constraints()/.subh制約確認、を取得することが出来る。

decoded_result[0].sample
----------------------------
<Output>
{'x[0][1]': 0,
 'x[0][0]': 0,
 'x[1][0]': 0,
 'x[2][1]': 0,
 'x[0][2]': 1,
 'x[2][2]': 0,
 'x[1][2]': 0,
 'x[2][0]': 1,
 'x[1][1]': 1}

decoded_result[0].energy
----------------------------
<Output>
-105.0

decoded_result[0].constraints()
----------------------------
<Output>
{'C2': (True, 0.0), 'C1': (True, 0.0)}

decoded_result[0].subh
----------------------------
<Output>
{'C1': 0.0, 'C2': 0.0}

スピン毎の結果を示す.sampleの出力は見難いので、行列で書き直す。

def ans(i, decoded_result):
  spin = decoded_result[i].sample
  ans = np.zeros(R*C).reshape(R,C) #結果を代入する3x3の空行列を作成
  for key,value in zip(spin.keys(), spin.values()):
      key = key.replace("]","").split("[") # x[0][1]→["x","0","1"]に変換
      key = [int(_v) for _v in key[1:]] # ["x","0","1"]→ [0,1]に変換
      ans[key[0]][key[1]]=value #結果の代入
  return ans
ans(1,decoded_result)
----------------------------
<Output>
array([[0., 1., 0.],
       [0., 0., 1.],
       [1., 0., 0.]])

制約確認も使用しやすく整理する。

def judge(i,decoded_result):
  judge_result=[]
  judge = decoded_result[i].constraints()
  for key in reversed(judge.keys()): #Constraints()は逆順のリストを返すため、reversedで逆転させて降順にする
    _temp = 0
    if judge[key][0] == True: #制約を満たす場合、つまりTrueが返される場合、「1」を返す。そうでない場合は「0」
      _temp = 1
    judge_result.append(_temp)
  return judge_result
judge(1,decoded_result)
----------------------------
<Output>
[1, 1]

必要な情報をひとまとめにする関数を定義する。

def summary(i, decoded_result):
  return [abs(decoded_result[i].energy), judge(i,decoded_result), ans(i,decoded_result)]
summary(1,decoded_result)
----------------------------
<Output>
[105.0, [1, 1], array([[0., 1., 0.],[0., 0., 1.],[1., 0., 0.]])]

3.2 検証1: 制約条件を満たさない場合の挙動

$\lambda$の値を変えながら、decodeの挙動を確認してみる。

###λの調整
feed_dict = {"lambda_e1":1,"lambda_c1":100,"lambda_c2":10}

###Qubo行列の取り出し
model = H.compile()
Qubo,offset = model.to_qubo(feed_dict=feed_dict)

###Dwaveでの計算
sampler_dw = EmbeddingComposite(dw_sampler)
result = sampler_dw.sample_qubo(Qubo, num_reads=10)

###Decode
decoded_result = model.decode_sampleset(sampleset=result,feed_dict=feed_dict)

for i in range(len(result)):
  print(summary(i,decoded_result))
----------------------------
<Output>
[125.0, [1, 0], array([[0., 0., 0.],[0., 1., 0.],[1., 0., 1.]])]
[125.0, [1, 0], array([[0., 0., 0.],[1., 0., 0.],[0., 1., 1.]])]
[105.0, [1, 0], array([[0., 0., 0.],[1., 1., 0.],[0., 0., 1.]])]
[105.0, [1, 0], array([[0., 1., 0.],[0., 0., 0.],[1., 0., 1.]])]
[105.0, [1, 1], array([[0., 1., 0.],[1., 0., 0.],[0., 0., 1.]])]
[105.0, [1, 1], array([[1., 0., 0.],[0., 1., 0.],[0., 0., 1.]])]

制約条件が不成立時には確かにFalseを返していることが確認できた。

3.3 検証2: エネルギー式が複数ある場合の挙動

エネルギー式が複数存在する場合でも、PyquboのConstraintで制約式とラベリングしたもの以外がエネルギー式として認識されるはずである。上の定式化のE(x)を単純に複製して、エネルギー値が2倍なっているか確認してみる。

### ハミルトニアンの作成
lambda_e2 = Placeholder("lambda_e2")
H = lambda_e1*E + lambda_c1*Constraint(C1,"C1") + lambda_c2*Constraint(C2,"C2") + lambda_e2*E #λのみ新たに定義してエネルギー式を複製

###λの調整
feed_dict = {"lambda_e1":1,"lambda_c1":100,"lambda_c2":100, "lambda_e2":1}

###Qubo行列の取り出し
model = H.compile()
Qubo,offset = model.to_qubo(feed_dict=feed_dict)

###Dwaveでの計算
sampler_dw = EmbeddingComposite(dw_sampler)
result = sampler_dw.sample_qubo(Qubo, num_reads=10)

###Decode
decoded_result = model.decode_sampleset(sampleset=result,feed_dict=feed_dict)

for i in range(len(result)):
  print(summary(i,decoded_result))
----------------------------
<Output>
[210.0, [1, 1], array([[0., 1., 0.],[0., 0., 1.],[1., 0., 0.]])]
[210.0, [1, 1], array([[0., 1., 0.],[1., 0., 0.],[0., 0., 1.]])]

確かにエネルギー式が2倍になっていることが確認できた。

3.4 検証3: 制約条件が「-1」でない場合の挙動

例題では$(\sum_{}x_{rc}-1)^2$としたが、$(\sum_{}x_{rc}-2)^2$と2カ所変更する場合の.constraints()の挙動を確認してみる。

### 第1項の定式化
E = 0
for r in range(R):
  for c in range(C):
    E = E - S[r][c]*x[r][c]

### 第2項の定式化
C1 = 0
for c in range(C):
  _c = 0
  for r in range(R):
    _c = _c + x[r][c]
  C1 = C1 + (_c-2)**2 # 「-1」から「-2」に変更した

### 第3項の定式化
C2 = 0
for r in range(R):
  _c = 0
  for c in range(C):
    _c = _c + x[r][c] 
  C2 = C2 + (_c-2)**2 # 「-1」から「-2」に変更した

### lambdaの設定
from pyqubo import Placeholder
lambda_e1,lambda_c1,lambda_c2 = Placeholder("lambda_e1"),Placeholder("lambda_c1"),Placeholder("lambda_c2")

### ハミルトニアンの作成
H = lambda_e1*E + lambda_c1*Constraint(C1,"C1") + lambda_c2*Constraint(C2,"C2")

###λの調整
feed_dict = {"lambda_e1":1,"lambda_c1":10,"lambda_c2":10} #制約式のλに小さい値10を入れる

###Qubo行列の取り出し
model = H.compile()
Qubo,offset = model.to_qubo(feed_dict=feed_dict)

###Dwaveでの計算
sampler_dw = EmbeddingComposite(dw_sampler)
result = sampler_dw.sample_qubo(Qubo, num_reads=10)

###Decode
decoded_result = model.decode_sampleset(sampleset=result,feed_dict=feed_dict)

for i in range(len(result)):
  print(summary(i,decoded_result),decoded_result[i].constraints()) #制約式の成立可否の生データを追加した。
----------------------------
<Output>
[265.0, [0, 0], array([[0., 1., 1.],[1., 1., 1.],[1., 1., 1.]])]
{'C2': (False, 2.0), 'C1': (False, 2.0)}

制約条件の数値が変わっても.constraints()が機能していることを確認。一方で(False, 2.0)の(2.0)が何を示すかが不透明。

3.5 一連の流れの関数化

###Decode
decoded_result = model.decode_sampleset(sampleset=result,feed_dict=feed_dict)

###結果の可視化
def ans(i, decoded_result):
  spin = decoded_result[i].sample
  ans = np.zeros(R*C).reshape(R,C) #結果を代入する3x3の空行列を作成
  for key,value in zip(spin.keys(), spin.values()):
      key = key.replace("]","").split("[") # x[0][1]→["x","0","1"]に変換
      key = [int(_v) for _v in key[1:]] # ["x","0","1"]→ [0,1]に変換
      ans[key[0]][key[1]]=value #結果の代入
  return ans

###制約達成の可視化
def judge(i,decoded_result):
  judge_result=[]
  judge = decoded_result[i].constraints()
  for key in reversed(judge.keys()): #Constraints()は逆順のリストを返すため、reversedで逆転させて降順にする
    _temp = 0
    if judge[key][0] == True: #制約を満たす場合、つまりTrueが返される場合、「1」を返す。そうでない場合は「0」
      _temp = 1
    judge_result.append(_temp)
  return judge_result

###結果の表示
def summary(i, decoded_result):
  return [abs(decoded_result[i].energy), judge(i,decoded_result), ans(i,decoded_result)]

おわりに

これまでPyquboのDecodeを使用できていなかったが、制約達成の確認、エネルギー値の取得は極めて便利なことがよくわかった。一方、制約達成をどのように確認しているかは不透明。時間ある時にメカニズムを確認してみたい。。。

3
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?