Quantum subspace expansion

• 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.

McClean et al. state the shape of the guarantee rather than a bound: variational approaches "fit in a more general hierarchy of measurement and classical computation that allows one to obtain increasingly accurate solutions with additional classical resources", and they "demonstrate numerically on a sample electronic system that this method both allows for the accurate determination of excited electronic states as well as reduces the impact of decoherence, without using any additional quantum coherence time or formal error correction codes". The decoherence claim is separate and is argued from an "exactly solvable channel model of variational state preparation". What the abstract does NOT give is the mechanism at the level a reader would need to implement it — which operators span the expansion, and which matrix elements have to be measured — and neither is invented here; both sit in the body of the paper, and the honest summary of the abstract is the trade it names, not a construction it does not state.

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.

  The route enters this slot with an ordinary VQE loop. A parameterised $|\Psi(\vec{\theta})\rangle$ is prepared and its parameters driven to the Rayleigh-Ritz minimum $\min_{\vec{\theta}} \langle H \rangle(\vec{\theta}) = \langle \Psi(\vec{\theta})|H|\Psi(\vec{\theta})\rangle / \langle \Psi(\vec{\theta})|\Psi(\vec{\theta})\rangle$, the energy read out by repeated state preparation and partial tomography as the contraction $\langle H \rangle = \sum_{ik} h_{ik}\,{}^1D^i_k + \sum_{ijkl} h_{ijkl}\,{}^2D^{ij}_{lk}$. McClean et al. model the preparation as a channel rather than a unitary: minimising $\mathrm{Tr}[(\sum_i K_i |\Psi\rangle\langle\Psi| K_i^\dagger) H]$ is equivalent to solving $H'|\Psi\rangle = E|\Psi\rangle$ on the transformed $H' = \sum_i K_i^\dagger H K_i$ for the lowest eigenvalue and eigenvector pair. assumption: the model prepares an arbitrary pure state, so it fixes the optimal possible performance separately from ansatz choice or experimental protocol Symmetry can be imposed by a penalty $H \to H + \sum_i \lambda_i (O_i - o_i)^2$ approximation: exact only as $\lambda_i \to \infty$, and only if the parameterisation can satisfy the symmetry at all; the authors call it prohibitively expensive and warn it may not generalise to excited states of a different symmetry. More than the energy leaves this slot: the 1- and 2-RDM already measured are the reference the next hop expands about, which then buys the 3- and 4-RDM with a polynomial number of additional measurements.

  approximationassumption

• Estimate an observable

  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.

  Here the route stops optimising. It expands about $\rho$, the state at the exit of the channel, in the vectors $a_i^\dagger a_j |\Psi\rangle$. Partial tomography supplies the elements within that expansion, $O^{ij}_{kl} = \mathrm{Tr}[(a_i^\dagger a_j)^\dagger O\, a_k^\dagger a_l\, \rho]$ with metric $S^{ij}_{kl} = \delta_{ik}\,{}^1D^j_l - 2\,{}^2D^{jk}_{li}$, and those coupling it to the reference, $O^{ij}_g = \mathrm{Tr}[(a_i^\dagger a_j)^\dagger O \rho]$ with $S^{ij}_g = {}^1D^j_i$; that $g$ row and column border the matrices, so the reference sits inside the diagonalised space. A one-body operator reaches at most the 3-RDM, a two-body $V = \sum_{pqrs} a_p^\dagger a_q^\dagger a_r a_s$ the 4-RDM through a $-24\,{}^4D^{jkpq}_{srli}$ term. The classical step solves $H_{\mathrm{LR}} C = S_{\mathrm{LR}} C E$ for the ground and excited states; symmetry is imposed here with no $\lambda_i \to \infty$ limit, by building $O_{\mathrm{LR}}$ in the same subspace and projecting onto its eigenvectors of the wanted eigenvalue. approximation: truncating the hierarchy $\mathcal{B}^k_f$ at $k=1$ spans only the linear-response space, justified as the dominant contribution in a linear response theory of local time-dependent perturbations approximation: the commutator form $H^{ij}_{kl} = \langle \Psi|(a_i^\dagger a_j)^\dagger [H, a_k^\dagger a_l]|\Psi\rangle + E_g \langle \Psi|(a_i^\dagger a_j)^\dagger a_k^\dagger a_l|\Psi\rangle$ removes the 4-RDM at $O(M^6)$ measured terms and is exact only when $|\Psi\rangle$ is exact; ZC further sets the 3-particle cumulant ${}^3\Delta = 0$ approximation: ZA instead starts from the un-reduced 4-RDM and reconstructs it from the 2-RDM alone, setting ${}^4\Delta = {}^3\Delta = 0$

  approximation

