変分量子固有値ソルバー (VQE)

このモジュールでは、 Python 環境が動作し、以下のパッケージの最新版がインストールされている必要があります：

qiskit
qiskit_ibm_runtime
qiskit-aer
qiskit.visualization
numpy
pylatexenc

これらのパッケージをセットアップしてインストールするには、 Qiskitのインストールガイドをご覧ください。実際の量子コンピュータでジョブを実行するには、 IBM Cloud アカウントの設定ガイドの手順に従って、 IBM Cloud アカウントを設定する必要があります。

このモジュールはテストされ、約8分のQPU時間を使用した。これは概算であり、実際の使用量は異なる場合があります。

# Uncomment and modify this line as needed to install dependencies
#!pip install 'qiskit>=2.1.0' 'qiskit-ibm-runtime>=0.40.1' 'qiskit-aer>=0.17.0' 'numpy' 'pylatexenc'

概要

20世紀初頭に量子力学的モデルが開発されて以来、科学者たちは、電子は原子核の周りの固定された経路をたどらず、むしろ軌道と呼ばれる確率の領域に存在することを理解してきた。これらの軌道は、電子が占有できる特定の個別のエネルギー準位に対応している。電子は当然、基底状態として知られる、利用可能な最低のエネルギー準位に存在する。しかし、電子が十分なエネルギーを吸収すれば、より高いエネルギー準位に飛び、励起状態に入ることができる。この励起状態は一時的なもので、電子は最終的に低いエネルギー準位に戻り、吸収したエネルギーを放出する。エネルギーの吸収と放出というこの基本的なプロセスは、原子がどのように相互作用し、結合を形成するかを理解する上で重要である。

原子が集まって分子を形成するとき、その原子軌道は結合して分子軌道を形成する。これらの分子軌道内の電子の配置とエネルギー準位が、得られる分子の性質と化学結合の強さを決定する。例えば、2個の水素原子から水素分子（ $H_2$ ）が形成されるとき、それぞれの原子の電子が原子軌道を占める。原子が互いに近づくと、これらの原子軌道が重なり合い、結合して新しい分子軌道が形成される。1つはエネルギーが低く（結合軌道）、もう1つはエネルギーが高い（反結合軌道）。各水素原子から1個ずつ、合計2個の電子が、よりエネルギーの低い結合軌道を優先的に占有し、 $H_2$ 分子をつなぎとめる安定した共有結合を形成する。分離した原子と形成された分子との間のエネルギー差、特に分子軌道の電子のエネルギーが、結合の安定性と特性を決定する。

以下のセクションでは、 $H_2$ 分子を中心に、この分子形成の過程を探っていく。私たちは、古典的な最適化技術と本物の量子コンピューターを組み合わせて、この単純だが基本的なプロセスのエネルギーを求める。この実験は、量子計算が計算化学の問題解決にどのように応用できるかを実践的に示すものであり、電子エネルギーの役割についての洞察を提供するものである。

VQE - 固有値問題のための変分量子アルゴリズム

化学における近似手法 - 変分原理と基底関数

エルヴィン・シュレーディンガーの量子力学への貢献は、新しい電子モデルの導入にとどまらない。基本的には、有名な時間依存シュレーディンガー方程式を開発することによって、波動力学を確立したのである：

i\hbar \frac{d}{dt}|\psi\rangle = \hat{H}|\psi\rangle

ここで、 $\hat{H}$ は系の全エネルギーを表すハミルトニアン演算子であり、 $|\psi\rangle$ は系の量子状態に関するすべての情報を含む波動関数である。 (注: $\frac{d}{dt}$ は全時間微分であり、ここではエネルギー固有値 $E$ は明示的に含まない)

しかし、原子や分子の許容エネルギー準位の決定など、多くの実用的なアプリケーションでは、代わりに、定常状態を仮定することによって時間依存型から導かれる時間非依存型のシュレーディンガー方程式（エネルギー固有値方程式）を用いる。定常状態とは、ある粒子が空間上のある点に存在する確率密度が時間とともに変化しない量子状態のことである。

$\hat{H}|\psi\rangle = E|\psi\rangle$

この形式では、 $E$ は量子状態 $|\psi\rangle$ に対応するエネルギー固有値を表す。ハミルトニアンには、電子と原子核の運動エネルギー、電子と原子核間の引力、電子間の反発力など、さまざまなエネルギー寄与が含まれる。

エネルギー固有値方程式を解くことで、原子や分子系の量子化されたエネルギー準位を計算することができる。しかし分子の場合、電子の空間分布を記述する波動関数 $\Psi$ が複雑で高次元であるため、これを正確に解くことは難しい。

その結果、科学者たちは実用的で正確な解を得るために近似的なテクニックを用いる。この作品では、2つの重要な手法に焦点を当てる：

変分原理

この方法は波動関数を近似し、目標エネルギー（通常は系の基底状態エネルギー）にできるだけ近づくように調整する。変分原理の背後にある重要な考え方は単純だ：
- 波動関数 $\Psi_\text{trial}$ （「試行関数」）を推測すると、そこから計算されるエネルギーは常に系の基底状態エネルギー( $E_0$ )と等しいか、それ以上になる。 $E_\text{approx} = \frac{\langle \Psi_\text{trial}|\hat{H}|\Psi_\text{trial}\rangle}{\langle \Psi_\text{trial}|\Psi_\text{trial}\rangle} \geq E_0$
- 試行関数 $|\Psi_\text{trial}(\theta)\rangle$ のパラメータ $\theta$ を調整することで、基底状態のエネルギーをより良く近似することができる。
- その精度は、試行波動関数 $\Psi_\text{trial}$ の選択に大きく依存する。試行関数の選択を誤ると、正確さとは程遠いエネルギー推定値になる可能性がある。
基底セット近似

第二の近似法は、波動関数を構築する段階で登場する。量子化学では、分子のシュレーディンガー方程式を正確に解くことはほとんど不可能である。その代わりに、複雑な多電子の波動関数を、より単純であらかじめ定義された数学的関数から構築して近似する。基底セットとは、基本的に、分子内の原子を中心とした既知の数学関数の集まりであり、系内の電子の形状と挙動を表現するためのビルディングブロックとして使用される。標準的なレゴブロックのコレクションだけを使って、詳細な彫刻を再現しようとするようなものだと考えてほしい。ブロックの種類やサイズが多ければ多いほど（基本セットが大きければ大きいほど）、元の形に正確に近づけることができる。

これらの基底関数は、水素原子のような単純系の解析解から着想を得て、ガウス関数やスレーター型関数のような形をとることが多いが、それでも近似であることに変わりはない。理論的には "厳密 "だが扱いにくい完全な分子軌道を扱う代わりに、これらの基底関数の線形結合（係数の和）として表現する。この方法は、基底関数が原子軌道に似ている場合、原子軌道の線形結合（LCAO）アプローチとして知られている。この線形結合の係数を最適化することで、選択された基底セットの制限内で、最良の近似波動関数とエネルギーを見つけることができる。
- 基底セットに含まれる関数の数が多ければ多いほど、近似度は向上するが、その代償として計算量が増える。
- 小さな基底セットは大まかな推定値を与えるが、大きな基底セットはより正確な結果を与える。

要約すると、計算を実行可能にし、計算コストを削減するために、波動関数を近似することによって変分原理を使用し、計算の複雑さを軽減し、エネルギーを最小化するための反復最適化を可能にする。一方、基底セット法は、連続波動関数を直接解くのではなく、あらかじめ定義された関数の組み合わせとして原子軌道を表現することで、計算を単純化する。

理解度チェック

$\Psi_\text{trial}(\alpha,x) = Ae^{- \alpha x^2}$ $A$ は正規化定数、 $\alpha$ は調整可能なパラメータである。

(a) 試行波動関数を正規化するため、 $\int_{-\infty}^{\infty} |\Psi_\text{trial}|^2 dx = 1$ となるような $A$ を決定する。

