Repository: PennyLaneAI/pennylane
Pull request #5179 (2024-02-08 13:11:58)

Add grouping functionality to `Sum`
**Context:**
This PR adds support for grouping with `Sum`.

**Description of the Change:**
* Added `grouping_type` and `method` arguments to `qml.dot`, `qml.sum`, and `qml.ops.Sum` constructor.
* Added `Sum.grouping_indices` attribute which stores grouping indices in the same manner as `qml.Hamiltonian`.
* Added `Sum.compute_grouping` which updates `Sum.grouping_indices` based on the requested grouping type and method. The grouping indices refer to the list of ops returned by `Sum.terms`, not `Sum.operands`. I added a note to the docstring of the function, class, and constructor functions about this.

**Benefits:**
One step closer to `Hamiltonian`-`Sum` feature parity.

**Possible Drawbacks:**

**Related GitHub Issues:**


TODO: Add changelog entry with example.


github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
mudit2812: [sc-55328]
codecov-commenter: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5179?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
Attention: `23 lines` in your changes are missing coverage. Please review.
> Comparison is base [(`5c47649`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/5c47649381b94bd8ec12cff5d2edaa1f7c62981d?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.68% compared to head [(`0abf3e6`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5179?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.61%.

| [Files](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5179?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Patch % | Lines |
|---|---|---|
| [pennylane/ops/op\_math/sum.py](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5179?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9vcF9tYXRoL3N1bS5weQ==) | 20.68% | [23 Missing :warning: ](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5179?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) |

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff              @@
##           sumterms    #5179      +/-   ##
============================================
- Coverage     99.68%   99.61%   -0.07%     
============================================
  Files           394      394              
  Lines         36052    36081      +29     
============================================
+ Hits          35937    35943       +6     
- Misses          115      138      +23     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5179?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).


Patch:
diff --git a/doc/code/qml_pauli.rst b/doc/code/qml_pauli.rst
index 9183ec2f0..1dfe4da7b 100644
--- a/doc/code/qml_pauli.rst
+++ b/doc/code/qml_pauli.rst
@@ -106,6 +106,7 @@ We can also displace the Hamiltonian by an arbitrary amount. Here, for example,
 
 >>> H = H - np.min(np.linalg.eigvalsh(H.to_mat()))
 
+.. _graph_colouring:
 
 Graph colouring
 ---------------
@@ -147,3 +148,6 @@ the groups of observables:
  [PauliY(wires=[0]), PauliZ(wires=[1])]]
 >>> coeffs_groupings
 [[4.21], [1.43, 0.97]]
+
+For a larger example of how grouping can be used with PennyLane, check out the
+`Measurement Optimization demo <https://pennylane.ai/qml/demos/tutorial_measurement_optimize/>`_.
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 696d4d07a..e87ed6fb0 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -4,6 +4,51 @@
 
 <h3>New features since last release</h3>
 
+* `qml.ops.Sum` now supports storing grouping information. Grouping type and method can be
+  specified during construction using the `grouping_type` and `method` keyword arguments of
+  `qml.dot`, `qml.sum`, or `qml.ops.Sum`. The grouping indices are stored in `Sum.grouping_indices`.
+  [(#5179)](https://github.com/PennyLaneAI/pennylane/pull/5179)
+
+  ```python
+  import pennylane as qml
+
+  a = qml.X(0)
+  b = qml.prod(qml.X(0), qml.X(1))
+  c = qml.Z(0)
+  obs = [a, b, c]
+  coeffs = [1.0, 2.0, 3.0]
+
+  op = qml.dot(coeffs, obs, grouping_type="qwc")
+  ```
+  ```pycon
+  >>> op.grouping_indices
+  ((2,), (0, 1))
+  ```
+
+  Additionally, grouping type and method can be set or changed after construction using
+  `Sum.compute_grouping()`:
+
+  ```python
+  import pennylane as qml
+
+  a = qml.X(0)
+  b = qml.prod(qml.X(0), qml.X(1))
+  c = qml.Z(0)
+  obs = [a, b, c]
+  coeffs = [1.0, 2.0, 3.0]
+
+  op = qml.dot(coeffs, obs)
+  ```
+  ```pycon
+  >>> op.grouping_indices is None
+  True
+  >>> op.compute_grouping(grouping_type="qwc")
+  >>> op.grouping_indices
+  ((2,), (0, 1))
+  ```
+
+  Note that the grouping indices refer to the lists returned by `Sum.terms()`, not `Sum.operands`.
+
 * Added new `SpectralNormError` class to the new error tracking functionality.
   [(#5154)](https://github.com/PennyLaneAI/pennylane/pull/5154)
 
@@ -87,6 +132,7 @@
 
 This release contains contributions from (in alphabetical order):
 
+Korbinian Kottmann,
 Guillermo Alonso,
 Astral Cai,
 Amintor Dusko,
@@ -94,4 +140,5 @@ Pietropaolo Frisoni,
 Soran Jahangiri,
 Korbinian Kottmann,
 Christina Lee,
+Mudit Pandey,
 Matthew Silverman.
diff --git a/pennylane/ops/functions/dot.py b/pennylane/ops/functions/dot.py
index 951e040a0..743d868ff 100644
--- a/pennylane/ops/functions/dot.py
+++ b/pennylane/ops/functions/dot.py
@@ -29,6 +29,8 @@ def dot(
     coeffs: Sequence[Union[float, Callable]],
     ops: Sequence[Union[Operator, PauliWord, PauliSentence]],
     pauli=False,
+    grouping_type=None,
+    method="rlf",
 ) -> Union[Operator, ParametrizedHamiltonian, PauliSentence]:
     r"""Returns the dot product between the ``coeffs`` vector and the ``ops`` list of operators.
 
@@ -42,6 +44,12 @@ def dot(
             operator is used to represent the linear combination. If False, a :class:`Sum` operator
             is returned. Defaults to ``False``. Note that when ``ops`` consists solely of ``PauliWord``
             and ``PauliSentence`` instances, the function still returns a PennyLane operator when ``pauli=False``.
+        grouping_type (str): The type of binary relation between Pauli words used to compute
+            the grouping. Can be ``'qwc'``, ``'commuting'``, or ``'anticommuting'``. Note that if
+            ``pauli=True``, the grouping will be ignored.
+        method (str): The graph coloring heuristic to use in solving minimum clique cover for
+            grouping, which can be ``'lf'`` (Largest First) or ``'rlf'`` (Recursive Largest
+            First). This keyword argument is ignored if ``grouping_type`` is ``None``.
 
     Raises:
         ValueError: if the number of coefficients and operators does not match or if they are empty
@@ -49,6 +57,12 @@ def dot(
     Returns:
         Operator or ParametrizedHamiltonian: operator describing the linear combination
 
+    .. note::
+
+        If grouping is requested, the computed groupings are stored as a list of list of indices
+        in ``Sum.grouping_indices``. The indices refer to the operators and coefficients returned
+        by ``Sum.terms()``, not ``Sum.operands``, as these are not guaranteed to be equivalent.
+
     **Example**
 
     >>> coeffs = np.array([1.1, 2.2])
@@ -90,6 +104,31 @@ def dot(
         <lambda>(params_0, t) * X(0)
       + <lambda>(params_1, t) * Y(0)
     )
+
+    .. details::
+        :title: Grouping
+
+        Grouping information can be collected during construction using the ``grouping_type`` and ``method``
+        keyword arguments. For example:
+
+        .. code-block:: python
+
+            import pennylane as qml
+
+            a = qml.X(0)
+            b = qml.prod(qml.X(0), qml.X(1))
+            c = qml.Z(0)
+            obs = [a, b, c]
+            coeffs = [1.0, 2.0, 3.0]
+
+            op = qml.dot(coeffs, obs, grouping_type="qwc")
+
+        >>> op.grouping_indices
+        ((2,), (0, 1))
+
+        ``grouping_type`` can be ``"qwc"`` (qubit-wise commuting), ``"commuting"``, or ``"anticommuting"``, and
+        ``method`` can be ``"rlf"`` or ``"lf"``. To see more details about how these affect grouping, see
+        :ref:`Pauli Graph Colouring<graph_colouring>` and :func:`~pennylane.pauli.group_observables`.
     """
 
     if len(coeffs) != len(ops):
@@ -116,7 +155,11 @@ def dot(
     ops = (convert_to_opmath(op) for op in ops)
 
     operands = [op if coeff == 1 else qml.s_prod(coeff, op) for coeff, op in zip(coeffs, ops)]
-    return operands[0] if len(operands) == 1 else qml.sum(*operands)
+    return (
+        operands[0]
+        if len(operands) == 1
+        else qml.sum(*operands, grouping_type=grouping_type, method=method)
+    )
 
 
 def _dot_with_ops_and_paulis(coeffs: Sequence[float], ops: Sequence[Operator]):
diff --git a/pennylane/ops/functions/equal.py b/pennylane/ops/functions/equal.py
index 9fa19f2aa..920caed10 100644
--- a/pennylane/ops/functions/equal.py
+++ b/pennylane/ops/functions/equal.py
@@ -14,7 +14,7 @@
 """
 This module contains the qml.equal function.
 """
-# pylint: disable=too-many-arguments,too-many-return-statements
+# pylint: disable=too-many-arguments,too-many-return-statements,too-many-branches
 from collections.abc import Iterable
 from functools import singledispatch
 from typing import Union
@@ -234,6 +234,12 @@ def _equal_operators(
     ):  # clarifies cases involving PauliX/Y/Z (Observable/Operation)
         return False
 
+    if isinstance(op1, qml.Identity):
+        # All Identities are equivalent, independent of wires.
+        # We already know op1 and op2 are of the same type, so no need to check
+        # that op2 is also an Identity
+        return True
+
     if op1.arithmetic_depth != op2.arithmetic_depth:
         return False
 
diff --git a/pennylane/ops/op_math/sum.py b/pennylane/ops/op_math/sum.py
index 5874e6785..1686f1c8c 100644
--- a/pennylane/ops/op_math/sum.py
+++ b/pennylane/ops/op_math/sum.py
@@ -29,7 +29,7 @@ from pennylane.queuing import QueuingManager
 from .composite import CompositeOp
 
 
-def sum(*summands, id=None, lazy=True):
+def sum(*summands, grouping_type=None, method="rlf", id=None, lazy=True):
     r"""Construct an operator which is the sum of the given operators.
 
     Args:
@@ -39,6 +39,11 @@ def sum(*summands, id=None, lazy=True):
         id (str or None): id for the Sum operator. Default is None.
         lazy=True (bool): If ``lazy=False``, a simplification will be performed such that when any
             of the operators is already a sum operator, its operands (summands) will be used instead.
+        grouping_type (str): The type of binary relation between Pauli words used to compute
+            the grouping. Can be ``'qwc'``, ``'commuting'``, or ``'anticommuting'``.
+        method (str): The graph coloring heuristic to use in solving minimum clique cover for
+            grouping, which can be ``'lf'`` (Largest First) or ``'rlf'`` (Recursive Largest
+            First). This keyword argument is ignored if ``grouping_type`` is ``None``.
 
     Returns:
         ~ops.op_math.Sum: The operator representing the sum of summands.
@@ -56,6 +61,12 @@ def sum(*summands, id=None, lazy=True):
         >>> op = qml.sum(np.array([qml.RX(0.4, 0), qml.RZ(0.3, 0)]), qml.Z(0))
         AttributeError: 'numpy.ndarray' object has no attribute 'wires'
 
+    .. note::
+
+        If grouping is requested, the computed groupings are stored as a list of list of indices
+        in ``Sum.grouping_indices``. The indices refer to the operators and coefficients returned
+        by ``Sum.terms()``, not ``Sum.operands``, as these are not guaranteed to be equivalent.
+
     .. seealso:: :class:`~.ops.op_math.Sum`
 
     **Example**
@@ -66,13 +77,38 @@ def sum(*summands, id=None, lazy=True):
     >>> summed_op.matrix()
     array([[ 1,  1],
            [ 1, -1]])
+
+    .. details::
+        :title: Grouping
+
+        Grouping information can be collected during construction using the ``grouping_type`` and ``method``
+        keyword arguments. For example:
+
+        .. code-block:: python
+
+            import pennylane as qml
+
+            a = qml.s_prod(1.0, qml.X(0))
+            b = qml.s_prod(2.0, qml.prod(qml.X(0), qml.X(1)))
+            c = qml.s_prod(3.0, qml.Z(0))
+
+            op = qml.sum(a, b, c, grouping_type="qwc")
+
+        >>> op.grouping_indices
+        ((2,), (0, 1))
+
+        ``grouping_type`` can be ``"qwc"`` (qubit-wise commuting), ``"commuting"``, or ``"anticommuting"``, and
+        ``method`` can be ``"rlf"`` or ``"lf"``. To see more details about how these affect grouping, see
+        :ref:`Pauli Graph Colouring<graph_colouring>` and :func:`~pennylane.pauli.group_observables`.
     """
     summands = tuple(convert_to_opmath(op) for op in summands)
     if lazy:
-        return Sum(*summands, id=id)
+        return Sum(*summands, grouping_type=grouping_type, method=method, id=id)
 
     summands_simp = Sum(
         *itertools.chain.from_iterable([op if isinstance(op, Sum) else [op] for op in summands]),
+        grouping_type=grouping_type,
+        method=method,
         id=id,
     )
 
@@ -89,6 +125,11 @@ class Sum(CompositeOp):
         *summands (tuple[~.operation.Operator]): a tuple of operators which will be summed together.
 
     Keyword Args:
+        grouping_type (str): The type of binary relation between Pauli words used to compute
+            the grouping. Can be ``'qwc'``, ``'commuting'``, or ``'anticommuting'``.
+        method (str): The graph coloring heuristic to use in solving minimum clique cover for
+            grouping, which can be ``'lf'`` (Largest First) or ``'rlf'`` (Recursive Largest
+            First). This keyword argument is ignored if ``grouping_type`` is ``None``.
         id (str or None): id for the sum operator. Default is None.
 
     .. note::
@@ -106,6 +147,12 @@ class Sum(CompositeOp):
         >>> op = qml.sum(np.array([qml.RX(0.4, 0), qml.RZ(0.3, 0)]), qml.Z(0))
         AttributeError: 'numpy.ndarray' object has no attribute 'wires'
 
+    .. note::
+
+        If grouping is requested, the computed groupings are stored as a list of list of indices
+        in ``Sum.grouping_indices``. The indices refer to the operators and coefficients returned
+        by ``Sum.terms()``, not ``Sum.operands``, as these are not guaranteed to be equivalent.
+
     .. seealso:: :func:`~.ops.op_math.sum`
 
     **Example**
@@ -161,11 +208,39 @@ class Sum(CompositeOp):
     _op_symbol = "+"
     _math_op = math.sum
 
+    def _flatten(self):
+        return tuple(self.operands), (self.grouping_indices,)
+
+    @classmethod
+    def _unflatten(cls, data, metadata):
+        # pylint: disable=protected-access
+        new_op = cls(*data)
+        new_op._grouping_indices = metadata[0]
+        return new_op
+
+    def __init__(
+        self, *operands: Operator, grouping_type=None, method="rlf", id=None, _pauli_rep=None
+    ):
+        super().__init__(*operands, id=id, _pauli_rep=_pauli_rep)
+
+        self._grouping_indices = None
+        if grouping_type is not None:
+            self.compute_grouping(grouping_type=grouping_type, method=method)
+
     @property
     def hash(self):
         # Since addition is always commutative, we do not need to sort
         return hash(("Sum", frozenset(o.hash for o in self.operands)))
 
+    @property
+    def grouping_indices(self):
+        """Return the grouping indices attribute.
+
+        Returns:
+            list[list[int]]: indices needed to form groups of commuting observables
+        """
+        return self._grouping_indices
+
     def __str__(self):
         """String representation of the Sum."""
         ops = self.operands
@@ -351,6 +426,53 @@ class Sum(CompositeOp):
                 ops.append(factor)
         return coeffs, ops
 
+    def compute_grouping(self, grouping_type="qwc", method="rlf"):
+        """
+        Compute groups of operators and coefficients corresponding to commuting
+        observables of this Sum.
+
+        .. note::
+
+            If grouping is requested, the computed groupings are stored as a list of list of indices
+            in ``Sum.grouping_indices``. The indices refer to operators and coefficients returned
+            by ``Sum.terms()``, not ``Sum.operands``, as these are not guaranteed to be equivalent.
+
+        Args:
+            grouping_type (str): The type of binary relation between Pauli words used to compute
+                the grouping. Can be ``'qwc'``, ``'commuting'``, or ``'anticommuting'``.
+            method (str): The graph coloring heuristic to use in solving minimum clique cover for
+                grouping, which can be ``'lf'`` (Largest First) or ``'rlf'`` (Recursive Largest
+                First).
+
+        **Example**
+
+        .. code-block:: python
+
+            import pennylane as qml
+
+            a = qml.X(0)
+            b = qml.prod(qml.X(0), qml.X(1))
+            c = qml.Z(0)
+            obs = [a, b, c]
+            coeffs = [1.0, 2.0, 3.0]
+
+            op = qml.dot(coeffs, obs)
+
+        >>> op.grouping_indices is None
+        True
+        >>> op.compute_grouping(grouping_type="qwc")
+        >>> op.grouping_indices
+        ((2,), (0, 1))
+        """
+        if not self.pauli_rep:
+            raise ValueError("Cannot compute grouping for Sums containing non-Pauli operators.")
+
+        _, ops = self.terms()
+
+        with qml.QueuingManager.stop_recording():
+            op_groups = qml.pauli.group_observables(ops, grouping_type=grouping_type, method=method)
+        self._grouping_indices = tuple(tuple(ops.index(o) for o in group) for group in op_groups)
+
     @property
     def coeffs(self):
         r"""
@@ -360,7 +482,8 @@ class Sum(CompositeOp):
 
         .. seealso:: :attr:`~Sum.ops`, :class:`~Sum.pauli_rep`"""
         warnings.warn(
-            "Sum.coeffs is deprecated and will be removed in future releases. You can access both (coeffs, ops) via op.terms(). Also consider op.operands.",
+            "Sum.coeffs is deprecated and will be removed in future releases. You can access both "
+            "(coeffs, ops) via op.terms(). Also consider op.operands.",
             qml.PennyLaneDeprecationWarning,
         )
         coeffs, _ = self.terms()
@@ -375,7 +498,8 @@ class Sum(CompositeOp):
 
         .. seealso:: :attr:`~Sum.coeffs`, :class:`~Sum.pauli_rep`"""
         warnings.warn(
-            "Sum.ops is deprecated and will be removed in future releases. You can access both (coeffs, ops) via op.terms(). Also consider op.operands.",
+            "Sum.ops is deprecated and will be removed in future releases. You can access both "
+            "(coeffs, ops) via op.terms(). Also consider op.operands.",
             qml.PennyLaneDeprecationWarning,
         )
         _, ops = self.terms()


Test patch:
diff --git a/tests/ops/functions/test_equal.py b/tests/ops/functions/test_equal.py
index 8df050a1f..0fb184dfe 100644
--- a/tests/ops/functions/test_equal.py
+++ b/tests/ops/functions/test_equal.py
@@ -1313,6 +1313,20 @@ class TestMeasurementsEqual:
 class TestObservablesComparisons:
     """Tests comparisons between Hamiltonians, Tensors and PauliX/Y/Z operators"""
 
+    def test_identity_equal(self):
+        """Test that comparing two Identities always returns True regardless of wires"""
+        I1 = qml.Identity()
+        I2 = qml.Identity(wires=[-1])
+        I3 = qml.Identity(wires=[0, 1, 2, 3])
+        I4 = qml.Identity(wires=["a", "b"])
+
+        assert qml.equal(I1, I2)
+        assert qml.equal(I1, I3)
+        assert qml.equal(I1, I4)
+        assert qml.equal(I2, I3)
+        assert qml.equal(I2, I4)
+        assert qml.equal(I3, I4)
+
     @pytest.mark.parametrize(("H1", "H2", "res"), equal_hamiltonians)
     def test_hamiltonian_equal(self, H1, H2, res):
         """Tests that equality can be checked between Hamiltonians"""
diff --git a/tests/ops/op_math/test_sum.py b/tests/ops/op_math/test_sum.py
index f3282d5f7..f503b5a80 100644
--- a/tests/ops/op_math/test_sum.py
+++ b/tests/ops/op_math/test_sum.py
@@ -814,6 +814,38 @@ class TestProperties:
         """Test that the pauli rep gives the expected result."""
         assert op.pauli_rep == rep
 
+    def test_flatten_unflatten(self):
+        """Test that we can flatten/unflatten a Sum correctly."""
+        # pylint: disable=protected-access
+        op = Sum(qml.PauliX(0), qml.PauliZ(0))
+        new_op = Sum._unflatten(*op._flatten())
+
+        assert op == new_op
+        assert new_op.grouping_indices is None
+
+    @pytest.mark.parametrize("grouping_type", ("qwc", "commuting", "anticommuting"))
+    @pytest.mark.parametrize("method", ("lf", "rlf"))
+    def test_flatten_unflatten_with_groups(self, grouping_type, method):
+        """Test that we can flatten/unflatten a Sum correctly when grouping indices are available."""
+        # pylint: disable=protected-access
+        op = Sum(
+            qml.PauliX(0),
+            qml.s_prod(2.0, qml.PauliX(1)),
+            qml.s_prod(3.0, qml.PauliZ(0)),
+            grouping_type=grouping_type,
+            method=method,
+        )
+        new_op = Sum._unflatten(*op._flatten())
+
+        assert new_op == op
+        assert new_op.grouping_indices == op.grouping_indices
+
+        old_coeffs, old_ops = op.terms()
+        new_coeffs, new_ops = new_op.terms()
+
+        assert old_coeffs == new_coeffs
+        assert old_ops == new_ops
+
 
 class TestSimplify:
     """Test Sum simplify method and depth property."""
@@ -1272,6 +1304,156 @@ class TestArithmetic:
             assert s1.arithmetic_depth == s2.arithmetic_depth
 
 
+class TestGrouping:
+    """Test grouping functionality of Sum"""
+
+    def test_non_pauli_error(self):
+        """Test that grouping non-Pauli observables is not supported."""
+        op = Sum(qml.PauliX(0), Prod(qml.PauliZ(0), qml.PauliX(1)), qml.Hadamard(2))
+
+        with pytest.raises(
+            ValueError, match="Cannot compute grouping for Sums containing non-Pauli"
+        ):
+            op.compute_grouping()
+
+    def test_grouping_is_correct_kwarg(self):
+        """Basic test checking that grouping with a kwarg works as expected"""
+        a = qml.PauliX(0)
+        b = qml.PauliX(1)
+        c = qml.PauliZ(0)
+        obs = [a, b, c]
+        coeffs = [1.0, 2.0, 3.0]
+
+        # check with qml.dot
+        op1 = qml.dot(coeffs, obs, grouping_type="qwc")
+        assert op1.grouping_indices == ((0, 1), (2,))
+
+        # check with qml.sum
+        sprods = [qml.s_prod(c, o) for c, o in zip(coeffs, obs)]
+        op2 = qml.sum(*sprods, grouping_type="qwc")
+        assert op2.grouping_indices == ((0, 1), (2,))
+
+        # check with Sum directly
+        op3 = Sum(*sprods, grouping_type="qwc")
+        assert op3.grouping_indices == ((0, 1), (2,))
+
+    def test_grouping_is_correct_compute_grouping(self):
+        """Basic test checking that grouping with compute_grouping works as expected"""
+        a = qml.PauliX(0)
+        b = qml.PauliX(1)
+        c = qml.PauliZ(0)
+        obs = [a, b, c]
+        coeffs = [1.0, 2.0, 3.0]
+
+        op = qml.dot(coeffs, obs, grouping_type=None)
+        assert op.grouping_indices is None
+
+        op.compute_grouping(grouping_type="qwc")
+        assert op.grouping_indices == ((0, 1), (2,))
+
+    def test_grouping_does_not_alter_queue(self):
+        """Tests that grouping is invisible to the queue."""
+        a = qml.PauliX(0)
+        b = qml.PauliX(1)
+        c = qml.PauliZ(0)
+        obs = [a, b, c]
+        coeffs = [1.0, 2.0, 3.0]
+
+        with qml.queuing.AnnotatedQueue() as q:
+            op = qml.dot(coeffs, obs, grouping_type="qwc")
+
+        assert q.queue == [op]
+
+    def test_grouping_for_non_groupable_sums(self):
+        """Test that grouping is computed correctly, even if no observables commute"""
+        a = qml.PauliX(0)
+        b = qml.PauliY(0)
+        c = qml.PauliZ(0)
+        obs = [a, b, c]
+        coeffs = [1.0, 2.0, 3.0]
+
+        op = qml.dot(coeffs, obs, grouping_type="qwc")
+        assert op.grouping_indices == ((0,), (1,), (2,))
+
+    def test_grouping_method_can_be_set(self):
+        """Tests that the grouping method can be controlled by kwargs.
+        This is done by changing from default to 'rlf' and checking the result."""
+        a = qml.PauliX(0)
+        b = qml.PauliX(1)
+        c = qml.PauliZ(0)
+        obs = [a, b, c]
+        coeffs = [1.0, 2.0, 3.0]
+
+        # compute grouping during construction with qml.dot
+        op1 = qml.dot(coeffs, obs, grouping_type="qwc", method="lf")
+        assert op1.grouping_indices == ((2, 1), (0,))
+
+        # compute grouping during construction with qml.sum
+        sprods = [qml.s_prod(c, o) for c, o in zip(coeffs, obs)]
+        op2 = qml.sum(*sprods, grouping_type="qwc", method="lf")
+        assert op2.grouping_indices == ((2, 1), (0,))
+
+        # compute grouping during construction with Sum
+        op3 = Sum(*sprods, grouping_type="qwc", method="lf")
+        assert op3.grouping_indices == ((2, 1), (0,))
+
+        # compute grouping separately
+        op4 = qml.dot(coeffs, obs, grouping_type=None)
+        op4.compute_grouping(method="lf")
+        assert op4.grouping_indices == ((2, 1), (0,))
+
+    @pytest.mark.parametrize(
+        "grouping_type, grouping_indices",
+        [("commuting", ((0, 1), (2,))), ("anticommuting", ((1,), (0, 2)))],
+    )
+    def test_grouping_type_can_be_set(self, grouping_type, grouping_indices):
+        """Tests that the grouping type can be controlled by kwargs.
+        This is done by changing from default to 'commuting' or 'anticommuting'
+        and checking the result."""
+        a = qml.PauliX(0)
+        b = qml.PauliX(1)
+        c = qml.PauliZ(0)
+        obs = [a, b, c]
+        coeffs = [1.0, 2.0, 3.0]
+
+        # compute grouping during construction with qml.dot
+        op1 = qml.dot(coeffs, obs, grouping_type=grouping_type)
+        assert op1.grouping_indices == grouping_indices
+
+        # compute grouping during construction with qml.sum
+        sprods = [qml.s_prod(c, o) for c, o in zip(coeffs, obs)]
+        op2 = qml.sum(*sprods, grouping_type=grouping_type)
+        assert op2.grouping_indices == grouping_indices
+
+        # compute grouping during construction with Sum
+        op3 = Sum(*sprods, grouping_type=grouping_type)
+        assert op3.grouping_indices == grouping_indices
+
+        # compute grouping separately
+        op4 = qml.dot(coeffs, obs, grouping_type=None)
+        op4.compute_grouping(grouping_type=grouping_type)
+        assert op4.grouping_indices == grouping_indices
+
+    @pytest.mark.parametrize("shots", [None, 1000])
+    def test_grouping_integration(self, shots):
+        """Test that grouping does not impact the results of a circuit."""
+        dev = qml.device("default.qubit", shots=shots)
+
+        @qml.qnode(dev)
+        def qnode(grouping_type):
+            H = qml.dot(
+                [1.0, 2.0, 3.0],
+                [qml.X(0), qml.prod(qml.X(0), qml.X(1)), qml.prod(qml.X(0), qml.X(1), qml.X(2))],
+                grouping_type=grouping_type,
+            )
+            for i in range(3):
+                qml.Hadamard(i)
+            return qml.expval(H)
+
+        assert np.allclose(qnode("qwc"), 6.0)
+        assert np.allclose(qnode(None), 6.0)
+
+
 class TestSupportsBroadcasting:
     """Test that the Sum operator supports broadcasting if its operands support broadcasting."""