AlgorithmAdaptVQE and AlgorithmIQEB

AlgorithmAdaptVQE and AlgorithmIQEB are algorithms that construct an ansatz iteratively from a collection or pool of excitation operators. [8] [9] This differs from algorithms like AlgorithmVQE in which the ansatz is fixed, i.e. parameters of a fixed set of operators are optimized. At each ADAPT step, in AlgorithmAdaptVQE and AlgorithmIQEB, excitation operators are selected and their exponentials appended to the ansatz iteratively until convergence criteria are met. Practically, this usually results in more circuits measurements required for these algorithms compared to VQE, however the resulting ansatz will (usually) have fewer terms compared to a straight-forward application of a fixed ansatz e.g. UCCSD optimized in VQE for the same accuracy. Hence, compared to typical VQE, these algorithms trade-off number of measurements with circuit depth.

After \(N\) iterations of the algorithm, the resulting (ADAPT/IQEB) ansatz will have the form

(4)\[\hat{U}^{(N)}_{\text{ADAPT/IQEB}} = \left( e^{\theta_{N} \hat{A}_{N}} \right) \left( e^{\theta_{N-1} \hat{A}_{N-1}} \right) ... \left( e^{\theta_{\lambda} \hat{A}_{\lambda}} \right) ... \left( e^{\theta_{1} \hat{A}_{1}} \right)\]
(5)\[|\Psi^{(N)}_{\text{ADAPT/IQEB}}\rangle = \hat{U}_{\text{ADAPT/IQEB}} |\text{HF}\rangle\]

Where the excitation operators \(\hat{A}\) are chosen by the user, and the corresponding parameters \(\theta\) are optimized, by AlgorithmAdaptVQE or AlgorithmIQEB. In each iteration, the current ansatz is applied to a reference state; a Hartree-Fock state is a common choice. The two major differences between AlgorithmAdaptVQE and AlgorithmIQEB are i) the space in which the operators act on, and ii) the method to select the operators to append to the ansatz.

In the ADAPT (Adaptive Derivative-Assembled Pseudo-Trotter ansatz)-VQE algorithm [8], the operators \(\hat{A}_{\lambda}\) in the pool consist of all possible spin complemented anti-hermitian operators within unitary coupled cluster ansatz (UCC), which act in fermionic space. While these are fermionic operators, they must be transformed to qubit operators via a fermion-to-qubit mapping to be accepted by AlgorithmAdaptVQE. The operators to be appended to the ansatz at the \(n^{th}\) iteration of the algorithm are chosen by calculating the following gradient of the total energy \(E\)

(6)\[\frac{\partial E^{(n)}}{\partial \theta_{\lambda}} = \langle \Psi^{(n)}_{\text{ADAPT}} | [\hat{H}, \hat{A}_{\lambda}] | \Psi^{(n)}_{\text{ADAPT}} \rangle\]

for each excitation operator. The \(\hat{A}_{\lambda}\) which yields the largest gradient is appended to the ansatz, and a regular VQE calculation is performed to determine the optimal ansatz parameters at this iteration. In the next iteration, the gradients are recalculated as before but with the \(e^{\theta_{\lambda} \hat{A}_{\lambda}}\) from the previous iteration appended. This is repeated until all gradients are below a tolerance threshold.

The following shows an example of how to run AlgorithmAdaptVQE in which the operators of the pool are restricted to UCCSD (re-using the fermion space, qubit-encoded H2 Hamiltonian, and qubit mapping of the AlgorithmVQE example).

from inquanto.algorithms import AlgorithmAdaptVQE
from inquanto.spaces import FermionSpace
from inquanto.states import QubitState, FermionState
from inquanto.mappings import QubitMappingJordanWigner
from inquanto.ansatzes import FermionSpaceAnsatzUCCSD
from inquanto.minimizers import MinimizerScipy
from inquanto.express import get_system
from inquanto.protocols import SparseStatevectorProtocol
from pytket.extensions.qiskit import AerStateBackend

