👣

量子アニーリングを使って迷路を作ってみた

2024/04/07に公開

Python

はじめに

オンライン公開伴走型量子コンピューティング講座QA4U2卒業生のluna_moonlightです。この記事では、QA4U2発の論文"Individual subject evaluated difficulty of adjustable mazes generated using quantum annealing"で提案された量子アニーリングを用いた迷路生成について解説していきます。

量子アニーリングとは

量子アニーリングとは、量子ゆらぎと呼ばれる現象を用いた組合せ最適化問題^[1]を解く手法です。

組合せ最適化問題を解く古典的な手法のひとつとして、焼きなまし法というアルゴリズムがあります。焼きなまし法では、温度というパラメータを導入し、最初は高温に設定して局所解を脱しながら解を探索していき、徐々に温度を低くしていくことで最適解へ遷移していきます。量子アニーリングは、温度の代わりに量子ゆらぎを導入して最適解を探索するアルゴリズムとして、1998年に"Quantum annealing in the transverse Ising model"という論文で門脇正史と西森秀稔により提案されました。

迷路生成アルゴリズム

迷路生成は古くから研究され、さまざまなアルゴリズムが考案されています。有名な迷路生成アルゴリズムとして、棒倒し法、壁伸ばし法、穴掘り法などがあります。今回はその中でも棒倒し法のアルゴリズムを使用して迷路を生成していきます。

棒倒し法

棒倒し法は以下の手順で迷路を生成します。

迷路の外周を壁にして、その他は通路にする。
外周の内側のx,yが共に奇数の位置に壁(棒)を設置する。
すべての棒を、上下左右のいずれか一方に倒して壁とする。(ただし倒し方に以下の制約がある)
- 左から1列目の棒以外は左に倒してはいけない。
- すでに棒が倒れている方向に重ねて倒してはいけない。

棒倒し法で迷路を生成している様子がこちらです。

壁伸ばし法

壁伸ばし法は以下の手順で迷路を生成します。

迷路の外周を壁にして、その他は通路にする。
外周の内側のx,yが共に奇数の位置を壁伸ばしの起点候補としてリストにまとめる。
すべての起点候補が壁になるまで、通路である起点候補からランダムに起点とし、4.を実行する。
選んだ起点から壁を伸ばす。
i. 選んだ起点を仮壁にする。
ii. 上下左右のいずれか一方をランダムで選択し、状況に合わせて以下を実行する。
- 2マス先が通路の場合、1マス先を仮壁にして2マス先を起点とし、4.を再帰的に実行する。
- 2マス先が壁の場合、1マス先を仮壁にしてすべての仮壁を壁に変更し、3.に戻る。
- 2マス先が仮壁の場合、仮壁をすべて削除^[2]し、3に戻る。

壁伸ばし法で迷路を生成している様子がこちらです。

穴掘り法

穴掘り法は以下の手順で迷路を生成します。

迷路の盤面すべてを壁にする。
外周の内側のx,yが共に偶数の位置を穴掘りの起点候補としてリストにまとめる。
起点候補からランダムに起点を選び、通路にする。
すべての起点候補が通路になるまで、選んだ起点から穴を掘る。
i. 選んだ起点を通路にする。
ii. 上下左右のいずれか一方をランダムで選択し、状況に合わせて以下を実行する。
- 2マス先が通路の場合、通路である起点候補からランダムに起点とし、4.を再帰的に実行する。
- 2マス先が壁の場合、1マス先を通路にして2マス先を起点とし、4.を再帰的に実行する。

穴掘り法で迷路を生成している様子がこちらです。

量子アニーリングで迷路を生成してみる

それでは実際に量子アニーリングを用いて棒倒し法のアルゴリズムを実装していきます。

定式化

量子アニーリングで解けるようにするには、非制約二値変数二次最適化(QUBO: Quadratic Unconstrained Binary Optimization)の形に落とし込む必要があります。QUBOでは、 $0$ と $1$ の2つの値をとる $M$ 個のbinary変数 $x_1, x_2, ..., x_M$ に対して、二次形式の評価関数 $f:\{0, 1\}^M \to \Reals$ を最小にするように変数の組み合わせ $\bm{x}=(x_1, x_2, ..., x_M)$ を求めていきます。すなわち、

