Variational quantum deflation

What it fills

• Estimate an excited-state energy

  Given a Hamiltonian you can query and a statement of which state above the lowest one is wanted, return an estimate of that state's energy — or of its distance from the ground state, which is the quantity an experiment usually measures. Whether the answer comes back as an absolute energy or as an excitation energy is a property of the method, and the two are not interchangeable.

When it applies

Higgott et al. put the mechanism in one sentence — "overlap estimation can be used to deflate eigenstates once they are found, enabling the calculation of excited state energies and their degeneracies" — and price it against the ground-state run it reuses: an implementation "requires the same number of qubits as VQE and at most twice the circuit depth". They also state what survives noise: the method "is robust to control errors, is compatible with error-mitigation strategies and can be implemented on near-term quantum computers". The ordering is the binding condition and it is structural rather than a limitation of any implementation: deflation is defined against the states already found, so the k-th needs all k below it, and an error in one of those is inherited by every state above it.

Requires

These do not move the route along. The method needs each of them alongside its own work, and the cost of getting them is part of what the method costs.

• Estimate a Hamiltonian's ground-state energy

  Given a Hamiltonian you can query and some way of preparing trial states, return an estimate of its lowest eigenvalue to a stated additive error. The state achieving that energy may or may not come back with the number; almost every application wants the number.

  VQD's zeroth stage is plain VQE: the deflation sum in $F(\lambda_k)$ is empty at $k=0$, so $\lambda_0$ is fixed by minimising $E(\lambda)=\langle\psi(\lambda)|H|\psi(\lambda)\rangle=\sum_j c_j\langle\psi(\lambda)|P_j|\psi(\lambda)\rangle$ (Eq. 1). What this route takes from the slot is not only $E_0$ but the argument $\lambda_0$: the overlap against the ground state is later evaluated as $|\langle 0|R(\lambda_0)^\dagger R(\lambda_k)|0\rangle|^2$ (in general $R(\lambda_i)^\dagger$ for $i<k$), so the classical parameters must be retained. That inverse is not merely gate-reversed but refit — fixing the optimal $\lambda_0^*$ and varying $\lambda_0$ to maximise $|\langle 0|R(\lambda_0)^\dagger R(\lambda_0^*)|0\rangle|^2$ — which is what keeps VQE's robustness to control errors. VQE on $-H$ re-enters the same slot for $E_{d-1}$ ($d$ eigenvectors), giving $\Delta:=E_{d-1}-E_0\ge E_k-E_i$. assumption: this stage returns the exact ground state; at $|\langle\tilde{\psi}_0|\psi_0\rangle|^2=1-\varepsilon_0$ the paper bounds only the next state, $\varepsilon_1'=O(\varepsilon_0\beta_0+2\varepsilon_1\|H\|)$, and only for $\beta_0>\frac{E_1-E_0}{1-\varepsilon_0}$; accumulation is measured negligible for $k<4$ and substantial at $k=4,5$

  assumption

Example

given  H = sum_j c_j P_j on L Pauli terms; a parameterised ansatz family
       with state-preparation circuit R(lambda) such that R(lambda)|0> =
       |psi(lambda)>; the index k of the eigenstate wanted (k=0 is the
       ground state, d the total number of eigenvectors of H)

requires  states 0..k-1 already found by this SAME procedure -- their
    optimal parameters lambda_0, ..., lambda_(k-1) and circuits
    R(lambda_0), ..., R(lambda_(k-1)) retained. Deflation is defined
    against them, so evaluating the objective below needs
    lambda_0..lambda_(k-1) already known, and an iterative procedure is
    required to reach the k-th eigenvalue                        (Sec. 2)

# index convention, matching the paper: j indexes Pauli terms, i indexes an
#   already-found state being deflated against, m the state being found now

# --- stage 0: the ground state is plain VQE, the deflation sum is empty ---
lambda_0 = argmin_lambda  E(lambda)
    where  E(lambda) := <psi(lambda)| H |psi(lambda)>
                       = sum_j c_j <psi(lambda)| P_j |psi(lambda)>    (Eq. 1)

# --- stage m = 1 .. k: one full optimisation pass per state, in order -----
for m = 1 .. k:

    # weight selection -- the condition is on a GAP, not an arbitrary size
    choose beta_0, ..., beta_(m-1)
        # beta_i > E_m - E_i suffices to put F's minimum at E_m rather than
        #   at the lower state i -- but only "if the ansatz is sufficiently
        #   powerful"; expressiveness is a proviso, not a guarantee (Sec. 3)
        # an always-usable choice when H is given as a Pauli sum:
        #   beta_i = 2 * sum_j |c_j|, via Delta := E_(d-1) - E_0 >= E_m - E_i
        #   together with Delta <= 2*||H|| <= 2 * sum_j |c_j|         (Sec. 3)
        #   -- at (1+m)^2 times VQE's own sample bound                (App. A)
        # self-correcting fallback: an under-set beta_i = gamma - E_i makes
        #   the minimum land at F = gamma, which is how the mistake is
        #   caught; doubling gamma until a strictly lower energy is found
        #   costs O(log(E_m - E_0)) further runs of the algorithm     (Sec. 3)

    lambda_m = argmin_lambda  F(lambda)
        where  F(lambda) := <psi(lambda)| H |psi(lambda)>
                    + sum_{i=0}^{m-1} beta_i * overlap(lambda, i)      (Eq. 2)
            # equivalently: the ground state of the effective Hamiltonian
            #   H_m := H + sum_{i<m} beta_i |i><i|, i.e. VQE's own landscape
            #   with the first m levels lifted by beta_i               (Eq. 4)
            # the |i> are ASSUMED to be true eigenstates of H, with possibly
            #   non-distinct energies                     (Sec. 3, footnote 1)
        # each overlap term below is bought with real shots on a real
        #   circuit -- it is not an abstract inner product

# --- overlap: read off a circuit, no ancilla ------------------------------
def overlap(lambda_trial, i):                                        (Sec. 4)
    prepare  R(lambda_i)^dagger R(lambda_trial) |0>
        # trial preparation circuit, THEN the inverse of state i's own
        #   preparation circuit, on the same register -- that inverse is
        #   the one built by refit_inverse(i) below
    measure in the computational basis, O(1 / eps^2) times
    return  fraction of all-zero bitstrings
        # = |<0| R(lambda_i)^dagger R(lambda_trial) |0>|^2
        #     = |<psi(lambda_i)|psi(lambda_trial)>|^2
    # alternative: the destructive SWAP test -- almost VQE's own circuit
    #   depth but twice the qubits; 2N qubits, depth O(1) (Sec. 4; App. B)
    # given greater gate depth, alpha-QPE reduces the total RUNTIME of
    #   overlap estimation from O(1/eps^2) to O(1/eps); the paper states
    #   this as a runtime, not as a shot count                       (Sec. 4)

# --- the inverse circuit is REFIT, not gate-reversed ----------------------
def refit_inverse(i):                                                (Sec. 4)
    fix lambda_i*  # the optimum already found for state i
    vary trial lambda_i to maximise |<0| R(lambda_i)^dagger R(lambda_i*) |0>|^2
    # inverting the gates of the original decomposition is known in theory,
    #   but device errors can make that implementation inaccurate in
    #   practice; refitting is what lets VQD keep VQE's own robustness to
    #   control errors

return  lambda_0, ..., lambda_k and their circuits R(lambda_0), ..., R(lambda_k)
    # each E_m is then read off exactly as in stage 0, via the same
    #   Pauli-term measurements VQE already needs                    (Eq. 1)

# sample cost of one objective evaluation at stage k, states 0..k-1 ASSUMED
#   perfectly prepared -- the loose main-text bound:
#     M^(k) <= (1/eps^2) * ( sum_j |c_j| + (1/2) sum_i beta_i )^2      (Eq. 3)
#   (Appendix A carries a tighter, variance-weighted form; this is not it)
# same number of qubits as VQE, at most twice the circuit depth for this
#   (non-SWAP-test) implementation                         (Sec. 4; Sec. 8)

