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

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

この教材は、Qiskit textbookの”Simulating Molecules using VQE”を参考にしています。

\(\newcommand{\ket}[1]{| #1 \rangle}\) \(\newcommand{\bra}[1]{\langle #1 |}\) \(\newcommand{\braket}[2]{\langle #1 | #2 \rangle}\) \(\newcommand{\expval}[3]{\langle #1 | #2 | #3 \rangle}\)

はじめに

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

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

\[ \lambda_{min} \le \lambda_{\theta} \equiv \expval{ \psi(\theta)}{H}{\psi(\theta) } \]

を満たす、できるだけ小さい\(\lambda_{\theta}\)を求めることに対応します。ここで\(\ket{\psi(\theta)}\)は近似解\(\lambda_{\theta}\)に対応する固有状態で、\(\theta\)はパラメータです。つまり、適当な初期状態\(\ket{\psi}\)にユニタリー\(U(\theta)\)で表現されるパラメータ化された回路を適用することで、\(\ket{\psi_{min}}\)を近似する状態\(\ket{\psi(\theta)} \equiv U(\theta)\ket{\psi}\)を得ようというアイデアです。最適なパラメータ\(\theta\)の値は、期待値 \(\expval{\psi(\theta)}{H}{\psi(\theta)}\)が最小になるように古典計算を繰り返しながら求めていくことになります。

量子力学における変分法

背景

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

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

\[ H = \sum_{i = 1}^{N} \lambda_i \ket{\psi_i} \bra{ \psi_i } \]

ここで、各\(\lambda_i\)は固有ベクトル\(\ket{\psi_i}\)に対応する固有値です。任意の量子状態\(\ket{\psi}\)に対して観測量\(H\)を測定した時の期待値は、以下の式で与えられます。

\[ \langle H \rangle_{\psi} \equiv \expval{ \psi }{ H }{ \psi } \]

上式の\(H\)を期待値の式に代入すると

\[\begin{split} \begin{aligned} \langle H \rangle_{\psi} = \expval{ \psi }{ H }{ \psi } &= \bra{ \psi } \left(\sum_{i = 1}^{N} \lambda_i \ket{\psi_i} \bra{ \psi_i }\right) \ket{\psi}\\ &= \sum_{i = 1}^{N} \lambda_i \braket{ \psi }{ \psi_i} \braket{ \psi_i }{ \psi} \\ &= \sum_{i = 1}^{N} \lambda_i | \braket{ \psi_i }{ \psi} |^2 \end{aligned} \end{split}\]

になります。最後の式は、任意の状態\(\ket{\psi}\)に対する\(H\)の期待値は、\(\lambda_i\)を重みとした固有ベクトル\(\ket{\psi_i}\)\(\ket{\psi}\)の内積(の絶対値二乗)の線形結合として与えられることを示しています。この式から、\(| \braket{ \psi_i }{ \psi} |^2 \ge 0\) であるために

\[ \lambda_{min} \le \langle H \rangle_{\psi} = \expval{ \psi }{ H }{ \psi } = \sum_{i = 1}^{N} \lambda_i | \braket{ \psi_i }{ \psi} |^2 \]

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

基底状態の近似

系のハミルトニアンがエルミート行列\(H\)で表現されている場合、系の基底状態のエネルギーは\(H\)の最小固有値になります。まず\(\ket{\psi_{min}}\)の初期推定として任意の波動関数\(\ket{\psi }\)Ansatzと呼ばれる)を選び、その状態での期待値\(\langle H \rangle_{\psi}\)を計算します。変分法の鍵は、この期待値が小さくなるように波動関数を更新しながら計算を繰り返し、ハミルトニアンの基底状態エネルギーに近づけていくところにあります。

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

変分量子回路

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

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

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

単純な変分フォーム

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

ここでは、まず\(n=1\)の場合を考えます。\(U3\)ゲートは3つのパラメータ\(\theta\)\(\phi\)\(\lambda\)を使って以下の変換を表現します:

\[\begin{split} U3(\theta, \phi, \lambda) = \begin{pmatrix}\cos\frac{\theta}{2} & -e^{i\lambda}\sin\frac{\theta}{2} \\ e^{i\phi}\sin\frac{\theta}{2} & e^{i\lambda + i\phi}\cos\frac{\theta}{2} \end{pmatrix} \end{split}\]

系全体にかかるグローバルな位相を除けば、3つのパラメータを適切に設定して実装すれば任意の単一量子ビット状態への変換が行えます。そういう意味で、この変分フォームはユニバーサルな変換が可能で、かつ3つしかパラメータがないため効率的に最適化できるという特徴があります。ただユニバーサルに任意の状態を生成できるということは、この変分フォームが生成する状態を使ってあるハミルトニアン\(H\)の期待値を計算しようとした場合、その固有状態を近似する状態だけでなく、それ以外のさまざまな状態も含むということになります。つまり、こういう不要な状態もたくさん作るということは、VQEで最小固有値を効率的に決めるためには、古典計算でいかにうまくパラメータを最適化できるかが重要ということでもあるわけです。

パラメータの最適化

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

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

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

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

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

ではここで、単一量子ビットの変分フォームを利用してパラメータ最適化の例を実行してみましょう。例として、ランダムな確率分布のベクトル\(\vec{x}\)(要素数は2)を入力として与えた時、出力の確率分布が\(\vec{x}\)に近くなるように単一量子ビットの変分フォームを決定するという問題を考えます(2つの確率分布の近さはL1距離によって定義します)。

vqe_u3

最初に、Pythonでランダムな確率分布のベクトルを作成します。

# Tested with python 3.7.9, qiskit 0.23.5, numpy 1.20.1
import numpy as np
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, Aer, execute
from qiskit.aqua.components.optimizers import COBYLA
nq = 1  # 量子ビットの数

npar = 3*nq  # パラメータの数

np.random.seed(999999)
target_distr = np.random.rand(2**nq)
target_distr /= sum(target_distr)

次に、単一の\(U3\)変分フォームの3つのパラメータを引数として受け取り、対応する量子回路を返す関数を定義します。

def get_var_form(params):
    qr = QuantumRegister(nq, name="q")
    cr = ClassicalRegister(nq, name='c')
    qc = QuantumCircuit(qr, cr)

    for i in range(nq):
        qc.u(params[3*i], params[3*i+1], params[3*i+2], qr[i])

    for i in range(nq):
        qc.measure(qr[i], cr[i])
    return qc

get_var_form(np.random.rand(npar)).draw('mpl')

変分フォームのパラメータのリストを入力とし、パラメータに対応したコストを計算する目的関数を定義します。アルゴリズムを実行するバックエンドとして、QASMシミュレータを使用します。

backend = Aer.get_backend("qasm_simulator")
NUM_SHOTS = 10000  # 測定する回数

# 出力されるビット列の確率分布を計算
def get_probability_distribution(counts):
    output_distr = []
    for i in range(2**nq):
        match = False
        for (k,v) in counts.items():
            if i == int(k,2):
                output_distr.append(v/NUM_SHOTS)
                match = True
        if not match:
            output_distr.append(0)

    if len(output_distr) == 1:
        output_distr.append(0)
    return output_distr

# コストを計算する目的関数を定義
def objective_function(params):
    qc = get_var_form(params)
    result = execute(qc, backend, shots=NUM_SHOTS).result()
    output_distr = get_probability_distribution(result.get_counts(qc))
    cost = sum([np.abs(output_distr[i] - target_distr[i]) for i in range(2**nq])])
    return cost

