From 2147813c93a487c992cdee13bb8b09a2b64ced5a Mon Sep 17 00:00:00 2001
From: Pietropaolo Frisoni <pietropaolo.frisoni@xanadu.ai>
Date: Tue, 14 Jan 2025 16:41:59 -0500
Subject: [PATCH 1/5] Final sync to master (#6831)

Final manual sync PR from the `rc` branch to `master`

---------

Co-authored-by: Christina Lee <christina@xanadu.ai>
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: Mudit Pandey <mudit.pandey@xanadu.ai>
Co-authored-by: ringo-but-quantum <github-ringo-but-quantum@xanadu.ai>
Co-authored-by: Andrija Paurevic <46359773+andrijapau@users.noreply.github.com>
Co-authored-by: Austin Huang <65315367+austingmhuang@users.noreply.github.com>
Co-authored-by: lillian542 <38584660+lillian542@users.noreply.github.com>
Co-authored-by: Astral Cai <astral.cai@xanadu.ai>
Co-authored-by: Yushao Chen (Jerry) <chenys13@outlook.com>
Co-authored-by: Alex Preciado <alex.preciado@xanadu.ai>
Co-authored-by: Korbinian Kottmann <43949391+Qottmann@users.noreply.github.com>
Co-authored-by: dwierichs <david.wierichs@xanadu.ai>
Co-authored-by: Isaac De Vlugt <34751083+isaacdevlugt@users.noreply.github.com>
Co-authored-by: ANTH0NY <39093564+AntonNI8@users.noreply.github.com>
Co-authored-by: AntonNI8 <anton.naimibrahim@xanadu.ai>
Co-authored-by: Diego <diego_guala@hotmail.com>
Co-authored-by: Diego <67476785+DSGuala@users.noreply.github.com>
Co-authored-by: Jay Soni <jbsoni@uwaterloo.ca>
Co-authored-by: qottmann <korbinian.kottmann@gmail.com>
Co-authored-by: Josh Izaac <josh146@gmail.com>
Co-authored-by: Ashish Kanwar Singh <104938869+ashishks0522@users.noreply.github.com>
---
 doc/releases/changelog-0.40.0.md | 2 +-
 pennylane/compiler/compiler.py   | 2 +-
 setup.py                         | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/releases/changelog-0.40.0.md b/doc/releases/changelog-0.40.0.md
index 39e5b6ecc92..bb9cb7cd3bc 100644
--- a/doc/releases/changelog-0.40.0.md
+++ b/doc/releases/changelog-0.40.0.md
@@ -1170,4 +1170,4 @@ Anton Naim Ibrahim,
 Andrija Paurevic,
 Justin Pickering,
 Jay Soni,
-David Wierichs.
+David Wierichs.
\ No newline at end of file
diff --git a/pennylane/compiler/compiler.py b/pennylane/compiler/compiler.py
index cbb38243d83..6c38bb3d203 100644
--- a/pennylane/compiler/compiler.py
+++ b/pennylane/compiler/compiler.py
@@ -23,7 +23,7 @@
 
 from packaging.version import Version
 
-PL_CATALYST_MIN_VERSION = Version("0.9.0")
+PL_CATALYST_MIN_VERSION = Version("0.10.0")
 
 
 class CompileError(Exception):
diff --git a/setup.py b/setup.py
index a6ce2fb7a08..47d6484df65 100644
--- a/setup.py
+++ b/setup.py
@@ -30,7 +30,7 @@
     "appdirs",
     "autoray>=0.6.11",
     "cachetools",
-    "pennylane-lightning>=0.39",
+    "pennylane-lightning>=0.40",
     "requests",
     "typing_extensions",
     "packaging",

From 4fc894e0598d6f59194776039da64224a0adf0cc Mon Sep 17 00:00:00 2001
From: ringo-but-quantum <github-ringo-but-quantum@xanadu.ai>
Date: Wed, 15 Jan 2025 09:51:51 +0000
Subject: [PATCH 2/5] [no ci] bump nightly version

---
 pennylane/_version.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pennylane/_version.py b/pennylane/_version.py
index e13f6bf4d72..1b991cac3ad 100644
--- a/pennylane/_version.py
+++ b/pennylane/_version.py
@@ -16,4 +16,4 @@
 Version number (major.minor.patch[-label])
 """
 
-__version__ = "0.41.0-dev6"
+__version__ = "0.41.0-dev7"

From 3ce6690546d89e56d3ec941329c4d0b5f124fbd6 Mon Sep 17 00:00:00 2001
From: "Yushao Chen (Jerry)" <chenys13@outlook.com>
Date: Wed, 15 Jan 2025 11:59:28 -0500
Subject: [PATCH 3/5] Remove qsvt legacy (#6827)

**Context:**
The ``qml.qsvt_legacy`` function has been removed.
Instead, use ``qml.qsvt``. The new functionality takes an input
polynomial instead of angles.

**Description of the Change:**
removed
 - `qsvt_legacy` along with its tests
- `_qsp_to_qsvt` which is a private helper strictly only used by
`qsvt_legacy`

**Benefits:**

**Possible Drawbacks:**

**Related GitHub Issues:**
[sc-82148]
---
 doc/development/deprecations.rst              |  12 +-
 doc/releases/changelog-dev.md                 |   7 +-
 pennylane/templates/subroutines/__init__.py   |   2 +-
 pennylane/templates/subroutines/qsvt.py       | 159 +----------
 tests/templates/test_subroutines/test_qsvt.py | 255 ------------------
 5 files changed, 14 insertions(+), 421 deletions(-)

diff --git a/doc/development/deprecations.rst b/doc/development/deprecations.rst
index fd1c5051e25..2c4269de6c4 100644
--- a/doc/development/deprecations.rst
+++ b/doc/development/deprecations.rst
@@ -9,12 +9,6 @@ deprecations are listed below.
 Pending deprecations
 --------------------
 
-* The ``qsvt_legacy`` function has been deprecated.
-  Instead, use ``qml.qsvt``. The new functionality takes an input polynomial instead of angles.
-
-  - Deprecated in v0.40
-  - Will be removed in v0.41
-
 * The ``tape`` and ``qtape`` properties of ``QNode`` have been deprecated. 
   Instead, use the ``qml.workflow.construct_tape`` function.
 
@@ -90,6 +84,12 @@ for details on how to port your legacy code to the new system. The following fun
 Completed deprecation cycles
 ----------------------------
 
+* The ``qml.qsvt_legacy`` function has been removed.
+  Instead, use ``qml.qsvt``. The new functionality takes an input polynomial instead of angles.
+
+  - Deprecated in v0.40
+  - Removed in v0.41
+
 * The ``qml.qinfo`` module has been removed. Please see the respective functions in the ``qml.math`` and ``qml.measurements``
   modules instead.
 
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 0f9e46664b9..d99bdb3383d 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -8,6 +8,9 @@
 
 <h3>Breaking changes 💔</h3>
 
+* Removed method `qsvt_legacy` along with its private helper `_qsp_to_qsvt`
+  [(#6827)](https://github.com/PennyLaneAI/pennylane/pull/6827)
+
 <h3>Deprecations 👋</h3>
 
 <h3>Documentation 📝</h3>
@@ -20,4 +23,6 @@
 <h3>Contributors ✍️</h3>
 
 This release contains contributions from (in alphabetical order):
-Diksha Dhawan
\ No newline at end of file
+
+Yushao Chen,
+Diksha Dhawan,
diff --git a/pennylane/templates/subroutines/__init__.py b/pennylane/templates/subroutines/__init__.py
index 517c86fe003..b3960b17341 100644
--- a/pennylane/templates/subroutines/__init__.py
+++ b/pennylane/templates/subroutines/__init__.py
@@ -33,7 +33,7 @@
 from .hilbert_schmidt import HilbertSchmidt, LocalHilbertSchmidt
 from .flip_sign import FlipSign
 from .basis_rotation import BasisRotation
-from .qsvt import poly_to_angles, QSVT, qsvt, qsvt_legacy, transform_angles
+from .qsvt import poly_to_angles, QSVT, qsvt, transform_angles
 from .select import Select
 from .qdrift import QDrift
 from .controlled_sequence import ControlledSequence
diff --git a/pennylane/templates/subroutines/qsvt.py b/pennylane/templates/subroutines/qsvt.py
index 129e6d930a1..2898de207f0 100644
--- a/pennylane/templates/subroutines/qsvt.py
+++ b/pennylane/templates/subroutines/qsvt.py
@@ -16,7 +16,6 @@
 """
 # pylint: disable=too-many-arguments
 import copy
-import warnings
 from typing import Literal
 
 import numpy as np
@@ -24,147 +23,12 @@
 
 import pennylane as qml
 from pennylane.operation import AnyWires, Operation
-from pennylane.ops import BlockEncode, PCPhase
 from pennylane.ops.op_math import adjoint
 from pennylane.queuing import QueuingManager
 from pennylane.wires import Wires
 
 
-def qsvt_legacy(A, angles, wires, convention=None):
-    r"""Implements the
-    `quantum singular value transformation <https://arxiv.org/abs/1806.01838>`__ (QSVT) circuit.
-
-    .. warning::
-
-        The ``qsvt_legacy`` function has been deprecated.
-        Instead, use :func:`~pennylane.qsvt`. The new functionality takes an input polynomial instead of angles.
-
-    .. note ::
-
-        :class:`~.BlockEncode` and :class:`~.PCPhase` used in this implementation of QSVT
-        are matrix-based operators and well-suited for simulators.
-        To implement QSVT with user-defined circuits for the block encoding and
-        projector-controlled phase shifts, use the :class:`~.QSVT` template.
-
-    Given a matrix :math:`A`, and a list of angles :math:`\vec{\phi}`, this function applies a
-    circuit for the quantum singular value transformation using :class:`~.BlockEncode` and
-    :class:`~.PCPhase`.
-
-    When the number of angles is even (:math:`d` is odd), the QSVT circuit is defined as:
-
-    .. math::
-
-        U_{QSVT} = \tilde{\Pi}_{\phi_1}U\left[\prod^{(d-1)/2}_{k=1}\Pi_{\phi_{2k}}U^\dagger
-        \tilde{\Pi}_{\phi_{2k+1}}U\right]\Pi_{\phi_{d+1}},
-
-
-    and when the number of angles is odd (:math:`d` is even):
-
-    .. math::
-
-        U_{QSVT} = \left[\prod^{d/2}_{k=1}\Pi_{\phi_{2k-1}}U^\dagger\tilde{\Pi}_{\phi_{2k}}U\right]
-        \Pi_{\phi_{d+1}}.
-
-    Here, :math:`U` denotes a block encoding of :math:`A` via :class:`~.BlockEncode` and
-    :math:`\Pi_\phi` denotes a projector-controlled phase shift with angle :math:`\phi`
-    via :class:`~.PCPhase`.
-
-    This circuit applies a polynomial transformation (:math:`Poly^{SV}`) to the singular values of
-    the block encoded matrix:
-
-    .. math::
-
-        \begin{align}
-             U_{QSVT}(A, \phi) &=
-             \begin{bmatrix}
-                Poly^{SV}(A) & \cdot \\
-                \cdot & \cdot
-            \end{bmatrix}.
-        \end{align}
-
-    The polynomial transformation is determined by a combination of the block encoding and choice of angles,
-    :math:`\vec{\phi}`. The convention used by :class:`~.BlockEncode` is commonly refered to as the
-    reflection convention or :math:`R` convention. Another equivalent convention for the block encoding is
-    the :math:`Wx` or rotation convention.
-
-    Depending on the choice of convention for blockencoding, the same phase angles will produce different
-    polynomial transformations. We provide the functionality to swap between blockencoding conventions and
-    to transform the phase angles accordingly using the :code:`convention` keyword argument.
-
-    .. seealso::
-
-        :class:`~.QSVT` and `A Grand Unification of Quantum Algorithms <https://arxiv.org/pdf/2105.02859.pdf>`_.
-
-    Args:
-        A (tensor_like): the general :math:`(n \times m)` matrix to be encoded
-        angles (tensor_like): a list of angles by which to shift to obtain the desired polynomial
-        wires (Iterable[int, str], Wires): the wires the template acts on
-        convention (string): can be set to ``"Wx"`` to convert quantum signal processing angles in the
-            `Wx` convention to QSVT angles.
-
-    **Example**
-
-    To implement QSVT in a circuit, we can use the following method:
-
-    >>> dev = qml.device("default.qubit", wires=2)
-    >>> A = np.array([[0.1, 0.2], [0.3, 0.4]])
-    >>> angles = np.array([0.1, 0.2, 0.3])
-    >>> @qml.qnode(dev)
-    ... def example_circuit(A):
-    ...     qml.qsvt_legacy(A, angles, wires=[0, 1])
-    ...     return qml.expval(qml.Z(0))
-
-    The resulting circuit implements QSVT.
-
-    >>> print(qml.draw(example_circuit)(A))
-    0: ─╭QSVT─┤  <Z>
-    1: ─╰QSVT─┤
-
-    To see the implementation details, we can expand the circuit:
-
-    >>> q_script = qml.tape.QuantumScript(ops=[qml.qsvt_legacy(A, angles, wires=[0, 1])])
-    >>> print(q_script.expand().draw(decimals=2))
-    0: ─╭∏_ϕ(0.30)─╭BlockEncode(M0)─╭∏_ϕ(0.20)─╭BlockEncode(M0)†─╭∏_ϕ(0.10)─┤
-    1: ─╰∏_ϕ(0.30)─╰BlockEncode(M0)─╰∏_ϕ(0.20)─╰BlockEncode(M0)†─╰∏_ϕ(0.10)─┤
-    """
-
-    warnings.warn(
-        "`qml.qsvt_legacy` is deprecated and will be removed in version 0.41. "
-        "Instead, please use `qml.qsvt` functionality.",
-        qml.PennyLaneDeprecationWarning,
-    )
-
-    if qml.math.shape(A) == () or qml.math.shape(A) == (1,):
-        A = qml.math.reshape(A, [1, 1])
-
-    c, r = qml.math.shape(A)
-    global_phase, global_phase_op = None, None
-
-    with qml.QueuingManager.stop_recording():
-        UA = BlockEncode(A, wires=wires)
-    projectors = []
-
-    if convention == "Wx":
-        angles = _qsp_to_qsvt(angles)
-        global_phase = (len(angles) - 1) % 4
-
-        if global_phase:
-            with qml.QueuingManager.stop_recording():
-                global_phase_op = qml.GlobalPhase(-0.5 * np.pi * (4 - global_phase), wires=wires)
-
-    for idx, phi in enumerate(angles):
-        dim = c if idx % 2 else r
-        with qml.QueuingManager.stop_recording():
-            projectors.append(PCPhase(phi, dim=dim, wires=wires))
-
-    projectors = projectors[::-1]  # reverse order to match equation
-
-    if convention == "Wx" and global_phase:
-        return qml.prod(global_phase_op, QSVT(UA, projectors))
-    return QSVT(UA, projectors)
-
-
-# pylint: disable=too-many-branches
+# pylint: disable=too-many-branches, unused-argument
 def qsvt(A, poly, encoding_wires=None, block_encoding=None, **kwargs):
     r"""
     Implements the Quantum Singular Value Transformation (QSVT) for a matrix or Hamiltonian ``A``,
@@ -329,14 +193,6 @@ def circuit():
 
     """
 
-    if encoding_wires is None or block_encoding is None or "wires" in kwargs:
-        warnings.warn(
-            "You may be trying to use the old `qsvt` functionality (now `qml.qsvt_legacy`).\n"
-            "Make sure you pass a polynomial instead of angles.\n"
-            "Set a value for `block_encoding` to silence this warning.\n",
-            qml.PennyLaneDeprecationWarning,
-        )
-
     projectors = []
 
     # If the input A is a Hamiltonian
@@ -768,19 +624,6 @@ def compute_matrix(*args, **kwargs):
         return mat
 
 
-def _qsp_to_qsvt(angles):
-    r"""Converts qsp angles to qsvt angles."""
-    num_angles = len(angles)
-    update_vals = np.empty(num_angles)
-
-    update_vals[0] = 3 * np.pi / 4
-    update_vals[1:-1] = np.pi / 2
-    update_vals[-1] = -np.pi / 4
-    update_vals = qml.math.convert_like(update_vals, angles)
-
-    return angles + update_vals
-
-
 def _complementary_poly(poly_coeffs):
     r"""
     Computes the complementary polynomial Q given a polynomial P.
diff --git a/tests/templates/test_subroutines/test_qsvt.py b/tests/templates/test_subroutines/test_qsvt.py
index 3e16541c2b1..fe0efc8aec5 100644
--- a/tests/templates/test_subroutines/test_qsvt.py
+++ b/tests/templates/test_subroutines/test_qsvt.py
@@ -399,227 +399,6 @@ def test_copy(self):
         assert all(p1 is not p2 for p1, p2 in zip(orig_projectors, copy_projectors))
 
 
-class Testqsvt_legacy:
-    """Test the qml.qsvt_legacy function."""
-
-    def test_qsvt_legacy_deprecated(self):
-        """Test that my_feature is deprecated."""
-        with pytest.warns(qml.PennyLaneDeprecationWarning, match="`qml.qsvt_legacy` is deprecated"):
-            _ = qml.qsvt_legacy(0.3, [0.1, 0.2], [0])
-
-    @pytest.mark.parametrize(
-        ("A", "phis", "wires", "true_mat"),
-        [
-            (
-                [[0.1, 0.2], [0.3, 0.4]],
-                [0.2, 0.3],
-                [0, 1],
-                # mathematical order of gates:
-                qml.matrix(qml.PCPhase(0.2, dim=2, wires=[0, 1]))
-                @ qml.matrix(qml.BlockEncode([[0.1, 0.2], [0.3, 0.4]], wires=[0, 1]))
-                @ qml.matrix(qml.PCPhase(0.3, dim=2, wires=[0, 1])),
-            ),
-            (
-                [[0.3, 0.1], [0.2, 0.9]],
-                [0.1, 0.2, 0.3],
-                [0, 1],
-                # mathematical order of gates:
-                qml.matrix(qml.PCPhase(0.1, dim=2, wires=[0, 1]))
-                @ qml.matrix(qml.adjoint(qml.BlockEncode([[0.3, 0.1], [0.2, 0.9]], wires=[0, 1])))
-                @ qml.matrix(qml.PCPhase(0.2, dim=2, wires=[0, 1]))
-                @ qml.matrix(qml.BlockEncode([[0.3, 0.1], [0.2, 0.9]], wires=[0, 1]))
-                @ qml.matrix(qml.PCPhase(0.3, dim=2, wires=[0, 1])),
-            ),
-        ],
-    )
-    def test_output(self, A, phis, wires, true_mat):
-        """Test that qml.qsvt_legacy produces the correct output."""
-        dev = qml.device("default.qubit", wires=len(wires))
-
-        @qml.qnode(dev)
-        def circuit():
-            qml.qsvt_legacy(A, phis, wires)
-            return qml.expval(qml.PauliZ(wires=0))
-
-        observable_mat = np.kron(qml.matrix(qml.PauliZ(0)), np.eye(2))
-        true_expval = (np.conj(true_mat).T @ observable_mat @ true_mat)[0, 0]
-
-        with pytest.warns(qml.PennyLaneDeprecationWarning, match="`qml.qsvt_legacy` is deprecated"):
-            assert np.isclose(circuit(), true_expval)
-            assert np.allclose(qml.matrix(circuit)(), true_mat)
-
-    @pytest.mark.parametrize(
-        ("A", "phis", "wires", "result"),
-        [
-            (
-                [[0.1, 0.2], [0.3, 0.4]],
-                [-1.520692517929803, 0.05010380886509347],
-                [0, 1],
-                0.01,
-            ),  # angles from pyqsp give 0.1*x
-            (
-                0.3,
-                [-0.8104500678299933, 1.520692517929803, 0.7603462589648997],
-                [0],
-                0.009,
-            ),  # angles from pyqsp give 0.1*x**2
-            (
-                -1,
-                [-1.164, 0.3836, 0.383, 0.406],
-                [0],
-                -1,
-            ),  # angles from pyqsp give 0.5 * (5 * x**3 - 3 * x)
-        ],
-    )
-    def test_output_wx(self, A, phis, wires, result):
-        """Test that qml.qsvt_legacy produces the correct output."""
-        dev = qml.device("default.qubit", wires=len(wires))
-
-        @qml.qnode(dev)
-        def circuit():
-            qml.qsvt_legacy(A, phis, wires, convention="Wx")
-            return qml.expval(qml.PauliZ(wires=0))
-
-        with pytest.warns(qml.PennyLaneDeprecationWarning, match="`qml.qsvt_legacy` is deprecated"):
-            assert np.isclose(np.real(qml.matrix(circuit)())[0][0], result, rtol=1e-3)
-
-    @pytest.mark.parametrize(
-        ("A", "phis", "wires", "result"),
-        [
-            (
-                [[0.1, 0.2], [0.3, 0.4]],
-                [-1.520692517929803, 0.05010380886509347],
-                [0, 1],
-                0.01,
-            ),  # angles from pyqsp give 0.1*x
-            (
-                0.3,
-                [-0.8104500678299933, 1.520692517929803, 0.7603462589648997],
-                [0],
-                0.009,
-            ),  # angles from pyqsp give 0.1*x**2
-            (
-                -1,
-                [-1.164, 0.3836, 0.383, 0.406],
-                [0],
-                -1,
-            ),  # angles from pyqsp give 0.5 * (5 * x**3 - 3 * x)
-        ],
-    )
-    def test_matrix_wx(self, A, phis, wires, result):
-        """Assert that the matrix method produces the expected result using both call signatures."""
-
-        with pytest.warns(qml.PennyLaneDeprecationWarning, match="`qml.qsvt_legacy` is deprecated"):
-            m1 = qml.matrix(qml.qsvt_legacy(A, phis, wires, convention="Wx"))
-            m2 = qml.matrix(qml.qsvt_legacy, wire_order=wires)(A, phis, wires, convention="Wx")
-
-            assert np.isclose(np.real(m1[0, 0]), result, rtol=1e-3)
-            assert np.allclose(m1, m2)
-
-    @pytest.mark.torch
-    @pytest.mark.parametrize(
-        ("input_matrix", "angles", "wires"),
-        [([[0.1, 0.2], [0.3, 0.4]], [0.1, 0.2], [0, 1])],
-    )
-    def test_qsvt_torch(self, input_matrix, angles, wires):
-        """Test that the qsvt_legacy function matrix is correct for torch."""
-        import torch
-
-        with pytest.warns(qml.PennyLaneDeprecationWarning, match="`qml.qsvt_legacy` is deprecated"):
-            default_matrix = qml.matrix(qml.qsvt_legacy(input_matrix, angles, wires))
-
-            input_matrix = torch.tensor(input_matrix, dtype=float)
-            angles = torch.tensor(angles, dtype=float)
-
-            op = qml.qsvt_legacy(input_matrix, angles, wires)
-
-            assert np.allclose(qml.matrix(op), default_matrix)
-            assert qml.math.get_interface(qml.matrix(op)) == "torch"
-
-    @pytest.mark.jax
-    @pytest.mark.parametrize(
-        ("input_matrix", "angles", "wires"),
-        [([[0.1, 0.2], [0.3, 0.4]], [0.1, 0.2], [0, 1])],
-    )
-    def test_qsvt_jax(self, input_matrix, angles, wires):
-        """Test that the qsvt_legacy function matrix is correct for jax."""
-        import jax.numpy as jnp
-
-        with pytest.warns(qml.PennyLaneDeprecationWarning, match="`qml.qsvt_legacy` is deprecated"):
-
-            default_matrix = qml.matrix(qml.qsvt_legacy(input_matrix, angles, wires))
-
-            input_matrix = jnp.array(input_matrix)
-            angles = jnp.array(angles)
-
-            op = qml.qsvt_legacy(input_matrix, angles, wires)
-
-            assert np.allclose(qml.matrix(op), default_matrix)
-            assert qml.math.get_interface(qml.matrix(op)) == "jax"
-
-    @pytest.mark.tf
-    @pytest.mark.parametrize(
-        ("input_matrix", "angles", "wires"),
-        [([[0.1, 0.2], [0.3, 0.4]], [0.1, 0.2], [0, 1])],
-    )
-    def test_qsvt_tensorflow(self, input_matrix, angles, wires):
-        """Test that the qsvt_legacy function matrix is correct for tensorflow."""
-        import tensorflow as tf
-
-        with pytest.warns(qml.PennyLaneDeprecationWarning, match="`qml.qsvt_legacy` is deprecated"):
-
-            default_matrix = qml.matrix(qml.qsvt_legacy(input_matrix, angles, wires))
-
-            input_matrix = tf.Variable(input_matrix)
-            angles = tf.Variable(angles)
-
-            op = qml.qsvt_legacy(input_matrix, angles, wires)
-
-            assert np.allclose(qml.matrix(op), default_matrix)
-            assert qml.math.get_interface(qml.matrix(op)) == "tensorflow"
-
-    def test_qsvt_grad(self):
-        """Test that qml.grad results are the same as finite difference results"""
-
-        @qml.qnode(qml.device("default.qubit", wires=2))
-        def circuit(A, phis):
-            qml.qsvt_legacy(
-                A,
-                phis,
-                wires=[0, 1],
-            )
-            return qml.expval(qml.PauliZ(wires=0))
-
-        with pytest.warns(qml.PennyLaneDeprecationWarning, match="`qml.qsvt_legacy` is deprecated"):
-
-            A = np.array([[0.1, 0.2], [0.3, 0.4]], dtype=complex, requires_grad=True)
-            phis = np.array([0.1, 0.2, 0.3], dtype=complex, requires_grad=True)
-            y = circuit(A, phis)
-
-            mat_grad_results, phi_grad_results = qml.grad(circuit)(A, phis)
-
-            diff = 1e-8
-
-            manual_mat_results = [
-                (circuit(A + np.array([[diff, 0], [0, 0]]), phis) - y) / diff,
-                (circuit(A + np.array([[0, diff], [0, 0]]), phis) - y) / diff,
-                (circuit(A + np.array([[0, 0], [diff, 0]]), phis) - y) / diff,
-                (circuit(A + np.array([[0, 0], [0, diff]]), phis) - y) / diff,
-            ]
-
-            for idx, result in enumerate(manual_mat_results):
-                assert np.isclose(result, np.real(mat_grad_results.flatten()[idx]), atol=1e-6)
-
-            manual_phi_results = [
-                (circuit(A, phis + np.array([diff, 0, 0])) - y) / diff,
-                (circuit(A, phis + np.array([0, diff, 0])) - y) / diff,
-                (circuit(A, phis + np.array([0, 0, diff])) - y) / diff,
-            ]
-
-            for idx, result in enumerate(manual_phi_results):
-                assert np.isclose(result, np.real(phi_grad_results[idx]), atol=1e-6)
-
-
 phase_angle_data = (
     (
         [0, 0, 0],
@@ -632,43 +411,9 @@ def circuit(A, phis):
 )
 
 
-@pytest.mark.jax
-@pytest.mark.parametrize("initial_angles, expected_angles", phase_angle_data)
-def test_private_qsp_to_qsvt_jax(initial_angles, expected_angles):
-    """Test that the _qsp_to_qsvt function is jax compatible"""
-    import jax.numpy as jnp
-
-    from pennylane.templates.subroutines.qsvt import _qsp_to_qsvt
-
-    initial_angles = jnp.array(initial_angles)
-    expected_angles = jnp.array(expected_angles)
-
-    computed_angles = _qsp_to_qsvt(initial_angles)
-    jnp.allclose(computed_angles, expected_angles)
-
-
-def test_global_phase_not_alway_applied():
-    """Test that the global phase is not applied if it is 0"""
-
-    with pytest.warns(qml.PennyLaneDeprecationWarning, match="`qml.qsvt_legacy` is deprecated"):
-
-        decomposition = qml.qsvt_legacy(
-            [1], [0, 1, 2, 3, 4], wires=[0], convention="Wx"
-        ).decomposition()
-        for op in decomposition:
-            assert not isinstance(op, qml.GlobalPhase)
-
-
 class Testqsvt:
     """Test the qml.qsvt function."""
 
-    def test_qsvt_warning(self):
-        """Test that qsvt through warning."""
-        with pytest.warns(
-            qml.PennyLaneDeprecationWarning, match="You may be trying to use the old `qsvt`"
-        ):
-            qml.qsvt([[0.1, 0.2], [0.2, -0.1]], [0.1, 0, 0.1], [0, 1, 2])
-
     @pytest.mark.parametrize(
         ("A", "poly", "block_encoding", "encoding_wires"),
         [

From 14c61a6d2b252dbf19d9771d61ca5498f8ce03fb Mon Sep 17 00:00:00 2001
From: Christina Lee <christina@xanadu.ai>
Date: Wed, 15 Jan 2025 12:17:53 -0500
Subject: [PATCH 4/5] Support differentiable coefficients for observables
 (#6598)

**Context:**

With legacy operator arithmetic, we supported trainable coefficients
using the `hadamard_grad` transform. This doesn't really work anymore
with generic operator arithmetic.

**Description of the Change:**

Uses the `split_to_single_terms` transform in gradient preprocessing if
any of the observables have trainable coefficients.

In order to get this to pass tests, various bugs in
`split_non_commuting` needed fixing.

**Benefits:**

**Possible Drawbacks:**

**Related GitHub Issues:**

[sc-71490]

---------

Co-authored-by: Astral Cai <astral.cai@xanadu.ai>
Co-authored-by: Andrija Paurevic <46359773+andrijapau@users.noreply.github.com>
---
 doc/releases/changelog-dev.md                 |   4 +
 pennylane/gradients/parameter_shift.py        |  26 +++-
 pennylane/ops/functions/dot.py                |  10 +-
 pennylane/ops/op_math/sprod.py                |   2 -
 pennylane/transforms/split_non_commuting.py   | 111 +++++++++---------
 pennylane/transforms/split_to_single_terms.py |  17 +--
 pennylane/workflow/execution.py               |   1 +
 .../parameter_shift/test_parameter_shift.py   |  47 ++------
 .../test_parameter_shift_shot_vec.py          |  19 ++-
 .../interfaces/execute/test_autograd.py       |  37 ++++--
 tests/workflow/interfaces/execute/test_jax.py |  28 +++--
 .../interfaces/execute/test_tensorflow.py     |  24 ++--
 .../workflow/interfaces/execute/test_torch.py |  28 +++--
 13 files changed, 182 insertions(+), 172 deletions(-)

diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index d99bdb3383d..25c75208671 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -6,6 +6,9 @@
 
 <h3>Improvements 🛠</h3>
 
+* The coefficients of observables now have improved differentiability.
+  [(#6598)](https://github.com/PennyLaneAI/pennylane/pull/6598)
+
 <h3>Breaking changes 💔</h3>
 
 * Removed method `qsvt_legacy` along with its private helper `_qsp_to_qsvt`
@@ -26,3 +29,4 @@ This release contains contributions from (in alphabetical order):
 
 Yushao Chen,
 Diksha Dhawan,
+Christina Lee,
diff --git a/pennylane/gradients/parameter_shift.py b/pennylane/gradients/parameter_shift.py
index 5c25237db58..6fcdd17df19 100644
--- a/pennylane/gradients/parameter_shift.py
+++ b/pennylane/gradients/parameter_shift.py
@@ -752,6 +752,12 @@ def _param_shift_stopping_condition(op) -> bool:
     return True
 
 
+def _inplace_set_trainable_params(tape):
+    """Update all the trainable params in place."""
+    params = tape.get_parameters(trainable_only=False)
+    tape.trainable_params = qml.math.get_trainable_indices(params)
+
+
 def _expand_transform_param_shift(
     tape: QuantumScript,
     argnum=None,
@@ -769,11 +775,21 @@ def _expand_transform_param_shift(
         name="param_shift",
         error=qml.operation.DecompositionUndefinedError,
     )
-    if new_tape is tape:
-        return [tape], postprocessing
-    params = new_tape.get_parameters(trainable_only=False)
-    new_tape.trainable_params = qml.math.get_trainable_indices(params)
-    return [new_tape], postprocessing
+    if any(
+        qml.math.requires_grad(d) for mp in tape.measurements for d in getattr(mp.obs, "data", [])
+    ):
+        try:
+            batch, postprocessing = qml.transforms.split_to_single_terms(new_tape)
+        except RuntimeError as e:
+            raise ValueError(
+                "Can only differentiate Hamiltonian "
+                f"coefficients for expectations, not {tape.measurements}."
+            ) from e
+    else:
+        batch = [new_tape]
+    if len(batch) > 1 or batch[0] is not tape:
+        _ = [_inplace_set_trainable_params(t) for t in batch]
+    return batch, postprocessing
 
 
 @partial(
diff --git a/pennylane/ops/functions/dot.py b/pennylane/ops/functions/dot.py
index de4c2f76225..e2e2d11a470 100644
--- a/pennylane/ops/functions/dot.py
+++ b/pennylane/ops/functions/dot.py
@@ -142,10 +142,12 @@ def dot(
                 f"ops must be an Iterable of {t.__name__}'s, not a {t.__name__} itself."
             )
 
-    if len(coeffs) != len(ops):
-        raise ValueError("Number of coefficients and operators does not match.")
-    if len(coeffs) == 0 and len(ops) == 0:
-        raise ValueError("Cannot compute the dot product of an empty sequence.")
+    # tensorflow variables have no len
+    if qml.math.get_interface(coeffs) != "tensorflow":
+        if len(coeffs) != len(ops):
+            raise ValueError("Number of coefficients and operators does not match.")
+        if len(coeffs) == 0 and len(ops) == 0:
+            raise ValueError("Cannot compute the dot product of an empty sequence.")
 
     for t in (Operator, PauliWord, PauliSentence):
         if isinstance(ops, t):
diff --git a/pennylane/ops/op_math/sprod.py b/pennylane/ops/op_math/sprod.py
index bd8e6a3a3e5..eca3bee5cf9 100644
--- a/pennylane/ops/op_math/sprod.py
+++ b/pennylane/ops/op_math/sprod.py
@@ -15,7 +15,6 @@
 This file contains the implementation of the SProd class which contains logic for
 computing the scalar product of operations.
 """
-from copy import copy
 from typing import Union
 
 import pennylane as qml
@@ -148,7 +147,6 @@ def __init__(
         elif (base_pauli_rep := getattr(self.base, "pauli_rep", None)) and (
             self.batch_size is None
         ):
-            scalar = copy(self.scalar)
 
             pr = {pw: qnp.dot(coeff, scalar) for pw, coeff in base_pauli_rep.items()}
             self._pauli_rep = qml.pauli.PauliSentence(pr)
diff --git a/pennylane/transforms/split_non_commuting.py b/pennylane/transforms/split_non_commuting.py
index 1faab45de73..ba282962233 100644
--- a/pennylane/transforms/split_non_commuting.py
+++ b/pennylane/transforms/split_non_commuting.py
@@ -18,11 +18,11 @@
 
 # pylint: disable=too-many-arguments,too-many-boolean-expressions
 
-from functools import partial
+from functools import partial, wraps
 from typing import Optional
 
 import pennylane as qml
-from pennylane.measurements import ExpectationMP, MeasurementProcess, Shots, StateMP
+from pennylane.measurements import ExpectationMP, MeasurementProcess, StateMP
 from pennylane.ops import Prod, SProd, Sum
 from pennylane.tape import QuantumScript, QuantumScriptBatch
 from pennylane.transforms import transform
@@ -36,6 +36,16 @@ def null_postprocessing(results):
     return results[0]
 
 
+def shot_vector_support(initial_postprocessing: PostprocessingFn) -> PostprocessingFn:
+    """Convert a postprocessing function to one with shot vector support."""
+
+    @wraps(initial_postprocessing)
+    def shot_vector_postprocessing(results):
+        return tuple(initial_postprocessing(r) for r in zip(*results))
+
+    return shot_vector_postprocessing
+
+
 @transform
 def split_non_commuting(
     tape: QuantumScript, grouping_strategy: Optional[str] = "default"
@@ -280,13 +290,15 @@ def circuit(x):
     if grouping_strategy is None:
         measurements = list(single_term_obs_mps.keys())
         tapes = [tape.copy(measurements=[m]) for m in measurements]
-        return tapes, partial(
+        fn = partial(
             _processing_fn_no_grouping,
             single_term_obs_mps=single_term_obs_mps,
             offsets=offsets,
-            shots=tape.shots,
             batch_size=tape.batch_size,
         )
+        if tape.shots.has_partitioned_shots:
+            fn = shot_vector_support(fn)
+        return tapes, fn
 
     if grouping_strategy == "wires" or any(
         m.obs is not None and not qml.pauli.is_pauli_word(m.obs) for m in single_term_obs_mps
@@ -360,14 +372,16 @@ def _split_ham_with_grouping(tape: qml.tape.QuantumScript):
             group_sizes.append(group_size)
 
     tapes = [tape.copy(measurements=mps) for mps in mp_groups]
-    return tapes, partial(
+    fn = partial(
         _processing_fn_with_grouping,
         single_term_obs_mps=single_term_obs_mps,
         offsets=[offset],
         group_sizes=group_sizes,
-        shots=tape.shots,
         batch_size=tape.batch_size,
     )
+    if tape.shots.has_partitioned_shots:
+        fn = shot_vector_support(fn)
+    return tapes, fn
 
 
 def _split_using_qwc_grouping(
@@ -424,16 +438,17 @@ def _split_using_qwc_grouping(
             0,
         )
         group_sizes.append(1)
-
     tapes = [tape.copy(measurements=mps) for mps in mp_groups]
-    return tapes, partial(
+    fn = partial(
         _processing_fn_with_grouping,
         single_term_obs_mps=single_term_obs_mps_grouped,
         offsets=offsets,
         group_sizes=group_sizes,
-        shots=tape.shots,
         batch_size=tape.batch_size,
     )
+    if tape.shots.has_partitioned_shots:
+        fn = shot_vector_support(fn)
+    return tapes, fn
 
 
 def _split_using_wires_grouping(
@@ -497,14 +512,16 @@ def _split_using_wires_grouping(
             num_groups += 1
 
     tapes = [tape.copy(measurements=mps) for mps in mp_groups]
-    return tapes, partial(
+    fn = partial(
         _processing_fn_with_grouping,
         single_term_obs_mps=single_term_obs_mps_grouped,
         offsets=offsets,
         group_sizes=group_sizes,
-        shots=tape.shots,
         batch_size=tape.batch_size,
     )
+    if tape.shots.has_partitioned_shots:
+        fn = shot_vector_support(fn)
+    return tapes, fn
 
 
 def _split_all_multi_term_obs_mps(tape: qml.tape.QuantumScript):
@@ -572,8 +589,7 @@ def _processing_fn_no_grouping(
     res: ResultBatch,
     single_term_obs_mps: dict[MeasurementProcess, tuple[list[int], list[Union[float, TensorLike]]]],
     offsets: list[Union[float, TensorLike]],
-    shots: Shots,
-    batch_size: int,
+    batch_size: Union[None, int],
 ):
     """Postprocessing function for the split_non_commuting transform without grouping.
 
@@ -592,12 +608,22 @@ def _processing_fn_no_grouping(
     coeffs_for_each_mp = [[] for _ in offsets]
 
     for smp_idx, (_, (mp_indices, coeffs)) in enumerate(single_term_obs_mps.items()):
-
         for mp_idx, coeff in zip(mp_indices, coeffs):
             res_batch_for_each_mp[mp_idx].append(res[smp_idx])
             coeffs_for_each_mp[mp_idx].append(coeff)
 
-    return _res_for_each_mp(res_batch_for_each_mp, coeffs_for_each_mp, offsets, shots, batch_size)
+    result_shape = (batch_size,) if batch_size and batch_size > 1 else ()
+    # Sum up the results for each original measurement
+
+    res_for_each_mp = [
+        _sum_terms(_sub_res, coeffs, offset, result_shape)
+        for _sub_res, coeffs, offset in zip(res_batch_for_each_mp, coeffs_for_each_mp, offsets)
+    ]
+    # res_for_each_mp should have shape (n_mps, [,n_shots] [,batch_size])
+    if len(res_for_each_mp) == 1:
+        return res_for_each_mp[0]
+
+    return tuple(res_for_each_mp)
 
 
 def _processing_fn_with_grouping(
@@ -605,9 +631,8 @@ def _processing_fn_with_grouping(
     single_term_obs_mps: dict[
         MeasurementProcess, tuple[list[int], list[Union[float, TensorLike]], int, int]
     ],
-    offsets: list[Union[float, TensorLike]],
+    offsets: list[TensorLike],
     group_sizes: list[int],
-    shots: Shots,
     batch_size: int,
 ):
     """Postprocessing function for the split_non_commuting transform with grouping.
@@ -636,26 +661,16 @@ def _processing_fn_with_grouping(
         res_group = res[group_idx]  # ([n_shots] [,n_mps] [,batch_size])
         group_size = group_sizes[group_idx]
 
-        if group_size > 1 and shots.has_partitioned_shots:
-            # Each result should have shape ([n_shots] [,batch_size])
-            sub_res = [_res[mp_idx_in_group] for _res in res_group]
-        else:
-            # If there is only one term in the group, the n_mps dimension would have
-            # been squeezed out, use the entire result directly.
-            sub_res = res_group if group_size == 1 else res_group[mp_idx_in_group]
+        # If there is only one term in the group, the n_mps dimension would have
+        # been squeezed out, use the entire result directly.
+        sub_res = res_group if group_size == 1 else res_group[mp_idx_in_group]
 
         # Add this result to the result batch for the corresponding original measurement
         for mp_idx, coeff in zip(mp_indices, coeffs):
             res_batch_for_each_mp[mp_idx].append(sub_res)
             coeffs_for_each_mp[mp_idx].append(coeff)
 
-    return _res_for_each_mp(res_batch_for_each_mp, coeffs_for_each_mp, offsets, shots, batch_size)
-
-
-def _res_for_each_mp(res_batch_for_each_mp, coeffs_for_each_mp, offsets, shots, batch_size):
-    """Helper function that combines a result batch into results for each mp"""
-
-    result_shape = _infer_result_shape(shots, batch_size)
+    result_shape = (batch_size,) if batch_size and batch_size > 1 else ()
 
     # Sum up the results for each original measurement
     res_for_each_mp = [
@@ -667,14 +682,6 @@ def _res_for_each_mp(res_batch_for_each_mp, coeffs_for_each_mp, offsets, shots,
     if len(res_for_each_mp) == 1:
         return res_for_each_mp[0]
 
-    if shots.has_partitioned_shots:
-        # If the shot vector dimension exists, it should be moved to the first axis
-        # Basically, the shape becomes (n_shots, n_mps, [,batch_size])
-        res_for_each_mp = [
-            tuple(res_for_each_mp[j][i] for j in range(len(res_for_each_mp)))
-            for i in range(shots.num_copies)
-        ]
-
     return tuple(res_for_each_mp)
 
 
@@ -685,9 +692,13 @@ def _sum_terms(
     shape: tuple,
 ) -> Result:
     """Sum results from measurements of multiple terms in a multi-term observable."""
-
-    # Trivially return the original result
-    if coeffs == [1] and offset == 0:
+    if (
+        coeffs
+        and not qml.math.is_abstract(coeffs[0])
+        and not qml.math.is_abstract(offset)
+        and coeffs == [1]
+        and offset == 0
+    ):
         return res[0]
 
     # The shape of res at this point is (n_terms, [,n_shots] [,batch_size])
@@ -695,10 +706,11 @@ def _sum_terms(
     for c, r in zip(coeffs, res):
         if qml.math.get_interface(r) == "autograd":
             r = qml.math.array(r)
-        dot_products.append(qml.math.dot(qml.math.squeeze(r), c))
+        if isinstance(r, (list, tuple)):
+            r = qml.math.stack(r)
+        dot_products.append(qml.math.dot(c, qml.math.squeeze(r)))
     if len(dot_products) == 0:
         return qml.math.ones(shape) * offset
-
     summed_dot_products = qml.math.sum(qml.math.stack(dot_products), axis=0)
     if qml.math.get_interface(offset) == "autograd" and qml.math.requires_grad(summed_dot_products):
         offset = qml.math.array(offset)
@@ -718,14 +730,3 @@ def _mp_to_obs(mp: MeasurementProcess, tape: qml.tape.QuantumScript) -> qml.oper
 
     obs_wires = mp.wires if mp.wires else tape.wires
     return qml.prod(*(qml.Z(wire) for wire in obs_wires))
-
-
-def _infer_result_shape(shots: Shots, batch_size: int) -> tuple:
-    """Based on the result, infer the ([,n_shots] [,batch_size]) shape of the result."""
-
-    shape = ()
-    if shots.has_partitioned_shots:
-        shape += (shots.num_copies,)
-    if batch_size and batch_size > 1:
-        shape += (batch_size,)
-    return shape
diff --git a/pennylane/transforms/split_to_single_terms.py b/pennylane/transforms/split_to_single_terms.py
index 3bbe2cbb0f9..f794d8424ec 100644
--- a/pennylane/transforms/split_to_single_terms.py
+++ b/pennylane/transforms/split_to_single_terms.py
@@ -24,6 +24,7 @@
 from pennylane.transforms.split_non_commuting import (
     _processing_fn_no_grouping,
     _split_all_multi_term_obs_mps,
+    shot_vector_support,
 )
 
 
@@ -162,23 +163,13 @@ def post_processing_split_sums(res):
             _processing_fn_no_grouping,
             single_term_obs_mps=single_term_obs_mps,
             offsets=offsets,
-            shots=tape.shots,
             batch_size=tape.batch_size,
         )
 
-        if len(new_tape.measurements) == 1:
-            return process(res)
-
         # we go from ((mp1_res, mp2_res, mp3_res),) as result output
         # to (mp1_res, mp2_res, mp3_res) as expected by _processing_fn_no_grouping
-        res = res[0]
-        if tape.shots.has_partitioned_shots:
-            # swap dimension order of mps vs shot copies for _processing_fn_no_grouping
-            res = [
-                tuple(res[j][i] for j in range(tape.shots.num_copies))
-                for i in range(len(new_tape.measurements))
-            ]
-
-        return process(res)
+        return process(res if len(new_tape.measurements) == 1 else res[0])
 
+    if tape.shots.has_partitioned_shots:
+        return (new_tape,), shot_vector_support(post_processing_split_sums)
     return (new_tape,), post_processing_split_sums
diff --git a/pennylane/workflow/execution.py b/pennylane/workflow/execution.py
index 7c9fb93b001..74fe6949c96 100644
--- a/pennylane/workflow/execution.py
+++ b/pennylane/workflow/execution.py
@@ -44,6 +44,7 @@ def execute(
     device: Union["qml.devices.LegacyDevice", "qml.devices.Device"],
     diff_method: Optional[Union[Callable, str, qml.transforms.core.TransformDispatcher]] = None,
     interface: Optional[Union[str, Interface]] = Interface.AUTO,
+    *,
     transform_program=None,
     inner_transform=None,
     config=None,
diff --git a/tests/gradients/parameter_shift/test_parameter_shift.py b/tests/gradients/parameter_shift/test_parameter_shift.py
index 499fbdf2ad5..f35a4b4fc95 100644
--- a/tests/gradients/parameter_shift/test_parameter_shift.py
+++ b/tests/gradients/parameter_shift/test_parameter_shift.py
@@ -3363,14 +3363,14 @@ def test_not_expval_error(self, broadcast):
         tape = qml.tape.QuantumScript.from_queue(q)
         tape.trainable_params = {2, 3, 4}
 
-        with pytest.raises(ValueError, match="for expectations, not var"):
+        with pytest.raises(ValueError, match="for expectations, not"):
             qml.gradients.param_shift(tape, broadcast=broadcast)
 
     def test_not_expval_pass_if_not_trainable_hamiltonian(self, broadcast):
         """Test that if the variance of a non-trainable Hamiltonian is requested,
         no error is raised"""
         obs = [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliX(1), qml.PauliY(0)]
-        coeffs = np.array([0.1, 0.2, 0.3])
+        coeffs = np.array([0.1, 0.2, 0.3], requires_grad=False)
         H = qml.Hamiltonian(coeffs, obs)
 
         weights = np.array([0.4, 0.5])
@@ -3393,7 +3393,7 @@ def test_no_trainable_coeffs(self, mocker, tol, broadcast):
         spy = mocker.spy(qml.gradients, "hamiltonian_grad")
 
         obs = [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliX(1), qml.PauliY(0)]
-        coeffs = np.array([0.1, 0.2, 0.3])
+        coeffs = np.array([0.1, 0.2, 0.3], requires_grad=False)
         H = qml.Hamiltonian(coeffs, obs)
 
         weights = np.array([0.4, 0.5])
@@ -3433,10 +3433,9 @@ def test_no_trainable_coeffs(self, mocker, tol, broadcast):
         assert np.allclose(res[0], expected[0], atol=tol, rtol=0)
         assert np.allclose(res[1], expected[1], atol=tol, rtol=0)
 
-    def test_trainable_coeffs(self, mocker, tol, broadcast):
+    def test_trainable_coeffs(self, tol, broadcast):
         """Test trainable Hamiltonian coefficients"""
         dev = qml.device("default.qubit", wires=2)
-        spy = mocker.spy(qml.gradients, "hamiltonian_grad")
 
         obs = [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliX(1), qml.PauliY(0)]
         coeffs = np.array([0.1, 0.2, 0.3])
@@ -3462,33 +3461,21 @@ def test_trainable_coeffs(self, mocker, tol, broadcast):
         tapes, fn = qml.gradients.param_shift(tape, broadcast=broadcast)
         # two (broadcasted if broadcast=True) shifts per rotation gate
         # one circuit per trainable H term
-        assert len(tapes) == (2 + 2 if broadcast else 2 * 2 + 2)
-        assert [t.batch_size for t in tapes] == ([2, 2, None, None] if broadcast else [None] * 6)
-        spy.assert_called()
+        assert len(tapes) == (2 if broadcast else 2 * 2)
+        assert [t.batch_size for t in tapes] == ([2, 2] if broadcast else [None] * 4)
 
         res = fn(dev.execute(tapes))
-        assert isinstance(res, tuple)
-        assert len(res) == 4
-        assert res[0].shape == ()
-        assert res[1].shape == ()
-        assert res[2].shape == ()
-        assert res[3].shape == ()
 
         expected = [
             -c * np.cos(x) * np.sin(y) - np.sin(x) * (a + b * np.sin(y)),
             b * np.cos(x) * np.cos(y) - c * np.cos(y) * np.sin(x),
-            np.cos(x),
-            -(np.sin(x) * np.sin(y)),
         ]
         assert np.allclose(res[0], expected[0], atol=tol, rtol=0)
         assert np.allclose(res[1], expected[1], atol=tol, rtol=0)
-        assert np.allclose(res[2], expected[2], atol=tol, rtol=0)
-        assert np.allclose(res[3], expected[3], atol=tol, rtol=0)
 
-    def test_multiple_hamiltonians(self, mocker, tol, broadcast):
+    def test_multiple_hamiltonians(self, tol, broadcast):
         """Test multiple trainable Hamiltonian coefficients"""
         dev = qml.device("default.qubit", wires=2)
-        spy = mocker.spy(qml.gradients, "hamiltonian_grad")
 
         obs = [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliX(1), qml.PauliY(0)]
         coeffs = np.array([0.1, 0.2, 0.3])
@@ -3520,24 +3507,20 @@ def test_multiple_hamiltonians(self, mocker, tol, broadcast):
         tapes, fn = qml.gradients.param_shift(tape, broadcast=broadcast)
         # two shifts per rotation gate (in one batched tape if broadcasting),
         # one circuit per trainable H term
-        assert len(tapes) == 2 * (1 if broadcast else 2) + 3
-        spy.assert_called()
+        assert len(tapes) == 2 * (1 if broadcast else 2)
 
         res = fn(dev.execute(tapes))
         assert isinstance(res, tuple)
         assert len(res) == 2
-        assert len(res[0]) == 5
-        assert len(res[1]) == 5
+        assert len(res[0]) == 2
+        assert len(res[1]) == 2
 
         expected = [
             [
                 -c * np.cos(x) * np.sin(y) - np.sin(x) * (a + b * np.sin(y)),
                 b * np.cos(x) * np.cos(y) - c * np.cos(y) * np.sin(x),
-                np.cos(x),
-                -(np.sin(x) * np.sin(y)),
-                0,
             ],
-            [-d * np.sin(x), 0, 0, 0, np.cos(x)],
+            [-d * np.sin(x), 0],
         ]
 
         assert np.allclose(np.stack(res), expected, atol=tol, rtol=0)
@@ -3574,12 +3557,8 @@ def cost_fn_expected(weights, coeffs1, coeffs2):
             [
                 -c * np.cos(x) * np.sin(y) - np.sin(x) * (a + b * np.sin(y)),
                 b * np.cos(x) * np.cos(y) - c * np.cos(y) * np.sin(x),
-                np.cos(x),
-                np.cos(x) * np.sin(y),
-                -(np.sin(x) * np.sin(y)),
-                0,
             ],
-            [-d * np.sin(x), 0, 0, 0, 0, np.cos(x)],
+            [-d * np.sin(x), 0],
         ]
 
     @pytest.mark.autograd
@@ -3670,7 +3649,7 @@ def test_jax(self, tol, broadcast):
 
         res = self.cost_fn(weights, coeffs1, coeffs2, dev, broadcast)
         expected = self.cost_fn_expected(weights, coeffs1, coeffs2)
-        assert np.allclose(res, np.array(expected), atol=tol, rtol=0)
+        assert np.allclose(np.array(res)[:, :2], np.array(expected), atol=tol, rtol=0)
 
         # TODO: test when Hessians are supported with the new return types
         # second derivative wrt to Hamiltonian coefficients should be zero
diff --git a/tests/gradients/parameter_shift/test_parameter_shift_shot_vec.py b/tests/gradients/parameter_shift/test_parameter_shift_shot_vec.py
index 1967a44b20a..34478ccecd0 100644
--- a/tests/gradients/parameter_shift/test_parameter_shift_shot_vec.py
+++ b/tests/gradients/parameter_shift/test_parameter_shift_shot_vec.py
@@ -2032,7 +2032,7 @@ def test_not_expval_error(self, broadcast):
         tape = qml.tape.QuantumScript.from_queue(q, shots=shot_vec)
         tape.trainable_params = {2, 3, 4}
 
-        with pytest.raises(ValueError, match="for expectations, not var"):
+        with pytest.raises(ValueError, match="for expectations, not"):
             qml.gradients.param_shift(tape, broadcast=broadcast)
 
     def test_no_trainable_coeffs(self, mocker, broadcast, tol):
@@ -2078,7 +2078,6 @@ def test_no_trainable_coeffs(self, mocker, broadcast, tol):
             b * np.cos(x) * np.cos(y) - c * np.cos(y) * np.sin(x),
         ]
         for res in all_res:
-            assert isinstance(res, tuple)
             assert len(res) == 2
             assert res[0].shape == ()
             assert res[1].shape == ()
@@ -2086,17 +2085,16 @@ def test_no_trainable_coeffs(self, mocker, broadcast, tol):
             assert np.allclose(res[0], expected[0], atol=tol, rtol=0)
             assert np.allclose(res[1], expected[1], atol=tol, rtol=0)
 
-    def test_trainable_coeffs(self, mocker, broadcast, tol):
+    def test_trainable_coeffs(self, broadcast, tol):
         """Test trainable Hamiltonian coefficients"""
         shot_vec = many_shots_shot_vector
         dev = qml.device("default.qubit", wires=2, shots=shot_vec)
-        spy = mocker.spy(qml.gradients, "hamiltonian_grad")
 
         obs = [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliX(1), qml.PauliY(0)]
-        coeffs = np.array([0.1, 0.2, 0.3])
+        coeffs = qml.numpy.array([0.1, 0.2, 0.3])
         H = qml.Hamiltonian(coeffs, obs)
 
-        weights = np.array([0.4, 0.5])
+        weights = qml.numpy.array([0.4, 0.5])
 
         with qml.queuing.AnnotatedQueue() as q:
             qml.RX(weights[0], wires=0)
@@ -2116,19 +2114,16 @@ def test_trainable_coeffs(self, mocker, broadcast, tol):
         tapes, fn = qml.gradients.param_shift(tape, broadcast=broadcast)
         # two (broadcasted if broadcast=True) shifts per rotation gate
         # one circuit per trainable H term
-        assert len(tapes) == (2 + 2 if broadcast else 2 * 2 + 2)
-        assert [t.batch_size for t in tapes] == ([2, 2, None, None] if broadcast else [None] * 6)
-        spy.assert_called()
+        assert len(tapes) == (2 if broadcast else 2 * 2)
+        assert [t.batch_size for t in tapes] == ([2, 2] if broadcast else [None] * 4)
 
         res = fn(dev.execute(tapes))
         assert isinstance(res, tuple)
-        assert qml.math.shape(res) == (3, 4)
+        assert qml.math.shape(res) == (3, 2)
 
         expected = [
             -c * np.cos(x) * np.sin(y) - np.sin(x) * (a + b * np.sin(y)),
             b * np.cos(x) * np.cos(y) - c * np.cos(y) * np.sin(x),
-            np.cos(x),
-            -(np.sin(x) * np.sin(y)),
         ]
         for r in res:
             assert qml.math.allclose(r, expected, atol=shot_vec_tol)
diff --git a/tests/workflow/interfaces/execute/test_autograd.py b/tests/workflow/interfaces/execute/test_autograd.py
index 83a3e161654..7bab72e0c81 100644
--- a/tests/workflow/interfaces/execute/test_autograd.py
+++ b/tests/workflow/interfaces/execute/test_autograd.py
@@ -774,24 +774,32 @@ def cost_fn(x):
 
 
 @pytest.mark.parametrize("execute_kwargs, shots, device_name", test_matrix)
+@pytest.mark.parametrize("constructor", (qml.Hamiltonian, qml.dot, "dunders"))
 class TestHamiltonianWorkflows:
     """Test that tapes ending with expectations
     of Hamiltonians provide correct results and gradients"""
 
+    # pylint: disable=too-many-arguments, too-many-positional-arguments
     @pytest.fixture
-    def cost_fn(self, execute_kwargs, shots, device_name, seed):
+    def cost_fn(self, execute_kwargs, shots, device_name, seed, constructor):
         """Cost function for gradient tests"""
 
         device = get_device(device_name, seed=seed)
 
         def _cost_fn(weights, coeffs1, coeffs2):
-            obs1 = [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliX(1), qml.PauliY(0)]
-            H1 = qml.Hamiltonian(coeffs1, obs1)
-            H1 = qml.pauli.pauli_sentence(H1).operation()
 
-            obs2 = [qml.PauliZ(0)]
-            H2 = qml.Hamiltonian(coeffs2, obs2)
-            H2 = qml.pauli.pauli_sentence(H2).operation()
+            if constructor == "dunders":
+                H1 = (
+                    coeffs1[0] * qml.Z(0) + coeffs1[1] * qml.Z(0) @ qml.X(1) + coeffs1[2] * qml.Y(0)
+                )
+                H2 = coeffs2[0] * qml.Z(0)
+            else:
+
+                obs1 = [qml.Z(0), qml.Z(0) @ qml.X(1), qml.Y(0)]
+                H1 = constructor(coeffs1, obs1)
+
+                obs2 = [qml.PauliZ(0)]
+                H2 = constructor(coeffs2, obs2)
 
             with qml.queuing.AnnotatedQueue() as q:
                 qml.RX(weights[0], wires=0)
@@ -856,12 +864,16 @@ def test_multiple_hamiltonians_not_trainable(self, cost_fn, shots):
         else:
             assert np.allclose(res, expected, atol=atol_for_shots(shots), rtol=0)
 
-    def test_multiple_hamiltonians_trainable(self, execute_kwargs, cost_fn, shots):
+    def test_multiple_hamiltonians_trainable(self, execute_kwargs, cost_fn, shots, constructor):
         """Test hamiltonian with trainable parameters."""
         if execute_kwargs["diff_method"] == "adjoint":
             pytest.skip("trainable hamiltonians not supported with adjoint")
-        if execute_kwargs["diff_method"] != "backprop":
-            pytest.xfail(reason="parameter shift derivatives do not yet support sums.")
+        if constructor == "dunders":
+            pytest.xfail("autograd does not like constructing an sprod via dunder.")
+        if shots.has_partitioned_shots:
+            pytest.xfail(
+                "multiple hamiltonians with shot vectors does not seem to be differentiable."
+            )
 
         coeffs1 = pnp.array([0.1, 0.2, 0.3], requires_grad=True)
         coeffs2 = pnp.array([0.7], requires_grad=True)
@@ -878,8 +890,7 @@ def test_multiple_hamiltonians_trainable(self, execute_kwargs, cost_fn, shots):
         res = pnp.hstack(qml.jacobian(cost_fn)(weights, coeffs1, coeffs2))
         expected = self.cost_fn_jacobian(weights, coeffs1, coeffs2)
         if shots.has_partitioned_shots:
-            pytest.xfail(
-                "multiple hamiltonians with shot vectors does not seem to be differentiable."
-            )
+            assert qml.math.allclose(res[:2, :], expected, atol=atol_for_shots(shots), rtol=0)
+            assert qml.math.allclose(res[2:, :], expected, atol=atol_for_shots(shots), rtol=0)
         else:
             assert qml.math.allclose(res, expected, atol=atol_for_shots(shots), rtol=0)
diff --git a/tests/workflow/interfaces/execute/test_jax.py b/tests/workflow/interfaces/execute/test_jax.py
index e9d07ccabf7..7a7c9468ca5 100644
--- a/tests/workflow/interfaces/execute/test_jax.py
+++ b/tests/workflow/interfaces/execute/test_jax.py
@@ -718,24 +718,31 @@ def cost_fn(x):
 
 
 @pytest.mark.parametrize("execute_kwargs, shots, device_name", test_matrix)
+@pytest.mark.parametrize("constructor", (qml.Hamiltonian, qml.dot, "dunders"))
 class TestHamiltonianWorkflows:
     """Test that tapes ending with expectations
     of Hamiltonians provide correct results and gradients"""
 
+    # pylint: disable=too-many-arguments, too-many-positional-arguments
     @pytest.fixture
-    def cost_fn(self, execute_kwargs, shots, device_name, seed):
+    def cost_fn(self, execute_kwargs, shots, device_name, seed, constructor):
         """Cost function for gradient tests"""
 
         device = get_device(device_name, seed)
 
         def _cost_fn(weights, coeffs1, coeffs2):
-            obs1 = [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliX(1), qml.PauliY(0)]
-            H1 = qml.Hamiltonian(coeffs1, obs1)
-            H1 = qml.pauli.pauli_sentence(H1).operation()
+            if constructor == "dunders":
+                H1 = (
+                    coeffs1[0] * qml.Z(0) + coeffs1[1] * qml.Z(0) @ qml.X(1) + coeffs1[2] * qml.Y(0)
+                )
+                H2 = coeffs2[0] * qml.Z(0)
+            else:
 
-            obs2 = [qml.PauliZ(0)]
-            H2 = qml.Hamiltonian(coeffs2, obs2)
-            H2 = qml.pauli.pauli_sentence(H2).operation()
+                obs1 = [qml.Z(0), qml.Z(0) @ qml.X(1), qml.Y(0)]
+                H1 = constructor(coeffs1, obs1)
+
+                obs2 = [qml.PauliZ(0)]
+                H2 = constructor(coeffs2, obs2)
 
             with qml.queuing.AnnotatedQueue() as q:
                 qml.RX(weights[0], wires=0)
@@ -807,8 +814,6 @@ def test_multiple_hamiltonians_trainable(self, execute_kwargs, cost_fn, shots):
         """Test hamiltonian with trainable parameters."""
         if execute_kwargs["diff_method"] == "adjoint":
             pytest.xfail("trainable hamiltonians not supported with adjoint")
-        if execute_kwargs["diff_method"] != "backprop":
-            pytest.xfail(reason="parameter shift derivatives do not yet support sums.")
 
         coeffs1 = jnp.array([0.1, 0.2, 0.3])
         coeffs2 = jnp.array([0.7])
@@ -825,8 +830,7 @@ def test_multiple_hamiltonians_trainable(self, execute_kwargs, cost_fn, shots):
         res = jnp.hstack(jax.jacobian(cost_fn, argnums=[0, 1, 2])(weights, coeffs1, coeffs2))
         expected = self.cost_fn_jacobian(weights, coeffs1, coeffs2)
         if shots.has_partitioned_shots:
-            pytest.xfail(
-                "multiple hamiltonians with shot vectors does not seem to be differentiable."
-            )
+            assert np.allclose(res[:2, :], expected, atol=atol_for_shots(shots), rtol=0)
+            assert np.allclose(res[2:, :], expected, atol=atol_for_shots(shots), rtol=0)
         else:
             assert np.allclose(res, expected, atol=atol_for_shots(shots), rtol=0)
diff --git a/tests/workflow/interfaces/execute/test_tensorflow.py b/tests/workflow/interfaces/execute/test_tensorflow.py
index 5fe780f0f5a..a9f93f42c77 100644
--- a/tests/workflow/interfaces/execute/test_tensorflow.py
+++ b/tests/workflow/interfaces/execute/test_tensorflow.py
@@ -729,24 +729,31 @@ def cost_fn(x):
 
 
 @pytest.mark.parametrize("execute_kwargs, shots, device_name", test_matrix)
+@pytest.mark.parametrize("constructor", (qml.Hamiltonian, qml.dot, "dunders"))
 class TestHamiltonianWorkflows:
     """Test that tapes ending with expectations
     of Hamiltonians provide correct results and gradients"""
 
+    # pylint: disable=too-many-arguments, too-many-positional-arguments
     @pytest.fixture
-    def cost_fn(self, execute_kwargs, shots, device_name, seed):
+    def cost_fn(self, execute_kwargs, shots, device_name, seed, constructor):
         """Cost function for gradient tests"""
 
         device = qml.device(device_name, seed=seed)
 
         def _cost_fn(weights, coeffs1, coeffs2):
-            obs1 = [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliX(1), qml.PauliY(0)]
-            H1 = qml.Hamiltonian(coeffs1, obs1)
-            H1 = qml.pauli.pauli_sentence(H1).operation()
+            if constructor == "dunders":
+                H1 = (
+                    coeffs1[0] * qml.Z(0) + coeffs1[1] * qml.Z(0) @ qml.X(1) + coeffs1[2] * qml.Y(0)
+                )
+                H2 = coeffs2[0] * qml.Z(0)
+            else:
 
-            obs2 = [qml.PauliZ(0)]
-            H2 = qml.Hamiltonian(coeffs2, obs2)
-            H2 = qml.pauli.pauli_sentence(H2).operation()
+                obs1 = [qml.Z(0), qml.Z(0) @ qml.X(1), qml.Y(0)]
+                H1 = constructor(coeffs1, obs1)
+
+                obs2 = [qml.PauliZ(0)]
+                H2 = constructor(coeffs2, obs2)
 
             with qml.queuing.AnnotatedQueue() as q:
                 qml.RX(weights[0], wires=0)
@@ -811,9 +818,6 @@ def test_multiple_hamiltonians_trainable(self, cost_fn, execute_kwargs, shots):
         """Test hamiltonian with trainable parameters."""
         if execute_kwargs["diff_method"] == "adjoint":
             pytest.skip("trainable hamiltonians not supported with adjoint")
-        if execute_kwargs["diff_method"] != "backprop":
-            pytest.xfail(reason="parameter shift derivatives do not yet support sums.")
-
         coeffs1 = tf.Variable([0.1, 0.2, 0.3], dtype=tf.float64)
         coeffs2 = tf.Variable([0.7], dtype=tf.float64)
         weights = tf.Variable([0.4, 0.5], dtype=tf.float64)
diff --git a/tests/workflow/interfaces/execute/test_torch.py b/tests/workflow/interfaces/execute/test_torch.py
index 10b07951ac3..8d0890368f6 100644
--- a/tests/workflow/interfaces/execute/test_torch.py
+++ b/tests/workflow/interfaces/execute/test_torch.py
@@ -761,23 +761,30 @@ def cost_fn(x):
 
 
 @pytest.mark.parametrize("execute_kwargs, shots, device_name", test_matrix)
+@pytest.mark.parametrize("constructor", (qml.Hamiltonian, qml.dot, "dunders"))
 class TestHamiltonianWorkflows:
     """Test that tapes ending with expectations
     of Hamiltonians provide correct results and gradients"""
 
+    # pylint: disable=too-many-arguments, too-many-positional-arguments
     @pytest.fixture
-    def cost_fn(self, execute_kwargs, shots, device_name, seed):
+    def cost_fn(self, execute_kwargs, shots, device_name, seed, constructor):
         """Cost function for gradient tests"""
         device = get_device(device_name, seed)
 
         def _cost_fn(weights, coeffs1, coeffs2):
-            obs1 = [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliX(1), qml.PauliY(0)]
-            H1 = qml.Hamiltonian(coeffs1, obs1)
-            H1 = qml.pauli.pauli_sentence(H1).operation()
+            if constructor == "dunders":
+                H1 = (
+                    coeffs1[0] * qml.Z(0) + coeffs1[1] * qml.Z(0) @ qml.X(1) + coeffs1[2] * qml.Y(0)
+                )
+                H2 = coeffs2[0] * qml.Z(0)
+            else:
+
+                obs1 = [qml.Z(0), qml.Z(0) @ qml.X(1), qml.Y(0)]
+                H1 = constructor(coeffs1, obs1)
 
-            obs2 = [qml.PauliZ(0)]
-            H2 = qml.Hamiltonian(coeffs2, obs2)
-            H2 = qml.pauli.pauli_sentence(H2).operation()
+                obs2 = [qml.PauliZ(0)]
+                H2 = constructor(coeffs2, obs2)
 
             with qml.queuing.AnnotatedQueue() as q:
                 qml.RX(weights[0], wires=0)
@@ -854,8 +861,6 @@ def test_multiple_hamiltonians_trainable(self, execute_kwargs, cost_fn, shots):
         """Test hamiltonian with trainable parameters."""
         if execute_kwargs["diff_method"] == "adjoint":
             pytest.skip("trainable hamiltonians not supported with adjoint")
-        if execute_kwargs["diff_method"] != "backprop":
-            pytest.xfail(reason="parameter shift derivatives do not yet support sums.")
 
         coeffs1 = torch.tensor([0.1, 0.2, 0.3], requires_grad=True)
         coeffs2 = torch.tensor([0.7], requires_grad=True)
@@ -872,8 +877,7 @@ def test_multiple_hamiltonians_trainable(self, execute_kwargs, cost_fn, shots):
         res = torch.hstack(torch.autograd.functional.jacobian(cost_fn, (weights, coeffs1, coeffs2)))
         expected = self.cost_fn_jacobian(weights, coeffs1, coeffs2)
         if shots.has_partitioned_shots:
-            pytest.xfail(
-                "multiple hamiltonians with shot vectors does not seem to be differentiable."
-            )
+            assert torch.allclose(res[:2, :], expected, atol=atol_for_shots(shots), rtol=0)
+            assert torch.allclose(res[2:, :], expected, atol=atol_for_shots(shots), rtol=0)
         else:
             assert torch.allclose(res, expected, atol=atol_for_shots(shots), rtol=0)

From 3bec1cd2328658ff400358670287f09588e04941 Mon Sep 17 00:00:00 2001
From: "Yushao Chen (Jerry)" <chenys13@outlook.com>
Date: Wed, 15 Jan 2025 13:39:02 -0500
Subject: [PATCH 5/5] Remove `output_dim` property of `qml.tape.QuantumScript`
 (#6829)

**Context:**
The property `output_dim` of `QuantumScript` has been deprecated in
previous round, [sc-74594]
https://github.com/PennyLaneAI/pennylane/pull/6577. Now we would love to
remove it entirely.

**Description of the Change:**
 - Remove property `output_dim` along with its tests, warnings, etc.
 - Remove the detached attribute `_output_dim` after previous removal.
- Remove the detached private method `_update_output_dim` which used to
be strictly only helping update property `output_dim`.

**Benefits:**
Cleaner codebase with less redundancies.

**Possible Drawbacks:**
Not known. Perhaps some nuances hidden in some uncovered cornercase
tests. Anyways potential drawback won't affect any current feature in
principle.

**Related GitHub Issues:**
[sc-82150]

---------

Co-authored-by: lillian542 <38584660+lillian542@users.noreply.github.com>
---
 doc/development/deprecations.rst    | 10 ++--
 doc/releases/changelog-dev.md       |  3 +
 pennylane/devices/_legacy_device.py |  1 -
 pennylane/tape/qscript.py           | 51 +---------------
 pennylane/tape/tape.py              |  2 -
 tests/tape/test_qscript.py          | 91 +----------------------------
 6 files changed, 11 insertions(+), 147 deletions(-)

diff --git a/doc/development/deprecations.rst b/doc/development/deprecations.rst
index 2c4269de6c4..e4b8810b38d 100644
--- a/doc/development/deprecations.rst
+++ b/doc/development/deprecations.rst
@@ -25,11 +25,6 @@ Pending deprecations
   - Deprecated in v0.40
   - Will be removed in v0.41
 
-* The ``output_dim`` property of ``qml.tape.QuantumScript`` has been deprecated. Instead, use method ``shape`` of ``QuantumScript`` or ``MeasurementProcess`` to get the same information.
-
-  - Deprecated in v0.40
-  - Will be removed in v0.41
-
 * The ``QNode.get_best_method`` and ``QNode.best_method_str`` methods have been deprecated. 
   Instead, use the ``qml.workflow.get_best_diff_method`` function. 
 
@@ -84,6 +79,11 @@ for details on how to port your legacy code to the new system. The following fun
 Completed deprecation cycles
 ----------------------------
 
+* The ``output_dim`` property of ``qml.tape.QuantumScript`` has been removed. Instead, use method ``shape`` of ``QuantumScript`` or ``MeasurementProcess`` to get the same information.
+
+  - Deprecated in v0.40
+  - Removed in v0.41
+
 * The ``qml.qsvt_legacy`` function has been removed.
   Instead, use ``qml.qsvt``. The new functionality takes an input polynomial instead of angles.
 
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 25c75208671..ca233cdf4f0 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -11,6 +11,9 @@
 
 <h3>Breaking changes 💔</h3>
 
+* The `output_dim` property of `qml.tape.QuantumScript` has been removed. Instead, use method `shape` of `QuantumScript` or `MeasurementProcess` to get the same information.
+  [(#6829)](https://github.com/PennyLaneAI/pennylane/pull/6829)
+
 * Removed method `qsvt_legacy` along with its private helper `_qsp_to_qsvt`
   [(#6827)](https://github.com/PennyLaneAI/pennylane/pull/6827)
 
diff --git a/pennylane/devices/_legacy_device.py b/pennylane/devices/_legacy_device.py
index 33400c37695..b4f37a654b1 100644
--- a/pennylane/devices/_legacy_device.py
+++ b/pennylane/devices/_legacy_device.py
@@ -93,7 +93,6 @@ def _local_tape_expand(tape, depth, stop_at):
 
     # Update circuit info
     new_tape._batch_size = tape._batch_size
-    new_tape._output_dim = tape._output_dim
     return new_tape
 
 
diff --git a/pennylane/tape/qscript.py b/pennylane/tape/qscript.py
index 38cdf5ff6bc..f7f0bd2c918 100644
--- a/pennylane/tape/qscript.py
+++ b/pennylane/tape/qscript.py
@@ -25,7 +25,7 @@
 from typing import Any, Optional, TypeVar, Union
 
 import pennylane as qml
-from pennylane.measurements import MeasurementProcess, ProbabilityMP, StateMP
+from pennylane.measurements import MeasurementProcess
 from pennylane.measurements.shots import Shots, ShotsLike
 from pennylane.operation import _UNSET_BATCH_SIZE, Observable, Operation, Operator
 from pennylane.pytrees import register_pytree
@@ -183,7 +183,6 @@ def __init__(
         self._trainable_params = trainable_params
         self._graph = None
         self._specs = None
-        self._output_dim = None
         self._batch_size = _UNSET_BATCH_SIZE
 
         self._obs_sharing_wires = None
@@ -322,29 +321,6 @@ def batch_size(self) -> Optional[int]:
             self._update_batch_size()
         return self._batch_size
 
-    @property
-    def output_dim(self) -> int:
-        """The (inferred) output dimension of the quantum script.
-
-        .. warning::
-
-            ``QuantumScript.output_dim`` is being deprecated. Instead, considering
-            using method ``shape`` of ``QuantumScript`` or ``MeasurementProcess``
-            to get the same information. See ``qml.gradients.parameter_shift_cv.py::_get_output_dim``
-            for an example.
-
-        """
-        # pylint: disable=import-outside-toplevel
-        import warnings
-
-        warnings.warn(
-            "The 'output_dim' property is deprecated and will be removed in version 0.41",
-            qml.PennyLaneDeprecationWarning,
-        )
-        if self._output_dim is None:
-            self._update_output_dim()  # this will set _batch_size if it isn't already
-        return self._output_dim
-
     @property
     def diagonalizing_gates(self) -> list[Operation]:
         """Returns the gates that diagonalize the measured wires such that they
@@ -566,28 +542,6 @@ def _update_batch_size(self):
 
         self._batch_size = candidate
 
-    def _update_output_dim(self):
-        """Update the dimension of the output of the quantum script.
-
-        Sets:
-            self._output_dim (int): Size of the quantum script output (when flattened)
-
-        This method makes use of `self.batch_size`, so that `self._batch_size`
-        needs to be up to date when calling it.
-        Call `_update_batch_size` before `_update_output_dim`
-        """
-        self._output_dim = 0
-        for m in self.measurements:
-            # attempt to infer the output dimension
-            if isinstance(m, ProbabilityMP):
-                # TODO: what if we had a CV device here? Having the base as
-                # 2 would have to be swapped to the cutoff value
-                self._output_dim += 2 ** len(m.wires)
-            elif not isinstance(m, StateMP):
-                self._output_dim += 1
-        if self.batch_size:
-            self._output_dim *= self.batch_size
-
     # ========================================================
     # Parameter handling
     # ========================================================
@@ -984,9 +938,6 @@ def copy(self, copy_operations: bool = False, **update) -> "QuantumScript":
             # obs may change if measurements were updated
             new_qscript._obs_sharing_wires = self._obs_sharing_wires
             new_qscript._obs_sharing_wires_id = self._obs_sharing_wires_id
-        if not (update.get("measurements") or update.get("operations")):
-            # output_dim may change if either measurements or operations were updated
-            new_qscript._output_dim = self._output_dim
         return new_qscript
 
     def __copy__(self) -> "QuantumScript":
diff --git a/pennylane/tape/tape.py b/pennylane/tape/tape.py
index 4d08187b147..80bff4c29e5 100644
--- a/pennylane/tape/tape.py
+++ b/pennylane/tape/tape.py
@@ -294,7 +294,6 @@ def stop_at(obj):  # pylint: disable=unused-argument
 
     # Update circuit info
     new_tape._batch_size = tape._batch_size
-    new_tape._output_dim = tape._output_dim
     return new_tape
 
 
@@ -343,7 +342,6 @@ def expand_tape_state_prep(tape, skip_first=True):
 
     # Update circuit info
     new_tape._batch_size = tape._batch_size
-    new_tape._output_dim = tape._output_dim
     return new_tape
 
 
diff --git a/tests/tape/test_qscript.py b/tests/tape/test_qscript.py
index 7bad7342677..07ba6ae19ab 100644
--- a/tests/tape/test_qscript.py
+++ b/tests/tape/test_qscript.py
@@ -318,63 +318,21 @@ def test_error_inconsistent_batch_sizes(self, x, rot, y):
         ):
             _ = tape.batch_size
 
-    @pytest.mark.parametrize(
-        "m, output_dim",
-        [
-            ([qml.expval(qml.PauliX(0))], 1),
-            ([qml.expval(qml.PauliX(0)), qml.var(qml.PauliY(1))], 2),
-            ([qml.probs(wires=(0, 1))], 4),
-            ([qml.state()], 0),
-            ([qml.probs((0, 1)), qml.expval(qml.PauliX(0))], 5),
-        ],
-    )
-    @pytest.mark.parametrize("ops, factor", [([], 1), ([qml.RX([1.2, 2.3, 3.4], wires=0)], 3)])
-    def test_update_output_dim(self, m, output_dim, ops, factor):
-        """Test setting the output_dim property."""
-        qs = QuantumScript(ops, m)
-
-        with pytest.warns(
-            qml.PennyLaneDeprecationWarning, match="The 'output_dim' property is deprecated"
-        ):
-            assert qs.output_dim == output_dim * factor
-
-    def test_lazy_batch_size_and_output_dim(self):
-        """Test that batch_size and output_dim are computed lazily."""
+    def test_lazy_batch_size(self):
+        """Test that batch_size is computed lazily."""
         qs = QuantumScript([qml.RX([1.1, 2.2], 0)], [qml.expval(qml.PauliZ(0))])
         copied = qs.copy()
         assert qs._batch_size is _UNSET_BATCH_SIZE
-        assert qs._output_dim is None
         # copying did not evaluate them either
         assert copied._batch_size is _UNSET_BATCH_SIZE
-        assert copied._output_dim is None
 
         # now evaluate it
         assert qs.batch_size == 2
-        assert qs._output_dim is None  # setting batch_size didn't set output_dim
 
-        with pytest.warns(
-            qml.PennyLaneDeprecationWarning, match="The 'output_dim' property is deprecated"
-        ):
-            assert qs.output_dim == 2
         copied = qs.copy()
         assert qs._batch_size == 2
-        assert qs._output_dim == 2
         # copied tape has it pre-evaluated
         assert copied._batch_size == 2
-        assert copied._output_dim == 2
-
-    def test_lazy_setting_output_dim_sets_batch_size(self):
-        """Test that setting the output_dim also sets the batch_size."""
-        qs = QuantumScript([qml.RX([1.1, 2.2], 0)], [qml.expval(qml.PauliZ(0))])
-        assert qs._batch_size is _UNSET_BATCH_SIZE
-        assert qs._output_dim is None
-
-        with pytest.warns(
-            qml.PennyLaneDeprecationWarning, match="The 'output_dim' property is deprecated"
-        ):
-            assert qs.output_dim == 2  # getting this sets both _output_dim and _batch_size
-        assert qs._output_dim == 2
-        assert qs._batch_size == 2
 
 
 class TestIteration:
@@ -583,12 +541,6 @@ def test_shallow_copy(self):
         assert qs.data == copied_qs.data
         assert qs.shots is copied_qs.shots
 
-        # check that the output dim is identical
-        with pytest.warns(
-            qml.PennyLaneDeprecationWarning, match="The 'output_dim' property is deprecated"
-        ):
-            assert qs.output_dim == copied_qs.output_dim
-
     # pylint: disable=unnecessary-lambda
     @pytest.mark.parametrize(
         "copy_fn", [lambda tape: tape.copy(copy_operations=True), lambda tape: copy.copy(tape)]
@@ -621,12 +573,6 @@ def test_shallow_copy_with_operations(self, copy_fn):
         assert qs.data == copied_qs.data
         assert qs.shots is copied_qs.shots
 
-        # check that the output dim is identical
-        with pytest.warns(
-            qml.PennyLaneDeprecationWarning, match="The 'output_dim' property is deprecated"
-        ):
-            assert qs.output_dim == copied_qs.output_dim
-
     def test_deep_copy(self):
         """Test that deep copying a tape works, and copies all constituent data except parameters"""
         prep = [qml.BasisState(np.array([1, 0]), wires=(0, 1))]
@@ -644,12 +590,6 @@ def test_deep_copy(self):
         assert all(m1 is not m2 for m1, m2 in zip(copied_qs.measurements, qs.measurements))
         assert copied_qs.shots is qs.shots
 
-        # check that the output dim is identical
-        with pytest.warns(
-            qml.PennyLaneDeprecationWarning, match="The 'output_dim' property is deprecated"
-        ):
-            assert qs.output_dim == copied_qs.output_dim
-
         # The underlying operation data has also been copied
         assert copied_qs.operations[0].wires is not qs.operations[0].wires
 
@@ -763,10 +703,6 @@ def test_cached_properties_when_updating_operations(self):
         tape = QuantumScript(ops, measurements=[qml.counts()], shots=2500, trainable_params=[1])
 
         assert tape.batch_size == 2
-        with pytest.warns(
-            qml.PennyLaneDeprecationWarning, match="The 'output_dim' property is deprecated"
-        ):
-            assert tape.output_dim == 2
         assert tape.trainable_params == [1]
 
         new_ops = [qml.RX([1.2, 2.3, 3.4], 0)]
@@ -776,10 +712,6 @@ def test_cached_properties_when_updating_operations(self):
         assert new_tape.operations == new_ops
 
         assert new_tape.batch_size == 3
-        with pytest.warns(
-            qml.PennyLaneDeprecationWarning, match="The 'output_dim' property is deprecated"
-        ):
-            assert new_tape.output_dim == 3
         assert new_tape.trainable_params == [0]
 
     def test_cached_properties_when_updating_measurements(self):
@@ -797,10 +729,6 @@ def test_cached_properties_when_updating_measurements(self):
 
         assert tape.obs_sharing_wires == []
         assert tape.obs_sharing_wires_id == []
-        with pytest.warns(
-            qml.PennyLaneDeprecationWarning, match="The 'output_dim' property is deprecated"
-        ):
-            assert tape.output_dim == 2
         assert tape.trainable_params == [1]
 
         new_measurements = [qml.expval(qml.X(0)), qml.var(qml.Y(0))]
@@ -809,10 +737,6 @@ def test_cached_properties_when_updating_measurements(self):
         assert tape.measurements == measurements
         assert new_tape.measurements == new_measurements
 
-        with pytest.warns(
-            qml.PennyLaneDeprecationWarning, match="The 'output_dim' property is deprecated"
-        ):
-            assert new_tape.output_dim == 4
         assert new_tape.obs_sharing_wires == [qml.X(0), qml.Y(0)]
         assert new_tape.obs_sharing_wires_id == [0, 1]
         assert new_tape.trainable_params == [0, 1]
@@ -1734,14 +1658,3 @@ def test_jax_pytree_integration(qscript_type):
     assert data[3] == 3.4
     assert data[4] == 2.0
     assert qml.math.allclose(data[5], eye_mat)
-
-
-@pytest.mark.parametrize("qscript_type", (QuantumScript, qml.tape.QuantumTape))
-@pytest.mark.parametrize("shots", [None, 1, 10])
-def test_output_dim_is_deprecated(qscript_type, shots):
-    """Test that the output_dim property is deprecated."""
-    with pytest.warns(
-        qml.PennyLaneDeprecationWarning, match="The 'output_dim' property is deprecated"
-    ):
-        qscript = qscript_type([], [], shots=shots)
-        _ = qscript.output_dim