変分量子固有値ソルバー（VQE）

このレッスンでは、変分量子固有値ソルバーを紹介し、量子コンピュータの基礎アルゴリズムとしての重要性を説明し、その長所と短所を探る。 VQEだけでは、補強的な手法を用いなければ、現代の実用的な規模の量子計算には不十分であろう。とはいえ、古典と量子のハイブリッド法の原型として重要であり、より高度な多くのアルゴリズムが構築される重要な基礎となっている。

このビデオでは、VQEの概要と、その効率に影響を与える要因について説明する。以下の文章では、さらに詳細を補足し、Qiskitを使用してVQEを実装している。

1. VQEとは何ですか？

変分量子固有値ソルバーとは、古典計算と量子計算を組み合わせてタスクを実行するアルゴリズムである。 VQEの算出には、主に4つの要素があります：

演算子 ：多くの場合、ハミルトニアン（ここでは $H$ と呼ぶ）で、最適化したいシステムの特性を記述する。別の言い方をすれば、最小固有値に対応するこの作用素の固有ベクトルを求めていることになる。その固有ベクトルを「基底状態」と呼ぶことが多い。
ansatz」 （ドイツ語で「アプローチ」を意味する）：これは、求める固有ベクトルに近似した量子状態を準備する量子回路である。というのも、いくつかのゲートはパラメータ化されており、パラメータを変化させることができるからだ。この量子回路のファミリーは、基底状態に近似した量子状態のファミリーを準備することができる。
推定量 ：現在の変分量子状態における演算子 $H$ の期待値を推定する手段。時に、私たちが本当に重視しているのは、単にこの期待値であり、これをコスト関数と呼ぶ。場合によっては、1つ以上の期待値から導き出すことができる、より複雑な関数に関心を持つこともあります。
古典的なオプティマイザ ：パラメータを変化させてコスト関数を最小化しようとするアルゴリズム。

それぞれの構成要素をさらに詳しく見てみよう。

1.1 演算子（ハミルトニアン）

VQE問題の核心は、関心のあるシステムを記述する演算子である。ここでは、この演算子の最小固有値と対応する固有ベクトルが、何らかの科学的あるいはビジネス上の目的に有用であると仮定する。例えば、分子を記述する化学的ハミルトニアンでは、作用素の最小固有値が分子の基底状態エネルギーに対応し、対応する固有状態が分子の形状や電子配置を記述する。あるいは、オペレータは最適化されるべき特定のプロセスのコストを記述することができ、固有状態はルートまたはプラクティスに対応することができる。物理学のようないくつかの分野では、「ハミルトニアン」はほとんどの場合、物理系のエネルギーを記述する作用素を指す。しかし、量子コンピューティングでは、ビジネスやロジスティクスの問題を記述する量子演算子も「ハミルトニアン」と呼ばれるのが一般的だ。ここではその慣例を採用する。

原子軌道を示す図と、多数のノードとそれらを結ぶ接続からなるネットワークを示す図。

物理的な問題や最適化問題を量子ビットにマッピングすることは、一般的に自明なことではないが、その詳細はこのコースの焦点ではない。量子演算子への問題のマッピングに関する一般的な議論は、 Quantum Computing in Practiceに掲載されている。化学の問題を量子作用素にマッピングする方法については、『 Quantum Chemistry with VQE 』に詳しい。

このコースでは、ハミルトニアンの形が既知であると仮定する。例えば、単純な水素分子のハミルトニアン（ある活性空間の仮定の下で、ヨルダン・ウィグナー・マッパーを使用）は次のようになる：

from qiskit.quantum_info import SparsePauliOp

hamiltonian = SparsePauliOp(
    [
        "IIII",
        "IIIZ",
        "IZII",
        "IIZI",
        "ZIII",
        "IZIZ",
        "IIZZ",
        "ZIIZ",
        "IZZI",
        "ZZII",
        "ZIZI",
        "YYYY",
        "XXYY",
        "YYXX",
        "XXXX",
    ],
    coeffs=[
        -0.09820182 + 0.0j,
        -0.1740751 + 0.0j,
        -0.1740751 + 0.0j,
        0.2242933 + 0.0j,
        0.2242933 + 0.0j,
        0.16891402 + 0.0j,
        0.1210099 + 0.0j,
        0.16631441 + 0.0j,
        0.16631441 + 0.0j,
        0.1210099 + 0.0j,
        0.17504456 + 0.0j,
        0.04530451 + 0.0j,
        0.04530451 + 0.0j,
        0.04530451 + 0.0j,
        0.04530451 + 0.0j,
    ],
)

上記のハミルトニアンには、 ZZII や YYYY のように、互いに通約しない項があることに注意されたい。つまり、 ZZII を評価するためには、（他の測定の中でも）量子ビット3のパウリZ演算子を測定する必要がある。しかし、 YYYY を評価するためには、同じ量子ビットである量子ビット3上のパウリY演算子を測定する必要がある。同じ量子ビット上のY演算子とZ演算子の間には不確定性関係があり、その両方を同時に測定することはできない。この点については後述する。上記のハミルトニアンは $16\times 16$ 行列演算子である。この演算子を対角化して最小エネルギー固有値を求めることは難しくない。

import numpy as np

A = np.array(hamiltonian)
eigenvalues, eigenvectors = np.linalg.eigh(A)
print("The ground state energy is ", min(eigenvalues), "hartrees")

Output:

The ground state energy is  -1.1459778447627311 hartrees

古典的な固有値解法では、薬やタンパク質のような非常に大きな原子系のエネルギーや形状を記述することはできない。 VQEは、この問題に量子コンピューティングを活用する初期の試みのひとつである。

このレッスンでは、上記のハミルトニアンよりもはるかに大きなハミルトニアンに出会うことになる。しかし、このコースの後半で、VQEを補強したり置き換えたりできる、より高度なツールを紹介する前に、VQEでできることの限界に挑戦するのは無駄なことだろう。

1.2 アンツァッツ

ansatz」はドイツ語で「アプローチ」を意味する。ドイツ語の正しい複数形は "ansätze "だが、"ansatzes "や "ansatze "をよく見かける。 VQEの文脈では、アンサッツとは、研究している系の基底状態に最も近い多量子ビットの波動関数を作成するために使用する量子回路のことで、その結果、演算子の期待値が最も小さくなります。この量子回路には変分パラメータ（多くの場合、変数ベクトル $\vec{\theta}$ ）が含まれる。

