[new opmath 2] New LinearCombination class to succeed `qml.Hamilton…

…ian` (#5216)

Branching: #5269 >
#5216 >
#5322 >
#5335

The basic idea is to have a clone of `qml.Hamiltonian` with a better
name, `LinearCombination` - that inherits from `CompositeOp` and plays
nice with new opmath. The motivation is that sometimes it is still handy
to have an operator class for which you _know_ there is no funny nesting
and that has `coeffs` and `ops` separately on demand without processing.

ToDo

- [x] basic init
- [x] `tests/ops/op_math/test_linear_combination.py` pass
- [x] update map_wires
- [x] all tests pass
- [x] add pauli_rep
- [x] add operands attribute and make iterable (?)
- [x] update string repr
- [x] upgrade to internally use new opmath
- [x] update tests with default new opmath
- [x] dunder math methods support with new opmath
- [x] simplify and speed-up `simplify()`, also make sure to not act
in-place
- [x] Could we make the matrix generation workflow more general like
Sum?
- [x] Also, deferring to the pauli rep method for matrix generation
leads to substantial performance improvements.
- [x] toggle `__use_new_opmath`
- [x] tests_linear_combination pass
- [x] Integration tests
- [x] get all linear combination tests to pass locally
- [x] utilize super().sparse_matrix
- [x] use grouping functionality from Sum that Mudit is adding in
#5179
- [x] update matmul to handle other LinearCombination instances
- [x] remove top level import
- [x] remove unnecessary copies
- [x] add diagonalizing gates
- [x] bugfix trivial case of simplify empty LinearCombination
- [x] add trivial case for is_hermitian
- [ ] torch differentiation with simplify=True in constructor
- [x] changelog

Optional features (likely to be included in a follow-up)
- [x] `qml.Hamiltonian` points to either old Hamiltonian or
LinearCombination depending on `__use_new_opmath`
- [x] support Hermitian
- [ ] matrix method (optional)
- [ ] adjoint method (optional)
- [ ] batching support (optional)
- [ ] qml.dot points to LinearCombination
- [ ] take care of xfails (mostly about raiseing errors in other parts
of the codebase)
- [ ] update logic of adjoint differentiation to catch attempt to
differentiate lincomb coeffs
- [x] support in qml.equal
- [x]  add special matmul with other product or non-composite ops
- [x] add special matmul with other LinearCombination


[sc-56704]

---------

Co-authored-by: lillian542 <38584660+lillian542@users.noreply.github.com>
Co-authored-by: Astral Cai <astral.cai@xanadu.ai>
Co-authored-by: Mudit Pandey <mudit.pandey@xanadu.ai>
Co-authored-by: Christina Lee <christina@xanadu.ai>
Co-authored-by: Pietropaolo Frisoni <pietropaolo.frisoni@xanadu.ai>
Co-authored-by: albi3ro <chrissie.c.l@gmail.com>
Co-authored-by: Utkarsh <utkarshazad98@gmail.com>
Co-authored-by: lillian542 <Lillian.frederiksen@xanadu.ai>
Co-authored-by: Alex Preciado <alex.preciado@xanadu.ai>
Co-authored-by: Thomas R. Bromley <49409390+trbromley@users.noreply.github.com>
Co-authored-by: Vincent Michaud-Rioux <vincentm@nanoacademic.com>
Co-authored-by: Josh Izaac <josh146@gmail.com>
Co-authored-by: Nathan Killoran <co9olguy@users.noreply.github.com>
Co-authored-by: Matthew Silverman <matthews@xanadu.ai>
Co-authored-by: Mikhail Andrenkov <mikhail@xanadu.ai>
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

doc/development/deprecations.rst

@@ -34,6 +34,20 @@ Pending deprecations
 
   - Added and deprecated for ``Sum`` and ``Prod`` instances in v0.35
 
+* Accessing ``qml.ops.Hamiltonian`` with new operator arithmetic enabled is deprecated. Using ``qml.Hamiltonian``
+  with new operator arithmetic enabled now returns a ``LinearCombination`` instance. Some functionality
+  may not work as expected, and use of the Hamiltonian class with the new operator arithmetic will not
+  be supported in future releases of PennyLane.
+
+  You can update your code to the new operator arithmetic by using ``qml.Hamiltonian`` instead of importing
+  the Hamiltonian class directly or via ``qml.ops.Hamiltonian``. When the new operator arithmetic is enabled, 
+  ``qml.Hamiltonian`` will access the new corresponding implementation. 
+
+  Alternatively, to continue accessing the legacy functionality, you can use 
+  ``qml.operation.disable_new_opmath()``.
+
+  - Deprecated in v0.36
+
 Completed deprecation cycles
 ----------------------------
 

doc/releases/changelog-dev.md

@@ -64,6 +64,14 @@
 * Added new function `qml.operation.convert_to_legacy_H` to convert `Sum`, `SProd`, and `Prod` to `Hamiltonian` instances.
   [(#5309)](https://github.com/PennyLaneAI/pennylane/pull/5309)
 
+<h3>Improvements 🛠</h3>
+
+* The `qml.is_commuting` function now accepts `Sum`, `SProd`, and `Prod` instances.
+  [(#5351)](https://github.com/PennyLaneAI/pennylane/pull/5351)
+
+* Operators can now be left multiplied `x * op` by numpy arrays.
+  [(#5361)](https://github.com/PennyLaneAI/pennylane/pull/5361)
+
 * Create the `qml.Reflection` operator, useful for amplitude amplification and its variants.
   [(#5159)](https://github.com/PennyLaneAI/pennylane/pull/5159)
 
@@ -123,6 +131,10 @@
 
   ```
 
+* A new class `qml.ops.LinearCombination` is introduced. In essence, this class is an updated equivalent of `qml.ops.Hamiltonian`
+  but for usage with new operator arithmetic.
+  [(#5216)](https://github.com/PennyLaneAI/pennylane/pull/5216)
+
 * The `qml.TrotterProduct` operator now supports error estimation functionality. 
   [(#5384)](https://github.com/PennyLaneAI/pennylane/pull/5384)
 
@@ -147,9 +159,14 @@
 * The `molecular_hamiltonian` function calls `PySCF` directly when `method='pyscf'` is selected.
   [(#5118)](https://github.com/PennyLaneAI/pennylane/pull/5118)
 
-* All generators in the source code (except those in the `qchem` module) no longer return
-  `Hamiltonian` or `Tensor` instances. Wherever possible, these return `Sum`, `SProd`, and `Prod` instances.
+* The generators in the source code return operators consistent with the global setting for 
+  `qml.operator.active_new_opmath()` wherever possible. `Sum`, `SProd` and `Prod` instances 
+  will be returned even after disabling the new operator arithmetic in cases where they offer 
+  additional functionality not available using legacy operators.
   [(#5253)](https://github.com/PennyLaneAI/pennylane/pull/5253)
+  [(#5410)](https://github.com/PennyLaneAI/pennylane/pull/5410)
+  [(#5411)](https://github.com/PennyLaneAI/pennylane/pull/5411) 
+  [(#5421)](https://github.com/PennyLaneAI/pennylane/pull/5421)
 
 * Upgraded `null.qubit` to the new device API. Also, added support for all measurements and various modes of differentiation.
   [(#5211)](https://github.com/PennyLaneAI/pennylane/pull/5211)
@@ -160,6 +177,9 @@
 * `Hamiltonian.pauli_rep` is now defined if the hamiltonian is a linear combination of paulis.
   [(#5377)](https://github.com/PennyLaneAI/pennylane/pull/5377)
 
+* `Prod.eigvals()` is now compatible with Qudit operators.
+  [(#5400)](https://github.com/PennyLaneAI/pennylane/pull/5400)
+
 * Obtaining classical shadows using the `default.clifford` device is now compatible with
   [stim](https://github.com/quantumlib/Stim) `v1.13.0`.
   [(#5409)](https://github.com/PennyLaneAI/pennylane/pull/5409)
@@ -226,6 +246,10 @@
 * Attempting to multiply `PauliWord` and `PauliSentence` with `*` will raise an error. Instead, use `@` to conform with the PennyLane convention.
   [(#5341)](https://github.com/PennyLaneAI/pennylane/pull/5341)
 
+* When new operator arithmetic is enabled, `qml.Hamiltonian` is now an alias for `qml.ops.LinearCombination`.
+  `Hamiltonian` will still be accessible as `qml.ops.Hamiltonian`.
+  [(#5393)](https://github.com/PennyLaneAI/pennylane/pull/5393)
+
 * Since `default.mixed` does not support snapshots with measurements, attempting to do so will result in a `DeviceError` instead of getting the density matrix.
   [(#5416)](https://github.com/PennyLaneAI/pennylane/pull/5416)
 
@@ -245,6 +269,11 @@
   ...     circuit = qml.from_qasm(f.read())
   ```
 
+* Accessing `qml.ops.Hamiltonian` with new operator arithmetic is deprecated. Using `qml.Hamiltonian` with new operator arithmetic enabled now
+  returns a `LinearCombination` instance. Some functionality may not work as expected. To continue using the `Hamiltonian` class, you can use
+  `qml.operation.disable_new_opmath()` to disable the new operator arithmetic.
+  [(#5393)](https://github.com/PennyLaneAI/pennylane/pull/5393)
+
 <h3>Documentation 📝</h3>
 
 * Removed some redundant documentation for the `evolve` function.

pennylane/__init__.py

@@ -145,6 +145,18 @@ class PennyLaneDeprecationWarning(UserWarning):
     """Warning raised when a PennyLane feature is being deprecated."""
 
 
+del globals()["Hamiltonian"]
+
+
+def __getattr__(name):
+    if name == "Hamiltonian":
+        if pennylane.operation.active_new_opmath():
+            return pennylane.ops.LinearCombination
+        return pennylane.ops.Hamiltonian
+
+    raise AttributeError(f"module 'pennylane' has no attribute '{name}'")
+
+
 def _get_device_entrypoints():
     """Returns a dictionary mapping the device short name to the
     loadable entrypoint"""

pennylane/_device.py

@@ -42,7 +42,7 @@
 )
 
 from pennylane.operation import Observable, Operation, Tensor, Operator, StatePrepBase
-from pennylane.ops import Hamiltonian, Sum
+from pennylane.ops import Hamiltonian, Sum, LinearCombination
 from pennylane.tape import QuantumScript, QuantumTape, expand_tape_state_prep
 from pennylane.wires import WireError, Wires
 from pennylane.queuing import QueuingManager
@@ -678,9 +678,8 @@ def default_expand_fn(self, circuit, max_expansion=10):
         )
         obs_on_same_wire = len(circuit._obs_sharing_wires) > 0 or comp_basis_sampled_multi_measure
         obs_on_same_wire &= not any(
-            isinstance(o, qml.Hamiltonian) for o in circuit._obs_sharing_wires
+            isinstance(o, (Hamiltonian, LinearCombination)) for o in circuit._obs_sharing_wires
         )
-
         ops_not_supported = not all(self.stopping_condition(op) for op in circuit.operations)
 
         if obs_on_same_wire:
@@ -745,17 +744,20 @@ def batch_transform(self, circuit: QuantumTape):
             to be applied to the list of evaluated circuit results.
         """
         supports_hamiltonian = self.supports_observable("Hamiltonian")
+
         supports_sum = self.supports_observable("Sum")
         finite_shots = self.shots is not None
         grouping_known = all(
             obs.grouping_indices is not None
             for obs in circuit.observables
-            if isinstance(obs, Hamiltonian)
+            if isinstance(obs, (Hamiltonian, LinearCombination))
         )
         # device property present in braket plugin
         use_grouping = getattr(self, "use_grouping", True)
 
-        hamiltonian_in_obs = any(isinstance(obs, Hamiltonian) for obs in circuit.observables)
+        hamiltonian_in_obs = any(
+            isinstance(obs, (Hamiltonian, LinearCombination)) for obs in circuit.observables
+        )
         expval_sum_in_obs = any(
             isinstance(m.obs, Sum) and isinstance(m, ExpectationMP) for m in circuit.measurements
         )

pennylane/_qubit_device.py

@@ -1631,7 +1631,7 @@ def adjoint_jacobian(
                     f" measurement {m.__class__.__name__}"
                 )
 
-            if m.obs.name == "Hamiltonian":
+            if not m.obs.has_matrix:
                 raise qml.QuantumFunctionError(
                     "Adjoint differentiation method does not support Hamiltonian observables."
                 )

pennylane/data/attributes/operator/operator.py

@@ -61,7 +61,13 @@ def consumes_types(cls) -> FrozenSet[Type[Operator]]:
                 qml.QubitCarry,
                 qml.QubitSum,
                 # pennylane/ops/qubit/hamiltonian.py
-                qml.Hamiltonian,
+                qml.ops.Hamiltonian,
+                # pennylane/ops/op_math/linear_combination.py
+                qml.ops.LinearCombination,
+                # pennylane/ops/op_math - prod.py, s_prod.py, sum.py
+                qml.ops.Prod,
+                qml.ops.SProd,
+                qml.ops.Sum,
                 # pennylane/ops/qubit/matrix_qml.py
                 qml.QubitUnitary,
                 qml.DiagonalQubitUnitary,
@@ -206,6 +212,8 @@ def _ops_to_hdf5(
         op_class_names = []
         for i, op in enumerate(value):
             op_key = f"op_{i}"
+            if isinstance(op, (qml.ops.Prod, qml.ops.SProd, qml.ops.Sum)):
+                op = op.simplify()
             if type(op) not in self.consumes_types():
                 raise TypeError(
                     f"Serialization of operator type '{type(op).__name__}' is not supported."
@@ -214,11 +222,19 @@ def _ops_to_hdf5(
             if isinstance(op, Tensor):
                 self._ops_to_hdf5(bind, op_key, op.obs)
                 op_wire_labels.append("null")
-            elif isinstance(op, qml.Hamiltonian):
+            elif isinstance(op, (qml.ops.Hamiltonian, qml.ops.LinearCombination)):
                 coeffs, ops = op.terms()
                 ham_grp = self._ops_to_hdf5(bind, op_key, ops)
                 ham_grp["hamiltonian_coeffs"] = coeffs
                 op_wire_labels.append("null")
+            elif isinstance(op, (qml.ops.Prod, qml.ops.Sum)):
+                self._ops_to_hdf5(bind, op_key, op.operands)
+                op_wire_labels.append("null")
+            elif isinstance(op, qml.ops.SProd):
+                coeffs, ops = op.terms()
+                sprod_grp = self._ops_to_hdf5(bind, op_key, ops)
+                sprod_grp["sprod_scalar"] = coeffs
+                op_wire_labels.append("null")
             else:
                 bind[op_key] = op.data if len(op.data) else h5py.Empty("f")
                 op_wire_labels.append(wires_to_json(op.wires))
@@ -238,21 +254,29 @@ def _hdf5_to_ops(self, bind: HDF5Group) -> List[Operator]:
         wires_bind = bind["op_wire_labels"]
         op_class_names = [] if names_bind.shape == (0,) else names_bind.asstr()
         op_wire_labels = [] if wires_bind.shape == (0,) else wires_bind.asstr()
-
         with qml.QueuingManager.stop_recording():
             for i, op_class_name in enumerate(op_class_names):
                 op_key = f"op_{i}"
 
                 op_cls = self._supported_ops_dict()[op_class_name]
                 if op_cls is Tensor:
                     ops.append(Tensor(*self._hdf5_to_ops(bind[op_key])))
-                elif op_cls is qml.Hamiltonian:
+                elif op_cls in (qml.ops.Hamiltonian, qml.ops.LinearCombination):
                     ops.append(
                         qml.Hamiltonian(
                             coeffs=list(bind[op_key]["hamiltonian_coeffs"]),
                             observables=self._hdf5_to_ops(bind[op_key]),
                         )
                     )
+                elif op_cls in (qml.ops.Prod, qml.ops.Sum):
+                    ops.append(op_cls(*self._hdf5_to_ops(bind[op_key])))
+                elif op_cls is qml.ops.SProd:
+                    ops.append(
+                        qml.ops.s_prod(
+                            scalar=bind[op_key]["sprod_scalar"][0],
+                            operator=self._hdf5_to_ops(bind[op_key])[0],
+                        )
+                    )
                 else:
                     wire_labels = json.loads(op_wire_labels[i])
                     op_data = bind[op_key]

pennylane/devices/default_clifford.py

@@ -75,6 +75,7 @@
     "Identity",
     "Projector",
     "Hamiltonian",
+    "LinearCombination",
     "Sum",
     "SProd",
     "Prod",
@@ -1017,9 +1018,9 @@ def _measure_single_sample(stim_ct, meas_ops, meas_idx, meas_wire):
         """Sample a single qubit Pauli measurement from a stim circuit"""
         stim_sm = stim.TableauSimulator()
         stim_sm.do_circuit(stim_ct)
-        return stim_sm.measure_observable(
-            stim.PauliString([0] * meas_idx + meas_ops + [0] * (meas_wire - meas_idx - 1))
-        )
+        res = [0] * meas_idx + meas_ops + [0] * (meas_wire - meas_idx - 1)
+        res = [int(r) for r in res]
+        return stim_sm.measure_observable(stim.PauliString(res))
 
     def _sample_classical_shadow(self, meas, stim_circuit, shots, seed):
         """Measures classical shadows from the state of simulator device"""

pennylane/devices/default_qubit.py

@@ -68,6 +68,7 @@
     "Projector",
     "SparseHamiltonian",
     "Hamiltonian",
+    "LinearCombination",
     "Sum",
     "SProd",
     "Prod",

pennylane/devices/default_qubit_legacy.py

@@ -196,6 +196,7 @@ class DefaultQubitLegacy(QubitDevice):
         "Projector",
         "SparseHamiltonian",
         "Hamiltonian",
+        "LinearCombination",
         "Sum",
         "SProd",
         "Prod",
@@ -597,7 +598,7 @@ def expval(self, observable, shot_range=None, bin_size=None):
         # intercept other Hamiltonians
         # TODO: Ideally, this logic should not live in the Device, but be moved
         # to a component that can be re-used by devices as needed.
-        if observable.name not in ("Hamiltonian", "SparseHamiltonian"):
+        if observable.name not in ("Hamiltonian", "SparseHamiltonian", "LinearCombination"):
             return super().expval(observable, shot_range=shot_range, bin_size=bin_size)
 
         assert self.shots is None, f"{observable.name} must be used with shots=None"
@@ -606,7 +607,7 @@ def expval(self, observable, shot_range=None, bin_size=None):
         backprop_mode = (
             not isinstance(self.state, np.ndarray)
             or any(not isinstance(d, (float, np.ndarray)) for d in observable.data)
-        ) and observable.name == "Hamiltonian"
+        ) and observable.name in ["Hamiltonian", "LinearCombination"]
 
         if backprop_mode:
             # TODO[dwierichs]: This branch is not adapted to broadcasting yet
@@ -666,7 +667,7 @@ def expval(self, observable, shot_range=None, bin_size=None):
                     csr_matrix.dot(Hmat, csr_matrix(state[..., None])),
                 ).toarray()[0]
 
-        if observable.name == "Hamiltonian":
+        if observable.name in ["Hamiltonian", "LinearCombination"]:
             res = qml.math.squeeze(res)
 
         return self._real(res)
@@ -1086,6 +1087,7 @@ def _get_diagonalizing_gates(self, circuit: qml.tape.QuantumTape) -> List[Operat
         meas_filtered = [
             m
             for m in circuit.measurements
-            if m.obs is None or not isinstance(m.obs, qml.Hamiltonian)
+            if m.obs is None
+            or not isinstance(m.obs, (qml.ops.Hamiltonian, qml.ops.LinearCombination))
         ]
         return super()._get_diagonalizing_gates(qml.tape.QuantumScript(measurements=meas_filtered))

pennylane/devices/default_qutrit.py

@@ -33,6 +33,7 @@
 OMEGA = qml.math.exp(2 * np.pi * 1j / 3)
 
 
+# pylint: disable=too-many-arguments
 class DefaultQutrit(QutritDevice):
     """Default qutrit device for PennyLane.
 
@@ -87,11 +88,7 @@ class DefaultQutrit(QutritDevice):
 
     # Identity is supported as an observable for qml.state() to work correctly. However, any
     # measurement types that rely on eigenvalue decomposition will not work with qml.Identity
-    observables = {
-        "THermitian",
-        "GellMann",
-        "Identity",
-    }
+    observables = {"THermitian", "GellMann", "Identity", "Prod"}
 
     # Static methods to use qml.math to allow for backprop differentiation
     _reshape = staticmethod(qml.math.reshape)

pennylane/devices/qubit/measure.py

@@ -19,7 +19,7 @@
 from scipy.sparse import csr_matrix
 
 from pennylane import math
-from pennylane.ops import Sum, Hamiltonian
+from pennylane.ops import Sum, Hamiltonian, LinearCombination
 from pennylane.measurements import (
     StateMeasurement,
     MeasurementProcess,
@@ -194,7 +194,7 @@ def get_measurement_function(
                 return full_dot_products
 
             backprop_mode = math.get_interface(state, *measurementprocess.obs.data) != "numpy"
-            if isinstance(measurementprocess.obs, Hamiltonian):
+            if isinstance(measurementprocess.obs, (Hamiltonian, LinearCombination)):
                 # need to work out thresholds for when its faster to use "backprop mode" measurements
                 return sum_of_terms_method if backprop_mode else csr_dot_products
 

pennylane/devices/qubit/sampling.py

@@ -16,7 +16,7 @@
 
 import numpy as np
 import pennylane as qml
-from pennylane.ops import Sum, Hamiltonian, SProd, Prod
+from pennylane.ops import Sum, Hamiltonian, SProd, Prod, LinearCombination
 from pennylane.measurements import (
     SampleMeasurement,
     Shots,
@@ -59,7 +59,7 @@ def _group_measurements(mps: List[Union[SampleMeasurement, ClassicalShadowMP, Sh
         elif mp.obs is None:
             mp_no_obs.append(mp)
             mp_no_obs_indices.append(i)
-        elif isinstance(mp.obs, (Sum, Hamiltonian, SProd, Prod)):
+        elif isinstance(mp.obs, (Hamiltonian, Sum, SProd, Prod)):
             # Sum, Hamiltonian, SProd, and Prod are treated as valid Pauli words, but
             # aren't accepted in qml.pauli.group_observables
             mp_other_obs.append([mp])
@@ -108,7 +108,9 @@ def get_num_shots_and_executions(tape: qml.tape.QuantumTape) -> Tuple[int, int]:
     num_executions = 0
     num_shots = 0
     for group in groups:
-        if isinstance(group[0], ExpectationMP) and isinstance(group[0].obs, qml.Hamiltonian):
+        if isinstance(group[0], ExpectationMP) and isinstance(
+            group[0].obs, (qml.ops.Hamiltonian, qml.ops.LinearCombination)
+        ):
             indices = group[0].obs.grouping_indices
             H_executions = len(indices) if indices else len(group[0].obs.ops)
             num_executions += H_executions
@@ -184,7 +186,9 @@ def measure_with_samples(
 
     all_res = []
     for group in groups:
-        if isinstance(group[0], ExpectationMP) and isinstance(group[0].obs, Hamiltonian):
+        if isinstance(group[0], ExpectationMP) and isinstance(
+            group[0].obs, (Hamiltonian, LinearCombination)
+        ):
             measure_fn = _measure_hamiltonian_with_samples
         elif isinstance(group[0], ExpectationMP) and isinstance(group[0].obs, Sum):
             measure_fn = _measure_sum_with_samples