Variational Quantum Deflation for excited states

In this file, we will introduce the methodology for finding electronic excited state energies using the Variational Quantum Deflation (VQD) approach in InQuanto. For the original publication by Higgot, Wang and Brierley (2019) and more technical details, please go here.

In VQE, one finds a minimum of the energy by classically optimizing the function below with respect to the wavefunction parameters, \(\{\lambda\}\):

\(E(\lambda) = \langle \psi(\lambda) | H | \psi(\lambda)\rangle = \sum_j c_j \langle \psi (\lambda) | P_j | \psi (\lambda) \rangle\),

where \(P_j\) are terms in the qubit Hamiltonian, and \(c_j\) are classically pre-computed coefficients. In contrast, when executing a VQD simulation, the objective function is modified to include a penalty term, multiplied by a weight \(\beta\),

\(F(\lambda_k) = \langle \psi(\lambda_k) | H | \psi(\lambda_k)\rangle + \sum_{i=0}^{k-1} \beta_i |\langle \psi (\lambda_k) | P_j | \psi (\lambda_i) \rangle|^2\).

This enforces the requirement that each \(|\psi(\lambda_k)\rangle\) be orthogonal to the other \(|\psi(\lambda_0)\rangle, \ldots |\psi(\lambda_{k-1})\rangle\) found by previous optimizations of the objective function, \(F(\lambda_k)\).

In order to run a VQD calculation, one needs several things:

  1. A molecular, electronic Hamiltonian operator, \(H\),

  2. A mapping object for constructing the qubit Hamiltonian,

  3. An ansatz for computing the ground state energy, \(|\psi(\lambda)\rangle\),

  4. A classical minimizer

  5. A complete VQE experiment,

  6. A second ansatz for describing the excited states, \(\{|\psi(\lambda_k)\rangle\}\),

  7. A expression for evaluating the weights, \(\{\beta_i\}\),

  8. The number of excited states of interest, \(k\).

We are concerned with first finding the ground state energy of the second-quantized molecular electronic Hamiltonian using the Variational Quantum Eigensolver (please see tutorials 1, and 2 for a more in-depth explanation). So let’s proceed with points 1-4 in the list above.

First, we need to import a backend, and the appropriate space and state objects. Since we are looking at fermions, these are inquanto.spaces.FermionSpace and inquanto.states.FermionState. We will use the AerBackend available through the pytket qiskit extension to simulate the quantum hardware.

[1]:

from pytket.extensions.qiskit import AerBackend
from inquanto.spaces import FermionSpace
from inquanto.states import FermionState

We’re going to simulate the dihydrogen molecule in the STO-3G basis. There are 4 spin orbitals – two of which are occupied – and the reference state lives in 4-dimensional Fock space. The corresponding objects can be constructed as below:

[2]:

from inquanto.express import load_h5
from inquanto.mappings import QubitMappingJordanWigner
h2 = load_h5("h2_sto3g.h5")
fermion_hamiltonian = h2["hamiltonian_operator"]
qubit_hamiltonian = QubitMappingJordanWigner().operator_map(fermion_hamiltonian)

space = FermionSpace(4)
state = FermionState([1, 1, 0, 0])

print(fermion_hamiltonian.df())

    Coefficients            Terms
0       0.743018
1      -1.270293          F0^ F0
2       0.340031  F1^ F0^ F0  F1
3       0.340031  F1^ F0^ F0  F1
4       0.089834  F1^ F0^ F2  F3
5       0.089834  F1^ F0^ F2  F3
6      -1.270293          F1^ F1
7      -0.089834  F2^ F0^ F0  F2
8      -0.089834  F2^ F0^ F0  F2
9       0.334289  F2^ F0^ F0  F2
10      0.334289  F2^ F0^ F0  F2
11     -0.089834  F2^ F1^ F0  F3
12     -0.089834  F2^ F1^ F0  F3
13      0.334289  F2^ F1^ F1  F2
14      0.334289  F2^ F1^ F1  F2
15     -0.456807          F2^ F2
16      0.334289  F3^ F0^ F0  F3
17      0.334289  F3^ F0^ F0  F3
18     -0.089834  F3^ F0^ F1  F2
19     -0.089834  F3^ F0^ F1  F2
20     -0.089834  F3^ F1^ F1  F3
21     -0.089834  F3^ F1^ F1  F3
22      0.334289  F3^ F1^ F1  F3
23      0.334289  F3^ F1^ F1  F3
24      0.089834  F3^ F2^ F0  F1
25      0.089834  F3^ F2^ F0  F1
26      0.351407  F3^ F2^ F2  F3
27      0.351407  F3^ F2^ F2  F3
28     -0.456807          F3^ F3

where we have loaded in the molecular Hamiltonian using InQuanto’s express module, and mapped it to a qubit Hamiltonian using the Jordan-Wigner transformation.

Now, for points 3 and 4 in the list, we need to construct an ansatz for our ground state calculation and choose a classical minimizer. For this example we will use the k-UpCCGSD inquanto.ansatzes.FermionSpaceAnsatzkUpCCGSD ansatz and the COBYLA minimizer available through the inquanto.minimizers.MinimizerScipy object.

[3]:

from inquanto.ansatzes import FermionSpaceAnsatzkUpCCGSD
from inquanto.minimizers import MinimizerScipy

ansatz = FermionSpaceAnsatzkUpCCGSD(space, state, k_input=2)
minimizer = MinimizerScipy(method="COBYLA")

We’re now in a position to address item 5 in the list - running a complete VQE calculation. We know from the first equation in this notebook that the objective function is the expectation value of the qubit Hamiltonian.

[4]:

from inquanto.computables import ExpectationValue
from inquanto.algorithms import AlgorithmVQE

expectation_value = ExpectationValue(ansatz, qubit_hamiltonian)
vqe = AlgorithmVQE(
    objective_expression=expectation_value,
    minimizer=minimizer,
    initial_parameters=ansatz.state_symbols.construct_random(seed=0)
)

We have passed in some random \(\{\lambda\}\) using the state_symbols attribute of the ansatz object as our starting guess.

We now choose our measurement protocol – in this case, a direct measurement by operator averaging, so we choose inquanto.protocols.PauliAveraging – and the number of shots we wish to simulate in each iteration. Then, we build the algorithm object and execute.

[5]:

from inquanto.protocols import PauliAveraging

vqe.build(protocol_objective=PauliAveraging(AerBackend(), shots_per_circuit=10000))

vqe.run()

# VQE Energy:     -1.1354204303965678
# VQE Parameters: [ 1.301 -1.538  0.287  1.699 -1.339 -0.344]
print("VQE Energy:    ", vqe.final_value)
print("VQE Parameters:", vqe.final_parameters.to_array())

# TIMER BLOCK-0 BEGINS AT 2024-07-19 12:51:15.909454
# TIMER BLOCK-0 ENDS - DURATION (s): 156.5863703 [0:02:36.586370]
VQE Energy:     -1.127936169552505
VQE Parameters: [ 1.396 -1.573  0.407  1.687 -1.386 -0.431]

According to point 6, we now need a second, deflationary ansatz we can use to describe the excited states. To do this, we’ll use the same ansatz structure and just make a copy of the first ansatz. We update our symbols from those used in the ground state to some other symbols. Consider this as constructing the symbols in \(\{\lambda_k\}\) using those in \(\{\lambda\}\) as a reference.

[6]:

ansatz_2 = ansatz.subs("{}_2")

It is almost time to construct, build and execute our VQD algorithm. First, we need to write expressions corresponding to the terms in the functional

\(F(\lambda_k) = \langle \psi(\lambda_k) | H | \psi(\lambda_k)\rangle + \sum_{i=0}^{k-1} \beta_i |\langle \psi (\lambda_k) | P_j | \psi (\lambda_i) \rangle|^2\).

We will refer to the leading term as expectation_value, the weights as weight_expression, and the overlap term in the penalty as overlap_expression.

We will also select the weights as the expectation value of the deflationary ansatz with respect to the sign-flipped qubit Hamiltonian to ensure it is sufficiently large to act as a constraint rather than a weak penalty.

[7]:

from inquanto.computables import OverlapSquared

expectation_value = ExpectationValue(ansatz_2, qubit_hamiltonian)

weight_expression = ExpectationValue(ansatz_2, -1 * qubit_hamiltonian)

overlap_expression = OverlapSquared(ansatz, ansatz_2)

As was the case with the previous VQE experiment, we must choose protocols for measuring the overlaps. For this we choose to use the vacuum test, available through the inquanto.protocols.ComputeUncompute object.

We must also choose the number of excited states we wish to find. For this calculation, we’ll choose to find 3 and pass this into the inquanto.algorithms.AlgorithmVQD constructor in the n_vectors argument.

[8]:

from inquanto.algorithms import AlgorithmVQD
from inquanto.protocols import ComputeUncompute
# instantiate VQD object
vqd = AlgorithmVQD(
    objective_expression=expectation_value,
    overlap_expression=overlap_expression,
    weight_expression=weight_expression,
    minimizer=MinimizerScipy(method="COBYLA"),
    initial_parameters=ansatz_2.state_symbols.construct_random(seed=0),
    vqe_value=vqe.final_value,
    vqe_parameters=vqe.final_parameters,
    n_vectors=3,
)
# build object
backend=AerBackend()
vqd.build(
    #small number of shots for demonstration purposes leads to large stochastic error
    objective_protocol=PauliAveraging(backend, shots_per_circuit=1000),
    weight_protocol=PauliAveraging(backend, shots_per_circuit=1000),
    overlap_protocol=ComputeUncompute(backend, n_shots=1000),
)

# execute
vqd.run()

# print results
# VQD excited state energies:    [-1.1354204303965678, -0.4949467734066755, -0.11981345040527547, 0.5478565335949009]
print("VQD excited state energies:   ", vqd.final_values)

# TIMER BLOCK-1 BEGINS AT 2024-07-19 12:53:52.571887
# TIMER BLOCK-1 ENDS - DURATION (s): 79.0145929 [0:01:19.014593]
# TIMER BLOCK-2 BEGINS AT 2024-07-19 12:55:11.586617
# TIMER BLOCK-2 ENDS - DURATION (s): 142.2693380 [0:02:22.269338]
# TIMER BLOCK-3 BEGINS AT 2024-07-19 12:57:33.856439
# TIMER BLOCK-3 ENDS - DURATION (s): 147.8532623 [0:02:27.853262]
VQD excited state energies:    [-1.127936169552505, -0.47952989815415065, -0.16244947039556473, 0.5461336969835955]