変分パラメータの初期値セット $\vec{\theta_0}$ 。回路上のansatzのユニタリー演算を $U_{\text{var}}(\vec{\theta_0})$ と呼ぶことにする。デフォルトでは、 IBM® 量子コンピュータのすべての量子ビットは、 $|0\rangle$ の状態に初期化されている。回路が実行されると、量子ビットの状態は次のようになる

|\psi(\vec{\theta_0})\rangle=U_{\text{var}}(\vec{\theta_0})|0\rangle^{\otimes N}

最も低いエネルギー（物理系の言葉を使う）さえあればいいのであれば、エネルギーを何度も測定し、最も低いものを選ぶだけで推定できる。しかし、私たちは通常、最も低いエネルギーや固有値をもたらすコンフィギュレーションも求める。そこで次のステップは、ハミルトニアンの期待値の推定であり、これは量子測定によって達成される。そのためには多くのことが必要だ。しかし、あるエネルギー（ $E_j$ ）を測定する確率（ $P_j$ ）が、期待値と次のように関連していることに注目すれば、このプロセスを定性的に理解することができる：

\langle \psi(\vec{\theta_0}) |H|\psi (\vec{\theta_0}) \rangle

確率 $P_j$ は、固有状態 $|\phi_j\rangle$ とシステムの現在の状態 $|\psi(\vec{\theta_0})\rangle$ との重なりにも関係する：

P_j=|\langle \phi_j|\psi(\vec{\theta_0})\rangle|^2 = |\langle \phi_j|U_{\text{var}}(\vec{\theta_0})|0\rangle^{\otimes N}|^2

そこで、ハミルトニアンを構成するパウリ作用素を何度も測定することで、系の現在の状態（ $|\psi(\vec{\theta_0})\rangle$ ）におけるハミルトニアンの期待値を推定することができる。次のステップは、パラメータ $\vec{\theta}$ を変化させ、系の最低エネルギー（基底）状態に近づけようとすることである。アナザッツには変分パラメータがあるため、変分形式と呼ばれることもある。

変分プロセスに移る前に、"良い推測 "の状態から状態をスタートさせることがしばしば有効であることに注意しよう。 $|0\rangle^{\otimes N}$ 例えば、化学的なアプリケーションでは、量子ビットをハートリーフォック状態に初期化するのが一般的です。変分パラメータを含まないこの開始推測を参照状態と呼ぶ。参照状態を生成するための量子回路を $U_{ref}$ と呼ぶことにする。参照状態をアサッツの残りの部分と区別することが重要になる場合は、いつでも使用します： $U_{\text{ansatz}}(\vec{\theta}) =U_{\text{var}}(\vec{\theta})U_{\text{ref}}.$ 等価的に

\begin{aligned} |\psi_{\text{ref}}\rangle&=U_{\text{ref}}|0\rangle^{\otimes N}\\ |\psi_{\text{ansatz}}(\vec{\theta})\rangle&=U_{var}(\vec{\theta})|\psi_{\text{ref}}\rangle = U_{\text{var}}(\vec{\theta})U_{\text{ref}}|0\rangle^{\otimes N}. \end{aligned}

1.3 見積もり担当者

特定の変分状態 $|\psi(\vec{\theta})\rangle$ におけるハミルトニアンの期待値を推定する方法が必要である。もし作用素全体 $H$ を直接測定できれば、多くの（例えば $N$ ）測定を行い、測定値を平均化するだけでよい：

\langle \psi(\vec{\theta})|H|\psi(\vec{\theta})\rangle _N \approx \frac{1}{N}\sum_{j=1}^N {E_j}

ここで、 $\approx$ の記号は、この期待値が $N\rightarrow \infty$ のような極限においてのみ正確に正しいことを思い出させる。しかし、回路上で何千もの測定が行われているため、期待値のサンプリング誤差はかなり小さい。また、非常に精密な計算をする場合には、ノイズが問題になる。

しかし、一般的に、 $H$ を一度にすべて測定することは不可能です。 $H$ には、互いに可換でない複数のパウリ演算子X、Y、Zが含まれている可能性があります。したがって、ハミルトニアンは、同時に測定可能な演算子のグループに分割され、それぞれのグループについて個別に推定を行い、その結果を組み合わせることで期待値を求める必要がある。この点については、次のレッスンで、古典的アプローチと量子的アプローチのスケーリングについて取り上げる際に、さらに詳しく検討します。測定におけるこうした複雑さが、このような推定を行うために高効率なコードが必要とされる理由の一つである。このレッスンおよびそれ以降のレッスンでは、この目的のために IBM Quantum のプリミティブ・エスティメータを使用します。

1.4 古典的最適化アルゴリズム

古典的オプティマイザとは、対象関数の極値（通常は最小値）を求めるように設計された古典的アルゴリズムのことです。可能なパラメータ空間の中から、目的の関数を最小化するセットを探すのである。これらは、勾配情報を利用する勾配ベースの手法と、ブラックボックス最適化として動作する勾配なしの手法に大別できる。古典的なオプティマイザの選択は、特に量子ハードウェアにノイズが存在する場合、アルゴリズムの性能に大きな影響を与える可能性がある。この分野で人気のあるオプティマイザには、Adam、AMSGrad、SPSAなどがあり、ノイズの多い環境でも有望な結果を示している。より伝統的なオプティマイザには、COBYLAやSLSQPがある。

一般的なワークフロー（セクション 3.3 で示される）は、scipyの minimize 関数のような最小化器の中のメソッドとして、これらのアルゴリズムのいずれかを使用することである。これは引数として受け取る：

