dwavebinarycsp
量子アニーリングマシンD-Waveの dwavebinarycsp というライブラリを使うと、制約式充足問題を二次形式最小化問題すなちイジングモデル（もしくはQUBO）に変換できます。たとえば、
a, b, c \in \{-1, 1\} \\
abc = 1 \\
a+b+c < 3
といった制約式充足問題があったとするなら、
import dwavebinarycsp
csp = dwavebinarycsp.ConstraintSatisfactionProblem(dwavebinarycsp.SPIN)
csp.add_constraint(lambda a, b, c: a * b * c == 1, ['a', 'b', 'c'])
csp.add_constraint(lambda a, b, c: a + b + c < 3, ['a', 'b', 'c'])
bqm = dwavebinarycsp.stitch(csp)
print(bqm)
BinaryQuadraticModel({'a': -1.0, 'aux0': -2.0, 'aux1': 0.0, 'b': 2.0, 'c': 1.0, 'aux2': -1.0}, {('aux0', 'a'): 1.0, ('aux1', 'a'): 1.0, ('aux1', 'aux0'): -1.0, ('b', 'a'): -1.0, ('b', 'aux0'): -1.0, ('b', 'aux1'): -1.0, ('c', 'a'): 0.0, ('c', 'aux0'): -1.0, ('c', 'aux1'): -1.0, ('c', 'b'): 1.0, ('aux2', 'a'): 1.0, ('b', 'aux2'): -1.0, ('c', 'aux2'): 1.0}, 0.0, Vartype.SPIN)
となるので、補助変数 aux0 を $x$ 、 aux1 を $y$ 、 aux2 を $z$ 、と表すことにすると
\newcommand{\argmin}{\mathop{\rm arg~min}\limits}
\argmin_{a, b, c \in \{-1, 1\}} -a+2b+c-2x-z-ab+bc+ax+ay+az-bx-by-bz-cx-cy+cz-xy
というイジングモデルに変換できることが分かります。やりましたね。
さて、このライブラリはどうやってこの変換を計算しているのでしょうか。
注目すべきは add_constraint に任意のラムダ式を食わせられるところです。ということは制約式を二乗してN体を分解するといった、代数的な変換操作をしているわけではなさそうですね？
ということで add_constraint した後に制約式をどう保持しているか覗いてみると、
print(csp.constraints)
[Constraint.from_configurations(frozenset({(1, -1, -1), (1, 1, 1), (-1, -1, 1), (-1, 1, -1)}), ('a', 'b', 'c'), Vartype.SPIN, name='Constraint'), Constraint.from_configurations(frozenset({(1, 1, -1), (-1, 1, 1), (1, -1, -1), (-1, -1, 1), (1, -1, 1), (-1, 1, -1), (-1, -1, -1)}), ('a', 'b', 'c'), Vartype.SPIN, name='Constraint')]
と食わせたラムダ関数は消失していて、代わりに解集合が保持されていることが見てとれますね。
実際 add_constraint にラムダ式を食わせると、下記の関数が呼ばれて解集合に変換されています。
class Constraint(Sized):
...
def from_func(cls, func, variables, vartype, name=None):
...
configurations = frozenset(config
for config in itertools.product(vartype.value, repeat=len(variables))
if func(*config))
dwavebinarycsp/dwavebinarycsp/core/constraint.py
ということで、次は同じ解集合を持つ最小化問題を見つければいいという話になって、それをしている入り口の関数が stitch です。抜粋します。
def stitch(csp, min_classical_gap=2.0, max_graph_size=8):
...
bqm = dimod.BinaryQuadraticModel.empty(csp.vartype)
...
for const in csp.constraints:
...
pmodel = None
...
for G in iter_complete_graphs(const.variables, max_graph_size + 1, aux):
# construct a specification
spec = pm.Specification(
graph=G,
decision_variables=const.variables,
feasible_configurations=configurations,
vartype=csp.vartype
)
# try to use the penaltymodel ecosystem
try:
pmodel = pm.get_penalty_model(spec)
except pm.ImpossiblePenaltyModel:
# hopefully adding more variables will make it possible
continue
if pmodel.classical_gap >= min_classical_gap:
break
...
bqm.update(pmodel.model)
return bqm
dwavebinarycsp/dwavebinarycsp/compilers/stitcher.py
すなわち。制約式ごとに get_penalty_model で最小化問題への変換を試みて、解を満たすときと満たさないときの最低エネルギーギャップが閾値を超える最小化問題が見つかるまで、補助変数を追加します。そして、それらを足し合わせたものが答というわけです。ぱちぱち。
penaltymodel
それで肝心の get_penalty_model の実装ですが、これは penaltymodel という別のD-Waveライブラリを呼びだしています。
penaltymodel はちょっとややこしくて、インストールパッケージがコアの penaltymodel のほか penaltymodel-cache penaltymodel-maxgap penaltymodel-mip と用意されていて、何をインストールするかによって get_penalty_model の実装を切り替えられるようになっています。
どういうことかと penaltymodel-cache を覗いてみると、 penaltymodel_factory(100) と一番高い優先度が振られていて、これはまぁキャッシュから引いているだけですね。
キャッシュヒットしなかったら MissingPenaltyModel という例外が投げられて、もっと優先度の低いものが呼ばれます。
@pm.interface.penaltymodel_factory(100)
def get_penalty_model(specification, database=None):
"""Factory function for penaltymodel_cache.
...
Raises:
:class:`penaltymodel.MissingPenaltyModel`: If the penalty model is not in the
cache.
penaltymodel/penaltymodel_cache/penaltymodel/cache/interface.py
つまり、たとえばこの penaltymodel-maxgap です。
@pm.penaltymodel_factory(-100) # set the priority to low
def get_penalty_model(specification):
"""Factory function for penaltymodel_maxgap.
...
Raises:
:class:`penaltymodel.ImpossiblePenaltyModel`: If the penalty cannot be built.
penaltymodel/penaltymodel_maxgap/penaltymodel/maxgap/interface.py
さらにこの実装を追っていくと、コメントで論文参照となっています。
"""
.. [DO] Bian et al., "Discrete optimization using quantum annealing on sparse Ising models",
https://www.frontiersin.org/articles/10.3389/fphy.2014.00056/full
.. [MC] Z. Bian, F. Chudak, R. Israel, B. Lackey, W. G. Macready, and A. Roy
"Mapping constrained optimization problems to quantum annealing with application to fault diagnosis"
https://arxiv.org/pdf/1603.03111.pdf
"""
penaltymodel/penaltymodel_maxgap/penaltymodel/maxgap/generation.py
……なにかしらの最適化問題を解いているということのようですね。
量子アニーリングマシンは組合せ最適化問題が得意ではあるのですが、量子アニーリングマシンに解かせるための過程で古典マシンで最適化問題を解かなければならないことは、量子アニーリングあるあるみたいです。人生そんなもんです。それでも全体的には計算効率が良くなることを祈りましょう。
さて penaltymodel-maxgap がインストールされていない場合、優先度最低な penaltymodel-mip の get_penalty_model が使われます。
@pm.penaltymodel_factory(-150)
def get_penalty_model(specification):
"""Factory function for penaltymodel-mip.
...
Raises:
:class:`penaltymodel.ImpossiblePenaltyModel`: If the penalty cannot be built.
penaltymodel/penaltymodel_mip/penaltymodel/mip/interface.py
こちらの実装は Google Optimization Tools を使って、最適化問題を解いているようです。
dimod
さて、手に入れたイジングモデルが正しいかどうか確認してみましょう。
いざD-Waveの実機に投げ……てもいいのですが、問題サイズが小さいので古典的なソルバーで見てしまいましょう。
やはりD-Waveのライブラリ dimod を使うと、
import dimod
import itertools
res = dimod.ExactSolver().sample(bqm)
min_energy = next(res.data(sorted_by='energy')).energy
for sample in itertools.takewhile(lambda sample: sample.energy == min_energy, res.data(sorted_by='energy')):
print(csp.check(sample.sample), sample.energy, sample.sample)
True -7.0 {'a': 1, 'aux0': -1, 'aux1': -1, 'aux2': -1, 'b': -1, 'c': -1}
True -7.0 {'a': 1, 'aux0': -1, 'aux1': -1, 'aux2': 1, 'b': -1, 'c': -1}
True -7.0 {'a': -1, 'aux0': 1, 'aux1': 1, 'aux2': 1, 'b': -1, 'c': 1}
True -7.0 {'a': -1, 'aux0': 1, 'aux1': 1, 'aux2': -1, 'b': -1, 'c': 1}
True -7.0 {'a': -1, 'aux0': 1, 'aux1': 1, 'aux2': 1, 'b': 1, 'c': -1}
と制約式を満たす解が５つ見つかりました。
あれ３つなはずでは？と思うかもしれませんが、それは補助変数のバリエーションが出てしまっているからで、 (a,b,c) のパターンは３つですね。
エネルギーギャップ
ところでこの問題に補助変数３つもいるの？と思ったりするわけですが、実は stich を呼ぶときに min_classical_gap を小さく指定してやると、
bqm = dwavebinarycsp.stitch(csp, min_classical_gap=.1)
print(bqm)
BinaryQuadraticModel({'a': -0.5, 'aux0': -1.0, 'b': 1.5, 'c': 0.5, 'aux1': -1.0}, {('aux0', 'a'): 1.0, ('b', 'a'): -0.5, ('b', 'aux0'): -1.0, ('c', 'a'): 0.5, ('c', 'aux0'): -1.0, ('c', 'b'): 0.5, ('aux1', 'a'): 1.0, ('b', 'aux1'): -1.0, ('c', 'aux1'): 1.0}, 0.0, Vartype.SPIN)
補助変数が２つで済むようになりました。
ただし最低エネルギーギャップが低いということは、おそらく最適解以外に落ちる可能性が上がってしまうのと、制約付き最小化問題なんかを扱うとき制約式を破りやすくなってしまうので注意が必要です。
そしてエネルギーギャップの大きさなんて二次形式の係数次第なんじゃね？という気もしますが、係数の範囲は ising_linear_ranges ising_quadratic_ranges で指定することができて、デフォルト値は [-2, 2] [-1, 1] となっているようです。
class Specification(object):
...
def __init__(self, graph, decision_variables, feasible_configurations, vartype,
ising_linear_ranges=None,
ising_quadratic_ranges=None):
...
linear_ranges[v] = [-2, 2]
...
quad_ranges[u][v] = quad_ranges[v][u] = [-1, 1]
penaltymodel/penaltymodel_core/penaltymodel/core/classes/specification.py
実行環境
$ python --version
Python 3.7.0
$ pip list | grep -e dwavebinarycsp -e penaltymodel -e dimod
dimod 0.7.9
dwavebinarycsp 0.0.9
penaltymodel 0.15.3
penaltymodel-cache 0.3.2
penaltymodel-maxgap 0.4.1
penaltymodel-mip 0.1.3