Implemented PVQD algorithm with primitives.

qiskit/algorithms/__init__.py

@@ -113,8 +113,6 @@
     RealEvolver
     ImaginaryEvolver
     TrotterQRTE
-    PVQD
-    PVQDResult
     EvolutionResult
     EvolutionProblem
 
@@ -132,6 +130,8 @@
 
     RealTimeEvolver
     ImaginaryTimeEvolver
+    PVQD
+    PVQDResult
     TimeEvolutionResult
     TimeEvolutionProblem
 
@@ -310,7 +310,7 @@
 from .observables_evaluator import estimate_observables
 from .evolvers.trotterization import TrotterQRTE
 
-from .evolvers.pvqd import PVQD, PVQDResult
+from .time_evolvers.pvqd import PVQD, PVQDResult
 
 __all__ = [
     "AlgorithmJob",

qiskit/algorithms/state_fidelities/base_state_fidelity.py

@@ -19,7 +19,7 @@
 import numpy as np
 
 from qiskit import QuantumCircuit
-from qiskit.algorithms import AlgorithmJob
+from qiskit.algorithms.algorithm_job import AlgorithmJob
 from qiskit.circuit import ParameterVector
 from .state_fidelity_result import StateFidelityResult
 

qiskit/algorithms/evolvers/pvqd/pvqd.py → qiskit/algorithms/time_evolvers/pvqd/pvqd.py

@@ -11,33 +11,32 @@
 # that they have been altered from the originals.
 
 """The projected Variational Quantum Dynamics Algorithm."""
-
-from typing import Optional, Union, List, Tuple, Callable
+from __future__ import annotations
 
 import logging
+from typing import Callable
+
 import numpy as np
 
-from qiskit import QiskitError
-from qiskit.algorithms.optimizers import Optimizer, Minimizer
-from qiskit.circuit import QuantumCircuit, ParameterVector
+from qiskit.circuit import QuantumCircuit, ParameterVector, Parameter
 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.opflow import PauliSumOp
+from qiskit.primitives import BaseEstimator
+from qiskit.quantum_info.operators.base_operator import BaseOperator
 from qiskit.synthesis import EvolutionSynthesis, LieTrotter
-from qiskit.utils import QuantumInstance
-
+from ...exceptions import AlgorithmError, QiskitError
 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
+from ..time_evolution_problem import TimeEvolutionProblem
+from ..time_evolution_result import TimeEvolutionResult
+from ..real_time_evolver import RealTimeEvolver
+from ...state_fidelities.base_state_fidelity import BaseStateFidelity
+from ...optimizers import Optimizer, Minimizer
 
 logger = logging.getLogger(__name__)
 
 
-class PVQD(RealEvolver):
+class PVQD(RealTimeEvolver):
     """The projected Variational Quantum Dynamics (p-VQD) Algorithm.
 
     In each timestep, this algorithm computes the next state with a Trotter formula
@@ -52,7 +51,6 @@ class PVQD(RealEvolver):
 
         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
@@ -73,14 +71,19 @@ class PVQD(RealEvolver):
 
             import numpy as np
 
+            from qiskit.algorithms.state_fidelities import ComputeUncompute
+            from qiskit.algorithms.time_evolvers.pvqd import PVQD
+            from qiskit.primitives import Estimator
             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
+            from qiskit.quantum_info import Pauli, SparsePauliOp
+            from qiskit.algorithms.optimizers import L_BFGS_B
+
+            sampler = Sampler()
+            fidelity = ComputeUncompute(sampler)
+            estimator = Estimator()
+            hamiltonian = 0.1 * SparsePauliOp([Pauli("ZZ"), Pauli("IX"), Pauli("XI")])
+            observable = Pauli("ZZ")
             ansatz = EfficientSU2(2, reps=1)
             initial_parameters = np.zeros(ansatz.num_parameters)
 
@@ -89,12 +92,12 @@ class PVQD(RealEvolver):
 
             # setup the algorithm
             pvqd = PVQD(
+                fidelity,
                 ansatz,
+                estimator,
                 initial_parameters,
                 num_timesteps=100,
                 optimizer=optimizer,
-                quantum_instance=backend,
-                expectation=expectation
             )
 
             # specify the evolution problem
@@ -114,30 +117,32 @@ class PVQD(RealEvolver):
 
     def __init__(
         self,
+        fidelity: BaseStateFidelity,
         ansatz: QuantumCircuit,
         initial_parameters: np.ndarray,
-        expectation: ExpectationBase,
-        optimizer: Optional[Union[Optimizer, Minimizer]] = None,
-        num_timesteps: Optional[int] = None,
-        evolution: Optional[EvolutionSynthesis] = None,
+        estimator: BaseEstimator | None = None,
+        optimizer: Optimizer | Minimizer | None = None,
+        num_timesteps: int | None = None,
+        evolution: EvolutionSynthesis | None = None,
         use_parameter_shift: bool = True,
-        initial_guess: Optional[np.ndarray] = None,
-        quantum_instance: Optional[Union[Backend, QuantumInstance]] = None,
+        initial_guess: np.ndarray | None = None,
     ) -> None:
         """
         Args:
+            fidelity: A fidelity primitive used by the algorithm.
             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.
+            estimator: An estimator primitive used for calculating expected values of auxiliary
+                operators (if provided via the problem).
             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.
+            num_timesteps: 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.
@@ -147,7 +152,6 @@ def __init__(
             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.
         """
         super().__init__()
         if evolution is None:
@@ -158,36 +162,19 @@ def __init__(
         self.num_timesteps = num_timesteps
         self.optimizer = optimizer
         self.initial_guess = initial_guess
-        self.expectation = expectation
+        self.estimator = estimator
+        self.fidelity_primitive = fidelity
         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,
+        hamiltonian: BaseOperator | PauliSumOp,
         ansatz: QuantumCircuit,
         theta: np.ndarray,
         dt: float,
         initial_guess: np.ndarray,
-    ) -> Tuple[np.ndarray, float]:
+    ) -> tuple[np.ndarray, float]:
         """Perform a single time step.
 
         Args:
@@ -223,11 +210,11 @@ def step(
 
     def get_loss(
         self,
-        hamiltonian: OperatorBase,
+        hamiltonian: BaseOperator | PauliSumOp,
         ansatz: QuantumCircuit,
         dt: float,
         current_parameters: np.ndarray,
-    ) -> Tuple[Callable[[np.ndarray], float], Optional[Callable[[np.ndarray], np.ndarray]]]:
+    ) -> tuple[Callable[[np.ndarray], float], Callable[[np.ndarray], np.ndarray]] | None:
 
         """Get a function to evaluate the infidelity between Trotter step and ansatz.
 
@@ -247,42 +234,50 @@ def get_loss(
         # 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)
+        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]]:
+        def evaluate_loss(displacement: np.ndarray | list[np.ndarray]) -> 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.
+
+            Raises:
+                AlgorithmError: If a primitive job fails.
             """
             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
+            param_dicts = self._transpose_param_dicts(value_dict)
+            num_of_param_sets = len(param_dicts)
+            states1 = [trotterized] * num_of_param_sets
+            states2 = [shifted] * num_of_param_sets
+            param_dicts2 = [list(param_dict.values()) for param_dict in param_dicts]
+            # the first state does not have free parameters so values_1 will be None by default
+            try:
+                job = self.fidelity_primitive.run(states1, states2, values_2=param_dicts2)
+                fidelities = job.result().fidelities
+            except Exception as exc:
+                raise AlgorithmError("The primitive job failed!") from exc
+
+            if len(fidelities) == 1:
+                fidelities = fidelities[0]
+
+            # 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
+            return 1 - np.abs(fidelities) ** 2
 
         if _is_gradient_supported(ansatz) and self.use_parameter_shift:
 
@@ -314,8 +309,25 @@ def evaluate_gradient(displacement: np.ndarray) -> np.ndarray:
 
         return evaluate_loss, evaluate_gradient
 
-    def evolve(self, evolution_problem: EvolutionProblem) -> EvolutionResult:
-        """
+    def _transpose_param_dicts(self, params: dict) -> list[dict[Parameter, float]]:
+        p_0 = list(params.values())[0]
+        if isinstance(p_0, (list, np.ndarray)):
+            num_parameterizations = len(p_0)
+            param_bindings = [
+                {param: value_list[i] for param, value_list in params.items()}  # type: ignore
+                for i in range(num_parameterizations)
+            ]
+        else:
+            param_bindings = [params]
+
+        return param_bindings
+
+    def evolve(self, evolution_problem: TimeEvolutionProblem) -> TimeEvolutionResult:
+        r"""Perform real time evolution :math:`\exp(-i t H)|\Psi\rangle`.
+
+        Evolves an initial state :math:`|\Psi\rangle` for a time :math:`t`
+        under a Hamiltonian  :math:`H`, as provided in the ``evolution_problem``.
+
         Args:
             evolution_problem: The evolution problem containing the hamiltonian, total evolution
                 time and observables to evaluate.
@@ -324,7 +336,8 @@ def evolve(self, evolution_problem: EvolutionProblem) -> EvolutionResult:
             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.
+            ValueError: If ``aux_operators`` provided in the time evolution problem but no estimator
+                provided to the algorithm.
             NotImplementedError: If the evolution problem contains an initial state.
         """
         self._validate_setup()
@@ -346,8 +359,12 @@ def evolve(self, evolution_problem: EvolutionProblem) -> EvolutionResult:
 
         # get the function to evaluate the observables for a given set of ansatz parameters
         if observables is not None:
+            if self.estimator is None:
+                raise ValueError(
+                    "The evolution problem contained aux_operators but no estimator was provided. "
+                )
             evaluate_observables = _get_observable_evaluator(
-                self.ansatz, observables, self.expectation, self._sampler
+                self.ansatz, observables, self.estimator
             )
             observable_values = [evaluate_observables(self.initial_parameters)]
 
@@ -392,7 +409,7 @@ def _validate_setup(self, skip=None):
         if skip is None:
             skip = {}
 
-        required_attributes = {"quantum_instance", "optimizer"}.difference(skip)
+        required_attributes = {"optimizer"}.difference(skip)
 
         for attr in required_attributes:
             if getattr(self, attr, None) is None: