Repository: PennyLaneAI/pennylane
Pull request #7046 (2025-03-05 04:51:19)

HT gradient - additional modes
**Context:** Hadamard gradients method variants

**Description of the Change:** added support for different modes in hadamard gradient and implemented direct HT, reversed HT and reversed_direct HT, whose details could be found at https://arxiv.org/pdf/2408.05406

**Benefits:** more variants of HT + measurement grouping generates better efficiency for gradient estimation

**Possible Drawbacks:** for reversed methods, that is reversed HT and reversed_direct HT, the current version does not support circuits with multiple measurements.  One will need to modify the processing_fn in how it aggregates the results so support this. For non-reversed methods, we aggregate for each measurement and the code work because the shape is N_G*M  where N_G is the number of terms in the generator and M is the number of measurements. However, in reversed method, the shape is [N_1, ..., N_M], where for the i-th measurement, there are N_i gradient tapes. Different measurement has different number of tapes, and all the tape has only one measurement, so one does not need to aggregate across dimensions, just group neighboring ones would do.

**Related GitHub Issues:**


albi3ro: So for the failing test, we keep a hardcoded list of potential gradient keyword arguments here:

https://github.com/PennyLaneAI/pennylane/blob/7bdf71b44318176e937206e61ccdb47a78e38ce7/pennylane/gradients/gradient_transform.py#L24

This was done as a patch to add validation to the `QNode(qfunc, dev, ... , **gradient_kwargs)` way of passing gradient keyword arguments and catching mistakes.  Gladly, we will be getting rid of that in April, but for now, we will need to add `"mode"` to that hardcoded set.
albi3ro: The documentation check warnings are:
```
/github/workspace/pennylane/workflow/qnode.py:docstring of pennylane.workflow.qnode.QNode:70: WARNING: Bullet list ends without a blank line; unexpected unindent.
/github/workspace/pennylane/workflow/qnode.py:docstring of pennylane.workflow.qnode.QNode:74: WARNING: Bullet list ends without a blank line; unexpected unindent.
/github/workspace/pennylane/workflow/qnode.py:docstring of pennylane.workflow.qnode.QNode:82: WARNING: Bullet list ends without a blank line; unexpected unindent.
/github/workspace/pennylane/workflow/qnode.py:docstring of pennylane.workflow.qnode.QNode:88: WARNING: Bullet list ends without a blank line; unexpected unindent.
/github/workspace/pennylane/workflow/qnode.py:docstring of pennylane.workflow.qnode.qnode:70: WARNING: Bullet list ends without a blank line; unexpected unindent.
/github/workspace/pennylane/workflow/qnode.py:docstring of pennylane.workflow.qnode.qnode:74: WARNING: Bullet list ends without a blank line; unexpected unindent.
/github/workspace/pennylane/workflow/qnode.py:docstring of pennylane.workflow.qnode.qnode:82: WARNING: Bullet list ends without a blank line; unexpected unindent.
/github/workspace/pennylane/workflow/qnode.py:docstring of pennylane.workflow.qnode.qnode:88: WARNING: Bullet list ends without a blank line; unexpected unindent.
```

so something is happening with the qnode docstring.  

I think we should avoid changing the `QNode` docstring for the moment, and update the docstring for `hadamard_grad` with how to select a mode with the qnode:
```
@qml.qnode(dev, diff_method="hadamard", gradient_kwargs={"mode": "reversed"})
def circuit(x):
   ...
```

Then I can do a follow up where add `diff_method="hadamard-reversed"` (etc.) as a qnode option.
albi3ro: just need to add `from typing import Literal` up at the top of the file.
JerryChen97: Yeah we definitely need to `from typing import Literal` on the top to have it
JerryChen97: Sphinx stopped complaining finally :partying_face: 
JerryChen97: @Alex-Preciado @PietropaoloFrisoni I added this PR to 0.41 as well

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index dc5ffb9e6..ff64b9c34 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -4,6 +4,9 @@
 
 <h3>New features since last release</h3>
 
