Repository: PennyLaneAI/pennylane (python)

Issue:
Re-think operation generators
Currently in PennyLane, qubit operations define their generator in the following form:

```python
operation.generator = [G, coeff]
```

where `G` can be either

* a PennyLane operation *class*, e.g., `qml.PauliX`,  or
* a numeric NumPy matrix,

and `coeff` is the coefficient in the matrix exponential `exp(1j * coeff * G * param)`.

The generators are only used for metric tensor computations, and this approach served us well back when we only had single-qubit rotations. However, our gate set has grown to include multi-qubit parametrized gates that can be applied to an *arbitrary* number of wires. This has caused two issues:

1. How do we define the generator for a multi-parameter gate such as `qml.Rot`?
2. How do we best define the generator for a multi-wire gate such as `qml.MultiRZ`, `qml.PauliRot`, `qml.CRot`, etc?

(1) is probably easy to support, we just need to think about the matrix exponential form of these gates, and generalize the generator property to support linear combinations.

(2) we currently get around by defining the generator as a NumPy matrix (see #963), but this is not ideal and scales pretty badly. Instead, rather than returning an operation class or a matrix, we should instead allow returning an operation *instance* or a matrix.

That way we could support returning operators defined over the operation wires, which avoids creating a large matrices. For example, for `qml.PauliRot("XY", wires=[0, 2])`, we should support:

```python
operation.generator = [qml.PauliX(wires=self.wires[0]) @ qml.PauliY(wires=self.wires[1]), coeff]
```


Pull request description:
Update the `Operation.generator` property
**Context:**

The `Operation.generator` property has long been badly designed, and a blocker for extending PennyLane. In its current form, it returns a 2-tuple `(generator, coeff)`, where

- `generator` is an operation _class_ or NumPy array
- `coeff` is a scaling factor

An unfortunate side effect of `generator` being a class is that it does not allow us to scale beyond simple observables, and further, does not allow us to encode wires.

**Description of the Change:**

- The generator property has been updated to an instance method, `Operation.generator()`.

- It now returns an instantiated operation, representing the generator of the instantiated operator.

- Various operators have been updated to specify the generator as either an `Observable`, `Tensor`, `Hamiltonian`, `SparseHamiltonian`, or `Hermitian`.

- A temporary utility function `get_generator` has been added to the `utils` module, to automate:

  * Extracting the matrix representation
  * Extracting the 'coefficient' if possible (only occurs if the generator is a single Pauli word)
  * Converting a Hamiltonian to a sparse matrix if there are more than 1 Pauli word present.
  * Negating the coefficient/taking the adjoint of the matrix if the operation was inverted

  This utility logic is currently needed because:

  * Extracting the matrix representation is not supported natively on `Hamiltonians` and `SparseHamiltonians`.
  * By default, calling `op.generator()` does not take into account `op.inverse()`.
  * If the generator is a single Pauli word, it is convenient to have access to both the coefficient and the observable separately.

**Benefits:** Cleaner and more scalable representation.

**Possible Drawbacks:**

- I noticed that `self` was only ever used for `self.wires`. Should we have a `compute_generator(*params, wires, **hyperparameters)` and `generator(self)` split like with other methods?

- Splitting the methods like above would also solve the inverse problem; we could have the inverse logic in `generator(self)`.

- Currently, it is not possible to negate `SparseHamiltonian` or `Hermitian` natively

- We need to standardize the operator API; calling `op.generator().matrix()` should work regardless of what generator was returned.

- Less directly related, but we should add support for measuring probabilities in the basis of a Hamiltonian, e.g., `qml.probs(op=H)`. On a state simulator, this should be straightforward, and we can simply re-use the functionality added to support `expval(H)`. For hardware, this is non-trivial to do, and likely would require a `batch_transform` to manage multiple executions under-the-hood (one per commuting group).

**Related GitHub Issues:** closes #970 


Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index c376ec705..d4849d741 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -26,9 +26,9 @@
 
 The Operator class has undergone a major refactor with the following changes:
 
-* The `diagonalizing_gates()` representation has been moved to the highest-level 
-  `Operator` class and is therefore available to all subclasses. A condition 
-  `qml.operation.defines_diagonalizing_gates` has been added, which can be used 
+* The `diagonalizing_gates()` representation has been moved to the highest-level
+  `Operator` class and is therefore available to all subclasses. A condition
+  `qml.operation.defines_diagonalizing_gates` has been added, which can be used
   in tape contexts without queueing.
   [(#1985)](https://github.com/PennyLaneAI/pennylane/pull/1985)
 
@@ -38,9 +38,34 @@ The Operator class has undergone a major refactor with the following changes:
 * The `string_for_inverse` attribute is removed.
   [(#2021)](https://github.com/PennyLaneAI/pennylane/pull/2021)
 
+* The generator property has been updated to an instance method,
+  `Operation.generator()`. It now returns an instantiated operation,
+  representing the generator of the instantiated operator.
+  [(#2030)](https://github.com/PennyLaneAI/pennylane/pull/2030)
+
+  Various operators have been updated to specify the generator as either
+  an `Observable`, `Tensor`, `Hamiltonian`, `SparseHamiltonian`, or `Hermitian`
+  operator.
+
+  In addition, a temporary utility function get_generator has been added
+  to the utils module, to automate:
+
+  - Extracting the matrix representation
+  - Extracting the 'coefficient' if possible (only occurs if the generator is a single Pauli word)
+  - Converting a Hamiltonian to a sparse matrix if there are more than 1 Pauli word present.
+  - Negating the coefficient/taking the adjoint of the matrix if the operation was inverted
+
+  This utility logic is currently needed because:
+
+  - Extracting the matrix representation is not supported natively on
+    Hamiltonians and SparseHamiltonians.
+  - By default, calling `op.generator()` does not take into account `op.inverse()`.
+  - If the generator is a single Pauli word, it is convenient to have access to
+    both the coefficient and the observable separately.
+
 <h3>Contributors</h3>
 
 This release contains contributions from (in alphabetical order):
 
-Olivia Di Matteo, Maria Schuld, David Wierichs
+Olivia Di Matteo, Josh Izaac, Christina Lee, Maria Schuld, David Wierichs.
 
diff --git a/pennylane/fourier/utils.py b/pennylane/fourier/utils.py
index 9adf430ea..0905b3e69 100644
--- a/pennylane/fourier/utils.py
+++ b/pennylane/fourier/utils.py
@@ -18,6 +18,9 @@ from itertools import product, combinations
 import numpy as np
 
 
+from pennylane.utils import get_generator
+
+
 def format_nvec(nvec):
     """Nice strings representing tuples of integers."""
 
@@ -60,23 +63,9 @@ def get_spectrum(op, decimals):
     Returns:
         set[float]: non-negative frequencies contributed by this input-encoding gate
     """
-    no_generator = False
-    if hasattr(op, "generator"):
-        g, coeff = op.generator
-
-        if isinstance(g, np.ndarray):
-            matrix = g
-        elif hasattr(g, "matrix"):
-            matrix = g.matrix
-        else:
-            no_generator = True
-    else:
-        no_generator = True
-
-    if no_generator:
-        raise ValueError(f"Generator of operation {op} is not defined.")
-
+    matrix, coeff = get_generator(op, return_matrix=True)
     matrix = coeff * matrix
+
     # todo: use qml.math.linalg once it is tested properly
     evals = np.linalg.eigvalsh(matrix)
 
diff --git a/pennylane/operation.py b/pennylane/operation.py
index ae95c5e13..8c525424b 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -723,30 +723,27 @@ class Operation(Operator):
         param_shift = default_param_shift if recipe is None else recipe
         return param_shift
 
-    @property
     def generator(self):
-        r"""Generator of the operation.
-
-        A length-2 list ``[generator, scaling_factor]``, where
+        r"""list[.Operation] or None: Generator of an operation
+        with a single trainable parameter.
 
-        * ``generator`` is an existing PennyLane
-          operation class or :math:`2\times 2` Hermitian array
-          that acts as the generator of the current operation
+        For example, for operator
 
-        * ``scaling_factor`` represents a scaling factor applied
-          to the generator operation
+        .. math::
 
-        For example, if :math:`U(\theta)=e^{i0.7\theta \sigma_x}`, then
-        :math:`\sigma_x`, with scaling factor :math:`s`, is the generator
-        of operator :math:`U(\theta)`:
+            U(\phi) = e^{i\phi (0.5 Y + Z\otimes X)}
 
-        .. code-block:: python
+        >>> U.generator()
+          (0.5) [Y0]
+        + (1.0) [Z0 X1]
 
-            generator = [PauliX, 0.7]
+        The generator may also be provided in the form of a dense or sparse Hamiltonian
+        (using :class:`.Hermitian` and :class:`.SparseHamiltonian` respectively).
 
-        Default is ``[None, 1]``, indicating the operation has no generator.
+        The default value to return is ``None``, indicating that the operation has
+        no defined generator.
         """
-        return [None, 1]
+        return None
 
     @property
     def inverse(self):
@@ -1796,25 +1793,7 @@ def operation_derivative(operation) -> np.ndarray:
         ValueError: if the operation does not have a generator or is not composed of a single
             trainable parameter
     """
-    generator, prefactor = operation.generator
-
-    if generator is None:
-        raise ValueError(f"Operation {operation.name} does not have a generator")
-    if operation.num_params != 1:
-        # Note, this case should already be caught by the previous raise since we haven't worked out
-        # how to have an operator for multiple parameters. It is added here in case of a future
-        # change
-        raise ValueError(
-            f"Operation {operation.name} is not written in terms of a single parameter"
-        )
-
-    if not isinstance(generator, np.ndarray):
-        generator = generator.matrix
-
-    if operation.inverse:
-        prefactor *= -1
-        generator = qml.math.conj(qml.math.T(generator))
-
+    generator, prefactor = qml.utils.get_generator(operation, return_matrix=True)
     return 1j * prefactor * generator @ operation.matrix
 
 
@@ -1827,7 +1806,7 @@ def not_tape(obj):
 @qml.BooleanFn
 def has_gen(obj):
     """Returns ``True`` if an operator has a generator defined."""
-    return hasattr(obj, "generator") and obj.generator[0] is not None
+    return hasattr(obj, "generator") and obj.generator() is not None
 
 
 @qml.BooleanFn
diff --git a/pennylane/ops/qubit/observables.py b/pennylane/ops/qubit/observables.py
index e1d516cab..17b632331 100644
--- a/pennylane/ops/qubit/observables.py
+++ b/pennylane/ops/qubit/observables.py
@@ -145,6 +145,9 @@ class SparseHamiltonian(Observable):
     num_wires = AllWires
     grad_method = None
 
+    def __init__(self, H, wires=None, do_queue=True, id=None):
+        super().__init__(H, wires=wires, do_queue=do_queue, id=id)
+
     @property
     def num_params(self):
         return 1
@@ -247,6 +250,14 @@ class Projector(Observable):
         w[idx] = 1
         return w
 
+    @classmethod
+    def _matrix(cls, *params):
+        basis_state = params[0]
+        m = np.zeros((2 ** len(basis_state), 2 ** len(basis_state)))
+        idx = int("".join(str(i) for i in basis_state), 2)
+        m[idx, idx] = 1
+        return m
+
     def diagonalizing_gates(self):
         """Return the gate set that diagonalizes a circuit according to the
         specified Projector observable.
diff --git a/pennylane/ops/qubit/parametric_ops.py b/pennylane/ops/qubit/parametric_ops.py
index 8c64a9718..cc8e752dd 100644
--- a/pennylane/ops/qubit/parametric_ops.py
+++ b/pennylane/ops/qubit/parametric_ops.py
@@ -18,6 +18,8 @@ core parameterized gates.
 # pylint:disable=abstract-method,arguments-differ,protected-access
 import functools
 import math
+from operator import matmul
+
 import numpy as np
 
 import pennylane as qml
@@ -53,7 +55,9 @@ class RX(Operation):
     num_wires = 1
     basis = "X"
     grad_method = "A"
-    generator = [PauliX, -1 / 2]
+
+    def generator(self):
+        return -0.5 * PauliX(wires=self.wires)
 
     @property
     def num_params(self):
@@ -110,7 +114,9 @@ class RY(Operation):
     num_wires = 1
     basis = "Y"
     grad_method = "A"
-    generator = [PauliY, -1 / 2]
+
+    def generator(self):
+        return -0.5 * PauliY(wires=self.wires)
 
     @property
     def num_params(self):
@@ -161,7 +167,9 @@ class RZ(Operation):
     num_wires = 1
     basis = "Z"
     grad_method = "A"
-    generator = [PauliZ, -1 / 2]
+
+    def generator(self):
+        return -0.5 * PauliZ(wires=self.wires)
 
     @property
     def num_params(self):
@@ -223,7 +231,9 @@ class PhaseShift(Operation):
     num_wires = 1
     basis = "Z"
     grad_method = "A"
-    generator = [np.array([[0, 0], [0, 1]]), 1]
+
+    def generator(self):
+        return qml.Projector(np.array([1]), wires=self.wires)
 
     @property
     def num_params(self):
@@ -297,7 +307,9 @@ class ControlledPhaseShift(Operation):
     num_wires = 2
     basis = "Z"
     grad_method = "A"
-    generator = [np.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1]]), 1]
+
+    def generator(self):
+        return qml.Projector(np.array([1, 1]), wires=self.wires)
 
     @property
     def num_params(self):
@@ -489,13 +501,8 @@ class MultiRZ(Operation):
 
         return multi_Z_rot_matrix
 
-    _generator = None
-
-    @property
     def generator(self):
-        if self._generator is None:
-            self._generator = [np.diag(pauli_eigs(len(self.wires))), -1 / 2]
-        return self._generator
+        return -0.5 * functools.reduce(matmul, [qml.PauliZ(w) for w in self.wires])
 
     @property
     def matrix(self):
@@ -590,7 +597,6 @@ class PauliRot(Operation):
 
     def __init__(self, theta, pauli_word, wires=None, do_queue=True):
         super().__init__(theta, pauli_word, wires=wires, do_queue=do_queue)
-
         if not PauliRot._check_pauli_word(pauli_word):
             raise ValueError(
                 f'The given Pauli word "{pauli_word}" contains characters that are not allowed.'
@@ -706,46 +712,9 @@ class PauliRot(Operation):
             list(range(len(pauli_word))),
         )
 
-    _generator = None
-
-    @property
     def generator(self):
-        if self._generator is None:
-            pauli_word = self.parameters[1]
-
-            # Simplest case is if the Pauli is the identity matrix
-            if pauli_word == "I" * len(pauli_word):
-                self._generator = [np.eye(2 ** len(pauli_word)), -1 / 2]
-                return self._generator
-
-            # We first generate the matrix excluding the identity parts and expand it afterwards.
-            # To this end, we have to store on which wires the non-identity parts act
-            non_identity_wires, non_identity_gates = zip(
-                *[(wire, gate) for wire, gate in enumerate(pauli_word) if gate != "I"]
-            )
-
-            # get MultiRZ's generator
-            multi_Z_rot_generator = qml.math.diag(pauli_eigs(len(non_identity_gates)))
-
-            # now we conjugate with Hadamard and RX to create the Pauli string
-            conjugation_matrix = functools.reduce(
-                qml.math.kron,
-                [PauliRot._PAULI_CONJUGATION_MATRICES[gate] for gate in non_identity_gates],
-            )
-
-            self._generator = [
-                expand(
-                    qml.math.dot(
-                        qml.math.conj(qml.math.T(conjugation_matrix)),
-                        qml.math.dot(multi_Z_rot_generator, conjugation_matrix),
-                    ),
-                    non_identity_wires,
-                    list(range(len(pauli_word))),
-                ),
-                -1 / 2,
-            ]
-
-        return self._generator
+        pauli_word = self.parameters[1]
+        return -0.5 * qml.grouping.string_to_pauli_word(pauli_word)
 
     @classmethod
     def _eigvals(cls, theta, pauli_word):
@@ -841,10 +810,8 @@ class CRX(Operation):
     grad_method = "A"
     grad_recipe = four_term_grad_recipe
 
-    generator = [
-        np.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]),
-        -1 / 2,
-    ]
+    def generator(self):
+        return -0.5 * qml.Projector(np.array([1]), wires=self.wires[0]) @ qml.PauliX(self.wires[1])
 
     @property
     def num_params(self):
@@ -942,10 +909,8 @@ class CRY(Operation):
     grad_method = "A"
     grad_recipe = four_term_grad_recipe
 
-    generator = [
-        np.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, -1j], [0, 0, 1j, 0]]),
-        -1 / 2,
-    ]
+    def generator(self):
+        return -0.5 * qml.Projector(np.array([1]), wires=self.wires[0]) @ qml.PauliY(self.wires[1])
 
     @property
     def num_params(self):
