はじめに

この記事は Brainpad Advent Calender 2019 23日目の記事になります。
ブレインパッドでは、来たるべき量子コンピュータ時代に向けて量子コンピュータの勉強会を週1回行っています。この勉強会では、量子コンピュータに関連する論文を読んだり、論文を読むために必要な基礎の部分を勉強したりしています。そんな有志の集まりで細々と行っている勉強会ですが、今回はアドベントカレンダーの記事を書くぞ、ということで「現状我々が利用できる量子コンピュータで実用的な問題が解けるのだろうか？」ということにチャレンジすることになりました。

量子コンピュータといえば、通常のコンピュータ(量子コンピュータ界隈では古典コンピュータと呼ばれているやつです)では解くことが困難な問題も解くことができる、と言われているのを聞いたことがある人も多いと思います。実際に一部の問題に対しては古典コンピュータで利用されているアルゴリズムより高速なアルゴリズムが見つかっています。では、現実にそのような問題を解くことができるのか、と言われると、実際できていません。

主な原因としては実現している量子ビットの少なさが理由としてあげられます。ということで、問題が解けるようになる量子ビットが実現するまでじっと待ってましょう！というのも寂しいので、実現している量子ビットで古典コンピュータより高速に問題を解く方法や問題が色々と考えられています。(参考：量子コンピュータの現在とこれから)

そんな流れもあり、我々も「とりあえず使ってみようぜ！」となるべく実用的な問題で、現状でもなんとか量子コンピュータで解けそうな問題を考えました。

議論の結果、今回解く問題は**「秋葉原の特徴の異なるメイド喫茶を効率よく回る」**という問題になりました。
既に目前に迫る2020年、東京オリンピックが開催され、東京には例年以上に外国人観光客が訪れることが予想されます。そういった外国人観光客に向けて観光地のレコメンドをし、Japanese Kawaii Culture に触れてもらうということはとても有益な問題と考えられます。

実用的な問題を量子コンピュータで解く、僕たちが夢見た2020年がここにあることを感じます(ブレードランナーとか2019年の世界ですし)。

データ収集

候補となるメイド喫茶は、ネットのまとめ記事なども参考にしつつ、Google検索した結果から都内にあるレビュー件数100以上の23店舗を集めました。

Googleの検索結果から店舗のクチコミページを開くと、ありがたいことに自動でキーワードタグがレビュー件数とともに付与されています。これらのキーワード情報を店舗情報とともにコピペして、以下のCSV形式のテキストファイルに保存しました。

key,value
店名,キュアメイドカフェ
所在地,東京都千代田区外神田３丁目１５－５ ジーストア・アキバ6F
評価,4
レビュー件数,441
コラボ,31
アニメ,25
紅茶,20
喫茶店,18
萌え,10
メイド服,8
エレベーター,7
画像,7
ラブライブ,6
喫煙,5

このCSVを以下のPythonコードに読み込んで、キーワードの店舗ごとの出現頻度（Term Frequency、TF）と、付与店舗数（Document Frequency、DF）を集計します。

import numpy as np
import pandas as pd
import csv, re
from collections import defaultdict

input_path = 'maid_cafe_tokyo.csv'

cafe_list = []
cafe_addresses = []
cafe_ratings = []
cafe_reviews = []
cafe_TF_dict = {}
cafe_DF_dict = defaultdict(int)

with open(input_path, encoding='cp932') as f:
    reader = csv.reader(f)
    count = 0
    temp_dict = {}
    for row in reader:
        if row[0] == '店名':
            cafe_list.append(row[1])
            count += 1                
        elif row[0] == '所在地':
            cafe_addresses.append(row[1])
        elif row[0] == '評価':
            cafe_ratings.append(row[1])
        elif row[0] == 'レビュー件数':
            cafe_reviews.append(int(re.sub(',', '', row[1])))
        elif row[0] == 'key':
            next
        else:
            # 店舗番号とキーワードをキーとしてTF（クチコミ数）を取得
            cafe_TF_dict[str(count - 1) + '_' + row[0]] = int(row[1])
            # キーワードのDF（付与店舗数）を集計
            cafe_DF_dict[row[0]] += 1
