[BUG] Incorrect optimal eigenvalues (sector) in qubit tapering

[Jonas-Jaeger]

Issue: Bug Report

I found a bug in the implementation of qubit tapering [2] where the optimal eigenvalues (aka eigenvalue sector) is determined incorrectly. I assume that the implementation follows the procedure as outlined in [3]. This means that after tapering, the smalles eigenvalue that can be obtained is higher than the smallest eigenvalue of the original hamiltonian! This is regarding the implementation of tangelo.toolboxes.operators.z2_tapering.get_eigenvalues.

This bug became apparent when using the HeH+ example as in the pennylane tutorial/demo on qubit tapering [1]. Most interestingly and importantly, this exact issue is present in the pennylane implementation (as one can see in the tutorial that the smallest eigenvalue becomes higher through tapering). Therefore, I have high doubts about the correctness of this entire approach of determining the optimal sector according to Setia et al., 2020 [3]. I am interested in investigating this concern in more details, and feel free to reach out to me to join.

Expected Behavior When qubit tapering tapers qubits on which the hamiltonian has a non-trivial impact (i.e., Pauli-X instead of Identity), the eigenvalues (+1 or -1) of these qubits still have to be determined as they induce the tapered hamiltonian. An efficient implementation for this subroutine (to the best of my knowledge following the approach outlined in [3]) is part of Tangelo. It is expected that these eigenvalues are determined in a way that the tapered hamiltonian still has the same smallest eigenvalue as the original hamiltonian. In the HeH+ example, this means that the two tapered qubits should either have eigenvalues -1, +1 or +1, -1, determined by tangelo.toolboxes.operators.z2_tapering.get_eigenvalues. The expected lowest eigenvalue of the (tapered) hamiltonian would be -3.2666.

Current Behavior The current implementation of tangelo.toolboxes.operators.z2_tapering.get_eigenvalues does not seem to guarantee the computation of the correct/optimal eigenvalues. In the HeH+ example, the eigenvalues of -1, -1 are returned, which gives a lowest eigenvalue of the tapered hamiltonian of -2.8150.

Steps to Reproduce (minimal example) The minimal example including the scenario as in the Pennylane tutorial [1] can be reproduced with the following code snippet. Note that it relies on the Psi4 backend as I could not make the SecondQuantizedMolecule initialization work using PySCF.

import numpy as np
from tangelo import SecondQuantizedMolecule
from tangelo.toolboxes.molecular_computation import IntegralSolverPsi4
from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
from tangelo.toolboxes.operators.operators import qubitop_to_qubitham
from tangelo.toolboxes.operators.operators import count_qubits
from openfermion.linalg import get_sparse_operator
from tangelo.toolboxes.operators.taper_qubits import QubitTapering
from tangelo.toolboxes.operators import MultiformOperator

symbols = ["He", "H"]
geometry = np.array([[0.00000000, 0.00000000, -0.87818361],
                     [0.00000000, 0.00000000,  0.87818362]])
xyz = list(zip(symbols, map(lambda x: tuple(x.tolist()), geometry)))
print(xyz)

mol = SecondQuantizedMolecule(xyz, q=+1, spin=0, solver=IntegralSolverPsi4(), frozen_orbitals=0)
qubit_h = qubitop_to_qubitham(fermion_to_qubit_mapping(mol.fermionic_hamiltonian, "jw", up_then_down=False), mapping="jw", up_then_down=False)
qubit_h_dense = get_sparse_operator(qubit_h).todense()
print("Qubits", count_qubits(qubit_h), "Terms", qubit_h.n_terms, "Constant", qubit_h.constant)
print("Eigenvalues:", np.linalg.eigvalsh(qubit_h_dense).tolist())

# ###### 🐛 BUG: Chooses incorrect sector/eigenvalues 🐛 ######
# 
# *** This leads to a higher smallest eigenvalue than for the untapered Hamiltonain ***

print("\n\nIncorrect current implementation to select correct eigenvalues for sector:")

taper = QubitTapering(qubit_h, n_qubits=4, n_electrons=2, spin=0, up_then_down=False)

taper_h = taper.z2_tapered_op
print("Tapered Hamiltonian:", taper.z2_tapered_op)
print(count_qubits(taper_h), taper_h.n_terms, taper_h.constant)
print("Sector:", taper.z2_properties["eigenvalues"])
taper_h_dense = get_sparse_operator(taper_h).todense()
print("Eigenvalues (tapered):", np.linalg.eigvalsh(taper_h_dense).tolist())

# ###### 🐛 Verify with correct optimal sector/eigenvalues 🐛 ######
# 
# *** (Manually) Choosing the eigenvalues to be -1, 1 (also applies to 1, -1) yields the matching correct smallest eigenvalue ***

print("\n\nVerify with manually chosen correct eigenvalues for sector:")

hybrid_op = MultiformOperator.from_qubitop(qubit_h, count_qubits(qubit_h))
sector = [-1, 1]
print("Sector:", sector)
taper_h = taper.z2_taper(hybrid_op, eigenvalues=sector)
taper_h_dense = get_sparse_operator(taper_h).todense()
print("Eigenvalues (tapered):", np.linalg.eigvalsh(taper_h_dense).tolist())

Environment Please provide at least information about your OS, as well as the branch and version of Tangelo you are using. List other packages used and their version if relevant.

Possible Solution

Summary The computation of the optimal eigenvalues/sector in qubit tapering leads to sub-optimal lowest eigenvalues, i.e., the tapered hamiltonian might has a higher lowest eigenvalue as the lowest one in the un-tapered hamiltonian. This is not the expected behavior. However, since the same problem occurs in the Pennylane implementation, I highly suspect that there might be a deeper issue with the approach of determining the optimal eigenvalues/sector for the non-trivially tapered qubits as outlined in Setia et al. 2020 [3].

References

1. https://pennylane.ai/qml/demos/tutorial_qubit_tapering/
2. Bravyi, Sergey, Jay M. Gambetta, Antonio Mezzacapo, and Kristan Temme. “Tapering off Qubits to Simulate Fermionic Hamiltonians.” arXiv, January 27, 2017. https://doi.org/10.48550/arXiv.1701.08213.
3. Setia, Kanav, Richard Chen, Julia E. Rice, Antonio Mezzacapo, Marco Pistoia, and James D. Whitfield. “Reducing Qubit Requirements for Quantum Simulations Using Molecular Point Group Symmetries.” Journal of Chemical Theory and Computation 16, no. 10 (October 13, 2020): 6091–97. https://doi.org/10.1021/acs.jctc.0c00113.

[alexfleury-sb]

Hello @Jonas-Jaeger, thank you very much for the detailed explanation. I believe I understand the main point here.

I don't think we should expect the -3.2666 eigenvalue upon diagonalizing the tapered Hamiltonian with the q=+1 and spin=0 settings for the molecule. The ground state of this charged molecule should mostly be the configuration where two electrons occupy the $\sigma$ orbital. The lower eigenvalue mosltly corresponds to a configuration where two electrons occupy the $\sigma$ orbital and one electron resides in the $\sigma^*$ orbital. Typically, ions are less stable in gas phase than their neutral counterparts, and this configuration also aligns with the observed degeneracy of 2 for the -3.2666 eigenvalue.

The number of electrons (N) and the spin (Sz) are z2 symmetries, and we are in fact selecting an Hilbert space where N and Sz stays consistent with the initial Hartree-Fock state when performing qubit tapering.

As a side note, if we modify line 36 with:

taper = QubitTapering(qubit_h, n_qubits=4, n_electrons=3, spin=1, up_then_down=False)

the -3.2666 eigenvalue is kept upon diagonalization of the tapered Hamiltonian.

Did I understand your point? Does this resolve the issue? Happy to look into it further if not.

[Jonas-Jaeger]

Hi @AlexandreF-1qbit

Thank you a lot for taking your time and providing such a detailed explanation. I understand your point and hat definitely helped me out immensely. And I am happy that it is not a bug! 😄

Thanks again, Jonas