@@ -1037,10 +1002,8 @@ class CRZ(Operation):
     grad_method = "A"
     grad_recipe = four_term_grad_recipe
 
-    generator = [
-        np.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]]),
-        -1 / 2,
-    ]
+    def generator(self):
+        return -0.5 * qml.Projector(np.array([1]), wires=self.wires[0]) @ qml.PauliZ(self.wires[1])
 
     @property
     def num_params(self):
@@ -1223,7 +1186,9 @@ class U1(Operation):
     """
     num_wires = 1
     grad_method = "A"
-    generator = [np.array([[0, 0], [0, 1]]), 1]
+
+    def generator(self):
+        return qml.Projector(np.array([1]), wires=self.wires)
 
     @property
     def num_params(self):
@@ -1432,10 +1397,8 @@ class IsingXX(Operation):
     num_wires = 2
     grad_method = "A"
 
-    generator = [
-        np.array([[0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0], [1, 0, 0, 0]]),
-        -1 / 2,
-    ]
+    def generator(self):
+        return -0.5 * PauliX(wires=self.wires[0]) @ PauliX(wires=self.wires[1])
 
     @property
     def num_params(self):
@@ -1495,10 +1458,9 @@ class IsingYY(Operation):
     """
     num_wires = 2
     grad_method = "A"
