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")
#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))
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")
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
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!
Hi @Julien!
Just going through your code, I have a couple of suggestions:
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")
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
hi @whatsis
Thanks for looking at my code ! It’s still mostly based on the tutorial demo. I will bear your comment 2) in mind adjusting the learning rate when experimenting with different molecules.
I have tried your suggestions listed in 1) but these didn’t work for me.
To execute your examole above I changed
qml.BasisState(jnp.array([1, 1, 1, 1, 0, 0, 0, 0, 0, 0]), wires=range(num_wires))
to
qml.BasisState(jnp.array([1, 1, 1, 1, 0, 0, 0, 0]), wires=range(num_wires))
because my tapered Hamiltonian has 8 qubits only.
However execution then fails with the following traceback error messages
Thanks for reporting back! You might have to do this to fix this error:
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)
tap_wires = sorted(H_tapered.wires)
H_tapered = H_tapered.map_wires(dict(zip(tap_wires, range(len(tap_wires)))))
print(H_tapered)
I have opened a PR, to ensure this won’t cause any trouble in the future.
Hi @whatsis -
Yes that worked, I confirm I now get the same expectation values for the two Hamiltonians with a HF state occupancy.
However I think there is still a problem using the tapered Hamiltonian for a VQE simulation given the output below
singles, doubles = qml.qchem.excitations(n_electrons, len(H_tapered.wires))
tapered_doubles = [
qml.taper_operation(qml.DoubleExcitation, generators, paulixops, paulix_sector,
wire_order=H_tapered.wires, op_wires=double) for double in doubles
]
tapered_singles = [
qml.taper_operation(qml.SingleExcitation, generators, paulixops, paulix_sector,
wire_order=H_tapered.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)
qml.BasisState(jnp.array([1, 1, 1, 1, 0, 0, 0, 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: -4.92758152 Ha, Params: [-1.45322403e-18 0.00000000e+00 0.00000000e+00 2.74994521e-18
0.00000000e+00 2.02661812e-18 0.00000000e+00 0.00000000e+00
-1.36555428e-18 -3.06287114e-18 0.00000000e+00 0.00000000e+00
-2.16194043e-17 0.00000000e+00 -7.31294365e-18 0.00000000e+00
0.00000000e+00 -5.32089055e-18 7.77026132e-18 0.00000000e+00
2.91178160e-19 0.00000000e+00 3.02709247e-17 0.00000000e+00
-9.18048930e-18 0.00000000e+00]
n: 10, E: -4.92758152 Ha, Params: [-2.53510881e-18 0.00000000e+00 0.00000000e+00 8.29713203e-18
0.00000000e+00 4.27007667e-18 0.00000000e+00 0.00000000e+00
-2.68774048e-18 -5.18790740e-18 0.00000000e+00 0.00000000e+00
-4.17642936e-17 0.00000000e+00 -7.55417864e-18 0.00000000e+00
0.00000000e+00 -2.31510841e-18 9.64246282e-18 0.00000000e+00
-4.58477765e-19 0.00000000e+00 5.75060832e-17 0.00000000e+00
-1.89681318e-17 0.00000000e+00]
n: 15, E: -4.92758152 Ha, Params: [-3.22125980e-18 0.00000000e+00 0.00000000e+00 1.61645115e-17
0.00000000e+00 5.81964583e-18 0.00000000e+00 0.00000000e+00
-5.89643846e-18 -8.09356922e-18 0.00000000e+00 0.00000000e+00
-6.17357105e-17 0.00000000e+00 -9.37834879e-18 0.00000000e+00
0.00000000e+00 -3.10403387e-18 1.29024431e-17 0.00000000e+00
-1.67299605e-19 0.00000000e+00 8.47412418e-17 0.00000000e+00
-2.84088296e-17 0.00000000e+00]
n: 20, E: -4.92758152 Ha, Params: [-2.84760316e-18 0.00000000e+00 0.00000000e+00 2.20586430e-17
0.00000000e+00 7.41258308e-18 0.00000000e+00 0.00000000e+00
-8.99671622e-18 -1.02619736e-17 0.00000000e+00 0.00000000e+00
-8.28780658e-17 0.00000000e+00 -1.33844758e-17 0.00000000e+00
0.00000000e+00 -7.47082651e-18 2.06727044e-17 0.00000000e+00
-1.36329966e-19 0.00000000e+00 1.10848830e-16 0.00000000e+00
-3.72423742e-17 0.00000000e+00]
Hi @Julien ,
I think the issue is in how you set the BasisState within tapered_circuit(params). You have wires=H_tapered.wires while it should be wires=range(num_wires).
Let me know if this solves the issue for you!
No unfortunately I haven’t solved the issue. I think I understand what map_wires() did to H_tapered but then after updating the code to use a consecutive range(0,8) for the wires argument triggers a runtime error that suggests some tapered operations are still acting on wire 8.
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")
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')
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(6) @ Z(7), paulix_op: X(7)
generator 2: Z(8) @ Z(9), paulix_op: X(9)
generator 3: Z(0) @ Z(2) @ Z(4) @ Z(6) @ Z(8) @ Z(10), paulix_op: X(10)
generator 4: Z(1) @ Z(3) @ Z(5) @ Z(6) @ Z(8) @ Z(11), paulix_op: X(11)
paulix_sector = qml.qchem.optimal_sector(H, generators, n_electrons)
print(paulix_sector)
[1, 1, 1, 1]
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)
tap_wires = sorted(H_tapered.wires)
print (tap_wires)
print (H_tapered[0:3])
sorted_wires = range(len(tap_wires))
print (sorted_wires)
wire_map = dict(zip(tap_wires, sorted_wires))
print (wire_map)
H_tapered = H_tapered.map_wires(wire_map)
print(H_tapered[0:3])
[0, 1, 2, 3, 4, 5, 6, 8]
(-3.9776266333429193 * I(), -0.3857829818098049 * (Z(1) @ Z(3) @ Z(5) @ Z(6) @ Z(8)), -0.3857829818098048 * (Z(6) @ Z(8) @ Z(0) @ Z(2) @ Z(4)))
range(0, 8)
{0: 0, 1: 1, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6, 8: 7}
(-3.9776266333429193 * I(), -0.3857829818098049 * (Z(1) @ Z(3) @ Z(5) @ Z(6) @ Z(7)), -0.3857829818098048 * (Z(6) @ Z(7) @ Z(0) @ Z(2) @ Z(4)))
H_sparse = qml.SparseHamiltonian(H.sparse_matrix(), wires=H.wires)
num_wires = len(H_tapered.wires)
H_tapered_sparse = qml.SparseHamiltonian(H_tapered.sparse_matrix(), wires=sorted_wires)
print (H_sparse)
print("Eigenvalues of H:\n", qml.eigvals(H_sparse, k=16))
print (H_tapered_sparse)
print("\nEigenvalues of H_tapered:\n", qml.eigvals(H_tapered_sparse, k=4))
SparseHamiltonian(<Compressed Sparse Row sparse matrix of dtype 'complex128'
with 102400 stored elements and shape (4096, 4096)>, wires=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
Eigenvalues of H:
[-7.88241061 -7.80635148 -7.80635148 -7.74919247 -7.76638866 -7.76638866
-7.76638866 -7.72622282 -7.72622282 -7.72622282 -7.72622282 -7.71643788
-7.71643788 -7.71643788 -7.71643788 -7.71643788]
SparseHamiltonian(<Compressed Sparse Row sparse matrix of dtype 'complex128'
with 6464 stored elements and shape (256, 256)>, wires=[0, 1, 2, 3, 4, 5, 6, 7])
Eigenvalues of H_tapered:
[-7.88241061 -7.76638866 -7.74919247 -7.48243879]
state_tapered = qml.qchem.taper_hf(generators, paulixops, paulix_sector,
num_electrons=n_electrons, num_wires=len(sorted_wires))
print(state_tapered)
[1 1 1 1 0 0 0 0]
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")
# Tapered Hamiltonian
num_wires = len(sorted_wires)
dev = qml.device("lightning.qubit", wires=sorted_wires)
@qml.qnode(dev, interface="jax")
def circuit():
qml.BasisState(jnp.array([1, 1, 1, 1, 0, 0, 0, 0 ]), wires=sorted_wires)
return qml.state()
qubit_state = circuit()
HF_energy = qubit_state.T @ H_tapered.sparse_matrix(wire_order=sorted_wires).toarray() @ qubit_state
print(f"HF energy (tapered): {jnp.real(HF_energy):.8f} Ha")
HF energy: -7.86204208 Ha
HF energy (tapered): -7.86204208 Ha
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=sorted_wires)
@qml.qnode(dev, interface="jax")
def tapered_circuit(params):
qml.BasisState(jnp.array([1, 1, 1, 1, 0, 0, 0, 0]), wires=sorted_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}")
---------------------------------------------------------------------------
WireError Traceback (most recent call last)
Cell In[15], line 17
14 params = init_params
16 for i in range(1, 21):
---> 17 params, opt_state = update_step(i, params, opt_state)
18 energy = tapered_circuit(params)
19 if not i % 5:
Cell In[15], line 6, in update_step(i, params, opt_state)
4 def update_step(i, params, opt_state):
5 """Perform a single gradient update step"""
----> 6 grads = catalyst.grad(tapered_circuit)(params)
7 updates, opt_state = opt.update(grads, opt_state)
8 params = optax.apply_updates(params, updates)
File ~/miniforge3/envs/pennylane/lib/python3.12/site-packages/catalyst/api_extensions/differentiation.py:705, in GradCallable.__call__(self, *args, **kwargs)
703 results = jax.value_and_grad(self.fn, argnums=argnums)(*args, **kwargs)
704 else:
--> 705 results = jax.grad(self.fn, argnums=argnums)(*args, **kwargs)
706 else:
707 assert (
708 not self.grad_params.with_value
709 ), "value_and_grad cannot be used with a Jacobian"
[... skipping hidden 10 frame]
File ~/miniforge3/envs/pennylane/lib/python3.12/site-packages/pennylane/workflow/qnode.py:905, in QNode.__call__(self, *args, **kwargs)
903 if qml.capture.enabled():
904 return capture_qnode(self, *args, **kwargs)
--> 905 return self._impl_call(*args, **kwargs)
File ~/miniforge3/envs/pennylane/lib/python3.12/site-packages/pennylane/workflow/qnode.py:881, in QNode._impl_call(self, *args, **kwargs)
878 # Calculate the classical jacobians if necessary
879 self._transform_program.set_classical_component(self, args, kwargs)
--> 881 res = qml.execute(
882 (tape,),
883 device=self.device,
884 diff_method=self.diff_method,
885 interface=interface,
886 transform_program=self._transform_program,
887 gradient_kwargs=self.gradient_kwargs,
888 **self.execute_kwargs,
889 )
890 res = res[0]
892 # convert result to the interface in case the qfunc has no parameters
File ~/miniforge3/envs/pennylane/lib/python3.12/site-packages/pennylane/workflow/execution.py:227, in execute(tapes, device, diff_method, interface, transform_program, inner_transform, config, grad_on_execution, gradient_kwargs, cache, cachesize, max_diff, device_vjp, mcm_config, gradient_fn)
222 transform_program, inner_transform = _setup_transform_program(
223 transform_program, device, config, cache, cachesize
224 )
226 #### Executing the configured setup #####
--> 227 tapes, post_processing = transform_program(tapes)
229 if transform_program.is_informative:
230 return post_processing(tapes)
File ~/miniforge3/envs/pennylane/lib/python3.12/site-packages/pennylane/transforms/core/transform_program.py:580, in TransformProgram.__call__(self, tapes)
578 if argnums is not None:
579 tape.trainable_params = argnums[j]
--> 580 new_tapes, fn = transform(tape, *targs, **tkwargs)
581 execution_tapes.extend(new_tapes)
583 fns.append(fn)
File ~/miniforge3/envs/pennylane/lib/python3.12/site-packages/pennylane/devices/preprocess.py:148, in validate_device_wires(tape, wires, name)
142 raise WireError(
143 f"Cannot run circuit(s) on {name} as abstract wires are present in the device: {wires}. "
144 f"Abstract wires are not yet supported."
145 )
147 if extra_wires := set(tape.wires) - set(wires):
--> 148 raise WireError(
149 f"Cannot run circuit(s) on {name} as they contain wires "
150 f"not found on the device: {extra_wires}"
151 )
153 modified = False
154 new_ops = None
WireError: Cannot run circuit(s) on lightning.qubit as they contain wires not found on the device: {8}
Hi @Julien ,
The issue here is that you’re creating the single and double excitations based on the untapered Hamiltonian, so you’re creating excitations based on wires that don’t exist in the tapered one. You can fix this by changing the line below from len(H.wires) to len(H_tapered.wires).
singles, doubles = qml.qchem.excitations(n_electrons, len(H_tapered.wires))
In fact you had it correctly in the code you shared just before my last post (this one). You changed it back in your last code which is why it stopped working.
Yes I did make that change because I didn’t observe the behavior I was expecting when passing the tapered Hamiltonian wires. I saw instead the energy steadily increase with the number of iterations. Are you able to successfully optimise the ground state energy when running the code on your end ?
Hi @Julien , that’s probably because you’re using stochastic gradient descent as the optimizer. If you use a different optimizer, for instance opt = optax.adam(learning_rate=0.8), you’ll see that you start at a worse initial value but then get progressively better results.
