Parameter handling in SparsePauliOp

qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py

@@ -46,6 +46,31 @@ class SparsePauliOp(LinearOp):
     using the :attr:`~SparsePauliOp.paulis` attribute. The coefficients
     are stored as a complex Numpy array vector and can be accessed using
     the :attr:`~SparsePauliOp.coeffs` attribute.
+
+    .. rubric:: Data type of coefficients
+
+    The default ``dtype`` of the internal ``coeffs`` Numpy array is ``complex128``.  Users can
+    configure this by passing ``np.ndarray`` with a different dtype.  For example, a parameterized
+    :class:`SparsePauliOp` can be made as follows:
+
+    .. code-block:: python
+
+        >>> import numpy as np
+        >>> from qiskit.circuit import ParameterVector
+        >>> from qiskit.quantum_info import SparsePauliOp
+
+        >>> SparsePauliOp(["II", "XZ"], np.array(ParameterVector("a", 2)))
+        SparsePauliOp(['II', 'XZ'],
+              coeffs=[ParameterExpression(1.0*a[0]), ParameterExpression(1.0*a[1])])
+
+    .. note::
+
+      Parameterized :class:`SparsePauliOp` does not support the following methods:
+
+      - ``to_matrix(sparse=True)`` since ``scipy.sparse`` cannot have objects as elements.
+      - ``to_operator()`` since :class:`~.quantum_info.Operator` does not support objects.
+      - ``sort``, ``argsort`` since :class:`.ParameterExpression` does not support comparison.
+      - ``equiv`` since :class:`.ParameterExpression`. cannot be converted into complex.
     """
 
     def __init__(self, data, coeffs=None, *, ignore_pauli_phase=False, copy=True):
@@ -86,10 +111,12 @@ def __init__(self, data, coeffs=None, *, ignore_pauli_phase=False, copy=True):
 
         pauli_list = PauliList(data.copy() if copy and hasattr(data, "copy") else data)
 
+        dtype = coeffs.dtype if isinstance(coeffs, np.ndarray) else complex
+
         if coeffs is None:
-            coeffs = np.ones(pauli_list.size, dtype=complex)
+            coeffs = np.ones(pauli_list.size, dtype=dtype)
         else:
-            coeffs = np.array(coeffs, copy=copy, dtype=complex)
+            coeffs = np.array(coeffs, copy=copy, dtype=dtype)
 
         if ignore_pauli_phase:
             # Fast path used in copy operations, where the phase of the PauliList is already known
@@ -101,7 +128,7 @@ def __init__(self, data, coeffs=None, *, ignore_pauli_phase=False, copy=True):
             # move the phase of `pauli_list` to `self._coeffs`
             phase = pauli_list._phase
             count_y = pauli_list._count_y()
-            self._coeffs = np.asarray((-1j) ** (phase - count_y) * coeffs, dtype=complex)
+            self._coeffs = np.asarray((-1j) ** (phase - count_y) * coeffs, dtype=coeffs.dtype)
             pauli_list._phase = np.mod(count_y, 4)
             self._pauli_list = pauli_list
 
@@ -132,9 +159,14 @@ def __eq__(self, other):
         """Entrywise comparison of two SparsePauliOp operators"""
         return (
             super().__eq__(other)
+            and self.coeffs.dtype == other.coeffs.dtype
             and self.coeffs.shape == other.coeffs.shape
-            and np.allclose(self.coeffs, other.coeffs)
             and self.paulis == other.paulis
+            and (
+                np.allclose(self.coeffs, other.coeffs)
+                if self.coeffs.dtype != object
+                else (self.coeffs == other.coeffs).all()
+            )
         )
 
     def equiv(self, other, atol: Optional[float] = None):
@@ -385,7 +417,19 @@ def simplify(self, atol=None, rtol=None):
             rtol = self.rtol
 
         # Filter non-zero coefficients
-        non_zero = np.logical_not(np.isclose(self.coeffs, 0, atol=atol, rtol=rtol))
+        if self.coeffs.dtype == object:
+
+            def to_complex(coeff):
+                if not hasattr(coeff, "sympify"):
+                    return coeff
+                sympified = coeff.sympify()
+                return complex(sympified) if sympified.is_Number else np.nan
+
+            non_zero = np.logical_not(
+                np.isclose([to_complex(x) for x in self.coeffs], 0, atol=atol, rtol=rtol)
+            )
+        else:
+            non_zero = np.logical_not(np.isclose(self.coeffs, 0, atol=atol, rtol=rtol))
         paulis_x = self.paulis.x[non_zero]
         paulis_z = self.paulis.z[non_zero]
         nz_coeffs = self.coeffs[non_zero]
@@ -398,16 +442,21 @@ def simplify(self, atol=None, rtol=None):
             # No zero operator or duplicate operator
             return self.copy()
 
-        coeffs = np.zeros(indexes.shape[0], dtype=complex)
+        coeffs = np.zeros(indexes.shape[0], dtype=self.coeffs.dtype)
         np.add.at(coeffs, inverses, nz_coeffs)
         # Delete zero coefficient rows
-        is_zero = np.isclose(coeffs, 0, atol=atol, rtol=rtol)
+        if self.coeffs.dtype == object:
+            is_zero = np.array(
+                [np.isclose(to_complex(coeff), 0, atol=atol, rtol=rtol) for coeff in coeffs]
+            )
+        else:
+            is_zero = np.isclose(coeffs, 0, atol=atol, rtol=rtol)
         # Check edge case that we deleted all Paulis
         # In this case we return an identity Pauli with a zero coefficient
         if np.all(is_zero):
             x = np.zeros((1, self.num_qubits), dtype=bool)
             z = np.zeros((1, self.num_qubits), dtype=bool)
-            coeffs = np.array([0j], dtype=complex)
+            coeffs = np.array([0j], dtype=self.coeffs.dtype)
         else:
             non_zero = np.logical_not(is_zero)
             non_zero_indexes = indexes[non_zero]
@@ -659,7 +708,7 @@ def from_operator(obj, atol=None, rtol=None):
         return SparsePauliOp(paulis, coeffs, copy=False)
 
     @staticmethod
-    def from_list(obj):
+    def from_list(obj, dtype=complex):
         """Construct from a list of Pauli strings and coefficients.
 
         For example, the 5-qubit Hamiltonian