-    generator = [
-        np.array([[0, 0, 0, -1], [0, 0, 1, 0], [0, 1, 0, 0], [-1, 0, 0, 0]]),
-        -1 / 2,
-    ]
+
+    def generator(self):
+        return -0.5 * PauliY(wires=self.wires[0]) @ PauliY(wires=self.wires[1])
 
     @property
     def num_params(self):
@@ -1556,10 +1518,9 @@ class IsingZZ(Operation):
     """
     num_wires = 2
     grad_method = "A"
-    generator = [
-        np.array([[1, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]]),
-        -1 / 2,
-    ]
+
+    def generator(self):
+        return -0.5 * PauliZ(wires=self.wires[0]) @ PauliZ(wires=self.wires[1])
 
     @property
     def num_params(self):
diff --git a/pennylane/ops/qubit/qchem_ops.py b/pennylane/ops/qubit/qchem_ops.py
index 7a5e66bd3..fa5ecb368 100644
--- a/pennylane/ops/qubit/qchem_ops.py
+++ b/pennylane/ops/qubit/qchem_ops.py
@@ -18,6 +18,7 @@ from quantum chemistry applications.
 # pylint:disable=abstract-method,arguments-differ,protected-access
 import math
 import numpy as np
+from scipy.sparse import coo_matrix
 
 import pennylane as qml
 from pennylane.operation import Operation
@@ -80,10 +81,10 @@ class SingleExcitation(Operation):
     num_wires = 2
     grad_method = "A"
     grad_recipe = four_term_grad_recipe
-    generator = [
-        np.array([[0, 0, 0, 0], [0, 0, -1j, 0], [0, 1j, 0, 0], [0, 0, 0, 0]]),
-        -1 / 2,
-    ]
+
+    def generator(self):
+        w1, w2 = self.wires
+        return 0.25 * qml.PauliX(w1) @ qml.PauliY(w2) - 0.25 * qml.PauliY(w1) @ qml.PauliX(w2)
 
     @property
     def num_params(self):
@@ -142,10 +143,15 @@ class SingleExcitationMinus(Operation):
     """
     num_wires = 2
     grad_method = "A"
