MethodLayer 0
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.
A Hermitian H reachable as a sum of terms that can be measured separately; a statement of which state is wanted — an index k, a symmetry sector, or a target energy to sit nearest; for most of the methods here, the ground state already solved, and for the deflation route every lower state as well; a target additive error and a shot budget.
A scalar estimate — of that eigenvalue, or of the gap between it and the ground state — together with the run budget it consumed. Some methods return a whole low-energy subspace at once and others return one state per run, and that difference is a cost, not a convenience.
Same contract as the slot it fills.
This one, drawn
From Hamiltonian whose excited state is wanted to Number with an error bar
A circle is an object you are holding. This method is drawn heavier, opened into its own steps; the other lines between the same two ends are the alternatives recorded for the same slot. Circles are named on hover, and each one is a link.
7 lines here have ways through that this figure does not open. The map opens them in place. See it on the map
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
Kuroiwa and Nakagawa analyse the two penalty forms already in use and separate them: one "works properly in that eigenstates obtained by the VQE with the penalty term reside in the desired symmetry sector", and they "further give a convenient formula to determine the magnitude of the penalty term, which may lead to the faster convergence of the VQE". The other is not a weaker option but a wrong one — they "prove that the other type of penalty terms does not work for obtaining the target state with the desired symmetry in a rigorous sense and even gives completely wrong results in some cases". Two conditions follow. The target symmetry sector has to be known in advance, since it is what the penalty is written against; and this route reaches a symmetry-resolved spectrum rather than the k-th state as such, so it answers "the lowest state with these quantum numbers" and not "the third state overall". Their own framing is that the work "lay[s] the theoretical foundation for the use of the VQE with the penalty terms", validated by numerical simulation.
Requires
Every step this method names moves its route along, so there is nothing it needs alongside them.
Example
given Hamiltonian H on n qubits and a conserved-quantity observable C with
[H, C] = 0 (particle number N, S^2, S_z, ...) (Sec. II C)
target symmetry sector: an eigenvalue c of C, fixed before the run --
c determines the sector, it is an input here, not an output (Sec. II C)
ordinary ansatz |psi(theta)> = V(theta)|ref>,
V(theta) = V_l(theta_l) ... V_1(theta_1) (Eq. 1-2)
# |ref> is the paper's |psi_0> of Eq. 1, "some reference state" -- one
# FIXED state (|00001111> in the H4 runs, Sec. IV) reused unchanged by
# every round. Renamed here because Sec. II B writes |psi_0> for a
# different object: the converged ground state found earlier.
# no symmetry structure is built into V(theta); the guarantees below
# are then derived assuming this ansatz can reach an ARBITRARY
# n-qubit state, |psi(theta)> = sum_i a_i |i>, because the paper's
# purpose is to study the two cost functions themselves, not a
# specific circuit's reach (Sec. III, Eq. 8)
for an excited state: eigenstates |psi_0>, ..., |psi_{k-1}> already
obtained by earlier rounds of this same procedure (Sec. II B)
requires the target sector known in advance as the scalar c -- this route
resolves a symmetry sector, so it answers "the lowest state with C = c",
not "the k-th state overall" without one; denote that state's energy
E_{i0}, meaning C_i != c for every i < i0 under E_0 <= ... <= E_{2^n-1} (Sec. III)
# --- cost function: penalize deviation from the sector, don't build it in ---------
if excited state k:
L(theta) = LVQD(theta) = <psi(theta)|H|psi(theta)>
+ sum_{i=0}^{k-1} beta_i |<psi_i|psi(theta)>|^2 (Eq. 4)
# equivalently plain VQE against H' = H + sum_i beta_i |psi_i><psi_i|,
# beta_i > 0 chosen only to enforce orthogonality against the
# already-found lower states (Sec. II B)
else:
L(theta) = LVQE(theta) = <psi(theta)|H|psi(theta)> (Eq. 3)
F(theta) = L(theta) + mu_C * <psi(theta)|(C - c)^2|psi(theta)> = F^(1) (Eq. 5)
# the rival form F^(2) = L + mu_C*(<psi|C|psi> - c)^2 is NOT used here:
# proved to miss the target by O(1/mu_C) at every finite mu_C even in
# the best case, and when (c, E_{i0}) sits interior to the (C,E) convex
# envelope, no mu_C -- however large -- recovers it at all (Sec. III B, Eq. 18-23)
# --- set the penalty weight --------------------------------------------------------
# proved sufficient, but needs the full spectrum so is not directly usable:
# mu_C >= max_{i<i0} (E_{i0} - E_i) / (C_i - c)^2 (Eq. 11)
C_min = smallest gap between distinct eigenvalues of C
# 1 for the particle-number operator N; 3/4 for S^2; 1/2 for S_z (Sec. III A)
mu_C_simple = (E_{i0} - E_0) / C_min^2
# a looser but simpler bound than Eq. 11, obtained from E_i >= E_0
# for all i < i0 -- still proved sufficient (Eq. 12-13)
# E_{i0} - E_0 is itself unknown before the run; two ways to fill it in:
mu_C_ce = (E_{i0}_classical - E_0_classical) / C_min^2 # e.g. Hartree-Fock
# or Moeller-Plesset estimates of the two energies (Eq. 14)
mu_C_rough = (2 / C_min^2) * sum_j |c_j| # H = sum_j c_j P_j (Pauli decomposition)
# from the rigorous bound E_{i0} - E_0 <= 2||H|| <= 2 sum_j |c_j| (Eq. 15)
mu_C := mu_C_rough # pick exactly ONE of the three; F^(1) takes a single scalar
# overestimating E_{i0} - E_0 still finds the target, only slower;
# UNDERestimating can move the minimum off the target -- so
# mu_C_rough (never an underestimate) is the safe default (Sec. III A)
# but the paper also says mu_C_rough "may be too large for the fast
# convergence", and that a large mu_C slows the optimizer, so
# mu_C := mu_C_ce is the faster choice where the classical estimate
# can be trusted not to underestimate the gap (Sec. III A, Eq. 15)
# several conserved quantities at once: one term and one mu_C^(l) per
# C^(l), each set by the same formula with its own C_min^(l) (Eq. 16-17)
# --- optimize -----------------------------------------------------------------------
run a classical optimizer on F^(1)(theta) with that mu_C, theta -> theta*
# the paper's own numerics use BFGS on noiseless simulation with exact
# expectation values -- no shot noise, no device noise (Sec. IV)
while optimizing, track <psi(theta)|(C - c)^2|psi(theta)>
# if this has not fallen near 0 by convergence, mu_C was too small and
# the run missed the sector -- raise mu_C and rerun (Sec. III A)
return |psi(theta*)>, F^(1)(theta*) as the estimate of the sector's E_{i0}
# under a global depolarizing channel this minimizer is UNCHANGED -- the noisy
# F^(1) equals (1-p) times F^(1) plus a theta-independent constant, so the
# same theta* minimizes it; F^(2)'s minimizer instead shifts at O(p) (Sec. III C, Eq. 26)Cost, as the source states it
Kuroiwa and Nakagawa price penalty-constrained VQE in Pauli measurements — (cost-function evaluations) (Pauli operators per evaluation) — and bound neither factor: the working penalty measures once per conserved quantity, more Pauli operators than alone but never counted, plus overlap terms for the -th state by VQD. Their figures are measured, not proved — 8-qubit triplet, depth-12 symmetry-preserving ansatz, BFGS, noiseless, averaged over ten initial parameter sets. takes Pauli measurements at , rising to at where the residual degrades from to Ha. That cheapest accurate point still exceeds the the unguaranteed needs at , though the ordering reverses in their VQD run. On with BFGS the classically estimated needs evaluations against for the loose bound, which converged less accurately. Every guarantee assumes an ansatz spanning arbitrary -qubit states.
Implementations
3
The constrained VQE and VQD runs on and in Sec. IV
Section IV is the authors' own check on their theorems rather than an application: "In this section, we numerically simulate the constrained VQE/VQD for two molecules, H2 and H4 , to validate our results presented in the previous section." The problem it was built against is stated in the abstract — "systematic analysis of the effect of the penalty term has been lacking, and the use of the penalty term in the VQE has not been justified rigorously" — so these runs exist to separate the two penalty forms and to test the formula for the penalty coefficient, not to produce a chemistry number.
Two targets, each named by its symmetry sector rather than by an index: the T1 state, "the ground state in the triplet sector" at and , computed by constrained VQE, and the S1 state, "the first excited state in the singlet sector" at and , computed by constrained VQD; since "the ground state is the spin singlet ( and ) state for both the H2 and H4 molecules", both targets are excited states of the full Hamiltonian. Which conserved quantities are penalized varies by molecule and target (their TABLE I), and for no penalty is placed on the electron number because the symmetry-preserving ansatz already preserves the electron number of the reference state, which they take as . The ansatz is the hardware-efficient type at depth for the constrained VQE and for the constrained VQD; uses "the real-valued symmetry-preserving type ansatz of the depth ". Section IV then runs two separate sets of simulations, and their hyperparameter choices do not carry over from one to the other. The Sec. IV A set, which tests the coefficient formulas, fixes the penalty coefficient two ways, from CISD energies computed in PySCF and from the rough bound on the spectral norm; takes the deflation weight as "based on the original VQD proposal" and estimates it the same two ways; and optimizes with Powell, conjugate gradient and BFGS "implemented in SciPy, a numerical library of Python" with default parameters. The Sec. IV B set, which compares the two cost functions, instead sweeps over 0.01, 0.1, 1, 10 and 100; constrains only in its VQD run, so that "the S1 state is obtained as the fourth excited state of the constrained VQD ()" with the four hyperparameters "all set to 3.0"; and reports that "for both simulations, we use the BFGS optimizer, which shows the best performance in Sec. IV A". Every tabulated figure in both sets is "the average for the ten trials" over ten different initial values.
Two molecules only: "a hydrogen molecule H2 at bond distance Å and a hydrogen chain H4 where four hydrogen atoms are aligned in line with the identical bond distance Å", both in the STO-3G minimal basis with a restricted Hartree-Fock reference, giving 4 and 8 qubits. The electronic Hamiltonian is not hand-built: it is prepared "using PySCF and Openfermion" and mapped to qubits by the Jordan-Wigner transformation, and the conserved quantities , and "are also transformed into operators on qubits by Jordan-Wigner transformation".
"All simulations are performed by the high-speed quantum circuit simulator Qulacs", with the classical optimizers taken from SciPy. The circuits are noiseless and the expectation values exact — "We do not include any noise for quantum circuit simulations, and exact expectation values are used in numerical simulations" — so nothing here is a shot-noise or device measurement. The arXiv v2 read here names no repository and carries no code-availability statement, so this entry records a run and the stack it ran on, not a released artefact.
TABLE V and TABLE VI each put both cost functions on a single target across the same five penalty coefficients, and TABLE VI is the sharper of the two: the S1 energy of by constrained VQD under the single constraint , exact value Ha. reaches it — residuals and Ha at and , then Ha at , and — while does not converge toward it at all: , , , and Ha across the same five coefficients. The paper reads this as the case where no coefficient can help, "the target eigenstate is inside the convex envelope explained in Sec. III B and we simply cannot obtain that state by optimizing the cost function ". Here is also the cheaper of the two, to Pauli measurements against to for , which is the reverse of the ordering in their triplet run in TABLE V, where 's residual instead shrinks as grows, from Ha at to Ha at .
Tangelo's `penalty_terms` module and the `penalty_terms` option on `VQESolver`
Tangelo describes itself as "A python package for exploring end-to-end chemistry workflows on quantum computers and simulators", and the module that concerns this method states its own purpose in its docstring: "This module defines the penatly terms that can be added to the target fermionic Hamiltonian, providing the ability to restrict the Hilbert space of solutions using VQE" — the misspelling is the file's. The docstring then says "For example usages see" and names two papers, Ryabinkin, Genin and Izmaylov's constrained VQE and Greene-Diniz and Muñoz Ramo's generalized unitary coupled cluster. Both sit inside the group of earlier studies Kuroiwa and Nakagawa cite as having suggested penalty terms, and the first is the one they cite for , the form they prove does not work: Ryabinkin, Genin and Izmaylov's own Eq. (12) adds to the energy, squaring the expectation value, and their paper defines no other penalty functional. The package's own paper draws the feature as one of the optional inputs to `VQESolver`, labelled "Penalty terms" beside the initial variational parameters and the qubit-encoding options.
Three builders return a `FermionOperator`, each documented as the square of an operator: `number_operator_penalty` (line 30) returns "mu*(\hat{N}-n_electrons)^2", `spin_operator_penalty` (line 54) returns "mu*(\hat{Sz}-sz)^2" and `spin2_operator_penalty` (line 76) returns "mu*(\hat{S}^2-s2)^2"; `combined_penalty` (line 102) sums whichever of the keys "N", "Sz" and "S^2" carry a positive prefactor and raises `KeyError` on any other key. What is squared is the operator and not its expectation value: each builder prepends the scalar term `[(), -c]` to the operator's term list and passes it to `squared_normal_ordered`, whose body is `fe_op = list_to_fermionoperator(all_terms)`, then `fe_op *= fe_op`, then `return normal_ordered(fe_op)` (`tangelo/toolboxes/operators/operators.py`, line 352). `VQESolver.build` maps that operator through the same fermion-to-qubit mapping as the Hamiltonian and folds it in with `self.qubit_hamiltonian += pen_qubit` (line 203), so the quantity minimized is the energy plus — the form for which Kuroiwa and Nakagawa prove that minimization "can yield the desired state/energy when the penalty coefficient is larger than a certain threshold", not the of the paper the docstring points at. The threshold is not part of the feature: each of the three builders declares `mu=1` as its default and documents it only as "Positive number in front of penalty term", so the coefficient is the caller's to choose. The same option is threaded through `sa_vqe_solver.py`, `sa_oo_vqe_solver.py`, `iqcc_solver.py`, `iqcc_ilc_solver.py` and the two phase-estimation solvers, and the same three builders are reached by a second route in `tangelo/toolboxes/ansatz_generator/_qubit_mf.py`, whose `penalize_mf_ham` (line 188) adds them to the mean-field Hamiltonian of the QMF ansatz at `qmf.py` line 128; `qmf.py` lines 121 to 122 write a per-term default coefficient of 1.5, but the loop meant to apply it tests `if params is None` against the tuples it has just built, so it never fires and the coefficient reaches `penalize_mf_ham` exactly as the caller supplied it. The `deflation_circuits` and `deflation_coeff` options declared beside it on the same solver are a different construction and belong with variational quantum deflation.
GitHub, `sandbox-quantum/Tangelo` — the older `goodchemistryco/Tangelo` path redirects there — `main` at `4b478a10`, 2026-08-05, read 2026-09-03; Python; the repository is archived. `tangelo/toolboxes/ansatz_generator/penalty_terms.py` is 150 lines and was last touched by `ae8895b4`, 2024-05-21, the commit that brought release v0.4.3 into `main`; `tangelo/_version.py` ends `__version__ = "0.4.3"` and `setup.py` names the distribution `tangelo-gc`. Licence is Apache 2.0: `LICENSE` opens "# Copyright SandboxAQ 2021-2024" and "# Licensed under the Apache License, Version 2.0", and `penalty_terms.py` repeats that header. The consumer is `tangelo/algorithms/variational/vqe_solver.py` (720 lines), where `penalty_terms` is popped from the options dict at line 109 and consumed at lines 195 to 203.
What the tests assert is narrower than the method. `tangelo/toolboxes/ansatz_generator/tests/test_penalty_terms.py` (233 lines) has five tests — `test_number_penalty_terms_hf_ref`, `test_spin_penalty_terms_hf_ref`, `test_spin2_penalty_terms_hf_ref`, `test_doublet` and `test_triplet` — each of which builds the penalty operator at `mu=100` for `mol_H2_sto3g` and `mol_H4_sto3g`, evaluates it on a reference circuit already in the target sector, and checks `self.assertAlmostEqual(penval, 0.0, delta=1e-6)`. The mapping coverage is not uniform: the three Hartree-Fock-reference tests loop over jw, bk and scbk, while `test_doublet` and `test_triplet` run under jw alone, the first of them carrying the comment "may add tests for other mappings later". All of that exercises the operator, not the optimization: no test drives a constrained VQE to an excited state, and the VQE solver's own test file, `tangelo/algorithms/variational/tests/test_vqe_solver.py` (682 lines), never mentions the option.
Quket's `constraint_lambda` spin penalty
Quket heads its README "Quket: Quantum Unified Kernel for Emulator Toolbox", "Official Version 0.9.4", and describes itself as a suite that "simulates quantum computing for electronic Hamiltonian". It is worth recording here because it wires a penalty on a symmetry sector straight into a VQE cost function, and because its option list carries three separate penalty coefficients — `constraint_lambda`, `constraint_lambda_Sz` and `constraint_lambda_S2_expectation` — of which only the first is documented. That documentation is a single line of the option table, "`constraint_lambda` : The lambda value for S**2 penalty" (README line 748).
The penalty enters a single one of Quket's cost functions, `cost_exp` in `quket/ansatze/exp.py`, which `quket/vqe/vqe.py` dispatches to for the plain UCC-type ansatze — its final `else` branch names "hf", "uccsd", "sauccsd", "uccgsd", "uccd", "sauccd", "uccgd", "sauccgsd", "sauccgd", "uccgsdt" and "uccgsdtq" — and only when `Quket.multi.nstates == 0`; the projected, UHF, BCS and k-UpCCGSD paths reach `cost_proj`, `cost_uhf`, `cost_bcs` and `cost_upccgsd` instead, and `constraint_lambda` occurs nowhere else under `quket/ansatze/`, so the one documented penalty option is inert on those ansatze and on any multi-state run. Inside `cost_exp` the penalty follows immediately after the deflation term `cost += orthogonal_constraint(Quket, state)` at line 161, so a sector penalty and an orthogonality constraint against `Quket.lower_states` can be carried in the same run. Lines 163 to 167 read `if Quket.constraint_lambda > 0:`, `s = (Quket.spin - 1)/2`, `S4 = Quket.qulacs.S4.get_expectation_value(state)`, `penalty = Quket.constraint_lambda*(S4 - S2*(s*(s+1) + (s*(s+1))**2))`, `cost += penalty`; the target sector is thus set by the `spin` multiplicity keyword, and exists as an operator only when the penalty is on, built as the product `S4 = S2*S2` in `quket/quket_data/quket_data.py` and left `None` otherwise. No expectation value is squared in that expression — it is linear in and — so it is not Kuroiwa and Nakagawa's ; neither is it in general their penalty , which expands to , because line 166 puts on and carries no constant term at all; the two coincide only for a singlet target, where and both reduce to . The analytic Jacobian `jac_mpi_ana` builds a third expression again, applying to and adding the constant at `quket/utils/deriv.py` lines 178 and 180. Two further coefficients, `constraint_lambda_Sz` and `constraint_lambda_S2_expectation`, are declared alongside at `quket/quket_data/quket_data.py` lines 766 to 768 and mapped from the input keywords `lambda_sz` and `lambda_s2` in `quket/fileio/read.py` lines 72 to 74, but neither is referenced in the cost function: both are read only where `Quket.adapt.mode == 'pauli'`, in `quket/vqe/adapt.py` and in the analytic Jacobian at `quket/utils/deriv.py` lines 189 to 215, and `constraint_lambda_Sz` alone also in a finite-difference gradient at `quket/opelib/opelib.py` line 511. Neither of the two squares an expectation value either, and at `quket/vqe/adapt.py` line 392, inside the block guarded by `constraint_lambda_S2_expectation`, the coefficient actually multiplied is `constraint_lambda`.
GitHub, `quket/quket`, `main` at `b343f0fb`, 2022-09-16, read 2026-09-03; Python over Qulacs and OpenFermion; version 0.9.4. Licence is Apache 2.0 — `LICENSE` is the stock text opening "Apache License" and "Version 2.0, January 2004", and the per-file headers read "Copyright 2022 The Quket Developers" then "Licensed under the Apache License, Version 2.0". The organisation repository was created 2022-07-27 and has not been pushed to since 2022-09-16; the lead author's separate, unforked `ttsuchimochi/quket` stops at `fcfdbda8`, 2022-07-27, whose commit subject is "version 0.9", so it holds the earlier release. The files a reader opens are `quket/ansatze/exp.py` (350 lines), `quket/utils/deriv.py` (670), `quket/vqe/adapt.py` (640), `quket/vqe/vqe.py` (691), `quket/quket_data/quket_data.py` (2797) and `quket/fileio/read.py` (1195).
Nothing in the repository exercises the penalty, so there is no published run to report from it. Outside the `quket/` package the three coefficients appear only in that one README line: none of the sample inputs under `samples/` — which include an `ADAPT/`, a `Tapering/` and a `PHF/` directory — nor `test.py` sets any of them.
What it needs
3
- 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.
- 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.
- Estimate an observable 4 methods
Given the ability to prepare and a description of an observable , return a classical scalar within of at confidence . The state is never returned; only the number is.
Other ways to fill the same slot
6
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 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.
- 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
1
- Penalty-based excited-state VQE
Orthogonality or symmetry penalties augment the energy objective to exclude previously identified sectors.