f.close()

類似度の算出

キーワードの店舗ごとの出現頻度（TF）と付与店舗数（DF）が集計できたので、これらを使って店舗ベクトルを構成することにしました。キーワードの DF からキーワード重要度 IDF が計算できるので、いわゆる TFIDF ベクトルとして店舗ベクトルを定義して、それらのコサイン類似度を計算します。類似度計算の Python コードは以下になります。

import math

# 23店舗に付与されてるキーワードをリスト化
word_list = [k for k, v in cafe_DF_dict.items()]

# キーワードのDF値からIDF値を計算
cafe_IDF_dict = {}
for k, v in cafe_DF_dict.items():
    cafe_IDF_dict[k] = math.log(23/v)

# 店舗ベクトルの配列を用意
num_cafes = 23
cafe_vecs = np.zeros([num_cafes, len(word_list)])

# TF値の辞書からキー（店舗番号＋キーワード）を取得
key_list = cafe_TF_dict.keys()

# TF値にIDF値をかけて店舗のTFIDFベクトルを計算
for i in range(num_cafes):
    for j in range(len(word_list)):
        key = str(i) + '_' + word_list[j]
        if key in key_list:
            cafe_vecs[i, j] = cafe_TF_dict[key] * cafe_IDF_dict[word_list[j]]

# ベクトル間のコサイン類似度を関数として定義
def cossim(a, b):
    return np.dot(a, b)/ np.sqrt(np.dot(a, a) * np.dot(b, b))

# 店舗TFIDFベクトルのコサイン類似度を計算
cossim_dist = np.zeros([num_cafes, num_cafes])
const = 0.5

for i in range(num_cafes):
    for j in range(i, num_cafes):
        cossim_dist[i, j] = cossim(cafe_vecs[i], cafe_vecs[j])
        cossim_dist[j, i] = cossim_dist[i, j]

類似度について閾値を設けることで、店舗間の類似度が閾値より大きいペアは似ている、そうでないペアは似ていないとして区別できます。閾値として、23店舗の類似度分布のメジアンを採用しました。こうすれば類似店舗ペアと非類似店舗ペアが同数になるからです。

店舗間距離

巡回セールスマン問題を設定するために、対象となる23店舗を結ぶエッジの重みを定義します。この重みは単純に店舗間の距離を利用します。本来であれば道のりに沿って距離を求めるべきですが、(結構面倒くさいので)今回は直線距離で近似します。距離の算出は以下の手順で行いました。

各店舗の住所を抜き出す
住所から緯度経度を算出
- こちらの記事を参考にさせていただきました：Google スプレッドシートで住所→緯度経度変換
緯度差、経度差から距離を算出
- 計算式は以下の通りです： $$1000 \times \sqrt{\left(\frac{\mbox{緯度}1-\mbox{緯度}2}{0.0111}\right)^2\left(\frac{\mbox{経度}1-\mbox{経度}2}{0.0091}\right)^2} [m]$$
- こちらの計算式を参考にさせていただきました：https://gist.github.com/tilfin/c35d3ee0545b16059860c849da60a114

これによって、スプレッドシート上で一気に距離を計算します。
できあがったものがこんな感じになります。

最適化アプローチ

問題設定

さて、解きたいのは「秋葉原の特徴の異なるメイド喫茶を効率よく回る」という問題でした。これには、以下の2つの思いがありました。

色々な種類のメイド喫茶を知ってもらいたい
ご存知の通り、メイド喫茶は店舗ごとにコンセプトがちがい、個性があります。せっかく、はるばる日本に来られるのですから、色々な個性を体験してもらいたいのです。
できるだけ少ない移動距離で回ってもらいたい
メイドカフェで最高の体験をしてもらっても、カフェ間の移動で疲れてしまっては意味がありません。また、オリンピックで人が増えることを考えると、無駄な移動は渋滞の原因ともなりうるので、避けなければなりません。

