Edited at

IBMのQiskit Aquaで量子コンピュータ向けアルゴリズムを実装(maxcut)


はじめに

量子コンピュータのアプリケーションを1から作るのは大変です。最近では主要なアルゴリズムはライブラリという形で実装済みのものが配布されています。すべて無料で手に入ります。

IBMが提供する量子コンピュータはIBMQと呼ばれ、Qiskitと呼ばれるツールが配布されています。今回はQiskit向けのアルゴリズムライブラリであるQiskit Aquaを使ってみます。今回は概要を確認します。


Qiskit Aqua

QiskitはIBMが提供する無料の量子コンピュータ向けのSDKです。

「IBMの量子コンピュータSDKのQiskitをつかって簡単な回路を作ってみる」

https://qiita.com/YuichiroMinato/items/6882e43eafd2fd167e75

今回はそれに合わせて多くのアルゴリズムを簡単に使えるAquaをみます。

Qiskit Aqua

Algorithms for near-term quantum applications

https://qiskit.org/aqua

near-termと言うのは近年のと言う意味で、最近の量子コンピュータを指します。アルゴリズムには汎用型とハイブリッド型のように用途とエラーの量に応じてちょっと使い方が変わります。Aquaには両方収録されているので勉強になります。


インストール

インストールはpipから行います。python3.5以上対応です。

pip install qiskit-aqua

簡単ですね。ドキュメントなどは、qiskit aquaのページがあります。

https://qiskit.org/aqua


Algorithms

今回はアルゴリズムを見てみます。

https://qiskit.org/documentation/aqua/algorithms.html

2019年8月現在このようなアルゴリズムが実装されています。

Variational Quantum Eigensolver (VQE)

Quantum Approximate Optimization Algorithm (QAOA)

Evolution of Hamiltonian (EOH)

Quantum Phase Estimation (QPE)

Iterative Quantum Phase Estimation (IQPE)

Amplitude Estimation

Quantum Grover Search

Deutsch-Jozsa

Bernstein-Vazirani

Simon

Quantum Support Vector Machine (QSVM)

Variational Quantum Classifier (VQC)

HHL algorithm for solving linear systems (HHL)

Shor’s Factory Algorithm (Shor)

Quantum Generative Adversarial Network(qGAN)

今回はIBMのチュートリアルを参考にしながらちょっと改造して今流行りの最適化問題のmaxcutをqiskitシミュレータで解いてみたいと思います。


maxcut

maxcutを解くためには、問題をイジングモデルと呼ばれる物理モデルに落とし込み、そこからアルゴリズムにかけます。今回は今までやってきたようなXとかHとかの回路の組み合わせを使う必要はありません。

まずはmaxcutの問題を準備します。maxcutは点と点を結んでできる線を一筆書きでできるだけたくさん横切ると言う問題です。分類問題に利用できます。

Figure_1.png

4ノードのmaxcut問題を見てみます。

import numpy as np

import networkx as nx
import matplotlib.pyplot as plt

上記をnetworkxでみてみます。

#4ノード作ります。

n=4

G=nx.Graph()
G.add_nodes_from(np.arange(0,n,1))

#エッジijに対して、(i,j,結合荷重)で設定します。
elist=[(0,1,1.0),(0,2,1.0),(0,3,1.0),(1,2,1.0),(2,3,1.0)]

G.add_weighted_edges_from(elist)
colors = ['r' for node in G.nodes()]
pos = nx.spring_layout(G)
default_axes = plt.axes(frameon=True)
nx.draw_networkx(G, node_color=colors, node_size=600, alpha=.8, ax=default_axes, pos=pos)

これで上記のグラフがかけます。ここはQiskitとは関係ないのでわからない人は詳しくはみなくても大丈夫です。


次にノード間の接続のグラフを作ります

次に上記の繋がりを行列の形で表現します。

上記のグラフの繋がっているところは、結合荷重を入れ、繋がってないところには0を入れます。

w = np.zeros([n,n])

for i in range(n):
for j in range(n):
temp = G.get_edge_data(i,j,default=0)
if temp != 0:
w[i,j] = temp['weight']
print(w)

ここまでくると、行列の形でノードの関係が把握できます。縦横共に0から3番までのノードの関係性が行列かされています。例えば、1行目の3列目は0番目のノードと2番目のノードの結合を表します。

[[0. 1. 1. 1.]

[1. 0. 1. 0.]
[1. 1. 0. 1.]
[1. 0. 1. 0.]]


