VQE implementation with estimator primitive (#8702)

* Add minimally working VQE with estimator primitive implementation.

No gradients, no tests etc.

* Revert from dataclass results to original result classes.

* Enforce positional and keyword VQE arguments.

* Move aux op eval logic to function

* Update docstring.

* Remove max_evals_grouped. Force to set directly on optimizer.

* Remove validate min import.

* Make note that eval_observables will be used to eval aux ops.

* Add initial vqe tests.

* Have VQE inherit from VariationalAlgorithm.

* Move energy evaluation to unnested function.

* Construct h2_op using SparsePauliOp

* Add gradient with primitives support.

* Update docstrings

* update broadcast handling

* update eval_energy output for batching

* add incomplete QNSPA test

* fix batch evaluation of QNSPSA test

* remove vqe callback

* move estimator to first arg

* remove usused imports

* add minimum eigensolvers test init file

* add aux ops tests and prepare for new eval_operators

* no longer support account for Nones in aux_ops

* correct typing for MinimumEigensolver

* Compute default initial point using ansatz bounds.

* Add NumPyMinimumEigensolverResult

* Fix type hints

* Fix type hints

* Formatting

* Do not store NumPyMES result inside the algo.

* Provide default values for ansatz and estimator

* Formatting

* fix old and new batching

* Add tests for NumpyMES and import in module.

* Use lazy formatting in log messages

* Use lazy formatting in log messages

* Add back callback to VQE.

* minor renaming

* raise algorithm error if primitive jobs fail

* Add return documentation

* Improve var names and docstrings.

* Apply suggestions from code review

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>
Co-authored-by: dlasecki <dal@zurich.ibm.com>

* minor formatting

* minor formatting

* Ensure evaluate_energy/gradient match

Co-authored-by: Max Rossmannek <max.rossmannek@uzh.ch>

* return bounds logic; fix some docstrings and minor refactor

* Force keyword arguments in vqe.

* Use estimate_observables function

* break up numpy mes tests with subTest

* formatting

* remove redundant eval_aux_ops

* add typehints to docstring attributes

* remove usused imports

* remove default ansatz

* remove usused imports

* update typehints

* avoid changing the original ansatz

* avoid changing the original ansatz

* create separate function to build vqe result

* Correct aux operator eignvalue type hint

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>

* add variance to callback

* use std_dev in callback rather than variance

* formatting

* return variance and shots in callback

* return full metadata in callback

* Move validation functions to algorithms/utils

* correct the callback attribute documentation

* correct the callback attribute typehint docstring

* update VQE docstring

* release note and pending-depreciate old algs

* update vqe class docstring

* Apply suggestions from code review

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>

* Do not copy ansatz

* Note pending depreciation of old algs

* fix docstrings and imports

* Fix issues with building docs

* Include OptimizerResult in VQEResult

* Remove trailing whitespace

* Fix math notation in docstring

* estimate_obervables to return metadata @ElePT +VQE

* add example in release note

* Update evaluate_observables docstring

* Fix observables_evaluator tests.

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>

* fix trotter_qrte tests and remove depreciation

* formatting

* remove unused import

* Apply suggestions to docstring from code review

Co-authored-by: Steve Wood <40241007+woodsp-ibm@users.noreply.github.com>

* Update VQE docstring

* Remove printing

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>

* Add parameters to estimate observables

* Fix estimate obs. unit test

* Update arg description

* Update arg description

* keep equation part of sentence

* dict -> dict[str, Any]

* Update qiskit/algorithms/optimizers/spsa.py

Apply suggestion Imamichi-san

* introduce FilterType and aux_operator_eigenvalues -> aux_operators_evaluated

* Correct typehint

Co-authored-by: Ikko Hamamura <ikkoham@users.noreply.github.com>

* update vqe docstring and use old typing for type alias

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>
Co-authored-by: dlasecki <dal@zurich.ibm.com>
Co-authored-by: Max Rossmannek <max.rossmannek@uzh.ch>
Co-authored-by: Steve Wood <40241007+woodsp-ibm@users.noreply.github.com>
Co-authored-by: ElePT <epenatap@gmail.com>
Co-authored-by: Ikko Hamamura <ikkoham@users.noreply.github.com>

qiskit/algorithms/__init__.py

@@ -183,10 +183,12 @@ Algorithms to solve linear systems of equations.
    linear_solvers
 
 
-Minimum Eigensolvers
---------------------
+Minimum Eigen Solvers
+---------------------
 
 Algorithms that can find the minimum eigenvalue of an operator.
+These algorithms are pending depreciation. One should instead make use of the
+Minimum Eigensolver classes in the section below, which leverage Runtime primitives.
 
 .. autosummary::
    :toctree: ../stubs/
@@ -203,6 +205,16 @@ Algorithms that can find the minimum eigenvalue of an operator.
    QAOA
    VQE
 
+Minimum Eigensolvers
+--------------------
+
+Algorithms that can find the minimum eigenvalue of an operator and leverage primitives.
+
+.. autosummary::
+   :toctree: ../stubs/
+
+   minimum_eigensolvers
+
 
 Optimizers
 ----------

qiskit/algorithms/minimum_eigensolvers/__init__.py

@@ -0,0 +1,48 @@
+# 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.
+
+"""
+============================================================================
+Minimum Eigensolvers Package (:mod:`qiskit.algorithms.minimum_eigensolvers`)
+============================================================================
+
+.. currentmodule:: qiskit.algorithms.minimum_eigensolvers
+
+Minimum Eigensolvers
+====================
+.. autosummary::
+   :toctree: ../stubs/
+
+   MinimumEigensolver
+   NumPyMinimumEigensolver
+   VQE
+
+.. autosummary::
+   :toctree: ../stubs/
+
+   MinimumEigensolverResult
+   NumPyMinimumEigensolverResult
+   VQEResult
+"""
+
+from .minimum_eigensolver import MinimumEigensolver, MinimumEigensolverResult
+from .numpy_minimum_eigensolver import NumPyMinimumEigensolver, NumPyMinimumEigensolverResult
+from .vqe import VQE, VQEResult
+
+__all__ = [
+    "MinimumEigensolver",
+    "MinimumEigensolverResult",
+    "NumPyMinimumEigensolver",
+    "NumPyMinimumEigensolverResult",
+    "VQE",
+    "VQEResult",
+]

qiskit/algorithms/minimum_eigensolvers/minimum_eigensolver.py

@@ -0,0 +1,99 @@
+# 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 minimum eigensolver interface and result."""
+
+from __future__ import annotations
+
+from abc import ABC, abstractmethod
+from typing import Any
+
+from qiskit.opflow import PauliSumOp
+from qiskit.quantum_info.operators.base_operator import BaseOperator
+
+from ..algorithm_result import AlgorithmResult
+from ..list_or_dict import ListOrDict
+
+
+class MinimumEigensolver(ABC):
+    """The minimum eigensolver interface.
+
+    Algorithms that can compute a minimum eigenvalue for an operator may implement this interface to
+    allow different algorithms to be used interchangeably.
+    """
+
+    @abstractmethod
+    def compute_minimum_eigenvalue(
+        self,
+        operator: BaseOperator | PauliSumOp,
+        aux_operators: ListOrDict[BaseOperator | PauliSumOp] | None = None,
+    ) -> "MinimumEigensolverResult":
+        """
+        Computes the minimum eigenvalue. The ``operator`` and ``aux_operators`` can be supplied here
+        and if not ``None`` will override any already set into algorithm so it can be reused with
+        different operators. While an ``operator`` is required by algorithms, ``aux_operators`` are
+        optional.
+
+        Args:
+            operator: Qubit operator of the observable.
+            aux_operators: Optional list of auxiliary operators to be evaluated with the
+                parameters of the minimum eigenvalue main result and their expectation values
+                returned. For instance in chemistry these can be dipole operators and total particle
+                count operators, so we can get values for these at the ground state.
+
+        Returns:
+            A minimum eigensolver result.
+        """
+        return MinimumEigensolverResult()
+
+    @classmethod
+    def supports_aux_operators(cls) -> bool:
+        """Whether computing the expectation value of auxiliary operators is supported.
+
+        If the minimum eigensolver computes an eigenvalue of the main ``operator`` then it can
+        compute the expectation value of the ``aux_operators`` for that state. Otherwise they will
+        be ignored.
+
+        Returns:
+            True if aux_operator expectations can be evaluated, False otherwise
+        """
+        return False
+
+
+class MinimumEigensolverResult(AlgorithmResult):
+    """Minimum eigensolver result."""
+
+    def __init__(self) -> None:
+        super().__init__()
+        self._eigenvalue = None
+        self._aux_operators_evaluated = None
+
+    @property
+    def eigenvalue(self) -> complex | None:
+        """The computed minimum eigenvalue."""
+        return self._eigenvalue
+
+    @eigenvalue.setter
+    def eigenvalue(self, value: complex) -> None:
+        self._eigenvalue = value
+
+    @property
+    def aux_operators_evaluated(self) -> ListOrDict[tuple[complex, dict[str, Any]]] | None:
+        """The aux operator expectation values.
+
+        These values are in fact tuples formatted as (mean, (variance, shots)).
+        """
+        return self._aux_operators_evaluated
+
+    @aux_operators_evaluated.setter
+    def aux_operators_evaluated(self, value: ListOrDict[tuple[complex, dict[str, Any]]]) -> None:
+        self._aux_operators_evaluated = value

qiskit/algorithms/minimum_eigensolvers/numpy_minimum_eigensolver.py

@@ -0,0 +1,106 @@
+# 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 NumPy minimum eigensolver algorithm and result."""
+
+from __future__ import annotations
+
+from typing import Callable, List, Union, Optional
+import logging
+import numpy as np
+
+from qiskit.opflow import PauliSumOp
+from qiskit.quantum_info.operators.base_operator import BaseOperator
+
+# TODO this path will need updating when VQD is merged
+from ..eigen_solvers.numpy_eigen_solver import NumPyEigensolver
+from .minimum_eigensolver import MinimumEigensolver, MinimumEigensolverResult
+from ..list_or_dict import ListOrDict
+
+logger = logging.getLogger(__name__)
+
+# future type annotations not supported in type aliases in 3.8
+FilterType = Callable[[Union[List, np.ndarray], float, Optional[ListOrDict[float]]], bool]
+
+
+class NumPyMinimumEigensolver(MinimumEigensolver):
+    """
+    The NumPy minimum eigensolver algorithm.
+    """
+
+    def __init__(
+        self,
+        filter_criterion: FilterType | None = None,
+    ) -> None:
+        """
+        Args:
+            filter_criterion: callable that allows to filter eigenvalues/eigenstates. The minimum
+                eigensolver is only searching over feasible states and returns an eigenstate that
+                has the smallest eigenvalue among feasible states. The callable has the signature
+                ``filter(eigenstate, eigenvalue, aux_values)`` and must return a boolean to indicate
+                whether to consider this value or not. If there is no feasible element, the result
+                can even be empty.
+        """
+        self._eigensolver = NumPyEigensolver(filter_criterion=filter_criterion)
+
+    @property
+    def filter_criterion(
+        self,
+    ) -> FilterType | None:
+        """Returns the criterion for filtering eigenstates/eigenvalues."""
+        return self._eigensolver.filter_criterion
+
+    @filter_criterion.setter
+    def filter_criterion(
+        self,
+        filter_criterion: FilterType,
+    ) -> None:
+        self._eigensolver.filter_criterion = filter_criterion
+
+    @classmethod
+    def supports_aux_operators(cls) -> bool:
+        return NumPyEigensolver.supports_aux_operators()
+
+    def compute_minimum_eigenvalue(
+        self,
+        operator: BaseOperator | PauliSumOp,
+        aux_operators: ListOrDict[BaseOperator | PauliSumOp] | None = None,
+    ) -> NumPyMinimumEigensolverResult:
+        super().compute_minimum_eigenvalue(operator, aux_operators)
+        eigensolver_result = self._eigensolver.compute_eigenvalues(operator, aux_operators)
+        result = NumPyMinimumEigensolverResult()
+        if eigensolver_result.eigenvalues is not None and len(eigensolver_result.eigenvalues) > 0:
+            result.eigenvalue = eigensolver_result.eigenvalues[0]
+            result.eigenstate = eigensolver_result.eigenstates[0]
+            if eigensolver_result.aux_operator_eigenvalues:
+                result.aux_operators_evaluated = eigensolver_result.aux_operator_eigenvalues[0]
+
+        logger.debug("NumPy minimum eigensolver result: %s", result)
+
+        return result
+
+
+class NumPyMinimumEigensolverResult(MinimumEigensolverResult):
+    """NumPy minimum eigensolver result."""
+
+    def __init__(self) -> None:
+        super().__init__()
+        self._eigenstate = None
+
+    @property
+    def eigenstate(self) -> np.ndarray | None:
+        """Returns the eigenstate corresponding to the computed minimum eigenvalue."""
+        return self._eigenstate
+
+    @eigenstate.setter
+    def eigenstate(self, value: np.ndarray) -> None:
+        self._eigenstate = value

qiskit/algorithms/minimum_eigensolvers/vqe.py

@@ -0,0 +1,341 @@
+# 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 variational quantum eigensolver algorithm."""
+
+from __future__ import annotations
+
+import logging
+from time import time
+from collections.abc import Callable, Sequence
+from typing import Any
+
+import numpy as np
+
+from qiskit.algorithms.gradients import BaseEstimatorGradient
+from qiskit.circuit import QuantumCircuit
+from qiskit.opflow import PauliSumOp
+from qiskit.primitives import BaseEstimator
+from qiskit.quantum_info.operators.base_operator import BaseOperator
+
+from ..exceptions import AlgorithmError
+from ..list_or_dict import ListOrDict
+from ..optimizers import Optimizer, Minimizer, OptimizerResult
+from ..variational_algorithm import VariationalAlgorithm, VariationalResult
+from .minimum_eigensolver import MinimumEigensolver, MinimumEigensolverResult
+from ..observables_evaluator import estimate_observables
+from ..utils import validate_initial_point, validate_bounds
+
+logger = logging.getLogger(__name__)
+
+
+class VQE(VariationalAlgorithm, MinimumEigensolver):
+    r"""The variational quantum eigensolver (VQE) algorithm.
+
+    VQE is a hybrid quantum-classical algorithm that uses a variational technique to find the
+    minimum eigenvalue of a given Hamiltonian operator :math:`H`.
+
+    The ``VQE`` algorithm is executed using an :attr:`estimator` primitive, which computes
+    expectation values of operators (observables).
+
+    An instance of ``VQE`` also requires an :attr:`ansatz`, a parameterized
+    :class:`.QuantumCircuit`, to prepare the trial state :math:`|\psi(\vec\theta)\rangle`. It also
+    needs a classical :attr:`optimizer` which varies the circuit parameters :math:`\vec\theta` such
+    that the expectation value of the operator on the corresponding state approaches a minimum,
+
+    .. math::
+
+        \min_{\vec\theta} \langle\psi(\vec\theta)|H|\psi(\vec\theta)\rangle.
+
+    The :attr:`estimator` is used to compute this expectation value for every optimization step.
+
+    The optimizer can either be one of Qiskit's optimizers, such as
+    :class:`~qiskit.algorithms.optimizers.SPSA` or a callable with the following signature:
+
+    .. code-block:: python
+
+        from qiskit.algorithms.optimizers import OptimizerResult
+
+        def my_minimizer(fun, x0, jac=None, bounds=None) -> OptimizerResult:
+            # Note that the callable *must* have these argument names!
+            # Args:
+            #     fun (callable): the function to minimize
+            #     x0 (np.ndarray): the initial point for the optimization
+            #     jac (callable, optional): the gradient of the objective function
+            #     bounds (list, optional): a list of tuples specifying the parameter bounds
+
+            result = OptimizerResult()
+            result.x = # optimal parameters
+            result.fun = # optimal function value
+            return result
+
+    The above signature also allows one to use any SciPy minimizer, for instance as
+
+    .. code-block:: python
+
+        from functools import partial
+        from scipy.optimize import minimize
+
+        optimizer = partial(minimize, method="L-BFGS-B")
+
+    The following attributes can be set via the initializer but can also be read and updated once
+    the VQE object has been constructed.
+
+    Attributes:
+        estimator (BaseEstimator): The estimator primitive to compute the expectation value of the
+            Hamiltonian operator.
+        ansatz (QuantumCircuit): A parameterized quantum circuit to prepare the trial state.
+        optimizer (Optimizer | Minimizer): A classical optimizer to find the minimum energy. This
+            can either be a Qiskit :class:`.Optimizer` or a callable implementing the
+            :class:`.Minimizer` protocol.
+        gradient (BaseEstimatorGradient | None): An optional estimator gradient to be used with the
+            optimizer.
+        callback (Callable[[int, np.ndarray, float, dict[str, Any]], None] | None): A callback that
+            can access the intermediate data at each optimization step. These data are: the
+            evaluation count, the optimizer parameters for the ansatz, the evaluated mean, and the
+            metadata dictionary.
+
+    References:
+        [1] Peruzzo et al, "A variational eigenvalue solver on a quantum processor"
+            `arXiv:1304.3061 <https://arxiv.org/abs/1304.3061>`__
+    """
+
+    def __init__(
+        self,
+        estimator: BaseEstimator,
+        ansatz: QuantumCircuit,
+        optimizer: Optimizer | Minimizer,
+        *,
+        gradient: BaseEstimatorGradient | None = None,
+        initial_point: Sequence[float] | None = None,
+        callback: Callable[[int, np.ndarray, float, dict[str, Any]], None] | None = None,
+    ) -> None:
+        r"""
+        Args:
+            estimator: The estimator primitive to compute the expectation value of the
+                Hamiltonian operator.
+            ansatz: A parameterized quantum circuit to prepare the trial state.
+            optimizer: A classical optimizer to find the minimum energy. This
+                can either be a Qiskit :class:`.Optimizer` or a callable implementing the
+                :class:`.Minimizer` protocol.
+            gradient: An optional estimator gradient to be used with the optimizer.
+            initial_point: An optional initial point (i.e. initial parameter values) for the
+                optimizer. The length of the initial point must match the number of :attr:`ansatz`
+                parameters. If ``None``, a random point will be generated within certain parameter
+                bounds. ``VQE`` will look to the ansatz for these bounds. If the ansatz does not
+                specify bounds, bounds of :math:`-2\pi`, :math:`2\pi` will be used.
+            callback: A callback that can access the intermediate data at each optimization step.
+                These data are: the evaluation count, the optimizer parameters for the ansatz, the
+                estimated value, and the metadata dictionary.
+        """
+        super().__init__()
+
+        self.estimator = estimator
+        self.ansatz = ansatz
+        self.optimizer = optimizer
+        self.gradient = gradient
+        # this has to go via getters and setters due to the VariationalAlgorithm interface
+        self.initial_point = initial_point
+        self.callback = callback
+
+    @property
+    def initial_point(self) -> Sequence[float] | None:
+        return self._initial_point
+
+    @initial_point.setter
+    def initial_point(self, value: Sequence[float] | None) -> None:
+        self._initial_point = value
+
+    def compute_minimum_eigenvalue(
+        self,
+        operator: BaseOperator | PauliSumOp,
+        aux_operators: ListOrDict[BaseOperator | PauliSumOp] | None = None,
+    ) -> VQEResult:
+        self._check_operator_ansatz(operator)
+
+        initial_point = validate_initial_point(self.initial_point, self.ansatz)
+
+        bounds = validate_bounds(self.ansatz)
+
+        start_time = time()
+
+        evaluate_energy = self._get_evaluate_energy(self.ansatz, operator)
+
+        if self.gradient is not None:
+            evaluate_gradient = self._get_evaluate_gradient(self.ansatz, operator)
+        else:
+            evaluate_gradient = None
+
+        # perform optimization
+        if callable(self.optimizer):
+            optimizer_result = self.optimizer(
+                fun=evaluate_energy, x0=initial_point, jac=evaluate_gradient, bounds=bounds
+            )
+        else:
+            optimizer_result = self.optimizer.minimize(
+                fun=evaluate_energy, x0=initial_point, jac=evaluate_gradient, bounds=bounds
+            )
+
+        optimizer_time = time() - start_time
+
+        logger.info(
+            "Optimization complete in %s seconds.\nFound optimal point %s",
+            optimizer_time,
+            optimizer_result.x,
+        )
+
+        if aux_operators is not None:
+            aux_operators_evaluated = estimate_observables(
+                self.estimator, self.ansatz, aux_operators, optimizer_result.x
+            )
+        else:
+            aux_operators_evaluated = None
+
+        return self._build_vqe_result(optimizer_result, aux_operators_evaluated, optimizer_time)
+
+    @classmethod
+    def supports_aux_operators(cls) -> bool:
+        return True
+
+    def _get_evaluate_energy(
+        self,
+        ansatz: QuantumCircuit,
+        operator: BaseOperator | PauliSumOp,
+    ) -> Callable[[np.ndarray], np.ndarray | float]:
+        """Returns a function handle to evaluate the energy at given parameters for the ansatz.
+        This is the objective function to be passed to the optimizer that is used for evaluation.
+
+        Args:
+            ansatz: The ansatz preparing the quantum state.
+            operator: The operator whose energy to evaluate.
+
+        Returns:
+            A callable that computes and returns the energy of the hamiltonian of each parameter.
+
+        Raises:
+            AlgorithmError: If the primitive job to evaluate the energy fails.
+        """
+        num_parameters = ansatz.num_parameters
+
+        # avoid creating an instance variable to remain stateless regarding results
+        eval_count = 0
+
+        def evaluate_energy(parameters: np.ndarray) -> np.ndarray | float:
+            nonlocal eval_count
+
+            # handle broadcasting: ensure parameters is of shape [array, array, ...]
+            parameters = np.reshape(parameters, (-1, num_parameters)).tolist()
+            batch_size = len(parameters)
+
+            try:
+                job = self.estimator.run(batch_size * [ansatz], batch_size * [operator], parameters)
+                estimator_result = job.result()
+            except Exception as exc:
+                raise AlgorithmError("The primitive job to evaluate the energy failed!") from exc
+
+            values = estimator_result.values
+
+            if self.callback is not None:
+                metadata = estimator_result.metadata
+                for params, value, meta in zip(parameters, values, metadata):
+                    eval_count += 1
+                    self.callback(eval_count, params, value, meta)
+
+            energy = values[0] if len(values) == 1 else values
+
+            return energy
+
+        return evaluate_energy
+
+    def _get_evaluate_gradient(
+        self,
+        ansatz: QuantumCircuit,
+        operator: BaseOperator | PauliSumOp,
+    ) -> Callable[[np.ndarray], np.ndarray]:
+        """Get a function handle to evaluate the gradient at given parameters for the ansatz.
+
+        Args:
+            ansatz: The ansatz preparing the quantum state.
+            operator: The operator whose energy to evaluate.
+
+        Returns:
+            A function handle to evaluate the gradient at given parameters for the ansatz.
+
+        Raises:
+            AlgorithmError: If the primitive job to evaluate the gradient fails.
+        """
+
+        def evaluate_gradient(parameters: np.ndarray) -> np.ndarray:
+            # broadcasting not required for the estimator gradients
+            try:
+                job = self.gradient.run([ansatz], [operator], [parameters])
+                gradients = job.result().gradients
+            except Exception as exc:
+                raise AlgorithmError("The primitive job to evaluate the gradient failed!") from exc
+
+            return gradients[0]
+
+        return evaluate_gradient
+
+    def _check_operator_ansatz(self, operator: BaseOperator | PauliSumOp):
+        """Check that the number of qubits of operator and ansatz match and that the ansatz is
+        parameterized.
+        """
+        if operator.num_qubits != self.ansatz.num_qubits:
+            try:
+                logger.info(
+                    "Trying to resize ansatz to match operator on %s qubits.", operator.num_qubits
+                )
+                self.ansatz.num_qubits = operator.num_qubits
+            except AttributeError as error:
+                raise AlgorithmError(
+                    "The number of qubits of the ansatz does not match the "
+                    "operator, and the ansatz does not allow setting the "
+                    "number of qubits using `num_qubits`."
+                ) from error
+
+        if self.ansatz.num_parameters == 0:
+            raise AlgorithmError("The ansatz must be parameterized, but has no free parameters.")
+
+    def _build_vqe_result(
+        self,
+        optimizer_result: OptimizerResult,
+        aux_operators_evaluated: ListOrDict[tuple[complex, tuple[complex, int]]],
+        optimizer_time: float,
+    ) -> VQEResult:
+        result = VQEResult()
+        result.eigenvalue = optimizer_result.fun
+        result.cost_function_evals = optimizer_result.nfev
+        result.optimal_point = optimizer_result.x
+        result.optimal_parameters = dict(zip(self.ansatz.parameters, optimizer_result.x))
+        result.optimal_value = optimizer_result.fun
+        result.optimizer_time = optimizer_time
+        result.aux_operators_evaluated = aux_operators_evaluated
+        result.optimizer_result = optimizer_result
+        return result
+
+
+class VQEResult(VariationalResult, MinimumEigensolverResult):
+    """Variational quantum eigensolver result."""
+
+    def __init__(self) -> None:
+        super().__init__()
+        self._cost_function_evals = None
+
+    @property
+    def cost_function_evals(self) -> int | None:
+        """The number of cost optimizer evaluations."""
+        return self._cost_function_evals
+
+    @cost_function_evals.setter
+    def cost_function_evals(self, value: int) -> None:
+        self._cost_function_evals = value

qiskit/algorithms/observables_evaluator.py

@@ -9,8 +9,12 @@
 # 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.
+
 """Evaluator of observables for algorithms."""
+
 from __future__ import annotations
+from collections.abc import Sequence
+from typing import Any
 
 import numpy as np
 
@@ -18,7 +22,7 @@ from qiskit import QuantumCircuit
 from qiskit.opflow import PauliSumOp
 from .exceptions import AlgorithmError
 from .list_or_dict import ListOrDict
-from ..primitives import EstimatorResult, BaseEstimator
+from ..primitives import BaseEstimator
 from ..quantum_info.operators.base_operator import BaseOperator
 
 
@@ -26,38 +30,31 @@ def estimate_observables(
     estimator: BaseEstimator,
     quantum_state: QuantumCircuit,
     observables: ListOrDict[BaseOperator | PauliSumOp],
+    parameter_values: Sequence[float] | None = None,
     threshold: float = 1e-12,
-) -> ListOrDict[tuple[complex, tuple[complex, int]]]:
+) -> ListOrDict[tuple[complex, dict[str, Any]]]:
     """
     Accepts a sequence of operators and calculates their expectation values - means
-    and standard deviations. They are calculated with respect to a quantum state provided. A user
+    and metadata. They are calculated with respect to a quantum state provided. A user
     can optionally provide a threshold value which filters mean values falling below the threshold.
 
     Args:
         estimator: An estimator primitive used for calculations.
-        quantum_state: An unparametrized quantum circuit representing a quantum state that
-            expectation values are computed against.
+        quantum_state: A (parameterized) quantum circuit preparing a quantum state that expectation
+            values are computed against.
         observables: A list or a dictionary of operators whose expectation values are to be
             calculated.
+        parameter_values: Optional list of parameters values to evaluate the quantum circuit on.
         threshold: A threshold value that defines which mean values should be neglected (helpful for
             ignoring numerical instabilities close to 0).
 
     Returns:
-        A list or a dictionary of tuples (mean, (variance, shots)).
+        A list or a dictionary of tuples (mean, metadata).
 
     Raises:
-        ValueError: If a ``quantum_state`` with free parameters is provided.
         AlgorithmError: If a primitive job is not successful.
     """
 
