Use eigh, eigsh in NumPyEigensolver

qiskit/algorithms/eigen_solvers/numpy_eigen_solver.py

@@ -13,15 +13,16 @@
 """The Eigensolver algorithm."""
 
 import logging
-from typing import List, Optional, Union, Tuple, Callable
+from typing import Callable, List, Optional, Tuple, Union
 
 import numpy as np
 from scipy import sparse as scisparse
 
-from qiskit.opflow import OperatorBase, I, StateFn, ListOp
+from qiskit.opflow import I, ListOp, OperatorBase, StateFn
 from qiskit.utils.validation import validate_min
-from .eigen_solver import Eigensolver, EigensolverResult
+
 from ..exceptions import AlgorithmError
+from .eigen_solver import Eigensolver, EigensolverResult
 
 logger = logging.getLogger(__name__)
 
@@ -125,11 +126,19 @@ def _solve(self, operator: OperatorBase) -> None:
         else:
             if self._k >= 2 ** operator.num_qubits - 1:
                 logger.debug("SciPy doesn't support to get all eigenvalues, using NumPy instead.")
-                eigval, eigvec = np.linalg.eig(operator.to_matrix())
+                if operator.is_hermitian():
+                    eigval, eigvec = np.linalg.eigh(operator.to_matrix())
+                else:
+                    eigval, eigvec = np.linalg.eig(operator.to_matrix())
             else:
-                eigval, eigvec = scisparse.linalg.eigs(
-                    operator.to_spmatrix(), k=self._k, which="SR"
-                )
+                if operator.is_hermitian():
+                    eigval, eigvec = scisparse.linalg.eigsh(
+                        operator.to_spmatrix(), k=self._k, which="SA"
+                    )
+                else:
+                    eigval, eigvec = scisparse.linalg.eigs(
+                        operator.to_spmatrix(), k=self._k, which="SR"
+                    )
             indices = np.argsort(eigval)[: self._k]
             eigval = eigval[indices]
             eigvec = eigvec[:, indices]

qiskit/opflow/operator_base.py

