QAP.04：D-Wave QUBO用行列を直接作る【量子コンピュータ/アニーリング@Python/D-Wave】

【はじめに】

前回は、D-Wave Leapの使い方をざっとまとめた。
簡単なプログラムとして「Binaryオブジェクト(dimod版)」を生成し、「目的関数に対するBQMオブジェクト」を作成した後、「Sampler」にBQMを渡すことで量子アニーリング計算をした。

※再掲：D-Waveの全体の流れ

その他に、D-Waveでは直接「QUBO用のオブジェクト(※)」を作成して目的関数とすることで、量子アニーリング計算をさせることもできる。（※実体はpythonのdictオブジェクト）

今回はこのQUBOについてざっとまとめて、簡単なプログラムを作ってみる。

【QUBOについて】

QUBOという略語の正式名称とざっくりした説明をすると以下のような感じになる。

■QUBO
Quadratic Unconstrained Binary Optimization：二次制約なし二値最適化

■特徴
・決定変数は 0 or 1 のもの
・目的関数の次数が2次まで（＋多項式）
・明示的な制約条件はない

　　　　　　　　　　　↓
『「QUBO問題」とは上記特徴を持つ数式を最適化（最小化）する』ということ。

■目的関数の例
2xy + 2xz + 2yz - x - y - z + 1
→　(a, b, c)どれかが１・他が0の時に、最小値：0をとってくれる

■QUBO用行列

先ほど書いた「QUBO問題」に対して、「所定の形式に従った行列（QUBO用行列）」を用意すれば、D-Waveがいい感じに量子アニーリング計算をしてくれる。

「QUBO用の行列」は具体的にいうと「目的関数（数式）」の「係数の値を使った上三角行列」である。
例えば、例に挙げている数式の場合は、次のような行列を作ればよい。

■ 2xy + 2xz + 2yz - x - y - z + 1のQUBO用行列をつくる

▲係数に対する上三角行列（定数「+1」は気にしなくてOK）

与えられた式の末尾の「+1」部分を無視しているが、これは「+1」部分が定数だからである。
→　目的関数をできるだけ小さくするための「x,y,z」の組み合わせがどうであれ常に変わらず「+1」なので気にしなくていい、ということ。

なお
・対角線上の成分を「linear」
・それ以外の成分を「quadratic」
と呼ぶ。

▲quadratic部分とlinear部分を表現する

【ここまでのまとめ】
長くなってしまったが、結局のところ、

「目的関数（数式）」に対して「係数の値を使った上三角行列」を作ることができれば、D-Waveが量子アニーリング計算してくれる、

ということ。

これを踏まえて簡単なプログラムを作ってみる。

【例題】

「2ab+2ac+2bc - a - b - c + 1 」を最小にする(a, b, c)の組み合わせを探す

■考え方
(a, b, c)の組み合わせの違いで値の変動が出るのは
「2ab+2ac+2bc - a - b - c」
部分だけ。

つまり、「+1」部分を除いた「2ab+2ac+2bc - a - b - c」を最小にする組み合わせを考える問題だと考える。
（※「+1」部分は定数なので、a,b,cの値の組み合わせがどうであれ変わらないから無視してよい）

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

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

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

【２】トークン設定

※再掲：
トークンを使ったD-Waveへのアクセスについては以下の通り

Configuring Access to Leap’s Solvers — Ocean Documentation 4.4.0 documentation

上記ドキュメント記載の通り、本来はアクセストークン漏洩防止など、セキュリティの観点から「コンフィグファイル」や「環境変数」に仕込む。

【注意】
今回は簡単なプログラムをちょっと試してみるだけなので、「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')]

【３】QUBO用行列を作成する（目的関数の作成）

QUBO用の行列は係数に対する上三角行列を作る。

「D-Wave」では「dict」＋「tuple」を使って表現する。

my_qubo = {
    ('x0','x1'):2,
    ('x0','x2'):2,
    ('x1','x2'):2,
    ('x0','x0'):-1, # tuple型で2値key(縦横の組み合わせ)を表現
    ('x1','x1'):-1,
    ('x2','x2'):-1,
}
print(my_qubo)

【ここまでの実行結果】
{('x0', 'x0'): -1,
 ('x0', 'x1'): 2,
 ('x0', 'x2'): 2,
 ('x1', 'x1'): -1,
 ('x1', 'x2'): 2,
 ('x2', 'x2'): -1}

▲QUBO用行列にあわせたdictオブジェクトの生成

■dictオブジェクト

組み込み型 — Python 3.10.0b2 ドキュメント

■tupleオブジェクト

組み込み型 — Python 3.10.0b2 ドキュメント

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

今回は「Sampler」として「DWaveSampler」を使用してみる。
作成したQUBO用行列フォーマットのデータは「sample.qubo()」に渡してコールすればよい。

from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite

# samplerオブジェクト生成
sampler = EmbeddingComposite(DWaveSampler(token=MY_DWAVE_TOKEN))

# qubo用行列フォーマットを渡して、アニーリングサイクル1000回で試行する
response = sampler.sample_qubo(my_qubo, num_reads=1000)

■sample.qubo()

dwave.system.samplers.DWaveSampler.sample_qubo — dwave-system 1.10.0 documentation

▲リファレンスの記載にあるが、実際は「BQM」に変換して「sample()」をコールしている。