-    if (
-        isinstance(quantum_state, QuantumCircuit)  # State cannot be parametrized
-        and len(quantum_state.parameters) > 0
-    ):
-        raise ValueError(
-            "A parametrized representation of a quantum_state was provided. It is not "
-            "allowed - it cannot have free parameters."
-        )
     if isinstance(observables, dict):
         observables_list = list(observables.values())
     else:
@@ -65,18 +62,19 @@ def estimate_observables(
 
     observables_list = _handle_zero_ops(observables_list)
     quantum_state = [quantum_state] * len(observables)
+    if parameter_values is not None:
+        parameter_values = [parameter_values] * len(observables)
     try:
-        estimator_job = estimator.run(quantum_state, observables_list)
+        estimator_job = estimator.run(quantum_state, observables_list, parameter_values)
         expectation_values = estimator_job.result().values
     except Exception as exc:
         raise AlgorithmError("The primitive job failed!") from exc
 
-    variance_and_shots = _prep_variance_and_shots(estimator_job, len(expectation_values))
-
+    metadata = estimator_job.result().metadata
     # Discard values below threshold
     observables_means = expectation_values * (np.abs(expectation_values) > threshold)
-    # zip means and standard deviations into tuples
-    observables_results = list(zip(observables_means, variance_and_shots))
+    # zip means and metadata into tuples
+    observables_results = list(zip(observables_means, metadata))
 
     return _prepare_result(observables_results, observables)
 
@@ -95,20 +93,20 @@ def _handle_zero_ops(
 
 
 def _prepare_result(
-    observables_results: list[tuple[complex, tuple[complex, int]]],
+    observables_results: list[tuple[complex, dict]],
     observables: ListOrDict[BaseOperator | PauliSumOp],
-) -> ListOrDict[tuple[complex, tuple[complex, int]]]:
+) -> ListOrDict[tuple[complex, dict[str, Any]]]:
     """
-    Prepares a list of tuples of eigenvalues and (variance, shots) tuples from
+    Prepares a list of tuples of eigenvalues and metadata tuples from
     ``observables_results`` and ``observables``.
 
     Args:
-        observables_results: A list of tuples (mean, (variance, shots)).
+        observables_results: A list of tuples (mean, metadata).
         observables: A list or a dictionary of operators whose expectation values are to be
             calculated.
 
     Returns:
-        A list or a dictionary of tuples (mean, (variance, shots)).
+        A list or a dictionary of tuples (mean, metadata).
     """
 
     if isinstance(observables, list):
@@ -122,36 +120,3 @@ def _prepare_result(
     for key, value in key_value_iterator:
         observables_eigenvalues[key] = value
     return observables_eigenvalues
-
-
-def _prep_variance_and_shots(
-    estimator_result: EstimatorResult,
-    results_length: int,
-) -> list[tuple[complex, int]]:
-    """
-    Prepares a list of tuples with variances and shots from results provided by expectation values
-    calculations. If there is no variance or shots data available from a primitive, the values will
-    be set to ``0``.
-
-    Args:
-        estimator_result: An estimator result.
-        results_length: Number of expectation values calculated.
-
-    Returns:
-        A list of tuples of the form (variance, shots).
-    """
-    if not estimator_result.metadata:
-        return [(0, 0)] * results_length
-
-    results = []
-    for metadata in estimator_result.metadata:
-        variance, shots = 0.0, 0
-        if metadata:
-            if "variance" in metadata.keys():
-                variance = metadata["variance"]
-            if "shots" in metadata.keys():
-                shots = metadata["shots"]
-
-        results.append((variance, shots))
-
-    return results

qiskit/algorithms/optimizers/qnspsa.py

@@ -187,7 +187,9 @@ class QNSPSA(SPSA):
         self._nfev += 6
 
         loss_values = _batch_evaluate(loss, loss_points, self._max_evals_grouped)
-        fidelity_values = _batch_evaluate(self.fidelity, fidelity_points, self._max_evals_grouped)
+        fidelity_values = _batch_evaluate(
+            self.fidelity, fidelity_points, self._max_evals_grouped, unpack_points=True
+        )
 
         # compute the gradient approximation and additionally return the loss function evaluations
         gradient_estimate = (loss_values[0] - loss_values[1]) / (2 * eps) * delta1
@@ -268,11 +270,13 @@ class QNSPSA(SPSA):
             sampler = CircuitSampler(backend)
 
             def fidelity(values_x, values_y=None):
-                if values_y is not None:  # no batches
+                # no batches
+                if isinstance(values_x, np.ndarray) and isinstance(values_y, np.ndarray):
                     value_dict = dict(
                         zip(params_x[:] + params_y[:], values_x.tolist() + values_y.tolist())
                     )
-                else:
+                # legacy batching -- remove once QNSPSA.get_fidelity is only supported with sampler
+                elif values_y is None:
                     value_dict = {p: [] for p in params_x[:] + params_y[:]}
                     for values_xy in values_x:
                         for value_x, param_x in zip(values_xy[0, :], params_x):
@@ -280,6 +284,14 @@ class QNSPSA(SPSA):
 
                         for value_y, param_y in zip(values_xy[1, :], params_y):
                             value_dict[param_y].append(value_y)
+                else:
+                    value_dict = {p: [] for p in params_x[:] + params_y[:]}
+                    for values_i_x, values_i_y in zip(values_x, values_y):
+                        for value_x, param_x in zip(values_i_x, params_x):
+                            value_dict[param_x].append(value_x)
+
+                        for value_y, param_y in zip(values_i_y, params_y):
+                            value_dict[param_y].append(value_y)
 
                 return np.abs(sampler.convert(expression, params=value_dict).eval()) ** 2
 

qiskit/algorithms/optimizers/spsa.py

@@ -711,7 +711,7 @@ def constant(eta=0.01):
         yield eta
 
 
-def _batch_evaluate(function, points, max_evals_grouped):
+def _batch_evaluate(function, points, max_evals_grouped, unpack_points=False):
     # if the function cannot handle lists of points as input, cover this case immediately
     if max_evals_grouped == 1:
         # support functions with multiple arguments where the points are given in a tuple
@@ -731,11 +731,32 @@ def _batch_evaluate(function, points, max_evals_grouped):
 
     results = []
     for batch in batched_points:
-        results += function(batch).tolist()
+        if unpack_points:
+            batch = _repack_points(batch)
+            results += _as_list(function(*batch))
+        else:
+            results += _as_list(function(batch))
 
     return results
 
 
+def _as_list(obj):
+    """Convert a list or numpy array into a list."""
+    return obj.tolist() if isinstance(obj, np.ndarray) else obj
+
+
+def _repack_points(points):
+    """Turn a list of tuples of points into a tuple of lists of points.
+    E.g. turns
+        [(a1, a2, a3), (b1, b2, b3)]
+    into
+        ([a1, b1], [a2, b2], [a3, b3])
+    where all elements are np.ndarray.
+    """
+    num_sets = len(points[0])  # length of (a1, a2, a3)
+    return ([x[i] for x in points] for i in range(num_sets))
+
+
 def _make_spd(matrix, bias=0.01):
     identity = np.identity(matrix.shape[0])
     psd = scipy.linalg.sqrtm(matrix.dot(matrix))

qiskit/algorithms/time_evolvers/trotterization/trotter_qrte.py

@@ -168,6 +168,7 @@ class TrotterQRTE(RealTimeEvolver):
                 self.estimator,
                 evolved_state,
                 evolution_problem.aux_operators,
+                None,
                 evolution_problem.truncation_threshold,
             )
 

qiskit/algorithms/utils/__init__.py

@@ -0,0 +1,21 @@
+# 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.
+
+"""Common Qiskit algorithms utility functions."""
+
+from .validate_initial_point import validate_initial_point
+from .validate_bounds import validate_bounds
+
+__all__ = [
+    "validate_initial_point",
+    "validate_bounds",
+]

qiskit/algorithms/utils/validate_bounds.py

@@ -0,0 +1,44 @@
+# 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.
+
+"""Validate parameter bounds."""
+
+from __future__ import annotations
+
+from qiskit.circuit import QuantumCircuit
+
+
+def validate_bounds(circuit: QuantumCircuit) -> list[tuple(float | None, float | None)]:
+    """
+    Validate the bounds provided by a quantum circuit against its number of parameters.
+    If no bounds are obtained, return ``None`` for all lower and upper bounds.
+
+    Args:
+        circuit: A parameterized quantum circuit.
+
+    Returns:
+        A list of tuples (lower_bound, upper_bound)).
+
+    Raises:
+        ValueError: If the number of bounds does not the match the number of circuit parameters.
+    """
+    if hasattr(circuit, "parameter_bounds") and circuit.parameter_bounds is not None:
+        bounds = circuit.parameter_bounds
+        if len(bounds) != circuit.num_parameters:
+            raise ValueError(
+                f"The number of bounds ({len(bounds)}) does not match the number of "
+                f"parameters in the circuit ({circuit.num_parameters})."
+            )
+    else:
+        bounds = [(None, None)] * circuit.num_parameters
+
+    return bounds

qiskit/algorithms/utils/validate_initial_point.py

@@ -0,0 +1,69 @@
+# 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.
+
+"""Validate an initial point."""
+
+from __future__ import annotations
+
+from collections.abc import Sequence
+
+import numpy as np
+
+from qiskit.circuit import QuantumCircuit
+from qiskit.utils import algorithm_globals
+
+
+def validate_initial_point(
+    point: Sequence[float] | None, circuit: QuantumCircuit
+) -> Sequence[float]:
+    r"""
+    Validate a choice of initial point against a choice of circuit. If no point is provided, a
+    random point will be generated within certain parameter bounds. It will first look to the
+    circuit for these bounds. If the circuit does not specify bounds, bounds of :math:`-2\pi`,
+    :math:`2\pi` will be used.
+
+    Args:
+        point: An initial point.
+        circuit: A parameterized quantum circuit.
+
+    Returns:
+        A validated initial point.
+
+    Raises:
+        ValueError: If the dimension of the initial point does not match the number of circuit
+        parameters.
+    """
+    expected_size = circuit.num_parameters
+
+    if point is None:
+        # get bounds if circuit has them set, otherwise use [-2pi, 2pi] for each parameter
+        bounds = getattr(circuit, "parameter_bounds", None)
+        if bounds is None:
+            bounds = [(-2 * np.pi, 2 * np.pi)] * expected_size
+
+        # replace all Nones by [-2pi, 2pi]
+        lower_bounds = []
+        upper_bounds = []
+        for lower, upper in bounds:
+            lower_bounds.append(lower if lower is not None else -2 * np.pi)
+            upper_bounds.append(upper if upper is not None else 2 * np.pi)
+
+        # sample from within bounds
+        point = algorithm_globals.random.uniform(lower_bounds, upper_bounds)
+
+    elif len(point) != expected_size:
+        raise ValueError(
+            f"The dimension of the initial point ({len(point)}) does not match the "
+            f"number of parameters in the circuit ({expected_size})."
+        )
+
+    return point

qiskit/algorithms/variational_algorithm.py

@@ -31,6 +31,7 @@ from abc import ABC, abstractmethod
 import numpy as np
 
 from .algorithm_result import AlgorithmResult
+from .optimizers import OptimizerResult
 
 
 class VariationalAlgorithm(ABC):
@@ -109,3 +110,13 @@ class VariationalResult(AlgorithmResult):
     def optimal_parameters(self, value: Dict) -> None:
         """Sets optimal parameters"""
         self._optimal_parameters = value
+
+    @property
+    def optimizer_result(self) -> Optional[OptimizerResult]:
+        """Returns the optimizer result"""
+        return self._optimizer_result
+
+    @optimizer_result.setter
+    def optimizer_result(self, value: OptimizerResult) -> None:
+        """Sets optimizer result"""
+        self._optimizer_result = value

releasenotes/notes/0.19/add-getters-and-setters-for-vqe-edc753591b368980.yaml

@@ -3,8 +3,10 @@ features:
   - |
     Every attribute of the :class:`~qiskit.algorithms.VQE` class that is set at
     the initialization is now accessible with getters and setters. Further, the
-    default values of the VQE attributes :attr:`~.VQE.ansatz` and
-    :attr:`~.VQE.optimizer` can be reset by assigning ``None`` to them::
+    default values of the VQE attributes 
+    :attr:`~qiskit.algorithms.minimimum_eigen_solvers.VQE.ansatz` and
+    :attr:`~qiskit.algorithms.minimimum_eigen_solvers.VQE.optimizer` can be 
+    reset by assigning ``None`` to them::
     
         vqe = VQE(my_ansatz, my_optimizer)
         vqe.ansatz = None   # reset to default: RealAmplitudes ansatz

releasenotes/notes/0.19/support-dict-for-aux-operators-c3c9ad380c208afd.yaml

@@ -1,9 +1,11 @@
 ---
 features:
   - |
-    The :obj:`.Eigensolver` and :obj:`.MinimumEigensolver` interfaces now support the type
+    The :obj:`.Eigensolver` and :obj:`~qiskit.algorithms.minimimum_eigen_solvers.MinimumEigensolver`
+    interfaces now support the type
     ``Dict[str, Optional[OperatorBase]]`` for the ``aux_operators`` parameter in their respective
-    :meth:`~.Eigensolver.compute_eigenvalues` and :meth:`~.MinimumEigensolver.compute_minimum_eigenvalue` methods.
+    :meth:`~.Eigensolver.compute_eigenvalues` and 
+    :meth:`~qiskit.algorithms.minimimum_eigen_solvers.MinimumEigensolver.compute_minimum_eigenvalue` methods.
     In this case, the auxiliary eigenvalues are also stored in a dictionary under the same keys
     provided by the ``aux_operators`` dictionary. Keys that correspond to an operator that does not commute
     with the main operator are dropped.

releasenotes/notes/0.19/taper_empty_operator_fix-53ce20e5d2b68fd6.yaml

@@ -1,11 +1,11 @@
 ---
 fixes:
   - |
-    When tapering an empty zero operator in :mod:`qiskit.opflow`, the code, on detecting it was zero, logged a
-    warning and returned the original operator.  Such operators are commonly found in
-    the auxiliary operators, when using Qiskit Nature, and the above behavior caused :obj:`.VQE`
-    to throw an exception as tapered non-zero operators were a different number of qubits
-    from the tapered zero operators (since taper has returned the input operator unchanged).
-    The code will now correctly taper a zero operator such that the number of qubits is
-    reduced as expected and matches to tapered non-zero operators e.g ```0*"IIII"``` when we are
-    tapering by 3 qubits will become ``0*"I"``.
+    When tapering an empty zero operator in :mod:`qiskit.opflow`, the code, on detecting it was
+    zero, logged a warning and returned the original operator.  Such operators are commonly found in
+    the auxiliary operators, when using Qiskit Nature, and the above behavior caused
+    :obj:`~qiskit.algorithms.minimimum_eigen_solvers.VQE` to throw an exception as tapered non-zero
+    operators were a different number of qubits from the tapered zero operators (since taper has
+    returned the input operator unchanged). The code will now correctly taper a zero operator such
+    that the number of qubits is reduced as expected and matches to tapered non-zero operators e.g
+    ```0*"IIII"``` when we are tapering by 3 qubits will become ``0*"I"``.

releasenotes/notes/vqe-with-estimator-primitive-7cbcc462ad4dc593.yaml

@@ -0,0 +1,34 @@
+---
+features:
+  - |
+    Added the :class:`qiskit.algorithms.minimum_eigensolvers` package to include interfaces for
+    primitive-enabled algorithms. :class:`~qiskit.algorithms.minimum_eigensolvers.VQE` has been
+    refactored in this implementation to leverage primitives.
+
+    To use the new implementation with a reference primitive, one can do, for example:
+
+    .. code-block:: python
+
+        from qiskit.algorithms.minimum_eigensolvers import VQE
+        from qiskit.algorithms.optimizers import SLSQP
+        from qiskit.circuit.library import TwoLocal
+        from qiskit.primitives import Estimator
+
+        h2_op = SparsePauliOp(
+            ["II", "IZ", "ZI", "ZZ", "XX"],
+            coeffs=[
+                -1.052373245772859,
+                0.39793742484318045,
+                -0.39793742484318045,
+                -0.01128010425623538,
+                0.18093119978423156,
+            ],
+        )
+
+        estimator = Estimator()
+        ansatz = TwoLocal(rotation_blocks=["ry", "rz"], entanglement_blocks="cz")
+        optimizer = SLSQP()
+
+        vqe = VQE(estimator, ansatz, optimizer)
+        result = vqe.compute_minimum_eigenvalue(h2_op)
+        eigenvalue = result.eigenvalue        

test/python/algorithms/minimum_eigensolvers/__init__.py

@@ -0,0 +1,11 @@
+# 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.

test/python/algorithms/minimum_eigensolvers/test_numpy_minimum_eigensolver.py

@@ -0,0 +1,241 @@
+# 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.
+
+""" Test NumPy minimum eigensolver """
+
+import unittest
+from test.python.algorithms import QiskitAlgorithmsTestCase
+
+import numpy as np
+from ddt import ddt, data
+
+from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver
+from qiskit.opflow import PauliSumOp
+
+
+@ddt
+class TestNumPyMinimumEigensolver(QiskitAlgorithmsTestCase):
+    """Test NumPy minimum eigensolver"""
+
+    def setUp(self):
+        super().setUp()
+        self.qubit_op = PauliSumOp.from_list(
+            [
+                ("II", -1.052373245772859),
+                ("ZI", 0.39793742484318045),
+                ("IZ", -0.39793742484318045),
+                ("ZZ", -0.01128010425623538),
+                ("XX", 0.18093119978423156),
+            ]
+        )
+
+        aux_op1 = PauliSumOp.from_list([("II", 2.0)])
+        aux_op2 = PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
+        self.aux_ops_list = [aux_op1, aux_op2]
+        self.aux_ops_dict = {"aux_op1": aux_op1, "aux_op2": aux_op2}
+
+    def test_cme(self):
+        """Basic test"""
+        algo = NumPyMinimumEigensolver()
+        result = algo.compute_minimum_eigenvalue(
+            operator=self.qubit_op, aux_operators=self.aux_ops_list
+        )
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+        self.assertEqual(len(result.aux_operators_evaluated), 2)
+        np.testing.assert_array_almost_equal(result.aux_operators_evaluated[0], [2, 0])
+        np.testing.assert_array_almost_equal(result.aux_operators_evaluated[1], [0, 0])
+
+    def test_cme_reuse(self):
+        """Test reuse"""
+        algo = NumPyMinimumEigensolver()
+
+        with self.subTest("Test with no operator or aux_operators, give via compute method"):
+            result = algo.compute_minimum_eigenvalue(operator=self.qubit_op)
+            self.assertEqual(result.eigenvalue.dtype, np.float64)
+            self.assertAlmostEqual(result.eigenvalue, -1.85727503)
+            self.assertIsNone(result.aux_operators_evaluated)
+
+        with self.subTest("Test with added aux_operators"):
+            result = algo.compute_minimum_eigenvalue(
+                operator=self.qubit_op, aux_operators=self.aux_ops_list
+            )
+            self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+            self.assertEqual(len(result.aux_operators_evaluated), 2)
+            np.testing.assert_array_almost_equal(result.aux_operators_evaluated[0], [2, 0])
+            np.testing.assert_array_almost_equal(result.aux_operators_evaluated[1], [0, 0])
+
+        with self.subTest("Test with aux_operators removed"):
+            result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=[])
+            self.assertEqual(result.eigenvalue.dtype, np.float64)
+            self.assertAlmostEqual(result.eigenvalue, -1.85727503)
+            self.assertIsNone(result.aux_operators_evaluated)
+
+        with self.subTest("Test with aux_operators set again"):
+            result = algo.compute_minimum_eigenvalue(
+                operator=self.qubit_op, aux_operators=self.aux_ops_list
+            )
+            self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+            self.assertEqual(len(result.aux_operators_evaluated), 2)
+            np.testing.assert_array_almost_equal(result.aux_operators_evaluated[0], [2, 0])
+            np.testing.assert_array_almost_equal(result.aux_operators_evaluated[1], [0, 0])
+
+        with self.subTest("Test after setting first aux_operators as main operator"):
+            result = algo.compute_minimum_eigenvalue(
+                operator=self.aux_ops_list[0], aux_operators=[]
+            )
+            self.assertAlmostEqual(result.eigenvalue, 2 + 0j)
+            self.assertIsNone(result.aux_operators_evaluated)
+
+    def test_cme_filter(self):
+        """Basic test"""
+
+        # define filter criterion
+        # pylint: disable=unused-argument
+        def criterion(x, v, a_v):
+            return v >= -0.5
+
+        algo = NumPyMinimumEigensolver(filter_criterion=criterion)
+        result = algo.compute_minimum_eigenvalue(
+            operator=self.qubit_op, aux_operators=self.aux_ops_list
+        )
+        self.assertAlmostEqual(result.eigenvalue, -0.22491125 + 0j)
+        self.assertEqual(len(result.aux_operators_evaluated), 2)
+        np.testing.assert_array_almost_equal(result.aux_operators_evaluated[0], [2, 0])
+        np.testing.assert_array_almost_equal(result.aux_operators_evaluated[1], [0, 0])
+
+    def test_cme_filter_empty(self):
+        """Test with filter always returning False"""
+
+        # define filter criterion
+        # pylint: disable=unused-argument
+        def criterion(x, v, a_v):
+            return False
+
+        algo = NumPyMinimumEigensolver(filter_criterion=criterion)
+        result = algo.compute_minimum_eigenvalue(
+            operator=self.qubit_op, aux_operators=self.aux_ops_list
+        )
+        self.assertEqual(result.eigenvalue, None)
+        self.assertEqual(result.eigenstate, None)
+        self.assertEqual(result.aux_operators_evaluated, None)
+
+    @data("X", "Y", "Z")
+    def test_cme_1q(self, op):
+        """Test for 1 qubit operator"""
+        algo = NumPyMinimumEigensolver()
+        operator = PauliSumOp.from_list([(op, 1.0)])
+        result = algo.compute_minimum_eigenvalue(operator=operator)
+        self.assertAlmostEqual(result.eigenvalue, -1)
+
+    def test_cme_aux_ops_dict(self):
+        """Test dictionary compatibility of aux_operators"""
+        # Start with an empty dictionary
+        algo = NumPyMinimumEigensolver()
+
+        with self.subTest("Test with an empty dictionary."):
+            result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators={})
+            self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+            self.assertIsNone(result.aux_operators_evaluated)
+
+        with self.subTest("Test with two auxiliary operators."):
+            result = algo.compute_minimum_eigenvalue(
+                operator=self.qubit_op, aux_operators=self.aux_ops_dict
+            )
+            self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+            self.assertEqual(len(result.aux_operators_evaluated), 2)
+            np.testing.assert_array_almost_equal(result.aux_operators_evaluated["aux_op1"], [2, 0])
+            np.testing.assert_array_almost_equal(result.aux_operators_evaluated["aux_op2"], [0, 0])
+
+        with self.subTest("Test with additional zero and None operators."):
+            extra_ops = {"None_op": None, "zero_op": 0, **self.aux_ops_dict}
+            result = algo.compute_minimum_eigenvalue(
+                operator=self.qubit_op, aux_operators=extra_ops
+            )
+            self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+            self.assertEqual(len(result.aux_operators_evaluated), 3)
+            np.testing.assert_array_almost_equal(result.aux_operators_evaluated["aux_op1"], [2, 0])
+            np.testing.assert_array_almost_equal(result.aux_operators_evaluated["aux_op2"], [0, 0])
+            self.assertEqual(result.aux_operators_evaluated["zero_op"], (0.0, 0))
+
+    def test_aux_operators_list(self):
+        """Test list-based aux_operators."""
+        aux_op1 = PauliSumOp.from_list([("II", 2.0)])
+        aux_op2 = PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
+        aux_ops = [aux_op1, aux_op2]
+        algo = NumPyMinimumEigensolver()
+
+        with self.subTest("Test with two auxiliary operators."):
+            result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=aux_ops)
+            self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+            self.assertEqual(len(result.aux_operators_evaluated), 2)
+            # expectation values
+            self.assertAlmostEqual(result.aux_operators_evaluated[0][0], 2, places=6)
+            self.assertAlmostEqual(result.aux_operators_evaluated[1][0], 0, places=6)
+            # standard deviations
+            self.assertAlmostEqual(result.aux_operators_evaluated[0][1], 0.0)
+            self.assertAlmostEqual(result.aux_operators_evaluated[1][1], 0.0)
+
+        with self.subTest("Test with additional zero and None operators."):
+            extra_ops = [*aux_ops, None, 0]
+            result = algo.compute_minimum_eigenvalue(
+                operator=self.qubit_op, aux_operators=extra_ops
+            )
+            self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+            self.assertEqual(len(result.aux_operators_evaluated), 4)
+            # expectation values
+            self.assertAlmostEqual(result.aux_operators_evaluated[0][0], 2, places=6)
+            self.assertAlmostEqual(result.aux_operators_evaluated[1][0], 0, places=6)
+            self.assertIsNone(result.aux_operators_evaluated[2], None)
+            self.assertEqual(result.aux_operators_evaluated[3][0], 0.0)
+            # standard deviations
+            self.assertAlmostEqual(result.aux_operators_evaluated[0][1], 0.0)
+            self.assertAlmostEqual(result.aux_operators_evaluated[1][1], 0.0)
+            self.assertEqual(result.aux_operators_evaluated[3][1], 0.0)
+
+    def test_aux_operators_dict(self):
+        """Test dict-based aux_operators."""
+        aux_op1 = PauliSumOp.from_list([("II", 2.0)])
+        aux_op2 = PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
+        aux_ops = {"aux_op1": aux_op1, "aux_op2": aux_op2}
+        algo = NumPyMinimumEigensolver()
+
+        with self.subTest("Test with two auxiliary operators."):
+            result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=aux_ops)
+            self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+            self.assertEqual(len(result.aux_operators_evaluated), 2)
+            # expectation values
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op1"][0], 2, places=6)
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op2"][0], 0, places=6)
+            # standard deviations
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op1"][1], 0.0)
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op2"][1], 0.0)
+
+        with self.subTest("Test with additional zero and None operators."):
+            extra_ops = {**aux_ops, "None_operator": None, "zero_operator": 0}
+            result = algo.compute_minimum_eigenvalue(
+                operator=self.qubit_op, aux_operators=extra_ops
+            )
+            self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+            self.assertEqual(len(result.aux_operators_evaluated), 3)
+            # expectation values
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op1"][0], 2, places=6)
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op2"][0], 0, places=6)
+            self.assertEqual(result.aux_operators_evaluated["zero_operator"][0], 0.0)
+            self.assertTrue("None_operator" not in result.aux_operators_evaluated.keys())
+            # standard deviations
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op1"][1], 0.0)
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op2"][1], 0.0)
+            self.assertAlmostEqual(result.aux_operators_evaluated["zero_operator"][1], 0.0)
+
+
+if __name__ == "__main__":
+    unittest.main()

