注釈

このページは tutorials/algorithms/10_pvqd.ipynb から生成されました。

投影変分量子ダイナミクス

投影変分量子ダイナミクス(p-VQD)アルゴリズムは、実時間発展のアルゴリズムための量子アルゴリズムです。 Trotterizationで計算されたように、 \(t + \Delta_t\) 時の状態をパラメータ化された量子回路に投影する変分アルゴリズムです。

パラメーター化された量子回路 \(U(\theta)\) と Hamiltonian \(H\) で構成された量子状態 \(|\phi(\theta)\rangle = U(\theta)|0\rangle\) の更新ルールは次のように書き表せます。

\[\theta_{n+1} = \theta_n + \arg\min_{\delta\theta} 1 - |\langle\phi(\theta_n + \delta\theta)|e^{-i\Delta_t H}|\phi(\theta_n)\rangle|^2,\]

ここで \(e^{-i\Delta_t H}\) はトロッター展開で計算されます 。(例えば、Qiskitの `PauliEvolutionGate <https://qiskit.org/documentation/stubs/qiskit.circuit.library.PauliEvolutionGate.html>`__ を使って！)

以下のチュートリアルでは、 qiskit.algorithms.time_evolvers.PVQD として実装され利用可能な、Qiskit の p-VQD アルゴリズムについて説明します。アルゴリズムの詳細については、元の論文である Barison et al. Quantum 5, 512 (2021) を参照してください。

例えば、以下のハミルトニアンのもとでの状態 \(|00\rangle\) の時間発展について見ていきます。

\[H = 0.1 Z_1 Z_2 + X_1 + X_2\]

という2 つの隣接するスピン上のイジング・ハミルトニアンで、\(T=1\) を observableとして経過を追いたいとします。

from qiskit.quantum_info import SparsePauliOp

final_time = 1
hamiltonian = SparsePauliOp.from_sparse_list([
    ("ZZ", [0, 1], 0.1), ("X", [0], 1), ("X", [1], 1),
], num_qubits=2)
observable = SparsePauliOp(["ZI", "IZ"])

Hamiltonianとobservableを定義した後、更新ステップを投影するパラメータ化されたansatzを選択する必要があります。 ここでは様々な選択肢がありますが、時間発展するハミルトニアンの構成要素を含むansatzは実時間発展のために通常非常にうまく機能します。

from qiskit.circuit import QuantumCircuit, ParameterVector

theta = ParameterVector("th", 5)
ansatz = QuantumCircuit(2)
ansatz.rx(theta[0], 0)
ansatz.rx(theta[1], 1)
ansatz.rzz(theta[2], 0, 1)
ansatz.rx(theta[3], 0)
ansatz.rx(theta[4], 1)

# you can try different circuits, like:
# from qiskit.circuit.library import EfficientSU2
# ansatz = EfficientSU2(2, reps=1)

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

この ansatz では、すべてのパラメーターが 0である場合、 \(|00\rangle\) 状態ベクトルになります。 したがって、最初のパラメーターは \(\theta_0 = 0\) に設定します。:

import numpy as np

initial_parameters = np.zeros(ansatz.num_parameters)

P-VQDアルゴリズムを実行する前に、バックエンドと期待値の計算方法を選択する必要があります。 ここでは、 qiskit.primitives にあるprimitiveの実装を参照し、厳密な状態ベクトルシミュレーション(我々が研究しているのは2量子ビットのため、非常に高速です。) を実行します。

from qiskit.primitives import Sampler, Estimator
from qiskit.algorithms.state_fidelities import ComputeUncompute

# the fidelity is used to evaluate the objective: the overlap of the variational form and the trotter step
sampler = Sampler()
fidelity = ComputeUncompute(sampler)

# the estimator is used to evaluate the observables
estimator = Estimator()

Since p-VQD performs a classical optimization in each timestep to determine the best parameters for the projection, we also have to specify the classical optimizer. As a first example we’re using BFGS, which typically works well in statevector simulations, but later we can switch to gradient descent.

from qiskit.algorithms.optimizers import L_BFGS_B

bfgs = L_BFGS_B()

Now we can define p-VQD and execute it!

from qiskit.algorithms.time_evolvers.pvqd import PVQD

pvqd = PVQD(
    fidelity,
    ansatz,
    initial_parameters,
    estimator=estimator,
    num_timesteps=100,
    optimizer=bfgs
)

The p-VQD implementation follows Qiskit’s time evolution interface, thus we pack all information of the evolution problem into an input class: the hamiltonian under which we evolve the state, the final_time of the evolution and the observables (aux_operators) we keep track of.

from qiskit.algorithms.time_evolvers.time_evolution_problem import TimeEvolutionProblem

problem = TimeEvolutionProblem(hamiltonian, time=final_time, aux_operators=[hamiltonian, observable])

And then run the algorithm!

result = pvqd.evolve(problem)

Now we can have a look at the results, which are stored in a PVQDResult object. This class has the fields

• evolved_state: The quantum circuit with the parameters at the final evolution time.

• times: The timesteps of the time integration. At these times we have the parameter values and evaluated the observables.

• parameters: The parameter values at each timestep.

