Fix NumPyEigensolver with SparsePauliOp (#9101)

* Fix NumPyEigensolver for SparsePauliOp

* Use to_matrix instead for non-Operator BaseOperator children.
* Make NumPyEigensolver stateless w.r.t. results.
* Build operators as SparsePauliOps rather than PauliSumOps.

* Test NumPyMinimumEigensolver with different ops

Test NumPyMinimumEigensolver with:

* SparsePauliOp
* Operator
* PauliSumOp

* add releasenote

* Use sparse matrix computation for non-Operator

* Use sparse matmul for non-Operator

* improve structure

* try dense matrix for all BaseOperators or raise ex

* remove unused import

* test Operator.to_matrix()

* use instance to access sparse

* raise error for invalid num_qubits. pauliop/scalarop tests

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>
Co-authored-by: Julien Gacon <gaconju@gmail.com>
Co-authored-by: mergify[bot] <37929162+mergify[bot]@users.noreply.github.com>

Author: 4 people

qiskit/algorithms/eigensolvers/numpy_eigensolver.py

@@ -20,11 +20,12 @@
 from scipy import sparse as scisparse
 
 from qiskit.opflow import PauliSumOp
+from qiskit.quantum_info import SparsePauliOp, Statevector
 from qiskit.quantum_info.operators.base_operator import BaseOperator
-from qiskit.quantum_info import Statevector
 from qiskit.utils.validation import validate_min
 
 from .eigensolver import Eigensolver, EigensolverResult
+from ..exceptions import AlgorithmError
 from ..list_or_dict import ListOrDict
 
 logger = logging.getLogger(__name__)
@@ -53,12 +54,12 @@ def __init__(
     ) -> None:
         """
         Args:
-            k: number of eigenvalues are to be computed, with a min. value of 1.
-            filter_criterion: callable that allows to filter eigenvalues/eigenstates. Only feasible
+            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
+                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)
@@ -69,8 +70,6 @@ def __init__(
 
         self._filter_criterion = filter_criterion
 
-        self._ret = NumPyEigensolverResult()
-
     @property
     def k(self) -> int:
         """Return k (number of eigenvalues requested)."""
@@ -109,69 +108,64 @@ def _check_set_k(self, operator: BaseOperator | PauliSumOp) -> None:
             else:
                 self._k = self._in_k
 
-    def _solve(self, operator: BaseOperator | PauliSumOp) -> None:
+    def _solve(self, operator: BaseOperator | PauliSumOp) -> tuple[np.ndarray, np.ndarray]:
         if isinstance(operator, PauliSumOp):
-            sp_mat = operator.to_spmatrix()
+            op_matrix = operator.to_spmatrix()
+        else:
+            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(sp_mat.diagonal()).nnz == sp_mat.nnz:
-                diag = sp_mat.diagonal()
+            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((sp_mat.shape[0], self._k))
+                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."
                     )
-                    if operator.is_hermitian():
-                        eigval, eigvec = np.linalg.eigh(operator.to_matrix())
-                    else:
-                        eigval, eigvec = np.linalg.eig(operator.to_matrix())
+                    eigval, eigvec = self._solve_dense(operator.to_matrix())
                 else:
-                    if operator.is_hermitian():
-                        eigval, eigvec = scisparse.linalg.eigsh(sp_mat, k=self._k, which="SA")
-                    else:
-                        eigval, eigvec = scisparse.linalg.eigs(sp_mat, k=self._k, which="SR")
-
-                indices = np.argsort(eigval)[: self._k]
-                eigval = eigval[indices]
-                eigvec = eigvec[:, indices]
+                    eigval, eigvec = self._solve_sparse(op_matrix, self._k)
         else:
-            logger.debug("SciPy not supported, using NumPy instead.")
-
-            if operator.data.all() == operator.data.conj().T.all():
-                eigval, eigvec = np.linalg.eigh(operator.data)
-            else:
-                eigval, eigvec = np.linalg.eig(operator.data)
-
-            indices = np.argsort(eigval)[: self._k]
-            eigval = eigval[indices]
-            eigvec = eigvec[:, indices]
-
-        self._ret.eigenvalues = eigval
-        self._ret.eigenstates = eigvec.T
+            # Sparse SciPy matrix not supported, use dense NumPy computation.
+            eigval, eigvec = self._solve_dense(operator.to_matrix())
 
-    def _get_ground_state_energy(self, operator: BaseOperator | PauliSumOp) -> None:
-        if self._ret.eigenvalues is None or self._ret.eigenstates is None:
-            self._solve(operator)
+        indices = np.argsort(eigval)[: self._k]
+        eigval = eigval[indices]
+        eigvec = eigvec[:, indices]
+        return eigval, eigvec.T
 
-    def _get_energies(
-        self,
-        operator: BaseOperator | PauliSumOp,
-        aux_operators: ListOrDict[BaseOperator | PauliSumOp] | None,
-    ) -> None:
-        if self._ret.eigenvalues is None or self._ret.eigenstates is None:
-            self._solve(operator)
+    @staticmethod
+    def _solve_sparse(op_matrix: scisparse.csr_matrix, k: int) -> tuple[np.ndarray, np.ndarray]:
+        if (op_matrix != op_matrix.H).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")
 
-        if aux_operators is not None:
-            aux_op_vals = []
-            for i in range(self._k):
-                aux_op_vals.append(
-                    self._eval_aux_operators(aux_operators, self._ret.eigenstates[i])
-                )
-            self._ret.aux_operators_evaluated = aux_op_vals
+    @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 np.linalg.eigh(op_matrix)
+        else:
+            return np.linalg.eig(op_matrix)
 
     @staticmethod
     def _eval_aux_operators(
@@ -190,25 +184,42 @@ def _eval_aux_operators(
         else:
             values = {}
             key_op_iterator = aux_operators.items()
+
         for key, operator in key_op_iterator:
             if operator is None:
                 continue
-            value = 0.0
+
+            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
             if isinstance(operator, PauliSumOp):
                 if operator.coeff != 0:
-                    mat = operator.to_spmatrix()
-                    # Terra doesn't support sparse yet, so do the matmul directly if so
-                    # This is necessary for the particle_hole and other chemistry tests because the
-                    # pauli conversions are 2^12th large and will OOM error if not sparse.
-                    if isinstance(mat, scisparse.spmatrix):
-                        value = mat.dot(wavefn).dot(np.conj(wavefn))
-                    else:
-                        value = (
-                            Statevector(wavefn).expectation_value(operator.primitive)
-                            * operator.coeff
-                        )
+                    op_matrix = operator.to_spmatrix()
             else:
+                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).
@@ -225,8 +236,12 @@ def compute_eigenvalues(
 
         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 = PauliSumOp.from_list([("I", 1)]).tensorpower(operator.num_qubits) * 0.0
+
+        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]
@@ -244,49 +259,51 @@ def compute_eigenvalues(
             # need to consider all elements if a filter is set
             self._k = 2**operator.num_qubits
 
-        self._ret = NumPyEigensolverResult()
-        self._solve(operator)
+        eigvals, eigvecs = self._solve(operator)
 
         # compute energies before filtering, as this also evaluates the aux operators
-        self._get_energies(operator, 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:
-
-            eigvecs = []
-            eigvals = []
-            aux_ops = []
-            cnt = 0
-            for i in range(len(self._ret.eigenvalues)):
-                eigvec = self._ret.eigenstates[i]
-                eigval = self._ret.eigenvalues[i]
-                if self._ret.aux_operators_evaluated is not None:
-                    aux_op = self._ret.aux_operators_evaluated[i]
+            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 = None
-                if self._filter_criterion(eigvec, eigval, aux_op):
-                    cnt += 1
-                    eigvecs += [eigvec]
-                    eigvals += [eigval]
-                    if self._ret.aux_operators_evaluated is not None:
-                        aux_ops += [aux_op]
-                if cnt == k_orig:
+                    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
 
-            self._ret.eigenstates = np.array(eigvecs)
-            self._ret.eigenvalues = np.array(eigvals)
-            # conversion to np.array breaks in case of aux_ops
-            self._ret.aux_operators_evaluated = aux_ops
+            eigvals = np.array(filt_eigvals)
+            eigvecs = np.array(filt_eigvecs)
+            aux_op_vals = filt_aux_op_vals
 
             self._k = k_orig
 
-        # evaluate ground state after filtering (in case a filter is set)
-        self._get_ground_state_energy(operator)
-        if self._ret.eigenstates is not None:
-            self._ret.eigenstates = [Statevector(vec) for vec in self._ret.eigenstates]
+        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", self._ret)
-        return self._ret
+        logger.debug("NumpyEigensolverResult:\n%s", result)
+        return result
 
 
 class NumPyEigensolverResult(EigensolverResult):

qiskit/algorithms/minimum_eigensolvers/numpy_minimum_eigensolver.py

@@ -43,7 +43,7 @@ def __init__(
     ) -> None:
         """
         Args:
-            filter_criterion: callable that allows to filter eigenvalues/eigenstates. The minimum
+            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

qiskit/quantum_info/operators/operator.py

@@ -470,6 +470,10 @@ def reverse_qargs(self):
         ret._op_shape = self._op_shape.reverse()
         return ret
 
+    def to_matrix(self):
+        """Convert operator to NumPy matrix."""
+        return self.data
+
     @classmethod
     def _einsum_matmul(cls, tensor, mat, indices, shift=0, right_mul=False):
         """Perform a contraction using Numpy.einsum

releasenotes/notes/fix-numpy-eigensolver-sparse-pauli-op-b09a9ac8fb93fe4a.yaml

@@ -0,0 +1,10 @@
+---
+fixes:
+  - |
+    Fixed an issue with
+    :class:`~qiskit.algorithms.eigensolvers.NumPyEigensolver` and by extension
+    :class:`~qiskit.algorithms.minimum_eigensolvers.NumPyMinimumEigensolver`
+    where solving for
+    :class:`~qiskit.quantum_info.operators.base_operator.BaseOperator`
+    subclasses other than :class:`~qiskit.quantum_info.operators.Operator`
+    would cause an error.