与えられた試行波動関数を正規化するには：
$\int_{-\infty}^{\infty} |\Psi_\text{trial}|^2 dx = \int_{-\infty}^{\infty} A^2 e^{-2 \alpha x^2} dx = 1$
ガウス積分を使う：
$\int_{-\infty}^{\infty} e^{-a x^2} dx = \sqrt{\frac{\pi}{a}} \text{, for } a>0$
set $a = 2\alpha$ then get： $A^2\sqrt{\frac{\pi}{a}} = 1$ $\therefore A = (\frac{2\alpha}{\pi})^{1/4}$

(b) ハミルトニアン $\hat{H}$ の期待値を計算せよ。これは、 $\hat{H} = -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + V(x)$ で与えられ、ここで $V(x) = \frac{1}{2}m\omega^2x^2$ であり、これは単純な調和振動子のポテンシャルに対応する。

調和振動子のハミルトニアンは次の通りである：
$\hat{H} = -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + \frac{1}{2} m \omega^2 x^2$
運動エネルギーの期待値
$\langle T \rangle = -\frac{\hbar^2}{2m} \int_{-\infty}^{\infty} \Psi_\text{trial}^* \frac{d^2}{dx^2} \Psi_\text{trial} dx$
二次微分を取る：
$\frac{d}{dx} \Psi_\text{trial} = -2\alpha x A e^{-\alpha x^2}$
$\frac{d^2}{dx^2} \Psi_\text{trial} = A e^{-\alpha x^2} (4\alpha^2 x^2 - 2\alpha)$
したがって、
$T = -\frac{\hbar^2}{2m} \int_{-\infty}^{\infty} A^2 e^{-2\alpha x^2} (4\alpha^2 x^2 - 2\alpha) dx$
標準的なガウス積分の結果を使用：
$\langle T \rangle = \frac{\hbar^2 \alpha}{2m}$
位置エネルギーの期待値
$\langle V \rangle = \frac{1}{2} m \omega^2 \int_{-\infty}^{\infty} x^2 |\Psi_\text{trial}|^2 dx$
使っている：
$\int_{-\infty}^{\infty} x^2 e^{-a x^2} dx = \frac{\sqrt{\pi}}{2a^{3/2}}$
となる：
$\langle V \rangle = \frac{m \omega^2}{4\alpha}$
総エネルギーの期待値
$\therefore E_\text{approx}(\alpha) = \frac{\hbar^2 \alpha}{2m} + \frac{m \omega^2}{4\alpha}$

$\alpha$ を最小エネルギーとなるよう最適化する
区別する
$\frac{d}{d\alpha} \left( \frac{\hbar^2 \alpha}{2m} + \frac{m \omega^2}{4\alpha} \right) = 0$
解決する：
$\frac{\hbar^2}{2m} - \frac{m \omega^2}{4\alpha^2} = 0$
$\alpha_\text{opt} = \frac{m\omega}{2\hbar}$
$\alpha_\text{opt}$ を $E_\text{approx}$ に代入する：
$\therefore E_\text{approx} = \frac{\hbar \omega}{2}$
これは量子調和振動子の基底状態のエネルギーと正確に一致する。

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

変分量子固有値ソルバー（VQE）は、 $H+H = H_2$ プロセスを探求するために使用する主な手法である。ここでは、VQEとは何か、どのように機能するのかを見ていこう。しかし、まず一歩立ち止まって、チェックインの質問を通して非常に重要なことを一つ見てみよう。

理解度チェック

化学の問題に対して、すでにそれだけの戦略があるのなら、なぜ量子コンピューターが必要なのか？量子コンピュータと古典コンピュータを併用する目的は何ですか？

量子コンピュータは、量子状態の指数関数的なスケーリングのために古典的なコンピュータが苦手とする問題に取り組むことで、化学に革命を起こすチャンスがある。リチャード・ファインマンは、自然をシミュレートするためには、計算も量子的でなければならないと指摘した（[参考文献1] ）。
例えば、最も単純な基底セット( STO-3G )でカフェインをシミュレートする場合、 $10^{48}$ ビットが必要となり、観測可能な宇宙に存在する星の総数( $10^{24}$ )よりも遥かに大きい( [参考文献2] )。量子コンピューターは160量子ビットでカフェインの電子軌道を記述できる。
量子コンピューターは、重ね合わせともつれを使って量子相互作用を自然に処理するため、正確な分子シミュレーションを可能にする有望な方法となる。さらに、量子コンピュータ（電子シミュレーション）と古典コンピュータ（データ前後処理、アルゴリズムプロセス管理、最適化など）の両方の利点を組み合わせることができる。これらにより、材料探索、薬物設計、反応予測を強化し、コストのかかる試行錯誤の実験を減らすことが期待される。 [レフ3][ref 4]
量子コンピューターがなぜ化学の問題に必要なのか、なぜ量子コンピューターと古典コンピューターの両方のリソースを使う必要があるのかを知りたい方は、以下の記事をご覧ください：

さて、VQEに話を戻そう。

VQEは量子コンピュータと古典コンピュータのパワーを組み合わせたもので、基本的には変分原理を使って系の基底状態エネルギーを求める。 VQEを理解するには、まず3つの部分に分けて考える：

(量子) 観測可能量：分子ハミルトニアン（分子のエネルギー）

VQEでは、分子／原子のハミルトニアンは観測可能であり、実験によってその値を測定することができる。我々の目標は、分子の可能な限り低いエネルギー（基底状態エネルギー）を見つけることである。そのために、パラメータ化された量子回路（ansatz）によって生成された試行量子状態を用いる。観測可能なものを測定し、可能な限り低いエネルギーに達するまで量子状態を最適化する。

分子ハミルトニアンに使用される基底セットは、必要な量子ビット数を決定し、VQEの精度に直接影響する。適切な基底セットを選択することは、効率と精度のバランスをとる上で非常に重要である。基底セットを変えずに計算を簡略化するには、対称性の付与やアクティブ空間の縮小といった戦略を用いることができる。多くの分子は対称的な形をしており（蝶や雪の結晶のような）、これはいくつかの部分が同じように振る舞うことを意味している。すべてを別々に計算する代わりに、ユニークな部分だけに集中することで、量子リソースを節約し、対称性を活用することができる。活性空間の縮小では、すべての電子が分子エネルギーに大きな影響を与えるわけではないので、重要な軌道のみを考慮する。原子核に近い電子はほとんど変化しないが、他の電子は結合に影響を与える。これらの方法を適用することで、精度を維持しながらVQEをより効率的にすることができる。

上記の適切な基底セットと戦略を使って分子ハミルトニアンを得たら、このハミルトニアンを量子コンピューターに適したものに変換する必要がある。問題をパウリ作用素にマッピングするのは非常に複雑な場合がある。これは特に、区別できない粒子（電子）を扱う量子化学において当てはまり、量子ビットは区別可能だからである。マッピングの詳細についてはここでは触れないが、以下の資料を参照されたい。量子演算子への問題のマッピングに関する一般的な議論は、 Quantum Computing in Practiceに掲載されている。化学の問題を量子演算子にマッピングすることに関するより詳細な議論は、 VQEによる量子化学に掲載されている。

このモジュールでは、 $H$ と $H_2$ の適切な（1量子ビットの）ハミルトニアンを提供し、量子コンピュータの使用に集中できるようにする。これらの 1 量子ビットのハミルトニアンは、 STO-6G 基底関数セットと Jordan-Wigner マッピングを使用して作成されます。これは、1 つのスピン軌道の占有を 1 量子ビットの占有にマッピングするため、最も直接的で物理的な解釈が最も単純です。また、ハミルトニアンの対称性を利用した量子ビット削減手法も使用しました。これは、スピン占有の動作パターンを使用して量子ビットの数を削減するものです。 $H_2$ 分子については、2つの水素原子間の距離を 0.735 $\mathring A$ と仮定する。

(量子) アントザス：試行波動関数（量子回路を用いて自明な量子状態を構築する方法）

VQEの場合、アンサッツ（複数形：ansätze）は2つの重要な要素で構成されている。 1つ目は初期状態の準備で、変分パラメータを持たない量子ゲートを適用して量子ビットの状態を設定する。つ目のコンポーネントはパラメータ化量子回路で、ラジオのダイヤルのようにパラメータを調整できる特殊な量子回路である。これらのパラメータは、最後の部分である古典的オプティマイザーで使用され、可能な限り最良の基底状態に到達する助けとなる。