f(\bm{x})=\bm{x}Q\bm{x}^T=\sum_{i=1}^M\sum_{j=1}^M Q_{i,j}x_ix_j

を最小にする $\bm{x}=(x_1, x_2, ..., x_M)$ を求めます。ここで、行列 $Q \in \Reals ^ {M \times M}$ は、対角成分は $x_i$ の重みを、非対角成分は $x_i, x_j$ 成分の相互作用を表しており、QUBO行列と呼ばれます。

ここで、棒倒し法のアルゴリズムについて考えてみましょう。 $x_{i,j,d}$ は、上から $i$ 行目、左から $j$ 列目の棒を方向 $d$ (0:上、1:右、2:下、3:左)に倒すか否かを表すこととします。例えば、上から2行目、左から1列目の棒を左に倒す場合は $x_{2,1,3}=1$ に、上から1行目、左から3列目の棒を下に倒さない場合は $x_{1,3,2}=0$ になります。棒の数が $H \times W$ 本の迷路を生成する場合の変数の組み合わせは、

\bm{x}=(x_{1,1,0}, x_{1,1,1}, ..., x_{1,1,3}, x_{2,1,0}, ..., x_{1,2,0}, x_{1,2,1}, x_{1,2,2}, x_{2,2,0}, ..., x_{i,j,0}, x_{i,j,1}, x_{i,j,2}, ..., x_{H,W,2})

となります。量子アニーリングで解けるような形に落とし込むために、この組み合わせ $\bm{x}$ が正しい迷路を生成するときは小さな値を、間違った迷路を生成するときは大きな値をとるように評価関数 $f$ (QUBO行列 $Q$ )を決定していきます。

1.すでに棒が倒れている方向に重ねて倒してはいけない

棒倒し法のアルゴリズムでは、すでに棒が倒れている方向に重ねて倒してはいけないという制約があります。すでに棒が倒れている方向に重ねて倒れてしまうと、経路が複数存在する迷路ができてしまいます。この制約は以下のように表すことができます。

f_1(\bm{x})=\sum_{i_1, i_2}\sum_{j_1, j_2}\sum_{d_1, d_2} Q_{(i_1,j_1,d_1),(i_2,j_2,d_2)}x_{i_1,j_1,d_1}x_{i_2,j_2,d_2}

ただし、

