add support for Parameters to SparseLabelOp (qiskit-community#891)

* add support for Parameters to SparseLabelOp

* fix isinstance for python < 3.10

* improve docs

* chop zero Parameter

* Update test/second_q/operators/test_sparse_label_op.py

Co-authored-by: Max Rossmannek <[email protected]>

* fix test

* add is_parameterized and check it for unsupported operations

* add assign_parameters

* fix chop doc

* lint

* add parameters method

* use dict instead of Mapping for data

* return self instead of None for inplace

* Revert "use dict instead of Mapping for data"

This reverts commit 3480c79.

* remove inplace option for assign_parameters

* address comments

* Update qiskit_nature/second_q/operators/sparse_label_op.py

* fix equiv

Co-authored-by: Max Rossmannek <[email protected]>

qiskit_nature/second_q/mappers/qubit_mapper.py

@@ -17,6 +17,7 @@
 from abc import ABC, abstractmethod
 from functools import lru_cache
 
+import numpy as np
 from qiskit.opflow import PauliSumOp
 from qiskit.quantum_info.operators import Pauli, SparsePauliOp
 
@@ -182,7 +183,7 @@ def mode_based_mapping(
         else:
             for terms, coeff in second_q_op.terms():
                 # 1. Initialize an operator list with the identity scaled by the `coeff`
-                ret_op = SparsePauliOp("I" * nmodes, coeffs=[coeff])
+                ret_op = SparsePauliOp("I" * nmodes, coeffs=np.array([coeff]))
 
                 # Go through the label and replace the fermion operators by their qubit-equivalent, then
                 # save the respective Pauli string in the pauli_str list.
@@ -202,4 +203,4 @@ def mode_based_mapping(
                         )
                 ret_op_list.append(ret_op)
 
-        return PauliSumOp(SparsePauliOp.sum(ret_op_list).simplify().chop())
+        return PauliSumOp(SparsePauliOp.sum(ret_op_list).simplify())

qiskit_nature/second_q/operators/fermionic_op.py

@@ -27,7 +27,7 @@
 
 from ._bits_container import _BitsContainer
 from .polynomial_tensor import PolynomialTensor
-from .sparse_label_op import SparseLabelOp
+from .sparse_label_op import _TCoeff, SparseLabelOp, _to_number
 
 
 class FermionicOp(SparseLabelOp):
@@ -149,6 +149,14 @@ class FermionicOp(SparseLabelOp):
             considered a lower bound, which means that mathematical operations acting on two or more
             operators will result in a new operator with the maximum number of spin orbitals of any
             of the involved operators.
+
+    .. note::
+
+        A FermionicOp can contain :class:`qiskit.circuit.ParameterExpression` objects as coefficients.
+        However, a FermionicOp containing parameters does not support the following methods:
+
+        - ``is_hermitian``
+        - ``to_matrix``
     """
 
     _OPERATION_REGEX = re.compile(r"([\+\-]_\d+\s)*[\+\-]_\d+")
@@ -285,7 +293,7 @@ def __str__(self) -> str:
         )
         return pre + ret
 
-    def terms(self) -> Iterator[tuple[list[tuple[str, int]], complex]]:
+    def terms(self) -> Iterator[tuple[list[tuple[str, int]], _TCoeff]]:
         """Provides an iterator analogous to :meth:`items` but with the labels already split into
         pairs of operation characters and indices.
 
@@ -353,7 +361,12 @@ def to_matrix(self, sparse: bool | None = True) -> csc_matrix | np.ndarray:
 
         Returns:
             The matrix of the operator in the Fock basis
+
+        Raises:
+            ValueError: Operator contains parameters.
         """
+        if self.is_parameterized():
+            raise ValueError("to_matrix is not supported for operators containing parameters.")
 
         csc_data, csc_col, csc_row = [], [], []
 
@@ -445,11 +458,11 @@ def normal_order(self) -> FermionicOp:
             {
                 label: coeff
                 for label, coeff in ordered_op.items()
-                if not np.isclose(coeff, 0.0, atol=self.atol)
+                if not np.isclose(_to_number(coeff), 0.0, atol=self.atol)
             }
         )
 
-    def _normal_order(self, terms: list[tuple[str, int]], coeff: complex) -> FermionicOp:
+    def _normal_order(self, terms: list[tuple[str, int]], coeff: _TCoeff) -> FermionicOp:
         if not terms:
             return self._new_instance({"": coeff})
 
@@ -550,25 +563,32 @@ def is_hermitian(self, atol: float | None = None) -> bool:
 
         Returns:
             True if the operator is hermitian up to numerical tolerance, False otherwise.
+
+        Raises:
+            ValueError: Operator contains parameters.
         """
+        if self.is_parameterized():
+            raise ValueError("is_hermitian is not supported for operators containing parameters.")
         atol = self.atol if atol is None else atol
         diff = (self - self.adjoint()).normal_order().simplify(atol=atol)
         return all(np.isclose(coeff, 0.0, atol=atol) for coeff in diff.values())
 
     def simplify(self, atol: float | None = None) -> FermionicOp:
         atol = self.atol if atol is None else atol
 
-        data = defaultdict(complex)  # type: dict[str, complex]
+        data = defaultdict(complex)  # type: dict[str, _TCoeff]
         # TODO: use parallel_map to make this more efficient (?)
         for label, coeff in self.items():
             label, coeff = self._simplify_label(label, coeff)
             data[label] += coeff
         simplified_data = {
-            label: coeff for label, coeff in data.items() if not np.isclose(coeff, 0.0, atol=atol)
+            label: coeff
+            for label, coeff in data.items()
+            if not np.isclose(_to_number(coeff), 0.0, atol=atol)
         }
         return self._new_instance(simplified_data)
 
-    def _simplify_label(self, label: str, coeff: complex) -> tuple[str, complex]:
+    def _simplify_label(self, label: str, coeff: _TCoeff) -> tuple[str, _TCoeff]:
         bits = _BitsContainer[int]()
 
         # Since Python 3.7, dictionaries are guaranteed to be insert-order preserving. We use this

qiskit_nature/second_q/operators/sparse_label_op.py

@@ -16,12 +16,12 @@
 
 from abc import ABC, abstractmethod
 from collections.abc import Collection, Mapping
-from numbers import Number
-from typing import Iterator, Sequence
+from numbers import Complex
+from typing import Iterator, Sequence, Union
 
 import cmath
 import numpy as np
-
+from qiskit.circuit import ParameterExpression
 from qiskit.quantum_info.operators.mixins import (
     AdjointMixin,
     GroupMixin,
@@ -32,6 +32,16 @@
 from .polynomial_tensor import PolynomialTensor
 
 
+_TCoeff = Union[complex, ParameterExpression]
+
+
+def _to_number(a: _TCoeff) -> complex:
+    if isinstance(a, ParameterExpression):
+        sympified = a.sympify()
+        return complex(sympified) if sympified.is_Number else np.nan
+    return a
+
+
 class SparseLabelOp(LinearMixin, AdjointMixin, GroupMixin, TolerancesMixin, ABC, Mapping):
     """The base class for sparse second-quantized operators.
 
@@ -48,11 +58,19 @@ class SparseLabelOp(LinearMixin, AdjointMixin, GroupMixin, TolerancesMixin, ABC,
     - equality and equivalence (using the :attr:`atol` and :attr:`rtol` tolerances) comparisons
 
     Furthermore, several general utility methods exist which are documented below.
+
+    .. note::
+
+        A SparseLabelOp can contain :class:`qiskit.circuit.ParameterExpression` objects as coefficients.
+        However, a SparseLabelOp containing parameters does not support the following methods:
+
+        - ``equiv``
+        - ``induced_norm``
     """
 
     def __init__(
         self,
-        data: Mapping[str, complex],
+        data: Mapping[str, _TCoeff],
         *,
         copy: bool = True,
         validate: bool = True,
@@ -72,7 +90,7 @@ def __init__(
         Raises:
             QiskitNatureError: when an invalid key is encountered during validation.
         """
-        self._data: Mapping[str, complex] = {}
+        self._data: Mapping[str, _TCoeff] = {}
         if copy:
             if validate:
                 self._validate_keys(data.keys())
@@ -191,7 +209,7 @@ def _add(self, other: SparseLabelOp, qargs: None = None) -> SparseLabelOp:
 
         return self._new_instance(new_data, other=other)
 
-    def _multiply(self, other: complex) -> SparseLabelOp:
+    def _multiply(self, other: _TCoeff) -> SparseLabelOp:
         """Return scalar multiplication of self and other.
 
         Args:
@@ -203,7 +221,7 @@ def _multiply(self, other: complex) -> SparseLabelOp:
         Raises:
             TypeError: if ``other`` is not compatible type (int, float or complex)
         """
-        if not isinstance(other, Number):
+        if not isinstance(other, (Complex, ParameterExpression)):
             raise TypeError(
                 f"Unsupported operand type(s) for *: 'SparseLabelOp' and '{type(other).__name__}'"
             )
@@ -296,20 +314,23 @@ def equiv(
     ) -> bool:
         """Check equivalence of two ``SparseLabelOp`` instances up to an accepted tolerance.
 
-        The absolute and relative tolerances can be changed via the `atol` and `rtol` attributes,
-        respectively.
-
         Args:
             other: the second ``SparseLabelOp`` to compare with this instance.
             atol: Absolute numerical tolerance. The default behavior is to use ``self.atol``.
             rtol: Relative numerical tolerance. The default behavior is to use ``self.rtol``.
 
         Returns:
             True if operators are equivalent, False if not.
+
+        Raises:
+            ValueError: Raised if either operator contains parameters
         """
         if not isinstance(other, self.__class__):
             return False
 
+        if self.is_parameterized() or other.is_parameterized():
+            raise ValueError("Cannot compare an operator that contains parameters.")
+
         if self._data.keys() != other._data.keys():
             return False
 
@@ -335,7 +356,7 @@ def __eq__(self, other: object) -> bool:
 
         return self._data == other._data
 
-    def __getitem__(self, __k: str) -> complex:
+    def __getitem__(self, __k: str) -> _TCoeff:
         """Get the requested element of the ``SparseLabelOp``."""
         return self._data.__getitem__(__k)
 
@@ -348,7 +369,7 @@ def __iter__(self) -> Iterator[str]:
         return self._data.__iter__()
 
     @abstractmethod
-    def terms(self) -> Iterator[tuple[list[tuple[str, int]], complex]]:
+    def terms(self) -> Iterator[tuple[list[tuple[str, int]], _TCoeff]]:
         """Provides an iterator analogous to :meth:`items` but with the labels already split into
         pairs of operation characters and indices.
 
@@ -394,10 +415,10 @@ def sort(self, *, weight: bool = False) -> SparseLabelOp:
         return self._new_instance({ind: self[ind] for ind in indices})
 
     def chop(self, atol: float | None = None) -> SparseLabelOp:
-        """Chops the real and imaginary phases of the operator coefficients.
+        """Chops the real and imaginary parts of the operator coefficients.
 
-        This function separately chops the real and imaginary phase of all coefficients to the
-        provided tolerance.
+        This function separately chops the real and imaginary parts of all coefficients to the
+        provided tolerance. Parameters are chopped only if they are exactly zero.
 
         Args:
             atol: the tolerance to which to chop. If ``None``, :attr:`atol` will be used.
@@ -409,16 +430,18 @@ def chop(self, atol: float | None = None) -> SparseLabelOp:
 
         new_data = {}
         for key, value in self.items():
-            zero_real = cmath.isclose(value.real, 0.0, abs_tol=atol)
-            zero_imag = cmath.isclose(value.imag, 0.0, abs_tol=atol)
-            if zero_real and zero_imag:
+            if _to_number(value) == 0:
                 continue
-            if zero_imag:
-                new_data[key] = value.real
-            elif zero_real:
-                new_data[key] = value.imag * 1j
-            else:
-                new_data[key] = value
+            if not isinstance(value, ParameterExpression):
+                zero_real = cmath.isclose(value.real, 0.0, abs_tol=atol)
+                zero_imag = cmath.isclose(value.imag, 0.0, abs_tol=atol)
+                if zero_real and zero_imag:
+                    continue
+                if zero_imag:
+                    value = value.real
+                elif zero_real:
+                    value = value.imag * 1j
+            new_data[key] = value
 
         return self._new_instance(new_data)
 
@@ -463,9 +486,33 @@ def induced_norm(self, order: int = 1) -> float:
 
         .. _https://en.wikipedia.org/wiki/Norm_(mathematics)#p-norm:
             https://en.wikipedia.org/wiki/Norm_(mathematics)#p-norm
+
+        Raises:
+            ValueError: Operator contains parameters.
         """
+        if self.is_parameterized():
+            raise ValueError("Cannot compute norm of an operator that contains parameters.")
         return sum(abs(coeff) ** order for coeff in self.values()) ** (1 / order)
 
+    def is_parameterized(self) -> bool:
+        """Returns whether the operator contains any parameters."""
+        return any(isinstance(coeff, ParameterExpression) for coeff in self.values())
+
+    def assign_parameters(self, parameters: Mapping[ParameterExpression, _TCoeff]) -> SparseLabelOp:
+        """Assign parameters to new parameters or values.
+
+        Args:
+            parameters: The mapping from parameters to new parameters or values.
+
+        Returns:
+            A new operator with the parameters assigned.
+        """
+        data = {
+            key: parameters[value] if value in parameters else value
+            for key, value in self._data.items()
+        }
+        return self._new_instance(data, other=self)
+
     def round(self, decimals: int = 0) -> SparseLabelOp:
         """Rounds the operator coefficients to a specified number of decimal places.
 
@@ -494,3 +541,11 @@ def is_zero(self, tol: int | None = None) -> bool:
             return True
         tol = tol if tol is not None else self.atol
         return all(np.isclose(val, 0, atol=tol) for val in self._data.values())
+
+    def parameters(self) -> list[ParameterExpression]:
+        """Returns a list of the parameters in the operator.
+
+        Returns:
+            A list of the parameters in the operator.
+        """
+        return [coeff for coeff in self.values() if isinstance(coeff, ParameterExpression)]

test/second_q/mappers/test_jordan_wigner_mapper.py

@@ -15,6 +15,8 @@
 import unittest
 from test import QiskitNatureTestCase
 
+from qiskit.circuit import Parameter
+from qiskit.quantum_info import SparsePauliOp
 from qiskit.opflow import I, PauliSumOp, X, Y, Z
 
 import qiskit_nature.optionals as _optionals
@@ -92,6 +94,12 @@ def test_mapping_for_single_op(self):
             expected = PauliSumOp.from_list([("I", 1)])
             self.assertEqual(JordanWignerMapper().map(op), expected)
 
+        with self.subTest("test parameters"):
+            a = Parameter("a")
+            op = FermionicOp({"+_0": a})
+            expected = SparsePauliOp.from_list([("X", 0.5 * a), ("Y", -0.5j * a)], dtype=object)
+            self.assertEqual(JordanWignerMapper().map(op).primitive, expected)
+
 
 if __name__ == "__main__":
     unittest.main()