変分原理のセクションでは、試行状態の質が変分アルゴリズムの結果の質に影響することを学んだ。つまり、VQEでは良いアサッツを選ぶことが重要なのだ。もう一度言うが、これは豊かで複雑なテーマだ。ここでは、さまざまなタイプのアサッツやその起源については触れない。パラメータ化された量子回路とansatzについてもっと学びたい方は、変分アルゴリズム設計コースの ansatzと変分形式のレッスンで、ansätzeの詳細な説明と例を見ることができます。

このモジュールでは1量子ビットのハミルトニアンを使うので、アンサッツとして1量子ビットのパラメタライズされた量子回路が必要である。次のセクションでは、3つのタイプの1量子ビット解析を見ていく。私たちはそれらを比較し、アサッツを選択する際の重要な考慮点について議論する。

(古典)オプティマイザー：量子回路の微調整

量子コンピュータがアナザッツから観測値のエネルギーを測定すると、アナザッツのパラメータとエネルギー値が古典的オプティマイザに送られ、チューニングが行われる。この最適化プロセスは、通常、 SciPy のような汎用科学パッケージを使って、古典的なコンピュータ上で実行される。

古典的なオプティマイザは、測定されたエネルギーをコスト関数として扱う。最適化問題において、コスト関数（目的関数と呼ばれることもある）とは、特定の解がどの程度「良い」かを測る数学的関数である。オプティマイザの目標は、このコスト関数を最小化するパラメータセットを見つけることである。分子の基底状態エネルギーを求めるという文脈では、エネルギーそのものがコスト関数として機能する。私たちは、可能な限り低いエネルギーをもたらす量子回路のパラメータ（私たちの「解」）を見つけたいのだ。古典的なオプティマイザは、この測定されたエネルギー値（コスト）を用いて、量子ansatzの次の最適化パラメータセットを決定する。更新されたパラメータは量子回路に送られ、このプロセスが繰り返される。古典的なオプティマイザは、反復のたびに、あらかじめ定義された収束基準が満たされるまで、エネルギーを削減しようと（コスト関数を最小化しようと）パラメータを調整し、理想的には、可能な限り低いエネルギー（その結合距離と基底セットに対する分子の基底状態に対応する）が見つかるようにする。

SciPy のような科学的パッケージが提供する最適化戦略は数多くある。詳細は、変分アルゴリズム設計コースの最適化ループのレッスンでご覧いただけます。ここでは、複雑なエネルギー・ランドスケープに適した最適化アルゴリズムであるCOBYLA（Constrained Optimization BY Linear Approximations）を使用する。特に、COBYLAは研究対象の関数の勾配を計算しようとはしません。これを勾配なしオプティマイザと呼びます。目を閉じて山脈の最高峰を見つけようとしていると想像してほしい。全体を見渡すことができないので、上か下かを確認しながら、さまざまな方向に少しずつ歩を進める。 COBYLAも同じように動作する。パラメータ空間を移動し、さまざまな値をテストし、最適なものを見つけるまで徐々に結果を改善する。

これでVQE計算の準備は整った。そのために、全体的なプロセスをまとめた以下のチェックイン問題を試してみてほしい。

理解度チェック

VQEプロセスの概要を完成させるため、空欄に正しい用語を記入し、クリックして答えを確認してください。

VQEは変分量子アルゴリズムであり、分子の(2)を見つけるために使われる、(1)と(2)の力を組み合わせたものである。このプロセスは、系の全エネルギーを表し、量子測定の観測量として働く(3)Ⓐを定義することから始まる。次に、分子の試行波動関数を表すパラメータを調整可能な量子回路である(4)㊟を用意する。これらのパラメータは、測定されたエネルギーを最小化するようにパラメータを反復的に調整する古典的アルゴリズムである(5)Ⓐを用いて最適化される。上記の議論では、微分計算を必要とせずにansatzパラメータを改良する(6)オプティマイザを使用した。このプロセスは(7)に達するまで続けられ、分子の可能な限り低いエネルギーを見つけたことになる。

ワードバンク

古典オプティマイザ
基底状態エネルギー
ハードウェア効率
アナザッツ
分子ハミルトニアン
コビラ
量子コンピューティング
集中性

1 → 量子コンピューティング
2 →基底状態のエネルギー
3 → 分子ハミルトニアン
4 → アンサッツ
5 → 古典的オプティマイザー
6 → コビラ
7 → 収束

VQEを用いて水素原子の基底状態エネルギーを計算する

では、学んだことを使って水素原子の基底状態エネルギーを計算してみよう。このモジュールでは、「Qiskitパターン」として知られる量子コンピューティングのフレームワークを使用する：

ステップ1：古典的入力を量子問題にマップする
ステップ2：量子実行のための問題の最適化
ステップ 3: IBM Quantum プリミティブを使用して実行する
ステップ4：後処理と古典的分析

私たちは通常、以下のステップを踏む。

まずは、 IBM Quantum のプリミティブを含む、必要なパッケージをいくつか読み込んでみましょう。また、利用可能な量子コンピュータの中から、稼働率が最も低いものを選定します。

初回使用時に認証情報を保存するためのコードが以下にあります。ノートブックを自分の環境に保存した後、必ずこの情報をノートブックから削除してください。そうすれば、ノートブックを共有するときにあなたの認証情報が誤って共有されることはありません。詳しいガイダンスについては、 IBM Cloud アカウントの設定および信頼できない環境でのサービスの初期化を参照してください。

# Load IBM Quantum Compute Service
from qiskit_ibm_runtime import QiskitRuntimeService

# Load the Runtime primitive and session
from qiskit_ibm_runtime import EstimatorV2 as Estimator

# Syntax for first saving your token.  Delete these lines after saving your credentials.
# QiskitRuntimeService.save_account(channel='ibm_quantum_platform',
# instance = '<YOUR_IBM_INSTANCE_CRN>', token='<YOUR-API_KEY>', overwrite=True, set_as_default=True)
# service = QiskitRuntimeService(channel='ibm_quantum_platform')

# Load saved credentials
service = QiskitRuntimeService()

# Use the least busy backend, or uncomment the loading of a specific backend like "ibm_brisbane".
backend = service.least_busy(operational=True, simulator=False, min_num_qubits=127)
# backend = service.backend("ibm_brisbane")
print(backend.name)

Output:

ibm_brisbane

下のセルは、ノートブック全体を通して、シミュレーターを使うか、実際のハードウェアを使うかを切り替えることができます。今すぐ実行することをお勧めする：

# Load the Aer simulator and generate a noise model based on the currently-selected backend.
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel

# Alternatively, load a fake backend with generic properties and define a simulator.


noise_model = NoiseModel.from_backend(backend)

# Define a simulator using Aer, and use it in Sampler.
backend_sim = AerSimulator(noise_model=noise_model)

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

VQEの計算は、特定の結合距離における水素分子のハミルトニアン（ $H_2$ ）を定義することから始めます。このハミルトニアンは、標準的な手順に従って分子系から導出・変換されたものであり、量子ビット演算子を用いて系の総エネルギーを表しています。その手順は以下の通りです：1) STO-6G 基底セット（電子軌道を近似するために用いられる特定の数学関数の集合）を採用し、 2) ジョーダン・ウィグナー写像（電子を記述するフェルミオン演算子を量子ビット演算子に変換する手法）を適用し、3) 問題を簡略化するためにハミルトニアンのパリティを用いて量子ビットの縮約を行う。

先に説明したように、計算された基底状態エネルギーは、基底セットの選択と分子形状（結合距離など）に大きく依存する。この特定の構成で、これらの変換を行った場合、結果として得られる量子ビットのハミルトニアンは単純である：

$\hat{H} = -0.2355 I + 0.2355 Z$

ここで、 $I$ は恒等演算子、 $Z$ は単一量子ビットに作用するパウリ-Z演算子を表す。係数は、この特定の結合距離において、 STO-6G の基底セットを用いて計算された積分から、適切な変換を施して導き出される。