+* Added additional Hadamard gradient modes and `"reversed"`, `"direct"`, and `"reversed-direct"` modes are now available for use with the hadamard gradient.
+  [(#7046)](https://github.com/PennyLaneAI/pennylane/pull/7046)
+
 <h4>Resource-efficient decompositions 🔎</h4>
 
 A new, experimental graph-based decomposition system is now available in PennyLane under the `qml.decomposition` 
@@ -1121,6 +1124,7 @@ Soran Jahangiri,
 Korbinian Kottmann,
 Christina Lee,
 Joseph Lee,
+Dantong Li,
 William Maxwell,
 Anton Naim Ibrahim,
 Lee J. O'Riordan,
diff --git a/pennylane/gradients/gradient_transform.py b/pennylane/gradients/gradient_transform.py
index d939de723..8c9e9d5b1 100644
--- a/pennylane/gradients/gradient_transform.py
+++ b/pennylane/gradients/gradient_transform.py
@@ -34,6 +34,7 @@ SUPPORTED_GRADIENT_KWARGS = {
     "force_order2",
     "gradient_recipes",
     "h",
+    "mode",
     "n",
     "num_directions",
     "num_split_times",
@@ -81,6 +82,21 @@ def assert_no_variance(measurements, transform_name):
         )
 
 
+def assert_no_probability(measurements, transform_name):
+    """Check whether a set of measurements contains a probability measurement
+    raise an error if this is the case.
+
+    Args:
+        measurements (list[MeasurementProcess]): measurements to analyze
+        transform_name (str): Name of the gradient transform that queries the measurements
+    """
+    if any(isinstance(m, ProbabilityMP) for m in measurements):
+        raise ValueError(
+            f"Computing the gradient of probabilities with the {transform_name} "
+            "gradient transform is not supported."
+        )
+
+
 def assert_no_trainable_tape_batching(tape, transform_name):
     """Check whether a tape is broadcasted and raise an error if this is the case.
 
diff --git a/pennylane/gradients/hadamard_gradient.py b/pennylane/gradients/hadamard_gradient.py
old mode 100644
new mode 100755
index fc125d1b4..93392192b
--- a/pennylane/gradients/hadamard_gradient.py
+++ b/pennylane/gradients/hadamard_gradient.py
@@ -17,6 +17,7 @@ of a qubit-based quantum tape.
 """
 from functools import partial
 from itertools import islice
+from typing import Literal
 
 import numpy as np
 
@@ -29,6 +30,7 @@ from .gradient_transform import (
     _contract_qjac_with_cjac,
     _no_trainable_grad,
     _try_zero_grad_from_graph_or_get_grad_method,
+    assert_no_probability,
     assert_no_state_returns,
     assert_no_trainable_tape_batching,
     assert_no_variance,
@@ -59,6 +61,7 @@ def _expand_transform_hadamard(
     argnum=None,
     aux_wire=None,
     device_wires=None,
+    mode: Literal["standard", "reversed", "direct", "reversed-direct"] = "standard",
 ) -> tuple[QuantumScriptBatch, PostprocessingFn]:
     """Expand function to be applied before hadamard gradient."""
     batch, postprocessing = qml.devices.preprocess.decompose(
@@ -94,6 +97,7 @@ def hadamard_grad(
     argnum=None,
     aux_wire=None,
     device_wires=None,
+    mode: Literal["standard", "reversed", "direct", "reversed-direct"] = "standard",
 ) -> tuple[QuantumScriptBatch, PostprocessingFn]:
     r"""Transform a circuit to compute the Hadamard test gradient of all gates
     with respect to their inputs.
@@ -109,6 +113,9 @@ def hadamard_grad(
             the original circuit and ``device_wires``.
         device_wires (pennylane.wires.Wires): Wires of the device that are going to be used for the
             gradient. Facilitates finding a default for ``aux_wire`` if ``aux_wire`` is ``None``.
+        mode (str): Specifies the gradient computation mode. Accepted values are
+            ``"standard"``, ``"reversed"``, ``"direct"``, and ``"reversed-direct"``. The default
+            is ``"standard"``.
 
     Returns:
         qnode (QNode) or tuple[List[QuantumTape], function]:
@@ -270,7 +277,17 @@ def hadamard_grad(
         Array([-0.3875172 , -0.18884787, -0.38355704], dtype=float64)
     """
 
-    transform_name = "Hadamard test"
+    modes = {
+        "standard": ("Hadamard test", _hadamard_test),
+        "reversed": ("Reversed hadamard test", _reversed_hadamard_test),
+        "direct": ("Direct hadamard test", _direct_hadamard_test),
+        "reversed-direct": ("Reversed direct hadamard test", _reversed_direct_hadamard_test),
+    }
+    try:
+        transform_name, gradient_method = modes[mode]
+    except KeyError as exc:
+        raise ValueError(f"Invalid mode: {mode}") from exc
+
     assert_no_state_returns(tape.measurements, transform_name)
     assert_no_variance(tape.measurements, transform_name)
     assert_no_trainable_tape_batching(tape, transform_name)
@@ -278,6 +295,13 @@ def hadamard_grad(
         raise NotImplementedError(
             "hadamard gradient does not support multiple measurements with partitioned shots."
         )
+    if mode in ["reversed", "direct", "reversed-direct"]:
+        assert_no_probability(tape.measurements, transform_name)
+    if mode in ["reversed", "reversed-direct"]:
+        if len(tape.measurements) > 1:
+            raise NotImplementedError(
+                "Reversed hadamard gradient does not support multiple measurements."
+            )
 
     if argnum is None and not tape.trainable_params:
         return _no_trainable_grad(tape)
@@ -285,7 +309,11 @@ def hadamard_grad(
     trainable_param_indices = choose_trainable_param_indices(tape, argnum)
 
     # Validate or get default for aux_wire
-    aux_wire = _get_aux_wire(aux_wire, tape, device_wires)
+    # unless using direct or reversed-direct modes
+
+    aux_wire = (
+        _get_aux_wire(aux_wire, tape, device_wires) if mode in ["standard", "reversed"] else None
+    )
 
     g_tapes = []
     coeffs = []
@@ -303,7 +331,7 @@ def hadamard_grad(
         else:
             # can dispatch between different algorithms here in the future
             # hadamard test, direct hadamard test, reversed, reversed direct, and flexible
-            batch, new_coeffs = _hadamard_test(tape, trainable_param_idx, aux_wire)
+            batch, new_coeffs = gradient_method(tape, trainable_param_idx, aux_wire)
             g_tapes += batch
             coeffs += new_coeffs
             generators_per_parameter.append(len(batch))
@@ -336,27 +364,136 @@ def _hadamard_test(tape, trainable_param_idx, aux_wire) -> tuple[list, list]:
     return new_batch, sub_coeffs
 
 
+def _direct_hadamard_test(tape, trainable_param_idx, aux_wire) -> tuple[list, list]:
+
+    trainable_op, idx, _ = tape.get_operation(trainable_param_idx)
+
+    ops_to_trainable_op = tape.operations[: idx + 1]
+    ops_after_trainable_op = tape.operations[idx + 1 :]
+
+    # Get a generator and coefficients
+    sub_coeffs, generators = _get_pauli_generators(trainable_op)
+
+    measurements = tape.measurements
+
+    new_batch = []
+    new_coeffs = []
+    for idx, gen in enumerate(generators):
+        pos_rot = [qml.evolve(gen, np.pi / 4)]
+        neg_rot = [qml.evolve(gen, -np.pi / 4)]
+        pos_ops = ops_to_trainable_op + pos_rot + ops_after_trainable_op
+        neg_ops = ops_to_trainable_op + neg_rot + ops_after_trainable_op
+
+        pos_tape = qml.tape.QuantumScript(pos_ops, measurements, shots=tape.shots)
+        neg_tape = qml.tape.QuantumScript(neg_ops, measurements, shots=tape.shots)
+        new_batch.append(pos_tape)
+        new_batch.append(neg_tape)
+        new_coeffs.append(-1 / 2 * sub_coeffs[idx])
+        new_coeffs.append(1 / 2 * sub_coeffs[idx])
+    return new_batch, new_coeffs
+
+
+def _reversed_hadamard_test(tape, trainable_param_idx, aux_wire) -> tuple[list, list]:
+
+    trainable_op, idx, _ = tape.get_operation(trainable_param_idx)
+
+    ops_before_trainable_op = tape.operations[:]
+    ops_after_trainable_op = [qml.adjoint(op) for op in reversed(tape.operations[idx + 1 :])]
+
+    # Create measurement with gate generators
+    mp = qml.expval(trainable_op.generator() @ qml.Y(aux_wire))
+    measurements = [mp]
+
+    # Get the observable from tape measurement
+    # Assume there's only one observable in the tape ################ processing function aggregation
+    coeffs, observables = _get_pauli_terms(tape.measurements[0].obs)
+    coeffs = [-c for c in coeffs]
+
+    new_batch = []
+    for obs in observables:
+        ctrl_obs = [qml.ctrl(obs, control=aux_wire)]
+        hadamard = [qml.Hadamard(wires=aux_wire)]
+        ops = ops_before_trainable_op + hadamard + ctrl_obs + hadamard + ops_after_trainable_op
+
+        new_tape = qml.tape.QuantumScript(ops, measurements, shots=tape.shots)
+        new_batch.append(new_tape)
+    return new_batch, coeffs
+
+
+def _reversed_direct_hadamard_test(tape, trainable_param_idx, aux_wire) -> tuple[list, list]:
+
+    trainable_op, idx, _ = tape.get_operation(trainable_param_idx)
+
+    ops_before_trainable_op = tape.operations[:]
+    ops_after_trainable_op = [qml.adjoint(op) for op in reversed(tape.operations[idx + 1 :])]
+
+    # Create measurement with gate generators
+    measurements = [qml.expval(trainable_op.generator())]
+
+    # Get the observable from tape measurement
+    # Assume there's only one observable in the tape ################ processing function aggregation
+    coeffs, observables = _get_pauli_terms(tape.measurements[0].obs)
+
+    new_batch = []
+    new_coeffs = []
+    for idx, obs in enumerate(observables):
+        pos_rot = [qml.evolve(obs, np.pi / 4)]
+        neg_rot = [qml.evolve(obs, -np.pi / 4)]
+        pos_ops = ops_before_trainable_op + pos_rot + ops_after_trainable_op
+        neg_ops = ops_before_trainable_op + neg_rot + ops_after_trainable_op
+
+        pos_tape = qml.tape.QuantumScript(pos_ops, measurements, shots=tape.shots)
+        neg_tape = qml.tape.QuantumScript(neg_ops, measurements, shots=tape.shots)
+        new_batch.append(pos_tape)
+        new_batch.append(neg_tape)
+        new_coeffs.append(1 / 2 * coeffs[idx])
+        new_coeffs.append(-1 / 2 * coeffs[idx])
+    return new_batch, new_coeffs
+
+
 def _new_measurement(mp, aux_wire, all_wires: qml.wires.Wires):
     obs = mp.obs or qml.prod(*(qml.Z(w) for w in mp.wires or all_wires))
     new_obs = qml.simplify(obs @ qml.Y(aux_wire))
     return type(mp)(obs=new_obs)
 
 
+def _get_pauli_terms(op):
+    """Extract the Pauli terms (generators) and their coefficients for an operator.
+
+    If the operator has no pre-computed pauli_rep, the function computes the matrix
+    and performs a Pauli decomposition.
+
+    Parameters:
+        op: The operator (e.g., a Hamiltonian) for which to extract the Pauli terms.
+
+    Returns:
+        The Pauli terms (generators) and their coefficients.
+    """
+    if op.pauli_rep is None:
+        mat = qml.matrix(op, wire_order=op.wires)
+        pauli_rep = qml.pauli_decompose(mat, wire_order=op.wires, pauli=True)
+    else:
+        pauli_rep = op.pauli_rep
+
+    # Remove identity term if present.
+    id_pw = qml.pauli.PauliWord({})
+    if id_pw in pauli_rep:
+        del pauli_rep[qml.pauli.PauliWord({})]
+
+    # qml.PauliZ has no defined terms() behavior
+    return (
+        pauli_rep.operation().terms()
+        if isinstance(pauli_rep.operation(), qml.ops.op_math.Sum)
+        else (1 * pauli_rep.operation()).terms()
+    )
+
+
 def _get_pauli_generators(trainable_op):
     """From a trainable operation, extract the unitary generators and their coefficients.
     Any operator with a generator is supported.
     """
     generator = trainable_op.generator()
-    if generator.pauli_rep is None:
-        mat = qml.matrix(generator, wire_order=generator.wires)
-        pauli_rep = qml.pauli_decompose(mat, wire_order=generator.wires, pauli=True)
-    else:
-        pauli_rep = generator.pauli_rep
-    id_pw = qml.pauli.PauliWord({})
-    if id_pw in pauli_rep:
-        del pauli_rep[qml.pauli.PauliWord({})]  # remove identity term
-    sum_op = pauli_rep.operation()
-    return sum_op.terms()
+    return _get_pauli_terms(generator)
 
 
 def _postprocess_probs(res, measurement, tape):
diff --git a/pennylane/workflow/qnode.py b/pennylane/workflow/qnode.py
index 76275ae19..9c5c248bb 100644
--- a/pennylane/workflow/qnode.py
+++ b/pennylane/workflow/qnode.py
@@ -273,9 +273,10 @@ class QNode:
               rule for all supported quantum operation arguments, with finite-difference
               as a fallback.
 
-            * ``"hadamard"``: Use the analytic hadamard gradient test
-              rule for all supported quantum operation arguments. More info is in the documentation
-              :func:`qml.gradients.hadamard_grad <.gradients.hadamard_grad>`.
+            * ``"hadamard"``: Use the standard analytic hadamard gradient test rule for
+              all supported quantum operation arguments. More info is in the documentation
+              for :func:`qml.gradients.hadamard_grad <.gradients.hadamard_grad>`. Reversed,
+              direct, and reversed-direct modes can be selected via a ``"mode"`` in ``gradient_kwargs``.
 
             * ``"finite-diff"``: Uses numerical finite-differences for all quantum operation
               arguments.


Test patch:
diff --git a/tests/gradients/core/test_hadamard_gradient.py b/tests/gradients/core/test_hadamard_gradient.py
old mode 100644
new mode 100755
index 2474ff4bf..5326789ee
--- a/tests/gradients/core/test_hadamard_gradient.py
+++ b/tests/gradients/core/test_hadamard_gradient.py
@@ -14,7 +14,7 @@
 """
 Tests for the gradients.hadamard_gradient module.
 """
-
+# pylint: disable=use-implicit-booleaness-not-comparison
 import pytest
 
 import pennylane as qml
@@ -105,27 +105,315 @@ def cost12(x):
     return qml.probs(op=qml.Hadamard(0) @ qml.Y(1) @ qml.Y(2))
 
 
-# pylint: disable=too-many-public-methods
-class TestHadamardGrad:
-    """Unit tests for the hadamard_grad function"""
+class TestHadamardValidation:
+    """Test validation of edge cases with the hadamard gradient."""
+
+    @pytest.mark.parametrize("mode", ["invalid-mode"])
+    def test_invalid_mode(self, mode):
+        """Test that a ValueError is raised if an invalid mode is provided."""
+        tape = qml.tape.QuantumScript([qml.RX(0.543, 0), qml.RY(-0.654, 0)], [qml.expval(qml.Z(0))])
+        _match = r"Invalid mode"
+        with pytest.raises(ValueError, match=_match):
+            qml.gradients.hadamard_grad(tape, mode=mode)
 
-    def test_trainable_batched_tape_raises(self):
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
+    def test_trainable_batched_tape_raises(self, mode):
         """Test that an error is raised for a broadcasted/batched tape if the broadcasted
         parameter is differentiated."""
         tape = qml.tape.QuantumScript([qml.RX([0.4, 0.2], 0)], [qml.expval(qml.PauliZ(0))])
-        _match = r"Computing the gradient of broadcasted tapes .* using the Hadamard test gradient"
+        _match = r"Computing the gradient of broadcasted tapes"
         with pytest.raises(NotImplementedError, match=_match):
-            qml.gradients.hadamard_grad(tape)
+            qml.gradients.hadamard_grad(tape, mode=mode)
 
-    def test_error_trainable_obs_probs(self):
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
+    def test_error_trainable_obs_probs(self, mode):
         """Test that there's an error if the observable is trainable with a probs."""
 
         H = qml.numpy.array(2.0) * (qml.X(0) + qml.Y(0))
         tape = qml.tape.QuantumScript([], [qml.probs(op=H)])
         with pytest.raises(ValueError, match="Can only differentiate Hamiltonian coefficients"):
-            qml.gradients.hadamard_grad(tape)
+            qml.gradients.hadamard_grad(tape, mode=mode)
+
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
+    def test_tape_with_partitioned_shots_multiple_measurements_raises(self, mode):
+        """Test that an error is raised with multiple measurements and partitioned shots."""
+        tape = qml.tape.QuantumScript(
+            [qml.RX(0.1, wires=0)],
+            [qml.expval(qml.PauliZ(0)), qml.expval(qml.PauliY(0))],
+            shots=(1000, 10000),
+        )
+        with pytest.raises(NotImplementedError):
+            qml.gradients.hadamard_grad(tape, mode=mode)
+
+    @pytest.mark.parametrize("mode", ["standard", "reversed"])
+    @pytest.mark.parametrize("aux_wire", [qml.wires.Wires(0), qml.wires.Wires(1)])
+    @pytest.mark.parametrize("device_wires", [qml.wires.Wires([0, 1, "aux"])])
+    def test_aux_wire_already_used_wires(self, aux_wire, device_wires, mode):
+        """Test the aux wire is available."""
+        dev = qml.device("default.qubit", wires=device_wires)
+        x = 0.543
+        y = -0.654
+
+        with qml.queuing.AnnotatedQueue() as q:
+            qml.RX(x, wires=[0])
+            qml.RY(y, wires=[1])
+            qml.CNOT(wires=[0, 1])
+            qml.expval(qml.PauliZ(0) @ qml.PauliX(1))
+
+        tape = qml.tape.QuantumScript.from_queue(q)
+
+        _match = "The requested auxiliary wire is already in use by the circuit"
+        with pytest.raises(qml.wires.WireError, match=_match):
+            qml.gradients.hadamard_grad(tape, aux_wire=aux_wire, device_wires=dev.wires, mode=mode)
+
+    @pytest.mark.parametrize("mode", ["standard", "reversed"])
+    @pytest.mark.parametrize("device_wires", [qml.wires.Wires([0, 1, 2])])
+    def test_requested_wire_not_exist(self, device_wires, mode):
+        """Test if the aux wire is not on the device an error is raised."""
+        aux_wire = qml.wires.Wires("aux")
+        dev = qml.device("default.qubit", wires=device_wires)
+        x = 0.543
+        y = -0.654
+
+        with qml.queuing.AnnotatedQueue() as q:
+            qml.RX(x, wires=[0])
+            qml.RY(y, wires=[1])
+            qml.CNOT(wires=[0, 1])
+            qml.expval(qml.PauliZ(0) @ qml.PauliX(1))
+
+        tape = qml.tape.QuantumScript.from_queue(q)
+        _match = "The requested auxiliary wire does not exist on the used device"
+        with pytest.raises(qml.wires.WireError, match=_match):
+            qml.gradients.hadamard_grad(tape, aux_wire=aux_wire, device_wires=dev.wires, mode=mode)
+
+    @pytest.mark.parametrize("mode", ["standard", "reversed"])
+    @pytest.mark.parametrize("aux_wire", [None, qml.wires.Wires(0), qml.wires.Wires(1)])
+    def test_device_not_enough_wires(self, aux_wire, mode):
+        """Test that an error is raised when the device cannot accept an auxiliary wire
+        because it is full."""
+        dev = qml.device("default.qubit", wires=2)
+
+        m = qml.expval(qml.PauliZ(0) @ qml.PauliX(1))
+        tape = qml.tape.QuantumScript([qml.RX(0.543, wires=[0]), qml.RY(-0.654, wires=[1])], [m])
+
+        if aux_wire is None:
+            _match = "The device has no free wire for the auxiliary wire."
+        else:
+            _match = "The requested auxiliary wire is already in use by the circuit."
+        with pytest.raises(qml.wires.WireError, match=_match):
+            qml.gradients.hadamard_grad(tape, aux_wire=aux_wire, device_wires=dev.wires, mode=mode)
+
+    @pytest.mark.parametrize("mode", ["standard", "reversed"])
+    def test_device_wire_does_not_exist(self, mode):
+        """Test that an error is raised when the device cannot accept an auxiliary wire
+        because it does not exist on the device."""
+        aux_wire = qml.wires.Wires("aux")
+        dev = qml.device("default.qubit", wires=2)
+        m = qml.expval(qml.PauliZ(0) @ qml.PauliX(1))
+        tape = qml.tape.QuantumScript([qml.RX(0.543, wires=[0]), qml.RY(-0.654, wires=[1])], [m])
+
+        _match = "The requested auxiliary wire does not exist on the used device."
+        with pytest.raises(qml.wires.WireError, match=_match):
+            qml.gradients.hadamard_grad(tape, aux_wire=aux_wire, device_wires=dev.wires, mode=mode)
+
+    @pytest.mark.parametrize("measurement", [qml.state(), qml.var(qml.Z(0))])
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
+    def test_state_non_differentiable_error(self, mode, measurement):
+        """Test error raised if attempting to differentiate with
+        respect to a state"""
+
+        tape = qml.tape.QuantumScript([qml.RX(0.543, 0), qml.RY(-0.654, 0)], [measurement])
+        _match = r"is not supported"
+        with pytest.raises(ValueError, match=_match):
+            qml.gradients.hadamard_grad(tape, mode=mode)
+
+    @pytest.mark.parametrize("measurement", [qml.probs(wires=[0])])
+    @pytest.mark.parametrize("mode", ["reversed", "direct", "reversed-direct"])
+    def test_no_probs_with_new_modes(self, mode, measurement):
+        """Test that a ValueError is raised if qml.probs is used with reversed, direct, or reversed direct"""
+
+        tape = qml.tape.QuantumScript([qml.RX(0.543, 0), qml.RY(-0.654, 0)], [measurement])
+        _match = r"is not supported"
+        with pytest.raises(ValueError, match=_match):
+            qml.gradients.hadamard_grad(tape, mode=mode)
+
+    @pytest.mark.parametrize("mode", ("reversed", "reversed-direct"))
+    def test_at_most_one_measurement_with_reversed(self, mode):
+        """Test that a ValueError is raised if more than one expectation value is used
+        with reversed and reversed-direct methods."""
+
+        tape = qml.tape.QuantumScript(
+            [qml.RX(0.543, 0), qml.RY(-0.654, 0)], [qml.expval(qml.Z(0)), qml.expval(qml.X(0))]
+        )
+        with pytest.raises(NotImplementedError):
+            qml.gradients.hadamard_grad(tape, mode=mode)
+
+
+class TestDifferentModes:
+    """Unit test the different modes of gradient."""
+
+    def test_reversed_mode(self):
+        """Test that reversed mode hadamard gradient constructs the correct batch of tapes."""
+
+        tape = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1))], [qml.expval(qml.X(0) + qml.Y(0))], shots=50
+        )
+        batch, fn = qml.gradients.hadamard_grad(tape, mode="reversed")
+        assert len(batch) == 2
+
+        expected_H = qml.IsingXY(0.5, wires=(0, 1)).generator() @ qml.Y(2)
+        expected0 = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1)), qml.H(2), qml.CNOT((2, 0)), qml.H(2)],
+            [qml.expval(expected_H)],
+            shots=50,
+        )
+        qml.assert_equal(batch[0], expected0)
+        expected1 = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1)), qml.H(2), qml.CY((2, 0)), qml.H(2)],
+            [qml.expval(expected_H)],
+            shots=50,
+        )
+        qml.assert_equal(batch[1], expected1)
+
+        assert qml.math.allclose(fn((1.0, 2.0)), -6.0)
+
+    def test_direct_mode(self):
+        """Directly test the batch output of applying the direct mode hadamard gradient."""
+        tape = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1))], [qml.expval(qml.X(0) + qml.Y(0))], shots=50
+        )
+        batch, fn = qml.gradients.hadamard_grad(tape, mode="direct")
+        assert len(batch) == 4
+
+        measurement = tape.measurements
+
+        expected0 = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1)), qml.evolve(qml.X(0) @ qml.X(1), np.pi / 4)],
+            measurement,
+            shots=50,
+        )
+        qml.assert_equal(batch[0], expected0)
+        expected1 = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1)), qml.evolve(qml.X(0) @ qml.X(1), -np.pi / 4)],
+            measurement,
+            shots=50,
+        )
+        qml.assert_equal(batch[1], expected1)
+        expected2 = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1)), qml.evolve(qml.Y(0) @ qml.Y(1), np.pi / 4)],
+            measurement,
+            shots=50,
+        )
+        qml.assert_equal(batch[2], expected2)
+        expected3 = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1)), qml.evolve(qml.Y(0) @ qml.Y(1), -np.pi / 4)],
+            measurement,
+            shots=50,
+        )
+        qml.assert_equal(batch[3], expected3)
+
+        assert qml.math.allclose(fn((1.0, 2.0, 3.0, 4.0)), 0.5)
+
+    def test_reversed_direct_mode(self):
+        """Directly test tht batch output of applying the reversed direct mode of hadamard gradient."""
+        tape = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1))], [qml.expval(qml.X(0) + qml.Y(0))], shots=50
+        )
+        batch, fn = qml.gradients.hadamard_grad(tape, mode="reversed-direct")
+        assert len(batch) == 4
+
+        expected_H = qml.IsingXY(0.5, wires=(0, 1)).generator()
+        expected0 = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1)), qml.evolve(qml.X(0), np.pi / 4)],
+            [qml.expval(expected_H)],
+            shots=50,
+        )
+        qml.assert_equal(batch[0], expected0)
+
+        expected1 = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1)), qml.evolve(qml.X(0), -np.pi / 4)],
+            [qml.expval(expected_H)],
+            shots=50,
+        )
+        qml.assert_equal(batch[1], expected1)
+
+        expected2 = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1)), qml.evolve(qml.Y(0), np.pi / 4)],
+            [qml.expval(expected_H)],
+            shots=50,
+        )
+        qml.assert_equal(batch[2], expected2)
 
