D-waveの量子アニーリングのクラウドサービスLeapで数独を解いてみる

はじめに

AIREV株式会社の森岡です。この記事では量子アニーリングとは何か、どのように使っていくのか、について解説していきたいと思います。また、D-Wave社が提供する量子アニーリングのクラウドサービスである Leap とそのSDKである Ocean を使って実際に数独を量子アニーリングで解いていきます。量子アニーリングの細かい理論についてはこの記事では扱いません。量子アニーリングを使って例題を解いてみる、という点にフォーカスした記事となりますこと、ご留意ください。

量子アニーリングとは

量子アニーリングとは、量子ゆらぎという現象を用いた、組合せ最適化問題の解法のことを指します。組合せ最適化問題とは、膨大な候補のなかからある制約のもとで最も良い解を選ぶ問題で、身近な例では出発地点から目的地点までの経路探索や、決められたサイズのリュックにおやつを詰め込んでおやつの金額を最大にするナップザック問題等が挙げられます。

組合せ最適化問題を解く方法については古くから研究されており、局所探索法という、今いる解の近くを探索してより良い解へと更新していく枠組みがあります。局所探索法では、その周辺では良い解のように見えるけれど全体で見ると最適ではない解、局所最適解が求まることがありますが、局所解にはまらずに最適解を探索しようとするアルゴリズムがいくつか研究されています。そのうちの一つに焼きなまし法という、金属材料を高温で熱した後に徐々に冷やすという再結晶化を行うことで材料組織を均質化し欠陥を減らしていく、という金属工学の焼きなましを模倣したアルゴリズムがあります。焼きなまし法では温度というパラメータを導入して、最初は高温に設定して解を探索し、徐々に温度を低くしていくことで局所解を脱し最適解へと遷移していきます。東工大の門脇、西森は、温度による遷移に代わり、量子ゆらぎによる遷移を導入して解を探索するアルゴリズム、量子アニーリングを1998年に "Quantum annealing in the transverse Ising model" という論文で提案しました。特定の条件下で、量子アニーリングは従来の焼きなまし法に比べてより良い解が求まる確率が高い、という結果が得られています。

このアルゴリズムをハードウェア実装したのがD-Wave社の量子アニーリングマシンになります。2011年に超電導量子ビットを用いた量子アニーリングマシンを開発し、2018年にはクラウドで量子アニーリングマシンにアクセスすることができる Leap の提供を開始しました。

Leap、Ocean SDKについて

Leap とは、D-Wave社が開発している量子アニーリングマシンにアクセスして実際に問題を解くことができる量子クラウドサービスで、D-Waveのマシンを使用するためのSDKが Ocean になります。
Ocean SDKには量子アニーリングのほか、タブー法や焼きなまし法など、古典コンピュータ上で動作する組合せ最適化のアルゴリズムも実装されています。

Leap、Ocean SDKを使って数独を解いてみる

数独とは $N^{2} \times N^{2}$ （$N$ は整数）正方形の各マスに $1 \sim N^{2}$ の数字を入れるパズルです。ここでは、 $N = 2$ 、すなわち $4 \times 4$ の数独の問題を量子アニーリングを使って解いていきたいと思います。

数独のルール

まず、数独のルール（制約）についておさらいしておきます。

1. $N^{2} \times N^{2}$ の正方形の各マスには $1 \sim N^{2}$ のいずれかの値が入ります
2. 各行には、同じ数字が重複して入ってはいけません
3. 各列には、同じ数字が重複して入ってはいけません
4. $N \times N$ のブロックに分割した時に、各ブロックには同じ数字が重複して入ってはいけません

数独の問題を定式化

量子アニーリングで解けるようするには、数独の問題を組合せ最適化の一つである二次制約なし二値最適化（QUBO; Quadratic Unconstrained Binary Optimization）の形に落とし込む必要があります。QUBOでは、$0$と$1$の2値をとる$M$個の変数 $x_{1}, x_{2}, \cdots, x_{M}$ に対し、2次形式の評価関数 $f: \left\{ 0, 1\right\}^{M} \rightarrow \mathbb{R}$ を最小にするような変数の組合せ $\mathbf{x} = \left(x_{1}, x_{2}, \cdots, x_{M} \right)$ を求めていきます。