-    generator = [
-        np.array([[1, 0, 0, 0], [0, 0, -1j, 0], [0, 1j, 0, 0], [0, 0, 0, 1]]),
-        -1 / 2,
-    ]
+
+    def generator(self):
+        w1, w2 = self.wires
+        return (
+            -0.25 * qml.Identity(w1)
+            + 0.25 * qml.PauliX(w1) @ qml.PauliY(w2)
+            - 0.25 * qml.PauliY(w1) @ qml.PauliX(w2)
+            - 0.25 * qml.PauliZ(w1) @ qml.PauliZ(w2)
+        )
 
     @property
     def num_params(self):
@@ -218,10 +224,15 @@ class SingleExcitationPlus(Operation):
     """
     num_wires = 2
     grad_method = "A"
-    generator = [
-        np.array([[-1, 0, 0, 0], [0, 0, -1j, 0], [0, 1j, 0, 0], [0, 0, 0, -1]]),
-        -1 / 2,
-    ]
+
+    def generator(self):
+        w1, w2 = self.wires
+        return (
+            0.25 * qml.Identity(w1)
+            + 0.25 * qml.PauliX(w1) @ qml.PauliY(w2)
+            - 0.25 * qml.PauliY(w1) @ qml.PauliX(w2)
+            + 0.25 * qml.PauliZ(w1) @ qml.PauliZ(w2)
+        )
 
     @property
     def num_params(self):
@@ -320,10 +331,20 @@ class DoubleExcitation(Operation):
     grad_method = "A"
     grad_recipe = four_term_grad_recipe
 
-    G = np.zeros((16, 16), dtype=np.complex64)
-    G[3, 12] = -1j  # 3 (dec) = 0011 (bin)
-    G[12, 3] = 1j  # 12 (dec) = 1100 (bin)
-    generator = [G, -1 / 2]
+    def generator(self):
+        w0, w1, w2, w3 = self.wires
+        coeffs = [0.0625, 0.0625, -0.0625, 0.0625, -0.0625, 0.0625, -0.0625, -0.0625]
+        obs = [
+            qml.PauliX(w0) @ qml.PauliX(w1) @ qml.PauliX(w2) @ qml.PauliY(w3),
+            qml.PauliX(w0) @ qml.PauliX(w1) @ qml.PauliY(w2) @ qml.PauliX(w3),
+            qml.PauliX(w0) @ qml.PauliY(w1) @ qml.PauliX(w2) @ qml.PauliX(w3),
+            qml.PauliX(w0) @ qml.PauliY(w1) @ qml.PauliY(w2) @ qml.PauliY(w3),
+            qml.PauliY(w0) @ qml.PauliX(w1) @ qml.PauliX(w2) @ qml.PauliX(w3),
+            qml.PauliY(w0) @ qml.PauliX(w1) @ qml.PauliY(w2) @ qml.PauliY(w3),
+            qml.PauliY(w0) @ qml.PauliY(w1) @ qml.PauliX(w2) @ qml.PauliY(w3),
+            qml.PauliY(w0) @ qml.PauliY(w1) @ qml.PauliY(w2) @ qml.PauliX(w3),
+        ]
+        return qml.Hamiltonian(coeffs, obs)
 
     @property
     def num_params(self):
@@ -415,12 +436,13 @@ class DoubleExcitationPlus(Operation):
     num_wires = 4
     grad_method = "A"
 
-    G = -1 * np.eye(16, dtype=np.complex64)
-    G[3, 3] = 0
-    G[12, 12] = 0
-    G[3, 12] = -1j  # 3 (dec) = 0011 (bin)
-    G[12, 3] = 1j  # 12 (dec) = 1100 (bin)
-    generator = [G, -1 / 2]
+    def generator(self):
+        G = -1 * np.eye(16, dtype=np.complex64)
+        G[3, 3] = G[12, 12] = 0
+        G[3, 12] = -1j  # 3 (dec) = 0011 (bin)
+        G[12, 3] = 1j  # 12 (dec) = 1100 (bin)
+        H = coo_matrix(-0.5 * G)
+        return qml.SparseHamiltonian(H, wires=self.wires)
 
     @property
     def num_params(self):
@@ -487,12 +509,14 @@ class DoubleExcitationMinus(Operation):
     num_wires = 4
     grad_method = "A"
 
-    G = np.eye(16, dtype=np.complex64)
-    G[3, 3] = 0
-    G[12, 12] = 0
-    G[3, 12] = -1j  # 3 (dec) = 0011 (bin)
-    G[12, 3] = 1j  # 12 (dec) = 1100 (bin)
-    generator = [G, -1 / 2]
+    def generator(self):
+        G = np.eye(16, dtype=np.complex64)
+        G[3, 3] = 0
+        G[12, 12] = 0
+        G[3, 12] = -1j  # 3 (dec) = 0011 (bin)
+        G[12, 3] = 1j  # 12 (dec) = 1100 (bin)
+        H = coo_matrix(-0.5 * G)
+        return qml.SparseHamiltonian(H, wires=self.wires)
 
     @property
     def num_params(self):
@@ -581,29 +605,15 @@ class OrbitalRotation(Operation):
     num_wires = 4
     grad_method = "A"
     grad_recipe = four_term_grad_recipe
-    generator = [
-        qml.math.array(
-            [
-                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, -1j, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, 0, 0, 0, 0, -1j, 0, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, 0, 0, -1j, 0, 0, -1j, 0, 0, 0, 0, 0, 0],
-                [0, 1j, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 1j, 0, 0, 0, 0, 0, 0, 0, 0, -1j, 0, 0, 0],
-                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1j, 0, 0],
-                [0, 0, 1j, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 1j, 0, 0, 0, 0, 0, 0, 0, 0, -1j, 0, 0, 0],
-                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1j, 0],
-                [0, 0, 0, 0, 0, 0, 1j, 0, 0, 1j, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, 0, 0, 0, 1j, 0, 0, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1j, 0, 0, 0, 0],
-                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
-            ]
-        ),
-        -1 / 2,
-    ]
+
+    def generator(self):
+        w0, w1, w2, w3 = self.wires
+        return (
+            0.25 * qml.PauliX(w0) @ qml.PauliY(w2)
+            - 0.25 * qml.PauliY(w0) @ qml.PauliX(w2)
+            + 0.25 * qml.PauliX(w1) @ qml.PauliY(w3)
+            - 0.25 * qml.PauliY(w1) @ qml.PauliX(w3)
+        )
 
     @property
     def num_params(self):
diff --git a/pennylane/tape/reversible.py b/pennylane/tape/reversible.py
index d61afec6c..82f6b3cee 100644
--- a/pennylane/tape/reversible.py
+++ b/pennylane/tape/reversible.py
@@ -190,10 +190,10 @@ class ReversibleTape(JacobianTape):
 
         if op.name == "Rot":
             decomp = op.decompose()
-            generator, multiplier = decomp[p_idx].generator
+            generator, multiplier = qml.utils.get_generator(decomp[p_idx])
             between_ops = decomp[p_idx + 1 :] + between_ops
         else:
-            generator, multiplier = op.generator
+            generator, multiplier = qml.utils.get_generator(op)
 
         # construct circuit to compute differentiated state
         between_ops_inverse = [copy.copy(op) for op in between_ops[::-1]]
@@ -207,7 +207,7 @@ class ReversibleTape(JacobianTape):
                 op.queue().inv()
 
             # apply generator needed for differentiation
-            generator(wires)
+            qml.apply(generator)
 
             # evolve forwards again
             for op in between_ops:
diff --git a/pennylane/transforms/metric_tensor.py b/pennylane/transforms/metric_tensor.py
index 0bb985639..78fc1bdb3 100644
--- a/pennylane/transforms/metric_tensor.py
+++ b/pennylane/transforms/metric_tensor.py
@@ -24,25 +24,6 @@ from pennylane.circuit_graph import LayerData
 
 from .batch_transform import batch_transform
 
-# This dictionary maps generators to their controlled gate equivalent.
-# It is equivalent to using ``controlled`` on the generators, but this
-# approach via registering specific generators allows us to guarantee
-# unitary circuits for ``allow_nonunitary=False``
-_GEN_TO_CGEN = {
-    qml.PauliX: qml.CNOT,
-    qml.PauliY: qml.CY,
-    qml.PauliZ: qml.CZ,
-}
-
-# In contrast to the above, this dictionary maps *operators* to controlled
-# generator gates, skipping over the step of retrieving the generator from
-# the original operation and identifying it as a gate. It replaces matrix
-# arithmetics that otherwise would be required at this point.
-_OP_TO_CGEN = {
-    # PhaseShift is the same as RZ up to a global phase
-    qml.PhaseShift: qml.CZ,
-}
-
 
 def expand_fn(tape, approx=None, allow_nonunitary=True, aux_wire=None, device_wires=None):
     """Set the metric tensor based on whether non-unitary gates are allowed."""
@@ -378,27 +359,10 @@ def _metric_tensor_cov_matrix(tape, diag_approx):
 
         # for each operation in the layer, get the generator
         for op in curr_ops:
-            gen, s = op.generator
-            if op.inverse:
-                s = -s
-            w = op.wires
+            obs, s = qml.utils.get_generator(op)
+            obs_list[-1].append(obs)
             coeffs_list[-1].append(s)
 
-            # get the observable corresponding to the generator of the current operation
-            if isinstance(gen, np.ndarray):
-                # generator is a Hermitian matrix
-                obs_list[-1].append(qml.Hermitian(gen, w))
-
-            elif issubclass(gen, qml.operation.Observable):
-                # generator is an existing PennyLane operation
-                obs_list[-1].append(gen(w))
-
-            else:
-                raise qml.QuantumFunctionError(
-                    f"Can't generate metric tensor, generator {gen}"
-                    "has no corresponding observable"
-                )
-
         # Create a quantum tape with all operations
         # prior to the parametrized layer, and the rotations
         # to measure in the basis of the parametrized layer generators.
@@ -452,21 +416,22 @@ def _get_gen_op(op, allow_nonunitary, aux_wire):
     generator but ``allow_nonunitary=False``, the operation ``op`` should have been decomposed
     before, leading to a ``ValueError``.
     """