@@ -677,6 +726,7 @@ def from_list(obj):
 
         Args:
             obj (Iterable[Tuple[str, complex]]): The list of 2-tuples specifying the Pauli terms.
+            dtype (type): The dtype of coeffs (Default complex).
 
         Returns:
             SparsePauliOp: The SparsePauliOp representation of the Pauli terms.
@@ -693,7 +743,7 @@ def from_list(obj):
         # determine the number of qubits
         num_qubits = len(obj[0][0])
 
-        coeffs = np.zeros(size, dtype=complex)
+        coeffs = np.zeros(size, dtype=dtype)
         labels = np.zeros(size, dtype=f"<U{num_qubits}")
         for i, item in enumerate(obj):
             labels[i] = item[0]
@@ -703,7 +753,7 @@ def from_list(obj):
         return SparsePauliOp(paulis, coeffs, copy=False)
 
     @staticmethod
-    def from_sparse_list(obj, num_qubits, do_checks=True):
+    def from_sparse_list(obj, num_qubits, do_checks=True, dtype=complex):
         """Construct from a list of local Pauli strings and coefficients.
 
         Each list element is a 3-tuple of a local Pauli string, indices where to apply it,
@@ -729,6 +779,7 @@ def from_sparse_list(obj, num_qubits, do_checks=True):
             obj (Iterable[Tuple[str, List[int], complex]]): The list 3-tuples specifying the Paulis.
             num_qubits (int): The number of qubits of the operator.
             do_checks (bool): The flag of checking if the input indices are not duplicated.
+            dtype (type): The dtype of coeffs (Default complex).
 
         Returns:
             SparsePauliOp: The SparsePauliOp representation of the Pauli terms.
@@ -744,7 +795,7 @@ def from_sparse_list(obj, num_qubits, do_checks=True):
         if size == 0:
             raise QiskitError("Input Pauli list is empty.")
 
-        coeffs = np.zeros(size, dtype=complex)
+        coeffs = np.zeros(size, dtype=dtype)
         labels = np.zeros(size, dtype=f"<U{num_qubits}")
 
         for i, (paulis, indices, coeff) in enumerate(obj):
@@ -781,7 +832,9 @@ def to_list(self, array=False):
         """
         # Dtype for a structured array with string labels and complex coeffs
         pauli_labels = self.paulis.to_labels(array=True)
-        labels = np.zeros(self.size, dtype=[("labels", pauli_labels.dtype), ("coeffs", "c16")])
+        labels = np.zeros(
+            self.size, dtype=[("labels", pauli_labels.dtype), ("coeffs", self.coeffs.dtype)]
+        )
         labels["labels"] = pauli_labels
         labels["coeffs"] = self.coeffs
         if array:

releasenotes/notes/operator-parameters-c81b7c05bffb740b.yaml

@@ -0,0 +1,25 @@
+---
+features:
+  - |
+    :class:`.SparsePauliOp`\ s can now be constructed with coefficient arrays
+    that are general Python objects.  This is purely intended for use with Terra's
+    :class:`.ParameterExpression` objects; other objects may work, but do not
+    have first-class support.  Some :class:`.SparsePauliOp` methods (such as
+    conversion to other class representations) may not work when using ``object``
+    arrays, if the desired target cannot represent these general arrays.
+
+    For example, a :class:`.ParameterExpression` :class:`.SparsePauliOp` could
+    be constructed by::
+
+        import numpy as np
+        from qiskit.circuit import Parameter
+        from qiskit.quantum_info import SparsePauliOp
+
+        print(SparsePauliOp(["II", "XZ"], np.array([Parameter("a"), Parameter("b")])))
+
+    which gives
+
+    .. code-block:: text
+
+        SparsePauliOp(['II', 'XZ'],
+              coeffs=[ParameterExpression(1.0*a), ParameterExpression(1.0*b)])