D-Waveの制約式をイジングモデルに変換するライブラリの仕組みを追ってみた

# dwavebinarycsp

量子アニーリングマシンD-Waveの [dwavebinarycsp](https://github.com/dwavesystems/dwavebinarycsp) というライブラリを使うと、制約式充足問題を二次形式最小化問題すなちイジングモデル（もしくはQUBO）に変換できます。たとえば、

```math
a, b, c \in \{-1, 1\} \\
abc = 1 \\
a+b+c < 3
```

といった制約式充足問題があったとするなら、

```python
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)
```

```python
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$ 、と表すことにすると

```math
\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` した後に制約式をどう保持しているか覗いてみると、

```python
print(csp.constraints)
```

```python
[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` にラムダ式を食わせると、下記の関数が呼ばれて解集合に変換されています。

```python
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](https://github.com/dwavesystems/dwavebinarycsp/blob/19068bdf8e9941a51b4d35911b69c0fe2eed836e/dwavebinarycsp/core/constraint.py)

ということで、次は同じ解集合を持つ最小化問題を見つければいいという話になって、それをしている入り口の関数が `stitch` です。抜粋します。

```python
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](https://github.com/dwavesystems/dwavebinarycsp/blob/f30f2f0d69c7794670930b8b9b415921221fe431/dwavebinarycsp/compilers/stitcher.py)

すなわち。制約式ごとに `get_penalty_model` で最小化問題への変換を試みて、解を満たすときと満たさないときの最低エネルギーギャップが閾値を超える最小化問題が見つかるまで、補助変数を追加します。そして、それらを足し合わせたものが答というわけです。ぱちぱち。

# penaltymodel

それで肝心の `get_penalty_model` の実装ですが、これは [penaltymodel](https://github.com/dwavesystems/penaltymodel) という別のD-Waveライブラリを呼びだしています。
penaltymodel はちょっとややこしくて、インストールパッケージがコアの `penaltymodel` のほか `penaltymodel-cache` `penaltymodel-maxgap` `penaltymodel-mip` と用意されていて、何をインストールするかによって `get_penalty_model` の実装を切り替えられるようになっています。

どういうことかと `penaltymodel-cache` を覗いてみると、 `penaltymodel_factory(100)` と一番高い優先度が振られていて、これはまぁキャッシュから引いているだけですね。
キャッシュヒットしなかったら `MissingPenaltyModel` という例外が投げられて、もっと優先度の低いものが呼ばれます。

```python
@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](https://github.com/dwavesystems/penaltymodel/blob/79533236e4c389297c972a5113c99400de2ea580/penaltymodel_cache/penaltymodel/cache/interface.py)

つまり、たとえばこの `penaltymodel-maxgap` です。

```python
@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](https://github.com/dwavesystems/penaltymodel/blob/002f2c1117be03263d06ed93d8c8326f6ffab484/penaltymodel_maxgap/penaltymodel/maxgap/interface.py)

さらにこの実装を追っていくと、コメントで論文参照となっています。

```python
"""
.. [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](https://github.com/dwavesystems/penaltymodel/blob/5a07feaae1aff9bcc762fd2684068a2e245a017b/penaltymodel_maxgap/penaltymodel/maxgap/generation.py)

……なにかしらの最適化問題を解いているということのようですね。
量子アニーリングマシンは組合せ最適化問題が得意ではあるのですが、量子アニーリングマシンに解かせるための過程で古典マシンで最適化問題を解かなければならないことは、量子アニーリングあるあるみたいです。人生そんなもんです。それでも全体的には計算効率が良くなることを祈りましょう。

さて `penaltymodel-maxgap` がインストールされていない場合、優先度最低な `penaltymodel-mip` の `get_penalty_model` が使われます。

```python
@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](https://github.com/dwavesystems/penaltymodel/blob/176ed2e155d364ed0672f7b5fda70a2503aaf944/penaltymodel_mip/penaltymodel/mip/interface.py)

こちらの実装は [Google Optimization Tools](https://developers.google.com/optimization/) を使って、最適化問題を解いているようです。

# dimod

さて、手に入れたイジングモデルが正しいかどうか確認してみましょう。
いざD-Waveの実機に投げ……てもいいのですが、問題サイズが小さいので古典的なソルバーで見てしまいましょう。
やはりD-Waveのライブラリ [dimod](https://github.com/dwavesystems/dimod) を使うと、

```python
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)
```

```python
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` を小さく指定してやると、

```python
bqm = dwavebinarycsp.stitch(csp, min_classical_gap=.1)

print(bqm)
```

```python
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]` となっているようです。

```python
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](https://github.com/dwavesystems/penaltymodel/blob/b1951fbad8faa528c25a7099bfaa2668d02e3ab2/penaltymodel_core/penaltymodel/core/classes/specification.py)

# 実行環境

```sh
$ 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     
```