test/python/algorithms/minimum_eigensolvers/test_vqe.py

@@ -0,0 +1,444 @@
+# 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.
+
+"""Test the variational quantum eigensolver algorithm."""
+
+import unittest
+from test.python.algorithms import QiskitAlgorithmsTestCase
+
+from functools import partial
+import numpy as np
+from scipy.optimize import minimize as scipy_minimize
+from ddt import data, ddt
+
+from qiskit import QuantumCircuit
+from qiskit.algorithms import AlgorithmError
+from qiskit.algorithms.gradients import ParamShiftEstimatorGradient
+from qiskit.algorithms.minimum_eigensolvers import VQE
+from qiskit.algorithms.optimizers import (
+    CG,
+    COBYLA,
+    GradientDescent,
+    L_BFGS_B,
+    OptimizerResult,
+    P_BFGS,
+    QNSPSA,
+    SLSQP,
+    SPSA,
+    TNC,
+)
+from qiskit.algorithms.state_fidelities import ComputeUncompute
+from qiskit.circuit.library import RealAmplitudes, TwoLocal
+from qiskit.opflow import PauliSumOp, TwoQubitReduction
+from qiskit.quantum_info import SparsePauliOp, Operator, Pauli
+from qiskit.primitives import Estimator, Sampler
+from qiskit.utils import algorithm_globals
+
+# pylint: disable=invalid-name
+def _mock_optimizer(fun, x0, jac=None, bounds=None, inputs=None) -> OptimizerResult:
+    """A mock of a callable that can be used as minimizer in the VQE."""
+    result = OptimizerResult()
+    result.x = np.zeros_like(x0)
+    result.fun = fun(result.x)
+    result.nit = 0
+
+    if inputs is not None:
+        inputs.update({"fun": fun, "x0": x0, "jac": jac, "bounds": bounds})
+    return result
+
+
+@ddt
+class TestVQE(QiskitAlgorithmsTestCase):
+    """Test VQE"""
+
+    def setUp(self):
+        super().setUp()
+        self.seed = 50
+        algorithm_globals.random_seed = self.seed
+        self.h2_op = SparsePauliOp(
+            ["II", "IZ", "ZI", "ZZ", "XX"],
+            coeffs=[
+                -1.052373245772859,
+                0.39793742484318045,
+                -0.39793742484318045,
+                -0.01128010425623538,
+                0.18093119978423156,
+            ],
+        )
+        self.h2_energy = -1.85727503
+
+        self.ryrz_wavefunction = TwoLocal(rotation_blocks=["ry", "rz"], entanglement_blocks="cz")
+        self.ry_wavefunction = TwoLocal(rotation_blocks="ry", entanglement_blocks="cz")
+
+    @data(L_BFGS_B(), COBYLA())
+    def test_basic_aer_statevector(self, estimator):
+        """Test VQE using reference Estimator."""
+        vqe = VQE(Estimator(), self.ryrz_wavefunction, estimator)
+
+        result = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
+
+        with self.subTest(msg="test eigenvalue"):
+            self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=5)
+
+        with self.subTest(msg="test optimal_value"):
+            self.assertAlmostEqual(result.optimal_value, self.h2_energy)
+
+        with self.subTest(msg="test dimension of optimal point"):
+            self.assertEqual(len(result.optimal_point), 16)
+
+        with self.subTest(msg="assert cost_function_evals is set"):
+            self.assertIsNotNone(result.cost_function_evals)
+
+        with self.subTest(msg="assert optimizer_time is set"):
+            self.assertIsNotNone(result.optimizer_time)
+
+        with self.subTest(msg="assert optimizer_result is set"):
+            self.assertIsNotNone(result.optimizer_result)
+
+        with self.subTest(msg="assert optimizer_result."):
+            self.assertAlmostEqual(result.optimizer_result.fun, self.h2_energy, places=5)
+
+    def test_invalid_initial_point(self):
+        """Test the proper error is raised when the initial point has the wrong size."""
+        ansatz = self.ryrz_wavefunction
+        initial_point = np.array([1])
+
+        vqe = VQE(
+            Estimator(),
+            ansatz,
+            SLSQP(),
+            initial_point=initial_point,
+        )
+
+        with self.assertRaises(ValueError):
+            _ = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
+
+    def test_ansatz_resize(self):
+        """Test the ansatz is properly resized if it's a blueprint circuit."""
+        ansatz = RealAmplitudes(1, reps=1)
+        vqe = VQE(Estimator(), ansatz, SLSQP())
+        result = vqe.compute_minimum_eigenvalue(self.h2_op)
+        self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=5)
+
+    def test_invalid_ansatz_size(self):
+        """Test an error is raised if the ansatz has the wrong number of qubits."""
+        ansatz = QuantumCircuit(1)
+        ansatz.compose(RealAmplitudes(1, reps=2))
+        vqe = VQE(Estimator(), ansatz, SLSQP())
+
+        with self.assertRaises(AlgorithmError):
+            _ = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
+
+    def test_missing_ansatz_params(self):
+        """Test specifying an ansatz with no parameters raises an error."""
+        ansatz = QuantumCircuit(self.h2_op.num_qubits)
+        vqe = VQE(Estimator(), ansatz, SLSQP())
+        with self.assertRaises(AlgorithmError):
+            vqe.compute_minimum_eigenvalue(operator=self.h2_op)
+
+    def test_max_evals_grouped(self):
+        """Test with SLSQP with max_evals_grouped."""
+        optimizer = SLSQP(maxiter=50, max_evals_grouped=5)
+        vqe = VQE(
+            Estimator(),
+            self.ryrz_wavefunction,
+            optimizer,
+        )
+        result = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
+        self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=5)
+
+    @data(
+        CG(),
+        L_BFGS_B(),
+        P_BFGS(),
+        SLSQP(),
+        TNC(),
+    )
+    def test_with_gradient(self, optimizer):
+        """Test VQE using gradient primitive."""
+        estimator = Estimator()
+        vqe = VQE(
+            estimator,
+            self.ry_wavefunction,
+            optimizer,
+            gradient=ParamShiftEstimatorGradient(estimator),
+        )
+        result = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
+        self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=5)
+
+    def test_gradient_passed(self):
+        """Test the gradient is properly passed into the optimizer."""
+        inputs = {}
+        estimator = Estimator()
+        vqe = VQE(
+            estimator,
+            RealAmplitudes(),
+            partial(_mock_optimizer, inputs=inputs),
+            gradient=ParamShiftEstimatorGradient(estimator),
+        )
+        _ = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
+
+        self.assertIsNotNone(inputs["jac"])
+
+    def test_gradient_run(self):
+        """Test using the gradient to calculate the minimum."""
+        estimator = Estimator()
+        vqe = VQE(
+            estimator,
+            RealAmplitudes(),
+            GradientDescent(maxiter=200, learning_rate=0.1),
+            gradient=ParamShiftEstimatorGradient(estimator),
+        )
+        result = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
+        self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=5)
+
+    def test_with_two_qubit_reduction(self):
+        """Test the VQE using TwoQubitReduction."""
+        qubit_op = PauliSumOp.from_list(
+            [
+                ("IIII", -0.8105479805373266),
+                ("IIIZ", 0.17218393261915552),
+                ("IIZZ", -0.22575349222402472),
+                ("IZZI", 0.1721839326191556),
+                ("ZZII", -0.22575349222402466),
+                ("IIZI", 0.1209126326177663),
+                ("IZZZ", 0.16892753870087912),
+                ("IXZX", -0.045232799946057854),
+                ("ZXIX", 0.045232799946057854),
+                ("IXIX", 0.045232799946057854),
+                ("ZXZX", -0.045232799946057854),
+                ("ZZIZ", 0.16614543256382414),
+                ("IZIZ", 0.16614543256382414),
+                ("ZZZZ", 0.17464343068300453),
+                ("ZIZI", 0.1209126326177663),
+            ]
+        )
+        tapered_qubit_op = TwoQubitReduction(num_particles=2).convert(qubit_op)
+        vqe = VQE(
+            Estimator(),
+            self.ry_wavefunction,
+            SPSA(maxiter=300, last_avg=5),
+        )
+        result = vqe.compute_minimum_eigenvalue(tapered_qubit_op)
+        self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=2)
+
+    def test_callback(self):
+        """Test the callback on VQE."""
+        history = {"eval_count": [], "parameters": [], "mean": [], "metadata": []}
+        # expected_shots = options["shots"] if "shots" in options else 0
+
+        def store_intermediate_result(eval_count, parameters, mean, metadata):
+            history["eval_count"].append(eval_count)
+            history["parameters"].append(parameters)
+            history["mean"].append(mean)
+            history["metadata"].append(metadata)
+
+        optimizer = COBYLA(maxiter=3)
+        wavefunction = self.ry_wavefunction
+
+        estimator = Estimator()
+
+        vqe = VQE(
+            estimator,
+            wavefunction,
+            optimizer,
+            callback=store_intermediate_result,
+        )
+        vqe.compute_minimum_eigenvalue(operator=self.h2_op)
+
+        self.assertTrue(all(isinstance(count, int) for count in history["eval_count"]))
+        self.assertTrue(all(isinstance(mean, float) for mean in history["mean"]))
+        self.assertTrue(all(isinstance(metadata, dict) for metadata in history["metadata"]))
+        for params in history["parameters"]:
+            self.assertTrue(all(isinstance(param, float) for param in params))
+
+    def test_reuse(self):
+        """Test re-using a VQE algorithm instance."""
+        ansatz = TwoLocal(rotation_blocks=["ry", "rz"], entanglement_blocks="cz")
+        vqe = VQE(Estimator(), ansatz, SLSQP(maxiter=300))
+        with self.subTest(msg="assert VQE works once all info is available"):
+            result = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
+            self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=5)
+
+        operator = Operator(np.array([[1, 0, 0, 0], [0, -1, 0, 0], [0, 0, 2, 0], [0, 0, 0, 3]]))
+
+        with self.subTest(msg="assert vqe works on re-use."):
+            result = vqe.compute_minimum_eigenvalue(operator=operator)
+            self.assertAlmostEqual(result.eigenvalue.real, -1.0, places=5)
+
+    def test_vqe_optimizer_reuse(self):
+        """Test running same VQE twice to re-use optimizer, then switch optimizer"""
+        vqe = VQE(
+            Estimator(),
+            self.ryrz_wavefunction,
+            SLSQP(),
+        )
+
+        def run_check():
+            result = vqe.compute_minimum_eigenvalue(operator=self.h2_op)
+            self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=5)
+
+        run_check()
+
+        with self.subTest("Optimizer re-use."):
+            run_check()
+
+        with self.subTest("Optimizer replace."):
+            vqe.optimizer = L_BFGS_B()
+            run_check()
+
+    def test_batch_evaluate_with_qnspsa(self):
+        """Test batch evaluating with QNSPSA works."""
+        ansatz = TwoLocal(2, rotation_blocks=["ry", "rz"], entanglement_blocks="cz")
+
+        wrapped_sampler = Sampler()
+        inner_sampler = Sampler()
+
+        wrapped_estimator = Estimator()
+        inner_estimator = Estimator()
+
+        callcount = {"sampler": 0, "estimator": 0}
+
+        def wrapped_estimator_run(*args, **kwargs):
+            kwargs["callcount"]["estimator"] += 1
+            return inner_estimator.run(*args, **kwargs)
+
+        def wrapped_sampler_run(*args, **kwargs):
+            kwargs["callcount"]["sampler"] += 1
+            return inner_sampler.run(*args, **kwargs)
+
+        wrapped_estimator.run = partial(wrapped_estimator_run, callcount=callcount)
+        wrapped_sampler.run = partial(wrapped_sampler_run, callcount=callcount)
+
+        fidelity = ComputeUncompute(wrapped_sampler)
+
+        def fidelity_callable(left, right):
+            batchsize = np.asarray(left).shape[0]
+            job = fidelity.run(batchsize * [ansatz], batchsize * [ansatz], left, right)
+            return job.result().fidelities
+
+        qnspsa = QNSPSA(fidelity_callable, maxiter=5)
+        qnspsa.set_max_evals_grouped(100)
+
+        vqe = VQE(
+            wrapped_estimator,
+            ansatz,
+            qnspsa,
+        )
+        _ = vqe.compute_minimum_eigenvalue(Pauli("ZZ"))
+
+        # 5 (fidelity)
+        expected_sampler_runs = 5
+        # 1 calibration + 1 stddev estimation + 1 initial blocking
+        # + 5 (1 loss + 1 blocking) + 1 return loss
+        expected_estimator_runs = 1 + 1 + 1 + 5 * 2 + 1
+
+        self.assertEqual(callcount["sampler"], expected_sampler_runs)
+        self.assertEqual(callcount["estimator"], expected_estimator_runs)
+
+    def test_optimizer_scipy_callable(self):
+        """Test passing a SciPy optimizer directly as callable."""
+        vqe = VQE(
+            Estimator(),
+            self.ryrz_wavefunction,
+            partial(scipy_minimize, method="L-BFGS-B", options={"maxiter": 10}),
+        )
+        result = vqe.compute_minimum_eigenvalue(self.h2_op)
+        self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=2)
+
+    def test_optimizer_callable(self):
+        """Test passing a optimizer directly as callable."""
+        ansatz = RealAmplitudes(1, reps=1)
+        vqe = VQE(Estimator(), ansatz, _mock_optimizer)
+        result = vqe.compute_minimum_eigenvalue(SparsePauliOp("Z"))
+        self.assertTrue(np.all(result.optimal_point == np.zeros(ansatz.num_parameters)))
+
+    def test_aux_operators_list(self):
+        """Test list-based aux_operators."""
+        vqe = VQE(Estimator(), self.ry_wavefunction, SLSQP(maxiter=300))
+
+        with self.subTest("Test with an empty list."):
+            result = vqe.compute_minimum_eigenvalue(self.h2_op, aux_operators=[])
+            self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=6)
+            self.assertIsInstance(result.aux_operators_evaluated, list)
+            self.assertEqual(len(result.aux_operators_evaluated), 0)
+
+        with self.subTest("Test with two auxiliary operators."):
+            aux_op1 = PauliSumOp.from_list([("II", 2.0)])
+            aux_op2 = PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
+            aux_ops = [aux_op1, aux_op2]
+            result = vqe.compute_minimum_eigenvalue(self.h2_op, aux_operators=aux_ops)
+            self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=5)
+            self.assertEqual(len(result.aux_operators_evaluated), 2)
+            # expectation values
+            self.assertAlmostEqual(result.aux_operators_evaluated[0][0], 2.0, places=6)
+            self.assertAlmostEqual(result.aux_operators_evaluated[1][0], 0.0, places=6)
+            # metadata
+            self.assertIsInstance(result.aux_operators_evaluated[0][1], dict)
+            self.assertIsInstance(result.aux_operators_evaluated[1][1], dict)
+
+        with self.subTest("Test with additional zero operator."):
+            extra_ops = [*aux_ops, 0]
+            result = vqe.compute_minimum_eigenvalue(self.h2_op, aux_operators=extra_ops)
+            self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=5)
+            self.assertEqual(len(result.aux_operators_evaluated), 3)
+            # expectation values
+            self.assertAlmostEqual(result.aux_operators_evaluated[0][0], 2.0, places=6)
+            self.assertAlmostEqual(result.aux_operators_evaluated[1][0], 0.0, places=6)
+            self.assertAlmostEqual(result.aux_operators_evaluated[2][0], 0.0)
+            # metadata
+            self.assertIsInstance(result.aux_operators_evaluated[0][1], dict)
+            self.assertIsInstance(result.aux_operators_evaluated[1][1], dict)
+            self.assertIsInstance(result.aux_operators_evaluated[2][1], dict)
+
+    def test_aux_operators_dict(self):
+        """Test dictionary compatibility of aux_operators"""
+        vqe = VQE(Estimator(), self.ry_wavefunction, SLSQP(maxiter=300))
+
+        with self.subTest("Test with an empty dictionary."):
+            result = vqe.compute_minimum_eigenvalue(self.h2_op, aux_operators={})
+            self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=6)
+            self.assertIsInstance(result.aux_operators_evaluated, dict)
+            self.assertEqual(len(result.aux_operators_evaluated), 0)
+
+        with self.subTest("Test with two auxiliary operators."):
+            aux_op1 = PauliSumOp.from_list([("II", 2.0)])
+            aux_op2 = PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
+            aux_ops = {"aux_op1": aux_op1, "aux_op2": aux_op2}
+            result = vqe.compute_minimum_eigenvalue(self.h2_op, aux_operators=aux_ops)
+            self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=6)
+            self.assertEqual(len(result.aux_operators_evaluated), 2)
+
+            # expectation values
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op1"][0], 2.0, places=5)
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op2"][0], 0.0, places=5)
+            # metadata
+            self.assertIsInstance(result.aux_operators_evaluated["aux_op1"][1], dict)
+            self.assertIsInstance(result.aux_operators_evaluated["aux_op2"][1], dict)
+
+        with self.subTest("Test with additional zero operator."):
+            extra_ops = {**aux_ops, "zero_operator": 0}
+            result = vqe.compute_minimum_eigenvalue(self.h2_op, aux_operators=extra_ops)
+            self.assertAlmostEqual(result.eigenvalue.real, self.h2_energy, places=6)
+            self.assertEqual(len(result.aux_operators_evaluated), 3)
+            # expectation values
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op1"][0], 2.0, places=5)
+            self.assertAlmostEqual(result.aux_operators_evaluated["aux_op2"][0], 0.0, places=5)
+            self.assertAlmostEqual(result.aux_operators_evaluated["zero_operator"][0], 0.0)
+            # metadata
+            self.assertIsInstance(result.aux_operators_evaluated["aux_op1"][1], dict)
+            self.assertIsInstance(result.aux_operators_evaluated["aux_op2"][1], dict)
+            self.assertIsInstance(result.aux_operators_evaluated["zero_operator"][1], dict)
+
+
+if __name__ == "__main__":
+    unittest.main()

