変分法と変分量子固有値ソルバー法を学習する

変分法と変分量子固有値ソルバー法を学習する#

この実習では、変分法の基本的な考え方と、その方法に基づいた変分量子アルゴリズムと呼ばれる量子計算の手法を学びます。特に、量子計算と古典計算を組み合わせた「量子・古典ハイブリッドアルゴリズム」としての変分量子アルゴリズムに着目します。この手法を用いて、近似的な固有値計算を可能にする変分量子固有値ソルバー法と呼ばれる方法へ拡張していきます。

はじめに #

行列で表現されるある物理系に対して、その最も小さい固有値を見つけるという操作は、多くのアプリケーションで必要となる重要な技術です。例えば化学の計算では、分子を特徴づけるエルミート行列の最小固有値はそのシステムの最もエネルギーの低い状態（基底状態）のエネルギーになります。最小固有値を見つけるには「量子位相推定」と呼ばれる手法（この課題を参照）を使うことができますが、この手法を使って実用的な問題を解こうとすると、そのために必要な量子回路はNISQコンピュータでは実現できないほど長くなることが知られています。そのために、短い量子回路を利用して分子の基底状態エネルギーを推定する手法として、変分量子固有値ソルバー法（Variational Quantum Eigensolver, VQE）が提案されました [PMS+14]。

まず、VQEの元になる関係を形式的に表現してみましょう。何か分からない最小固有値 $λ_{m i n}$ とその固有状態 $| ψ_{m i n} ⟩$ をもったエルミート行列 $H$ が与えられたとして、VQEはその系のエネルギーの下限である $λ_{m i n}$ の近似解 $λ_{θ}$ を求める手法です。つまり

λ_{m i n} \leq λ_{θ} \equiv ⟨ ψ (θ) | H | ψ (θ) ⟩

を満たす、できるだけ小さい $λ_{θ}$ を求めることに対応します。ここで $| ψ (θ) ⟩$ は近似解 $λ_{θ}$ に対応する状態で、 $θ$ はパラメータです。つまり、適当な初期状態 $| ψ ⟩$ にユニタリー $U (θ)$ で表現されるパラメータ化された回路を適用することで、 $| ψ_{m i n} ⟩$ を近似する状態 $| ψ (θ) ⟩ \equiv U (θ) | ψ ⟩$ を得ようというアイデアです。最適なパラメータ $θ$ の値は、期待値 $⟨ ψ (θ) | H | ψ (θ) ⟩$ が最小になるように古典計算を繰り返しながら求めていくことになります。

量子力学における変分法 #

背景 #

VQEは量子力学の変分法を応用した手法です。変分法をより良く理解するために、基礎的な数学的背景を見てみましょう。

行列 $A$ の固有ベクトル $| ψ_{i} ⟩$ とその固有値 $λ_{i}$ は、 $A | ψ_{i} ⟩ = λ_{i} | ψ_{i} ⟩$ という関係を持っていますね。行列 $H$ がエルミート行列 $H = H^{†}$ の場合、スペクトル定理から $H$ の固有値は実数になります（ $λ_{i} = λ_{i}^{*}$ ）。実際に実験で測定できる量は実数である必要があるため、量子系のハミルトニアンを記述するためにはエルミート行列が適切です。さらに、 $H$ は以下のように表現することもできます。

H = \sum_{i = 1}^{N} λ_{i} | ψ_{i} ⟩ ⟨ ψ_{i} |

ここで、各 $λ_{i}$ は固有ベクトル $| ψ_{i} ⟩$ に対応する固有値です。任意の量子状態 $| ψ ⟩$ に対して観測量 $H$ を測定した時の期待値は、以下の式で与えられます。

⟨ H ⟩_{ψ} \equiv ⟨ ψ | H | ψ ⟩

上式の $H$ を期待値の式に代入すると

\begin{array}{r} \begin{aligned} ⟨ H ⟩_{ψ} = ⟨ ψ | H | ψ ⟩ & = ⟨ ψ | (\sum_{i = 1}^{N} λ_{i} | ψ_{i} ⟩ ⟨ ψ_{i} |) | ψ ⟩ \\ = \sum_{i = 1}^{N} λ_{i} ⟨ ψ | ψ_{i} ⟩ ⟨ ψ_{i} | ψ ⟩ \\ = \sum_{i = 1}^{N} λ_{i} | ⟨ ψ_{i} | ψ ⟩ |^{2} \end{aligned} \end{array}