【５】計算結果を確認する

計算結果は「SampleSetオブジェクト」として返ってくる。今回は「data()+list化」でイテレータを回して少し見やすく表示した。

for sample, energy, num_occurrences, chain_break_fraction in list(response.data()):
  print(sample, "Energy: ", energy, "Occurrences: ", num_occurrences)

【実行結果例】
{'x0': 0, 'x1': 1, 'x2': 0} Energy: -1.0 Occurrences: 372
{'x0': 0, 'x1': 0, 'x2': 1} Energy: -1.0 Occurrences: 294
{'x0': 1, 'x1': 0, 'x2': 0} Energy: -1.0 Occurrences: 333
{'x0': 1, 'x1': 0, 'x2': 1} Energy: 0.0 Occurrences: 1

※実行結果はやるたびに変わる

▲1000回の試行で
・(0,1,0)：372回（正解パターン）
・(0,0,1)：294回（正解パターン）
・(1,0,0)：333回（正解パターン）
・(1,0,1)：1回　（※）
の結果となった。（正解の組み合わせが大多数観測されている）

※何回も解探索をするので、「最適ではない解にたどり着くパターン」も観測されることもある。これが「DWaveSampler」を使った量子アニーリングの計算の特徴とも言える。

■SampleSet

Samples — Ocean Documentation 4.4.0 documentation

【６】全体コード

【例題】
「2ab+2ac+2bc - a - b - c + 1 」を最小にする(a, b, c)の組み合わせを探す
　　　　↓
今回の問題設定上では「2ab+2ac+2bc - a - b - c 」を最小にする(a,b,c)の組み合わせを探せばOK。

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

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

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

■ここからが全体コード

from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite


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


from dwave.cloud import Client

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


my_qubo = {
    ('x0','x1'):2,
    ('x0','x2'):2,
    ('x1','x2'):2,
    ('x0','x0'):-1, # tuple型で2値key(縦横の組み合わせ)を表現
    ('x1','x1'):-1,
    ('x2','x2'):-1,
}
print(my_qubo)




# samplerオブジェクト生成
sampler = EmbeddingComposite(DWaveSampler(token=MY_DWAVE_TOKEN))

# qubo用行列フォーマットを渡して、アニーリングサイクル1000回で試行する
response = sampler.sample_qubo(my_qubo, num_reads=1000)

# 結果の確認
for sample, energy, num_occurrences, chain_break_fraction in list(response.data()):
  print(sample, "Energy: ", energy, "Occurrences: ", num_occurrences)

【実行結果例】
{'x0': 0, 'x1': 1, 'x2': 0} Energy: -1.0 Occurrences: 372
{'x0': 0, 'x1': 0, 'x2': 1} Energy: -1.0 Occurrences: 294
{'x0': 1, 'x1': 0, 'x2': 0} Energy: -1.0 Occurrences: 333
{'x0': 1, 'x1': 0, 'x2': 1} Energy: 0.0 Occurrences: 1

【おまけ】BQMをつくって後からlinearやquadraticを積んでいく書き方

BQMを先に生成して、そこにlinear成分やquadratic成分を積んでいく方法もある。

■Binary Quadratic Models

Quadratic Models: Unconstrained — Ocean Documentation 4.4.0 documentation

from dimod import BinaryQuadraticModel
from dwave.system.samplers import LeapHybridBQMSampler # LeapHybridSamplerのエイリアス


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


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


bqm = BinaryQuadraticModel(3, 'BINARY')
print(bqm) 
# 出力例：変数のラベルは0,1,2...とついていく
# BinaryQuadraticModel({0: 0.0, 1: 0.0, 2: 0.0}, {}, 0.0, 'BINARY')



# BQMではオフセットも指定可能
bqm.offset = 1.0
bqm

# linear部分
bqm.add_linear(0,-1)
bqm.add_linear(1,-1)
bqm.add_linear(2,-1)

# quadratic部分
bqm.add_quadratic(0,1,2)
bqm.add_quadratic(0,2,2)
bqm.add_quadratic(1,2,2)

# QUBOフォーマットを出力してみる
print(bqm.to_qubo())

# LeapHybridBQMSamplerの生成
sampler = LeapHybridBQMSampler(token = MY_DWAVE_TOKEN)

# Solverに投げて計算実施
sampleset = sampler.sample(bqm) #LeapHybridの場合はnum_reads指定はできない

# 結果を確認
print(sampleset)

【実行結果例】
※bqm.to_qubo()部分
({(0, 0): -1.0,
(1, 0): 2.0,
(1, 1): -1.0,
(2, 0): 2.0,
(2, 1): 2.0,
(2, 2): -1.0},
1.0)

※print(sampleset)部分の出力
SampleSet(rec.array([([1, 0, 0], 0., 1)],
dtype=[('sample', '<i4', (3,)), ('energy', '<f8'), ('num_occurrences', '<i8')]), Variables([0, 1, 2]), {'qpu_access_time': 51842, 'charge_time': 2999713, 'run_time': 2999713, 'problem_id': 'XXXXXXXX-XXXX-XXXX-XXXX-XXXXXXXXXXXX'}, 'BINARY')

※note掲載上、problem_idはXXXでマスキング

▲今回は定数「+1」部分も含めて、(1,0,0)で最小値０をとる組み合わせが答えとして返ってきた（正解）