D-Wave SystemのOceanを使った足し算

D-Wave System社が提供する量子アニーリングでイジングやQUBO問題を解くPython製パッケージであるOceanを説明して、簡単な足し算のプログラミングを例に解説します。

1. 量子コンピューターのSDK

量子コンピュータのSDKは実機別に公開されてますが、それらをまとめて汎用的なもので提供してくれようとしてるものもあります。

主な量子コンピューターSDKの一覧

方式	名称	提供元	開発言語
量子アニーリング	Ocean	D-Wave System	Python
量子アニーリング	Wildqat	MDR社	Python
量子ゲート	QISKit	IBM	Python
量子ゲート	QDK	Microsoft	Q#
量子ゲート	Cirq	Google	Python
量子ゲート	ProjectQ	チューリッヒ工科大学	Python
量子ゲート	Forest	Rigetti Computing	Python
量子ゲート	Blueqat	MDR社	Python

2. イジングモデルとは、QUBOとは

量子コンピューターには主に量子ゲート方式と量子アニーリング方式の２種類があり、量子ゲート方式は汎用の量子計算が可能であるが、量子アニーリング方式はイジング・QUBO問題に見られる組み合わせ最適化問題に特化した方式です。
現在の量子アニーリング方式では、D-WaveSystemsが2048量子ビットまで実現できている。
どちらの方式もハミルトニアンを構築するまでは同じだが、量子ゲート方式では、解きたい問題のハミルトニアンを量子アルゴリズムに従い、ビット反転演算Xや位相反転演算Z、制御NOTゲートCNOTなどの量子ゲートを使って組み立てる必要がある。一方、量子アニーリング方式は、解きたい問題のハミルトニアンをイジングモデルやQUBOモデルのバイナリ2次モデルに変換すれば、あとは量子アニーラーに任せるだけで、そのハミルトニアンの最小エネルギーを求めるので複雑な量子ゲートを組み立てる必要はなく、非常に簡単で短い量子コードで作成することができます。
イジングモデルやQUBOモデルの詳しい説明は以下をご参考ください。
量子アニーリング、イジングモデルとフレームワーク
QUBOとイジングモデルは、具体的にどう対応するのか

3. Oceanとは

Oceanは、D-Wave Systemが提供する量子アニーリングでイジングやQUBO問題を解くPython製パッケージです。
ライセンスは、Apache License 2.0で提供されています。

D-Wave SystemのGPUを使うとお金がかかってしまうのですが、Oceanを使えば、D-Wave SystemのGPU（量子コンピューター）とローカルマシン上のCPU/GPU（古典的コンピューター）を簡単な切り替えで使うことができるので、どんなものか試したかったり、初期段階の開発や試験運用では無料でできます。これに比べて、同じ量子アニーリングでも富士通のデジタルアニーラーは全て有料の商用サービスなので、こういったアプローチがとりずらく、この結果、興味のある人々へまずは試してもらうという機会を逃してしまうことになります。

4. Oceanのインストール

早速、Oceanをインストールしてみましょう。
Pythonパッケージなので、pipコマンドを使えば簡単にインストールされます。
Macの場合はターミナル、Windowsの場合はAnacondaプロンプトで以下のコマンドを実行するだけです。

$ pip install dwave-ocean-sdk
Collecting dwave-ocean-sdk
  Downloading https://files.pythonhosted.org/pac
   ：
uccessfully installed dimod-0.7.9 dwave-cloud-client-0.4.16 dwave-neal-0.4.3 dwave-networkx-0.6.6 dwave-ocean-sdk-1.0.3 dwave-qbsolv-0.2.9 dwave-system-0.5.4 dwavebinarycsp-0.0.9 homebase-1.0.1 minorminer-0.1.6 ortools-6.7.4973 penaltymodel-0.15.3 penaltymodel-cache-0.3.2 penaltymodel-mip-0.1.3

念の為にインストールできたかも確認しましょう。

$ python -c 'import dimod; print(dimod.ExactSolver())'
<dimod.reference.samplers.exact_solver.ExactSolver object at 0x107abbc50>

詳細については本家のマニュアルを参照してください。
https://docs.ocean.dwavesys.com/en/latest/overview/install.html

5. 足し算（1+2）をQUBOモデルにしてみよう

今回は$ 1 + 2 = x $となるようなxを求めます。