# NOT part of this method: the ansatz itself. The paper's own H2/STO-3G run
#   used a generalised UCC ansatz (UCCGSD), |psi> = exp(T - T^dagger)|HF>,
#   with T1, T2 indexing ANY orbital pair/quadruple regardless of occupancy
#   in the reference (unlike standard UCCSD) -- 6 and 3 independent
#   parameters for T1, T2 -- but the paper states many other ansatz
#   circuits can be used instead                  (Sec. 5; App. C; Sec. 8)
# NOT part of this method: symmetry. Electron-number conservation falls out
#   of T automatically; spin or other symmetries instead need SQUARED
#   penalty terms
#   F_C(lambda_k) = F(lambda_k)
#       + sum_i mu_i * ( <psi(lambda_k)|C_i|psi(lambda_k)> - c_i )^2
#   on symmetry operators C_i (e.g. N_e, S^2, S_z), outside the core
#   deflation objective above                             (App. E, Eq. 28)
# NOT accounted for above: imperfect states 0..k-1. A fidelity error eps_0
#   in state 0 forces beta_0 > (E_1 - E_0)/(1 - eps_0) even to bound state
#   1, and the resulting energy error stays linear in eps_0 -- measured
#   negligible for k < 4 and substantially higher at k = 4, 5 in the
#   paper's own H2 simulation                              (Sec. 6; App. D)

Cost, as the source states it

Higgott, Wang and Brierley: the $k$-th eigenstate ($k=0$ the ground state) takes $k+1$ sequential optimisation runs, because deflation adds $k$ overlap terms $\beta_i|\langle\psi(\lambda_k)|\psi(\lambda_i)\rangle|^2$ to VQE's $L=O(N^4)$ Pauli terms ($N$ spin orbitals) and needs $\lambda_0,\dots,\lambda_{k-1}$ first — plus $O(\log(E_k-E_0))$ further runs if a large enough weight must be found by doubling. One objective evaluation costs $M^{(k)}\le\varepsilon^{-2}(\sum_j|c_j|+\frac{1}{2}\sum_i\beta_i)^2$ samples (Eq. 3), assuming states $0,\dots,k-1$ are prepared perfectly, with $\varepsilon$ one standard deviation (Eq. 8 defines $\varepsilon^2$ as the variance), no confidence attached. At the always-safe $\beta_i=2\sum_j|c_j|$ that is $(1+k)^2$ times VQE's bound — a bound-to-bound ratio the authors flag as such and expect to be far smaller in practice. Each overlap takes $O(1/\varepsilon^2)$ computational-basis measurements; $\alpha$-QPE cuts overlap-estimation runtime, not shot count, to $O(1/\varepsilon)$ at greater depth. Same qubits as VQE, at most twice the circuit depth. The $\beta$ guarantee holds only for a sufficiently expressive ansatz, and nothing bounds the optimiser iterations.

Implementations

