common.title
Cloud support

Nobisuke

Dekisugi

RAG

autoQAOA
RAG for dev
Fortune telling app
Annealing
DEEPSCORE
Translation
Overview
Service overview
Terms of service

Privacy policy

Contact
Research
Sign in
Sign up
common.title

量子コンピュータデータサイエンティスト講座：VQE

Yuichiro Minato

2021/02/21 07:50

#量子ゲート

0.はじめに

量子コンピュータを始める人に向け、実践的な量子古典ハイブリッド計算の勉強会をしています。

0-1.インストール

今回はblueqatを利用します。

0-2.概要

量子アルゴリズムには将来的な誤り訂正が搭載されたときに利用される「汎用アルゴリズム」と、誤り訂正なしで利用される「変分アルゴリズム」があります。

・汎用アルゴリズム（グローバー、ショア、位相推定、量子フーリエ変換、HHL、QSVMなど）
・変分アルゴリズム（VQE,QAOAなどの量子古典ハイブリッドアルゴリズム）

ここでは実用性を重視して、量子古典ハイブリッドアルゴリズムを学びたいと思います。量子古典ハイブリッドアルゴリズムの代表格VQE(Variational Quantum Eigensolver)は、「位相推定」アルゴリズムの代替として2013年に当時ハーバード大学（現在はトロント大学）のアラン・アスプル・グジック教授のチームによって開発されました。

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

vqe1

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

1.理論編

最初に理論を確認します。最終的にはどんどんツール化されていますが、内容を確認します。

1-1.固有値・固有ベクトル問題

ある行列に対する固有値と固有ベクトルを見つけることができれば様々な問題を解くことができます。与えられた行列HH（ハミルトニアンと呼ばれます）に対して、固有値をE0E_0、固有ベクトルを∣ψ⟩\mid \psi \rangleとしたときに、

H∣ψ⟩=E0∣ψ⟩H \mid \psi \rangle = E_0 \mid \psi \rangle

を満たすような、E0E_0と∣ψ⟩\mid \psi \rangleを見つけるのが目的です。∣ψ⟩\mid \psi \rangleは状態ベクトル（波動関数）と呼ばれます。

1-2.量子状態と状態ベクトル

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

∣0⟩=[10],∣1⟩=[01]\mid 0 \rangle = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \mid 1 \rangle = \begin{bmatrix} 0 \\ 1 \end{bmatrix}

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

∣ψ⟩=α∣0⟩+β∣1⟩=α[10]+β[01]=[αβ]\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}

また、下記条件を満たします。

∣α∣2+∣β∣2=1|\alpha|^2+|\beta|^2 = 1

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

1-3.量子状態の極座標表示

ブロッホ球は二つの直交する純粋状態の重ね合わせで表現できる量子状態を単位球面上に表す表記法。
https://ja.wikipedia.org/wiki/%E3%83%96%E3%83%AD%E3%83%83%E3%83%9B%E7%90%83

vqe2

∣ψ⟩=[cosθ−sinθe−iϕsinθeiϕcosθ][10]=[cosθsinθeiϕ]\mid \psi \rangle = \begin{bmatrix} cos\theta & -sin\theta e^{-i\phi}\\ sin\theta e^{i\phi} & cos\theta \end{bmatrix} \begin{bmatrix} 1\\0 \end{bmatrix} = \begin{bmatrix} cos\theta \\ sin\theta e^{i\phi} \end{bmatrix}

のように角度を利用したパラメータ表記もできます。

1-4.状態ベクトルの操作と量子ゲート

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

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

X=[0110],Y=[0−ii0],Z=[100−1]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(θ)=[cos⁡(θ2)−isin⁡(θ2)−isin⁡(θ2)cos⁡(θ2)],Ry(θ)=[cos⁡(θ2)−sin⁡(θ2)sin⁡(θ2)cos⁡(θ2)],Rz(θ)=[e−iθ200eiθ2]Rx(\theta) = \begin{bmatrix} \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{bmatrix}, Ry(\theta) = \begin{bmatrix} \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{bmatrix}, Rz(\theta) = \begin{bmatrix} e^{-i\frac{\theta}{2}} & 0\\ 0 & e^{i\frac{\theta}{2}} \end{bmatrix}

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

X[10]=[0110][10]=[01]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}

1-5.複数量子ビットの状態ベクトルとテンソル積

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

[ab]⊗[cd]=[a∗[cd]b∗[cd]]=[acadbcbd]\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}

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