Q = \left\{ \begin{array}{ll} 1 & (i_1=i_2+1, j_1=j_2, d_1=2, d_2=0)\\ 1 & (i_1=i_2-1, j_1=j_2, d_1=0, d_2=2)\\ 0 & (otherwise) \end{array} \right.

です。棒が重なる場合、例えば $i_1=1, j_1=3, d_1=2, i_2=2, j_2=3, d_2=0$ のとき。すなわち、上から1行目、左から3列目の棒が下に、上から2行目、左から3列目の棒が上に倒れる場合、棒が重なってしまうため $Q_{(1,3,2),(2,3,0)}x_{1,3,2}x_{2,3,0}=1$ をとります。逆の場合も同様に、 $Q_{(2,3,0),(1,3,2)}x_{2,3,0}x_{1,3,2}=1$ をとります。それ以外の棒が重ならない場合、例えば $i_1=3, j_1=1, d_1=0, i_2=4, j_2=1, d_2=0$ のとき。すなわち、上から3行目、左から1列目の棒と上から4行目、左から1列目の棒がともに上に倒れる場合、棒は重ならないため $Q_{(3,1,0),(4,1,0)}x_{3,1,0}x_{4,1,0}=0$ をとります。

2.棒は1方向にしか倒してはいけない

棒倒し法のアルゴリズムでは、棒は1方向にしか倒してはいけないという制約が暗にあります。棒が倒れない場合は経路が複数存在する迷路が、棒が2方向以上に倒れる場合は辿り着けない場所が存在する迷路ができてしまいます。この制約は以下のように表すことができます。

f_2(\bm{x})=\sum_{i}\sum_{j}\left(\sum_{d} x_{i,j,d} - 1\right)^2

棒が倒れない場合、例えば上から1行目、左から3列目について考えると、 $\left(\sum_{d} x_{1,3,d} - 1\right)^2=(x_{1,3,0}+x_{1,3,1}+x_{1,3,2} - 1)^2=1$ をとります。
棒が全方向に倒れる場合、例えば2行目、左から1列目について考えると、 $\left(\sum_{d} x_{2,1,d} - 1\right)^2=(x_{2,1,0}+x_{2,1,1}+x_{2,1,2}+x_{2,1,3} - 1)^2=9$ をとります。
ちょうど1方向にのみ倒れる場合、例えば3行目、左から2列目について考えると、 $\left(\sum_{d} x_{3,2,d} - 1\right)^2=(x_{3,2,0}+x_{3,2,1}+x_{3,2,2} - 1)^2=0$ をとります。

迷路生成の評価関数

ここまでで、制約ひとつひとつに対して評価関数を作ってきました。これらを足し合わせることで、迷路生成に必要な評価関数を作ることができます。具体的には以下のようになります。

\begin{split} f(\bm{x}) &= f_1(\bm{x})+ \lambda f_2(\bm{x})\\ &= \sum_{i_1, i_2}\sum_{j_1, j_2}\sum_{d_1, d_2} Q_{(i_1,j_1,d_1),(i_2,j_2,d_2)}x_{i_1,j_1,d_1}x_{i_2,j_2,d_2} + \lambda\sum_{i}\sum_{j}\left(\sum_{d} x_{i,j,d} - 1\right)^2 \end{split}

実装

ここまでで、迷路生成の評価関数をQUBOの形で定式化することができました。それでは早速実装していきましょう。

まず、D-Waveの量子アニーリングマシンを使うために、D-Wave Leapでアカウントを作成します。

次に、必要な外部ライブラリをインストールしましょう。量子アニーリングを使うためにdwave-ocean-sdkと迷路の可視化のためにmatplotlibをインストールします。

pip install dwave-ocean-sdk
pip install matplotlib

必要な外部ライブラリをインストールしたところで、さっそくプログラムを書いていきましょう。コードの全体像はこちらです。

コードの全体像

from collections import defaultdict
import matplotlib.pyplot as plt
from dwave.system import DWaveSampler, EmbeddingComposite


'''-----パラメータ-----'''
STICK_NUM_H = 5 # 縦の棒の本数
STICK_NUM_W = 7 # 横の棒の本数
l = 1           # 罰金項の強さ
num_reads = 10  # アニーリングを実行する回数
token = 'XXXXX' # API token(個人のものを使用)
solver='Advantage_system6.4' # 量子アニーリングマシン
'''------------------'''

Q = defaultdict(lambda: 0) # Q_i,j,d   (i, j)における棒を倒す方向d
DIRECTION = [(-1, 0), (0, 1), (1, 0), (0, -1)] # 棒の倒れる方向

# x_i,j,dの通し番号を記録
ijd_to_idx = {}
idx = 0
for i in range(1, STICK_NUM_H + 1):
  for j in range(1, STICK_NUM_W + 1):
    # 左から1列目の場合
    if j == 1:
      for d in range(4):
        ijd_to_idx[(i, j, d)] = idx
        idx += 1
    # 左から2列目以降の場合
    else:
      for d in range(3):
        ijd_to_idx[(i, j, d)] = idx
        idx += 1

# 迷路の初期盤面を作成
field = [[0 for i in range(2 * STICK_NUM_W + 3)] for j in range(2 * STICK_NUM_H + 3)]
for j in range(2 * STICK_NUM_W + 3):
  field[0][j] = 1
  field[-1][j] = 1
for i in range(2 * STICK_NUM_H + 3):
  field[i][0] = 1
  field[i][-1] = 1
for i in range(2, 2 * STICK_NUM_H + 3, 2):
  for j in range(2, 2 * STICK_NUM_W + 3, 2):
    field[i][j] = 1



# 倒れる方向が重ならないように制限
for i in range(1, STICK_NUM_H + 1):
  for j in range(1, STICK_NUM_W + 1):
    # 左から1列目の場合
    if i != STICK_NUM_H:
      Q[(ijd_to_idx[i, j, 2], ijd_to_idx[i+1, j, 0])] += 1
    if i != 1:
      Q[(ijd_to_idx[i, j, 0], ijd_to_idx[i-1, j, 2])] += 1
        
# 倒れる方向をひとつに制限
for i in range(1, STICK_NUM_H + 1):
  for j in range(1, STICK_NUM_W + 1):
    # 左から1列目の場合
    if j == 1:
      for d1 in range(4):
        for d2 in range(4):
          Q[(ijd_to_idx[(i, j, d1)], ijd_to_idx[(i, j, d2)])] += l
        Q[(ijd_to_idx[(i, j, d1)], ijd_to_idx[(i, j, d1)])] -= 2 * l
    # 左から2列目以降の場合
    else:
      for d1 in range(3):
        for d2 in range(3):
          Q[(ijd_to_idx[(i, j, d1)], ijd_to_idx[(i, j, d2)])] += l
        Q[(ijd_to_idx[(i, j, d1)], ijd_to_idx[(i, j, d1)])] -= 2 * l



# 量子アニーリングの実行
endpoint = 'https://cloud.dwavesys.com/sapi/'
dw_sampler = DWaveSampler(solver=solver, token=token, endpoint=endpoint)
sampler = EmbeddingComposite(dw_sampler)
sampleset = sampler.sample_qubo(Q, num_reads=num_reads)



# 結果に基づいて迷路を生成
result = [i for i in sampleset.first[0].values()]
for i in range(1, STICK_NUM_H + 1):
  for j in range(1, STICK_NUM_W + 1):
    # 左から1列目の場合
    if j == 1:
      for d in range(4):
        if result[ijd_to_idx[(i, j, d)]] == 1:
          field[2 * i + DIRECTION[d][0]][2 * j + DIRECTION[d][1]] = 1
    # 左から2列目の場合
    else:
      for d in range(3):
        if result[ijd_to_idx[(i, j, d)]] == 1:
          field[2 * i + DIRECTION[d][0]][2 * j + DIRECTION[d][1]] = 1

# 迷路の表示
plt.imshow(field, cmap='gray_r')
plt.axis('off')
plt.show()

ここからは、このコードを分けて解説していきます。

初期化

以下のコードは、必要なライブラリのインポートや変数、定数を定義しています。パラメータと書かれている部分を変更することで、迷路の大きさや使う量子アニーリングマシンなどを変更することができます。

from collections import defaultdict
import matplotlib.pyplot as plt
from dwave.system import DWaveSampler, EmbeddingComposite


'''-----パラメータ-----'''
STICK_NUM_H = 5 # 縦の棒の本数
STICK_NUM_W = 7 # 横の棒の本数
l = 1           # 罰金項の強さ
num_reads = 10  # アニーリングを実行する回数
token = 'XXXXX' # API token(個人のものを使用)
solver='Advantage_system6.4' # 量子アニーリングマシン
'''------------------'''

Q = defaultdict(lambda: 0) # Q_i,j,d   (i, j)における棒を倒す方向d
DIRECTION = [(-1, 0), (0, 1), (1, 0), (0, -1)] # 棒の倒れる方向

# x_i,j,dの通し番号を記録
ijd_to_idx = {}
idx = 0
for i in range(1, STICK_NUM_H + 1):
  for j in range(1, STICK_NUM_W + 1):
    # 左から1列目の場合
    if j == 1:
      for d in range(4):
        ijd_to_idx[(i, j, d)] = idx
        idx += 1
    # 左から2列目以降の場合
    else:
      for d in range(3):
        ijd_to_idx[(i, j, d)] = idx
        idx += 1

# 迷路の初期盤面を作成
field = [[0 for i in range(2 * STICK_NUM_W + 3)] for j in range(2 * STICK_NUM_H + 3)]
for j in range(2 * STICK_NUM_W + 3):
  field[0][j] = 1
  field[-1][j] = 1
for i in range(2 * STICK_NUM_H + 3):
  field[i][0] = 1
  field[i][-1] = 1
for i in range(2, 2 * STICK_NUM_H + 3, 2):
  for j in range(2, 2 * STICK_NUM_W + 3, 2):
    field[i][j] = 1

QUBO行列 $Q$ の設定

ここでは、定式化した迷路生成の評価関数をプログラムに載せています。

1.すでに棒が倒れている方向に重ねて倒してはいけない

以下のコードは、 $f_1(\bm{x})=\sum_{i_1, i_2}\sum_{j_1, j_2}\sum_{d_1, d_2} Q_{(i_1,j_1,d_1),(i_2,j_2,d_2)}x_{i_1,j_1,d_1}x_{i_2,j_2,d_2}$ に対応しています。

# 倒れる方向が重ならないように制限
for i in range(1, STICK_NUM_H + 1):
  for j in range(1, STICK_NUM_W + 1):
    # 左から1列目の場合
    if i != STICK_NUM_H:
      Q[(ijd_to_idx[i, j, 2], ijd_to_idx[i+1, j, 0])] += 1
    if i != 1:
      Q[(ijd_to_idx[i, j, 0], ijd_to_idx[i-1, j, 2])] += 1

2.棒は1方向にしか倒してはいけない

以下のコードは、 $f_2(\bm{x})=\sum_{i}\sum_{j}\left(\sum_{d} x_{i,j,d} - 1\right)^2$ に対応しています。

# 倒れる方向をひとつに制限
for i in range(1, STICK_NUM_H + 1):
  for j in range(1, STICK_NUM_W + 1):
    # 左から1列目の場合
    if j == 1:
      for d1 in range(4):
        for d2 in range(4):
          Q[(ijd_to_idx[(i, j, d1)], ijd_to_idx[(i, j, d2)])] += l
        Q[(ijd_to_idx[(i, j, d1)], ijd_to_idx[(i, j, d1)])] -= 2 * l
    # 左から2列目以降の場合
    else:
      for d1 in range(3):
        for d2 in range(3):
          Q[(ijd_to_idx[(i, j, d1)], ijd_to_idx[(i, j, d2)])] += l
        Q[(ijd_to_idx[(i, j, d1)], ijd_to_idx[(i, j, d1)])] -= 2 * l

量子アニーリングの実行

以下のコードは、solver(今回はAdvantage_system6.4)という量子アニーリングマシンで、QUBO行列 $Q$ に基づいて、num_reads(今回は10)回量子アニーリングを実行しています。

# 量子アニーリングの実行
endpoint = 'https://cloud.dwavesys.com/sapi/'
dw_sampler = DWaveSampler(solver=solver, token=token, endpoint=endpoint)
sampler = EmbeddingComposite(dw_sampler)
sampleset = sampler.sample_qubo(Q, num_reads=num_reads)

迷路の可視化

以下のコードは、num_reads回のアニーリングの中で最も上手くいった結果に基づいて迷路を生成し、それを可視化しています。

# 結果に基づいて迷路を生成
result = [i for i in sampleset.first[0].values()]
for i in range(1, STICK_NUM_H + 1):
  for j in range(1, STICK_NUM_W + 1):
    # 左から1列目の場合
    if j == 1:
      for d in range(4):
        if result[ijd_to_idx[(i, j, d)]] == 1:
          field[2 * i + DIRECTION[d][0]][2 * j + DIRECTION[d][1]] = 1
    # 左から2列目の場合
    else:
      for d in range(3):
        if result[ijd_to_idx[(i, j, d)]] == 1:
          field[2 * i + DIRECTION[d][0]][2 * j + DIRECTION[d][1]] = 1

# 迷路の表示
plt.imshow(field, cmap='gray_r')
plt.axis('off')
plt.show()

結果

上のコードを実行すると以下のような出力が得られます。

このように、迷路を生成した結果が得られました。

まとめ

この記事では、量子アニーリングを用いて迷路を生成する方法について解説しました。参考にした論文"Individual subject evaluated difficulty of adjustable mazes generated using quantum annealing"では、今回紹介した迷路生成に加えて、スタートとゴールを決定する方法や個人に対して迷路の難易度を上昇させる方法についても述べられています。これらについても解説していると非常に長くなってしまうため、今回は最小限の迷路生成について解説しました。

脚注

組合せ最適化問題とは、多くの選択肢からある制約をもとに最も良い解を求める問題のこと。ナップサック問題やSAT問題などが挙げられる。 ↩︎
4.の操作が続けられる場所まで戻る方法もあるが、今回は簡単のためすべて削除してやり直す方法をとっている。 ↩︎

はじめに

量子アニーリングとは

迷路生成アルゴリズム

棒倒し法

壁伸ばし法

穴掘り法

量子アニーリングで迷路を生成してみる

定式化

1.すでに棒が倒れている方向に重ねて倒してはいけない

2.棒は1方向にしか倒してはいけない

迷路生成の評価関数

実装

初期化

QUBO行列Qの設定

1.すでに棒が倒れている方向に重ねて倒してはいけない

2.棒は1方向にしか倒してはいけない

量子アニーリングの実行

迷路の可視化

結果

まとめ

Discussion

QUBO行列 $Q$ の設定