diff --git a/qiskit/algorithms/__init__.py b/qiskit/algorithms/__init__.py index 322acba21dac..a57691d60811 100644 --- a/qiskit/algorithms/__init__.py +++ b/qiskit/algorithms/__init__.py @@ -92,6 +92,7 @@ NumPyEigensolver VQD + SSVQE Evolvers @@ -299,7 +300,15 @@ MaximumLikelihoodAmplitudeEstimationResult, EstimationProblem, ) -from .eigen_solvers import NumPyEigensolver, Eigensolver, EigensolverResult, VQD, VQDResult +from .eigen_solvers import ( + NumPyEigensolver, + Eigensolver, + EigensolverResult, + VQD, + VQDResult, + SSVQE, + SSVQEResult, +) from .factorizers import Shor, ShorResult from .linear_solvers import HHL, LinearSolver, NumPyLinearSolver, LinearSolverResult from .minimum_eigen_solvers import ( @@ -360,6 +369,8 @@ "EigensolverResult", "Shor", "ShorResult", + "SSVQE", + "SSVQEResult", "VQE", "VQEResult", "QAOA", diff --git a/qiskit/algorithms/eigen_solvers/__init__.py b/qiskit/algorithms/eigen_solvers/__init__.py index c567ab4f5695..f33c3cf8dfba 100644 --- a/qiskit/algorithms/eigen_solvers/__init__.py +++ b/qiskit/algorithms/eigen_solvers/__init__.py @@ -15,5 +15,14 @@ from .numpy_eigen_solver import NumPyEigensolver from .eigen_solver import Eigensolver, EigensolverResult from .vqd import VQD, VQDResult +from .ssvqe import SSVQE, SSVQEResult -__all__ = ["NumPyEigensolver", "Eigensolver", "EigensolverResult", "VQD", "VQDResult"] +__all__ = [ + "NumPyEigensolver", + "Eigensolver", + "EigensolverResult", + "VQD", + "VQDResult", + "SSVQE", + "SSVQEResult", +] diff --git a/qiskit/algorithms/eigen_solvers/ssvqe.py b/qiskit/algorithms/eigen_solvers/ssvqe.py new file mode 100644 index 000000000000..6375912550f3 --- /dev/null +++ b/qiskit/algorithms/eigen_solvers/ssvqe.py @@ -0,0 +1,953 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 2022. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""The Subspace Search Variational Quantum Eigensolver algorithm. +See https://arxiv.org/abs/1810.09434 +""" + +from typing import Optional, List, Callable, Union, Dict, Tuple +import logging +import warnings +from time import time +import numpy as np +import scipy + +from qiskit.circuit import QuantumCircuit, Parameter +from qiskit.circuit.library import RealAmplitudes +from qiskit.providers import Backend +from qiskit.opflow import ( + OperatorBase, + ExpectationBase, + ExpectationFactory, + StateFn, + CircuitStateFn, + ListOp, + CircuitSampler, + PauliSumOp, +) +from qiskit.opflow.gradients import GradientBase +from qiskit.utils.validation import validate_min +from qiskit.utils.backend_utils import is_aer_provider +from qiskit.utils import QuantumInstance, algorithm_globals +from qiskit.quantum_info import Statevector +from ..list_or_dict import ListOrDict +from ..optimizers import Optimizer, SLSQP, OptimizerResult +from ..variational_algorithm import VariationalAlgorithm, VariationalResult +from .eigen_solver import Eigensolver, EigensolverResult +from ..exceptions import AlgorithmError +from ..aux_ops_evaluator import eval_observables + +logger = logging.getLogger(__name__) + + +OBJECTIVE = Callable[[np.ndarray], float] +GRADIENT = Callable[[np.ndarray], np.ndarray] +RESULT = Union[scipy.optimize.OptimizeResult, OptimizerResult] +BOUNDS = List[Tuple[float, float]] + +MINIMIZER = Callable[ + [ + OBJECTIVE, # the objective function to minimize (the energy in the case of the VQE) + np.ndarray, # the initial point for the optimization + Optional[GRADIENT], # the gradient of the objective function + Optional[BOUNDS], # parameters bounds for the optimization + ], + RESULT, # a result object (either SciPy's or Qiskit's) +] + + +class SSVQE(VariationalAlgorithm, Eigensolver): + r"""The Subspace Search Variational Quantum Eigensolver algorithm. + `SSVQE `__ is a quantum algorithm that uses a + variational technique to find + the low-lying eigenvalues of the Hamiltonian :math:`H` of a given system. + SSVQE can be seen as a natural generalization of VQE. Instead of + minimizing the expectation value of :math:`H`, SSVQE takes a set of + mutually orthogonal input states, applies the same parameterized + ansatz to all of them, then minimizes a weighted sum + of the expectation values of :math:`H` with respect to these states. + An instance of SSVQE requires defining two algorithmic sub-components: + a trial state (a.k.a. ansatz) which is a :class:`QuantumCircuit`, and one of the classical + :mod:`~qiskit.algorithms.optimizers`. The ansatz is varied, via its set of parameters, by the + optimizer, such that it works towards a set of mutually orthogonal states, as determined by the + parameters applied to the ansatz, that will result in the minimum weighted sum of expectation values + being measured of the input operator (Hamiltonian) with respect to these states. The weights + given to this list of expectation values is given by *weight_vector*. + An optional array of parameter values, via the *initial_point*, may be provided as the + starting point for the search of the low-lying eigenvalues. This feature is particularly useful + such as when there are reasons to believe that the solution point is close to a particular + point. As an example, when building the dissociation profile of a molecule, + it is likely that using the previous computed optimal solution as the starting + initial point for the next interatomic distance is going to reduce the number of iterations + necessary for the variational algorithm to converge. It provides an + `initial point tutorial `__ detailing this use case. + The length of the *initial_point* list value must match the number of the parameters + expected by the ansatz being used. If the *initial_point* is left at the default + of ``None``, then SSVQE will look to the ansatz for a preferred value, based on its + given initial state. If the ansatz returns ``None``, + then a random point will be generated within the parameter bounds set, as per above. + If the ansatz provides ``None`` as the lower bound, then SSVQE + will default it to :math:`-2\pi`; similarly, if the ansatz returns ``None`` + as the upper bound, the default value will be :math:`2\pi`. + An optional list of initial states, via the *initial_states*, may also be provided. Choosing + these states appropriately is a critical part of the algorithm. They must be mutually orthogonal + as this is how the algorithm enforces the mutual orthogonality of the solution states. If + the *initial_states* is left as ``None``, then SSVQE will automatically generate a list of + computational basis states and use these as the initial states. For many physically-motivated + problems, it is advised to not rely on these default values as doing so can easily result in + an unphysical solution being returned. For example, if one wishes to find the low-lying + excited states of a molecular Hamiltonian, then we expect the output states to belong to + a particular particle-number subspace. If an ansatz that preserves particle number such as + :class:`UCCSD` is used, then states belonging to the incorrect particle number subspace + will be returned if the *initial_states* are not in the correct particle number subspace. + A similar statement can often be made for the spin-magnetization quantum number. + The optimizer can either be one of Qiskit's optimizers, such as + :class:`~qiskit.algorithms.optimizers.SPSA` or a callable with the following signature: + .. note:: + The callable _must_ have the argument names ``fun, x0, jac, bounds`` as indicated + in the following code block. + .. 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 to directly pass 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") + """ + + def __init__( + self, + num_states: Optional[int] = 2, + ansatz: Optional[QuantumCircuit] = None, + optimizer: Optional[Union[Optimizer, MINIMIZER]] = None, + initial_point: Optional[np.ndarray] = None, + initial_states: Optional[ + Union[List[QuantumCircuit], List[Statevector]] + ] = None, # Set of initial orthogonal states expressed as a list of QuantumCircuit objects. + weight_vector: Optional[ + Union[np.ndarray, List] + ] = None, # set of weight factors to be used in the cost function + gradient: Optional[Union[GradientBase, Callable]] = None, + expectation: Optional[ExpectationBase] = None, + include_custom: bool = False, + max_evals_grouped: int = 1, + callback: Optional[Callable[[int, np.ndarray, float, float], None]] = None, + quantum_instance: Optional[Union[QuantumInstance, Backend]] = None, + ) -> None: + """ + Args: + num_states: The number of eigenstates that the algorithm will attempt to find. + ansatz: A parameterized circuit used as Ansatz for the wave function. + optimizer: A classical optimizer. Can either be a Qiskit optimizer or a callable + that takes an array as input and returns a Qiskit or SciPy optimization result. + initial_point: An optional initial point (i.e. initial parameter values) + for the optimizer. If ``None`` then VQE will look to the ansatz for a preferred + point and if not will simply compute a random one. + initial_states: An optional list of mutually orthogonal initial states. + If ``None``, then SSVQE will set these to be a list of mutually orthogonal + computational basis states. + weight_vector: An optional list or array of real positive numbers with length + equal to the value of *num_states* to be used in the weighted energy summation + objective function. This fixes the ordering of the returned eigenstate/eigenvalue + pairs. If ``None``, then SSVQE will default to [n, n-1, ..., 1] for `num_states` = n. + gradient: An optional gradient function or operator for optimizer. + expectation: The Expectation converter for taking the average value of the + Observable over the ansatz state function. When ``None`` (the default) an + :class:`~qiskit.opflow.expectations.ExpectationFactory` is used to select + an appropriate expectation based on the operator and backend. When using Aer + qasm_simulator backend, with paulis, it is however much faster to leverage custom + Aer function for the computation but, although VQE performs much faster + with it, the outcome is ideal, with no shot noise, like using a state vector + simulator. If you are just looking for the quickest performance when choosing Aer + qasm_simulator and the lack of shot noise is not an issue then set `include_custom` + parameter here to ``True`` (defaults to ``False``). + include_custom: When `expectation` parameter here is None setting this to ``True`` will + allow the factory to include the custom Aer pauli expectation. + max_evals_grouped: Max number of evaluations performed simultaneously. Signals the + given optimizer that more than one set of parameters can be supplied so that + potentially the expectation values can be computed in parallel. Typically this is + possible when a finite difference gradient is used by the optimizer such that + multiple points to compute the gradient can be passed and if computed in parallel + improve overall execution time. Deprecated if a gradient operator or function is + given. + callback: a callback that can access the intermediate data during the optimization. + Four parameter values are passed to the callback as follows during each evaluation + by the optimizer for its current set of parameters as it works towards the minimum. + These are: the evaluation count, the optimizer parameters for the + ansatz, the evaluated mean and the evaluated standard deviation.` + quantum_instance: Quantum Instance or Backend + """ + validate_min("max_evals_grouped", max_evals_grouped, 1) + + self._num_states = num_states + + self._initial_states = None + self.initial_states = initial_states + self._weight_vector = None + self.weight_vector = weight_vector + self._initialized_ansatz_list = None + + super().__init__() + + self._max_evals_grouped = max_evals_grouped + + self._expectation = None + self.expectation = expectation + self._include_custom = include_custom + + self._ansatz = None + self.ansatz = ansatz + + self._optimizer = None + self.optimizer = optimizer + + self._initial_point = None + self.initial_point = initial_point + self._gradient = None + self.gradient = gradient + self._quantum_instance = None + + if quantum_instance is not None: + self.quantum_instance = quantum_instance + + self._eval_time = None + self._eval_count = 0 + self._callback = None + self.callback = callback + + logger.info(self.print_settings()) + + # TODO remove this once the stateful methods are deleted + self._ret = None + + @property + def weight_vector(self) -> Union[np.ndarray, List]: + """Returns the weight vector used to order the eigenvalues.""" + return self._weight_vector + + @weight_vector.setter + def weight_vector(self, vec: Union[np.ndarray, List]): + """Sets the weight vector used to order the eigenvalues. + Args: + vec: The weight vector used in the objective function. If + None is passed, then the weight vector is set to an array + of the form [num_states, num_states - 1, .... ,1]. + """ + + if vec is not None: + self._weight_vector = vec + else: + self._weight_vector = np.asarray( + [self._num_states - n for n in range(self._num_states)] + ) + + @property + def num_states(self) -> int: + """Returns the number of low-lying eigenstates which the algorithm will attempt to find.""" + return self._num_states + + @num_states.setter + def num_states(self, num: int) -> None: + """Sets the number of low-lying eigenstates which the algorithm will attempt to find. + Args: + num: The number of eigenstates to find. + """ + self._num_states = num + + @property + def ansatz(self) -> QuantumCircuit: + """Returns the ansatz.""" + return self._ansatz + + @ansatz.setter + def ansatz(self, ansatz: Optional[QuantumCircuit]) -> None: + """Sets the ansatz. + Args: + ansatz: The parameterized circuit used as an ansatz. + If None is passed, RealAmplitudes is used by default. + """ + if ansatz is None: + ansatz = RealAmplitudes(reps=6) + + self._ansatz = ansatz + + @property + def initial_states(self) -> List[Union[QuantumCircuit, Statevector]]: + """Returns the initial states.""" + + return self._initial_states + + @initial_states.setter + def initial_states(self, states: List[Union[QuantumCircuit, Statevector]]) -> None: + """Sets the initial states. + Args: + states: The initial states to be used by the algorithm.""" + + self._initial_states = states + + @property + def initialized_ansatz_list(self) -> List[QuantumCircuit]: + """Returns a list of ansatz circuits, where the nth element + in the list is a QuantumCircuit consisting of the ansatz + initialized with the nth element of the list + of initial states.""" + + return self._initialized_ansatz_list + + @initialized_ansatz_list.setter + def initialized_ansatz_list(self, initialized_states: List[QuantumCircuit]) -> None: + """Sets the list of ansatz circuits initialized in the set of initial states. + Args: initial_states: The list of orthogonal initial states.""" + + self._initialized_ansatz_list = initialized_states + + @property + def gradient(self) -> Optional[Union[GradientBase, Callable]]: + """Returns the gradient.""" + return self._gradient + + @gradient.setter + def gradient(self, gradient: Optional[Union[GradientBase, Callable]]): + """Sets the gradient.""" + self._gradient = gradient + + @property + def quantum_instance(self) -> Optional[QuantumInstance]: + """Returns quantum instance.""" + return self._quantum_instance + + @quantum_instance.setter + def quantum_instance(self, quantum_instance: Union[QuantumInstance, Backend]) -> None: + """Sets quantum_instance""" + if not isinstance(quantum_instance, QuantumInstance): + quantum_instance = QuantumInstance(quantum_instance) + + self._quantum_instance = quantum_instance + self._circuit_sampler_list = [ + CircuitSampler(quantum_instance, param_qobj=is_aer_provider(quantum_instance.backend)) + for n in range(self.num_states) + ] + + @property + def initial_point(self) -> Optional[np.ndarray]: + """Returns initial point""" + return self._initial_point + + @initial_point.setter + def initial_point(self, initial_point: np.ndarray): + """Sets initial point""" + self._initial_point = initial_point + + @property + def max_evals_grouped(self) -> int: + """Returns max_evals_grouped""" + return self._max_evals_grouped + + @max_evals_grouped.setter + def max_evals_grouped(self, max_evals_grouped: int): + """Sets max_evals_grouped""" + self._max_evals_grouped = max_evals_grouped + self.optimizer.set_max_evals_grouped(max_evals_grouped) + + @property + def include_custom(self) -> bool: + """Returns include_custom""" + return self._include_custom + + @include_custom.setter + def include_custom(self, include_custom: bool) -> None: + """Sets include_custom. If set to another value than the one that was previsously set, + the expectation attribute is reset to None. + """ + if include_custom != self._include_custom: + self._include_custom = include_custom + self.expectation = None + + @property + def callback(self) -> Optional[Callable[[int, np.ndarray, float, float], None]]: + """Returns callback""" + return self._callback + + @callback.setter + def callback(self, callback: Optional[Callable[[int, np.ndarray, float, float], None]]): + """Sets callback""" + self._callback = callback + + @property + def expectation(self) -> Optional[ExpectationBase]: + """The expectation value algorithm used to construct the expectation measurement from + the observable.""" + return self._expectation + + @expectation.setter + def expectation(self, exp: Optional[ExpectationBase]) -> None: + self._expectation = exp + + def _check_operator_ansatz(self, operator: OperatorBase): + """Check that the number of qubits of operator and ansatz match.""" + if operator is not None and self.ansatz is not None: + if operator.num_qubits != self.ansatz.num_qubits: + # try to set the number of qubits on the ansatz, if possible + try: + self.ansatz.num_qubits = operator.num_qubits + except AttributeError as ex: + 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 ex + + @property + def optimizer(self) -> Optimizer: + """Returns optimizer""" + return self._optimizer + + @optimizer.setter + def optimizer(self, optimizer: Optional[Optimizer]): + """Sets the optimizer attribute. + Args: + optimizer: The optimizer to be used. If None is passed, SLSQP is used by default. + """ + if optimizer is None: + optimizer = SLSQP() + + if isinstance(optimizer, Optimizer): + optimizer.set_max_evals_grouped(self.max_evals_grouped) + + self._optimizer = optimizer + + @property + def setting(self): + """Prepare the setting of VQE as a string.""" + ret = f"Algorithm: {self.__class__.__name__}\n" + params = "" + for key, value in self.__dict__.items(): + if key[0] == "_": + if "initial_point" in key and value is None: + params += "-- {}: {}\n".format(key[1:], "Random seed") + else: + params += f"-- {key[1:]}: {value}\n" + ret += f"{params}" + return ret + + def print_settings(self): + """ + Preparing the setting of SSVQE into a string. + Returns: + str: the formatted setting of SSVQE + """ + ret = "\n" + ret += "==================== Setting of {} ============================\n".format( + self.__class__.__name__ + ) + ret += f"{self.setting}" + ret += "===============================================================\n" + if self.ansatz is not None: + ret += "{}".format(self.ansatz.draw(output="text")) + else: + ret += "ansatz has not been set" + ret += "===============================================================\n" + if callable(self.optimizer): + ret += "Optimizer is custom callable\n" + else: + ret += f"{self._optimizer.setting}" + ret += "===============================================================\n" + return ret + + def construct_expectation_list( + self, + parameter: Union[List[float], List[Parameter], np.ndarray], + operator: OperatorBase, + return_expectation: bool = False, + ) -> List[Union[OperatorBase, Tuple[OperatorBase, ExpectationBase]]]: + r""" + Generate a list of the initialized ansatz circuits and + expectation value measurements, and return their + runnable compositions. + Args: + parameter: Parameters for the ansatz circuit. + operator: Qubit operator of the Observable + return_expectation: If True, return the ``ExpectationBase`` expectation converter used + in the construction of the expectation value. Useful e.g. to compute the standard + deviation of the expectation value. + Returns: + A list of the Operators equalling the measurement of the initialized ansatzes + :class:`StateFn` (where the ansatzes are initialized in the initial orthogonal states) + by the Observable's expectation :class:`StateFn`, + and, optionally, the expectation converters. + Raises: + AlgorithmError: If no operator has been provided. + AlgorithmError: If no expectation is passed and None could be inferred via the + ExpectationFactory. + """ + if operator is None: + raise AlgorithmError("The operator was never provided.") + + self._check_operator_ansatz(operator) + + # if expectation was never created, try to create one + if self.expectation is None: + expectation = ExpectationFactory.build( + operator=operator, + backend=self.quantum_instance, + include_custom=self._include_custom, + ) + else: + expectation = self.expectation + + wave_function_list = [ + self.initialized_ansatz_list[n].assign_parameters(parameter) + for n in range(self.num_states) + ] + + observable_meas = expectation.convert(StateFn(operator, is_measurement=True)) + ansatz_circuit_op_list = [ + CircuitStateFn(wave_function_list[n]) for n in range(self.num_states) + ] + expect_op_list = [ + observable_meas.compose(ansatz_circuit_op_list[n]).reduce() + for n in range(self.num_states) + ] + + if return_expectation: + return expect_op_list, expectation + + return expect_op_list + + def construct_circuits( + self, + parameter: Union[List[float], List[Parameter], np.ndarray], + operator: OperatorBase, + ) -> List[QuantumCircuit]: + """Return the circuits used to compute the expectation values with respect to + each state involved in the weighted summation of expectation values. + Args: + parameter: Parameters for the ansatz circuit. + operator: Qubit operator of the Observable + Returns: + A list of lists of the circuits used to compute the expectation value, + where the nth list is a list of the circuits used to compute the expectation + value of the nth state. + """ + expect_op_list = self.construct_expectation_list(parameter, operator) + + for n in range(self.num_states): + expect_op_list[n] = expect_op_list[n].to_circuit_op() + + circuits = [] + list_of_circuits = [] + + def extract_circuits(op): + if isinstance(op, CircuitStateFn): + circuits.append(op.primitive) + elif isinstance(op, ListOp): + for op_i in op.oplist: + extract_circuits(op_i) + + for n in range(self.num_states): + extract_circuits(expect_op_list[n]) + list_of_circuits.append(circuits) + circuits = [] + + return list_of_circuits + + @classmethod + def supports_aux_operators(cls) -> bool: + return True + + def compute_eigenvalues( + self, operator: OperatorBase, aux_operators: Optional[ListOrDict[OperatorBase]] = None + ) -> EigensolverResult: + super().compute_eigenvalues(operator, aux_operators) + + if self.quantum_instance is None: + raise AlgorithmError( + "A QuantumInstance or Backend must be supplied to run the quantum algorithm." + ) + self.quantum_instance.circuit_summary = True + + # this sets the size of the ansatz, so it must be called before the initial point + # validation + self._check_operator_ansatz(operator) + + # set an expectation for this algorithm run (will be reset to None at the end) + initial_point = _validate_initial_point(self.initial_point, self.ansatz) + + bounds = _validate_bounds(self.ansatz) + # We need to handle the array entries being zero or Optional i.e. having value None + + if self.initial_states is None: + self.initial_states = [ + QuantumCircuit(self.ansatz.num_qubits) for n in range(self.num_states) + ] + for n in range(self.num_states): + # if no initial states are provided, set them to be the first + # n computational basis states. + print(self.ansatz.num_qubits) + self.initial_states[n].initialize( + Statevector.from_int(n, 2**self.ansatz.num_qubits) + ) + + for n in range(self.num_states): + # if any element in the list of initial states is :class:`Statevector`, + # convert it to :class:`QuantumCircuit` + if isinstance(self.initial_states[n], Statevector): + temporary_circ = QuantumCircuit(self.ansatz.num_qubits) + temporary_circ.initialize(self.initial_states[n]) + self.initial_states[n] = temporary_circ + + if self.initialized_ansatz_list is None: + self.initialized_ansatz_list = [ + self.initial_states[n].compose(self.ansatz) for n in range(self.num_states) + ] + + if aux_operators: + zero_op = PauliSumOp.from_list([("I" * self.ansatz.num_qubits, 0)]) + + # Convert the None and zero values when aux_operators is a list. + # Drop None and convert zero values when aux_operators is a dict. + if isinstance(aux_operators, list): + key_op_iterator = enumerate(aux_operators) + converted = [zero_op] * len(aux_operators) + else: + key_op_iterator = aux_operators.items() + converted = {} + for key, op in key_op_iterator: + if op is not None: + converted[key] = zero_op if op == 0 else op + + aux_operators = converted + + else: + aux_operators = None + + if isinstance(self._gradient, GradientBase): + + gradient_op = 0 + for n in range(self.num_states): + gradient_op += float(self.weight_vector[n]) * ( + ~StateFn(operator) @ StateFn(self.initialized_ansatz_list[n]) + ) + + gradient = self._gradient.gradient_wrapper( + gradient_op, + bind_params=list(self.ansatz.parameters), + backend=self._quantum_instance, + ) + else: + gradient = self._gradient + + self._eval_count = 0 + weighted_energy_sum_evaluation, expectation = self.get_weighted_energy_sum_evaluation( + operator, return_expectation=True + ) + + start_time = time() + + if callable(self.optimizer): + opt_result = self.optimizer( # pylint: disable=not-callable + fun=weighted_energy_sum_evaluation, x0=initial_point, jac=gradient, bounds=bounds + ) + else: + # keep this until Optimizer.optimize is removed + try: + opt_result = self.optimizer.minimize( + fun=weighted_energy_sum_evaluation, + x0=initial_point, + jac=gradient, + bounds=bounds, + ) + except AttributeError: + # self.optimizer is an optimizer with the deprecated interface that uses + # ``optimize`` instead of ``minimize``` + warnings.warn( + "Using an optimizer that is run with the ``optimize`` method is " + "deprecated as of Qiskit Terra 0.19.0 and will be unsupported no " + "sooner than 3 months after the release date. Instead use an optimizer " + "providing ``minimize`` (see qiskit.algorithms.optimizers.Optimizer).", + DeprecationWarning, + stacklevel=2, + ) + + opt_result = self.optimizer.optimize( + len(initial_point), + weighted_energy_sum_evaluation, + gradient, + bounds, + initial_point, + ) + + eval_time = time() - start_time + + result = SSVQEResult() + result.optimal_point = opt_result.x + result.optimal_parameters = dict(zip(self.ansatz.parameters, opt_result.x)) + result.optimal_value = opt_result.fun + result.cost_function_evals = opt_result.nfev + result.optimizer_time = eval_time + result.eigenvalues = self._get_eigenvalues(result.optimal_parameters, operator) + result.eigenstates = self._get_eigenstates(result.optimal_parameters) + + logger.info( + "Optimization complete in %s seconds.\nFound opt_params %s in %s evals", + eval_time, + result.optimal_point, + self._eval_count, + ) + + # TODO delete as soon as get_optimal_vector etc are removed + self._ret = result + + if aux_operators is not None: + bound_ansatz_list = [ + self.initialized_ansatz_list[n].bind_parameters(result.optimal_point) + for n in range(self.num_states) + ] + + aux_values_list = [ + eval_observables( + self.quantum_instance, + bound_ansatz_list[n], + aux_operators, + expectation=expectation, + ) + for n in range(self.num_states) + ] + result.aux_operator_eigenvalues = aux_values_list + + return result + + def get_weighted_energy_sum_evaluation( + self, + operator: OperatorBase, + return_expectation: bool = False, + ) -> Callable[[np.ndarray], Union[float, List[float]]]: + """Returns a function handle to evaluates the weighted energy sum at given parameters + for the ansatz. This is the objective function to be passed to the optimizer + that is used for evaluation. + Args: + operator: The operator whose energy levels to evaluate. + return_expectation: If True, return the ``ExpectationBase`` expectation converter used + in the construction of the expectation value. Useful e.g. to evaluate other + operators with the same expectation value converter. + Returns: + Weighted energy sum of the hamiltonian of each parameter, and, optionally, the expectation + converter. + Raises: + RuntimeError: If the circuit is not parameterized (i.e. has 0 free parameters). + """ + num_parameters = self.ansatz.num_parameters + if num_parameters == 0: + raise RuntimeError("The ansatz must be parameterized, but has 0 free parameters.") + + ansatz_params = self.ansatz.parameters + expect_op_list, expectation = self.construct_expectation_list( + ansatz_params, operator, return_expectation=True + ) + + def weighted_energy_sum_evaluation(parameters): + parameter_sets = np.reshape(parameters, (-1, num_parameters)) + # Create dict associating each parameter with the lists of parameterization values for it + param_bindings = dict(zip(ansatz_params, parameter_sets.transpose().tolist())) + + start_time = time() + sampled_expect_op_list = [ + self._circuit_sampler_list[n].convert(expect_op_list[n], params=param_bindings) + for n in range(self.num_states) + ] + list_of_means = np.asarray( + [np.real(sampled_expect_op_list[n].eval()) for n in range(self.num_states)] + ) + + energy_values = np.copy(list_of_means) + + for n in range(self.num_states): + list_of_means[n, :] *= self.weight_vector[n] + + weighted_sum_means = np.sum(list_of_means, axis=0) + + if self._callback is not None: + variance_list = [ + np.real(expectation.compute_variance(sampled_expect_op_list[n])) + for n in range(self.num_states) + ] + estimator_error_list = [ + np.sqrt(variance_list[n] / self.quantum_instance.run_config.shots) + for n in range(self.num_states) + ] + for i, param_set in enumerate(parameter_sets): + self._eval_count += 1 + self._callback( + self._eval_count, + param_set, + [energy_values[n][i] for n in range(self.num_states)], + weighted_sum_means[i], + [estimator_error_list[n][i] for n in range(self.num_states)], + ) + else: + self._eval_count += len(weighted_sum_means) + + end_time = time() + logger.info( + "Energies evaluation returned %s - %.5f (ms), eval count: %s", + list_of_means, + (end_time - start_time) * 1000, + self._eval_count, + ) + + return weighted_sum_means if len(weighted_sum_means) > 1 else weighted_sum_means[0] + + if return_expectation: + return weighted_energy_sum_evaluation, expectation + + return weighted_energy_sum_evaluation + + def _get_eigenstates(self, optimal_parameters) -> Union[List[float], Dict[str, int]]: + """Get the simulation outcome of the ansatz, provided with parameters.""" + optimal_circuits_list = [ + self.initialized_ansatz_list[n].bind_parameters(optimal_parameters) + for n in range(self.num_states) + ] + state_fns_list = [ + self._circuit_sampler_list[n].convert(StateFn(optimal_circuits_list[n])).eval() + for n in range(self.num_states) + ] + if self.quantum_instance.is_statevector: + list_of_states = [ + state_fns_list[n].primitive.data for n in range(self.num_states) + ] # VectorStateFn -> Statevector -> np.array + else: + list_of_states = [ + state_fns_list[n].to_dict_fn().primitive for n in range(self.num_states) + ] # SparseVectorStateFn -> DictStateFn -> dict + + return list_of_states + + def _get_eigenvalues(self, optimal_parameters, operator) -> List[float]: + """Get the eigenvalue results at the optimal parameters.""" + + if self.expectation is None: + expectation = ExpectationFactory.build( + operator=operator, + backend=self.quantum_instance, + include_custom=self._include_custom, + ) + else: + expectation = self.expectation + + optimal_circuits_list = [ + self.initialized_ansatz_list[n].assign_parameters(optimal_parameters) + for n in range(self.num_states) + ] + + observable_meas = expectation.convert(StateFn(operator, is_measurement=True)) + ansatz_circuit_op_list = [ + CircuitStateFn(optimal_circuits_list[n]) for n in range(self.num_states) + ] + expect_op_list = [ + observable_meas.compose(ansatz_circuit_op_list[n]).reduce() + for n in range(self.num_states) + ] + sampled_expect_op_list = [ + self._circuit_sampler_list[n].convert(expect_op_list[n]) for n in range(self.num_states) + ] + list_of_means = np.asarray( + [sampled_expect_op_list[n].eval() for n in range(self.num_states)] + ) + return list_of_means + + +class SSVQEResult(VariationalResult, EigensolverResult): + """SSVQE Result.""" + + def __init__(self) -> None: + super().__init__() + self._cost_function_evals = None + + @property + def cost_function_evals(self) -> Optional[int]: + """Returns number of cost optimizer evaluations""" + return self._cost_function_evals + + @cost_function_evals.setter + def cost_function_evals(self, value: int) -> None: + """Sets number of cost function evaluations""" + self._cost_function_evals = value + + @property + def eigenstates(self) -> Optional[np.ndarray]: + """return eigen states""" + return self._eigenstates + + @eigenstates.setter + def eigenstates(self, value: np.ndarray) -> None: + """set eigen state""" + self._eigenstates = value + + +def _validate_initial_point(point, ansatz): + expected_size = ansatz.num_parameters + + # try getting the initial point from the ansatz + if point is None and hasattr(ansatz, "preferred_init_points"): + point = ansatz.preferred_init_points + # if the point is None choose a random initial point + + if point is None: + # get bounds if ansatz has them set, otherwise use [-2pi, 2pi] for each parameter + bounds = getattr(ansatz, "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 + + +def _validate_bounds(ansatz): + if hasattr(ansatz, "parameter_bounds") and ansatz.parameter_bounds is not None: + bounds = ansatz.parameter_bounds + if len(bounds) != ansatz.num_parameters: + raise ValueError( + f"The number of bounds ({len(bounds)}) does not match the number of " + f"parameters in the circuit ({ansatz.num_parameters})." + ) + else: + bounds = [(None, None)] * ansatz.num_parameters + + return bounds diff --git a/releasenotes/notes/adds-SSVQE-algorithm-c0f339ea56b2e04b.yaml b/releasenotes/notes/adds-SSVQE-algorithm-c0f339ea56b2e04b.yaml new file mode 100644 index 000000000000..6fd00665117a --- /dev/null +++ b/releasenotes/notes/adds-SSVQE-algorithm-c0f339ea56b2e04b.yaml @@ -0,0 +1,40 @@ +--- +features: + - | + Adds the Subspace Search Variational Quantum Eigensolver (:class:`SSVQE`) + algorithm. SSVQE is a generalization of the Variational Quantum Eigensolver + (VQE) which aims to find the low-lying eigenvalue/eigenvector pairs of a + Hermitian operator :math:`H`. It takes a set of mutually orthogonal input + states, applies a parameterized ansatz (:class:`QuantumCircuit`) to all of + them, then minimizes a weighted sum of the expectation values of :math:`H` + with respect to these parameterized states using a classical optimizer + (:class:`qiskit.algorithms.optimizers.Optimizer`). + For example: + + .. code-block:: python + + from qiskit import QuantumCircuit + from qiskit.opflow import Z + from qiskit.circuit.library import RealAmplitudes + from qiskit.utils import QuantumInstance + from qiskit.providers.aer import AerSimulator + from qiskit.algorithms.optimizers import SPSA + from ssvqe import SSVQE + + operator = Z^Z + input_states = [QuantumCircuit(2), QuantumCircuit(2)] + input_states[0].x(0) + input_states[1].x(1) + quantum_instance = QuantumInstance(backend=AerSimulator()) + + ssvqe_instance = SSVQE(num_states=2, + optimizer=SPSA(), + ansatz=RealAmplitudes(2), + initial_states=input_states, + quantum_instance=quantum_instance) + + result = ssvqe_instance.compute_eigenvalues(operator) + + + + diff --git a/test/python/algorithms/test_ssvqe.py b/test/python/algorithms/test_ssvqe.py new file mode 100644 index 000000000000..d8a307776f24 --- /dev/null +++ b/test/python/algorithms/test_ssvqe.py @@ -0,0 +1,636 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2018, 2021. +# +# 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 SSVQE """ + +import unittest +from test.python.algorithms import QiskitAlgorithmsTestCase + +import numpy as np +from ddt import data, ddt, unpack + +from qiskit import BasicAer, QuantumCircuit +from qiskit.algorithms import AlgorithmError +from qiskit.algorithms.eigen_solvers import SSVQE, SSVQEResult +from qiskit.algorithms.optimizers import ( + COBYLA, + L_BFGS_B, + SLSQP, +) +from qiskit.circuit.library import EfficientSU2, RealAmplitudes, TwoLocal +from qiskit.exceptions import MissingOptionalLibraryError +from qiskit.opflow import ( + AerPauliExpectation, + I, + MatrixExpectation, + MatrixOp, + PauliExpectation, + PauliSumOp, + PrimitiveOp, + X, + Z, +) + +from qiskit.utils import QuantumInstance, algorithm_globals, has_aer + + +if has_aer(): + from qiskit import Aer + + +@ddt +class TestSSVQE(QiskitAlgorithmsTestCase): + """Test SSVQE""" + + def setUp(self): + super().setUp() + self.seed = 50 + algorithm_globals.random_seed = self.seed + self.h2_op = ( + -1.052373245772859 * (I ^ I) + + 0.39793742484318045 * (I ^ Z) + - 0.39793742484318045 * (Z ^ I) + - 0.01128010425623538 * (Z ^ Z) + + 0.18093119978423156 * (X ^ X) + ) + self.h2_energy = -1.85727503 + self.h2_energy_excited = [-1.85727503, -1.24458455] + + self.test_op = MatrixOp(np.diagflat([3, 5, -1, 0.8, 0.2, 2, 1, -3])).to_pauli_op() + self.test_results = [-3, -1] + + self.ryrz_wavefunction = TwoLocal( + rotation_blocks=["ry", "rz"], entanglement_blocks="cz", reps=1 + ) + self.ry_wavefunction = TwoLocal(rotation_blocks="ry", entanglement_blocks="cz") + + self.qasm_simulator = QuantumInstance( + BasicAer.get_backend("qasm_simulator"), + shots=2048, + seed_simulator=self.seed, + seed_transpiler=self.seed, + ) + self.statevector_simulator = QuantumInstance( + BasicAer.get_backend("statevector_simulator"), + shots=1, + seed_simulator=self.seed, + seed_transpiler=self.seed, + ) + + def test_basic_aer_statevector(self): + """Test SSVQE on BasicAer's statevector simulator.""" + wavefunction = self.ryrz_wavefunction + ssvqe = SSVQE( + num_states=2, + ansatz=wavefunction, + weight_vector=[2, 1], + optimizer=COBYLA(), + quantum_instance=QuantumInstance( + BasicAer.get_backend("statevector_simulator"), + basis_gates=["u1", "u2", "u3", "cx", "id"], + coupling_map=[[0, 1]], + seed_simulator=algorithm_globals.random_seed, + seed_transpiler=algorithm_globals.random_seed, + ), + ) + + result = ssvqe.compute_eigenvalues(operator=self.h2_op) + + with self.subTest(msg="test eigenvalue"): + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=1 + ) + + with self.subTest(msg="test dimension of optimal point"): + self.assertEqual(len(result.optimal_point), 8) + + 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) + + def test_mismatching_num_qubits(self): + """Ensuring circuit and operator mismatch is caught""" + wavefunction = QuantumCircuit(1) + optimizer = SLSQP(maxiter=50) + ssvqe = SSVQE( + num_states=1, + ansatz=wavefunction, + weight_vector=[1], + optimizer=optimizer, + quantum_instance=self.statevector_simulator, + ) + with self.assertRaises(AlgorithmError): + _ = ssvqe.compute_eigenvalues(operator=self.h2_op) + + @data( + (MatrixExpectation(), 1), + (AerPauliExpectation(), 1), + (PauliExpectation(), 2), + ) + @unpack + def test_construct_circuits(self, expectation, num_circuits): + """Test construct circuits returns List[QuantumCircuits] and the right number of them.""" + try: + wavefunction = EfficientSU2(2, reps=1) + ssvqe = SSVQE( + num_states=2, + optimizer=COBYLA(maxiter=1), + ansatz=wavefunction, + weight_vector=[2, 1], + expectation=expectation, + quantum_instance=self.statevector_simulator, + ) + ssvqe.compute_eigenvalues(operator=self.h2_op) + params = [0] * wavefunction.num_parameters + circuits = ssvqe.construct_circuits(parameter=params, operator=self.h2_op) + + for circuits_list in circuits: + self.assertEqual(len(circuits_list), num_circuits) + for circ_list in circuits: + for circ in circ_list: + self.assertIsInstance(circ, QuantumCircuit) + except MissingOptionalLibraryError as ex: + self.skipTest(str(ex)) + return + + def test_missing_varform_params(self): + """Test specifying a variational form with no parameters raises an error.""" + circuit = QuantumCircuit(self.h2_op.num_qubits) + ssvqe = SSVQE( + num_states=1, + ansatz=circuit, + quantum_instance=BasicAer.get_backend("statevector_simulator"), + ) + with self.assertRaises(RuntimeError): + ssvqe.compute_eigenvalues(operator=self.h2_op) + + def test_basic_aer_qasm(self): + """Test SSVQE on BasicAer's QASM simulator.""" + optimizer = COBYLA(maxiter=1000) + wavefunction = self.ry_wavefunction + + ssvqe = SSVQE( + num_states=2, + weight_vector=[2, 1], + ansatz=wavefunction, + optimizer=optimizer, + max_evals_grouped=1, + quantum_instance=self.qasm_simulator, + ) + + result = ssvqe.compute_eigenvalues(operator=self.h2_op) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=1 + ) + + @unittest.skipUnless(has_aer(), "qiskit-aer doesn't appear to be installed.") + def test_with_aer_statevector(self): + """Test SSVQE with Aer's statevector_simulator.""" + backend = Aer.get_backend("aer_simulator_statevector") + wavefunction = self.ry_wavefunction + optimizer = L_BFGS_B() + + quantum_instance = QuantumInstance( + backend, + seed_simulator=algorithm_globals.random_seed, + seed_transpiler=algorithm_globals.random_seed, + ) + ssvqe = SSVQE( + num_states=2, + weight_vector=[2, 1], + ansatz=wavefunction, + optimizer=optimizer, + max_evals_grouped=1, + quantum_instance=quantum_instance, + ) + + result = ssvqe.compute_eigenvalues(operator=self.h2_op) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + + @unittest.skipUnless(has_aer(), "qiskit-aer doesn't appear to be installed.") + def test_with_aer_qasm(self): + """Test SSVQE with Aer's qasm_simulator.""" + backend = Aer.get_backend("aer_simulator") + optimizer = COBYLA(maxiter=1000) + wavefunction = self.ry_wavefunction + + quantum_instance = QuantumInstance( + backend, + seed_simulator=algorithm_globals.random_seed, + seed_transpiler=algorithm_globals.random_seed, + ) + + ssvqe = SSVQE( + num_states=2, + weight_vector=[2, 1], + ansatz=wavefunction, + optimizer=optimizer, + expectation=PauliExpectation(), + quantum_instance=quantum_instance, + ) + + result = ssvqe.compute_eigenvalues(operator=self.h2_op) + + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=1 + ) + + @unittest.skipUnless(has_aer(), "qiskit-aer doesn't appear to be installed.") + def test_with_aer_qasm_snapshot_mode(self): + """Test SSVQE using Aer's qasm_simulator snapshot mode.""" + + backend = Aer.get_backend("aer_simulator") + optimizer = COBYLA(maxiter=400) + wavefunction = self.ryrz_wavefunction + + quantum_instance = QuantumInstance( + backend, + shots=100, + seed_simulator=algorithm_globals.random_seed, + seed_transpiler=algorithm_globals.random_seed, + ) + ssvqe = SSVQE( + num_states=2, + weight_vector=[2, 1], + ansatz=wavefunction, + optimizer=optimizer, + expectation=AerPauliExpectation(), + quantum_instance=quantum_instance, + ) + + result = ssvqe.compute_eigenvalues(operator=self.test_op) + np.testing.assert_array_almost_equal(result.eigenvalues.real, self.test_results, decimal=1) + + def test_callback(self): + """Test the callback on SSVQE.""" + history = { + "eval_count": [], + "parameters": [], + "mean_list": [], + "weighted_sum_mean_list": [], + "std_list": [], + } + + def store_intermediate_result(eval_count, parameters, mean, weighted_sum_mean, std): + history["eval_count"].append(eval_count) + history["parameters"].append(parameters) + history["mean_list"].append(mean) + history["weighted_sum_mean_list"].append(weighted_sum_mean) + history["std_list"].append(std) + + optimizer = COBYLA(maxiter=3) + wavefunction = self.ry_wavefunction + + ssvqe = SSVQE( + num_states=2, + weight_vector=[2, 1], + ansatz=wavefunction, + optimizer=optimizer, + callback=store_intermediate_result, + quantum_instance=self.qasm_simulator, + ) + ssvqe.compute_eigenvalues(operator=self.h2_op) + + self.assertTrue(all(isinstance(count, int) for count in history["eval_count"])) + self.assertTrue( + all(isinstance(mean, float) for mean_list in history["mean_list"] for mean in mean_list) + ) + self.assertTrue(all(isinstance(mean, float) for mean in history["weighted_sum_mean_list"])) + self.assertTrue( + all(isinstance(std, float) for std_list in history["std_list"] for std in std_list) + ) + for params in history["parameters"]: + self.assertTrue(all(isinstance(param, float) for param in params)) + + ref_eval_count = [1, 2, 3] + ref_mean = [[-1.06, -1.43], [-1.46, -1.06], [-1.36, -0.94]] + ref_weighted_sum_mean = [-3.55, -3.98, -3.66] + ref_std = [[0.01, 0.01], [0.01, 0.01], [0.01, 0.01]] + + np.testing.assert_array_almost_equal(history["eval_count"], ref_eval_count, decimal=0) + np.testing.assert_array_almost_equal(history["mean_list"], ref_mean, decimal=2) + np.testing.assert_array_almost_equal(history["std_list"], ref_std, decimal=2) + np.testing.assert_array_almost_equal( + history["weighted_sum_mean_list"], ref_weighted_sum_mean, decimal=2 + ) + + def test_reuse(self): + """Test re-using an SSVQE algorithm instance.""" + ssvqe = SSVQE(num_states=1, weight_vector=[1]) + with self.subTest(msg="assert running empty raises AlgorithmError"): + with self.assertRaises(AlgorithmError): + _ = ssvqe.compute_eigenvalues(operator=self.h2_op) + + ansatz = TwoLocal(rotation_blocks=["ry", "rz"], entanglement_blocks="cz") + ssvqe.ansatz = ansatz + with self.subTest(msg="assert missing operator raises AlgorithmError"): + with self.assertRaises(AlgorithmError): + _ = ssvqe.compute_eigenvalues(operator=self.h2_op) + + ssvqe.expectation = MatrixExpectation() + ssvqe.quantum_instance = self.statevector_simulator + with self.subTest(msg="assert SSVQE works once all info is available"): + result = ssvqe.compute_eigenvalues(operator=self.h2_op) + np.testing.assert_array_almost_equal(result.eigenvalues.real, self.h2_energy, decimal=2) + + operator = PrimitiveOp(np.array([[1, 0, 0, 0], [0, -1, 0, 0], [0, 0, 2, 0], [0, 0, 0, 3]])) + + with self.subTest(msg="assert minimum eigensolver interface works"): + result = ssvqe.compute_eigenvalues(operator=operator) + self.assertAlmostEqual(result.eigenvalues.real[0], -1.0, places=5) + + def test_ssvqe_optimizer(self): + """Test running same SSVQE twice to re-use optimizer, then switch optimizer""" + ssvqe = SSVQE( + num_states=2, + weight_vector=[2, 1], + optimizer=SLSQP(), + quantum_instance=QuantumInstance(BasicAer.get_backend("statevector_simulator")), + ) + + def run_check(): + result = ssvqe.compute_eigenvalues(operator=self.h2_op) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=3 + ) + + run_check() + + with self.subTest("Optimizer re-use"): + run_check() + + with self.subTest("Optimizer replace"): + ssvqe.optimizer = L_BFGS_B() + run_check() + + @data(MatrixExpectation(), None) + def test_backend_change(self, user_expectation): + """Test that SSVQE works when backend changes.""" + ssvqe = SSVQE( + num_states=1, + weight_vector=[1], + ansatz=TwoLocal(rotation_blocks=["ry", "rz"], entanglement_blocks="cz"), + optimizer=SLSQP(maxiter=2), + expectation=user_expectation, + quantum_instance=BasicAer.get_backend("statevector_simulator"), + ) + result0 = ssvqe.compute_eigenvalues(operator=self.h2_op) + if user_expectation is not None: + with self.subTest("User expectation kept."): + self.assertEqual(ssvqe.expectation, user_expectation) + + ssvqe.quantum_instance = BasicAer.get_backend("qasm_simulator") + + # works also if no expectation is set, since it will be determined automatically + result1 = ssvqe.compute_eigenvalues(operator=self.h2_op) + + if user_expectation is not None: + with self.subTest("Change backend with user expectation, it is kept."): + self.assertEqual(ssvqe.expectation, user_expectation) + + with self.subTest("Check results."): + self.assertEqual(len(result0.optimal_point), len(result1.optimal_point)) + + def test_set_ansatz_to_none(self): + """Tests that setting the ansatz to None results in the default behavior""" + ssvqe = SSVQE( + num_states=1, + weight_vector=[1], + ansatz=self.ryrz_wavefunction, + optimizer=L_BFGS_B(), + quantum_instance=self.statevector_simulator, + ) + ssvqe.ansatz = None + self.assertIsInstance(ssvqe.ansatz, RealAmplitudes) + + def test_set_optimizer_to_none(self): + """Tests that setting the optimizer to None results in the default behavior""" + ssvqe = SSVQE( + num_states=1, + weight_vector=[1], + ansatz=self.ryrz_wavefunction, + optimizer=L_BFGS_B(), + quantum_instance=self.statevector_simulator, + ) + ssvqe.optimizer = None + self.assertIsInstance(ssvqe.optimizer, SLSQP) + + def test_aux_operators_list(self): + """Test list-based aux_operators.""" + wavefunction = self.ry_wavefunction + ssvqe = SSVQE( + num_states=2, + weight_vector=[2, 1], + ansatz=wavefunction, + quantum_instance=self.statevector_simulator, + ) + + # Start with an empty list + result = ssvqe.compute_eigenvalues(self.h2_op, aux_operators=[]) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + self.assertIsNone(result.aux_operator_eigenvalues) + + # Go again 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 = ssvqe.compute_eigenvalues(self.h2_op, aux_operators=aux_ops) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + self.assertEqual(len(result.aux_operator_eigenvalues), 2) + # expectation values + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0][0], 2, places=2) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][1][0], 0, places=2) + # standard deviations + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][1][1], 0.0) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][1][1], 0.0) + + # Go again with additional None and zero operators + extra_ops = [*aux_ops, None, 0] + result = ssvqe.compute_eigenvalues(self.h2_op, aux_operators=extra_ops) + np.testing.assert_array_almost_equal( + np.sort(result.eigenvalues.real), self.h2_energy_excited, decimal=2 + ) + self.assertEqual(len(result.aux_operator_eigenvalues), 2) + # expectation values + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0][0], 2, places=2) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][1][0], 0, places=2) + self.assertEqual(result.aux_operator_eigenvalues[0][2][0], 0.0) + self.assertEqual(result.aux_operator_eigenvalues[0][3][0], 0.0) + # standard deviations + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0][1], 0.0) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][1][1], 0.0) + self.assertEqual(result.aux_operator_eigenvalues[0][3][1], 0.0) + + def test_aux_operators_dict(self): + """Test dictionary compatibility of aux_operators""" + wavefunction = self.ry_wavefunction + ssvqe = SSVQE(ansatz=wavefunction, quantum_instance=self.statevector_simulator) + + # Start with an empty dictionary + result = ssvqe.compute_eigenvalues(self.h2_op, aux_operators={}) + np.testing.assert_array_almost_equal( + result.eigenvalues.real, self.h2_energy_excited, decimal=2 + ) + self.assertIsNone(result.aux_operator_eigenvalues) + + # Go again 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 = ssvqe.compute_eigenvalues(self.h2_op, aux_operators=aux_ops) + self.assertEqual(len(result.eigenvalues), 2) + self.assertEqual(len(result.eigenstates), 2) + self.assertEqual(result.eigenvalues.dtype, np.complex128) + self.assertAlmostEqual(result.eigenvalues[0], -1.85727503) + self.assertEqual(len(result.aux_operator_eigenvalues), 2) + self.assertEqual(len(result.aux_operator_eigenvalues[0]), 2) + # expectation values + self.assertAlmostEqual(result.aux_operator_eigenvalues[0]["aux_op1"][0], 2, places=6) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0]["aux_op2"][0], 0, places=1) + # standard deviations + self.assertAlmostEqual(result.aux_operator_eigenvalues[0]["aux_op1"][1], 0.0) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0]["aux_op2"][1], 0.0) + + # Go again with additional None and zero operators + extra_ops = {**aux_ops, "None_operator": None, "zero_operator": 0} + result = ssvqe.compute_eigenvalues(self.h2_op, aux_operators=extra_ops) + self.assertEqual(len(result.eigenvalues), 2) + self.assertEqual(len(result.eigenstates), 2) + self.assertEqual(result.eigenvalues.dtype, np.complex128) + self.assertAlmostEqual(np.sort(result.eigenvalues)[0], -1.85727503) + self.assertEqual(len(result.aux_operator_eigenvalues), 2) + self.assertEqual(len(result.aux_operator_eigenvalues[0]), 3) + # expectation values + self.assertAlmostEqual(result.aux_operator_eigenvalues[0]["aux_op1"][0], 2, places=6) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0]["aux_op2"][0], 0, places=6) + self.assertEqual(result.aux_operator_eigenvalues[0]["zero_operator"][0], 0.0) + self.assertTrue("None_operator" not in result.aux_operator_eigenvalues[0].keys()) + # standard deviations + self.assertAlmostEqual(result.aux_operator_eigenvalues[0]["aux_op1"][1], 0.0) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0]["aux_op2"][1], 0.0) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0]["zero_operator"][1], 0.0) + + def test_aux_operator_std_dev_pauli(self): + """Test non-zero standard deviations of aux operators with PauliExpectation.""" + wavefunction = self.ry_wavefunction + ssvqe = SSVQE( + ansatz=wavefunction, + expectation=PauliExpectation(), + initial_point=[ + 1.70256666, + -5.34843975, + -0.39542903, + 5.99477786, + -2.74374986, + -4.85284669, + 0.2442925, + -1.51638917, + ], + optimizer=COBYLA(maxiter=0), + quantum_instance=self.qasm_simulator, + ) + + # Go again 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 = ssvqe.compute_eigenvalues(self.h2_op, aux_operators=aux_ops) + self.assertEqual(len(result.aux_operator_eigenvalues), 2) + # expectation values + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0][0], 2.0, places=1) + self.assertAlmostEqual( + result.aux_operator_eigenvalues[0][1][0], 0.0019531249999999445, places=1 + ) + # standard deviations + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0][1], 0.0) + self.assertAlmostEqual( + result.aux_operator_eigenvalues[0][1][1], 0.015183867579396111, places=1 + ) + + # Go again with additional None and zero operators + aux_ops = [*aux_ops, None, 0] + result = ssvqe.compute_eigenvalues(self.h2_op, aux_operators=aux_ops) + self.assertEqual(len(result.aux_operator_eigenvalues[0]), 4) + # expectation values + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0][0], 2.0, places=1) + self.assertAlmostEqual( + result.aux_operator_eigenvalues[0][1][0], 0.0019531249999999445, places=1 + ) + self.assertEqual(result.aux_operator_eigenvalues[0][2][0], 0.0) + self.assertEqual(result.aux_operator_eigenvalues[0][3][0], 0.0) + # # standard deviations + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0][1], 0.0) + self.assertAlmostEqual( + result.aux_operator_eigenvalues[0][1][1], 0.01548658094658011, places=1 + ) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][2][1], 0.0) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][3][1], 0.0) + + @unittest.skipUnless(has_aer(), "qiskit-aer doesn't appear to be installed.") + def test_aux_operator_std_dev_aer_pauli(self): + """Test non-zero standard deviations of aux operators with AerPauliExpectation.""" + wavefunction = self.ry_wavefunction + ssvqe = SSVQE( + ansatz=wavefunction, + expectation=AerPauliExpectation(), + optimizer=COBYLA(maxiter=0), + quantum_instance=QuantumInstance( + backend=Aer.get_backend("qasm_simulator"), + shots=1, + seed_simulator=algorithm_globals.random_seed, + seed_transpiler=algorithm_globals.random_seed, + ), + ) + + # Go again 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 = ssvqe.compute_eigenvalues(self.h2_op, aux_operators=aux_ops) + self.assertEqual(len(result.aux_operator_eigenvalues), 2) + # expectation values + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0][0], 2.0, places=1) + self.assertAlmostEqual( + result.aux_operator_eigenvalues[0][1][0], 0.6698863565455391, places=1 + ) + # standard deviations + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0][1], 0.0) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][1][1], 0.0, places=6) + + # Go again with additional None and zero operators + aux_ops = [*aux_ops, None, 0] + result = ssvqe.compute_eigenvalues(self.h2_op, aux_operators=aux_ops) + self.assertEqual(len(result.aux_operator_eigenvalues[-1]), 4) + # expectation values + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0][0], 2.0, places=6) + self.assertAlmostEqual( + result.aux_operator_eigenvalues[0][1][0], 0.6036400943063891, places=6 + ) + self.assertEqual(result.aux_operator_eigenvalues[0][2][0], 0.0) + self.assertEqual(result.aux_operator_eigenvalues[0][3][0], 0.0) + # standard deviations + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0][1], 0.0) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][1][1], 0.0, places=6) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][2][1], 0.0) + self.assertAlmostEqual(result.aux_operator_eigenvalues[0][3][1], 0.0) + + +if __name__ == "__main__": + unittest.main()