fermion_hamiltonian, fermion_space, fermion_state = get_system("h2_sto3g.h5")

jw = QubitMappingJordanWigner()
qubit_hamiltonian = jw.operator_map(fermion_hamiltonian)


space = FermionSpace(4)
fermion_state = FermionState([1, 1, 0, 0])
qubit_state = QubitState([1, 1, 0, 0])
jw_map = QubitMappingJordanWigner()

pool = space.construct_single_ucc_operators(fermion_state)
pool += space.construct_double_ucc_operators(fermion_state)
pool = jw_map.operator_map(pool)

scipy_minimizer = MinimizerScipy(method="L-BFGS-B", disp=False)

adapt = AlgorithmAdaptVQE(
    pool,
    qubit_state,
    qubit_hamiltonian,
    scipy_minimizer,
    tolerance=1.0e-3
)

protocol = SparseStatevectorProtocol(AerStateBackend())

adapt.build(
    protocol,
    protocol,
    protocol
)

adapt.run()

results = adapt.generate_report()

print("Minimum Energy: {}".format(results["final_value"]))
param_report = results["final_parameters"]
for i in range(len(param_report)):
    print(param_report[i]["symbol"], ":", param_report[i]["value"])

# TIMER BLOCK-0 BEGINS AT 2024-10-30 18:12:31.346622

# TIMER BLOCK-0 ENDS - DURATION (s):  0.2259321 [0:00:00.225932]
Minimum Energy: -1.1368465754720527
d0 : -0.10723347230091601

Here, the pool of operators was generated by the construct_single_ucc_operators() and construct_double_ucc_operators() methods of the FermionSpace class. Like the Hamiltonian, these operators were qubit-encoded before passing into AlgorithmAdaptVQE. The tolerance parameter sets the gradient threshold for convergence of the algorithm. The protocol SparseStatevectorProtocol is needed to calculate the gradients defined above.

Alternatively, one can use fermionic states and operators as arguments to the algorithm class AlgorithmFermionicAdaptVQE, which inherits from AlgorithmAdaptVQE, and performs the qubit mapping internally. Below is an example of AlgorithmFermionicAdaptVQE where again the fermionic space has been re-used (also the scipy_minimizer), and this time we use the fermionic H2 Hamiltonian (defined in the AlgorithmFermionicAdaptVQE example).

from inquanto.algorithms import AlgorithmFermionicAdaptVQE

state = FermionState([1, 1, 0, 0])

pool = space.construct_single_ucc_operators(state)
pool += space.construct_double_ucc_operators(state)

fermionic_adapt = AlgorithmFermionicAdaptVQE(
    pool,
    state,
    fermion_hamiltonian,
    scipy_minimizer,
    tolerance=1.0e-3
)

protocol = SparseStatevectorProtocol(AerStateBackend())

fermionic_adapt.build(
    protocol,
    protocol,
    protocol
)

fermionic_adapt.run()

results = fermionic_adapt.generate_report()

print("Minimum Energy: {}".format(results["final_value"]))
param_report = results["final_parameters"]
for i in range(len(param_report)):
    print(param_report[i]["symbol"], ":", param_report[i]["value"])

# TIMER BLOCK-1 BEGINS AT 2024-10-30 18:12:31.668526

# TIMER BLOCK-1 ENDS - DURATION (s):  0.2131633 [0:00:00.213163]
Minimum Energy: -1.1368465754720527
d0 : -0.10723347230091601

Note that in this case the reference state and pool have not been qubit encoded before passing into AlgorithmFermionicAdaptVQE. Jordan-Wigner encoding is performed by default.

In the IQEB (Iterative Qubit-Excitation Based)-VQE algorithm [9], the operators in the pool correspond to qubit excitations. Qubit excitation operators are generated from the following ladder operators which obey the so-called parafermionic [10] commutation relations