-    def test_nontrainable_batched_tape(self):
+        expected3 = qml.tape.QuantumScript(
+            [qml.IsingXY(0.5, wires=(0, 1)), qml.evolve(qml.Y(0), -np.pi / 4)],
+            [qml.expval(expected_H)],
+            shots=50,
+        )
+        qml.assert_equal(batch[3], expected3)
+
+        assert qml.math.allclose(fn((1.0, 2.0, 3.0, 4.0)), -2.0)
+
+    def test_reversed_mode_hadamard_obs(self):
+        """Test that reversed mode can handle measuring an observable that does not have a pauli rep."""
+
+        tape = qml.tape.QuantumScript([qml.RX(0.5, 0)], [qml.expval(qml.H(0))])
+        batch, fn = qml.gradients.hadamard_grad(tape, mode="reversed")
+
+        assert len(batch) == 2
+
+        obs_terms = [qml.X(0), qml.Z(0)]
+        H = qml.RX(0.5, 0).generator() @ qml.Y(1)
+        for i, ob in enumerate(obs_terms):
+            expected = qml.tape.QuantumScript(
+                [qml.RX(0.5, 0), qml.H(1), qml.ctrl(ob, 1), qml.H(1)], [qml.expval(H)]
+            )
+            qml.assert_equal(expected.measurements[0].obs, batch[i].measurements[0].obs)
+            qml.assert_equal(batch[i], expected)
+
+        out = fn((1.0, 2.0))
+        expected = -2 / np.sqrt(2) * (1.0 + 2.0)
+        assert qml.math.allclose(out, expected)
+
+    def test_reversed_direct_mode_hadamard_obs(self):
+        """Test that reversed-direct mode can measure an observable that does not have a pauli rep."""
+        # similar to the above test, we just want to directly check the four outputted tapes for
+        # using reversed direct when taking the expectation value of a hadamard
+        tape = qml.tape.QuantumScript([qml.RX(0.5, 0)], [qml.expval(qml.H(0))])
+        batch, fn = qml.gradients.hadamard_grad(tape, mode="reversed-direct")
+
+        assert len(batch) == 4
+
+        obs_terms = [qml.X(0), qml.Z(0)]
+        H = qml.RX(0.5, 0).generator()
+        for i, ob in enumerate(obs_terms):
+            expected0 = qml.tape.QuantumScript(
+                [qml.RX(0.5, 0), qml.evolve(ob, np.pi / 4)], [qml.expval(H)]
+            )
+            qml.assert_equal(expected0.measurements[0].obs, batch[2 * i].measurements[0].obs)
+            qml.assert_equal(batch[2 * i], expected0)
+            expected1 = qml.tape.QuantumScript(
+                [qml.RX(0.5, 0), qml.evolve(ob, -np.pi / 4)], [qml.expval(H)]
+            )
+            qml.assert_equal(expected1.measurements[0].obs, batch[2 * i + 1].measurements[0].obs)
+            qml.assert_equal(batch[2 * i + 1], expected1)
+
+            out = fn((1.0, 2.0, 3.0, 4.0))
+            expected = 1 / np.sqrt(2) * (1.0 - 2.0) + 1 / np.sqrt(2) * (3.0 - 4.0)
+            assert qml.math.allclose(out, expected)
+
+    @pytest.mark.parametrize("mode", ["direct", "reversed-direct"])
+    def test_no_available_work_wire_direct_methods(self, mode):
+        """Test that direct and reversed direct work with no available work wires."""
+        tape = qml.tape.QuantumScript([qml.RX(0.5, 0)], [qml.expval(qml.Z(0))])
+        _ = qml.gradients.hadamard_grad(tape, mode=mode, device_wires=qml.wires.Wires(0))
+
+
+# pylint: disable=too-many-public-methods
+class TestHadamardGrad:
+    """Unit tests for the hadamard_grad function"""
+
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
+    def test_nontrainable_batched_tape(self, mode):
         """Test that no error is raised for a broadcasted/batched tape if the broadcasted
         parameter is not differentiated, and that the results correspond to the stacked
         results of the single-tape derivatives."""
