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

量子コンピュータVariational Quantum Eigensolver(VQE)、QAOA4周目まとめ

はじめに

2020年の初頭現在で実用的に量子コンピュータで解ける問題は「変分アルゴリズム」と呼ばれるVQEやQAOAと言った量子コンピュータと古典コンピュータをハイブリッドで利用したものが主流になっています。ここでは実用性を重視して、これらの量子古典ハイブリッドアルゴリズムの使い方を見てみたいと思います。

概要

量子古典ハイブリッドアルゴリズムのVQE(Variational Quantum Eigensolver)は、現在の量子コンピュータはエラーがまだ多いため、従来利用されると思われていた「位相推定」アルゴリズムの代替として開発され、2013年に当時ハーバード大学(現在はトロント大学)のアラン・アスプル・グジック教授のチームによって開発されました。

現在の量子コンピュータはエラーが多く、長い量子回路を組むとエラーが蓄積し正しい答えを得ることができません。想定されていたアルゴリズムは全て長い回路が基本となっていました。それに対して、短い量子回路をたくさん計算し、それを集計し、最適化をかけることで同じ機能を実現するVQEアルゴリズムが現在の主流となっています。

スクリーンショット 2020-02-29 14.42.40.png

引用:https://arxiv.org/pdf/1304.3061.pdf

固有値・固有ベクトル問題

VQEで解ける問題は、「固有値・固有ベクトル問題」と呼ばれるもので、ある行列に対する固有値と固有ベクトルを見つけることができれば様々な問題を解くことができます。ある与えられた行列$H$に対して、固有値を$E_0$、固有ベクトルを$\mid \psi \rangle$としたときに、

H \mid \psi \rangle = E_0 \mid \psi \rangle

を満たすような、$E_0$と$\psi$を見つけるのが目的です。$H$はハミルトニアン、$\mid \psi \rangle$は状態ベクトル(波動関数)と呼ばれ、状態ベクトルの固有ベクトルのことを固有状態と呼びます。

実用的に解ける問題

VQEは汎用的なアルゴリズムでハミルトニアン$H$が与えられれば本来どのような固有値・固有ベクトルも探せるはずですが、実際には最適化計算で探索するパラーメータの範囲がとても広いため、与えられたハミルトニアン$H$と求める固有ベクトルである$\mid \psi \rangle$の間にある程度の関連性がないと実質的に解くことができません。

2020年初頭現在で解くことができると言われているのが、「量子化学計算」と「組合せ最適化計算」と呼ばれる2種類で、それぞれ与えられたハミルトニアン$H$と求めたい固有状態$\mid \psi \rangle$の間に関連性を見出すことができ、それらをある程度定式化することができています。

上記の固有状態を作り出すためには量子回路を組んで、その結果として固有状態$\mid \psi \rangle$を得ますが、この量子回路のことを特に、「ansatz(アンザッツ)」と呼びます。量子化学では、UCCansatzと呼ばれるUCC(Unitary Coupled Cluster)理論を利用しansatzを構成し、組合せ最適化では、QAOAansatzと呼ばれるQAOA(Quantum Approximate Optimization Algorithm)理論を利用してansatzを組み立てます。