そこで、今回はこの問題を以下の2つの部分問題として解くことにしました。

店舗抽出問題
与えられたメイド喫茶リストから、できるだけ個性の違う店舗を抽出する
最短経路問題
与えられたメイド喫茶リストを全て巡る最短の経路をみつける

店舗抽出問題

まず、店舗抽出問題は以下のように定式化しました。

$$
{\rm minimize} \sum_{ij}^{N} S_{ij}x_{i}x_{j} \\
\begin{eqnarray}
S_{ij} &=& \cos(i, j) -C_0, \\
x_i &\in& \{0, -1\}
\end{eqnarray}
$$

ここで$x_{i}$は店舗$i$を抽出するか否かを表すフラグ、$\cos(i, j)$は店舗間のコサイン類似度、$C_0$はしきい値です。普通は一日で巡れる店舗数についての制約を陽に加えることが多いですが、今回はしきい値$C_0$を定義して、それを調整することで抽出される店舗の数を調整することにしました。
これは、一日に巡れる店舗の数が人によって異なっているのと、大雑把に言って$C_0$以上であれば似ていて、$C_0$以下であれば似ていない店舗であるという意味で、解釈がしやすいという点、また量子コンピュータを使うのであればできるだけ問題をシンプルにしておきたいという点を考慮に入れた結果です。

最短経路問題

次に最短経路問題ですが、これはいわゆる巡回セールスマン問題と呼ばれるもので、以下のように定式化できます。

$$
{\rm minimize}\ \sum_{t,i,j}d_{ij}y_{ti}y_{(t+1)j} \\
\begin{eqnarray}
\sum_{t}y_{ti} &=& 1 \\
\sum_{i}y_{ti} &=& 1 \\
y_{ti} &\in& \{0, 1\}
\end{eqnarray}
$$

ここで、$d_{ij}$は店舗$i$と店舗$j$の間の距離、$y_{ti}$は$t$番目に店舗$i$を訪れる場合のみ$1$となるフラグです。定式化の詳細はこちらのブログをご覧ください。

古典的アプローチ

上記の問題は、規模の小さい問題であれば、量子コンピュータを使わずとも0-1整数計画問題として定式化して解くことができます。

抽出問題

前出のままの形だと、2次の項$S_{ij}x_{i}x_{j}$があり、線形でないため多少やっかいです。そこで簡単な変数変換をして2次の項をなくすことを考えます。

まず、新しい変数 $z_{ij} = x_ix_j$ を定義すると、

$$
z_{ij} = \begin{cases}
0 & (x_i = 0 \ {\rm or} \ x_j = 0) \\
1 & {\rm others}
\end{cases}
$$

となります。そこで、$z_{ij} \in \{0, 1\}$が上の条件を満たすような制約を考えてみます。具体的には、

$$
\begin{eqnarray}
z_{ij} &\leq& x_i \\
z_{ij} &\leq& x_j \\
x_i + x_j - 1 &\leq& z_{ij} \\
z_{ij} &\in& \{0, 1\}
\end{eqnarray}
$$

という制約を考えます。$z_{ij} \in \{0, 1\}$という定義を一旦忘れて、この制約を満たす$z_{ij}$を全ての$x_{i}$の組合せについて列挙すると以下のようになります。

$x_i$	$x_j$	$x_i + x_j - 1$	可能な$z_{ij}$の値
0	0	-1	0
0	1	0	0
1	0	0	0
1	1	1	1

上図の可能な$z_{ij}$の値は、まさに $x_ix_j$ と同じです。つまり、$z_{ij}=x_{i}x_{j}$と定義する代わりに、$z_{ij}$に上の制約を科せばよい、ということになります。まとめると、この問題は

