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

【ほぼ決定版】量子コンピュータVQE/QAOAセミナーまとめ

はじめに

量子コンピュータ勉強会で、変分アルゴリズムを勉強して2周終わりました。さらに内容をアップデートして、3周目をまとめます。変分アルゴリズムは2020年現在で利用できる量子コンピュータ向けのアルゴリズムを理解する上でとても大事です。ここでは理論から応用までをできるだけ幅広く丁寧にフォローしたいと思います。

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

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

2-1-1 VQE/QAOAとは?

VQE(Variational Quantum Eigensolver)は、2020年現在の量子コンピュータはエラーが多く長い回路の実行が厳しいため従来想定された位相推定アルゴリズムの実行が厳しく、その代替としてエルミート行列の固有値の期待値を求める量子古典ハイブリッドアルゴリズムとして、2013年に開発されました。様々な問題に応用でき、量子化学計算や組合せ最適化問題を解くことができます。

また、定式化と量子回路を特に組合せ最適化問題に特化した形で発展したのがQAOA(Quantum Approximate Optimization Algorithm)で、基本の形はVQEを踏襲しながらも独自の機能特化をしています。

両方のアルゴリズムは変分法を使うアルゴリズムとして量子コンピュータの基本を構成しています。

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

量子コンピュータで利用される量子ビットは|0>状態と|1>状態を重ね合わせて量子状態を表現します。

\mid \psi \rangle = \alpha \mid 0 \rangle + \beta \mid 1 \rangle

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

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

を満たします。$|\alpha|^2$や$|\beta|^2$は「測定」をすることで取得できます。例えば、100回計算して、0が45回、1が55回でたら、$|\alpha|^2$の期待値は0.45、$|\beta|^2$の期待値は0.55と計算できます。

今回、状態ベクトルは縦ベクトルで表現をします。

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

より、任意の1量子ビットの量子状態は、

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

と表されます。

2-1-3 下準備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ゲートです。アダマールゲートは、

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-4 下準備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

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

VQEやQAOAなどの変分アルゴリズムの計算手順は、

1、パラメータ化された量子回路(ansatz)を準備(量子計算)
2、量子ビットの測定値からハミルトニアン(固有値を求めたい行列)の期待値を計算(古典計算)
3、古典最適化計算により、ハミルトニアンの期待値を採用にするパラメータを探索(古典計算)

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

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

ansatz(アンザッツ)と呼ばれる、状態ベクトルを作るために任意回転ゲートと角度でパラメータ化された量子回路を組み立て、ハミルトニアンの期待値を求めるための測定値を提供します。後ほど確認しますが、ansatz回路はときたい問題によって量子回路の形が大きく変わり、変分アルゴリズムで大事な要素の1つです。

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

ハミルトニアンと呼ばれる演算子を準備しますが、これは基本のパウリゲートXYZと単位行列Iから構成される行列です。今回の変分アルゴリズムは、このハミルトニアンの固有状態である状態ベクトルから導き出される固有値の期待値を求めることを目的としているため、ハミルトニアンの構築も大事な要素の1つです。

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

ハミルトニアンの期待値は、状態ベクトルとハミルトニアンが準備されているとき人、計算によってもとまります。ある状態ベクトル$\mid \psi \rangle$とハミルトニアン$H$が与えられた時のハミルトニアンの期待値は、

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

で表され、単純な計算によってもとまります。

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

ansatzを使って任意の状態ベクトルが作れたとして、求めたいハミルトニアンが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が出る確率に対応しますので、最小になるのは、$|\beta|$=1となるときの-1になります。

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

ハミルトニアンがXの場合には、

\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 ハミルトニアンの変形

ハミルトニアンの要素が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 = 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

とすることができます。これらは、実際に状態ベクトルを組み込んだゲート操作を測定の直前に行うことで解を得ます。

2-1-12 線形結合

ハミルトニアンは下記のような場合、分解できます。

\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 計算はツールがやってくれる

多くの理論の計算は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))

変分アルゴリズム#2-2(実装初級編)

学んだ内容を元にプログラムに落とします。ソフトウェアの導入からです。

2-2-1 インストール

blueqatsdkを使います。

https://github.com/Blueqat/Blueqat

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

インストールは、

pip3 install blueqat

で完了します。quantumcomputing.comというサービスを利用するには登録が必要です。

https://app.quantumcomputing.com/

2020年初頭現在はβ運用中で今年製品版が出ます。無料です。こちらは最初からblueqatが入っているのでインストール不要です。

2-2-2 コード

