Qubit tapering / LiH

Hi,

I have succesfully reproduced the qubit tapering demo at Qubit tapering | PennyLane Demos

using the provided example HeH+

However when I replaced HeH+ by LiH

I noticed a large difference in HF energy between the tapered Hamiltonian and the reference Hamiltonian

The subsequent VQE calculation converges to E ca. -7.120 Ha which is far from the exact ground state energies

Am I doing something wrong when tapering the reference state?

Hi @Julien,

Can you please copy your full code to see if I can replicate it? You should be able to paste it as code here. Make sure to include the code for creating H_tapered.

Hi @CatalinaAlbornoz - here is a mardkown export of my run without the output cells to stick under word count limit for posts on the forum

import pennylane as qml
from jax import numpy as jnp
import jax
jax.config.update('jax_enable_x64', True)
jax.config.update("jax_platform_name", "cpu")

Tapering the molecular Hamiltonian

#symbols = ['He','H']
#geometry = jnp.array([[0.0000, 0.0000, 0.0],
#                      [0.0000, 0.0000, 1.40]])
#n_electrons = 2
#charge = 1
symbols = ['Li','H']
geometry = jnp.array([[0.0000, 0.0000, -1.5065],
                      [0.0000, 0.0000, 1.5065]])
n_electrons = 4
charge = 0
molecule = qml.qchem.Molecule(symbols, geometry, charge=charge)
H, qubits = qml.qchem.molecular_hamiltonian(molecule, mapping='jordan_wigner')#, active_electrons=2)# active_electrons=2, active_orbitals=3, mapping='parity')
H
print (qubits, len(H))
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]}")

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

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)

print (len(H_tapered))

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))

Tapering the reference state

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

dev = qml.device("default.qubit", wires=H.wires)
@qml.qnode(dev, interface="jax")
def circuit():
    qml.BasisState(jnp.array([1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 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, 1, 1, 0, 0, 0, 0]), 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")

VQE simulation

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)

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, 101):
    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}")

Hi @Julien ,

I was able to run your code and here are my conclusions.

The eigenvalues for the tapered Hamiltonian look ok.

Eigenvalues of H:
 [-7.88241061 -7.80635148 -7.80635148 -7.74919247 -7.76638866 -7.76638866
 -7.76638866 -7.72622282 -7.72622282 -7.71643788 -7.72622282 -7.72622282
 -7.71643788 -7.71643788 -7.71643788 -7.71643788]

Eigenvalues of H_tapered:
 [-7.88241061 -7.76638866 -7.74919247 -7.48243879]

However when I print H_tapered_sparse I see that the wires are in the wrong order SparseHamiltonian(<256x256 sparse matrix of type '<class 'numpy.complex128'>' with 6464 stored elements in Compressed Sparse Row format>, wires=[1, 3, 5, 6, 8, 0, 2, 4])

I’m thinking this might be the cause for the discrepancy.

The way I found to fix this was to match the basis state with the wire order, when using the tapered Hamiltonian.

So instead of having

qml.BasisState(jnp.array([1, 1, 1, 1, 0, 0, 0, 0]), wires=H_tapered.wires)

You can write:

qml.BasisState(jnp.array([1, 1, 0, 0, 0, 1, 1, 0]), wires=H_tapered.wires)

This is definitely confusing and not ideal, so I’ll forward it to the team. But for now you have a solution!

I’ve created this bug report for this issue. Let us know if you encounter further issues with this or other features.

Hi @CatalinaAlbornoz - that’s great thank you ! I wouldn’t have figured this out.
I confirm that if I update the BasisState in my code I get the same energies for the reference HF and tapered HF states.
However the VQE simulation still doesn’t work. The energy keeps going up with the code below.
Could it be that the code in qml.taper_operations to generate tapered_doubles and tapered_singles excitations also needs updating ?

Here is a markdown transcript of the VQE simulation section of my script

VQE simulation

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):
    # see https://github.com/PennyLaneAI/pennylane/issues/6934
    # qml.BasisState(state_tapered, wires=H_tapered.wires)
    qml.BasisState(jnp.array([1, 1, 0, 0, 0, 1, 1, 0]), wires=H_tapered.wires)
    for idx, tapered_op in enumerate(tapered_doubles + tapered_singles):
        tapered_op(params[idx])
    return qml.expval(H_tapered)
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, 21):
    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: -7.87935652 Ha, Params: [ 2.46492422e-02  0.00000000e+00  0.00000000e+00  6.55901686e-03
  0.00000000e+00  0.00000000e+00 -6.55667685e-03  1.31656184e-02
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  1.31631049e-02  0.00000000e+00  0.00000000e+00  3.37953434e-02
  0.00000000e+00  0.00000000e+00 -9.71598204e-04  0.00000000e+00
  0.00000000e+00  0.00000000e+00 -4.15302663e-04  0.00000000e+00
  0.00000000e+00 -7.23695329e-04  0.00000000e+00  0.00000000e+00
  1.69429837e-03  5.70135731e-03  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  5.70126746e-03  0.00000000e+00
  0.00000000e+00 -8.33430563e-04  4.06219740e-04  0.00000000e+00
  0.00000000e+00  1.69762442e-03  0.00000000e+00  0.00000000e+00
 -7.25831981e-04 -5.70444802e-03  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00 -5.70441375e-03  0.00000000e+00
  0.00000000e+00  8.32574746e-04  0.00000000e+00  0.00000000e+00
 -9.70812452e-04  0.00000000e+00  0.00000000e+00  0.00000000e+00
  3.48040293e-02  0.00000000e+00  0.00000000e+00  8.21948600e-02
  0.00000000e+00  0.00000000e+00 -8.22174936e-02  4.72213244e-02
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  4.72284006e-02  0.00000000e+00  0.00000000e+00  2.18852544e-01
  1.68591837e-05  0.00000000e+00  0.00000000e+00 -1.66550919e-05
  1.93953635e-05  0.00000000e+00  0.00000000e+00  7.08757148e-05
 -6.32481222e-03  0.00000000e+00  0.00000000e+00  1.57798858e-03
  6.23143187e-03  0.00000000e+00  0.00000000e+00 -1.82676403e-03]
n: 10, E: -7.67830945 Ha, Params: [-0.03118367  0.          0.         -0.04921943  0.          0.
  0.04916761 -0.0250059   0.          0.          0.          0.
 -0.02500058  0.          0.         -0.34165131  0.          0.
 -0.00085892  0.          0.          0.          0.0013031   0.
  0.          0.00130729  0.          0.         -0.00047911  0.00678703
  0.          0.          0.          0.          0.00678736  0.
  0.          0.01277258 -0.00135846  0.          0.         -0.00047649
  0.          0.          0.00134594 -0.00680488  0.          0.
  0.          0.         -0.00680526  0.          0.         -0.01272762
  0.          0.         -0.00083867  0.          0.          0.
  0.04829753  0.          0.          0.10051709  0.          0.
 -0.10052185  0.05293173  0.          0.          0.          0.
  0.05294715  0.          0.          0.23077611 -0.00327594  0.
  0.         -0.00724606  0.00354432  0.          0.          0.0072207
 -0.0202453   0.          0.          0.00416127  0.01988543  0.
  0.         -0.00448381]
n: 15, E: -5.74740593 Ha, Params: [ 3.42301895e-02  0.00000000e+00  0.00000000e+00  9.39111696e-02
  0.00000000e+00  0.00000000e+00 -9.32016095e-02  2.35288130e-02
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  2.35105474e-02  0.00000000e+00  0.00000000e+00  1.21263304e+00
  0.00000000e+00  0.00000000e+00 -1.90042338e-03  0.00000000e+00
  0.00000000e+00  0.00000000e+00 -4.50195872e-04  0.00000000e+00
  0.00000000e+00 -4.53322932e-03  0.00000000e+00  0.00000000e+00
  6.58253330e-03  5.58413456e-03  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  5.58413946e-03  0.00000000e+00
  0.00000000e+00 -3.83959985e-02  5.38354613e-04  0.00000000e+00
  0.00000000e+00  6.62677421e-03  0.00000000e+00  0.00000000e+00
 -4.73063180e-03 -5.70589006e-03  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00 -5.70595640e-03  0.00000000e+00
  0.00000000e+00  3.80898796e-02  0.00000000e+00  0.00000000e+00
 -2.04658878e-03  0.00000000e+00  0.00000000e+00  0.00000000e+00
  5.61029128e-02  0.00000000e+00  0.00000000e+00  1.07079555e-01
  0.00000000e+00  0.00000000e+00 -1.07251698e-01  6.78489032e-02
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  6.78754348e-02  0.00000000e+00  0.00000000e+00  2.34990594e-01
  4.48768232e-03  0.00000000e+00  0.00000000e+00  4.26932623e-02
 -3.20485551e-04  0.00000000e+00  0.00000000e+00 -4.34416401e-02
 -6.21571705e-02  0.00000000e+00  0.00000000e+00  1.30569851e-03
  6.10879286e-02  0.00000000e+00  0.00000000e+00 -5.30038488e-03]
