T-QARD Harbor

N-クイーン問題 をD-Waveマシンで解く

N-クイーン問題とは

導入

5 × 5マスのチェス盤に5個のクイーンをどのように置いたら、お互いの効きが当たらないか考えてみよう。クイーンは、縦横斜めに他の駒によって遮られない限り進める。

5×5マスのチェス盤と5個のクイーン
5×5マスのチェス盤と5個のクイーン

これを一般化した問題はN-クイーン問題 (N-Queens Problem)と呼ばれ、$N \times N$マスのチェス盤にN個のクイーンをお互いの効きが当たらないように配置する問題である。

さて、先ほどの場合の答えは分かっただろうか。もし分からなくても、愚直にクイーンの配置を試していけば見つかるかもしれない。しかし、クイーンの可能な配置は${}_{25} C _5=53,130$通り存在し、これには時間と忍耐力を要する。実際には、以下の2つの基本解が存在する。

5×5マスの場合の基本解
5×5マスの場合の基本解

また、基本解の回転と反転によるバリエーション解がある。左の基本解に対して回転と反転により8パターンの変形、右の基本解は点対称であるため2パターンの変形があり、合わせて10個のバリエーション解が見つかっている。$N=5$の場合であれば大したことはないが、Nが大きくなると解を見つけていくことが困難になることは想像に難くない。このように膨大な組み合わせが存在する問題については、動画「『フカシギの数え方』 おねえさんといっしょ! みんなで数えてみよう!」が面白い。

これまでに見つかった解

人類がすべての解を見つけることができた最大の場合は、$N=27$である。

N基本解バリエーション解可能な配置数
1111
2006
30084
4121820
521053130
6141947792
764085900584
812924426165368
946352260887834350
109272417310309456440
1134126801276749965026536
12178714200103619293824707388
139233737129176358300744339432
1445752365596880530516383349192480
15285053227918491005567811177478095440
1618469551477251210078751602022313874633200
1711977939958151041190739044344491048895397910
1883263591666090624149482492334195165714038760136
19621012754496805784819870867053543756004133247695400
204878666808390291888842788360983670896737872851072994080
2139333324973314666222712411887396336567398822620727355402190
22336376244042269100870164463887407766986865702182544710470036560
2330292426582102423393768444010381758958529585222885358558747563185920
24284392729569342275141719737361763783520005146433232953016554504214270600
252759866837434342207893435808352312690620414907617326451186807195415244296900
2627897124665102892231769961636404457746226578042013138408988185727715132037050952
272936379196767819923490796715412252811091107763254898773425731705373527055193637625824
表: これまでに見つかっている問題の大きさと解の数

歴史

N-クイーン問題の歴史は150年以上あり、チェスにおけるパズルとして始まった。

1848年チェスコンポーザー(チェスに関する問題を作成する人)のMax Bezzelが8-クイーン問題を発表した。
1850年Franz Nauckが8-クイーン問題の初めての解を発見。N-クイーン問題に拡張した。

これ以降、Carl Friedrich Gaussをはじめとする数学者がこれらの問題を研究の対象とし始めた。

1874年S.Guntherは解を求めるために、行列式を用いる方法を提案し、後にJWL Glaisherがその方法を改良した。
1972年Edsger Dijkstraは、この問題を題材として構造化プログラミングや深さ優先探索、バックトラックアルゴリズムを説明した。
2009年ドレスデン工科大学で、26-クイーン問題が計算された。
2016年Q27 Projectによって、27-クイーン問題が計算された。

本記事で行うこと

今回の記事では、

  • 問題の定式化を変えることで、D-Waveマシンへの埋め込みや得られる結果にどのような違いが見られるか実験する
  • 得られる解とエネルギー分布から、D-Waveマシンの特性を調べる

を主な目標とする。

D-Waveマシンによる実験

QUBO表現による問題の定式化

変数:$q_{i,j} \in \{0, 1\}$ (i行j列にクイーンが配置される場合に$1$, そうでなければ$0$)