になります。最後の式は、任意の状態 $| ψ ⟩$ に対する $H$ の期待値は、 $λ_{i}$ を重みとした固有ベクトル $| ψ_{i} ⟩$ と $| ψ ⟩$ の内積（の絶対値二乗）の線形結合として与えられることを示しています。この式から、 $| ⟨ ψ_{i} | ψ ⟩ |^{2} \geq 0$ であるために

λ_{m i n} \leq ⟨ H ⟩_{ψ} = ⟨ ψ | H | ψ ⟩ = \sum_{i = 1}^{N} λ_{i} | ⟨ ψ_{i} | ψ ⟩ |^{2}

が成り立つことは明らかです。上記の式が変分法と呼ばれるもの（テキストによっては変分原理と呼ぶ）で、波動関数を「うまく取る」ことで、ハミルトニアン $H$ の期待値の下限として最小固有値を近似的に求めることができることを表しています。この式から、 $| ψ_{m i n} ⟩$ 状態の期待値は $⟨ ψ_{m i n} | H | ψ_{m i n} ⟩ = ⟨ ψ_{m i n} | λ_{m i n} | ψ_{m i n} ⟩ = λ_{m i n}$ になることも分かるでしょう。

基底状態の近似 #

系のハミルトニアンがエルミート行列 $H$ で表現されている場合、系の基底状態のエネルギーは $H$ の最小固有値になります。まず $| ψ_{m i n} ⟩$ の初期推定としてある波動関数 $| ψ ⟩$ （Ansatzと呼ばれる）を選び、その状態での期待値 $⟨ H ⟩_{ψ}$ を計算します。変分法の鍵は、この期待値が小さくなるように波動関数を更新しながら計算を繰り返し、ハミルトニアンの基底状態エネルギーに近づけていくところにあります。

変分量子アルゴリズム #

まず、VQEの元になっている変分量子アルゴリズム（Variational Quantum Algorithm, VQA）と呼ばれる手法について見ていきます。

変分量子回路 #

量子コンピューター上で変分法を実装するには、Ansatzを更新する仕組みが必要です。量子状態の更新には量子ゲートが使えることを、私たちは知っていますね。VQAも量子ゲートを使いますが、VQAは決まった構造を持つパラメータ化された量子回路（変分量子回路と呼ばれる）を使って行います。この回路は変分フォーム（variational form）と呼ばれる場合もあり、回路をひとまとめにしてユニタリー変換 $U (θ)$ と書くことも多いです（ $θ$ はパラメータで、複数ある場合はベクトルになります）。

変分フォームを初期状態 $| ψ ⟩$ （例えば標準状態 $| 0 ⟩$ ）に適用すると、出力として $| ψ (θ) ⟩ \equiv U (θ) | ψ ⟩$ が生成されます。この状態の元で期待値 $⟨ ψ (θ) | H | ψ (θ) ⟩$ が $λ_{m i n}$ に近付くように、 $| ψ (θ) ⟩$ に対してパラメータ $θ$ の最適化を行うのがVQEです。パラメータの最適化は古典計算で実行することを想定しており、その意味でVQEは典型的な量子・古典ハイブリッドアルゴリズムの一つです。

変分フォームの決め方ですが、解きたい問題のドメインに応じて特定の構造を持つ変分フォームを導入することがあります。そうではなく、幅広い問題への応用ができるようにドメインに依存しない形の変分フォーム（例えば $R_{X}$ や $R_{Y}$ などの回転ゲート）を使うこともあります。後で高エネルギー実験へのVQEの応用を課題として考えますが、そこでは $R_{Y}$ と制御 $Z$ ゲートを使った変分フォームを実装します。

単純な変分フォーム #

変分フォームを決める時には、2つの相反する目的に対するバランスを考える必要があります。 $n$ 量子ビットの変分フォームは、パラメータの数を増やせば実自由度 $2^{n + 1} - 2$ の任意の状態ベクトル $| ψ ⟩$ を生成できるでしょう。しかし、パラメータを最適化することを考えれば、できるだけ少ないパラメータで変分フォームを構築したいですよね。回転角をパラメータとする回転ゲートの数が増えれば、量子コンピュータで動かす場合はそれだけノイズの影響を受けやすくなります。なので、できるだけ少ないパラメータ（やゲート）で求める状態を生成できればベストでしょう。

ここでは、まず $n = 1$ の場合を考えます。Qiskitの $U$ ゲート（上の $U (θ)$ ノーテーションと紛らわしいですが、こちらは単一のゲート）は3つのパラメータ $θ$ 、 $ϕ$ 、 $λ$ を使って以下の変換を表現します:

\begin{array}{r} U (θ, ϕ, λ) = (\begin{array}{c} \cos \frac{θ}{2} & - e^{i λ} \sin \frac{θ}{2} \\ e^{i ϕ} \sin \frac{θ}{2} & e^{i λ + i ϕ} \cos \frac{θ}{2} \end{array}) \end{array}

変分フォームの初期状態を $| 0 ⟩$ に取るならば、上の行列の第一列のみが状態生成に寄与し、 $θ$ と $ϕ$ の2つのパラメータで任意の単一量子ビット状態を表現できます。そういう意味で、この変分フォームはユニバーサルであると言います。ただ、ユニバーサルに任意の状態を生成できるということは、この変分フォームが生成する状態を使ってあるハミルトニアン $H$ の期待値を計算しようとした場合、その固有状態を近似する状態だけでなく、それ以外のさまざまな状態も含むということになります。つまり、VQEで最小固有値を効率的に決められるかどうかは、こういう不要な状態を避けつつ、古典計算でいかに適切にパラメータを最適化できるかにかかっているわけです。

パラメータの最適化 #

パラメータ化された変分フォームを選択したら、所望の量子状態を作ることができるように、変分法に従ってパラメータを最適化する必要があります。パラメータの最適化のプロセスには様々な課題があります。例えば量子ハードウェアには様々なタイプのノイズがあるため、その状態でエネルギーを測定しても正しい答えが返ってくるという保証はありません。そのために、パラメータの最適化に使う目的関数の評価が実際のエネルギーの値からずれてしまい、正しくパラメータの更新ができない可能性があります。また、最適化の手法（オプティマイザー）によっては、パラメータの数に依って目的関数を評価する回数が増えることがあり、ノイズの影響をさらに受けやすくなります。つまり、アプリケーションの要求を考慮しながら、オプティマイザーの選択にも気を配る必要があります。

最も一般的な最適化手法は、エネルギーの減少が極大になるような方向に各パラメータを更新する勾配降下法です。各パラメータごとに勾配を計算するため、最適化すべきパラメータの数に応じて目的関数を評価する回数は増えます。また、この性質から探索空間の中で局所的な最適値を素早く発見することは可能ですが、逆に探索が局所的な最小値に留まってしまうことがあります。勾配降下法は直感的で理解しやすい最適化の手法ですが、少なくとも現在のNISQコンピュータでは精度良く実行することは難しいと考えられています（VQEを実装する時に、勾配計算を使った最適化の方法を紹介します）。

ノイズのある量子コンピュータで目的関数を最適化する適切なオプティマイザーとして、Simultaneous Perturbation Stochastic Approximation（SPSA）[BPP13]があります。SPSAは2回の測定だけで目的関数の勾配を近似できるという特徴があります。勾配降下法では各パラメータを独立に変化させるのに対して、SPSAでは全てのパラメータを同時にランダムに変化させます。以上のことから、VQEを利用する場合のオプティマイザーとしてSPSAが推奨されることがあります。

ノイズがない量子コンピュータで目的関数を評価する場合（例えば状態ベクトルシミュレータで実行する場合など）は、PythonのSciPyパッケージで提供されているオプティマイザーなど、様々な選択肢があります。この実習では、Qiskitでサポートされているオプティマイザーの中で、特にConstrained Optimization by Linear Approximation（COBYLA）と呼ばれるオプティマイザーも使用します。COBYLAは目的関数の評価を1回しか実行しない（つまり評価の回数がパラメータの数に依存しない）ため、ノイズがない状態でかつ評価の回数を少なくしたい場合にはCOBYLAの利用が推奨されているようです。いずれにしろ、どのオプティマイザーがベストかはVQEアルゴリズムの実装形式や実行環境によって変わるため、ある程度経験によって決める必要があると考えられます。

変分フォームを使った実例 #

ではここで、 $U$ ゲート一つからなる単一量子ビットの変分フォームを利用してパラメータ最適化の例を実行してみましょう。ランダムに量子状態 $| ψ_{0} ⟩$ を選び、 $| ψ (θ, ϕ) ⟩ := U (θ, ϕ, 0) | 0 ⟩$ でそれを近似するという問題を考えます[1]。1量子ビットの状態は観測量 $X, Y, Z$ の期待値 $⟨ X ⟩, ⟨ Y ⟩, ⟨ Z ⟩$ の値がわかれば（全体位相を除いて）完全に決まるので、 $| ψ (θ, ϕ) ⟩$ による $X, Y, Z$ の期待値 $⟨ X ⟩_{θ, ϕ}, ⟨ Y ⟩_{θ, ϕ}, ⟨ Z ⟩_{θ, ϕ}$ が $| ψ_{0} ⟩$ による対応する期待値 $⟨ X ⟩_{0}, ⟨ Y ⟩_{0}, ⟨ Z ⟩_{0}$ に等しくなるように $θ, ϕ$ を選ぶことになります。したがって、問題は目的関数

L (θ, ϕ) = [⟨ X ⟩_{θ, ϕ} - ⟨ X ⟩_{0}]^{2} + [⟨ Y ⟩_{θ, ϕ} - ⟨ Y ⟩_{0}]^{2} + [⟨ Z ⟩_{θ, ϕ} - ⟨ Z ⟩_{0}]^{2}

の最小化となります。

import numpy as np
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, transpile
from qiskit.circuit import Parameter, ParameterVector
from qiskit.primitives import BackendEstimator
from qiskit.quantum_info import Statevector, Operator, SparsePauliOp
from qiskit_algorithms.optimizers import SPSA, COBYLA
from qiskit_aer import AerSimulator

最初に、ターゲットとなる量子状態ベクトルをランダムに生成する関数と、状態ベクトルから $X, Y, Z$ の期待値を計算する関数を定義します。状態ベクトルはQiskitのStatevectorというクラスで表現し、パウリ演算子にはSparsePauliOpを使用します。

rng = np.random.default_rng(999999)

# 与えられた量子ビット数（nq）に応じたランダムな状態ベクトルを作る関数
def random_statevector(nq):
    # 2^nqの複素数をランダムに生成
    data = rng.random(2 ** nq) + 1.j * rng.random(2 ** nq)
    # 正規化
    data /= np.sqrt(np.sum(np.square(np.abs(data))))

    return Statevector(data)

# 例：U(π/3, π/6, 0)|0>
statevector = Statevector(np.array([np.cos(np.pi / 6.), np.exp(1.j * np.pi / 6.) * np.sin(np.pi / 6.)]))
for pauli in ['X', 'Y', 'Z']:
    op = SparsePauliOp(pauli)
    print(f'<{pauli}> = {statevector.expectation_value(op).real:.5f}')

次に、変分フォーム回路を定義します。このとき、 $U$ ゲートの回転角として、具体的な数値を設定せず、QiskitのParameterというオブジェクトを利用します。Parameterはあとから数値を代入できる名前のついた箱として使えます。

theta = Parameter('θ')
phi = Parameter('φ')

ansatz_1q = QuantumCircuit(1)
ansatz_1q.u(theta, phi, 0., 0)

# Parameterの値は未定
ansatz_1q.draw('mpl')

Parameterに値を代入するには、回路のassign_parametersメソッドを利用します。

# thetaとphiにπ/3とπ/6を代入
ansatz_1q.assign_parameters({theta: np.pi / 3., phi: np.pi / 6.}, inplace=False).draw('mpl')

変分フォーム回路が作る状態における $X, Y, Z$ の期待値を測定するための回路を定義します。

circuits = {}

# <X>を測るにはHゲートで基底を変換する
circuits['X'] = ansatz_1q.copy()
circuits['X'].h(0)
circuits['X'].measure_all()

# <Y>を測るにはSdg, Hゲートで基底を変換する
circuits['Y'] = ansatz_1q.copy()
circuits['Y'].sdg(0)
circuits['Y'].h(0)
circuits['Y'].measure_all()

# <Z>はそのままの回路で測れる
circuits['Z'] = ansatz_1q.copy()
circuits['Z'].measure_all()

それぞれの回路を通常通りバックエンドのrun()メソッドで実行し、結果から期待値を計算する関数を定義します。

backend = AerSimulator()

def circuit_expval(circuit, param_vals):
    bound_circuit = circuit.assign_parameters({theta: param_vals[0], phi: param_vals[1]}, inplace=False)

    bound_circuit_tr = transpile(bound_circuit, backend=backend)
    # shotsは関数の外で定義する
    job = backend.run(bound_circuit_tr, shots=shots)
    counts = job.result().get_counts()

    return (counts.get('0', 0) - counts.get('1', 0)) / shots

# 例：U(π/3, π/6, 0)|0>
shots = 10000
param_vals = [np.pi / 3., np.pi / 6.]
for pauli in ['X', 'Y', 'Z']:
    print(f'<{pauli}> = {circuit_expval(circuits[pauli], param_vals)}')

最小化する目的関数を定義します。

def objective_function(param_vals):
    loss = 0.
    for pauli in ['X', 'Y', 'Z']:
        # target_state_1qは関数の外で定義する
        op = SparsePauliOp(pauli)
        target = target_state_1q.expectation_value(op).real
        current = circuit_expval(circuits[pauli], param_vals)
        loss += (target - current) ** 2

    return loss

# 最適化の1ステップごとに呼び出される関数。目的関数の値をリストに記録しておく
def callback_function(param_vals):
    # lossesは関数の外で定義する
    losses.append(objective_function(param_vals))

最適化には使用しませんが、解を得たあとで変分フォームの終状態とターゲット状態とのフィデリティ $| ⟨ ψ_{0} | ψ (θ, ϕ) ⟩ |^{2}$ を計算する関数も定義しておきます。厳密に最適化が成功すれば、この関数の返り値は1になります。

def fidelity(ansatz, param_vals, target_state):
    # 量子回路のパラメータのリストはcircuit.parametersで取得できる
    parameters = ansatz.parameters

    param_binding = dict(zip(parameters, param_vals))
    opt_ansatz = ansatz.assign_parameters(param_binding, inplace=False)

    # Statevectorは回路からも生成可能（回路を|0>に対して適用した終状態になる）
    circuit_state = Statevector(opt_ansatz)

    return np.square(np.abs(target_state.inner(circuit_state)))

最後にCOBYLAオプティマイザーのインスタンスを作成し、アルゴリズムを実行します。

# COBYLAの最大ステップ数
maxiter = 300
# COBYLAの収束条件（小さいほどよい近似を目指す）
tol = 0.0001
# バックエンドでのショット数
shots = 1000

# オプティマイザーのインスタンス生成
optimizer = COBYLA(maxiter=maxiter, tol=tol, callback=callback_function)

# ターゲット状態
target_state_1q = random_statevector(1)

# thetaを[0, π), phiを[0, 2π)からランダムに選ぶ
init = [rng.uniform(0., np.pi), rng.uniform(0., 2. * np.pi)]

# 最適化を実行
losses = []
min_result = optimizer.minimize(objective_function, x0=init)

最適化プロセスにおけるロス（目的関数の返り値）の推移をプロットします。

plt.plot(losses);

fidelity(ansatz_1q, min_result.x, target_state_1q)

Estimatorの利用 #

VQEを含む変分量子アルゴリズムでは、上のように変分フォームにパラメータ値を代入し複数の観測量の期待値を計算するという手順の繰り返しが頻出します。そのため、これを自動化し、かつ（今は利用しませんが）様々なエラー補正なども適応してくれるEstimatorというクラスを使用することが推奨されています。特に、ここではBackendEstimatorという、特定のバックエンドを利用して計算をするタイプのEstimatorを利用します。

# BackendEstimatorインスタンスの生成
estimator = BackendEstimator(backend)

# 観測量はSparsePauliOpオブジェクトで表現
observables = [SparsePauliOp('X'), SparsePauliOp('Y'), SparsePauliOp('Z')]

param_vals = [np.pi / 3., np.pi / 6.]

# 変分フォーム、観測量、パラメータ値をrun()に渡す
# 観測量が3つあるので、ansatz_1qとparam_valuesも3つずつ
job = estimator.run([ansatz_1q] * 3, observables, [param_vals] * 3, shots=10000)
result = job.result()
result.values

Estimatorを使った目的関数を定義します。

observables_1q = [SparsePauliOp('X'), SparsePauliOp('Y'), SparsePauliOp('Z')]

def objective_function_estimator(param_vals):
    target = np.array(list(target_state_1q.expectation_value(op).real for op in observables_1q))

    job = estimator.run([ansatz_1q] * len(observables_1q), observables_1q, [param_vals] * len(observables_1q), shots=shots)
    current = np.array(job.result().values)

    return np.sum(np.square(target - current))

def callback_function_estimator(param_vals):
    # lossesは関数の外で定義する
    losses.append(objective_function_estimator(param_vals))

上の目的関数を最適化します。

# COBYLAの最大ステップ数
maxiter = 300
# COBYLAの収束条件（小さいほどよい近似を目指す）
tol = 0.0001
# バックエンドでのショット数
shots = 1000

# オプティマイザーのインスタンス生成
optimizer = COBYLA(maxiter=maxiter, tol=tol, callback=callback_function_estimator)

# ターゲット状態
target_state_1q = random_statevector(1)

# thetaを[0, π), phiを[0, 2π)からランダムに選ぶ
init = [rng.uniform(0., np.pi), rng.uniform(0., 2. * np.pi)]

# 最適化を実行
losses = []
min_result = optimizer.minimize(objective_function_estimator, x0=init)

fidelity(ansatz_1q, min_result.x, target_state_1q)

エンタングルメントの導入 #

では次に、この問題を2量子ビットに拡張してやってみましょう。2量子ビットの純粋状態は6個の実自由度を持ちますが、ここでは最も一般的に2量子ビット状態を決定する15個の観測量の期待値

⟨ O_{1} O_{2} ⟩ (O_{1}, O_{2} = I, X, Y, Z; O_{1} O_{2} \neq I I)

を測定します。ここで $I$ は恒等演算子です。

ターゲット状態に関する関数random_statevectorとpauli_expvalはそのまま利用できます。まず変分フォームとして2つの量子ビットに $U$ ゲートが一つずつかかっているものを考えて、最小化すべき目的関数を定義します。

# パラメータ数4なので、4要素のパラメータベクトルを作る
params = ParameterVector('params', 4)

ansatz_2q = QuantumCircuit(2)
ansatz_2q.u(params[0], params[1], 0., 0)
ansatz_2q.u(params[2], params[3], 0., 1)

paulis_1q = ['I', 'X', 'Y', 'Z']
paulis_2q = list(f'{op1}{op2}' for op1 in paulis_1q for op2 in paulis_1q if (op1, op2) != ('I', 'I'))
observables_2q = list(SparsePauliOp(pauli) for pauli in paulis_2q)

def objective_function_2q(param_vals):
    # target_state_2qは関数の外で定義
    target = np.array(list(target_state_2q.expectation_value(op).real for op in observables_2q))

    job = estimator.run([ansatz_2q] * len(observables_2q), observables_2q, [param_vals] * len(observables_2q), shots=shots)
    current = np.array(job.result().values)

    return np.sum(np.square(target - current))

def callback_function_2q(param_vals):
    # lossesは関数の外で定義する
    losses.append(objective_function_2q(param_vals))

# COBYLAの最大ステップ数
maxiter = 300
# COBYLAの収束条件（小さいほどよい近似を目指す）
tol = 0.0001
# バックエンドでのショット数
shots = 1000

# オプティマイザーのインスタンス生成
optimizer = COBYLA(maxiter=maxiter, tol=tol, callback=callback_function_2q)

# ターゲット状態
target_state_2q = random_statevector(2)

# パラメータの初期値
init = rng.uniform(0., 2. * np.pi, size=4)

# 最適化を実行
losses = []
min_result = optimizer.minimize(objective_function_2q, x0=init)

fidelity(ansatz_2q, min_result.x, target_state_2q)

やってみると分かりますが、結果は1量子ビットの場合と比べて良くないですね。どうすれば良くなるでしょうか？（やり方は複数あると思います）

一つの解決策：変分フォームにエンタングルメントを導入する

ansatz_2q = QuantumCircuit(2)
ansatz_2q.u(params[0], params[1], 0., 0)
ansatz_2q.u(params[2], params[3], 0., 1)
ansatz_2q.cx(0, 1)

どうなるか確かめてください。

2量子ビットの一般の状態では2つのビットがエンタングルしているので、変分フォームに2量子ビットゲートを入れると近似精度が良くなるのはある意味当然です。例えば、ベル状態（CHSH不等式の破れを確認するを参照）を再現したいときにこの状況をクリアに見ることができます。上で

target_state_2q = random_statevector(2)

を

target_state_2q = Statevector(np.array([1., 0., 0., 1.], dtype=complex) / np.sqrt(2.))

として実行するとどうなるでしょうか。エンタングルさせる場合とさせない場合で大きな違いが見えるでしょう。

問題を3量子ビットに拡張し、GHZ状態（単純な量子回路をゼロから書くを参照）

target_state_3q = Statevector(np.array([1.] + [0.] * 6 + [1.], dtype=complex) / np.sqrt(2.))

をターゲットにするなどして、遊んでみてください。

変分量子固有値ソルバー法 #

それでは、ここからVQEの簡単な実装例を見ていくことにします。

パラメータシフト法 #

VQEを実装する前に、目的関数の勾配を使った最適化の方法（パラメータシフト法と呼ばれる）について説明します。あるタイプの変分量子回路では、目的関数の勾配を厳密に計算できることが知られています。その勾配の値を使って、古典計算でパラメータを最適化することができます。

パラメータシフト法を使って勾配を導出するために、まず一般的なパラメータ付きユニタリー演算 $U (θ) = \prod_{j = 1}^{L} U_{j} (θ_{j})$ を考えます。 $U_{j} (θ_{j})$ はパラメータ $θ_{j}$ を持つユニタリーで、例えば $θ_{j}$ を回転角とする回転ゲート $U_{j} (θ_{j}) = \exp (- i θ_{j} P_{j} / 2)$ を考えてみます（ $P_{j} \in {X, Y, Z}$ はパウリ演算子です）。初期状態 $ρ$ を $U (θ)$ で時間発展させた状態に対して、物理量 $M$ を測定したときの期待値 $⟨ M (θ) ⟩$ は

