Help us understand the problem. What is going on with this article?

量子アニーリングで組合せ最適化

More than 3 years have passed since last update.

これはNextremer Advent Calendarの23日目の記事です。

現在、株式会社Nextremerでは早稲田大学の田中宗様と量子アニーリングに関する共同研究を行っております。

はじめに

量子アニーリングは組合せ最適化問題を解くための有効なアルゴリズムと言われています。

本記事では、量子モンテカルロ法による量子アニーリングを用いて、代表的な組合せ最適化問題であるTSP(巡回セールスマン問題)に対して解を求めたいと思います。

量子論の基礎

TSPの議論に入る前に、少しだけ量子論の基礎的なことに触れておきます。

物理量(エネルギーや運動量など)は、量子の世界では演算子というものに対応します。例えば、エネルギーという物理量に対応するものはハミルトニアン($\hat{H}$と表します)という演算子です。

このハミルトニアン $\hat{H}$ という演算子は、ある量子状態 $\psi$ に作用すると、

\hat{H}\psi = E\psi

となり、エネルギー $E$ が固有値という形で目に見える(観測できる)という構造を持っています。

このように、演算子は、量子の世界を知るための"呪文"とも言えます。
と同時に、演算子は量子論を記述する上で不可欠なツールです。

TSPの定式化

TSPは、n都市の集合と各2都市間の距離が与えられた時、全都市を1度ずつ経由して出発地に戻るときの総距離が最小となる経路を求める問題です。

スクリーンショット 2017-01-04 14.25.10.png

今 $n$ 個の都市 $1,2,\cdots,n$ があるとします。

都市 $i,j$ 間の距離を $d_{ij}$、 都市$1$を出発して $a$ ステップ目に都市 $i$ にいる場合に$1$, いない場合に $0$ を表す変数 $n_{i,a}$ を用いると、総距離 $L$ は、

L = \sum_{i,j,a=1}^{n}d_{i,j} n_{i,a} n_{j,a+1} \tag{1}

で表せます。(出発地にいるステップを1ステップ目としています。)

ただし、1都市には1度しか止まれないという条件と、巡回するという条件から、

\sum_{i=1}^{n}n_{i,a}=\sum_{a=1}^{n}n_{i,a}=1 \\
n_{i,n+1}=n_{i,1}  \ \ \ \ \ \ (1\leq i \leq n)

を満たします。$n_{i,a}$ を図で表すと、

スクリーンショット 2016-12-23 18.54.56.png

という感じで、0と1の組合せで表されます。

今回の目的は、この $L$ を最小にする $n_{i,a}$ の組を求めることです。

イジングモデル

電子はスピンという物理量を持っています。先ほども述べましたが、物理量は量子論における演算子というものに対応していますから、スピンに対応する演算子が存在し、Pauli行列と呼ばれるものです。

Pauli行列には3つあり、どれを用いても構わないですが、一般的には $\hat{\sigma}^{z}$ (z方向のスピン角運動量を表す演算子)を用いて、2つの固有値$1,-1$ を観測値として得ます。

スクリーンショット 2017-01-04 14.03.48.png

スピンを用いれば2値の変数を表現できるため、組合せ最適化における2値変数として利用します。

また、実際に組合せ最適化を解くためには、2値の変数をたくさん用意しなければなりません。

そこで、スピンを格子点上に並べたモデル(イジングモデル)を用います。

スクリーンショット 2017-01-04 14.14.40.png

このモデルにおける系のハミルトニアンは、

\hat{H}=-\sum_{i,j}J_{i,j}\hat{\sigma}^{z}_{i}\hat{\sigma}^{z}_{j} \tag{2}

と記述されます。($J_{i,j}$はスピン間の相互作用の強さ、$\hat{\sigma}_{i}$は格子点 $i$ にあるスピンの変数を表します)

よって、このハミルトニアンに対応するエネルギーは

H = -\sum_{i,j}J_{i,j}\sigma^{z}_{i}\sigma^{z}_{j}

と表せます。($\sigma_{i}^{z}$はスピンの変数で、$1,-1$を表す)

