Qubit tapering

Tapering the molecular Hamiltonian

In PennyLane, a molecular Hamiltonian can be created by specifying the atomic symbols and coordinates.

import pennylane as qml
from jax import numpy as jnp
import jax

jax.config.update("jax_enable_x64", True)
symbols = ["He", "H"]
geometry = jnp.array([[0.00000000, 0.00000000, -0.87818361],
                     [0.00000000, 0.00000000,  0.87818362]])

molecule = qml.qchem.Molecule(symbols, geometry, charge=1)
H, qubits = qml.qchem.molecular_hamiltonian(molecule)
H

(
    -1.7101641607058429 * I([0, 1, 2, 3])
  + 0.7306880119481007 * Z(0)
  + 0.05804316091277385 * (Y(0) @ Z(1) @ Y(2))
  + 0.05804316091277385 * (X(0) @ Z(1) @ X(2))
  + 0.21459548889590133 * Z(2)
  + 0.11793755332554709 * (Z(0) @ Z(2))
  + 0.7306880119481007 * Z(1)
  + 0.2363564449774022 * (Z(0) @ Z(1))
  + 0.04383965196887748 * (Y(0) @ Y(2))
  + 0.04383965196887748 * (X(0) @ X(2))
  + 0.05804316091277386 * (Y(1) @ Z(2) @ Y(3))
  + 0.04383965196887748 * (Z(0) @ Y(1) @ Z(2) @ Y(3))
  + 0.05804316091277386 * (X(1) @ Z(2) @ X(3))
  + 0.04383965196887748 * (Z(0) @ X(1) @ Z(2) @ X(3))
  + 0.031707550351754266 * (Y(0) @ X(1) @ X(2) @ Y(3))
  + -0.031707550351754266 * (Y(0) @ Y(1) @ X(2) @ X(3))
  + -0.031707550351754266 * (X(0) @ X(1) @ Y(2) @ Y(3))
  + 0.031707550351754266 * (X(0) @ Y(1) @ Y(2) @ X(3))
  + 0.2145954888959013 * Z(3)
  + 0.14964510367730136 * (Z(0) @ Z(3))
  + -0.01420350937754952 * (Y(0) @ Z(1) @ Y(2) @ Z(3))
  + -0.01420350937754952 * (X(0) @ Z(1) @ X(2) @ Z(3))
  + 0.11793755332554709 * (Z(1) @ Z(3))
  + 0.14964510367730136 * (Z(1) @ Z(2))
  + -0.014203509377549522 * (Y(1) @ Y(3))
  + -0.014203509377549522 * (X(1) @ X(3))
  + 0.18678942625805933 * (Z(2) @ Z(3))
)

This Hamiltonian contains 27 terms where each term acts on up to four qubits.

We can now obtain the symmetry generators and the \(X^{j}\) operators that are used to construct the unitary \(U\) operator that transforms the \(\textrm{HeH}^+\) Hamiltonian. In PennyLane, these are constructed by using the symmetry_generators() and paulix_ops() functions.

generators = qml.symmetry_generators(H)
paulixops = qml.paulix_ops(generators, qubits)

for idx, generator in enumerate(generators):
    print(f"generator {idx+1}: {generator}, paulix_op: {paulixops[idx]}")

generator 1: Z(0) @ Z(2), paulix_op: X(2)
generator 2: Z(1) @ Z(3), paulix_op: X(3)

Once the operator \(U\) is applied, each of the Hamiltonian terms will act on the qubits \(q_2, q_3\) either with the identity or with a Pauli-X operator. For each of these qubits, we can simply replace the Pauli-X operator with one of its eigenvalues \(+1\) or \(-1.\) This results in a total number of \(2^k\) Hamiltonians, where \(k\) is the number of tapered-off qubits and each Hamiltonian corresponds to one eigenvalue sector. The optimal sector corresponding to the ground-state energy of the molecule can be obtained by using the optimal_sector() function.

n_electrons = 2
paulix_sector = qml.qchem.optimal_sector(H, generators, n_electrons)
print(paulix_sector)

[-1, -1]

The optimal eigenvalues are \(-1, -1\) for qubits \(q_2, q_3,\) respectively. We can now build the tapered Hamiltonian with the taper() function which constructs the operator \(U,\) applies it to the Hamiltonian and finally tapers off the qubits \(q_2, q_3\) by replacing the Pauli-X operators acting on those qubits with the optimal eigenvalues.