このハミルトニアンが定義されたので、VQEを使って基底状態のエネルギーを計算することができる。我々の計算した基底状態のエネルギーを期待値と比較することは有益である。単一の孤立した水素原子（H）の場合、基底状態のエネルギーは（相対論的効果がない場合）正確には -0.5 Hartreeである。上で定義した特定の量子ビット・ハミルトニアンの正確な基底状態エネルギーを計算し、関連する既知の値と比較してみよう。

from qiskit.quantum_info import SparsePauliOp
import numpy as np

# Qubit Hamiltonian of the hydrogen atom generated by using STO-3G basis set and parity mapping
Hamiltonian = SparsePauliOp.from_list([("I", -0.2355), ("Z", 0.2355)])

# exact ground state energy of Hamiltonian

A = np.array(Hamiltonian)
eigenvalues, eigenvectors = np.linalg.eig(A)
print(
    "The exact ground state energy of the Hamiltonian is ",
    min(eigenvalues).real,
    "hartree",
)
h = min(eigenvalues.real)

Output:

The exact ground state energy of the Hamiltonian is  -0.471 hartree

次に、基底状態の波動関数（ $\Psi_\text{trial}$ ）の試行を準備するために、パラメータ化された量子回路（ansatz）が必要である。目標は、エネルギー期待値 $\langle\psi(\theta)|\hat{H}|\psi(\theta)\rangle$ を最小化するパラメータ $\theta$ を見つけることである。ansatzの選択は、回路が準備できる量子状態の集合を決定するため、非常に重要である。良い "ansatzとは、我々が研究しているハミルトニアンの真の基底状態に非常に近い状態を表現するのに十分柔軟でありながら、現在の量子コンピュータにとって複雑すぎるパラメータや深すぎる回路を必要としないものである。

ここでは、3つの異なる1量子ビット法を試して、1つの量子ビットがとりうる量子状態をよりよく "カバー "するのはどれかを見てみよう。カバレッジ」とは、アンサッツ回路がパラメーターを変化させることで生成できる量子状態の範囲のことである。

異なる1量子ビット回転ゲートの組み合わせに基づく3つの解法を用いる：

1軸回転ゲートアサッツ：このアサッツでは、1つの軸（ $R_x(\theta)$ ）だけを回転させる。ブロッホ球面上では、これは特定の円に沿ってのみ動くことに相当する。これは最も柔軟性が低く、限られた州をカバーする。
つの2軸回転ゲート：これらのゲートでは、2つの異なる軸（ $R_x(\theta_1) R_z(\theta_2)$ と $R_x(\theta_1) R_z(\theta_2) R_x(\theta_3)$ ）回りの回転が組み合わされる。これにより、1軸回転に比べてブロッホ球のより広い部分に到達することができる。

これら3つのansätzeで得られたVQEの結果を比較することで、ansatzの柔軟性と状態空間カバレッジが、簡略化したハミルトニアンの真の基底状態エネルギーを求める能力にどのような影響を与えるかを見ることができる。より柔軟なansatzは、より良い近似を見つける可能性があるが、古典的なオプティマイザにとっては難しくなるかもしれない。

from qiskit import QuantumCircuit
from qiskit.circuit import Parameter
from qiskit.quantum_info import Statevector, DensityMatrix, Pauli

theta = Parameter("θ")
phi = Parameter("φ")
lam = Parameter("λ")

ansatz1 = QuantumCircuit(1)
ansatz1.rx(theta, 0)

ansatz2 = QuantumCircuit(1)
ansatz2.rx(theta, 0)
ansatz2.rz(phi, 0)

ansatz3 = QuantumCircuit(1)
ansatz3.rx(theta, 0)
ansatz3.rz(phi, 0)
ansatz3.rx(lam, 0)

Output:

<qiskit.circuit.instructionset.InstructionSet at 0x1059def80>

では、各パラメータについて5000個の乱数を発生させ、これらの乱数パラメータを用いて3つのansätzeが生成するランダム量子状態の分布をプロットしてみよう。これらのパラメータは、球面上の異なる軸の周りの回転のように考えることができる。量子状態の分布を見るために、単一量子ビットの状態を示す 3 次元球であるブロッホ球を使用します。北極と南極は古典的な "0 "と "1 "のようなものだが、量子ビットはその中間にあることもでき、重ね合わせのような特別な量子特性を示す。まず、 3D ブロッホ球をプロットするのに必要な関数を用意し、5000個のランダム・パラメータを用意する。

import matplotlib.pyplot as plt


def plot_bloch(bloch_vectors):
    # Extract X, Y, Z coordinates for 3D projection
    X_coords = bloch_vectors[:, 0]
    Z_coords = bloch_vectors[:, 2]

    # Compute Y coordinates from X and Z to approximate the full Bloch sphere projection
    Y_coords = bloch_vectors[:, 1]

    # Create 3D plot
    fig = plt.figure(figsize=(8, 8))
    ax = fig.add_subplot(111, projection="3d")
    ax.scatter(X_coords, Y_coords, Z_coords, color="blue", alpha=0.6)

    # Labels and title
    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    ax.set_zlabel("Z")
    ax.set_title("Parameterized 1-Qubit Circuit on 3D Bloch Sphere")

    # Set axis limits and make them equal
    ax.set_xlim([-1, 1])
    ax.set_ylim([-1, 1])
    ax.set_zlim([-1, 1])

    # Ensure equal aspect ratio for all axes
    ax.set_box_aspect([1, 1, 1])  # Equal scaling for x, y, z axes

    # Show grid
    ax.grid(True)

    plt.show()


num_samples = 5000  # Number of random states
theta_vals = np.random.uniform(0, 2 * np.pi, num_samples)
phi_vals = np.random.uniform(0, 2 * np.pi, num_samples)
lam_vals = np.random.uniform(0, 2 * np.pi, num_samples)

最初のアナザッツがどのように機能するか見てみよう。

# List to store Bloch Sphere XZ coordinates
bloch_vectors = []

# Generate quantum states and extract Bloch vectors
for i in range(num_samples):
    # Create a circuit and bind parameters
    qc = ansatz1
    bound_qc = qc.assign_parameters({theta: theta_vals[i]})  # , lam: lam_vals[i]})
    state = Statevector.from_instruction(bound_qc)
    rho = DensityMatrix(state)

    X = rho.expectation_value(Pauli("X")).real
    Y = rho.expectation_value(Pauli("Y")).real
    Z = rho.expectation_value(Pauli("Z")).real
    bloch_vectors.append([X, Y, Z])  # Store X, Z components

# Convert to a numpy array for plotting
bloch_vectors = np.array(bloch_vectors)

plot_bloch(bloch_vectors)

Output:

最初のアナザッツは、ブロッホ球のリング状の分散量子状態を返すことがわかる。これは理にかなっている。というのも、アサッツには回転パラメーターを1つしか与えていないからだ。したがって、1つの軸を中心に回転した状態しか作り出すことができない。 $(0,0,1)$ 、1つの軸を中心に回転させると、必ずリングができる。次に、直交する2つの回転ゲート（ Rx と Rz ）を持つ2番目のアナザッツをチェックしてみよう。

bloch_vectors = []

# Generate quantum states and extract Bloch vectors
for i in range(num_samples):
    # Create circuit and bind parameters
    qc = ansatz2
    bound_qc = qc.assign_parameters(
        {theta: theta_vals[i], phi: phi_vals[i]}
    )  # , lam: lam_vals[i]})
    state = Statevector.from_instruction(bound_qc)
    rho = DensityMatrix(state)

    X = rho.expectation_value(Pauli("X")).real
    Y = rho.expectation_value(Pauli("Y")).real
    Z = rho.expectation_value(Pauli("Z")).real
    bloch_vectors.append([X, Y, Z])  # Store X, Z components

# Convert to numpy array for plotting
bloch_vectors = np.array(bloch_vectors)

plot_bloch(bloch_vectors)

Output:

ここで、私たちの2番目のアナザッツがブロッホ球のより広い部分をカバーしていることがわかるが、点は極付近に集中し、赤道付近ではより広がっていることに注意されたい。さて、最後のansatzを確認しよう。

bloch_vectors = []