@@ -134,7 +422,7 @@ class TestHadamardGrad:
         tape = qml.tape.QuantumScript(
             [qml.RY(0.6, 0), qml.RX(x, 0)], [qml.expval(qml.PauliZ(0))], trainable_params=[0]
         )
-        batched_tapes, batched_fn = qml.gradients.hadamard_grad(tape)
+        batched_tapes, batched_fn = qml.gradients.hadamard_grad(tape, mode=mode)
         batched_grad = batched_fn(dev.execute(batched_tapes))
         separate_tapes = [
             qml.tape.QuantumScript(
@@ -142,20 +430,10 @@ class TestHadamardGrad:
             )
             for _x in x
         ]
-        separate_tapes_and_fns = [qml.gradients.hadamard_grad(t) for t in separate_tapes]
+        separate_tapes_and_fns = [qml.gradients.hadamard_grad(t, mode=mode) for t in separate_tapes]
         separate_grad = [_fn(dev.execute(_tapes)) for _tapes, _fn in separate_tapes_and_fns]
         assert np.allclose(batched_grad, separate_grad)
 
-    def test_tape_with_partitioned_shots_multiple_measurements_raises(self):
-        """Test that an error is raised with multiple measurements and partitioned shots."""
-        tape = qml.tape.QuantumScript(
-            [qml.RX(0.1, wires=0)],
-            [qml.expval(qml.PauliZ(0)), qml.expval(qml.PauliY(0))],
-            shots=(1000, 10000),
-        )
-        with pytest.raises(NotImplementedError):
-            qml.gradients.hadamard_grad(tape)
-
     @pytest.mark.parametrize("theta", np.linspace(-2 * np.pi, 2 * np.pi, 7))
     @pytest.mark.parametrize("G", [qml.RX, qml.RY, qml.RZ, qml.PhaseShift, qml.U1])
     def test_pauli_rotation_gradient(self, G, theta, tol):