• Numerical simulation of VQD on $\mathrm{H}_2$ in the STO-3G basis

  • Variational Quantum Computation of Excited States

  Higgott, Wang and Brierley open on the application — "The calculation of excited state energies of electronic structure Hamiltonians has many important applications, such as the calculation of optical spectra and reaction rates" — and on the gap: "While low-depth quantum algorithms, such as the variational quantum eigenvalue solver (VQE), have been used to determine ground state energies, methods for calculating excited states currently involve the implementation of high-depth controlled-unitaries or a large number of additional samples." Their own runs answer two different questions about the deflation objective they construct. Section 5 asks whether it converges to the right levels at all, by sweeping a molecule whose spectrum is known exactly: "We simulated VQD on H2 in the STO-3G basis for a range of internuclear separations and compared it to exact diagonalisation, as shown in Fig. 2." Section 6 asks the question the ordering forces — whether the error in state $i$ is inherited by every state above it — and answers it by simulation rather than by the Appendix D bound alone: "We analysed this accumulation of errors further through numerical simulations of VQD in the presence of sampling error, shown in Fig. 3." Degeneracy is what the authors call "another key advantage of our method", and they name the alternatives it separates them from: "the folded spectrum and WAVES methods rely on the energies of states to differentiate between them and have no systematic way of determining the degeneracy of the eigenvalues."

  The Fig. 2 sweep sets "$\beta_i = 3$ Ha for all $i$" and uses a generalised rather than the standard coupled-cluster ansatz: "We instead use a generalised unitary coupled cluster ansatz (UCCGSD) with $|\psi_0\rangle = |\mathrm{HF}\rangle$ set to the Hartree-Fock state", the cluster operators re-indexed so that "$p, q, r, s$ can now index any orbital (irrespective of its occupation in the reference state)". Appendix C then counts what the fermionic anticommutation relations leave free, for an instance of four spin orbitals — the examples it gives index orbitals 0 to 3: "it can be directly verified that there are only 6 and 3 independent parameters for $T_1$ ... and $T_2$ ... respectively". Optimisation is classical and derivative-free — "A tolerance of $10^{-2}$ was used with a Nelder-Mead optimiser (xatol=fatol=$10^{-2}$, as implemented in the scipy Python scientific library), and the best of two consecutive (randomly initialised) runs was used for each bond length and energy level" — with no warm start: "While we initialised the UCCGSD parameters randomly in this work". The Fig. 3 study is a different design on the same algorithm: "We analysed three different sampling rates: $M = 10^6$, $10^7$ and $10^8$ samples per Hamiltonian subterm and overlap term, running 225 simulations of VQD (with random initial parameters) for each of these three scenarios." Runs that came out in the wrong order were dropped, and the fraction is reported rather than buried: "Of these runs, $\sim 20\%$ of the simulations found the eigenstates in the incorrect order and were discarded for consistency in the analysis", leaving "the remaining $\sim 180$ runs for each state $k$". A control arm was run beside it — "we also simulated 130 runs (dashed lines) using 'exact' states $i < k$ in $H_k$ ($< 10^{-7}$ energy error in each state $i < k$)" — which is what makes the accumulated error separable from the sampling error. Error bars are "$1\sigma$ standard errors for the median estimates (calculated using bootstrap resampling)".

  Molecular hydrogen in a minimal basis, twice over. Fig. 2 shows "All ground and excited state energy levels of H2 in the STO-3G basis, calculated using exact diagonalisation (blue dotted line) and our variational quantum deflation (VQD) method (red filled circles) over a range of internuclear separations"; Fig. 3 fixes the geometry, reporting "each energy level $k$ in the spectrum of H2 at bond distance (0.7414 Å) in the STO-3G basis". Neither Section 5 nor Section 6 writes the Hamiltonian down. The form the paper does write is Eq. (27), introduced in Appendix E as an illustration rather than as the operator simulated — "For example, consider the electronic structure Hamiltonian in second quantised form" — $H = \sum_{ij} h_{ij} a_i^\dagger a_j + \sum_{ijkl} h_{ijkl} a_i^\dagger a_j^\dagger a_k a_l$, "where the coefficients $h_{ij}$ and $h_{ijkl}$ denote the one- and two-electron integrals, respectively"; no provenance is stated for the integrals behind the H2 runs. Electron number is pinned by the ansatz rather than by a penalty term: the "fermionic unitary coupled cluster ansatz we use in Section 5 conserves the desired number of electrons ($\eta = 2$) of neutral molecular Hydrogen for all input parameters".

  What the paper names is the stack the Fig. 2 simulation ran on: "The results in Fig. 2 were simulated using ProjectQ and FermiLib", citing [25, 37] — McClean et al., "OpenFermion: the electronic structure package for quantum computers", and Steiger, Häner and Troyer, "ProjectQ: an open source software framework for quantum computing" (Quantum 2, 49 (2018)). ProjectQ is a Python framework at https://github.com/ProjectQ-Framework/ProjectQ; FermiLib is a Python library at https://github.com/ProjectQ-Framework/FermiLib, described by its own README as "designed as a library on top of ProjectQ". The classical optimiser is scipy's Nelder-Mead, taken "as implemented in the scipy Python scientific library".

  Every number here comes from classical simulation; the paper reports no run of its own on a quantum device. Against exact diagonalisation of the same Hamiltonian, "the median error of our method relative to exact diagonalisation is less than $4\times10^{-6}$ Ha for all energy levels, significantly better than chemical accuracy of $1.6 \times 10^{-3}$ Ha", and the degenerate levels are resolved rather than sidestepped: "Our method finds all 6 eigenstates systematically, including all those in the 3-dimensional degenerate subspace spanned by the 1st, 2nd and 3rd excited states." Those two figures are the Fig. 2 sweep, and the simulator behind Fig. 2 is ProjectQ with FermiLib. The sampling study that follows is where the ordering condition shows its price, and it shows it as a threshold rather than a slope — but the paper names no software at all for it, so the machine behind every number in the rest of this section is not stated. Below the fourth level the inheritance is invisible: "the median errors for all states $k < 4$ are very similar (for a given $M$) to the errors when using an exact $H_k$, and are all below chemical accuracy, demonstrating that error accumulation is negligible for these states". At and above it, it is not: "For $k = 4$ and $k = 5$ the accumulated error is substantially higher than the error using an exact $H_k$, however, showing that VQD is most effective for low-lying states", and the extra cost is quoted in samples — "Achieving chemical accuracy for $k = 5$ requires $10^7$ samples, instead of $10^6$ for $0 < k < 4$."

