@yuichirominato 2018.07.01更新

D-WaveとWildqatで巡回セールスマン問題とmaxcut問題を解いてみた

D-Wave QUBO wildqat 量子アニーリング

はじめに

さまざまなイジングの問題を解いてきましたが、巡回セールスマンは条件も多くて難しい部類の問題です。プログラムでチャチャっとやってしまいたいところですが、一応順を追って見てみます。

巡回セールスマン問題を解いてみる

巡回セールスマン問題は下記の記事の中に一般式がありますが、よりわかりやすく例題を元に行って見たいと思います。

「NP問題のイジング」
http://blog.mdrft.com/post/42

D-Waveマシンとは?

カナダのベンチャー企業で、量子アニーリングと言われる量子効果を活用したアルゴリズムをハードウェア実装したマシンです。

企業情報:
本社所在地:カナダブリティッシュコロンビア州バーナビー
設立:1999年
事業内容:量子コンピュータのハードウェアとソフトウェアを商用製品として提供
従業員数:160名(Ph.D 55名)


引用元:https://www.dwavesys.com/resources/media-resources


引用元:https://www.dwavesys.com/resources/media-resources

詳細はカナダのホームページを参照してくださいませ。
https://www.dwavesys.com/

日本語サイトもあるようです。
http://dwavejapan.com/

参考

下記の記事など重複する知識が多いので参照します。

「量子アニーリング、イジングモデルとフレームワーク」
http://blog.mdrft.com/post/6

「NP問題のイジング(量子アニーリングなど向け)」
http://blog.mdrft.com/post/42

例題

今回は4都市をAからDを一筆書きで最短距離で行って見ます。なんとなく恣意的な感じもしますが、都市間の距離も設定しておきます。LUCASさんの上記のNP問題のイジングに出る際には結合がない都市間はペナルティを与えていますが、今回は全結合で考えて見ます。

この際、イジングで使用する量子ビット数は16量子ビットになります。かつ、量子ビットはバイナリ値{0,1}で考えます。考え方は、都市を何番目に回るかを16量子ビットで表現をします。下記はその関係を示して見た例です。


横の列に何番目という順番を割り当て、縦の行に都市を割り当てると、上記の表現は、D > B > C > Aの順番で都市を回るということになります。

上記で大事なのは、1都市は一度だけ出現し、1つの順番は一度だけ出現するという制約条件です。

$
H = \sum_{v=1}^n\left( 1-\sum_{j=1}^N x_{v,j} \right)^2 + \sum_{j=1}^N\left(1-\sum_{v=1}^Nx_{v,j} \right)^2 + B\sum_{(u,v)\in E}W_{u,v}\sum_{j=1}^N x_{u,j} x_{v,j}
$

前の2つの項がそれぞれ上記の条件を示しています。これだとよくわからないので具体的にやってみます。

制約条件を考える

制約条件の考え方は簡単です。1都市は一度だけ出現は、横の行で考えると、1行に1量子ビットだけ+1になるという条件が付与されます。たとえば、A1+A2+A3+A4=1です。これをすべての4行で考えます。

同様に、1つの順番は1度だけ出現となると、縦の列の中にも1量子ビットだけ+1になります。たとえばA1+B1+C1+D1=1です。これを4列で考えます。

左上から順番に下記のように量子ビットを割り当てると、

上記条件は、

$
(1-q_0-q_1-q_2-q_3)^2+(1-q_4-q_5-q_6-q_7)^2+(1-q_8-q_9-q_{10}-q_{11})^2+(1-q_{12}-q_{13}-q_{14}-q_{15})^2
$

$
(1-q_0-q_4-q_8-q_{12})^2+(1-q_1-q_5-q_9-q_{13})^2+(1-q_2-q_6-q_{10}-q_{14})^2+(1-q_3-q_7-q_{11}-q_{15})^2
$

足し合わせて展開すると、

$
2q_0q_1 + 2q_0q_{12} + 2q_0q_2 + 2q_0q_3 + 2q_0q_4 + 2q_0q_8 – 2q_0 + 2q_1q_{13} + 2q_1q_2 + 2q_1q_3 + 2q_1q_5 + 2q_1q_9 – 2q_1 + 2q_{10}q_{11} + 2q_{10}q_{14} + 2q_{10}q_2 + 2q_{10}q_6 + 2q_{10}q_8 + 2q_{10}q_9 – 2q_{10} + 2q_{11}q_{15} + 2q_{11}q_3 + 2q_{11}q_7 + 2q_{11}q_8 + 2q_{11}q_9 – 2q_{11} + 2q_{12}q_{13} + 2q_{12}q_{14} + 2q_{12}q_{15} + 2q_{12}q_4 + 2q_{12}q_8 – 2q_{12} + 2q_{13}q_{14}+ 2q_{13}q_{15} + 2q_{13}q_5 + 2q_{13}q_9 – 2q_{13} + 2q_{14}q_{15} + 2q_{14}q_2 + 2q_{14}q_6 – 2q_{14} + 2q_{15}q_3 + 2q_{15}q_7 – 2q_{15} + 2q_2q_3 + 2q_2q_6 – 2q_2 + 2q_3q_7 – 2q_3 + 2q_4q_5 + 2q_4q_6 + 2q_4q_7 + 2q_4q_8 – 2q_4 + 2q_5q_6 + 2q_5q_7 + 2q_5q_9 – 2q_5 + 2q_6q_7 – 2q_6 – 2q_7 + 2q_8q_9 – 2q_8 – 2q_9 + 8
$

距離を考える

距離は都市間の距離が出ていますので、それに対応する量子ビット間にその距離をセットします。例えば、$q_0$と$q_5$の量子ビットは都市Aと都市Bの移動に対応するので、距離「2」を相互作用項にセットします。そのように全ての都市間の距離をセットします。

$
2q_0q_5+q_0q_9+3q_0q_{13} + 2q_4q_1+4q_4q_9+2q_4q_{13} + q_8q_1+4q_8q_5+2q_8q_{13} + 3q_{12}q_1 + 2q_{12}q_5 + 2q_{12}q_9\\
+2q_1q_6+q_1q_{10}+3q_1q_{14} + 2q_5q_2+4q_5q_{10}+2q_5q_{14} + q_9q_2+4q_9q_6+2q_9q_{14} + 3q_{13}q_2+2q_{13}q_6+2q_{13}q_{10}\\
+2q_2q_7+q_2q_{11}+3q_2q_{15} + 2q_6q_3+4q_6q_{11}+2q_6q_{15} + q_{10}q_3+4q_{10}q_7+2q_{10}q_{15} + 3q_{14}q_3+2q_{14}q_7+2q_{14}q_{11}\\
+2q_3q_4+q_3q_8+3q_3q_{12} + 2q_7q_0+4q_7q_8+2q_7q_{12} + q_{11}q_0+4q_{11}q_4+2q_{11}q_{12} + 3q_{15}q_0 + 2q_{15}q_4 + 2q_{15}q_8\\
$

制約条件とコスト関数を組み合わせる

上記制約条件と距離のコスト関数は等しい割合で参入ができません。結合する際にはどちらかに係数をかけてその係数を調整する必要があります。今回は、距離のコスト関数側にBという係数をかけて足し合わせて見ます。

まずWildqcatで

まずWildqatでQUBOを使いながら、B=0.2で制約条件に影響を与えすぎないようにコスト関数を設定しました。

import wildqat as wq
a = wq.opt()

J1 = np.array([ 
[-2,2,2,2,2,0,0,0,2,0,0,0,2,0,0,0], 
[0,-2,2,2,0,2,0,0,0,2,0,0,0,2,0,0], 
[0,0,-2,2,0,0,2,0,0,0,2,0,0,0,2,0], 
[0,0,0,-2,0,0,0,2,0,0,0,2,0,0,0,2], 
[0,0,0,0,-2,2,2,2,2,0,0,0,2,0,0,0], 
[0,0,0,0,0,-2,2,2,0,2,0,0,0,2,0,0], 
[0,0,0,0,0,0,-2,2,0,0,2,0,0,0,2,0], 
[0,0,0,0,0,0,0,-2,0,0,0,2,0,0,0,2], 
[0,0,0,0,0,0,0,0,-2,2,2,2,2,0,0,0], 
[0,0,0,0,0,0,0,0,0,-2,2,2,0,2,0,0], 
[0,0,0,0,0,0,0,0,0,0,-2,2,0,0,2,0], 
[0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,2], 
[0,0,0,0,0,0,0,0,0,0,0,0,-2,2,2,2], 
[0,0,0,0,0,0,0,0,0,0,0,0,0,-2,2,2], 
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,2], 
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2], 
]) 
 