• observables: The observable values at each timestep.

• fidelities: The fidelity of projecting the Trotter timestep onto the variational form at each timestep.

• estimated_error: The estimated error as product of all fidelities.

The energy should be constant in a real time evolution. However, we are projecting the time-evolved state onto a variational form, which might violate this rule. Ideally the energy is still more or less constant. In this evolution here we observe shifts of ~5% of the energy.

import matplotlib.pyplot as plt

energies = np.real(result.observables)[:, 0]

plt.plot(result.times, energies, color="royalblue")
plt.xlabel("time $t$")
plt.ylabel("energy $E$")
plt.title("Energy over time")

Text(0.5, 1.0, 'Energy over time')

Since we also kept track of the total magnetization of the system, we can plot that quantity too. However let’s first compute exact reference values to verify our algorithm results.

import scipy as sc

def exact(final_time, timestep, hamiltonian, initial_state):
    """Get the exact values for energy and the observable."""
    O = observable.to_matrix()
    H = hamiltonian.to_matrix()

    energ, magn = [], []  # list of energies and magnetizations evaluated at timesteps timestep
    times = []  # list of timepoints at which energy/obs are evaluated
    time = 0
    while time <= final_time:
        # get exact state at time t
        exact_state = initial_state.evolve(sc.linalg.expm(-1j * time * H))
        # store observables and time
        times.append(time)
        energ.append(exact_state.expectation_value(H).real)
        magn.append(exact_state.expectation_value(observable).real)

        # next timestep
        time += timestep

    return times, energ, magn

from qiskit.quantum_info import Statevector

initial_state = Statevector(ansatz.bind_parameters(initial_parameters))
exact_times, exact_energies, exact_magnetizations = exact(final_time, 0.01, hamiltonian, initial_state)

magnetizations = np.real(result.observables)[:, 1]

plt.plot(result.times, magnetizations.real, color="crimson", label="PVQD")
plt.plot(exact_times, exact_magnetizations, ":", color="k", label="Exact")
plt.xlabel("time $t$")
plt.ylabel(r"magnetization $\langle Z_1 Z_2 \rangle$")
plt.title("Magnetization over time")
plt.legend(loc="best")

<matplotlib.legend.Legend at 0x7fbc98539480>

Looks pretty good!

Gradient-based optimizations

The PVQD class also implements parameter-shift gradients for the loss function and we can use a gradient descent optimization routine

\[\theta_{k+1} = \theta_{k} - \eta_k \nabla\ell(\theta_k).\]

Here we’re using a learning rate of

\[\eta_k = 0.1 k^{-0.602}\]

and 80 optimization steps in each timestep.

from qiskit.algorithms.optimizers import GradientDescent

maxiter = 80
learning_rate = 0.1 * np.arange(1, maxiter + 1) ** (-0.602)
gd = GradientDescent(maxiter, lambda: iter(learning_rate))

pvqd.optimizer = gd

The following cell would take a few minutes to run for 100 timesteps, so we reduce them here.

n = 10
pvqd.num_timesteps = n
problem.time = 0.1

result_gd = pvqd.evolve(problem)

energies_gd = np.real(result_gd.observables)[:, 0]

plt.plot(result.times[:n + 1], energies[:n + 1], "-", color="royalblue", label="BFGS")
plt.plot(result_gd.times, energies_gd, "--", color="royalblue", label="Gradient descent")
plt.plot(exact_times[:n + 1], exact_energies[:n + 1], ":", color="k", label="Exact")
plt.legend(loc="best")
plt.xlabel("time $t$")
plt.ylabel("energy $E$")
plt.title("Energy over time")

Text(0.5, 1.0, 'Energy over time')

We can observe here, that the energy does vary quite a bit! But as we mentioned before, p-VQD does not preserve the energy.

magnetizations_gd = np.real(result_gd.observables)[:, 1]

plt.plot(result.times[:n + 1], magnetizations[:n + 1], "-", color="crimson", label="BFGS")
plt.plot(result_gd.times, magnetizations_gd, "--", color="crimson", label="Gradient descent")
plt.plot(exact_times[:n + 1], exact_magnetizations[:n + 1], ":", color="k", label="Exact")
plt.legend(loc="best")
plt.xlabel("time $t$")
plt.ylabel(r"magnetization $\langle Z_1 + Z_2 \rangle$")
plt.title("Magnetization over time")

Text(0.5, 1.0, 'Magnetization over time')

The magnetization, however, is computed very precisely.

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit Software	Version
qiskit-terra	0.24.0.dev0+d814ad4
qiskit-aer	0.11.1
qiskit-ignis	0.7.1
qiskit-ibmq-provider	0.19.2
qiskit-nature	0.6.0
qiskit-finance	0.3.4
qiskit-optimization	0.6.0
qiskit-machine-learning	0.6.0
System information
Python version	3.10.4
Python compiler	Clang 12.0.0
Python build	main, Mar 31 2022 03:38:35
OS	Darwin
CPUs	4
Memory (Gb)	32.0
Tue May 09 13:16:11 2023 CEST

This code is a part of Qiskit

© Copyright IBM 2017, 2023.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.