n: 20, E: -5.78581997 Ha, Params: [-9.46279277e-03  0.00000000e+00  0.00000000e+00 -7.80011958e-02
  0.00000000e+00  0.00000000e+00  7.66902936e-02 -7.92613509e-03
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
 -7.91298433e-03  0.00000000e+00  0.00000000e+00 -1.20082698e+00
  0.00000000e+00  0.00000000e+00  5.20635391e-04  0.00000000e+00
  0.00000000e+00  0.00000000e+00  1.40394832e-03  0.00000000e+00
  0.00000000e+00  2.89277508e-03  0.00000000e+00  0.00000000e+00
 -3.63781755e-03  6.16988482e-03  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  6.17058162e-03  0.00000000e+00
  0.00000000e+00  4.17530200e-02 -1.62449347e-03  0.00000000e+00
  0.00000000e+00 -3.73815158e-03  0.00000000e+00  0.00000000e+00
  3.27053673e-03 -6.73505752e-03  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00 -6.73584356e-03  0.00000000e+00
  0.00000000e+00 -4.19311000e-02  0.00000000e+00  0.00000000e+00
  7.00433749e-04  0.00000000e+00  0.00000000e+00  0.00000000e+00
  5.76136626e-02  0.00000000e+00  0.00000000e+00  1.15437532e-01
  0.00000000e+00  0.00000000e+00 -1.14992514e-01  8.34750191e-02
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  8.35272181e-02  0.00000000e+00  0.00000000e+00  2.25481894e-01
 -7.89427205e-04  0.00000000e+00  0.00000000e+00  3.35269298e-02
  6.15510769e-03  0.00000000e+00  0.00000000e+00 -3.39722792e-02
 -1.39341141e-01  0.00000000e+00  0.00000000e+00  6.60354817e-03
  1.35634038e-01  0.00000000e+00  0.00000000e+00 -1.06372547e-02]


I think you’re right @Julien , you may need to reorder some things here too.

I’m not sure however if the reordering needs to happen at the level of the operations, excitations, or the actual state.

One of my colleagues is taking a look at this so we’ll get back to you on this by early next week.

Thanks again for helping us discover this bug!

1 Like

Hi @Julien!

Just going through your code, I have a couple of suggestions:

  1. First, try defining your device and HF state in the usual sorted order and also providing wire_order when you compute the sparse_matrix method to adjust for the different wire order via H_tapered.wires:
num_wires = len(H_tapered.wires)
dev = qml.device("lightning.qubit", wires=num_wires) # device wires are [0, 1, ..., num_wires-1] 

@qml.qnode(dev, interface="jax")
def circuit():
    qml.BasisState(jnp.array([1, 1, 1, 1, 0, 0, 0, 0, 0, 0]), wires=range(num_wires))
    return qml.state()
qubit_state = circuit()

HF_energy = qubit_state.T @ H_tapered.sparse_matrix(wire_order=range(num_wires)).toarray() @ qubit_state
print(f"HF energy (tapered): {jnp.real(HF_energy):.8f} Ha")
  1. When you optimize, implement the circuit as usual once again. However, do remember to tweak your learning rate according to your problem at hand. Larger molecules might have more complex energy landscapes to optimize over, leading to need for smaller learning rates to not overshoot during subsequent iterations.
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
]
state_tapered = qml.qchem.taper_hf(generators, paulixops, paulix_sector,
                                   num_electrons=n_electrons, num_wires=len(H.wires))

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

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

opt = optax.sgd(learning_rate=0.1) # adjust the learning rate