# Generate quantum states and extract Bloch vectors
for i in range(num_samples):
    # Create circuit and bind parameters
    qc = ansatz3
    bound_qc = qc.assign_parameters(
        {theta: theta_vals[i], phi: phi_vals[i], lam: lam_vals[i]}
    )
    state = Statevector.from_instruction(bound_qc)
    rho = DensityMatrix(state)

    X = rho.expectation_value(Pauli("X")).real
    Y = rho.expectation_value(Pauli("Y")).real
    Z = rho.expectation_value(Pauli("Z")).real
    bloch_vectors.append([X, Y, Z])  # Store X, Z components

# Convert to numpy array for plotting
bloch_vectors = np.array(bloch_vectors)

plot_bloch(bloch_vectors)

Output:

ここでは、前回のアナザッツによって生成された、より均等に分布した量子状態を見ることができる。

前述のように、求める基底状態についての知識を得て、その基底状態に近い状態をプローブするのに適したアナザッツを使うのがベストである。例えば、基底状態が極の近くにあるとわかっていれば、ansatz 2を選ぶかもしれない。簡単のため、ブロッホ球全体を一様に探索するansatz 3にこだわる。

さて、アサッツを選択したので、回路を描いてみよう。

# Pre-defined ansatz circuit and operator class for Hamiltonian

ansatz = ansatz3

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

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

Output:

This circuit has  3 parameters

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

実際の量子コンピューターで計算を行う場合、量子回路のロジックだけを気にするわけではない。また、その量子コンピュータでどのような演算が可能か、量子コンピュータのどこに量子ビットがあるのかといったことも気になる。隣同士ですか？距離は離れていますか？従って、次のステップは、私たちが使う量子コンピューターにとって自然なゲートを使い、量子ビットのレイアウトを考慮して回路を書き直すことだ。これは transpilation 。このプロセスの後、単純なansatzが異なるゲートのセットに変換され、抽象的な量子ビットが実際の量子コンピューター上の物理的な量子ビットにマッピングされるのを見ることができる。

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

config = backend.configuration()

print("Backend: {config.backend_name}")
print("Native gates: ", config.supported_instructions, ",")


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:

Backend: {config.backend_name}
Native gates:  ['ecr', 'id', 'delay', 'measure', 'reset', 'rz', 'sx', 'x'] ,

我々のアンサッツの rx, rz ゲートが、我々のバックエンドのネイティブゲートである rz, sx ゲートのシリーズに変換されたのがわかるだろう。また、 q0 、5番目の物理量子ビットにマッピングされていることがわかる。また、以下のコードのように、これらの変更に従ってハミルトニアンをマッピングする必要がある：

Hamiltonian_isa = Hamiltonian.apply_layout(layout=ansatz_isa.layout)

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

さて、いよいよ本物のQPUでVQEを実行する時が来た。そのためにはまず、最適化プロセスのためのコスト関数が必要である。このコスト関数は、量子状態を持つハミルトニアンの期待値を評価するもので、ansatzによって生成される。心配無用です! すべてを自分でコーディングする必要はない。そのための関数を用意しましたので、以下のセルを実行するだけです。

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

最後に、アサッツとその最適化プロセスのための初期パラメータを準備する。単純にゼロやランダムな値を使うこともできる。以下に初期パラメータを選択したが、0から $2\pi$ まで一様に、ランダムにパラメータをサンプリングするために、セル内の行を自由にコメントまたはコメント解除してください。

# x0 = np.random.uniform(0, 2*pi, 3)
x0 = [1, 1, 0]

# QPU Est. 2min for ibm_brisbane

from scipy.optimize import minimize
from qiskit_ibm_runtime import Batch

batch = Batch(backend=backend)

cost_history_dict = {
    "prev_vector": None,
    "iters": 0,
    "cost_history": [],
}
estimator = Estimator(mode=batch)
estimator.options.default_shots = 10000

res = minimize(
    cost_func,
    x0,
    args=(ansatz_isa, Hamiltonian_isa, estimator),
    method="cobyla",
    options={"maxiter": 10, "tol": 0.01},
)

batch.close()

Output:

Iters. done: 1 [Current cost: -0.3361517318448143]
Iters. done: 2 [Current cost: -0.4682546422099432]
Iters. done: 3 [Current cost: -0.38985802144149584]
Iters. done: 4 [Current cost: -0.38319217316749354]
Iters. done: 5 [Current cost: -0.4628720756579032]
Iters. done: 6 [Current cost: -0.4683301936226905]
Iters. done: 7 [Current cost: -0.45480498699294747]
Iters. done: 8 [Current cost: -0.4690533242050814]
Iters. done: 9 [Current cost: -0.465867415110354]
Iters. done: 10 [Current cost: -0.4606882723137227]

h_vqe = res.fun
print("The reference ground state energy is ", min(eigenvalues))
print("The computed ground state energy is ", h_vqe)

Output:

The reference ground state energy is  (-0.471+0j)
The computed ground state energy is  -0.4690533242050814

おめでとうございます! あなたは今、初めての量子化学実験を成功裏に終えたところです。ハミルトニアンの正確な基底状態エネルギーと我々のものとの間に違いが見られるが、我々はデフォルトのエラー緩和技術（読み出しエラーを修正する）を使用したため、その差はわずかである。これは非常に良いスタートだ！

注：を使用してエラー軽減のレベル resilience_levelを設定すると、より良い結果が得られます。デフォルト値は 1 です。これより大きい値に設定すると、QPU の処理時間は長くなりますが、より良い結果が得られる可能性があります。

ステップ4: 後処理

古典的なオプティマイザーがどのように機能したかを見てみよう。下のセルを動かして収束パターンを見てみよう。

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

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

Output:

私たちはかなり良い初期値から始めたので、わずか10ステップで良い最終値を得ることができた。大小のピークが見えるが、これはCOBYLAオプティマイザーの典型的な特徴で、まるで景色が見えないかのように空間を探索し、測定のたびにステップサイズを調整する。

理解度チェック

どうお考えですか？より理論値に近い結果、あるいはハミルトニアンの正確な基底状態エネルギーに近い結果を得るために、上記のプロセスのどの部分に改善の余地があるか？そのために考慮すべきことは？

まず考えなければならないのは、分子のハミルトニアンを計算する際に使われる塩基のセットの変化である。先に述べたように、よく知られているように、H原子の基底状態エネルギーは -0.5 Hartreeであり、我々が選んだ STO-6G の基底ではこの値を正確に導き出すには不十分である。
より複雑な基底を選択すると、ハミルトニアンが使用する量子ビットの数が増える。したがって、化学の問題では、より複雑で適切なアサッツを選択する必要がある。
次に最適化されるべきは、QPUにおけるノイズの管理である。より高度なエラー軽減技術はより良い結果をもたらすが、使用には時間がかかるかもしれない。また、 shot_number 。
最後に、異なるオプティマイザを試すことで、より優れた収束性能を達成することもできる。

VQEを用いて水素分子の基底状態エネルギーを計算する

ここまで、 $H$ 原子を使ったVQEの全体的なプロセスを見てきたが、ここからは、 $H_2$ 分子の基底状態エネルギーをより迅速に計算してみよう。

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

ここでは、 STO-6G の基底とJordan-Wigner変換を用いた1量子ビットのハミルトニアンも提供する。ハミルトニアンの対称性を利用することで、量子ビットの削減が可能である。なお、ここでは2つの水素原子間の原子間距離を 0.735 $\mathring A$ とした。

水素原子1個の計算( $H$ )とは異なり、水素分子の基底状態( $H_2$ )を計算するためには、電子軌道に関連するエネルギーに加えて、2個の水素原子の原子核間に働く反発力も考慮しなければならない。このステップでは、この値を定数として与え、チェックイン問題で実際にこの値を計算する。 $\hat{H} = -1.04886 I + -0.79674 Z + 0.18122 X$

h2_hamiltonian = SparsePauliOp.from_list(
    [("I", -1.04886087), ("Z", -0.7967368), ("X", 0.18121804)]
)

# exact ground state energy of hamiltonian
nuclear_repulsion = 0.71997
A = np.array(h2_hamiltonian)
eigenvalues, eigenvectors = np.linalg.eig(A)
print("Electronic ground state energy (Hartree): ", min(eigenvalues).real)
print("Nuclear repulsion energy (Hartree): ", nuclear_repulsion)
print(
    "Total ground state energy (Hartree): ", min(eigenvalues).real + nuclear_repulsion
)
h2 = min(eigenvalues).real + nuclear_repulsion