test/python/algorithms/test_observables_evaluator.py

@@ -9,7 +9,9 @@
 # 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 evaluator of auxiliary operators for algorithms."""
+
 from __future__ import annotations
 import unittest
 from typing import Tuple
@@ -53,7 +55,7 @@ class TestObservablesEvaluator(QiskitAlgorithmsTestCase):
         # the exact value is a list of (mean, (variance, shots)) where we expect 0 variance and
         # 0 shots
         exact = [
-            (Statevector(ansatz).expectation_value(observable), (0, 0))
+            (Statevector(ansatz).expectation_value(observable), {})
             for observable in observables_list
         ]
 
@@ -70,7 +72,7 @@ class TestObservablesEvaluator(QiskitAlgorithmsTestCase):
         observables: ListOrDict[BaseOperator | PauliSumOp],
         estimator: Estimator,
     ):
-        result = estimate_observables(estimator, quantum_state, observables, self.threshold)
+        result = estimate_observables(estimator, quantum_state, observables, None, self.threshold)
 
         if isinstance(observables, dict):
             np.testing.assert_equal(list(result.keys()), list(expected_result.keys()))
@@ -142,8 +144,8 @@ class TestObservablesEvaluator(QiskitAlgorithmsTestCase):
         state = bound_ansatz
         estimator = Estimator()
         observables = [SparsePauliOp(["XX", "YY"]), 0]