$$
\begin{eqnarray}
{\rm minimize} \ H_{\rm ext} &=& \sum_{ij} S_{ij} z_{ij} \\
S_{ij} &=& \cos_(i, j) -C_0 \\
z_{ij} &\leq& x_i \\
z_{ij} &\leq& x_j \\
x_i + x_j - 1 &\leq& z_{ij} \
x_i &\in& \{0, 1\} \\
z_{ij} &\in& \{0, 1\}
\end{eqnarray}
$$

という0-1整数（線形）計画問題に帰着できます。

では、実際に解いてみましょう。

import pulp


median = 0.013

prob = pulp.LpProblem()

x = pulp.LpVariable.dicts('x', range(n_nodes), cat='Binary')
z = pulp.LpVariable.dicts('z', (range(n_nodes), range(n_nodes)), cat='Binary')


for i in range(n_nodes):
    for j in range(i + 1, n_nodes):
        prob += z[i][j] <= x[i]
        prob += z[i][j] <= x[j]
        prob += z[i][j] >= x[i] + x[j] - 1

prob += pulp.lpSum([
    (sim_reduced[i, j] - median)*z[i][j]
    for i in range(n_nodes)
    for j in range(i + 1, n_nodes)
])

%time print(pulp.LpStatus[tsp.solve()])

for i in range(n_nodes):
    v = y[i].value()
    if v > 0:
        print(caffe_list_reduced[i])

問題を解くのにかかった時間は118msでした

Optimal
CPU times: user 6.25 ms, sys: 7.71 ms, total: 14 ms
Wall time: 118 ms
@home cafe ドンキ店（＠ほぉ?むカフェ）
カフェ メイリッシュ
私設図書館 シャッツキステ
めいどりーみんアイドル通り店
AKIBAスモーカーズ
鉄道居酒屋

12店舗のうち、6店舗が選ばれたようです。確かに、ダイバーシティのある選択になっています。

巡回セールスマン問題

つづいて巡回セールスマン問題を解いていきましょう。巡回セールスマンについても、最小化したい関数に2次の項が含まれていました。先程と同様に線形な式に変換することもできますが、今回はいわゆる MTZ 方式の定式化をしました。

$$
{\rm minimize.} \ \sum_{ij}d_{ij}x_{ij} \\
\begin{eqnarray}
\sum_{i}x_{ij} &=& 1 \\
\sum_{j}x_{ij} &=& 1 \\
u_{i} + 1 - N(1 - x_{ij}) &\leq& u_j \quad {\rm for} \ 1 \leq \ i,j \\
\\
u_{i} &\in& \{1, 2, \dots, N \} \\
x &\in& \{0, 1\}
\end{eqnarray}
$$

ここで、$x_{ij}$は、店舗$i$から店舗$j$への経路を通るか否かのフラグ、$u_{i}$ が回る順番を表しています。こちらもさらっと実装できます。

tsp = pulp.LpProblem('tsp')

n_nodes = dist_to_go.shape[0]

x = pulp.LpVariable.dicts('x', (range(n_nodes), range(n_nodes)), cat='Binary')
u = pulp.LpVariable.dicts('u', range(n_nodes), cat='Integer', lowBound=1, upBound=n_nodes - 1)

for j in range(n_nodes):
    tsp += pulp.lpSum(x[i][j] for i in range(n_nodes) if i != j) == 1

for i in range(n_nodes):
    tsp += pulp.lpSum(x[i][j] for j in range(n_nodes) if i != j) == 1

for i in range(n_nodes):
    for j in range(1, n_nodes):
        if i == j:
            continue
        tsp += u[i] + 1 - n_nodes*(1 - x[i][j]) <= u[j]

tsp += pulp.lpSum(dist.values[i, j]*x[i][j] for i in range(n_nodes) for j in range(n_nodes) if i != j)

%time print(pulp.LpStatus[tsp.solve()])

sorted((u[i].value(), caffe_to_go[i]) for i in range(n_nodes))