QUBOにおいて足し算は上記右辺から左辺を引いてその最小値を求めることに相当しますので、$ x - ( 1 + 2 ) = 0$の目的関数を最小にするような変数xを求める問題になります。

$\begin{align} f(x) &= \{x-(1+2)\}^2 \\ &= (x-3)^2 \end{align}$

QUBOモデルでは、目的関数$f(x)$には、バイナリ変数$x\in\{0,1\}$が使われ、スカラー係数$ a_i,b_{i,j} $を用いたスカラー形式にできる。

$\begin{align} f(x) &= \sum_{i}{a_i}{x_i} + \sum_{i<j}b_{ij}x_{i}x_{j} \end{align}$

QUBOは、行列の対角要素$ Q_{i,j} $を一次係数とし、非対角係数$ Q_{i,j} $を二次係数とした行列形式で表せれる。

$\begin{align} f(x) &= \sum_{i}{Q_i}{x_i} + \sum_{i<j}Q_{ij}x_{i}x_{j} = \sum_{i<j}x_{i}Q_{ij}x_{j} \end{align}$

一方、xは量子ビットを使って、
$\begin{align} x = x_0 + 2x_1 \end{align}$
という二進数表記ができますので、さらにこれを上記のコスト関数に代入すると、

$\begin{align}
f(x) &= \bigl(x_0 + 2x_1\bigl)^2 - 6\bigl(x_0+2x_1\bigl) + 9 \\
&= x_0^2 + 4{x_0}{x_1} + 4x_1^2 - 6x_0 -12x_1 + 9 \end{align}$

と展開されます。ここで、QUBOモデルの変数$x\in\{0,1\}$はバイナリ値{0,1}を取るので、二乗の項はすべて指数がとれます。

$\begin{align} x_0^2 = x_0 ,\quad x_1^2 = x_1 \end{align}$

これにより、目的関数はオフセット値9とスカラー形式の目的関数

$\begin{align} f(x) &= \sum_{i}{a_i}{x_i} + \sum_{i<j}b_{ij}x_{i}x_{j}+9 \\ &=-5x_{0}-8x_{1}-4x_{0}x_{1}+9 \end{align}$

となります。これをQUBOの行列表記にすると、下記のようになります。

$\begin{align} f(x) &= \sum_{i<j}{x_i}Q_{ij}x_{j} + 9 \\ &=
\begin{bmatrix}
x_0 & x_1
\end{bmatrix}
\begin{bmatrix}
-5 & 4 \\
0 & -8
\end{bmatrix}
\begin{bmatrix}
x_0 \\
x_1
\end{bmatrix}+9
\end{align}$

6. Oceanパッケージのサンプラー種類とコード例

ここから先はOceanパッケージの各サンプラーを紹介し、コード例も紹介します。
問題が初歩過ぎてつまらないかもしれませんが、どの辺が変わるかだけでも参考にしていただければ幸いです。

種類	ツール	サンプラー	概要
古典的	dimod	ExactSovler	簡単な問題で正確な解を求めるためのサンプラー
古典的	dimod	RandomSampler	ランダムに探索するサンプラー
古典的	dimod	SimulatedAnnealingSampler	シミュレーテッド・アニーリング（焼きなまし法）で求めるサンプラー
古典的	neal	SimulatedAnnealingSampler	シミュレーテッド・アニーリングで、dimodよりう詳細な設定が行えるサンプラー
量子実機	dwave.system	DWaveSampler	D-Wave System社量子コンピューターのハードウェア実機（QPU）を使って解を求めるサンプラー。有償利用となります。
量子実機	dwave.cloud	Solver	上記より細かい制御が可能なもの。こちらも有償利用となります。
古典的/量子実機	dimod	QBSolv	大きなQUBO問題を解くソルバー。デフォルトは古典的ソルバーで探索するが、QPUへの切り替えもできる。

6.1. ExractSolver

簡単な問題で正確な解を求めたい時はこちらを利用します。

QUBOモデルの行列Qの値の$\begin{align}\begin{bmatrix}
-5 & 4 \\
0 & -8
\end{bmatrix}\end{align}
$と、オフセット値の9を用いて、以下のように設定します。

from dimod import *
Q = {(0, 0):-5,
     (0, 1): 4,
     (1, 1): -8
    } #QUBOモデル行列Q