制約条件は4つからなり、クイーンの総数に関する制約条件を異なる形で加えた場合を含む。

各行に配置するクイーンは1個のみ

$$ H_\mathrm{row} = \lambda_\mathrm{rowcol} \sum_{i=0}^{N-1} \left( 1-\sum_{j=0}^{N-1} q_{ij} \right)^2 \tag{1} $$

各列に配置するクイーンは1個のみ

$$ H_\mathrm{column} = \lambda_\mathrm{rowcol} \sum_{j=0}^{N-1} \left( 1-\sum_{i=0}^{N-1} q_{ij} \right)^2 \tag{2} $$

行と列の制約項については、他の組合せ最適化問題にもよく表れる形のため説明は割愛する。
(例えば、「宮城県の市区町村を塗り分けろ – グラフ彩色問題 をD-Waveマシンで解く」を参照)

斜めに配置するクイーンは0または1個のみ

$$ \begin{aligned} & H_\mathrm{diagonal} \\ &= \lambda_\mathrm{diagonal} \left( \sum_{s=1}^{2N-3} \left(\frac{1}{2}-\sum_{i,j;i+j = s} q_{ij} \right)^2 + \sum_{s=-\lfloor \frac{2N-3}{2} \rfloor}^{\lfloor \frac{2N-3}{2} \rfloor} \left(\frac{1}{2}-\sum_{i,j;i-j = s} q_{ij} \right)^2 \right) \tag{3} \end{aligned} $$

それぞれの斜めに配置するクイーンの数が0または1のときに限り、式(3)の2乗の部分が最小値$\frac{1}{4}$となるようにしている。

N個のクイーンを配置する

[1]クイーンの総数をなるべくNに近づける
素朴に2乗の項を用いて、クイーンの総数がNのときに最小値0となるようにする。

$$ H_\mathrm{queens} = \lambda_\mathrm{queens} \left(N-\sum_{i=0}^{N-1} \sum_{j=0}^{N-1} q_{ij} \right)^2 \tag{4-1} $$

[2]クイーンの総数をなるべく多くする
他の制約が存在するため、このように単純にしても上手くいくと考えられる。

$$ H_\mathrm{queens} = -\lambda_\mathrm{queens} \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} q_{ij} \tag{4-2} $$

[3]制約項を加えず、行と列の制約を強くする
行と列の制約を満たす場合クイーンの総数はNとなるため、$H_{row}$と$H_{column}$の効きを強くする($\lambda_{rowcol}$の値を大きくする)。

$$ H_\mathrm{queens} = 0 \tag{4-3} $$

注意として、今回考える場合以外にも様々な制約の与え方が存在する。なお、行と列は対称であるため、(1)と(2)の項には同じパラメータを設定する。同様に、2方向の斜めに関する部分も同じパラメータとする。
よって、最小化する目的関数は以下のようになる。

目的関数:

$$ H = H_\mathrm{row} + H_\mathrm{column} + H_\mathrm{diagonal} + H_\mathrm{queens} \tag{5} $$

各制約を満たす場合の各項の値は、$H_{row}=0$、$H_{column}=0$、$H_{diagonal}=\frac{1}{2} \times \left( 2N-3\right) = N-\frac{3}{2}$、式(4-1),(4-3)のとき$H_{queens}=0$、式(4-2)のとき$H_{queens}=-N$である。よって、すべての制約を満たすときの目的関数のエネルギーは、式(4-1),(4-3)のとき$H=N-\frac{3}{2}$、式(4-2)のとき$H=-\frac{3}{2}$となる。

QUBO行列

D-Wave社が提供している埋め込みアルゴリズムを利用すると、$N=8$の場合における隣接行列は以下のようになった。式(4-1)の場合は、式(4-1)の$H_{queens}$からも分かるように全結合であり、必要な物理量子ビットの数が多くなる。式(4-2),(4-3)の場合、先ほどよりスパースであるが、対角成分付近にある行の制約部分が密な結合となっている。

式(4-1)による制約の場合の隣接行列
式(4-1)による制約の場合の隣接行列

