Repository: PennyLaneAI/pennylane
Pull request #1974 (2021-12-02 16:12:16)

Add symmetry generation functions for qubit tapering
**Context:**
This PR adds functions to calculate the possible symmetries for a Hamiltonian, that can be used to transform the Hamiltonian such that it acts trivially on a subset of qubits, which can be then tapered or removed from the Hamiltonian.

**Description of the Change:**
1. Contains functions for computing row echelon form `_row_echelon_form` and kernel `_kernel` of a binary matrix in a binary field `Z2`. 
2. Contains function `generate_taus` for computing symmetry generators `taus` using the kernel of the binary matrix.
3. Contains function `generate_paulis` for computing PauliX operators that will be used to calculate the Clifford unitaries `U`. 

**Benefits:**

**Possible Drawbacks:**

**Related GitHub Issues:**


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.
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1974?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#1974](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1974?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (c374685) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/d7fe61be52c2e5f28f58da68b65cea0d3c1149c2?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (d7fe61b) will **decrease** coverage by `0.44%`.
> The diff coverage is `0.00%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1974/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1974?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@            Coverage Diff             @@
##           master    #1974      +/-   ##
==========================================
- Coverage   98.81%   98.37%   -0.45%     
==========================================
  Files         225      226       +1     
  Lines       17312    17390      +78     
==========================================
  Hits        17107    17107              
- Misses        205      283      +78     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1974?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/hf/tapering.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1974/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2hmL3RhcGVyaW5nLnB5) | `0.00% <0.00%> (ø)` | |

------

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


Patch:
diff --git a/pennylane/hf/__init__.py b/pennylane/hf/__init__.py
index bb0353074..de79245cc 100644
--- a/pennylane/hf/__init__.py
+++ b/pennylane/hf/__init__.py
@@ -44,3 +44,8 @@ from .matrices import (
     molecular_density_matrix,
 )
 from .molecule import Molecule
+from .tapering import (
+    get_generators,
+    generate_paulis,
+    generate_symmetries,
+)
diff --git a/pennylane/hf/tapering.py b/pennylane/hf/tapering.py
new file mode 100644
index 000000000..c978ea60d
--- /dev/null
+++ b/pennylane/hf/tapering.py
@@ -0,0 +1,303 @@
+# 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.
+"""
+This module contains the functions needed for tapering qubits using symmetries.
+"""
+import numpy as np
+import pennylane as qml
+
+
+def _binary_matrix(terms, num_qubits):
+    r"""Get a binary matrix representation of the Hamiltonian where each row corresponds to a
+    Pauli term, which is represented by a concatenation of Z and X vectors.
+
+    Args:
+        terms (Iterable[pennylane.Observable]): operators defining the Hamiltonian
+        num_qubits (int): number of wires required to define the Hamiltonian
+    Returns:
+        array[int]: binary matrix representation of the Hamiltonian of shape
+        :math:`len(terms) \times 2*num_qubits`
+
+    **Example**
+
+    .. code-block:: python
+
+        >>> terms = [qml.PauliZ(wires=[0]) @ qml.PauliX(wires=[1]),
+        ...          qml.PauliZ(wires=[0]) @ qml.PauliY(wires=[2]),
+        ...          qml.PauliX(wires=[0]) @ qml.PauliY(wires=[3])]
+        >>> _binary_matrix(terms, 4)
+         array([[1, 0, 0, 0, 0, 1, 0, 0],
+                [1, 0, 1, 0, 0, 0, 1, 0],
+                [0, 0, 0, 1, 1, 0, 0, 1]]))
+
+    """
+
+    binary_matrix = np.zeros((len(terms), 2 * num_qubits), dtype=int)
+    for idx, term in enumerate(terms):
+        ops, wires = term.name, term.wires
+        if len(wires) == 1:
+            ops = [ops]
+        for op, wire in zip(ops, wires):
+            if op in ["PauliX", "PauliY"]:
+                binary_matrix[idx][wire + num_qubits] = 1
+            if op in ["PauliZ", "PauliY"]:
+                binary_matrix[idx][wire] = 1
+
+    return binary_matrix
+
+
+def _reduced_row_echelon(binary_matrix):
+    r"""Returns the reduced row echelon form (RREF) of a matrix in a binary finite field :math:`\mathbb{Z}_2`.
+
+    Args:
+        binary_matrix (array[int]): binary matrix representation of the Hamiltonian
+    Returns:
+        array[int]: reduced row-echelon form of the given `binary_matrix`
+
+    **Example**
+
+    .. code-block:: python
+
+        >>> binary_matrix = np.array([[1, 0, 0, 0, 0, 1, 0, 0],
+        ...                           [1, 0, 1, 0, 0, 0, 1, 0],
+        ...                           [0, 0, 0, 1, 1, 0, 0, 1]])
+        >>> _reduced_row_echelon(binary_matrix)
+         array([[1, 0, 0, 0, 0, 1, 0, 0],
+                [0, 0, 1, 0, 0, 1, 1, 0],
+                [0, 0, 0, 1, 1, 0, 0, 1]])
+
+    """
+
+    rref_mat = binary_matrix.copy()
+    shape = rref_mat.shape
+    icol = 0
+
+    for irow in range(shape[0]):
+
+        while icol < shape[1] and not rref_mat[irow][icol]:
+
+            # get the nonzero indices in the remainder of column icol
+            non_zero_idx = rref_mat[irow:, icol].nonzero()[0]
+
+            if len(non_zero_idx) == 0:  # if remainder of column icol is all zero
+                icol += 1
+            else:
+                # find value and index of largest element in remainder of column icol
+                krow = irow + non_zero_idx[0]
+
+                # swap rows krow and irow
+                rref_mat[irow, icol:], rref_mat[krow, icol:] = (
+                    rref_mat[krow, icol:].copy(),
+                    rref_mat[irow, icol:].copy(),
+                )
+        if icol < shape[1] and rref_mat[irow][icol]:
+
+            # store remainder right hand side columns of the pivot row irow
+            rpvt_cols = rref_mat[irow, icol:].copy()
+
+            # get the column icol and set its irow element to 0 to avoid XORing pivot row with itself
+            currcol = rref_mat[:, icol].copy()
+            currcol[irow] = 0
+
+            # XOR the right hand side of the pivot row irow with all of the other rows
+            rref_mat[:, icol:] ^= np.outer(currcol, rpvt_cols)
+            icol += 1
+
+    return rref_mat.astype(int)
+
+
+def _kernel(binary_matrix):
+    r"""Computes the kernel of a binary matrix on the binary finite field :math:`\mathbb{Z}_2`.
+
+    Args:
+        binary_matrix (array[int]): binary matrix representation of the Hamiltonian
+    Returns:
+        array[int]: nullspace of the `binary_matrix` where each row corresponds to a
+        basis vector in the nullspace
+
+    **Example**
+
+    .. code-block:: python
+
+        >>> binary_matrix = np.array([[1, 0, 0, 0, 0, 1, 0, 0],
+        ...                          [0, 0, 1, 1, 1, 1, 1, 1],
+        ...                          [0, 0, 0, 1, 1, 0, 0, 1]])
+        >>> _kernel(binary_matrix)
+         array([[0, 1, 0, 0, 0, 0, 0, 0],
+                [0, 0, 1, 1, 1, 0, 0, 0],
+                [1, 0, 1, 0, 0, 1, 0, 0],
+                [0, 0, 1, 0, 0, 0, 1, 0],
+                [0, 0, 1, 1, 0, 0, 0, 1]])
+
+    """
+
+    # Get the columns with and without pivots
+    pivots = (binary_matrix.T != 0).argmax(axis=0)
+    nonpivots = np.setdiff1d(range(len(binary_matrix[0])), pivots)
+
+    # Initialize the nullspace
+    null_vector = np.zeros((binary_matrix.shape[1], len(nonpivots)), dtype=int)
+    null_vector[nonpivots, np.arange(len(nonpivots))] = 1
+
+    # Fill up the nullspace vectors from the binary matrix
+    null_vector_indices = np.ix_(pivots, np.arange(len(nonpivots)))
+    binary_vector_indices = np.ix_(np.arange(len(pivots)), nonpivots)
+    null_vector[null_vector_indices] = -binary_matrix[binary_vector_indices] % 2
+
+    nullspace = null_vector.T
+    return nullspace
+
+
+def get_generators(nullspace, num_qubits):
+    r"""Compute the generators :math:`\{\tau_1, \ldots, \tau_k\}` from the nullspace of
+    the binary matrix form of a Hamiltonian over the binary field :math:`\mathbb{Z}_2`.
+    These correspond to the generator set of the :math:`\mathbb{Z}_2`-symmetries present
+    in the Hamiltonian as given in `arXiv:1910.14644 <https://arxiv.org/abs/1910.14644>`_.
+
+    Args:
+        nullspace (array[int]): kernel of the binary matrix corresponding to the Hamiltonian
+        num_qubits (int): number of wires required to define the Hamiltonian
+
+    Returns:
+        list[pennylane.Hamiltonian]: list of generators of symmetries, taus, for the Hamiltonian
+
+    **Example**
+
+    .. code-block:: python
+
+        >>> kernel = np.array([[0, 1, 0, 0, 0, 0, 0, 0],
+        ...                    [0, 0, 1, 1, 1, 0, 0, 0],
+        ...                    [1, 0, 1, 0, 0, 1, 0, 0],
+        ...                    [0, 0, 1, 0, 0, 0, 1, 0],
+        ...                    [0, 0, 1, 1, 0, 0, 0, 1]])
+        >>> generate_taus(kernel, 4)
+         [(1.0) [X1], (1.0) [Z0 X2 X3], (1.0) [X0 Z1 X2], (1.0) [Y2], (1.0) [X2 Y3]]
+
+    """
+
+    generators = []
+    pauli_map = {"00": qml.Identity, "10": qml.PauliX, "11": qml.PauliY, "01": qml.PauliZ}
+
+    for null_vector in nullspace:
+        tau = qml.Identity(0)
+        for idx, op in enumerate(zip(null_vector[:num_qubits], null_vector[num_qubits:])):
+            x, z = op
+            tau @= pauli_map[f"{x}{z}"](idx)
+
+        ham = qml.Hamiltonian([1.0], [tau], simplify=True)
+        generators.append(ham)
+
+    return generators
+
+
+def generate_paulis(generators, num_qubits):
+    r"""Generate the single qubit Pauli-X operators :math:`\sigma^{x}_{i}` for each symmetry :math:`\tau_j`,
+    such that it anti-commutes with :math:`\tau_j` and commutes with all others symmetries :math:`\tau_{k\neq j}`.
+    These are required to obtain the Clifford operators :math:`U` for the Hamiltonian :math:`H`.
+
+    Args:
+        generators (list[pennylane.Hamiltonian]): list of generators of symmetries, taus, for the Hamiltonian
+        num_qubits (int): number of wires required to define the Hamiltonian
+    Return:
+        list[pennylane.Observable]: list of single-qubit Pauli-X operators which will be used to build the
+        Clifford operators :math:`U`.
+
+    **Example**
+
+    .. code-block:: python
+
+        >>> generators = [qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(1)]),
+        ...               qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(2)]),
+        ...               qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(3)])]
+        >>> generate_paulis(generators, qubits)
+         [PauliX(wires=[1]), PauliX(wires=[2]), PauliX(wires=[3])]
+
+    """
+
+    ops_generator = [g.ops[0] if isinstance(g.ops, list) else g.ops for g in generators]
+    bmat = _binary_matrix(ops_generator, num_qubits)
+
+    pauli_x_ops = []
+    for row in range(bmat.shape[0]):
+        bmatrow = bmat[row]
+        bmatrest = np.delete(bmat, row, axis=0)
+        for col in range(bmat.shape[1] // 2):
+            # Anti-commutes with the (row) and commutes with all other symmetries.
+            if bmatrow[col] and np.array_equal(
+                bmatrest[:, col], np.zeros(bmat.shape[0] - 1, dtype=int)
+            ):
+                pauli_x_ops.append(qml.PauliX(col))
+                break
+
+    return pauli_x_ops
+
+
+def generate_symmetries(qubit_op, num_qubits):
+    r"""Compute the generator set of the symmetries :math:`\mathbf{\tau}` and the corresponding single-qubit
+    set of the Pauli-X operators :math:`\mathbf{\sigma^x}` that are used to build the Clifford operators
+    :math:`U`, according to the following relation:
+
+    .. math::
+
+        U_i = \frac{1}{\sqrt{2}}(\tau_i+\sigma^{x}_{q}).
+
+    Here, :math:`\sigma^{x}_{q}` is the Pauli-X operator acting on qubit :math:`q`. These :math:`U_i` can be
+    used to transform the Hamiltonian :math:`H` in such a way that it acts trivially or at most with one
+    Pauli-gate on a subset of qubits, which allows us to taper off those qubits from the simulation
+    using :func:`~.transform_hamiltonian`.
+
+    Args:
+        qubit_op (pennylane.Hamiltonian): Hamiltonian for which symmetries are to be generated to perform tapering
+        num_qubits (int): number of wires required to define the Hamiltonian
+
+    Returns:
+        tuple (list[pennylane.Hamiltonian], list[pennylane.operation.Observable]):
+
+            * list[pennylane.Hamiltonian]: list of generators of symmetries, :math:`\mathbf{\tau}`,
+              for the Hamiltonian.
+            * list[pennylane.operation.Operation]: list of single-qubit Pauli-X operators which will be used
+              to build the Clifford operators :math:`U`.
+
+    .. code-block:: python
+
+        >>> symbols = ['H', 'H']
+        >>> coordinates = np.array([0., 0., -0.66140414, 0., 0., 0.66140414])
+        >>> mol = qml.hf.Molecule(symbols, coordinates)
+        >>> H, qubits = qml.hf.generate_hamiltonian(mol)(), 4
+        >>> generators, pauli_x = generate_symmetries(H, qubits)
+        >>> generators, pauli_x
+         ([(1.0) [Z0 Z1], (1.0) [Z0 Z2], (1.0) [Z0 Z3]],
+          [PauliX(wires=[1]), PauliX(wires=[2]), PauliX(wires=[3])])
+
+    """
+
+    # Generate binary matrix for qubit_op
+    binary_matrix = _binary_matrix(qubit_op.ops, num_qubits)
+
+    # Get reduced row echelon form of binary matrix
+    rref_binary_matrix = _reduced_row_echelon(binary_matrix)
+    rref_binary_matrix_red = rref_binary_matrix[
+        ~np.all(rref_binary_matrix == 0, axis=1)
+    ]  # remove all-zero rows
+
+    # Get kernel (i.e., nullspace) for trimmed binary matrix using gaussian elimination
+    nullspace = _kernel(rref_binary_matrix_red)
+
+    # Get generators tau from the calculated nullspace
+    generators = get_generators(nullspace, num_qubits)
+
+    # Get unitaries from the calculated nullspace
+    pauli_x_ops = generate_paulis(generators, num_qubits)
+
+    return generators, pauli_x_ops


Test patch:
diff --git a/tests/hf/test_tapering.py b/tests/hf/test_tapering.py
new file mode 100644
index 000000000..0d77b5ea4
--- /dev/null
+++ b/tests/hf/test_tapering.py
@@ -0,0 +1,357 @@
+# 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 functions needed for qubit tapering.
+"""
+
+import pytest
+import pennylane as qml
+from pennylane import numpy as np
+from pennylane.hf.tapering import (
+    _binary_matrix,
+    _reduced_row_echelon,
+    _kernel,
+    get_generators,
+    generate_paulis,
+    generate_symmetries,
+)
+
+
+@pytest.mark.parametrize(
+    ("terms", "num_qubits", "result"),
+    [
+        (
+            [
+                qml.Identity(wires=[0]),
+                qml.PauliZ(wires=[0]),
+                qml.PauliZ(wires=[1]),
+                qml.PauliZ(wires=[2]),
+                qml.PauliZ(wires=[3]),
+                qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[1]),
+                qml.PauliY(wires=[0])
+                @ qml.PauliX(wires=[1])
+                @ qml.PauliX(wires=[2])
+                @ qml.PauliY(wires=[3]),
+                qml.PauliY(wires=[0])
+                @ qml.PauliY(wires=[1])
+                @ qml.PauliX(wires=[2])
+                @ qml.PauliX(wires=[3]),
+                qml.PauliX(wires=[0])
+                @ qml.PauliX(wires=[1])
+                @ qml.PauliY(wires=[2])
+                @ qml.PauliY(wires=[3]),
+                qml.PauliX(wires=[0])
+                @ qml.PauliY(wires=[1])
+                @ qml.PauliY(wires=[2])
+                @ qml.PauliX(wires=[3]),
+                qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[2]),
+                qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[3]),
+                qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[2]),
+                qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[3]),
+                qml.PauliZ(wires=[2]) @ qml.PauliZ(wires=[3]),
+            ],
+            4,
+            np.array(
+                [
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [1, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 1, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 1, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 1, 0, 0, 0, 0],
+                    [1, 1, 0, 0, 0, 0, 0, 0],
+                    [1, 0, 0, 1, 1, 1, 1, 1],
+                    [1, 1, 0, 0, 1, 1, 1, 1],
+                    [0, 0, 1, 1, 1, 1, 1, 1],
+                    [0, 1, 1, 0, 1, 1, 1, 1],
+                    [1, 0, 1, 0, 0, 0, 0, 0],
+                    [1, 0, 0, 1, 0, 0, 0, 0],
+                    [0, 1, 1, 0, 0, 0, 0, 0],
+                    [0, 1, 0, 1, 0, 0, 0, 0],
+                    [0, 0, 1, 1, 0, 0, 0, 0],
+                ]
+            ),
+        ),
+    ],
+)
+def test_binary_matrix(terms, num_qubits, result):
+    r"""Test that _binary_matrix returns the correct result."""
+    binary_matrix = _binary_matrix(terms, num_qubits)
+    assert (binary_matrix == result).all()
+
+
+@pytest.mark.parametrize(
+    ("binary_matrix", "result"),
+    [
+        (
+            np.array(
+                [
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [1, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 1, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 1, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 1, 0, 0, 0, 0],
+                    [1, 1, 0, 0, 0, 0, 0, 0],
+                    [1, 0, 0, 1, 1, 1, 1, 1],
+                    [1, 1, 0, 0, 1, 1, 1, 1],
+                    [0, 0, 1, 1, 1, 1, 1, 1],
+                    [0, 1, 1, 0, 1, 1, 1, 1],
+                    [1, 0, 1, 0, 0, 0, 0, 0],
+                    [1, 0, 0, 1, 0, 0, 0, 0],
+                    [0, 1, 1, 0, 0, 0, 0, 0],
+                    [0, 1, 0, 1, 0, 0, 0, 0],
+                    [0, 0, 1, 1, 0, 0, 0, 0],
+                ]
+            ),
+            np.array(
+                [
+                    [1, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 1, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 1, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 1, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 1, 1, 1, 1],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                ]
+            ),
+        ),
+    ],
+)
+def test_reduced_row_echelon(binary_matrix, result):
+    r"""Test that _reduced_row_echelon returns the correct result."""
+
+    # build row echelon form of the matrix
+    shape = binary_matrix.shape
+    for irow in range(shape[0]):
+        pivot_index = 0
+        if np.count_nonzero(binary_matrix[irow, :]):
+            pivot_index = np.nonzero(binary_matrix[irow, :])[0][0]
+
+        for jrow in range(shape[0]):
+            if jrow != irow and binary_matrix[jrow, pivot_index]:
+                binary_matrix[jrow, :] = (binary_matrix[jrow, :] + binary_matrix[irow, :]) % 2
+
+    indices = [
+        irow
+        for irow in range(shape[0] - 1)
+        if np.array_equal(binary_matrix[irow, :], np.zeros(shape[1]))
+    ]
+
+    temp_row_echelon_matrix = binary_matrix.copy()
+    for row in indices[::-1]:
+        temp_row_echelon_matrix = np.delete(temp_row_echelon_matrix, row, axis=0)
+
+    row_echelon_matrix = np.zeros(shape, dtype=int)
+    row_echelon_matrix[: shape[0] - len(indices), :] = temp_row_echelon_matrix
+
+    # build reduced row echelon form of the matrix from row echelon form
+    for idx in range(len(row_echelon_matrix))[:0:-1]:
+        nonzeros = np.nonzero(row_echelon_matrix[idx])[0]
+        if len(nonzeros) > 0:
+            redrow = (row_echelon_matrix[idx, :] % 2).reshape(1, -1)
+            coeffs = (
+                (-row_echelon_matrix[:idx, nonzeros[0]] / row_echelon_matrix[idx, nonzeros[0]]) % 2
+            ).reshape(1, -1)
+            row_echelon_matrix[:idx, :] = (
+                row_echelon_matrix[:idx, :] + (coeffs.T * redrow) % 2
+            ) % 2
+
+    # get reduced row echelon form from the _reduced_row_echelon function
+    rref_bin_mat = _reduced_row_echelon(binary_matrix)
+
+    assert (rref_bin_mat == row_echelon_matrix).all()
+    assert (rref_bin_mat == result).all()
+
+
+@pytest.mark.parametrize(
+    ("binary_matrix", "result"),
+    [
+        (
+            np.array(
+                [
+                    [1, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 1, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 1, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 1, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 1, 1, 1, 1],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                    [0, 0, 0, 0, 0, 0, 0, 0],
+                ]
+            ),
+            np.array(
+                [[0, 0, 0, 0, 1, 1, 0, 0], [0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 1, 0, 0, 1]]
+            ),
+        ),
+    ],
+)
+def test_kernel(binary_matrix, result):
+    r"""Test that _kernel returns the correct result."""
+
+    # get the kernel from the gaussian elimination.
+    pivots = (binary_matrix.T != 0).argmax(axis=0)
+    nonpivots = np.setdiff1d(range(len(binary_matrix[0])), pivots)
+
+    kernel = []
+    for col in nonpivots:
+        col_vector = binary_matrix[:, col]
+        null_vector = np.zeros((binary_matrix.shape[1]), dtype=int)
+        null_vector[col] = 1
+        for i in pivots:
+            first_entry = np.where(binary_matrix[:, i] == 1)[0][0]
+            if col_vector[first_entry] == 1:
+                null_vector[i] = 1
+        kernel.append(null_vector.tolist())
+
+    # get the nullspace from the _kernel function.
+    nullspace = _kernel(binary_matrix)
+
+    for nullvec in kernel:
+        assert nullvec in nullspace.tolist()
+
+    assert (nullspace == result).all()
+
+
+@pytest.mark.parametrize(
+    ("nullspace", "num_qubits", "result"),
+    [
+        (
+            np.array(
+                [[0, 0, 0, 0, 1, 1, 0, 0], [0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 1, 0, 0, 1]]
+            ),
+            4,
+            [  # Correct generators as given by Bravyi et al. (arXiv:1701.08213).
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(1)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(2)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(3)]),
+            ],
+        ),
+    ],
+)
+def test_get_generators(nullspace, num_qubits, result):
+    r"""Test that get_generators returns the correct result."""
+
+    generators = get_generators(nullspace, num_qubits)
+    for g1, g2 in zip(generators, result):
+        assert g1.compare(g2)
+
+
+@pytest.mark.parametrize(
+    ("generators", "num_qubits", "result"),
+    [
+        (
+            [
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(1)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(2)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(3)]),
+            ],
+            4,
+            [qml.PauliX(wires=[1]), qml.PauliX(wires=[2]), qml.PauliX(wires=[3])],
+        ),
+    ],
+)
+def test_generate_paulis(generators, num_qubits, result):
+    r"""Test that generate_paulis returns the correct result."""
+    pauli_ops = generate_paulis(generators, num_qubits)
+    for p1, p2 in zip(pauli_ops, result):
+        assert p1.compare(p2)
+
+
+@pytest.mark.parametrize(
+    ("hamiltonian", "num_qubits", "res_generators", "res_pauli_ops"),
+    [
+        (
+            qml.Hamiltonian(
+                np.array(
+                    [
+                        -0.09886397,
+                        0.17119775,
+                        0.17119775,
+                        -0.22278593,
+                        -0.22278593,
+                        0.16862219,
+                        0.0453222,
+                        -0.0453222,
+                        -0.0453222,
+                        0.0453222,
+                        0.12054482,
+                        0.16586702,
+                        0.16586702,
+                        0.12054482,
+                        0.17434844,
+                    ]
+                ),
+                [
+                    qml.Identity(wires=[0]),
+                    qml.PauliZ(wires=[0]),
+                    qml.PauliZ(wires=[1]),
+                    qml.PauliZ(wires=[2]),
+                    qml.PauliZ(wires=[3]),
+                    qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[1]),
+                    qml.PauliY(wires=[0])
+                    @ qml.PauliX(wires=[1])
+                    @ qml.PauliX(wires=[2])
+                    @ qml.PauliY(wires=[3]),
+                    qml.PauliY(wires=[0])
+                    @ qml.PauliY(wires=[1])
+                    @ qml.PauliX(wires=[2])
+                    @ qml.PauliX(wires=[3]),
+                    qml.PauliX(wires=[0])
+                    @ qml.PauliX(wires=[1])
+                    @ qml.PauliY(wires=[2])
+                    @ qml.PauliY(wires=[3]),
+                    qml.PauliX(wires=[0])
+                    @ qml.PauliY(wires=[1])
+                    @ qml.PauliY(wires=[2])
+                    @ qml.PauliX(wires=[3]),
+                    qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[2]),
+                    qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[3]),
+                    qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[2]),
+                    qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[3]),
+                    qml.PauliZ(wires=[2]) @ qml.PauliZ(wires=[3]),
+                ],
+            ),
+            4,
+            [
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(1)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(2)]),
+                qml.Hamiltonian([1.0], [qml.PauliZ(0) @ qml.PauliZ(3)]),
+            ],
+            [qml.PauliX(wires=[1]), qml.PauliX(wires=[2]), qml.PauliX(wires=[3])],
+        ),
+    ],
+)
+def test_generate_symmetries(hamiltonian, num_qubits, res_generators, res_pauli_ops):
+    r"""Test that generate_symmetries returns the correct result."""
+
+    generators, pauli_ops = generate_symmetries(hamiltonian, num_qubits)
+    for g1, g2 in zip(generators, res_generators):
+        assert g1.compare(g2)
+
+    for p1, p2 in zip(pauli_ops, res_pauli_ops):
+        assert p1.compare(p2)