• The VQD eigensolver in Qiskit Algorithms

  A released library implementation rather than a study. The module header states its subject and its source in two lines — "The Variational Quantum Deflation Algorithm for computing higher energy states. See https://arxiv.org/abs/1805.08138." — and the class docstring restates the method: "VQD is a quantum algorithm that uses a variational technique to find the k lowest eigenvalues of the Hamiltonian $H$ of a given system", which it computes "by optimizing over a modified cost function where each successive eigenvalue is calculated iteratively by introducing an overlap term with all the previously computed eigenstates that must be minimised, thus ensuring higher energy eigenstates are found". The file carries "(C) Copyright IBM 2022, 2026". The package README carries a warning block above its installation instructions: "Qiskit Algorithms is no longer officially supported by IBM".

  Energies come from a Qiskit Estimator primitive (`BaseEstimatorV2`) and the overlaps from a `BaseStateFidelity` the caller supplies. The fidelity used throughout the shipped tests is `ComputeUncompute`, whose docstring defines it as "the state fidelity of two quantum circuits following the compute-uncompute method", giving $|\langle\psi(x)|\phi(y)\rangle|^2$, and cites Havlíček et al., "Supervised learning with quantum-enhanced feature spaces" — the same reference Higgott et al. give for the low-depth overlap $|\langle 0|R(\lambda_i)^\dagger R(\lambda_k)|0\rangle|^2$ of their Section 4, rather than for the destructive SWAP test of their Appendix B. The deflation weights are the paper's own, exposed as `betas` and documented as "Beta parameters in the VQD paper". Left unset they are computed, not guessed: `upper_bound = sum(np.abs(operator.coeffs))`, then `betas = np.asarray([upper_bound * 10] * self.k)` — ten times $\sum_j |c_j|$, where the always-safe choice in the paper is $\beta_i = 2\sum_j |c_j|$. In 0.4.0 the docstring beside that argument describes a different rule, "a default value computed as the mean square sum of the coefficients of the observable", which is not what those two lines compute. The class also carries a guard of its own, `convergence_threshold`, documented as "the maximal average fidelity an eigenstate is allowed to have with the previous eigenstates". It is checked once a step's optimisation has returned, against the beta-weighted mean fidelity at that step's converged optimum — `current_optimal_point["total_fidelity"][0] / (step - 1)` — and if the threshold is exceeded it raises an `AlgorithmError`, so the remaining eigenvalues are never computed.

  None of its own: the Hamiltonian is whatever the caller passes. The shipped test module supplies a two-qubit $\mathrm{H}_2$ operator as a literal, `SparsePauliOp.from_list([("II", -1.052373245772859), ("IZ", 0.39793742484318045), ("ZI", -0.39793742484318045), ("ZZ", -0.01128010425623538), ("XX", 0.18093119978423156)])`, together with the reference levels it is checked against, `self.h2_energy_excited = [-1.85727503, -1.24458455, -0.88272215, -0.22491125]`. The file gives no source for either set of numbers.

  `qiskit_algorithms.eigensolvers.VQD`, Python, in `qiskit_algorithms/eigensolvers/vqd.py` at https://github.com/qiskit-community/qiskit-algorithms, under the Apache License 2.0 — "This code is licensed under the Apache License, Version 2.0". It is distributed on PyPI as `qiskit-algorithms` ("Qiskit Algorithms: A library of quantum computing algorithms"), version 0.4.0 at the time of reading, with its tests in `test/eigensolvers/test_vqd.py` in the same repository.

  What shows that it runs is the repository's own test suite, and the machine is a statevector simulator. `test/eigensolvers/test_vqd.py` builds the algorithm on `StatevectorEstimator(seed=self.seed)` with `ComputeUncompute(StatevectorSampler(seed=self.seed, default_shots=10_000))` at `self.seed = 50`, a `COBYLA` optimiser and a two-qubit ansatz `n_local(2, rotation_blocks=["ry", "rz"], entanglement_blocks="cz", reps=1)`. `test_full_spectrum` asks for the whole spectrum of the two-qubit $\mathrm{H}_2$ operator with `k=4, betas=[3, 3, 3]` and asserts the returned eigenvalues against `[-1.85727503, -1.24458455, -0.88272215, -0.22491125]` to `decimal=2`; the basic test asks for the lowest two and asserts to `decimal=1`. The transpiler path is exercised against `GenericBackendV2(num_qubits=3, coupling_map=[[0, 1], [1, 2]], seed=54)` — a fake backend from `qiskit.providers.fake_provider`, not a device. The test module names no quantum hardware.