-        result = estimate_observables(estimator, state, observables, self.threshold)
-        expected_result = [(0.015607318055509564, (0, 0)), (0.0, (0, 0))]
+        result = estimate_observables(estimator, state, observables, None, self.threshold)
+        expected_result = [(0.015607318055509564, {}), (0.0, {})]
         means = [element[0] for element in result]
         expected_means = [element[0] for element in expected_result]
         np.testing.assert_array_almost_equal(means, expected_means, decimal=0.01)
@@ -152,6 +154,32 @@ class TestObservablesEvaluator(QiskitAlgorithmsTestCase):
         expected_vars_and_shots = [element[1] for element in expected_result]
         np.testing.assert_array_equal(vars_and_shots, expected_vars_and_shots)
 
+    def test_estimate_observables_shots(self):
+        """Tests that variances and shots are returned properly."""
+        ansatz = EfficientSU2(2)
+        parameters = np.array(
+            [1.2, 4.2, 1.4, 2.0, 1.2, 4.2, 1.4, 2.0, 1.2, 4.2, 1.4, 2.0, 1.2, 4.2, 1.4, 2.0],
+            dtype=float,
+        )
+
+        bound_ansatz = ansatz.bind_parameters(parameters)
+        state = bound_ansatz
+        estimator = Estimator(options={"shots": 2048})
+        observables = [PauliSumOp.from_list([("ZZ", 2.0)])]
+        result = estimate_observables(estimator, state, observables, None, self.threshold)
+        exact_result = self.get_exact_expectation(bound_ansatz, observables)
+        expected_result = [(exact_result[0][0], {"variance": 1.0898, "shots": 2048})]
+
+        means = [element[0] for element in result]
+        expected_means = [element[0] for element in expected_result]
+        np.testing.assert_array_almost_equal(means, expected_means, decimal=0.01)
+
+        vars_and_shots = [element[1] for element in result]
+        expected_vars_and_shots = [element[1] for element in expected_result]
+        for computed, expected in zip(vars_and_shots, expected_vars_and_shots):
+            self.assertAlmostEqual(computed.pop("variance"), expected.pop("variance"), 2)
+            self.assertEqual(computed.pop("shots"), expected.pop("shots"))
+
 
 if __name__ == "__main__":
     unittest.main()