@@ -540,16 +818,19 @@ class TestHadamardGrad:
         (cost9, (2,), tuple),
     ]
 
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
     @pytest.mark.parametrize("cost, exp_shape, exp_type", costs_and_expected_expval_scalar)
-    def test_output_shape_matches_qnode_expval_scalar(self, cost, exp_shape, exp_type):
+    def test_output_shape_matches_qnode_expval_scalar(self, cost, exp_shape, exp_type, mode):
         """Test that the transform output shape matches that of the QNode for
         expectation values and a scalar parameter."""
+        if mode in {"reversed", "reversed-direct"} and exp_shape == (2,):
+            pytest.xfail("reversed modes do not work with more than one output value.")
         dev = qml.device("default.qubit", wires=4)
 
         x = np.array(0.419)
         circuit = qml.QNode(cost, dev)
 
-        res_hadamard = qml.gradients.hadamard_grad(circuit)(x)
+        res_hadamard = qml.gradients.hadamard_grad(circuit, mode=mode)(x)
 
         assert isinstance(res_hadamard, exp_type)
         assert np.array(res_hadamard).shape == exp_shape
@@ -560,16 +841,20 @@ class TestHadamardGrad:
         (cost3, [2, 3], tuple),
     ]
 
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
     @pytest.mark.parametrize("cost, exp_shape, exp_type", costs_and_expected_expval_array)
-    def test_output_shape_matches_qnode_expval_array(self, cost, exp_shape, exp_type):
+    def test_output_shape_matches_qnode_expval_array(self, cost, exp_shape, exp_type, mode):
         """Test that the transform output shape matches that of the QNode for
         expectation values and an array-valued parameter."""
+        if mode in {"reversed", "reversed-direct"} and exp_shape == [2, 3]:
+            pytest.xfail("reversed modes do not work with more than one output value.")
+
         dev = qml.device("default.qubit", wires=4)
 
         x = np.random.rand(3)
         circuit = qml.QNode(cost, dev)
 
-        res_hadamard = qml.gradients.hadamard_grad(circuit)(x)
+        res_hadamard = qml.gradients.hadamard_grad(circuit, mode=mode)(x)
 
         assert isinstance(res_hadamard, exp_type)
         if len(res_hadamard) == 1:
@@ -590,15 +875,18 @@ class TestHadamardGrad:
         (cost12, [8, 3], np.ndarray),
     ]
 
+    @pytest.mark.parametrize(
+        "mode", ["standard"]
+    )  # other modes should not work with probs at the moment
     @pytest.mark.parametrize("cost, exp_shape, exp_type", costs_and_expected_probs)
-    def test_output_shape_matches_qnode_probs(self, cost, exp_shape, exp_type):
+    def test_output_shape_matches_qnode_probs(self, cost, exp_shape, exp_type, mode):
         """Test that the transform output shape matches that of the QNode."""
         dev = qml.device("default.qubit", wires=5)
 
         x = np.random.rand(3)
         circuit = qml.QNode(cost, dev)
 
-        res_hadamard = qml.gradients.hadamard_grad(circuit)(x)
+        res_hadamard = qml.gradients.hadamard_grad(circuit, mode=mode)(x)
         assert isinstance(res_hadamard, exp_type)
         if len(res_hadamard) == 1:
             res_hadamard = res_hadamard[0]
@@ -606,7 +894,7 @@ class TestHadamardGrad:
 
         # Also check on the tape level
         tape = qml.workflow.construct_tape(circuit)(x)
