Quantum Subspace Expansion AlgorithmQSE

The Quantum Subspace Expansion (QSE) obtains extra correlation and energetics of excited states, on top of the VQE ground state \(|\Psi_{\mathrm{VQE}}\rangle\) by doing a linear excitation expansion as [15]

(16)\[|\Psi_{\mathrm{QSE}} \rangle = \sum_{i > j} c_{ij} |\psi_{ij}\rangle\]


(17)\[|\psi_{ij}\rangle = a^\dagger_i a_j |\Psi_{\mathrm{VQE}} \rangle\]

and \(c_{ij}\) are the complex coefficients.

If \(|\Psi_{\mathrm{VQE}}\rangle\) is a correlated state, the basis functions of the subspace generated by applying the excitation operators are not orthogonal in general. The optimal solution for the ground and excited states within this subspace can be found by solving the generalized eigenvalue problem expressed as

(18)\[HC = SCE\]

where \(H\) is the Hamiltonian matrix and \(S\) is the overlap matrix in this subspace. \(E\) is the diagonal matrix of eigenvalues and \(C\) is the matrix of eigenvectors. The matrix elements of \(H\) and S are evaluated by quantum computers as

(19)\[H_{ij,kl} = \langle \psi_{ij} | \hat{H} | \psi_{kl} \rangle\]
(20)\[S_{ij,kl} = \langle \psi_{ij} | \psi_{kl} \rangle\]

And the eigenstates \(C\) and \(E\) are obtained with classical diagonalization.

AlgorithmQSE is a high level interface of InQuanto to run the QSE algorithm for given qubit operator and ansatz with optimal parameters to evaluate the eigenstates. We take the minimal basis hydrogen molecule as an example to calculate the spin-adapted excited states from the VQE state with UCCSD ansatz.

The first step is constructing the molecular Hamiltonian and the Fock state in fermion objects.

from inquanto.express import get_system
hamiltonian, fermion_space, fermion_state = get_system("h2_sto3g.h5")

And then map them into the qubit objects. Here we use the UCCSD ansatz and then address the excited states with QSE, but one can use a cheaper ansatz for the VQE and obtain the extra correlation at the later QSE step, through the linear expansion and matrix diagonalization.

from inquanto.mappings import QubitMappingJordanWigner
from inquanto.ansatzes import FermionSpaceAnsatzUCCSD

jw = QubitMappingJordanWigner()
qubit_hamiltonian = jw.operator_map(hamiltonian)
ansatz = FermionSpaceAnsatzUCCSD(fermion_space, fermion_state, qubit_mapping=jw)

Subsequently, we need to define the subspace in which to solve the generalized eigenvalue equations; namely the set of expansion operators to be applied to our ground state.

expansion_operators = jw.operator_map(

Note that we use singlet single excitation operators \(E_{ij} = \sum_{\sigma}^{\mathrm{spin}}a_{i\sigma}^{\dagger}a_{j\sigma}\) to span the spin-adapted subspace. We then define a new Computable object for the \(H\) and \(S\) matrices in Eq. (18), using the QSEMatricesComputable class. This computable and the values of the parameters of the ansatz are then passed to AlgorithmQSE. In general we first need to perform a VQE calculation to obtain the optimized values of the parameters of the ansatz. In this specific example, we already know them from the previous VQE example (see in AlgorithmVQE), so we provide the values as an array.

from inquanto.computables.composite import QSEMatricesComputable
from inquanto.algorithms import AlgorithmQSE

vqe_parameters = ansatz.state_symbols.construct_from_array(
    [0.0, 0.0, -0.107233493519281]

computable = QSEMatricesComputable(

algorithm = AlgorithmQSE(

The following procedure is basically the same as the examples above, but we use a shot-based protocol. No solver is required, as it is a non-variational method.

from inquanto.protocols import PauliAveraging
from pytket.extensions.qiskit import AerBackend

protocol_expression = PauliAveraging(AerBackend(), shots_per_circuit=2000)


# TIMER Measure and calculate H and S matrix element BEGINS AT 2024-06-21 13:05:23.557277
# TIMER Measure and calculate H and S matrix element ENDS - DURATION (s):  1.1041073 [0:00:01.104107]
# TIMER Solve HC=CSE BEGINS AT 2024-06-21 13:05:24.661504
# TIMER Solve HC=CSE ENDS - DURATION (s):  0.0006954 [0:00:00.000695]
{'final_value': array([-1.335, -1.103, -0.128]), 'final_states': array([[-0.427+0.116j, -0.012-0.479j,  0.008-0.002j],
       [-8.422-0.176j, -2.037-0.018j,  0.055+0.009j],
       [ 0.017-0.126j,  0.009-0.064j, -0.616+0.343j],
       [-0.369-4.439j, -0.089-1.07j ,  0.033+0.015j]])}

The eigenvalues and eigenvectors are stored as final_value and the final_states, respectively. Numerical instabilities in the matrix diagonalization are handled by removing near linear dependencies in the basis.