Repository: PennyLaneAI/pennylane
Pull request #1181 (2021-03-31 13:47:19)

Incorporate Pauli group functions into existing grouping module.
**Context:** Iterating over the Pauli group, and manipulation of Pauli operations at the level of words are common tasks. Pauli operators as observables are available in PennyLane, however more advanced functionality is not fully implemented.

**Description of the Change:** Adds the function pauli_group to construct the Pauli group, and pauli_mult/pauli_mult_with_phase to compute products of Pauli words. Adds other miscellaneous utility functions.

**Benefits:** Users can now quickly iterate over the full Pauli group, convert back and forth between Pauli observables and their string representation, and multiply Paulis together at the level of their words without full matrix multiplication.

**Possible Drawbacks:** None

**Related GitHub Issues:** None


codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1181?src=pr&el=h1) Report
> Merging [#1181](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1181?src=pr&el=desc) (ee3477c) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/27896b40033badf1e7f0538d9a566e20b554ba42?el=desc) (27896b4) will **increase** coverage by `0.01%`.
> The diff coverage is `99.23%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1181/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1181?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master    #1181      +/-   ##
==========================================
+ Coverage   98.12%   98.14%   +0.01%     
==========================================
  Files         144      145       +1     
  Lines       10846    10977     +131     
==========================================
+ Hits        10643    10773     +130     
- Misses        203      204       +1     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1181?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [pennylane/grouping/utils.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1181/diff?src=pr&el=tree#diff-cGVubnlsYW5lL2dyb3VwaW5nL3V0aWxzLnB5) | `99.50% <98.57%> (-0.50%)` | :arrow_down: |
| [pennylane/grouping/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1181/diff?src=pr&el=tree#diff-cGVubnlsYW5lL2dyb3VwaW5nL19faW5pdF9fLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/grouping/pauli.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1181/diff?src=pr&el=tree#diff-cGVubnlsYW5lL2dyb3VwaW5nL3BhdWxpLnB5) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1181?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1181?src=pr&el=footer). Last update [27896b4...ee3477c](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1181?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

antalszava: @glassnotes nice! :1st_place_medal: overall looks good, it seems that codecov found a line that is not covered `pennylane/grouping/utils.py#L450`. Was that something newly added or was just missed before and rearranged now?
glassnotes: > @glassnotes nice! 1st_place_medal overall looks good, it seems that codecov found a line that is not covered `pennylane/grouping/utils.py#L450`. Was that something newly added or was just missed before and rearranged now?

I'm not sure, but I've updated the tests to get it incorporated now. Thanks for the heads up!

Patch:
diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index cd9dcb07c..a86ce0254 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -2,6 +2,60 @@
 
 <h3>New features since last release</h3>
 
+* Added functionality for constructing and manipulating the Pauli group
+  [(#1181)](https://github.com/PennyLaneAI/pennylane/pull/1181).
+  The function `pennylane.grouping.pauli_group` provides a generator to
+  easily loop over the group, or construct and store it in its entirety.
+  For example, we can construct the 3-qubit Pauli group like so:
+
+  ```pycon
+  >>> from pennylane.grouping import pauli_group
+  >>> pauli_group_1_qubit = list(pauli_group(1))
+  >>> pauli_group_1_qubit
+  [Identity(wires=[0]), PauliZ(wires=[0]), PauliX(wires=[0]), PauliY(wires=[0])]
+  ```
+
+  We can multiply together its members at the level of Pauli words
+  using the `pauli_mult` and `pauli_multi_with_phase` functions.
+  This can be done on arbitrarily-labeled wires as well, by defining a wire map.
+
+  ```pycon
+  >>> from pennylane.grouping import pauli_group, pauli_mult
+  >>> wire_map = {'a' : 0, 'b' : 1, 'c' : 2}
+  >>> pg = list(pauli_group(3, wire_map=wire_map))
+  >>> pg[3]
+  PauliZ(wires=['b']) @ PauliZ(wires=['c'])
+  >>> pg[55]
+  PauliY(wires=['a']) @ PauliY(wires=['b']) @ PauliZ(wires=['c'])
+  >>> pauli_mult(pg[3], pg[55], wire_map=wire_map)
+  PauliY(wires=['a']) @ PauliX(wires=['b'])
+  ```
+
+  Functions for conversion of Pauli observables to strings (and back),
+  are included.
+
+  ```pycon
+  >>> from pennylane.grouping import pauli_word_to_string, string_to_pauli_word
+  >>> pauli_word_to_string(pg[55], wire_map=wire_map)
+  'YYZ'
+  >>> string_to_pauli_word('ZXY', wire_map=wire_map)
+  PauliZ(wires=['a']) @ PauliX(wires=['b']) @ PauliY(wires=['c'])
+  ```
+
+  Calculation of the matrix representation for arbitrary Paulis and wire maps is now
+  also supported.
+
+  ```pycon
+  >>> from pennylane.grouping import pauli_word_to_matrix
+  >>> wire_map = {'a' : 0, 'b' : 1}
+  >>> pauli_word = qml.PauliZ('b')  # corresponds to Pauli 'IZ'
+  >>> pauli_word_to_matrix(pauli_word, wire_map=wire_map)
+  array([[ 1.,  0.,  0.,  0.],
+         [ 0., -1.,  0., -0.],
+         [ 0.,  0.,  1.,  0.],
+         [ 0., -0.,  0., -1.]])
+  ```
+
 * TBD: Add `qml.qaoa.cycle` (when structure is finalised)
   [(#1207)](https://github.com/PennyLaneAI/pennylane/pull/1207)
   [(#1209)](https://github.com/PennyLaneAI/pennylane/pull/1209)
diff --git a/doc/code/qml_grouping.rst b/doc/code/qml_grouping.rst
index 1b441cffd..f1699c49d 100644
--- a/doc/code/qml_grouping.rst
+++ b/doc/code/qml_grouping.rst
@@ -1,13 +1,69 @@
 qml.grouping
 ============
-
-This subpackage defines functions and classes for Pauli-word partitioning
+This subpackage defines functions and classes for generating and manipulating
+elements of the Pauli group. It also contains Pauli-word partitioning
 functionality used in measurement optimization.
 
-A Pauli word is defined as :math:`P_J = \prod_{i=1}^{N}\sigma_i^{(J)}`, where
-:math:`\sigma_i^{(J)}` is one of the Pauli operators (:class:`~pennylane.PauliX`,
-:class:`~pennylane.PauliY`, :class:`~pennylane.PauliZ`) or identity
-(:class:`~pennylane.Identity`) acting on the :math:`i^{th}` qubit.
+
+Pauli group
+-----------
+
+The single-qubit Pauli group consists of the four single-qubit Pauli operations
+:class:`~pennylane.Identity`, :class:`~pennylane.PauliX`,
+:class:`~pennylane.PauliY` , and :class:`~pennylane.PauliZ`. The :math:`n`-qubit
+Pauli group is constructed by taking all possible :math:`N`-fold tensor products
+of these elements. Elements of the :math:`n`-qubit Pauli group are often known
+as Pauli words, and have the form :math:`P_J = \otimes_{i=1}^{n}\sigma_i^{(J)}`,
+where :math:`\sigma_i^{(J)}` is one of the Pauli operators
+(:class:`~pennylane.PauliX`, :class:`~pennylane.PauliY`,
+:class:`~pennylane.PauliZ`) or identity (:class:`~pennylane.Identity`) acting on
+the :math:`i^{th}` qubit. The full :math:`n`-qubit Pauli group has size
+:math:`4^n` (neglecting the four possible global phases that may arise from
+multiplication of its elements).
+
+The Pauli group can be constructed using the :func:`~.pauli_group`
+function. To construct the group, it is recommended to provide a wire map in
+order to indicate the names and indices of the wires. (If none is provided, a
+default mapping of integers will be used.)
+
+>>> from pennylane.grouping import pauli_group
+>>> pg_3 = list(pauli_group(3))
+
+Multiplication of Pauli group elements can be performed using
+:func:`~.pauli_mult` or
+:func:`~.pauli_mult_with_phase`:
+
+>>> from pennylane.grouping import pauli_mult
+>>> wire_map = {'a' : 0, 'b' : 1, 'c' : 2}
+>>> pg = list(pauli_group(3, wire_map=wire_map))
+>>> pg[3]
+PauliZ(wires=['b']) @ PauliZ(wires=['c'])
+>>> pg[55]
+PauliY(wires=['a']) @ PauliY(wires=['b']) @ PauliZ(wires=['c'])
+>>> pauli_mult(pg[3], pg[55], wire_map=wire_map)
+PauliY(wires=['a']) @ PauliX(wires=['b'])
+
+Pauli observables can be converted to strings (and vice versa):
+
+>>> from pennylane.grouping import pauli_word_to_string, string_to_pauli_word
+>>> pauli_word_to_string(pg[55], wire_map=wire_map)
+'YYZ'
+>>> string_to_pauli_word('ZXY', wire_map=wire_map)
+PauliZ(wires=['a']) @ PauliX(wires=['b']) @ PauliY(wires=['c'])
+
+The matrix representation for arbitrary Paulis and wire maps can also be performed.
+
+>>> pennylane.grouping import pauli_word_to_matrix
+>>> wire_map = {'a' : 0, 'b' : 1}
+>>> pauli_word = qml.PauliZ('b')  # corresponds to Pauli 'IZ'
+>>> pauli_word_to_matrix(pauli_word, wire_map=wire_map)
+array([[ 1.,  0.,  0.,  0.],
+       [ 0., -1.,  0., -0.],
+       [ 0.,  0.,  1.,  0.],
+       [ 0., -0.,  0., -1.]])
+
+Grouping observables
+--------------------
 
 Pauli words can be used for expressing a qubit :class:`~pennylane.Hamiltonian`.
 A qubit Hamiltonian has the form :math:`H_{q} = \sum_{J} C_J P_J` where
diff --git a/pennylane/grouping/__init__.py b/pennylane/grouping/__init__.py
index 458dd94c3..36cac90ab 100644
--- a/pennylane/grouping/__init__.py
+++ b/pennylane/grouping/__init__.py
@@ -30,7 +30,12 @@ from .utils import (
     are_identical_pauli_words,
     pauli_to_binary,
     binary_to_pauli,
+    pauli_word_to_string,
+    string_to_pauli_word,
+    pauli_word_to_matrix,
+    is_commuting,
     is_qwc,
     observables_to_binary_matrix,
     qwc_complement_adj_matrix,
 )
+from .pauli import pauli_group, pauli_mult, pauli_mult_with_phase
diff --git a/pennylane/grouping/pauli.py b/pennylane/grouping/pauli.py
new file mode 100644
index 000000000..f314c4a09
--- /dev/null
+++ b/pennylane/grouping/pauli.py
@@ -0,0 +1,252 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Functions for constructing the :math:`n`-qubit Pauli group, and performing the
+group operation (multiplication).
+"""
+
+import numpy as np
+
+from pennylane import Identity
+from pennylane.grouping.utils import (
+    binary_to_pauli,
+    pauli_to_binary,
+    pauli_word_to_string,
+    are_identical_pauli_words,
+    _wire_map_from_pauli_pair,
+)
+
+
+def _pauli_group_generator(n_qubits, wire_map=None):
+    """Generator function for the Pauli group.
+
+    This function is called by ``pauli_group`` in order to actually generate the
+    group elements. They are split so that the outer function can handle input
+    validation, while this generator is responsible for performing the actual
+    operations.
+
+    Args:
+        n_qubits (int): The number of qubits for which to create the group.
+        wire_map (dict[Union[str, int], int]): dictionary containing all wire labels
+            used in the Pauli word as keys, and unique integer labels as their values.
+            If no wire map is provided, wires will be labeled by consecutive integers between :math:`0` and ``n_qubits``.
+
+    Returns:
+        .Operation: The next Pauli word in the group.
+    """
+
+    element_idx = 0
+
+    if not wire_map:
+        wire_map = {wire_idx: wire_idx for wire_idx in range(n_qubits)}
+
+    while element_idx < 4 ** n_qubits:
+        binary_string = format(element_idx, f"#0{2*n_qubits+2}b")[2:]
+        binary_vector = [float(b) for b in binary_string]
+        yield binary_to_pauli(binary_vector, wire_map=wire_map)
+        element_idx += 1
+
+
+def pauli_group(n_qubits, wire_map=None):
+    """Generate the :math:`n`-qubit Pauli group.
+
+    This function enables the construction of the :math:`n`-qubit Pauli group with no
+    storage involved.  The :math:`n`-qubit Pauli group has size :math:`4^n`,
+    thus it may not be desirable to construct it in full and store.
+
+    The order of iteration is based on the binary symplectic representation of
+    the Pauli group as :math:`2n`-bit strings. Ordering is done by converting
+    the integers :math:`0` to :math:`2^{2n}` to binary strings, and converting those
+    strings to Pauli operators using the :func:`~.binary_to_pauli` method.
+
+    Args:
+        n_qubits (int): The number of qubits for which to create the group.
+        wire_map (dict[Union[str, int], int]): dictionary containing all wire labels
+            used in the Pauli word as keys, and unique integer labels as their values.
+            If no wire map is provided, wires will be labeled by integers between 0 and ``n_qubits``.
+
+    Returns:
+        .Operation: The next Pauli word in the group.
+
+    **Example**
+
+    The ``pauli_group`` generator can be used to loop over the Pauli group as follows.
+    (Note: in the example below, we display only the first 5 elements for brevity.)
+
+    >>> from pennylane.grouping import pauli_group
+    >>> n_qubits = 3
+    >>> for p in pauli_group(n_qubits):
+    ...     print(p)
+    ...
+    Identity(wires=[0])
+    PauliZ(wires=[2])
+    PauliZ(wires=[1])
+    PauliZ(wires=[1]) @ PauliZ(wires=[2])
+    PauliZ(wires=[0])
+
+    The full Pauli group can then be obtained like so:
+
+    >>> full_pg = list(pauli_group(n_qubits))
+
+    The group can also be created using a custom wire map; if no map is
+    specified, a default map of label :math:`i` to wire ``i`` as in the example
+    above will be created. (Note: in the example below, we display only the first
+    5 elements for brevity.)
+
+    >>> wire_map = {'a' : 0, 'b' : 1, 'c' : 2}
+    >>> for p in pauli_group(n_qubits, wire_map=wire_map):
+    ...     print(p)
+    ...
+    Identity(wires=['a'])
+    PauliZ(wires=['c'])
+    PauliZ(wires=['b'])
+    PauliZ(wires=['b']) @ PauliZ(wires=['c'])
+    PauliZ(wires=['a'])
+
+    """
+    # Cover the case where n_qubits may be passed as a float
+    if isinstance(n_qubits, float):
+        if n_qubits.is_integer():
+            n_qubits = int(n_qubits)
+
+    # If not an int, or a float representing a int, raise an error
+    if not isinstance(n_qubits, int):
+        raise TypeError("Must specify an integer number of qubits to construct the Pauli group.")
+
+    if n_qubits <= 0:
+        raise ValueError("Number of qubits must be at least 1 to construct Pauli group.")
+
+    return _pauli_group_generator(n_qubits, wire_map=wire_map)
+
+
+def pauli_mult(pauli_1, pauli_2, wire_map=None):
+    """Multiply two Pauli words together and return the product as a Pauli word.
+
+    Two Pauli operations can be multiplied together by taking the additive
+    OR of their binary symplectic representations.
+
+    Args:
+        pauli_1 (.Operation): A Pauli word.
+        pauli_2 (.Operation): A Pauli word to multiply with the first one.
+        wire_map (dict[Union[str, int], int]): dictionary containing all wire labels used in the Pauli
+            word as keys, and unique integer labels as their values. If no wire map is
+            provided, the map will be constructed from the set of wires acted on
+            by the input Pauli words.
+
+    Returns:
+        .Operation: The product of pauli_1 and pauli_2 as a Pauli word
+        (ignoring the global phase).
+
+    **Example**
+
+    This function enables multiplication of Pauli group elements at the level of
+    Pauli words, rather than matrices. For example,
+
+    >>> from pennylane.grouping import pauli_mult
+    >>> pauli_1 = qml.PauliX(0) @ qml.PauliZ(1)
+    >>> pauli_2 = qml.PauliY(0) @ qml.PauliZ(1)
+    >>> product = pauli_mult(pauli_1, pauli_2)
+    >>> print(product)
+    PauliZ(wires=[0])
+    """
+
+    if wire_map is None:
+        wire_map = _wire_map_from_pauli_pair(pauli_1, pauli_2)
+
+    # Check if pauli_1 and pauli_2 are the same; if so, the result is the Identity
+    if are_identical_pauli_words(pauli_1, pauli_2):
+        first_wire = list(pauli_1.wires)[0]
+        return Identity(first_wire)
+
+    # Compute binary symplectic representations
+    pauli_1_binary = pauli_to_binary(pauli_1, wire_map=wire_map)
+    pauli_2_binary = pauli_to_binary(pauli_2, wire_map=wire_map)
+
+    bin_symp_1 = np.array([int(x) for x in pauli_1_binary])
+    bin_symp_2 = np.array([int(x) for x in pauli_2_binary])
+
+    # Shorthand for bitwise XOR of numpy arrays
+    pauli_product = bin_symp_1 ^ bin_symp_2
+
+    return binary_to_pauli(pauli_product, wire_map=wire_map)
+
+
+def pauli_mult_with_phase(pauli_1, pauli_2, wire_map=None):
+    r"""Multiply two Pauli words together, and return both their product as a Pauli word
+    and the global phase.
+
+    Two Pauli operations can be multiplied together by taking the additive
+    OR of their binary symplectic representations. The phase is computed by
+    looking at the number of times we have the products  :math:`XY, YZ`, or :math:`ZX` (adds a
+    phase of :math:`i`), or :math:`YX, ZY, XZ` (adds a phase of :math:`-i`).
+
+    Args:
+        pauli_1 (.Operation): A Pauli word.
+        pauli_2 (.Operation): A Pauli word to multiply with the first one.
+        wire_map  (dict[Union[str, int], int]): dictionary containing all wire labels used in the Pauli
+            word as keys, and unique integer labels as their values. If no wire map is
+            provided, the map will be constructed from the set of wires acted on
+            by the input Pauli words.
+
+    Returns:
+        tuple[.Operation, complex]: The product of ``pauli_1`` and ``pauli_2``, and the
+        global phase.
+
+    **Example**
+
+    This function works the same as :func:`~.pauli_mult` but also returns the global
+    phase accumulated as a result of the ordering of Paulis in the product (e.g., :math:`XY = iZ`,
+    and :math:`YX = -iZ`).
+
+    >>> from pennylane.grouping import pauli_mult_with_phase
+    >>> pauli_1 = qml.PauliX(0) @ qml.PauliZ(1)
+    >>> pauli_2 = qml.PauliY(0) @ qml.PauliZ(1)
+    >>> product, phase = pauli_mult_with_phase(pauli_1, pauli_2)
+    >>> product
+    PauliZ(wires=[0])
+    >>> phase
+    1j
+    """
+
+    if wire_map is None:
+        wire_map = _wire_map_from_pauli_pair(pauli_1, pauli_2)
+
+    # Get the product; use our earlier function
+    pauli_product = pauli_mult(pauli_1, pauli_2, wire_map)
+
+    # Get the names of the operations; in cases where only one single-qubit Pauli
+    # is present, the operation name is stored as a string rather than a list, so convert it
+    pauli_1_string = pauli_word_to_string(pauli_1, wire_map=wire_map)
+    pauli_2_string = pauli_word_to_string(pauli_2, wire_map=wire_map)
+
+    pos_phases = (("X", "Y"), ("Y", "Z"), ("Z", "X"))
+
+    phase = 1
+
+    for qubit_1_char, qubit_2_char in zip(pauli_1_string, pauli_2_string):
+        # If we have identities anywhere we don't pick up a phase
+        if qubit_1_char == "I" or qubit_2_char == "I":
+            continue
+
+        # Likewise, no additional phase if the Paulis are the same
+        if qubit_1_char == qubit_2_char:
+            continue
+
+        # Use Pauli commutation rules to determine the phase
+        if (qubit_1_char, qubit_2_char) in pos_phases:
+            phase *= 1j
+        else:
+            phase *= -1j
+
+    return pauli_product, phase
diff --git a/pennylane/grouping/utils.py b/pennylane/grouping/utils.py
index fef471784..13890a1ef 100644
--- a/pennylane/grouping/utils.py
+++ b/pennylane/grouping/utils.py
@@ -20,12 +20,33 @@ representation of Pauli words and applications, see:
 * `arXiv:1701.08213 <https://arxiv.org/abs/1701.08213>`_
 * `arXiv:1907.09386 <https://arxiv.org/abs/1907.09386>`_
 """
+from functools import reduce
 
+import pennylane as qml
 from pennylane import PauliX, PauliY, PauliZ, Identity
 from pennylane.operation import Observable, Tensor
 from pennylane.wires import Wires
 import numpy as np
 
+# To make this quicker later on
+ID_MAT = np.eye(2)
+
+
+def _wire_map_from_pauli_pair(pauli_word_1, pauli_word_2):
+    """Generate a wire map from the union of wires of two Paulis.
+
+    Args:
+        pauli_word_1 (.Operation): A Pauli word.
+        pauli_word_2 (.Operation): A second Pauli word.
+
+    Returns:
+        dict[Union[str, int], int]): dictionary containing all wire labels used
+        in the Pauli word as keys, and unique integer labels as their values.
+    """
+    wire_labels = Wires.all_wires([pauli_word_1.wires, pauli_word_2.wires]).labels
+    wire_map = {label: i for i, label in enumerate(wire_labels)}
+    return wire_map
+
 
 def is_pauli_word(observable):
     """
@@ -333,6 +354,240 @@ def binary_to_pauli(binary_vector, wire_map=None):  # pylint: disable=too-many-b
     return pauli_word
 
 
+def pauli_word_to_string(pauli_word, wire_map=None):
+    """Convert a Pauli word from a tensor to a string.
+
+    Given a Pauli in observable form, convert it into string of
+    characters from ``['I', 'X', 'Y', 'Z']``. This representation is required for
+    functions such as :class:`.PauliRot`.
+
+    Args:
+        pauli_word (Observable): an observable, either a :class:`~.Tensor` instance or
+            single-qubit observable representing a Pauli group element.
+        wire_map (dict[Union[str, int], int]): dictionary containing all wire labels used in
+            the Pauli word as keys, and unique integer labels as their values
+
+    Returns:
+        str: The string representation of the observable in terms of ``'I'``, ``'X'``, ``'Y'``,
+        and/or ``'Z'``.
+
+    Raises:
+        TypeError: if the input observable is not a proper Pauli word.
+
+    **Example**
+
+    >>> wire_map = {'a' : 0, 'b' : 1, 'c' : 2}
+    >>> pauli_word = qml.PauliX('a') @ qml.PauliY('c')
+    >>> pauli_word_to_string(pauli_word, wire_map=wire_map)
+    'XIY'
+    """
+
+    if not is_pauli_word(pauli_word):
+        raise TypeError(f"Expected Pauli word observables, instead got {pauli_word}")
+
+    character_map = {"Identity": "I", "PauliX": "X", "PauliY": "Y", "PauliZ": "Z"}
+
+    # If there is no wire map, we must infer from the structure of Paulis
+    if wire_map is None:
+        wire_map = {pauli_word.wires.labels[i]: i for i in range(len(pauli_word.wires))}
+
+    n_qubits = len(wire_map)
+
+    # Set default value of all characters to identity
+    pauli_string = ["I"] * n_qubits
+
+    # Special case is when there is a single Pauli term
+    if not isinstance(pauli_word.name, list):
+        if pauli_word.name != "Identity":
+            wire_idx = wire_map[pauli_word.wires[0]]
+            pauli_string[wire_idx] = character_map[pauli_word.name]
+        return "".join(pauli_string)
+
+    for name, wire_label in zip(pauli_word.name, pauli_word.wires):
+        wire_idx = wire_map[wire_label]
+        pauli_string[wire_idx] = character_map[name]
+
+    return "".join(pauli_string)
+
+
+def string_to_pauli_word(pauli_string, wire_map=None):
+    """Convert a string in terms of ``'I'``, ``'X'``, ``'Y'``, and ``'Z'`` into a Pauli word
+    for the given wire map.
+
+    Args:
+        pauli_string (str): A string of characters consisting of ``'I'``, ``'X'``, ``'Y'``, and ``'Z'``
+            indicating a Pauli word.
+        wire_map (dict[Union[str, int], int]): dictionary containing all wire labels used in
+            the Pauli word as keys, and unique integer labels as their values
+
+    Returns:
+        .Observable: The Pauli word representing of ``pauli_string`` on the wires
+        enumerated in the wire map.
+
+    **Example**
+
+    >>> wire_map = {'a' : 0, 'b' : 1, 'c' : 2}
+    >>> string_to_pauli_word('XIY', wire_map=wire_map)
+    PauliX(wires=['a']) @ PauliY(wires=['c'])
+    """
+    character_map = {"I": Identity, "X": PauliX, "Y": PauliY, "Z": PauliZ}
+
+    if not isinstance(pauli_string, str):
+        raise TypeError(f"Input to string_to_pauli_word must be string, obtained {pauli_string}")
+
+    # String can only consist of I, X, Y, Z
+    if any(char not in character_map.keys() for char in pauli_string):
+        raise ValueError(
+            "Invalid characters encountered in string_to_pauli_word "
+            f"string {pauli_string}. Permitted characters are 'I', 'X', 'Y', and 'Z'"
+        )
+
+    # If no wire map is provided, construct one using integers based on the length of the string
+    if wire_map is None:
+        wire_map = {x: x for x in range(len(pauli_string))}
+
+    if len(pauli_string) != len(wire_map):
+        raise ValueError(
+            "Wire map and pauli_string must have the same length to convert "
+            "from string to Pauli word."
+        )
+
+    # Special case: all-identity Pauli
+    if pauli_string == "I" * len(wire_map):
+        first_wire = list(wire_map.keys())[0]
+        return Identity(wire_map[first_wire])
+
+    pauli_word = None
+
+    for wire_name, wire_idx in wire_map.items():
+        pauli_char = pauli_string[wire_idx]
+
+        # Don't care about the identity
+        if pauli_char == "I":
+            continue
+
+        if pauli_word is not None:
+            pauli_word = pauli_word @ character_map[pauli_char](wire_name)
+        else:
+            pauli_word = character_map[pauli_char](wire_name)
+
+    return pauli_word
+
+
+def pauli_word_to_matrix(pauli_word, wire_map=None):
+    """Convert a Pauli word from a tensor to its matrix representation.
+
+    The matrix representation of a Pauli word has dimension :math:`2^n \\times 2^n`,
+    where :math:`n` is the number of qubits provided in ``wire_map``. For wires
+    that the Pauli word does not act on, identities must be inserted into the tensor
+    product at the correct positions.
+
+    Args:
+        pauli_word (Observable): an observable, either a :class:`~.Tensor` instance or
+            single-qubit observable representing a Pauli group element.
+        wire_map (dict[Union[str, int], int]): dictionary containing all wire labels used in
+            the Pauli word as keys, and unique integer labels as their values
+
+    Returns:
+        array[complex]: The matrix representation of the multi-qubit Pauli over the
+        specified wire map.
+
+    Raises:
+        TypeError: if the input observable is not a proper Pauli word.
+
+    **Example**
+
+    >>> wire_map = {'a' : 0, 'b' : 1}
+    >>> pauli_word = qml.PauliX('a') @ qml.PauliY('b')
+    >>> pauli_word_to_matrix(pauli_word, wire_map=wire_map)
+    array([[0.+0.j, 0.-0.j, 0.+0.j, 0.-1.j],
+           [0.+0.j, 0.+0.j, 0.+1.j, 0.+0.j],
+           [0.+0.j, 0.-1.j, 0.+0.j, 0.-0.j],
+           [0.+1.j, 0.+0.j, 0.+0.j, 0.+0.j]])
+    """
+    if not is_pauli_word(pauli_word):
+        raise TypeError(f"Expected Pauli word observables, instead got {pauli_word}")
+
+    # If there is no wire map, we must infer from the structure of Paulis
+    if wire_map is None:
+        wire_map = {pauli_word.wires.labels[i]: i for i in range(len(pauli_word.wires))}
+
+    n_qubits = len(wire_map)
+
+    # If there is only a single qubit, we can return the matrix directly
+    if n_qubits == 1:
+        return pauli_word.matrix
+
+    # There may be more than one qubit in the Pauli but still only
+    # one of them with anything acting on it, so take that into account
+    pauli_names = [pauli_word.name] if isinstance(pauli_word.name, str) else pauli_word.name
+
+    # Special case: the identity Pauli
+    if pauli_names == ["Identity"]:
+        return np.eye(2 ** n_qubits)
+
+    # If there is more than one qubit, we must go through the wire map wire
+    # by wire and pick out the relevant matrices
+    pauli_mats = [ID_MAT for x in range(n_qubits)]
+
+    for wire_label, wire_idx in wire_map.items():
+        if wire_label in pauli_word.wires.labels:
+            op_idx = pauli_word.wires.labels.index(wire_label)
+            pauli_mats[wire_idx] = getattr(qml, pauli_names[op_idx]).matrix
+
+    return reduce(np.kron, pauli_mats)
+
+
+def is_commuting(pauli_word_1, pauli_word_2, wire_map=None):
+    r"""Checks if two Pauli words commute.
+
+    To determine if two Pauli words commute, we can check the value of the
+    symplectic inner product of their binary vector representations.
+    For two binary vectors representing Pauli words, :math:`p_1 = [x_1, z_1]`
+    and :math:`p_2 = [x_2, z_2],` the symplectic inner product is defined as
+    :math:`\langle p_1, p_2 \rangle_{symp} = z_1 x_2^T + z_2 x_1^T`. If the symplectic
+    product is :math:`0` they commute, while if it is :math:`1`, they don't commute.
+
+    Args:
+        pauli_word_1 (Observable): first Pauli word in commutator
+        pauli_word_2 (Observable): second Pauli word in commutator
+        wire_map (dict[Union[str, int], int]): dictionary containing all wire labels used in
+            the Pauli word as keys, and unique integer labels as their values
+
+    Returns:
+        bool: returns True if the input Pauli commute, False otherwise
+
+    Raises:
+        TypeError: if either of the Pauli words is not valid.
+
+    **Example**
+
+    >>> wire_map = {'a' : 0, 'b' : 1, 'c' : 2}
+    >>> pauli_word_1 = qml.PauliX('a') @ qml.PauliY('b')
+    >>> pauli_word_2 = qml.PauliZ('a') @ qml.PauliZ('c')
+    >>> is_commuting(pauli_word_1, pauli_word_2, wire_map=wire_map)
+    False
+    """
+
+    if not (is_pauli_word(pauli_word_1) and is_pauli_word(pauli_word_2)):
+        raise TypeError(
+            f"Expected Pauli word observables, instead got {pauli_word_1} and {pauli_word_2}"
+        )
+
+    if wire_map is None:
+        wire_map = _wire_map_from_pauli_pair(pauli_word_1, pauli_word_2)
+
+    n_qubits = len(wire_map)
+
+    pauli_vec_1 = pauli_to_binary(pauli_word_1, n_qubits=n_qubits, wire_map=wire_map)
+    pauli_vec_2 = pauli_to_binary(pauli_word_2, n_qubits=n_qubits, wire_map=wire_map)
+
+    x1, z1 = pauli_vec_1[:n_qubits], pauli_vec_1[n_qubits:]
+    x2, z2 = pauli_vec_2[:n_qubits], pauli_vec_2[n_qubits:]
+
+    return (np.dot(z1, x2) + np.dot(z2, x1)) % 2 == 0
+
+
 def is_qwc(pauli_vec_1, pauli_vec_2):
     """Checks if two Pauli words in the binary vector representation are qubit-wise commutative.
 


Test patch:
diff --git a/tests/grouping/test_graph_colouring.py b/tests/grouping/test_graph_colouring.py
index 359d42a75..b7adb7228 100644
--- a/tests/grouping/test_graph_colouring.py
+++ b/tests/grouping/test_graph_colouring.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2020 Xanadu Quantum Technologies Inc.
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
 
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
diff --git a/tests/grouping/test_group_observables.py b/tests/grouping/test_group_observables.py
index b756727f2..72156c8dd 100644
--- a/tests/grouping/test_group_observables.py
+++ b/tests/grouping/test_group_observables.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2020 Xanadu Quantum Technologies Inc.
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
 
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
diff --git a/tests/grouping/test_grouping_utils.py b/tests/grouping/test_grouping_utils.py
index 797a293c4..0906a419e 100644
--- a/tests/grouping/test_grouping_utils.py
+++ b/tests/grouping/test_grouping_utils.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2020 Xanadu Quantum Technologies Inc.
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
 
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -16,6 +16,7 @@ Unit tests for the :mod:`grouping` utility functions in ``grouping/utils.py``.
 """
 import pytest
 import numpy as np
+import pennylane as qml
 from pennylane import Identity, PauliX, PauliY, PauliZ, Hadamard, Hermitian, U3
 from pennylane.operation import Tensor
 from pennylane.wires import Wires
@@ -24,6 +25,10 @@ from pennylane.grouping.utils import (
     are_identical_pauli_words,
     pauli_to_binary,
     binary_to_pauli,
+    pauli_word_to_string,
+    string_to_pauli_word,
+    pauli_word_to_matrix,
+    is_commuting,
     is_qwc,
     observables_to_binary_matrix,
     qwc_complement_adj_matrix,
@@ -294,3 +299,178 @@ class TestGroupingUtils:
 
         with pytest.raises(ValueError, match="Expected a binary array, instead got"):
             qwc_complement_adj_matrix(not_binary_observables)
+
+    @pytest.mark.parametrize(
+        "pauli_word,wire_map,expected_string",
+        [
+            (PauliX(0), {0: 0}, "X"),
+            (Identity(0), {0: 0}, "I"),
+            (PauliZ(0) @ PauliY(1), {0: 0, 1: 1}, "ZY"),
+            (PauliX(1), {0: 0, 1: 1}, "IX"),
+            (PauliX(1), None, "X"),
+            (PauliX(1), {1: 0, 0: 1}, "XI"),
+            (PauliZ("a") @ PauliY("b") @ PauliZ("d"), {"a": 0, "b": 1, "c": 2, "d": 3}, "ZYIZ"),
+            (PauliZ("a") @ PauliY("b") @ PauliZ("d"), None, "ZYZ"),
+            (PauliX("a") @ PauliY("b") @ PauliZ("d"), {"d": 0, "c": 1, "b": 2, "a": 3}, "ZIYX"),
+        ],
+    )
+    def test_pauli_word_to_string(self, pauli_word, wire_map, expected_string):
+        """Test that Pauli words are correctly converted into strings."""
+        obtained_string = pauli_word_to_string(pauli_word, wire_map)
+        assert obtained_string == expected_string
+
+    @pytest.mark.parametrize("non_pauli_word", non_pauli_words)
+    def test_pauli_word_to_string_invalid_input(self, non_pauli_word):
+        """Ensure invalid inputs are handled properly when converting Pauli words to strings."""
+        with pytest.raises(TypeError):
+            pauli_word_to_string(non_pauli_word)
+
+    @pytest.mark.parametrize(
+        "pauli_string,wire_map,expected_pauli",
+        [
+            ("I", {"a": 0}, Identity("a")),
+            ("X", {0: 0}, PauliX(0)),
+            ("XI", {1: 0, 0: 1}, PauliX(1)),
+            ("II", {0: 0, 1: 1}, Identity(0)),
+            ("ZYIZ", {"a": 0, "b": 1, "c": 2, "d": 3}, PauliZ("a") @ PauliY("b") @ PauliZ("d")),
+            ("ZYZ", None, PauliZ(0) @ PauliY(1) @ PauliZ(2)),
+            ("ZIYX", {"d": 0, "c": 1, "b": 2, "a": 3}, PauliZ("d") @ PauliY("b") @ PauliX("a")),
+        ],
+    )
+    def test_string_to_pauli_word(self, pauli_string, wire_map, expected_pauli):
+        """Test that valid strings are correctly converted into Pauli words."""
+        obtained_pauli = string_to_pauli_word(pauli_string, wire_map)
+        assert obtained_pauli.compare(expected_pauli)
+
+    @pytest.mark.parametrize(
+        "non_pauli_string,wire_map,error_type,error_message",
+        [
+            (Identity("a"), None, TypeError, "must be string"),
+            ("XAYZ", None, ValueError, "Invalid characters encountered"),
+            ("XYYZ", {0: 0, 1: 1, 2: 2}, ValueError, "must have the same length"),
+        ],
+    )
+    def test_string_to_pauli_word_invalid_input(
+        self, non_pauli_string, wire_map, error_type, error_message
+    ):
+        """Ensure invalid inputs are handled properly when converting strings to Pauli words."""
+        with pytest.raises(error_type, match=error_message):
+            string_to_pauli_word(non_pauli_string, wire_map)
+
+    @pytest.mark.parametrize(
+        "pauli_word,wire_map,expected_matrix",
+        [
+            (PauliX(0), {0: 0}, PauliX.matrix),
+            (Identity(0), {0: 0}, np.eye(2)),
+            (
+                PauliZ(0) @ PauliY(1),
+                {0: 0, 1: 1},
+                np.array([[0, -1j, 0, 0], [1j, 0, 0, 0], [0, 0, 0, 1j], [0, 0, -1j, 0]]),
+            ),
+            (
+                PauliY(1) @ PauliZ(0),
+                {0: 0, 1: 1},
+                np.array([[0, -1j, 0, 0], [1j, 0, 0, 0], [0, 0, 0, 1j], [0, 0, -1j, 0]]),
+            ),
+            (
+                PauliY(1) @ PauliZ(0),
+                {1: 0, 0: 1},
+                np.array([[0, 0, -1j, 0], [0, 0, 0, 1j], [1j, 0, 0, 0], [0, -1j, 0, 0]]),
+            ),
+            (Identity(0), {0: 0, 1: 1}, np.eye(4)),
+            (PauliX(2), None, PauliX.matrix),
+            (
+                PauliX(2),
+                {0: 0, 1: 1, 2: 2},
+                np.array(
+                    [
+                        [0, 1, 0, 0, 0, 0, 0, 0],
+                        [1, 0, 0, 0, 0, 0, 0, 0],
+                        [0, 0, 0, 1, 0, 0, 0, 0],
+                        [0, 0, 1, 0, 0, 0, 0, 0],
+                        [0, 0, 0, 0, 0, 1, 0, 0],
+                        [0, 0, 0, 0, 1, 0, 0, 0],
+                        [0, 0, 0, 0, 0, 0, 0, 1],
+                        [0, 0, 0, 0, 0, 0, 1, 0],
+                    ]
+                ),
+            ),
+            (
+                PauliZ("a") @ PauliX(2),
+                {"a": 0, 1: 1, 2: 2},
+                np.array(
+                    [
+                        [0, 1, 0, 0, 0, 0, 0, 0],
+                        [1, 0, 0, 0, 0, 0, 0, 0],
+                        [0, 0, 0, 1, 0, 0, 0, 0],
+                        [0, 0, 1, 0, 0, 0, 0, 0],
+                        [0, 0, 0, 0, 0, -1, 0, 0],
+                        [0, 0, 0, 0, -1, 0, 0, 0],
+                        [0, 0, 0, 0, 0, 0, 0, -1],
+                        [0, 0, 0, 0, 0, 0, -1, 0],
+                    ]
+                ),
+            ),
+            (
+                PauliX(2) @ PauliZ("a"),
+                {"a": 0, 1: 1, 2: 2},
+                np.array(
+                    [
+                        [0, 1, 0, 0, 0, 0, 0, 0],
+                        [1, 0, 0, 0, 0, 0, 0, 0],
+                        [0, 0, 0, 1, 0, 0, 0, 0],
+                        [0, 0, 1, 0, 0, 0, 0, 0],
+                        [0, 0, 0, 0, 0, -1, 0, 0],
+                        [0, 0, 0, 0, -1, 0, 0, 0],
+                        [0, 0, 0, 0, 0, 0, 0, -1],
+                        [0, 0, 0, 0, 0, 0, -1, 0],
+                    ]
+                ),
+            ),
+        ],
+    )
+    def test_pauli_word_to_matrix(self, pauli_word, wire_map, expected_matrix):
+        """Test that Pauli words are correctly converted into matrices."""
+        obtained_matrix = pauli_word_to_matrix(pauli_word, wire_map)
+        assert np.allclose(obtained_matrix, expected_matrix)
+
+    @pytest.mark.parametrize("non_pauli_word", non_pauli_words)
+    def test_pauli_word_to_matrix_invalid_input(self, non_pauli_word):
+        """Ensure invalid inputs are handled properly when converting Pauli words to matrices."""
+        with pytest.raises(TypeError):
+            pauli_word_to_matrix(non_pauli_word)
+
+    @pytest.mark.parametrize(
+        "pauli_word_1,pauli_word_2,wire_map,commute_status",
+        [
+            (Identity(0), PauliZ(0), {0: 0}, True),
+            (PauliY(0), PauliZ(0), {0: 0}, False),
+            (PauliX(0), PauliX(1), {0: 0, 1: 1}, True),
+            (PauliY("x"), PauliX("y"), None, True),
+            (
+                PauliZ("a") @ PauliY("b") @ PauliZ("d"),
+                PauliX("a") @ PauliZ("c") @ PauliY("d"),
+                {"a": 0, "b": 1, "c": 2, "d": 3},
+                True,
+            ),
+            (
+                PauliX("a") @ PauliY("b") @ PauliZ("d"),
+                PauliX("a") @ PauliZ("c") @ PauliY("d"),
+                {"a": 0, "b": 1, "c": 2, "d": 3},
+                False,
+            ),
+        ],
+    )
+    def test_is_commuting(self, pauli_word_1, pauli_word_2, wire_map, commute_status):
+        """Test that (non)-commuting Pauli words are correctly identified."""
+        do_they_commute = is_commuting(pauli_word_1, pauli_word_2, wire_map=wire_map)
+        assert do_they_commute == commute_status
+
+    @pytest.mark.parametrize(
+        "pauli_word_1,pauli_word_2",
+        [(non_pauli_words[0], PauliX(0) @ PauliY(2)), (PauliX(0) @ PauliY(2), non_pauli_words[0])],
+    )
+    def test_is_commuting_invalid_input(self, pauli_word_1, pauli_word_2):
+        """Ensure invalid inputs are handled properly when determining commutativity."""
+        with pytest.raises(TypeError):
+            is_commuting(pauli_word_1, pauli_word_2)
diff --git a/tests/grouping/test_optimize_measurements.py b/tests/grouping/test_optimize_measurements.py
index 9b24af5c5..1f0184f82 100644
--- a/tests/grouping/test_optimize_measurements.py
+++ b/tests/grouping/test_optimize_measurements.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2020 Xanadu Quantum Technologies Inc.
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
 
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
diff --git a/tests/grouping/test_pauli_group.py b/tests/grouping/test_pauli_group.py
new file mode 100644
index 000000000..22967f278
--- /dev/null
+++ b/tests/grouping/test_pauli_group.py
@@ -0,0 +1,198 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Unit tests for the :mod:`pauli_group`  functions in ``grouping/pauli.py``.
+"""
+import pytest
+from pennylane import Identity, PauliX, PauliY, PauliZ
+
+from pennylane.grouping.pauli import pauli_group, pauli_mult, pauli_mult_with_phase
+
+
+class TestPauliGroup:
+    """Testing for Pauli group construction and manipulation functions."""
+
+    def test_pauli_group_size(self):
+        """Test that the size of the returned Pauli group is correct."""
+        for n_qubits in range(1, 5):
+            pg = list(pauli_group(n_qubits))
+            assert len(pg) == 4 ** n_qubits
+
+    def test_pauli_group_invalid_input(self):
+        """Test that invalid inputs to the Pauli group are handled correctly."""
+        with pytest.raises(TypeError, match="Must specify an integer number"):
+            pauli_group("3")
+
+        with pytest.raises(ValueError, match="Number of qubits must be at least 1"):
+            pauli_group(-1)
+
+    def test_one_qubit_pauli_group_integer_wire_map(self):
+        """Test that the single-qubit Pauli group is constructed correctly with a wire
+        labeled by an integer."""
+
+        expected_pg_1 = [Identity(0), PauliZ(0), PauliX(0), PauliY(0)]
+        pg_1 = list(pauli_group(1))
+        assert all([expected.compare(obtained) for expected, obtained in zip(expected_pg_1, pg_1)])
+
+    def test_one_qubit_pauli_group_valid_float_input(self):
+        """Test that the single-qubit Pauli group is constructed correctly when a float
+        that represents an integer is passed."""
+
+        expected_pg_1 = [Identity(0), PauliZ(0), PauliX(0), PauliY(0)]
+        pg_1 = list(pauli_group(1.0))
+        assert all([expected.compare(obtained) for expected, obtained in zip(expected_pg_1, pg_1)])
+
+    def test_one_qubit_pauli_group_string_wire_map(self):
+        """Test that the single-qubit Pauli group is constructed correctly with a wire
+        labeled by a string."""
+
+        wire_map = {"qubit": 0}
+        expected_pg_1_wires = [
+            Identity("qubit"),
+            PauliZ("qubit"),
+            PauliX("qubit"),
+            PauliY("qubit"),
+        ]
+        pg_1_wires = list(pauli_group(1, wire_map=wire_map))
+        assert all(
+            [
+                expected.compare(obtained)
+                for expected, obtained in zip(expected_pg_1_wires, pg_1_wires)
+            ]
+        )
+
+    def test_two_qubit_pauli_group(self):
+        """Test that the two-qubit Pauli group is constructed correctly."""
+        # With no wire map; ordering is based on construction from binary representation
+        wire_map = {"a": 0, "b": 1}
+
+        expected_pg_2 = [
+            Identity("a"),
+            PauliZ("b"),
+            PauliZ("a"),
+            PauliZ("a") @ PauliZ("b"),
+            PauliX("b"),
+            PauliY("b"),
+            PauliZ("a") @ PauliX("b"),
+            PauliZ("a") @ PauliY("b"),
+            PauliX("a"),
+            PauliX("a") @ PauliZ("b"),
+            PauliY("a"),
+            PauliY("a") @ PauliZ("b"),
+            PauliX("a") @ PauliX("b"),
+            PauliX("a") @ PauliY("b"),
+            PauliY("a") @ PauliX("b"),
+            PauliY("a") @ PauliY("b"),
+        ]
+
+        pg_2 = list(pauli_group(2, wire_map=wire_map))
+        assert all([expected.compare(obtained) for expected, obtained in zip(expected_pg_2, pg_2)])
+
+    @pytest.mark.parametrize(
+        "pauli_word_1,pauli_word_2,wire_map,expected_product",
+        [
+            (PauliX(0), Identity(0), {0: 0}, PauliX(0)),
+            (PauliZ(0), PauliY(0), {0: 0}, PauliX(0)),
+            (PauliZ(0), PauliZ(0), {0: 0}, Identity(0)),
+            (Identity("a"), Identity("b"), None, Identity("a")),
+            (PauliZ("b") @ PauliY("a"), PauliZ("b") @ PauliY("a"), None, Identity("b")),
+            (
+                PauliZ("b") @ PauliY("a"),
+                PauliZ("b") @ PauliY("a"),
+                {"b": 0, "a": 1},
+                Identity("b"),
+            ),
+            (
+                PauliZ(0) @ PauliY(1),
+                PauliX(0) @ PauliZ(1),
+                {0: 0, 1: 1},
+                PauliY(0) @ PauliX(1),
+            ),
+            (PauliZ(0) @ PauliY(1), PauliX(1) @ PauliY(0), {0: 0, 1: 1}, PauliX(0) @ PauliZ(1)),
+            (
+                PauliZ(0) @ PauliY(3) @ PauliZ(1),
+                PauliX(1) @ PauliX(2) @ PauliY(0),
+                None,
+                PauliX(0) @ PauliY(1) @ PauliX(2) @ PauliY(3),
+            ),
+            (
+                PauliX(0) @ PauliX(2),
+                PauliX(0) @ PauliZ(2),
+                None,
+                PauliY(2),
+            ),
+            (PauliZ("a"), PauliX("b"), {"a": 0, "b": 1}, PauliZ("a") @ PauliX("b")),
+            (
+                PauliZ("a"),
+                PauliX("e"),
+                {"a": 0, "b": 1, "c": 2, "d": 3, "e": 4},
+                PauliZ("a") @ PauliX("e"),
+            ),
+            (PauliZ("a"), PauliY("e"), None, PauliZ("a") @ PauliY("e")),
+            (
+                PauliZ(0) @ PauliZ(2) @ PauliZ(4),
+                PauliZ(0) @ PauliY(1) @ PauliX(3),
+                None,
+                PauliZ(2) @ PauliY(1) @ PauliZ(4) @ PauliX(3),
+            ),
+            (
+                PauliZ(0) @ PauliZ(4) @ PauliZ(2),
+                PauliZ(0) @ PauliY(1) @ PauliX(3),
+                {0: 0, 2: 1, 4: 2, 1: 3, 3: 4},
+                PauliZ(2) @ PauliY(1) @ PauliZ(4) @ PauliX(3),
+            ),
+            (
+                PauliZ(0) @ PauliY(3) @ PauliZ(1),
+                PauliX(1) @ PauliX(2) @ PauliY(0),
+                None,
+                PauliX(0) @ PauliY(3) @ PauliY(1) @ PauliX(2),
+            ),
+        ],
+    )
+    def test_pauli_mult(self, pauli_word_1, pauli_word_2, wire_map, expected_product):
+        """Test that Pauli words are multiplied together correctly."""
+        obtained_product = pauli_mult(pauli_word_1, pauli_word_2, wire_map=wire_map)
+        assert obtained_product.compare(expected_product)
+
+    @pytest.mark.parametrize(
+        "pauli_word_1,pauli_word_2,wire_map,expected_phase",
+        [
+            (PauliX(0), Identity(0), {0: 0}, 1),
+            (PauliZ(0), PauliY(0), {0: 0}, -1j),
+            (PauliZ(0), PauliZ(0), {0: 0}, 1),
+            (Identity("a"), Identity("b"), None, 1),
+            (PauliZ("b") @ PauliY("a"), PauliZ("b") @ PauliY("a"), None, 1),
+            (PauliZ(0), PauliY("b"), None, 1),
+            (
+                PauliZ("a") @ PauliY("b"),
+                PauliX("a") @ PauliZ("b"),
+                {"a": 0, "b": 1},
+                -1,
+            ),
+            (PauliX(0) @ PauliX(2), PauliX(0) @ PauliZ(2), None, -1j),
+            (
+                PauliX(0) @ PauliY(1) @ PauliZ(2),
+                PauliY(0) @ PauliY(1),
+                {0: 0, 1: 1, 2: 2},
+                1j,
+            ),
+            (PauliZ(0) @ PauliY(1), PauliX(1) @ PauliY(0), {0: 0, 1: 1}, -1),
+            (PauliZ(0) @ PauliY(1), PauliX(1) @ PauliY(0), {1: 0, 0: 1}, -1),
+            (PauliZ(0) @ PauliY(3) @ PauliZ(1), PauliX(1) @ PauliX(2) @ PauliY(0), None, 1),
+        ],
+    )
+    def test_pauli_mult_with_phase(self, pauli_word_1, pauli_word_2, wire_map, expected_phase):
+        """Test that multiplication including phases works as expected."""
+        _, obtained_phase = pauli_mult_with_phase(pauli_word_1, pauli_word_2, wire_map=wire_map)
+        assert obtained_phase == expected_phase