-        tapes, fn = qml.gradients.hadamard_grad(tape)
+        tapes, fn = qml.gradients.hadamard_grad(tape, mode=mode)
         res_hadamard_tape = qml.math.moveaxis(qml.math.stack(fn(dev.execute(tapes))), -2, -1)
 
         for res in [res_hadamard, res_hadamard_tape]:
@@ -642,182 +930,9 @@ class TestHadamardGrad:
 
         assert all(new_tape.shots == tape.shots for new_tape in tapes)
 
-
-class TestHadamardGradEdgeCases:
-    """Test the Hadamard gradient transform and edge cases such as non diff parameters, auxiliary wires, etc..."""
-
-    # pylint:disable=too-many-public-methods
-
-    device_wires = [qml.wires.Wires([0, 1, "aux"])]
-    device_wires_no_aux = [qml.wires.Wires([0, 1, 2])]
-
-    working_wires = [None, qml.wires.Wires("aux"), "aux"]
-    already_used_wires = [qml.wires.Wires(0), qml.wires.Wires(1)]
-
-    @pytest.mark.parametrize("aux_wire", working_wires)
-    @pytest.mark.parametrize("device_wires", device_wires)
-    def test_aux_wire(self, aux_wire, device_wires):
-        """Test the aux wire is available."""
-        # One qubit is added to the device 2 + 1
-        dev = qml.device("default.qubit", wires=device_wires)
-        x = 0.543
-        y = -0.654
-
-        with qml.queuing.AnnotatedQueue() as q:
-            qml.RX(x, wires=[0])
-            qml.RY(y, wires=[1])
-            qml.CNOT(wires=[0, 1])
-            qml.expval(qml.PauliZ(0) @ qml.PauliX(1))
-
-        tape = qml.tape.QuantumScript.from_queue(q)
-
-        tapes, _ = qml.gradients.hadamard_grad(tape, aux_wire=aux_wire, device_wires=dev.wires)
-        assert len(tapes) == 2
-        tapes, _ = qml.gradients.hadamard_grad(tape, aux_wire=aux_wire)
-        assert len(tapes) == 2
-
-    @pytest.mark.parametrize("aux_wire", already_used_wires)
-    @pytest.mark.parametrize("device_wires", device_wires)
-    def test_aux_wire_already_used_wires(self, aux_wire, device_wires):
-        """Test the aux wire is available."""
-        dev = qml.device("default.qubit", wires=device_wires)
-        x = 0.543
-        y = -0.654
-
-        with qml.queuing.AnnotatedQueue() as q:
-            qml.RX(x, wires=[0])
-            qml.RY(y, wires=[1])
-            qml.CNOT(wires=[0, 1])
-            qml.expval(qml.PauliZ(0) @ qml.PauliX(1))
-
-        tape = qml.tape.QuantumScript.from_queue(q)
-
-        _match = "The requested auxiliary wire is already in use by the circuit"
-        with pytest.raises(qml.wires.WireError, match=_match):
-            qml.gradients.hadamard_grad(tape, aux_wire=aux_wire, device_wires=dev.wires)
-
-    @pytest.mark.parametrize("device_wires", device_wires_no_aux)
-    def test_requested_wire_not_exist(self, device_wires):
-        """Test if the aux wire is not on the device an error is raised."""
-        aux_wire = qml.wires.Wires("aux")
-        dev = qml.device("default.qubit", wires=device_wires)
-        x = 0.543
-        y = -0.654
-
-        with qml.queuing.AnnotatedQueue() as q:
-            qml.RX(x, wires=[0])
-            qml.RY(y, wires=[1])
-            qml.CNOT(wires=[0, 1])
-            qml.expval(qml.PauliZ(0) @ qml.PauliX(1))
-
-        tape = qml.tape.QuantumScript.from_queue(q)
-        _match = "The requested auxiliary wire does not exist on the used device"
-        with pytest.raises(qml.wires.WireError, match=_match):
-            qml.gradients.hadamard_grad(tape, aux_wire=aux_wire, device_wires=dev.wires)
-
-    @pytest.mark.parametrize("aux_wire", [None] + already_used_wires)
-    def test_device_not_enough_wires(self, aux_wire):
-        """Test that an error is raised when the device cannot accept an auxiliary wire
-        because it is full."""
-        dev = qml.device("default.qubit", wires=2)
-
-        m = qml.expval(qml.PauliZ(0) @ qml.PauliX(1))
-        tape = qml.tape.QuantumScript([qml.RX(0.543, wires=[0]), qml.RY(-0.654, wires=[1])], [m])
-
-        if aux_wire is None:
-            _match = "The device has no free wire for the auxiliary wire."
-        else:
-            _match = "The requested auxiliary wire is already in use by the circuit."
-        with pytest.raises(qml.wires.WireError, match=_match):
-            qml.gradients.hadamard_grad(tape, aux_wire=aux_wire, device_wires=dev.wires)
-
-    def test_device_wire_does_not_exist(self):
-        """Test that an error is raised when the device cannot accept an auxiliary wire
-        because it does not exist on the device."""
-        aux_wire = qml.wires.Wires("aux")
-        dev = qml.device("default.qubit", wires=2)
-        m = qml.expval(qml.PauliZ(0) @ qml.PauliX(1))
-        tape = qml.tape.QuantumScript([qml.RX(0.543, wires=[0]), qml.RY(-0.654, wires=[1])], [m])
-
-        _match = "The requested auxiliary wire does not exist on the used device."
-        with pytest.raises(qml.wires.WireError, match=_match):
-            qml.gradients.hadamard_grad(tape, aux_wire=aux_wire, device_wires=dev.wires)
-
-    def test_empty_circuit(self):
-        """Test that an empty circuit works correctly"""
-        with qml.queuing.AnnotatedQueue() as q:
-            qml.expval(qml.PauliZ(0))
-
-        tape = qml.tape.QuantumScript.from_queue(q)
-        with pytest.warns(UserWarning, match="gradient of a tape with no trainable parameters"):
-            tapes, _ = qml.gradients.hadamard_grad(tape)
-        assert not tapes
-
-    def test_all_parameters_independent(self):
-        """Test that a circuit where all parameters do not affect the output"""
-        with qml.queuing.AnnotatedQueue() as q:
-            qml.RX(0.4, wires=0)
-            qml.expval(qml.PauliZ(1))
-
-        tape = qml.tape.QuantumScript.from_queue(q)
-        tapes, _ = qml.gradients.hadamard_grad(tape)
-        assert not tapes
-
-    def test_state_non_differentiable_error(self):
-        """Test error raised if attempting to differentiate with
-        respect to a state"""
-
-        with qml.queuing.AnnotatedQueue() as q:
-            qml.RX(0.543, wires=[0])
-            qml.RY(-0.654, wires=[1])
-            qml.state()
-
-        tape = qml.tape.QuantumScript.from_queue(q)
-        _match = r"return the state with the Hadamard test gradient transform"
-        with pytest.raises(ValueError, match=_match):
-            qml.gradients.hadamard_grad(tape)
-
-    def test_variance_non_differentiable_error(self):
-        """Test error raised if attempting to differentiate with respect to a variance"""
-
-        with qml.queuing.AnnotatedQueue() as q:
-            qml.RX(0.543, wires=[0])
-            qml.RY(-0.654, wires=[1])
-            qml.var(qml.PauliZ(wires=0))
-
-        tape = qml.tape.QuantumScript.from_queue(q)
-        with pytest.raises(
-            ValueError,
-            match=(
-                r"Computing the gradient of variances with the Hadamard test "
-                "gradient transform is not supported."
-            ),
-        ):
-            qml.gradients.hadamard_grad(tape)
-
-    def test_independent_parameter(self):
-        """Test that an independent parameter is skipped
-        during the Jacobian computation."""
-
-        with qml.queuing.AnnotatedQueue() as q:
-            qml.RX(0.543, wires=[0])
-            qml.RY(-0.654, wires=[1])  # does not have any impact on the expval
-            qml.expval(qml.PauliZ(0))
-
-        tape = qml.tape.QuantumScript.from_queue(q)
-        dev = qml.device("default.qubit", wires=3)
-
-        res_hadamard, tapes = grad_fn(tape, dev)
-
-        assert len(tapes) == 1
-
-        assert isinstance(res_hadamard, (tuple, list))
-        assert len(res_hadamard) == 2
-        assert res_hadamard[0].shape == ()
-        assert res_hadamard[1].shape == ()
-
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
     @pytest.mark.autograd