参考までに、式(4-2),(4-3)の場合における隣接行列を以下に示す。ただし、クイーンの総数に関する制約項$H_{queens}$は対角成分のみである。

式(4-2),(4-3)による制約の場合の隣接行列
式(4-2),(4-3)による制約の場合の隣接行列

チェス盤の行と列は対称であるものの、上の図では行と列の成分に違いが見られる。これは、i行j列を表すイジング変数$x_{ij}$を$x_{i \times N + j}$のようにして、一次元で表しているためである。もちろん、行と列を入れ替えて$x_{j \times N + i}$としてもよく、隣接行列も入れ替わる。

埋め込み

$N=5$から$N=8$まで場合で、それぞれの定式化に対して必要な物理量子ビットの数を以下に示す。D-Waveが提供している埋め込み関数はヒューリスティックなアルゴリズムを用いているため、ここでは10回の平均値を算出した。

表: 埋め込みに必要な物理量子ビットの数(10回の平均値)

N5678
式(4-1)による定式化の場合2515319881710
式(4-2),(4-3)による定式化の場合2044539171557

各項のパラメータ

D-Wave 2000Qをソルバーとし、SAPI(D-WaveのAPI)を用いてそれぞれの定式化における解を求めた。アニーリング回数は100回と1000回で行なった。また、各項のパラメータを変えながら基底状態が多く得られる場合を探し、以下のように設定した。ただし、今回試した中で最も良かった値を設定しており、最適な値とは限らない。

表: 設定した各項のパラメータ

各項のパラメータ$\lambda_{rowcol}$$\lambda_{diagonal}$$\lambda_{queens}$
式(4-1)の場合221
式(4-2)の場合1.41.21
式(4-3)の場合1.71

結果

成功例・失敗例

ひとまず、得られた結果を見てみよう。式(4-2)の制約項による場合において100回アニーリングを行なった場合の、基本解(a)と(b)、失敗例(c)を示す。失敗例(c)は、図中8行目にクイーンが配置されておらず、6行目に2個のクイーンが配置されているため、行に関する制約を満たしていない。

8×8マスの基本解(a),(b)と失敗例(c)
8×8マスの基本解(a),(b)と失敗例(c)

また、得られた結果の中にバリエーション解が存在した。上記の図(a)のバリエーション解を以下に示す。図(a)の反転と90度の回転によって図(d)が得られ、図(e)がそのバリエーション解であることが確認できる。

8×8マスのバリエーション解
8×8マスのバリエーション解

いくつか解は得られたのだが、そもそも量子アニーリングに基づいているD-Waveマシンは組合せ最適化問題の近似解法であり、この問題のすべての基本解やバリエーション解を厳密に得ることを目的としている訳ではない。上記の表にあるような膨大な可能性の中から、古典コンピュータにおけるアルゴリズムよりも速く最適解を得られることを期待し、D-Waveマシンを使用している。

エネルギーの分布

先ほど示したように、D-Waveマシンから返ってくる結果にはすべての制約を満たすものとそうでないものが含まれる。また、それぞれの解に対するエネルギーのや同じエネルギーの解がいくつ含まれるかといった、様々な情報が得られる。
そこで、得られた結果のエネルギー分布を見てみよう。エネルギーの値を見ることによっても、各制約を満たしているかどうかが分かる。以下の図は、D-Waveマシンから直接返されたエネルギーのヒストグラムである。式(4-3)の場合を除き、アニーリング回数を増やすことでエネルギーの低い度数が多くなっており、エネルギーのより低い解を得られていることが分かる。

式(4-1)の場合のエネルギー分布
式(4-1)の場合のエネルギー分布
式(4-2)の場合のエネルギー分布
式(4-2)の場合のエネルギー分布
式(4-3)の場合のエネルギー分布
式(4-3)の場合のエネルギー分布