b = BinaryQuadraticModel.from_qubo(Q, 9.0)
r = SimulatedAnnealingSampler().sample(b)
for s, e in r.data(['sample','energy']):
    print(s,'\t E = ', e)

[実行結果]

{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0

最小エネルギー$E=0$を与えるバイナリ変数は、$x_0=1,x_1=1$なので次のように計算ができます。

$x=x_0+2x_1=1+2 \times 1=3$

6.2. RandomSampler

これはランダムに探索するサンプラーです。開発時の試験運営に利用します。

from dimod import *
Q = {(0, 0):-5,
     (0, 1): 4,
     (1, 1): -8
    }
b = BinaryQuadraticModel.from_qubo(Q, 9.0)
r = RandomSampler().sample(b) # ここだけ修正
for s, e in r.data(['sample','energy']):
    print(s,'\t E = ', e)

[実行結果]

{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 0, 1: 1}     E =  1.0
{0: 1, 1: 0}     E =  4.0
{0: 1, 1: 0}     E =  4.0
{0: 0, 1: 0}     E =  9.0
{0: 0, 1: 0}     E =  9.0
{0: 0, 1: 0}     E =  9.0
{0: 0, 1: 0}     E =  9.0

一番エネルギーの小さい（$E=0$）バイナリ変数は、$x_0=1,x_1=1$なのであってます。

6.3. SimulatedAnnealingSampler

これは焼きなまし法で近似解を求めるサンプラーです。開発時の試験運営に利用します。

from dimod import *
Q = {(0, 0):-5,
     (0, 1): 4,
     (1, 1): -8
    }
b = BinaryQuadraticModel.from_qubo(Q, 9.0)
r = SimulatedAnnealingSampler().sample(b) # ここだけ修正
for s, e in r.data(['sample','energy']):
    print(s,'\t E = ', e)

[実行結果]

問題が簡単なので全部あってますね。

{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0

6.4. neal.SimulatedAnnealingSampler

これは焼きなまし法で解を求めるサンプラーです。開発時の試験運営に利用します。
Simodツールに比べると細かいパラメータの設定が可能です。

線形

from dimod import *
import neal
Q = {(0, 0):-5,
     (0, 1): 4,
     (1, 1): -8
    }
b = BinaryQuadraticModel.from_qubo(Q, 9.0)
r = neal.SimulatedAnnealingSampler().sample(b,
     seed=1111,
     sweeps=20,
     beta_schedule_type='linear') # ここだけ修正
for s, e in r.data(['sample','energy']):
    print(s,'\t E = ', e)

[実行結果]

{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0

非線形

from dimod import *
import neal
Q = {(0, 0):-5,
     (0, 1): 4,
     (1, 1): -8
    }
b = BinaryQuadraticModel.from_qubo(Q, 9.0)
r = neal.SimulatedAnnealingSampler().sample(b,
     seed=1111,
     sweeps=20,
     beta_schedule_type='geometric') # ここだけ修正
for s, e in r.data(['sample','energy']):
    print(s,'\t E = ', e)

[実行結果]

{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0
{0: 1, 1: 1}     E =  0.0

6.5 DWaveSampler（有料）

D-Wave Systemの量子アニーリングの実機を使う方法です。
ただし、日本からは個人で使わせてもらえる手段がわかりませんのでコード例だけ記述します。

from dimod import *
from dwave.system.samplers import *
Q = {(0, 0):-5,
     (0, 1): 4,
     (1, 1): -8
    }
b = BinaryQuadraticModel.from_qubo(Q, 9.0)
r = DWaveSampler().sample(b) # ここだけ修正
for s, e in r.data(['sample','energy']):
    print(s,'\t E = ', e)

[実行結果]

（有料機能を使える方は結果を教えてくれたら嬉しいです。）

6.6 D-Wave Cloud Solve（有料）

D-Wave Systemの量子アニーリングの実機を使う方法です。
ただし、日本からは個人で使わせてもらえる手段がわかりませんのでコード例だけ記述します。

from dimod import *
from dwave.system.samplers import *
from dwave.system.composites import EmbeddingComposite

Q = {(0, 0):-5,
     (0, 1): 4,
     (1, 1): -8
    }
b = BinaryQuadraticModel.from_qubo(Q, 9.0)

token = 'ここにあなたのTokenを記載' # ここだけ修正
r = EmbeddingComposite(DWaveSampler(token=token)).sample(b)
for s, e in r.data(['sample','energy']):
    print(s,'\t E = ', e)

[実行結果]

{0: 1, 1: 1}     E =  0.0

*2019/6/29 有料機能（D-wave 実機）にて検証済み

6.7. QBSolv

大きな問題をサブ問題に分割するソルバーです。
デフォルトはローカルでの古典的で実行されますが、D-WaveのQPUを利用することも可能です（その場合は有償）。

from dimod import *
from dwave_qbsolv import QBSolv
Q = {(0, 0):-5,
     (0, 1): 4,
     (1, 1): -8
    }
b = BinaryQuadraticModel.from_qubo(Q, 9.0)
r = QBSolv().sample(b) # ここだけ修正
for s, e in r.data(['sample','energy']):
    print(s,'\t E = ', e)

[実行結果]

{0: 1, 1: 1}     E =  0.0

Oceanはサンプラーを変えるだけで、すぐに切り替えができるのでとても便利ですね。

7. （参考）Wildqat SDK

日本を代表する量子コンピューター企業にMDR社があります。MDR社は量子アニーリングのSDKとしてWildqatを公開しています。また、ゲートモデル用ではBlueqatも公開してます。WildqatでもQUBOを使って解を求めることができるか確認してみました。

WildqatもPython製なのでpipにより簡単にインストールができます。

$pip install wildqat

コード例は以下になります。オフセットの設定が見当たりませんでしたが、最小コストを求めるのには変わらないのでQUBOだけ設定して解を求めます。

import wildqat as wq
a = wq.opt()
a.qubo = [[-5,4],[0,-8]]
a.sa()

[実行結果]
正しい結果（$x_0,x_1$）が出力されました。

[1, 1]

8. 便利なユーティリティ関数

Oceanには便利なユーティリティ関数が存在するので紹介します。

ユーティリティー関数	概要
ising_to_qubo(h,J,offset=0.0)	イジング問題をQUBO問題へ変換
qubo_to_ising(Q,offset=0.0)	QUBO問題をイジング問題へ変換
ising_energy(sample,h,J,offset=0.0)	イジングモデルのサンプルからエネルギーを算出
qubo_energy(sample,Q,offset=0.0)	QUBOのサンプルからエネルギーを算出

QUBOをイジングモデルに変換して実行した例

Q = {(0, 0): -5.0, (0, 1): 4.0, (1, 1): -8.0}
print(qubo_to_ising(Q, 9.0))

イジングモデルへの変換結果

({0: -1.5, 1: -3.0}, {(0, 1): 1.0}, 3.5)

イジングモデルからバイナリ二次モデルを生成して、解を求める。

b = BinaryQuadraticModel({0: -1.5, 1: -3.0}, {(0, 1): 1.0}, 3.5,SPIN)
r = ExactSolver().sample(b)
for s, e in r.data(['sample','energy']):
    print(s,'\t E = ', e)

あたり前ですが、実行結果はQUBOの時と同じです。

{0: 1, 1: 1}     E =  0.0
{0: -1, 1: 1}    E =  1.0
{0: 1, 1: -1}    E =  4.0
{0: -1, 1: -1}   E =  9.0

OceanのBinaryQuadraticModelクラス（バイナリ二次モデル）にも、便利なメソッドが揃ってましたので、こちらも紹介します。

メソッド	概要
from_ising(h,offset])	イジングモデルからBQMへ変換
from_qubo(Q[,offset])	QUBOモデルからBQMへ変換
from_numpy_vectors()	NumPyベクトルからBQMへ変換
from_numpy_matrix(mat)	NumPy配列からBQMへ変換
from_pandas_dataframe(df)	pandasからBQMへ変換
from_coo(obj[,type])	座標からBQMへ変換する
to_ising()	BQMからイジングモデルへ変換
to_qubo()	BQMからQUBOモデルへ変換
to_numpy_vectors(])	BQMからNumPyベクトルへ変換
to_numpy_matrix([])	BQMからNumPy配列へ変換
to_pandas_dataframe()	BQMからpandasデータフレームへ変換
to_coo([fp,typer])	BQMから座標フォーマットへ変換
to_networkx_graph([])	BQMからNetworkxグラフへ変換

参考文献

https://www.slideshare.net/masayukiminato/jan-2018-85514793
Ocean量子アニーリング入門 中山茂