次にイジングに直す

次に式をイジングに直します。こちらもQiskit Aquaに収録されている簡単にイジングを変換するツールを使います。

from qiskit.aqua.translators.ising import max_cut

from qiskit.aqua.input import EnergyInput

qubitOp, offset = max_cut.get_max_cut_qubitops(w)
algo_input = EnergyInput(qubitOp)

これでできました。念の為中身を見てみます。

vars(qubitOp) 

こちらには接続情報がパウリ演算子の形で入っています。Zゲートを使っています。

{'_paulis': [[0.5,

Pauli(z=[True, True, False, False], x=[False, False, False, False])],
[0.5, Pauli(z=[True, False, True, False], x=[False, False, False, False])],
[0.5, Pauli(z=[False, True, True, False], x=[False, False, False, False])],
[0.5, Pauli(z=[True, False, False, True], x=[False, False, False, False])],
[0.5, Pauli(z=[False, False, True, True], x=[False, False, False, False])]],
'_coloring': 'largest-degree',
'_grouped_paulis': None,
'_matrix': None,
'_dia_matrix': None,
'_paulis_table': {'IIZZ': 0, 'IZIZ': 1, 'IZZI': 2, 'ZIIZ': 3, 'ZZII': 4},
'_summarize_circuits': False}

offsetにはエネルギー値のオフセットが入っていました。上記maxcutをイジングに直すと数字が{-1,1}になるので、変換が必要になります。

vars(algo_input)

こちらは、エネルギーに関する書式が規定されていました。

{'_configuration': {'name': 'EnergyInput',

'description': 'Energy problem input',
'input_schema': {'$schema': 'http://json-schema.org/schema#',
'id': 'energy_state_schema',
'type': 'object',
'properties': {'qubit_op': {'type': 'object', 'default': {}},
'aux_ops': {'type': ['array', 'null'], 'default': None}},
'additionalProperties': False},
'problems': ['energy', 'excited_states', 'eoh', 'ising']},
'_qubit_op': <qiskit.aqua.operator.Operator at 0x10df98b90>,
'_aux_ops': []}

上記のオブジェクトの中身は知らなくてもできますので気にしないようにしましょう。


まずは状態ベクトルで解いてみる

Qiskitには様々なモードがあります。状態ベクトルを使うモードや、実際に何回かサンプリングをして計算をする方法です。まずは状態ベクトルと呼ばれる解の出る確率の平方根に対応する確率振幅を求める方をやってみましょう。

さらにツールを読み込みます。

from qiskit import BasicAer

from qiskit.aqua.algorithms import VQE
from qiskit.aqua.components.optimizers import SPSA
from qiskit.aqua.components.variational_forms import RY
from qiskit.aqua import QuantumInstance

seed = 10598

spsa = SPSA(max_trials=300)
ry = RY(qubitOp.num_qubits, depth=5, entanglement='linear')
vqe = VQE(qubitOp, ry, spsa, 'matrix')

backend = BasicAer.get_backend('statevector_simulator')
quantum_instance = QuantumInstance(backend, seed=seed, seed_transpiler=seed)

result = vqe.run(quantum_instance)

これで計算ができました。早速答えを可視化してみます。計算結果をnetworkxの形に落とし込みます。

x = max_cut.sample_most_likely(result['eigvecs'][0])

print('energy:', result['energy'])
print('time:', result['eval_time'])
print('max-cut objective:', result['energy'] + offset)
print('solution:', max_cut.get_graph_solution(x))
print('solution objective:', max_cut.max_cut_value(x, w))

colors = ['r' if max_cut.get_graph_solution(x)[i] == 0 else 'b' for i in range(n)]
nx.draw_networkx(G, node_color=colors, node_size=600, alpha = .8, pos=pos)

無事書くことができました。

Figure_1.png

青と赤の間のエッジが切られることで最大化されます。


コードを全部まとめると、

import numpy as np

import networkx as nx
import matplotlib.pyplot as plt

from qiskit.aqua.translators.ising import max_cut
from qiskit.aqua.input import EnergyInput

from qiskit import BasicAer
from qiskit.aqua.algorithms import VQE
from qiskit.aqua.components.optimizers import SPSA
from qiskit.aqua.components.variational_forms import RY
from qiskit.aqua import QuantumInstance

n=4

G=nx.Graph()
G.add_nodes_from(np.arange(0,n,1))