test/python/algorithms/time_evolvers/test_trotter_qrte.py

@@ -88,7 +88,10 @@ class TestTrotterQRTE(QiskitAlgorithmsTestCase):
         )
 
         aux_ops_result = evolution_result.aux_ops_evaluated
-        expected_aux_ops_result = [(0.078073, (0.0, 0.0)), (0.268286, (0.0, 0.0))]
+        expected_aux_ops_result = [
+            (0.078073, {"variance": 0, "shots": 0}),
+            (0.268286, {"variance": 0, "shots": 0}),
+        ]
 
         means = [element[0] for element in aux_ops_result]
         expected_means = [element[0] for element in expected_aux_ops_result]
@@ -96,7 +99,10 @@ class TestTrotterQRTE(QiskitAlgorithmsTestCase):
 
         vars_and_shots = [element[1] for element in aux_ops_result]
         expected_vars_and_shots = [element[1] for element in expected_aux_ops_result]
-        np.testing.assert_array_equal(vars_and_shots, expected_vars_and_shots)
+
+        for computed, expected in zip(vars_and_shots, expected_vars_and_shots):
+            self.assertAlmostEqual(computed.pop("variance", 0), expected["variance"], 2)
+            self.assertEqual(computed.pop("shots", 0), expected["shots"])
 
     @data(
         (

test/python/algorithms/utils/__init__.py

@@ -0,0 +1,11 @@
+# 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.

test/python/algorithms/utils/test_validate_bounds.py

@@ -0,0 +1,53 @@
+# 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.
+
+"""Test validate bounds."""
+
+from test.python.algorithms import QiskitAlgorithmsTestCase
+
+from unittest.mock import Mock
+
+import numpy as np
+
+from qiskit.algorithms.utils import validate_bounds
+from qiskit.utils import algorithm_globals
+
+
+class TestValidateBounds(QiskitAlgorithmsTestCase):
+    """Test the ``validate_bounds`` utility function."""
+
+    def setUp(self):
+        super().setUp()
+        algorithm_globals.random_seed = 0
+        self.bounds = [(-np.pi / 2, np.pi / 2)]
+        self.ansatz = Mock()
+
+    def test_with_no_ansatz_bounds(self):
+        """Test with no ansatz bounds."""
+        self.ansatz.num_parameters = 1
+        self.ansatz.parameter_bounds = None
+        bounds = validate_bounds(self.ansatz)
+        self.assertEqual(bounds, [(None, None)])
+
+    def test_with_ansatz_bounds(self):
+        """Test with ansatz bounds."""
+        self.ansatz.num_parameters = 1
+        self.ansatz.parameter_bounds = self.bounds
+        bounds = validate_bounds(self.ansatz)
+        self.assertEqual(bounds, self.bounds)
+
+    def test_with_mismatched_num_params(self):
+        """Test with a mismatched number of parameters and bounds"""
+        self.ansatz.num_parameters = 2
+        self.ansatz.parameter_bounds = self.bounds
+        with self.assertRaises(ValueError):
+            _ = validate_bounds(self.ansatz)

test/python/algorithms/utils/test_validate_initial_point.py

@@ -0,0 +1,50 @@
+# 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.
+
+"""Test validate initial point."""
+
+from test.python.algorithms import QiskitAlgorithmsTestCase
+
+from unittest.mock import Mock
+
+import numpy as np
+
+from qiskit.algorithms.utils import validate_initial_point
+from qiskit.utils import algorithm_globals
+
+
+class TestValidateInitialPoint(QiskitAlgorithmsTestCase):
+    """Test the ``validate_initial_point`` utility function."""
+
+    def setUp(self):
+        super().setUp()
+        algorithm_globals.random_seed = 0
+        self.ansatz = Mock()
+        self.ansatz.num_parameters = 1
+
+    def test_with_no_initial_point_or_bounds(self):
+        """Test with no user-defined initial point and no ansatz bounds."""
+        self.ansatz.parameter_bounds = None
+        initial_point = validate_initial_point(None, self.ansatz)
+        np.testing.assert_array_almost_equal(initial_point, [1.721111])
+
+    def test_with_no_initial_point(self):
+        """Test with no user-defined initial point with ansatz bounds."""
+        self.ansatz.parameter_bounds = [(-np.pi / 2, np.pi / 2)]
+        initial_point = validate_initial_point(None, self.ansatz)
+        np.testing.assert_array_almost_equal(initial_point, [0.430278])
+
+    def test_with_mismatched_params(self):
+        """Test with mistmatched parameters and bounds.."""
+        self.ansatz.parameter_bounds = None
+        with self.assertRaises(ValueError):
+            _ = validate_initial_point([1.0, 2.0], self.ansatz)