最小化すべき機能これは多くの場合、エネルギー期待値である。しかし、これらは一般的に「コスト関数」と呼ばれている。
検索を開始するパラメータのセット。しばしば $x_0$ または $\theta_0$ と呼ばれる。
引数（コスト関数の引数を含む）。 Qiskitを用いた量子コンピューティングでは、これらの引数にはアンザッツ、ハミルトニアン、および「Estimator」プリミティブが含まれます。後者については、次の小節で詳しく説明します。
最小化の「方法」。これは、パラメータ空間の探索に使用される特定のアルゴリズムを指す。ここで、例えばCOBYLAやSLSQPを指定する。
オプション利用できるオプションは方法によって異なる場合がある。しかし、実質的にすべてのメソッドが含むであろう例は、探索を終了するまでのオプティマイザの最大反復回数である：maxiter'である。

エネルギーを表す曲線と、その最小値を求めるために値が検証されているいくつかの点が示された図。

各反復ステップで、ハミルトニアンの期待値は、多くの測定を行うことによって推定される。この推定エネルギーはコスト関数によって返され、最小化器はエネルギーランドスケープに関する情報を更新する。次のステップを選択するためにオプティマイザが何をするかは、方法によって異なる。勾配を利用し、最も急な下り坂の方向を選択するものもある。また、ノイズを考慮し、真のエネルギーがその方向に沿って減少することを受け入れる前に、コストが大きく減少することを要求する場合もある。

# Example syntax for minimization
# from scipy.optimize import minimize
# res = minimize(cost_func, x0, args=(ansatz, hamiltonian, estimator), method="cobyla",
# options={'maxiter': 200})

1.5 変分原理

この文脈では、変分原理が非常に重要である。変分原理は、いかなる変分波動関数も、基底状態の波動関数がもたらす期待値よりも低いエネルギー（またはコスト）をもたらすことはできないと述べている。数学的に、

E_\text{var}=\langle \psi_\text{var}|H|\psi_\text{var}\rangle \geq E_\text{min}=\langle \psi_\text{0}|H|\psi_\text{0}\rangle

これは、 $H$ の全固有状態の集合 $\{|\psi_0\rangle, |\psi_1\rangle, |\psi_2\rangle, ...|\psi_n \rangle\}$ がヒルベルト空間の完全な基底を形成していることに注意すれば、簡単に検証できる。言い換えれば、どのような状態も、特に $|\psi_\text{var}\rangle$ は、 $H$ のこれらの固有状態の重み付けされた（正規化された）和として書くことができる：

|\psi_\text{var}\rangle=\sum_{i=0}^n c_i |\psi_i\rangle

ここで、 $c_i$ は決定すべき定数であり、 $\sum_{i=0} |c_i|^2 = 1$ 。これは読者への練習として残しておく。最もエネルギーの低い期待値を生み出す変分状態が*、* 真の基底状態の最良の推定値である。

理解度チェック

任意の変分状態 $|\psi_\text{var}\rangle$ について、 $E_\text{var}\geq E_0$ が成り立つことを数学的に証明せよ。

与えられた変分状態をエネルギー固有状態の観点から展開する、
$|\psi_\text{var}\rangle=\sum_{i=0}^n c_i |\psi_i\rangle,$
変分エネルギー期待値は次のように書ける
$\begin{aligned} E_\text{var}&=\langle \psi_\text{var}|H|\psi_\text{var}\rangle =\left(\sum_{i=0}^n c^*_i \langle \psi_i|\right)H\left(\sum_{j=0}^n c_j |\psi_j\rangle\right)\\ &=\left(\sum_{i=0}^n c^*_i \langle \psi_i|\right)\left(\sum_{j=0}^n c_j E_j|\psi_j\rangle\right)\\ &=\sum_{i,j=0}^n c^*_i c_j E_j \langle \psi_i|\psi_j\rangle\\ &=\sum_{i,j=0}^n c^*_i c_j E_j \delta_{i,j}\\ &=\sum_{i=0}^n |c_i|^2 E_i. \end{aligned}$
すべての係数について $0\leq|c_i|^2\leq 1$ 。と書くことができる
$\begin{aligned} E_\text{var}&=\sum_{i=0}^n |c_i|^2 E_i\geq \sum_{i=0}^n |c_i|^2 E_0 = E_0 \sum_{i=0}^n |c_i|^2 = E_0(1) \\ E_\text{var}&\geq E_0 \end{aligned}$

2. 従来のワークフローとの比較

N行N列の行列に興味があるとしよう。行列が非常に大きく、厳密な対角化が不可能だとしよう。さらに、ターゲットとなる固有状態の全体的な構造についてある程度推測できるほど、問題について十分な知識があり、コストやエネルギーをさらに下げられるかどうかを確認するために、最初に推測した状態に似た状態を調べたいとする。これは変分法であり、厳密な対角化が不可能な場合に用いられる方法のひとつである。

2.1 従来のワークフロー

古典的なコンピュータを使えば、次のようになる：

$\vec{\theta}_i$ 、いくつかのパラメーターを変化させながら、状態を推測する。 $|\psi(\vec{\theta}_i)\rangle$ 。この最初の推測はランダムでもよいが、それは好ましくない。私たちは、目の前の問題についての知識を使って、可能な限り私たちの推測を調整したい。
システムがその状態にあるときの演算子の期待値を計算する： $\langle\psi(\vec{\theta}_i)|H|\psi(\vec{\theta}_i)\rangle$
変分パラメータを変更し、 $\vec{\theta}_i\rightarrow \vec{\theta}_{i+1}$ を繰り返す。
変分部分空間における可能な状態のランドスケープに関する蓄積された情報を使って、より良い推測を行い、目標状態に近づいていく。変分原理は、われわれの変分状態が、目標とする基底状態の固有値よりも低い固有値を得ることができないことを保証している。つまり、期待値が低ければ低いほど、基底状態の近似度は高くなる：

\min_{\vec{\theta}} \{ E_{\text{var},i} = \langle\psi(\vec{\theta_i})|H|\psi(\vec{\theta_i})\rangle \} \geq E_0

このアプローチの各ステップの難易度を検証してみよう。パラメーターの設定や更新は計算上簡単で、難しいのは有用で物理的に動機づけられた初期パラメーターを選択することである。以前の反復から蓄積された情報を使って、基底状態に近づくようにパラメータを更新することは、自明なことではない。しかし、これを非常に効率的に行う古典的な最適化アルゴリズムが存在する。この古典的な最適化が高価なのは、多くの反復を必要とする可能性があるからにほかならない。最悪の場合、反復の回数はNに対して指数関数的に増加する可能性がある。最も計算量の多いステップは、与えられた状態 $|\psi(\vec{\theta_i})\rangle$ を用いて行列の期待値を計算することである： $\langle\psi(\vec{\theta_i})|H|\psi(\vec{\theta_i})\rangle.$

