Lowering qubit requirements using binary codes

View on QuantumAI Run in Google Colab View source on GitHub Download notebook

Setup

Install the OpenFermion package:

try:
    import openfermion
except ImportError:
    !pip install git+https://github.com/quantumlib/OpenFermion.git@master#egg=openfermion

Introduction

Molecular Hamiltonians are known to have certain symmetries that are not taken into account by mappings like the Jordan-Wigner or Bravyi-Kitaev transform. The most notable of such symmetries is the conservation of the total number of particles in the system. Since those symmetries effectively reduce the degrees of freedom of the system, one is able to reduce the number of qubits required for simulation by utilizing binary codes (arXiv:1712.07067).

We can represent the symmetry-reduced Fermion basis by binary vectors of a set \(\mathcal{V} \ni \boldsymbol{\nu}\), with \( \boldsymbol{\nu} = (\nu_0, \, \nu_1, \dots, \, \nu_{N-1} ) \), where every component \(\nu_i \in \lbrace 0, 1 \rbrace \) and \(N\) is the total number of Fermion modes. These binary vectors \( \boldsymbol{\nu}\) are related to the actual basis states by: \( \left[\prod_{i=0}^{N-1} (a_i^{\dagger})^{\nu_i} \right] \left|{\text{vac} }\right\rangle \, , \) where \( (a_i^\dagger)^{0}=1\). The qubit basis, on the other hand, can be characterized by length-\(n\) binary vectors \(\boldsymbol{\omega}=(\omega_0, \, \dots , \, \omega_{n-1})\), that represent an \(n\)-qubit basis state by:

\[ \left|{\omega_0}\right\rangle \otimes \left|\omega_1\right\rangle \otimes \dots \otimes \left|{\omega_{n-1} }\right\rangle \, . \]

Since \(\mathcal{V}\) is a mere subset of the \(N\)-fold binary space, but the set of the vectors \(\boldsymbol{\omega}\) spans the entire \(n\)-fold binary space we can assign every vector \(\boldsymbol{\nu}\) to a vector \( \boldsymbol{\omega}\), such that \(n<N\). This reduces the amount of qubits required by \((N-n)\). The mapping can be done by a binary code, a classical object that consists of an encoder function \(\boldsymbol{e}\) and a decoder function \(\boldsymbol{d}\). These functions relate the binary vectors \(\boldsymbol{e}(\boldsymbol{\nu})=\boldsymbol{\omega}\), \(\boldsymbol{d}(\boldsymbol{\omega})=\boldsymbol{\nu}\), such that \(\boldsymbol{d}(\boldsymbol{e}(\boldsymbol{\nu}))=\boldsymbol{\nu}\). In OpenFermion, we, at the moment, allow for non-linear decoders \(\boldsymbol{d}\) and linear encoders \(\boldsymbol{e}(\boldsymbol{\nu})=A \boldsymbol{\nu}\), where the matrix multiplication with the \((n\times N)\)-binary matrix \(A\) is \((\text{mod 2})\) in every component.

Symbolic binary functions

The non-linear binary functions for the components of the decoder are here modeled by the \(\text{BinaryPolynomial}\) class in openfermion.ops. For initialization we can conveniently use strings ('w0 w1 + w1 +1' for the binary function \(\boldsymbol{\omega} \to \omega_0 \omega_1 + \omega_1 + 1 \;\text{mod 2}\)), the native data structure or symbolic addition and multiplication.

from openfermion.ops import BinaryPolynomial

binary_1 = BinaryPolynomial('w0 w1 + w1 + 1')

print("These three expressions are equivalent: \n", binary_1)
print(BinaryPolynomial('w0') * BinaryPolynomial('w1 + 1') + BinaryPolynomial('1'))
print(BinaryPolynomial([(1, 0), (1, ), ('one', )]))

print('The native data type structure can be seen here:')
print(binary_1.terms)
print('We can always evaluate the expression for instance by the vector (w0, w1, w2) = (1, 0, 0):',
      binary_1.evaluate('100'))
These three expressions are equivalent: 
 [W0 W1] + [W1] + [1]
[W0 W1] + [W0] + [1]
[W0 W1] + [W1] + [1]
The native data type structure can be seen here:
[(0, 1), (1,), ('one',)]
We can always evaluate the expression for instance by the vector (w0, w1, w2) = (1, 0, 0): 1

Binary codes