J2 = np.array([ 
[0,0,0,0,0,2,0,2,0,1,0,1,0,3,0,3], 
[0,0,0,0,2,0,2,0,1,0,1,0,3,0,3,0], 
[0,0,0,0,0,2,0,2,0,1,0,1,0,3,0,3], 
[0,0,0,0,2,0,2,0,1,0,1,0,3,0,3,0], 
[0,0,0,0,0,0,0,0,0,4,0,4,0,2,0,2], 
[0,0,0,0,0,0,0,0,4,0,4,0,2,0,2,0], 
[0,0,0,0,0,0,0,0,0,4,0,4,0,2,0,2], 
[0,0,0,0,0,0,0,0,4,0,4,0,2,0,2,0], 
[0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2], 
[0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0], 
[0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2], 
[0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0], 
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 
]) 

B = 0.2
a.qubo = J1+J2*B
a.sa()
[0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0]

これはD>C>A>Bの順に巡回し、コストの高いルートを回避して巡回しているのがわかりました。

D-Waveで実行

続いて、D-Waveの実機での実行を行いますが、
注意すべきポイントは、

1、キメラグラフを考慮する
2、パラメータのレンジをJijを{-1,1}に局所場hiを{-2,2}の範囲に収める。

解けました。

経路はB>D>C>Aでコストの高い順路を避けることができています。
キメラグラフへの実装は完全結合を実行することで、
今回は80量子ビットを使用しました。(もう少し減らすことはできます。)

続いてmaxcut

頂点を2つのグループに分けるような辺の切り方のうち、一番最大数のものを探します。イジング問題で解くときには、隣り合う頂点同士がなるべく異なる符号に落ちるようにエネルギー関数の最小値を探していきます。

引用:https://ja.wikipedia.org/wiki/%E3%82%AB%E3%83%83%E3%83%88_(%E3%82%B0%E3%83%A9%E3%83%95%E7%90%86%E8%AB%96)

今回、ノード間のエッジの重みは1として固定し、maxcutをときます。

maxcut例題

せっかくなので、上記wikipediaにのっている5ノード、6エッジのグラフを解いて見ます。

頂点に今回は量子ビットを割り当てますが、それぞれ0から4の名前をつけます。頂点自体に局所場やバイアスのような値はかけません。ノード間のエッジの重みは1とします。

maxcutのコスト関数一般式は、頂点の量子ビットが{-1,1}をとりうるとして、

$
H = -\sum_{i,j} \frac{1}{2}(1-q_iq_j)
$

ノードが等しいとコストが高くなるようになっています。今回は5ノードあるので、

$
H = -\frac{1}{2}\bigl[(1-q_0q_1)+(1-q_0q_3)+(1-q_1q_2)+(1-q_2q_3)+(1-q_3q_4)+(1-q_2q_4) \bigr]\\
=\frac{1}{2}(q_0q_1+q_0q_3+q_1q_2+q_2q_3+q_3q_4+q_2q_4)-3
$

これはイジング式なので、イジングmatrixに落としてみると、

$
\begin{array}
-&q_0&q_1&q_2&q_3&q_4\\
q_0& 0 & 0.5 & 0 & 0.5 & 0\\
q_1& & 0 & 0.5 & 0 & 0\\
q_2& & & 0 & 0.5 & 0.5\\
q_3& & & & 0 & 0.5\\
q_4& & & & & 0
\end{array}
$

こちらをD-Waveに入力しますがキメラグラフを考慮するとそのままでは入らないので、6量子ビット使います。

こんな感じで、結合をキメラグラフに落とします。q4のクローンを作って対応しました。

maxcutをD-Waveで計算した結果

入力して見た結果、結果として数種類もどってきました。コストが同一なので色々な解があるようです。

1つ見てみると、

確かにきちんととけています。+1と-1の間のエッジをカットすると読めるので、

Wildqatでも

手前味噌ですが、Wildqatでもやってみます。そのままイジングmatrixを入れると計算結果が得られます。

import wildqat as wq                                                                                                                
a = wq.opt()                                                                                                                                                                                                                                               
a.J = [[0,0.5,0,0.5,0],[0,0,0.5,0,0],[0,0,0,0.5,0.5],[0,0,0,0,0.5],[0,0,0,0,0]]                                                     
a.sa()                                                                                                                              
[1, 0, 1, 0, 0]

無事結果が出ました。

あわせて読みたい