H_tapered = qml.taper(H, generators, paulixops, paulix_sector)
H_tapered_coeffs, H_tapered_ops = H_tapered.terms()
H_tapered = qml.Hamiltonian(jnp.real(jnp.array(H_tapered_coeffs)), H_tapered_ops)
print(H_tapered)

-1.9460392673569362 * I() + 0.5160925230521993 * Z(1) + 0.5160925230521993 * Z(0) + 0.12385566388085878 * (Z(0) @ Z(1)) + -0.11608632182554765 * (X(0) @ Z(1)) + -0.11608632269285395 * X(0) + 0.11608632269285396 * X(1) + 0.11608632182554766 * (Z(0) @ X(1)) + -0.126830201407017 * (Y(0) @ Y(1))

The new Hamiltonian has only 9 non-zero terms acting on only 2 qubits! We can verify that the original and the tapered Hamiltonian both give the correct ground state energy of the \(\textrm{HeH}^+\) cation, which is \(-2.862595242378\) Ha computed with the full configuration interaction (FCI) method. In PennyLane, it’s possible to build a sparse matrix representation of Hamiltonians. This allows us to directly diagonalize them to obtain exact values of the ground-state energies.

H_sparse = qml.SparseHamiltonian(H.sparse_matrix(), wires=H.wires)
H_tapered_sparse = qml.SparseHamiltonian(H_tapered.sparse_matrix(), wires=H_tapered.wires)

print("Eigenvalues of H:\n", qml.eigvals(H_sparse, k=16))
print("\nEigenvalues of H_tapered:\n", qml.eigvals(H_tapered_sparse, k=4))

Eigenvalues of H:
 [-3.12541987 -3.12541987 -2.86259524 -2.64241998 -2.19672513 -2.19672513
 -2.19672513 -2.18547545 -2.18547545 -2.02874709 -1.35709798 -1.35709798
 -0.69608961 -0.17266334 -0.17266334  1.13871403]

Eigenvalues of H_tapered:
 [-2.86259524 -2.19672513 -2.02874709 -0.69608961]

Note that a second-quantized Hamiltonian is independent of the number of electrons and its eigenspectrum contains the energies of the neutral and charged molecules. Therefore, the smallest eigenvalue returned by eigvals() for a molecular Hamiltonian might correspond to the neutral or charged molecule. While in the case of \(\textrm{HeH}^+,\) qubit tapering allows specifying the optimal sector of the eigenvectors corresponding only to the correct number of electrons, it is generally guaranteed that the optimal sector covers all eigenvectors with the correct number of electrons, but may contain additional eigenvectors of different charge. Therefore, the ground-state energy of the \(\textrm{HeH}^+\) cation is the smallest eigenvalue of the tapered Hamiltonian.

Tapering the reference state

The ground state Hartree-Fock energy of \(\textrm{HeH}^+\) can be computed by directly applying the Hamiltonians to the Hartree-Fock state. For the tapered Hamiltonian, this requires transforming the Hartree-Fock state with the same symmetries obtained for the original Hamiltonian. This reduces the number of qubits in the Hartree-Fock state to match that of the tapered Hamiltonian. It can be done with the taper_hf() function.

state_tapered = qml.qchem.taper_hf(generators, paulixops, paulix_sector,
                                   num_electrons=n_electrons, num_wires=len(H.wires))
print(state_tapered)

[1 1]

Recall that the original Hartree-Fock state for the \(\textrm{HeH}^+\) cation is \([1 1 0 0].\) We can now generate the qubit representation of these states and compute the Hartree-Fock energies for each Hamiltonian.

dev = qml.device("default.qubit", wires=H.wires)
@qml.qnode(dev, interface="jax")
def circuit():
    qml.BasisState(jnp.array([1, 1, 0, 0]), wires=H.wires)
    return qml.state()

qubit_state = circuit()
HF_energy = qubit_state.T @ H.sparse_matrix().toarray() @ qubit_state
print(f"HF energy: {jnp.real(HF_energy):.8f} Ha")

dev = qml.device("lightning.qubit", wires=H_tapered.wires)
@qml.qnode(dev, interface="jax")
def circuit():
    qml.BasisState(jnp.array([1, 1]), wires=H_tapered.wires)
    return qml.state()