[ab]⊗[cd]⊗[ef]=[ac∗[ef]ad∗[ef]bc∗[ef]bd∗[ef]]=[aceacfadeadfbcebcfbdebdf]\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}

量子ゲートもテンソル積をとれます。ハミルトニアンの計算の際に利用します。ゲートがない場合には単位行列Iを組み込みます。

[abcd]⊗[efgh]=[a∗[efgh]b∗[efgh]c∗[efgh]d∗[efgh]]=[aeafbebfagahbgbhcecfdedfcgchdgdh]\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}

1-6.ハミルトニアンとその期待値

ハミルトニアンHHはパウリゲートXYZと単位行列Iから構成されるエルミート行列です。ある量子状態∣ψ⟩\mid \psi\rangleが与えられた場合の、ハミルトニアンHHの期待値EEは、上記の式の両方に⟨ψ∣\langle \psi \midを左からかけてみると、

⟨ψ∣H∣ψ⟩=⟨ψ∣E∣ψ⟩⟨ψ∣H∣ψ⟩=E⟨ψ∣ψ⟩⟨ψ∣H∣ψ⟩=E\langle \psi \mid H \mid \psi \rangle = \langle \psi \mid E \mid \psi \rangle\\ \langle \psi \mid H \mid \psi \rangle = E \langle \psi \mid \psi \rangle\\ \langle \psi \mid H \mid \psi \rangle = E

のように計算できます。1量子ビットの例として、任意の量子状態∣ψ⟩\mid \psi \rangleと任意のハミルトニアンHHから、

⟨ψ∣H∣ψ⟩=[α∗β∗][abcd][αβ]\langle \psi \mid H \mid \psi \rangle = \begin{bmatrix} \alpha^* & \beta^* \end{bmatrix} \begin{bmatrix} a&b\\ c&d \end{bmatrix} \begin{bmatrix} \alpha\\ \beta \end{bmatrix}

のように計算できます。

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

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

⟨ψ∣H∣ψ⟩=[α∗β∗][100−1][αβ]=∣α∣2−∣β∣2\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

となります。∣α∣2|\alpha|^2は0がでる確率、∣β∣2|\beta|^2は1が出る確率に対応します。確かめてみましょう。

numpyを導入して、状態ベクトルを設定します。

import numpy as np
v_a = np.array([0.9387913 +0.j, -0.23971277+0.j, 0. +0.06120872j, 0. -0.23971277j])
v_a
array([ 0.9387913 +0.j        , -0.23971277+0.j        ,

        0.        +0.06120872j,  0.        -0.23971277j])

複素共役を取ります。

v_b = np.conjugate(v_a)
v_b
array([ 0.9387913 -0.j        , -0.23971277-0.j        ,

        0.        -0.06120872j,  0.        +0.23971277j])

求めたいのは、v_b@Z@v_aになりますので、素直に計算します。ZとIを準備して、

Z = np.array([[1,0],[0,-1]])
I = np.eye(2)

ZとIのテンソル積をとった行列をbとaで挟んで計算します。

v_b@np.kron(Z,I)@v_a
(0.8775825975516516+0j)

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

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

⟨ψ∣H∣ψ⟩=[α∗β∗][0110][αβ]=α∗β+αβ∗\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^*

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

1-9.ハミルトニアンの変形

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

X=HZHX = HZH

より、

⟨ψ∣X∣ψ⟩=⟨ψ∣HZH∣ψ⟩=⟨ψ′∣Z∣ψ′⟩\langle \psi \mid X \mid \psi \rangle \\ = \langle \psi \mid HZH \mid \psi \rangle\\ = \langle \psi' \mid Z \mid \psi' \rangle

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

Y=RX(−π/2)ZRX(π/2)Y = RX(-\pi/2) Z RX(\pi/2)

と、なり同様に、

⟨ψ∣Y∣ψ⟩=⟨ψ∣RX(−π/2)ZRX(π/2)∣ψ⟩=⟨ψ′′∣Z∣ψ′′⟩\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ゲートを入れて測定を行うことで期待値を計算することができます。

1-10.線形結合

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

⟨ψ∣aH1+bH2∣ψ⟩=⟨ψ∣aH1∣ψ⟩+⟨ψ∣bH2∣ψ⟩=a⟨ψ∣H1∣ψ⟩+b⟨ψ∣H2∣ψ⟩\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.2X0Z2+2.5Z0X1Y2−3.4Z2X1H = 1.2 X_0 Z_2 + 2.5 Z_0 X_1 Y_2 - 3.4 Z_2 X_1

期待値は、