Optimal
CPU times: user 2.51 ms, sys: 4.94 ms, total: 7.45 ms
Wall time: 52.9 ms
[(1.0, '@home cafe ドンキ店（＠ほぉ?むカフェ）'),
 (2.0, 'カフェ メイリッシュ'),
 (3.0, '私設図書館 シャッツキステ'),
 (4.0, 'AKIBAスモーカーズ'),
 (5.0, 'めいどりーみんアイドル通り店'),
 (6.0, '鉄道居酒屋')]

52.9ms かかったようです。さて、量子コンピュータを使うともっと高速になるのでしょうか？

量子的アプローチ

古典的アプローチでは、最適化問題の目的関数と制約条件は、$x_i \in \{0, 1\}$ といった2値変数により記述されました。ただし、目的関数が2値変数の2次式のままでは 0-1整数計画問題として解くことができず、$x_i x_j$ と等価な別の2値変数 $z_{ij}$ により目的関数や制約条件を書き直す必要がありました。

一方、物理学の問題では、エネルギー関数（ハミルトニアン）が変位や速度や2次式として表される例が多数あります。その代表例として、強磁性体のイジングモデルがあります。このモデルでは、強磁性体は格子状に並んだ原子磁石の相互作用としてエネルギー関数を定義します。ハミルトニアンは、スピン変数 $\sigma^{(i)}$ の2次式として以下のように記述されます。

$$
H = -\sum_{i,j} J_{ij} ,\sigma^{(i)} \sigma^{(j)}
-\sum_i h_i ,\sigma^{(i)}
$$

原子や分子といったミクロな対象は、量子力学で記述されます。量子力学では、すべての物理量は行列（または演算子）であり、その観測値は行列の固有値に限定されます。スピン変数は、2行2列のスピン行列:

\sigma^{(i)} = \left(\begin{array}{cc}
 1 & 0 \\ 0 & -1 \end{array}\right)

として表されるので、その取り得る値は2つの固有ベクトル $\left({0 \atop 1}\right),$ $\left({1 \atop 0}\right)$ に対応する2つの固有値 $-1,$ $1$ に限定されます。

\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}

\sigma^{(i)} \ket{s_i} = s_i\ket{s_i}, \quad s_i \in \{-1,1\}

これは、古典的アプローチの2値変数 $x_i \in \{0,1\}$ と同様に振る舞います。

そこで量子的アプローチでは、古典的アプローチの2値変数をスピン行列に変換して、スピン系ハミルトニアンの基底状態（最低エネルギー状態）として、最適解を見つけます。実際、店舗抽出問題の目的関数で以下の変換を施すと、イジングモデルのハミルトニアンに変換することが確認できます。

\begin{eqnarray}
x_i & \rightarrow & {\sigma^{(i)} + I \over 2}, \quad
I = \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right)\\
S_{ij} & \rightarrow & -4J_{ij}, \quad
\sum_i S_{ij} \rightarrow\, -2 h_j
\end{eqnarray}

量子コンピュータによる解法

イジングハミルトニアンの基底状態を求める方法として、厳密解法と近似解法があります。

厳密解法とは、ハミルトニアンを表現する行列の固有値問題を直接解いて、最小固有値に対応する固有ベクトル（すなわち基底状態）を求める方法です。スピン変数の数を $N$ とすると、ハミルトニアン行列の次元数は $2^N$ となるので、$N$ が100のオーダーより大きな規模になると古典コンピュータでは現実的な時間内で解くことはかないません。

理想的な量子コンピュータを使えば、こうした大規模な問題も現実的な時間で解くことができると期待されています。しかし、現状の技術では理想からほど遠く、以下のような近似解法が提案されています。この近似解法では、スピン行列の2量子状態 $N$ 個分を掛け算（直積）して、ハミルトニアンの $2^N$ 個の量子状態を重ね合わせた状態を用意します。近似解法では、この重ね合わせた量子状態について、その制御パラメータ $\lambda$ を断熱的に更新しながら基底状態に近づけていきます。