実質的に行っていることは似通っていますが、特に後者はQAOAという名前で呼ばれることが通常です。ここでは、VQEとQAOAを中心に理論の下準備と実際の活用をblueqatSDK( https://github.com/blueqat/blueqat )を利用して解説していきます。

変分アルゴリズム#2-1(理論編)

手計算を中心に理論を確認します。主に変分アルゴリズムと呼ばれるVQEとQAOAはほぼ同じ手順で実行をすることが可能です。苦手な方はこの章は飛ばしても結構です。主に学ぶ内容は、

・ハミルトニアンについて
・ansatzについて

の二つです。最後に簡単な例題を実行してみます。

blueqatのインストール

今後例題にはblueqatを使います。

https://github.com/Blueqat/Blueqat

手元にpythonやpipのインストールが完了している、
もしくはgoogle colab ( https://colab.research.google.com/ )
などを使う場合にインストールは、

pip3 install blueqat

で完了します。

2-1-1 下準備1:量子ビットと状態ベクトル

量子コンピュータで利用される量子ビットは|0>状態と|1>状態を持っていてベクトルで表現できます。

\mid 0 \rangle =
\begin{bmatrix}
1 \\
0
\end{bmatrix},
\mid 1 \rangle =
\begin{bmatrix}
0 \\
1
\end{bmatrix}

また、任意の量子ビットは重ね合わせて量子状態を表現できます。

\mid \psi \rangle = \alpha \mid 0 \rangle + \beta \mid 1 \rangle =
\alpha 
\begin{bmatrix}
1 \\
0
\end{bmatrix}
+ \beta
\begin{bmatrix}
0 \\
1
\end{bmatrix}
=
\begin{bmatrix}
\alpha \\
\beta
\end{bmatrix}

$\mid \psi \rangle$は状態ベクトル(波動関数)と呼ばれ量子ビットの量子状態を表します。$\alpha$と$\beta$は複素数で「確率振幅」と呼ばれ、二乗するとその対応する計算基底の出現確率になり、$|\alpha|^2$や$|\beta|^2$は「測定」をすることで取得できます。また、

|\alpha|^2+|\beta|^2 = 1

を満たします。例えば、100回計算して、0が45回、1が55回でたら、$|\alpha|^2$の期待値は0.45、$|\beta|^2$の期待値は0.55と計算できます。実際の状態ベクトルは直接測定をすることはできないため、通常は何度も計算を実行して統計的なサンプルを取ります(シミュレータでは直接状態ベクトルを取る機能もありますが、サイズの大きいものは取れません)。

2-1-2 下準備2:状態ベクトルの操作と量子ゲート

状態ベクトルを操作するには、量子ゲートを利用します。変分アルゴリズムでは基本のパウリゲート、任意回転ゲート、アダマールゲートを主に利用します。

パウリゲートはXYZの3種類があり、それぞれブロッホ球でのXYZ軸周りでの180度の回転を意味し、対応する行列は、

X=
\begin{bmatrix}
0&1 \\
1&0
\end{bmatrix},
Y=
\begin{bmatrix}
0&-i \\
i&0
\end{bmatrix},
Z=
\begin{bmatrix}
1&0 \\
0&-1
\end{bmatrix}

となります。また、XYZ軸周りに任意回転にしたのが、RX,RY,RZゲートで、

Rx(\theta) = \left( \begin{array}{cc} \cos\left(\frac{\theta}{2}\right) &
-i\sin\left(\frac{\theta}{2}\right)\\
-i\sin\left(\frac{\theta}{2}\right) &
\cos\left(\frac{\theta}{2}\right) \end{array} \right),
Ry(\theta) = \left( \begin{array}{cc} \cos\left(\frac{\theta}{2}\right) &
-\sin\left(\frac{\theta}{2}\right)\\
\sin\left(\frac{\theta}{2}\right) &
\cos\left(\frac{\theta}{2}\right) \end{array} \right),
Rz(\theta) = \left( \begin{array}{cc} e^{-i\frac{\theta}{2}} & 0\\ 0 & e^{i\frac{\theta}{2}} \end{array} \right)

となります。アダマールゲートは、

H =
\frac{1}{\sqrt{2}} 
\begin{bmatrix}
1&1 \\
1&-1
\end{bmatrix}

となります。主にX軸とZ軸の間で変換を行います。状態ベクトルにこれらの量子ゲートを作用させることで量子状態を変化させることができます。

X
\begin{bmatrix}
1 \\
0
\end{bmatrix}
=
\begin{bmatrix}
0&1 \\
1&0
\end{bmatrix}
\begin{bmatrix}
1 \\
0
\end{bmatrix}
=
\begin{bmatrix}
0 \\
1
\end{bmatrix}

2-1-3 下準備3:複数量子ビットの状態ベクトルとテンソル積

上記の状態ベクトルは1量子ビットの時のものです。2量子ビット以上の状態ベクトルはテンソル積を使って表現します。

\begin{bmatrix}
a \\
b
\end{bmatrix}
\otimes
\begin{bmatrix}
c \\
d
\end{bmatrix}
=
\begin{bmatrix}
a*
\begin{bmatrix}
c\\
d
\end{bmatrix}
 \\
b*
\begin{bmatrix}
c\\
d
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
ac \\
ad \\
bc \\
bd
\end{bmatrix}

3量子ビット以上の場合も同様です。

\begin{bmatrix}
a \\
b
\end{bmatrix}
\otimes
\begin{bmatrix}
c \\
d
\end{bmatrix}
\otimes
\begin{bmatrix}
e \\
f
\end{bmatrix}
=

\begin{bmatrix}
ac *
\begin{bmatrix}
e \\
f
\end{bmatrix} \\
ad*
\begin{bmatrix}
e \\
f
\end{bmatrix} \\
bc*
\begin{bmatrix}
e \\
f
\end{bmatrix} \\
bd*
\begin{bmatrix}
e \\
f
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
ace\\
acf\\
ade\\
adf\\
bce\\
bcf\\
bde\\
bdf
\end{bmatrix}

量子ゲートもテンソル積をとれます。ハミルトニアンの計算の際に利用します。

\begin{bmatrix}
a&b \\
c&d
\end{bmatrix}
\otimes
\begin{bmatrix}
e&f \\
g&h
\end{bmatrix} \\
=
\begin{bmatrix}
a*
\begin{bmatrix}
e&f \\
g&h
\end{bmatrix}
&b*
\begin{bmatrix}
e&f \\
g&h
\end{bmatrix}
\\
c*
\begin{bmatrix}
e&f \\
g&h
\end{bmatrix}
&d*
\begin{bmatrix}
e&f \\
g&h
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
ae&af&be&bf\\
ag&ah&bg&bh\\
ce&cf&de&df\\
cg&ch&dg&dh
\end{bmatrix}

2-1-4 量子変分原理

今回固有値を直接求めるのは厳しいので、間接的に固有値の期待値にできるだけ近づける最適化計算に落とし込みます。量子変分原理より、任意の状態ベクトルにおけるハミルトニアンの期待値は基底状態を上回るというのがあります。

\langle \Psi (\theta) \mid H \mid \Psi (\theta) \rangle \geq E_0

これは、与えられたハミルトニアン$H$に対して任意の状態ベクトルは常に最小の固有値の期待値を上回るということを示しています。ハミルトニアン$H$は与えられた式なので、私たちにできるのは状態ベクトル$\mid \Psi (\theta) \rangle$をできるだけ固有状態に近づけるということです。これにより、できるだけ固有状態に近い状態ベクトルを角度の最適化計算として探していくのが今回の目的です。

2-1-5 変分アルゴリズムの計算手順

変分アルゴリズムの計算手順は、

1、パラメータ$\theta$を導入した量子回路(ansatz)を準備(量子計算)
2、量子ビットの測定値からハミルトニアン$H$の期待値$\langle \Psi (\theta) \mid H \mid \Psi (\theta) \rangle$を計算(古典計算)
3、古典最適化計算により、ハミルトニアン$H$の期待値を最小にするパラメータ$\theta$を探索(古典計算)

のように、量子計算と古典計算のハイブリッドで成り立っています。

スクリーンショット 2020-02-29 14.42.40.png

引用:https://arxiv.org/pdf/1304.3061.pdf

注意点は、VQEは量子回路はハミルトニアン$H$がユニタリ行列の場合の期待値しか計算できません。ハミルトニアン$H$は通常エルミート行列で与えられるので、ユニタリ行列の和に分解する必要があります。

2-1-6 ansatz(アンザッツ)

ansatz(アンザッツ)は、状態ベクトルを作るための、主に任意回転ゲートでパラメータ化された量子回路で、ハミルトニアン$H$の期待値を求めるための測定値を提供します。

2-1-7 ハミルトニアン

ハミルトニアン$H$はパウリゲートXYZと単位行列Iから構成される行列です。今回の変分アルゴリズムは、このハミルトニアン$H$の固有状態である状態ベクトルから導き出される固有値の期待値を求めることを目的としています。

2-1-8 ハミルトニアンの期待値

ある状態ベクトル$\mid \psi \rangle$とハミルトニアン$H$が与えられた時のハミルトニアン$H$の期待値は、

\langle \psi (\theta) \mid H \mid \psi (\theta) \rangle

の簡単な計算によってもとまります。

2-1-9 例題:1量子ビットのansatzを使ってハミルトニアンZの期待値を求める

ansatzを使った任意の状態ベクトルにおけるハミルトニアンが$H = Z$の場合の期待値は、

\langle \psi \mid H \mid \psi \rangle = 

\begin{bmatrix}
\alpha^* & \beta^*
\end{bmatrix}

\begin{bmatrix}
1&0\\
0&-1
\end{bmatrix}

\begin{bmatrix}
\alpha\\
\beta
\end{bmatrix}

= |\alpha|^2 - |\beta|^2

となります。ここで、$|\alpha|^2$は0がでる確率、$|\beta|^2$は1が出る確率に対応します。

2-1-10 例題2:1量子ビットのansatzを使ってハミルトニアンXの期待値を求める

ハミルトニアンが$H = X$の場合のハミルトニアン$H$の期待値は、

\langle \psi \mid H \mid \psi \rangle = 

\begin{bmatrix}
\alpha^* & \beta^*
\end{bmatrix}

\begin{bmatrix}
0&1\\
1&0
\end{bmatrix}

\begin{bmatrix}
\alpha\\
\beta
\end{bmatrix}

= \alpha^* \beta + \alpha \beta^*

となりますが、通常の測定ではこの値を求めることができませんので、ハミルトニアンの変形を通じて解を得ます。

2-1-11 ハミルトニアンの変形

ハミルトニアン$H$が$Z$の場合には、期待値の計算ができますが、XやYの場合には求めることができません。そこでハミルトニアンの変形を利用します。

X = HZH

より、

\langle \psi \mid X \mid \psi \rangle \\
= \langle \psi \mid HZH \mid \psi \rangle\\
= \langle \psi' \mid Z \mid \psi' \rangle

のように状態ベクトルを変形することで対応できます。Yの場合は、

Y = RX(-\pi/2) Z RX(\pi/2)

と、なり同様に、

\langle \psi \mid Y \mid \psi \rangle \\
= \langle \psi \mid RX(-\pi/2) Z RX(\pi/2) \mid \psi \rangle\\
= \langle \psi'' \mid Z \mid \psi'' \rangle

とすることができます。これらは、量子回路の測定の直前に上記のHやRXゲートを入れて測定を行うことで期待値を計算することができます。

2-1-12 線形結合

実際にはハミルトニアン$H$はエルミート行列で与えられ、ユニタリ行列の和の形で与えられ、下記のように分解できます。

\langle \psi \mid aH_1 + bH_2 \mid \psi \rangle \\ = \langle \psi \mid aH_1 \mid \psi \rangle + \langle \psi \mid bH_2 \mid \psi \rangle \\ = a\langle \psi \mid H_1 \mid \psi \rangle + b\langle \psi \mid H_2 \mid \psi \rangle

これを使えば、複雑なハミルトニアンの場合でも、

H = 1.2 X_0 Z_2 + 2.5 Z_0 X_1 Y_2 - 3.4 Z_2 X_1

期待値は、

\langle \psi \mid 1.2 X_0 Z_2 + 2.5 Z_0 X_1 Y_2 - 3.4 Z_2 X_1 \mid \psi \rangle\\
= 1.2\langle \psi \mid X_0 Z_2 \mid \psi \rangle + 2.5 \langle \psi \mid Z_0 X_1 Y_2\mid \psi \rangle - 3.4 \langle \psi \mid Z_2 X_1 \mid \psi \rangle

のように個別に計算をして和を求めることで計算できます。上記の具体的な行列の確認の仕方は、次の実装初級編で触れます。

2-1-13 例:任意のansatzで実装

ハミルトニアン$H$に対して対応したansatzを使わないと太平洋から一円玉を見つけるような作業になってしまいますが、ここでは例題として1量子ビットのハミルトニアン$H$の固有値を求めてみます。求めたいハミルトニアン$H$は、

from blueqat.pauli import X, Y, Z, I

h = 1.23 * I - 4.56 * X(0) + 2.45 * Y(0) + 2.34 * Z(0)

のように作ってみました。読み込まれたのはパウリゲートで行列が対応します。パウリゲートの後ろの数字はそれぞれを適用する量子ビットの通し番号です。上記のハミルトニアン$H$は参考として行列表記は、

[[ 3.57,   -4.56-2.45j]
 [-4.56+2.45j, -1.11  ]]

となります。そして、ansatzはこれから探しますが探し方はパラメータを割り振ってそれを探します。1量子ビットの任意の状態ベクトルは極座標で2つの角度があれば表現できます。

Circuit().ry(a)[0].rz(b)[0]

としてみます。上記は0番目の量子ビットに、Ryゲートを角度aで、Rzゲートを角度bで適用します。aとbでいい値が見つかれば、上記のansatzはハミルトニアンの期待値に対していい答えを返します。

全体の実装は、

import numpy as np
from blueqat import Circuit
from blueqat.pauli import X, Y, Z, I
from blueqat.vqe import AnsatzBase, Vqe

class OneQubitAnsatz(AnsatzBase):
    def __init__(self, hamiltonian):
        super().__init__(hamiltonian.to_expr(), 2)
        self.step = 1

    def get_circuit(self, params):
        a, b = params
        return Circuit().ry(a)[0].rz(b)[0] #パラメータでansatzを作る


# ハミルトニアン
h = 1.23 * I - 4.56 * X(0) + 2.45 * Y(0) + 2.34 * Z(0)
runner = Vqe(OneQubitAnsatz(h))
result = runner.run()

print('Result by VQE')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))