(7)\[\{\hat{Q}_i, \hat{Q}_i^\dagger\} = I,\]
(8)\[[\hat{Q}_i, \hat{Q}_j^\dagger] = 0 \ \ \ \ (i \neq j),\]
(9)\[[\hat{Q}_i, \hat{Q}_j] = [\hat{Q}_i^\dagger, \hat{Q}_j^\dagger] = 0 \ \ \ \ \forall i, j\]

where \(\hat{Q}_i\) (\(\hat{Q}_i^\dagger\)) is a qubit annihilation (creation) operator which changes the occupation of spin orbital \(i\) (assuming a Jordan-Wigner encoding of the Hamiltonian and reference state), and which can be represented in terms of Pauli gates

(10)\[\hat{Q}_i = \frac{1}{2} \left( X_i + iY_i \right)\]
(11)\[\hat{Q}_i^\dagger = \frac{1}{2} \left( X_i - iY_i \right)\]

The pool in AlgorithmIQEB consists of one- and two-body qubit excitation operators, built from these parafermionic operators, and acting on qubit (or spin orbital) indexes \(i, j, k, l\) (where the set of indexes is unique to the \(\lambda\)-th operator in the pool).

(12)\[\hat{A}_{\lambda}^{(ik)} = \hat{Q}_{i_\lambda}^\dagger\hat{Q}_{k_\lambda} - \hat{Q}_{k_\lambda}^\dagger\hat{Q}_{i_\lambda} = \frac{1}{2} \left( X_{i_\lambda}Y_{k_\lambda} - Y_{i_\lambda}X_{k_\lambda} \right)\]
(13)\[\hat{A}_{\lambda}^{(ijkl)} = \hat{Q}_{i_\lambda}^\dagger\hat{Q}_{j_\lambda}^\dagger\hat{Q}_{k_\lambda}\hat{Q}_{l_\lambda} - \hat{Q}_{k_\lambda}^\dagger\hat{Q}_{l_\lambda}^\dagger\hat{Q}_{i_\lambda}\hat{Q}_{j_\lambda}\]
(14)\[\begin{split}= \frac{1}{8} \left( X_{i_\lambda}Y_{j_\lambda}X_{k_\lambda}X_{l_\lambda} + Y_{i_\lambda}X_{j_\lambda}X_{k_\lambda}X_{l_\lambda} + Y_{i_\lambda}Y_{j_\lambda}Y_{k_\lambda}X_{l_\lambda} + Y_{i_\lambda}Y_{j_\lambda}X_{k_\lambda}Y_{l_\lambda} \right. \\ \left.- X_{i_\lambda}X_{j_\lambda}Y_{k_\lambda}X_{l_\lambda} - X_{i_\lambda}X_{j_\lambda}X_{k_\lambda}Y_{l_\lambda} - Y_{i_\lambda}X_{j_\lambda}Y_{k_\lambda}Y_{l_\lambda} - X_{i_\lambda}Y_{j_\lambda}Y_{k_\lambda}Y_{l_\lambda} \right)\end{split}\]

Note that the AlgorithmIQEB class inherits from AlgorithmAdaptVQE. While gradients are also used in the selection process of AlgorithmIQEB, their purpose here is to narrow down the candidate operators from the total IQEB pool. Hence the convergence of AlgorithmIQEB is not evaluated directly by gradients as in AlgorithmAdaptVQE. Instead, AlgorithmIQEB checks the total energy difference between iterations, and convergence is achieved when the decrease of energy between iterations is less than a threshold (the energy_tolerance parameter in the code block below).