⟨ M (θ) ⟩ = Tr [M U (θ) ρ U (θ)^{†}] = Tr [M U_{L : 1} ρ U_{L : 1}^{†}]

になります（ここで $U_{l : m} := \prod_{j = m}^{l} U_{j} (θ_{j})$ と略記することにします）。この期待値のパラメータ $θ_{j}$ に対する勾配は

\frac{\partial}{\partial θ_{j}} ⟨ M (θ) ⟩ = Tr [M \frac{\partial U_{L : 1}}{\partial θ_{j}} ρ U_{L : 1}^{†}] + Tr [M U_{L : 1} ρ \frac{\partial U_{L : 1}^{†}}{\partial θ_{j}}]

です。 $P_{j}^{†} = P_{j}$ なので

\begin{array}{r} \begin{aligned} \frac{\partial U_{L : 1}}{\partial θ_{j}} & = U_{L} \dots U_{j + 1} \frac{\partial U_{j}}{\partial θ_{j}} U_{j - 1} \dots U_{1} = - \frac{i}{2} U_{L : j} P_{j} U_{j - 1 : 1} \\ \frac{\partial U_{L : 1}^{†}}{\partial θ_{j}} & = \frac{i}{2} U_{j - 1 : 1}^{†} P_{j} U_{L : j}^{†} \end{aligned} \end{array}

という関係が成り立ちますが、この関係から

\frac{\partial}{\partial θ_{j}} ⟨ M (θ) ⟩ = - \frac{i}{2} Tr [M U_{L : j} [P_{j}, U_{j - 1 : 1} ρ U_{j - 1 : 1}^{†}] U_{L : j}^{†}]

が得られます。ここで $U_{j} (θ_{j})$ の定義に戻ると、 $U_{j} (θ_{j}) = \exp (- i θ_{j} P_{j} / 2) = \cos (θ_{j} / 2) I - i \sin (θ_{j} / 2) P_{j}$ （ $I$ は恒等演算子）なので、 $U (\pm π / 2) = (1 / \sqrt{2}) (I \mp i P_{j})$ になります。従って、

U_{j} (\frac{π}{2}) ρ U_{j}^{†} (\frac{π}{2}) - U_{j} (- \frac{π}{2}) ρ U_{j}^{†} (- \frac{π}{2}) = \frac{1}{2} (I - i P_{j}) ρ (I + i P_{j}^{†}) - \frac{1}{2} (I + i P_{j}) ρ (I - i P_{j}^{†}) = - i [P_{j}, ρ]

と書くことができます。この結果を上の $\partial ⟨ M (θ) ⟩ / \partial θ_{j}$ の式に適用すると、

\begin{array}{r} \begin{aligned} \frac{\partial}{\partial θ_{j}} ⟨ M (θ) ⟩ & = - \frac{i}{2} Tr [M U_{L : j} [P_{j}, U_{j - 1 : 1} ρ U_{j - 1 : 1}^{†}] U_{L : j}^{†}] \\ = \frac{1}{2} Tr [M U_{L : j + 1} U_{j} (θ_{j} + \frac{π}{2}) U_{j - 1 : 1} ρ U_{j - 1 : 1}^{†} U_{j}^{†} (θ_{j} + \frac{π}{2}) U_{L : j + 1}^{†} - M U_{L : j + 1} U_{j} (θ_{j} - \frac{π}{2}) U_{j - 1 : 1} ρ U_{j - 1 : 1}^{†} U_{j}^{†} (θ_{j} - \frac{π}{2}) U_{L : j + 1}^{†})] \\ = \frac{1}{2} [⟨ M (θ + \frac{π}{2} e_{j}) ⟩ - ⟨ M (θ - \frac{π}{2} e_{j}) ⟩] \end{aligned} \end{array}

が得られます。ここで $e_{j}$ は $j$ 番目の要素だけが1でそれ以外は0のベクトルです。

この式から、期待値 $⟨ M (θ) ⟩$ のパラメータ $θ_{j}$ に対する勾配は、 $θ_{j}$ を $\pm π / 2$ だけ増減させて求めた期待値の差として求めることができることが分かりました。これがパラメータシフト法と呼ばれるものです。

VQEの実装 #

では、パラメータシフト法を使って簡単なVQEの例を実装してみます。ある観測量の期待値が最小になるように、VQEを使ってAnsatzを更新する回路パラメータを決定するという問題を考えてみます。