What it needs

1. Choose a parameterised trial state 13 methods

   Fix the gate structure of a circuit family and leave its rotation angles open. What comes back is not a circuit but the set of states the later optimisation is allowed to search — which is why this is a slot of its own and not a paragraph in one method's write-up.

2. Minimise the objective over the parameters 6 methods

   Search the parameters of a fixed circuit family for the ones that minimise a chosen objective, given that every evaluation of that objective is a noisy estimate someone paid shots for. The output is a routine that prepares one particular state — the family collapsed to a member.

3. Estimate an observable 4 methods

   Given the ability to prepare $|ψ⟩$ and a description of an observable $O$, return a classical scalar within $\varepsilon$ of $⟨O⟩$ at confidence $1−δ$. The state is never returned; only the number is.

Other ways to fill the same slot

Different approaches

• Subspace-search variational eigensolver

  Send several mutually orthogonal input states through one parameterised circuit and minimise their energies together. A unitary keeps orthogonal inputs orthogonal, so the whole low-energy subspace comes out of a single optimisation — no earlier state to deflate, and no ancilla to test overlaps with.

• Quantum subspace expansion

  Stop optimising and start measuring: take the state the ground-state run already produced, measure matrix elements in a small space of operators applied to it, and let a classical solver diagonalise the little matrix that results. The extra accuracy is bought with classical work and more measurements, not with more coherence.

• Quantum equation of motion

  Ask for the gaps directly rather than for two energies to subtract. Measure the matrix elements the classical equation-of-motion formalism needs between excitation operators on the ground state, and solve its generalised eigenvalue problem classically; what comes back is an excitation energy, which is the quantity a spectrum is made of.

• Folded-spectrum variational eigensolver

  Point the same search somewhere other than the bottom. Minimising the variance around a chosen energy makes every eigenstate a minimum and the one nearest that energy the reachable one — so a state can be found without knowing its index, and the bill arrives as a squared Hamiltonian with far more terms to measure.

• Penalty-constrained variational eigensolver

  Add a term to the objective that punishes the trial state for leaving the symmetry sector you asked for, and the ordinary ground-state search returns that sector's lowest state — an excited state of the whole Hamiltonian whenever the sector is not the one the ground state lives in. Which penalty is used matters: one common form is proved not to work.

• Multistate contracted variational eigensolver

  Optimise one circuit for several states at once and read off the transitions between them — both the energy of each transition and the oscillator strength that says how strongly light drives it. The answer is a spectrum with intensities, which is what an absorption experiment actually produces.

In the Atlas

• Variational quantum deflation

  Overlap penalties against previously found states turn excited-state search into a sequence of VQE objectives.

Sources

• Variational Quantum Computation of Excited States↗Oscar Higgott, Daochen Wang, Stephen Brierley · 2018theory yes · simulation yes · hardware no