よく見ると先の(1)式に酷似していますね。

よって、TSPをうまくイジングモデルに変換すれば、TSPで最適な組合せ $n_{i,a}$ を求めることと、イジングモデルにおいてエネルギーが最小となるスピンの変数の組 $\sigma_{i,a}$ を求めることとが等価になることがわかります。

TSPをイジングモデルに変換

TSPで求めたいものは $n_{i,a}$ の組ですが、これは $0,1$ をとり、イジングモデルにおけるスピンの変数 $1,-1$ に変換する必要があります。

そこで、 $n_{i,a}$ を

n_{i,a}\equiv\frac{\sigma_{i,a}+1}{2}

とおくと、$(1)$式は、

L=\frac{1}{4}\sum_{a,i,j=1}^{n} d_{i,j} \sigma_{i,a}\sigma_{j,a+1}+(const.) \tag{3}

と表せます。(ただし、第1項の定数倍($\frac{1}{4}$)と第2項の定数は、最適化には影響しない項なので無視します。)

よって、Lが$1,-1$ の2値をとる $\sigma_{i,a}$ によって書き直せました。

また、このLを物理量として得るためのハミルトニアンは、

\hat{H}=\sum_{a,i,j=1}^{n} d_{i,j} \hat{\sigma}^{z}_{i,a}\hat{\sigma}^{z}_{j,a+1} \tag{4}

となり、2次元イジングモデルと等価なハミルトニアンであることがわかります。

横磁場

(4)の固有値として得られるエネルギーが最小となる $\sigma_{i,a}$ の組を求めればいいという問題設定に落とし込めました。

ここで、いよいよ量子アニーリングの手法を使っていきます。

(3)のままだといわゆるシミュレーテッドアニーリングの手法で求められますが、
ここに量子効果として、横磁場項 $\color{red}{-\Gamma\sum_{i,a=1}^{n}\hat{\sigma}_{i,a}^{x}}$ を加えます。

\hat{H}=\sum_{a,i,j=1}^{n} d_{i,j} \hat{\sigma}^{z}_{i,a}\hat{\sigma}^{z}_{j,a+1}\color{red}{-\Gamma\sum_{i,a=1}^{n}\hat{\sigma}_{i,a}^{x}} \tag{5}

$\hat{\sigma}_{i,a}^{x}$はPauli行列のx成分、$\Gamma$はアニーリング係数と呼ばれるもので、後で再び出てきます。

分配関数

(5)を使って量子モンテカルロ法に必要な、分配関数 $Z$ というものを求めましょう。その定義は、

Z\equiv tr e^{-\beta\hat{H}}=\sum_{\sigma}<\sigma|e^{-\beta\hat{H}}|\sigma>

で与えられます。($\beta$は温度の逆数、$tr$は全てのスピンの量子状態に関して和を取ります。)

今、横磁場の項があるため非対角化項が出てきてしまい、素直に計算できません。

そこで、鈴木トロッタ分解という手法を用いて、分配関数を近似してやります。
十分大きな $m$ を用いて、

Z=tre^{-\beta\hat{H}}\simeq tr(e^{-\frac{\beta}{m}\hat{H}_{c}}e^{-\frac{\beta}{m}\hat{H}_{q}})^m

とします。($\hat{H}_c,\hat{H}_q$は(5)式の第1項、第2項を表します)

それに伴い、スピン系$|\sigma>$の次元にも新たな次元(トロッタ次元)が加わって、$n\times n\times m$個の3次元スピン系$|\sigma_{i,a,k}>$ (kはトロッタ次元におけるスピンの番号)を考えます。

ここから少々長い計算を行えば、分配関数は、

Z \simeq \biggl[\frac{1}{2}sinh\biggl(\frac{2\beta\Gamma}{m}\biggr)\biggr]^{\frac{mn^2}{2}}\sum_{\sigma_{i,a,k}^z=\pm 1}exp\biggl[-\frac{\beta}{m} \sum_{a=1}^{n}\sum_{i,j=1}^{n}\sum_{k=1}^{m}d_{i,j}\sigma_{i,a,k}^{z}\sigma_{j,a+1,k}^{z} +\frac{1}{2}\log{coth\biggl(\frac{\beta \Gamma}{m}\biggr)}\sum_{i,a=1}^{n}\sum_{k=1}^{m}\sigma_{i,a,k}^{z}\sigma_{i,a,k+1}^{z}\biggr] \tag{6}