Output:

Electronic ground state energy (Hartree):  -1.8659468547627318
Nuclear repulsion energy (Hartree):  0.71997
Total ground state energy (Hartree):  -1.1459768547627318

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

以前のVQEとハミルトニアンで使用された量子ビット数は、実行に使用されるバックエンドと同じであるため、既存のansatzとその最適化された形式を使用する。

h2_hamiltonian_isa = h2_hamiltonian.apply_layout(layout=ansatz_isa.layout)

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

さて、次は実際のQPUで計算する番だ。ほとんどすべてが同じだが、ハミルトニアンに適合させるために適切な初期点を使用する。また、反復計算部分では、QPUのansatzのハミルトニアンの期待値を計算するために使用される Estimator の設定の一部が、以前の計算とは若干異なる設定になります。この変更については、チェックインの質問でさらに説明する。

x0 = [2, 0, 0]

# QPU time 4min for ibm_brisbane
batch = Batch(backend=backend)

cost_history_dict = {
    "prev_vector": None,
    "iters": 0,
    "cost_history": [],
}
estimator = Estimator(mode=batch)
estimator.options.default_shots = 10000

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

batch.close()

Output:

Iters. done: 1 [Current cost: -0.710621837568328]
Iters. done: 2 [Current cost: -0.2603208441168329]
Iters. done: 3 [Current cost: -0.25548711201326424]
Iters. done: 4 [Current cost: -0.581129450619904]
Iters. done: 5 [Current cost: -1.722920997605439]
Iters. done: 6 [Current cost: -1.6633324849371915]
Iters. done: 7 [Current cost: -1.8066989598929164]
Iters. done: 8 [Current cost: -1.8051093803839542]
Iters. done: 9 [Current cost: -1.802692217571555]
Iters. done: 10 [Current cost: -1.8233585485263144]
Iters. done: 11 [Current cost: -1.6904116652617205]
Iters. done: 12 [Current cost: -1.8245120321245392]
Iters. done: 13 [Current cost: -1.6837021361383608]
Iters. done: 14 [Current cost: -1.8166632606115467]
Iters. done: 15 [Current cost: -1.863446212658907]

h2_vqe = res.fun + nuclear_repulsion
print(
    "The reference ground state energy is ", min(eigenvalues).real + nuclear_repulsion
)
print("The computed ground state energy is ", h2_vqe)

Output:

The reference ground state energy is  -1.1459768547627318
The computed ground state energy is  -1.143476212658907

VQEは理論的には真の基底状態エネルギーの上限を与えるにもかかわらず、実際の量子ハードウェアやノイズの多いシミュレーション量子ハードウェアでの実用的な実装、およびハミルトニアンの準備（基底セットや量子ビットの削減など）で行われた近似により、誤差が生じることがあり、その結果、測定されたエネルギーが正確な理論値や特定の数値参照値よりもわずかに低くなることがある。多少の誤差はあるものの、特にステップ数が少ないことを考えれば、結果は満足のいくものと思われる。さて、オプティマイザーの働きを見て、このVQE計算を終えよう。

ステップ4: 後処理

fig, ax = plt.subplots()
x = np.linspace(0, 5, 15)

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

Output:

理解度チェック

$H_2$ 、定数値（ 0.71997 Hartree）として入れた分子の核反発エネルギーを計算してみよう。

クーロンの法則と原子単位を用いて、値がとなることを確認してください Hartree。

どちらの水素原子核もプラスに帯電しているため、静電気力によって反発し合う。この斥力はクーロンの法則で説明される：
$E_{repulsive} = \frac{e^2}{4\pi\epsilon_0R}$
ここで、 $e$ は陽子の電荷、 $\epsilon_0$ は真空の誘電率、 $R$ は2つの原子核間の距離で、単位はメートルまたはボーア半径、単位はジュール(J)。
このエネルギーをハートリー単位で計算するには、上記の式を原子単位（AU）システムに変換する必要がある。 AUでは、 $e^2 = 1$ 、 $4\pi\epsilon_0=1$ 、ボーア半径( $a_0$ )が1となり、AUの基本的な長さスケールになる。このように単純化すると、クーロンの法則は次のようになる：
$E_{repulsion} = \frac{1}{R}$
ここで、 $R$ はボーア半径( $a_0$ )で測定しなければならない。
与えられた核分離（ $\r{A}$ ）を $a_0$ に変換するには、この変換関係が必要である：
$1\r{A} = 1.88973 a_0$
だから $0.735\r{A}$ は $0.735 * 1.88973 = 1.38895 a_0$ になる。
したがって、与えられた $H_2$ の核反発エネルギーは次のようになる
$E_{repulsion} = \frac{1}{R} = \frac{1}{1.38895} = 0.71997 Hartree$

$H + H = H_2$ の反応エネルギーを計算する

では、手に入れたものを使ってみよう！あなたはVQE（変分量子固有値ソルバー）を使って、 $H$ 原子と $H_2$ 分子の基底状態エネルギーを計算した。あとは計算値を使って、 $H+H=H_2$ プロセスの反応エネルギーを求めるだけだ。

反応エネルギーとは、物質が反応して新しい物質を形成するときに起こるエネルギー変化のことである。ブロックを積み上げるように）エネルギーを注ぎ込む必要があるときもあれば、（坂道を転がるボールのように）エネルギーが放出されるときもある。化学では、反応はエネルギーを吸収するか（吸熱）、エネルギーを放出する（発熱）。

$H+H = H_2$ プロセスの反応エネルギーは以下の式で計算できる：

$E_{reaction} = E_{H_2} - (E_H + E_H)$

下のセルを実行して、これを視覚的に見てみよう。ここでは、各ハミルトニアンの厳密な基底状態の値を用い、厳密解とVQEの結果の反応エネルギーを比較する。

# Theoretical values
E_H_theo = h.real
E_H2_theo = h2

# Experimental values
E_H_exp = h_vqe
E_H2_exp = h2_vqe

# Calculate reaction energies
E_reaction_theo = E_H2_theo - (2 * E_H_theo)
E_reaction_exp = E_H2_exp - (2 * E_H_exp)

# Set up the plot
fig, ax = plt.subplots(figsize=(8, 6))
ax.set_xlim(0, 3)
ax.set_ylim(-1.16, -0.93)  # Adjust y-axis range to highlight differences
ax.set_xticks([])
ax.set_ylabel("Energy (Hartree)")
ax.set_title("H + H → H₂ Reaction Energy Diagram")

# Plot theoretical energy levels
ax.hlines(
    y=2 * E_H_theo, xmin=0.5, xmax=1.3, linewidth=2, color="r", label="2H (Exact)"
)
ax.hlines(y=E_H2_theo, xmin=1.3, xmax=2, linewidth=2, color="b", label="H₂ (Exact)")

# Plot experimental energy levels
ax.hlines(
    y=2 * E_H_exp,
    xmin=0.5,
    xmax=1.5,
    linewidth=2,
    color="r",
    linestyle="dashed",
    label="2H (VQE)",
)
ax.hlines(
    y=E_H2_exp,
    xmin=1.5,
    xmax=2.5,
    linewidth=2,
    color="b",
    linestyle="dashed",
    label="H₂ (VQE)",
)

# Add labels
ax.text(
    1,
    2 * E_H_theo,
    f"2H: {2*E_H_theo:.4f}",
    verticalalignment="top",
    horizontalalignment="left",
)
ax.text(
    2,
    E_H2_theo,
    f"H₂: {E_H2_theo:.4f}",
    verticalalignment="top",
    horizontalalignment="left",
)
ax.text(
    1,
    2 * E_H_exp,
    f"2H_VQE: {2*E_H_exp:.4f}",
    verticalalignment="bottom",
    horizontalalignment="right",
)
ax.text(
    2,
    E_H2_exp,
    f"H₂_VQE: {E_H2_exp:.4f}",
    verticalalignment="bottom",
    horizontalalignment="right",
)

