[new opmath 1] Toggle __use_new_opmath

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)
@@ -230,6 +250,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)
 
@@ -249,6 +273,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",
@@ -428,10 +429,11 @@ def _apply_sx(self, state, axes, inverse=False):
         Returns:
             array[complex]: output state
         """
-        if inverse:
-            return 0.5 * ((1 - 1j) * state + (1 + 1j) * self._apply_x(state, axes))
-
-        return 0.5 * ((1 + 1j) * state + (1 - 1j) * self._apply_x(state, axes))
+        return (
+            0.5 * ((1 - 1j) * state + (1 + 1j) * self._apply_x(state, axes))
+            if inverse
+            else 0.5 * ((1 + 1j) * state + (1 - 1j) * self._apply_x(state, axes))
+        )
 
     def _apply_cnot(self, state, axes, **kwargs):
         """Applies a CNOT gate by slicing along the first axis specified in ``axes`` and then
@@ -597,7 +599,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 +608,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 +668,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 +1088,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