下記のコードを入れてみます。今回は1量子ビットの行列の固有値を求めてみます。

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.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]


# Calculate by VQE
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))

# Hamiltonian to matrix
mat = h.to_matrix()

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

上記は、hで決められたハミルトニアンの期待値を求め、その答え合わせをnumpyで確かめをしています。

上記のポイントは2カ所で、

  1. Circuit().ry(a)[0].rz(b)[0]のansatz
  2. h = 1.23 * I - 4.56 * X(0) + 2.45 * Y(0) + 2.34 * Z(0)のハミルトニアン

です。

2-2-3 ハミルトニアンの定式化

ハミルトニアンの行列はパウリゲートを使って自分で作る必要があります。blueqatの場合、量子化学用のUCCansatzと呼ばれるものと、組合せ最適化問題のQAOAは簡単に利用できます。今回はどちらも使わず手作業でやってみます。

適当に係数とパウリゲートを組み合わせて作ってみました。

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

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

上記は、

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

の例題からですが、IとXとYとZというゲートを使って2*2の行列を作っています。XやYのあとの0という数字は何番目の量子ビットを使うかという通し番号になっています。量子ビットの通し番号は通常0からスタートします。

上記の行列は参考として、

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

となっています。

2-2-4 Ansatzの準備

ansatzの準備は、

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]

この部分で、Circuit().ry(a)[0].rz(b)[0]がansatzです。0番目の量子ビットに、RyゲートとRzゲートを利用して状態ベクトルが表現されており、aとbが角度パラメータです。今回は主にこれを計算します。

2-2-5 結果

量子回路が実行され、ハミルトニアンの期待値がもとまり、最適化が進むと結果は出力されます。numpyでの念のための確かめもできます。

Result by VQE
-4.450591481563331
Result by numpy
-4.450818602983201

期待値なので、ほぼ近い値が出ています。

2-2-6 quantumcomputing.com

今回のコードの共有リンクを準備しました。

https://app.quantumcomputing.com/public/projects/udFDCw

実行すると、

式の途中にパラメータが使われていて、上記で設定したansatzのパラメータが最適化されて答えが出ています。

VQE#2-3(量子化学計算H2編)

量子化学計算はハミルトニアンの作成をパウリ行列を使って行い、VQEを利用して解を得ることができます。

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 量子コンピュータ向け変換

ブラビキタエフ変換やジョルダンウィグナー変換がありますが、今回は前者の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]

このようにきちんと変換されました。ちょっと複雑な形をしていますが、これこそがこれまで学んできたパウリ行列の和で書かれたハミルトニアンです。そして、このhをVQEにかけます。

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: [1.00803084 0.88310007] val: 0.10544430487391798
params: [-1.61000316  0.88310007] val: 0.48593858308326665
params: [0.00803084 0.88310007] val: -0.3685674394653079
params: [-0.61000314  0.88310007] val: -0.24619628358384726
params: [0.38999684 0.88310007] val: -0.27563514593504546
params: [-0.07664849  0.88310007] val: -0.371744879738929
params: [-0.07664849  0.88310007] val: -0.371744879738929
params: [-0.07664849  1.88310007] val: 1.0541144286439048
params: [-0.07664849 -0.73493393] val: -0.44860542849648677
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [-0.16132782 -0.82515946] val: -0.32853423118712466
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [0.92335151 0.0289703 ] val: -0.3234253290431143
params: [-1.69468249  0.0289703 ] val: 0.7219628224155241
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [-0.69468246  0.0289703 ] val: -0.5631470828007092
params: [0.30531751 0.0289703 ] val: -0.8360012977529356
params: [-0.0024937  0.0289703] val: -0.903736928997239
params: [-0.0024937  0.0289703] val: -0.903736928997239
params: [-0.0024937  1.0289703] val: -0.19293944513988606
params: [-0.0024937 -1.5890637] val: 0.5730186565494813
params: [-0.0024937  0.0289703] val: -0.903736928997239
params: [-0.0024937  -0.58906367] val: -0.6556036819150325
params: [-0.0024937  0.4109363] val: -0.7807595396151468
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(組合せ最適化イジング+QUBO編)

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

1、まずプログラムを覚える
2、壁に当たったら理論に戻る

という順でやることを想定します。疑問があったらslackの量子コンピュータコミュニティでいつでも質問できます。

https://join.slack.com/t/blueqat/shared_invite/enQtNzc3NjgxNjA5OTA2LTNmNjg0YzAzNGUxYmYwNDIxNWFlZGVkODI0MTE2NDMwM2I3OTRjNGExZjMxZWFmZDI3ZTRlNDU3NTcwZDUyOTk