The \(\text{BinaryCode}\) class bundles a decoder - a list of decoder components, which are instances of \(\text{BinaryPolynomial}\) - and an encoder - the matrix \(A\) as sparse numpy array - as a binary code. The constructor however admits (dense) numpy arrays, nested lists or tuples as input for \(A\), and arrays, lists or tuples of \(\text{BinaryPolynomial}\) objects - or valid inputs for \(\text{BinaryPolynomial}\) constructors - as input for \(\boldsymbol{d}\). An instance of the \(\text{BinaryCode}\) class knows about the number of qubits and the number of modes in the mapping.

from openfermion.ops import BinaryCode

code_1 = BinaryCode([[1, 0, 0], [0, 1, 0]], ['w0', 'w1', 'w0 + w1 + 1' ])

print(code_1)
print('number of qubits: ', code_1.n_qubits, '  number of Fermion modes: ', code_1.n_modes )
print('encoding matrix: \n', code_1.encoder.toarray())
print('decoder: ', code_1.decoder)
[[[1, 0, 0], [0, 1, 0]], '[[W0],[W1],[W0] + [W1] + [1]]']
number of qubits:  2   number of Fermion modes:  3
encoding matrix: 
 [[1 0 0]
 [0 1 0]]
decoder:  [[W0], [W1], [W0] + [W1] + [1]]

The code used in the example above, is in fact the (odd) Checksum code, and is implemented already - along with a few other examples from arxiv:1712.07067. In addition to the \(\text{checksum_code}\) the functions \(\text{weight_one_segment_code}\), \(\text{weight_two_segment_code}\), that output a subcode each, as well as \(\text{weight_one_binary_addressing_code}\) can be found under openfermion.transforms._code_transform_functions.

There are two other ways to construct new codes from the ones given - both of them can be done conveniently with symbolic operations between two code objects \((\boldsymbol{e}, \boldsymbol{d})\) and \((\boldsymbol{e^\prime}, \boldsymbol{d^\prime})\) to yield a new code \((\boldsymbol{e^{\prime\prime} }, \boldsymbol{d^{\prime\prime} })\):

Appendage

Input and output vectors of two codes are appended to each other such that:

\[ e^{\prime\prime}(\boldsymbol{\nu} \oplus \boldsymbol{\nu^{\prime} })=\boldsymbol{e}(\boldsymbol{\nu}) \oplus \boldsymbol{e^\prime}(\boldsymbol{\nu^\prime})\, , \qquad d^{\prime\prime}(\boldsymbol{\omega} \oplus \boldsymbol{\omega^{\prime} })=\boldsymbol{d}(\boldsymbol{\omega}) \oplus \boldsymbol{d^\prime}(\boldsymbol{\omega^\prime}) \, . \]

This is implemented with symbolic addition of two \(\text{BinaryCode}\) objects (using + or += ) or, for appending several instances of the same code at once, multiplication of the \(\text{BinaryCode}\) with an integer. Appending codes is useful when we want to obtain a segment code, or a segmented transform.

Concatenation

Two codes can (if the corresponding vectors match in size) be applied consecutively, in the sense that the output of the encoder of the first code is input to the encoder of the second code. This defines an entirely new encoder, and the corresponding decoder is defined to undo this operation.

\[ \boldsymbol{e^{\prime\prime} }(\boldsymbol{\nu^{\prime\prime} })=\boldsymbol{e^\prime}\left(\boldsymbol{e}(\boldsymbol{\nu^{\prime\prime} }) \right) \, , \qquad \boldsymbol{d^{\prime\prime} }(\boldsymbol{\omega^{\prime\prime} })=\boldsymbol{d}\left(\boldsymbol{d^\prime}(\boldsymbol{\omega^{\prime\prime} }) \right) \]

This is done by symbolic multiplication of two \(\text{BinaryCode}\) instances (with * or *= ). One can concatenate the codes with each other such that additional qubits can be saved (e.g. checksum code * segment code ), or to modify the resulting gates after transform (e.g. checksum code * Bravyi-Kitaev code).

A broad palette of codes is provided to help construct codes symbolically. The \(\text{jordan_wigner_code}\) can be appended to every code to fill the number of modes, concatenating the \(\text{bravyi_kitaev_code}\) or \(\text{parity_code}\) will modify the appearance of gates after the transform. The \(\text{interleaved_code}\) is useful to concatenate appended codes with if in Hamiltonians, Fermion operators are ordered by spin indexing even-odd (up-down-up-down-up ...) . This particular instance is used in the demonstration below.

Before we turn to describe the transformation, a word of warning has to be spoken here. Controlled gates that occur in the Hamiltonian by using non-linear codes are decomposed into Pauli strings, e.g. \(\text{CPHASE}(1,2)=\frac{1}{2}(1+Z_1+Z_2-Z_1Z_2)\). In that way the amount of terms in a Hamiltonian might rise exponentially, if one chooses to use strongly non-linear codes.