と得られます。(この計算結果があってるかどうかは保証できませんので、ご注意ください。)

量子モンテカルロシミュレーション

(6)式をよく見ると、次のエネルギー関数を持つ、3次元のイジングモデルの分配関数と等価であることがわかります。

E = \frac{1}{m} \sum_{a=1}^{n}\sum_{i,j=1}^{n}\sum_{k=1}^{m}d_{i,j}\sigma_{i,a,k}^{z}\sigma_{j,a+1,k}^{z} -\frac{1}{2\beta}\log{coth\biggl(\frac{\beta \Gamma}{m}\biggr)}\sum_{i,a=1}^{n}\sum_{k=1}^{m}\sigma_{i,a,k}^{z}\sigma_{i,a,k+1}^{z} \tag{7}

ただし、$\sigma_{i,a,k}^{z}\in(1,-1)$。
この3次元イジングモデルに対して量子モンテカルロ法を用いた量子アニーリングのアルゴリズムを適用します。

そのアルゴリズムは、

  1. 各トロッタ次元において、各行と各列が 1 が 1 個だけで他は全て-1 となる条件を満たすように、1,-1をランダムに初期化する。
  2. ランダムにトロッタ次元 c とステップ番号 $a 、 b (a\neq b, a,b\geq2)$を選び、aとbの各行について要素が 1 となる列番号をそれぞれp と q とする。行a,bと列p,qが交差する位置にある 4 つの spin を flip させる前後のエネルギー差を∆E とすると、確率 $min(1,e^{-\beta \Delta E})$で flip を受容する。
  3. 2.のステップをモンテカルロステップ回繰り返した後、アニーリングパラメータ $\Gamma$を小さくする。
  4. 2.と3.のステップを繰り返し、 $\Gamma$ を0に近づける。

ここで、$\Delta{E}$は、

\Delta E_{1}=\frac{2}{m}\sum_{i=1}^{n}(\sigma_{i,a-1,c}^{z}+\sigma_{i,a+1,c}^{z}-\sigma_{i,b-1,c}^{z}-\sigma_{i,b+1,c}^{z})(d_{q,i}-d_{p,i}) \\
\Delta E_{2}=\frac{1}{\beta}\log{coth\biggl(\frac{\beta \Gamma}{m}\biggr)}\{\sigma_{p,a,c}^{z}(\sigma_{p,a,c-1}^{z}+\sigma_{p,a,c+1}^{z})+\sigma_{q,a,c}^{z}(\sigma_{q,a,c-1}^{z}+\sigma_{q,a,c+1}^{z})\} \\
\Delta E_{3}=\frac{1}{\beta}\log{coth\biggl(\frac{\beta \Gamma}{m}\biggr)}\{\sigma_{p,b,c}^{z}(\sigma_{p,b,c-1}^{z}+\sigma_{p,b,c+1}^{z})+\sigma_{q,b,c}^{z}(\sigma_{q,b,c-1}^{z}+\sigma_{q,b,c+1}^{z})\}

を用いて、

\Delta E=\Delta E_{1}+\Delta E_{2}+\Delta E_{3}

と表せます。
このアルゴリズムに従い、元のコスト関数(1)によるコスト値が最小となるトロッタ次元の配位が、求めたい組合せです。

あとは、-1,1の変数を0,1の変数に直してやると、元の欲しい最適な組合せになり、当初の目的が達成されます。

実行

ベンチマークとして、こちらのサイトにある、" Djibouti "(ジブチ共和国)の都市データを用いました。

シミュレーションに使ったパラメータは以下です。

  • 初期逆温度 $\beta$ : 37
  • トロッタ数 $m$ : 10
  • アニーリングパラメータの初期値 $\Gamma_{init}$ : 1.0
  • 量子モンテカルロステップ数:13320
  • アニーリングパラメータの減衰率:0.99