The following example shows how to run AlgorithmIQEB (using the previously defined H2 qubit Hamiltonian and scipy_minimizer (see here for minimizers).

from inquanto.algorithms import AlgorithmIQEB
from inquanto.spaces import ParaFermionSpace

space = ParaFermionSpace(4)
state = QubitState([1, 1, 0, 0])

pool = space.construct_single_qubit_excitation_operators()
pool += space.construct_double_qubit_excitation_operators()

iqeb = AlgorithmIQEB(
    pool,
    state,
    qubit_hamiltonian,
    scipy_minimizer,
    n_grads=3,
    energy_tolerance=1.0e-10
)

protocol = SparseStatevectorProtocol(AerStateBackend())

iqeb.build(
    protocol,
    protocol,
    protocol,
)

iqeb.run()

results = iqeb.generate_report()

print("Minimum Energy: {}".format(results["final_value"]))
param_report = results["final_parameters"]
for i in range(len(param_report)):
    print(param_report[i]["symbol"], ":", param_report[i]["value"])

System has zero net spin -> will append spin-complementary exponents.
# TIMER BLOCK-2 BEGINS AT 2024-10-30 18:12:32.004505

# TIMER BLOCK-2 ENDS - DURATION (s):  0.2165020 [0:00:00.216502]
# TIMER BLOCK-3 BEGINS AT 2024-10-30 18:12:32.221612

# TIMER BLOCK-3 ENDS - DURATION (s):  0.0180810 [0:00:00.018081]
# TIMER BLOCK-4 BEGINS AT 2024-10-30 18:12:32.240230
# TIMER BLOCK-4 ENDS - DURATION (s):  0.0169789 [0:00:00.016979]

# TIMER BLOCK-5 BEGINS AT 2024-10-30 18:12:32.450963


# TIMER BLOCK-5 ENDS - DURATION (s):  0.0739638 [0:00:00.073964]
# TIMER BLOCK-6 BEGINS AT 2024-10-30 18:12:32.525579
# TIMER BLOCK-6 ENDS - DURATION (s):  0.0531882 [0:00:00.053188]
# TIMER BLOCK-7 BEGINS AT 2024-10-30 18:12:32.579392
# TIMER BLOCK-7 ENDS - DURATION (s):  0.0531284 [0:00:00.053128]

CONVERGED!!!
Final ansatz elements after 2 iteration(s):
r_1_1     [(0.125j, X0 Y1 X2 X3), (0.125j, Y0 X1 X2 X3), (0.125j, Y0 Y1 Y2 X3), (0.125j, Y0 Y1 X2 Y3), (-0.125j, X0 X1 Y2 X3), (-0.125j, X0 X1 X2 Y3), (-0.125j, X0 Y1 Y2 Y3), (-0.125j, Y0 X1 Y2 Y3)] 

Minimum Energy: -1.1368465754720527
r_1_1 : -0.10723347230091601

Notice that six separate VQE calculations have been performed (one for each TIMER BLOCK in the output log). This is due to two reasons. i) Our choice of n_grads=3, which tells AlgorithmIQEB that we want to narrow down the pool to those terms which have the three largest gradients, and a VQE calculation for each term will be run. Of these three, the term which has the largest effect on the energy will be appended to the ansatz in this iteration. ii) Since AlgorithmIQEB establishes convergence by comparing the energy difference between iterations, a second iteration is performed, again with three separate terms (appended to the previously found term). In this case, convergence is found at the second iteration, which means the resulting ansatz will have the form of the first iteration. Note that the internal VQE initial parameter coefficients are set to be all zeros.

As in the case of AlgorithmAdaptVQE, we define a pool of operators. However here we employ the ParaFermionSpace class to handle the parafermionic operator algebra. The operators in the IQEB pool consist of all unique permutations of qubit indexes for one- and two-body terms, obtained by the construct_single_qubit_excitation_operators() and construct_double_qubit_excitation_operators() methods of ParaFermionSpace() (which do not need a reference state). This results in an asymptotically larger pool than AlgorithmAdaptVQE [9]. However, the advantage of AlgorithmIQEB is that excitation operators act directly in parafermionic space, hence the strings of Pauli-Z operators resulting from Jordan-Wigner encoding, in order to maintain fermionic exchange symmetry, are not required. Therefore each qubit excitation of AlgorithmIQEB acts on a fixed number of qubits, independent of the system size.