Add Bravyi-Kitaev mapping (#5390)

**Context:**
The PR adds Bravyi-Kitaev mapping

**Description of the Change:**

**Benefits:**
New mapping scheme from fermionic to qubit Hamiltonian

**Possible Drawbacks:**

**Related GitHub Issues:**

---------

Co-authored-by: soranjh <40344468+soranjh@users.noreply.github.com>
Co-authored-by: Austin Huang <65315367+austingmhuang@users.noreply.github.com>
Co-authored-by: Thomas R. Bromley <49409390+trbromley@users.noreply.github.com>

doc/releases/changelog-dev.md

@@ -134,6 +134,22 @@
 
   ```
 
+
+* Added new function `qml.bravyi_kitaev` to map fermionic Hamiltonians to qubit Hamiltonians.
+  [(#5390)](https://github.com/PennyLaneAI/pennylane/pull/5390)
+
+  ```python
+  import pennylane as qml
+  fermi_ham = qml.fermi.from_string('0+ 1-')
+
+  qubit_ham = qml.bravyi_kitaev(fermi_ham, n=6)
+  ```
+
+  ```pycon
+  >>> print(qubit_ham)
+  -0.25j * Y(0.0) + (-0.25+0j) * X(0) @ Z(1.0) + (0.25+0j) * X(0.0) + 0.25j * Y(0) @ Z(1.0)
+  ```
+
 * A new class `qml.ops.LinearCombination` is introduced. In essence, this class is an updated equivalent of `qml.ops.Hamiltonian`
   but for usage with new operator arithmetic.
   [(#5216)](https://github.com/PennyLaneAI/pennylane/pull/5216)
@@ -381,6 +397,7 @@ Mikhail Andrenkov,
 Utkarsh Azad,
 Gabriel Bottrill,
 Astral Cai,
+Diksha Dhawan,
 Isaac De Vlugt,
 Amintor Dusko,
 Pietropaolo Frisoni,

pennylane/__init__.py

@@ -36,7 +36,13 @@
 from pennylane.resource import specs
 import pennylane.resource
 import pennylane.qchem
-from pennylane.fermi import FermiC, FermiA, jordan_wigner
+from pennylane.fermi import (
+    FermiC,
+    FermiA,
+    jordan_wigner,
+    parity_transform,
+    bravyi_kitaev,
+)
 from pennylane.qchem import (
     taper,
     symmetry_generators,

pennylane/fermi/__init__.py

@@ -15,5 +15,5 @@
 fermionic operators. """
 
 
-from .conversion import jordan_wigner, parity_transform
+from .conversion import jordan_wigner, parity_transform, bravyi_kitaev
 from .fermionic import FermiWord, FermiC, FermiA, FermiSentence, from_string

pennylane/fermi/conversion.py

@@ -16,6 +16,7 @@
 from functools import singledispatch
 from typing import Union
 
+import numpy as np
 import pennylane as qml
 from pennylane.operation import Operator
 from pennylane.pauli import PauliSentence, PauliWord
@@ -62,7 +63,7 @@ def jordan_wigner(
 
     **Example**
 
-    >>> w = FermiWord({(0, 0) : '+', (1, 1) : '-'})
+    >>> w = qml.fermi.from_string('0+ 1-')
     >>> jordan_wigner(w)
     (
         -0.25j * (Y(0) @ X(1))
@@ -145,6 +146,8 @@ def _(fermi_operator: FermiSentence, ps=False, wire_map=None, tol=None):
             if tol is not None and abs(qml.math.imag(qubit_operator[pw])) <= tol:
                 qubit_operator[pw] = qml.math.real(qubit_operator[pw])
 
+    qubit_operator.simplify(tol=1e-16)
+
     if not ps:
         qubit_operator = qubit_operator.operation(wire_order=[identity_wire])
 
@@ -203,7 +206,7 @@ def parity_transform(
 
     **Example**
 
-    >>> w = FermiWord({(0, 0) : '+', (1, 1) : '-'})
+    >>> w = qml.fermi.from_string('0+ 1-')
     >>> parity_transform(w, n=6)
     (
         -0.25j * Y(0)
@@ -285,7 +288,277 @@ def _(fermi_operator: FermiSentence, n, ps=False, wire_map=None, tol=None):
         fermi_word_as_ps = parity_transform(fw, n, ps=True)
 
         for pw in fermi_word_as_ps:
-            qubit_operator[pw] += fermi_word_as_ps[pw] * coeff
+            qubit_operator[pw] = qubit_operator[pw] + fermi_word_as_ps[pw] * coeff
+
+            if tol is not None and abs(qml.math.imag(qubit_operator[pw])) <= tol:
+                qubit_operator[pw] = qml.math.real(qubit_operator[pw])
+
+    qubit_operator.simplify(tol=1e-16)
+
+    if not ps:
+        qubit_operator = qubit_operator.operation(wire_order=[identity_wire])
+
+    if wire_map:
+        return qubit_operator.map_wires(wire_map)
+
+    return qubit_operator
+
+
+def bravyi_kitaev(
+    fermi_operator: Union[FermiWord, FermiSentence],
+    n: int,
+    ps: bool = False,
+    wire_map: dict = None,
+    tol: float = None,
+) -> Union[Operator, PauliSentence]:
+    r"""Convert a fermionic operator to a qubit operator using the Bravyi-Kitaev mapping.
+
+    .. note::
+
+        Hamiltonians created with this mapping should be used with operators and states that are
+        compatible with the Bravyi-Kitaev basis.
+
+    In the Bravyi-Kitaev mapping, both occupation number and parity of the orbitals are stored non-locally.
+    In comparison, :func:`~.jordan_wigner` stores the occupation number locally while storing the parity
+    non-locally and vice-versa for :func:`~.parity_transform`. In the Bravyi-Kitaev mapping, the
+    fermionic creation and annihilation operators for even-labelled orbitals are mapped to the Pauli operators as
+
+    .. math::
+
+        \begin{align*}
+           a^{\dagger}_0 &= \frac{1}{2} \left ( X_0  -iY_{0} \right ), \\\\
+           a^{\dagger}_n &= \frac{1}{2} \left ( X_{U(n)} \otimes X_n \otimes Z_{P(n)} -iX_{U(n)} \otimes Y_{n} \otimes Z_{P(n)}\right ), \\\\
+        \end{align*}
+
+    and
+
+    .. math::
+        \begin{align*}
+           a_0 &= \frac{1}{2} \left ( X_0  + iY_{0} \right ), \\\\
+           a_n &= \frac{1}{2} \left ( X_{U(n)} \otimes X_n \otimes Z_{P(n)} +iX_{U(n)} \otimes Y_{n} \otimes Z_{P(n)}\right ). \\\\
+        \end{align*}
+
+    Similarly, the fermionic creation and annihilation operators for odd-labelled orbitals are mapped to the Pauli operators as
+
+    .. math::
+        \begin{align*}
+           a^{\dagger}_n &= \frac{1}{2} \left ( X_{U(n)} \otimes X_n \otimes Z_{P(n)} -iX_{U(n)} \otimes Y_{n} \otimes Z_{R(n)}\right ), \\\\
+        \end{align*}
+
+    and
+
+    .. math::
+        \begin{align*}
+           a_n &= \frac{1}{2} \left ( X_{U(n)} \otimes X_n \otimes Z_{P(n)} +iX_{U(n)} \otimes Y_{n} \otimes Z_{R(n)}\right ), \\\\
+        \end{align*}
+
+    where :math:`X`, :math:`Y`, and :math:`Z` are the Pauli operators, and :math:`U(n)`, :math:`P(n)` and :math:`R(n)`
+    represent the update, parity and remainder sets, respectively [`arXiv:1812.02233 <https://arxiv.org/abs/1812.02233>`_].
+
+    Args:
+        fermi_operator(FermiWord, FermiSentence): the fermionic operator
+        n (int): number of qubits, i.e., spin orbitals in the system
+        ps (bool): whether to return the result as a PauliSentence instead of an
+            Operator. Defaults to False.
+        wire_map (dict): a dictionary defining how to map the orbitals of
+            the Fermi operator to qubit wires. If None, the integers used to
+            order the orbitals will be used as wire labels. Defaults to None.
+        tol (float): tolerance for discarding the imaginary part of the coefficients
+
+    Returns:
+        Union[PauliSentence, Operator]: a linear combination of qubit operators
+
+    **Example**
+
+    >>> w = qml.fermi.from_string('0+ 1-')
+    >>> bravyi_kitaev(w, n=6)
+    (
+    -0.25j * Y(0)
+    + (-0.25+0j) * (X(0) @ Z(1))
+    + (0.25+0j) * X(0)
+    + 0.25j * (Y(0) @ Z(1))
+    )
+
+    >>> bravyi_kitaev(w, n=6, ps=True)
+    -0.25j * Y(0)
+    + (-0.25+0j) * X(0) @ Z(1)
+    + (0.25+0j) * X(0)
+    + 0.25j * Y(0) @ Z(1)
+
+    >>> bravyi_kitaev(w, n=6, ps=True, wire_map={0: 2, 1: 3})
+    -0.25j * Y(2)
+    + (-0.25+0j) * X(2) @ Z(3)
+    + (0.25+0j) * X(2)
+    + 0.25j * Y(2) @ Z(3)
+
+    """
+    return _bravyi_kitaev_dispatch(fermi_operator, n, ps, wire_map, tol)
+
+
+def _update_set(j, bin_range, n):
+    """
+    Computes the update set of the j-th orbital in n qubits.
+
+    Args:
+        j (int) : the orbital index
+        bin_range (int) : smallest power of 2 equal to or greater than
+                          given number of qubits, e.g., 8 for 5 qubits
+        n (int) : number of qubits
+
+    Returns:
+        numpy.ndarray: Array containing the update set
+    """
+
+    indices = np.array([], dtype=int)
+    midpoint = int(bin_range / 2)
+    if bin_range % 2 != 0:
+        return indices
+
+    if j < midpoint:
+        indices = np.append(indices, np.append(bin_range - 1, _update_set(j, midpoint, n)))
+    else:
+        indices = np.append(indices, _update_set(j - midpoint, midpoint, n) + midpoint)
+
+    indices = np.array([u for u in indices if u < n])
+    return indices
+
+
+def _parity_set(j, bin_range):
+    """
+    Computes the parity set of the j-th orbital in n qubits.
+
+    Args:
+        j (int) : the orbital index
+        bin_range (int) : smallest power of 2 equal to or greater than
+                          given number of qubits, e.g., 8 for 5 qubits
+
+    Returns:
+        numpy.ndarray: Array of qubits which determine the parity of qubit j
+    """
+
+    indices = np.array([], dtype=int)
+    midpoint = int(bin_range / 2)
+    if bin_range % 2 != 0:
+        return indices
+
+    if j < midpoint:
+        indices = np.append(indices, _parity_set(j, midpoint))
+    else:
+        indices = np.append(
+            indices,
+            np.append(
+                _parity_set(j - midpoint, midpoint) + midpoint,
+                midpoint - 1,
+            ),
+        )
+
+    return indices
+
+
+def _flip_set(j, bin_range):
+    """
+    Computes the flip set of the j-th orbital in n qubits.
+
+    Args:
+        j (int) : the orbital index
+        bin_range (int) : smallest power of 2 equal to or greater than
+                          given number of qubits, e.g., 8 for 5 qubits
+
+    Returns:
+        numpy.ndarray: Array containing information if the phase of orbital j is same as qubit j.
+    """
+
+    indices = np.array([])
+    midpoint = int(bin_range / 2)
+    if bin_range % 2 != 0:
+        return indices
+
+    if j < midpoint:
+        indices = np.append(indices, _flip_set(j, midpoint))
+    elif midpoint <= j < bin_range - 1:
+        indices = np.append(indices, _flip_set(j - midpoint, midpoint) + midpoint)
+    else:
+        indices = np.append(
+            np.append(indices, _flip_set(j - midpoint, midpoint) + midpoint),
+            midpoint - 1,
+        )
+    return indices
+
+
+@singledispatch
+def _bravyi_kitaev_dispatch(fermi_operator, n, ps, wire_map, tol):
+    """Dispatches to appropriate function if fermi_operator is a FermiWord or FermiSentence."""
+    raise ValueError(f"fermi_operator must be a FermiWord or FermiSentence, got: {fermi_operator}")
+
+
+@_bravyi_kitaev_dispatch.register
+def _(fermi_operator: FermiWord, n, ps=False, wire_map=None, tol=None):
+    wires = list(fermi_operator.wires) or [0]
+    identity_wire = wires[0]
+
+    coeffs = {"+": -0.5j, "-": 0.5j}
+    qubit_operator = PauliSentence({PauliWord({}): 1.0})  # Identity PS to multiply PSs with
+
+    bin_range = int(2 ** np.ceil(np.log2(n)))
+
+    for (_, wire), sign in fermi_operator.items():
+        if wire >= n:
+            raise ValueError(
+                f"Can't create or annihilate a particle on qubit index {wire} for a system with only {n} qubits."
+            )
+
+        u_set = _update_set(wire, bin_range, n)
+        update_string = dict(zip(u_set, ["X"] * len(u_set)))
+
+        p_set = _parity_set(wire, bin_range)
+        parity_string = dict(zip(p_set, ["Z"] * len(p_set)))
+
+        if wire % 2 == 0:
+            qubit_operator @= PauliSentence(
+                {
+                    PauliWord({**parity_string, **{wire: "X"}, **update_string}): 0.5,
+                    PauliWord({**parity_string, **{wire: "Y"}, **update_string}): coeffs[sign],
+                }
+            )
+        else:
+            f_set = _flip_set(wire, bin_range)
+
+            r_set = np.setdiff1d(p_set, f_set)
+            remainder_string = dict(zip(r_set, ["Z"] * len(r_set)))
+
+            qubit_operator @= PauliSentence(
+                {
+                    PauliWord({**parity_string, **{wire: "X"}, **update_string}): 0.5,
+                    PauliWord({**remainder_string, **{wire: "Y"}, **update_string}): coeffs[sign],
+                }
+            )
+
+    for pw in qubit_operator:
+        if tol is not None and abs(qml.math.imag(qubit_operator[pw])) <= tol:
+            qubit_operator[pw] = qml.math.real(qubit_operator[pw])
+
+    if not ps:
+        # wire_order specifies wires to use for Identity (PauliWord({}))
+        qubit_operator = qubit_operator.operation(wire_order=[identity_wire])
+
+    if wire_map:
+        return qubit_operator.map_wires(wire_map)
+
+    return qubit_operator
+
+
+@_bravyi_kitaev_dispatch.register
+def _(fermi_operator: FermiSentence, n, ps=False, wire_map=None, tol=None):
+    wires = list(fermi_operator.wires) or [0]
+    identity_wire = wires[0]
+
+    qubit_operator = PauliSentence()  # Empty PS as 0 operator to add Pws to
+
+    for fw, coeff in fermi_operator.items():
+        fermi_word_as_ps = parity_transform(fw, n, ps=True)
+
+        for pw in fermi_word_as_ps:
+            qubit_operator[pw] = qubit_operator[pw] + fermi_word_as_ps[pw] * coeff
 
             if tol is not None and abs(qml.math.imag(qubit_operator[pw])) <= tol:
                 qubit_operator[pw] = qml.math.real(qubit_operator[pw])