$$\begin{aligned} f(\mathbf{x}) & = \mathbf{x}^{\mathrm{T}} \mathbf{Q} \mathbf{x} \\ & = \sum_{i=1}^{M} \sum_{j=i}^{M} q_{ij} x_{i} x_{j} \end{aligned}$$

ただし、 $\mathbf{Q} \in \mathbb{R}^{M \times M}$ は上三角行列で、対角成分は成分$i$の重み、非対角成分は成分$i$と成分$j$の相互関係を表しています。

ここで、数独の $i$ 行目、 $j$ 列目のマスの値が $k$ であるかどうか $x_{ijk}$ が、QUBOの2値変数に対応します。上の例では、$1$行目$2$列目のセルの値は$1$なので、$x_{121} = 1$ 、$2$行目$1$列目のセルの値は$3$ではないので、$x_{213} = 0$ になります。
最終的に求めたいのは、 $N = 2$ のときの数独に対し、どのマスにどの値が入るかの組合せ $\mathbf{x} = \left(x_{111}, x_{112}, \cdots, x_{114}, x_{121}, \cdots, x_{124}, \cdots, x_{ij1}, \cdots, x_{ijk}, \cdots, x_{ijN^{2}}, \cdots, x_{441}, \cdots, x_{444} \right)$ になります。
これをQUBOの枠組みで求めるには、マスに入る値の組合せ $\mathbf{x}$ に対して、数独の条件を満たさなかったときに大きい値（反則点）をとり、全ての条件を満たした時に $0$ になるような評価関数 $f$ を自分で設計する必要があります。上記で挙げた数独のそれぞれのルールを評価関数に落とし込んでいきましょう。

1. $N^{2} \times N^{2}$ の正方形の各マスには $1 \sim N^{2}$ のいずれかの値が入る

$i$ 行目、 $j$ 列目のマスには $1 \sim N^{2}$ のいずれか1つの値が入り、2つ以上の値が同時に入ることはありません。このルールを次のような式で表現することができます。

$$\left( \sum_{k=1}^{N^{2}} x_{ijk} - 1 \right)^{2}$$

数独のルールを満たすとき、すなわち $x_{ij1}, x_{ij2}, \cdots, x_{ijN^{2}}$ のいずれか1つの値が $1$ で、他が $0$ になるとき、上式の評価値は $0$ となります。数独のルールを満たさないとき、すなわち同行に値の重複が $m (\geq 2)$ 個存在する場合、 $x_{ij1}, x_{ij2}, \cdots, x_{ijN^{2}}$ のうち $m$ 個の値が $1$ で、残りが $0$ となるので、上式の評価値は $(m - 1)^{2}$ となります。また、 $i$ 行目 $j$ 列目のマスにはどんな値も入らない場合は $x_{ij1} = x_{ij2} = \cdots = x_{ijN^{2}} = 0$ となり、上式の評価値は $1$ となります。

全てのマスに対して同様の議論ができるので、「$N^{2} \times N^{2}$ の正方形の各マスには $1 \sim N^{2}$ のいずれかの値が入る」というルールは以下の関数 $f_{1}$ で表現することができます。

$$f_{1}\left( \mathbf{x} \right) = \sum_{i=1}^{N^{2}} \sum_{j=1}^{N^{2}} \left( \sum_{k=1}^{N^{2}} x_{ijk} - 1 \right)^{2}$$

2. 各行には、同じ数字が重複して入ってはいけない

$i$ 行目にあるマスに $k$ という値は1回しかとることができない、というルールを次のような式で表現することができます。

$$\left( \sum_{j=1}^{N^{2}} x_{ijk} - 1 \right)^{2}$$

数独のルールを満たすとき、すなわち $x_{i1k}, x_{i2k}, \cdots, x_{iN^{2}k}$ のいずれか1つの値が $1$ で、他が $0$ になるとき、上式の評価値は $0$ となります。数独のルールを満たさないとき、すなわち同行に値の重複が $m (\geq 2)$ 個存在する場合、 $x_{i1k}, x_{i2k}, \cdots, x_{iN^{2}k}$ のうち $m$ 個の値が $1$ で、残りが $0$ となるので、上式の評価値は $(m - 1)^{2}$ となります。また、 $i$ 行目にあるどのマスも $k$ という値を取らない場合は $x_{i1k} = x_{i2k} = \cdots = x_{iN^{2}k} = 0$ となり、上式の評価値は $1$ となります。以下の図は、2行目に適当な数値を入れたときの評価値の例になります。