# これ以降は確認のために
mat = h.to_matrix()

# Calculate by numpy
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])

これを計算すると、

Result by VQE
-4.450591481563331
Result by numpy
-4.450818602983201

となり、ほぼ近い値が固有値として求まっているのがわかります。実際には、上記のように任意でハミルトニアン$H$を書くことは多くないので定式化にルールを定めます。今後は全てblueqatで自動でやってくれますので、理解だけすれば実装は下記のように非常に簡単に実装できます。

from blueqat import vqe
from blueqat.pauli import Z

H = Z(0)*Z(1) #ハミルトニアン
step = 2

result = vqe.Vqe(vqe.QaoaAnsatz(H, step)).run()
print(result.most_common(12))

これ以降はansatzとハミルトニアンの関連に触れながら実問題を見ていきます。実問題でこれからみるのは、

1、量子化学計算
2、組合せ最適化問題

になります。

2-3 VQEで量子化学計算

量子化学計算はハミルトニアンの作成をパウリ行列を使って行い、VQEを利用して解を得ることができます。量子化学計算はハミルトニアンに対応するansatzの研究が進んでおり、課題として見通しが他の問題よりもよくなっています。

2-3-1 追加のインストール

追加でツールのインストールが必要です。
今回は、blueqatに追加で、openfermionblueqatとopenfermionを追加します。

pip install blueqat openfermionblueqat openfermion

準備OKです。

2-3-2 コードをざっくり

今回はopenfermionに準備されたコードを使います。

from blueqat import *
from openfermion import *
from openfermionblueqat import*
import numpy as np

def get_molecule(bond_len):
  geometry = [('H',(0.,0.,0.)),('H',(0.,0.,bond_len))]

  description = format(bond_len)
  molecule = MolecularData(geometry, "sto-3g",1,description=description)

  molecule.load()
  return molecule

x = [];e=[];fullci=[]
for bond_len in np.arange(0.2,2.5,0.1):
  m = get_molecule("{:.2}".format(bond_len))
  h = bravyi_kitaev(get_fermion_operator(m.get_molecular_hamiltonian()))
  runner = vqe.Vqe(UCCAnsatz(h,6,Circuit().x[0]))
  result = runner.run()
  x.append(bond_len)
  e.append(runner.ansatz.get_energy(result.circuit,runner.sampler))
  fullci.append(m.fci_energy)

%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(x,fullci)
plt.plot(x,e,"o")

こちらが正しく実行されると下記のような画像が出力されます。

ダウンロード.png

今回は水素分子の原子の結合間の距離に応じて、各状態の安定エネルギーをVQEを用いて求めています。

2-3-3 中身確認

MoluculearDataというのでOpenFermionに最初から入っているデータを使えるようです。H2でSto-3g系を読み込んでみます。

def get_molecule(bond_len):
  geometry = [('H',(0.,0.,0.)),('H',(0.,0.,bond_len))]

  description = format(bond_len)
  molecule = MolecularData(geometry, "sto-3g",1,description=description)

  molecule.load()
  return molecule

原子間距離が0.4のものを読み込んで、

