Advanced VQE: Quantum Subspace Expansion¶

References: https://arxiv.org/abs/1603.05681, https://arxiv.org/abs/1707.06408, https://arxiv.org/abs/1807.10050

The VQE algorithm exhibits a "natural" robustness against errors, especially regarding $\vec{\theta}^\star$, the optimal value of the parameter. Unfortunately, the energy evalutation (i.e. mean-value measurement) can still suffer from important errors.

McClean et al. drew inspiration from Linear Response Theory to design an extension to the VQE, the Quantum Subspace Expansion (QSE). The core idea is to expand the Hamiltonian post-VQE on a well-chosen subspace (i.e. where an improved, lower, energy lies) and solve classically the associated generalized eigenvalue problem with the hope of getting an improved value for the ground state energy.

More precisely, the QSE can be split into different steps:

  1. Choice of qubit operators;
  2. Expansion of the Hamiltonian on the subspace defined by the two previous choices; Construction of the overlap matrix;
  3. Resolution of the generalized eigenvalue problem.

Thus, the $n$-qubit QSE using $G$ as the chosen set of $n$-qubit operators, is associated with the following state subspace: $$ \{ \hat{\sigma}|\psi^\star\rangle, \qquad \hat{\sigma} \in G \} $$ where $|\psi^\star\rangle = |\mathrm{UCC}(\vec{\theta}^\star)\rangle$ is the output of the VQE.

The expanded Hamiltonian and overlap matrices, $(H_{i, j})$ and $(S_{i, j})$, are then measured via a quantum computer, i.e. $$ H_{i, j} = \langle \psi^\star | \hat{\sigma}_i^\dagger \hat{H} \hat{\sigma}_j | \psi^\star\rangle \qquad S_{i, j} = \langle \psi^\star | \hat{\sigma}_i^\dagger \hat{\sigma}_j | \psi^\star\rangle $$

Finally, the associated generalized eigenvalue problem is solved classically and the minimal solution is extracted, i.e. $$ E_{\mathrm{QSE}} = \min\{E, \qquad H \vec{x} = E S \vec{x}\} $$

Part 1: Problem definition and UCC preparation¶

In [1]:
import numpy as np

n_electrons = 2
one_body_integrals = np.array([[-1.25246357, 0], [0, -0.475948715]])
two_body_integrals = np.array(
    [
        [[[0.674488766, 0], [0, 0.181288808]], [[0, 0.181288808], [0.663468096, 0]]],
        [[[0, 0.663468096], [0.181288808, 0]], [[0.181288808, 0], [0, 0.697393767]]],
    ]
)
orbital_energies = np.array([-0.57797481, 0.66969867])
nuclear_repulsion = 0.7137539936876182

## The natural-orbital occupation numbers (NOONs) are computed from 1-RDM (computed in CISD here)
noons = np.array([1.9745399697399246, 0.025460030260075376])

Note :¶

If you have installed the pySCF module, you can use the following lines of code to perform the quantum-chemistry part of the computation

import numpy as np
from qat.fermion.chemistry.pyscf_tools import perform_pyscf_computation

geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))]
basis = 'sto-3g'

rdm1, orbital_energies, nuclear_repulsion,\
nels, one_body_integrals, two_body_integrals, info = perform_pyscf_computation(geometry=geometry, basis=basis, spin=0, charge=0, verbose=True)

# Get NOONs from 1-RDM (computed in CISD)
noons = list(reversed(sorted(np.linalg.eigvalsh(rdm1))))

UCC preparation¶

In [2]:
from qat.fermion.chemistry import MolecularHamiltonian, MoleculeInfo
from qat.fermion.chemistry.ucc import (
    guess_init_params,
    get_hf_ket,
    get_cluster_ops,
)

# Wrap the hamiltonian data into the `MolecularHamiltonian` class.
mol_h = MolecularHamiltonian(one_body_integrals, two_body_integrals, nuclear_repulsion)

# Use the MoleculeInfo class to then restrict the active space
molecule = MoleculeInfo(hamiltonian=mol_h, n_electrons=n_electrons, noons=noons, orbital_energies=orbital_energies)

# Computation of the initial parameters
theta_init = guess_init_params(
    molecule.two_body_integrals,
    molecule.n_electrons,
    molecule.orbital_energies,
)

# Define the initial Hartree-Fock state
ket_hf_init = get_hf_ket(molecule.n_electrons, nqbits=molecule.nqbits)

# Compute the cluster operators
cluster_ops = get_cluster_ops(molecule.n_electrons, noons=molecule.noons)

# Get the ElectronicStructureHamiltonian
H_active = molecule.hamiltonian.get_electronic_hamiltonian()

Transformation to qubit space¶

In [3]:
from qat.fermion.transforms import transform_to_parity_basis, get_parity_code, recode_integer

transformation, code = transform_to_parity_basis, get_parity_code

H_active_sp = transformation(H_active)
nqbits = H_active_sp.nbqbits

# Express the cluster operator in spin terms
cluster_ops_sp = [transformation(t_o) for t_o in cluster_ops]

# Encoding the initial state to new encoding
hf_init_sp = recode_integer(ket_hf_init, code(nqbits))

# Finally: build_uccsd
from qat.fermion.chemistry.ucc import construct_ucc_ansatz

prog = construct_ucc_ansatz(cluster_ops_sp, hf_init_sp)
circ = prog.to_circ()

circ.display()
No description has been provided for this image

Part 2: VQE optimization¶

Let us first compute the exact ground state energy:

In [4]:
eigvals = np.linalg.eigvalsh(H_active_sp.get_matrix())
E_min = min(eigvals)
print("E_min (exact diagonalization) = %s" % (E_min))

E_min (exact diagonalization) = -1.137270167926503

Minimizing with COBYLA¶

In [5]:
from qat.qpus import get_default_qpu
from qat.plugins import ScipyMinimizePlugin

qpu = get_default_qpu()

optimizer_scipy = ScipyMinimizePlugin(method="COBYLA", tol=1e-3, options={"maxiter": 1000}, x0=theta_init)
stack = optimizer_scipy | qpu
res = stack.submit(circ.to_job(job_type="OBS", observable=H_active_sp))
theta_VQE = res.meta_data["parameters"]
print("Optimal theta (VQE): %s" % theta_VQE)
print("E (VQE) = %s (err = %s %%)" % (res.value, 100 * abs((res.value - E_min) / E_min)))

Optimal theta (VQE): [np.float64(-0.0002372411065788997), np.float64(0.0001072258641688005), np.float64(0.11248152980500957)]
E (VQE) = -1.137269567465197 (err = 5.279847506291602e-05 %)

Minimizing with SPSA¶

Let us try with an alternative optimizer:

In [6]:
from qat.plugins import SPSAMinimizePlugin

optimizer_spsa = ScipyMinimizePlugin(x0=theta_init)
stack = optimizer_spsa | qpu

res = stack.submit(circ.to_job(job_type="OBS", observable=H_active_sp))
theta_VQE = res.meta_data["parameters"]
print("Optimal theta (VQE): %s" % theta_VQE)
print("E (VQE) = %s (err = %s %%)" % (res.value, 100 * abs((res.value - E_min) / E_min)))
e_vqe = res.value

Optimal theta (VQE): [np.float64(-4.590068224901082e-09), np.float64(-4.590068224901082e-09), np.float64(0.11306820868992833)]
E (VQE) = -1.1372701679264878 (err = 1.327655781425361e-12 %)

Part 3: Quantum subspace expansion¶

In [7]:
from qat.fermion import SpinHamiltonian
from qat.core import Term

expansion_operators = [
    SpinHamiltonian(nqbits, [], 1.0),
    SpinHamiltonian(nqbits, [Term(1.0, "ZZ", [0, 1])]),
]

from qat.fermion.chemistry.qse import apply_quantum_subspace_expansion

# we use the optimal parameters found by VQE
opt_circ = circ.bind_variables(eval(res.meta_data["parameter_map"]))

e_qse = apply_quantum_subspace_expansion(H_active_sp, opt_circ, expansion_operators, qpu, return_matrices=False)
print("E(QSE) = %s (err = %s %%)" % (e_qse, 100 * abs((e_qse - E_min) / E_min)))

E(QSE) = -1.1372701679264972 (err = 5.076330928979322e-13 %)