-    def test_no_trainable_params_qnode_autograd(self):
+    def test_no_trainable_params_qnode_autograd(self, mode):
         """Test that the correct ouput and warning is generated in the absence of any trainable
         parameters"""
         dev = qml.device("default.qubit", wires=2)
@@ -830,11 +945,12 @@ class TestHadamardGradEdgeCases:
 
         weights = [0.1, 0.2]
         with pytest.raises(qml.QuantumFunctionError, match="No trainable parameters."):
-            qml.gradients.hadamard_grad(circuit)(weights)
+            qml.gradients.hadamard_grad(circuit, mode=mode)(weights)
 
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
     @pytest.mark.torch
-    def test_no_trainable_params_qnode_torch(self):
-        """Test that the correct ouput and warning is generated in the absence of any trainable
+    def test_no_trainable_params_qnode_torch(self, mode):
+        """Test that the correct output and warning is generated in the absence of any trainable
         parameters"""
         dev = qml.device("default.qubit", wires=2)
 
@@ -846,10 +962,11 @@ class TestHadamardGradEdgeCases:
 
         weights = [0.1, 0.2]
         with pytest.raises(qml.QuantumFunctionError, match="No trainable parameters."):
-            qml.gradients.hadamard_grad(circuit)(weights)
+            qml.gradients.hadamard_grad(circuit, mode=mode)(weights)
 
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
     @pytest.mark.tf
-    def test_no_trainable_params_qnode_tf(self):
+    def test_no_trainable_params_qnode_tf(self, mode):
         """Test that the correct ouput and warning is generated in the absence of any trainable
         parameters"""
         dev = qml.device("default.qubit", wires=2)
@@ -862,10 +979,11 @@ class TestHadamardGradEdgeCases:
 
         weights = [0.1, 0.2]
         with pytest.raises(qml.QuantumFunctionError, match="No trainable parameters."):
-            qml.gradients.hadamard_grad(circuit)(weights)
+            qml.gradients.hadamard_grad(circuit, mode=mode)(weights)
 
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
     @pytest.mark.jax
-    def test_no_trainable_params_qnode_jax(self):
+    def test_no_trainable_params_qnode_jax(self, mode):
         """Test that the correct ouput and warning is generated in the absence of any trainable
         parameters"""
         dev = qml.device("default.qubit", wires=2)
@@ -878,9 +996,10 @@ class TestHadamardGradEdgeCases:
 
         weights = [0.1, 0.2]
         with pytest.raises(qml.QuantumFunctionError, match="No trainable parameters."):
-            qml.gradients.hadamard_grad(circuit)(weights)
+            qml.gradients.hadamard_grad(circuit, mode=mode)(weights)
 
-    def test_no_trainable_params_tape(self):
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
+    def test_no_trainable_params_tape(self, mode):
         """Test that the correct ouput and warning is generated in the absence of any trainable
         parameters"""
         dev = qml.device("default.qubit", wires=2)
@@ -896,13 +1015,93 @@ class TestHadamardGradEdgeCases:
         tape.trainable_params = []
 
         with pytest.warns(UserWarning, match="gradient of a tape with no trainable parameters"):
-            g_tapes, post_processing = qml.gradients.hadamard_grad(tape)
+            g_tapes, post_processing = qml.gradients.hadamard_grad(tape, mode=mode)
         res_hadamard = post_processing(qml.execute(g_tapes, dev, None))
 
         assert g_tapes == []
         assert isinstance(res_hadamard, np.ndarray)
         assert res_hadamard.shape == (0,)
 
+
+class TestHadamardGradEdgeCases:
+    """Test the Hadamard gradient transform and edge cases such as non diff parameters, auxiliary wires, etc..."""
+
+    # pylint:disable=too-many-public-methods
+
+    device_wires = [qml.wires.Wires([0, 1, "aux"])]
+    device_wires_no_aux = [qml.wires.Wires([0, 1, 2])]
+
+    working_wires = [None, qml.wires.Wires("aux"), "aux"]
+    already_used_wires = [qml.wires.Wires(0), qml.wires.Wires(1)]
+
+    @pytest.mark.parametrize("mode", ["standard"])  # test is fairly standard specific
+    @pytest.mark.parametrize("aux_wire", working_wires)
+    @pytest.mark.parametrize("device_wires", device_wires)
+    def test_aux_wire(self, aux_wire, device_wires, mode):
+        """Test the aux wire is available."""
+        # One qubit is added to the device 2 + 1
+        dev = qml.device("default.qubit", wires=device_wires)
+        x = 0.543
+        y = -0.654
+
+        with qml.queuing.AnnotatedQueue() as q:
+            qml.RX(x, wires=[0])
+            qml.RY(y, wires=[1])
+            qml.CNOT(wires=[0, 1])
+            qml.expval(qml.PauliZ(0) @ qml.PauliX(1))
+
+        tape = qml.tape.QuantumScript.from_queue(q)
+
+        tapes, _ = qml.gradients.hadamard_grad(
+            tape, aux_wire=aux_wire, device_wires=dev.wires, mode=mode
+        )
+        assert len(tapes) == 2
+        tapes, _ = qml.gradients.hadamard_grad(tape, aux_wire=aux_wire, mode=mode)
+        assert len(tapes) == 2
+
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
+    def test_empty_circuit(self, mode):
+        """Test that an empty circuit works correctly"""
+        with qml.queuing.AnnotatedQueue() as q:
+            qml.expval(qml.PauliZ(0))
+
+        tape = qml.tape.QuantumScript.from_queue(q)
+        with pytest.warns(UserWarning, match="gradient of a tape with no trainable parameters"):
+            tapes, _ = qml.gradients.hadamard_grad(tape, mode=mode)
+        assert not tapes
+
+    @pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
+    def test_all_parameters_independent(self, mode):
+        """Test that a circuit where all parameters do not affect the output"""
+        with qml.queuing.AnnotatedQueue() as q:
+            qml.RX(0.4, wires=0)
+            qml.expval(qml.PauliZ(1))
+
+        tape = qml.tape.QuantumScript.from_queue(q)
+        tapes, _ = qml.gradients.hadamard_grad(tape, mode=mode)
+        assert not tapes
+
+    def test_independent_parameter(self):
+        """Test that an independent parameter is skipped
+        during the Jacobian computation."""
+
+        with qml.queuing.AnnotatedQueue() as q:
+            qml.RX(0.543, wires=[0])
+            qml.RY(-0.654, wires=[1])  # does not have any impact on the expval
+            qml.expval(qml.PauliZ(0))
+
+        tape = qml.tape.QuantumScript.from_queue(q)
+        dev = qml.device("default.qubit", wires=3)
+
+        res_hadamard, tapes = grad_fn(tape, dev)
+
+        assert len(tapes) == 1
+
+        assert isinstance(res_hadamard, (tuple, list))
+        assert len(res_hadamard) == 2
+        assert res_hadamard[0].shape == ()
+        assert res_hadamard[1].shape == ()
+
     def test_no_trainable_params_multiple_return_tape(self):
         """Test that the correct ouput and warning is generated in the absence of any trainable
         parameters with multiple returns."""
@@ -1013,8 +1212,9 @@ class TestHadamardGradEdgeCases:
         assert res_hadamard[1][2].shape == (4,)
         assert np.allclose(res_hadamard[1][2], 0)
 
+    @pytest.mark.parametrize("mode", ["standard"])
     @pytest.mark.parametrize("prefactor", [1.0, 2.0])
-    def test_all_zero_diff_methods(self, prefactor):
+    def test_all_zero_diff_methods(self, prefactor, mode):
         """Test that the transform works correctly when the diff method for every parameter is
         identified to be 0, and that no tapes were generated."""
         dev = qml.device("default.qubit", wires=4)
@@ -1026,19 +1226,19 @@ class TestHadamardGradEdgeCases:
 
         params = np.array([0.5, 0.5, 0.5], requires_grad=True)
 
-        result = qml.gradients.hadamard_grad(circuit)(params)
+        result = qml.gradients.hadamard_grad(circuit, mode=mode)(params)
 
         assert isinstance(result, np.ndarray)
-        print(result)
         assert result.shape == (4, 3)
         assert np.allclose(result, 0)
 
 
+@pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
 class TestHadamardTestGradDiff:
     """Test that the transform is differentiable"""
 
     @pytest.mark.autograd
-    def test_autograd(self):
+    def test_autograd(self, mode):
         """Tests that the output of the hadamard gradient transform
         can be differentiated using autograd, yielding second derivatives."""
         dev = qml.device("default.qubit", wires=3)
@@ -1053,7 +1253,7 @@ class TestHadamardTestGradDiff:
 
             tape = qml.tape.QuantumScript.from_queue(q)
             tape.trainable_params = {0, 1}
-            tapes, fn = qml.gradients.hadamard_grad(tape)
+            tapes, fn = qml.gradients.hadamard_grad(tape, mode=mode)
             jac = fn(dev.execute(tapes))
             return qml.math.stack(jac)
 
@@ -1075,7 +1275,7 @@ class TestHadamardTestGradDiff:
         assert np.allclose(res_hadamard, res_param_shift)
 
     @pytest.mark.tf
-    def test_tf(self):
+    def test_tf(self, mode):
         """Tests that the output of the hadamard gradient transform
         can be differentiated using TF, yielding second derivatives."""
         import tensorflow as tf
@@ -1092,7 +1292,7 @@ class TestHadamardTestGradDiff:
 
             tape = qml.tape.QuantumScript.from_queue(q)
             tape.trainable_params = {0, 1}
-            tapes, fn = qml.gradients.hadamard_grad(tape)
+            tapes, fn = qml.gradients.hadamard_grad(tape, mode=mode)
             jac_h = fn(dev.execute(tapes))
             jac_h = qml.math.stack(jac_h)
 
@@ -1115,7 +1315,7 @@ class TestHadamardTestGradDiff:
         assert np.allclose(res_hadamard, res_param_shift)
 
     @pytest.mark.torch
-    def test_torch(self):
+    def test_torch(self, mode):
         """Tests that the output of the hadamard gradient transform
         can be differentiated using Torch, yielding second derivatives."""
         import torch
@@ -1131,7 +1331,7 @@ class TestHadamardTestGradDiff:
                 qml.expval(qml.PauliZ(0) @ qml.PauliX(1))
 
             tape = qml.tape.QuantumScript.from_queue(q)
-            tapes, fn = qml.gradients.hadamard_grad(tape)
+            tapes, fn = qml.gradients.hadamard_grad(tape, mode=mode)
 
             return tuple(fn(dev.execute(tapes)))
 
@@ -1154,7 +1354,7 @@ class TestHadamardTestGradDiff:
         assert np.allclose(res_hadamard[1].detach(), res_param_shift[1].detach())
 
     @pytest.mark.jax
-    def test_jax(self):
+    def test_jax(self, mode):
         """Tests that the output of the hadamard gradient transform
         can be differentiated using JAX, yielding second derivatives."""
         import jax
@@ -1172,7 +1372,7 @@ class TestHadamardTestGradDiff:
 
             tape = qml.tape.QuantumScript.from_queue(q)
             tape.trainable_params = {0, 1}
-            tapes, fn = qml.gradients.hadamard_grad(tape)
+            tapes, fn = qml.gradients.hadamard_grad(tape, mode=mode)
 
             return fn(dev.execute(tapes))
 
@@ -1185,7 +1385,7 @@ class TestHadamardTestGradDiff:
 
             tape = qml.tape.QuantumScript.from_queue(q)
             tape.trainable_params = {0, 1}
-            tapes, fn = qml.gradients.hadamard_grad(tape)
+            tapes, fn = qml.gradients.hadamard_grad(tape, mode=mode)
 
             return fn(dev.execute(tapes))
 
@@ -1196,6 +1396,7 @@ class TestHadamardTestGradDiff:
 
 @pytest.mark.parametrize("argnums", [[0], [1], [0, 1]])
 @pytest.mark.parametrize("interface", ["jax"])
+@pytest.mark.parametrize("mode", ["standard", "reversed", "direct", "reversed-direct"])
 @pytest.mark.jax
 class TestJaxArgnums:
     """Class to test the integration of argnums (Jax) and the hadamard grad transform."""
@@ -1203,7 +1404,7 @@ class TestJaxArgnums:
     expected_jacs = []
     interfaces = ["auto", "jax"]
 
-    def test_argnum_error(self, argnums, interface):
+    def test_argnum_error(self, argnums, interface, mode):
         """Test that giving argnum to Jax, raises an error."""
         import jax
 
@@ -1223,9 +1424,9 @@ class TestJaxArgnums:
             qml.QuantumFunctionError,
             match="argnum does not work with the Jax interface. You should use argnums instead.",
         ):
-            qml.gradients.hadamard_grad(circuit, argnum=argnums)(x, y)
+            qml.gradients.hadamard_grad(circuit, argnum=argnums, mode=mode)(x, y)
 
-    def test_single_expectation_value(self, argnums, interface):
+    def test_single_expectation_value(self, argnums, interface, mode):
         """Test for single expectation value."""
         import jax
 
@@ -1241,7 +1442,7 @@ class TestJaxArgnums:
         x = jax.numpy.array([0.543, 0.2])
         y = jax.numpy.array(-0.654)
 
-        res = qml.gradients.hadamard_grad(circuit, argnums=argnums)(x, y)
+        res = qml.gradients.hadamard_grad(circuit, argnums=argnums, mode=mode)(x, y)
 
         expected_0 = np.array([-np.sin(y) * np.sin(x[0]), 0])
         expected_1 = np.array(np.cos(y) * np.cos(x[0]))
@@ -1254,8 +1455,10 @@ class TestJaxArgnums:
             assert np.allclose(res[0], expected_0)
             assert np.allclose(res[1], expected_1)
 
-    def test_multi_expectation_values(self, argnums, interface):
+    def test_multi_expectation_values(self, argnums, interface, mode):
         """Test for multiple expectation values."""
+        if mode in {"reversed", "reversed-direct"}:
+            pytest.skip("these methods do not work with more than one observable.")
         import jax
 
         dev = qml.device("default.qubit", wires=3)
@@ -1270,7 +1473,7 @@ class TestJaxArgnums:
         x = jax.numpy.array([0.543, 0.2])
         y = jax.numpy.array(-0.654)
 
-        res = qml.gradients.hadamard_grad(circuit, argnums=argnums)(x, y)
+        res = qml.gradients.hadamard_grad(circuit, argnums=argnums, mode=mode)(x, y)
 
         expected_0 = np.array([[-np.sin(x[0]), 0.0], [0.0, 0.0]])
         expected_1 = np.array([0, np.cos(y)])
@@ -1286,7 +1489,7 @@ class TestJaxArgnums:
             assert np.allclose(res[1][0], expected_1[0])
             assert np.allclose(res[1][1], expected_1[1])
 
-    def test_hessian(self, argnums, interface):
+    def test_hessian(self, argnums, interface, mode):
         """Test for hessian."""
         import jax
 
@@ -1302,9 +1505,9 @@ class TestJaxArgnums:
         x = jax.numpy.array([0.543, -0.654])
         y = jax.numpy.array(-0.123)
 
-        res = jax.jacobian(qml.gradients.hadamard_grad(circuit, argnums=argnums), argnums=argnums)(
-            x, y
-        )
+        res = jax.jacobian(
+            qml.gradients.hadamard_grad(circuit, argnums=argnums, mode=mode), argnums=argnums
+        )(x, y)
         res_expected = jax.hessian(circuit, argnums=argnums)(x, y)
 
         if len(argnums) == 1: