2023/5/23追記:本記事の内容は大幅にアップデートされています。以下の記事でまとめてあります。
量子アニーリングのQUBOで設定可能な条件式まとめ(保存版)
やること
これまでは「n個の量子ビットからm個を1にする」という条件設定ができればいろいろな問題が解ける!と思っていましたが、実際にいろいろなパズルを解こうとすると、それだけでは足りないことがわかりました。
ここではD-waveマシンを使いこなすための多彩な条件設定を学びます。
実行環境
WinPython3.6をおすすめしています。
Google Colaboratoryが利用可能です。
pip, import
必要なパッケージをインストールし、
pip install pyqubo
pip install dwave-ocean-sdk
インポートします。
from pyqubo import Binary
from dwave.system.composites import EmbeddingComposite
from dwave.system.samplers import DWaveSampler
基本となる共通コード
n個の量子ビットを用意する
3つの量子ビットを用意するコマンドは次のとおりでした。
#3つの量子ビットを用意する
q0, q1, q2 = Binary('q0'), Binary('q1'), Binary('q2')
D-waveでサンプリングしまくる
その後、条件式を設定してサンプリングしまくる部分のコードはこちらです。
#条件を記述する
H = (条件式)
#コンパイル
model = H.compile()
#print('model\n{}'.format(model))
qubo, offset = model.to_qubo()
#print('qubo\n{}'.format(qubo))
#D-waveの設定
sampler = EmbeddingComposite(DWaveSampler(endpoint='https://cloud.dwavesys.com/sapi',
token='your_token', #あなたのAPIトークンを記入
solver='DW_2000Q_5')) #利用可能なソルバーを指定
#サンプリングしまくる
response = sampler.sample_qubo(qubo, num_reads=1000)
#レスポンスの確認
print('response')
for (sample, energy, occurrence, _) in response.data():
print('Sample:{} Energy:{} Occurrence:{}'.format(list(sample.values()), energy, occurrence))
おさらい
n個の量子ビットからm個を1にする
「3個の量子ビットから2個を1にする」は次のとおりでした。
#条件を記述する
H = (q0 + q1 + q2 - 2)**2
Sample:[1, 1, 0] Energy:-4.0 Occurrence:244
Sample:[0, 1, 1] Energy:-4.0 Occurrence:350
Sample:[1, 0, 1] Energy:-4.0 Occurrence:403
Sample:[1, 0, 1] Energy:-4.0 Occurrence:1
Sample:[1, 0, 0] Energy:-3.0 Occurrence:2
3C2 =3通りが出力されています。どうしてこの設定でうまくいくのかを考えます。
H(ハミルトニアン)は量子ビット全体のエネルギー値みたいなものです。q0~q2はそれぞれ0か1をとります。D-waveは、Hがもっとも低くなるq0~q2の組み合わせを探すためのマシンです。
q0~q2にそれぞれ0か1を適当に代入してみるとH = 0, 1, 4を取りうることがわかります。
- H = 0になるのは、q0~q2のうち2個が1のとき
- H = 1になるのは、q0~q2のうち1個または3個が1のとき
- H = 4になるのは、q0~q2のうち0個が1のとき
ですから、上記の設定では「3個の量子ビットから2個を1にする」組み合わせが得られるのです。
「x = q0 + q1 + q2」と置くと、xはまさに「いくつを1にしたいか」を表しています。条件式は「H = (x – 2)**2」と書けます。
これをプロットしてみると、x = 2が最適解であることが明確にわかります。よって「3個の量子ビットから2個を1にする」が得られるというわけです。
n個の量子ビットからm個またはm+1個を1にする
「3個の量子ビットから1個または2個を1にする」は次のとおりでした。
#条件を記述する
H = (q0 + q1 + q2 - 1.5)**2
Sample:[0, 1, 0] Energy:-2.0 Occurrence:167
Sample:[1, 0, 0] Energy:-2.0 Occurrence:185
Sample:[0, 1, 1] Energy:-2.0 Occurrence:154
Sample:[0, 0, 1] Energy:-2.0 Occurrence:166
Sample:[1, 0, 1] Energy:-2.0 Occurrence:180
Sample:[1, 1, 0] Energy:-2.0 Occurrence:148
3C1 + 3C2 =6通りが出力されています。なぜこうなるのかを理解するために、H = (x – 1.5)**2のグラフを見てみます。
xは0, 1, 2, 3の離散値を取りうるため、x = 1.5はとれず、x = 1, 2が最適解であることがわかります。よって「3個の量子ビットから1個または2個を1にする」ができるのです。
n個の量子ビットを全部1または全部0にする
「3個の量子ビットから0個または3個を1にする」はどうすればよいでしょうか。さっきの式にマイナスを付けてみます。
#条件を記述する
H = - (q0 + q1 + q2 - 1.5)**2
Sample:[1, 1, 1] Energy:0.0 Occurrence:610
Sample:[0, 0, 0] Energy:0.0 Occurrence:390
できました。なぜこれでうまくいくかというと、H = – (x – 1.5)**2のグラフを見ると、
x = 0, 3におけるエネルギー値がもっとも低いからです。
n個の量子ビットをm個またはm+2個を1にする(失敗編)
ここからが本番です。「3個の量子ビットから1個または3個を1にする」はどうすればよいでしょうか。
明らかにダメな例
一瞬、「1と3の中間である2を引けばいいんだろう」と思ったのですが、
#条件を記述する
H = (q0 + q1 + q2 - 2)**2
これは「2個を1にする」やつなので、実行するまでもなくダメです。
足してみる
「1個を1にする」と「3個を1にする」を足してみたらどうでしょうか?
#条件を記述する
H = (q0 + q1 + q2 - 1)**2 + (q0 + q1 + q2 - 3)**2
Sample:[1, 0, 1] Energy:-8.0 Occurrence:365
Sample:[1, 1, 0] Energy:-8.0 Occurrence:285
Sample:[0, 1, 1] Energy:-8.0 Occurrence:350
これも「2個を1にする」になってしまいました。なぜかというと、H = (x – 1)**2 + (x – 3)**2のグラフはx = 2で落ち込むからです。
2次関数はどんなに足しても2次関数ですので、足す方法では不可能です。
かけてみる
かけてみたらどうでしょうか?理論上はこれでいけるはずです。
#条件を記述する
H = (q0 + q1 + q2 - 1)**2 * (q0 + q1 + q2 - 3)**2
Sample:[1, 0, 0, 1] Energy:-13.0 Occurrence:31
Sample:[1, 0, 0, 1] Energy:-13.0 Occurrence:12
Sample:[1, 1, 1, 1] Energy:-9.0 Occurrence:1
Sample:[1, 1, 0, 1] Energy:-8.0 Occurrence:616
Sample:[1, 1, 0, 1] Energy:-8.0 Occurrence:243
Sample:[0, 1, 1, 1] Energy:-5.0 Occurrence:97
あれ、いろいろおかしいです。量子ビットは3個しか用意していないのに4個に増えています。解も実行するたびに滅茶苦茶な結果になります。
ちなみに、上の式を展開した
#条件を記述する
H = (q0+q1+q2)**4 - 8*(q0+q1+q2)**3 + 22*(q0+q1+q2)**2 - 24*(q0+q1+q2) + 9
#条件を記述する
H = q0**4 + 4*q0**3*q1 + 4*q0**3*q2 - 8*q0**3 + 6*q0**2*q1**2 + 12*q0**2*q1*q2 - 24*q0**2*q1 \
+ 6*q0**2*q2**2 - 24*q0**2*q2 + 22*q0**2 + 4*q0*q1**3 + 12*q0*q1**2*q2 - 24*q0*q1**2 \
+ 12*q0*q1*q2**2 - 48*q0*q1*q2 + 44*q0*q1 + 4*q0*q2**3 - 24*q0*q2**2 + 44*q0*q2 \
- 24*q0 + q1**4 + 4*q1**3*q2 - 8*q1**3 + 6*q1**2*q2**2 - 24*q1**2*q2 + 22*q1**2 \
+ 4*q1*q2**3 - 24*q1*q2**2 + 44*q1*q2 - 24*q1 + q2**4 - 8*q2**3 + 22*q2**2 - 24*q2 + 9
も同様にダメでした。
関数の形状はちゃんとx = 1, 3で落ち込んでいるのに、なぜうまく行かないのでしょうか。
そこで、マシンに投げる前のQUBO(読み方:きゅーぼ)(量子ビット間の重み)を見てみると、
qubo
{('q0', 'q1'): 10.0, ('q0', 'q2'): 5.0, ('q0', 'q0*q2'): -10.0, ('q1', 'q2'): 10.0, ('q0*q2', 'q1'): -12.0, ('q0*q2', 'q2'): -10.0, ('q0', 'q0'): -9.0, ('q1', 'q1'): -9.0, ('q2', 'q2'): -9.0, ('q0*q2', 'q0*q2'): 25.0}
「q0*q2」という謎の補助ビットを用意してくれていることがわかりました。ということは、おそらく[1, 0, 0, 1] = [q0, q0*q2, q1, q2]なので[q0, q1, q2] = [1, 0, 1]しか得られなかったということです。(実行するたびに変わります)
「1個または3個」とは程遠い結果になりました。D-waveマシンはハードウェアの都合上、2量子ビットよりも多くの相互作用を扱うことができません。その場合には補助ビットを用意する必要がありますが、単にHに条件式を入れる方式ではうまいことやってくれないということが分かりました。
n個の量子ビットをm個またはm+2個を1にする(解決編)
2023/5/23追記:この章の内容は大幅にアップデートされているので読まなくて大丈夫です。以下の記事でまとめてあります。
量子アニーリングのQUBOで設定可能な条件式まとめ(保存版)
3個の量子ビットから1個または3個を1にする
泣く泣く、日本最大の(もしかしたら世界最大の?)量子コンピュータのオンラインコミュニティであるSlack「Blueqat(クリックで参加)」で有識者に助けを求めました。
こちらの記事の方法でできるかもよ、と教えていただきました。
結論から書くと、
#5つの量子ビットを用意する
q0, q1, q2, q3, q4 = Binary('q0'), Binary('q1'), Binary('q2'), Binary('q3'), Binary('q4')
#条件を記述する
aux0 = q3 #補助ビット
aux1 = q4 #補助ビット
H = (4.0*aux1)+(6.0*q0)+(-2.0*q1)+(6.0*q2)+(-4.0*aux0*aux1)+(-4.0*aux0*q0)+(4.0*aux0*q1)+(-4.0*aux0*q2)+(-4.0*aux1*q0)+(4.0*aux1*q1)+(-4.0*aux1*q2)+(-4.0*q0*q1)+(4.0*q0*q2)+(-4.0*q1*q2)+(2.0)
Sample:[0, 1, 0, 0, 0] Energy:-2.0 Occurrence:62
Sample:[1, 1, 1, 1, 1] Energy:-2.0 Occurrence:23
Sample:[1, 0, 0, 1, 1] Energy:-2.0 Occurrence:175
Sample:[0, 0, 1, 1, 1] Energy:-2.0 Occurrence:56
Sample:[0, 0, 1, 1, 1] Energy:-2.0 Occurrence:1
Sample:[1, 0, 1, 1, 1] Energy:0.0 Occurrence:683
これで「3個の量子ビットから1個または3個を1にする」ができました。左側の3つがq0, q1, q2で、右側の2つが補助ビットです。条件式は、2つの補助ビットを巧みに使って2次より大きな相互作用が現れないようにしてあります。
係数の求め方
この魔法のような係数の求め方は、先程の記事を参考にさせていただきました。
以下はGoogle Colaboratoryでの実行を強く推奨します。ローカルの実行環境にpipしようとすると既存のパッケージ環境がなんか破壊されます(破壊されました)。
!pip install dwavebinarycsp
!pip install dwavebinarycsp[maxgap]
怒られるので、途中にある「RESTART RUNTIME」を押します。
import dwavebinarycsp
csp = dwavebinarycsp.ConstraintSatisfactionProblem(dwavebinarycsp.BINARY)
csp.add_constraint(lambda a, b, c: (a + b + c == 1) or (a + b + c == 3), ['q0', 'q1', 'q2'])
bqm = dwavebinarycsp.stitch(csp)
print(bqm)
BinaryQuadraticModel({'aux0': 0.0, 'aux1': 4.0, 'q0': 6.0, 'q1': -2.0, 'q2': 6.0}, {('aux0', 'aux1'): -4.0, ('aux0', 'q0'): -4.0, ('aux0', 'q1'): 4.0, ('aux0', 'q2'): -4.0, ('aux1', 'q0'): -4.0, ('aux1', 'q1'): 4.0, ('aux1', 'q2'): -4.0, ('q0', 'q1'): -4.0, ('q0', 'q2'): 4.0, ('q1', 'q2'): -4.0}, 2.0, Vartype.BINARY)
これが係数です。空気を読んで対応させてください。どんな計算を行っているかはまったく分かりません。。。
感謝の思いを込めて、スペシャルなコードを公開させていただきます。↑の「空気を読んで対応させてください」という面倒な作業を自動的にやってくれるコードです。
#Qビット名を取得
names = []
for name in bqm:
names.append(name)
print(names)
#単体の係数
for name in names:
print(name, bqm.linear[name])
#2体の係数
for i in range(len(names)):
for j in range(i + 1, len(names)):
try:
print(names[i], names[j], bqm.adj[names[i]][names[j]])
except:
pass
#オフセット
print(bqm.offset)
#H=の形式で出力
command = 'H = 0 '
#単体の係数
for name in names:
#ゼロでなければ式に追加
if bqm.linear[name] != 0.0:
command += '+({}*{})'.format(bqm.linear[name], name)
#2体の係数
for i in range(len(names)):
for j in range(i + 1, len(names)):
try:
#ゼロでなければ式に追加
if bqm.adj[names[i]][names[j]] != 0.0:
command += '+({}*{}*{})'.format(bqm.adj[names[i]][names[j]], names[i], names[j])
except:
pass
#オフセット
command += '+({})'.format(bqm.offset)
print(command)
H = 0 +(4.0*aux1)+(6.0*q0)+(-2.0*q1)+(6.0*q2)+(-4.0*aux0*aux1)+(-4.0*aux0*q0)+(4.0*aux0*q1)+(-4.0*aux0*q2)+(-4.0*aux1*q0)+(4.0*aux1*q1)+(-4.0*aux1*q2)+(-4.0*q0*q1)+(4.0*q0*q2)+(-4.0*q1*q2)+(2.0)
結果をコピペすればそのまま条件式として使えます。(やったね)
3個の量子ビットから0個または2個を1にする
Colabの方で係数を求め、
csp.add_constraint(lambda a, b, c: (a + b + c == 0) or (a + b + c == 2), ['q0', 'q1', 'q2'])
H = 0 +(-4.0*aux0)+(10.0*q0)+(-6.0*q1)+(-6.0*q2)+(-4.0*aux0*aux1)+(-4.0*aux0*q0)+(4.0*aux0*q1)+(4.0*aux0*q2)+(-4.0*aux1*q0)+(4.0*aux1*q1)+(4.0*aux1*q2)+(-4.0*q0*q1)+(-4.0*q0*q2)+(4.0*q1*q2)+(8.0)
D-waveマシンに投げます。
#条件を記述する
aux0 = q3 #補助ビット
aux1 = q4 #補助ビット
H = (-4.0*aux0)+(10.0*q0)+(-6.0*q1)+(-6.0*q2)+(-4.0*aux0*aux1)+(-4.0*aux0*q0)+(4.0*aux0*q1)+(4.0*aux0*q2)+(-4.0*aux1*q0)+(4.0*aux1*q1)+(4.0*aux1*q2)+(-4.0*q0*q1)+(-4.0*q0*q2)+(4.0*q1*q2)+(8.0)
Sample:[1, 1, 0, 1, 1] Energy:-8.0 Occurrence:392
Sample:[0, 0, 0, 1, 1] Energy:-8.0 Occurrence:235
Sample:[1, 0, 1, 1, 1] Energy:-8.0 Occurrence:158
Sample:[0, 1, 1, 0, 0] Energy:-8.0 Occurrence:209
Sample:[1, 1, 1, 1, 0] Energy:-6.0 Occurrence:6
・・・できています、魔法です。
4個の量子ビットから1個または3個を1にする
同様に係数を求めました。
#4つの量子ビットを用意する
q0, q1, q2, q3, q4, q5 = Binary('q0'), Binary('q1'), Binary('q2'), Binary('q3'), Binary('q4'), Binary('q5')
#条件を記述する
aux0 = q4 #補助ビット
aux1 = q5 #補助ビット
H = (-6.0*q0)+(-6.0*q1)+(-2.0*q2)+(-2.0*q3)+(-8.0*aux1)+(4.0*q0*q1)+(4.0*q0*aux0)+(4.0*q0*aux1)+(4.0*q1*aux0)+(4.0*q1*aux1)+(4.0*q2*q3)+(-4.0*q2*aux0)+(4.0*q2*aux1)+(-4.0*q3*aux0)+(4.0*q3*aux1)+(10.0)
Sample:[0, 0, 1, 0, 1, 1] Energy:-10.0 Occurrence:57
Sample:[1, 1, 0, 1, 0, 0] Energy:-10.0 Occurrence:59
Sample:[1, 1, 0, 1, 0, 0] Energy:-10.0 Occurrence:48
Sample:[0, 0, 0, 1, 1, 1] Energy:-10.0 Occurrence:29
Sample:[0, 1, 1, 1, 1, 0] Energy:-10.0 Occurrence:24
Sample:[1, 1, 1, 0, 0, 0] Energy:-10.0 Occurrence:77
Sample:[0, 1, 0, 0, 0, 1] Energy:-10.0 Occurrence:102
Sample:[0, 1, 1, 1, 1, 0] Energy:-10.0 Occurrence:55
Sample:[1, 1, 1, 0, 0, 0] Energy:-10.0 Occurrence:137
Sample:[1, 0, 1, 1, 1, 0] Energy:-10.0 Occurrence:32
Sample:[1, 0, 1, 1, 1, 0] Energy:-10.0 Occurrence:53
Sample:[1, 0, 0, 0, 0, 1] Energy:-10.0 Occurrence:123
Sample:[0, 1, 0, 1, 1, 1] Energy:-8.0 Occurrence:17
Sample:[0, 1, 0, 1, 0, 1] Energy:-8.0 Occurrence:66
Sample:[0, 1, 1, 0, 1, 1] Energy:-8.0 Occurrence:46
Sample:[1, 0, 0, 1, 0, 1] Energy:-8.0 Occurrence:75
やたらに重複が多いですが、4C1 + 4C3 =8通りが出てきます。まじかよ!
4個の量子ビットから0個または3個を1にする
応答エラー
4個の量子ビットから1個または4個を1にする
応答エラー
4個の量子ビットから0個または2個または4個を1にする
期待した出力が得られませんでした。
まとめ
n個の量子ビットについて、できる(と思われる)ことはこちらです。
- n個の量子ビットからm個を1にする
- n個の量子ビットからm個またはm+1個を1にする
- n個の量子ビットを全部1または全部0にする
一般に、以下はできない(と思われる)。
n個の量子ビットからm個またはk個を1にする(m<k, k≠m+1)
ただし、次に限って成功しました。
3個の量子ビットから1個または3個を1にする3個の量子ビットから0個または2個を1にする4個の量子ビットから1個または3個を1にする
これらの限られた手札を使って、様々な最適化問題やパズル問題にチャレンジしていきたいと思います。
2023/5/23追記:本記事の内容は大幅にアップデートされています。以下の記事でまとめてあります。
量子アニーリングのQUBOで設定可能な条件式まとめ(保存版)