Repository: PennyLaneAI/pennylane
Pull request #1966 (2021-12-01 22:47:29)

Add Hamiltonian transform functions for qubit tapering
**Context:**
This PR adds functions to transform a Hamiltonian, based on molecular symmetries, such that it acts trivially on a subset of qubits. These qubits can then be removed from the Hamiltonian.

**Description of the Change:**
The function transform_hamiltonian is added.

**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/1966?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#1966](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1966?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (9a717aa) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/eb3a7b38eaa6724c5841614668e829a77607774f?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (eb3a7b3) will **decrease** coverage by `0.31%`.
> The diff coverage is `0.00%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1966/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/1966?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@            Coverage Diff             @@
##           master    #1966      +/-   ##
==========================================
- Coverage   98.81%   98.50%   -0.32%     
==========================================
  Files         225      226       +1     
  Lines       17299    17354      +55     
==========================================
  Hits        17094    17094              
- Misses        205      260      +55     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1966?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/1966/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/1966?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/1966?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [eb3a7b3...9a717aa](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1966?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).

josh146: @soranjh I have a quick minor question: would this functionality be better suited for `pennylane.transforms` rather than `pennylane.hf`? The functionality seems somewhat separate from HF logic.
soranjh: > @soranjh I have a quick minor question: would this functionality be better suited for `pennylane.transforms` rather than `pennylane.hf`? The functionality seems somewhat separate from HF logic.

Hi @josh146. The tapering scheme we are using here is conceptually based on molecular symmetries. That is the main reason we initially put it in `qml.hf`. The `tapering.py` module is self-contained and can be relocated anytime.

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index f0cef4ccb..0b61df0c6 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -4,6 +4,30 @@
 
 <h3>New features since last release</h3>
 
+* Functions for tapering qubits based on molecular symmetries is added.
+  [(#1966)](https://github.com/PennyLaneAI/pennylane/pull/1966)
+* [(#1974)](https://github.com/PennyLaneAI/pennylane/pull/1974)
+
+  With this functionality, a molecular Hamiltonian can be transformed to a new Hamiltonian that acts
+  on a reduced number of qubits.
+
+  ```python
+  symbols = ["H", "H"]
+  geometry = np.array([[0.0, 0.0, -0.69440367], [0.0, 0.0, 0.69440367]])
+  mol = qml.hf.Molecule(symbols, geometry)
+  H = qml.hf.generate_hamiltonian(mol)(geometry)
+  generators, paulix_ops = qml.hf.generate_symmetries(H, len(H.wires))
+  paulix_sector = [1, -1, -1]
+  H_tapered = qml.hf.transform_hamiltonian(H, generators, paulix_ops, paulix_sector)
+  ```
+
+  ```pycon
+  >>> print(H_tapered)
+    ((-0.321034397355719+0j)) [I0]
+  + ((0.1809270275619743+0j)) [X0]
+  + ((0.7959678503870796+0j)) [Z0]
+  ```
+
 * Added the adjoint method for the metric tensor.
   [(#1992)](https://github.com/PennyLaneAI/pennylane/pull/1992)
 
diff --git a/pennylane/hf/__init__.py b/pennylane/hf/__init__.py
index de79245cc..0fdda9ec6 100644
--- a/pennylane/hf/__init__.py
+++ b/pennylane/hf/__init__.py
@@ -45,7 +45,9 @@ from .matrices import (
 )
 from .molecule import Molecule
 from .tapering import (
-    get_generators,
+    clifford,
     generate_paulis,
     generate_symmetries,
+    get_generators,
+    transform_hamiltonian,
 )
diff --git a/pennylane/hf/tapering.py b/pennylane/hf/tapering.py
index c978ea60d..11af063ff 100644
--- a/pennylane/hf/tapering.py
+++ b/pennylane/hf/tapering.py
@@ -14,8 +14,12 @@
 """
 This module contains the functions needed for tapering qubits using symmetries.
 """
-import numpy as np
+# pylint: disable=unnecessary-lambda
+import functools
+
+import autograd.numpy as anp
 import pennylane as qml
+from pennylane import numpy as np
 
 
 def _binary_matrix(terms, num_qubits):
@@ -23,7 +27,7 @@ def _binary_matrix(terms, num_qubits):
     Pauli term, which is represented by a concatenation of Z and X vectors.
 
     Args:
-        terms (Iterable[pennylane.Observable]): operators defining the Hamiltonian
+        terms (Iterable[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
@@ -31,18 +35,14 @@ def _binary_matrix(terms, 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]]))
-
+    >>> 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
@@ -67,18 +67,14 @@ def _reduced_row_echelon(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]])
-
+    >>> 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
@@ -128,20 +124,16 @@ def _kernel(binary_matrix):
 
     **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]])
-
+    >>> 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)
@@ -170,22 +162,18 @@ def get_generators(nullspace, num_qubits):
         num_qubits (int): number of wires required to define the Hamiltonian
 
     Returns:
-        list[pennylane.Hamiltonian]: list of generators of symmetries, taus, for the Hamiltonian
+        list[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]]
-
+    >>> 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]])
+    >>> get_generators(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}
 
@@ -207,28 +195,24 @@ def generate_paulis(generators, num_qubits):
     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
+        generators (list[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
+        list[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])]
-
+    >>> 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, 4)
+    [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 = []
+    paulix_ops = []
     for row in range(bmat.shape[0]):
         bmatrow = bmat[row]
         bmatrest = np.delete(bmat, row, axis=0)
@@ -237,10 +221,10 @@ def generate_paulis(generators, num_qubits):
             if bmatrow[col] and np.array_equal(
                 bmatrest[:, col], np.zeros(bmat.shape[0] - 1, dtype=int)
             ):
-                pauli_x_ops.append(qml.PauliX(col))
+                paulix_ops.append(qml.PauliX(col))
                 break
 
-    return pauli_x_ops
+    return paulix_ops
 
 
 def generate_symmetries(qubit_op, num_qubits):
@@ -258,30 +242,28 @@ def generate_symmetries(qubit_op, num_qubits):
     using :func:`~.transform_hamiltonian`.
 
     Args:
-        qubit_op (pennylane.Hamiltonian): Hamiltonian for which symmetries are to be generated to perform tapering
+        qubit_op (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]):
+        tuple (list[Hamiltonian], list[Observable]):
 
-            * list[pennylane.Hamiltonian]: list of generators of symmetries, :math:`\mathbf{\tau}`,
+            * list[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
+            * list[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])])
+    **Example**
 
+    >>> symbols = ['H', 'H']
+    >>> geometry = np.array([[0., 0., -0.66140414], [0., 0., 0.66140414]])
+    >>> mol = qml.hf.Molecule(symbols, geometry)
+    >>> H, qubits = qml.hf.generate_hamiltonian(mol)(geometry), 4
+    >>> generators, paulix_ops = qml.hf.tapering.generate_symmetries(H, qubits)
+    >>> generators, paulix_ops
+    ([(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)
 
@@ -298,6 +280,187 @@ def generate_symmetries(qubit_op, num_qubits):
     generators = get_generators(nullspace, num_qubits)
 
     # Get unitaries from the calculated nullspace
-    pauli_x_ops = generate_paulis(generators, num_qubits)
+    paulix_ops = generate_paulis(generators, num_qubits)
+
+    return generators, paulix_ops
+
+
+def _observable_mult(obs_a, obs_b):
+    r"""Multiply two PennyLane observables together.
+
+    Each observable should be a linear combination of Pauli words, e.g.,
+    :math:`\sum_{k=0}^{N} c_k P_k`, and represented as a PennyLane Hamiltonian.
+
+    Args:
+        obs_a (Hamiltonian): first observable
+        obs_b (Hamiltonian): second observable
+
+    Returns:
+        qml.Hamiltonian: Observable expressed as a PennyLane Hamiltonian
+
+    **Example**
+
+    >>> c = np.array([0.5, 0.5])
+    >>> obs_a = qml.Hamiltonian(c, [qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliZ(1)])
+    >>> obs_b = qml.Hamiltonian(c, [qml.PauliX(0) @ qml.PauliX(1), qml.PauliZ(0) @ qml.PauliZ(1)])
+    >>> print(_observable_mult(obs_a, obs_b))
+      (-0.25j) [Z1]
+    + (-0.25j) [Y0]
+    + ( 0.25j) [Y1]
+    + ((0.25+0j)) [Y0 X1]
+    """
+    o = []
+    c = []
+    for i in range(len(obs_a.terms[0])):
+        for j in range(len(obs_b.terms[0])):
+            op, phase = qml.grouping.pauli_mult_with_phase(obs_a.terms[1][i], obs_b.terms[1][j])
+            o.append(op)
+            c.append(phase * obs_a.terms[0][i] * obs_b.terms[0][j])
+
+    return _simplify(qml.Hamiltonian(qml.math.stack(c), o))
 
-    return generators, pauli_x_ops
+
+def _simplify(h, cutoff=1.0e-12):
+    r"""Add together identical terms in the Hamiltonian.
+
+    The Hamiltonian terms with identical Pauli words are added together and eliminated if the
+    overall coefficient is zero.
+
+    Args:
+        h (Hamiltonian): PennyLane Hamiltonian
+        cutoff (float): cutoff value for discarding the negligible terms
+
+    Returns:
+        Hamiltonian: Simplified PennyLane Hamiltonian
+
+    **Example**
+
+    >>> c = np.array([0.5, 0.5])
+    >>> h = qml.Hamiltonian(c, [qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliY(1)])
+    >>> print(_simplify(h))
+    (1.0) [X0 Y1]
+    """
+    s = []
+    wiremap = dict(zip(h.wires, range(len(h.wires) + 1)))
+    for term in h.terms[1]:
+        term = qml.operation.Tensor(term).prune()
+        s.append(qml.grouping.pauli_word_to_string(term, wire_map=wiremap))
+
+    o = list(set(s))
+    c = [0.0] * len(o)
+    for i, item in enumerate(s):
+        c[o.index(item)] += h.terms[0][i]
+    c = qml.math.stack(c)
+
+    coeffs = []
+    ops = []
+    nonzero_ind = np.argwhere(abs(c) > cutoff).flatten()
+    for i in nonzero_ind:
+        coeffs.append(c[i])
+        ops.append(qml.grouping.string_to_pauli_word(o[i], wire_map=wiremap))
+    try:
+        coeffs = qml.math.stack(coeffs)
+    except ValueError:
+        pass
+
+    return qml.Hamiltonian(coeffs, ops)
+
+
+def clifford(generators, paulix_ops):
+    r"""Compute a Clifford operator from a set of generators and Pauli-X operators.
+
+    This function computes :math:`U = U_0U_1...U_k` for a set of :math:`k` generators and
+    :math:`k` Pauli-X operators.
+
+    Args:
+        generators (list[Hamiltonian]): generators expressed as PennyLane Hamiltonians
+        paulix_ops (list[Operation]): list of single-qubit Pauli-X operators
+
+    Returns:
+        Hamiltonian: Clifford operator expressed as a PennyLane Hamiltonian
+
+    **Example**
+
+    >>> t1 = qml.Hamiltonian([1.0], [qml.grouping.string_to_pauli_word('ZZII')])
+    >>> t2 = qml.Hamiltonian([1.0], [qml.grouping.string_to_pauli_word('ZIZI')])
+    >>> t3 = qml.Hamiltonian([1.0], [qml.grouping.string_to_pauli_word('ZIIZ')])
+    >>> generators = [t1, t2, t3]
+    >>> paulix_ops = [qml.PauliX(1), qml.PauliX(2), qml.PauliX(3)]
+    >>> u = clifford(generators, paulix_ops)
+    >>> print(u)
+      (0.3535533905932737) [Z1 Z2 X3]
+    + (0.3535533905932737) [X1 X2 X3]
+    + (0.3535533905932737) [Z1 X2 Z3]
+    + (0.3535533905932737) [X1 Z2 Z3]
+    + (0.3535533905932737) [Z0 X1 X2 Z3]
+    + (0.3535533905932737) [Z0 Z1 Z2 Z3]
+    + (0.3535533905932737) [Z0 X1 Z2 X3]
+    + (0.3535533905932737) [Z0 Z1 X2 X3]
+    """
+    cliff = []
+    for i, t in enumerate(generators):
+        cliff.append(1 / 2 ** 0.5 * (paulix_ops[i] + t))
+
+    u = functools.reduce(lambda i, j: _observable_mult(i, j), cliff)
+
+    return u
+
+
+def transform_hamiltonian(h, generators, paulix_ops, paulix_sector):
+    r"""Transform a Hamiltonian with a Clifford operator and taper qubits.
+
+    The Hamiltonian is transformed as :math:`H' = U^{\dagger} H U` where :math:`U` is a Clifford
+    operator. The transformed Hamiltonian acts trivially on some qubits which are then replaced
+    with the eigenvalues of their corresponding Pauli-X operator. The list of these
+    eigenvalues is defined as the Pauli sector.
+
+    Args:
+        h (Hamiltonian): PennyLane Hamiltonian
+        generators (list[Hamiltonian]): generators expressed as PennyLane Hamiltonians
+        paulix_ops (list[Operation]): list of single-qubit Pauli-X operators
+        paulix_sector (llist[int]): eigenvalues of the Pauli-X operators
+
+    Returns:
+        Hamiltonian: the tapered Hamiltonian
+
+    **Example**
+
+    >>> symbols = ["H", "H"]
+    >>> geometry = np.array([[0.0, 0.0, -0.69440367], [0.0, 0.0, 0.69440367]])
+    >>> mol = qml.hf.Molecule(symbols, geometry)
+    >>> H = qml.hf.generate_hamiltonian(mol)(geometry)
+    >>> generators, paulix_ops = qml.hf.generate_symmetries(H, len(H.wires))
+    >>> paulix_sector = [1, -1, -1]
+    >>> H_tapered = transform_hamiltonian(H, generators, paulix_ops, paulix_sector)
+    >>> print(H_tapered)
+      ((-0.321034397355757+0j)) [I0]
+    + ((0.1809270275619003+0j)) [X0]
+    + ((0.7959678503869626+0j)) [Z0]
+    """
+    u = clifford(generators, paulix_ops)
+    h = _observable_mult(_observable_mult(u, h), u)
+
+    val = np.ones(len(h.terms[0])) * complex(1.0)
+
+    wiremap = dict(zip(h.wires, range(len(h.wires) + 1)))
+    paulix_wires = [x.wires[0] for x in paulix_ops]
+    for idx, w in enumerate(paulix_wires):
+        for i in range(len(h.terms[0])):
+            s = qml.grouping.pauli_word_to_string(h.terms[1][i], wire_map=wiremap)
+            if s[w] == "X":
+                val[i] *= paulix_sector[idx]
+
+    o = []
+    wires_tap = [h.wires[i] for i in range(len(h.wires)) if i not in paulix_wires]
+    wiremap_tap = dict(zip(wires_tap, range(len(wires_tap) + 1)))
+    for i in range(len(h.terms[0])):
+        s = qml.grouping.pauli_word_to_string(h.terms[1][i], wire_map=wiremap)
+        wires = [x for x in h.wires if x not in paulix_wires]
+        o.append(
+            qml.grouping.string_to_pauli_word("".join([s[i] for i in wires]), wire_map=wiremap_tap)
+        )
+
+    c = anp.multiply(val, h.terms[0])
+    c = qml.math.stack(c)
+
+    return _simplify(qml.Hamiltonian(c, o))


Test patch:
diff --git a/tests/hf/test_tapering.py b/tests/hf/test_tapering.py
index 0d77b5ea4..6511e91a3 100644
--- a/tests/hf/test_tapering.py
+++ b/tests/hf/test_tapering.py
@@ -14,17 +14,20 @@
 """
 Unit tests for functions needed for qubit tapering.
 """
-
-import pytest
 import pennylane as qml
+import pytest
 from pennylane import numpy as np
 from pennylane.hf.tapering import (
     _binary_matrix,
-    _reduced_row_echelon,
     _kernel,
-    get_generators,
+    _observable_mult,
+    _reduced_row_echelon,
+    _simplify,
+    clifford,
     generate_paulis,
     generate_symmetries,
+    get_generators,
+    transform_hamiltonian,
 )
 
 
@@ -355,3 +358,144 @@ def test_generate_symmetries(hamiltonian, num_qubits, res_generators, res_pauli_
 
     for p1, p2 in zip(pauli_ops, res_pauli_ops):
         assert p1.compare(p2)
+
+
+@pytest.mark.parametrize(
+    ("obs_a", "obs_b", "result"),
+    [
+        (
+            qml.Hamiltonian(np.array([-1.0]), [qml.PauliX(0) @ qml.PauliY(1) @ qml.PauliX(2)]),
+            qml.Hamiltonian(np.array([-1.0]), [qml.PauliX(0) @ qml.PauliY(1) @ qml.PauliX(2)]),
+            qml.Hamiltonian(np.array([1.0]), [qml.Identity(0)]),
+        ),
+        (
+            qml.Hamiltonian(
+                np.array([0.5, 0.5]), [qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliZ(1)]
+            ),
+            qml.Hamiltonian(
+                np.array([0.5, 0.5]), [qml.PauliX(0) @ qml.PauliX(1), qml.PauliZ(0) @ qml.PauliZ(1)]
+            ),
+            qml.Hamiltonian(
+                np.array([-0.25j, 0.25j, -0.25j, 0.25]),
+                [qml.PauliY(0), qml.PauliY(1), qml.PauliZ(1), qml.PauliY(0) @ qml.PauliX(1)],
+            ),
+        ),
+    ],
+)
+def test_observable_mult(obs_a, obs_b, result):
+    r"""Test that observable_mult returns the correct result."""
+    o = _observable_mult(obs_a, obs_b)
+    assert o.compare(result)
+
+
+@pytest.mark.parametrize(
+    ("hamiltonian", "result"),
+    [
+        (
+            qml.Hamiltonian(
+                np.array([0.5, 0.5]), [qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliY(1)]
+            ),
+            qml.Hamiltonian(np.array([1.0]), [qml.PauliX(0) @ qml.PauliY(1)]),
+        ),
+        (
+            qml.Hamiltonian(
+                np.array([0.5, -0.5]),
+                [qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliY(1)],
+            ),
+            qml.Hamiltonian([], []),
+        ),
+        (
+            qml.Hamiltonian(
+                np.array([0.0, -0.5]),
+                [qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliZ(1)],
+            ),
+            qml.Hamiltonian(np.array([-0.5]), [qml.PauliX(0) @ qml.PauliZ(1)]),
+        ),
+        (
+            qml.Hamiltonian(
+                np.array([0.25, 0.25, 0.25, -0.25]),
+                [
+                    qml.PauliX(0) @ qml.PauliY(1),
+                    qml.PauliX(0) @ qml.PauliZ(1),
+                    qml.PauliX(0) @ qml.PauliY(1),
+                    qml.PauliX(0) @ qml.PauliY(1),
+                ],
+            ),
+            qml.Hamiltonian(
+                np.array([0.25, 0.25]),
+                [qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliZ(1)],
+            ),
+        ),
+    ],
+)
+def test_simplify(hamiltonian, result):
+    r"""Test that simplify returns the correct hamiltonian."""
+    h = _simplify(hamiltonian)
+    assert h.compare(result)
+
+
+@pytest.mark.parametrize(
+    ("generator", "paulix_ops", "result"),
+    [
+        (
+            [
+                qml.Hamiltonian(np.array([1.0]), [qml.PauliZ(0)]),
+                qml.Hamiltonian(np.array([1.0]), [qml.PauliZ(1)]),
+            ],
+            [qml.PauliX(0), qml.PauliX(1)],
+            qml.Hamiltonian(
+                np.array(
+                    [
+                        (1 / np.sqrt(2)) ** 2,
+                        (1 / np.sqrt(2)) ** 2,
+                        (1 / np.sqrt(2)) ** 2,
+                        (1 / np.sqrt(2)) ** 2,
+                    ]
+                ),
+                [
+                    qml.PauliZ(0) @ qml.PauliZ(1),
+                    qml.PauliZ(0) @ qml.PauliX(1),
+                    qml.PauliX(0) @ qml.PauliZ(1),
+                    qml.PauliX(0) @ qml.PauliX(1),
+                ],
+            ),
+        ),
+    ],
+)
+def test_cliford(generator, paulix_ops, result):
+    r"""Test that clifford returns the correct operator."""
+    u = clifford(generator, paulix_ops)
+    assert u.compare(result)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "generator", "paulix_ops", "paulix_sector", "ham_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, -0.69440367], [0.0, 0.0, 0.69440367]], requires_grad=False),
+            [
+                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(1), qml.PauliX(2), qml.PauliX(3)],
+            [1, -1, -1],
+            qml.Hamiltonian(
+                np.array([-0.3210344, 0.18092703, 0.79596785]),
+                [qml.Identity(0), qml.PauliX(0), qml.PauliZ(0)],
+            ),
+        ),
+    ],
+)
+def test_transform_hamiltonian(symbols, geometry, generator, paulix_ops, paulix_sector, ham_ref):
+    r"""Test that transform_hamiltonian returns the correct hamiltonian."""
+    mol = qml.hf.Molecule(symbols, geometry)
+    h = qml.hf.generate_hamiltonian(mol)()
+    ham_calc = transform_hamiltonian(h, generator, paulix_ops, paulix_sector)
+
+    # sort Hamiltonian terms and then compare with reference
+    sorted_terms = list(sorted(zip(ham_calc.terms[0], ham_calc.terms[1])))
+    for i, term in enumerate(sorted_terms):
+        assert np.allclose(term[0], ham_ref.terms[0][i])
+        assert term[1].compare(ham_ref.terms[1][i])