全ての取り得る行、値に対しても同様に議論できるので、「各行には、同じ数字が重複して入ってはいけない」というルールは以下の関数 $f_{2}$ で表現することができます。

$$f_{2}\left( \mathbf{x} \right) = \sum_{i=1}^{N^{2}} \sum_{k=1}^{N^{2}} \left( \sum_{j=1}^{N^{2}} x_{ijk} -1 \right)^{2}$$

3. 各列には、同じ数字が重複して入ってはいけない

「各行には、同じ数字が重複して入ってはいけない」と同様に、$j$ 列目にあるマスに $k$ という値は1回しかとることができない、というルールを次のような式で表現することができます。

$$\left( \sum_{i=1}^{N^{2}} x_{ijk} - 1 \right)^{2}$$

数独のルールを満たすとき、すなわち $x_{1jk}, x_{2jk}, \cdots, x_{N^{2}jk}$ のいずれか1つの値が $1$ で、他が $0$ になるとき、上式の評価値は $0$ となります。数独のルールを満たさないとき、すなわち同列に値の重複が $m (\geq 2)$ 個存在する場合、 $x_{1jk}, x_{2jk}, \cdots, x_{N^{2}jk}$ のうち $m$ 個の値が $1$ で、残りが $0$ となるので、上式の評価値は $(m - 1)^{2}$ となります。また、 $j$ 列目にあるどのマスも $k$ という値を取らない場合は $x_{1jk} = x_{2jk} = \cdots = x_{N^{2}jk} = 0$ となり、上式の評価値は $1$ となります。以下の図は、2列目に適当な数値を入れたときの評価値の例になります。

全ての取り得る列、値に対しても同様に議論できるので、「各列には、同じ数字が重複して入ってはいけない」というルールは以下の関数 $f_{3}$ で表現することができます。

$$f_{3}\left( \mathbf{x} \right) = \sum_{j=1}^{N^{2}} \sum_{k=1}^{N^{2}} \left( \sum_{i=1}^{N^{2}} x_{ijk} -1 \right)^{2}$$

4. $N \times N$ のブロックに分割した時に、各ブロックには同じ数字が重複して入ってはいけない

数独を $N \times N$ の $N^{2}$ 個のブロック$B_{1,1}, \cdots, B_{1,N}, \cdots, B_{r,c}, \cdots, B_{N,1}, \cdots, B_{N,N}$ に分割し、各ブロックには同じ数字が重複しないというルールを考えます。
$N = 2$のときに、ブロック $B_{r,c}$ の行番号、列番号は以下の通りになります。

ブロック	行番号、列番号
$B_{1,1}$	$\left\{ (1, 1), (1, 2), (2, 1), (2, 2) \right\}$
$B_{1,2}$	$\left\{ (1, 3), (1, 4), (2, 3), (2, 4) \right\}$
$B_{2,1}$	$\left\{ (3, 1), (3, 2), (4, 1), (4, 2) \right\}$
$B_{2,2}$	$\left\{ (3, 3), (3, 4), (4, 3), (4, 4) \right\}$

ブロック $B_{r,c}$ 内には $k$ という数字は1回しか現れてはいけない、というルールを次のような式で表現することできます。

$$\left(\sum_{(i, j) \in B_{r,c}} x_{ijk} - 1 \right)^{2}$$

数独のルールを満たすとき、すなわちブロック $B_{r,c}$ 内のいずれか1つの値が $1$ で、他が $0$ になるとき、上式の評価値は $0$ となります。数独のルールを満たさないとき、すなわちブロック内に値の重複が $m (\geq 2)$ 個存在する場合、 ブロック内の変数の $m$ 個の値が $1$ で、残りが $0$ となるので、上式の評価値は $(m - 1)^{2}$ となります。また、 ブロック内にあるマスのどれもが $k$ という値を取らない場合はブロック内の変数の値は全て $0$ となり、上式の評価値は $1$ となります。以下の図は、ブロック $B_{1,1}$ に適当な数値を入れたときの評価値の例になります。