Operator transform

The actual transform of Fermion operators into qubit operators is done with the routine \(\text{binary_code_transform}\), that takes a Hamiltonian and a suitable code as inputs, outputting a qubit Hamiltonian.
Let us consider the case of a molecule with 4 modes where, due to the absence of magnetic interactions, the set of valid modes is only \( \mathcal{V}=\lbrace (1,\, 1,\, 0,\, 0 ),\,(1,\, 0,\, 0,\, 1 ),\,(0,\, 1,\, 1,\, 0 ),\,(0,\, 0,\, 1,\, 1 )\rbrace \, .\) One can either use an (even weight) checksum code to save a single qubit, or use and (odd weight) checksum code on spin-up and -down modes each to save two qubits. Since the ordering is even-odd, however, this requires to concatenate the with the interleaved code, which switches the spin indexing of the qubits from even-odd ordering to up-then-down. Instead of using the interleaved code, we can also use the reorder function to apply up-then-down ordering on the hamiltonian.

from openfermion.transforms import *
from openfermion.chem import MolecularData
from openfermion.transforms import binary_code_transform
from openfermion.transforms import get_fermion_operator
from openfermion.linalg import eigenspectrum
from openfermion.transforms import normal_ordered, reorder
from openfermion.utils import up_then_down

def LiH_hamiltonian():
    geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., 1.45))]
    molecule = MolecularData(geometry, 'sto-3g', 1,
                             description="1.45")
    molecule.load()
    molecular_hamiltonian = molecule.get_molecular_hamiltonian(occupied_indices = [0], active_indices = [1,2])
    hamiltonian = normal_ordered(get_fermion_operator(molecular_hamiltonian))
    return hamiltonian

hamiltonian = LiH_hamiltonian()
print('Fermionic Hamiltonian')
print (hamiltonian)
print("The eigenspectrum")
print(eigenspectrum(hamiltonian))

print('\n-----\n')
jw = binary_code_transform(hamiltonian, jordan_wigner_code(4))
print('Jordan-Wigner transformed Hamiltonian')
print(jw)
print("the eigenspectrum of the transformed hamiltonian")
print(eigenspectrum(jw))

print('\n-----\n')
cksm_save_one = binary_code_transform(hamiltonian, checksum_code(4,0))
print('Even-weight checksum code')
print(cksm_save_one)
print("the eigenspectrum of the transformed hamiltonian")
print(eigenspectrum(cksm_save_one))

print('\n-----\n')
up_down_save_two = binary_code_transform(hamiltonian, interleaved_code(4)*(2*checksum_code(2,1)))
print('Double odd-weight checksum codes')
print(up_down_save_two )
print("the eigenspectrum of the transformed hamiltonian")
print(eigenspectrum(up_down_save_two ))

print('\n-----\n')
print('Instead of interleaving, we can apply up-then-down ordering using the reorder function:')
up_down_save_two = binary_code_transform(reorder(hamiltonian,up_then_down), 2*checksum_code(2,1))
print(up_down_save_two)
print("the eigenspectrum of the transformed hamiltonian")
print(eigenspectrum(up_down_save_two))
Fermionic Hamiltonian
-6.7698132180879735 [] +
-0.7952726864779313 [0^ 0] +
0.04614563473199314 [0^ 2] +
-0.4977908053255035 [1^ 0^ 1 0] +
-0.046145652803099894 [1^ 0^ 2 1] +
0.046145652803099894 [1^ 0^ 3 0] +
-0.011731985763800887 [1^ 0^ 3 2] +
-0.7952726864779313 [1^ 1] +
0.04614563473199314 [1^ 3] +
0.04614563473199324 [2^ 0] +
-0.21652178317319534 [2^ 0^ 2 0] +
-0.04614565280309991 [2^ 1^ 1 0] +
-0.22825376893699628 [2^ 1^ 2 1] +
0.011731985763800913 [2^ 1^ 3 0] +
0.005497504431583474 [2^ 1^ 3 2] +
-0.36549257026798354 [2^ 2] +
0.04614565280309991 [3^ 0^ 1 0] +
0.011731985763800913 [3^ 0^ 2 1] +
-0.22825376893699628 [3^ 0^ 3 0] +
-0.005497504431583474 [3^ 0^ 3 2] +
0.04614563473199324 [3^ 1] +
-0.21652178317319534 [3^ 1^ 3 1] +
-0.011731985763800917 [3^ 2^ 1 0] +
0.005497504431583467 [3^ 2^ 2 1] +
-0.005497504431583467 [3^ 2^ 3 0] +
-0.33918438174683924 [3^ 2^ 3 2] +
-0.36549257026798354 [3^ 3]
The eigenspectrum
[-7.86277316 -7.78339621 -7.78339621 -7.71405669 -7.71405669 -7.71405669
 -7.70047584 -7.56998474 -7.56998474 -7.51199971 -7.51199971 -7.36481744
 -7.15152548 -7.13040696 -7.13040696 -6.76981322]