given  a normalized reference state |Psi> already produced by a prior VQE run
       that minimized <H>(theta) over the ansatz parameters (Eq. 18) -- this
       method does not run that optimization, it consumes its output
       the reference's 1- and 2-RDM, 1D^i_k and 2D^ij_kl (Eq. 3-4), already
       measured as part of that VQE run
       M spin-orbitals (JW-encoded qubits), so excitation indices i,j range
       over [1, M]
       optionally, an operator O (a symmetry operator, or any observable
       whose value is wanted in the excited states) expressed as a sum of
       one- and two-body fermion operators, as H itself is (Eq. 10)

requires  |Psi> is the single state the subspace is built ABOUT, not a family
    of trial states -- the paper demonstrates the k=1 truncation's accuracy
    from an exact ground-state reference (Fig. 4) and separately from an
    approximate, channel-degraded one (Fig. 5, amplitude/phase damping); it
    states no general bound on truncation error for an arbitrary reference

# --- build the k=1 excitation subspace -------------------------------------
Bf1 = { a_i^dagger a_j |Psi>  :  i, j in [1, M] }                (Eq. 6, Eq. 8 at k=1)
    # justified only as the dominant contribution in a linear-response theory
    #   of local time-dependent perturbations to the system -- that result is
    #   cited to a separate reference, not derived in this paper       (ref [28])
    # Bf1 is level k=1 of a hierarchy Bf^k; k=Ne would span the whole
    #   Ne-fermion space and make the classical step exact diagonalization,
    #   but its cost then scales exponentially in system size      (text after Eq. 8)

# --- measure the subspace's Hamiltonian and overlap by partial tomography --
for each (i,j), (k,l) in [1,M] x [1,M]:
    H^ij_kl = <Psi| (a_i^dagger a_j)^dagger H a_k^dagger a_l |Psi>          (Eq. 45, O = H)
    S^ij_kl = delta_ik * 1D^j_l  -  2 * 2D^jk_li                            (Eq. 47)
    # S is the overlap/metric of the (non-orthogonal) excitation vectors --
    #   it is NOT the identity, which is exactly why the next step is a
    #   GENERALIZED eigenvalue problem, not a plain diagonalization
    # S alone is fixed by the 1- and 2-RDM (Eq. 47); H is NOT -- a one-body
    #   term of H reaches the 3-RDM, a two-body term the 4-RDM (matrix
    #   elements given explicitly, Eq. 48-51), and determining these costs a
    #   polynomial number of measurements beyond what the ground-state VQE
    #   run already took                                        (text after Eq. 5)

# --- solve the subspace problem classically ---------------------------------
solve  H_LR C = S_LR C E   for eigenpairs (C, E)                    (Eq. 7 / Eq. 21)
    # E's lowest eigenvalue approximates the ground state within Bf1; its
    #   higher eigenvalues are the excited-state estimates this method is for

# --- optional: enforce an exact symmetry in the subspace --------------------
if a symmetry operator O with target eigenvalue o is supplied:
    build  O^ij_kl  the same way as H^ij_kl above
    project onto the eigenvectors of O_LR carrying eigenvalue o
    diagonalize H_LR restricted to that projected subspace          (Supp. Sec. C)
    # exact for any finite subspace here -- unlike the penalty
    #   H -> H + sum_i lambda_i (O_i - o_i)^2 used to impose a symmetry on
    #   the reference state itself, which the paper notes is only exact as
    #   lambda_i -> infinity and "may not generalize well to excited states
    #   if they are of a different symmetry than the ground state"   (main text, Eq. 20)
    # the paper gives this projection only in prose, not as an explicit
    #   projector formula                                            (Supp. Sec. C)

