はじめに
D-Wave Systems社（以下、D-Wave社）のマシン等を使って何か始める人の参考になればと思い、とりあえずQUBO問題を自分で作成し、なんらかのソルバーを使って解けるようになるレベルを目指してその手順を書きました。誤りのご指摘はコメント欄にてお願いいたします。
QUBOとは
Quadratic Unconstrained Binary Optimization (QUBO)は2値の制約なし２次最適化問題のことです。次のように表されます。
$$ \min_{x\in{0,1}^n}f(x)$$
ただしここで
$$f(x)=\sum_{i\leq j}^nq_{i,j}x_ix_j$$
とします。言い換えると、「任意の係数 $q_{i,j}\in\mathbb R$($i,j$は$i\leq j$をみたす$1\sim n$の整数)が与えられた時に、関数$f(x)$を最小にするような２値のベクトル$x\in{0,1}^n$を求めなさい」ということです。
これらは一般にNP困難と呼ばれる難しい問題です[Murty & Kabadi, 1987]。
例題1
次の問題の解を全列挙によって求めましょう。
$$ \min_{x\in{0,1}^4}f(x)$$
ただしここで
$$
\begin{align}
f(x) &= x_1-1x_2-2x_3 \\
& +2x_1x_2-x_1x_3 \\
& -2x_2x_3-x_2x_4 \\
& +2x_3x_4\end{align}$$
とします。以下に定義した関数f(x)とallX(n)、numpyライブラリのapply_along_axis関数を利用します。
import numpy as np
def f(x):
'''
x: 長さ4で要素に0, 1のみを持つNumpy行列
'''
assert np.all((x == 0) | (x == 1))
Q = np.array([
[ 1, 2, -1, 0],
[ 0, -1, -2, -1],
[ 0, 0, -2, 2],
[ 0, 0, 0, 0]
])
return np.dot(x.T, np.dot(Q, x))
def allX(n):
'''
長さnの0,1ベクトルの全パターンを列挙して2次元のNumpy行列として返す。
'''
Xs = np.arange(2**n).reshape((-1,1))
mask = 1 << np.arange(n-1, -1, -1).reshape((1,-1))
return (Xs & mask) / mask
1次式による制約の付加
QUBOはその名の通り「制約なし」問題ですが、特定の種類の制約であれば簡単にQUBOに組み込むことができます。そのひとつが、1次式による制約です。
1次式による制約条件
$$ g(x) = 0$$
ただし
$$g(x) = c_0 + \sum_{i=1}^n c_i x_i$$
というものを考えます。この一時式が満たされるような$x$をQUBOの解として表現したい場合、$g(x)$を二乗した関数$g(x)^2$を考えれば、制約が満たされる時のみ最小値0をとることがわかります。
例えば、「$x_1, x_2, x_3$のうちひとつだけが1になる」という条件は$g(x) = x_1 + x_2 + x_3 - 1 = 0$という線形制約で書かれるので、左辺を二乗した関数
$$
\begin{align}
g(x)^2 &= (x_1 + x_2 + x_3 - 1)^2 \\
&= x_1 + x_2 + x_3 + 1 \\
&+ 2 x_1 x_2 + 2x_2x_3 + 2x_3x_1 \\
&-2x_1 - 2x_2 - 2x_3 \\
&= 1 - x_1 - x_2 - x_3 + 2x_1x_2 + 2x_2x_3 + 2x_3x_1
\end{align}
$$
を十分大きな係数$\lambda$をかけてもとの問題の$f(x)$に足してあげると、この制約を満たしつつ$f(x)$を最小化するものだけが解として出てくるようになります。
つまり、QUBO形式の問題
$$ \min_{x\in{0,1}^n}f(x) + \lambda g(x)^2$$
ただし$\lambda$は十分大きいものとする。は、制約付きの問題
$$\min_{x\in{0,1}^n|g(x)=0}f(x)$$
と等価であるということができます。
実際は$\lambda$はいくらでも大きい値にすればよいというわけではなく、これから紹介するようなソルバーを利用したときに良い解が出てくるように$f(x)$がとる値の大きさを考慮しながら必要最小限の値にすることが求められます。
例題2
例題1の問題に、「$x$の要素のうちちょうど3つが1になる」という制約を$\lambda=2$の重さで加えた関数f2(x)を定義し、例題1と同様の方法で解いてみましょう。
PyQUBOの利用
PyQUBOというPythonライブラリを用いると、QUBO問題の作成や制約の付与をより直感的に行うことができます。PyQUBOの利用スタイルの一番の特徴は、変数$x$の各要素$x_0, x_1,\cdots$をBinary型のオブジェクトとして扱い、それらに和・差・積といった演算を施すことによって関数$f(x)$を構成する点にあります。
例えば$f(x)=x_1 - 2x_2+3x_1x_2$という関数は次のように定義できます。
> import pyqubo
> x1 = pyqubo.Binary('x1')
> x2 = pyqubo.Binary('x2')
> f = x1 - 2 * x2 + 3 * x1 * x2
ただしこのfはまだ確定前のような状態で、compileメソッドを使ってModel型のオブジェクトに変換することで初めて関数値を評価できるようになります。例えば$$x_1=1, x_2=1$$を代入したときの値は次のようにして確認できます。
> Q = f.compile()
> Q.energy({'x1': 1, 'x2': 1}, 'BINARY')
2.0
変数の種類としては0/1の値を取るBinary型以外に、+1/-1の値をとるSpin型を使うこともできます（併用は不可）。詳しい利用方法は公式のドキュメントを参考にしてください。
PyQUBOドキュメント: https://pyqubo.readthedocs.io/en/latest/
例題3
課題1, 2でそれぞれ定義した関数f(x)とf2(x)を表現するPyQUBOオブジェクトを作成するコードを書いてください。制約条件を表現するにはConstraintクラスを利用することができます。
ソルバーの利用
問題のサイズが大きくなると、ここまでおこなってきた、すべての$x$に対してQUBOの目的関数を評価して確かめるというアプローチが現実的ではなくなってきます。探索空間が広かったり、探索に使える時間が限られたりしているような問題に対しては、近似的な手法で解を求めるための各種ソルバーを利用することになります。
イジングモデルとQUBOの等価性
統計物理の分野には、状態ベクトル変数$\sigma$の要素が$\pm 1$に値をもち、系のエネルギーがそれらの2次式として記述されるイジングモデルと呼ばれる物理モデルがあります。すなわち、系のエネルギーを記述するハミルトニアン$H(\sigma)$が、次のような形式でかかれます。
$$ H(\sigma) = \sum_i h_i\sigma_i + \sum_{i < j}^NJ_{i,j}\sigma_i\sigma_j $$
これは磁性体上の電子スピンの分布を説明するために考案されたモデルで、温度$T$において系の状態変数がボルツマン分布：
$$ p(\sigma) \propto \exp(-\frac{H(\sigma)}{k_BT})$$
に従うものとしてうまく記述できることが知られています。
式からは最も存在確率が高い状態が$H(\sigma)$を最小化するような状態$\sigma^\star = \mathrm{arg}\min_\sigma H(\sigma)$であることがわかります(ただし$\sigma^\star$の唯一性を仮定)。そして分布確率は$T\to 0$の極限ではその最適解のみに収束します。
$$p(\sigma)|_{T\to 0} = \begin{cases}
1 & (\mathrm{if}\ \sigma = \sigma^\star) \\
0 & (\mathrm{otherwise}) \\
\end{cases}$$
つまりこのような物理系を用意してゆっくりと温度を下げていけば、$H(\sigma)$を最小化するような状態$\sigma^\star$を求められるはずです。さらに$H(\sigma)$は先述のQUBOの目的関数$f(x)$と、変数が0/1をとるか$\pm 1$をとるかの違いを除いて一致しています。つまり適当な線形変換をかませることでQUBO問題を解くソルバーとして利用することもできるということです。
こうした架空の系を物理的に実現することは難しいですが、Markov Chain Monte Carlo (MCMC; マルコフ連鎖モンテカルロ)と呼ばれる確率過程で$T\to 0$に近づけていったときのボルツマン分布の収束をシミュレーションするSimulated Annealing (SA; 焼きなまし法)という手法があります。これまで扱ってきた問題をSAで解いてみましょう。
PythonでのSAの利用
SAの実装としては、D-Wave社が提供しているnealと呼ばれるPythonパッケージなどが利用しやすいでしょう。ここではこのSA実装を利用するとして、D-Wave 社が提供するQUBO問題の表現形式も合わせて解説します。
dimodおよびdwave-nealの二つのPythonライブラリをインストールします。
$ pip install dimod dwave-neal
前者はQUBOとイジングモデルを表現したり相互に変換できる画一的なPythonクラスBinaryQuadraticModel(略してBQM)などを定義しており、後者はそれを解くためのSAの実装です。インストールが完了したら、QUBO問題の係数をもとにBQMのインスタンスを作成してみましょう。
> from dimod import BinaryQuadraticModel as BQM
> Q = np.array([
[ 1, 2, -1, 0],
[ 0, -1, -2, -1],
[ 0, 0, -2, 2],
[ 0, 0, 0, 0]
])
> bqm = BQM.from_qubo(Q)
> bqm
BinaryQuadraticModel({0: 1.0, 1: -1.0, 2: -2.0, 3: 0.0}, {(0, 1): 2, (0, 2): -1, (1, 2): -2, (1, 3): -1, (2, 3): 2}, 0, 'BINARY')
次に、この問題をSAで解きます。neal.SimulatedAnnealingSamplerというクラスのインスタンスを作って、そのsampleメソッドにBQMのインスタンスを与えることで解(SampleSetクラス)が得られます。
> from neal import SimulatedAnnealingSampler
> sa_sampler = SimulatedAnnealingSampler()
> res = sa_sampler.sample(bqm)
> res
SampleSet(
rec.array( # res.recordでアクセス可能
[([0, 1, 1, 0], -5., 1)],
dtype=[ # レコードの型情報を記載
('sample', 'i1', (4,)),
('energy', '<f8'),
('num_occurrences', '<i8')
]
),
[0, 1, 2, 3], # res.variable: 変数の名前を表す。明示しなかったので、勝手に[0,1,2,3]と割り振られている
{ # res.info: アニーリングの詳細情報が格納される
'beta_range': [0.2772588722239781, 18.420680743952367],
'beta_schedule_type': 'geometric'
},
'BINARY' # res.vartype: 変数の種類。QUBO問題としてBQMを作ったのでBINARY
)
上記の出力にあるrec.arrayの内容は、SAを(1回)行った結果[0,1,1,0]というビットパターンがエネルギー-5を取り、それが1回出現したということを表しています。この解は例題1の解に一致しています。
さらに、BQMをPyQUBOからも作成して解いてみましょう。PyQUBOでコンパイルしたものから、さらにto_dimod_bqmメソッドを呼び出すとBQM型のオブジェクトが得られます。
> import pyqubo
> Q = np.array([
[ 1, 2, -1, 0],
[ 0, -1, -2, -1],
[ 0, 0, -2, 2],
[ 0, 0, 0, 0]
])
> x = pyqubo.Array.create('x', shape=(4), vartype='BINARY')
> xQx = x.dot(Q.dot(x))
> xQx_cons = xQx + 2 * pyqubo.Constraint((sum(x)-3)**2, label='sums to 3')
> bqm = xQx_cons.compile().to_dimod_bqm()
> sa_sampler = SimulatedAnnealingSampler()]
> res = sa_sampler.sample(bqm)
> res
SampleSet(rec.array([([0, 1, 1, 1], -4., 1)], ...), ...)
ここでは制約付きの問題を扱っていますが、きちんと解けていることがわかります。SimulatedAnnealingSamplerのsampleメソッドにはSAに関する様々なオプションがあります(ドキュメント)。問題の規模が大きくなって正しい解が得られないなどしたときは、これらのパラメータを調節すると改善するかもしれません。
D-Wave 量子アニーリングマシンの利用
D-Wave社は量子ビットを利用して、上述したイジングモデルに相当する系の物理的な実装を行いました。これを量子アニーリングマシンといい、先ほどのSAと同様にQUBO問題を解くことに使用できます。
量子アニーリングの直接的な利用
ここからはD-Waveのクラウドサービスに登録して発行される認証用のトークンが必要になります。D-WaveのマシンはWeb API経由で利用します。PythonでそのAPIをたたくライブラリも提供されているのでインストールしましょう。
$ pip install dwave-system
準備はこれだけです。実行例がこちらです。
> from dwave.system import DWaveSampler, AutoEmbeddingComposite
> dw = DWaveSampler(
endpoint = 'https://cloud.dwavesys.com/sapi',
token = 'SECRET TOKEN GOES HERE',
solver = 'DW_2000Q_6'
)
> dw_sampler = AutoEmbeddingComposite(dw)
> res = dw_sampler.sample(xQx_cons.compile().to_dimod_bqm())
> res
SampleSet(rec.array([([0, 1, 1, 1], -4., 1, 0.)], ... ), ...)
まずDWaveSamplerクラスのインスタンスを作成するときに、endpoint, token, solverという3つのオプションを指定しています。endpointはサーバーのURLで、自前でマシンを購入してホストしている場合を除いてこのままです。tokenは認証用のトークン文字列で、solverは用いるマシン名です。不定期のアップデートで使えるマシンが変わっていくので、その時使えるマシンをD-WaveのWebページ等で確認しておきましょう。
pipで一緒にインストールされる
dwave config createというコマンドを使ってこれらオプションのデフォルト値を設定しておくこともできます。ソースコードにトークンをハードコードするのは避けたほうがよいので、本格的に使用する場合はそちらで設定しておいてください。
名称からわかる通りDWaveSamplerのインスタンスが先ほどのSimulatedAnnealingSamplerのものに対応するのですが、ここではさらにAutoEmbeddingCompositeというクラスでラップしてから使用しています。このようにする理由は、D-Waveマシンが物理的な実装であるために、QUBOを指定する行列$Q$の任意の位置に要素を設定することができないという問題を隠蔽するためです[参考: D-Wave QPU Architecture: Topologies]。例えばキメラグラフという構造を持った量子アニーリングチップの上では、0番目と4番目のビット間の接続(カップラー)が存在するのでQ[0, 4]には0以外の値を設定しても良いですが、Q[0, 1]にはそれがないためできません。こういった制約を考慮しつつ問題を設計することは難しいので、隣接する複数のビットを仮想的な一つのビットとして扱うなど、グラフを勝手に変形して自動的にハードウェアにあった形に変換してくれるのがAutoEmbeddingCompositeの役割です。このような処理を「埋め込み(Embedding)」と呼び、デフォルトではminorminerと呼ばれるヒューリスティックなアルゴリズムが利用されます。
実行時間に関する注意
D-Waveのマシンには利用時間(QPU time)の上限が設定されていて、それを使い果たすと次に利用時間が補充されるまで使用できなくなってしまいます。D-Wave Leapのページで適宜時間の残高を確認するようにしましょう。
D-Waveマシンでサンプリングをするときに各種パラメータを設定することができますが、最も基本的な項目としてはnum_readsとannealing_timeというパラメータが挙げられるでしょう。sampleメソッドを呼び出すときに指定します。num_readsとはひとつのQUBO問題を投げたときに回路の上でアニーリングと解の読み出しを行う回数（デフォルトは1）で、annealing_timeは一回のアニーリングにかける時間[μs] (デフォルトは20)を指定します。
問題を投げて消費される時間の目安としては、
- 問題を回路に読み込む時間：約10000μs
- アニーリングにかかる時間：(annealing_time + 読み出し等のオーバーヘッド≒300μs) * num_reads
の合計で考えてください。annealing_timeは長くしすぎても熱的なノイズが混入して解の精度が落ちるだけなのでそこまで大きくすることはあまりないですが、より良い解を得るためにnum_readsを増やしたくなることがあるかと思います。式からわかる通り、num_readsを30くらいにするとアニーリングの時間が大体問題の読み込みにかかる時間10000μs = 10msと同じ程度になります。より簡略化して消費時間は
10 * (1 + num_reads / 30) [ms]
程度と概算できるでしょう。気を付けたほうがいいのは、num_readsを数百とかに設定してしまうと、意外とすぐに時間を使いきってしまうということです。サンプリング用途など特別の必要がない限りは、num_readsは1~30, 高々50くらいにとどめて利用するのがよいかと思います。
仮に1問のサンプリングが30[ms]消費する場合、月1時間の割り当てで3600 [s] / 30 [ms] = 120000 問の問題を解くことができることになります。利用時の参考にしてください。
時間に関する詳細はこちらにあります。
https://docs.dwavesys.com/docs/latest/timing_cloud_client_use.html
ハイブリッドソルバーの利用
量子アニーリングマシンは扱える問題の規模が埋め込みを利用しない場合であっても2000や5000ビット程度に限られています。D-Wave Hybrid Solversの枠組みでは、商業規模の問題への利用を促す目的もあってか、変数の多いQUBO問題を分割してそれぞれをD-Waveマシンで解いて結果を古典的なアルゴリズムで統合するという「ハイブリッド」なアプローチがとられています。
多層的で中身が不透明なきらいはありますが、性能が向上することは確認されています。
利用にはdwave-hybridというPythonパッケージをpipでインストールしてください。分割数の設定など自由度が高くなっていますので、利用法はドキュメントをご参照ください。githubのプロジェクトページにあるサンプルコードだけ、さらっと眺めておきます。
# https://github.com/dwavesystems/dwave-hybrid より引用
import dimod
import hybrid
# Construct a problem
bqm = dimod.BinaryQuadraticModel({}, {'ab': 1, 'bc': -1, 'ca': 1}, 0, dimod.SPIN)
# Define the workflow
iteration = hybrid.RacingBranches(
hybrid.InterruptableTabuSampler(),
hybrid.EnergyImpactDecomposer(size=2)
| hybrid.QPUSubproblemAutoEmbeddingSampler()
| hybrid.SplatComposer()
) | hybrid.ArgMin()
workflow = hybrid.LoopUntilNoImprovement(iteration, convergence=3)
# Solve the problem
init_state = hybrid.State.from_problem(bqm)
final_state = workflow.run(init_state).result()
# Print results
print("Solution: sample={.samples.first}".format(final_state))
この例ではhybrid.RacingBranchesでTabu探索という手法によるSamplerと(問題分割→量子アニーリングによるそれぞれの解の探索→結果の統合)によるSamplerとを競わせて、より良い解を出した方をhybrid.ArgMinで採用するという手順(workflow)を変数iterationに格納しています。これを収束するまで何度も適用するという内容を変数workflowに格納し、問題bqmを解かせています。
ただしbqmをそのまま渡すのではなく、Stateというクラスに変換して内部状態を持たせているあたりなども、これまでのSimulatedAnnealingSamplerやDWaveSamplerといったソルバーとは大きく異なる点です。
最後に
以上の通り、QUBO問題を作成してソルバーを利用して解くまでの手順を説明しました。各種ツールの開発者に感謝いたします。