全ての取りうる値、ブロックに対して同様に議論ができるので、「$N \times N$ のブロックに分割した時に、各ブロックには同じ数字が重複して入ってはいけない」というルールは以下の関数 $f_{4}$ で表現することができます。

$$f_{4}\left( \mathbf{x} \right) = \sum_{k=1}^{N^{2}} \sum_{r=1}^{N} \sum_{c=1}^{N} \left( \sum_{(i,j) \in B_{r,c}} x_{ijk} - 1 \right)^{2}$$

今回設計した数独の評価関数

上記で立式した評価関数を全て満たしてほしいので、最終的に用意する評価関数は以下のようになります。

$$\begin{aligned} f\left( \mathbf{x} \right) & = f_{1}\left( \mathbf{x} \right) + f_{2}\left( \mathbf{x} \right) + f_{3}\left( \mathbf{x} \right) + f_{4}\left( \mathbf{x} \right) \\ & = \sum_{i=1}^{N^{2}} \sum_{j=1}^{N^{2}} \left( \sum_{k=1}^{N^{2}} x_{ijk} - 1 \right)^{2} + \sum_{i=1}^{N^{2}} \sum_{k=1}^{N^{2}} \left( \sum_{j=1}^{N^{2}} x_{ijk} -1 \right)^{2} \\ & \quad + \sum_{j=1}^{N^{2}} \sum_{k=1}^{N^{2}} \left( \sum_{i=1}^{N^{2}} x_{ijk} -1 \right)^{2} + \sum_{k=1}^{N^{2}} \sum_{r=1}^{N} \sum_{c=1}^{N} \left( \sum_{(i,j) \in B_{r,c}} x_{ijk} - 1 \right)^{2} \end{aligned}$$

Ocean SDKで実装

数独の問題をQUBOで定式化できました。これをOcean SDKで実装していきます（実装は GitHub にあります ）。実装の手続きとしては以下のような流れになります。

1. 数独の問題を読み込みます
2. 読み込んだ数独の問題を QUBO の形式に変換します
3. モデル化した QUBO をソルバー（量子アニーリング、タブー法、焼きなまし法など）に送り、数独を解いてもらいます
4. 答えを出力し、正しいかどうか判定する

メインスクリプトは以下のようになっており、特に数独をQUBO形式に変換する関数 build_bqm と、実際に問題を解いている solve_sudoku 関数について解説していきたいと思います。

from __future__ import print_function

import typer

from puzzle import get_matrix, is_correct
from models import build_original_bqm, build_bqm
from solvers import solve_sudoku


def main(
        filename: str = "problem.txt",
        solver_name: str = "qpu",
):
    # 問題を読み込む
    mat = get_matrix(filename)
    # 数独の問題を QUBO 形式に変換する
    bqm = build_bqm(mat)
    # QUBO形式にした数独の問題を解く
    result, energy = solve_sudoku(bqm, mat, solver_name)

    for line in result:
        print(*line, sep=" ")   # Print list without commas or brackets

    print("bqm energy: ", energy)

    # Verify
    if is_correct(result):
        print("The solution is correct")
    else:
        print("The solution is incorrect")


if __name__ == "__main__":
    typer.run(main)

QUBO形式に変換する build_bqm

Ocean SDKでは、Binary Quadratic Model (BQM) と呼ばれるクラスを使って、QUBOの定式化をしていきます。
まず以下のように、モデル bqm を初期化していきます。

import math
import itertools

import dimod
from dimod.generators import combinations

# 行番号、列番号、値を引数にとり、変数の文字列表現に変換します
# （例） var_repr(1, 1, 3) -> '1,1_3'
from puzzle import var_repr


# 数独に関する基本情報
matrix = [
    [0, 1, 0, 0],
    [2, 0, 0, 0],
    [0, 0, 0, 0],
    [4, 0, 0, 1],
]
n = len(matrix)
m = int(math.sqrt(n))
digits = list(range(1, n+1)) # マスに入る数字一覧