m = get_molecule(0.4
m.get_molecular_hamiltonian()

ハミルトニアンと呼ばれる基本の式が取得できます。

() 1.322943021475
((0, 1), (0, 0)) -1.4820918858979102
((1, 1), (1, 0)) -1.4820918858979102
((2, 1), (2, 0)) -0.1187350527865787
((3, 1), (3, 0)) -0.1187350527865787
((0, 1), (0, 1), (0, 0), (0, 0)) 0.36843967630348756
((0, 1), (0, 1), (2, 0), (2, 0)) 0.08225771204699692
((0, 1), (1, 1), (1, 0), (0, 0)) 0.36843967630348756
((0, 1), (1, 1), (3, 0), (2, 0)) 0.08225771204699692
((0, 1), (2, 1), (0, 0), (2, 0)) 0.082257712046997
((0, 1), (2, 1), (2, 0), (0, 0)) 0.3626667179796745
((0, 1), (3, 1), (1, 0), (2, 0)) 0.082257712046997
((0, 1), (3, 1), (3, 0), (0, 0)) 0.3626667179796745
((1, 1), (0, 1), (0, 0), (1, 0)) 0.36843967630348756
((1, 1), (0, 1), (2, 0), (3, 0)) 0.08225771204699692
((1, 1), (1, 1), (1, 0), (1, 0)) 0.36843967630348756
((1, 1), (1, 1), (3, 0), (3, 0)) 0.08225771204699692
((1, 1), (2, 1), (0, 0), (3, 0)) 0.082257712046997
((1, 1), (2, 1), (2, 0), (1, 0)) 0.3626667179796745
((1, 1), (3, 1), (1, 0), (3, 0)) 0.082257712046997
((1, 1), (3, 1), (3, 0), (1, 0)) 0.3626667179796745
((2, 1), (0, 1), (0, 0), (2, 0)) 0.36266671797967454
((2, 1), (0, 1), (2, 0), (0, 0)) 0.08225771204699726
((2, 1), (1, 1), (1, 0), (2, 0)) 0.36266671797967454
((2, 1), (1, 1), (3, 0), (0, 0)) 0.08225771204699726
((2, 1), (2, 1), (0, 0), (0, 0)) 0.08225771204699728
((2, 1), (2, 1), (2, 0), (2, 0)) 0.38272169831413727
((2, 1), (3, 1), (1, 0), (0, 0)) 0.08225771204699728
((2, 1), (3, 1), (3, 0), (2, 0)) 0.38272169831413727
((3, 1), (0, 1), (0, 0), (3, 0)) 0.36266671797967454
((3, 1), (0, 1), (2, 0), (1, 0)) 0.08225771204699726
((3, 1), (1, 1), (1, 0), (3, 0)) 0.36266671797967454
((3, 1), (1, 1), (3, 0), (1, 0)) 0.08225771204699726
((3, 1), (2, 1), (0, 0), (1, 0)) 0.08225771204699728
((3, 1), (2, 1), (2, 0), (3, 0)) 0.38272169831413727
((3, 1), (3, 1), (1, 0), (1, 0)) 0.08225771204699728
((3, 1), (3, 1), (3, 0), (3, 0)) 0.38272169831413727

こちらはパウリ行列の表現ではないので、今回VQEで計算できるように変換をします。

2-3-4 量子コンピュータ向け変換

Bravi-Kitaev(ブラビキタエフ)変換やJordan–Wigner(ジョルダンウィグナー)変換がありますが、今回は前者のBK変換をしてみると、

h = bravyi_kitaev(get_fermion_operator(m.get_molecular_hamiltonian()))

確認してみると、

print(h)

(0.7407724940116754+0j)*I + (0.23528824284103544+0j)*Z[0] + (0.23528824284103542+0j)*Z[0]*Z[1] + (-0.45353118471995524+0j)*Z[2] + (-0.45353118471995524+0j)*Z[1]*Z[2]*Z[3] + (0.18421983815174378+0j)*Z[1] + (0.041128856023498556+0j)*Y[0]*Z[1]*Y[2]*Z[3] + (0.041128856023498556+0j)*X[0]*Z[1]*X[2] + (0.041128856023498556+0j)*X[0]*Z[1]*X[2]*Z[3] + (0.041128856023498556+0j)*Y[0]*Z[1]*Y[2] + (0.14020450296633868+0j)*Z[0]*Z[2] + (0.18133335898983727+0j)*Z[0]*Z[1]*Z[2]*Z[3] + (0.18133335898983727+0j)*Z[0]*Z[1]*Z[2] + (0.14020450296633868+0j)*Z[0]*Z[2]*Z[3] + (0.19136084915706864+0j)*Z[1]*Z[3]

このようにきちんと変換されました。ちょっと複雑な形をしていますが、これこそがこれまで学んできたパウリ行列の和で書かれたハミルトニアンです。これをVQEにかけますが、今回はUCC理論と呼ばれる上記のハミルトニアンに対応するansatzがあるため効率的に計算ができます。

runner = vqe.Vqe(UCCAnsatz(h,2,Circuit().x[0]))
result = runner.run(verbose = True)

そうすると、最適化計算が進み、

params: [0.00803084 0.88310007] val: -0.3685674394653079
params: [0.00803084 0.88310007] val: -0.3685674394653079
(中略)
params: [-0.0024937  -0.00256402] val: -0.9043519838381495
params: [-0.0024937  -0.16023561] val: -0.885812805585636

最適化のプロセスが出てきました。パラメータが最適化されながら最小基底状態が探索されています。今回はansatzは量子化学の理論に基づいて電子の配置をするUCCansatzと呼ばれるものが利用されました。通常ansatzの構成は難しいので、理論的な探究が求まれます。最小エネルギーは、

runner.ansatz.get_energy(result.circuit,runner.sampler)

-0.9043519838381495

原子間距離0.4の場合の最小値の期待値が求まりました。これを全体でやったのが最初のコードです。再掲します。

from blueqat import *
from openfermion import *
from openfermionblueqat import*
import numpy as np

def get_molecule(bond_len):
  geometry = [('H',(0.,0.,0.)),('H',(0.,0.,bond_len))]

  description = format(bond_len)
  molecule = MolecularData(geometry, "sto-3g",1,description=description)

  molecule.load()
  return molecule

x = [];e=[];fullci=[]
for bond_len in np.arange(0.2,2.5,0.1):
  m = get_molecule("{:.2}".format(bond_len))
  h = bravyi_kitaev(get_fermion_operator(m.get_molecular_hamiltonian()))
  runner = vqe.Vqe(UCCAnsatz(h,6,Circuit().x[0]))
  result = runner.run()
  x.append(bond_len)
  e.append(runner.ansatz.get_energy(result.circuit,runner.sampler))
  fullci.append(m.fci_energy)

%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(x,fullci)
plt.plot(x,e,"o")

実行結果は、

ダウンロード.png

2-4 組合せ最適化問題

変分アルゴリズムの応用として組合せ最適化問題があります。多くの人に馴染みのある問題なので、量子コンピュータは正直どっから手をつけていいかわからんと思うひとおすすめです。

組合せ最適化問題は、多くの選択肢からベストな答えを選ぶ問題で、社会問題をうまく定式化し、最小値問題を解くことで最適な組合せを得ることができます。主に選択肢を選んだ状態を1、選ばない状態を0としますが、様々な条件をつけながら問題を解きます。

2-4-1 例:VQEで組合せ最適化問題の定式化

定式化はハミルトニアンをコスト関数として扱い、最小値を返す量子ビットの組合せを解とします。ルールは、

・イジングモデルに落とし込む
・Zを使う
・最終的にはZの代わりにQUBOを使う

です。今回は通常VQEを組合せ最適化問題としてはあまり利用しませんが、例題のハミルトニアンを実行してみます。

h = -Z(0) - Z(0)*Z(1)

Zの後ろの数字は、量子ビットの通し番号を表します。0番目と1番目の量子ビットの2つを使っています。また、問題設定で大事なのは、Zの前の係数です。

Z(0)の前は-1のバイアス
Z(0)*Z(1)の前は-1のウェイト

が設定されています。Zは期待値として-1か+1のどちらかをとります。hはより小さい値を取ると正解になります。最終的な答えは、Z(0),Z(1)の値で場合わけすると、

Z(0) Z(1) h
-1 -1 0
-1 1 2
1 -1 0
1 1 -2

VQEは上記の表で最小となるZ(0)=1,Z(1)=1の時-2となるものを計算で求めてくれます。今回ansatzはa,b,c,dの4パラメータを利用した極座標表記のものを使ってみます。ansatzを含む全体のコードは、

import numpy as np
from blueqat import Circuit
from blueqat.pauli import X, Y, Z, I
from blueqat.vqe import AnsatzBase, Vqe

class OneQubitAnsatz(AnsatzBase):
    def __init__(self, hamiltonian):
        super().__init__(hamiltonian.to_expr(), 4)
        self.step = 1

    def get_circuit(self, params):
        a, b, c, d = params
        return Circuit().ry(a)[0].rz(b)[0].ry(c)[1].rz(d)[1]


# この定式化が大事
h = -Z(0) - Z(0)*Z(1)
runner = Vqe(OneQubitAnsatz(h))
result = runner.run()

print('Result by VQE')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))

としました。計算を実行すると、

Result by VQE
-1.9999943170511056

約-2が答えとして出てきましたので、正解です。今回は例題として任意の量子状態を使いましたが、通常はQAOAを使います。そして実際にはより大きな問題を解きます。上記問題は役にたたなさそうですが、Z(0)をAさん、Z(1)をBさんに見立てて、Aさんはグループ1に属し、BさんはAさんと同じグループに属するという条件をつけた分類問題と同じです。

ただ、このままでは毎回問題を解くのが大変なので様々な工夫が必要です、それをみていきましょう。

2-4-2 定式化は0と1のバイナリ値で

定式化は組合せ最適化問題を+1と-1の値の組合せで表される式に落とし込みます。ただ、課題があります。物理学で使われるイジングモデルは-1,+1を利用しますが、通常の産業用では0と1を計算として利用します(計算基底)。幸い-1と+1は自動的に変換ができますので、01を使っても組合せ最適化問題としては問題がありません。

01で定式化をするのをQUBOと呼びます。QUBOは量子ビットの01を利用できます。係数は様々な呼び方がありますが「バイアス」と「ウェイト(結合荷重)」と呼びます。定式化で作るのは「コスト関数」と呼びたいと思います。

-1と+1でかかれたイジング式を0と1で書かれたQUBO式に変換するには、イジングのZを下記のように変換するだけでできます。

q = \frac{Z + 1}{2}

これで、-1の時が0に、1の時は1のままで変換されます。定式化は社会問題をQUBO形式で、コスト関数を作ることで実現でき、コスト関数は01の値をとる量子ビットに設定するバイアスとウェイトを設定することで実現できます。この手法は量子コンピュータに限ったことではないので普通の計算に慣れている人でも受け入れられやすいでしょう。

2-4-3 QUBO式のプログラミング

別の問題をQUBOで解いてみるため、blueqatの機能を利用します。QUBOへの変換はZを代入して書き換えればいいだけでした。

今回のときたいコスト関数は、

h = -3*q(0)-3*q(1)-2*q(0)*q(1)

明らかにq(0)=1,q(1)=1の時に最小値-3-3-2=-8を取ります。上記のコスト関数は以前のように-1と+1ではなく、0と1で考えることができるので簡単です。これをVQEで解くことで、最小値の期待値-8が得られます。

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comを利用する時に必要
import numpy as np
from blueqat import Circuit
from blueqat.pauli import X, Y, Z, I
from blueqat.pauli import qubo_bit as q
from blueqat.vqe import AnsatzBase, Vqe

class QubitAnsatz(AnsatzBase):
    def __init__(self, hamiltonian):
        super().__init__(hamiltonian, 4)
        self.step = 1

    def get_circuit(self, params):
        a, b, c, d = params
        return Circuit().ry(a)[0].rz(b)[0].ry(c)[1].rz(d)[1]

h = -3*q(0)-3*q(1)-2*q(0)*q(1)
h = h.to_expr().simplify()
runner = Vqe(QubitAnsatz(h))
result = runner.run()