qubit_state = circuit()
HF_energy = qubit_state.T @ H_tapered.sparse_matrix().toarray() @ qubit_state
print(f"HF energy (tapered): {jnp.real(HF_energy):.8f} Ha")

HF energy: -2.85436865 Ha
HF energy (tapered): -2.85436865 Ha

These values are identical to the reference Hartree-Fock energy \(-2.8543686493\) Ha.

VQE simulation

Finally, we can use the tapered Hamiltonian and the tapered reference state to perform a VQE simulation and compute the ground-state energy of the \(\textrm{HeH}^+\) cation. We build a tapered variational ansatz [3] that prepares an entangled state by evolving the tapered Hartree-Fock state using the tapered particle-conserving gates, i.e., the SingleExcitation() and DoubleExcitation() operations tapered using taper_operation().

singles, doubles = qml.qchem.excitations(n_electrons, len(H.wires))
tapered_doubles = [
    qml.taper_operation(qml.DoubleExcitation, generators, paulixops, paulix_sector,
                        wire_order=H.wires, op_wires=double) for double in doubles
]
tapered_singles = [
    qml.taper_operation(qml.SingleExcitation, generators, paulixops, paulix_sector,
                        wire_order=H.wires, op_wires=single) for single in singles
]

dev = qml.device("lightning.qubit", wires=H_tapered.wires)

@qml.qnode(dev, interface="jax")
def tapered_circuit(params):
    qml.BasisState(state_tapered, wires=H_tapered.wires)
    for idx, tapered_op in enumerate(tapered_doubles + tapered_singles):
        tapered_op(params[idx])
    return qml.expval(H_tapered)

We define an optimizer and the initial values of the circuit parameters and optimize the circuit parameters with respect to the ground state energy.

import optax
import catalyst

opt = optax.sgd(learning_rate=0.8) # sgd stands for StochasticGradientDescent
init_params = jnp.zeros(len(doubles) + len(singles))

def update_step(i, params, opt_state):
    """Perform a single gradient update step"""
    grads = catalyst.grad(tapered_circuit)(params)
    updates, opt_state = opt.update(grads, opt_state)
    params = optax.apply_updates(params, updates)
    return (params, opt_state)

loss_history = []

opt_state = opt.init(init_params)
params = init_params

for i in range(1, 41):
    params, opt_state = update_step(i, params, opt_state)
    energy = tapered_circuit(params)
    if not i % 5:
        print(f"n: {i}, E: {energy:.8f} Ha, Params: {params}")

n: 5, E: -2.86256216 Ha, Params: [ 0.12622598  0.02570331 -0.02570331]
n: 10, E: -2.86259421 Ha, Params: [ 0.12846324  0.03336569 -0.03336569]
n: 15, E: -2.86259521 Ha, Params: [ 0.12885292  0.03472042 -0.03472042]
n: 20, E: -2.86259524 Ha, Params: [ 0.12892183  0.03496048 -0.03496048]
n: 25, E: -2.86259524 Ha, Params: [ 0.12893404  0.03500304 -0.03500304]
n: 30, E: -2.86259524 Ha, Params: [ 0.1289362   0.03501058 -0.03501058]
n: 35, E: -2.86259524 Ha, Params: [ 0.12893659  0.03501192 -0.03501192]
n: 40, E: -2.86259524 Ha, Params: [ 0.12893666  0.03501216 -0.03501216]

The computed energy matches the FCI energy, \(-2.862595242378\) Ha, while the number of qubits and the number of Hamiltonian terms are significantly reduced with respect to their original values.

Conclusions

Molecular Hamiltonians possess symmetries that can be leveraged to reduce the number of qubits required in quantum computing simulations. This tutorial introduces a PennyLane functionality that can be used for qubit tapering based on \(\mathbb{Z}_2\) symmetries. The procedure includes obtaining tapered Hamiltonians and tapered reference states that can be used in variational quantum algorithms such as VQE.

References

1(1,2)

Sergey Bravyi, Jay M. Gambetta, Antonio Mezzacapo, Kristan Temme, “Tapering off qubits to simulate fermionic Hamiltonians”. arXiv:1701.08213

2

Kanav Setia, Richard Chen, Julia E. Rice, Antonio Mezzacapo, Marco Pistoia, James Whitfield, “Reducing qubit requirements for quantum simulation using molecular point group symmetries”. arXiv:1910.14644