2-4-1 組合せ最適化問題

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

2-4-2 組合せ最適化問題の定式化

組合せ最適化問題の定式化はハミルトニアンをコスト関数として扱い、答えを解いて最小値を求めたときに最適な答えになるように調整ができます。ルールは、

・問題をイジングモデルと呼ばれるモデルに落とし込む
・Z演算子を使う

ということです。

今回は下記の例題を実行してみます。

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

ansatzを含む全体のコードは、

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.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が答えとして出てきました。ルールを説明します。

2-4-3 Zは1か-1をとる

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

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

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

が設定されています。

Zは期待値として-1か+1のどちらかをとります。
また、hはより小さい値を取ると正解になります。
上記VQEでは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パラメータ利用しました。

2-4-4 この問題を大きくするだけで社会問題がたくさん解ける

そして実際にはより大きな問題を解きます。上記問題は役にたたなさそうですが、Z(0)をAさん、Z(1)をBさんに見立てて、Aさんはグループ1に属し、BさんはAさんと同じグループに属するという条件をつけた分類問題と同じです。

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

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

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

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

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

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

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

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

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

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])

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

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が得られます。

Result by VQE
-7.981654210576551
Result by numpy
-8.0

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

変分アルゴリズム#2-5(QAOA編)

変分アルゴリズムでは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_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

試しに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ライブラリを実行します。

今回は簡単な定式化を解いてみましょう。

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))

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

上記で設定する必要があるのはコスト関数です。

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のウェイトを設定します。

これをQAOAで解くことで、最小値-8が得られます。

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

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

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

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

youtubeでは、
交通最適化
https://www.youtube.com/watch?v=JU2omOzuOQg

がわかりやすいです。記事は、

「D-WaveでVW社の交通最適化アプリケーションの実装を解く」
https://qiita.com/YuichiroMinato/items/19dd89e381f14ad24134

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

今回は、交通最適をときます。やることは、

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

一つ目は混雑度をコスト関数として、

(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))

この二つを繋げて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-4 実例:最小値問題における素因数分解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-5 Quantum Alternating Operator Ansatz

初期状態とmixerをswap時間発展へと変更することで、ハミルトニアンの制約条件を減らすことができます。前述の交通最適化問題などは3つのルートから1つのルートを選ぶという制約式が正答率を下げる原因でしたが、今回はその制約式を減らすテクニックを確認します。

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

XX+YYゲートを導入することで、01と10を入れ替える操作を通じてmixerを作れます。

from blueqat import Circuit, pauli, vqe
from blueqat.pauli import qubo_bit as q
from math import pi
import numpy as np

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

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

op = (cr(1) * an(0) + cr(0) * an(1)).to_expr().simplify()
print(op)
x = []
y00 = []
y01 = []
y10 = []
y11 = []
evos = [term.get_time_evolution() for term in op.terms]
for i in range(100):
    c = Circuit().h[0].cx[0, 1].x[1]
    for evo in evos:
        evo(c, i / 100 * pi)
    c.rz(i / 100 * pi)[1]
    for evo in evos:
        evo(c, i / 100 * pi)
    c.rz(i / 100 * pi)[1]
    state = c.run()
    _a, _b, _c, _d = list(np.abs(state) ** 2)
    x.append(i)
    y00.append(_a)
    y10.append(_b)
    y01.append(_c)
    y11.append(_d)

import matplotlib.pyplot as plt
plt.plot(x, y00, label='|00>')
plt.plot(x, y10, label='|10>')
plt.plot(x, y01, label='|01>')
plt.plot(x, y11, label='|11>')
plt.legend()
plt.show()

時間発展を見ても、入れ替え操作が見れます。

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

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

quantumcomputing.comのリンクも貼っておきます。

https://app.quantumcomputing.com/public/projects/GMsbjC

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%近くの確率で制約条件を満たしながら混雑の最も少ないルートが選択されています。
(1,0,0,1)

変分アルゴリズムセミナー#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])

上記では古典最適化のアルゴリズムを指定できています。

2-6-2 最適化アルゴリズムの選択

主に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-3 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-4 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を解いているのが分かります。どのツールでも基本的には使っている回路や最適化は大きくは変わらないので、安心して量子コンピュータを学びましょう。

最後に

一通り変分アルゴリズムを見てみました。この基本を抑えれば様々な問題が解けます。

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
ユーザーは見つかりませんでした