#エッジijに対して、(i,j,結合荷重)で設定します。
elist=[(0,1,1.0),(0,2,1.0),(0,3,1.0),(1,2,1.0),(2,3,1.0)]

G.add_weighted_edges_from(elist)
colors = ['r' for node in G.nodes()]
pos = nx.spring_layout(G)
default_axes = plt.axes(frameon=True)
nx.draw_networkx(G, node_color=colors, node_size=600, alpha=.8, ax=default_axes, pos=pos)

w = np.zeros([n,n])
for i in range(n):
for j in range(n):
temp = G.get_edge_data(i,j,default=0)
if temp != 0:
w[i,j] = temp['weight']
print(w)

qubitOp, offset = max_cut.get_max_cut_qubitops(w)
algo_input = EnergyInput(qubitOp)

seed = 10598

spsa = SPSA(max_trials=300)
ry = RY(qubitOp.num_qubits, depth=5, entanglement='linear')
vqe = VQE(qubitOp, ry, spsa, 'matrix')

backend = BasicAer.get_backend('statevector_simulator')
quantum_instance = QuantumInstance(backend, seed=seed, seed_transpiler=seed)

result = vqe.run(quantum_instance)

x = max_cut.sample_most_likely(result['eigvecs'][0])
print('energy:', result['energy'])
print('time:', result['eval_time'])
print('max-cut objective:', result['energy'] + offset)
print('solution:', max_cut.get_graph_solution(x))
print('solution objective:', max_cut.max_cut_value(x, w))

colors = ['r' if max_cut.get_graph_solution(x)[i] == 0 else 'b' for i in range(n)]
nx.draw_networkx(G, node_color=colors, node_size=600, alpha = .8, pos=pos)

これは機能のほんの一部なので、今後もいろんな問題を解いていきます!

Aquaも整理されてないものが多々あるので、今後はもっと使いやすくなると思います。


(追記)QAOA

QAOAを追記します。ほとんど同じですが、アルゴリズムでVQEの代わりにQAOAを使います。参考文献を探していたら、こちらを見つけました。内容はほぼ一緒かと思います。最適化アルゴリズムもPOWELLにしました。

「【量子コンピューター】QAOAでMaxcutを解こう! (最終回)」

https://base.terrasky.co.jp/articles/ONSHr

import numpy as np

import networkx as nx
import matplotlib.pyplot as plt

from qiskit.aqua.translators.ising import max_cut
from qiskit.aqua.input import EnergyInput

from qiskit import BasicAer
from qiskit.aqua.algorithms import QAOA
from qiskit.aqua.components.optimizers import POWELL
from qiskit.aqua import QuantumInstance

n=4

G=nx.Graph()
G.add_nodes_from(np.arange(0,n,1))

#エッジijに対して、(i,j,結合荷重)で設定します。
elist=[(0,1,1.0),(0,2,1.0),(0,3,1.0),(1,2,1.0),(2,3,1.0)]

G.add_weighted_edges_from(elist)
colors = ['r' for node in G.nodes()]
pos = nx.spring_layout(G)
default_axes = plt.axes(frameon=True)
nx.draw_networkx(G, node_color=colors, node_size=600, alpha=.8, ax=default_axes, pos=pos)

w = np.zeros([n,n])
for i in range(n):
for j in range(n):
temp = G.get_edge_data(i,j,default=0)
if temp != 0:
w[i,j] = temp['weight']
print(w)

qubitOp, offset = max_cut.get_max_cut_qubitops(w)
algo_input = EnergyInput(qubitOp)

seed = 10598

optimizer = POWELL()
qaoa = QAOA(qubitOp, optimizer, 1, operator_mode='matrix')

backend = BasicAer.get_backend('statevector_simulator')
quantum_instance = QuantumInstance(backend, seed=seed, seed_transpiler=seed)

result = qaoa.run(quantum_instance)

x = max_cut.sample_most_likely(result['eigvecs'][0])
print('energy:', result['energy'])
print('time:', result['eval_time'])
print('max-cut objective:', result['energy'] + offset)
print('solution:', max_cut.get_graph_solution(x))
print('solution objective:', max_cut.max_cut_value(x, w))

colors = ['r' if max_cut.get_graph_solution(x)[i] == 0 else 'b' for i in range(n)]
nx.draw_networkx(G, node_color=colors, node_size=600, alpha = .8, pos=pos)