量子回路として、 $R_{Y} R_{Z}$ ゲートを繰り返すシンプルなパラメータ回路を使い、観測量として、パウリ演算子のテンソル積 $Z X Y$ を使います。

パラメータシフト法の実装は、QiskitではParamShiftEstimatorGradientというAPIを使うことで一行で済んでしまいます（実際の勾配計算に興味がある人は、期待値から勾配を直接求めるコードを書いて、パラメータ毎に $\pm π / 2$ シフトさせた回路を走らせることで、このAPIと同じ勾配が得られることを確認してみて下さい）。パラメータの最適化は、勾配を使って勾配降下を行うConjugate Descent (CG)とGradient Descentの2つのオプティマイザーを使って行いますが、比較のためにCOBYLAも使うことにします。

最終的に、3通りのVQEを使って求めた最小エネルギーの近似解を、厳密対角化して求めた最小エネルギーの値と比較することにします。

from qiskit_algorithms.minimum_eigensolvers import VQE, NumPyMinimumEigensolver
from qiskit_algorithms.optimizers import CG, GradientDescent
from qiskit_algorithms.gradients import ParamShiftEstimatorGradient

# Ansatzの定義
num_qubits = 3   # 量子ビット数
num_layers = 2  # レイヤー数

ansatz = QuantumCircuit(num_qubits)

# 長さ0のパラメータ配列
theta = ParameterVector('θ')

# 配列に要素を一つ足して、新しく足されたパラメータを返す
def new_theta():
    theta.resize(len(theta) + 1)
    return theta[-1]

for _ in range(num_layers):
    for iq in range(num_qubits):
        ansatz.ry(new_theta(), iq)

    for iq in range(num_qubits):
        ansatz.rz(new_theta(), iq)

    #for iq in range(num_qubits - 1):
    #    ansatz.cx(iq, iq + 1)

ansatz.draw('mpl')

# 最小固有値を求める観測量
obs = SparsePauliOp('ZXY')

# パラメータの初期値
init = rng.uniform(0., 2. * np.pi, size=len(theta))

# Estimatorを使って観測量の勾配を計算するオブジェクト
grad = ParamShiftEstimatorGradient(estimator)

# Conjugate gradientを使ったVQE
optimizer_cg = CG(maxiter=200)
vqe_cg = VQE(estimator, ansatz, optimizer_cg, gradient=grad, initial_point=init)

# Gradient descentを使ったVQE
optimizer_gd = GradientDescent(maxiter=200)
vqe_gd = VQE(estimator, ansatz, optimizer_gd, gradient=grad, initial_point=init)

# COBYLAを使ったVQE
optimizer_cobyla = COBYLA(maxiter=300)
vqe_cobyla = VQE(estimator, ansatz, optimizer_cobyla, initial_point=init)

# 厳密解を計算するソルバー
ee = NumPyMinimumEigensolver()

result_vqe_cg = vqe_cg.compute_minimum_eigenvalue(obs)
result_vqe_gd = vqe_gd.compute_minimum_eigenvalue(obs)
result_vqe_cobyla = vqe_cobyla.compute_minimum_eigenvalue(obs)
result_ee = ee.compute_minimum_eigenvalue(obs)

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[34], line 5
      3 if os.getenv('JUPYTERBOOK_BUILD') == '1':
      4     with open('data/vqe_results.pkl', 'rb') as source:
----> 5         result_ee, result_vqe_cobyla, result_vqe_cg, result_vqe_gd = pickle.load(source)

File symengine_wrapper.pyx:831, in symengine.lib.symengine_wrapper.load_basic()

RuntimeError: SymEngine-0.13 was asked to deserialize an object created using SymEngine-0.11.

print('Result:')
print(f'  Exact      = {result_ee.eigenvalue}')
print(f'  VQE(COBYLA) = {result_vqe_cobyla.optimal_value}')
print(f'  VQE(CG)    = {result_vqe_cg.optimal_value}')
print(f'  VQE(GD)    = {result_vqe_gd.optimal_value}')

COBYLAを使ってVQEで求めた答えは、厳密解(=-1.0)に近くなっていると思います。勾配計算を使ったVQEも多くの場合正しい答えに近い答えを返しますが、パラメータの初期値によってかなり悪い結果を返してしまう場合があります。

Ansatzの回路構成を変えたり、観測量を変えてみるなどして色々試してみてください。

変分法と変分量子固有値ソルバー法を学習する

Contents

変分法と変分量子固有値ソルバー法を学習する#