\begin{eqnarray}
E(\lambda) &=& \braket{\lambda}{H_{\rm Ising}(\sigma^{(1)},\dots,\sigma^{(N)}) \middle|\lambda} \\
\ket{\lambda} &=& U(\lambda)\ket{s_1}\otimes\dots\otimes\ket{s_N}
\end{eqnarray}

QAOA（Quantum Approximate Optimization Algorithm）という手法では、スピン行列の固有ベクトルを量子ビットとみなして、$N$ 個の量子ビットの基底状態 $\ket{s_1}\otimes\dots\otimes\ket{s_N}$ に対する量子ゲート演算 $U(\lambda)$ により近似的な量子状態 $\ket{\lambda}$ を表現します。この量子状態のもとでイジングハミルトニアンの期待値 $E(\lambda)$ を求めたら、その最小化問題は古典コンピュータで解きます。この処理をエネルギー期待値が最小値に収束するまで繰り返すという方法です。詳しくは、以下の記事を参照してください。

QAOA - ゲート式量子コンピューターで最適化問題を解く近似アルゴリズム

IBM Q と Qiskit について

IBM Q では、複数の量子ビットから構成されるゲート型量子コンピュータを利用できます。ただし、量子コンピュータの実機は、量子ビット数に制限（最大で14個）があるので、古典コンピュータで量子コンピュータの動作を模倣する量子シミュレータも用意されています。

Python で IBM Q を動かすには Qiskit というライブラリが使えます。Max-Cut問題や巡回セールスマン問題を扱ったサンプルコードが Jupyter Notebook で用意されています。今回はこれを参考にしました。

Qiskitチュートリアルのサンプルコード

Python には0-1整数計画問題を解くためのライブラリとして docplex があります。これを使うと、Python コード上で2値変数を定義して、目的関数や制約条件を2値変数による定義式として書くだけで、最適化問題をコード化できます。Qiskit には、docplex で指定した目的関数や制約条件をイジングハミルトニアンに変換してくれるライブラリがあります。今回、これを使うことで、複雑な量子ゲート回路を構成しなくても問題を量子ハミルトニアンにエンコードすることができました。以下のコードです。

from docplex.mp.model import Model
from qiskit.aqua.translators.ising import docplex

# 店舗抽出問題のモデルと変数のインスタンスを生成
mdl = Model(name='extract_cafe')
x = {i: mdl.binary_var(name='x_{0}'.format(i)) for i in range(n)}

# 目的関数
object_func = mdl.sum((sim2[i,j] - 0.026/2) * x[i] * x[j] 
                      for i in range(n) for j in range(i))
mdl.minimize(object_func)

# モデルのインスタンスを量子計算用に変換
qubitOp_docplex, offset_docplex = docplex.get_qubitops(mdl)

対象店舗数が100件以上の大規模な問題では、状態数は $2^{100} \simeq 1.3 \times 10^{30}$という途方もなく大きな数になるので、古典的アプローチはもちろん、固有値問題を古典コンピュータで解くことも事実上不可能です。そこで、量子的アプローチが期待されるわけですが、現状では IBM Q でも最大で14量子ビット（状態数は $2^{14} \simeq 1.6万$）までが限界です。この規模なら、古典的アプローチの方が厳密解を短時間で返すことができてしまいます。

店舗抽出問題の量子シミュレータ解

Qiskit では、固有値問題を直接に解く厳密解法と、量子コンピュータの実機またはシミュレータを使って QAOA を解く近似解法の両方を試すことができます。

# イジングハミルトニアンの固有値問題を直接に解いて基底状態を見つける
ee = ExactEigensolver(qubitOp_docplex, k=1)
result = ee.run()

x = docplex.sample_most_likely(result['eigvecs'][0])
extract_idx = [i for i in range(len(x)) if x[i] > 0]
extract_cafes = [cafe_list2[c] for c in extract_idx]

print('objective:', result['energy'] + offset_docplex)
print('solution:',x)
print('extracted indexes', extract_idx)
print('extracted maid cafes', extract_cafes)

店舗抽出問題については、当初は23店舗を抽出対象としていました。この場合、古典的解法と同様、厳密解法で解くと個性的な8店舗を抽出できました。

