$$
\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}}
$$

はじめに

若干いまさら感はありますが、XY mixerを使ったQAOAアルゴリズムを一度動かしてみようと思いました。

XY mixerの利点について簡単に確認すると、

ある最適化問題において例えば$\sum_i^N x_i = 1$といったような制約条件があった場合
通常であればコスト関数(ハミルトニアン)に制約項を追加する代わりに
ansatzとしてXY mixerを導入することで
時間発展における探索空間自体を制約を満たす領域内に制限することができる。

と言えます。

今回はこちらの論文で扱っている内容からつまんでいきます。

問題設定

概要

Max-$\kappa$-Colorable-Subgraph problem と呼ばれる問題です。
グラフ彩色の一種で、$\kappa$種類の色でグラフの隣り合う頂点が同じ色にならない塗り分け方が可能な最大のsubgraphを探す問題です。

本記事では論文でも扱われてるプリズムを例とします。
3色で塗り分け可能なことが知られています。

https://arxiv.org/abs/1904.09314

関連する有名な問題としては四色定理、いわゆる地図の塗り分け問題があります。
地図を平面グラフとして扱えばグラフの塗り分け問題になり、平面グラフという条件付きですが4色で塗り分け可能と示されています。

また、地図を3色で塗り分け可能か判定する問題はNP完全問題とされているようです。

定式化

グラフ $G=(V, E)$があり、

頂点数：$V$
エッジ数：$m$
色の種類：$\kappa$
とします。

各頂点はそれぞれ$\kappa$種類の色に塗られる可能性があります。
よって、下記のように変数$x_{v, j}$を導入します。

$$
x_{v,j} = \begin{cases}
1 & (v\ が色\ j\ に塗られている) \\
0 & (その他)
\end{cases}
$$

全てのエッジの頂点$v,v'$について$x_{v,j} x_{v',j}=0$であれば、隣接する頂点が異なる色になります。

コスト関数：