$N\times N$ 行列は、 $N$ -要素ベクトルに作用しなければならない。これは最悪の場合、 $O(N^2)$ の乗算演算に相当する。これはパラメーターの反復ごとに行わなければならない。非常に大きな行列の場合、これは計算コストが高い。

2.2 量子ワークフローと通勤パウリ群

では、この計算部分を量子コンピューターに任せることを想像してみよう。この期待値を計算する代わりに、変分アナザッツを使って量子コンピュータ上で状態 $|\psi(\vec{\theta_i})\rangle$ を準備し、測定を行うことで期待値を推定する。

$H$ を測定するのは一般的に容易ではない。例えば、多くの非共通パウリX、Y、Z演算子で構成される。しかし、 $H$ は、 $h_\alpha$ という項の線形結合として書くことができ、各項は簡単に測定可能である（例えば、パウリ演算子や、量子ビットごとに通約するパウリ演算子のグループなど）。ある状態 $|\Psi\rangle$ に対する $H$ の期待値は、構成項の期待値の加重和 $h_\alpha$ である。この式はどの状態 $|\Psi⟩$ でも成り立つが、特に変分状態 $|\psi(\theta_i)\rangle$ でこれを使う。

H = \sum_{\alpha = 1}^T{c_\alpha h_\alpha}

ここで、 $h_\alpha$ は、 IZZX…XIYX のようなパウリ文字列、あるいは、そのような文字列がいくつかあり、互いに通約する。つまり、量子コンピューターでの測定の実態により近い期待値の記述は次のようになる

\langle \Psi |H|\Psi \rangle =\sum_{\alpha} c_\alpha \langle \Psi | h_\alpha|\Psi \rangle.

そして、変分波動関数の文脈で：

\langle \psi(\vec{\theta}_i) |H|\psi(\vec{\theta}_i) \rangle =\sum_{\alpha} c_\alpha \langle \psi(\vec{\theta}_i) | h_\alpha|\psi(\vec{\theta}_i) \rangle

各項 $h_\alpha$ は $M$ 回測定することができ、測定サンプル $s_{\alpha j}$ を $j=1…M$ と共にもたらし、期待値 $\mu_\alpha$ と標準偏差 $\sigma_\alpha$ を返す。これらの項を合計し、合計を通して誤差を伝播させることで、全体的な期待値 $\mu$ と標準偏差 $\sigma$ を得ることができる。

\begin{aligned} \langle \psi(\vec{\theta}_i) |h_\alpha|\psi(\vec{\theta}_i) \rangle &\simeq \mu _\alpha \pm \frac{\sigma_\alpha}{\sqrt{M}} &\qquad \mu_\alpha &=\frac{1}{M}\sum_j s_{\alpha,j} &\qquad \sigma^2_\alpha &=\frac{1}{M-1}\sum_j (s_{\alpha,j}-\mu_\alpha)^2\\ \langle \psi(\vec{\theta}_i) |H|\psi(\vec{\theta}_i) \rangle &\simeq \mu \pm \sigma &\qquad \mu &= \sum_\alpha c_\alpha \mu_\alpha &\qquad \sigma^2&=\sum_\alpha c^2_\alpha \frac{\sigma^2_\alpha }{M} \end{aligned}

これには大規模な乗算も、 $N^2$ のようなスケールを必要とするプロセスも必要ない。その代わり、量子コンピューター上で複数の測定を行う必要がある。もし、あまり数が必要ないのであれば、この方法は効率的かもしれない。それがVQEの量子の部分だ。

しかし、これが効率的でない理由を説明しよう。多くの測定を行う理由のひとつは、非常に精度の高い計算を行うために、推定値の統計的な不確かさを減らすためである。もう一つの理由は、マトリックス全体をカバーするのに必要なパウリ弦の数だ。パウリ行列（と恒等式：X、Y、Z、I）は、与えられた次元のすべての作用素の空間にまたがるので、前と同じように、目的の行列をパウリ作用素の加重和として書けることが保証されている。

H = \sum_{\alpha = 1}^T{c_\alpha h_\alpha}

ここで、 $h_\alpha$ は、 IZZX…XIYX のように、あなたのシステムを記述するすべての量子ビットに作用するパウリ文字列であり、あるいは、そのような文字列がいくつかあり、互いにコミュートしている。 Qiskitはリトルエンディアン表記を用いており、右からの $n^\text{th}$ Pauli演算子が $n^\text{th}$ qubitに作用することを思い出してほしい。つまり、一連のパウリ作用素を測定することで、我々の作用素を測定することができる。

しかし、これらすべてのパウリ演算子を同時に測定することはできない。同じ量子ビットに関連付けられている場合、パウリ演算子（Iを除く）は互いに可換ではない。例えば、第3の量子ビットについてはIとZを同時に測定でき、第1の量子ビットについてはIとXを同時に把握できるため、と ZZXZ を同時に測定 IZIZ することができる。しかし、ZとXは可換ではなく、どちらも0番目の量子ビットに作用するため、これらを ZZZX 同時に測定 ZZZZ することはできない。経験豊富な読者ならご存じかもしれませんが、個々の量子ビットの測定は互いに可換でないにもかかわらず、2つのパウリ演算子の集合としては可換である場合があります。推定器は、量子ビットごとに可換なグループ化演算子に対応する、テンソル積のパウリ測定（基底回転による）を前提としている。したがって、Estimatorを用いて2つの文字列（AとB）のパウリ演算子を同時に推定するには、AとBに含まれる各量子ビットのパウリ演算子が可換でなければならない。つまり、とを ZZXX 同時に測定 ZZZZ することもできないということになる。

そこで、私たちは行列 $H$ 、異なる量子ビットに作用するポーリの和に分解する。その合計のいくつかの要素は一度に測定できます。これを*、可換なパウリス群と*呼びます。非通学期間がどれだけあるかにもよるが、このようなグループが何個も必要になるかもしれない。このようなパウリ弦のグループの数を $N_\text{GCP}$ と呼ぶ。 $N_\text{GCP}$ が小さければ、これはうまくいくだろう。 $H$ に何百万ものグループがある場合、これは役に立たないだろう。

