diff --git a/qiskit/algorithms/__init__.py b/qiskit/algorithms/__init__.py
index 452f9a4f1220..63e44798a710 100644
--- a/qiskit/algorithms/__init__.py
+++ b/qiskit/algorithms/__init__.py
@@ -108,6 +108,8 @@
RealEvolver
ImaginaryEvolver
TrotterQRTE
+ PVQD
+ PVQDResult
EvolutionResult
EvolutionProblem
@@ -246,6 +248,7 @@
from .exceptions import AlgorithmError
from .aux_ops_evaluator import eval_observables
from .evolvers.trotterization import TrotterQRTE
+from .evolvers.pvqd import PVQD, PVQDResult
__all__ = [
"AlgorithmResult",
@@ -292,6 +295,8 @@
"PhaseEstimationScale",
"PhaseEstimation",
"PhaseEstimationResult",
+ "PVQD",
+ "PVQDResult",
"IterativePhaseEstimation",
"AlgorithmError",
"eval_observables",
diff --git a/qiskit/algorithms/evolvers/evolution_problem.py b/qiskit/algorithms/evolvers/evolution_problem.py
index b069effe1747..e0f9fe3063c6 100644
--- a/qiskit/algorithms/evolvers/evolution_problem.py
+++ b/qiskit/algorithms/evolvers/evolution_problem.py
@@ -31,7 +31,7 @@ def __init__(
self,
hamiltonian: OperatorBase,
time: float,
- initial_state: Union[StateFn, QuantumCircuit],
+ initial_state: Optional[Union[StateFn, QuantumCircuit]] = None,
aux_operators: Optional[ListOrDict[OperatorBase]] = None,
truncation_threshold: float = 1e-12,
t_param: Optional[Parameter] = None,
@@ -41,7 +41,9 @@ def __init__(
Args:
hamiltonian: The Hamiltonian under which to evolve the system.
time: Total time of evolution.
- initial_state: Quantum state to be evolved.
+ initial_state: The quantum state to be evolved for methods like Trotterization.
+ For variational time evolutions, where the evolution happens in an ansatz,
+ this argument is not required.
aux_operators: Optional list of auxiliary operators to be evaluated with the
evolved ``initial_state`` and their expectation values returned.
truncation_threshold: Defines a threshold under which values can be assumed to be 0.
diff --git a/qiskit/algorithms/evolvers/pvqd/__init__.py b/qiskit/algorithms/evolvers/pvqd/__init__.py
new file mode 100644
index 000000000000..9377ce631b4e
--- /dev/null
+++ b/qiskit/algorithms/evolvers/pvqd/__init__.py
@@ -0,0 +1,18 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022.
+#
+# 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.
+
+"""The projected Variational Quantum Dynamic (p-VQD) module."""
+
+from .pvqd_result import PVQDResult
+from .pvqd import PVQD
+
+__all__ = ["PVQD", "PVQDResult"]
diff --git a/qiskit/algorithms/evolvers/pvqd/pvqd.py b/qiskit/algorithms/evolvers/pvqd/pvqd.py
new file mode 100644
index 000000000000..e4a1d5893bb5
--- /dev/null
+++ b/qiskit/algorithms/evolvers/pvqd/pvqd.py
@@ -0,0 +1,414 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2019, 2022.
+#
+# 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.
+
+"""The projected Variational Quantum Dynamics Algorithm."""
+
+from typing import Optional, Union, List, Tuple, Callable
+
+import logging
+import numpy as np
+
+from qiskit import QiskitError
+from qiskit.algorithms.optimizers import Optimizer, Minimizer
+from qiskit.circuit import QuantumCircuit, ParameterVector
+from qiskit.circuit.library import PauliEvolutionGate
+from qiskit.extensions import HamiltonianGate
+from qiskit.providers import Backend
+from qiskit.opflow import OperatorBase, CircuitSampler, ExpectationBase, StateFn, MatrixOp
+from qiskit.synthesis import EvolutionSynthesis, LieTrotter
+from qiskit.utils import QuantumInstance
+
+from .pvqd_result import PVQDResult
+from .utils import _get_observable_evaluator, _is_gradient_supported
+
+from ..evolution_problem import EvolutionProblem
+from ..evolution_result import EvolutionResult
+from ..real_evolver import RealEvolver
+
+logger = logging.getLogger(__name__)
+
+
+class PVQD(RealEvolver):
+ """The projected Variational Quantum Dynamics (p-VQD) Algorithm.
+
+ In each timestep, this algorithm computes the next state with a Trotter formula
+ (specified by the ``evolution`` argument) and projects the timestep onto a variational form
+ (``ansatz``). The projection is determined by maximizing the fidelity of the Trotter-evolved
+ state and the ansatz, using a classical optimization routine. See Ref. [1] for details.
+
+ The following attributes can be set via the initializer but can also be read and
+ updated once the PVQD object has been constructed.
+
+ Attributes:
+
+ ansatz (QuantumCircuit): The parameterized circuit representing the time-evolved state.
+ initial_parameters (np.ndarray): The parameters of the ansatz at time 0.
+ expectation (ExpectationBase): The method to compute expectation values.
+ optimizer (Optional[Union[Optimizer, Minimizer]]): The classical optimization routine
+ used to maximize the fidelity of the Trotter step and ansatz.
+ num_timesteps (Optional[int]): The number of timesteps to take. If None, it is automatically
+ selected to achieve a timestep of approximately 0.01.
+ evolution (Optional[EvolutionSynthesis]): The method to perform the Trotter step.
+ Defaults to first-order Lie-Trotter evolution.
+ use_parameter_shift (bool): If True, use the parameter shift rule for loss function
+ gradients (if the ansatz supports).
+ initial_guess (Optional[np.ndarray]): The starting point for the first classical optimization
+ run, at time 0. Defaults to random values in :math:`[-0.01, 0.01]`.
+
+ Example:
+
+ This snippet computes the real time evolution of a quantum Ising model on two
+ neighboring sites and keeps track of the magnetization.
+
+ .. code-block:: python
+
+ import numpy as np
+
+ from qiskit import BasicAer
+ from qiskit.circuit.library import EfficientSU2
+ from qiskit.opflow import X, Z, I, MatrixExpectation
+
+ backend = BasicAer.get_backend("statevector_simulator")
+ expectation = MatrixExpectation()
+ hamiltonian = 0.1 * (Z ^ Z) + (I ^ X) + (X ^ I)
+ observable = Z ^ Z
+ ansatz = EfficientSU2(2, reps=1)
+ initial_parameters = np.zeros(ansatz.num_parameters)
+
+ time = 1
+ optimizer = L_BFGS_B()
+
+ # setup the algorithm
+ pvqd = PVQD(
+ ansatz,
+ initial_parameters,
+ num_timesteps=100,
+ optimizer=optimizer,
+ quantum_instance=backend,
+ expectation=expectation
+ )
+
+ # specify the evolution problem
+ problem = EvolutionProblem(
+ hamiltonian, time, aux_operators=[hamiltonian, observable]
+ )
+
+ # and evolve!
+ result = pvqd.evolve(problem)
+
+ References:
+
+ [1] Stefano Barison, Filippo Vicentini, and Giuseppe Carleo (2021), An efficient
+ quantum algorithm for the time evolution of parameterized circuits,
+ `Quantum 5, 512 `_.
+ """
+
+ def __init__(
+ self,
+ ansatz: QuantumCircuit,
+ initial_parameters: np.ndarray,
+ expectation: ExpectationBase,
+ optimizer: Optional[Union[Optimizer, Minimizer]] = None,
+ num_timesteps: Optional[int] = None,
+ evolution: Optional[EvolutionSynthesis] = None,
+ use_parameter_shift: bool = True,
+ initial_guess: Optional[np.ndarray] = None,
+ quantum_instance: Optional[Union[Backend, QuantumInstance]] = None,
+ ) -> None:
+ """
+ Args:
+ ansatz: A parameterized circuit preparing the variational ansatz to model the
+ time evolved quantum state.
+ initial_parameters: The initial parameters for the ansatz. Together with the ansatz,
+ these define the initial state of the time evolution.
+ expectation: The expectation converter to evaluate expectation values.
+ optimizer: The classical optimizers used to minimize the overlap between
+ Trotterization and ansatz. Can be either a :class:`.Optimizer` or a callable
+ using the :class:`.Minimizer` protocol. This argument is optional since it is
+ not required for :meth:`get_loss`, but it has to be set before :meth:`evolve`
+ is called.
+ num_timestep: The number of time steps. If ``None`` it will be set such that the timestep
+ is close to 0.01.
+ evolution: The evolution synthesis to use for the construction of the Trotter step.
+ Defaults to first-order Lie-Trotter decomposition, see also
+ :mod:`~qiskit.synthesis.evolution` for different options.
+ use_parameter_shift: If True, use the parameter shift rule to compute gradients.
+ If False, the optimizer will not be passed a gradient callable. In that case,
+ Qiskit optimizers will use a finite difference rule to approximate the gradients.
+ initial_guess: The initial guess for the first VQE optimization. Afterwards the
+ previous iteration result is used as initial guess. If None, this is set to
+ a random vector with elements in the interval :math:`[-0.01, 0.01]`.
+ quantum_instance: The backend or quantum instance used to evaluate the circuits.
+ """
+ if evolution is None:
+ evolution = LieTrotter()
+
+ self.ansatz = ansatz
+ self.initial_parameters = initial_parameters
+ self.num_timesteps = num_timesteps
+ self.optimizer = optimizer
+ self.initial_guess = initial_guess
+ self.expectation = expectation
+ self.evolution = evolution
+ self.use_parameter_shift = use_parameter_shift
+
+ self._sampler = None
+ self.quantum_instance = quantum_instance
+
+ @property
+ def quantum_instance(self) -> Optional[QuantumInstance]:
+ """Return the current quantum instance."""
+ return self._quantum_instance
+
+ @quantum_instance.setter
+ def quantum_instance(self, quantum_instance: Optional[Union[Backend, QuantumInstance]]) -> None:
+ """Set the quantum instance and circuit sampler."""
+ if quantum_instance is not None:
+ if not isinstance(quantum_instance, QuantumInstance):
+ quantum_instance = QuantumInstance(quantum_instance)
+ self._sampler = CircuitSampler(quantum_instance)
+
+ self._quantum_instance = quantum_instance
+
+ def step(
+ self,
+ hamiltonian: OperatorBase,
+ ansatz: QuantumCircuit,
+ theta: np.ndarray,
+ dt: float,
+ initial_guess: np.ndarray,
+ ) -> Tuple[np.ndarray, float]:
+ """Perform a single time step.
+
+ Args:
+ hamiltonian: The Hamiltonian under which to evolve.
+ ansatz: The parameterized quantum circuit which attempts to approximate the
+ time-evolved state.
+ theta: The current parameters.
+ dt: The time step.
+ initial_guess: The initial guess for the classical optimization of the
+ fidelity between the next variational state and the Trotter-evolved last state.
+ If None, this is set to a random vector with elements in the interval
+ :math:`[-0.01, 0.01]`.
+
+ Returns:
+ A tuple consisting of the next parameters and the fidelity of the optimization.
+ """
+ self._validate_setup()
+
+ loss, gradient = self.get_loss(hamiltonian, ansatz, dt, theta)
+
+ if initial_guess is None:
+ initial_guess = np.random.random(self.initial_parameters.size) * 0.01
+
+ if isinstance(self.optimizer, Optimizer):
+ optimizer_result = self.optimizer.minimize(loss, initial_guess, gradient)
+ else:
+ optimizer_result = self.optimizer(loss, initial_guess, gradient)
+
+ # clip the fidelity to [0, 1]
+ fidelity = np.clip(1 - optimizer_result.fun, 0, 1)
+
+ return theta + optimizer_result.x, fidelity
+
+ def get_loss(
+ self,
+ hamiltonian: OperatorBase,
+ ansatz: QuantumCircuit,
+ dt: float,
+ current_parameters: np.ndarray,
+ ) -> Tuple[Callable[[np.ndarray], float], Optional[Callable[[np.ndarray], np.ndarray]]]:
+
+ """Get a function to evaluate the infidelity between Trotter step and ansatz.
+
+ Args:
+ hamiltonian: The Hamiltonian under which to evolve.
+ ansatz: The parameterized quantum circuit which attempts to approximate the
+ time-evolved state.
+ dt: The time step.
+ current_parameters: The current parameters.
+
+ Returns:
+ A callable to evaluate the infidelity and, if gradients are supported and required,
+ a second callable to evaluate the gradient of the infidelity.
+ """
+ self._validate_setup(skip={"optimizer"})
+
+ # use Trotterization to evolve the current state
+ trotterized = ansatz.bind_parameters(current_parameters)
+
+ if isinstance(hamiltonian, MatrixOp):
+ evolution_gate = HamiltonianGate(hamiltonian.primitive, time=dt)
+ else:
+ evolution_gate = PauliEvolutionGate(hamiltonian, time=dt, synthesis=self.evolution)
+
+ trotterized.append(evolution_gate, ansatz.qubits)
+
+ # define the overlap of the Trotterized state and the ansatz
+ x = ParameterVector("w", ansatz.num_parameters)
+ shifted = ansatz.assign_parameters(current_parameters + x)
+ overlap = StateFn(trotterized).adjoint() @ StateFn(shifted)
+
+ converted = self.expectation.convert(overlap)
+
+ def evaluate_loss(
+ displacement: Union[np.ndarray, List[np.ndarray]]
+ ) -> Union[float, List[float]]:
+ """Evaluate the overlap of the ansatz with the Trotterized evolution.
+
+ Args:
+ displacement: The parameters for the ansatz.
+
+ Returns:
+ The fidelity of the ansatz with parameters ``theta`` and the Trotterized evolution.
+ """
+ if isinstance(displacement, list):
+ displacement = np.asarray(displacement)
+ value_dict = {x_i: displacement[:, i].tolist() for i, x_i in enumerate(x)}
+ else:
+ value_dict = dict(zip(x, displacement))
+
+ sampled = self._sampler.convert(converted, params=value_dict)
+
+ # in principle we could add different loss functions here, but we're currently
+ # not aware of a use-case for a different one than in the paper
+ return 1 - np.abs(sampled.eval()) ** 2
+
+ if _is_gradient_supported(ansatz) and self.use_parameter_shift:
+
+ def evaluate_gradient(displacement: np.ndarray) -> np.ndarray:
+ """Evaluate the gradient with the parameter-shift rule.
+
+ This is hardcoded here since the gradient framework does not support computing
+ gradients for overlaps.
+
+ Args:
+ displacement: The parameters for the ansatz.
+
+ Returns:
+ The gradient.
+ """
+ # construct lists where each element is shifted by plus (or minus) pi/2
+ dim = displacement.size
+ plus_shifts = (displacement + np.pi / 2 * np.identity(dim)).tolist()
+ minus_shifts = (displacement - np.pi / 2 * np.identity(dim)).tolist()
+
+ evaluated = evaluate_loss(plus_shifts + minus_shifts)
+
+ gradient = (evaluated[:dim] - evaluated[dim:]) / 2
+
+ return gradient
+
+ else:
+ evaluate_gradient = None
+
+ return evaluate_loss, evaluate_gradient
+
+ def evolve(self, evolution_problem: EvolutionProblem) -> EvolutionResult:
+ """
+ Args:
+ evolution_problem: The evolution problem containing the hamiltonian, total evolution
+ time and observables to evaluate.
+
+ Returns:
+ A result object containing the evolution information and evaluated observables.
+
+ Raises:
+ ValueError: If the evolution time is not positive or the timestep is too small.
+ NotImplementedError: If the evolution problem contains an initial state.
+ """
+ self._validate_setup()
+
+ time = evolution_problem.time
+ observables = evolution_problem.aux_operators
+ hamiltonian = evolution_problem.hamiltonian
+
+ # determine the number of timesteps and set the timestep
+ num_timesteps = (
+ int(np.ceil(time / 0.01)) if self.num_timesteps is None else self.num_timesteps
+ )
+ timestep = time / num_timesteps
+
+ if evolution_problem.initial_state is not None:
+ raise NotImplementedError(
+ "Setting an initial state for the evolution is not yet supported for PVQD."
+ )
+
+ # get the function to evaluate the observables for a given set of ansatz parameters
+ if observables is not None:
+ evaluate_observables = _get_observable_evaluator(
+ self.ansatz, observables, self.expectation, self._sampler
+ )
+ observable_values = [evaluate_observables(self.initial_parameters)]
+
+ fidelities = [1]
+ parameters = [self.initial_parameters]
+ times = np.linspace(0, time, num_timesteps + 1).tolist() # +1 to include initial time 0
+
+ initial_guess = self.initial_guess
+
+ for _ in range(num_timesteps):
+ # perform VQE to find the next parameters
+ next_parameters, fidelity = self.step(
+ hamiltonian, self.ansatz, parameters[-1], timestep, initial_guess
+ )
+
+ # set initial guess to last parameter update
+ initial_guess = next_parameters - parameters[-1]
+
+ parameters.append(next_parameters)
+ fidelities.append(fidelity)
+ if observables is not None:
+ observable_values.append(evaluate_observables(next_parameters))
+
+ evolved_state = self.ansatz.bind_parameters(parameters[-1])
+
+ result = PVQDResult(
+ evolved_state=evolved_state,
+ times=times,
+ parameters=parameters,
+ fidelities=fidelities,
+ estimated_error=1 - np.prod(fidelities),
+ )
+ if observables is not None:
+ result.observables = observable_values
+ result.aux_ops_evaluated = observable_values[-1]
+
+ return result
+
+ def _validate_setup(self, skip=None):
+ """Validate the current setup and raise an error if something misses to run."""
+
+ if skip is None:
+ skip = {}
+
+ required_attributes = {"quantum_instance", "optimizer"}.difference(skip)
+
+ for attr in required_attributes:
+ if getattr(self, attr, None) is None:
+ raise ValueError(f"The {attr} cannot be None.")
+
+ if self.num_timesteps is not None and self.num_timesteps <= 0:
+ raise ValueError(
+ f"The number of timesteps must be positive but is {self.num_timesteps}."
+ )
+
+ if self.ansatz.num_parameters == 0:
+ raise QiskitError(
+ "The ansatz cannot have 0 parameters, otherwise it cannot be trained."
+ )
+
+ if len(self.initial_parameters) != self.ansatz.num_parameters:
+ raise QiskitError(
+ f"Mismatching number of parameters in the ansatz ({self.ansatz.num_parameters}) "
+ f"and the initial parameters ({len(self.initial_parameters)})."
+ )
diff --git a/qiskit/algorithms/evolvers/pvqd/pvqd_result.py b/qiskit/algorithms/evolvers/pvqd/pvqd_result.py
new file mode 100644
index 000000000000..e14e7a0c9db5
--- /dev/null
+++ b/qiskit/algorithms/evolvers/pvqd/pvqd_result.py
@@ -0,0 +1,55 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022.
+#
+# 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.
+
+"""Result object for p-VQD."""
+
+from typing import Union, Optional, List, Tuple
+import numpy as np
+
+from qiskit.circuit import QuantumCircuit
+from qiskit.opflow import StateFn, OperatorBase
+
+from ..evolution_result import EvolutionResult
+
+
+class PVQDResult(EvolutionResult):
+ """The result object for the p-VQD algorithm."""
+
+ def __init__(
+ self,
+ evolved_state: Union[StateFn, QuantumCircuit, OperatorBase],
+ aux_ops_evaluated: Optional[List[Tuple[complex, complex]]] = None,
+ times: Optional[List[float]] = None,
+ parameters: Optional[List[np.ndarray]] = None,
+ fidelities: Optional[List[float]] = None,
+ estimated_error: Optional[float] = None,
+ observables: Optional[List[List[float]]] = None,
+ ):
+ """
+ Args:
+ evolved_state: An evolved quantum state.
+ aux_ops_evaluated: Optional list of observables for which expected values on an evolved
+ state are calculated. These values are in fact tuples formatted as (mean, standard
+ deviation).
+ times: The times evaluated during the time integration.
+ parameters: The parameter values at each evaluation time.
+ fidelities: The fidelity of the Trotter step and variational update at each iteration.
+ estimated_error: The overall estimated error evaluated as one minus the
+ product of all fidelities.
+ observables: The value of the observables evaluated at each iteration.
+ """
+ super().__init__(evolved_state, aux_ops_evaluated)
+ self.times = times
+ self.parameters = parameters
+ self.fidelities = fidelities
+ self.estimated_error = estimated_error
+ self.observables = observables
diff --git a/qiskit/algorithms/evolvers/pvqd/utils.py b/qiskit/algorithms/evolvers/pvqd/utils.py
new file mode 100644
index 000000000000..589f12005de8
--- /dev/null
+++ b/qiskit/algorithms/evolvers/pvqd/utils.py
@@ -0,0 +1,100 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022.
+#
+# 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.
+
+
+"""Utilities for p-VQD."""
+
+from typing import Union, List, Callable
+import logging
+
+import numpy as np
+
+from qiskit.circuit import QuantumCircuit, Parameter, ParameterExpression
+from qiskit.compiler import transpile
+from qiskit.exceptions import QiskitError
+from qiskit.opflow import ListOp, CircuitSampler, ExpectationBase, StateFn, OperatorBase
+from qiskit.opflow.gradients.circuit_gradients import ParamShift
+
+logger = logging.getLogger(__name__)
+
+
+def _is_gradient_supported(ansatz: QuantumCircuit) -> bool:
+ """Check whether we can apply a simple parameter shift rule to obtain gradients."""
+
+ # check whether the circuit can be unrolled to supported gates
+ try:
+ unrolled = transpile(ansatz, basis_gates=ParamShift.SUPPORTED_GATES, optimization_level=0)
+ except QiskitError:
+ # failed to map to supported basis
+ logger.log(
+ logging.INFO,
+ "No gradient support: Failed to unroll to gates supported by parameter-shift.",
+ )
+ return False
+
+ # check whether all parameters are unique and we do not need to apply the chain rule
+ # (since it's not implemented yet)
+ total_num_parameters = 0
+ for circuit_instruction in unrolled.data:
+ for param in circuit_instruction.operation.params:
+ if isinstance(param, ParameterExpression):
+ if isinstance(param, Parameter):
+ total_num_parameters += 1
+ else:
+ logger.log(
+ logging.INFO,
+ "No gradient support: Circuit is only allowed to have plain parameters, "
+ "as the chain rule is not yet implemented.",
+ )
+ return False
+
+ if total_num_parameters != ansatz.num_parameters:
+ logger.log(
+ logging.INFO,
+ "No gradient support: Circuit is only allowed to have unique parameters, "
+ "as the product rule is not yet implemented.",
+ )
+ return False
+
+ return True
+
+
+def _get_observable_evaluator(
+ ansatz: QuantumCircuit,
+ observables: Union[OperatorBase, List[OperatorBase]],
+ expectation: ExpectationBase,
+ sampler: CircuitSampler,
+) -> Callable[[np.ndarray], Union[float, List[float]]]:
+ """Get a callable to evaluate a (list of) observable(s) for given circuit parameters."""
+
+ if isinstance(observables, list):
+ observables = ListOp(observables)
+
+ expectation_value = StateFn(observables, is_measurement=True) @ StateFn(ansatz)
+ converted = expectation.convert(expectation_value)
+
+ ansatz_parameters = ansatz.parameters
+
+ def evaluate_observables(theta: np.ndarray) -> Union[float, List[float]]:
+ """Evaluate the observables for the ansatz parameters ``theta``.
+
+ Args:
+ theta: The ansatz parameters.
+
+ Returns:
+ The observables evaluated at the ansatz parameters.
+ """
+ value_dict = dict(zip(ansatz_parameters, theta))
+ sampled = sampler.convert(converted, params=value_dict)
+ return sampled.eval()
+
+ return evaluate_observables
diff --git a/releasenotes/notes/project-dynamics-2f848a5f89655429.yaml b/releasenotes/notes/project-dynamics-2f848a5f89655429.yaml
new file mode 100644
index 000000000000..a33fdfefed28
--- /dev/null
+++ b/releasenotes/notes/project-dynamics-2f848a5f89655429.yaml
@@ -0,0 +1,45 @@
+features:
+ - |
+ Added the :class:`PVQD` class to the time evolution framework. This class implements the
+ projected Variational Quantum Dynamics (p-VQD) algorithm as :class:`.PVQD` of
+ `Barison et al. `_.
+
+ In each timestep this algorithm computes the next state with a Trotter formula and projects it
+ onto a variational form. The projection is determined by maximizing the fidelity of the
+ Trotter-evolved state and the ansatz, using a classical optimization routine.
+
+ .. code-block:: python
+
+ import numpy as np
+
+ from qiskit import BasicAer
+ from qiskit.circuit.library import EfficientSU2
+ from qiskit.opflow import X, Z, I, MatrixExpectation
+
+ backend = BasicAer.get_backend("statevector_simulator")
+ expectation = MatrixExpectation()
+ hamiltonian = 0.1 * (Z ^ Z) + (I ^ X) + (X ^ I)
+ observable = Z ^ Z
+ ansatz = EfficientSU2(2, reps=1)
+ initial_parameters = np.zeros(ansatz.num_parameters)
+
+ time = 1
+ optimizer = L_BFGS_B()
+
+ # setup the algorithm
+ pvqd = PVQD(
+ ansatz,
+ initial_parameters,
+ num_timesteps=100,
+ optimizer=optimizer,
+ quantum_instance=backend,
+ expectation=expectation
+ )
+
+ # specify the evolution problem
+ problem = EvolutionProblem(
+ hamiltonian, time, aux_operators=[hamiltonian, observable]
+ )
+
+ # and evolve!
+ result = pvqd.evolve(problem)
diff --git a/test/python/algorithms/evolvers/test_pvqd.py b/test/python/algorithms/evolvers/test_pvqd.py
new file mode 100644
index 000000000000..5c450fa24822
--- /dev/null
+++ b/test/python/algorithms/evolvers/test_pvqd.py
@@ -0,0 +1,325 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2018, 2022.
+#
+# 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.
+
+"""Tests for PVQD."""
+
+from functools import partial
+from ddt import ddt, data, unpack
+import numpy as np
+
+from qiskit.test import QiskitTestCase
+
+from qiskit import BasicAer, QiskitError
+from qiskit.circuit import QuantumCircuit, Parameter, Gate
+from qiskit.algorithms.evolvers import EvolutionProblem
+from qiskit.algorithms.evolvers.pvqd import PVQD
+from qiskit.algorithms.optimizers import L_BFGS_B, GradientDescent, SPSA, OptimizerResult
+from qiskit.circuit.library import EfficientSU2
+from qiskit.opflow import X, Z, I, MatrixExpectation, PauliExpectation
+
+
+# pylint: disable=unused-argument, invalid-name
+def gradient_supplied(fun, x0, jac, info):
+ """A mock optimizer that checks whether the gradient is supported or not."""
+ result = OptimizerResult()
+ result.x = x0
+ result.fun = 0
+ info["has_gradient"] = jac is not None
+
+ return result
+
+
+class WhatAmI(Gate):
+ """An custom opaque gate that can be inverted but not decomposed."""
+
+ def __init__(self, angle):
+ super().__init__(name="whatami", num_qubits=2, params=[angle])
+
+ def inverse(self):
+ return WhatAmI(-self.params[0])
+
+
+@ddt
+class TestPVQD(QiskitTestCase):
+ """Tests for the pVQD algorithm."""
+
+ def setUp(self):
+ super().setUp()
+ self.sv_backend = BasicAer.get_backend("statevector_simulator")
+ self.qasm_backend = BasicAer.get_backend("qasm_simulator")
+ self.expectation = MatrixExpectation()
+ self.hamiltonian = 0.1 * (Z ^ Z) + (I ^ X) + (X ^ I)
+ self.observable = Z ^ Z
+ self.ansatz = EfficientSU2(2, reps=1)
+ self.initial_parameters = np.zeros(self.ansatz.num_parameters)
+
+ @data(
+ ("ising", MatrixExpectation, True, "sv", 2),
+ ("ising_matrix", MatrixExpectation, True, "sv", None),
+ ("ising", PauliExpectation, True, "qasm", 2),
+ ("pauli", PauliExpectation, False, "qasm", None),
+ )
+ @unpack
+ def test_pvqd(self, hamiltonian_type, expectation_cls, gradient, backend_type, num_timesteps):
+ """Test a simple evolution."""
+ time = 0.02
+
+ if hamiltonian_type == "ising":
+ hamiltonian = self.hamiltonian
+ elif hamiltonian_type == "ising_matrix":
+ hamiltonian = self.hamiltonian.to_matrix_op()
+ else: # hamiltonian_type == "pauli":
+ hamiltonian = X ^ X
+
+ # parse input arguments
+ if gradient:
+ optimizer = GradientDescent(maxiter=1)
+ else:
+ optimizer = L_BFGS_B(maxiter=1)
+
+ backend = self.sv_backend if backend_type == "sv" else self.qasm_backend
+ expectation = expectation_cls()
+
+ # run pVQD keeping track of the energy and the magnetization
+ pvqd = PVQD(
+ self.ansatz,
+ self.initial_parameters,
+ num_timesteps=num_timesteps,
+ optimizer=optimizer,
+ quantum_instance=backend,
+ expectation=expectation,
+ )
+ problem = EvolutionProblem(hamiltonian, time, aux_operators=[hamiltonian, self.observable])
+ result = pvqd.evolve(problem)
+
+ self.assertTrue(len(result.fidelities) == 3)
+ self.assertTrue(np.all(result.times == [0.0, 0.01, 0.02]))
+ self.assertTrue(np.asarray(result.observables).shape == (3, 2))
+ num_parameters = self.ansatz.num_parameters
+ self.assertTrue(
+ len(result.parameters) == 3
+ and np.all([len(params) == num_parameters for params in result.parameters])
+ )
+
+ def test_step(self):
+ """Test calling the step method directly."""
+
+ pvqd = PVQD(
+ self.ansatz,
+ self.initial_parameters,
+ optimizer=L_BFGS_B(maxiter=100),
+ quantum_instance=self.sv_backend,
+ expectation=MatrixExpectation(),
+ )
+
+ # perform optimization for a timestep of 0, then the optimal parameters are the current
+ # ones and the fidelity is 1
+ theta_next, fidelity = pvqd.step(
+ self.hamiltonian.to_matrix_op(),
+ self.ansatz,
+ self.initial_parameters,
+ dt=0.0,
+ initial_guess=np.zeros_like(self.initial_parameters),
+ )
+
+ self.assertTrue(np.allclose(theta_next, self.initial_parameters))
+ self.assertAlmostEqual(fidelity, 1)
+
+ def test_get_loss(self):
+ """Test getting the loss function directly."""
+
+ pvqd = PVQD(
+ self.ansatz,
+ self.initial_parameters,
+ quantum_instance=self.sv_backend,
+ expectation=MatrixExpectation(),
+ use_parameter_shift=False,
+ )
+
+ theta = np.ones(self.ansatz.num_parameters)
+ loss, gradient = pvqd.get_loss(
+ self.hamiltonian, self.ansatz, dt=0.0, current_parameters=theta
+ )
+
+ displacement = np.arange(self.ansatz.num_parameters)
+
+ with self.subTest(msg="check gradient is None"):
+ self.assertIsNone(gradient)
+
+ with self.subTest(msg="check loss works"):
+ self.assertGreater(loss(displacement), 0)
+ self.assertAlmostEqual(loss(np.zeros_like(theta)), 0)
+
+ def test_invalid_num_timestep(self):
+ """Test raises if the num_timestep is not positive."""
+ pvqd = PVQD(
+ self.ansatz,
+ self.initial_parameters,
+ num_timesteps=0,
+ optimizer=L_BFGS_B(),
+ quantum_instance=self.sv_backend,
+ expectation=self.expectation,
+ )
+ problem = EvolutionProblem(
+ self.hamiltonian, time=0.01, aux_operators=[self.hamiltonian, self.observable]
+ )
+
+ with self.assertRaises(ValueError):
+ _ = pvqd.evolve(problem)
+
+ def test_initial_guess_and_observables(self):
+ """Test doing no optimizations stays at initial guess."""
+ initial_guess = np.zeros(self.ansatz.num_parameters)
+
+ pvqd = PVQD(
+ self.ansatz,
+ self.initial_parameters,
+ num_timesteps=10,
+ optimizer=SPSA(maxiter=0, learning_rate=0.1, perturbation=0.01),
+ initial_guess=initial_guess,
+ quantum_instance=self.sv_backend,
+ expectation=self.expectation,
+ )
+ problem = EvolutionProblem(
+ self.hamiltonian, time=0.1, aux_operators=[self.hamiltonian, self.observable]
+ )
+
+ result = pvqd.evolve(problem)
+
+ observables = result.aux_ops_evaluated
+ self.assertEqual(observables[0], 0.1) # expected energy
+ self.assertEqual(observables[1], 1) # expected magnetization
+
+ def test_missing_attributesquantum_instance(self):
+ """Test appropriate error is raised if the quantum instance is missing."""
+ pvqd = PVQD(
+ self.ansatz,
+ self.initial_parameters,
+ optimizer=L_BFGS_B(maxiter=1),
+ expectation=self.expectation,
+ )
+ problem = EvolutionProblem(self.hamiltonian, time=0.01)
+
+ attrs_to_test = [
+ ("optimizer", L_BFGS_B(maxiter=1)),
+ ("quantum_instance", self.qasm_backend),
+ ]
+
+ for attr, value in attrs_to_test:
+ with self.subTest(msg=f"missing: {attr}"):
+ # set attribute to None to invalidate the setup
+ setattr(pvqd, attr, None)
+
+ with self.assertRaises(ValueError):
+ _ = pvqd.evolve(problem)
+
+ # set the correct value again
+ setattr(pvqd, attr, value)
+
+ with self.subTest(msg="all set again"):
+ result = pvqd.evolve(problem)
+ self.assertIsNotNone(result.evolved_state)
+
+ def test_zero_parameters(self):
+ """Test passing an ansatz with zero parameters raises an error."""
+ problem = EvolutionProblem(self.hamiltonian, time=0.02)
+
+ pvqd = PVQD(
+ QuantumCircuit(2),
+ np.array([]),
+ optimizer=SPSA(maxiter=10, learning_rate=0.1, perturbation=0.01),
+ quantum_instance=self.sv_backend,
+ expectation=self.expectation,
+ )
+
+ with self.assertRaises(QiskitError):
+ _ = pvqd.evolve(problem)
+
+ def test_initial_state_raises(self):
+ """Test passing an initial state raises an error for now."""
+ initial_state = QuantumCircuit(2)
+ initial_state.x(0)
+
+ problem = EvolutionProblem(
+ self.hamiltonian,
+ time=0.02,
+ initial_state=initial_state,
+ )
+
+ pvqd = PVQD(
+ self.ansatz,
+ self.initial_parameters,
+ optimizer=SPSA(maxiter=0, learning_rate=0.1, perturbation=0.01),
+ quantum_instance=self.sv_backend,
+ expectation=self.expectation,
+ )
+
+ with self.assertRaises(NotImplementedError):
+ _ = pvqd.evolve(problem)
+
+
+class TestPVQDUtils(QiskitTestCase):
+ """Test some utility functions for PVQD."""
+
+ def setUp(self):
+ super().setUp()
+ self.sv_backend = BasicAer.get_backend("statevector_simulator")
+ self.expectation = MatrixExpectation()
+ self.hamiltonian = 0.1 * (Z ^ Z) + (I ^ X) + (X ^ I)
+ self.ansatz = EfficientSU2(2, reps=1)
+
+ def test_gradient_supported(self):
+ """Test the gradient support is correctly determined."""
+ # gradient supported here
+ wrapped = EfficientSU2(2) # a circuit wrapped into a big instruction
+ plain = wrapped.decompose() # a plain circuit with already supported instructions
+
+ # gradients not supported on the following circuits
+ x = Parameter("x")
+ duplicated = QuantumCircuit(2)
+ duplicated.rx(x, 0)
+ duplicated.rx(x, 1)
+
+ needs_chainrule = QuantumCircuit(2)
+ needs_chainrule.rx(2 * x, 0)
+
+ custom_gate = WhatAmI(x)
+ unsupported = QuantumCircuit(2)
+ unsupported.append(custom_gate, [0, 1])
+
+ tests = [
+ (wrapped, True), # tuple: (circuit, gradient support)
+ (plain, True),
+ (duplicated, False),
+ (needs_chainrule, False),
+ (unsupported, False),
+ ]
+
+ # used to store the info if a gradient callable is passed into the
+ # optimizer of not
+ info = {"has_gradient": None}
+ optimizer = partial(gradient_supplied, info=info)
+
+ pvqd = PVQD(
+ ansatz=None,
+ initial_parameters=np.array([]),
+ optimizer=optimizer,
+ quantum_instance=self.sv_backend,
+ expectation=self.expectation,
+ )
+ problem = EvolutionProblem(self.hamiltonian, time=0.01)
+ for circuit, expected_support in tests:
+ with self.subTest(circuit=circuit, expected_support=expected_support):
+ pvqd.ansatz = circuit
+ pvqd.initial_parameters = np.zeros(circuit.num_parameters)
+ _ = pvqd.evolve(problem)
+ self.assertEqual(info["has_gradient"], expected_support)