⟨ψ∣1.2X0Z2+2.5Z0X1Y2−3.4Z2X1∣ψ⟩=1.2⟨ψ∣X0Z2∣ψ⟩+2.5⟨ψ∣Z0X1Y2∣ψ⟩−3.4⟨ψ∣Z2X1∣ψ⟩\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

のように個別に計算をして和を求めることで計算できます。

1-11.量子変分原理

次にハミルトニアンの期待値の計算がわかったところで、今度は計算手法を確認します。量子変分原理より、任意の状態ベクトルにおけるハミルトニアンの期待値は基底状態を上回るというのがあります。

⟨ψ(θ)∣H∣ψ(θ)⟩≥E0\langle \psi (\theta) \mid H \mid \psi (\theta) \rangle \geq E_0

固有値を直接求めるのは厳しいので、間接的に固有値の期待値にできるだけ近づける最適化計算に落とし込みます。ハミルトニアンHHは与えられた式なので、私たちにできるのは状態ベクトル∣ψ(θ)⟩\mid \psi (\theta) \rangleをできるだけ固有状態に近づける状態ベクトルを角度の最適化計算として探していくのが今回の目的です。状態ベクトルの探索はansatz（アンザッツ）と呼ばれる回路を作ることで実現します。

1-12.ansatz（アンザッツ）

上記の∣ψ⟩\mid \psi \rangleは量子回路によって作ることができます。その回路のことをansatz（アンザッツ）といいます。ansatzは通常極座標で量子状態を表現する角度パラメータを持った回路です。

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

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

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

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

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

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

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

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

vqe3.png

1-13.変分アルゴリズムの計算手順

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

１、パラメータθ\thetaを導入した量子回路（ansatz）を準備（量子計算） ２、量子ビットの測定値からハミルトニアンHHの期待値⟨ψ(θ)∣H∣ψ(θ)⟩\langle \psi (\theta) \mid H \mid \psi (\theta) \rangleを計算（古典計算） ３、古典最適化計算により、ハミルトニアンHHの期待値を最小にするパラメータθ\thetaを探索（古典計算）

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

vqe0.png

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

全体の実装は、

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().rx(a)[0].rz(b)[0]

# hamiltonian
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_sparse(result.circuit))

# This is for check
mat = h.to_matrix()
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])
Result by VQE

-4.450744388205187

Result by numpy

-4.450818602983201

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

2.VQEの実問題

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

１、量子化学計算
２、組合せ最適化問題

になります。

2-1.VQEで量子化学計算

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

追加のインストール

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

pip install blueqat openfermionblueqat openfermion

準備OKです。

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

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

#ツールの読み込み
from blueqat import *
from openfermion import *
from openfermionblueqat import*
import numpy as np

#水素分子の原子間の結合の距離を指定し、OpenFermionに入っているデータを読み込み、返します
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=[]

#0.2から2.5まで0.1刻みで原子間距離を入れる
for bond_len in np.arange(0.2,2.5,0.1):

    #分子の情報を得る
    m = get_molecule("{:.2}".format(bond_len))
    
    #ハミルトニアンを得て、第二量子化し、最後にbravi-kitaev変換でパウリ行列に変換
    h = bravyi_kitaev(get_fermion_operator(m.get_molecular_hamiltonian()))
    
    #あらかじめ用意したUCCAnsatzに入れて実行
    runner = vqe.Vqe(UCCAnsatz(h,6,Circuit().x[0]))
    result = runner.run()
    x.append(bond_len)
    e.append(runner.ansatz.get_energy_sparse(result.circuit))
    fullci.append(m.fci_energy)

#グラフを書く
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(x,fullci)
plt.plot(x,e,"o")
[<matplotlib.lines.Line2D at 0x7f07394e9dc0>]
<Figure size 432x288 with 1 Axes>output

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

2-3.中身確認

例題として、原子間距離が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で計算できるように変換をします。Bravi-Kitaev(ブラビキタエフ)変換やJordan–Wigner(ジョルダンウィグナー)変換がありますが、今回は前者のBK変換をしてみると、

h = bravyi_kitaev(get_fermion_operator(m.get_molecular_hamiltonian()))
print(h)
(0.7407724940116754+0j) [] +

(0.041128856023498556+0j) [X0 Z1 X2] +

(0.041128856023498556+0j) [X0 Z1 X2 Z3] +

(0.041128856023498556+0j) [Y0 Z1 Y2] +

(0.041128856023498556+0j) [Y0 Z1 Y2 Z3] +

(0.23528824284103544+0j) [Z0] +

(0.23528824284103542+0j) [Z0 Z1] +

(0.18133335898983727+0j) [Z0 Z1 Z2] +

(0.18133335898983727+0j) [Z0 Z1 Z2 Z3] +

(0.14020450296633868+0j) [Z0 Z2] +

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

runner = vqe.Vqe(UCCAnsatz(h,2,Circuit().x[0]))
result = runner.run(verbose = True)
params: [0.49026319 0.86349653] val: -0.16365997559518825

params: [0.49026319 0.86349653] val: -0.16365997559518825

params: [1.49026319 0.86349653] val: 0.6916397682896052

params: [-1.12777081  0.86349653] val: -0.09221605745375022

params: [0.49026319 0.86349653] val: -0.16365997559518825

params: [-0.12777079  0.86349653] val: -0.41575667374223013

params: [-0.12777079  0.86349653] val: -0.41575667374223013

params: [-0.12777079  1.86349653] val: 1.0948067779569997

params: [-0.12777079 -0.75453747] val: -0.3792627355192761

params: [-0.12777079  0.86349653] val: -0.41575667374223013

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

runner.ansatz.get_energy_sparse(result.circuit)
-0.9043560246890365

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

2-4.VQEで組合せ最適化問題

応用として組合せ最適化問題があります。多くの人に馴染みのある問題なので、量子コンピュータは正直どっから手をつけていいかわからんと思うひとおすすめです。組合せ最適化問題は、多くの選択肢からベストな答えを選ぶ問題で、社会問題をうまく定式化し、最小値問題を解くことで最適な組合せを得ることができます。

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

2-5.VQEで組合せ最適化問題の定式化

定式化は、答えはわからないが、条件はわかると言う状況に対して、最小値が出れば答えになるように、問題の条件を式の形に落としこみます。

・Z演算子を使う ・最終的にはZ演算子の代わりにQUBOを使う

です。まずは例題で、下記を行ってみましょう。定式化は、

from blueqat.pauli import X, Y, Z, I
h = -Z(0) - Z(0)*Z(1)

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

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

が設定されています。Zは期待値として<ψ|Z|ψ>=-1か<ψ|Z|ψ>1のどちらかをとります。hはより小さい値を取ると正解になります。Zの組合せで最終的な答えを判断できます。

VQEは最小となるZ(0)=1,Z(1)=1の時-2となるものを計算で求めてくれます。今回ansatzはa,bの2パラメータを利用した極座標表記のものを使ってみます。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)[1].cx[0,1]


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

print('Result by VQE')
print(runner.ansatz.get_energy_sparse(result.circuit))

# This is for check
mat = h.to_matrix()
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])
Result by VQE

-1.9753796902152434

Result by numpy

-2.0

約-2が答えとして出てきましたので正解です。上記問題は役にたたなさそうですが、

・Z(0)をAさん、Z(1)をBさんに見立てる ・Aさんはグループ１に属する（バイアスの設定） ・BさんはAさんと同じグループに属する（ウェイトの設定）

という条件をつけた(max-cut)分類問題と同じです。ただ、このままでは毎回問題を解くのが大変なので様々な工夫が必要です、それをみていきましょう。

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

定式化は+1と-1の値の組合せで定式化しました。ただ、産業では0と1を計算として利用します。そこで、01で定式化をするのをQUBOと呼びます。-1と+1でかかれた定式化を0と1で書かれたQUBO式に変換するには、Z演算子を下記のように変換するだけでできます。

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

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

QUBOで便利なルールとして、02=00^2=0そして12=11^2=1のように変形して字数を落とすことができます。

2-7.QUBO式のプログラミング

QUBOへの変換はZを代入して書き換えればいいのですが、毎回行うのは大変なのでblueqatで自動変換してくれます。今回のときたいコスト関数は、

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

にZ=2q−1Z = 2q-1を代入すると定式化は、

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

となりました。上記のコスト関数は以前のように-1と+1ではなく、0と1で考えることができるので簡単です。これをVQEで解くことで、最小値が得られます。

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, 2)
        self.step = 1

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

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

print('Result by VQE')
print(runner.ansatz.get_energy_sparse(result.circuit))

# This is for check
mat = h.to_matrix()
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])
Result by VQE

-1.999998044672769

Result by numpy

-2.0

vqe4.png

上記は定式化を変更しただけですので、量子回路自体は変更されていません。正しく得られていますが実際にはansatzを作るのが難しかったり、その他定式化に有利な回路を作るためには、VQEではなく、組合せ最適化問題用の回路を導入したQAOAを利用するのが無難です。

© 2024, blueqat Inc. All rights reserved