# --- optional variants that avoid measuring the 4-RDM -----------------------
# the fermionic QSE above is exact within Bf1 but, in the paper's own words,
#   "can still be prohibitively expensive in the number of measurements
#   required"; both variants below are approximations aimed at that   (Supp. Sec. I)
# ZC: H^ij_kl ~= <Psi|(a_i^dagger a_j)^dagger [H, a_k^dagger a_l]|Psi>
#                + Eg <Psi|(a_i^dagger a_j)^dagger a_k^dagger a_l|Psi>        (Eq. 59)
#     exact only when |Psi> is the exact ground state (Eg its eigenvalue);
#     the commutator drops the expression to at most the 3-RDM, O(M^6)
#     terms; setting the irreducible 3-particle cumulant 3Delta = 0 removes
#     even that, needing nothing beyond the ORIGINAL ground-state 2-RDM
# ZA: skips the commutator and reconstructs the 4-RDM from the 2-RDM alone,
#     assuming 4Delta = 3Delta = 0
# the two differ in ACCURACY, not in measurement cost: the paper states
#   "both methods require only the original measurements used for the ground
#   state", with ZC at "extremely high accuracy" while ZA is "qualitatively
#   correct in some cases but produces sub-variational solutions in others";
#   the extra structure of ZC is what buys the accuracy         (Fig. 7 caption)
# both accuracy comparisons are a single numerical demonstration on H2, not
#   a bound                                                            (Fig. 7)

return  the eigenpairs (C, E) -- E holds ground- and excited-state energy
    estimates within Bf1, C the coefficients of each state in the
    a_i^dagger a_j |Psi> basis
# the identical C, with O^ij_kl computed in place of H^ij_kl, reads off <O>
#   for any of these states with no new optimization -- the mechanism this
#   method is reused for whenever the wanted quantity is not the energy

McClean, Schwartz, Carter and de Jong buy excited states with measurements, not coherence. At the linear-response level $k=1$ of the hierarchy of Eq. (8), whose $k$ index pairs each range over $[1,M]$ spin-orbitals — the paper writes the index ranges, not the counts — that is $M^2$ vectors $a_i^\dagger a_j |\Psi\rangle$, $M^4$ labelled matrix elements $O^{ij}_{kl}$, and one generalized eigenvalue problem $H_{LR}C = S_{LR}CE$ (Eq. 7) solved classically. The 1- and 2-RDM come from the ground-state run; for the extra 3- and 4-RDM the authors claim only "a polynomial number of additional measurements", never an exponent. Their one explicit count is the commutator form of Eq. (59), which "depends at most on the 3-RDM requiring only $O(M^6)$ terms to be measured", and is exact only when $|\Psi\rangle$ is exact; the ZC and ZA cumulant approximations need nothing beyond the ground-state 2-RDM. No shot count or precision bound is given. Truncated levels the authors call "efficient"; at $k=N_e$ the classical diagonalization "scales exponentially in the size of the system".