-    gen, _ = op.generator
+    op_to_cgen = {
+        qml.RX: qml.CNOT,
+        qml.RY: qml.CY,
+        qml.RZ: qml.CZ,
+        qml.PhaseShift: qml.CZ,  # PhaseShift is the same as RZ up to a global phase
+    }
+
     try:
-        if isinstance(gen, np.ndarray):
-            cgen = _OP_TO_CGEN[op.__class__]
-        else:
-            cgen = _GEN_TO_CGEN.get(gen, None)
-            if cgen is None:
-                cgen = _OP_TO_CGEN[op.__class__]
+        cgen = op_to_cgen[op.__class__]
         return cgen(wires=[aux_wire, *op.wires])
 
     except KeyError as e:
         if allow_nonunitary:
-            if (not isinstance(gen, np.ndarray)) and issubclass(gen, qml.operation.Observable):
-                gen = gen.matrix
-            return qml.ControlledQubitUnitary(gen, control_wires=aux_wire, wires=op.wires)
+            gen, coeff = qml.utils.get_generator(op, return_matrix=True)
+            return qml.ControlledQubitUnitary(coeff * gen, control_wires=aux_wire, wires=op.wires)
+
         raise ValueError(
             f"Generator for operation {op} not known and non-unitary operations "
             "deactivated via allow_nonunitary=False."
diff --git a/pennylane/utils.py b/pennylane/utils.py
index 05a6e81b3..27aae5894 100644
--- a/pennylane/utils.py
+++ b/pennylane/utils.py
@@ -417,3 +417,57 @@ def expand_vector(vector, original_wires, expanded_wires):
     )
 
     return qml.math.reshape(expanded_tensor, 2 ** M)