しかし、IBM Q の実機は最大でも14量子ビットまでしか計算できません。23店舗からの抽出問題は、23量子ビットが必要となります。そこで、やむなく店舗数を12店舗に減らすことにしました。その際、店舗ベクトルの23店舗で計算した類似度をそのまま使うことにしました。ただし、閾値を類似度分布のメジアンの半分の値としました。その結果、厳密解法では古典計算と同じ6店舗を抽出できました。

seed = 10598

# 最適化エンジンを指定
spsa = SPSA(max_trials=300)
ry = RY(qubitOp_docplex.num_qubits, depth=5, entanglement='linear')

# VQE（量子変分固有値ソルバー）を指定 
vqe = VQE(qubitOp_docplex, ry, spsa)

# 最適化を実行するバックエンドを指定（実機またはシミュレータ）
backend = Aer.get_backend('qasm_simulator')
quantum_instance = QuantumInstance(backend, shots=1024, seed_simulator=seed,
                                   seed_transpiler=seed)

# VQEを実行
result = vqe.run(quantum_instance)

# 近似解を抽出
x = docplex.sample_most_likely(result['eigvecs'][0])
extract_idx = [i for i in range(len(x)) if x[i] > 0]
extract_cafes = [cafe_list2[c] for c in extract_idx]

# 結果の書き出し
print('time:', result['eval_time'])
print('objective:', result['energy'] + offset_docplex)
print('solution:', x)
print('extracted indexes', extract_idx)
print('extracted maid cafes', extract_cafes)

plot_histogram(result['eigvecs'][0])

IBM Q の実機で計算するには、14量子ビットを搭載した ibmq_16_melborne が必要です。今回、この実機にジョブを投げたのですが、エラーが出て受け付けてもらえませんでした。残念ですが量子シミュレータによる計算にとどめることにします。3.5時間にわたる量子シミュレータによる最適化計算の結果、以下の4店舗が抽出されました。

time: 12474.151748418808
objective: -0.037108848625413526
solution: [0 0 0 0 1 1 0 0 0 0 1 1]
extracted indexes [4, 5, 10, 11]
extracted maid cafes ['もののぷ', 'カフェ メイリッシュ', 'AKIBAスモーカーズ', '鉄道居酒屋']

巡回セールスマン問題の量子シミュレータ解

4店舗の巡回セールスマン問題を解くには、$4^2 = 16$個の量子ビットが必要です。IBM Q の実機は、最大でも14量子ビットまでなので、現状では巡回セールスマン問題は、実機では解けません。代わりに量子シミュレータで解くことにしました。

# 巡回セールスマン問題のモデルと変数のインスタンスを生成
mdl = Model(name='tsp')
x = {(i,p): mdl.binary_var(name='x_{0}_{1}'.format(i,p)) 
     for i in range(n) for p in range(n)}

# 目的関数
tsp_func = mdl.sum(dist[i,j] * x[(i,p)] * x[(j,(p+1)%n)] 
                   for i in range(n) for j in range(n) for p in range(n))
mdl.minimize(tsp_func)

# 制約条件
for i in range(n):
    mdl.add_constraint(mdl.sum(x[(i,p)] for p in range(n)) == 1)
for p in range(n):
    mdl.add_constraint(mdl.sum(x[(i,p)] for i in range(n)) == 1)

24時間にわたり計算しましたが、残念なことに制約条件を満たす近似解は得られませんでした。4店舗のうち「もののぷ」を除く3店舗を巡回する経路が得られました。興味深いことに、この近似解は厳密解から「もののぷ」を除いて、逆回りに巡回するパスと同じです。21世紀の量子シミュレータにとって、「もののぷ」のような中世のコスプレは受け入れられないようですね。

# 厳密解
objective: 569.0
solution: [0. 1. 0. 0. 0. 0. 0. 1. 1. 0. 0. 0. 0. 0. 1. 0.]
feasible: True
solution: [2, 0, 3, 1]
-> ['カフェ メイリッシュ','鉄道居酒屋','もののぷ','AKIBAスモーカーズ']
solution objective: 569.0