# Add arrows for reaction energy with ΔE label in the middle
mid_y_theo = (2 * E_H_theo + E_H2_theo) / 2
mid_y_exp = (2 * E_H_exp + E_H2_exp) / 2
ax.annotate(
    "",
    xy=(1.3, E_H2_theo),
    xytext=(1.3, 2 * E_H_theo),
    arrowprops=dict(arrowstyle="<->", color="g"),
)
ax.text(
    1.35, mid_y_theo, f"ΔE: {E_reaction_theo:.4f}", color="g", verticalalignment="top"
)

ax.annotate(
    "",
    xy=(1.5, E_H2_exp),
    xytext=(1.5, 2 * E_H_exp),
    arrowprops=dict(arrowstyle="<->", color="g", linestyle="dashed"),
)
ax.text(
    1.55,
    mid_y_exp,
    f"ΔE_VQE: {E_reaction_exp:.4f}",
    color="g",
    verticalalignment="center",
)

# Add legend
ax.legend()

plt.show()

Output:

図に示すように、多少の誤差はあるものの、ハミルトニアンの正確な基底状態エネルギーと、VQEの結果を用いて計算された反応エネルギーは、 -0.2 Hartreeに近く、類似している。

ここで注意しなければならないのは、このプロセスの反応エネルギーは負の値を持つということである。これは、プロセスを通じてエネルギーが放出され、結果として得られる分子は、2つの単原子よりもエネルギーが低いことを意味する。

おわりに

これまで学んだことをまとめてみよう。

まず、量子化学の問題を解くために必要な2つの重要な近似技術、すなわち変分原理と基底セットの選択について調べた。われわれは変分原理を手で探り、単純調和振動子の基底状態エネルギーを計算した。

次に、量子系の基底状態エネルギーを計算するアルゴリズムとして広く使われているVQEについて調べた。原子状水素( $H$ )と水素分子( $H_2$ )の基底状態エネルギーを計算するコードを実行した。特に、系の適切な分子ハミルトニアンを求め、それを量子コンピュータで実行可能な形に変換する必要があることを学んだ。また、VQE内で試行量子状態を準備するためには、パラメータ化された量子回路であるアンサッツが必要であることを確認し、適切なアンサッツ回路構造を選択することの重要性を議論した。私たちはまた、VQEが古典コンピュータを使った反復最適化プロセスに依存していることを学び、量子回路を導いて最低エネルギー状態を見つけ、そのプロセスがどのように収束するかを見た。

最後に、VQEによって得られた $H$ と $H_2$ の計算基底状態エネルギーを用いて、 $H + H \rightarrow H_2$ プロセスの反応エネルギーを計算した。

VQEは、近い将来強力な量子アルゴリズムとなるが、その限界を認識しておくことが重要である。 VQEの性能はアサッツの選択に大きく依存し、真の基底状態を正確に表現できる効率的に準備可能なアサッツを見つけることは、より大きく複雑な分子では困難となる。さらに、現在の量子ハードウェアはノイズの影響を受けやすく、特に回路が深い場合や量子ビット数が多い場合には、VQE結果の精度に影響を与える可能性がある。このような課題にもかかわらず、VQEは基礎的なアルゴリズムとして機能しており、現在進行中の研究では、より洗練された変分法とエラー緩和技術を探求し、近い将来の量子コンピュータで量子化学の可能性の限界を押し広げようとしている。例えば、サンプルベースの量子対角化（SQD）のようなアルゴリズムが開発されている。これは、エネルギー推定を改善し、VQEが直面するいくつかの制限、特に測定効率とノイズロバスト性に対処するために、部分空間における古典的対角化と組み合わせた量子回路から得られたサンプルを活用するものである。

レビューと質問

重要な概念：

変分量子アルゴリズムは、古典コンピュータと量子コンピュータが協力して問題を解決するコンピューティングパラダイムである。
VQEでは、システムのハミルトニアンから始め、それを量子コンピュータ上で実行するための量子ビットにマッピングする。パラメータ化された量子回路（アナザッツ）を選択し、アナザッツのパラメータを変化させながら、最低エネルギー値に達するまで測定を繰り返す。パラメータ空間の探索は、古典的なオプティマイザを用いて行われる。良い結果を得るためには、良いアナザッツと適切なオプティマイザを選択する必要がある。
反応エネルギーとは、化学反応における総エネルギー変化のことで、反応物と生成物のエネルギーの差によって決まる。

True/False

変分原理は、任意の試行波動関数のエネルギーの期待値は、常に真の基底状態のエネルギー以上であることを示す。
基底セットとは、量子波動関数を近似するための関数の集まりである。
VQEは、与えられたハミルトニアンに対するシュレーディンガー方程式を正確に解くために使用される量子アルゴリズムである。
VQEでは、パラメータ化された量子回路（アンサッツ）を用いて試行波動関数を準備する。
VQEにおけるオプティマイザの選択（例えば、COBYLA、SPSA、ADAM）は、結果の品質に影響を与えない。
VQEにおけるハミルトニアンの期待値を直接計算するために、Qiskitの Estimator 。

選択問題：

VQEにおけるハミルトニアンの目的は何ですか？

A) ランダムな量子状態を生成する
B) 量子状態のエネルギーを決定する
C)量子回路の最適化
D) エンタングルメントを起こす

VQEアルゴリズムの主な目的は何ですか？

A) ハミルトニアンの基底状態エネルギーを求めるには
B) 量子ビット間にエンタングルメントを作る
C)グローバーの検索を実行する
D) RSA暗号を解読する

このノートでは、ansatzを比較するために何個の量子状態が生成されますか？

A) 100
B) 1000
C)5000
D) 10,000

なぜVQEには古典的なオプティマイザが必要なのか？

A) 量子測定
B) エネルギーが最小になるようにansatzパラメータを更新する
C)量子ビットを絡ませる
D) 量子ランダム性の生成

なぜアサッツはパラメータ化されるように設計されているのですか？

A) 量子状態の準備を可能にする
B) 量子状態の広い空間を探索できるようにする
C)回路の複雑さを軽減する
D) 固有値を直接測定する

良いアンサッツの選択に関する次の記述のうち、最も正しいものはどれか？

A) アナザッツはブロッホ球に均等に分布した状態を生成しなければならない。
B) アナサッツは、基底状態に近い状態を生成できることを確認するために、あなたのシステムに合わせるべきである。
C)アナザッツは、変分パラメータを使ってランダムな状態を生成する。
D) より良いアナザッツは常に、より多くの変分パラメーターを持つ。

(オプション) 付録：仮説の複雑さによるオプティマイザのオーバーヘッド

VQEはいくつかのよく知られた課題ref [6に]直面しており、以下は上記で学んだことと関連している。

アンサッツ選択の課題

適切な変分アサッツを選択することには、固有の課題がある。 UCCSDのような）化学にインスパイアされた解析手法は物理的な正確さを提供するが、深い回路を必要とする。一方、ハードウェア的に効率的な解析手法は回路が浅いが、物理的な解釈可能性に欠ける可能性がある。また、多くのansätzeは、精度の向上にはほとんど寄与しないが、最適化の難易度を著しく高める過剰な変分パラメータを導入している。

最適化の難しさ

VQEの最適化ランドスケープには、勾配が指数関数的に消滅する領域（不毛のプラトー）が存在することがあり、古典的な最適化手法では変分パラメータを効率的に更新することが難しくなります。このため、研究者たちは、勾配ベースと勾配なしという異なるタイプのオプティマイザを使おうとしてきたが、どちらも課題を抱えている。勾配ベースのオプティマイザは不毛のプラトーに悩まされ、勾配なしの方法は大量の関数評価を必要とする。

オプティマイザのオーバーヘッド

さらによく知られた課題として、問題の規模に関係するオプティマイザのオーバーヘッドがある。 VQEに必要な量子回路は、問題サイズが大きくなるにつれて深さと複雑さが増し、その結果、最適化すべきパラメータの数も一般的に増加する。最適化プロセスは、パラメータの数が増えるにつれて難解になり、収束が遅くなり、最適解を見つけるのが困難になる。

ここでは、 $H_2$ 、2つの異なるタイプのansätzeを持つ分子に対してVQEを使用することで、これらの課題を見てみよう。

(注：この作業にはQPUの時間がかかるので、時間がない場合はシミュレーターを使ってください)