+
+
+def get_generator(op, return_matrix=False):
+    """Utility function to help return the generator of an operation.
+
+    This utility function is used to abstract away the API differences
+    between different operators that could be returned as
+    generators, including ``SparseHamiltonian``, ``Hamiltonian``,
+    ``Observable``, ``Tensor``, and ``Hermitian``.
+
+    This utility function should be removed once the aforementioned classes
+    no longer differ in behaviour.
+    """
+    gen = getattr(op, "generator", lambda: None)()
+
+    if gen is None:
+        raise ValueError(f"Operation {op.name} does not have a generator")
+
+    if op.num_params != 1:
+        raise ValueError(f"Operation {op.name} is not written in terms of a single parameter")
+
+    if isinstance(gen, (qml.Hermitian, qml.SparseHamiltonian)):
+        obs = gen
+        s = 1.0
+
+    elif isinstance(gen, qml.operation.Observable):
+        gen = 1.0 * gen  # convert to a qml.Hamiltonian
+
+        if len(gen.ops) == 1:
+            # case where the Hamiltonian is a single Pauli word
+            obs = gen.ops[0]
+            s = gen.coeffs[0]
+        else:
+            # otherwise, we convert to a sparse array
+            H = qml.utils.sparse_hamiltonian(gen)
+            obs = qml.SparseHamiltonian(H, wires=gen.wires)
+            s = 1.0
+
+    else:
+        raise qml.QuantumFunctionError(f"Generator {gen} is not an observable")
+
+    if op.inverse:
+        s *= -1.0
+
+    if return_matrix:
+        obs = obs.matrix
+
+        if not isinstance(obs, np.ndarray):
+            obs = obs.toarray()
+
+        if op.inverse:
+            obs = qml.math.conj(qml.math.T(obs))
+
+    return obs, s


Test patch:
diff --git a/tests/fourier/test_fourier_utils.py b/tests/fourier/test_fourier_utils.py
index 371710d3a..9e74bb327 100644
--- a/tests/fourier/test_fourier_utils.py
+++ b/tests/fourier/test_fourier_utils.py
@@ -71,9 +71,9 @@ def test_get_spectrum_complains_no_generator():
     """Test that an error is raised if the operator has no generator defined."""
 
     # Observables have no generator attribute