# 量子シミュレータ
objective: 30731.8
solution: [0. 0. 0. 1. 0. 1. 0. 0. 0. 0. 1. 0. 0. 0. 0. 1.]
feasible: False
solution: [0, 2, 1]
-> ['鉄道居酒屋', 'カフェ メイリッシュ', 'AKIBAスモーカーズ']
solution objective: 569.0

検証

古典と量子の比較

店舗抽出問題では、古典的アプローチが短時間で厳密解に至ったのに対して、量子的アプローチでは3.5時間かかって準最適な解しか得られませんでした。しかし、抽出結果を見ると量子的アプローチで得られた4店舗のうち「もののぷ」を除く4店舗は、古典的アプローチで抽出された6店舗に含まれていました。量子的アプローチはお財布に優しいようです。

お店を巡回してみた

実際に量子シミュレータの結果に従ってメイド喫茶を巡回してきました。経路としては、鉄道居酒屋 → カフェメイリッシュ → AKIBAスモーカーズ を回ることにします(量子コンピュータの指令にないもののぷは除外しました)。

秋葉原駅電気街口 → 鉄道居酒屋

開始地点は秋葉原駅電気街口とし、ここに戻ってくるまで時間を測ります。最初の目的地鉄道居酒屋 LittleTGVまでの経路は徒歩6分です。

そもそも居酒屋と言っているけれどメイド喫茶なのだろうか……？などと疑問を抱きつつ到着です。

名前の通り鉄道をコンセプトとしたお店で、店内には様々な鉄道関係の道具が並んでおり、入店時には切符をもらうなど、とても凝ったお店になっていました。

メイド喫茶なの？という問題ですが、駅員の格好をした女性店員さんが接客をしてくれるため、基本的な部分は同じなのかな？と思いました。

鉄道居酒屋 → カフェメイリッシュ

次にカフェメイリッシュに向かいます。鉄道居酒屋からは徒歩3分です。

すぐに到着しました。なんだか見覚えがあるな、と思っていたら、シュ○ゲに出てくるメイド喫茶のモデルになったお店でした。

こちらのお店はイメージ通りのメイド喫茶、という感じで、入店時に「お帰りなさいませ、旦那様！」を体験しました。

カフェメイリッシュ → AKIBAスモーカーズ

最後にAKIBAスモーカーズに向かいます。シーシャカフェという水タバコを楽しむメイドカフェというコンセプトになっており、メイドカフェ巡回を行っているメンバー3名は誰一人として喫煙者がいないのですが、量子コンピュータの指令に逆らうわけにはいかないため、ちゃんと行きます。カフェメイリッシュからは徒歩2分です。

到着です。割とビビりながら入店します。

店内は薄暗く甘い匂いが漂っており、水タバコを楽しむ場の雰囲気が出ています。

AKIBAスモーカーズ → 秋葉原駅電気街口

きちんと巡回をするために、最後に秋葉原駅電気街口まで戻ります。AKIBAスモーカーズからは徒歩7分です。

巡回にかかった時間を見てみると5時間21分でした！

量子コンピュータの命令に従ってメイド喫茶を巡回した感想は以下の通りです。

似ていないお店の抽出はよくできていると思いました
- イメージ通りのメイド喫茶1店、（メイド喫茶としては）尖り気味の2店という印象です
全体的にメイド喫茶同士が近い
- そもそも秋葉原の中でメイド喫茶の集まっている地域があるそうで、その辺りから抽出されているようです
- 店舗抽出の段階で距離を考慮できると良いかも、と思いました
時間を測る意味があまりない
- それは充分わかっています
- ちなみに、メイド喫茶の多くは時間制限付きだそうです（1時間半くらい）
コンピュータの指令に従って行動するという体験が楽しかったです
- ~~ディストピア~~SF体験ですね