from qiskit.circuit import ParameterVector

num_iter = 4
alpha = ParameterVector("alpha", 3)
beta = ParameterVector("beta", 3 * num_iter)

# step1: Map problem to quantum circuits and operators
hamiltonian = SparsePauliOp.from_list(
    [("I", -1.04886087), ("Z", -0.7967368), ("X", 0.18121804)]
)

ansatz_1 = ansatz3
ansatz_2 = QuantumCircuit(1)
for i in range(num_iter):
    ansatz_2.rx(beta[i * 3 + 0], 0)
    ansatz_2.rz(beta[i * 3 + 1], 0)
    ansatz_2.rx(beta[i * 3 + 2], 0)

ansatz_1.draw("mpl")

Output:

ansatz_2.draw("mpl")

Output:

# Step 2: Optimize for target hardware

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

ansatz_isa_1 = pm.run(ansatz_1)
ansatz_isa_2 = pm.run(ansatz_2)
hamiltonian_isa_1 = hamiltonian.apply_layout(layout=ansatz_isa_1.layout)
hamiltonian_isa_2 = hamiltonian.apply_layout(layout=ansatz_isa_2.layout)

では、初期点がすべて1であるVQEを最大20ステップで実行し、両者の収束性を比較してみよう。

# QPU time 3m 40s for ibm_brisbane
# Step 3: Execute on target hardware

from scipy.optimize import minimize

x0 = np.ones(ansatz_1.num_parameters)

batch = Batch(backend=backend)


cost_history_dict = {
    "prev_vector": None,
    "iters": 0,
    "cost_history": [],
}
estimator = Estimator(mode=batch)
estimator.options.default_shots = 2048

res = minimize(
    cost_func,
    x0,
    args=(ansatz_isa_1, hamiltonian_isa_1, estimator),
    method="cobyla",
    options={"maxiter": 20},
)

batch.close()

Output:

Iters. done: 1 [Current cost: -0.8782202668652658]
Iters. done: 2 [Current cost: -0.43473160695469165]
Iters. done: 3 [Current cost: -0.4076372093159749]
Iters. done: 4 [Current cost: -1.3587839859772106]
Iters. done: 5 [Current cost: -1.774529906754082]
Iters. done: 6 [Current cost: -1.541934983115727]
Iters. done: 7 [Current cost: -1.2732403113465345]
Iters. done: 8 [Current cost: -1.820842221085785]
Iters. done: 9 [Current cost: -1.8065762857059005]
Iters. done: 10 [Current cost: -1.8126394095981146]
Iters. done: 11 [Current cost: -1.8205831886180421]
Iters. done: 12 [Current cost: -1.8086715778994924]
Iters. done: 13 [Current cost: -1.8307676638629322]
Iters. done: 14 [Current cost: -1.8177328827556327]
Iters. done: 15 [Current cost: -1.8179426218088064]
Iters. done: 16 [Current cost: -1.8109239667991088]
Iters. done: 17 [Current cost: -1.824271872489647]
Iters. done: 18 [Current cost: -1.813167587671394]
Iters. done: 19 [Current cost: -1.824647343397313]
Iters. done: 20 [Current cost: -1.8219785311686143]

# Save Cost_history as a new list
ansatz_1_history = cost_history_dict["cost_history"]

# QPU time 3m 40s for ibm_brisbane

x0 = np.ones(ansatz_2.num_parameters)

batch = Batch(backend=backend)


cost_history_dict = {
    "prev_vector": None,
    "iters": 0,
    "cost_history": [],
}
estimator = Estimator(mode=batch)
estimator.options.default_shots = 2048

res = minimize(
    cost_func,
    x0,
    args=(ansatz_isa_2, hamiltonian_isa_2, estimator),
    method="cobyla",
    options={"maxiter": 20},
)

batch.close()

Output:

Iters. done: 1 [Current cost: -0.738191173881188]
Iters. done: 2 [Current cost: -0.42636037194506304]
Iters. done: 3 [Current cost: -1.3503788613797374]
Iters. done: 4 [Current cost: -0.9109204349776897]
Iters. done: 5 [Current cost: -0.9060873157510835]
Iters. done: 6 [Current cost: -0.7735065414083984]
Iters. done: 7 [Current cost: -1.586889197437709]
Iters. done: 8 [Current cost: -1.659215191584943]
Iters. done: 9 [Current cost: -1.245445981794618]
Iters. done: 10 [Current cost: -1.1608385766138023]
Iters. done: 11 [Current cost: -1.1551733876027737]
Iters. done: 12 [Current cost: -1.8143337768286332]
Iters. done: 13 [Current cost: -1.2510951563756598]
Iters. done: 14 [Current cost: -1.6918311531865413]
Iters. done: 15 [Current cost: -1.8163783305531838]
Iters. done: 16 [Current cost: -1.8434877732947152]
Iters. done: 17 [Current cost: -1.8461898233304472]
Iters. done: 18 [Current cost: -1.0346471214915485]
Iters. done: 19 [Current cost: -1.8322518854150687]
Iters. done: 20 [Current cost: -1.717144678705999]

ansatz_2_history = cost_history_dict["cost_history"]

fig, ax = plt.subplots()

# Define the constant function)
ax.plot(
    range(cost_history_dict["iters"]),
    ansatz_1_history,
    label="Ansatz with 3 parameters",
)
ax.plot(
    range(cost_history_dict["iters"]),
    ansatz_2_history,
    label="Ansatz with 12 parameters",
)
ax.set_xlabel("Iterations")
ax.set_ylabel("Cost (Hartree)")
plt.legend()
plt.draw()

Output:

上のグラフは、より多くの変数を持つansatzの最適化プロセスが、安定した収束に至るまでに多くの時間を要することを明確に示している。

単純な単一量子ビット回路と単純なansatzに頼るのではなく、より大きな量子回路とより複雑な構造のansätzeが必要になると、最適化の複雑さが増す。これは、オプティマイザのオーバーヘッドという、VQEにおけるよく知られた課題を浮き彫りにしている。

研究者らは化学問題に量子コンピュータを活用できる様々な先進的手法の開発を続けている。さまざまな教材をで閲覧できます IBM Quantum Learning。

参照

[参考文献1 ] Richard P. Feynman, Simulating Physics with Computers, International Journal of Theoretical Physics, 1982.
[ref 2] Marov, M.Y. (2015). 宇宙の構造で：現代宇宙物理学の基礎. シュプリンガー、ニューヨーク、ニューヨーク
[参考文献3 ] How to solve difficult chemical engineering problems with quantum computing, IBM Research Blog, 2023.
[文献4 ] Y. 曹、J. ロメロ、A. Aspuru-Guzik, "Potential of quantum computing for drug discovery," in IBM Journal of Research and Development, vol.62, no.6, pp.6:1-6:20, 1 Nov.-Dec.2018
[参考文献5 ] 分子構造計算の現状, REv. モジュールフィジーク。 32, 170, 1960
[参考文献6 ] Fedorov, D.A、 Peng, B、 Govind, N. et al. VQE法：簡単な調査と最近の進展。メーテル理論 6, 2 (2022)

このページは役に立ちましたか？

バグや誤字の報告、またはコンテンツの要求はGitHubで行ってください。

変分量子固有値ソルバー (VQE)

概要

VQE - 固有値問題のための変分量子アルゴリズム

化学における近似手法 - 変分原理と基底関数

理解度チェック

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

理解度チェック

(量子) 観測可能量：分子ハミルトニアン（分子のエネルギー）

(量子) アントザス：試行波動関数（量子回路を用いて自明な量子状態を構築する方法）

(古典)オプティマイザー：量子回路の微調整

理解度チェック

VQEを用いて水素原子の基底状態エネルギーを計算する

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

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

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

ステップ4: 後処理

理解度チェック

VQEを用いて水素分子の基底状態エネルギーを計算する

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

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

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

ステップ4: 後処理

理解度チェック

H+H=H2H + H = H_2H+H=H2​ の反応エネルギーを計算する

レビューと質問

重要な概念：

True/False

選択問題：

(オプション) 付録：仮説の複雑さによるオプティマイザのオーバーヘッド

参照

$H + H = H_2$ の反応エネルギーを計算する