得られた解、つまり$x_{ij}$の値を目的関数に代入することで、目的関数のエネルギーを計算することができる。ちなみに先ほど示したエネルギーは、D-Waveマシン上での埋め込みに対するエネルギーである。
以下の図は、それぞれの定式化における目的関数のエネルギーのヒストグラムを示している。なお、基底状態のエネルギーは、式(4-1)と式(4-3)のときに-6.5、式(4-2)のときに-1.5となる(各図における左端のエネルギー)。

式(4-1)の場合のエネルギー分布
式(4-1)の場合の目的関数のエネルギー分布
式(4-2)の場合のエネルギー分布
式(4-2)の場合の目的関数のエネルギー分布
式(4-3)の場合のエネルギー分布
式(4-3)の場合の目的関数のエネルギー分布

式(4-2)の場合ではエネルギーが様々な値をとっている一方で、式(4-1)と式(4-3)の場合ではエネルギーが特定の値にまとまっている。これはクイーンの総数に関する制約について、式(4-1)より式(4-2)の方が緩い制約であることによると考えられる。ここでは示していないが、異なるパラメータの設定で実験した場合では、同じエネルギーの値に対して異なる割合で分布していた。

これまで示した結果の中には、2つ以上の異なった解が同じエネルギーをとることがあった。物理学において、ハミルトニアン(先ほどの目的関数のこと)が同じエネルギー準位をとることを縮退(degeneracy)という。物理的状態が縮退している場合、その状態は対称性をもっていることが多く、N-クイーン問題においては基本解に加えて対称性によるバリエーション解がある。そこで、上の図のようなエネルギー分布に対して、クイーンがどのような配置となっているのか気になってくる。

基底状態のエネルギーと解

前節の結果において、基底状態のエネルギー分布に対するクイーンの配置を以下に示す。式(4-1)の場合、得られた解に多様さがほとんどなく、アニーリング回数を増やすとそれが顕著に表れた。

式(4-2)の場合における解とエネルギーの割合(100回アニーリングした場合)

式(4-1)の場合における解とエネルギーの割合(100回アニーリングした場合)
式(4-1)の場合における解とエネルギーの割合(1000回アニーリングした場合)

式(4-1)の場合における解とエネルギーの割合(1000回アニーリングした場合)

式(4-2)の場合、式(4-1)の場合よりも多種類の解が得られた。またアニーリング回数を増やすことでも、解の種類を増えているることが分かる。
アニーリング回数が100のとき、解(a)の方が解(b)よりも多く得られたが、アニーリング回数を1000にすると、解(b)の方が解(a)よりも多く得られた。アニーリング回数によって得られる各解の個数が変化することは興味深い。

式(4-2)の場合における解とエネルギーの割合(100回アニーリングした場合)

式(4-2)の場合における解とエネルギーの割合(100回アニーリングした場合)
式(4-2)の場合における解とエネルギーの割合(1000回アニーリングした場合)

式(4-2)の場合における解とエネルギーの割合(1000回アニーリングした場合)

式(4-3)の場合、式(4-2)よりもさらに多種類の解が得られた。同様にアニーリング回数を増やすことで、得られる解のバリエーションが増えている。そして、割合の小さい解は割愛しているが、得られた解の中で様々なバリエーション解が得られていた(少々分かりづらいが、下の図の中にもバリエーション解がある)。

式(4-2)の場合における解とエネルギーの割合(100回アニーリングした場合)

式(4-3)の場合における解とエネルギーの割合(100回アニーリングした場合)
式(4-2)の場合における解とエネルギーの割合(1000回アニーリングした場合)

式(4-3)の場合における解とエネルギーの割合(1000回アニーリングした場合)

ここまで、単に解にのみ注目するのではなく、そのエネルギー分布、エネルギーと解の関係を見てきた。
以降はおまけとして、問題の特徴を考慮する、そしてNがより大きな問題を解くことを試みる。

対称性を考慮する