print('Result by VQE')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))

# Hamiltonian to matrix
mat = h.to_matrix()

# Calculate by numpy
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])

結果は、

Result by VQE
-7.981654210576551
Result by numpy
-8.0

正しく得られています。実際にはansatzを作るのが難しかったり、その他定式化に有利な回路を作るためにQAOAを利用するのが無難です。

2-5 QAOA

VQEはansatzを高度化することで問題特化型のアルゴリズムとして発展をさせることができます。ここでは、VQEの変分法をベースに、ansatzを組合せ最適化問題に特化した形で発展をしたQAOA(Quantum Approximate Opitmization Alogirthm)をみてみます。

2-5-1 量子断熱計算

量子断熱計算は、状態ベクトルを断続的に変化させることで基底状態をキープしたまま時間発展を行うことができる理論です。

初期状態のハミルトニアンを$H_{start}$として、最終的に求めたい問題のハミルトニアンを$H_{fin}$としたときに、時間$t$と全体のスケジュール$T$から、

H_{temp} = (1-\frac{t}{T})H_{start} + \frac{t}{T}H_{fin}

としたときに$T\rightarrow\infty$とすれば、時間発展で変化させた状態ベクトルが、その瞬間瞬間のハミルトニアンに追従し、固有状態をとり固有値$\lambda$を持つようにすることができます。

H_{temp}\mid \psi \rangle = E_{0temp}\mid \psi \rangle

時間発展計算は、

\mid \psi_{t+1} \rangle = U \mid \psi_t \rangle = e^{-iHt}  \mid \psi_t \rangle

となります。課題は基底状態と第一励起状態が最接近する部分ですが、$E_1(t)-E_0(t)$のエネルギー差に注意して計算することによって、基底状態をキープできます。

2-5-2 QAOA

QAOAは上記の量子断熱計算の時間発展計算をansatzとして変分アルゴリズムに適用したものです。

https___qiita-user-contents.imgix.net_https%3A%2F%2Fqiita-image-store.s3.ap-northeast-1.amazonaws.com%2F0%2F218694%2Fe10f1843-cc16-cdfe-e4a6-e2fbaab6df9f.png_ixlib=rb-1.2.png

一番左のHは初期の固有状態を作っています。これに対応するハミルトニアンはXです。また、CX-Rz-CXは問題設定のハミルトニアンのウェイトに対応し、次のRzはハミルトニアンのバイアスに対応しています。Rxは上記のハミルトニアンXを時間発展させたものです。

試しにblueqatに用意されたQaoaAnsatzをみてみます。ハミルトニアンは今回Z演算子から構成されており、自動的にハミルトニアンから時間発展のansatzを構成しているのがみて取れます。

class QaoaAnsatz(AnsatzBase):
    def __init__(self, hamiltonian, step=1, init_circuit=None):
        super().__init__(hamiltonian, step * 2)
        self.hamiltonian = hamiltonian.to_expr().simplify()
        if not self.check_hamiltonian():
            raise ValueError("Hamiltonian terms are not commutable")

        self.step = step
        self.n_qubits = self.hamiltonian.max_n() + 1
        if init_circuit:
            self.init_circuit = init_circuit
            if init_circuit.n_qubits > self.n_qubits:
                self.n_qubits = init_circuit.n_qubits
        else:
            self.init_circuit = Circuit(self.n_qubits).h[:]
        self.init_circuit.make_cache()
        self.time_evolutions = [term.get_time_evolution() for term in self.hamiltonian]

    def check_hamiltonian(self):
        """Check hamiltonian is commutable. This condition is required for QaoaAnsatz"""
        return self.hamiltonian.is_all_terms_commutable()

    def get_circuit(self, params):
        c = self.init_circuit.copy()
        betas = params[:self.step]
        gammas = params[self.step:]
        for beta, gamma in zip(betas, gammas):
            beta *= np.pi
            gamma *= 2 * np.pi
            for evo in self.time_evolutions:
                evo(c, gamma)
            c.rx(beta)[:]
        return c

ライブラリ側で適切な計算をしてくれているので、これ以降は、ほぼSDK側で処理されたQAOAを使った計算をしてみます。今回はQUBO形式でQAOAライブラリを実行します。今回は簡単な定式化を解いてみましょう。

cost = -3*q(0)-3*q(1)-2*q(0)*q(1)

明らかにq(0)=1,q(1)=1の時に最小値-3-3-2=-8を取ります。上記のコスト関数は以前のように-1と+1ではなく、0と1で考えることができるので簡単です。

設定するのは、q(0)とq(1)に-3のバイアスを設定し、q(0)*q(1)に-2のウェイトを設定します。

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comの時だけ必要
from blueqat import vqe
from blueqat.pauli import qubo_bit as q

h = -3*q(0)-3*q(1)-2*q(0)*q(1)
step = 2

result = vqe.Vqe(vqe.QaoaAnsatz(h, step)).run()
print(result.most_common(12))

ハミルトニアンの設定と時間発展の分割を指定しました。回路はちょっと複雑になりました。

これをQAOAで解くことで、量子ビットが共に1の時、最小値-8が得られます。

(((1, 1), 0.9922852962030535), ((0, 0), 0.0047331664261775755), ((1, 0), 0.0014907686853841205), ((0, 1), 0.0014907686853841168))

最初の数字が答えの候補で、次の数字はそれが現れる確率です。(1,1)の組合せが最大確率になっているので正解です。01の量子ビットでの定式化をして、あとはツールに任せればいいというのが分かりました。最初はツールに任せて定式化を頑張りましょう!

2-5-3 実例:最小値問題における素因数分解VQF

素因数分解を最小値問題に置き換えて解いてみます。

Quantum factorization of 56153 with only 4 qubits

Nikesh S. Dattani,1,2,∗ Nathaniel Bryans3,†
1Quantum Chemistry Laboratory, Kyoto University, 606-8502, Kyoto, Japan, 2Physical & Theoretical Chemistry Laboratory,
Oxford University, OX1 3QZ, Oxford, UK, 3University of Calgary, T2N 4N1, Calgary, Canada. ∗
dattani.nike@gmail.com,

nbryans1@gmail.com

https://arxiv.org/pdf/1411.6758.pdf

Variational Quantum Factoring

Eric R. Anschuetz, Jonathan P. Olson, Alán Aspuru-Guzik, Yudong Cao
(Submitted on 27 Aug 2018)
https://arxiv.org/abs/1808.08927

Prime Factorization Using Quantum Annealing and Algebraic Geometry
https://1qbit.com/wp-content/uploads/2016/09/1QBit-Research-Paper-%E2%80%93-Prime-Factorization-Using-Quantum-Annealing-and-Algebraic-Geometry.pdf

具体的手順

手順はとてもシンプルです。これまでゲート方式はShorのアルゴリズムを使って問題を解いていましたが、そうではなくて最小値問題に落とし込みます。効率的な落とし込みかたがあり、前処理をいれて143の素因数分解は4量子ビットだけで解けます。

1、量子ビットを使って二進数で乗算を考える
2、桁数の特性を生かしてQUBO問題を簡略化する
3、単純化されたQUBO問題をQAOAアルゴリズムにかけて最小基底状態を求める。

1、量子ビットを使って二進数の乗算(143の素因数分解)

まず二進数の乗算を考えます。二進数の乗算と加算を学ぶ必要がありますが、乗算は筆算では10進数とほとんど同じです。

11.png

これは、m=p*qを考えますが、pを二進数で考えて、一番下と一番上の位は1になることはわかっています。一番上が1は自明で、素数なので一番下の位は0ではない(偶数ではない)です。pとqを二進数で表すと上記のようになります。

まずはこれを単純にかけます。$2^0$の位から10進数のようにかけて、ずらして足し合わせてきます。途中carriesというのは桁上がりになります。二進数で1+1をすると桁が上がるので、$Z_{12}$は$2^1$の位から$2^2$の位に上がると言う意味です。

2、桁数の特性を生かしてイジング問題を簡略化する次に各位を足し合わせます。足し合わせると連立方程式は、

\begin{eqnarray}p_1+q_1 &=& 1+2z_{12}\\ p_2+p_1q_1+q_2+z_{12} &=& 1+2z_{23}+4z_{24}\\ 1+p_2q_1+p_1q_2+1+z_{23} &=& 1+2z_{34}+4z_{35}\\ … q_2+p_2+z_{45}+z_{35} &=& 0+2z_{56}+4z_{57}\\ 1+z_{56}+z_{46} &=& 0+2z_{67}\\ z_{67}+z_{57} &=& 1\end{eqnarray}

となります。ここで上から見ていくと、、、まず$p_1+q_1=1+2z_{12}$において、1が確定しているため、$p_1$と$q_1$はどちらかが1であるため、2の位はでてこないため、$z_{12}=0$がわかります。

同様に、$p_2+p_1q_1+q_2+z_{12} = 1+2z_{23}+4z_{24}$はまず、上記より$p_1$と$q_1$はどちらかが1でどちらかが0であるため、$p_1q_1=0$がわかる。また、$z_{12}=0$だったため、$p_2+q_2 = 1+2z_{23}+4z_{24}$となる。ここで、前述と同様1が確定しているため、$z_{23}=0,z_{24}=0$がわかります。

また、$1+p_2q_1+p_1q_2+1+z_{23} = 1+2z_{34}+4z_{35}$は、$z_{23}=0$より、$1+p_2q_1+p_1q_2 = 2z_{34}+4z_{35}$となる。まず桁数が合わないので$z_{35}=0$となり、$1+p_2q_1+p_1q_2 = 2z_{34}$ここで、左辺に1があるため、$z_{34}$は0でないので1が確定する。

これらを進めていってまとめると、下記の式が得られます。

\begin{eqnarray}p_1+q_1-1 &=&0\\
p_2+q_2-1 &=& 0\\
p_2q_1+p_1q_2-1 &=&0
\end{eqnarray}

この3つの連立方程式を最小値問題に落とすためには、各方程式を二乗すると必ず0の時に最小になるので、連立方程式から下記のハミルトニアンを得ます。

H=(p_1+q_1-1)^2 + (p_2+q_2-1)^2+(p_2q_1+p_1q_2-1)^2\\=5-3p_1-p_2-q_1+2p_1q_1-3p_2q_1+2p_1p_2q_1-3q_2+p_1q_2+2p_2q_2+2p_2q_1q_2

これを確認のために、総当たりでコスト計算してみると、バイナリ値で$\mid 0000 \rangle$から$\mid 1111\rangle$までで、

$5,2,4,1,4,3,0,1,2,0,3,1,1,1,1,3$

となり、$\mid 6 \rangle$と$\mid 9 \rangle$の時に$H=0$となるので、先に対応する、$\mid 0110 \rangle$と$\mid 1001 \rangle$が答えとなり、

$p=13,q=11$と$p=11,q=13$が答えとなります。

3、問題をQAOAアルゴリズムにかけて最小基底状態を求める。

実際にQAOAを使って正しいか確認します。

from blueqat import vqe 
from blueqat.pauli import qubo_bit as q 

hamiltonian = -3*q(0)-q(1)-q(2)+2*q(0)*q(2)-3*q(1)*q(2)+2*q(0)*q(1)*q(2)-3*q(3)+q(0)*q(3)+2*q(1)*q(3)+2*q(1)*q(2)*q(3) 
step = 4 

result = vqe.Vqe(vqe.QaoaAnsatz(hamiltonian, step)).run() 
print(result.most_common(5)) 

計算結果は、

(((0, 1, 1, 0), 0.38282701807018904), ((1, 0, 0, 1), 0.17092057836615537), ((1, 1, 1, 0), 0.08073422675654193), ((0, 1, 1, 1), 0.08073422675654193), ((1, 0, 1, 1), 0.045786141364148956))

ここで、きちんと答えがえら得ていて、0110と1001が十分な確率振幅で得られている。このため量子ゲート方式で基底状態が求まり、素因数分解ができました。

ちなみに量子アニーリング方式でも解くことができるが、この場合には上記にある3量子ビットの多体問題の2体問題への分解が必要となる。量子ゲートでは3量子ビットの相互作用をハミルトニアンに入れたまま計算できるので、量子ゲートの方が前処理がとても簡単に済むと言う利点があります。

2-5-4 実例:交通最適化

交通最適化は度々出てくるテーマです。やることは、

1、1台の車に3つのルート選択を許容する(古典)
2、現在の選択ルートから混雑状況を計算(古典)
3、混雑度を最小にするようにルート選択を最小化(QAOA)

となります。参考を確認します。

「Quantum Computing at Volkswagen:
Traffic Flow Optimization using the D-Wave Quantum Annealer」
引用:https://www.dwavesys.com/sites/default/files/VW.pdf

スタートAからゴールBまで12の道路が通し番号0から11まで付けられています。自動車1と自動車2があり、それぞれ指定されたルートを取ります。ルートはそれぞれ3つずつ提案されており、そのうちの一つを取ります。今回は一番混雑度が低くなるようなルートを選択します。

一つ目は混雑度をコスト関数として実装します。上記からs0,s3,...のように、道路を指定して各出現回数を足し合わせて二乗します。上記提案ルートには、q0からq5まで通し番号を降っています。

(q_0+q_1+q_3+q_4)^2+(q_0+q_1+q_3+q_4)^2+(q_0+q_3)^2+(q_0+q_3)^2+(q_1+q_4)^2+(q_1+q_2+q_4+q_5)^2+(q_2+q_5)^2+(q_2+q_5)^2+(q_2+q_5)^2\\ 
=2(q_0+q_1+q_3+q_4)^2+2(q_0+q_3)^2+(q_1+q_4)^2+(q_1+q_2+q_4+q_5)^2+3(q_2+q_5)^2\\ 
=4q_0^2+4q_0q_1+8q_0q_3+4q_0q_4+4q_1^2+2q_1q_2+4q_1q_3+8q_1q_4+2q_1q_5+4q_2^2 +2q_2q_4+8q_2q_5+4q_3^2+4q_3q_4+4q_4^2+2q_4q_5+4q_5^2\\
=4q_0+4q_0q_1+8q_0q_3+4q_0q_4+4q_1+2q_1q_2+4q_1q_3+8q_1q_4+2q_1q_5+4q_2 +2q_2q_4+8q_2q_5+4q_3+4q_3q_4+4q_4+2q_4q_5+4q_5

上記は混雑度ですが、もう一つ設定する必要があります。3つのルートのうちに1つをえらぶ制約条件です。

K(q_0+q_1+q_2-1)^2+K(q_3+q_4+q_5-1)^2 = -Kq_0+2Kq_0q_1+2Kq_0q_2-Kq_1+2Kq_1q_2-Kq_2-Kq_3+2Kq_3q_4+2Kq_3q_5-Kq_4+2Kq_4q_5-Kq_5+2K

試しに制約条件だけを解いてみます。

from blueqat import vqe
from blueqat.pauli import qubo_bit as q

K=2
cost2 = -K*q(0)+2*K*q(0)*q(1)+2*K*q(0)*q(2)-K*q(1)+2*K*q(1)*q(2)-K*q(2)-K*q(3)+2*K*q(3)*q(4)+2*K*q(3)*q(5)-K*q(4)+2*K*q(4)*q(5)-K*q(5)

step = 4

result = vqe.Vqe(vqe.QaoaAnsatz(cost2, step)).run()
print(result.most_common(12))

答えは、

(((0, 0, 1, 1, 0, 0), 0.10912346345562347), ((1, 0, 0, 1, 0, 0), 0.10912346345562343), ((0, 0, 1, 0, 0, 1), 0.10912346345562343), ((1, 0, 0, 0, 1, 0), 0.10912346345562342), ((1, 0, 0, 0, 0, 1), 0.10912346345562342), ((0, 1, 0, 1, 0, 0), 0.1091234634556234), ((0, 1, 0, 0, 1, 0), 0.1091234634556234), ((0, 0, 1, 0, 1, 0), 0.10912346345562339), ((0, 1, 0, 0, 0, 1), 0.10912346345562336), ((1, 0, 0, 0, 0, 0), 0.0015359617223092625), ((0, 0, 0, 1, 0, 0), 0.0015359617223092618), ((0, 0, 0, 0, 0, 1), 0.0015359617223092577))

上手くできました。この二つを繋げてQUBO式とします。

そして、実行します。

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comを利用する時に必要

from blueqat import vqe
from blueqat.pauli import qubo_bit as q

K=10
cost1 = 4*q(0)+4*q(0)*q(1)+8*q(0)*q(3)+4*q(0)*q(4)+4*q(1)+2*q(1)*q(2)+4*q(1)*q(3)+8*q(1)*q(4)+2*q(1)*q(5)+4*q(2)+2*q(2)*q(4)+8*q(2)*q(5)+4*q(3)+4*q(3)*q(4)+4*q(4)+2*q(4)*q(5)+4*q(5)
cost2 = -K*q(0)+2*K*q(0)*q(1)+2*K*q(0)*q(2)-K*q(1)+2*K*q(1)*q(2)-K*q(2)-K*q(3)+2*K*q(3)*q(4)+2*K*q(3)*q(5)-K*q(4)+2*K*q(4)*q(5)-K*q(5)

step = 8

result = vqe.Vqe(vqe.QaoaAnsatz(cost1+cost2, step)).run()
print(result.most_common(12))

なかなか答えが出ませんが、、、、

(((0, 1, 0, 0, 0, 1), 0.12986970172628198), ((0, 0, 1, 0, 1, 0), 0.1298697017262819), ((0, 0, 0, 0, 0, 0), 0.08232027492914641), ((0, 1, 0, 0, 1, 0), 0.05463053577530196), ((0, 0, 0, 1, 0, 0), 0.04609411831952029), ((1, 0, 0, 0, 0, 0), 0.04609411831952023), ((0, 0, 1, 1, 0, 1), 0.037740176136034545), ((1, 0, 1, 0, 0, 1), 0.0377401761360345), ((1, 0, 0, 0, 1, 1), 0.03131986281902896), ((0, 1, 1, 1, 0, 0), 0.031319862819028946), ((0, 0, 1, 0, 0, 1), 0.03071276141604334), ((1, 0, 0, 1, 0, 0), 0.02140326399538612))

うまくいけば制約条件を満たす解を与えます。VQEセミナーの中ですでにQAOAに触れていましたので、今後は統合しようと思います。

2-5-5 Quantum Alternating Operator Ansatz

上記はなかなか答えが出ませんでした。ここでは、正答率を大きく上げる方法を確認します。ハミルトニアンの制約条件を定式化から外すできます。前述の交通最適化問題などは3つのルートから1つのルートを選ぶという制約式が正答率を下げる原因でしたが、今回はその制約式を減らすテクニックを確認します。

XX+YYゲートを導入することで、01と10を入れ替える操作をできます。時間発展を見ても、入れ替え操作が見れます。最初から00と11をださず、01と10のみで計算をします。


引用:https://qiita.com/gyu-don/items/c51a9e3d5d16a6d5baf6

これらの操作を導入して簡単な問題を解いてみます。

2-5-6 Quantum Alternating Operator Ansatzの例題

今回は初期状態に|10>と|01>のもつれ状態を採用し、途中のmixerにXYミキサーというXX+YYの時間発展を導入したものを利用してみます。

import numpy as np
from blueqat import Circuit
from blueqat.pauli import X, Y, Z, I
from blueqat.pauli import qubo_bit as q
from blueqat.vqe import AnsatzBase, Vqe

def an(index):
    return 0.5 * X[index] + 0.5j * Y[index]

def cr(index):
    return 0.5 * X[index] - 0.5j * Y[index]

op = (cr(1) * an(0) + cr(0) * an(1)).to_expr().simplify()

class QaoaQaoaAnsatz(AnsatzBase):
    def __init__(self, hamiltonian, step=1, init_circuit=None):
        super().__init__(hamiltonian, step * 2)
        self.hamiltonian = hamiltonian.to_expr().simplify()
        if not self.check_hamiltonian():
            raise ValueError("Hamiltonian terms are not commutable")

        self.step = step
        self.n_qubits = self.hamiltonian.max_n() + 1
        if init_circuit:
            self.init_circuit = init_circuit
            if init_circuit.n_qubits > self.n_qubits:
                self.n_qubits = init_circuit.n_qubits
        else:
            self.init_circuit = Circuit(self.n_qubits).h[:]
        self.init_circuit.make_cache()
        self.time_evolutions = [term.get_time_evolution() for term in self.hamiltonian]
        self.time_evolutions2 = [term.get_time_evolution() for term in op]
    def check_hamiltonian(self):
        """Check hamiltonian is commutable. This condition is required for QaoaAnsatz"""
        return self.hamiltonian.is_all_terms_commutable()

    def get_circuit(self, params):
        c = self.init_circuit.copy()
        betas = params[:self.step]
        gammas = params[self.step:]
        for beta, gamma in zip(betas, gammas):
            beta *= np.pi
            gamma *= 2 * np.pi
            for evo in self.time_evolutions:
                evo(c, gamma)
            for evo2 in self.time_evolutions2:
                evo2(c, beta)
        return c

h = -3*q(0)-3*q(1)-2*q(0)*q(1)
h = h.to_expr().simplify()
runner = Vqe(QaoaQaoaAnsatz(h,4,Circuit().h[0].cx[0,1].x[0]))
result = runner.run()

# get probability
print(result.most_common(12))

print('Result by QAOA')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))

# Hamiltonian to matrix
mat = h.to_matrix()

# Calculate by numpy
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])

通常、最小基底状態は-8を選びますが、上記は|01>と|10>のもつれ状態を利用しているため、制約条件をハード制約としてもち、|00>と|11>をとることができないため、QAOAで得られる最小値の期待値は-3となっています。

(((1, 0), 0.4999999999999993), ((0, 1), 0.4999999999999993), ((0, 0), 4.0059342843254495e-32), ((1, 1), 3.0814879110195774e-32))
Result by VQE
-2.999999999999996
Result by numpy
-8.0

このように、通常ソフト制約として出てくるはずのハード制約をかすことで、正答率を大幅に上げることができます。回路自体は長くなっています。

スクリーンショット 2020-02-25 0.39.32.png

2-5-7 量子もつれを使った制約条件項なしの交通最適化計算

ケーススタディとして交通最適を取り上げます。
量子アニーリングでは探索はほぼランダムで行われます。初期状態は|+>からスタートし、途中過程の探索は$\sigma_x$の横磁場を使います。

しかしここでは、ルート選択は1つだけを選ぶわけですから、$q_0$か$q_1$の片方は1で、もう片方は0となるのが条件としてついてきます。また、$q_2$と$q_3$も同様です。

ここで、量子回路でまず量子もつれを作ります。量子もつれは多くのパターンから限られた組み合わせだけが出るように調整できるのでした。簡単な例として、

from blueqat import Circuit
Circuit().h[0].cx[0,1].m[:].run(shots=100)

こうすると、本来4種類出るところから00と11のみ答えが出ます。

Counter({'00': 51, '11': 49})

ちょっと形を変えることで、01と10のもつれに変更できます。

from blueqat import Circuit
Circuit().h[0].cx[0,1].x[0].m[:].run(shots=100)

#=>Counter({'01': 41, '10': 59})

これを使います。つまり、自動車1と自動車2にそれぞれもつれを活用して、01と10に制限をかけるように初期状態を設定します。

今回4量子ビットの回路では、

from blueqat import Circuit
Circuit().h[0].cx[0,1].x[0].h[2].cx[2,3].x[2].m[:].run(shots=100)

#=>Counter({'0101': 27, '0110': 16, '1001': 29, '1010': 28})

この4通りに絞ります。

入れ替え作業

$\sigma_x$の横磁場での探索の代わりに明示的に(XX+YY)/2というゲートを利用することで、01と10を入れ替えることができます。

これにより、初期状態で準備した01と10をひたすら入れ替えながら、量子断熱計算で、ハミルトニアンの最小値を求めるQAOAと組み合わせることができます。

具体的コード

具体的な実行コードはこちらです。元々のQAOAのAnsatzを改造しています。詳しい内容は、

「【ほぼ決定版】量子コンピュータVQE/QAOAセミナーまとめ」
https://qiita.com/YuichiroMinato/items/3f68b9fed67ff3a01370

こちらをご覧ください。

import numpy as np
from blueqat import Circuit
from blueqat.pauli import X, Y, Z, I
from blueqat.pauli import qubo_bit as q
from blueqat.vqe import AnsatzBase, Vqe

def an(index):
    return 0.5 * X[index] + 0.5j * Y[index]

def cr(index):
    return 0.5 * X[index] - 0.5j * Y[index]

