diff --git a/.pylintdict b/.pylintdict index b667eec08..1dae4034b 100644 --- a/.pylintdict +++ b/.pylintdict @@ -6,7 +6,9 @@ aer's al annealers ansatz +ansatzes apidocs +apl applegate args arxiv @@ -15,7 +17,9 @@ autosummary backend backends barkoutsos +batchsize benchmarking +bfgs bitstring bitstrings bixby @@ -23,6 +27,7 @@ bool boolean boyd bravyi +callables catol cartan chu @@ -35,16 +40,20 @@ conv const cplex cplexoptimizer +crs cts currentmodule cvar +dataclass deterministically dic dicts +dimensionality disp dmitry docplex docplex's +docstring doi dok dtype @@ -55,22 +64,31 @@ eigen eigensolver eigensolvers eigenstate +eigenstates embeddings entangler enum eq et +eval +evals exponentiated +failsafe farhi fmin formatter +frac func +functools fred fval fx +f'spsa gambella geq getter +getters +globals glover goemans goldstone @@ -81,11 +99,13 @@ gurobi gurobioptimizer gurobipy gutmann +hadfield hamilton hamiltonian hamiltonians hastings hayashi +hessians hoyer https ibm @@ -99,7 +119,9 @@ iprint ising iter iteratively +jac july +kandala karimi kirkpatrick kwargs @@ -108,22 +130,28 @@ len leq lhs lin +linalg linearconstraint linexpr +linter lowerbound lp +lse lucas macos makefile +marecek masahito matplotlib maxcut +maxfev maxfun maxiter mdl milp minimizer minimumeigenoptimizer +modelspace mmp mpm multiset @@ -132,9 +160,13 @@ nannicini natively ndarray ndarrays +nones noop +nelder networkx neven +nfev +nft nosignatures np num @@ -147,18 +179,23 @@ optimality optimizationresult optimizationresultstatus optimizers +packagebut panchenko param +parameterizations params parikh +passmanager pauli paulis peleato pmm +polyfit pooya pos ppp pre +preconditioner preprint prepend presolver @@ -187,6 +224,8 @@ qubo readme repr representable +resamplings +rescaling rhobeg rhoend rhs @@ -196,22 +235,30 @@ robert ronagh rtype runtime +rustworkx ry rz sahar scipy sdp +serializable sherrington simonetto slsqp smode smoothen +spall +spedalieri spsa src statevector stdout stephen +steppable +stepsize str +subclassed +subclasses subcollection subgraph submodules @@ -222,6 +269,7 @@ summands tavernelli terra th +tnc toctree todok tol @@ -241,13 +289,17 @@ variational vartype vqe vqeresult +utils writelines +xatol xixj +xopt wavefunction wecker whitespace wiesner williamson +woerner xs ys zemlin diff --git a/qiskit_optimization/__init__.py b/qiskit_optimization/__init__.py index db2b7281a..f1d42f11c 100644 --- a/qiskit_optimization/__init__.py +++ b/qiskit_optimization/__init__.py @@ -85,9 +85,15 @@ """ -from .exceptions import QiskitOptimizationError +from .exceptions import QiskitOptimizationError, AlgorithmError from .infinity import INFINITY # must be at the top of the file from .problems.quadratic_program import QuadraticProgram from .version import __version__ -__all__ = ["__version__", "QuadraticProgram", "QiskitOptimizationError", "INFINITY"] +__all__ = [ + "__version__", + "QuadraticProgram", + "QiskitOptimizationError", + "AlgorithmError", + "INFINITY", +] diff --git a/qiskit_optimization/algorithm_job.py b/qiskit_optimization/algorithm_job.py new file mode 100644 index 000000000..abd6def46 --- /dev/null +++ b/qiskit_optimization/algorithm_job.py @@ -0,0 +1,45 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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. + +""" +AlgorithmJob class +""" +from qiskit.primitives.primitive_job import PrimitiveJob + + +class AlgorithmJob(PrimitiveJob): + """ + This class is introduced for typing purposes and provides no + additional function beyond that inherited from its parents. + + Update: :meth:`AlgorithmJob.submit()` method added. See its + documentation for more info. + """ + + def submit(self) -> None: + """ + Submit the job for execution. + + For V1 primitives, Qiskit ``PrimitiveJob`` subclassed JobV1 and defined ``submit()``. + ``PrimitiveJob`` was updated for V2 primitives, no longer subclasses ``JobV1``, and + now has a private ``_submit()`` method, with ``submit()`` being deprecated as of + Qiskit version 0.46. This maintains the ``submit()`` for ``AlgorithmJob`` here as + it's called in many places for such a job. An alternative could be to make + 0.46 the required minimum version and alter all algorithm's call sites to use + ``_submit()`` and make this an empty class again as it once was. For now this + way maintains compatibility with the current min version of 0.44. + """ + # TODO: Considering changing this in the future - see above docstring. + try: + super()._submit() + except AttributeError: + super().submit() # pylint: disable=no-member diff --git a/qiskit_optimization/algorithm_result.py b/qiskit_optimization/algorithm_result.py new file mode 100644 index 000000000..95b45d829 --- /dev/null +++ b/qiskit_optimization/algorithm_result.py @@ -0,0 +1,65 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2020, 2024. +# +# 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. + +""" +This module implements the abstract base class for algorithm results. +""" + +from abc import ABC +import inspect +import pprint + + +class AlgorithmResult(ABC): + """Abstract Base Class for algorithm results.""" + + def __str__(self) -> str: + result = {} + for name, value in inspect.getmembers(self): + if ( + not name.startswith("_") + and not inspect.ismethod(value) + and not inspect.isfunction(value) + and hasattr(self, name) + ): + + result[name] = value + + return pprint.pformat(result, indent=4) + + def combine(self, result: "AlgorithmResult") -> None: + """ + Any property from the argument that exists in the receiver is + updated. + Args: + result: Argument result with properties to be set. + Raises: + TypeError: Argument is None + """ + if result is None: + raise TypeError("Argument result expected.") + if result == self: + return + + # find any result public property that exists in the receiver + for name, value in inspect.getmembers(result): + if ( + not name.startswith("_") + and not inspect.ismethod(value) + and not inspect.isfunction(value) + and hasattr(self, name) + ): + try: + setattr(self, name, value) + except AttributeError: + # some attributes may be read only + pass diff --git a/qiskit_optimization/algorithms/minimum_eigen_optimizer.py b/qiskit_optimization/algorithms/minimum_eigen_optimizer.py index fc0b24262..67d72d0bf 100644 --- a/qiskit_optimization/algorithms/minimum_eigen_optimizer.py +++ b/qiskit_optimization/algorithms/minimum_eigen_optimizer.py @@ -15,7 +15,7 @@ import numpy as np from qiskit.quantum_info import SparsePauliOp -from qiskit_algorithms import ( +from ..minimum_eigensolvers import ( NumPyMinimumEigensolver, NumPyMinimumEigensolverResult, SamplingMinimumEigensolver, diff --git a/qiskit_optimization/eigensolvers/__init__.py b/qiskit_optimization/eigensolvers/__init__.py new file mode 100644 index 000000000..a60bd2520 --- /dev/null +++ b/qiskit_optimization/eigensolvers/__init__.py @@ -0,0 +1,23 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 Eigensolvers package""" + +from .numpy_eigensolver import NumPyEigensolver, NumPyEigensolverResult +from .eigensolver import Eigensolver, EigensolverResult + +__all__ = [ + "NumPyEigensolver", + "NumPyEigensolverResult", + "Eigensolver", + "EigensolverResult", +] diff --git a/qiskit_optimization/eigensolvers/eigensolver.py b/qiskit_optimization/eigensolvers/eigensolver.py new file mode 100644 index 000000000..74fcc7d4a --- /dev/null +++ b/qiskit_optimization/eigensolvers/eigensolver.py @@ -0,0 +1,103 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 eigensolver interface and result.""" + +from __future__ import annotations + +from abc import ABC, abstractmethod +from typing import Any +import numpy as np + +from qiskit.quantum_info.operators.base_operator import BaseOperator + +from ..algorithm_result import AlgorithmResult +from ..list_or_dict import ListOrDict + + +class Eigensolver(ABC): + """The eigensolver interface. + + Algorithms that can compute eigenvalues for an operator + may implement this interface to allow different algorithms to be + used interchangeably. + """ + + @abstractmethod + def compute_eigenvalues( + self, + operator: BaseOperator, + aux_operators: ListOrDict[BaseOperator] | None = None, + ) -> "EigensolverResult": + """ + Computes the minimum eigenvalue. The ``operator`` and ``aux_operators`` are supplied here. + 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 + eigenstate 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: + An eigensolver result. + """ + return EigensolverResult() + + @classmethod + def supports_aux_operators(cls) -> bool: + """Whether computing the expectation value of auxiliary operators is supported. + + If the eigensolver computes the eigenvalues 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 EigensolverResult(AlgorithmResult): + """Eigensolver result.""" + + def __init__(self) -> None: + super().__init__() + self._eigenvalues: np.ndarray | None = None + self._aux_operators_evaluated: list[ListOrDict[tuple[float, dict[str, Any]]]] | None = None + + @property + def eigenvalues(self) -> np.ndarray | None: + """Return the eigenvalues.""" + return self._eigenvalues + + @eigenvalues.setter + def eigenvalues(self, value: np.ndarray) -> None: + """Set the eigenvalues.""" + self._eigenvalues = value + + @property + def aux_operators_evaluated( + self, + ) -> list[ListOrDict[tuple[float, dict[str, Any]]]] | None: + """Return the aux operator expectation values. + + These values are in fact tuples formatted as (mean, metadata). + """ + return self._aux_operators_evaluated + + @aux_operators_evaluated.setter + def aux_operators_evaluated( + self, value: list[ListOrDict[tuple[float, dict[str, Any]]]] + ) -> None: + """Set the aux operator eigenvalues.""" + self._aux_operators_evaluated = value diff --git a/qiskit_optimization/eigensolvers/numpy_eigensolver.py b/qiskit_optimization/eigensolvers/numpy_eigensolver.py new file mode 100644 index 000000000..9ebaa1bfe --- /dev/null +++ b/qiskit_optimization/eigensolvers/numpy_eigensolver.py @@ -0,0 +1,320 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 eigensolver algorithm.""" + +from __future__ import annotations + +from collections.abc import Iterable +from typing import Callable, Union, Tuple, Dict, List, Optional, cast +import logging +import numpy as np +from scipy import sparse as scisparse + +from qiskit.quantum_info import SparsePauliOp, Statevector +from qiskit.quantum_info.operators.base_operator import BaseOperator + +from ..utils.validation import validate_min +from .eigensolver import Eigensolver, EigensolverResult +from ..exceptions import AlgorithmError +from ..list_or_dict import ListOrDict + +logger = logging.getLogger(__name__) + +FilterType = Callable[ + [Union[List, np.ndarray], float, Optional[ListOrDict[Tuple[float, Dict[str, float]]]]], bool +] + + +class NumPyEigensolver(Eigensolver): + r""" + The NumPy eigensolver algorithm. + + The NumPy Eigensolver computes up to the first :math:`k` eigenvalues of a complex-valued square + matrix of dimension :math:`n \times n`, with :math:`k \leq n`. + + Note: + Operators are automatically converted to SciPy's ``spmatrix`` + as needed and this conversion can be costly in terms of memory and performance as the + operator size, mostly in terms of number of qubits it represents, gets larger. + """ + + def __init__( + self, + k: int = 1, + filter_criterion: FilterType | None = None, + ) -> None: + """ + Args: + k: Number of eigenvalues are to be computed, with a minimum value of 1. + filter_criterion: Callable that allows to filter eigenvalues/eigenstates. Only feasible + eigenstates are returned in the results. The callable has the signature + ``filter(eigenstate, eigenvalue, aux_values)`` and must return a boolean to indicate + whether to keep this value in the final returned result or not. If the number of + elements that satisfies the criterion is smaller than ``k``, then the returned list will + have fewer elements and can even be empty. + """ + validate_min("k", k, 1) + super().__init__() + + self._in_k = k + self._k = k # pylint: disable=invalid-name + + self._filter_criterion = filter_criterion + + @property + def k(self) -> int: + """Return k (number of eigenvalues requested).""" + return self._in_k + + @k.setter + def k(self, k: int) -> None: + """Set k (number of eigenvalues requested).""" + validate_min("k", k, 1) + self._in_k = k + self._k = k + + @property + def filter_criterion( + self, + ) -> FilterType | None: + """Return the filter criterion if set.""" + return self._filter_criterion + + @filter_criterion.setter + def filter_criterion(self, filter_criterion: FilterType | None) -> None: + """Set the filter criterion.""" + self._filter_criterion = filter_criterion + + @classmethod + def supports_aux_operators(cls) -> bool: + return True + + def _check_set_k(self, operator: BaseOperator) -> None: + if operator is not None: + if self._in_k > 2**operator.num_qubits: + self._k = 2**operator.num_qubits + logger.debug( + "WARNING: Asked for %s eigenvalues but max possible is %s.", self._in_k, self._k + ) + else: + self._k = self._in_k + + def _solve(self, operator: BaseOperator) -> tuple[np.ndarray, np.ndarray]: + + try: + op_matrix = operator.to_matrix(sparse=True) + except TypeError: + logger.debug( + "WARNING: operator of type `%s` does not support sparse matrices. " + "Trying dense computation", + type(operator), + ) + try: + op_matrix = operator.to_matrix() + except AttributeError as ex: + raise AlgorithmError(f"Unsupported operator type `{type(operator)}`.") from ex + + if isinstance(op_matrix, scisparse.csr_matrix): + # If matrix is diagonal, the elements on the diagonal are the eigenvalues. Solve by sorting. + if scisparse.csr_matrix(op_matrix.diagonal()).nnz == op_matrix.nnz: + diag = op_matrix.diagonal() + indices = np.argsort(diag)[: self._k] + eigval = diag[indices] + eigvec = np.zeros((op_matrix.shape[0], self._k)) + for i, idx in enumerate(indices): + eigvec[idx, i] = 1.0 + else: + if self._k >= 2**operator.num_qubits - 1: + logger.debug( + "SciPy doesn't support to get all eigenvalues, using NumPy instead." + ) + eigval, eigvec = self._solve_dense(operator.to_matrix()) + else: + eigval, eigvec = self._solve_sparse(op_matrix, self._k) + else: + # Sparse SciPy matrix not supported, use dense NumPy computation. + eigval, eigvec = self._solve_dense(operator.to_matrix()) + + indices = np.argsort(eigval)[: self._k] + eigval = eigval[indices] + eigvec = eigvec[:, indices] + return eigval, eigvec.T + + @staticmethod + def _solve_sparse(op_matrix: scisparse.csr_matrix, k: int) -> tuple[np.ndarray, np.ndarray]: + if (op_matrix != op_matrix.getH()).nnz == 0: + # Operator is Hermitian + return scisparse.linalg.eigsh(op_matrix, k=k, which="SA") + else: + return scisparse.linalg.eigs(op_matrix, k=k, which="SR") + + @staticmethod + def _solve_dense(op_matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + if op_matrix.all() == op_matrix.conj().T.all(): + # Operator is Hermitian + return cast(Tuple[np.ndarray, np.ndarray], np.linalg.eigh(op_matrix)) + else: + return cast(Tuple[np.ndarray, np.ndarray], np.linalg.eig(op_matrix)) + + @staticmethod + def _eval_aux_operators( + aux_operators: ListOrDict[BaseOperator], + wavefn: np.ndarray, + threshold: float = 1e-12, + ) -> ListOrDict[tuple[float, dict[str, float]]]: + + values: ListOrDict[tuple[float, dict[str, float]]] + + # As a list, aux_operators can contain None operators for which None values are returned. + # As a dict, the None operators in aux_operators have been dropped in compute_eigenvalues. + key_op_iterator: Iterable[tuple[str | int, BaseOperator]] + if isinstance(aux_operators, list): + values = [None] * len(aux_operators) + key_op_iterator = enumerate(aux_operators) + else: + values = {} + key_op_iterator = aux_operators.items() + + for key, operator in key_op_iterator: + if operator is None: + continue + + if operator.num_qubits is None or operator.num_qubits < 1: + logger.info( + "The number of qubits of the %s operator must be greater than zero.", key + ) + continue + + op_matrix = None + try: + op_matrix = operator.to_matrix(sparse=True) + except TypeError: + logger.debug( + "WARNING: operator of type `%s` does not support sparse matrices. " + "Trying dense computation", + type(operator), + ) + try: + op_matrix = operator.to_matrix() + except AttributeError as ex: + raise AlgorithmError(f"Unsupported operator type {type(operator)}.") from ex + + if isinstance(op_matrix, scisparse.csr_matrix): + value = op_matrix.dot(wavefn).dot(np.conj(wavefn)) + elif isinstance(op_matrix, np.ndarray): + value = Statevector(wavefn).expectation_value(operator) + else: + value = 0.0 + + value = value if np.abs(value) > threshold else 0.0 + # The value gets wrapped into a tuple: (mean, metadata). + # The metadata includes variance (and, for other eigensolvers, shots). + # Since this is an exact computation, there are no shots + # and the variance is known to be zero. + values[key] = (value, {"variance": 0.0}) # type: ignore[index] + return values + + def compute_eigenvalues( + self, + operator: BaseOperator, + aux_operators: ListOrDict[BaseOperator] | None = None, + ) -> NumPyEigensolverResult: + + super().compute_eigenvalues(operator, aux_operators) + + if operator.num_qubits is None or operator.num_qubits < 1: + raise AlgorithmError("The number of qubits of the operator must be greater than zero.") + + self._check_set_k(operator) + + zero_op = SparsePauliOp(["I" * operator.num_qubits], coeffs=[0.0]) + if isinstance(aux_operators, list) and len(aux_operators) > 0: + # For some reason Chemistry passes aux_ops with 0 qubits and paulis sometimes. + aux_operators = [zero_op if op == 0 else op for op in aux_operators] + elif isinstance(aux_operators, dict) and len(aux_operators) > 0: + aux_operators = { + key: zero_op if op == 0 else op # Convert zero values to zero operators + for key, op in aux_operators.items() + if op is not None # Discard None values + } + else: + aux_operators = None + + k_orig = self._k + if self._filter_criterion: + # need to consider all elements if a filter is set + self._k = 2**operator.num_qubits + + eigvals, eigvecs = self._solve(operator) + + # compute energies before filtering, as this also evaluates the aux operators + if aux_operators is not None: + aux_op_vals = [ + self._eval_aux_operators(aux_operators, eigvecs[i]) for i in range(self._k) + ] + else: + aux_op_vals = None + + # if a filter is set, loop over the given values and only keep + if self._filter_criterion: + filt_eigvals = [] + filt_eigvecs = [] + filt_aux_op_vals = [] + count = 0 + for i, (eigval, eigvec) in enumerate(zip(eigvals, eigvecs)): + if aux_op_vals is not None: + aux_op_val = aux_op_vals[i] + else: + aux_op_val = None + + if self._filter_criterion(eigvec, eigval, aux_op_val): + count += 1 + filt_eigvecs.append(eigvec) + filt_eigvals.append(eigval) + if aux_op_vals is not None: + filt_aux_op_vals.append(aux_op_val) + + if count == k_orig: + break + + eigvals = np.array(filt_eigvals) + eigvecs = np.array(filt_eigvecs) + aux_op_vals = filt_aux_op_vals + + self._k = k_orig + + result = NumPyEigensolverResult() + result.eigenvalues = eigvals + result.eigenstates = [Statevector(vec) for vec in eigvecs] + result.aux_operators_evaluated = aux_op_vals + + logger.debug("NumpyEigensolverResult:\n%s", result) + return result + + +class NumPyEigensolverResult(EigensolverResult): + """NumPy eigensolver result.""" + + def __init__(self) -> None: + super().__init__() + self._eigenstates: list[Statevector] | None = None + + @property + def eigenstates(self) -> list[Statevector] | None: + """Return eigenstates.""" + return self._eigenstates + + @eigenstates.setter + def eigenstates(self, value: list[Statevector]) -> None: + """Set eigenstates.""" + self._eigenstates = value diff --git a/qiskit_optimization/exceptions.py b/qiskit_optimization/exceptions.py index 7b6c654f4..536bf4c4a 100644 --- a/qiskit_optimization/exceptions.py +++ b/qiskit_optimization/exceptions.py @@ -1,6 +1,6 @@ # This code is part of a Qiskit project. # -# (C) Copyright IBM 2019, 2023. +# (C) Copyright IBM 2019, 2024. # # 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 @@ -19,3 +19,9 @@ class QiskitOptimizationError(QiskitError): """Class for errors returned by Qiskit optimization module.""" pass + + +class AlgorithmError(QiskitError): + """For Algorithm specific errors.""" + + pass diff --git a/qiskit_optimization/list_or_dict.py b/qiskit_optimization/list_or_dict.py new file mode 100644 index 000000000..ead479bce --- /dev/null +++ b/qiskit_optimization/list_or_dict.py @@ -0,0 +1,18 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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. + +"""Introduced new type to maintain readability.""" + +from typing import TypeVar, List, Union, Optional, Dict + +_T = TypeVar("_T") # Pylint does not allow single character class names. +ListOrDict = Union[List[Optional[_T]], Dict[str, _T]] diff --git a/qiskit_optimization/minimum_eigensolvers/__init__.py b/qiskit_optimization/minimum_eigensolvers/__init__.py new file mode 100644 index 000000000..a331d64fe --- /dev/null +++ b/qiskit_optimization/minimum_eigensolvers/__init__.py @@ -0,0 +1,30 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 Eigensolvers package.""" + +from .sampling_mes import SamplingMinimumEigensolver, SamplingMinimumEigensolverResult +from .minimum_eigensolver import MinimumEigensolver, MinimumEigensolverResult +from .numpy_minimum_eigensolver import NumPyMinimumEigensolver, NumPyMinimumEigensolverResult +from .qaoa import QAOA +from .sampling_vqe import SamplingVQE + +__all__ = [ + "SamplingMinimumEigensolver", + "SamplingMinimumEigensolverResult", + "MinimumEigensolver", + "MinimumEigensolverResult", + "NumPyMinimumEigensolver", + "NumPyMinimumEigensolverResult", + "SamplingVQE", + "QAOA", +] diff --git a/qiskit_optimization/minimum_eigensolvers/diagonal_estimator.py b/qiskit_optimization/minimum_eigensolvers/diagonal_estimator.py new file mode 100644 index 000000000..e6b162350 --- /dev/null +++ b/qiskit_optimization/minimum_eigensolvers/diagonal_estimator.py @@ -0,0 +1,203 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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. + +"""Expectation value for a diagonal observable using a sampler primitive.""" + +from __future__ import annotations + +from collections.abc import Callable, Sequence, Mapping, Iterable, MappingView +from typing import Any + +from dataclasses import dataclass + +import numpy as np +from qiskit.circuit import QuantumCircuit +from qiskit.primitives import BaseSampler, BaseEstimator, EstimatorResult +from qiskit.primitives.utils import init_observable, _circuit_key +from qiskit.quantum_info import SparsePauliOp +from qiskit.quantum_info.operators.base_operator import BaseOperator + +from ..algorithm_job import AlgorithmJob + + +@dataclass(frozen=True) +class _DiagonalEstimatorResult(EstimatorResult): + """A result from an expectation of a diagonal observable.""" + + # TODO make each measurement a dataclass rather than a dict + best_measurements: Sequence[Mapping[str, Any]] | None = None + + +class _DiagonalEstimator(BaseEstimator): + """An estimator for diagonal observables.""" + + def __init__( + self, + sampler: BaseSampler, + aggregation: float | Callable[[Iterable[tuple[float, float]]], float] | None = None, + callback: Callable[[Sequence[Mapping[str, Any]]], None] | None = None, + **options, + ) -> None: + r"""Evaluate the expectation of quantum state with respect to a diagonal operator. + + Args: + sampler: The sampler used to evaluate the circuits. + aggregation: The aggregation function to aggregate the measurement outcomes. If a float + this specified the CVaR :math:`\alpha` parameter. + callback: A callback which is given the best measurements of all circuits in each + evaluation. + run_options: Options for the sampler. + + """ + super().__init__(options=options) + self._circuits: list[QuantumCircuit] = [] # See Qiskit pull request 11051 + self._parameters: list[MappingView] = [] + self._observables: list[SparsePauliOp] = [] + + self.sampler = sampler + if not callable(aggregation): + aggregation = _get_cvar_aggregation(aggregation) + + self.aggregation = aggregation + self.callback = callback + self._circuit_ids: dict[int, QuantumCircuit] = {} + self._observable_ids: dict[int, BaseOperator] = {} + + def _run( + self, + circuits: Sequence[QuantumCircuit], + observables: Sequence[BaseOperator], + parameter_values: Sequence[Sequence[float]], + **run_options, + ) -> AlgorithmJob: + circuit_indices = [] + for circuit in circuits: + key = _circuit_key(circuit) + index = self._circuit_ids.get(key) + if index is not None: + circuit_indices.append(index) + else: + circuit_indices.append(len(self._circuits)) + self._circuit_ids[key] = len(self._circuits) + self._circuits.append(circuit) + self._parameters.append(circuit.parameters) + observable_indices = [] + for observable in observables: + index = self._observable_ids.get(id(observable)) + if index is not None: + observable_indices.append(index) + else: + observable_indices.append(len(self._observables)) + self._observable_ids[id(observable)] = len(self._observables) + converted_observable = init_observable(observable) + _check_observable_is_diagonal(converted_observable) # check it's diagonal + self._observables.append(converted_observable) + job = AlgorithmJob( + self._call, circuit_indices, observable_indices, parameter_values, **run_options + ) + job.submit() + return job + + def _call( + self, + circuits: Sequence[int], + observables: Sequence[int], + parameter_values: Sequence[Sequence[float]], + **run_options, + ) -> _DiagonalEstimatorResult: + job = self.sampler.run( + [self._circuits[i] for i in circuits], + parameter_values, + **run_options, + ) + sampler_result = job.result() + samples = sampler_result.quasi_dists + + # a list of dictionaries containing: {state: (measurement probability, value)} + evaluations: list[dict[int, tuple[float, float]]] = [ + { + state: (probability, _evaluate_sparsepauli(state, self._observables[i])) + for state, probability in sampled.items() + } + for i, sampled in zip(observables, samples) + ] + + results = np.array([self.aggregation(evaluated.values()) for evaluated in evaluations]) + + # get the best measurements + best_measurements = [] + num_qubits = self._circuits[0].num_qubits + for evaluated in evaluations: + best_result = min(evaluated.items(), key=lambda x: x[1][1]) + best_measurements.append( + { + "state": best_result[0], + "bitstring": bin(best_result[0])[2:].zfill(num_qubits), + "value": best_result[1][1], + "probability": best_result[1][0], + } + ) + + if self.callback is not None: + self.callback(best_measurements) + + return _DiagonalEstimatorResult( + values=results, metadata=sampler_result.metadata, best_measurements=best_measurements + ) + + +def _get_cvar_aggregation(alpha: float | None) -> Callable[[Iterable[tuple[float, float]]], float]: + """Get the aggregation function for CVaR with confidence level ``alpha``.""" + if alpha is None: + alpha = 1 + elif not 0 <= alpha <= 1: + raise ValueError(f"alpha must be in [0, 1] but was {alpha}") + + # if alpha is close to 1 we can avoid the sorting + if np.isclose(alpha, 1): + + def aggregate(measurements: Iterable[tuple[float, float]]) -> float: + return sum(probability * value for probability, value in measurements) + + else: + + def aggregate(measurements: Iterable[tuple[float, float]]) -> float: + # sort by values + sorted_measurements = sorted(measurements, key=lambda x: x[1]) + + accumulated_percent = 0.0 # once alpha is reached, stop + cvar = 0.0 + for probability, value in sorted_measurements: + cvar += value * min(probability, alpha - accumulated_percent) + accumulated_percent += probability + if accumulated_percent >= alpha: + break + + return cvar / alpha + + return aggregate + + +_PARITY = np.array([-1 if bin(i).count("1") % 2 else 1 for i in range(256)], dtype=np.complex128) + + +def _evaluate_sparsepauli(state: int, observable: SparsePauliOp) -> float: + packed_uint8 = np.packbits(observable.paulis.z, axis=1, bitorder="little") + state_bytes = np.frombuffer(state.to_bytes(packed_uint8.shape[1], "little"), dtype=np.uint8) + reduced = np.bitwise_xor.reduce(packed_uint8 & state_bytes, axis=1) + return np.sum(observable.coeffs * _PARITY[reduced]) + + +def _check_observable_is_diagonal(observable: SparsePauliOp) -> None: + is_diagonal = not np.any(observable.paulis.x) + if not is_diagonal: + raise ValueError("The observable must be diagonal.") diff --git a/qiskit_optimization/minimum_eigensolvers/minimum_eigensolver.py b/qiskit_optimization/minimum_eigensolvers/minimum_eigensolver.py new file mode 100644 index 000000000..78c05d9c8 --- /dev/null +++ b/qiskit_optimization/minimum_eigensolvers/minimum_eigensolver.py @@ -0,0 +1,96 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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.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, + aux_operators: ListOrDict[BaseOperator] | None = None, + ) -> "MinimumEigensolverResult": + """ + Computes the minimum eigenvalue. The ``operator`` and ``aux_operators`` are supplied here. + 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: complex | None = None + self._aux_operators_evaluated: ListOrDict[tuple[complex, dict[str, Any]]] | None = 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 diff --git a/qiskit_optimization/minimum_eigensolvers/numpy_minimum_eigensolver.py b/qiskit_optimization/minimum_eigensolvers/numpy_minimum_eigensolver.py new file mode 100644 index 000000000..5458d5e50 --- /dev/null +++ b/qiskit_optimization/minimum_eigensolvers/numpy_minimum_eigensolver.py @@ -0,0 +1,109 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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, Union, Tuple, Dict, List, Optional +import logging +import numpy as np + +from qiskit.quantum_info import Statevector +from qiskit.quantum_info.operators.base_operator import BaseOperator + +from ..eigensolvers.numpy_eigensolver 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[Tuple[float, Dict[str, 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, + aux_operators: ListOrDict[BaseOperator] | 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_operators_evaluated: + result.aux_operators_evaluated = eigensolver_result.aux_operators_evaluated[ + 0 + ] # type: ignore[assignment] + + 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: Statevector | None = None + + @property + def eigenstate(self) -> Statevector | None: + """Returns the eigenstate corresponding to the computed minimum eigenvalue.""" + return self._eigenstate + + @eigenstate.setter + def eigenstate(self, value: Statevector) -> None: + self._eigenstate = value diff --git a/qiskit_optimization/minimum_eigensolvers/qaoa.py b/qiskit_optimization/minimum_eigensolvers/qaoa.py new file mode 100644 index 000000000..849ee5d96 --- /dev/null +++ b/qiskit_optimization/minimum_eigensolvers/qaoa.py @@ -0,0 +1,136 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 quantum approximate optimization algorithm.""" + +from __future__ import annotations + +from typing import Any, Callable + +import numpy as np +from qiskit.circuit import QuantumCircuit +from qiskit.circuit.library.n_local.qaoa_ansatz import QAOAAnsatz +from qiskit.passmanager import BasePassManager +from qiskit.primitives import BaseSamplerV1, BaseSamplerV2 +from qiskit.quantum_info.operators.base_operator import BaseOperator +from ..optimizers.optimizer import Minimizer, Optimizer +from ..utils.validation import validate_min +from .sampling_vqe import SamplingVQE + + +class QAOA(SamplingVQE): + r""" + The Quantum Approximate Optimization Algorithm (QAOA). + QAOA is a well-known algorithm for finding approximate solutions to combinatorial-optimization + problems [1]. + The QAOA implementation directly extends :class:`.SamplingVQE` and inherits its optimization + structure. However, unlike VQE, which can be configured with arbitrary ansatzes, QAOA uses its + own fine-tuned ansatz, which comprises :math:`p` parameterized global :math:`x` rotations and + :math:`p` different parameterizations of the problem hamiltonian. QAOA is thus principally + configured by the single integer parameter, ``reps``, which dictates the depth of the ansatz, + and thus affects the approximation quality. + An optional array of :math:`2p` parameter values, as the :attr:`initial_point`, may be provided + as the starting :math:`\beta` and :math:`\gamma` parameters for the QAOA ansatz [1]. + An operator or a parameterized quantum circuit may optionally also be provided as a custom + :attr:`mixer` Hamiltonian. This allows in the case of quantum annealing [2] and QAOA [3], to run + constrained optimization problems where the mixer constrains the evolution to a feasible + subspace of the full Hilbert space. + The following attributes can be set via the initializer but can also be read and updated once + the QAOA object has been constructed. + Attributes: + sampler (BaseSampler): The sampler primitive to sample the circuits. + optimizer (Optimizer | Minimizer): A classical optimizer to find the minimum energy. This + can either be an :class:`.Optimizer` or a callable implementing the + :class:`.Minimizer` protocol. + reps (int): The integer parameter :math:`p`. Has a minimum valid value of 1. + initial_state: An optional initial state to prepend the QAOA circuit with. + mixer (QuantumCircuit | BaseOperator): The mixer Hamiltonian to evolve with or + a custom quantum circuit. Allows support of optimizations in constrained subspaces [2, + 3] as well as warm-starting the optimization [4]. + aggregation (float | Callable[[list[float]], float] | None): A float or callable to specify + how the objective function evaluated on the basis states should be aggregated. If a + float, this specifies the :math:`\alpha \in [0,1]` parameter for a CVaR expectation + value. + 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 value, and + the metadata dictionary. + References: + [1]: Farhi, E., Goldstone, J., Gutmann, S., "A Quantum Approximate Optimization Algorithm" + `arXiv:1411.4028 `__ + [2]: Hen, I., Spedalieri, F. M., "Quantum Annealing for Constrained Optimization" + `PhysRevApplied.5.034007 `__ + [3]: Hadfield, S. et al, "From the Quantum Approximate Optimization Algorithm to a Quantum + Alternating Operator Ansatz" `arXiv:1709.03489 `__ + [4]: Egger, D. J., Marecek, J., Woerner, S., "Warm-starting quantum optimization" + `arXiv: 2009.10095 `__ + """ + + def __init__( + self, + sampler: BaseSamplerV1 | BaseSamplerV2, + optimizer: Optimizer | Minimizer, + *, + reps: int = 1, + initial_state: QuantumCircuit | None = None, + mixer: QuantumCircuit | BaseOperator = None, + initial_point: np.ndarray | None = None, + aggregation: float | Callable[[list[float]], float] | None = None, + callback: Callable[[int, np.ndarray, float, dict[str, Any]], None] | None = None, + passmanager: BasePassManager | None = None, + ) -> None: + r""" + Args: + sampler: The sampler primitive to sample the circuits. + optimizer: A classical optimizer to find the minimum energy. This can either be + an :class:`.Optimizer` or a callable implementing the :class:`.Minimizer` + protocol. + reps: The integer parameter :math:`p`. Has a minimum valid value of 1. + initial_state: An optional initial state to prepend the QAOA circuit with. + mixer: The mixer Hamiltonian to evolve with or a custom quantum circuit. Allows support + of optimizations in constrained subspaces [2, 3] as well as warm-starting the + optimization [4]. + 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. ``QAOA`` 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. + aggregation: A float or callable to specify how the objective function evaluated on the + basis states should be aggregated. If a float, this specifies the :math:`\alpha \in + [0,1]` parameter for a CVaR expectation value. + 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 + evaluated value, the metadata dictionary. + passmanager: A pass manager to transpile the circuits. + """ + validate_min("reps", reps, 1) + + self.reps = reps + self.mixer = mixer + self.initial_state = initial_state + self._cost_operator = None + + super().__init__( + sampler=sampler, + ansatz=None, + optimizer=optimizer, + initial_point=initial_point, + aggregation=aggregation, + callback=callback, + passmanager=passmanager, + ) + + def _check_operator_ansatz(self, operator: BaseOperator): + # Recreates a circuit based on operator parameter. + self.ansatz = QAOAAnsatz( + operator, self.reps, initial_state=self.initial_state, mixer_operator=self.mixer + ).decompose() # TODO remove decompose once #6674 is fixed <-- I don't know what this issue is diff --git a/qiskit_optimization/minimum_eigensolvers/sampling_mes.py b/qiskit_optimization/minimum_eigensolvers/sampling_mes.py new file mode 100644 index 000000000..1a6749d27 --- /dev/null +++ b/qiskit_optimization/minimum_eigensolvers/sampling_mes.py @@ -0,0 +1,125 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 Sampling Minimum Eigensolver interface.""" + +from __future__ import annotations +from abc import ABC, abstractmethod +from collections.abc import Mapping +from typing import Any + +from qiskit.quantum_info.operators.base_operator import BaseOperator +from qiskit.result import QuasiDistribution +from ..algorithm_result import AlgorithmResult +from ..list_or_dict import ListOrDict + + +class SamplingMinimumEigensolver(ABC): + """The Sampling Minimum Eigensolver Interface.""" + + @abstractmethod + def compute_minimum_eigenvalue( + self, + operator: BaseOperator, + aux_operators: ListOrDict[BaseOperator] | None = None, + ) -> "SamplingMinimumEigensolverResult": + """Compute the minimum eigenvalue of a diagonal operator. + + Args: + operator: Diagonal qubit operator. + aux_operators: Optional list of auxiliary operators to be evaluated with the + final state. + + Returns: + A :class:`~.SamplingMinimumEigensolverResult` containing the optimization result. + """ + pass + + @classmethod + def supports_aux_operators(cls) -> bool: + """Whether computing the expectation value of auxiliary operators is supported. + + If the minimum eigensolver computes an eigenstate 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 SamplingMinimumEigensolverResult(AlgorithmResult): + """Sampling Minimum Eigensolver Result. + + In contrast to the result of a :class:`~.MinimumEigenSolver`, this result also contains + the best measurement of the overall optimization and the samples of the final state. + """ + + def __init__(self) -> None: + super().__init__() + self._eigenvalue: complex | None = None + self._eigenstate: QuasiDistribution | None = None + self._aux_operator_values: ListOrDict[tuple[complex, dict[str, Any]]] | None = None + self._best_measurement: Mapping[str, Any] | None = None + + @property + def eigenvalue(self) -> complex | None: + """Return the approximation to the eigenvalue.""" + return self._eigenvalue + + @eigenvalue.setter + def eigenvalue(self, value: complex | None) -> None: + """Set the approximation to the eigenvalue.""" + self._eigenvalue = value + + @property + def eigenstate(self) -> QuasiDistribution | None: + """Return the quasi-distribution sampled from the final state. + + The ansatz is sampled when parameterized with the optimal parameters that where obtained + computing the minimum eigenvalue. The keys represent a measured classical value and the + value is a float for the quasi-probability of that result. + """ + return self._eigenstate + + @eigenstate.setter + def eigenstate(self, value: QuasiDistribution | None) -> None: + """Set the quasi-distribution sampled from the final state.""" + self._eigenstate = value + + @property + def aux_operators_evaluated(self) -> ListOrDict[tuple[complex, dict[str, Any]]] | None: + """Return aux operator expectation values and metadata. + + These are formatted as (mean, metadata). + """ + return self._aux_operator_values + + @aux_operators_evaluated.setter + def aux_operators_evaluated( + self, value: ListOrDict[tuple[complex, dict[str, Any]]] | None + ) -> None: + self._aux_operator_values = value + + @property + def best_measurement(self) -> Mapping[str, Any] | None: + """Return the best measurement over the entire optimization. + + Possesses keys: ``state``, ``bitstring``, ``value``, ``probability``. + """ + return self._best_measurement + + @best_measurement.setter + def best_measurement(self, value: Mapping[str, Any]) -> None: + """Set the best measurement over the entire optimization.""" + self._best_measurement = value diff --git a/qiskit_optimization/minimum_eigensolvers/sampling_vqe.py b/qiskit_optimization/minimum_eigensolvers/sampling_vqe.py new file mode 100644 index 000000000..fec4b5134 --- /dev/null +++ b/qiskit_optimization/minimum_eigensolvers/sampling_vqe.py @@ -0,0 +1,403 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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, optimized for diagonal Hamiltonians.""" + +from __future__ import annotations + +import logging +from collections.abc import Callable +from time import time +from typing import Any + +import numpy as np + +from qiskit.circuit import QuantumCircuit +from qiskit.passmanager import BasePassManager +from qiskit.primitives import BaseSamplerV1, BaseSamplerV2 +from qiskit.primitives.utils import init_observable +from qiskit.quantum_info.operators.base_operator import BaseOperator +from qiskit.result import QuasiDistribution + +from ..variational_algorithm import VariationalAlgorithm, VariationalResult +from ..exceptions import AlgorithmError +from ..list_or_dict import ListOrDict +from ..minimum_eigensolvers.sampling_mes import ( + SamplingMinimumEigensolver, + SamplingMinimumEigensolverResult, +) + +from ..observables_evaluator import estimate_observables +from ..optimizers.optimizer import Minimizer, Optimizer, OptimizerResult +from ..utils import validate_bounds, validate_initial_point + +# private function as we expect this to be updated in the next released +from ..utils.set_batching import _set_default_batchsize +from .diagonal_estimator import _DiagonalEstimator + +logger = logging.getLogger(__name__) + + +class SamplingVQE(VariationalAlgorithm, SamplingMinimumEigensolver): + r"""The Variational Quantum Eigensolver algorithm, optimized for diagonal Hamiltonians. + VQE is a hybrid quantum-classical algorithm that uses a variational technique to find the + minimum eigenvalue of a given diagonal Hamiltonian operator :math:`H_{\text{diag}}`. + In contrast to the :class:`~qiskit_algorithms.minimum_eigensolvers.VQE` class, the + ``SamplingVQE`` algorithm is executed using a :attr:`sampler` primitive. + An instance of ``SamplingVQE`` 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` to + minimize the objective function, which depends on the chosen :attr:`aggregation`. + 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 ``SamplingVQE`` object has been constructed. + Attributes: + sampler (BaseSamplerV1 or BaseSamplerV2): The sampler primitive to sample the circuits. + 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 an :class:`.Optimizer` or a callable implementing the + :class:`.Minimizer` protocol. + aggregation (float | Callable[[list[tuple[float, complex]], float] | None): + A float or callable to specify how the objective function evaluated on the basis states + should be aggregated. If a float, this specifies the :math:`\alpha \in [0,1]` parameter + for a CVaR expectation value [1]. If a callable, it takes a list of basis state + measurements specified as ``[(probability, objective_value)]`` and return an objective + value as float. If None, all an ordinary expectation value is calculated. + 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 value, and the + metadata dictionary. + References: + [1]: Barkoutsos, P. K., Nannicini, G., Robert, A., Tavernelli, I., and Woerner, S., + "Improving Variational Quantum Optimization using CVaR" + `arXiv:1907.04769 `_ + """ + + def __init__( + self, + sampler: BaseSamplerV1 | BaseSamplerV2, + ansatz: QuantumCircuit, + optimizer: Optimizer | Minimizer, + *, + initial_point: np.ndarray | None = None, + aggregation: float | Callable[[list[float]], float] | None = None, + callback: Callable[[int, np.ndarray, float, dict[str, Any]], None] | None = None, + passmanager: BasePassManager | None = None, + ) -> None: + r""" + Args: + sampler: The sampler primitive to sample the circuits. + ansatz: A parameterized quantum circuit to prepare the trial state. + optimizer: A classical optimizer to find the minimum energy. This can either be an + :class:`.Optimizer` or a callable implementing the :class:`.Minimizer` protocol. + 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. ``SamplingVQE`` 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. + aggregation: A float or callable to specify how the objective function evaluated on the + basis states should be aggregated. + 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.sampler = sampler + self.ansatz = ansatz + self.optimizer = optimizer + self.aggregation = aggregation + self.callback = callback + self.passmanager = passmanager + + # this has to go via getters and setters due to the VariationalAlgorithm interface + self._initial_point = initial_point + + @property + def initial_point(self) -> np.ndarray | None: + """Return the initial point.""" + return self._initial_point + + @initial_point.setter + def initial_point(self, value: np.ndarray | None) -> None: + """Set the initial point.""" + self._initial_point = value + + def _check_operator_ansatz(self, operator: BaseOperator): + """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.") + + @classmethod + def supports_aux_operators(cls) -> bool: + return True + + def compute_minimum_eigenvalue( + self, + operator: BaseOperator, + aux_operators: ListOrDict[BaseOperator] | None = None, + ) -> SamplingMinimumEigensolverResult: + # check that the number of qubits of operator and ansatz match, and resize if possible + self._check_operator_ansatz(operator) + + if len(self.ansatz.clbits) > 0: + self.ansatz.remove_final_measurements() + self.ansatz.measure_all() + + initial_point = validate_initial_point(self.initial_point, self.ansatz) + + bounds = validate_bounds(self.ansatz) + + if self.passmanager: + ansatz: QuantumCircuit = self.passmanager.run(self.ansatz) + layout = ansatz.layout + operator = init_observable(operator) + operator = operator.apply_layout(layout) + if aux_operators: + if isinstance(aux_operators, list): + aux_operators = [op.apply_layout(layout) for op in aux_operators] + else: + aux_operators = { + key: op.apply_layout(layout) for key, op in aux_operators.items() + } + else: + ansatz = self.ansatz + + # NOTE: we type ignore below because the `return_best_measurement=True` is guaranteed to + # return a tuple + evaluate_energy, best_measurement = self._get_evaluate_energy( # type: ignore[misc] + operator, ansatz, return_best_measurement=True + ) + + start_time = time() + + if callable(self.optimizer): + optimizer_result = self.optimizer( + fun=evaluate_energy, # type: ignore[arg-type] + x0=initial_point, + jac=None, + bounds=bounds, + ) + else: + # we always want to submit as many estimations per job as possible for minimal + # overhead on the hardware + was_updated = _set_default_batchsize(self.optimizer) + + optimizer_result = self.optimizer.minimize( + fun=evaluate_energy, # type: ignore[arg-type] + x0=initial_point, + bounds=bounds, + ) + + # reset to original value + if was_updated: + self.optimizer.set_max_evals_grouped(None) + + optimizer_time = time() - start_time + + logger.info( + "Optimization complete in %s seconds.\nFound opt_params %s.", + optimizer_time, + optimizer_result.x, + ) + + if isinstance(self.sampler, BaseSamplerV1): + final_state = self.sampler.run([ansatz], [optimizer_result.x]).result().quasi_dists[0] + else: + result = self.sampler.run([(ansatz, optimizer_result.x)]).result()[0] + creg = ansatz.cregs[0].name + counts = getattr(result.data, creg).get_counts() + shots = sum(counts.values()) + final_state = QuasiDistribution( + {key: val / shots for key, val in counts.items()}, shots=shots + ) + + if aux_operators is not None: + aux_operators_evaluated = estimate_observables( + _DiagonalEstimator(sampler=self.sampler), + ansatz, + aux_operators, + optimizer_result.x, # type: ignore[arg-type] + ) + else: + aux_operators_evaluated = None + + return self._build_sampling_vqe_result( + self.ansatz.copy(), + optimizer_result, + aux_operators_evaluated, # type: ignore[arg-type] + best_measurement, + final_state, + optimizer_time, + ) + + def _get_evaluate_energy( + self, + operator: BaseOperator, + ansatz: QuantumCircuit, + return_best_measurement: bool = False, + ) -> ( + Callable[[np.ndarray], np.ndarray | float] + | tuple[Callable[[np.ndarray], np.ndarray | float], dict[str, Any]] + ): + """Returns a function handle to evaluate the energy at given parameters. + This is the objective function to be passed to the optimizer that is used for evaluation. + Args: + operator: The operator whose energy to evaluate. + ansatz: The ansatz preparing the quantum state. + return_best_measurement: If True, a handle to a dictionary containing the best + measurement evaluated with the cost function. + Returns: + A tuple of a callable evaluating the energy and (optionally) a dictionary containing the + best measurement of the energy evaluation. + Raises: + AlgorithmError: If the circuit is not parameterized (i.e. has 0 free parameters). + """ + num_parameters = ansatz.num_parameters + if num_parameters == 0: + raise AlgorithmError("The ansatz must be parameterized, but has 0 free parameters.") + + # avoid creating an instance variable to remain stateless regarding results + eval_count = 0 + + best_measurement = {"best": None} + + def store_best_measurement(best): + for best_i in best: + if best_measurement["best"] is None or _compare_measurements( + best_i, best_measurement["best"] + ): + best_measurement["best"] = best_i + + estimator = _DiagonalEstimator( + sampler=self.sampler, + callback=store_best_measurement, + aggregation=self.aggregation, # type: ignore[arg-type] + ) + + 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) + + estimator_result = estimator.run( + batch_size * [ansatz], batch_size * [operator], parameters + ).result() + 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) + + result = values if len(values) > 1 else values[0] + return np.real(result) + + if return_best_measurement: + return evaluate_energy, best_measurement + + return evaluate_energy + + def _build_sampling_vqe_result( # pylint: disable=too-many-positional-arguments + self, + ansatz: QuantumCircuit, + optimizer_result: OptimizerResult, + aux_operators_evaluated: ListOrDict[tuple[complex, tuple[complex, int]]], + best_measurement: dict[str, Any], + final_state: QuasiDistribution, + optimizer_time: float, + ) -> SamplingVQEResult: + result = SamplingVQEResult() + result.eigenvalue = optimizer_result.fun + result.cost_function_evals = optimizer_result.nfev + result.optimal_point = optimizer_result.x # type: ignore[assignment] + result.optimal_parameters = dict( + zip(self.ansatz.parameters, optimizer_result.x) # type: ignore[arg-type] + ) + result.optimal_value = optimizer_result.fun + result.optimizer_time = optimizer_time + result.aux_operators_evaluated = aux_operators_evaluated # type: ignore[assignment] + result.optimizer_result = optimizer_result + result.best_measurement = best_measurement["best"] + result.eigenstate = final_state + result.optimal_circuit = ansatz + return result + + +class SamplingVQEResult(VariationalResult, SamplingMinimumEigensolverResult): + """The SamplingVQE Result.""" + + def __init__(self) -> None: + super().__init__() + self._cost_function_evals: int | None = None + + @property + def cost_function_evals(self) -> int | None: + """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 + + +def _compare_measurements(candidate, current_best): + """Compare two best measurements. Returns True if the candidate is better than current value. + + This compares the following two criteria, in this precedence: + + 1. The smaller objective value is better + 2. The higher probability for the objective value is better + + """ + if candidate["value"] < current_best["value"]: + return True + elif candidate["value"] == current_best["value"]: + return candidate["probability"] > current_best["probability"] + return False diff --git a/qiskit_optimization/observables_evaluator.py b/qiskit_optimization/observables_evaluator.py new file mode 100644 index 000000000..2e705c468 --- /dev/null +++ b/qiskit_optimization/observables_evaluator.py @@ -0,0 +1,129 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2021, 2024. +# +# 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. + +"""Evaluator of observables for algorithms.""" + +from __future__ import annotations +from collections.abc import Sequence +from typing import Any + +import numpy as np + +from qiskit import QuantumCircuit +from qiskit.quantum_info import SparsePauliOp +from qiskit.primitives import BaseEstimator +from qiskit.quantum_info.operators.base_operator import BaseOperator + +from .exceptions import AlgorithmError +from .list_or_dict import ListOrDict + + +def estimate_observables( + estimator: BaseEstimator, + quantum_state: QuantumCircuit, + observables: ListOrDict[BaseOperator], + parameter_values: Sequence[float] | None = None, + threshold: float = 1e-12, +) -> ListOrDict[tuple[float, dict[str, Any]]]: + """ + Accepts a sequence of operators and calculates their expectation values - means + 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: 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, metadata). + + Raises: + AlgorithmError: If a primitive job is not successful. + """ + + if isinstance(observables, dict): + observables_list = list(observables.values()) + else: + observables_list = observables + + if len(observables_list) > 0: + observables_list = _handle_zero_ops(observables_list) + quantum_state = [quantum_state] * len(observables) + parameter_values_: Sequence[float] | Sequence[Sequence[float]] | None = parameter_values + if parameter_values is not None: + parameter_values_ = [parameter_values] * len(observables) + try: + 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 + + metadata = estimator_job.result().metadata + # Discard values below threshold + observables_means = expectation_values * (np.abs(expectation_values) > threshold) + # zip means and metadata into tuples + observables_results = list(zip(observables_means, metadata)) + else: + observables_results = [] + + return _prepare_result(observables_results, observables) + + +def _handle_zero_ops( + observables_list: list[BaseOperator], +) -> list[BaseOperator]: + """Replaces all occurrence of operators equal to 0 in the list with an equivalent ``SparsePauliOp`` + operator.""" + if observables_list: + zero_op = SparsePauliOp.from_list([("I" * observables_list[0].num_qubits, 0)]) + for ind, observable in enumerate(observables_list): + if observable == 0: + observables_list[ind] = zero_op + return observables_list + + +def _prepare_result( + observables_results: list[tuple[float, dict]], + observables: ListOrDict[BaseOperator], +) -> ListOrDict[tuple[float, dict[str, Any]]]: + """ + Prepares a list of tuples of eigenvalues and metadata tuples from + ``observables_results`` and ``observables``. + + Args: + 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, metadata). + """ + + observables_eigenvalues: ListOrDict[tuple[float, dict]] + + if isinstance(observables, list): + observables_eigenvalues = [] + for value in observables_results: + observables_eigenvalues.append(value) + + else: + observables_eigenvalues = {} + for key, value in zip(observables.keys(), observables_results): + observables_eigenvalues[key] = value + + return observables_eigenvalues diff --git a/qiskit_optimization/optimizers/__init__.py b/qiskit_optimization/optimizers/__init__.py new file mode 100644 index 000000000..e64c6e57f --- /dev/null +++ b/qiskit_optimization/optimizers/__init__.py @@ -0,0 +1,134 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2018, 2024. +# +# 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. + +""" +Optimizers (:mod:`qiskit_algorithms.optimizers`) +================================================ +Classical Optimizers. + +This package contains a variety of classical optimizers and were designed for use by +qiskit_algorithm's quantum variational algorithms, such as :class:`~qiskit_algorithms.VQE`. +Logically, these optimizers can be divided into two categories: + +`Local Optimizers`_ + Given an optimization problem, a **local optimizer** is a function + that attempts to find an optimal value within the neighboring set of a candidate solution. + +`Global Optimizers`_ + Given an optimization problem, a **global optimizer** is a function + that attempts to find an optimal value among all possible solutions. + +.. currentmodule:: qiskit_algorithms.optimizers + +Optimizer Base Classes +---------------------- + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + + OptimizerResult + Optimizer + Minimizer + +Steppable Optimization +---------------------- + +.. autosummary:: + :toctree: ../stubs/ + + optimizer_utils + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + + SteppableOptimizer + AskData + TellData + OptimizerState + + +Local Optimizers +---------------- + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + + ADAM + AQGD + CG + COBYLA + L_BFGS_B + GSLS + GradientDescent + GradientDescentState + NELDER_MEAD + NFT + P_BFGS + POWELL + SLSQP + SPSA + QNSPSA + TNC + SciPyOptimizer + UMDA + +Qiskit also provides the following optimizers, which are built-out using the optimizers from +`scikit-quant `_. The ``scikit-quant`` package +is not installed by default but must be explicitly installed, if desired, by the user. The +optimizers therein are provided under various licenses, hence it has been made an optional install. +To install the ``scikit-quant`` dependent package you can use ``pip install scikit-quant``. + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + + BOBYQA + IMFIL + SNOBFIT + +Global Optimizers +----------------- +The global optimizers here all use `NLOpt `_ for their +core function and can only be used if the optional dependent ``NLOpt`` package is installed. +To install the ``NLOpt`` dependent package you can use ``pip install nlopt``. + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + + CRS + DIRECT_L + DIRECT_L_RAND + ESCH + ISRES + +""" + +from .optimizer import Minimizer, Optimizer, OptimizerResult, OptimizerSupportLevel +from .spsa import SPSA +from .cobyla import COBYLA +from .nelder_mead import NELDER_MEAD +from .scipy_optimizer import SciPyOptimizer + +__all__ = [ + "Optimizer", + "OptimizerSupportLevel", + "OptimizerResult", + "Minimizer", + "SPSA", + "COBYLA", + "NELDER_MEAD", + "SciPyOptimizer", +] diff --git a/qiskit_optimization/optimizers/cobyla.py b/qiskit_optimization/optimizers/cobyla.py new file mode 100644 index 000000000..fe4e4a570 --- /dev/null +++ b/qiskit_optimization/optimizers/cobyla.py @@ -0,0 +1,59 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2018, 2024. +# +# 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. + +"""Constrained Optimization By Linear Approximation optimizer.""" + +from __future__ import annotations + +from .scipy_optimizer import SciPyOptimizer + + +class COBYLA(SciPyOptimizer): + """ + Constrained Optimization By Linear Approximation optimizer. + + COBYLA is a numerical optimization method for constrained problems + where the derivative of the objective function is not known. + + Uses scipy.optimize.minimize COBYLA. + For further detail, please refer to + https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html + """ + + _OPTIONS = ["maxiter", "disp", "rhobeg"] + + # pylint: disable=unused-argument + def __init__( # pylint: disable=too-many-positional-arguments + self, + maxiter: int = 1000, + disp: bool = False, + rhobeg: float = 1.0, + tol: float | None = None, + options: dict | None = None, + **kwargs, + ) -> None: + """ + Args: + maxiter: Maximum number of function evaluations. + disp: Set to True to print convergence messages. + rhobeg: Reasonable initial changes to the variables. + tol: Final accuracy in the optimization (not precisely guaranteed). + This is a lower bound on the size of the trust region. + options: A dictionary of solver options. + kwargs: additional kwargs for scipy.optimize.minimize. + """ + if options is None: + options = {} + for k, v in list(locals().items()): + if k in self._OPTIONS: + options[k] = v + super().__init__(method="COBYLA", options=options, tol=tol, **kwargs) diff --git a/qiskit_optimization/optimizers/nelder_mead.py b/qiskit_optimization/optimizers/nelder_mead.py new file mode 100644 index 000000000..12b3bb191 --- /dev/null +++ b/qiskit_optimization/optimizers/nelder_mead.py @@ -0,0 +1,73 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2018, 2024. +# +# 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. + +"""Nelder-Mead optimizer.""" +from __future__ import annotations + + +from .scipy_optimizer import SciPyOptimizer + + +class NELDER_MEAD(SciPyOptimizer): # pylint: disable=invalid-name + """ + Nelder-Mead optimizer. + + The Nelder-Mead algorithm performs unconstrained optimization; it ignores bounds + or constraints. It is used to find the minimum or maximum of an objective function + in a multidimensional space. It is based on the Simplex algorithm. Nelder-Mead + is robust in many applications, especially when the first and second derivatives of the + objective function are not known. + + However, if the numerical computation of the derivatives can be trusted to be accurate, + other algorithms using the first and/or second derivatives information might be preferred to + Nelder-Mead for their better performance in the general case, especially in consideration of + the fact that the Nelder–Mead technique is a heuristic search method that can converge to + non-stationary points. + + Uses scipy.optimize.minimize Nelder-Mead. + For further detail, please refer to + See https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html + """ + + _OPTIONS = ["maxiter", "maxfev", "disp", "xatol", "adaptive"] + + # pylint: disable=unused-argument + def __init__( # pylint: disable=too-many-positional-arguments + self, + maxiter: int | None = None, + maxfev: int = 1000, + disp: bool = False, + xatol: float = 0.0001, + tol: float | None = None, + adaptive: bool = False, + options: dict | None = None, + **kwargs, + ) -> None: + """ + Args: + maxiter: Maximum allowed number of iterations. If both maxiter and maxfev are set, + minimization will stop at the first reached. + maxfev: Maximum allowed number of function evaluations. If both maxiter and + maxfev are set, minimization will stop at the first reached. + disp: Set to True to print convergence messages. + xatol: Absolute error in xopt between iterations that is acceptable for convergence. + tol: Tolerance for termination. + adaptive: Adapt algorithm parameters to dimensionality of problem. + options: A dictionary of solver options. + kwargs: additional kwargs for scipy.optimize.minimize. + """ + if options is None: + options = {} + for k, v in list(locals().items()): + if k in self._OPTIONS: + options[k] = v + super().__init__(method="Nelder-Mead", options=options, tol=tol, **kwargs) diff --git a/qiskit_optimization/optimizers/optimizer.py b/qiskit_optimization/optimizers/optimizer.py new file mode 100644 index 000000000..b9effe9b8 --- /dev/null +++ b/qiskit_optimization/optimizers/optimizer.py @@ -0,0 +1,389 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2018, 2024. +# +# 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. + +"""Optimizer interface""" + +from __future__ import annotations + +from abc import ABC, abstractmethod +from collections.abc import Callable +from enum import IntEnum +import logging +from typing import Any, Union, Protocol + +import numpy as np +import scipy + +from ..algorithm_result import AlgorithmResult + +logger = logging.getLogger(__name__) + +POINT = Union[float, np.ndarray] # pylint: disable=invalid-name + + +class OptimizerResult(AlgorithmResult): + """The result of an optimization routine.""" + + def __init__(self) -> None: + super().__init__() + self._x: POINT | None = None # pylint: disable=invalid-name + self._fun: float | None = None + self._jac: POINT | None = None + self._nfev: int | None = None + self._njev: int | None = None + self._nit: int | None = None + + @property + def x(self) -> POINT | None: + """The final point of the minimization.""" + return self._x + + @x.setter + def x(self, x: POINT | None) -> None: + """Set the final point of the minimization.""" + self._x = x + + @property + def fun(self) -> float | None: + """The final value of the minimization.""" + return self._fun + + @fun.setter + def fun(self, fun: float | None) -> None: + """Set the final value of the minimization.""" + self._fun = fun + + @property + def jac(self) -> POINT | None: + """The final gradient of the minimization.""" + return self._jac + + @jac.setter + def jac(self, jac: POINT | None) -> None: + """Set the final gradient of the minimization.""" + self._jac = jac + + @property + def nfev(self) -> int | None: + """The total number of function evaluations.""" + return self._nfev + + @nfev.setter + def nfev(self, nfev: int | None) -> None: + """Set the total number of function evaluations.""" + self._nfev = nfev + + @property + def njev(self) -> int | None: + """The total number of gradient evaluations.""" + return self._njev + + @njev.setter + def njev(self, njev: int | None) -> None: + """Set the total number of gradient evaluations.""" + self._njev = njev + + @property + def nit(self) -> int | None: + """The total number of iterations.""" + return self._nit + + @nit.setter + def nit(self, nit: int | None) -> None: + """Set the total number of iterations.""" + self._nit = nit + + +class Minimizer(Protocol): + """Callable Protocol for minimizer. + + This interface is based on `SciPy's optimize module + `__. + + This protocol defines a callable taking the following parameters: + + fun + The objective function to minimize (for example the energy in the case of the VQE). + x0 + The initial point for the optimization. + jac + The gradient of the objective function. + bounds + Parameters bounds for the optimization. Note that these might not be supported + by all optimizers. + + and which returns a minimization result object (either SciPy's or Qiskit's). + """ + + # pylint: disable=invalid-name + def __call__( + self, + fun: Callable[[np.ndarray], float], + x0: np.ndarray, + jac: Callable[[np.ndarray], np.ndarray] | None, + bounds: list[tuple[float, float]] | None, + ) -> scipy.optimize.OptimizeResult | OptimizerResult: + """Minimize the objective function. + + This interface is based on `SciPy's optimize module `__. + + Args: + fun: The objective function to minimize (for example the energy in the case of the VQE). + x0: The initial point for the optimization. + jac: The gradient of the objective function. + bounds: Parameters bounds for the optimization. Note that these might not be supported + by all optimizers. + + Returns: + The minimization result object (either SciPy's or Qiskit's). + """ + ... # pylint: disable=unnecessary-ellipsis + + +class OptimizerSupportLevel(IntEnum): + """Support Level enum for features such as bounds, gradient and initial point""" + + # pylint: disable=invalid-name + not_supported = 0 # Does not support the corresponding parameter in optimize() + ignored = 1 # Feature can be passed as non None but will be ignored + supported = 2 # Feature is supported + required = 3 # Feature is required and must be given, None is invalid + + +class Optimizer(ABC): + """Base class for optimization algorithm.""" + + @abstractmethod + def __init__(self): + """ + Initialize the optimization algorithm, setting the support + level for _gradient_support_level, _bound_support_level, + _initial_point_support_level, and empty options. + """ + self._gradient_support_level = self.get_support_level()["gradient"] + self._bounds_support_level = self.get_support_level()["bounds"] + self._initial_point_support_level = self.get_support_level()["initial_point"] + self._options = {} + self._max_evals_grouped = None + + @abstractmethod + def get_support_level(self): + """Return support level dictionary""" + raise NotImplementedError + + def set_options(self, **kwargs): + """ + Sets or updates values in the options dictionary. + + The options dictionary may be used internally by a given optimizer to + pass additional optional values for the underlying optimizer/optimization + function used. The options dictionary may be initially populated with + a set of key/values when the given optimizer is constructed. + + Args: + kwargs (dict): options, given as name=value. + """ + for name, value in kwargs.items(): + self._options[name] = value + logger.debug("options: %s", self._options) + + # pylint: disable=invalid-name + @staticmethod + def gradient_num_diff(x_center, f, epsilon, max_evals_grouped=None): + """ + We compute the gradient with the numeric differentiation in the parallel way, + around the point x_center. + + Args: + x_center (ndarray): point around which we compute the gradient + f (func): the function of which the gradient is to be computed. + epsilon (float): the epsilon used in the numeric differentiation. + max_evals_grouped (int): max evals grouped, defaults to 1 (i.e. no batching). + Returns: + grad: the gradient computed + + """ + if max_evals_grouped is None: # no batching by default + max_evals_grouped = 1 + + forig = f(*((x_center,))) + grad = [] + ei = np.zeros((len(x_center),), float) + todos = [] + for k in range(len(x_center)): + ei[k] = 1.0 + d = epsilon * ei + todos.append(x_center + d) + ei[k] = 0.0 + + counter = 0 + chunk = [] + chunks = [] + length = len(todos) + # split all points to chunks, where each chunk has batch_size points + for i in range(length): + x = todos[i] + chunk.append(x) + counter += 1 + # the last one does not have to reach batch_size + if counter == max_evals_grouped or i == length - 1: + chunks.append(chunk) + chunk = [] + counter = 0 + + for chunk in chunks: # eval the chunks in order + parallel_parameters = np.concatenate(chunk) + todos_results = f(parallel_parameters) # eval the points in a chunk (order preserved) + if isinstance(todos_results, float): + grad.append((todos_results - forig) / epsilon) + else: + for todor in todos_results: + grad.append((todor - forig) / epsilon) + + return np.array(grad) + + @staticmethod + def wrap_function(function, args): + """ + Wrap the function to implicitly inject the args at the call of the function. + + Args: + function (func): the target function + args (tuple): the args to be injected + Returns: + function_wrapper: wrapper + """ + + def function_wrapper(*wrapper_args): + return function(*(wrapper_args + args)) + + return function_wrapper + + @property + def setting(self): + """Return setting""" + ret = f"Optimizer: {self.__class__.__name__}\n" + params = "" + for key, value in self.__dict__.items(): + if key[0] == "_": + params += f"-- {key[1:]}: {value}\n" + ret += f"{params}" + return ret + + @property + def settings(self) -> dict[str, Any]: + """The optimizer settings in a dictionary format. + + The settings can for instance be used for JSON-serialization (if all settings are + serializable, which e.g. doesn't hold per default for callables), such that the + optimizer object can be reconstructed as + + .. code-block:: + + settings = optimizer.settings + # JSON serialize and send to another server + optimizer = OptimizerClass(**settings) + + """ + raise NotImplementedError("The settings method is not implemented per default.") + + @abstractmethod + def minimize( + self, + fun: Callable[[POINT], float], + x0: POINT, + jac: Callable[[POINT], POINT] | None = None, + bounds: list[tuple[float, float]] | None = None, + ) -> OptimizerResult: + """Minimize the scalar function. + + Args: + fun: The scalar function to minimize. + x0: The initial point for the minimization. + jac: The gradient of the scalar function ``fun``. + bounds: Bounds for the variables of ``fun``. This argument might be ignored if the + optimizer does not support bounds. + + Returns: + The result of the optimization, containing e.g. the result as attribute ``x``. + """ + raise NotImplementedError() + + @property + def gradient_support_level(self): + """Returns gradient support level""" + return self._gradient_support_level + + @property + def is_gradient_ignored(self): + """Returns is gradient ignored""" + return self._gradient_support_level == OptimizerSupportLevel.ignored + + @property + def is_gradient_supported(self): + """Returns is gradient supported""" + return self._gradient_support_level != OptimizerSupportLevel.not_supported + + @property + def is_gradient_required(self): + """Returns is gradient required""" + return self._gradient_support_level == OptimizerSupportLevel.required + + @property + def bounds_support_level(self): + """Returns bounds support level""" + return self._bounds_support_level + + @property + def is_bounds_ignored(self): + """Returns is bounds ignored""" + return self._bounds_support_level == OptimizerSupportLevel.ignored + + @property + def is_bounds_supported(self): + """Returns is bounds supported""" + return self._bounds_support_level != OptimizerSupportLevel.not_supported + + @property + def is_bounds_required(self): + """Returns is bounds required""" + return self._bounds_support_level == OptimizerSupportLevel.required + + @property + def initial_point_support_level(self): + """Returns initial point support level""" + return self._initial_point_support_level + + @property + def is_initial_point_ignored(self): + """Returns is initial point ignored""" + return self._initial_point_support_level == OptimizerSupportLevel.ignored + + @property + def is_initial_point_supported(self): + """Returns is initial point supported""" + return self._initial_point_support_level != OptimizerSupportLevel.not_supported + + @property + def is_initial_point_required(self): + """Returns is initial point required""" + return self._initial_point_support_level == OptimizerSupportLevel.required + + def print_options(self): + """Print algorithm-specific options.""" + for name in sorted(self._options): + logger.debug("%s = %s", name, str(self._options[name])) + + def set_max_evals_grouped(self, limit): + """Set max evals grouped""" + self._max_evals_grouped = limit diff --git a/qiskit_optimization/optimizers/scipy_optimizer.py b/qiskit_optimization/optimizers/scipy_optimizer.py new file mode 100644 index 000000000..d1f143d3f --- /dev/null +++ b/qiskit_optimization/optimizers/scipy_optimizer.py @@ -0,0 +1,191 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2018, 2024. +# +# 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. + +"""Wrapper class of scipy.optimize.minimize.""" +from __future__ import annotations + +from collections.abc import Callable +from typing import Any + +import numpy as np +from scipy.optimize import minimize + +from ..utils.validation import validate_min +from .optimizer import Optimizer, OptimizerSupportLevel, OptimizerResult, POINT + + +class SciPyOptimizer(Optimizer): + """A general Qiskit Optimizer wrapping scipy.optimize.minimize. + + For further detail, please refer to + https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html + """ + + _bounds_support_methods = {"l-bfgs-b", "tnc", "slsqp", "powell", "trust-constr"} + _gradient_support_methods = { + "cg", + "bfgs", + "newton-cg", + "l-bfgs-b", + "tnc", + "slsqp", + "dogleg", + "trust-ncg", + "trust-krylov", + "trust-exact", + "trust-constr", + } + + def __init__( + self, + method: str | Callable, + options: dict[str, Any] | None = None, + max_evals_grouped: int = 1, + **kwargs, + ): + """ + Args: + method: Type of solver. + options: A dictionary of solver options. + kwargs: additional kwargs for scipy.optimize.minimize. + max_evals_grouped: Max number of default gradient evaluations performed simultaneously. + """ + self._method = method.lower() if isinstance(method, str) else method + # Set support level + if self._method in self._bounds_support_methods: + self._bounds_support_level = OptimizerSupportLevel.supported + else: + self._bounds_support_level = OptimizerSupportLevel.ignored + if self._method in self._gradient_support_methods: + self._gradient_support_level = OptimizerSupportLevel.supported + else: + self._gradient_support_level = OptimizerSupportLevel.ignored + self._initial_point_support_level = OptimizerSupportLevel.required + + self._options = options if options is not None else {} + validate_min("max_evals_grouped", max_evals_grouped, 1) + self._max_evals_grouped = max_evals_grouped + self._kwargs = kwargs + + if "bounds" in self._kwargs: + raise RuntimeError( + "Optimizer bounds should be passed to SciPyOptimizer.minimize() and is not " + "supported in SciPyOptimizer constructor kwargs." + ) + if "bounds" in self._options: + raise RuntimeError( + "Optimizer bounds should be passed to SciPyOptimizer.minimize() and not as " + "options." + ) + + def get_support_level(self): + """Return support level dictionary""" + return { + "gradient": self._gradient_support_level, + "bounds": self._bounds_support_level, + "initial_point": self._initial_point_support_level, + } + + @property + def settings(self) -> dict[str, Any]: + options = self._options.copy() + if hasattr(self, "_OPTIONS"): + # all _OPTIONS should be keys in self._options, but add a failsafe here + attributes = [ + option + for option in self._OPTIONS # pylint: disable=no-member + if option in options.keys() + ] + + settings = {attr: options.pop(attr) for attr in attributes} + else: + settings = {} + + settings["max_evals_grouped"] = self._max_evals_grouped + settings["options"] = options + settings.update(self._kwargs) + + # the subclasses don't need the "method" key as the class type specifies the method + if self.__class__ == SciPyOptimizer: + settings["method"] = self._method + + return settings + + def minimize( + self, + fun: Callable[[POINT], float], + x0: POINT, + jac: Callable[[POINT], POINT] | None = None, + bounds: list[tuple[float, float]] | None = None, + ) -> OptimizerResult: + + # Remove ignored bounds to suppress the warning of scipy.optimize.minimize + if self.is_bounds_ignored: + bounds = None + + # Remove ignored gradient to suppress the warning of scipy.optimize.minimize + if self.is_gradient_ignored: + jac = None + + if self.is_gradient_supported and jac is None and self._max_evals_grouped > 1: + if "eps" in self._options: + epsilon = self._options["eps"] + else: + epsilon = ( + 1e-8 if self._method in {"l-bfgs-b", "tnc"} else np.sqrt(np.finfo(float).eps) + ) + jac = Optimizer.wrap_function( + Optimizer.gradient_num_diff, (fun, epsilon, self._max_evals_grouped) + ) + + # Workaround for L_BFGS_B because it does not accept np.ndarray. + # See https://github.com/Qiskit/qiskit/pull/6373. + if jac is not None and self._method == "l-bfgs-b": + jac = self._wrap_gradient(jac) + + # Starting in scipy 1.9.0 maxiter is deprecated and maxfun (added in 1.5.0) + # should be used instead + swapped_deprecated_args = False + if self._method == "tnc" and "maxiter" in self._options: + swapped_deprecated_args = True + self._options["maxfun"] = self._options.pop("maxiter") + + raw_result = minimize( + fun=fun, + x0=x0, + method=self._method, + jac=jac, + bounds=bounds, + options=self._options, + **self._kwargs, + ) + if swapped_deprecated_args: + self._options["maxiter"] = self._options.pop("maxfun") + + result = OptimizerResult() + result.x = raw_result.x + result.fun = raw_result.fun + result.nfev = raw_result.nfev + result.njev = raw_result.get("njev", None) + result.nit = raw_result.get("nit", None) + + return result + + @staticmethod + def _wrap_gradient(gradient_function): + def wrapped_gradient(x): + gradient = gradient_function(x) + if isinstance(gradient, np.ndarray): + return gradient.tolist() + return gradient + + return wrapped_gradient diff --git a/qiskit_optimization/optimizers/spsa.py b/qiskit_optimization/optimizers/spsa.py new file mode 100644 index 000000000..5e80f11e5 --- /dev/null +++ b/qiskit_optimization/optimizers/spsa.py @@ -0,0 +1,775 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2018, 2024. +# +# 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. + +"""Simultaneous Perturbation Stochastic Approximation (SPSA) optimizer. + +This implementation allows both standard first-order and second-order SPSA. +""" +from __future__ import annotations + +from collections import deque +from collections.abc import Iterator +from typing import Callable, Any, SupportsFloat +import logging +import warnings +from time import time + +import scipy +import numpy as np + +from ..utils import algorithm_globals + +from .optimizer import Optimizer, OptimizerSupportLevel, OptimizerResult, POINT + +# number of function evaluations, parameters, loss, stepsize, accepted +CALLBACK = Callable[[int, np.ndarray, float, SupportsFloat, bool], None] +TERMINATIONCHECKER = Callable[[int, np.ndarray, float, SupportsFloat, bool], bool] + +logger = logging.getLogger(__name__) + + +class SPSA(Optimizer): + """Simultaneous Perturbation Stochastic Approximation (SPSA) optimizer. + + SPSA [1] is an gradient descent method for optimizing systems with multiple unknown parameters. + As an optimization method, it is appropriately suited to large-scale population models, + adaptive modeling, and simulation optimization. + + .. see also:: + + Many examples are presented at the `SPSA Web site `__. + + The main feature of SPSA is the stochastic gradient approximation, which requires only two + measurements of the objective function, regardless of the dimension of the optimization + problem. + + Additionally, to standard first-order SPSA, where only gradient information is used, this + implementation also allows second-order SPSA (2-SPSA) [2]. In 2-SPSA we also estimate the + Hessian of the loss with a stochastic approximation and multiply the gradient with the + inverse Hessian to take local curvature into account and improve convergence. + Notably this Hessian estimate requires only a constant number of function evaluations + unlike an exact evaluation of the Hessian, which scales quadratically in the number of + function evaluations. + + .. note:: + + SPSA can be used in the presence of noise, and it is therefore indicated in situations + involving measurement uncertainty on a quantum computation when finding a minimum. + If you are executing a variational algorithm using a Quantum ASseMbly Language (QASM) + simulator or a real device, SPSA would be the most recommended choice among the optimizers + provided here. + + The optimization process can include a calibration phase if neither the ``learning_rate`` nor + ``perturbation`` is provided, which requires additional functional evaluations. + (Note that either both or none must be set.) For further details on the automatic calibration, + please refer to the supplementary information section IV. of [3]. + + .. note:: + + This component has some function that is normally random. If you want to reproduce behavior + then you should set the random number generator seed in the algorithm_globals + (``qiskit_algorithms.utils.algorithm_globals.random_seed = seed``). + + + Examples: + + This short example runs SPSA for the ground state calculation of the ``Z ^ Z`` + observable where the ansatz is a ``PauliTwoDesign`` circuit. + + .. code-block:: python + + import numpy as np + from qiskit_algorithms.optimizers import SPSA + from qiskit.circuit.library import PauliTwoDesign + from qiskit.primitives import Estimator + from qiskit.quantum_info import SparsePauliOp + + ansatz = PauliTwoDesign(2, reps=1, seed=2) + observable = SparsePauliOp("ZZ") + initial_point = np.random.random(ansatz.num_parameters) + estimator = Estimator() + + def loss(x): + job = estimator.run([ansatz], [observable], [x]) + return job.result().values[0] + + spsa = SPSA(maxiter=300) + result = spsa.minimize(loss, x0=initial_point) + + To use the Hessian information, i.e. 2-SPSA, you can add `second_order=True` to the + initializer of the `SPSA` class, the rest of the code remains the same. + + .. code-block:: python + + two_spsa = SPSA(maxiter=300, second_order=True) + result = two_spsa.minimize(loss, x0=initial_point) + + The `termination_checker` can be used to implement a custom termination criterion. + + .. code-block:: python + + import numpy as np + from qiskit_algorithms.optimizers import SPSA + + def objective(x): + return np.linalg.norm(x) + .04*np.random.rand(1) + + class TerminationChecker: + + def __init__(self, N : int): + self.N = N + self.values = [] + + def __call__(self, nfev, parameters, value, stepsize, accepted) -> bool: + self.values.append(value) + + if len(self.values) > self.N: + last_values = self.values[-self.N:] + pp = np.polyfit(range(self.N), last_values, 1) + slope = pp[0] / self.N + + if slope > 0: + return True + return False + + spsa = SPSA(maxiter=200, termination_checker=TerminationChecker(10)) + result = spsa.minimize(objective, x0=[0.5, 0.5]) + print(f'SPSA completed after {result.nit} iterations') + + References: + + [1]: J. C. Spall (1998). An Overview of the Simultaneous Perturbation Method for Efficient + Optimization, Johns Hopkins APL Technical Digest, 19(4), 482–492. + `Online at jhuapl.edu. `_ + + [2]: J. C. Spall (1997). Accelerated second-order stochastic optimization using only + function measurements, Proceedings of the 36th IEEE Conference on Decision and Control, + 1417-1424 vol.2. `Online at IEEE.org. `_ + + [3]: A. Kandala et al. (2017). Hardware-efficient Variational Quantum Eigensolver for + Small Molecules and Quantum Magnets. Nature 549, pages242–246(2017). + `arXiv:1704.05018v2 `_ + + """ + + def __init__( # pylint: disable=too-many-positional-arguments + self, + maxiter: int = 100, + blocking: bool = False, + allowed_increase: float | None = None, + trust_region: bool = False, + learning_rate: float | np.ndarray | Callable[[], Iterator] | None = None, + perturbation: float | np.ndarray | Callable[[], Iterator] | None = None, + last_avg: int = 1, + resamplings: int | dict[int, int] = 1, + perturbation_dims: int | None = None, + second_order: bool = False, + regularization: float | None = None, + hessian_delay: int = 0, + lse_solver: Callable[[np.ndarray, np.ndarray], np.ndarray] | None = None, + initial_hessian: np.ndarray | None = None, + callback: CALLBACK | None = None, + termination_checker: TERMINATIONCHECKER | None = None, + ) -> None: + r""" + Args: + maxiter: The maximum number of iterations. Note that this is not the maximal number + of function evaluations. + blocking: If True, only accepts updates that improve the loss (up to some allowed + increase, see next argument). + allowed_increase: If ``blocking`` is ``True``, this argument determines by how much + the loss can increase with the proposed parameters and still be accepted. + If ``None``, the allowed increases is calibrated automatically to be twice the + approximated standard deviation of the loss function. + trust_region: If ``True``, restricts the norm of the update step to be :math:`\leq 1`. + learning_rate: The update step is the learning rate is multiplied with the gradient. + If the learning rate is a float, it remains constant over the course of the + optimization. If a NumPy array, the :math:`i`-th element is the learning rate for + the :math:`i`-th iteration. It can also be a callable returning an iterator which + yields the learning rates for each optimization step. + If ``learning_rate`` is set ``perturbation`` must also be provided. + perturbation: Specifies the magnitude of the perturbation for the finite difference + approximation of the gradients. See ``learning_rate`` for the supported types. + If ``perturbation`` is set ``learning_rate`` must also be provided. + last_avg: Return the average of the ``last_avg`` parameters instead of just the + last parameter values. + resamplings: The number of times the gradient (and Hessian) is sampled using a random + direction to construct a gradient estimate. Per default the gradient is estimated + using only one random direction. If an integer, all iterations use the same number + of resamplings. If a dictionary, this is interpreted as + ``{iteration: number of resamplings per iteration}``. + perturbation_dims: The number of perturbed dimensions. Per default, all dimensions + are perturbed, but a smaller, fixed number can be perturbed. If set, the perturbed + dimensions are chosen uniformly at random. + second_order: If True, use 2-SPSA instead of SPSA. In 2-SPSA, the Hessian is estimated + additionally to the gradient, and the gradient is preconditioned with the inverse + of the Hessian to improve convergence. + regularization: To ensure the preconditioner is symmetric and positive definite, the + identity times a small coefficient is added to it. This generator yields that + coefficient. + hessian_delay: Start multiplying the gradient with the inverse Hessian only after a + certain number of iterations. The Hessian is still evaluated and therefore this + argument can be useful to first get a stable average over the last iterations before + using it as preconditioner. + lse_solver: The method to solve for the inverse of the Hessian. Per default an + exact LSE solver is used, but can e.g. be overwritten by a minimization routine. + initial_hessian: The initial guess for the Hessian. By default the identity matrix + is used. + callback: A callback function passed information in each iteration step. The + information is, in this order: the number of function evaluations, the parameters, + the function value, the stepsize, whether the step was accepted. + termination_checker: A callback function executed at the end of each iteration step. The + arguments are, in this order: the parameters, the function value, the number + of function evaluations, the stepsize, whether the step was accepted. If the callback + returns True, the optimization is terminated. + To prevent additional evaluations of the objective method, if the objective has not yet + been evaluated, the objective is estimated by taking the mean of the objective + evaluations used in the estimate of the gradient. + + + Raises: + ValueError: If ``learning_rate`` or ``perturbation`` is an array with less elements + than the number of iterations. + + + """ + super().__init__() + + # general optimizer arguments + self.maxiter = maxiter + self.trust_region = trust_region + self.callback = callback + self.termination_checker = termination_checker + + # if learning rate and perturbation are arrays, check they are sufficiently long + for attr, name in zip([learning_rate, perturbation], ["learning_rate", "perturbation"]): + if isinstance(attr, (list, np.ndarray)): + if len(attr) < maxiter: + raise ValueError(f"Length of {name} is smaller than maxiter ({maxiter}).") + + self.learning_rate = learning_rate + self.perturbation = perturbation + + # SPSA specific arguments + self.blocking = blocking + self.allowed_increase = allowed_increase + self.last_avg = last_avg + self.resamplings = resamplings + self.perturbation_dims = perturbation_dims + + # 2-SPSA specific arguments + if regularization is None: + regularization = 0.01 + + self.second_order = second_order + self.hessian_delay = hessian_delay + self.lse_solver = lse_solver + self.regularization = regularization + self.initial_hessian = initial_hessian + + # runtime arguments + self._nfev: int | None = None # the number of function evaluations + self._smoothed_hessian: np.ndarray | None = None # smoothed average of the Hessians + + @staticmethod + def calibrate( # pylint: disable=too-many-positional-arguments + loss: Callable[[np.ndarray], float], + initial_point: np.ndarray, + c: float = 0.2, + stability_constant: float = 0, + target_magnitude: float | None = None, # 2 pi / 10 + alpha: float = 0.602, + gamma: float = 0.101, + modelspace: bool = False, + max_evals_grouped: int = 1, + ) -> tuple[Callable, Callable]: + r"""Calibrate SPSA parameters with a power series as learning rate and perturbation coeffs. + + The power series are: + + .. math:: + + a_k = \frac{a}{(A + k + 1)^\alpha}, c_k = \frac{c}{(k + 1)^\gamma} + + Args: + loss: The loss function. + initial_point: The initial guess of the iteration. + c: The initial perturbation magnitude. + stability_constant: The value of `A`. + target_magnitude: The target magnitude for the first update step, defaults to + :math:`2\pi / 10`. + alpha: The exponent of the learning rate power series. + gamma: The exponent of the perturbation power series. + modelspace: Whether the target magnitude is the difference of parameter values + or function values (= model space). + max_evals_grouped: The number of grouped evaluations supported by the loss function. + Defaults to 1, i.e. no grouping. + + Returns: + tuple(generator, generator): A tuple of power series generators, the first one for the + learning rate and the second one for the perturbation. + """ + logger.info("SPSA: Starting calibration of learning rate and perturbation.") + if target_magnitude is None: + target_magnitude = 2 * np.pi / 10 + + dim = len(initial_point) + + # compute the average magnitude of the first step + steps = 25 + points = [] + for _ in range(steps): + # compute the random direction + pert = bernoulli_perturbation(dim) + points += [initial_point + c * pert, initial_point - c * pert] + + losses = _batch_evaluate(loss, points, max_evals_grouped) + + avg_magnitudes = 0.0 + for i in range(steps): + delta = losses[2 * i] - losses[2 * i + 1] + avg_magnitudes += np.abs(delta / (2 * c)) + + avg_magnitudes /= steps + + if modelspace: + a = target_magnitude / (avg_magnitudes**2) + else: + a = target_magnitude / avg_magnitudes + + # compute the rescaling factor for correct first learning rate + if a < 1e-10: + warnings.warn(f"Calibration failed, using {target_magnitude} for `a`") + a = target_magnitude + + logger.info("Finished calibration:") + logger.info( + " -- Learning rate: a / ((A + n) ^ alpha) with a = %s, A = %s, alpha = %s", + a, + stability_constant, + alpha, + ) + logger.info(" -- Perturbation: c / (n ^ gamma) with c = %s, gamma = %s", c, gamma) + + # set up the power series + def learning_rate(): + return powerseries(a, alpha, stability_constant) + + def perturbation(): + return powerseries(c, gamma) + + return learning_rate, perturbation + + @staticmethod + def estimate_stddev( + loss: Callable[[np.ndarray], float], + initial_point: np.ndarray, + avg: int = 25, + max_evals_grouped: int = 1, + ) -> float: + """Estimate the standard deviation of the loss function.""" + losses = _batch_evaluate(loss, avg * [initial_point], max_evals_grouped) + return np.std(losses) + + @property + def settings(self) -> dict[str, Any]: + # if learning rate or perturbation are custom iterators expand them + if callable(self.learning_rate): + iterator = self.learning_rate() + learning_rate = np.array([next(iterator) for _ in range(self.maxiter)]) + else: + learning_rate = self.learning_rate # type: ignore[assignment] + + if callable(self.perturbation): + iterator = self.perturbation() + perturbation = np.array([next(iterator) for _ in range(self.maxiter)]) + else: + perturbation = self.perturbation # type: ignore[assignment] + + return { + "maxiter": self.maxiter, + "learning_rate": learning_rate, + "perturbation": perturbation, + "trust_region": self.trust_region, + "blocking": self.blocking, + "allowed_increase": self.allowed_increase, + "resamplings": self.resamplings, + "perturbation_dims": self.perturbation_dims, + "second_order": self.second_order, + "hessian_delay": self.hessian_delay, + "regularization": self.regularization, + "lse_solver": self.lse_solver, + "initial_hessian": self.initial_hessian, + "callback": self.callback, + "termination_checker": self.termination_checker, + } + + def _point_sample( + self, loss, x, eps, delta1, delta2 + ): # pylint: disable=too-many-positional-arguments + """A single sample of the gradient at position ``x`` in direction ``delta``.""" + # points to evaluate + points = [x + eps * delta1, x - eps * delta1] + self._nfev += 2 + + if self.second_order: + points += [x + eps * (delta1 + delta2), x + eps * (-delta1 + delta2)] + self._nfev += 2 + + # batch evaluate the points (if possible) + values = _batch_evaluate(loss, points, self._max_evals_grouped) + + plus = values[0] + minus = values[1] + gradient_sample = (plus - minus) / (2 * eps) * delta1 + + hessian_sample = None + if self.second_order: + diff = (values[2] - plus) - (values[3] - minus) + diff /= 2 * eps**2 + + rank_one = np.outer(delta1, delta2) + hessian_sample = diff * (rank_one + rank_one.T) / 2 + + return np.mean(values), gradient_sample, hessian_sample + + def _point_estimate(self, loss, x, eps, num_samples): + """The gradient estimate at point x.""" + # set up variables to store averages + value_estimate = 0 + gradient_estimate = np.zeros(x.size) + hessian_estimate = np.zeros((x.size, x.size)) + + # iterate over the directions + deltas1 = [ + bernoulli_perturbation(x.size, self.perturbation_dims) for _ in range(num_samples) + ] + + if self.second_order: + deltas2 = [ + bernoulli_perturbation(x.size, self.perturbation_dims) for _ in range(num_samples) + ] + else: + deltas2 = None + + for i in range(num_samples): + delta1 = deltas1[i] + delta2 = deltas2[i] if self.second_order else None + + value_sample, gradient_sample, hessian_sample = self._point_sample( + loss, x, eps, delta1, delta2 + ) + value_estimate += value_sample + gradient_estimate += gradient_sample + + if self.second_order: + hessian_estimate += hessian_sample + + return ( + value_estimate / num_samples, + gradient_estimate / num_samples, + hessian_estimate / num_samples, + ) + + def _compute_update( + self, loss, x, k, eps, lse_solver + ): # pylint: disable=too-many-positional-arguments + # compute the perturbations + if isinstance(self.resamplings, dict): + num_samples = self.resamplings.get(k, 1) + else: + num_samples = self.resamplings + + # accumulate the number of samples + value, gradient, hessian = self._point_estimate(loss, x, eps, num_samples) + + # precondition gradient with inverse Hessian, if specified + if self.second_order: + smoothed = k / (k + 1) * self._smoothed_hessian + 1 / (k + 1) * hessian + self._smoothed_hessian = smoothed + + if k > self.hessian_delay: + spd_hessian = _make_spd(smoothed, self.regularization) + + # solve for the gradient update + gradient = np.real(lse_solver(spd_hessian, gradient)) + + return value, gradient + + def minimize( + self, + fun: Callable[[POINT], float], + x0: POINT, + jac: Callable[[POINT], POINT] | None = None, + bounds: list[tuple[float, float]] | None = None, + ) -> OptimizerResult: + # ensure learning rate and perturbation are correctly set: either none or both + # this happens only here because for the calibration the loss function is required + x0 = np.asarray(x0) + if self.learning_rate is None and self.perturbation is None: + get_eta, get_eps = self.calibrate(fun, x0, max_evals_grouped=self._max_evals_grouped) + else: + get_eta, get_eps = _validate_pert_and_learningrate( + self.perturbation, self.learning_rate + ) + eta, eps = get_eta(), get_eps() + + lse_solver = self.lse_solver + if self.lse_solver is None: + lse_solver = np.linalg.solve + + # prepare some initials + x = np.asarray(x0) + if self.initial_hessian is None: + self._smoothed_hessian = np.identity(x.size) + else: + self._smoothed_hessian = self.initial_hessian + + self._nfev = 0 + + # if blocking is enabled we need to keep track of the function values + if self.blocking: + fx = fun(x) # pylint: disable=invalid-name + + self._nfev += 1 + if self.allowed_increase is None: + self.allowed_increase = 2 * self.estimate_stddev( + fun, x, max_evals_grouped=self._max_evals_grouped + ) + + logger.info("SPSA: Starting optimization.") + start = time() + + # keep track of the last few steps to return their average + last_steps = deque([x]) + + # use a local variable and while loop to keep track of the number of iterations + # if the termination checker terminates early + k = 0 + while k < self.maxiter: + k += 1 + iteration_start = time() + # compute update + fx_estimate, update = self._compute_update(fun, x, k, next(eps), lse_solver) + + # trust region + if self.trust_region: + norm = np.linalg.norm(update) + if norm > 1: # stop from dividing by 0 + update = update / norm + + # compute next parameter value + update = update * next(eta) + x_next = x - update + fx_next = None + + # blocking + if self.blocking: + self._nfev += 1 + fx_next = fun(x_next) + + if fx + self.allowed_increase <= fx_next: # accept only if loss improved + if self.callback is not None: + self.callback( + self._nfev, # number of function evals + x_next, # next parameters + fx_next, # loss at next parameters + np.linalg.norm(update), # size of the update step + False, + ) # not accepted + + logger.info( + "Iteration %s/%s rejected in %s.", + k, + self.maxiter + 1, + time() - iteration_start, + ) + continue + fx = fx_next # pylint: disable=invalid-name + + logger.info( + "Iteration %s/%s done in %s.", k, self.maxiter + 1, time() - iteration_start + ) + + if self.callback is not None: + # if we didn't evaluate the function yet, do it now + if not self.blocking: + self._nfev += 1 + fx_next = fun(x_next) + + self.callback( + self._nfev, # number of function evals + x_next, # next parameters + fx_next, # loss at next parameters + np.linalg.norm(update), # size of the update step + True, + ) # accepted + + # update parameters + x = x_next + + # update the list of the last ``last_avg`` parameters + if self.last_avg > 1: + last_steps.append(x_next) + if len(last_steps) > self.last_avg: + last_steps.popleft() + + if self.termination_checker is not None: + fx_check = fx_estimate if fx_next is None else fx_next + if self.termination_checker( + self._nfev, x_next, fx_check, np.linalg.norm(update), True + ): + logger.info("terminated optimization at {k}/{self.maxiter} iterations") + break + + logger.info("SPSA: Finished in %s", time() - start) + + if self.last_avg > 1: + x = np.mean(np.asarray(last_steps), axis=0) + + result = OptimizerResult() + result.x = x + result.fun = fun(x) + result.nfev = self._nfev + result.nit = k + + return result + + def get_support_level(self): + """Get the support level dictionary.""" + return { + "gradient": OptimizerSupportLevel.ignored, + "bounds": OptimizerSupportLevel.ignored, + "initial_point": OptimizerSupportLevel.required, + } + + +def bernoulli_perturbation(dim, perturbation_dims=None): + """Get a Bernoulli random perturbation.""" + if perturbation_dims is None: + return 1 - 2 * algorithm_globals.random.binomial(1, 0.5, size=dim) + + pert = 1 - 2 * algorithm_globals.random.binomial(1, 0.5, size=perturbation_dims) + indices = algorithm_globals.random.choice( + list(range(dim)), size=perturbation_dims, replace=False + ) + result = np.zeros(dim) + result[indices] = pert + + return result + + +def powerseries(eta=0.01, power=2, offset=0): + """Yield a series decreasing by a power law.""" + + n = 1 + while True: + yield eta / ((n + offset) ** power) + n += 1 + + +def constant(eta=0.01): + """Yield a constant series.""" + + while True: + yield eta + + +def _batch_evaluate(function, points, max_evals_grouped, unpack_points=False): + """Evaluate a function on all points with batches of max_evals_grouped. + + The points are a list of inputs, as ``[in1, in2, in3, ...]``. If the individual + inputs are tuples (because the function takes multiple inputs), set ``unpack_points`` to ``True``. + """ + + # if the function cannot handle lists of points as input, cover this case immediately + if max_evals_grouped is None or max_evals_grouped == 1: + # support functions with multiple arguments where the points are given in a tuple + return [ + function(*point) if isinstance(point, tuple) else function(point) for point in points + ] + + num_points = len(points) + + # get the number of batches + num_batches = num_points // max_evals_grouped + if num_points % max_evals_grouped != 0: + num_batches += 1 + + # split the points + batched_points = np.array_split(np.asarray(points), num_batches) + + results = [] + for batch in batched_points: + 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)) + return psd + bias * identity + + +def _validate_pert_and_learningrate(perturbation, learning_rate): + if learning_rate is None or perturbation is None: + raise ValueError("If one of learning rate or perturbation is set, both must be set.") + + if isinstance(perturbation, float): + + def get_eps(): + return constant(perturbation) + + elif isinstance(perturbation, (list, np.ndarray)): + + def get_eps(): + return iter(perturbation) + + else: + get_eps = perturbation + + if isinstance(learning_rate, float): + + def get_eta(): + return constant(learning_rate) + + elif isinstance(learning_rate, (list, np.ndarray)): + + def get_eta(): + return iter(learning_rate) + + else: + get_eta = learning_rate + + return get_eta, get_eps diff --git a/qiskit_optimization/problems/quadratic_program.py b/qiskit_optimization/problems/quadratic_program.py index 3ab177c27..7e5c19e6d 100644 --- a/qiskit_optimization/problems/quadratic_program.py +++ b/qiskit_optimization/problems/quadratic_program.py @@ -364,7 +364,12 @@ def continuous_var_dict( # pylint: disable=too-many-positional-arguments nested substitution. """ return self._var_dict( - keys, lowerbound, upperbound, Variable.Type.CONTINUOUS, name, key_format + keys=keys, + lowerbound=lowerbound, + upperbound=upperbound, + vartype=Variable.Type.CONTINUOUS, + name=name, + key_format=key_format, ) def continuous_var_list( # pylint: disable=too-many-positional-arguments @@ -442,7 +447,14 @@ def binary_var_dict( QiskitOptimizationError: if `key_format` has more than one substitution or a nested substitution. """ - return self._var_dict(keys, 0, 1, Variable.Type.BINARY, name, key_format) + return self._var_dict( + keys=keys, + lowerbound=0, + upperbound=1, + vartype=Variable.Type.BINARY, + name=name, + key_format=key_format, + ) def binary_var_list( self, @@ -524,7 +536,14 @@ def integer_var_dict( # pylint: disable=too-many-positional-arguments QiskitOptimizationError: if `key_format` has more than one substitution or a nested substitution. """ - return self._var_dict(keys, lowerbound, upperbound, Variable.Type.INTEGER, name, key_format) + return self._var_dict( + keys=keys, + lowerbound=lowerbound, + upperbound=upperbound, + vartype=Variable.Type.INTEGER, + name=name, + key_format=key_format, + ) def integer_var_list( # pylint: disable=too-many-positional-arguments self, diff --git a/qiskit_optimization/utils/__init__.py b/qiskit_optimization/utils/__init__.py new file mode 100644 index 000000000..bb76dcc49 --- /dev/null +++ b/qiskit_optimization/utils/__init__.py @@ -0,0 +1,23 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 .algorithm_globals import algorithm_globals +from .validate_initial_point import validate_initial_point +from .validate_bounds import validate_bounds + +__all__ = [ + "algorithm_globals", + "validate_initial_point", + "validate_bounds", +] diff --git a/qiskit_optimization/utils/algorithm_globals.py b/qiskit_optimization/utils/algorithm_globals.py new file mode 100644 index 000000000..c762dba6e --- /dev/null +++ b/qiskit_optimization/utils/algorithm_globals.py @@ -0,0 +1,131 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2019, 2024. +# +# 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. + +""" +utils.algorithm_globals +======================= +Common (global) properties used across qiskit_algorithms. + +.. currentmodule:: qiskit_algorithms.utils.algorithm_globals + +Includes: + + * Random number generator and random seed. + + Algorithms can use the generator for random values, as needed, and it + can be seeded here for reproducible results when using such an algorithm. + This is often important, for example in unit tests, where the same + outcome is desired each time (reproducible) and not have it be variable + due to randomness. + +Attributes: + random_seed (int | None): Random generator seed (read/write). + random (np.random.Generator): Random generator (read-only) +""" + +from __future__ import annotations + +import warnings + +import numpy as np + + +class QiskitAlgorithmGlobals: + """Global properties for algorithms.""" + + # The code is done to work even after some future removal of algorithm_globals + # from Qiskit (qiskit.utils). All that is needed in the future, after that, if + # this is updated, is just the logic in the except blocks. + # + # If the Qiskit version exists this acts a redirect to that (it delegates the + # calls off to it). In the future when that does not exist this has similar code + # in the except blocks here, as noted above, that will take over. By delegating + # to the Qiskit instance it means that any existing code that uses that continues + # to work. Logic here in qiskit_algorithms though uses this instance and the + # random check here has logic to warn if the seed here is not the same as the Qiskit + # version so we can detect direct usage of the Qiskit version and alert the user to + # change their code to use this. So simply changing from: + # from qiskit.utils import algorithm_globals + # to + # from qiskit_algorithm.utils import algorithm_globals + + def __init__(self) -> None: + self._random_seed: int | None = None + self._random: np.random.Generator | None = None + + @property + def random_seed(self) -> int | None: + """Random seed property (getter/setter).""" + try: + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=DeprecationWarning) + + from qiskit.utils import algorithm_globals as qiskit_globals + + return qiskit_globals.random_seed + + except ImportError: + return self._random_seed + + @random_seed.setter + def random_seed(self, seed: int | None) -> None: + """Set the random generator seed. + + Args: + seed: If ``None`` then internally a random value is used as a seed + """ + try: + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=DeprecationWarning) + + from qiskit.utils import algorithm_globals as qiskit_globals + + qiskit_globals.random_seed = seed + # Mirror the seed here when set via this random_seed. If the seed is + # set on the qiskit.utils instance then we can detect it's different + self._random_seed = seed + + except ImportError: + self._random_seed = seed + self._random = None + + @property + def random(self) -> np.random.Generator: + """Return a numpy np.random.Generator (default_rng) using random_seed.""" + try: + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=DeprecationWarning) + + from qiskit.utils import algorithm_globals as qiskit_globals + + if self._random_seed != qiskit_globals.random_seed: + # If the seeds are different - likely this local is None and the qiskit.utils + # algorithms global was seeded directly then we will warn to use this here as + # the Qiskit version is planned to be removed in a future version of Qiskit. + warnings.warn( + "Using random that is seeded via qiskit.utils algorithm_globals is deprecated " + "since version 0.2.0. Instead set random_seed directly to " + "qiskit_algorithms.utils algorithm_globals.", + category=DeprecationWarning, + stacklevel=2, + ) + + return qiskit_globals.random + + except ImportError: + if self._random is None: + self._random = np.random.default_rng(self._random_seed) + return self._random + + +# Global instance to be used as the entry point for globals. +algorithm_globals = QiskitAlgorithmGlobals() diff --git a/qiskit_optimization/utils/set_batching.py b/qiskit_optimization/utils/set_batching.py new file mode 100644 index 000000000..1891b2f07 --- /dev/null +++ b/qiskit_optimization/utils/set_batching.py @@ -0,0 +1,27 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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. + +"""Set default batch sizes for the optimizers.""" + +from qiskit_optimization.optimizers import Optimizer, SPSA + + +def _set_default_batchsize(optimizer: Optimizer) -> bool: + """Set the default batchsize, if None is set and return whether it was updated or not.""" + if isinstance(optimizer, SPSA): + updated = optimizer._max_evals_grouped is None + if updated: + optimizer.set_max_evals_grouped(50) + else: # we only set a batchsize for SPSA + updated = False + + return updated diff --git a/qiskit_optimization/utils/validate_bounds.py b/qiskit_optimization/utils/validate_bounds.py new file mode 100644 index 000000000..2affd7b8c --- /dev/null +++ b/qiskit_optimization/utils/validate_bounds.py @@ -0,0 +1,44 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 diff --git a/qiskit_optimization/utils/validate_initial_point.py b/qiskit_optimization/utils/validate_initial_point.py new file mode 100644 index 000000000..f9a96173b --- /dev/null +++ b/qiskit_optimization/utils/validate_initial_point.py @@ -0,0 +1,65 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 + +import numpy as np + +from qiskit.circuit import QuantumCircuit +from qiskit_optimization.utils.algorithm_globals import algorithm_globals + + +def validate_initial_point(point: np.ndarray | None | None, circuit: QuantumCircuit) -> np.ndarray: + 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 diff --git a/qiskit_optimization/utils/validation.py b/qiskit_optimization/utils/validation.py new file mode 100644 index 000000000..edd8e06e1 --- /dev/null +++ b/qiskit_optimization/utils/validation.py @@ -0,0 +1,138 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2019, 2024. +# +# 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. + +""" +Validation module +""" + +from typing import Set + + +def validate_in_set(name: str, value: object, values: Set[object]) -> None: + """ + Args: + name: value name. + value: value to check. + values: set that should contain value. + Raises: + ValueError: invalid value + """ + if value not in values: + raise ValueError(f"{name} must be one of '{values}', was '{value}'.") + + +def validate_min(name: str, value: float, minimum: float) -> None: + """ + Args: + name: value name. + value: value to check. + minimum: minimum value allowed. + Raises: + ValueError: invalid value + """ + if value < minimum: + raise ValueError(f"{name} must have value >= {minimum}, was {value}") + + +def validate_min_exclusive(name: str, value: float, minimum: float) -> None: + """ + Args: + name: value name. + value: value to check. + minimum: minimum value allowed. + Raises: + ValueError: invalid value + """ + if value <= minimum: + raise ValueError(f"{name} must have value > {minimum}, was {value}") + + +def validate_max(name: str, value: float, maximum: float) -> None: + """ + Args: + name: value name. + value: value to check. + maximum: maximum value allowed. + Raises: + ValueError: invalid value + """ + if value > maximum: + raise ValueError(f"{name} must have value <= {maximum}, was {value}") + + +def validate_max_exclusive(name: str, value: float, maximum: float) -> None: + """ + Args: + name: value name. + value: value to check. + maximum: maximum value allowed. + Raises: + ValueError: invalid value + """ + if value >= maximum: + raise ValueError(f"{name} must have value < {maximum}, was {value}") + + +def validate_range(name: str, value: float, minimum: float, maximum: float) -> None: + """ + Args: + name: value name. + value: value to check. + minimum: minimum value allowed. + maximum: maximum value allowed. + Raises: + ValueError: invalid value + """ + if value < minimum or value > maximum: + raise ValueError(f"{name} must have value >= {minimum} and <= {maximum}, was {value}") + + +def validate_range_exclusive(name: str, value: float, minimum: float, maximum: float) -> None: + """ + Args: + name: value name. + value: value to check. + minimum: minimum value allowed. + maximum: maximum value allowed. + Raises: + ValueError: invalid value + """ + if value <= minimum or value >= maximum: + raise ValueError(f"{name} must have value > {minimum} and < {maximum}, was {value}") + + +def validate_range_exclusive_min(name: str, value: float, minimum: float, maximum: float) -> None: + """ + Args: + name: value name. + value: value to check. + minimum: minimum value allowed. + maximum: maximum value allowed. + Raises: + ValueError: invalid value + """ + if value <= minimum or value > maximum: + raise ValueError(f"{name} must have value > {minimum} and <= {maximum}, was {value}") + + +def validate_range_exclusive_max(name: str, value: float, minimum: float, maximum: float) -> None: + """ + Args: + name: value name. + value: value to check. + minimum: minimum value allowed. + maximum: maximum value allowed. + Raises: + ValueError: invalid value + """ + if value < minimum or value >= maximum: + raise ValueError(f"{name} must have value >= {minimum} and < {maximum}, was {value}") diff --git a/qiskit_optimization/variational_algorithm.py b/qiskit_optimization/variational_algorithm.py new file mode 100644 index 000000000..ad5ddd30a --- /dev/null +++ b/qiskit_optimization/variational_algorithm.py @@ -0,0 +1,137 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2019, 2024. +# +# 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 Algorithm Base Class. + +This class can be used an interface for working with Variation Algorithms, such as VQE, +QAOA, or QSVM, and also provides helper utilities for implementing new variational algorithms. +Writing a new variational algorithm is a simple as extending this class, implementing a cost +function for the new algorithm to pass to the optimizer, and running :meth:`find_minimum` method +of this class to carry out the optimization. Alternatively, all of the functions below can be +overridden to opt-out of this infrastructure but still meet the interface requirements. + +.. note:: + + This component has some function that is normally random. If you want to reproduce behavior + then you should set the random number generator seed in the algorithm_globals + (``qiskit_algorithms.utils.algorithm_globals.random_seed = seed``). +""" + +from __future__ import annotations +from abc import ABC, abstractmethod +import numpy as np + +from qiskit.circuit import QuantumCircuit + +from .algorithm_result import AlgorithmResult +from .optimizers import OptimizerResult + + +class VariationalAlgorithm(ABC): + """The Variational Algorithm Base Class.""" + + @property + @abstractmethod + def initial_point(self) -> np.ndarray | None: + """Returns initial point.""" + pass + + @initial_point.setter + @abstractmethod + def initial_point(self, initial_point: np.ndarray | None) -> None: + """Sets initial point.""" + pass + + +class VariationalResult(AlgorithmResult): + """Variation Algorithm Result.""" + + def __init__(self) -> None: + super().__init__() + self._optimizer_evals: int | None = None + self._optimizer_time: float | None = None + self._optimal_value: float | None = None + self._optimal_point: np.ndarray | None = None + self._optimal_parameters: dict | None = None + self._optimizer_result: OptimizerResult | None = None + self._optimal_circuit: QuantumCircuit | None = None + + @property + def optimizer_evals(self) -> int | None: + """Returns number of optimizer evaluations""" + return self._optimizer_evals + + @optimizer_evals.setter + def optimizer_evals(self, value: int) -> None: + """Sets number of optimizer evaluations""" + self._optimizer_evals = value + + @property + def optimizer_time(self) -> float | None: + """Returns time taken for optimization""" + return self._optimizer_time + + @optimizer_time.setter + def optimizer_time(self, value: float) -> None: + """Sets time taken for optimization""" + self._optimizer_time = value + + @property + def optimal_value(self) -> float | None: + """Returns optimal value""" + return self._optimal_value + + @optimal_value.setter + def optimal_value(self, value: int) -> None: + """Sets optimal value""" + self._optimal_value = value + + @property + def optimal_point(self) -> np.ndarray | None: + """Returns optimal point""" + return self._optimal_point + + @optimal_point.setter + def optimal_point(self, value: np.ndarray) -> None: + """Sets optimal point""" + self._optimal_point = value + + @property + def optimal_parameters(self) -> dict | None: + """Returns the optimal parameters in a dictionary""" + return self._optimal_parameters + + @optimal_parameters.setter + def optimal_parameters(self, value: dict) -> None: + """Sets optimal parameters""" + self._optimal_parameters = value + + @property + def optimizer_result(self) -> OptimizerResult | None: + """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 + + @property + def optimal_circuit(self) -> QuantumCircuit: + """The optimal circuits. Along with the optimal parameters, + these can be used to retrieve the minimum eigenstate. + """ + return self._optimal_circuit + + @optimal_circuit.setter + def optimal_circuit(self, optimal_circuit: QuantumCircuit) -> None: + self._optimal_circuit = optimal_circuit diff --git a/requirements-dev.txt b/requirements-dev.txt index fb91bae47..06cccd3ca 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -16,3 +16,4 @@ mypy>=0.981 mypy-extensions>=0.4.3 nbsphinx qiskit_sphinx_theme~=1.16.0 +rustworkx diff --git a/test/__init__.py b/test/__init__.py index 8ab19fd32..25b3fb172 100755 --- a/test/__init__.py +++ b/test/__init__.py @@ -1,6 +1,6 @@ # This code is part of a Qiskit project. # -# (C) Copyright IBM 2018, 2023. +# (C) Copyright IBM 2018, 2024. # # 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 @@ -13,5 +13,6 @@ """ Optimization test packages """ from .optimization_test_case import QiskitOptimizationTestCase +from .algorithms_test_case import QiskitAlgorithmsTestCase -__all__ = ["QiskitOptimizationTestCase"] +__all__ = ["QiskitOptimizationTestCase", "QiskitAlgorithmsTestCase"] diff --git a/test/algorithms_test_case.py b/test/algorithms_test_case.py new file mode 100644 index 000000000..684046af6 --- /dev/null +++ b/test/algorithms_test_case.py @@ -0,0 +1,88 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2018, 2024. +# +# 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. + +"""Algorithms Test Case""" + +from typing import Optional +from abc import ABC +import warnings +import inspect +import logging +import os +import unittest +import time + +from qiskit_optimization.utils import algorithm_globals + +# disable deprecation warnings that can cause log output overflow +# pylint: disable=unused-argument + + +def _noop(*args, **kargs): + pass + + +# disable warning messages +# warnings.warn = _noop + + +class QiskitAlgorithmsTestCase(unittest.TestCase, ABC): + """Optimization Test Case""" + + moduleName = None + log = None + + def setUp(self) -> None: + warnings.filterwarnings("default", category=DeprecationWarning) + self._started_at = time.time() + self._class_location = __file__ + + def tearDown(self) -> None: + algorithm_globals.random_seed = None + elapsed = time.time() - self._started_at + if elapsed > 5.0: + print(f"({round(elapsed, 2):.2f}s)", flush=True) + + @classmethod + def setUpClass(cls) -> None: + cls.moduleName = os.path.splitext(inspect.getfile(cls))[0] + cls.log = logging.getLogger(cls.__name__) + + # Set logging to file and stdout if the LOG_LEVEL environment variable + # is set. + if os.getenv("LOG_LEVEL"): + # Set up formatter. + log_fmt = f"{cls.__name__}.%(funcName)s:%(levelname)s:%(asctime)s:" " %(message)s" + formatter = logging.Formatter(log_fmt) + + # Set up the file handler. + log_file_name = f"{cls.moduleName}.log" + file_handler = logging.FileHandler(log_file_name) + file_handler.setFormatter(formatter) + cls.log.addHandler(file_handler) + + # Set the logging level from the environment variable, defaulting + # to INFO if it is not a valid level. + level = logging._nameToLevel.get(os.getenv("LOG_LEVEL"), logging.INFO) + cls.log.setLevel(level) + + def get_resource_path(self, filename: str, path: Optional[str] = None) -> str: + """Get the absolute path to a resource. + Args: + filename: filename or relative path to the resource. + path: path used as relative to the filename. + Returns: + str: the absolute path to the resource. + """ + root = os.path.dirname(self._class_location) + path = root if path is None else os.path.join(root, path) + return os.path.normpath(os.path.join(path, filename)) diff --git a/test/eigensolvers/__init__.py b/test/eigensolvers/__init__.py new file mode 100644 index 000000000..ce4a89ca8 --- /dev/null +++ b/test/eigensolvers/__init__.py @@ -0,0 +1,13 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Tests for the eigensolvers.""" diff --git a/test/eigensolvers/test_numpy_eigensolver.py b/test/eigensolvers/test_numpy_eigensolver.py new file mode 100644 index 000000000..8d57dbb93 --- /dev/null +++ b/test/eigensolvers/test_numpy_eigensolver.py @@ -0,0 +1,214 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2018, 2024. +# +# 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 NumPyEigensolver""" + +import unittest +from test import QiskitAlgorithmsTestCase + +import numpy as np +from ddt import data, ddt + +from qiskit.quantum_info import Operator, SparsePauliOp, Pauli, ScalarOp + +from qiskit_optimization.eigensolvers import NumPyEigensolver +from qiskit_optimization import AlgorithmError + +H2_SPARSE_PAULI = SparsePauliOp( + ["II", "ZI", "IZ", "ZZ", "XX"], + coeffs=[ + -1.052373245772859, + 0.39793742484318045, + -0.39793742484318045, + -0.01128010425623538, + 0.18093119978423156, + ], +) + +H2_OP = Operator(H2_SPARSE_PAULI.to_matrix()) + + +@ddt +class TestNumPyEigensolver(QiskitAlgorithmsTestCase): + """Test NumPy Eigen solver""" + + @data(H2_SPARSE_PAULI, H2_OP) + def test_ce(self, op): + """Test basics""" + algo = NumPyEigensolver() + result = algo.compute_eigenvalues(operator=op, aux_operators=[]) + self.assertEqual(len(result.eigenvalues), 1) + self.assertEqual(len(result.eigenstates), 1) + self.assertEqual(result.eigenvalues.dtype, np.float64) + self.assertAlmostEqual(result.eigenvalues[0], -1.85727503) + + @data(H2_SPARSE_PAULI, H2_OP) + def test_ce_k4(self, op): + """Test for k=4 eigenvalues""" + algo = NumPyEigensolver(k=4) + result = algo.compute_eigenvalues(operator=op, aux_operators=[]) + self.assertEqual(len(result.eigenvalues), 4) + self.assertEqual(len(result.eigenstates), 4) + self.assertEqual(result.eigenvalues.dtype, np.float64) + np.testing.assert_array_almost_equal( + result.eigenvalues, [-1.85727503, -1.24458455, -0.88272215, -0.22491125] + ) + + @data(H2_SPARSE_PAULI, H2_OP) + def test_ce_k4_filtered(self, op): + """Test for k=4 eigenvalues with filter""" + + # define filter criterion + # pylint: disable=unused-argument + def criterion(x, v, a_v): + return v >= -1 + + algo = NumPyEigensolver(k=4, filter_criterion=criterion) + result = algo.compute_eigenvalues(operator=op, aux_operators=[]) + self.assertEqual(len(result.eigenvalues), 2) + self.assertEqual(len(result.eigenstates), 2) + self.assertEqual(result.eigenvalues.dtype, np.float64) + np.testing.assert_array_almost_equal(result.eigenvalues, [-0.88272215, -0.22491125]) + + @data(H2_SPARSE_PAULI, H2_OP) + def test_ce_k4_filtered_empty(self, op): + """Test for k=4 eigenvalues with filter always returning False""" + + # define filter criterion + # pylint: disable=unused-argument + def criterion(x, v, a_v): + return False + + algo = NumPyEigensolver(k=4, filter_criterion=criterion) + result = algo.compute_eigenvalues(operator=op, aux_operators=[]) + self.assertEqual(len(result.eigenvalues), 0) + self.assertEqual(len(result.eigenstates), 0) + + @data( + SparsePauliOp(["X"], coeffs=[1.0]), + SparsePauliOp(["Y"], coeffs=[1.0]), + SparsePauliOp(["Z"], coeffs=[1.0]), + ) + def test_ce_k1_1q(self, op): + """Test for 1 qubit operator""" + algo = NumPyEigensolver(k=1) + result = algo.compute_eigenvalues(operator=op) + np.testing.assert_array_almost_equal(result.eigenvalues, [-1]) + + @data( + SparsePauliOp(["X"], coeffs=[1.0]), + SparsePauliOp(["Y"], coeffs=[1.0]), + SparsePauliOp(["Z"], coeffs=[1.0]), + ) + def test_ce_k2_1q(self, op): + """Test for 1 qubit operator""" + algo = NumPyEigensolver(k=2) + result = algo.compute_eigenvalues(operator=op) + np.testing.assert_array_almost_equal(result.eigenvalues, [-1, 1]) + + @data(H2_SPARSE_PAULI, H2_OP) + def test_aux_operators_list(self, op): + """Test list-based aux_operators.""" + aux_op1 = Operator(SparsePauliOp(["II"], coeffs=[2.0]).to_matrix()) + aux_op2 = SparsePauliOp(["II", "ZZ", "YY", "XX"], coeffs=[0.5, 0.5, 0.5, -0.5]) + aux_ops = [aux_op1, aux_op2] + algo = NumPyEigensolver() + result = algo.compute_eigenvalues(operator=op, aux_operators=aux_ops) + self.assertEqual(len(result.eigenvalues), 1) + self.assertEqual(len(result.eigenstates), 1) + self.assertEqual(result.eigenvalues.dtype, np.float64) + self.assertAlmostEqual(result.eigenvalues[0], -1.85727503) + self.assertEqual(len(result.aux_operators_evaluated), 1) + self.assertEqual(len(result.aux_operators_evaluated[0]), 2) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][0], 2, places=6) + self.assertAlmostEqual(result.aux_operators_evaluated[0][1][0], 0, places=6) + # metadata + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][1].pop("variance"), 0.0) + self.assertAlmostEqual(result.aux_operators_evaluated[0][1][1].pop("variance"), 0.0) + + # Go again with additional None and zero operators + extra_ops = [*aux_ops, None, 0] + result = algo.compute_eigenvalues(operator=op, aux_operators=extra_ops) + self.assertEqual(len(result.eigenvalues), 1) + self.assertEqual(len(result.eigenstates), 1) + self.assertEqual(result.eigenvalues.dtype, np.float64) + self.assertAlmostEqual(result.eigenvalues[0], -1.85727503) + self.assertEqual(len(result.aux_operators_evaluated), 1) + self.assertEqual(len(result.aux_operators_evaluated[0]), 4) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][0], 2, places=6) + self.assertAlmostEqual(result.aux_operators_evaluated[0][1][0], 0, places=6) + self.assertIsNone(result.aux_operators_evaluated[0][2], None) + self.assertEqual(result.aux_operators_evaluated[0][3][0], 0.0) + # metadata + self.assertAlmostEqual(result.aux_operators_evaluated[0][0][1].pop("variance"), 0.0) + self.assertAlmostEqual(result.aux_operators_evaluated[0][1][1].pop("variance"), 0.0) + self.assertEqual(result.aux_operators_evaluated[0][3][1].pop("variance"), 0.0) + + @data(H2_SPARSE_PAULI, H2_OP) + def test_aux_operators_dict(self, op): + """Test dict-based aux_operators.""" + aux_op1 = Operator(SparsePauliOp(["II"], coeffs=[2.0]).to_matrix()) + aux_op2 = SparsePauliOp(["II", "ZZ", "YY", "XX"], coeffs=[0.5, 0.5, 0.5, -0.5]) + aux_ops = {"aux_op1": aux_op1, "aux_op2": aux_op2} + algo = NumPyEigensolver() + result = algo.compute_eigenvalues(operator=op, aux_operators=aux_ops) + self.assertEqual(len(result.eigenvalues), 1) + self.assertEqual(len(result.eigenstates), 1) + self.assertEqual(result.eigenvalues.dtype, np.float64) + self.assertAlmostEqual(result.eigenvalues[0], -1.85727503) + self.assertEqual(len(result.aux_operators_evaluated), 1) + self.assertEqual(len(result.aux_operators_evaluated[0]), 2) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0]["aux_op1"][0], 2, places=6) + self.assertAlmostEqual(result.aux_operators_evaluated[0]["aux_op2"][0], 0, places=6) + # metadata + self.assertAlmostEqual(result.aux_operators_evaluated[0]["aux_op1"][1].pop("variance"), 0.0) + self.assertAlmostEqual(result.aux_operators_evaluated[0]["aux_op2"][1].pop("variance"), 0.0) + + # Go again with additional None and zero operators + extra_ops = {**aux_ops, "None_operator": None, "zero_operator": 0} + result = algo.compute_eigenvalues(operator=op, aux_operators=extra_ops) + self.assertEqual(len(result.eigenvalues), 1) + self.assertEqual(len(result.eigenstates), 1) + self.assertEqual(result.eigenvalues.dtype, np.float64) + self.assertAlmostEqual(result.eigenvalues[0], -1.85727503) + self.assertEqual(len(result.aux_operators_evaluated), 1) + self.assertEqual(len(result.aux_operators_evaluated[0]), 3) + # expectation values + self.assertAlmostEqual(result.aux_operators_evaluated[0]["aux_op1"][0], 2, places=6) + self.assertAlmostEqual(result.aux_operators_evaluated[0]["aux_op2"][0], 0, places=6) + self.assertEqual(result.aux_operators_evaluated[0]["zero_operator"][0], 0.0) + self.assertTrue("None_operator" not in result.aux_operators_evaluated[0].keys()) + # metadata + self.assertAlmostEqual(result.aux_operators_evaluated[0]["aux_op1"][1].pop("variance"), 0.0) + self.assertAlmostEqual(result.aux_operators_evaluated[0]["aux_op2"][1].pop("variance"), 0.0) + self.assertAlmostEqual( + result.aux_operators_evaluated[0]["zero_operator"][1].pop("variance"), 0.0 + ) + + def test_pauli_op(self): + """Test simple pauli operator""" + algo = NumPyEigensolver(k=1) + result = algo.compute_eigenvalues(operator=Pauli("X")) + np.testing.assert_array_almost_equal(result.eigenvalues, [-1]) + + def test_scalar_op(self): + """Test scalar operator""" + algo = NumPyEigensolver(k=1) + with self.assertRaises(AlgorithmError): + algo.compute_eigenvalues(operator=ScalarOp(1)) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/minimum_eigensolvers/__init__.py b/test/minimum_eigensolvers/__init__.py new file mode 100644 index 000000000..c870af8ba --- /dev/null +++ b/test/minimum_eigensolvers/__init__.py @@ -0,0 +1,11 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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. diff --git a/test/minimum_eigensolvers/test_numpy_minimum_eigensolver.py b/test/minimum_eigensolvers/test_numpy_minimum_eigensolver.py new file mode 100644 index 000000000..9e23fe837 --- /dev/null +++ b/test/minimum_eigensolvers/test_numpy_minimum_eigensolver.py @@ -0,0 +1,238 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 import QiskitAlgorithmsTestCase + +import numpy as np +from ddt import ddt, data + +from qiskit.quantum_info import Operator, SparsePauliOp + +from qiskit_optimization.minimum_eigensolvers import NumPyMinimumEigensolver + +H2_SPARSE_PAULI = SparsePauliOp( + ["II", "ZI", "IZ", "ZZ", "XX"], + coeffs=[ + -1.052373245772859, + 0.39793742484318045, + -0.39793742484318045, + -0.01128010425623538, + 0.18093119978423156, + ], +) + +H2_OP = Operator(H2_SPARSE_PAULI.to_matrix()) + + +@ddt +class TestNumPyMinimumEigensolver(QiskitAlgorithmsTestCase): + """Test NumPy minimum eigensolver""" + + def setUp(self): + super().setUp() + aux_op1 = Operator(SparsePauliOp(["II"], coeffs=[2.0]).to_matrix()) + aux_op2 = SparsePauliOp(["II", "ZZ", "YY", "XX"], coeffs=[0.5, 0.5, 0.5, -0.5]) + self.aux_ops_list = [aux_op1, aux_op2] + self.aux_ops_dict = {"aux_op1": aux_op1, "aux_op2": aux_op2} + + @data(H2_SPARSE_PAULI, H2_OP) + def test_cme(self, op): + """Basic test""" + algo = NumPyMinimumEigensolver() + result = algo.compute_minimum_eigenvalue(operator=op, aux_operators=self.aux_ops_list) + self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j) + self.assertEqual(len(result.aux_operators_evaluated), 2) + self.assertAlmostEqual(result.aux_operators_evaluated[0][0], 2) + self.assertAlmostEqual(result.aux_operators_evaluated[1][0], 0) + + @data(H2_SPARSE_PAULI, H2_OP) + def test_cme_reuse(self, op): + """Test reuse""" + algo = NumPyMinimumEigensolver() + + with self.subTest("Test with no operator or aux_operators, give via compute method"): + result = algo.compute_minimum_eigenvalue(operator=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=op, aux_operators=self.aux_ops_list) + self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j) + self.assertEqual(len(result.aux_operators_evaluated), 2) + self.assertAlmostEqual(result.aux_operators_evaluated[0][0], 2) + self.assertAlmostEqual(result.aux_operators_evaluated[1][0], 0) + + with self.subTest("Test with aux_operators removed"): + result = algo.compute_minimum_eigenvalue(operator=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=op, aux_operators=self.aux_ops_list) + self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j) + self.assertEqual(len(result.aux_operators_evaluated), 2) + self.assertAlmostEqual(result.aux_operators_evaluated[0][0], 2) + self.assertAlmostEqual(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) + + @data(H2_SPARSE_PAULI, H2_OP) + def test_cme_filter(self, op): + """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=op, aux_operators=self.aux_ops_list) + self.assertAlmostEqual(result.eigenvalue, -0.22491125 + 0j) + self.assertEqual(len(result.aux_operators_evaluated), 2) + self.assertAlmostEqual(result.aux_operators_evaluated[0][0], 2) + self.assertAlmostEqual(result.aux_operators_evaluated[1][0], 0) + + @data(H2_SPARSE_PAULI, H2_OP) + def test_cme_filter_empty(self, op): + """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=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 = SparsePauliOp([op], coeffs=1.0) + result = algo.compute_minimum_eigenvalue(operator=operator) + self.assertAlmostEqual(result.eigenvalue, -1) + + @data(H2_SPARSE_PAULI, H2_OP) + def test_cme_aux_ops_dict(self, op): + """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=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=op, aux_operators=self.aux_ops_dict) + self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j) + self.assertEqual(len(result.aux_operators_evaluated), 2) + self.assertAlmostEqual(result.aux_operators_evaluated["aux_op1"][0], 2) + self.assertAlmostEqual(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=op, aux_operators=extra_ops) + self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j) + self.assertEqual(len(result.aux_operators_evaluated), 3) + self.assertAlmostEqual(result.aux_operators_evaluated["aux_op1"][0], 2) + self.assertAlmostEqual(result.aux_operators_evaluated["aux_op2"][0], 0) + self.assertEqual(result.aux_operators_evaluated["zero_op"], (0.0, {"variance": 0})) + + @data(H2_SPARSE_PAULI, H2_OP) + def test_aux_operators_list(self, op): + """Test list-based aux_operators.""" + algo = NumPyMinimumEigensolver() + + with self.subTest("Test with two auxiliary operators."): + result = algo.compute_minimum_eigenvalue(operator=op, aux_operators=self.aux_ops_list) + 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].pop("variance"), 0.0) + self.assertAlmostEqual(result.aux_operators_evaluated[1][1].pop("variance"), 0.0) + + with self.subTest("Test with additional zero and None operators."): + extra_ops = [*self.aux_ops_list, None, 0] + result = algo.compute_minimum_eigenvalue(operator=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].pop("variance"), 0.0) + self.assertAlmostEqual(result.aux_operators_evaluated[1][1].pop("variance"), 0.0) + self.assertEqual(result.aux_operators_evaluated[3][1].pop("variance"), 0.0) + + @data(H2_SPARSE_PAULI, H2_OP) + def test_aux_operators_dict(self, op): + """Test dict-based aux_operators.""" + algo = NumPyMinimumEigensolver() + + with self.subTest("Test with two auxiliary operators."): + result = algo.compute_minimum_eigenvalue(operator=op, aux_operators=self.aux_ops_dict) + 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].pop("variance"), 0.0 + ) + self.assertAlmostEqual( + result.aux_operators_evaluated["aux_op2"][1].pop("variance"), 0.0 + ) + + with self.subTest("Test with additional zero and None operators."): + extra_ops = {**self.aux_ops_dict, "None_operator": None, "zero_operator": 0} + result = algo.compute_minimum_eigenvalue(operator=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].pop("variance"), 0.0 + ) + self.assertAlmostEqual( + result.aux_operators_evaluated["aux_op2"][1].pop("variance"), 0.0 + ) + self.assertAlmostEqual( + result.aux_operators_evaluated["zero_operator"][1].pop("variance"), 0.0 + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/minimum_eigensolvers/test_qaoa.py b/test/minimum_eigensolvers/test_qaoa.py new file mode 100644 index 000000000..ca6caa985 --- /dev/null +++ b/test/minimum_eigensolvers/test_qaoa.py @@ -0,0 +1,304 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 QAOA algorithm.""" + +import unittest +from functools import partial +from test import QiskitAlgorithmsTestCase + +import numpy as np +import rustworkx as rx +from ddt import ddt, idata, unpack +from scipy.optimize import minimize as scipy_minimize + +from qiskit import QuantumCircuit +from qiskit.circuit import Parameter +from qiskit.primitives import Sampler +from qiskit.quantum_info import Pauli, SparsePauliOp +from qiskit.result import QuasiDistribution + +from qiskit_optimization.minimum_eigensolvers import QAOA +from qiskit_optimization.optimizers import COBYLA, NELDER_MEAD +from qiskit_optimization.utils import algorithm_globals + +W1 = np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]]) +P1 = 1 +M1 = SparsePauliOp.from_list( + [ + ("IIIX", 1), + ("IIXI", 1), + ("IXII", 1), + ("XIII", 1), + ] +) +S1 = {"0101", "1010"} + + +W2 = np.array( + [ + [0.0, 8.0, -9.0, 0.0], + [8.0, 0.0, 7.0, 9.0], + [-9.0, 7.0, 0.0, -8.0], + [0.0, 9.0, -8.0, 0.0], + ] +) +P2 = 1 +M2 = None +S2 = {"1011", "0100"} + +CUSTOM_SUPERPOSITION = [1 / np.sqrt(15)] * 15 + [0] + + +@ddt +class TestQAOA(QiskitAlgorithmsTestCase): + """Test QAOA with MaxCut.""" + + def setUp(self): + super().setUp() + self.seed = 10598 + algorithm_globals.random_seed = self.seed + self.sampler = Sampler() + + @idata( + [ + [W1, P1, M1, S1], + [W2, P2, M2, S2], + ] + ) + @unpack + def test_qaoa(self, w, reps, mixer, solutions): + """QAOA test""" + self.log.debug("Testing %s-step QAOA with MaxCut on graph\n%s", reps, w) + + qubit_op, _ = self._get_operator(w) + + qaoa = QAOA(self.sampler, COBYLA(), reps=reps, mixer=mixer) + result = qaoa.compute_minimum_eigenvalue(operator=qubit_op) + + x = self._sample_most_likely(result.eigenstate) + graph_solution = self._get_graph_solution(x) + self.assertIn(graph_solution, solutions) + + @idata( + [ + [W1, P1, S1], + [W2, P2, S2], + ] + ) + @unpack + def test_qaoa_qc_mixer(self, w, prob, solutions): + """QAOA test with a mixer as a parameterized circuit""" + self.log.debug( + "Testing %s-step QAOA with MaxCut on graph with a mixer as a parameterized circuit\n%s", + prob, + w, + ) + + optimizer = COBYLA() + qubit_op, _ = self._get_operator(w) + + num_qubits = qubit_op.num_qubits + mixer = QuantumCircuit(num_qubits) + theta = Parameter("θ") + mixer.rx(theta, range(num_qubits)) + + qaoa = QAOA(self.sampler, optimizer, reps=prob, mixer=mixer) + result = qaoa.compute_minimum_eigenvalue(operator=qubit_op) + x = self._sample_most_likely(result.eigenstate) + graph_solution = self._get_graph_solution(x) + self.assertIn(graph_solution, solutions) + + def test_qaoa_qc_mixer_many_parameters(self): + """QAOA test with a mixer as a parameterized circuit with the num of parameters > 1.""" + optimizer = COBYLA() + qubit_op, _ = self._get_operator(W1) + + num_qubits = qubit_op.num_qubits + mixer = QuantumCircuit(num_qubits) + for i in range(num_qubits): + theta = Parameter("θ" + str(i)) + mixer.rx(theta, range(num_qubits)) + + qaoa = QAOA(self.sampler, optimizer, reps=2, mixer=mixer) + result = qaoa.compute_minimum_eigenvalue(operator=qubit_op) + x = self._sample_most_likely(result.eigenstate) + self.log.debug(x) + graph_solution = self._get_graph_solution(x) + self.assertIn(graph_solution, S1) + + def test_qaoa_qc_mixer_no_parameters(self): + """QAOA test with a mixer as a parameterized circuit with zero parameters.""" + qubit_op, _ = self._get_operator(W1) + + num_qubits = qubit_op.num_qubits + mixer = QuantumCircuit(num_qubits) + # just arbitrary circuit + mixer.rx(np.pi / 2, range(num_qubits)) + + qaoa = QAOA(self.sampler, COBYLA(), reps=1, mixer=mixer) + result = qaoa.compute_minimum_eigenvalue(operator=qubit_op) + # we just assert that we get a result, it is not meaningful. + self.assertIsNotNone(result.eigenstate) + + def test_change_operator_size(self): + """QAOA change operator size test""" + qubit_op, _ = self._get_operator( + np.array([[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]]) + ) + qaoa = QAOA(self.sampler, COBYLA(), reps=1) + result = qaoa.compute_minimum_eigenvalue(operator=qubit_op) + x = self._sample_most_likely(result.eigenstate) + graph_solution = self._get_graph_solution(x) + with self.subTest(msg="QAOA 4x4"): + self.assertIn(graph_solution, {"0101", "1010"}) + + qubit_op, _ = self._get_operator( + np.array( + [ + [0, 1, 0, 1, 0, 1], + [1, 0, 1, 0, 1, 0], + [0, 1, 0, 1, 0, 1], + [1, 0, 1, 0, 1, 0], + [0, 1, 0, 1, 0, 1], + [1, 0, 1, 0, 1, 0], + ] + ) + ) + result = qaoa.compute_minimum_eigenvalue(operator=qubit_op) + x = self._sample_most_likely(result.eigenstate) + graph_solution = self._get_graph_solution(x) + with self.subTest(msg="QAOA 6x6"): + self.assertIn(graph_solution, {"010101", "101010"}) + + @idata([[W2, S2, None], [W2, S2, [0.0, 0.0]], [W2, S2, [1.0, 0.8]]]) + @unpack + def test_qaoa_initial_point(self, w, solutions, init_pt): + """Check first parameter value used is initial point as expected""" + qubit_op, _ = self._get_operator(w) + + first_pt = [] + + def cb_callback(eval_count, parameters, mean, metadata): + nonlocal first_pt + if eval_count == 1: + first_pt = list(parameters) + + qaoa = QAOA( + self.sampler, + COBYLA(), + initial_point=init_pt, + callback=cb_callback, + ) + result = qaoa.compute_minimum_eigenvalue(operator=qubit_op) + + x = self._sample_most_likely(result.eigenstate) + graph_solution = self._get_graph_solution(x) + + with self.subTest("Initial Point"): + # If None the preferred random initial point of QAOA variational form + if init_pt is None: + self.assertLess(result.eigenvalue, -0.97) + else: + self.assertListEqual(init_pt, first_pt) + + with self.subTest("Solution"): + self.assertIn(graph_solution, solutions) + + def test_qaoa_random_initial_point(self): + """QAOA random initial point""" + # the function undirected_gnp_random_graph() does exist in + # rustworkx packagebut the linter can't see it + w = rx.adjacency_matrix( + rx.undirected_gnp_random_graph( # pylint: disable=no-member + 5, 0.5, seed=algorithm_globals.random_seed + ) + ) + qubit_op, _ = self._get_operator(w) + qaoa = QAOA(self.sampler, NELDER_MEAD(disp=True), reps=2) + result = qaoa.compute_minimum_eigenvalue(operator=qubit_op) + + self.assertLess(result.eigenvalue, -0.97) + + def test_optimizer_scipy_callable(self): + """Test passing a SciPy optimizer directly as callable.""" + w = rx.adjacency_matrix( + rx.undirected_gnp_random_graph( # pylint: disable=no-member + 5, 0.5, seed=algorithm_globals.random_seed + ) + ) + qubit_op, _ = self._get_operator(w) + qaoa = QAOA( + self.sampler, + partial(scipy_minimize, method="Nelder-Mead", options={"maxiter": 2}), + ) + result = qaoa.compute_minimum_eigenvalue(qubit_op) + self.assertEqual(result.cost_function_evals, 5) + + def _get_operator(self, weight_matrix): + """Generate Hamiltonian for the max-cut problem of a graph. + + Args: + weight_matrix (numpy.ndarray) : adjacency matrix. + + Returns: + PauliSumOp: operator for the Hamiltonian + float: a constant shift for the obj function. + + """ + num_nodes = weight_matrix.shape[0] + pauli_list = [] + shift = 0 + for i in range(num_nodes): + for j in range(i): + if weight_matrix[i, j] != 0: + x_p = np.zeros(num_nodes, dtype=bool) + z_p = np.zeros(num_nodes, dtype=bool) + z_p[i] = True + z_p[j] = True + pauli_list.append([0.5 * weight_matrix[i, j], Pauli((z_p, x_p))]) + shift -= 0.5 * weight_matrix[i, j] + lst = [(pauli[1].to_label(), pauli[0]) for pauli in pauli_list] + return SparsePauliOp.from_list(lst), shift + + def _get_graph_solution(self, x: np.ndarray) -> str: + """Get graph solution from binary string. + + Args: + x : binary string as numpy array. + + Returns: + a graph solution as string. + """ + + return "".join([str(int(i)) for i in 1 - x]) + + def _sample_most_likely(self, state_vector: QuasiDistribution) -> np.ndarray: + """Compute the most likely binary string from state vector. + Args: + state_vector: Quasi-distribution. + + Returns: + Binary string as numpy.ndarray of ints. + """ + values = list(state_vector.values()) + n = int(np.log2(len(values))) + k = np.argmax(np.abs(values)) + x = np.zeros(n) + for i in range(n): + x[i] = k % 2 + k >>= 1 + return x + + +if __name__ == "__main__": + unittest.main() diff --git a/test/minimum_eigensolvers/test_sampling_vqe.py b/test/minimum_eigensolvers/test_sampling_vqe.py new file mode 100644 index 000000000..c789c4b68 --- /dev/null +++ b/test/minimum_eigensolvers/test_sampling_vqe.py @@ -0,0 +1,101 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2018, 2024. +# +# 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 Sampler VQE.""" + + +import unittest +from functools import partial +from test import QiskitAlgorithmsTestCase + +import numpy as np +from ddt import ddt +from scipy.optimize import minimize as scipy_minimize + +from qiskit.circuit.library import RealAmplitudes +from qiskit.primitives import Sampler +from qiskit.quantum_info import Pauli + +from qiskit_optimization.minimum_eigensolvers import SamplingVQE +from qiskit_optimization.optimizers import OptimizerResult +from qiskit_optimization.utils import algorithm_globals + + +# pylint: disable=invalid-name +def _mock_optimizer(fun, x0, jac=None, bounds=None, inputs=None): + """A mock of a callable that can be used as minimizer in the VQE. + + If ``inputs`` is given as a dictionary, stores the inputs in that dictionary. + """ + 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 TestSamplerVQE(QiskitAlgorithmsTestCase): + """Test VQE""" + + def setUp(self): + super().setUp() + self.optimal_value = -1.38 + self.optimal_bitstring = "10" + algorithm_globals.random_seed = 42 + + def test_optimizer_scipy_callable(self): + """Test passing a SciPy optimizer directly as callable.""" + vqe = SamplingVQE( + Sampler(), + RealAmplitudes(), + partial(scipy_minimize, method="COBYLA", options={"maxiter": 2}), + ) + result = vqe.compute_minimum_eigenvalue(Pauli("Z")) + self.assertEqual(result.cost_function_evals, 2) + + def test_optimizer_callable(self): + """Test passing a optimizer directly as callable.""" + ansatz = RealAmplitudes(1, reps=1) + vqe = SamplingVQE(Sampler(), ansatz, _mock_optimizer) + result = vqe.compute_minimum_eigenvalue(Pauli("Z")) + self.assertTrue(np.all(result.optimal_point == np.zeros(ansatz.num_parameters))) + + def test_aggregation(self): + """Test the aggregation works.""" + + # test a custom aggregation that just uses the best measurement + def best_measurement(measurements): + res = min(measurements, key=lambda meas: meas[1])[1] + return res + + # test CVaR with alpha of 0.4 (i.e. 40% of the best measurements) + alpha = 0.4 + + ansatz = RealAmplitudes(1, reps=0) + ansatz.h(0) + + for aggregation in [alpha, best_measurement]: + with self.subTest(aggregation=aggregation): + vqe = SamplingVQE(Sampler(), ansatz, _mock_optimizer, aggregation=best_measurement) + result = vqe.compute_minimum_eigenvalue(Pauli("Z")) + + # evaluation at x0=0 samples -1 and 1 with 50% probability, and our aggregation + # takes the smallest value + self.assertAlmostEqual(result.optimal_value, -1) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/optimizers/__init__.py b/test/optimizers/__init__.py new file mode 100644 index 000000000..c9d7aa94a --- /dev/null +++ b/test/optimizers/__init__.py @@ -0,0 +1,13 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2021, 2024. +# +# 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. + +"""Qiskit's algorithm optimizer tests.""" diff --git a/test/optimizers/test_optimizers.py b/test/optimizers/test_optimizers.py new file mode 100644 index 000000000..1d45e9ece --- /dev/null +++ b/test/optimizers/test_optimizers.py @@ -0,0 +1,382 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2018, 2024. +# +# 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 Optimizers""" + +import unittest +from test import QiskitAlgorithmsTestCase + +from typing import Optional, List, Tuple +from ddt import ddt, data, unpack +import numpy as np +from scipy.optimize import rosen, rosen_der + +from qiskit_optimization.optimizers import ( + COBYLA, + NELDER_MEAD, + Optimizer, + SPSA, + SciPyOptimizer, + OptimizerSupportLevel, +) +from qiskit_optimization.utils import algorithm_globals + + +@ddt +class TestOptimizers(QiskitAlgorithmsTestCase): + """Test Optimizers""" + + def setUp(self): + super().setUp() + algorithm_globals.random_seed = 52 + self.optimizer = SPSA() + self.optimizer_2 = SPSA() + self.optimizer_2.set_options(tolerance=1e-6, maxiter=100, method="SPSA") + + def run_optimizer( + self, + optimizer: Optimizer, + max_nfev: int, + grad: bool = False, + bounds: Optional[List[Tuple[float, float]]] = None, + ): + """Test the optimizer. + + Args: + optimizer: The optimizer instance to test. + max_nfev: The maximal allowed number of function evaluations. + grad: Whether to pass the gradient function as input. + bounds: Optimizer bounds. + """ + x_0 = np.asarray([1.3, 0.7, 0.8, 1.9, 1.2]) + jac = rosen_der if grad else None + + res = optimizer.minimize(rosen, x_0, jac, bounds) + x_opt = res.x + nfev = res.nfev + + np.testing.assert_array_almost_equal(x_opt, [1.0] * len(x_0), decimal=2) + self.assertLessEqual(nfev, max_nfev) + + def test_cobyla(self): + """cobyla test""" + optimizer = COBYLA(maxiter=100000, tol=1e-06) + self.run_optimizer(optimizer, max_nfev=100000) + + def test_nelder_mead(self): + """nelder mead test""" + optimizer = NELDER_MEAD(maxfev=10000, tol=1e-06) + self.run_optimizer(optimizer, max_nfev=10000) + + @unittest.skip("Skipping SPSA as it does not do well on non-convex rozen") + def test_spsa(self): + """spsa test""" + optimizer = SPSA(maxiter=10000) + self.run_optimizer(optimizer, max_nfev=100000) + + def test_scipy_optimizer(self): + """scipy_optimizer test""" + optimizer = SciPyOptimizer("BFGS", options={"maxiter": 1000}) + self.run_optimizer(optimizer, max_nfev=10000) + + def test_scipy_optimizer_callback(self): + """scipy_optimizer callback test""" + values = [] + + def callback(x): + values.append(x) + + optimizer = SciPyOptimizer("BFGS", options={"maxiter": 1000}, callback=callback) + self.run_optimizer(optimizer, max_nfev=10000) + self.assertTrue(values) # Check the list is nonempty. + + def test_scipy_optimizer_parse_bounds(self): + """ + Test the parsing of bounds in SciPyOptimizer.minimize method. Verifies that the bounds are + correctly parsed and set within the optimizer object. + + Raises: + AssertionError: If any of the assertions fail. + AssertionError: If a TypeError is raised unexpectedly while parsing bounds. + + """ + try: + # Initialize SciPyOptimizer instance with SLSQP method + optimizer = SciPyOptimizer("SLSQP") + + # Call minimize method with a simple lambda function and bounds + optimizer.minimize(lambda x: -x, 1.0, bounds=[(0.0, 1.0)]) + + # Assert that "bounds" is not present in optimizer options and kwargs + self.assertFalse("bounds" in optimizer._options) + self.assertFalse("bounds" in optimizer._kwargs) + + except TypeError: + # This would give: https://github.com/qiskit-community/qiskit-machine-learning/issues/570 + self.fail( + "TypeError was raised unexpectedly when parsing bounds in SciPyOptimizer.minimize(...)." + ) + + # Finally, expect exceptions if bounds are parsed incorrectly, i.e. differently than as in Scipy + with self.assertRaises(RuntimeError): + _ = SciPyOptimizer("SLSQP", bounds=[(0.0, 1.0)]) + with self.assertRaises(RuntimeError): + _ = SciPyOptimizer("SLSQP", options={"bounds": [(0.0, 1.0)]}) + + def test_gradient_num_diff(self): + """Test the gradient_num_diff function.""" + + # Define a simple quadratic function and its gradient + def func(x): + return (x[0] - 2) ** 2 + (x[1] - 3) ** 2 + + def expected_gradient(x): + return np.array([2 * (x[0] - 2), 2 * (x[1] - 3)]) + + # Set the point around which we compute the gradient + x_center = np.array([1.0, 1.0]) + epsilon = 1e-5 # Small perturbation for numerical differentiation + + # Compute the numerical gradient using the optimizer method + numerical_gradient = self.optimizer.gradient_num_diff(x_center, func, epsilon) + + # Compute the expected gradient + expected_grad = expected_gradient(x_center) + + # Assert that the computed gradient is close to the expected gradient + np.testing.assert_allclose(numerical_gradient, expected_grad, rtol=1e-5, atol=1e-8) + + def test_set_options(self): + """Test the set_options method.""" + + # Define some options to set + options = {"max_iter": 100, "tolerance": 1e-6, "verbose": True} + + # Set options using the set_options method + self.optimizer.set_options(**options) + + # Assert that the options dictionary is updated correctly + for key, value in options.items(): + self.assertIn(key, self.optimizer._options) + self.assertEqual(self.optimizer._options[key], value) + + # Test updating an existing option + self.optimizer.set_options(max_iter=200) + self.assertEqual(self.optimizer._options["max_iter"], 200) + + def test_wrap_function(self): + """Test the wrap_function method.""" + + # Define a simple function to test + def simple_function(x, y): + return x + y + + # Wrap the function, injecting the argument (5,) + wrapped_function = self.optimizer.wrap_function(simple_function, (5,)) + + # Call the wrapped function with a single argument + result = wrapped_function(10) # Should compute 10 + 5 + + # Assert that the result is as expected + self.assertEqual(result, 15) + + def test_wrap_function_with_multiple_args(self): + """Test wrap_function with multiple injected args.""" + + # Define a simple function to test + def multiply_function(a, b, c): + return a * b * c + + # Wrap the function, injecting the arguments (2, 3) + wrapped_function = self.optimizer.wrap_function(multiply_function, (2, 3)) + + # Call the wrapped function with a single argument + result = wrapped_function(4) # Should compute 4 * 2 * 3 + + # Assert that the result is as expected + self.assertEqual(result, 24) + + def test_setting(self): + """Test the setting property.""" + + actual_output = self.optimizer.setting + + # Check if key parts are present in the settings output + self.assertIn("Optimizer: SPSA", actual_output) + self.assertIn("max_evals_grouped: None", actual_output) + + # Optional: check for specific support levels if required + self.assertIn("gradient_support_level", actual_output) + self.assertIn("bounds_support_level", actual_output) + self.assertIn("initial_point_support_level", actual_output) + + def test_gradient_support_level(self): + """Test for gradient support level property""" + self.assertEqual(self.optimizer.gradient_support_level, OptimizerSupportLevel.ignored) + + def test_is_gradient_ignored(self): + """Test for is_gradient_ignored property""" + self.assertTrue(self.optimizer.is_gradient_ignored) + + def test_is_gradient_supported(self): + """Test for is_gradient_supported property""" + self.assertTrue(self.optimizer.is_gradient_supported) + + def test_is_gradient_required(self): + """Test for is_gradient_required property""" + self.assertFalse(self.optimizer.is_gradient_required) + + def test_bounds_support_level(self): + """Test for bounds support level property""" + self.assertNotEqual(self.optimizer.bounds_support_level, OptimizerSupportLevel.supported) + + def test_is_bounds_ignored(self): + """Test for is_bounds_ignored property""" + self.assertTrue(self.optimizer.is_bounds_ignored) + + def test_is_bounds_supported(self): + """Test for is_bounds_supported property""" + self.assertTrue(self.optimizer.is_bounds_supported) + + def test_is_bounds_required(self): + """Test for is_bounds_required property""" + self.assertFalse(self.optimizer.is_bounds_required) + + def test_initial_point_support_level(self): + """Test for initial point support level property""" + self.assertEqual(self.optimizer.initial_point_support_level, OptimizerSupportLevel.required) + + def test_is_initial_point_ignored(self): + """Test for is_initial_point_ignored property""" + self.assertFalse(self.optimizer.is_initial_point_ignored) + + def test_is_initial_point_supported(self): + """Test for is_initial_point_supported property""" + self.assertTrue(self.optimizer.is_initial_point_supported) + + def test_is_initial_point_required(self): + """Test for is_initial_point_required property""" + self.assertTrue(self.optimizer.is_initial_point_required) + + def test_set_max_evals_grouped(self): + """Test for set_max_evals_grouped method""" + self.optimizer.set_max_evals_grouped(10) + self.assertEqual(self.optimizer._max_evals_grouped, 10) + + +@ddt +class TestOptimizerSerialization(QiskitAlgorithmsTestCase): + """Tests concerning the serialization of optimizers.""" + + @data( + ("COBYLA", {"maxiter": 10}), + ("NELDER_MEAD", {"maxiter": 0}), + ("dogleg", {"maxiter": 100}), + ("trust-constr", {"maxiter": 10}), + ("trust-ncg", {"maxiter": 100}), + ("trust-exact", {"maxiter": 120}), + ("trust-krylov", {"maxiter": 150}), + ) + @unpack + def test_scipy(self, method, options): + """Test the SciPyOptimizer is serializable.""" + + optimizer = SciPyOptimizer(method, options=options) + serialized = optimizer.settings + from_dict = SciPyOptimizer(**serialized) + + self.assertEqual(from_dict._method, method.lower()) + self.assertEqual(from_dict._options, options) + + def test_independent_reconstruction(self): + """Test the SciPyOptimizers don't reset all settings upon creating a new instance. + + COBYLA is used as representative example here.""" + + kwargs = {"coffee": "without sugar"} + options = {"tea": "with milk"} + optimizer = COBYLA(maxiter=1, options=options, **kwargs) + serialized = optimizer.settings + from_dict = COBYLA(**serialized) + + with self.subTest(msg="test attributes"): + self.assertEqual(from_dict.settings["maxiter"], 1) + + with self.subTest(msg="test options"): + # options should only contain values that are *not* already in the initializer + # (e.g. should not contain maxiter) + self.assertEqual(from_dict.settings["options"], {"tea": "with milk"}) + + with self.subTest(msg="test kwargs"): + self.assertEqual(from_dict.settings["coffee"], "without sugar") + + with self.subTest(msg="option ids differ"): + self.assertNotEqual(id(serialized["options"]), id(from_dict.settings["options"])) + + def test_spsa(self): + """Test SPSA optimizer is serializable.""" + options = { + "maxiter": 100, + "blocking": True, + "allowed_increase": 0.1, + "second_order": True, + "learning_rate": 0.02, + "perturbation": 0.05, + "regularization": 0.1, + "resamplings": 2, + "perturbation_dims": 5, + "trust_region": False, + "initial_hessian": None, + "lse_solver": None, + "hessian_delay": 0, + "callback": None, + "termination_checker": None, + } + spsa = SPSA(**options) + + self.assertDictEqual(spsa.settings, options) + + def test_spsa_custom_iterators(self): + """Test serialization works with custom iterators for learning rate and perturbation.""" + rate = 0.99 + + def powerlaw(): + n = 0 + while True: + yield rate**n + n += 1 + + def steps(): + n = 1 + divide_after = 20 + epsilon = 0.5 + while True: + yield epsilon + n += 1 + if n % divide_after == 0: + epsilon /= 2 + + learning_rate = powerlaw() + expected_learning_rate = np.array([next(learning_rate) for _ in range(200)]) + + perturbation = steps() + expected_perturbation = np.array([next(perturbation) for _ in range(200)]) + + spsa = SPSA(maxiter=200, learning_rate=powerlaw, perturbation=steps) + settings = spsa.settings + + self.assertTrue(np.allclose(settings["learning_rate"], expected_learning_rate)) + self.assertTrue(np.allclose(settings["perturbation"], expected_perturbation)) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/optimizers/test_spsa.py b/test/optimizers/test_spsa.py new file mode 100644 index 000000000..f4b1f01e8 --- /dev/null +++ b/test/optimizers/test_spsa.py @@ -0,0 +1,183 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2021, 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Tests for the SPSA optimizer.""" + +from test import QiskitAlgorithmsTestCase +from ddt import ddt, data + +import numpy as np + +from qiskit.circuit.library import PauliTwoDesign +from qiskit.quantum_info import SparsePauliOp, Statevector + +from qiskit_optimization.optimizers import SPSA +from qiskit_optimization.utils import algorithm_globals + + +@ddt +class TestSPSA(QiskitAlgorithmsTestCase): + """Tests for the SPSA optimizer.""" + + def setUp(self): + super().setUp() + np.random.seed(12) + algorithm_globals.random_seed = 12 + + # @slow_test + def test_pauli_two_design(self): + """Test SPSA on the Pauli two-design example.""" + circuit = PauliTwoDesign(3, reps=1, seed=1) + parameters = list(circuit.parameters) + obs = SparsePauliOp("ZZI") # Z^Z^I + + initial_point = np.array( + [0.82311034, 0.02611798, 0.21077064, 0.61842177, 0.09828447, 0.62013131] + ) + + def objective(x): + bound_circ = circuit.assign_parameters(dict(zip(parameters, x))) + return Statevector(bound_circ).expectation_value(obs).real + + settings = {"maxiter": 100, "blocking": True, "allowed_increase": 0} + + expected_nfev = settings["maxiter"] * 3 + 1 + + spsa = SPSA(**settings) + + result = spsa.minimize(objective, x0=initial_point) + + with self.subTest("check final accuracy"): + self.assertLess(result.fun, -0.95) # final loss + + with self.subTest("check number of function calls"): + self.assertEqual(result.nfev, expected_nfev) # function evaluations + + def test_recalibrate_at_optimize(self): + """Test SPSA calibrates anew upon each optimization run, if no auto-calibration is set.""" + + def objective(x): + return -(x**2) + + spsa = SPSA(maxiter=1) + _ = spsa.minimize(objective, x0=np.array([0.5])) + + self.assertIsNone(spsa.learning_rate) + self.assertIsNone(spsa.perturbation) + + def test_learning_rate_perturbation_as_iterators(self): + """Test the learning rate and perturbation can be callables returning iterators.""" + + def get_learning_rate(): + def learning_rate(): + x = 0.99 + while True: + x *= x + yield x + + return learning_rate + + def get_perturbation(): + def perturbation(): + x = 0.99 + while True: + x *= x + yield max(x, 0.01) + + return perturbation + + def objective(x): + return (np.linalg.norm(x) - 2) ** 2 + + spsa = SPSA(learning_rate=get_learning_rate(), perturbation=get_perturbation()) + result = spsa.minimize(objective, np.array([0.5, 0.5])) + + self.assertAlmostEqual(np.linalg.norm(result.x), 2, places=2) + + def test_learning_rate_perturbation_as_arrays(self): + """Test the learning rate and perturbation can be arrays.""" + + learning_rate = np.linspace(1, 0, num=100, endpoint=False) ** 2 + perturbation = 0.01 * np.ones(100) + + def objective(x): + return (np.linalg.norm(x) - 2) ** 2 + + spsa = SPSA(learning_rate=learning_rate, perturbation=perturbation) + result = spsa.minimize(objective, x0=np.array([0.5, 0.5])) + + self.assertAlmostEqual(np.linalg.norm(result.x), 2, places=2) + + def test_termination_checker(self): + """Test the termination_callback""" + + def objective(x): + return np.linalg.norm(x) + np.random.rand(1) + + class TerminationChecker: + """Example termination checker""" + + def __init__(self): + self.values = [] + + def __call__( # pylint: disable=too-many-positional-arguments + self, nfev, point, fvalue, stepsize, accepted + ) -> bool: + self.values.append(fvalue) + + if len(self.values) > 10: + return True + return False + + maxiter = 400 + spsa = SPSA(maxiter=maxiter, termination_checker=TerminationChecker()) + result = spsa.minimize(objective, x0=[0.5, 0.5]) + + self.assertLess(result.nit, maxiter) + + def test_callback(self): + """Test using the callback.""" + + def objective(x): + return (np.linalg.norm(x) - 2) ** 2 + + history = {"nfevs": [], "points": [], "fvals": [], "updates": [], "accepted": []} + + def callback(nfev, point, fval, update, accepted): + history["nfevs"].append(nfev) + history["points"].append(point) + history["fvals"].append(fval) + history["updates"].append(update) + history["accepted"].append(accepted) + + maxiter = 10 + spsa = SPSA(maxiter=maxiter, learning_rate=0.01, perturbation=0.01, callback=callback) + _ = spsa.minimize(objective, x0=np.array([0.5, 0.5])) + + expected_types = [int, np.ndarray, float, float, bool] + for i, (key, values) in enumerate(history.items()): + self.assertTrue(all(isinstance(value, expected_types[i]) for value in values)) + self.assertEqual(len(history[key]), maxiter) + + @data(1, 2, 3, 4) + def test_estimate_stddev(self, max_evals_grouped): + """Test the estimate_stddev + See https://github.com/Qiskit/qiskit-nature/issues/797""" + + def objective(x): + if len(x.shape) == 2: + return np.array([sum(x_i) for x_i in x]) + return sum(x) + + point = np.ones(5) + result = SPSA.estimate_stddev(objective, point, avg=10, max_evals_grouped=max_evals_grouped) + self.assertAlmostEqual(result, 0) diff --git a/test/test_algorithm_result.py b/test/test_algorithm_result.py new file mode 100644 index 000000000..823007b5c --- /dev/null +++ b/test/test_algorithm_result.py @@ -0,0 +1,70 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2019, 2024. +# +# 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 AlgorithmResult""" + +import unittest + +from test import QiskitAlgorithmsTestCase + +from qiskit_optimization.algorithm_result import AlgorithmResult + + +class TestAlgorithmResult(AlgorithmResult): + """Concrete subclass for testing purposes""" + + def __init__(self, data): + self.data = data + self.name = "Test Result" + + +class TestAlgorithmResultMethods(QiskitAlgorithmsTestCase): + """AlgorithmResult tests""" + + def setUp(self): + """Setting up initial test objects""" + super().setUp() + self.result1 = TestAlgorithmResult({"value1": 10, "value2": 20}) + self.result2 = TestAlgorithmResult({"value1": 100, "value2": 200}) + self.result3 = TestAlgorithmResult({"value3": 300}) + + def test_str_method(self): + """Test the __str__ method""" + expected_str = "{'data': {'value1': 10, 'value2': 20}, 'name': 'Test Result'}" + self.assertEqual(str(self.result1), expected_str) + + def test_combine_with_another_result(self): + """Test the combine method with another result""" + self.result1.combine(self.result2) + self.assertEqual(self.result1.data, {"value1": 100, "value2": 200}) + self.assertEqual(self.result1.name, "Test Result") + + def test_combine_with_new_property(self): + """Test combining with a result that has a new property""" + self.result1.combine(self.result3) + self.assertEqual(self.result1.data, {"value3": 300}) + self.assertEqual(self.result1.name, "Test Result") + + def test_combine_with_none_raises_error(self): + """Test that passing None to combine raises TypeError""" + with self.assertRaises(TypeError): + self.result1.combine(None) + + def test_combine_with_self_does_nothing(self): + """Test that combining with self doesn't change anything""" + original_data = self.result1.data.copy() + self.result1.combine(self.result1) + self.assertEqual(self.result1.data, original_data) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/test_observables_evaluator.py b/test/test_observables_evaluator.py new file mode 100644 index 000000000..c0d1564a6 --- /dev/null +++ b/test/test_observables_evaluator.py @@ -0,0 +1,133 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2019, 2024. +# +# 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 Observables Evaluators""" + +import unittest +from unittest.mock import MagicMock + +from test import QiskitAlgorithmsTestCase + +import numpy as np + +from qiskit import QuantumCircuit +from qiskit.quantum_info.operators import SparsePauliOp +from qiskit.primitives import BaseEstimator + +from qiskit_optimization.exceptions import AlgorithmError +from qiskit_optimization.observables_evaluator import ( + estimate_observables, + _handle_zero_ops, + _prepare_result, +) + + +class TestEstimateObservables(QiskitAlgorithmsTestCase): + """Observables Evaluators tests""" + + def setUp(self): + """Set up a basic quantum circuit and estimator for testing.""" + super().setUp() + self.estimator = MagicMock(spec=BaseEstimator) + self.quantum_state = QuantumCircuit(2) # Simple 2-qubit circuit + self.observable = SparsePauliOp.from_list([("Z", 1)]) + self.observable_2 = SparsePauliOp.from_list([("X", 1)]) + self.observable_3 = SparsePauliOp.from_list([("Z", 1)]) + + def test_estimate_observables_success_with_list(self): + """Test estimation with a list of observables and successful estimator.""" + self.estimator.run.return_value.result.return_value.values = np.array([1.0, 0.5]) + self.estimator.run.return_value.result.return_value.metadata = [{} for _ in range(2)] + observables = [self.observable, self.observable_2] + result = estimate_observables(self.estimator, self.quantum_state, observables) + # Verify results + expected_results = [(1.0, {}), (0.5, {})] + self.assertEqual(result, expected_results) + self.estimator.run.assert_called_once() + + def test_estimate_observables_success_with_dict(self): + """Test estimation with a dictionary of observables.""" + self.estimator.run = unittest.mock.MagicMock() + self.estimator.run.return_value.result.return_value.values = np.array([1.0, 0.5]) + self.estimator.run.return_value.result.return_value.metadata = [{} for _ in range(2)] + # Use valid BaseOperator instances for testing + observables = {"obs1": self.observable, "obs2": self.observable_2} + result = estimate_observables(self.estimator, self.quantum_state, observables) + # Verify results + expected_results = [(1.0, {}), (0.5, {})] + self.assertEqual( + result, _prepare_result(expected_results, observables) + ) # Adjust according to your expected output + # Assert that the estimator was called correctly + self.estimator.run.assert_called_once() + + def test_estimate_observables_below_threshold(self): + """Test estimation with values below threshold.""" + self.estimator.run.return_value.result.return_value.values = np.array([0.1, 0.0]) + self.estimator.run.return_value.result.return_value.metadata = [{} for _ in range(2)] + observables = [self.observable, self.observable_2] + result = estimate_observables( + self.estimator, self.quantum_state, observables, threshold=0.2 + ) + # Verify that the results are filtered below the threshold + expected_results = [(0.0, {}), (0.0, {})] # Both should be ignored (mean <= threshold) + self.assertEqual(result, expected_results) + self.estimator.run.assert_called_once() + + def test_estimate_observables_empty_observables(self): + """Test estimation with an empty list of observables.""" + observables = [] + result = estimate_observables(self.estimator, self.quantum_state, observables) + # Verify that the result is an empty list + self.assertEqual(result, []) + self.estimator.run.assert_not_called() + + def test_estimate_observables_algorithm_error(self): + """Test handling of AlgorithmError when estimator fails.""" + self.estimator.run.side_effect = Exception("Failed job") + observables = [self.observable, self.observable_2] + with self.assertRaises(AlgorithmError): + estimate_observables(self.estimator, self.quantum_state, observables) + + def test_handle_zero_ops(self): + """Test replacing zero operators with SparsePauliOp.""" + observables_list = [self.observable_3, 0, self.observable_3] + # Ensure num_qubits is accessible and valid + num_qubits = self.observable_3.num_qubits + self.assertIsInstance(num_qubits, int) + result = _handle_zero_ops(observables_list) + # Check if the zero operator was replaced correctly + zero_op = SparsePauliOp.from_list([("I" * num_qubits, 0)]) + # Validate that the zero operator was replaced + self.assertEqual(result[1], zero_op) + + def test_prepare_result_with_list(self): + """Test the _prepare_result function with a list of observables.""" + observables_results = [(1.0, {"meta1": "data1"}), (0.5, {"meta2": "data2"})] + observables = [self.observable, self.observable_2] + result = _prepare_result(observables_results, observables) + # Verify the output is as expected + expected_results = [(1.0, {"meta1": "data1"}), (0.5, {"meta2": "data2"})] + self.assertEqual(result, expected_results) + + def test_prepare_result_with_dict(self): + """Test the _prepare_result function with a dictionary of observables.""" + observables_results = [(1.0, {"meta1": "data1"}), (0.5, {"meta2": "data2"})] + observables = {"obs1": self.observable, "obs2": self.observable_2} + result = _prepare_result(observables_results, observables) + # Verify the output is as expected + expected_results = {"obs1": (1.0, {"meta1": "data1"}), "obs2": (0.5, {"meta2": "data2"})} + self.assertEqual(result, expected_results) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/test_validation.py b/test/test_validation.py new file mode 100644 index 000000000..70663ce19 --- /dev/null +++ b/test/test_validation.py @@ -0,0 +1,91 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2019, 2024. +# +# 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 Validation""" + +import unittest + +from test import QiskitAlgorithmsTestCase + +from qiskit_optimization.utils.validation import ( + validate_in_set, + validate_min, + validate_min_exclusive, + validate_max, + validate_max_exclusive, + validate_range, + validate_range_exclusive, + validate_range_exclusive_min, + validate_range_exclusive_max, +) + + +class TestValidation(QiskitAlgorithmsTestCase): + """Validation tests.""" + + def test_validate_in_set(self): + """validate in set test""" + test_value = "value1" + validate_in_set("test_value", test_value, {"value1", "value2"}) + with self.assertRaises(ValueError): + validate_in_set("test_value", test_value, {"value3", "value4"}) + + def test_validate_min(self): + """validate min test""" + test_value = 2.5 + validate_min("test_value", test_value, -1) + validate_min("test_value", test_value, 2.5) + with self.assertRaises(ValueError): + validate_min("test_value", test_value, 4) + validate_min_exclusive("test_value", test_value, -1) + with self.assertRaises(ValueError): + validate_min_exclusive("test_value", test_value, 2.5) + with self.assertRaises(ValueError): + validate_min_exclusive("test_value", test_value, 4) + + def test_validate_max(self): + """validate max test""" + test_value = 2.5 + with self.assertRaises(ValueError): + validate_max("test_value", test_value, -1) + validate_max("test_value", test_value, 2.5) + validate_max("test_value", test_value, 4) + with self.assertRaises(ValueError): + validate_max_exclusive("test_value", test_value, -1) + with self.assertRaises(ValueError): + validate_max_exclusive("test_value", test_value, 2.5) + validate_max_exclusive("test_value", test_value, 4) + + def test_validate_range(self): + """validate range test""" + test_value = 2.5 + with self.assertRaises(ValueError): + validate_range("test_value", test_value, 0, 2) + with self.assertRaises(ValueError): + validate_range("test_value", test_value, 3, 4) + validate_range("test_value", test_value, 2.5, 3) + validate_range_exclusive("test_value", test_value, 0, 3) + with self.assertRaises(ValueError): + validate_range_exclusive("test_value", test_value, 0, 2.5) + validate_range_exclusive("test_value", test_value, 2.5, 3) + validate_range_exclusive_min("test_value", test_value, 0, 3) + with self.assertRaises(ValueError): + validate_range_exclusive_min("test_value", test_value, 2.5, 3) + validate_range_exclusive_min("test_value", test_value, 0, 2.5) + validate_range_exclusive_max("test_value", test_value, 2.5, 3) + with self.assertRaises(ValueError): + validate_range_exclusive_max("test_value", test_value, 0, 2.5) + validate_range_exclusive_max("test_value", test_value, 2.5, 3) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/translators/test_prettyprint.py b/test/translators/test_prettyprint.py index 6fead2566..02232c494 100644 --- a/test/translators/test_prettyprint.py +++ b/test/translators/test_prettyprint.py @@ -1,6 +1,6 @@ # This code is part of a Qiskit project. # -# (C) Copyright IBM 2022, 2023. +# (C) Copyright IBM 2022, 2024. # # 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 @@ -112,11 +112,11 @@ def test_prettyprint(self): q_p.linear_constraint({"x": 1, "y": 2}, "<=", 1, "lin_leq") q_p.linear_constraint({"x": 1, "y": 2}, ">=", 1, "lin_geq") q_p.quadratic_constraint( - {"x": 1, "y": 1}, - {("x", "x"): 1, ("y", "z"): -1, ("z", "z"): 2}, - "==", - 1, - "quad_eq", + linear={"x": 1, "y": 1}, + quadratic={("x", "x"): 1, ("y", "z"): -1, ("z", "z"): 2}, + sense="==", + rhs=1, + name="quad_eq", ) q_p.quadratic_constraint( {"x": 1, "y": 1}, diff --git a/test/utils/__init__.py b/test/utils/__init__.py new file mode 100644 index 000000000..c870af8ba --- /dev/null +++ b/test/utils/__init__.py @@ -0,0 +1,11 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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. diff --git a/test/utils/test_set_batching.py b/test/utils/test_set_batching.py new file mode 100644 index 000000000..55e12c3bd --- /dev/null +++ b/test/utils/test_set_batching.py @@ -0,0 +1,55 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 set batching.""" + +from test import QiskitAlgorithmsTestCase + +from qiskit_optimization.optimizers import SPSA, COBYLA + +from qiskit_optimization.utils.set_batching import _set_default_batchsize + + +class TestSetBatching(QiskitAlgorithmsTestCase): + """Set Batching tests.""" + + def test_set_default_batchsize_updates(self): + """Test that the default batchsize is set when it is None.""" + # Create an instance of SPSA with _max_evals_grouped as None + optimizer = SPSA() + optimizer._max_evals_grouped = None # Directly set the private variable for testing + # Call the function + updated = _set_default_batchsize(optimizer) + # Check that the batch size was updated + self.assertTrue(updated) + self.assertEqual(optimizer._max_evals_grouped, 50) + + def test_set_default_batchsize_no_update(self): + """Test that the batchsize is not updated when it is already set.""" + # Create an instance of SPSA with _max_evals_grouped already set + optimizer = SPSA() + optimizer._max_evals_grouped = 10 # Already set to a value + # Call the function + updated = _set_default_batchsize(optimizer) + # Check that the batch size was not updated + self.assertFalse(updated) + self.assertEqual(optimizer._max_evals_grouped, 10) + + def test_set_default_batchsize_not_spsa(self): + """Test that the function does not update when not an SPSA instance.""" + # Create a mock optimizer that is not an instance of SPSA + optimizer = COBYLA() + optimizer._max_evals_grouped = None # COBYLA doesn't need the actual implementation + # Call the function + updated = _set_default_batchsize(optimizer) + # Check that the batch size was not updated + self.assertFalse(updated) diff --git a/test/utils/test_validate_bounds.py b/test/utils/test_validate_bounds.py new file mode 100644 index 000000000..6a8a8597e --- /dev/null +++ b/test/utils/test_validate_bounds.py @@ -0,0 +1,52 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 import QiskitAlgorithmsTestCase + +from unittest.mock import Mock + +import numpy as np + +from qiskit_optimization.utils import algorithm_globals, validate_bounds + + +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) diff --git a/test/utils/test_validate_initial_point.py b/test/utils/test_validate_initial_point.py new file mode 100644 index 000000000..2ca013577 --- /dev/null +++ b/test/utils/test_validate_initial_point.py @@ -0,0 +1,49 @@ +# This code is part of a Qiskit project. +# +# (C) Copyright IBM 2022, 2024. +# +# 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 import QiskitAlgorithmsTestCase + +from unittest.mock import Mock + +import numpy as np + +from qiskit_optimization.utils import algorithm_globals, validate_initial_point + + +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 mismatched parameters and bounds..""" + self.ansatz.parameter_bounds = None + with self.assertRaises(ValueError): + _ = validate_initial_point([1.0, 2.0], self.ansatz)