期待値の推定に必要な処理は、 IBM Quantum の「Estimator」というプリミティブにまとめられています。 Estimator の詳細については、『 IBM Quantum® 』のドキュメントにある API リファレンスをご覧ください。 Estimatorを直接使用することも可能ですが、Estimatorは単に最小エネルギーの固有値だけでなく、それ以上の情報を返します。例えば、アンサンブルの標準誤差に関する情報も返されます。したがって、最小化問題の文脈では、コスト関数の中に推定量が含まれていることがよく見られる。 Estimatorの入力と出力について詳しく知りたい場合は、 IBM Quantum のドキュメントにあるこのガイドをご覧ください。

$\vec{\theta_i}$ 、その状態で使用されるパラメータのセットに対する期待値（またはコスト関数）を記録し、パラメータを更新する。時間が経てば、アナザッツによってサンプリングされた状態の部分空間におけるコスト関数の勾配を近似するために、期待値や推定したコスト関数の値を使うことができる。勾配ベースと勾配なしの古典的オプティマイザが存在する。両者とも、複数の局所極小値や、 不毛のプラトーと呼ばれる勾配がゼロに近いパラメータ空間の大きな領域など、潜在的な学習性の問題に悩まされている。

最小値を持つ曲線の図が 2つ。一方では、最小値を探すために点をランダムにチェックし、もう一方では、隣接する2つの点間に線を引くことで勾配を推定します。

2.3 計算コストを決定する要因

VQEは量子化学の難問をすべて解決するわけではない。いいえしかし、すべての計算が優れていることが重要なのではない。私たちは、計算コストを決定するものをシフトさせた。

古典的アプローチと量子変分アプローチを比較した表。どちらの場合も、適切な初期推定値が必要です。従来、計算コストは行列の次元の二乗に比例しますが、量子論的なアプローチでは、交換するパウリ演算子のグループがいくつあるかによって決まります。

我々は、複雑さが行列の次元のみに依存するプロセスから、必要な精度と行列を構成する非共約パウリ作用素の数に依存するプロセスへとシフトした。最後のビットは、古典的なコンピューティングには類似していない。

これらの依存関係に基づき、疎な行列や、少数の非共和パウリ文字列を含む行列では、この処理が有効である可能性がある。これは例えば、相互作用するスピンの系の場合である。密な行列の場合は、あまり役に立たないかもしれない。例えば、化学系はしばしば数百、数千、数百万のパウリ文字列を含むハミルトニアンを持つことがわかっている。この用語数を減らすために興味深い研究がなされている。しかし、化学的なシステムは、このコースで説明する他のアルゴリズムに適しているかもしれない。

理解度チェック

項を含む4量子ビットのハミルトニアンを考える：

IIXX, IIXZ, IIZZ, IZXZ, IXXZ, ZZXZ, XZXZ, ZIXZ, ZZZZ, XXXX

これらの用語をグループに分類し、グループ内のすべての用語を同時に測定できるようにしたい。すべての項が説明できるような、そのようなグループの最小数はいくつですか？

4つのグループに分けて行うことができます。なお、このような解は通常、一意ではないことに注意してください。
IIXX, XXXX, IIZZ, ZZZZ
IIXZ, IZXZ, ZIXZ, ZZXZ
IXXZ
XZXZ

一般的に、VQEを用いた量子化学を困難にしているのは、ハミルトニアンの項の数と、適切な近似解を見つけることのどちらだと思いますか？

化学的なコンテクストに高度に最適化されたアンセッツがあることがわかった。ハミルトニアンの項の数、したがって必要な測定の数は、通常、より多くの問題を引き起こす。

3. ハミルトニアンの例

各ステップで何が起こっているかがわかるように、小さなハミルトニアン行列を使ってこのアルゴリズムを実践してみよう。我々はQiskitパターンのフレームワークを採用する：

ステップ1 : 量子回路と演算子への問題のマッピング - ステップ2 : ターゲットハードウェアの最適化 - ステップ3 : ターゲットハードウェア上での実行 - ステップ4 ：結果の後処理

3.1 ステップ1：問題を量子回路と演算子に写像する

ここでは化学の文脈で定義されたものを使用する。まずは一般的な輸入品から。

# General imports
import numpy as np

# SciPy minimizer routine
from scipy.optimize import minimize

# Plotting functions
import matplotlib.pyplot as plt

ここでも、目的のハミルトニアンは既知であると仮定する。このコースで説明する他の方法の方が、より大きな問題を解くのに効率的だからである。

from qiskit.quantum_info import SparsePauliOp
import numpy as np

hamiltonian = SparsePauliOp.from_list(
    [("YZ", 0.3980), ("ZI", -0.3980), ("ZZ", -0.0113), ("XX", 0.1810)]
)

A = np.array(hamiltonian)
eigenvalues, eigenvectors = np.linalg.eigh(A)
print("The ground state energy is ", min(eigenvalues))

Output:

The ground state energy is  -0.702930394459531

Qiskitにはプレハブのアンサッツの選択肢がたくさんある。 efficient_su2 を使用する。

# Pre-defined ansatz circuit and operator class for Hamiltonian
from qiskit.circuit.library import efficient_su2

# Note that it is more common to place initial 'h' gates outside the ansatz.
# Here we specifically wanted this layer structure.
ansatz = efficient_su2(
    hamiltonian.num_qubits, su2_gates=["h", "rz", "y"], entanglement="circular", reps=1
)

num_params = ansatz.num_parameters
print("This circuit has ", num_params, "parameters")

ansatz.decompose().draw("mpl", style="iqp")

Output:

This circuit has  4 parameters

異なるansätzeは異なるもつれ構造と異なる回転ゲートを持つ。ここに示したものは、エンタングルにはCNOTゲートを使用し、回転にはYゲートとパラメトリックRZゲートの両方を使用している。つまり、4つの変数（RZゲートのパラメータ）にわたってコスト関数を最小化しなければならない。これは規模を拡大することはできるが、無限ではない。 efficient_su2 、デフォルトの3レップを使用して4量子ビットで同様の問題を実行すると、16の変分パラメータが得られる。