@@ -170,6 +170,13 @@ def to_spmatrix(self) -> spmatrix:
         """
         return csr_matrix(self.to_matrix())
 
+    def is_hermitian(self) -> bool:
+        """Return True if the operator is hermitian.
+
+        Returns: Boolean value
+        """
+        return (self.to_spmatrix() != self.to_spmatrix().getH()).nnz == 0
+
     @staticmethod
     def _indent(lines: str, indentation: str = INDENTATION) -> str:
         """Indented representation to allow pretty representation of nested operators."""

qiskit/opflow/primitive_ops/pauli_sum_op.py

@@ -446,3 +446,6 @@ def is_zero(self) -> bool:
         op = self.reduce()
         primitive: SparsePauliOp = op.primitive
         return op.coeff == 1 and len(op) == 1 and primitive.coeffs[0] == 0
+
+    def is_hermitian(self):
+        return np.isreal(self.coeffs).all() and np.all(self.primitive.paulis.phase == 0)

releasenotes/notes/add-opflow-is-hermitian-method-6a461549e3c6b32c.yaml

@@ -0,0 +1,7 @@
+---
+features:
+  - |
+    The :meth:`~qiskit.opflow.OperatorBase.is_hermitian` method is added to :obj:`~qiskit.opflow.OperatorBase` to check whether
+    the operator is Hermitian or not. :class:`~qiskit.algorithms.NumPyEigensolver` and
+    :class:`~qiskit.algorithms.NumPyMinimumEigensolver` use eigh or eigsh to solve the
+    eigenvalue problem when the operator is Hermitian.

test/python/algorithms/test_numpy_eigen_solver.py

@@ -44,16 +44,18 @@ def test_ce(self):
         result = algo.compute_eigenvalues(operator=self.qubit_op, aux_operators=[])
         self.assertEqual(len(result.eigenvalues), 1)
         self.assertEqual(len(result.eigenstates), 1)
-        self.assertAlmostEqual(result.eigenvalues[0], -1.85727503 + 0j)
+        self.assertEqual(result.eigenvalues.dtype, np.float64)
+        self.assertAlmostEqual(result.eigenvalues[0], -1.85727503)
 
     def test_ce_k4(self):
         """Test for k=4 eigenvalues"""
         algo = NumPyEigensolver(k=4)
         result = algo.compute_eigenvalues(operator=self.qubit_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.real, [-1.85727503, -1.24458455, -0.88272215, -0.22491125]
+            result.eigenvalues, [-1.85727503, -1.24458455, -0.88272215, -0.22491125]
         )
 
     def test_ce_k4_filtered(self):
@@ -68,7 +70,8 @@ def criterion(x, v, a_v):
         result = algo.compute_eigenvalues(operator=self.qubit_op, aux_operators=[])
         self.assertEqual(len(result.eigenvalues), 2)
         self.assertEqual(len(result.eigenstates), 2)
-        np.testing.assert_array_almost_equal(result.eigenvalues.real, [-0.88272215, -0.22491125])
+        self.assertEqual(result.eigenvalues.dtype, np.float64)
+        np.testing.assert_array_almost_equal(result.eigenvalues, [-0.88272215, -0.22491125])
 
     def test_ce_k4_filtered_empty(self):
         """Test for k=4 eigenvalues with filter always returning False"""

test/python/algorithms/test_numpy_minimum_eigen_solver.py

@@ -46,7 +46,8 @@ def test_cme(self):
         """Basic test"""
         algo = NumPyMinimumEigensolver()
         result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=self.aux_ops)
-        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+        self.assertEqual(result.eigenvalue.dtype, np.float64)
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503)
         self.assertEqual(len(result.aux_operator_eigenvalues), 2)
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[0], [2, 0])
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[1], [0, 0])
@@ -56,24 +57,28 @@ def test_cme_reuse(self):
         # Start with no operator or aux_operators, give via compute method
         algo = NumPyMinimumEigensolver()
         result = algo.compute_minimum_eigenvalue(operator=self.qubit_op)
-        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+        self.assertEqual(result.eigenvalue.dtype, np.float64)
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503)
         self.assertIsNone(result.aux_operator_eigenvalues)
 
         # Add aux_operators and go again
         result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=self.aux_ops)
-        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+        self.assertEqual(result.eigenvalue.dtype, np.float64)
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503)
         self.assertEqual(len(result.aux_operator_eigenvalues), 2)
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[0], [2, 0])
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[1], [0, 0])
 
         # "Remove" aux_operators and go again
         result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=[])
-        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+        self.assertEqual(result.eigenvalue.dtype, np.float64)
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503)
         self.assertIsNone(result.aux_operator_eigenvalues)
 
         # Set aux_operators and go again
         result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=self.aux_ops)
-        self.assertAlmostEqual(result.eigenvalue, -1.85727503 + 0j)
+        self.assertEqual(result.eigenvalue.dtype, np.float64)
+        self.assertAlmostEqual(result.eigenvalue, -1.85727503)
         self.assertEqual(len(result.aux_operator_eigenvalues), 2)
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[0], [2, 0])
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[1], [0, 0])
@@ -93,7 +98,8 @@ def criterion(x, v, a_v):
 
         algo = NumPyMinimumEigensolver(filter_criterion=criterion)
         result = algo.compute_minimum_eigenvalue(operator=self.qubit_op, aux_operators=self.aux_ops)
-        self.assertAlmostEqual(result.eigenvalue, -0.22491125 + 0j)
+        self.assertEqual(result.eigenvalue.dtype, np.float64)
+        self.assertAlmostEqual(result.eigenvalue, -0.22491125)
         self.assertEqual(len(result.aux_operator_eigenvalues), 2)
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[0], [2, 0])
         np.testing.assert_array_almost_equal(result.aux_operator_eigenvalues[1], [0, 0])

test/python/opflow/test_op_construction.py

@@ -24,7 +24,13 @@
 from scipy.stats import unitary_group
 
 from qiskit import QiskitError
-from qiskit.circuit import Instruction, Parameter, ParameterVector, QuantumCircuit, QuantumRegister
+from qiskit.circuit import (
+    Instruction,
+    Parameter,
+    ParameterVector,
+    QuantumCircuit,
+    QuantumRegister,
+)
 from qiskit.circuit.library import CZGate, ZGate
 from qiskit.extensions.exceptions import ExtensionError
 from qiskit.opflow import (
@@ -1092,6 +1098,29 @@ def test_listop_num_qubits(self):
             with self.assertRaises(ValueError):
                 X @ op  # pylint: disable=pointless-statement
 
+    def test_is_hermitian(self):
+        """Test is_hermitian method."""
+        with self.subTest("I"):
+            self.assertTrue(I.is_hermitian())
+
+        with self.subTest("X"):
+            self.assertTrue(X.is_hermitian())
+
+        with self.subTest("Y"):
+            self.assertTrue(Y.is_hermitian())
+
+        with self.subTest("Z"):
+            self.assertTrue(Z.is_hermitian())
+
+        with self.subTest("XY"):
+            self.assertFalse((X @ Y).is_hermitian())
+
+        with self.subTest("CX"):
+            self.assertTrue(CX.is_hermitian())
+
+        with self.subTest("T"):
+            self.assertFalse(T.is_hermitian())
+
 
 @ddt
 class TestListOpMethods(QiskitOpflowTestCase):

test/python/opflow/test_pauli_sum_op.py

@@ -306,6 +306,32 @@ def test_matrix_iter_sparse(self):
                 np.array_equal(i.toarray(), coeff * coeffs[idx] * Pauli(labels[idx]).to_matrix())
             )
 
+    def test_is_hermitian(self):
+        """Test is_hermitian method"""
+        with self.subTest("True test"):
+            target = PauliSumOp.from_list(
+                [
+                    ("II", -1.052373245772859),
+                    ("IZ", 0.39793742484318045),
+                    ("ZI", -0.39793742484318045),
+                    ("ZZ", -0.01128010425623538),
+                    ("XX", 0.18093119978423156),
+                ]
+            )
+            self.assertTrue(target.is_hermitian())
+
+        with self.subTest("False test"):
+            target = PauliSumOp.from_list(
+                [
+                    ("II", -1.052373245772859),
+                    ("IZ", 0.39793742484318045j),
+                    ("ZI", -0.39793742484318045),
+                    ("ZZ", -0.01128010425623538),
+                    ("XX", 0.18093119978423156),
+                ]
+            )
+            self.assertFalse(target.is_hermitian())
+
 
 if __name__ == "__main__":
     unittest.main()