• The $\mathrm{H}_2$ numerical demonstration under the variational channel state model

  • Hybrid Quantum-Classical Hierarchy for Mitigation of Decoherence and Determination of Excited States

  The method paper's own numerics, and the authors state the target in one sentence: "We assess the performance of the QSE extension to the original hybrid quantum classical approach on the spectrum of a simple molecule, H2 in a minimal STO-6G basis [29] under the Jordan-Wigner (JW) qubit encoding [30], using the VCS model." One molecule carries both halves of the paper's title. The excited states are shown from an exact reference, the decoherence mitigation from a degraded one, and they are different figures on different reference states rather than one experiment that demonstrates both at once.

  Three calculations on the same molecule: the first two differ in what the reference state is, the third in how the subspace matrix elements are approximated. The first applies the fermionic linear-response expansion to an exact ground state: "First we examine the performance of the fermionic LR expansion in determining excited states on the exact ground state of H2. This allows one to understand properties of the method in situations where very good approximations to the ground state may be prepared." The second replaces that reference with the optimum of the paper's own noise model, the variational channel state preparation, defined as "the preparation of an arbitrary pure state followed by the action of a quantum channel defined by a set of Kraus operators" and intended "to allow one to the study the optimal possible performance of quantum-classical variational algorithms in experiments separate from the considerations of ansatz choice or experimental protocol"; the channel is amplitude and phase damping at $T_p/T_2 = 0.05$, which the supplement glosses as "a total gate sequence time corresponding to roughly 5% of an expected coherence time". The third holds the reference fixed and swaps the matrix elements for the two cumulant approximations ZC and ZA, which reach the excited states from the ground-state measurements alone. A fourth ingredient appears in Fig. 5 as a constraint rather than a run: the kink in the noisy curve is repaired by "enforcing the correct spin symmetry", the $S^2 = 0$ sector.

  $\mathrm{H}_2$ under the Jordan-Wigner encoding in "a minimal STO-6G basis", which is also the basis the Fig. 4 caption names. A different basis appears once elsewhere in the paper — the supplement's channel-count comparison of Fig. 6 is run on "the electronic ground state of H2 in a STO-3G basis" — but Fig. 6 plots channels against the exact solution and carries no QSE curve, so it is a separate calculation rather than a second answer for this one. Beyond the basis the inputs are not printed: there is no geometry list, no source for the electronic integrals and no table of Hamiltonian coefficients. The internuclear separation axes of Figs. 5 and 7 are labelled $R$ (Å), with printed ticks reaching 3.0 — Fig. 7 from 0.0, Fig. 5 from 0.5.

  There is no code artefact to point at. The paper names no repository, no package, no library and no code-availability statement anywhere in its text: a pattern search over the full text extraction for repository hosts, availability phrasings, and language and library names returns only prose uses of "implementation", and the paper's single software-adjacent citation is the STO-nG basis-set reference, cited for the basis and not for any code. What survives of this demonstration is its figures, so the honest reading is that the construction was published without the code that produced the plots, which is a different fact from the code having been lost.

  Curves, not error bars. From an exact reference the truncation costs nothing on this system: "In this case the linear response (LR) approximation from an exact reference is sufficient to capture exactly the excited states with the correct number symmetry, as the excitation operators conserve the number from the exact reference state", and the text adds that LR "confines one to the physical subspace of Ne = 2 particles exactly even though the Jordan-Wigner transformation encodes the unphysical space of Ne = 0 to the number of spin-orbital sites". Under the damping channel the claim is comparative and unquantified: "one sees that in the linear response subspace (AP LR), the solution quality in increased with respect to the optimal solution of the quantum channel model (AP)", with no number attached. For the two cheap approximations, "ZC achieves an extremely high accuracy, while ZA is qualitatively correct in some cases but produces sub-variational solutions in others". No simulator, shot count or sampling noise is named anywhere in the demonstration: these are calculations on exactly solved channel states, so nothing here measures what finite measurement statistics would do to the generalized eigenvalue problem.