(愚直な)コードは、こんな感じ。

QuantumAnnealing.py
#coding:utf-8
import time
import math
import numpy as np
import os
import random
import matplotlib.pyplot as plt

# パラメータの入力
TROTTER_DIM = int(input("Trotter dimension: "))
ANN_PARA =  float(input("initial annealing parameter: "))
ANN_STEP = int(input("Annealing Step: "))
MC_STEP = int(input("MC step: "))
BETA = float(input("inverse Temperature: "))
REDUC_PARA = 0.99

"""
都市についてのデータ(都市の番号とx,y座標)を取得
"""
FILE_NAME = 'FILE_NAME'

f = open(os.path.dirname(os.path.abspath(FILE_NAME))+FILE_NAME).read().split("\n")

POINT = []
for i in f:
    POINT.append(i.split(" "))

# 都市データ
NCITY = len(POINT)
TOTAL_TIME = NCITY
for i in range(NCITY):
    POINT[i].remove(POINT[i][0])
for i in range(NCITY):
    for j in range(2):
        POINT[i][j] = float(POINT[i][j])

## 2都市間の距離
def distance(point1, point2):
    return math.sqrt((point1[1]-point2[1])**2 + (point1[0]-point2[0])**2)


## 全体のspinの配位
## spinの各次元数は[TROTTER_DIM, TOTAL_TIME, NCITY]で表す
def getSpinConfig():

    ## あるトロッタ次元の、ある時刻におけるspinの配位
    def spin_config_at_a_time_in_a_TROTTER_DIM(tag):
        config = list(-np.ones(NCITY, dtype = np.int))
        config[tag] = 1
        return config


    ## あるトロッタ次元におけるspinの配位
    def spin_config_in_a_TROTTER_DIM(tag):
        spin = []
        spin.append(config_at_init_time)
        for i in xrange(TOTAL_TIME-1):
            spin.append(list(spin_config_at_a_time_in_a_TROTTER_DIM(tag[i])))
        return spin

    spin = []
    for i in xrange(TROTTER_DIM):
        tag = np.arange(1,NCITY)
        np.random.shuffle(tag)
        spin.append(spin_config_in_a_TROTTER_DIM(tag)) 
    return spin


# 最短距離となるTrotter次元を選び、その時のルートを出力
def getBestRoute(config):   
    length = []
    for i in xrange(TROTTER_DIM):
        route = []
        for j in xrange(TOTAL_TIME):
            route.append(config[i][j].index(1))
        length.append(getTotaldistance(route))

    min_Tro_dim = np.argmin(length)
    Best_Route = []
    for i in xrange(TOTAL_TIME):
        Best_Route.append(config[min_Tro_dim][i].index(1))
    return Best_Route


## あるルートに対するmax値で割った総距離
def getTotaldistance(route):
    Total_distance = 0
    for i in xrange(TOTAL_TIME):
        Total_distance += distance(POINT[route[i]],POINT[route[(i+1)%TOTAL_TIME]])/max_distance
    return Total_distance


## あるルートに対する総距離
def getRealTotaldistance(route):
    Total_distance = 0
    for i in xrange(TOTAL_TIME):
        Total_distance += distance(POINT[route[i]], POINT[route[(i+1)%TOTAL_TIME]])
    return Total_distance