\sum^{\kappa}_{j=1} \sum_{(v,v')\in E}x_{v,j} x_{v',j}

これをハミルトニアンに書き換えると$x \to (I-\sigma^z)/2$の置き換えを経て以下のようになります。

H_C = \kappa m I - \sum^{V}_{v=1}d_v \sum^{\kappa}_{j=1}\sigma^{z}_{v,j} + \sum_{(v,v')\in E}\sum^{\kappa}_{j=1}\sigma^{z}_{v,j}\sigma^{z}_{v',j}

$d_v$はある頂点$v$につながるエッジの数です。
このコストハミルトニアンの期待値を最小化することが目的です。

量子回路

qubitの割当て

頂点$x_{v,j}$それぞれにqubit $q_{v,j}$を割り当てます。$n \times \kappa$ qubitsが必要です。
$[q_{1,1}, q_{1,2}, q_{1,3},q_{2,1}, q_{2,2}, q_{2,3},\cdots,q_{6,3}]=[q_{1}, q_{2}, q_{3},q_{4}, q_{5}, q_{6},\cdots,q_{18}]$
というふうに割り当てます。

初期状態

色の種類$\kappa =3$とすると、各頂点$v$について$q_{v,1}, q_{v,2}, q_{v,3}$のどれか1つのみが$\ket{1}$で他が$\ket{0}$の状態が制約を満たします。

$$
\ket{100}\ \rm{or}\ \ket{010}\ \rm{or}\ \ket{001}
$$
よって初期状態として上の3状態が等しく重ね合わせれた状態が欲しくなります。
このような状態を作る量子回路は論文のAppendix Cで説明されています。
細かな説明は省きますが、N量子ビット中n量子ビットが$\ket{1}$で他が$\ket{0}$である状態の重ね合わせはW-stateと呼ばれており、補助量子ビット1つを用いる量子回路で実現可能なようです。
頂点が6つある場合は、6個のW-stateの重ね合わせが初期状態となります。

XY mixer ansatz

XY mixerは式で書くと以下の通りです。

$$
H_{XY,v,j,j'} = \frac{1}{2}(\sigma^x_{v,j}\sigma^x_{v,j'}+\sigma^y_{v,j}\sigma^y_{v,j'})
$$

$\ket{q_{v,j}q_{v,j'}}$が$\ket{01}$と$\ket{10}$の間でスワップされるような時間発展を実現します。
上記の初期状態に対して
Mixerハミルトニアンは今回以下のようにしました。

$$
\sum_{v} (H_{XY,v,1,2} + H_{XY,v,2,3})
$$

もっと色の種類が増えるような場合にどの程度ミキサーの結合を密にするか、といった議論もあります。
全結合に近くするのか、環状にローテーションするように組むかなどの違いが回路の長さやパフォーマンスに影響します。

実装

初期状態生成

頂点6個、色が3種類で18qubitが必要で、さらに初期状態として6個のW-state生成に補助量子ビットが6個必要なので計24量子ビット系です。
(補助量子ビットを回路途中でリセットしながら繰り返し使えれば1個で済みそうなのですが、現時点でサポートされてる？)

V = 6
K = 3
n = K * V
Xinit = Circuit(n + V)

# ユニタリ変換を定義
def Hamming1(Xinit, anscilla, t, theta):
    Xinit.cx[anscilla, t].x[t].cry(2 * theta)[t, anscilla].x[t].cx[anscilla, t]
    
# 定義したユニタリ変換を適用qubitを変えながら実行
for i in range(V):
    for j in range(K):
        theta = np.arcsin(1 / np.sqrt(K + 1 - (j + 1))) # 角度を計算
        Hamming1(Xinit, n+i, i*K+j, theta)

# 確認用
# Xinit.m[:]
# Xinit.run(shots = 100)

コストハミルトニアン生成

上で記載したコストハミルトニアンの通りに。

d = [3, 3, 3, 3, 3, 3]
Edges = [(0, 1),
        (0, 2),
        (0, 4),
        (1, 2),
        (1, 5),
        (2, 3),
        (3, 4),
        (3, 5),
        (4, 5)]

# Z(18)~Z(23)は本来は不要だが、どうもハミルトニアンの次元と量子回路の次元が合わないと後でエラーを吐くようだったので追加
h = Z(18) + Z(19) + Z(20) + Z(21) + Z(22) + Z(23)

for i in range(V):
    h += -d[i]*(Z(K*i)+Z(K*i+1)+Z(K*i+2))
    
for edge in Edges:
    v1 = edge[0]
    v2 = edge[1]
    for j in range(K):
        h += Z(K*v1+j) * Z(K*v2+j)
    
# 確認用
# print(h)

Mixerハミルトニアン生成

for i in range(6):
    if i == 0:
        XYmixer = 0.5*X[0]*X[1] + 0.5*Y[0]*Y[1] + 0.5*X[1]*X[2] + 0.5*Y[1]*Y[2]
    else:
        XYmixer += 0.5*X[i*3+0]*X[i*3+1] + 0.5*Y[i*3+0]*Y[i*3+1] + 0.5*X[i*3+1]*X[i*3+2] + 0.5*Y[i*3+1]*Y[i*3+2]

# 確認用
# print(XYmixer)

実行

ここがお手軽なのは助かりますね。
補助ビットも含めてハミルトニアン期待値最小化してしまってますけど、含めないでもできるのだろうか...。
初期状態回路生成の段階かVqeに投げる段階かで補助量子ビットを除きたいところです。

補助量子ビット含めた24量子ビットで実行するとかなり時間がかかり、1晩寝て起きたら終わってました。
(むしろ実行不可能なのではと思ってましたがいけました)

step = 2
result = vqe.Vqe(vqe.QaoaAnsatz(h, step, Xinit, XYmixer)).run()

結果の確認

最適化により得られたパラメータによる量子回路実行後の状態を確認します。
測定によって得られる可能性が高い状態の上位12個を見てみましょう。

print(result.most_common(12))

実行結果
最後尾の補助量子ビット6個を除き3qubitごとに"1"がひとつずつ含まれています。
制約は満たしていそうです。

(((0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1), 0.013232081518121674), ((0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1), 0.01323208151812166), ((0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1), 0.013232081518121655), ((1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1), 0.013232081518121651), ((1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1), 0.01323208151812165), ((1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1), 0.013232081518121644), ((0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1), 0.013232081518121644), ((0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1), 0.013232081518121643), ((0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1), 0.013232081518121643), ((0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1), 0.013232081518121639), ((1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1), 0.013232081518121636), ((0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1), 0.013232081518121636))

今度は上位30個についてコスト関数を計算します。

def cost_func(state, Edges):
    cost = 0
    K = 3
    for edge in Edges:
        for j in range(K):
            v1 = edge[0]
            v2 = edge[1]
            cost += state[v1*3+j] * state[v2*3+j]
    return cost

states = result.most_common(30)
for i in range(len(states)):
    state = states[i][0]
    print(state, cost_func(state, Edges))

実行結果
コスト関数0(塗り分け成功)の状態12個が上位に来ていて、コスト関数1の状態がそれに続きます。
最適化される方向には進んでいるので、あとはQAOA step数などあげればコスト関数の小さい状態がより高い確率で得られそうです。

(0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1) 0
(0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1) 0
(0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1) 0
(1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1) 0
(1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1) 0
(1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1) 0
(0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1) 0
(0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1) 0
(0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1) 0
(0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1) 0
(1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1) 0
(0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1) 0
(0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1) 1
(1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1) 1
(1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1) 1
(0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1) 1
(0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1) 1
(0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1) 1
(1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1) 1
(1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1) 1
(0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1) 1
(1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1) 1
(1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1) 1
(0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1) 1
(0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1) 1
(0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1) 1
(1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1) 1
(0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1) 1
(0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1) 1
(1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1) 1

まとめ

XY mixerとそれに合わせた初期条件を適切に生成することで、制約を満たす領域内のみを探索することができました。
初期状態生成用に補助ビットを初期化・再利用したり、最適化時に補助ビットを無視したりできるとさらに効率良く計算できるので今後の課題としたいです。

ただそもそもとして、補助ビット除いても18 qubits使ってやっと6頂点の3色塗り分け問題を解いています。
仮に50 qubits全結合で使っても日本の都道府県さえ塗り分けられないためまだちょっと厳しそうですね。

以上、読んでいただきありがとうございました。