-----

Jordan-Wigner transformed Hamiltonian
-7.498946902010707 [] +
-0.002932996440950227 [X0 X1 Y2 Y3] +
0.002932996440950227 [X0 Y1 Y2 X3] +
0.012910780273117489 [X0 Z1 X2] +
-0.0013743761078958677 [X0 Z1 X2 Z3] +
0.011536413200774975 [X0 X2] +
0.002932996440950227 [Y0 X1 X2 Y3] +
-0.002932996440950227 [Y0 Y1 X2 X3] +
0.012910780273117489 [Y0 Z1 Y2] +
-0.0013743761078958677 [Y0 Z1 Y2 Z3] +
0.011536413200774975 [Y0 Y2] +
0.1619947538800418 [Z0] +
0.011536413200774975 [Z0 X1 Z2 X3] +
0.011536413200774975 [Z0 Y1 Z2 Y3] +
0.12444770133137588 [Z0 Z1] +
0.054130445793298836 [Z0 Z2] +
0.05706344223424907 [Z0 Z3] +
0.012910780273117489 [X1 Z2 X3] +
-0.0013743761078958677 [X1 X3] +
0.012910780273117489 [Y1 Z2 Y3] +
-0.0013743761078958677 [Y1 Y3] +
0.1619947538800418 [Z1] +
0.05706344223424907 [Z1 Z2] +
0.054130445793298836 [Z1 Z3] +
-0.013243698330265952 [Z2] +
0.08479609543670981 [Z2 Z3] +
-0.013243698330265966 [Z3]
the eigenspectrum of the transformed hamiltonian
[-7.86277316 -7.78339621 -7.78339621 -7.71405669 -7.71405669 -7.71405669
 -7.70047584 -7.56998474 -7.56998474 -7.51199971 -7.51199971 -7.36481744
 -7.15152548 -7.13040696 -7.13040696 -6.76981322]

-----

Even-weight checksum code
-7.498946902010707 [] +
0.005865992881900454 [X0 Y1 Y2] +
0.012910780273117489 [X0 Z1 X2] +
0.012910789308670843 [X0 X2] +
-0.005865992881900454 [Y0 Y1 X2] +
0.012910780273117489 [Y0 Z1 Y2] +
0.012910789308670843 [Y0 Y2] +
0.1619947538800418 [Z0] +
-0.012910780273117489 [Z0 X1] +
0.012910789308670843 [Z0 X1 Z2] +
0.2092437967680857 [Z0 Z1] +
-0.013243698330265966 [Z0 Z1 Z2] +
0.10826089158659767 [Z0 Z2] +
-0.012910789308670843 [X1] +
0.012910780273117489 [X1 Z2] +
0.1619947538800418 [Z1] +
0.11412688446849814 [Z1 Z2] +
-0.013243698330265952 [Z2]
the eigenspectrum of the transformed hamiltonian
[-7.86277316 -7.71405669 -7.71405669 -7.71405669 -7.70047584 -7.36481744
 -7.15152548 -6.76981322]

-----

Double odd-weight checksum codes
-7.607207793597305 [] +
0.025821578617341686 [X0] +
0.025821560546234978 [X0 Z1] +
-0.011731985763800908 [Y0 Y1] +
0.1752384522103078 [Z0] +
-0.025821560546234978 [Z0 X1] +
0.09511691229958755 [Z0 Z1] +
-0.025821578617341686 [X1] +
0.1752384522103078 [Z1]
the eigenspectrum of the transformed hamiltonian
[-7.86277316 -7.71405669 -7.70047584 -7.15152548]

-----

Instead of interleaving, we can apply up-then-down ordering using the reorder function:
-7.607207793597305 [] +
0.025821560546234978 [X0] +
0.011731985763800908 [X0 X1] +
0.025821578617341686 [X0 Z1] +
0.1752384522103078 [Z0] +
0.025821578617341686 [Z0 X1] +
0.09511691229958755 [Z0 Z1] +
0.025821560546234978 [X1] +
0.1752384522103078 [Z1]
the eigenspectrum of the transformed hamiltonian
[-7.86277316 -7.71405669 -7.70047584 -7.15152548]