## 量子モンテカルロステップ
def QMC_move(config, ann_para):
    # 異なる2つの時刻a,bを選ぶ
    c = np.random.randint(0,TROTTER_DIM)
    a_ = range(1,TOTAL_TIME)
    a = np.random.choice(a_)
    a_.remove(a)
    b = np.random.choice(a_)

    # あるトロッタ数cで、時刻a,bにいる都市p,q
    p = config[c][a].index(1)
    q = config[c][b].index(1)

    # エネルギー差の値を初期化
    delta_cost = delta_costc = delta_costq_1 = delta_costq_2 = delta_costq_3 = delta_costq_4 = 0

    # (7)式の第一項についてspinをフリップする前後のエネルギーの差
    for j in range(NCITY):
        l_p_j = distance(POINT[p], POINT[j])/max_distance
        l_q_j = distance(POINT[q], POINT[j])/max_distance
        delta_costc += 2*(-l_p_j*config[c][a][p] - l_q_j*config[c][a][q])*(config[c][a-1][j]+config[c][(a+1)%TOTAL_TIME][j])+2*(-l_p_j*config[c][b][p] - l_q_j*config[c][b][q])*(config[c][b-1][j]+config[c][(b+1)%TOTAL_TIME][j])

    # (7)式の第二項についてspinをフリップする前後のエネルギー差
    para = (1/BETA)*math.log(math.cosh(BETA*ann_para/TROTTER_DIM)/math.sinh(BETA*ann_para/TROTTER_DIM))
    delta_costq_1 = config[c][a][p]*(config[(c-1)%TROTTER_DIM][a][p]+config[(c+1)%TROTTER_DIM][a][p])
    delta_costq_2 = config[c][a][q]*(config[(c-1)%TROTTER_DIM][a][q]+config[(c+1)%TROTTER_DIM][a][q])
    delta_costq_3 = config[c][b][p]*(config[(c-1)%TROTTER_DIM][b][p]+config[(c+1)%TROTTER_DIM][b][p])
    delta_costq_4 = config[c][b][q]*(config[(c-1)%TROTTER_DIM][b][q]+config[(c+1)%TROTTER_DIM][b][q])

    # (7)式についてspinをフリップする前後のエネルギー差
    delta_cost = delta_costc/TROTTER_DIM+para*(delta_costq_1+delta_costq_2+delta_costq_3+delta_costq_4)

    # 確率min(1,exp(-BETA*delta_cost))でflip
    if delta_cost <= 0:
        config[c][a][p]*=-1
        config[c][a][q]*=-1
        config[c][b][p]*=-1
        config[c][b][q]*=-1
    elif np.random.random() < np.exp(-BETA*delta_cost):
        config[c][a][p]*=-1
        config[c][a][q]*=-1
        config[c][b][p]*=-1
        config[c][b][q]*=-1

    return config


"""
量子アニーリングシミュレーション
"""
if __name__ == '__main__':

    # 2都市間の距離の最大値
    max_distance = 0
    for i in range(NCITY):
        for j in range(NCITY):
            if max_distance <= distance(POINT[i], POINT[j]):
                max_distance = distance(POINT[i], POINT[j])


    # 初期時刻におけるspinの配位(初期時刻では必ず都市0にいる)
    config_at_init_time = list(-np.ones(NCITY,dtype=np.int))
    config_at_init_time[0] = 1

    print "start..."
    t0 = time.clock()

    np.random.seed(100)
    spin = getSpinConfig()
    LengthList = []
    for t in range(ANN_STEP):
        for i in range(MC_STEP):
            con = QMC_move(spin, ANN_PARA)
            rou = getBestRoute(con)
            length = getRealTotaldistance(rou)
        LengthList.append(length)
        print "Step:{}, Annealing Parameter:{}, length:{}".format(t+1,ANN_PARA, length)
        ANN_PARA *= REDUC_PARA

    Route = getBestRoute(spin)
    Total_Length = getRealTotaldistance(Route)
    elapsed_time = time.clock()-t0

    print "最短ルートは:{}".format(Route)
    print "最短距離は{}".format(Total_Length)
    print "処理時間は{}s".format(elapsed_time)

    plt.plot(LengthList)
    plt.show()

結果

QAforTSP_1_37.png

縦軸は総距離、横軸はアニーリングステップです。

収束時の解は6737と得られました。
ベンチマークにおける最適解が6656とあるので、比較的良い近似解が得られました。

最後に

今回、TSPのベンチマークを使って、量子モンテカルロ法を用いた量子アニーリングのシミュレーションを行いました。

量子アニーリングで組合せ最適化問題を解けると巷では言われていますが、具体的な情報があまりなかったため、作成しました。何かのお役に立てれば幸いです。

feifake
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away