パズルゲームRulloを量子アニーリングマシンで解いてみる

2019.04.07 説明図を修正。

# 概要
[Rullo](https://appget.com/appli/view/70100/)は、マス内に書かれた数字の合計が指定された値（目標値）に等しくなるように、セルの取捨選択を行うゲームです。

![explain.png](https://qiita-image-store.s3.amazonaws.com/0/152794/f6904325-e296-0757-9089-248ca2fd1765.png)

図の例では、列に上から順に3,8,5,1と並んでおり、目標値には9が指定されています。8と1を選択状態にし、3と5を非選択状態にすれば目標値に等しくなります。上記リンク先のゲーム紹介を読んでいただくのが早いと思いますが、実際には各行・各列に対して目標値が定められているため、上記例のように特定の列（または行）だけで一意に選択状態とすべきマスが決まりません。そこで、このゲームを選択状態/非選択状態を最適化する組み合わせ最適化問題として捉え直し、量子アニーリングの手法を適用して解いてみたいと思います。

# 解き方
非常にシンプルなモデル化で簡単に解くことができます。~~まるで量子アニーリングで解くために生まれてきたようなパズルです。~~

各セルにインデックスを振ってみます。行（縦方向）には$i$を、列（横方向）には$j$を割り当てて、セルに書かれた数字を$V_{ij}$で表現します。また、各セルの選択状態/非選択状態を表すバイナリ値を$q_{ij}$とします。すなわち、この$q_{ij}$が最適化対象となりますので、QUBO変数に該当します。さらに、下図（4x4の例）のように目標値$R_{i}$および$C_{j}$を定義します。

![problem.png](https://qiita-image-store.s3.amazonaws.com/0/152794/b0cd35e0-1928-337d-d61f-ce59aab319de.png)

行$I$および列$J$について、

```math
\sum _{j} V_{Ij}q_{Ij} = R_{I}, \quad \sum _{i} V_{iJ}q_{iJ} = C_{J}
```
を満たすように制約条件を決めてあげれば問題を解くことができます。

# モデル
組み合わせ最適化問題として、QUBOに落とし込みます。先程の式を満たすようにハミルトニアンを構成すると、

```math
H = \sum _{i} \left( \sum _{j} V_{ij} q_{ij} - R_{i}\right) ^{2} + \sum _{j} \left( \sum _{i} V_{ij} q_{ij} - C_{j} \right) ^{2}, \quad q_{ij} \in \{ 0, 1\}
```
となります。上記ハミルトニアンの最小値は0ですので、計算結果が0になっていれば正しい解が得られたことを確認することができます。

## ハミルトニアンの展開
展開してQUBOパラメータを決めます。計算の途中で、$q_{ij}^{2} = q_{ij}$であることを使います。

```math
\begin{align}
H &= \sum _{i} \left( \sum _{j} V_{ij} q_{ij} - R_{i}\right) ^{2} + \sum _{j} \left( \sum _{i} V_{ij} q_{ij} - C_{j} \right) ^{2} \\
&= \sum _{i} \left[ \left( \sum _{j} V_{ij} q_{ij} \right) ^{2} - 2R_{i} \sum _{j} V_{ij} q_{ij} + R_{i} ^{2} \right] 
+ \sum _{j} \left[ \left( \sum _{i} V_{ij} q_{ij} \right) ^{2} - 2C_{j} \sum _{j} V_{ij} q_{ij} + C_{j} ^{2} \right]  \\
&=  2\sum _{i}\sum _{j} \left( V_{ij} ^{2} - R_{i} - C_{j} \right) q_{ij} + 2\sum _{i} \sum _{j < j^{\prime}} V_{ij} V_{ij^{\prime}} q_{ij} q_{ij^{\prime}}
+ 2\sum _{j} \sum _{i < i^{\prime}} V_{ij} V_{i^{\prime}j} q_{ij} q_{i^{\prime}j} + \sum _{i} R_{i} ^{2} + \sum _{j} C_{j} ^{2}. 
\end{align}
```
定数項とハミルトニアン全体の定数倍は無視してよいので、結果は以下のようになります。

```math
H = \sum _{i}\sum _{j} \left( V_{ij} ^{2} - R_{i} - C_{j} \right) q_{ij} + \sum _{i} \sum _{j < j^{\prime}} V_{ij} V_{ij^{\prime}} q_{ij} q_{ij^{\prime}}
+ \sum _{j} \sum _{i < i^{\prime}} V_{ij} V_{i^{\prime}j} q_{ij} q_{i^{\prime}j}. 
```

# 結果
4x4の問題に上記QUBOを適用し、Blueqatで実装して解いてみました。問題（マスの値と目標値）は、以下のような配列で用意します。

```python
import numpy as np

TBL_ROW = np.array([10, 20, 7, 6]) # 目標値（R）
TBL_COL = np.array([11, 9, 15, 8]) # 目標値（C）
TBL_CELL = np.array([[6, 3, 4, 2], [4, 8, 4, 8], [1, 5, 7, 1], [5, 1, 6, 3]]) # マスの値
```

実装例は以下の通りです。

```python
import blueqat
import numpy as np

class RulloSolver:
    # コンストラクタ
    def __init__(self, _size):
        self.size = _size   # the size of problem: size x size

    # 結果解析用
    def _checkResult(self, _result):
        ret = 0
        for i in range(self.size):
            sum = 0
            for j in range(self.size):
                sum += _result[i][j] * self.area[i][j]
            ret += (sum - self.row[i])**2
        for j in range(self.size):
            sum = 0
            for i in range(self.size):
                sum += _result[i][j] * self.area[i][j]
            ret += (sum - self.column[j])**2
        return (ret == 0)

    # 問題を入力する
    def setProblem(self, _area, _row, _column):
        self.area = _area
        self.row = _row
        self.column = _column

    # 問題を（Blueqatで）解く
    def solve(self, _iteration):
        a = blueqat.opt.opt()
        Q = np.array([[0.0 for i in range(self.size * self.size)] for j in range(self.size * self.size)])
        
        # h
        for i in range(self.size):
            for j in range(self.size):
                Q[self.size * i + j][self.size * i + j] = self.area[i][j] * (self.area[i][j] - self.row[i] - self.column[j])

        # J
        for i in range(self.size):
            for j1 in range(self.size-1):
                for j2 in range(j1+1, self.size):
                    Q[self.size * i + j1][self.size * i + j2] += self.area[i][j1] * self.area[i][j2]

        for j in range(self.size):
            for i1 in range(self.size-1):
                for i2 in range(i1+1, self.size):
                    Q[self.size * i1 + j][self.size * i2 + j] += self.area[i1][j] * self.area[i2][j]

        a.qubo = Q
        responses = a.sa(shots=_iteration, sampler="fast")
        for response in responses:
            ret = np.reshape(response, (self.size, self.size))
            validity = self._checkResult(ret)
            if validity != False:
                return ret
        return []
       
if __name__ == '__main__':
    solver = RulloSolver(4)

    # 問題（マスの値と目標値）
    TBL_ROW = np.array([10, 20, 7, 6]) # R
    TBL_COL = np.array([11, 9, 15, 8]) # C
    TBL_CELL = np.array([[6, 3, 4, 2], [4, 8, 4, 8], [1, 5, 7, 1], [5, 1, 6, 3]])

    solver.setProblem(TBL_CELL, TBL_ROW, TBL_COL)
    ret = solver.solve(30)
    print("ret={}".format(ret))

```

結果は以下を得ました。正しい結果になっています。

```console
ret=[[1 0 1 0]
 [0 1 1 1]
 [0 0 1 0]
 [1 1 0 0]]
```

# まとめ
Rulloを量子アニーリングの手法で解く方法を紹介しました。選択状態を0/1で表現でき、また複雑な制約条件がないため、シンプルなモデルで解くことができます。初めて量子アニーリングに取り組む方には良い題材なのではないかと思います。