-    with pytest.raises(ValueError, match="Generator of operation"):
+    with pytest.raises(ValueError, match="does not have a generator"):
         get_spectrum(qml.P(wires=0), decimals=10)
 
     # CNOT is an operation where generator is an abstract property
-    with pytest.raises(ValueError, match="Generator of operation"):
+    with pytest.raises(ValueError, match="does not have a generator"):
         get_spectrum(qml.CNOT(wires=[0, 1]), decimals=10)
diff --git a/tests/ops/qubit/test_parametric_ops.py b/tests/ops/qubit/test_parametric_ops.py
index 5229d15f9..c752fbf30 100644
--- a/tests/ops/qubit/test_parametric_ops.py
+++ b/tests/ops/qubit/test_parametric_ops.py
@@ -1221,7 +1221,7 @@ class TestPauliRot:
     def test_multirz_generator(self, pauli_word):
         """Test that the generator of the MultiRZ gate is correct."""
         op = qml.PauliRot(0.3, pauli_word, wires=range(len(pauli_word)))
-        gen = op.generator
+        gen = op.generator()
 
         if pauli_word[0] == "I":
             # this is the identity
@@ -1236,10 +1236,7 @@ class TestPauliRot:
             else:
                 expected_gen = expected_gen @ getattr(qml, f"Pauli{pauli}")(wires=i)
 
-        expected_gen_mat = expected_gen.matrix
-
-        assert np.allclose(gen[0], expected_gen_mat)
-        assert gen[1] == -0.5
+        assert gen.compare(-0.5 * expected_gen)
 
     @pytest.mark.gpu
     @pytest.mark.parametrize("theta", np.linspace(0, 2 * np.pi, 7))
@@ -1382,16 +1379,13 @@ class TestMultiRZ:
     def test_multirz_generator(self, qubits, mocker):
         """Test that the generator of the MultiRZ gate is correct."""
         op = qml.MultiRZ(0.3, wires=range(qubits))
-        gen = op.generator
+        gen = op.generator()
 
         expected_gen = qml.PauliZ(wires=0)
         for i in range(1, qubits):
             expected_gen = expected_gen @ qml.PauliZ(wires=i)
 
-        expected_gen_mat = expected_gen.matrix
-
-        assert np.allclose(gen[0], expected_gen_mat)
-        assert gen[1] == -0.5
+        assert gen.compare(-0.5 * expected_gen)
 
         spy = mocker.spy(qml.utils, "pauli_eigs")
 
diff --git a/tests/ops/qubit/test_qchem_ops.py b/tests/ops/qubit/test_qchem_ops.py
index c99065024..2f0268758 100644
--- a/tests/ops/qubit/test_qchem_ops.py
+++ b/tests/ops/qubit/test_qchem_ops.py
@@ -148,7 +148,7 @@ class TestSingleExcitation:
     def test_single_excitation_generator(self, phi):
         """Tests that the SingleExcitation operation calculates the correct generator"""
         op = qml.SingleExcitation(phi, wires=[0, 1])
-        g, a = op.generator
+        g, a = qml.utils.get_generator(op, return_matrix=True)
         res = expm(1j * a * g * phi)
         exp = SingleExcitation(phi)
         assert np.allclose(res, exp)
@@ -165,7 +165,7 @@ class TestSingleExcitation:
     def test_single_excitation_plus_generator(self, phi):
         """Tests that the SingleExcitationPlus operation calculates the correct generator"""
         op = qml.SingleExcitationPlus(phi, wires=[0, 1])
-        g, a = op.generator
+        g, a = qml.utils.get_generator(op, return_matrix=True)
         res = expm(1j * a * g * phi)
         exp = SingleExcitationPlus(phi)
         assert np.allclose(res, exp)
@@ -182,7 +182,7 @@ class TestSingleExcitation:
     def test_single_excitation_minus_generator(self, phi):
         """Tests that the SingleExcitationMinus operation calculates the correct generator"""
         op = qml.SingleExcitationMinus(phi, wires=[0, 1])
-        g, a = op.generator
+        g, a = qml.utils.get_generator(op, return_matrix=True)
         res = expm(1j * a * g * phi)
         exp = SingleExcitationMinus(phi)
         assert np.allclose(res, exp)
@@ -313,7 +313,7 @@ class TestDoubleExcitation:
     def test_double_excitation_generator(self, phi):
         """Tests that the DoubleExcitation operation calculates the correct generator"""
         op = qml.DoubleExcitation(phi, wires=[0, 1, 2, 3])
-        g, a = op.generator
+        g, a = qml.utils.get_generator(op, return_matrix=True)
 
         res = expm(1j * a * g * phi)
         exp = DoubleExcitation(phi)
@@ -379,7 +379,7 @@ class TestDoubleExcitation:
     def test_double_excitation_plus_generator(self, phi):
         """Tests that the DoubleExcitationPlus operation calculates the correct generator"""
         op = qml.DoubleExcitationPlus(phi, wires=[0, 1, 2, 3])
-        g, a = op.generator
+        g, a = qml.utils.get_generator(op, return_matrix=True)
 
         res = expm(1j * a * g * phi)
         exp = DoubleExcitationPlus(phi)
@@ -398,7 +398,7 @@ class TestDoubleExcitation:
     def test_double_excitation_minus_generator(self, phi):
         """Tests that the DoubleExcitationMinus operation calculates the correct generator"""
         op = qml.DoubleExcitationMinus(phi, wires=[0, 1, 2, 3])
-        g, a = op.generator
+        g, a = qml.utils.get_generator(op, return_matrix=True)
 
         res = expm(1j * a * g * phi)
         exp = DoubleExcitationMinus(phi)
@@ -579,7 +579,7 @@ class TestOrbitalRotation:
     def test_orbital_rotation_generator(self, phi):
         """Tests that the OrbitalRotation operation calculates the correct generator"""
         op = qml.OrbitalRotation(phi, wires=[0, 1, 2, 3])