op1 = (cr(1) * an(0) + cr(0) * an(1)).to_expr().simplify()
op2 = (cr(3) * an(2) + cr(2) * an(3)).to_expr().simplify()

class QaoaQaoaAnsatz(AnsatzBase):
    def __init__(self, hamiltonian, step=1, init_circuit=None):
        super().__init__(hamiltonian, step * 2)
        self.hamiltonian = hamiltonian.to_expr().simplify()
        if not self.check_hamiltonian():
            raise ValueError("Hamiltonian terms are not commutable")

        self.step = step
        self.n_qubits = self.hamiltonian.max_n() + 1
        if init_circuit:
            self.init_circuit = init_circuit
            if init_circuit.n_qubits > self.n_qubits:
                self.n_qubits = init_circuit.n_qubits
        else:
            self.init_circuit = Circuit(self.n_qubits).h[:]
        self.init_circuit.make_cache()
        self.time_evolutions = [term.get_time_evolution() for term in self.hamiltonian]
        self.time_evolutions1 = [term.get_time_evolution() for term in op1]
        self.time_evolutions2 = [term.get_time_evolution() for term in op2]
    def check_hamiltonian(self):
        """Check hamiltonian is commutable. This condition is required for QaoaAnsatz"""
        return self.hamiltonian.is_all_terms_commutable()

    def get_circuit(self, params):
        c = self.init_circuit.copy()
        betas = params[:self.step]
        gammas = params[self.step:]
        for beta, gamma in zip(betas, gammas):
            beta *= np.pi
            gamma *= 2 * np.pi
            for evo in self.time_evolutions:
                evo(c, gamma)
            for evo1 in self.time_evolutions1:
                evo1(c, beta)
            for evo2 in self.time_evolutions2:
                evo2(c, beta)
        return c

h = 4*q(0)+4*q(1)+4*q(2)+4*q(3)+4*q(0)*q(1)+4*q(0)*q(2)+8*q(1)*q(2)+2*q(1)*q(3)+2*q(2)*q(3)
h = h.to_expr().simplify()
runner = Vqe(QaoaQaoaAnsatz(h,4,Circuit().h[0].cx[0,1].x[0].h[2].cx[2,3].x[2]))
result = runner.run()

# get probability
print(result.most_common(12))

print('Result by QAOA')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))

# Hamiltonian to matrix
mat = h.to_matrix()

# Calculate by numpy
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])

これを実行すると、

(((1, 0, 0, 1), 0.9285092517501069), ((0, 1, 0, 1), 0.03727452567612003), ((1, 0, 1, 0), 0.019828254090636977), ((0, 1, 1, 0), 0.01438796848312899), ((1, 0, 0, 0), 4.949639957075196e-32), ((1, 1, 1, 0), 4.1497989458064e-32), ((1, 0, 1, 1), 3.9481563859938324e-32), ((0, 0, 1, 0), 2.7174584187380976e-32), ((1, 1, 0, 1), 9.7798953112461e-33), ((0, 0, 0, 1), 6.384745093566593e-33), ((0, 1, 0, 0), 6.162975822039154e-33), ((0, 1, 1, 1), 4.240055143190024e-33))

Result by QAOA
8.268965815579767

こうなりました。93%近くの確率で制約条件を満たしながら混雑の最も少ないルートが選択されています。

2-6 中級設定

変分アルゴリズムを使い続けていると、様々な設定をしたくなってきます。今回は古典最適化部分を見てみたいと思います。

2-6-1 古典最適化

変分アルゴリズムの内部構造は大きく分けて、

1、パラメータ付きの短い量子回路
2、量子回路を集計し、期待値を計算、次のパラメータの提案

です。次のトライアルのパラメータを渡すのは最適化アルゴリズムですが、それは古典で実行されます。

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comを利用する時に必要
import numpy as np
from blueqat import Circuit
from blueqat.pauli import X, Y, Z, I
from blueqat.pauli import qubo_bit as q
from blueqat.vqe import AnsatzBase, Vqe

class QubitAnsatz(AnsatzBase):
    def __init__(self, hamiltonian):
        super().__init__(hamiltonian, 4)
        self.step = 1

    def get_circuit(self, params):
        a, b, c, d = params
        return Circuit().ry(a)[0].rz(b)[0].ry(c)[1].rz(d)[1]

h = -3*q(0)-3*q(1)-2*q(0)*q(1)
h = h.to_expr().simplify()
minimizer=vqe.get_scipy_minimizer(method="COBYLA",options={"tol":5.0e-4})
runner = Vqe(QubitAnsatz(cost),minimizer=minimizer)
result = runner.run()

print('Result by VQE')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))

# Hamiltonian to matrix
mat = h.to_matrix()

# Calculate by numpy
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])

上記では古典最適化のアルゴリズムを指定できています。主にscipyから選べます。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

'Nelder-Mead', 'Powell', 'CG', 'BFGS', 'L-BFGS-B', 'TNC', 'COBYLA', 'SLSQP'

この辺りを見れば良さそうです。自分の好きなソルバーを組み込むことができます。

def hyperopt_minimizer(objective, n_params):
    from hyperopt import fmin, Trials, tpe, hp
    trials = Trials()
    best = fmin(objective, [hp.uniform(f'p{i}', 0., 2 * np.pi) for i in range(n_params)],
            algo=tpe.suggest, max_evals=10000, trials=trials, verbose=1)
    return list(best.values())

2-6-2 numbaバックエンドの採用

単純に高速化できます。元は、

from blueqat import vqe
from blueqat.pauli import qubo_bit as q

hamiltonian = q(0)*q(1)*q(2)*q(3)
step = 10

start = time.time()
result = vqe.Vqe(vqe.QaoaAnsatz(hamiltonian, step)).run()
end = time.time()

print(end - start)

こちらはハミルトニアンを作ってvqeでstep = 10で実行します。

32.705730676651
from blueqat import BlueqatGlobalSetting
BlueqatGlobalSetting.set_default_backend('numba')

上記を設定します。

hamiltonian = q(0)*q(1)*q(2)*q(3)
step = 10

start = time.time()
result = vqe.Vqe(vqe.QaoaAnsatz(hamiltonian, step)).run()
end = time.time()

print(end - start)

こちらは、

8.867603778839111

3、4倍ほど速くなりました。

詳しくはこちらをみてください。
「blueqatでいろいろベンチマークとってみる。最適化とかGPUとか。」
https://qiita.com/YuichiroMinato/items/f08420957896b96e6874

2-6-3 Qiskitでもみてみる

IBMのqiskitも参考に見てみます。インストールは、

pip install qiskit-aqua

Qiskitはちょっと長いです。

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)

ここで、上記はmaxcutを解いていますが、

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

のところで、最適化アルゴとansatzを準備してVQEを解いているのが分かります。どのツールでも基本的には使っている回路や最適化は大きくは変わらないので、安心して量子コンピュータを学びましょう。

参考文献

https://arxiv.org/pdf/1304.3061.pdf

「量子コンピュータの基礎(2)」
https://qiita.com/KeiichiroHiga/items/f5117efcf6a96c9d10b8

「1量子ビットVQEを世界一簡単に書く 2020年2月Ver.」
https://qiita.com/gyu-don/items/65f830ae3777534f9a04

「Quantum Computing at Volkswagen:
Traffic Flow Optimization using the D-Wave Quantum Annealer」
引用:https://www.dwavesys.com/sites/default/files/VW.pdf

「交通最適化問題を量子アニーリングで解く」
https://qiita.com/sychocola1/items/646d2c9504e34ed1589f

Variational Quantum Factoring
Eric R. Anschuetz, Jonathan P. Olson, Alán Aspuru-Guzik, Yudong Cao
(Submitted on 27 Aug 2018)
https://arxiv.org/abs/1808.08927

「[QAOA]SWAPの時間発展で|10>と|01>とを入れ替える」
https://qiita.com/gyu-don/items/c51a9e3d5d16a6d5baf6

「blueqatでいろいろベンチマークとってみる。最適化とかGPUとか。」
https://qiita.com/YuichiroMinato/items/f08420957896b96e6874

mdrft
量子コンピュータのアプリケーション、ミドルウェア、ハードウェアをフルスタックで開発
https://mdrft.com/
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
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  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
ユーザーは見つかりませんでした