最後にCOBYLAオプティマイザーのインスタンスを作成し、アルゴリズムを実行します。出力される確率分布は実行の度に異なり、ターゲットの確率分布と完全には同じにならないことに注意してください。出力の精度は量子計算の回数(ショット数=NUM_SHOTS)に依存するので、ショット数を増減させた時の一致具合を確認してみてください。

optimizer = COBYLA(maxiter=500, tol=0.0001)

params = np.random.rand(npar)
ret = optimizer.optimize(num_vars=npar, objective_function=objective_function, initial_point=params)

qc = get_var_form(ret[0])
counts = execute(qc, backend, shots=NUM_SHOTS).result().get_counts(qc)
output_distr = get_probability_distribution(counts)

print("Target Distribution:         ", np.round(target_distr,4))
print("Obtained Distribution:       ", np.round(np.array(output_distr),4))
print("Cost Value (L1-Distance):     {:.6f}".format(ret[1]))
print("Parameters Found:            ", np.round(ret[0],4))

では次に、この問題を2量子ビット(確率分布の要素数は4)に拡張してやってみましょう。上に戻って

nq = 2  # 量子ビットの数

として再度実行するとどういう結果が得られるでしょうか。量子回路とオプティマイザーの関係はこのようになってますね。

vqe_2q_u3

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

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

    for i in range(nq):
        qc.u(params[3*i], params[3*i+1], params[3*i+2], qr[i])
    qc.cx(qr[0],qr[1])

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

量子ビットをエンタングルさせることで相関のあるデータを表現しやすくなるという状況は、例えば、ベル状態(CHSH不等式の破れを確認するを参照)の確率分布を再現したいときにクリアに見ることができます。上で

target_distr = np.random.rand(2**nq)

# 00と11を測定する確率が50%、01と10の確率は0
target_distr = np.array([0.5,0.,0.,0.5])

として実行するとどうなるでしょうか。エンタングルさせる場合とさせない場合で大きな違いが見えるでしょう。3量子ビットのGHZ状態(単純な量子回路をゼロから書くを参照)

# 000と111を測定する確率が50%、それ以外の確率は0
target_distr = np.array([0.5,0.,0.,0.,0.,0.,0.,0.5])

に拡張してみるなどして、遊んでみてください。

参考文献

BPP13

S Bhatnagar, H L Prasad, and L A Prashanth. Stochastic Recursive Algorithms for Optimization. Springer, 2013. doi:10.1007/978-1-4471-4285-0.

PMS+14

Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O'Brien. A variational eigenvalue solver on a photonic quantum processor. Nature Communications, 5(1):4213, 2014. URL: https://doi.org/10.1038/ncomms5213, doi:10.1038/ncomms5213.