-        g, a = op.generator
+        g, a = qml.utils.get_generator(op, return_matrix=True)
 
         res = expm(1j * a * g * phi)
         exp = OrbitalRotation(phi)
diff --git a/tests/test_operation.py b/tests/test_operation.py
index 1102e8586..3011db952 100644
--- a/tests/test_operation.py
+++ b/tests/test_operation.py
@@ -1345,7 +1345,8 @@ class TestOperationDerivative:
         parameters"""
 
         class RotWithGen(qml.Rot):
-            generator = [np.zeros((2, 2)), 1]
+            def generator(self):
+                return qml.Hermitian(np.zeros((2, 2)), wires=self.wires)
 
         op = RotWithGen(0.1, 0.2, 0.3, wires=0)
 
diff --git a/tests/transforms/test_metric_tensor.py b/tests/transforms/test_metric_tensor.py
index 999a20db3..3fbbfd08d 100644
--- a/tests/transforms/test_metric_tensor.py
+++ b/tests/transforms/test_metric_tensor.py
@@ -262,7 +262,7 @@ class TestMetricTensor:
             qml.RX(a[1], wires=0)
             qml.RY(a[0], wires=0)
             qml.CNOT(wires=[0, 1])
-            qml.PhaseShift(b, wires=1)
+            qml.U1(b, wires=1)
             return qml.expval(qml.PauliX(0)), qml.expval(qml.PauliX(1))
 
         circuit = qml.QNode(circuit, dev)
@@ -577,13 +577,13 @@ class TestMetricTensor:
 
         tapes, proc_fn = qml.metric_tensor(tape, approx="block-diag")
         assert len(tapes) == 4
-        assert [len(tape.operations) for tape in tapes] == [2, 4, 5, 6]
+        assert [len(tape.operations) for tape in tapes] == [3, 5, 4, 5]
         assert [len(tape.measurements) for tape in tapes] == [1] * 4
         expected_ops = [
-            [qml.Hadamard, qml.QubitUnitary],
-            [qml.Hadamard, qml.Hadamard, qml.IsingXX, qml.QubitUnitary],
-            [qml.Hadamard, qml.Hadamard, qml.IsingXX, qml.IsingXX, qml.QubitUnitary],
-            [qml.Hadamard, qml.Hadamard, qml.IsingXX, qml.IsingXX, qml.IsingZZ, qml.QubitUnitary],
+            [qml.Hadamard, qml.Hadamard, qml.Hadamard],
+            [qml.Hadamard, qml.Hadamard, qml.IsingXX, qml.Hadamard, qml.Hadamard],
+            [qml.Hadamard, qml.Hadamard, qml.IsingXX, qml.IsingXX],
+            [qml.Hadamard, qml.Hadamard, qml.IsingXX, qml.IsingXX, qml.IsingZZ],
         ]
         assert [[type(op) for op in tape.operations] for tape in tapes] == expected_ops
 
@@ -1157,13 +1157,13 @@ def test_generator_no_expval(monkeypatch):
     """Test exception is raised if subcircuit contains an
     operation with generator object that is not an observable"""
     with monkeypatch.context() as m:
-        m.setattr("pennylane.RX.generator", [qml.RX, 1])
+        m.setattr("pennylane.RX.generator", lambda self: qml.RX(0.1, wires=0))
 
         with qml.tape.QuantumTape() as tape:
             qml.RX(np.array(0.5, requires_grad=True), wires=0)
             qml.expval(qml.PauliX(0))
 
-        with pytest.raises(qml.QuantumFunctionError, match="no corresponding observable"):
+        with pytest.raises(qml.QuantumFunctionError, match="is not an observable"):
             qml.metric_tensor(tape, approx="block-diag")
 
 
@@ -1236,16 +1236,19 @@ def test_error_generator_not_registered(allow_nonunitary, monkeypatch):
             with pytest.raises(ValueError, match="Generator for operation"):
                 qml.metric_tensor(circuit, approx=None, allow_nonunitary=allow_nonunitary)(x, z)
 
+    class RX(qml.RX):
+        def generator(self):
+            return qml.Hadamard(self.wires)
+
     @qml.qnode(dev)
     def circuit(x, z):
-        qml.RX(x, wires="wire1")
+        RX(x, wires="wire1")
         qml.RZ(z, wires="wire1")
         return qml.expval(qml.PauliZ("wire2"))
 
     with monkeypatch.context() as m:
         exp_fn = lambda tape, *args, **kwargs: tape
         m.setattr("pennylane.transforms.metric_tensor.expand_fn", exp_fn)
-        m.setattr("pennylane.RX.generator", [qml.Hadamard, 1.0])
 
         if allow_nonunitary:
             qml.metric_tensor(circuit, approx=None, allow_nonunitary=allow_nonunitary)(x, z)
diff --git a/tests/transforms/test_tape_expand.py b/tests/transforms/test_tape_expand.py
index 7f5e1d61d..ce0eb275c 100644
--- a/tests/transforms/test_tape_expand.py
+++ b/tests/transforms/test_tape_expand.py
@@ -147,7 +147,8 @@ class TestExpandMultipar:
         dev = qml.device("default.qubit", wires=3)
 
         class _CRX(qml.CRX):
-            generator = [None, 1]
+            def generator(self):
+                return None
 
         with qml.tape.QuantumTape() as tape:
             qml.RX(1.5, wires=0)
@@ -205,7 +206,8 @@ class TestExpandNonunitaryGen:
         unitary generators and non-parametric operations is not touched."""
 
         class _PhaseShift(qml.PhaseShift):
-            generator = [None, 1]
+            def generator(self):
+                return None
 
         with qml.tape.JacobianTape() as tape:
             qml.RX(0.2, wires=0)