• The $\mathrm{H}_2$ spectrum experiment on two superconducting transmon qubits

  • Robust determination of molecular spectra on a quantum processor

  The device run of this method, by a Berkeley and LBNL group. The stated goal is the whole spectrum: "We demonstrate a complete implementation of the Variational Quantum Eigensolver (VQE), augmented with a novel Quantum Subspace Expansion, to calculate the complete energy spectrum of the H2 molecule with near chemical accuracy." How much of the method paper's authorship is here depends on which version is read: the arXiv preprint cited above lists eight authors and carries three of the method paper's four — McClean, Carter and de Jong — while the published version, Phys. Rev. X 8, 011021 (2018), lists nine under the different title "Computation of Molecular Spectra on a Quantum Processor with an Error-Resilient Algorithm", the added author being M. E. Kimchi-Schwartz, who is the method paper's second author under the name Mollie E. Schwartz, so all four are on the published record. The second half of the paper's purpose is to test a conjecture rather than to compute: "it was conjectured in ref. [20] that the inclusion of specific measurement operators expanding the number of states under consideration, the QSE could improve the accuracy of the initial VQE ground state estimate", ref. 20 being the method paper, and the answer given is "we find experimental support for this conjecture" — support, not a bound.

  The construction is the paper's own restatement of the method. It measures "the matrix elements coupling the states generated by these operators", $H_{ij} = \langle\psi_0|O_i^\dagger H O_j|\psi_0\rangle$, together with "the corresponding identity operator in this space, also known as the overlap matrix", $S_{ij} = \langle\psi_0|O_i^\dagger O_j|\psi_0\rangle$, and then proceeds "by solving the generalized eigenvalue problem HC = SCE". The expansion set is Pauli rather than fermionic: "While S1 could serve as a subspace, a more natural choice when working with qubits involves the set of single Pauli flips", the set they call the linear response expansion. The reference state comes from a six-parameter circuit — two single-qubit rotations, in amplitude and phase, and a bSWAP entangler in length and phase — minimised by "a particle-swarm optimizer (PSO)", chosen because "it is likely to avoid getting trapped in local minima and it is more robust to noisy objective-function calls". Two choices depart from what a scaled-up run would do, and both are stated. Expectation values come from full two-qubit tomography, "a set of 32 tomographic measurements", from which the needed correlators are then calculated. And because "in future implementations of VQE on larger quantum systems, full tomography will be impossible" so that at scale "only the necessary correlators will be directly measured", the reconstruction deliberately "did not use any method such as maximum-likelihood estimation which enforces physicality (positivity and trace-normalization) on the result" — in the paper that future-scaling argument is the reason for dropping the physicality constraint, not the reason for doing tomography. The expansion is then applied not once but to every intermediate state: "we apply the linear-response QSE to each individually-reconstructed density matrix recovered during the minimization process", and the resulting swarm of energies is binned at 1.5 mHa, convolved with a Gaussian filter of standard deviation 7.5 mHa, and read off by peak-finding.

  $\mathrm{H}_2$ in the STO-3G basis reduced to two qubits, and the reduction is spelled out rather than cited: after the one- and two-electron integrals, "the Hamiltonian was projected into a particle conserving manifold defined by a determinantal configuration interaction expansion that does not flip spins", with the four determinants mapped explicitly onto the four computational basis states, giving the six-term $H_Q(R) = g_0(R) + g_1(R)\sigma_z^1 + g_2(R)\sigma_z^2 + g_3(R)\sigma_z^1\sigma_z^2 + g_4(R)\sigma_y^1\sigma_y^2 + g_5(R)\sigma_x^1\sigma_x^2$ of Eq. (S3). The coefficients are published, not described: "The exact coefficients used in this work are given in Table I of this SI." The sweep is "45 internuclear distances between 0.05 Å and 3.85 Å".

  One package is named and no repository is. "For the classical optimization routine, we used the Pyswarm package, a Python implementation of particle-swarm optimization"; the readout electronics are named down to the digitizer, an AlazarTech 9373 ADC, after which "the qubit state is determined in software". Nothing else about that software is stated, and the paper carries no code-availability statement, so the classical half of this run — the tomographic reconstruction, the subspace matrices, the generalized eigenvalue solve, the binning and peak-finding — is described but not released.

  Two transmon qubits whose coherence times the paper prints as 16 µs $T_1^A$, 13.5 µs $T_2^{*A}$, 12 µs $T_1^B$ and 3.5 µs $T_2^{*B}$, with single-qubit pulses of 50 and 70 ns at "fidelities of ∼99%", a bSWAP of up to 310 ns approaching "a fidelity of ∼96%", and a whole preparation-and-measurement sequence under about 1.5 µs. The run is small enough to price: "Each function evaluation consists of 10,000 acquisitions and lasts on the order of a minute, resulting in a total algorithm run time of approximately four hours." What is achieved is stated per state, not for the spectrum as a whole: "Chemical accuracy, the level required to make realistic chemical predictions, is achieved for the ground and highest excited state across a wide range of internuclear distances. Estimates of the second and third excited state energies are generally within an order of magnitude of this level", chemical accuracy being $1.6 \times 10^{-3}$ Ha. On the ground state the paper prices the expansion twice and not identically: the Fig. 4 caption claims "an improvement of more than an order of magnitude over the bare VQE estimate" for "the majority of internuclear distances computed", while the main text describing that same figure says the complete linear-response expansion "produces a reduction in the energy estimate error of almost two orders of magnitude over the entire range". A supplementary dataset qualifies even that: the Fig. S4 run was restarted mid-sweep for technical reasons, after which "the linear-response is no longer able to get below chemical accuracy", "most likely due to calibration errors in our tomographic reconstructions", and chemical accuracy is recovered there only by adding the $\sigma_x^1\sigma_x^2$ correlator to the expansion. The failure mode is reported too: "An unphysical spurious state emerges at internuclear distances greater than ∼ 1.2 Å due to uncorrected incoherent errors", removable by continuity in $R$ or by enlarging the operator set "at the cost of an increased tomographic measurement overhead", and "the exact conditions for the presence of a spurious state are currently being investigated".