3.2 ステップ2: 対象ハードウェア向けに最適化

アナザッツは使い慣れたゲートを使って書かれているが、我々の回路は、各量子コンピュータに実装可能な基底ゲートを利用するためにトランスパイルされなければならない。最も忙しくないバックエンドを選ぶ。

# runtime imports
from qiskit_ibm_runtime import QiskitRuntimeService, Session
from qiskit_ibm_runtime import EstimatorV2 as Estimator

# To run on hardware, select the backend with the fewest number of jobs in the queue
service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)

print(backend)

Output:

<IBMBackend('ibm_torino')>

このハードウェア用に回路をトランスパイルし、トランスパイルされたansatzを可視化することができる。

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

target = backend.target
pm = generate_preset_pass_manager(target=target, optimization_level=3)

ansatz_isa = pm.run(ansatz)

ansatz_isa.draw(output="mpl", idle_wires=False, style="iqp")

Output:

使用されているゲートが変更され、抽象回路の量子ビットが量子コンピューター上の異なる番号の量子ビットにマッピングされていることに注意されたい。結果を意味のあるものにするためには、ハミルトニアンを同一に写像しなければならない。

hamiltonian_isa = hamiltonian.apply_layout(layout=ansatz_isa.layout)

3.3 ステップ3: 対象ハードウェア上で実行する

3.3.1 値の報告

ここでは、これまでの手順で構築した構造、すなわちパラメータ、仮説、およびハミルトニアンを引数として受け取るコスト関数を定義する。また、まだ定義していないEstimatorも使用しています。収束の挙動を確認できるように、コスト関数の推移を追跡するコードを実装しています。