# QUBOによる定式化
linear = {} # 1次の項
quadratic = {} # 2次の項
offset = 0.
vartype = dimod.BINARY
bqm = dimod.BinaryQuadraticModel(linear, quadratic, offset, vartype)

bqm に項を追加するメソッドは add_linear 、 add_quadratic など様々なものが提供されていますが、今回は update メソッドを主に使っていきます。update メソッドでは2つのBQMを統合します。評価関数 $f_{1}, f_{2}, f_{3}, f_{4}$ に対応するBQMを作成し、それぞれのBQMを bqm.update で追加・統合していきます。
評価関数 $f_{1}, f_{2}, f_{3}, f_{4}$ はすべて $\left(\sum_{i} x_{i} - k \right)^{2}$ （ただし $k$ は実数）の形式の関数であり、この形式の関数をBQMに変換するメソッド dimod.generators.combinations を使うと、各評価関数のBQMを簡単に作成することができます。このメソッドの第1引数には取りうる変数のリストを、第2引数には $k$ の値を指定します。

1. 各マスには $1 \sim N^{2}$ の値が入る

以下のコードは、評価関数 $f_{1}(\mathbf{x}) = \sum_{i} \sum_{j} \left( \sum_{k} x_{ijk} - 1 \right)^{2}$ に対応しています。

for i, j in itertools.product(range(n), range(n)):
    node_digits = [var_repr(i, j, v) for v in digits]
    bqm.update(combinations(node_digits, 1))

2. 各行には、同じ数字が重複して入ってはいけない

以下のコードは、評価関数 $f_{2}(\mathbf{x}) = \sum_{i} \sum_{k} \left( \sum_{j} x_{ijk} - 1 \right)^{2}$ に対応しています。

for i, v in itertools.product(range(n), digits):
    row_nodes = [var_repr(i, j, v) for j in range(n)]
    bqm.update(combinations(row_nodes, 1))

3. 各列には、同じ数字が重複して入ってはいけない

以下のコードは、評価関数 $f_{3}(\mathbf{x}) = \sum_{j} \sum_{k} \left( \sum_{i} x_{ijk} - 1 \right)^{2}$ に対応しています。

for j, v in itertools.product(range(n), digits):
    col_nodes = [var_repr(i, j, v) for i in range(n)]
    bqm.update(combinations(col_nodes, 1))

4. $N \times N$ のブロックに分割した時に、各ブロックには同じ数字が重複して入ってはいけない

以下のコードは、評価関数 $f_{4}(\mathbf{x}) = \sum_{k} \sum_{r} \sum_{c} \left( \sum_{i,j \in B_{r,c}} x_{ijk} - 1 \right)^{2}$ に対応しています。

for digit in digits:
    for mi, mj in itertools.product(range(m), range(m)):
        ij = itertools.product(
            range(mi * m, (mi + 1) * m),
            range(mj * m, (mj + 1) * m),
        )
        square_nodes = [var_repr(i, j, digit) for i, j in ij]
        bqm.update(combinations(square_nodes, 1))

既に値が入っているマスに対応する変数は 1 で固定する

数独では問題が与えられたとき、既にいくつかのマスについては値がセットされています。すなわち、すでに正解（対応する変数 $x_{ijk}$ が $0$ をとるか $1$ をとるか）を知っていることになるので、値がセットされているマスに対応する変数 $x_{ijk}$ は 1 にしておきます。bqm の fix_variable あるいは fix_variables を用いて使って既知の変数の値を固定しておきます。

fixed_vars: set[tuple[str, int]] = set()
for i, row in enumerate(matrix):
    for j, val in enumerate(row):
        if val > 0:
            # マス(i,j) は既知なので 1 にする
            fixed_vars |= {(var_repr(i, j, val), 1)}
bqm.fix_variables(list(fixed_vars))

これで数独の問題をQUBOに変換できました。 Binary Quadratic Model bqm をソルバーに渡して解いていきます。

評価関数 $f(\mathbf{x})$ を最小化するようなマスの値を当てる solve_sudoku

Binary Quadratic Modelの形式になった組合せ最適化問題を解くためのソルバーが Ocean SDK には用意されています。ここでは4つのソルバーについて紹介していきます。