• The quantum subspace expansion routine in myQLM-fermion

  A released library function rather than a paper's supplement: `apply_quantum_subspace_expansion` in the chemistry subpackage of myQLM-fermion, the open-source many-body layer of the myQLM stack, whose package metadata gives the author as Atos Quantum Lab and whose documentation site is copyrighted to Bull SAS. Its tutorial notebook credits the construction to this method's authors — "McClean et al. drew inspiration from Linear Response Theory to design an extension to the VQE, the Quantum Subspace Expansion (QSE)" — and lists three references at the head of the page, https://arxiv.org/abs/1603.05681, https://arxiv.org/abs/1707.06408 and https://arxiv.org/abs/1807.10050. What the notebook says the routine is for is worth reading twice, because it is not the excited states: the aim is stated as expanding "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".

  The routine implements the qubit form of the expansion, not the fermionic one: the caller passes `expansion_operators`, a list of spin observables, and the docstring builds $H^{(s)}_{i,j} = \langle\Psi^\star|\hat{O}_i^\dagger \hat{H}\hat{O}_j|\Psi^\star\rangle$ and $S_{i,j} = \langle\Psi^\star|\hat{O}_i^\dagger\hat{O}_j|\Psi^\star\rangle$ on the VQE output state, then solves $H^{(s)}\vec{x} = E S \vec{x}$. Every element is a separate job on the QPU handler — the double loop forms `h_expanded = op_l.dag() * hamiltonian * op_r` and `overlap = op_l.dag() * op_r` and submits each as an observable on the same state-preparation circuit — with `nbshots` defaulting to 0, which the docstring glosses as "infinite number of shots", and elements below a `threshold` of $10^{-15}$ set to zero. A companion helper `build_linear_pauli_expansion` constructs the first-order set, "first-order all-qubit expansion from the listed Pauli gates". One detail decides what this artefact can be used for as shipped: the generalized problem is solved with `scipy.linalg.eig` and the result is then reduced by taking the minimum over the real parts of the returned eigenvalues, so the function returns the lowest eigenvalue only. The excited eigenvalues are computed and discarded; a caller who wants them must pass `return_matrices=True` and diagonalize $H^{(s)}$ and $S$ themselves.

  The tutorial's molecule is $\mathrm{H}_2$ with its integrals hard-coded into the notebook — a $2 \times 2$ one-body array, the two-body tensor, orbital energies, a nuclear repulsion of 0.7137539936876182 and natural-orbital occupation numbers noted as "computed from 1-RDM (computed in CISD here)" — with an optional cell that regenerates them through PySCF at the geometry `[('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))]` in the `sto-3g` basis. The state is a UCC ansatz in the parity encoding, optimised first by a `ScipyMinimizePlugin` running COBYLA and then by a second cell whose heading and import both say SPSA but whose code instantiates another `ScipyMinimizePlugin`, so no SPSA optimisation runs; it is that second cell's parameters the expansion is built on. The expansion set chosen for the demonstration is two operators: the identity and $Z_0Z_1$.

  Python and Apache-2.0, with the package's own metadata disagreeing with the file it ships. The repository is https://github.com/myQLM/myqlm-fermion, default branch `rc`, with a 201-line Apache-2.0 LICENSE at its root; the file is `qat/fermion/chemistry/qse.py`, which defines `apply_quantum_subspace_expansion`, the private `_build_quantum_subspace_expansion` that fills the two matrices, and `build_linear_pauli_expansion`, and `apply_quantum_subspace_expansion` is re-exported from `qat/fermion/chemistry/__init__.py`. The PyPI distribution `myqlm-fermion` 1.4.0 carries a free-text licence field reading "Atos myQLM EULA" while bundling that same Apache-2.0 text as its own LICENSE and declaring it as `License-File`, so the two statements sit inside one distribution rather than diverging between repository and index, and the wheel's copy of the module is byte-identical to the repository's at 239 lines. The notebook is a separate repository, https://github.com/myQLM/myqlm-notebooks on branch `master`, at `tutorials/digital_quantum_simulations_spin_fermion/qat_fermion_quantum_subspace_expansion_h2.ipynb`; that repository declares no licence.

  The notebook publishes none: every one of its code cells is stored with its output stripped, so its comparison of the VQE and QSE errors prints only when a reader runs it. The routine's own documentation does publish a pair, because the docstring example sits under a Sphinx `run-block` directive that executes at build time — on the two-qubit $XX + YY + ZZ$ Hamiltonian whose exact ground-state energy the example states as $-3$, the built page at https://myqlm.github.io/04_api_reference/module_qat/module_fermion/apply_quantum_subspace_expansion.html prints "E(VQE) = -3.0 (err = 0.0 %)" and "E(QSE) = -3.0 (err = 0.0 %)", the same value from both steps, so the expansion buys nothing on that toy. No hardware is involved anywhere: the example and the notebook both run on `get_default_qpu`, myQLM's local simulator.

• The AlgorithmQSE class in Quantinuum's InQuanto

  A commercial library's implementation, recorded from generated documentation because no source is published. Quantinuum's InQuanto exposes `AlgorithmQSE`, and its manual states the purpose in this method's own terms: "The Quantum Subspace Expansion (QSE) obtains extra correlation and energetics of excited states, on top of the VQE ground state" by "doing a linear excitation expansion", with the subspace vectors written $a^\dagger_i a_j |\Psi_{\mathrm{VQE}}\rangle$ and the generalized eigenvalue problem written $HC = SCE$. The page cites the method paper directly, as bibliography entry 15, in its published form: "J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong. Hybrid quantum-classical hierarchy for mitigation of decoherence and determination of excited states. Phys. Rev. A, 95():042308, 2017." Unlike the myQLM routine recorded beside it, this one keeps the excited states: the manual says the generalized problem gives "the optimal solution for the ground and excited states within this subspace", and that the eigenvalues and eigenvectors are returned as `final_value` and `final_states`.

  Two routes are documented and they differ in how much the caller controls. The explicit route builds a `QSEMatricesComputable` from an ansatz, a Hermitian operator and a list of expansion operators, evaluates it under a measurement protocol, and hands the resulting $H$ and $S$ to `pd_safe_eigh`. The black-box route wraps that same computable in `AlgorithmQSE`, which the manual says handles numerical instabilities in the diagonalization "by removing near linear dependencies in the basis", at the price that "the user has less control over the specific details of the process" — the emulator tutorial names launching, retrieving and PMSV noise mitigation as things it cannot do. The expansion set in both documented examples is spin-adapted rather than the raw fermionic set: `generate_subspace_singlet_singles` supplies singlet single excitations $E_{ij} = \sum_{\sigma}^{\mathrm{spin}} a_{i\sigma}^{\dagger} a_{j\sigma}$, mapped to qubits by `QubitMappingJordanWigner`. Expectation values are measured by the `PauliAveraging` protocol, and the emulator run partitions them with pytket's `CommutingSets` strategy and verifies symmetry with PMSV on the stabilizer $Z_0 Z_1 Z_2 Z_3 Z_4 Z_5$.

  Two molecules, both in the STO-3G basis. The manual's example is $\mathrm{H}_2$ loaded from a file shipped with the library, `get_system("h2_sto3g.h5")`, with a UCCSD ansatz whose parameters are supplied as the array 0.0, 0.0, -0.107233493519281 carried over from an earlier VQE example rather than re-optimised. The tutorial's molecule is methane, built from a Z-matrix with C-H bonds of 1.083000 and angles of 109.471000 and 120.000000, through the `ChemistryDriverPySCFMolecularRHF` driver with `frozen=[0, 1, 2, 3, 7, 8]` and point-group symmetry enabled, Jordan-Wigner encoded and compressed at `abs_tol=1e-6` into a six-qubit Hamiltonian whose terms the page prints in full; its reference state is a chemically aware UCCSD ansatz optimised by a noiseless VQE on the Qulacs simulator before any emulator work begins.

  No source. InQuanto is distributed by Quantinuum, is not on PyPI, and its own "Software License" page carries no licence text — it directs access inquiries to Quantinuum support. Everything recorded here is read off the generated documentation at https://docs.quantinuum.com/inquanto, which identifies itself as InQuanto 6.1.1 and is copyrighted 2025 to Quantinuum Ltd: the manual page `manual/algorithms/algorithms_qse.html` and the tutorial `tutorials/InQ_htut_qsys_qse.html`. The names those pages expose are `inquanto.algorithms.AlgorithmQSE`, `inquanto.computables.composite.QSEMatricesComputable`, `inquanto.core.pd_safe_eigh`, and `PauliAveraging` and `PMSV` in `inquanto.protocols`, with the PySCF driver in the separate `inquanto.extensions.pyscf`. There is no paper: the project asks to be cited as a `@misc` entry dated 2022 to Quantinuum. Nothing below the API surface was read, so no claim here is about the implementation's internals.

  Both documented runs print numbers and neither is on a device. The manual's $\mathrm{H}_2$ example runs `PauliAveraging` on an `AerBackend` at 2000 shots per circuit and reports `final_value` as the three eigenvalues -1.121, -0.141 and 0.227, with a timing block giving about one second to measure and calculate the $H$ and $S$ matrix elements and under a millisecond to solve the eigenvalue problem. The methane run goes to a Quantinuum emulator and not to a Quantinuum processor: the backend is `QuantinuumConfig(device_name="H2-1E")` reached through qnexus, which the page describes as designed "to emulate the associated devices" and as running "remotely on a server", at 6000 shots per circuit over 30 measurement circuits, and its printed eigenvalues are -39.855, -39.68, -38.889, -38.865, -38.488, -37.761 and -37.327. One qualifier the page makes easy to miss: those methane numbers come from the explicit `QSEMatricesComputable` and `pd_safe_eigh` route, because the `AlgorithmQSE` cell that follows has its `run` and report calls commented out, so no output from the class itself is published for methane. Neither run is given a reference energy, an error bar, or a comparison against a classical calculation.

Every step this method names is listed under Requires above. It walks its own span in one hop and calls out to the rest — that is a fact about the recorded route, not a claim that the span is simple.

Different approaches

• Variational quantum deflation

  Find the ground state first, then run the same variational search again with a term that punishes overlap with every state already found. Each state is reached by pushing the search off the ones below it, so they have to be found in order and each one costs another pass through the whole loop.

• 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 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.

• Quantum subspace expansion

  Measured response operators around a VQE state define a generalized eigenproblem for excitations and mitigation.

• Hybrid Quantum-Classical Hierarchy for Mitigation of Decoherence and Determination of Excited States↗Jarrod R. McClean, Mollie E. Schwartz, Jonathan Carter, Wibe A. de Jong · 2016theory yes · simulation yes · hardware no