def cost_func(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from Estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (EstimatorV2): Estimator primitive instance
        cost_history_dict: Dictionary for storing intermediate results

    Returns:
        float: Energy estimate
    """
    pub = (ansatz, [hamiltonian], [params])
    result = estimator.run(pubs=[pub]).result()
    energy = result[0].data.evs[0]

    cost_history_dict["iters"] += 1
    cost_history_dict["prev_vector"] = params
    cost_history_dict["cost_history"].append(energy)
    print(f"Iters. done: {cost_history_dict['iters']} [Current cost: {energy}]")

    return energy


cost_history_dict = {
    "prev_vector": None,
    "iters": 0,
    "cost_history": [],
}

手元の問題と目標状態の特徴に関する知識に基づいて、パラメータの初期値を選択できれば非常に有利である。そのような知識を前提とせず、ランダムな初期値を使用する。

x0 = 2 * np.pi * np.random.random(num_params)

# This required 13 min, 20 s QPU time on an Eagle processor, 28 min total time.
with Session(backend=backend) as session:
    estimator = Estimator(mode=session)
    estimator.options.default_shots = 10000

    res = minimize(
        cost_func,
        x0,
        args=(ansatz_isa, hamiltonian_isa, estimator),
        method="cobyla",
        options={"maxiter": 50},
    )

Output:

Iters. done: 1 [Current cost: 0.010575798722044727]
Iters. done: 2 [Current cost: 0.004040015974440895]
Iters. done: 3 [Current cost: 0.0020213258785942503]
Iters. done: 4 [Current cost: 0.18723082446726014]
Iters. done: 5 [Current cost: -0.2746792152068885]
Iters. done: 6 [Current cost: -0.3094547651648519]
Iters. done: 7 [Current cost: -0.05281985428356641]
Iters. done: 8 [Current cost: 0.00808560303514377]
Iters. done: 9 [Current cost: -0.0014821685303514388]
Iters. done: 10 [Current cost: -0.004759824281150161]
Iters. done: 11 [Current cost: 0.09942328705995292]
Iters. done: 12 [Current cost: 0.01092366214057508]
Iters. done: 13 [Current cost: 0.05017497496069776]
Iters. done: 14 [Current cost: 0.13028868414310696]
Iters. done: 15 [Current cost: 0.013747803514376994]
Iters. done: 16 [Current cost: 0.2583072432944498]
Iters. done: 17 [Current cost: -0.14422125655131562]
Iters. done: 18 [Current cost: -0.0004950150347678081]
Iters. done: 19 [Current cost: 0.00681082268370607]
Iters. done: 20 [Current cost: -0.0023377795527156544]
Iters. done: 21 [Current cost: 0.6027665591169237]
Iters. done: 22 [Current cost: 0.00596641373801917]
Iters. done: 23 [Current cost: -0.008318769968051117]
Iters. done: 24 [Current cost: -0.00026683306709265246]
Iters. done: 25 [Current cost: -0.007648222843450479]
Iters. done: 26 [Current cost: 0.004121086261980831]
Iters. done: 27 [Current cost: -0.004075019968051117]
Iters. done: 28 [Current cost: -0.004419369009584665]
Iters. done: 29 [Current cost: 0.213185460054037]
Iters. done: 30 [Current cost: -0.06505919572162797]
Iters. done: 31 [Current cost: -0.5334241316590271]
Iters. done: 32 [Current cost: 0.00218370607028754]
Iters. done: 33 [Current cost: 0.09579352143666908]
Iters. done: 34 [Current cost: -0.009274800319488819]
Iters. done: 35 [Current cost: -0.44395141360688106]
Iters. done: 36 [Current cost: 0.011747104632587858]
Iters. done: 37 [Current cost: -0.003344149361022364]
Iters. done: 38 [Current cost: 0.19138183916486304]
Iters. done: 39 [Current cost: 0.013513931813145209]

生の出力を見ることができる。

res

Output:

 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: -0.5334241316590271
       x: [ 1.024e+00  6.459e+00  3.625e+00  4.007e+00]
    nfev: 39
   maxcv: 0.0

3.4 ステップ4：結果の後処理

もしこの手順が正しく終了すれば、辞書の値はそれぞれ解ベクトルと関数評価総数に等しくなるはずである。これは簡単に確認できる：

cost_history_dict

Output:

{'prev_vector': array([1.02397956, 6.45886604, 3.62479262, 4.00744128]),
 'iters': 39,
 'cost_history': [np.float64(0.010575798722044727),
  np.float64(0.004040015974440895),
  np.float64(0.0020213258785942503),
  np.float64(0.18723082446726014),
  np.float64(-0.2746792152068885),
  np.float64(-0.3094547651648519),
  np.float64(-0.05281985428356641),
  np.float64(0.00808560303514377),
  np.float64(-0.0014821685303514388),
  np.float64(-0.004759824281150161),
  np.float64(0.09942328705995292),
  np.float64(0.01092366214057508),
  np.float64(0.05017497496069776),
  np.float64(0.13028868414310696),
  np.float64(0.013747803514376994),
  np.float64(0.2583072432944498),
  np.float64(-0.14422125655131562),
  np.float64(-0.0004950150347678081),
  np.float64(0.00681082268370607),
  np.float64(-0.0023377795527156544),
  np.float64(0.6027665591169237),
  np.float64(0.00596641373801917),
  np.float64(-0.008318769968051117),
  np.float64(-0.00026683306709265246),
  np.float64(-0.007648222843450479),
  np.float64(0.004121086261980831),
  np.float64(-0.004075019968051117),
  np.float64(-0.004419369009584665),
  np.float64(0.213185460054037),
  np.float64(-0.06505919572162797),
  np.float64(-0.5334241316590271),
  np.float64(0.00218370607028754),
  np.float64(0.09579352143666908),
  np.float64(-0.009274800319488819),
  np.float64(-0.44395141360688106),
  np.float64(0.011747104632587858),
  np.float64(-0.003344149361022364),
  np.float64(0.19138183916486304),
  np.float64(0.013513931813145209)]}

fig, ax = plt.subplots()
x = np.linspace(0, 10, 50)

# Define the constant function
constant = -0.7029
y_constant = np.full_like(x, constant)
ax.plot(
    range(cost_history_dict["iters"]), cost_history_dict["cost_history"], label="VQE"
)
ax.set_xlabel("Iterations")
ax.set_ylabel("Cost")
ax.plot(y_constant, label="Target")
plt.legend()
plt.draw()

Output:

IBM クァンタムには、VQEに関連する他のスキルアップ・サービスもある。 VQEを実践する準備ができたら、チュートリアルをご覧ください： VQEによるハイゼンベルグ鎖の基底状態エネルギー推定をご覧ください。分子ハミルトニアンの作成についてさらに詳しく知りたい場合は、 VQEによる量子化学コースのこのレッスンを参照してください。 VQEのような変分アルゴリズムがどのように機能するのかをより深く理解したい方には、変分アルゴリズム設計のコースをお勧めします。

理解度チェック

このセクションでは、ハミルトニアンから基底状態のエネルギーを計算した。これを例えば、分子の形状を決定することに応用するとしたら、どのように拡張するのだろうか？

原子間の間隔や結合間の角度の変数を導入する必要があるだろう。これらを変化させる必要があるだろう。エネルギーを記述する作用素は確かに幾何学に依存するので）これらのバリエーションごとに、新しいハミルトニアンを作り出すことになる。このようなハミルトニアンが生成され、量子ビットにマッピングされるごとに、上で行ったような最適化を行う必要がある。収束した最適化問題の中で、最も低いエネルギーを生み出す形状は、自然が採用したものだろう。これは上に示したものよりもかなり複雑だ。このような計算は、最も単純な分子、 $\text{H}_2$ 、ここで行われている。

4. VQEと他の手法との関係

このセクションでは、オリジナルのVQEアプローチの長所と短所をレビューし、他の最新のアルゴリズムとの関係を指摘する。

4.1 VQEの長所と短所

いくつかの長所はすでに指摘されている。 IBM iの柱は以下のとおりです。

最新のハードウェアへの適合性：量子アルゴリズムの中には、より低いエラーレートを必要とするものがあり、大規模なフォールトトレランスに近づいている。 VQEはそうではなく、現在の量子コンピュータで実装可能だ。
浅い回路：VQEは比較的浅い量子回路を採用することが多い。このため、VQEは累積ゲート・エラーの影響を受けにくく、多くのエラー緩和技術に適している。もちろん、回路は必ずしも浅いとは限らない。
汎用性：VQEは、（原理的には）固有値／固有ベクトルの問題としてキャストできるあらゆる問題に適用できる。 VQEには、いくつかの問題に対して実用的でなかったり、不利であったりする多くの注意点がある。以下にその一部を紹介する。

VQEの弱点や実用的でない問題点については、前述の通りである。次のものが含まれます。

ヒューリスティックな性質：VQEは正しい基底状態エネルギーへの収束を保証するものではなく、その性能はアナザッツと最適化手法の選択に依存するからである [1-2]。望む基底状態に必要なエンタングルメントを欠く貧弱なアナザッツが選ばれた場合、古典的オプティマイザはその基底状態に到達できない。
多数のパラメータの可能性：非常に表現力豊かなansatzは、最小化の反復に非常に時間がかかるほど多くのパラメータを持つ可能性がある。
計算負荷が大きい：VQEでは、ハミルトニアンの各項の期待値を推定するために推定器が用いられる。注目されるハミルトニアンの大半には、同時に推定できない項が含まれている。これにより、ハミルトニアンが複雑な大規模システムでは、VQEがリソースを大量に消費する可能性がある [1]。
ノイズの影響：古典的なオプティマイザが最小値を探索するとき、ノイズの多い計算はオプティマイザを混乱させ、真の最小値から遠ざけたり、収束を遅らせたりします。この解決策として考えられるのは、 IBM の最先端のエラー緩和・エラー抑制技術 [2-3] を活用することである。
不毛のプラトー：このような勾配が消失する領域 [2-3] は、ノイズがない場合でも存在するが、ノイズがあると、このような不毛な領域では、ノイズによる期待値の変化が、パラメータの更新による変化よりも大きくなる可能性があるため、より厄介になる。

4.2 他のアプローチとの関係

Adapt-VQE

ADAPT-VQE（Adaptive Derivative-Assembled Pseudo-Trotter Variational Quantum Eigensolver）アルゴリズムは、オリジナルのVQEアルゴリズムを改良したもので、特に量子化学における量子シミュレーションの効率、精度、スケーラビリティを向上させるように設計されています。

このレッスン全体を通して説明されるオリジナルのVQEアルゴリズムは、システムの基底状態を近似するために、あらかじめ定義された固定のansatzを使用する。私たちの場合、Y回転ゲートとRZ回転ゲートを使い、 efficient_su2、1回だけ繰り返した。 RZゲートのパラメータは変わったが、このansatzの構造と使用するゲートは変わらなかった。

ADAPT-VQEは、適応的なansatz構築によってVQEの限界に対処する。 ADAPT-VQEは、固定されたansatzから始めるのではなく、動的に反復的にansatzを構築する。各ステップで、あらかじめ定義されたプール（フェルミオン励起演算子など）から、エネルギーに対して最大の勾配を持つ演算子を選択する。これにより、最も影響力のある演算子のみが追加され、コンパクトで効率的なansatz [4-6] につながります。このアプローチはいくつかの有益な効果をもたらす：

回路深度の削減 ：ADAPT-VQEは、漸進的にansatzを成長させ、必要な演算子のみに焦点を当てることで、従来のVQEアプローチと比較してゲート演算を最小限に抑えます [5,7]。
精度の向上 ：適応的な性質により、ADAPT-VQEは各ステップでより多くの相関エネルギーを回復することができ、従来のVQEが苦戦するような相関の強いシステムに対して特に効果的です [8,9]。
スケーラビリティとノイズ耐性 ：コンパクトなansatzは、ゲートエラーの蓄積を減らし、計算オーバーヘッドを削減し、最小化しなければならない変分パラメータの数を制限します。

ADAPT-VQEはまだ完璧ではない。場合によっては、局所的な極小値にとらわれたり、速度が低下したりすることがあり、また、パラメータ化のしすぎに悩まされることもある。また、多くのゲート構造で勾配の計算とパラメーターの最適化を必要とするため、かなりリソースを消費する可能性がある。

量子位相推定（QPE）

QPEはVQEと目的は似ているが、実施方法は大きく異なる。 QPEは、一般的に深い量子回路と高レベルのコヒーレンスを必要とするため、耐障害性の高い量子コンピュータを必要とする。 QPEが実装されれば、VQEよりも精度が高くなるだろう。その違いを説明する一つの方法は、回路の深さの関数としての精度である。 QPEは、回路深さを $O(1/\epsilon)$ [10] のようにスケーリングしながら、精度 $\epsilon$ を達成する。 VQEは、同じ精度を達成するために $O(1/\epsilon^2)$ のサンプルを必要とする [10,11]。

クリロフ、SQD、QSCI、およびこのコースのその他の者

VQEは、量子コンピュータの操作だけでなく、アルゴリズムのかなりの部分を古典コンピュータに依存している量子アルゴリズムの確立に貢献した。このコースの残りの部分では、そのようなアルゴリズムに焦点を当てる。ここでは、単にVQEと比較対照するために、そのいくつかをざっと説明する。これらについては、この後のレッスンで詳しく説明する。

クリロフ量子対角化(KQD)

クリロフ部分空間法とは、行列を部分空間に投影してその次元を縮小し、最も重要な特徴を維持しながら扱いやすくする方法である。この方法の1つのコツは、これらの特徴を保持する部分空間を生成することである。この部分空間を生成することは、量子コンピュータで確立されたトロッター 化と呼ばれる方法と密接に関係していることが判明した。

量子クリロフ法にはいくつかのバリエーションがあるが、一般的には次のようなアプローチである：

量子コンピュータを使って、トロッター化による部分空間（クリロフ部分空間）を生成する
目的の行列をそのクリロフ部分空間に射影する
古典的なコンピュータを使って、新しい投影ハミルトニアンを対角化する

サンプリングに基づく量子対角化（SQD）

サンプリングに基づく量子対角化(SQD) は、クリロフ法に関連しており、主要な特徴を保持しながら対角化する行列の次元を縮小しようとするものである。 SQDはこれを次のように行っている：

基底状態の初期推測から始め、その基底状態で系を準備する。
この状態を構成するビット列をサンプリングするにはサンプラーを使う。
目的の行列を投影する部分空間として、サンプラーの計算基底状態のコレクションを使用する。
古典的なコンピュータを使って、投影された小さい行列を対角化する。

これは、アルゴリズムの実質的な構成要素に古典計算と量子計算を活用するという点で、VQEと関連している。両者とも、良い初期推測またはansatzを準備する必要があることも共通している。しかし、SQDにおける古典計算機と量子計算機の間の仕事の分配は、クリロフ法に似ている。

実際、クリロフ法とSQDは最近、サンプリングに基づくクリロフ量子対角化（SKQD）法として統合された [12]。

量子亜空間構成相互作用

量子選択配置相互作用（QSCI） [13] は、古典的対角化のための部分空間を生成するための重要な計算基底状態を識別するために、試行波動関数をサンプリングすることにより、ハミルトニアンの近似基底状態を生成するアルゴリズムである。 SQDもQSCIも、量子コンピュータを使って縮小部分空間を構築する。 QSCIのさらなる強みは、特に化学の問題における状態の準備にある。これは、時間発展した状態 [14] や化学にヒントを得た解法のセットを使用するなど、さまざまな戦略を活用している。効率的な状態準備に重点を置くことで、QSCIは化学ハミルトニアンの量子計算コストを削減すると同時に、高い忠実度を維持し、量子状態サンプリング技術から得られるノイズの頑健性を活用している [15]。 QSCIはまた、より良い結果を得るために、より多くのansätzeを提供する適応的な構築技術を提供します。

化学問題に対するQSCIのデフォルトのワークフローは以下の通りである：

お好みのソフトウェア（ SciPy など）を使用して分子ハミルトニアンを構築します。
QSCIアルゴリズムは、適切な初期状態と、あらかじめ選択されたパラメータセットを持つ化学に触発されたansatzを選択することによって準備します。
重要な基底状態をサンプリングし、古典コンピュータを使ってハミルトニアンを対角化し、基底状態のエネルギーを求める。
後処理手法として、構成回復 [16] や対称性ポストセレクション [15] がよく用いられる。
オプションとして、適応型QSCIのワークフローは、ランダムな初期状態を持つより多くのansätzeを使用することにより、 step2 から step3 までの追加の最適化ループを持つ。