行と列や各向きの斜めといった対称性を利用することで、得られる結果がどうなるか確かめる。具体的には、1個のクイーンの位置を固定する。固定する位置もいくつか考えられるが、ここでは1行2列目にクイーンを配置するようにする。このために横磁場項hに予め設定しておく。予めhの1行2列目だけ-5(それ以外は0)としてから、QUBO行列を作成した。式(4-2)の場合で100回アニーリングを行なったところ、1種類の解のみ得られた。また実験の条件を変えることで、クイーンの位置を固定した場合でも異なる解が得られるだろう。

1行2列目にクイーンを固定した場合の解

1行2列目にクイーンを固定した場合の解

このとき、目的関数のエネルギー分布は以下のようになった。固定しない場合と比べて、基底状態に達する割合が大きくなっていることが分かる。

クイーンの位置を固定した場合の目的関数のエネルギー分布

クイーンの位置を固定した場合の目的関数のエネルギー分布

上記のように、いくつかのクイーンを予め配置した上で残りのクイーンの配置を決定する問題(N-queens Complesion Problem)も考えられている。セント・アンドルーズ大学のIan Gent教授は、彼の論文(“Complexity of n-Queens Completion” byIan Gent, et al.)の中で、この問題とN-クイーン問題の困難さについて述べている。

qbsolvを用いてより大きな問題を解く

D-Wave社が開発しているソルバーqbsolvといった、大きな問題の分割を行なうソルバーを用いることでより大きな問題に対処できる。現在、qbsolvはD-Wave Ocean SDK(githubドキュメント)の一部として提供されている。
$N=12$において、目的関数とパラメータは上記の式(4-2)の場合とした。以下に得られた解の一部を示す。

qbsolvを用いた場合の12×12マスの解

qbsolvを用いた場合の12×12マスの解

応用例

ここでは、N-クイーン問題の応用例を取り上げる。IBMが開発している量子ゲート方式の量子コンピュータを用いたN-クイーン問題に関する論文(“A Novel Quantum N-Queens Solver Algorithm and its Simulation and Application to Satellite Communication Using IBM Quantum Experience” by Rounak Jha, et al.)では、$N=4$の場合におけるシミュレーションと人工衛星への応用が述べられている。
クイーンと人工衛星の反射器、チェス盤と人工衛星が周回する平面が対応し、情報をビームとしてやり取りが行われる。送信される情報の量と、反射器のために送信される地上面積を最大にすることを目的とする。2つのビームが同じ行や列、対角線に沿って伝送される場合、時間とともにランダムに変化する位相差を持つ時点で連続的に干渉する。それが信号波の重ね合わせにつながり、両方の反射器によって送られる情報が失われてしまう。N-クイーン問題の解のように反射器を配置することで、いずれのビームも同じ行、列または対角線に沿って情報を送信せず損失を抑えることができる。

まとめ

まとめと感想

N-クイーン問題の概要を紹介した後、異なる制約項の与え方に対して得られる解とエネルギーについて実験を行なった。制約項の与え方によって、エネルギー分布が異なり、基底状態の解にも違いがみられた。またアニーリング回数の変化や実験を行なうごとに、得られる解、つまりクイーンの配置が変わってくる。これは、ヒューリスティックな埋め込み、各アニーリングによると考えられる。

多くの問題において制約項の与え方はいくつか考えられ、どのような式の場合に望ましい解が多く得られるのか、その与え方によってどのような差異がみられるのかといった疑問を以前から感じていた。解くことができる問題の大きさは未いまだ小さいものの、この疑問に取り組む上での題材としてとても面白い問題であった。現在は量子コンピュータの黎明期であり、扱える問題の大きさや量子ビットのアーキテクチャによる制約を考慮して、問題に取り組むことが重要となる。本記事のように、実際にD-Waveマシンを用いて問題を解くことで、最適化問題のソルバーとしてD-Waveマシンを扱う上での知見を深めていきたい。

今後行いたいこと

  • 得られた結果をさらに分析し、基本解、バリエーション解について調べる
  • 各項のパラメータの最適な値を、効率的に探索できないかサーベイする、また実験を行う
  • D-Waveマシンの次世代機が登場したら、同様の実験を行う

参考文献・サイト

本記事の担当者

丸山尚貴