方式	概説
タブー法 InterruptableTabuSampler	一度探索した解を再探索しないように、最適化問題の解を探索していく手法。結果が悪くなる場合でも強制的に探索することで局所解で止まることを防いでいる。
シミュレーテッド・アニーリング InterruptableSimulatedAnnealingProblemSampler	古典コンピュータ上で実行される焼きなまし法。
量子アニーリング QPUSubproblemAutoEmbeddingSampler	D-Wave社の量子アニーリングマシンで実行される最適化手法。
ハイブリット法 KerberosSampler	タブー法とアニーリングと、量子アニーリング法のハイブリッド方式。3つの方式を並列で実行し、最も良い結果を採用する。

solver_name = "qpu"

solver = {
    'tabu': solve_by_tabu_search,
    'sa': solve_by_sa,
    'qpu': solve_by_qpu_subproblem_sampler,
    'kerberos': solve_by_kerberos,
}.get(solver_name, solve_by_qpu_subproblem_sampler)

solution = solver(bqm)

各変数が取る値は solution.first.sample に格納されています。辞書型のオブジェクト（変数名 → {0, 1}）で、値が 1 になっている変数が数独の答えになります。

import copy

# 数独の答えになるのは1の値を取る変数
best_solution = solution.first.sample
solution_list = [k for k, v in best_solution.items() if v == 1]

result = copy.deepcopy(matrix)
energy = solution.first.energy

for label in solution_list:
    row, col, val = parse_label(label)
    if result[row][col] > 0:
        continue
    result[row][col] = val

実行結果

数独の問題をQUBOで解けるように定式化ができましたので、実際に量子アニーリングで数独を解いてみましょう。実行環境を用意します。

アクセストークン

D-Waveの量子アニーリングマシンにクラウド経由でアクセスできるサービスleapでアカウントを作成し、アクセストークンを発行します。発行したアクセストークンは環境変数 DWAVE_API_TOKEN に設定してください。

export DWAVE_API_TOKEN=*************************************

python環境の準備

GitHub からコードを取得して、pipenv で依存するパッケージをインストールします。

$ git clone https://github.com/AIREVinc/quantum-annealing-leap-playground
$ cd quantum-annealing-leap-playground
$ pipenv --python 3.9.4
$ pipenv shell

数独を解くプログラムの実行

数独を解くプログラムは main.py になります。以下のコマンド引数を指定します。

引数	概説
--filename	数独の問題が記載されているテキストファイル。
--solver-name	使用するソルバー名。tabu、 sa、 qpu、 kerberos が指定できます。

実行すると以下のような出力が得られます。

(quantum-annealing-leap-playground) $ cd sudoku
(quantum-annealing-leap-playground) $ python main.py \
  --filename problems/4x4/problem-00.txt \
  --solver-name qpu

3 1 4 2
2 4 1 3
1 3 2 4
4 2 3 1
bqm energy:  0.0
The solution is correct

このように、数独を解いた結果が得られました。

まとめ

この記事では、量子アニーリングとはなにか、どう使っていくのかについて、D-Wave社が提供する量子アニーリングのクラウドサービスである Leap とそのSDKである Ocean を使い、数独の例題を通して解説しました。量子アニーリングを使ってを解くには、問題を QUBO の形式に落とし込む必要があり、ここの設計ができるかがポイントになるかと思います。今回は $N = 2$ の場合（$4 \times 4$ マス）の数独を解きましたが、 $N = 3$ のとき、すなわち $9 \times 9$ マスの数独も同じように適用することができます。数独の場合、最適化の対象となる変数の数は高々 （行数 $N^{2}$） $\times$ （列数 $N^{2}$） $\times$ （数字の種類 $N^{2}$） と $N^{6}$ 個で、 この記事で解説した数独の問題は $4 \times 4$ と、比較的小さいサイズの問題を扱いました。変数の数がより多い場合や、その他の実応用に対し、量子アニーリングができること、量子アニーリングで解くときの限界については今後も調査していきたいと思います。

参考文献・参考資料

• (PDF) Quantum annealing in the transverse Ising model
• 量子アニーリングの基礎
• Ocean SDK
• d-wave examples