diff --git a/source/common/references.bib b/source/common/references.bib
index 48287506..d70d1296 100644
--- a/source/common/references.bib
+++ b/source/common/references.bib
@@ -235,7 +235,7 @@ @incollection{cobyla
pages = {51--67}
}
-@article{PhysRevA.98.032309,
+@article{quantum_circuit_learning,
title = {Quantum circuit learning},
author = {Mitarai, K. and Negoro, M. and Kitagawa, M. and Fujii, K.},
journal = {Phys. Rev. A},
diff --git a/source/data/vqe_result_1q_estimator.pkl b/source/data/vqe_result_1q_estimator.pkl
new file mode 100644
index 00000000..603e61db
Binary files /dev/null and b/source/data/vqe_result_1q_estimator.pkl differ
diff --git a/source/data/vqe_result_2q.pkl b/source/data/vqe_result_2q.pkl
new file mode 100644
index 00000000..d9566de2
Binary files /dev/null and b/source/data/vqe_result_2q.pkl differ
diff --git a/source/data/vqe_results.pkl b/source/data/vqe_results.pkl
index 49992258..7569a148 100644
Binary files a/source/data/vqe_results.pkl and b/source/data/vqe_results.pkl differ
diff --git a/source/data/vqe_results_1q.pkl b/source/data/vqe_results_1q.pkl
new file mode 100644
index 00000000..a54879ea
Binary files /dev/null and b/source/data/vqe_results_1q.pkl differ
diff --git a/source/figs/vqe_u3.png b/source/figs/vqe_u3.png
index e79b9ec4..e983c21c 100644
Binary files a/source/figs/vqe_u3.png and b/source/figs/vqe_u3.png differ
diff --git a/source/ja/python_intro.md b/source/ja/python_intro.md
index 9ec4b6d2..148dfbc1 100644
--- a/source/ja/python_intro.md
+++ b/source/ja/python_intro.md
@@ -896,6 +896,19 @@ for key, value in d.items():
d.get('that', 39) # 39 (default value)
```
+リスト、タプル、辞書などを生成するときに時々使う構文を紹介します。
+
+```{code-cell} ipython3
+# list comprehension: 一行の中でforループを回す
+list(f'a{i}' for i in range(3)) # ['a0', 'a1', 'a2']
+
+# dictは（キー, 値）のタプルの配列からも生成できる
+dict((f'v{i}', i) for i in range(3)) # {'v0': 0, 'v1': 1, 'v2': 2}
+
+# dictのタプル配列からの生成とzipの組み合わせ
+dict(zip(['a', 'b', 'c'], [3, 4, 5])) # {'a': 3, 'b': 4, 'c': 5}
+```
+
### NumPy
Pythonで数値計算をする際に今や欠かせない存在になっているのが、NumPy（「なむぱい」もしくは「なむぴー」）ライブラリです。NumPyは基本的には数値の（多次元）配列に対して効率的に数値計算をすることを目的に書かれていますが、サポートされている計算オペレーションの多様さも支持を広げる要因になっています。
diff --git a/source/ja/vqe.md b/source/ja/vqe.md
index 4ada44d3..fbee89c4 100644
--- a/source/ja/vqe.md
+++ b/source/ja/vqe.md
@@ -28,8 +28,6 @@ language_info:
この実習では、変分法の基本的な考え方と、その方法に基づいた変分量子アルゴリズムと呼ばれる量子計算の手法を学びます。特に、量子計算と古典計算を組み合わせた「**量子・古典ハイブリッドアルゴリズム**」としての変分量子アルゴリズムに着目します。この手法を用いて、近似的な固有値計算を可能にする**変分量子固有値ソルバー法**と呼ばれる方法へ拡張していきます。
-この教材は、Qiskit textbookの["Simulating Molecules using VQE"](https://qiskit.org/textbook/ch-applications/vqe-molecules.html)を参考にしています。
-
```{contents} 目次
---
local: true
@@ -52,7 +50,7 @@ $$
\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)}$が最小になるように古典計算を繰り返しながら求めていくことになります。
+を満たす、できるだけ小さい$\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)}$が最小になるように古典計算を繰り返しながら求めていくことになります。
+++
@@ -102,6 +100,8 @@ $$
まず、VQEの元になっている**変分量子アルゴリズム**（*Variational Quantum Algorithm*, VQA）と呼ばれる手法について見ていきます。
++++
+
### 変分量子回路
量子コンピューター上で変分法を実装するには、*Ansatz*を更新する仕組みが必要です。量子状態の更新には量子ゲートが使えることを、私たちは知っていますね。VQAも量子ゲートを使いますが、VQAは決まった構造を持つパラメータ化された量子回路（**変分量子回路**と呼ばれる）を使って行います。この回路は**変分フォーム**（*variational form*）と呼ばれる場合もあり、回路をひとまとめにしてユニタリー変換$U(\theta)$と書くことも多いです（$\theta$はパラメータで、複数ある場合はベクトルになります）。
@@ -109,28 +109,40 @@ $$
変分フォームの決め方ですが、解きたい問題のドメインに応じて特定の構造を持つ変分フォームを導入することがあります。そうではなく、幅広い問題への応用ができるようにドメインに依存しない形の変分フォーム（例えば$R_X$や$R_Y$などの回転ゲート）を使うこともあります。後で高エネルギー実験へのVQEの応用を課題として考えますが、そこでは$R_Y$と制御$Z$ゲートを使った変分フォームを実装します。
++++
+
### 単純な変分フォーム
-変分フォームを決める時には、2つの相反する目的に対するバランスを考える必要があります。$n$量子ビットの変分フォームは、パラメータの数を増やせば$\ket{\psi} \in \mathbb{C}^N$（$N=2^n$）の任意の状態ベクトル$\ket{\psi}$を生成できるでしょう。しかし、パラメータを最適化することを考えれば、できるだけ少ないパラメータで変分フォームを構築したいですよね。回転角をパラメータとする回転ゲートの数が増えれば、量子コンピュータで動かす場合はそれだけノイズの影響を受けやすくなります。なので、できるだけ少ないパラメータ（やゲート）で求める状態を生成できればベストでしょう。
+変分フォームを決める時には、2つの相反する目的に対するバランスを考える必要があります。$n$量子ビットの変分フォームは、パラメータの数を増やせば実自由度$2^{n+1}-2$の任意の状態ベクトル$\ket{\psi}$を生成できるでしょう。しかし、パラメータを最適化することを考えれば、できるだけ少ないパラメータで変分フォームを構築したいですよね。回転角をパラメータとする回転ゲートの数が増えれば、量子コンピュータで動かす場合はそれだけノイズの影響を受けやすくなります。なので、できるだけ少ないパラメータ（やゲート）で求める状態を生成できればベストでしょう。
-ここでは、まず$n=1$の場合を考えます。$U3$ゲートは3つのパラメータ$\theta$、$\phi$、$\lambda$を使って以下の変換を表現します:
+ここでは、まず$n=1$の場合を考えます。Qiskitの$U$ゲート（上の$U(\theta)$ノーテーションと紛らわしいですが、こちらは単一のゲート）は3つのパラメータ$\theta$、$\phi$、$\lambda$を使って以下の変換を表現します:
$$
-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}
+U(\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}
$$
-系全体にかかるグローバルな位相を除けば、3つのパラメータを適切に設定して実装すれば任意の単一量子ビット状態への変換が行えます。そういう意味で、この変分フォームは**ユニバーサル**な変換が可能で、かつ3つしかパラメータがないため効率的に最適化できるという特徴があります。ただユニバーサルに任意の状態を生成できるということは、この変分フォームが生成する状態を使ってあるハミルトニアン$H$の期待値を計算しようとした場合、その固有状態を近似する状態だけでなく、それ以外のさまざまな状態も含むということになります。つまり、VQEで最小固有値を効率的に決めるためには、こういう不要な状態を避けつつ、古典計算でいかに適切にパラメータを最適化できるかにかかっているわけです。
+変分フォームの初期状態を$\ket{0}$に取るならば、上の行列の第一列のみが状態生成に寄与し、$\theta$と$\phi$の2つのパラメータで任意の単一量子ビット状態を表現できます。そういう意味で、この変分フォームは**ユニバーサル**であると言います。ただ、ユニバーサルに任意の状態を生成できるということは、この変分フォームが生成する状態を使ってあるハミルトニアン$H$の期待値を計算しようとした場合、その固有状態を近似する状態だけでなく、それ以外のさまざまな状態も含むということになります。つまり、VQEで最小固有値を効率的に決められるかどうかは、こういう不要な状態を避けつつ、古典計算でいかに適切にパラメータを最適化できるかにかかっているわけです。
+
++++
### パラメータの最適化
パラメータ化された変分フォームを選択したら、所望の量子状態を作ることができるように、変分法に従ってパラメータを最適化する必要があります。パラメータの最適化のプロセスには様々な課題があります。例えば量子ハードウェアには様々なタイプのノイズがあるため、その状態でエネルギーを測定しても正しい答えが返ってくるという保証はありません。そのために、パラメータの最適化に使う目的関数の評価が実際のエネルギーの値からずれてしまい、正しくパラメータの更新ができない可能性があります。また、最適化の手法（**オプティマイザー**）によっては、パラメータの数に依って目的関数を評価する回数が増えることがあり、ノイズの影響をさらに受けやすくなります。つまり、アプリケーションの要求を考慮しながら、オプティマイザーの選択にも気を配る必要があります。
最も一般的な最適化手法は、エネルギーの減少が極大になるような方向に各パラメータを更新する**勾配降下法**です。各パラメータごとに勾配を計算するため、最適化すべきパラメータの数に応じて目的関数を評価する回数は増えます。また、この性質から探索空間の中で局所的な最適値を素早く発見することは可能ですが、逆に探索が局所的な最小値に留まってしまうことがあります。勾配降下法は直感的で理解しやすい最適化の手法ですが、少なくとも現在のNISQコンピュータでは精度良く実行することは難しいと考えられています（VQEを実装する時に、勾配計算を使った最適化の方法を紹介します）。
-ノイズのある量子コンピュータで目的関数を最適化する適切なオプティマイザーとして、*Simultaneous Perturbation Stochastic Approximation*（**SPSA**）{cite}`bhatnagar_optimization`があります。SPSAは２回の測定だけで目的関数の勾配を近似できるという特徴があります。勾配降下法では各パラメータを独立に変化させるのに対して、SPSAでは全てのパラメータを同時にランダムに変化させます。以上のことから、VQEを利用する場合のオプティマイザーとしてSPSAが推奨されることがあります。
+ノイズのある量子コンピュータで目的関数を最適化する適切なオプティマイザーとして、*Simultaneous Perturbation Stochastic Approximation*（**SPSA**）{cite}`bhatnagar_optimization`があります。SPSAは2回の測定だけで目的関数の勾配を近似できるという特徴があります。勾配降下法では各パラメータを独立に変化させるのに対して、SPSAでは全てのパラメータを同時にランダムに変化させます。以上のことから、VQEを利用する場合のオプティマイザーとしてSPSAが推奨されることがあります。
-ノイズがない量子コンピュータで目的関数を評価する場合（例えば状態ベクトルシミュレータで実行する場合など）は、Pythonの[SciPy](https://www.scipy.org/scipylib/index.html)パッケージで提供されているオプティマイザーなど、様々な選択肢があります。この実習では、Qiskitでサポートされているオプティマイザーの中で、特に*Constrained Optimization by Linear Approximation*（**COBYLA**）と呼ばれるオプティマイザーも使用します。COBYLAは目的関数の評価を1回しか実行しない（つまり評価の回数がパラメータの数に依存しない）ため、ノイズがない状態でかつ評価の回数を少なくしたい場合にはCOBYLAの利用が推奨されているようです。いずれにしろ、どのオプティマイザーがベストかはVQEアルゴリズムの実装形式や実行環境によって変わるため、ある程度経験によって決める必要があると考えられます。
+ノイズがない量子コンピュータで目的関数を評価する場合（例えば状態ベクトルシミュレータで実行する場合など）は、PythonのSciPyパッケージで提供されているオプティマイザーなど、様々な選択肢があります。この実習では、Qiskitでサポートされているオプティマイザーの中で、特に*Constrained Optimization by Linear Approximation*（**COBYLA**）と呼ばれるオプティマイザーも使用します。COBYLAは目的関数の評価を1回しか実行しない（つまり評価の回数がパラメータの数に依存しない）ため、ノイズがない状態でかつ評価の回数を少なくしたい場合にはCOBYLAの利用が推奨されているようです。いずれにしろ、どのオプティマイザーがベストかはVQEアルゴリズムの実装形式や実行環境によって変わるため、ある程度経験によって決める必要があると考えられます。
+
++++
### 変分フォームを使った実例
-ではここで、単一量子ビットの変分フォームを利用してパラメータ最適化の例を実行してみましょう。例として、ランダムな確率分布のベクトル$\vec{x}$（要素数は2）を入力として与えた時、出力の確率分布が$\vec{x}$に近くなるように単一量子ビットの変分フォームを決定するという問題を考えます（2つの確率分布の近さはL1距離によって定義します）。
+ではここで、$U$ゲート一つからなる単一量子ビットの変分フォームを利用してパラメータ最適化の例を実行してみましょう。ランダムに量子状態$\ket{\psi_0}$を選び、$\ket{\psi(\theta, \phi)} := U(\theta, \phi, 0)\ket{0}$でそれを近似するという問題を考えます[^actually_exact]。1量子ビットの状態は観測量$X, Y, Z$の期待値$\langle X \rangle, \langle Y \rangle, \langle Z \rangle$の値がわかれば（全体位相を除いて）完全に決まるので、$\ket{\psi(\theta, \phi)}$による$X, Y, Z$の期待値$\langle X \rangle_{\theta, \phi}, \langle Y \rangle_{\theta, \phi}, \langle Z \rangle_{\theta, \phi}$が$\ket{\psi_0}$による対応する期待値$\langle X \rangle_0, \langle Y \rangle_0, \langle Z \rangle_0$に等しくなるように$\theta, \phi$を選ぶことになります。したがって、問題は目的関数
+
+$$
+L(\theta, \phi) = [\langle X \rangle_{\theta, \phi} - \langle X \rangle_0]^2 + [\langle Y \rangle_{\theta, \phi} - \langle Y \rangle_0]^2 + [\langle Z \rangle_{\theta, \phi} - \langle Z \rangle_0]^2
+$$
+
+の最小化となります。
```{image} figs/vqe_u3.png
:alt: vqe_u3
@@ -138,7 +150,7 @@ $$
:align: center
```
-最初に、Pythonでランダムな確率分布のベクトルを作成します。
+[^actually_exact]: $U(\theta, \phi, 0)$は単一量子ビットについてユニバーサルなので、原理的には近似ではなく厳密に一致させることができます。
```{code-cell} ipython3
---
@@ -149,14 +161,18 @@ pycharm:
'
---
-# Tested with python 3.8.12, qiskit 0.34.2, numpy 1.22.2
import numpy as np
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, transpile
+from qiskit.circuit import Parameter, ParameterVector
+from qiskit.primitives import BackendEstimator
+from qiskit.quantum_info import Operator, SparsePauliOp
from qiskit.algorithms.optimizers import SPSA, COBYLA
from qiskit_aer import AerSimulator
```
+最初に、ターゲットとなる量子状態ベクトルをランダムに生成する関数と、状態ベクトルから$X, Y, Z$の期待値を計算する関数を定義します。
+
```{code-cell} ipython3
---
jupyter:
@@ -166,161 +182,444 @@ pycharm:
'
---
-nq = 1 # 量子ビットの数
+rng = np.random.default_rng(999999)
+
+# 与えられた量子ビット数（nq）に応じたランダムな状態ベクトルを作る関数
+def random_statevector(nq):
+ # 2^nqの複素数をランダムに生成
+ target_statevector = rng.random(2 ** nq) + 1.j * rng.random(2 ** nq)
+ # 正規化
+ target_statevector /= np.sqrt(np.sum(np.square(np.abs(target_statevector))))
+
+ return target_statevector
+
+# 与えられた状態ベクトルにおけるパウリ演算子（X, Y, Z）の期待値を計算する関数
+# opは"X"や"XZZ"など
+def pauli_expval(statevector, op):
+ full_op = np.array(1., dtype=complex)
+ for qubit_op in op:
+ if qubit_op.upper() == 'I':
+ full_op = np.kron(full_op, np.eye(2, dtype=complex))
+ elif qubit_op.upper() == 'X':
+ full_op = np.kron(full_op, np.array([[0., 1.], [1., 0.]], dtype=complex))
+ elif qubit_op.upper() == 'Y':
+ full_op = np.kron(full_op, np.array([[0., -1.j], [1.j, 0.]]))
+ elif qubit_op.upper() == 'Z':
+ full_op = np.kron(full_op, np.array([[1., 0.], [0., -1.]], dtype=complex))
+
+ return (statevector.conjugate() @ full_op @ statevector).real
+
+# 例：U(π/3, π/6, 0)|0>
+statevector = np.array([np.cos(np.pi / 6.), np.exp(1.j * np.pi / 6.) * np.sin(np.pi / 6.)])
+print(pauli_expval(statevector, 'X'), pauli_expval(statevector, 'Y'), pauli_expval(statevector, 'Z'))
+```
-npar = 3*nq # パラメータの数
+次に、変分フォーム回路を定義します。このとき、$U$ゲートの回転角として、具体的な数値を設定せず、QiskitのParameterというオブジェクトを利用します。Parameterはあとから数値を代入できる名前のついた箱として使えます。
-np.random.seed(999999)
-target_distr = np.random.rand(2**nq)
-target_distr /= sum(target_distr)
+```{code-cell} ipython3
+:tags: [remove-output]
+
+theta = Parameter('θ')
+phi = Parameter('φ')
+
+ansatz_1q = QuantumCircuit(1)
+ansatz_1q.u(theta, phi, 0., 0)
```
-次に、単一の$U3$変分フォームの3つのパラメータを引数として受け取り、対応する量子回路を返す関数を定義します。
+Parameterに値を代入するには、回路の`bind_parameters`メソッドを利用します。
```{code-cell} ipython3
-def get_var_form(params):
- qr = QuantumRegister(nq, name="q")
- cr = ClassicalRegister(nq, name='c')
- qc = QuantumCircuit(qr, cr)
+# Parameterの値は未定
+ansatz_1q.draw('mpl')
+```
- for i in range(nq):
- qc.u(params[3*i], params[3*i+1], params[3*i+2], qr[i])
+```{code-cell} ipython3
+# thetaとphiにπ/3とπ/6を代入
+ansatz_1q.bind_parameters({theta: np.pi / 3., phi: np.pi / 6.}).draw('mpl')
+```
- for i in range(nq):
- qc.measure(qr[i], cr[i])
- return qc
+変分フォーム回路が作る状態における$X, Y, Z$の期待値を測定するための回路を定義します。
-get_var_form(np.random.rand(npar)).draw('mpl')
+```{code-cell} ipython3
+# を測るにはHゲートで基底を変換する
+x_circuit = ansatz_1q.copy()
+x_circuit.h(0)
+x_circuit.measure_all()
+
+# を測るにはSdg, Hゲートで基底を変換する
+y_circuit = ansatz_1q.copy()
+y_circuit.sdg(0)
+y_circuit.h(0)
+y_circuit.measure_all()
+
+# はそのままの回路で測れる
+z_circuit = ansatz_1q.copy()
+z_circuit.measure_all()
```
-変分フォームのパラメータのリストを入力とし、パラメータに対応したコストを計算する目的関数を定義します。アルゴリズムを実行するバックエンドとして、**QASMシミュレータ**を使用します。
+それぞれの回路を通常通りバックエンドの`run()`メソッドで実行し、結果から期待値を計算する関数を定義します。
```{code-cell} ipython3
backend = AerSimulator()
-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)
+def circuit_expval(circuit, param_vals):
+ bound_circuit = circuit.bind_parameters({theta: param_vals[0], phi: param_vals[1]})
- if len(output_distr) == 1:
- output_distr.append(0)
- return output_distr
+ bound_circuit_tr = transpile(bound_circuit, backend=backend)
+ # shotsは関数の外で定義する
+ job = backend.run(bound_circuit_tr, shots=shots)
+ counts = job.result().get_counts()
-# コストを計算する目的関数を定義
-def objective_function(params):
- qc = get_var_form(params)
- qc = transpile(qc, backend=backend)
- result = backend.run(qc, 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
+ 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.]
+print(circuit_expval(x_circuit, param_vals), circuit_expval(y_circuit, param_vals), circuit_expval(z_circuit, param_vals))
```
-最後にCOBYLAオプティマイザーのインスタンスを作成し、アルゴリズムを実行します。出力される確率分布は実行の度に異なり、ターゲットの確率分布と完全には同じにならないことに注意してください。出力の精度は量子計算の回数（ショット数＝NUM_SHOTS）に依存するので、ショット数を増減させた時の一致具合を確認してみてください。
+最小化する目的関数を定義します。
```{code-cell} ipython3
-optimizer = COBYLA(maxiter=500, tol=0.0001)
+def objective_function(param_vals):
+ loss = 0.
+ for op, circuit in [('X', x_circuit), ('Y', y_circuit), ('Z', z_circuit)]:
+ # target_state_1qは関数の外で定義する
+ target = pauli_expval(target_state_1q, op)
+ current = circuit_expval(circuit, param_vals)
+ loss += (target - current) ** 2
+
+ return loss
+
+# 最適化の1ステップごとに呼び出される関数。目的関数の値をリストに記録しておく
+def callback_function(param_vals):
+ # lossesは関数の外で定義する
+ losses.append(objective_function(param_vals))
+```
-params = np.random.rand(npar)
-ret = optimizer.optimize(num_vars=npar, objective_function=objective_function, initial_point=params)
+最適化には使用しませんが、解を得たあとで変分フォームの終状態とターゲット状態とのフィデリティ$|\langle \psi_0 | \psi(\theta, \phi) \rangle|^2$を計算する関数も定義しておきます。厳密に最適化が成功すれば、この関数の返り値は1になります。
+
+```{code-cell} ipython3
+def fidelity(ansatz, param_vals, target_state):
+ # 量子回路のパラメータのリストはcircuit.parametersで取得できる
+ parameters = ansatz.parameters
-qc = get_var_form(ret[0])
-qc = transpile(qc, backend=backend)
-counts = backend.run(qc, shots=NUM_SHOTS).result().get_counts(qc)
-output_distr = get_probability_distribution(counts)
+ param_binding = dict(zip(parameters, param_vals))
+ opt_ansatz = ansatz.bind_parameters(param_binding)
-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))
+ # QuantumCircuitからOperatorオブジェクトを作ると、AerSimulatorを利用しなくても状態ベクトルを得られる
+ circuit_state = Operator(opt_ansatz).data[:, 0]
+
+ return np.square(np.abs(target_state.conjugate() @ circuit_state))
```
-では次に、この問題を2量子ビット（確率分布の要素数は4）に拡張してやってみましょう。上に戻って
+最後にCOBYLAオプティマイザーのインスタンスを作成し、アルゴリズムを実行します。
```{code-cell} ipython3
----
-jupyter:
- outputs_hidden: false
-pycharm:
- name: '#%%
+# COBYLAの最大ステップ数
+maxiter = 500
+# COBYLAの収束条件（小さいほどよい近似を目指す）
+tol = 0.0001
+# バックエンドでのショット数
+shots = 1000
+
+# オプティマイザーのインスタンス生成
+optimizer = COBYLA(maxiter=maxiter, tol=tol, callback=callback_function)
+```
- '
----
-nq = 2 # 量子ビットの数
+```{code-cell} ipython3
+:tags: [remove-input]
+
+# テキスト作成用のセル
+
+import os
+if os.getenv('JUPYTERBOOK_BUILD') == '1':
+ del optimizer
```
-として再度実行するとどういう結果が得られるでしょうか。量子回路とオプティマイザーの関係はこのようになってますね。
+```{code-cell} ipython3
+:tags: [raises-exception, remove-output]
+
+# ターゲット状態
+target_state_1q = random_statevector(1)
-```{image} figs/vqe_2q_u3.png
-:alt: vqe_2q_u3
-:width: 500px
-:align: center
+# thetaを[0, π), phiを[0, 2π)からランダムに選ぶ
+init = [rng.uniform(0., np.pi), rng.uniform(0., 2. * np.pi)]
+
+# 最適化を実行
+losses = list()
+min_result = optimizer.minimize(objective_function, x0=init)
```
-やってみると分かりますが、結果は1量子ビットの場合と比べて良くないですね。どうすれば良くなるでしょうか？（やり方は複数あると思います）
+```{code-cell} ipython3
+:tags: [remove-input]
-+++
+# テキスト作成用のセル
-**一つの解決策：変分フォームにエンタングルメントを導入する**
+import pickle
+if os.getenv('JUPYTERBOOK_BUILD') == '1':
+ with open('data/vqe_results_1q.pkl', 'rb') as source:
+ min_result, losses = pickle.load(source)
+```
-```python
- 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])
+最適化プロセスにおけるロス（目的関数の返り値）の推移をプロットします。
+
+```{code-cell} ipython3
+plt.plot(losses);
```
-どうなるか確かめてください。
+```{raw-cell}
+`optimizer.minimize()`の返り値の`min_result`から最適化過程の様々な情報（目的関数の呼び出し回数や最適化に要したステップ数など）にアクセスできます。特に、最適化されたパラメータ値が`min_result.x`から得られるので、それを使ってフィデリティを計算してみます。
+```
-+++
+```{code-cell} ipython3
+fidelity(ansatz_1q, min_result.x, target_state_1q)
+```
+
+```{raw-cell}
+ショット数が有限なので統計誤差の影響から最適パラメータは厳密解とは一致しません。ショット数やステップ数を変えて、解の一致具合を確認してみてください。
+```
-量子ビットをエンタングルさせることで相関のあるデータを表現しやすくなるという状況は、例えば、ベル状態（[CHSH不等式の破れを確認する](https://utokyo-icepp.github.io/qc-workbook/chsh_inequality.html#id14)を参照）の確率分布を再現したいときにクリアに見ることができます。上で
+#### Estimatorの利用
+
+VQEを含む変分量子アルゴリズムでは、上のように変分フォームにパラメータ値を代入し複数の観測量の期待値を計算するという手順の繰り返しが頻出します。そのため、これを自動化し、かつ（今は利用しませんが）様々なエラー補正なども適応してくれるEstimatorというクラスを使用することが推奨されています。特に、ここではBackendEstimatorという、特定のバックエンドを利用して計算をするタイプのEstimatorを利用します。
+
+```{code-cell} ipython3
+# 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を使った目的関数を定義します。
+
+```{code-cell} ipython3
+def objective_function_estimator(param_vals):
+ target = list(pauli_expval(target_state_1q, op) for op in ['X', 'Y', 'Z'])
+
+ observables = [SparsePauliOp('X'), SparsePauliOp('Y'), SparsePauliOp('Z')]
+ job = estimator.run([ansatz_1q] * 3, observables, [param_vals] * 3, shots=shots)
+ current = job.result().values
+
+ return np.sum(np.square(np.array(target) - np.array(current)))
+
+def callback_function_estimator(param_vals):
+ # lossesは関数の外で定義する
+ losses.append(objective_function_estimator(param_vals))
+```
+
+上の目的関数を最適化します。
+
+```{code-cell} ipython3
+# COBYLAの最大ステップ数
+maxiter = 500
+# COBYLAの収束条件（小さいほどよい近似を目指す）
+tol = 0.0001
+# バックエンドでのショット数
+shots = 1000
+
+# オプティマイザーのインスタンス生成
+optimizer = COBYLA(maxiter=maxiter, tol=tol, callback=callback_function_estimator)
+```
+
+```{code-cell} ipython3
+:tags: [remove-input]
+
+# テキスト作成用のセル
+
+if os.getenv('JUPYTERBOOK_BUILD') == '1':
+ del optimizer
+```
+
+```{code-cell} ipython3
+:tags: [raises-exception, remove-output]
+
+# ターゲット状態
+target_state_1q = random_statevector(1)
+
+# thetaを[0, π), phiを[0, 2π)からランダムに選ぶ
+init = [rng.uniform(0., np.pi), rng.uniform(0., 2. * np.pi)]
+
+# 最適化を実行
+losses = list()
+min_result = optimizer.minimize(objective_function_estimator, x0=init)
+```
+
+```{code-cell} ipython3
+:tags: [remove-input]
+
+# テキスト作成用のセル
+
+if os.getenv('JUPYTERBOOK_BUILD') == '1':
+ with open('data/vqe_result_1q_estimator.pkl', 'rb') as source:
+ min_result = pickle.load(source)
+```
+
+```{code-cell} ipython3
+fidelity(ansatz_1q, min_result.x, target_state_1q)
+```
+
+### エンタングルメントの導入
+
+では次に、この問題を2量子ビットに拡張してやってみましょう。2量子ビットの純粋状態は6個の実自由度を持ちますが、ここでは最も一般的に2量子ビット状態を決定する15個の観測量の期待値
+
+$$
+\langle O_1 O_2 \rangle \quad (O_1, O_2 = I, X, Y, Z; O_1 O_2 \neq II)
+$$
+
+を測定します。ここで$I$は恒等演算子です。
+
+ターゲット状態に関する関数`random_statevector`と`pauli_expval`はそのまま利用できます。まず変分フォームとして2つの量子ビットに$U$ゲートが一つずつかかっているものを考えて、最小化すべき目的関数を定義します。
+
+```{code-cell} ipython3
+:tags: [remove-output]
+
+# パラメータ数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)
+```
+
+```{code-cell} ipython3
+ops_1q = ['I', 'X', 'Y', 'Z']
+ops = list(f'{op1}{op2}' for op1 in ops_1q for op2 in ops_1q if (op1, op2) != ('I', 'I'))
+
+def objective_function_2q(param_vals):
+ # target_state_2qは関数の外で定義
+ target = list(pauli_expval(target_state_2q, op) for op in ops)
+
+ observables = list(SparsePauliOp(op) for op in ops)
+ job = estimator.run([ansatz_2q] * len(ops), observables, [param_vals] * len(ops), shots=shots)
+ current = job.result().values
+
+ return np.sum(np.square(np.array(target) - np.array(current)))
+
+def callback_function_2q(param_vals):
+ # lossesは関数の外で定義する
+ losses.append(objective_function_2q(param_vals))
+```
```{code-cell} ipython3
---
+jupyter:
+ outputs_hidden: false
pycharm:
name: '#%%
'
---
-target_distr = np.random.rand(2**nq)
+# COBYLAの最大ステップ数
+maxiter = 500
+# 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)
```
-を
+```{code-cell} ipython3
+:tags: [remove-input]
+
+# テキスト作成用のセル
+
+if os.getenv('JUPYTERBOOK_BUILD') == '1':
+ del optimizer
+```
```{code-cell} ipython3
---
+jupyter:
+ outputs_hidden: false
pycharm:
name: '#%%
'
+tags: [raises-exception, remove-output]
---
-# 00と11を測定する確率が50%、01と10の確率は0
-target_distr = np.array([0.5,0.,0.,0.5])
+# 最適化を実行
+losses = list()
+min_result = optimizer.minimize(objective_function_2q, x0=init)
```
-として実行するとどうなるでしょうか。エンタングルさせる場合とさせない場合で大きな違いが見えるでしょう。3量子ビットのGHZ状態（[単純な量子回路をゼロから書く](https://utokyo-icepp.github.io/qc-workbook/circuit_from_scratch.html#ghz)を参照）
+```{code-cell} ipython3
+:tags: [remove-input]
+
+# テキスト作成用のセル
+
+if os.getenv('JUPYTERBOOK_BUILD') == '1':
+ with open('data/vqe_result_2q.pkl', 'rb') as source:
+ min_result = pickle.load(source)
+```
```{code-cell} ipython3
---
+jupyter:
+ outputs_hidden: false
pycharm:
name: '#%%
'
---
-# 000と111を測定する確率が50%、それ以外の確率は0
-target_distr = np.array([0.5,0.,0.,0.,0.,0.,0.,0.5])
+fidelity(ansatz_2q, min_result.x, target_state_2q)
+```
+
+やってみると分かりますが、結果は1量子ビットの場合と比べて良くないですね。どうすれば良くなるでしょうか？（やり方は複数あると思います）
+
++++
+
+**一つの解決策：変分フォームにエンタングルメントを導入する**
+
+```python
+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不等式の破れを確認する](https://utokyo-icepp.github.io/qc-workbook/chsh_inequality.html#id14)を参照）を再現したいときにこの状況をクリアに見ることができます。上で
+
+```python
+target_state_2q = random_statevector(2)
```
-に拡張してみるなどして、遊んでみてください。
+を
+
+```python
+target_state_2q = np.array([1., 0., 0., 1.], dtype=complex) / np.sqrt(2.)
+```
+
+として実行するとどうなるでしょうか。エンタングルさせる場合とさせない場合で大きな違いが見えるでしょう。
+
+問題を3量子ビットに拡張し、GHZ状態（[単純な量子回路をゼロから書く](https://utokyo-icepp.github.io/qc-workbook/circuit_from_scratch.html#ghz)を参照）
+
+```python
+target_state_3q = np.array([1.] + [0.] * 6 + [1.], dtype=complex) / np.sqrt(2.)
+```
+
+をターゲットにするなどして、遊んでみてください。
+++ {"pycharm": {"name": "#%% md\n"}}
@@ -386,75 +685,68 @@ $$
### VQEの実装
では、パラメータシフト法を使って簡単なVQEの例を実装してみます。ある観測量の期待値が最小になるように、VQEを使ってAnsatzを更新する回路パラメータを決定するという問題を考えてみます。
-量子回路として、$R_YR_Z$ゲートを繰り返すシンプルなパラメータ回路を使うことにします。
+量子回路として、$R_YR_Z$ゲートを繰り返すシンプルなパラメータ回路を使い、観測量として、パウリ演算子のテンソル積$ZXY$を使います。
-```{code-cell} ipython3
----
-jupyter:
- outputs_hidden: false
-pycharm:
- name: '#%%
+パラメータシフト法の実装は、QiskitではParamShiftEstimatorGradientというAPIを使うことで一行で済んでしまいます（実際の勾配計算に興味がある人は、期待値から勾配を直接求めるコードを書いて、パラメータ毎に$\pm\pi/2$シフトさせた回路を走らせることで、このAPIと同じ勾配が得られることを確認してみて下さい）。パラメータの最適化は、勾配を使って勾配降下を行うConjugate Descent (CG)とGradient Descentの2つのオプティマイザーを使って行いますが、比較のためにCOBYLAも使うことにします。
- '
----
-# Tested with python 3.8.12, qiskit 0.34.2, numpy 1.22.2
-from qiskit.circuit import ParameterVector
-from qiskit.utils import QuantumInstance
-from qiskit.opflow import I, X, Y, Z, StateFn, CircuitStateFn
-from qiskit.opflow.gradients import Gradient
-from qiskit.algorithms import VQE, NumPyMinimumEigensolver
+最終的に、3通りのVQEを使って求めた最小エネルギーの近似解を、厳密対角化して求めた最小エネルギーの値と比較することにします。
+
+```{code-cell} ipython3
+from qiskit.algorithms.minimum_eigensolvers import VQE, NumPyMinimumEigensolver
from qiskit.algorithms.optimizers import CG, GradientDescent
+from qiskit.algorithms.gradients import ParamShiftEstimatorGradient
```
```{code-cell} ipython3
----
-jupyter:
- outputs_hidden: false
-pycharm:
- name: '#%%
+# Ansatzの定義
- '
----
-n = 3 # 量子ビット数
-nl = 2 # レイヤー数
-npar = n*2*nl # パラメータ数
-
-qc = QuantumCircuit(n)
-param_list = ParameterVector('param_list',npar)
-for i in range(nl):
- qc.ry(param_list[6*i], 0)
- qc.ry(param_list[6*i+1], 1)
- qc.ry(param_list[6*i+2], 2)
- qc.rz(param_list[6*i+3], 0)
- qc.rz(param_list[6*i+4], 1)
- qc.rz(param_list[6*i+5], 2)
- #qc.cnot(0, 1)
- #qc.cnot(1, 2)
-```
+num_qubits = 3 # 量子ビット数
+num_layers = 2 # レイヤー数
+num_params = num_qubits * 2 * num_layers # パラメータ数
-+++ {"pycharm": {"name": "#%% md\n"}}
+ansatz = QuantumCircuit(num_qubits)
-観測量として、パウリ演算子のテンソル積$ZXY$を使います。
+param_list = ParameterVector('param_list', num_params)
+iparam = 0
+for _ in range(num_layers):
+ for iq in range(num_qubits):
+ ansatz.ry(param_list[iparam], iq)
+ iparam += 1
-パラメータシフト法の実装は、QiskitではGradientというAPIを使うことで一行で済んでしまいます（実際の勾配計算に興味がある人は、期待値から勾配を直接求めるコードを書いて、パラメータ毎に$\pm\pi/2$シフトさせた回路を走らせることで、このAPIと同じ勾配が得られることを確認してみて下さい）。パラメータの最適化は、勾配を使って勾配降下を行うConjugate Descent (CG)とGradient Descentの2つのオプティマイザーを使って行いますが、比較のためにCOBYLAも使うことにします。
+ for iq in range(num_qubits):
+ ansatz.rz(param_list[iparam], iq)
+ iparam += 1
-最終的に、3通りのVQEを使って求めた最小エネルギーの近似解を、厳密対角化して求めた最小エネルギーの値と比較することにします。
+ #for iq in range(num_qubits - 1):
+ # ansatz.cx(iq, iq + 1)
+
+ansatz.draw('mpl')
+```
```{code-cell} ipython3
----
-jupyter:
- outputs_hidden: false
-pycharm:
- name: '#%%
+# 最小固有値を求める観測量
+obs = SparsePauliOp('ZXY')
- '
----
-obs = Z ^ X ^ Y
+# パラメータの初期値
+init = rng.uniform(0., 2. * np.pi, size=num_params)
+
+# 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)
-grad = Gradient(grad_method="param_shift")
-cg = CG(maxiter=200)
-gd = GradientDescent(maxiter=200)
-cobyla = COBYLA(maxiter=300)
+# 厳密解を計算するソルバー
+ee = NumPyMinimumEigensolver()
```
```{code-cell} ipython3
@@ -462,40 +754,16 @@ cobyla = COBYLA(maxiter=300)
# テキスト作成用のセル
-import os
if os.getenv('JUPYTERBOOK_BUILD') == '1':
- del qc
+ del obs
```
```{code-cell} ipython3
----
-jupyter:
- outputs_hidden: false
-pycharm:
- name: '#%%
+:tags: [raises-exception, remove-output]
- '
-tags: [raises-exception, remove-output]
----
-random_seed = 10598
-backend = Aer.get_backend('qasm_simulator')
-quantum_instance = QuantumInstance(backend=backend,
- shots=1024,
- seed_simulator=random_seed,
- seed_transpiler=random_seed,
- skip_qobj_validation=True)
-# VQEアルゴリズムの実装
-vqe_gfree = VQE(ansatz=qc, optimizer=cobyla, quantum_instance=quantum_instance)
-result_vqe_gfree = vqe_gfree.compute_minimum_eigenvalue(obs)
-
-vqe_cg = VQE(ansatz=qc, optimizer=cg, gradient=grad, quantum_instance=quantum_instance)
result_vqe_cg = vqe_cg.compute_minimum_eigenvalue(obs)
-
-vqe_gd = VQE(ansatz=qc, optimizer=gd, gradient=grad, quantum_instance=quantum_instance)
result_vqe_gd = vqe_gd.compute_minimum_eigenvalue(obs)
-
-# 厳密解
-ee = NumPyMinimumEigensolver()
+result_vqe_cobyla = vqe_cobyla.compute_minimum_eigenvalue(obs)
result_ee = ee.compute_minimum_eigenvalue(obs)
```
@@ -504,9 +772,8 @@ result_ee = ee.compute_minimum_eigenvalue(obs)
# テキスト作成用のセルなので無視してよい
-import pickle
with open('data/vqe_results.pkl', 'rb') as source:
- result_ee, result_vqe_gfree, result_vqe_cg, result_vqe_gd = pickle.load(source)
+ result_ee, result_vqe_cobyla, result_vqe_cg, result_vqe_gd = pickle.load(source)
```
```{code-cell} ipython3
@@ -520,7 +787,7 @@ pycharm:
---
print('Result:')
print(f' Exact = {result_ee.eigenvalue}')
-print(f' VQE(GFree) = {result_vqe_gfree.optimal_value}')
+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}')
```