GCiM

GCiM, the generator-coordinate-inspired method, computes the smallest eigenvalue of the projected generalized eigenproblem \(Hf=ESf\), where \(E\) is the projected eigenvalue, \(f\) its coefficient vector, \(H_{ij}=\langle\phi_i|A|\phi_j\rangle\) and \(S_{ij}=\langle\phi_i|\phi_j\rangle\) for normalized basis states \(\phi_i\) and the Hermitian operator \(A\) (Zheng et al., Phys. Rev. Research 5, 023200 (2023), doi:10.1103/PhysRevResearch.5.023200, the projected Hill–Wheeler equation, Eq. (13) of arXiv:2212.09205v1). result.eigenvalue is the smallest projected value after the overlap directions whose eigenvalues are at or below a cutoff are removed, and it does not establish the full-space ground energy. FixedGCIM uses basis states you supply. ADAPT grows the basis from a reference state and a pool of anti-Hermitian generators, by default adding in each round the unselected generator \(A_j\) with the largest absolute gradient \(\langle[A,A_j]\rangle\) on the current product state (Zheng et al., npj Quantum Information 10, 127 (2024), arXiv:2312.07691v3, main text p. 4 and Methods p. 9). Here the product state after \(k\) selections is \(|\chi_k\rangle=\exp(\theta_k A_{j_k})\cdots\exp(\theta_1 A_{j_1})|\phi_{\rm ref}\rangle\), where \(|\phi_{\rm ref}\rangle\) is the normalized reference, \(j_1,\ldots,j_k\) are the selected pool indices in order and \(\theta_\ell\) is the real angle of the \(\ell\)th selected generator.

Use FixedGCIM when you have trial states and want the lowest energy in their span, and ADAPT when you have a reference state and a generator pool, such as the spin-adapted fermionic pool for chemistry, and want the basis chosen for you. Chebyshev Lanczos instead builds its trial space from one initial state and powers of the operator, and QPE estimates an eigenphase or energy from controlled time evolution of a prepared state.

from nwqlib import Eigenproblem, solve
from nwqlib.algorithms.gcim import FixedGCIM

problem = Eigenproblem(A=[[1, 0], [0, -1]])
result = solve(problem, method=FixedGCIM(basis=([1, 1], [1, 1j])), seed=7)
print(result.eigenvalue)

-0.9999999999999983

The two basis states span the whole two-dimensional space, so the projected value is the lowest eigenvalue \(-1\) of \(A=\operatorname{diag}(1,-1)\), up to roundoff. Replace A and the basis columns to treat your own operator. The eigenvalue notebook (examples/gcim_lanczos_qpe_eigenvalue_intro.ipynb) applies ADAPT-GCiM to a stretched H4 chain and compares Lanczos and quantum complex exponential least-squares (QCELS) on the same problem, and the API lists every option.

Inputs and matrix access

Eigenproblem holds the Hermitian operator. FixedGCIM holds the basis columns and uses each one normalized, and ADAPT holds the reference state, the generator pool and the stopping controls. Each supplied column is read in the original coordinates of \(A\). Classical FixedGCIM keeps dense, CSR and CSC inputs and all basis columns in their original dimension, and quantum dense and sparse conversions fill any added diagonal entries with \(\operatorname{tr}(A)/d\), where \(d\) is the original dimension of \(A\).

Quantum FixedGCIM and ADAPT use the same input_conversion="auto"|"dense_pauli" and max_conversion_work=100_000_000 options as Lanczos. Automatic conversion accepts only explicit dense matrices of original dimension at most 16, and this limit does not apply to compact Pauli inputs. Choose dense_pauli explicitly for a sparse quantum matrix or a sparse custom ADAPT generator. For example, ADAPT(initial_state=[1, 1], pool=(csr_matrix([[0., 1.], [-1., 0.]]),), input_conversion="dense_pauli") accepts this finite custom generator under the same conversion-work and byte limits.

Classical ADAPT applies a sparse Hamiltonian in sparse form and never converts it to a dense matrix. When the original Hamiltonian dimension \(d\) is not a power of two, classical ADAPT pads the matrix to the next power of two with dummy diagonal entries \(\operatorname{tr}(A)/d\), keeps CSR and CSC input sparse, and zero-pads the reference state. A quantum ADAPT Plan, the construction and its costs computed before any circuit exists, counts the exact synthesis of the dense unitaries of a supplied reference circuit, and Qiskit's control of them in the Hadamard-test and joint-state circuits, against max_products and max_bytes during planning (exact synthesis).

Within one classical evaluation, a state object that appears more than once in the basis has its vector and (A-cI)@vector computed once and shared. Distinct state objects are not merged, even when they agree in shape, declared name or circuit, and saving and loading keep this sharing. The known preparation products are counted against the projected-work limit before any state is computed. The simulation work of a custom circuit whose cost the library cannot count remains unknown.

Read the result

result.eigenvalue is the smallest obtained projected value and does not establish the full-space ground energy. result.analyze(overlap_cutoff=...) uses the same obtained H and S with a new absolute cutoff, whose default is 1e-12. result.pencil contains the complex Gram matrix, the projected Hamiltonian, the kept spectrum and the numerical diagnostics below. Both matrices hold the obtained upper-triangle entries, completed by complex conjugation, except that quantum execution sets every diagonal overlap entry to one under the normalized-basis convention. Classical FixedGCIM first averages its two independently computed triangles. No diagonal jitter or positive-semidefinite repair is applied. These fields describe the projected problem and can be examined after loading a saved Result, without obtaining its matrix elements again.

Field	Meaning
overlap_eigenvalues, overlap_cutoff, kept_rank	Full Gram spectrum before the cutoff, absolute cutoff and number of positive directions used in the solve.
eigenvalues, offset_free_eigenvalues	Ascending Ritz values of the kept subspace, and the lowest one or two of them as the solve computed them from \(H_0\) and \(S\), before \(c\) is added back. The first one or two eigenvalues are these values plus \(c\), rounded once each, so two Ritz values that round to the same physical value stay distinct in offset_free_eigenvalues, which is empty when the solve failed. Pilot and main shot allocation and the score of bounded Ritz growth use these offset-free values.
overlap_condition_number	Largest divided by smallest Gram eigenvalue among the directions above the cutoff.
projected_residual	Computed \(\lVert Hf-ESf\rVert_2\) for the lowest returned projected eigenpair, evaluated on the shifted pencil as \(\lVert H_0f-(E-c)Sf\rVert_2\), which gives the same value in exact arithmetic.
projected_backward_error	Computed residual divided by \((\lVert H\rVert_F+\lvert E\rvert\,\lVert S\rVert_F)\lVert f\rVert_2\), with the physical \(H=H_0+cS\) and \(E\). This is dimensionless and unchanged by a rescaled energy unit, but adding an identity term to the operator changes it. Lanczos states its value on a normalized pencil instead (Projected quantities). The reported value is zero when the computed residual is zero.
overlap_normalization_error	Computed \(\lvert f^\dagger Sf-1\rvert\) for the lowest returned coefficient vector.

Here \(H\) and \(S\) are the projected matrices, \(f\) is the corresponding coefficient vector and \(c\) is the identity coefficient in \(A=cI+A_0\), so \(H_0=H-cS\). The expression \(E-c\) denotes the shifted Ritz value computed before \(c\) is added back. \(\lVert\cdot\rVert_2\) is the Euclidean vector norm and \(\lVert\cdot\rVert_F\) the Frobenius matrix norm. The diagnostics are floating-point evaluations. Their sizes help explain the solve and its conditioning, but they do not by themselves bound trial-space approximation or full-space ground-energy error. A zero Hermiticity defect in a matrix assembled by conjugate completion describes that representation, not an independent validation of its obtained entries. Projected verification reports chosen saved diagnostics as numerical facts, the error and reference evidence attached to a Result, without obtaining the entries again or running another eigensolve.

For sampled quantum FixedGCIM and ADAPT, the current pencil also records the processed Pauli L1 spectral enclosure (sampled_enclosure), the Ritz estimate's distance outside it (sampled_enclosure_violation), whether that distance lies within the relative numerical comparison window (sampled_enclosure_passed) and the window itself (sampled_enclosure_rtol). A noisy pencil need not obey the exact operator's Rayleigh–Ritz enclosure even when its own eigensolve has small backward error. This is a physical plausibility diagnostic, not a joint sampling bound, a proven accuracy bound or a ground-state test. The raw estimate is kept, and no clipping, resampling or additional full-space eigensolve is performed. An enclosure that is not representable stays unavailable. The summary describes the current pencil's enclosure diagnostic, including after cutoff reanalysis. Solver failure and ADAPT's stop reason remain separate, and an enclosure excursion does not change generator selection, stopping or previous-energy updates.

A numerical error bound must identify both the quantity it bounds and its assumptions. For example, if \(S\) is exactly the identity and two Hermitian projected Hamiltonians have a justified difference bound \(\lVert\hat H-H_\star\rVert_F\le\delta\), their minimum eigenvalues differ by at most \(\delta\) by the Hermitian eigenvalue perturbation bound. This bounds the matrix-perturbation contribution to the projected energy. It does not include error in the subsequent eigensolve or error from the chosen trial space. With a perturbed or ill-conditioned \(S\), a bound on \(H\) alone is insufficient.

FixedGCIM's default total-error model declares six unresolved sources (native_preparation, native_simulation, sampling, projected_solve, ground_identification and physical_model), and its analysis does not automatically produce the corresponding physical-output error bounds. The default total result.assess(absolute_tolerance=...) therefore remains INCONCLUSIVE even for a full trial basis. Additional compatible facts= can supply established bounds, but a reference energy alone does not fill the missing contributions. The absolute-error criterion does not require a known ground energy. Limitations and open work tracks the work on these bounds.

How matrix elements are obtained

Exact quantum FixedGCIM and ADAPT obtain \(H_{ij}\) and \(S_{ij}\) once per basis pair. By default each pair's exact saved state is reduced, as soon as it is produced, to the offset-free Hamiltonian entry and the overlap entry. Sampled FixedGCIM and ADAPT both measure them with grouped joint-state counts unless ADAPT's acquisition=PencilAcquisition(...) option selects termwise Hadamard tests. Every path then uses the shared generalized Hermitian eigensolver.

• Off-diagonal pair, exact. One joint preparation holds \(\phi_i\) and \(\phi_j\) on the two branches of an ancilla, with their relative phase. The pair computation applies the offset-free Hamiltonian \(H_0\) once to the right branch with the grouped Pauli routine and contracts with the left branch, and it computes the overlap directly from the two branches. The Hamiltonian entry's error combines preparation error, grouped-action rounding and complex-dot rounding. The overlap entry has only preparation and dot error.
• Diagonal pair, exact. One system preparation returns its raw squared norm as the obtained overlap. Pencil assembly uses unit diagonal overlap under the normalized-basis convention, while the raw obtained scalars and the padded coordinates remain recorded with the pair computation that produced them.
• Sampled, grouped. shots=N gives N shots per measurement setting, and a setting measures the system qubits in the common eigenbasis of one qubit-wise commuting group of Hamiltonian terms. On an off-diagonal pair it prepares the same joint state as the exact computation and also measures the ancilla, in the X basis for the real parts of the pair's matrix elements or in the Y basis for their imaginary parts. The ancilla outcome alone in the first group's settings gives the overlap. A group that supplies the overlap also records its overlap–Hamiltonian cross-moment, so the variance of the physical Hamiltonian includes the covariance of their shared shots. A diagonal pair prepares its state alone, with one system setting per group.
• Sampled ADAPT, termwise. For comparison with the grouped measurement, acquisition=PencilAcquisition(matrix_measurement="termwise", pair_preparation="controlled") measures the matrix elements with Hadamard tests, each of which controls both complete preparations and the Pauli term on the ancilla and therefore needs pair_preparation="controlled". Each pair uses N shots for every Hadamard test it needs, one test per Pauli term or overlap and per real or imaginary part. FixedGCIM measures grouped joint states only and raises ValueError for this option, and ADAPT planning raises it for this option with shots and bounded Ritz growth.
• ADAPT, joint states that control only the reference. Write \(\phi_i=Q_iB|0\rangle\), where \(B\) is the reference preparation, \(B|0\rangle=|\phi_{\rm ref}\rangle\), \(|0\rangle\) is the all-zero state of the system, the vacuum with no occupied spin orbital, and \(Q_i\) is the generator chain of basis state \(i\), the product of its generator exponentials in the order they act. By default (pair_preparation="prefix") an ADAPT joint state runs \(B\) and the generators at equal angles that both chains start with once without control and controls the remaining generators of each chain on the ancilla, and "controlled" controls both complete preparations. When the exponential of every pool generator maps \(|0\rangle\) to itself up to a known phase, so that \(Q_i|0\rangle=e^{i\gamma_i}|0\rangle\) with the vacuum phase \(\gamma_i\), two further values control only \(B\) and run every generator without control. They adapt the controlled reference preparation of Yoshioka et al., arXiv:2407.14431v3, Fig. 1b and pp. 2–3, to generator chains. "vacuum_reference" prepares the pair of the reference and a state \(\phi_j\) with a Hadamard gate on the ancilla, \(B\) controlled on ancilla value 1, \(Q_j\) and \(B\) controlled on ancilla value 0. The zero branch keeps \(|0\rangle\), with the phase \(e^{i\gamma_j}\), until the second \(B\) prepares the reference on it. Every other pair is prepared as under "prefix". "vacuum_restore" orders the two chains of each pair lexicographically by their (pool index, angle) items, as Python compares tuples, which can reverse their basis order. Write \(Q_L\) for the chain that comes first and \(Q_R\) for the other, with the states \(\phi_L\), \(\phi_R\) and the vacuum phases \(\gamma_L\), \(\gamma_R\). "vacuum_restore" then applies the Hadamard gate, \(B\) controlled on 1, \(Q_R\), \(Q_L^\dagger\), \(B\) controlled on 0 and \(Q_L\). Under \(Q_R\) and \(Q_L^\dagger\) the zero branch stays at \(|0\rangle\) and gains the phase \(e^{i(\gamma_R-\gamma_L)}\), then \(B\) prepares the reference on it and \(Q_L\) makes it \(\phi_L\). The other branch holds \(Q_L^\dagger\phi_R\) until the final \(Q_L\) restores \(\phi_R\), which is why \(Q_L^\dagger\) is the adjoint of the circuit already built for \(Q_L\), its gates inverted in reverse order. When the tuple order reverses the basis order, the pair's entries of \(H\) and \(S\) are conjugated back to basis order. Under both values, phase gates on the ancilla, one per component of the vacuum phases, then give the other branch the phase that the zero branch gained, so that only a global phase remains. They act before any measurement basis change of the ancilla. Both values therefore prepare the same joint state as "controlled", for complex and real-only measurement alike, with the same settings and shots, and \(H\) and \(S\) follow from it as above. Under real-only measurement every vacuum phase is zero and no phase gate is added, because a real generator has neither an identity term nor a term of only I and Z letters. Without cancellations or qubit routing, a pair of the reference costs two controlled \(B\) and one uncontrolled \(Q_j\), and a restored pair also applies \(Q_L\) twice, once as its adjoint. A restored pair whose chains \(Q_L\) and \(Q_R\) hold \(m_L\) and \(m_R\) generators therefore runs \(m_R+2m_L\) uncontrolled generator exponentials, even where "prefix" would run their shared leading generators once for both branches. For example, with \(g_0\) and \(g_1\) the items of pool generators 0 and 1 at one angle, in the basis of the chains \(()\), \((g_0)\), \((g_1)\) and \((g_0,g_1)\), the pair at basis positions 2 and 3, counted from 0, takes \((g_0,g_1)\) as \(Q_L\) and \((g_1)\) as \(Q_R\) and runs 5 exponentials. "vacuum_reference" changes only the \(b-1\) pairs of the reference in a basis of \(b\) states. Whether this costs less or adds less noise than controlling the generators depends on the circuits and the device. ADAPTResult.pair_preparations records the construction of each pair.
• ADAPT, pools that allow controlling only the reference. In quantum execution, exact or with grouped shots, planning accepts "vacuum_reference" and "vacuum_restore" only when it shows, from the Pauli terms of every pool generator after the anti-Hermitian correction that Generator pools describes, that the generator's exponential maps \(|0\rangle\) to itself up to a phase. The whole pool is checked, not only the generators a run selects. Yoshioka et al. base their construction on particle-number conservation, which keeps \(|0\rangle\) up to a phase because the vacuum is the only state with no particle. The library shows directly, with the following vacuum rule at planning and the checks on each built circuit described below, that the complete generator circuits it builds, rounded angles included, keep \(|0\rangle\) up to a phase. It reads no pool name or conservation label, so the rule decides for built-in and supplied pools alike. A generator whose terms commute qualifies when, for each coefficient magnitude and each nonempty set of X and Y positions, its terms with that magnitude and those positions cancel on \(|0\rangle\). For example, the two terms of \(ic(XX+YY)\) map \(|00\rangle\) to \(|11\rangle\) with amplitudes \(1\) and \(-1\), so it qualifies, while \(icX\) does not. Its vacuum phase then comes from its identity term and its terms of only I and Z letters. A double excitation whose terms do not commute, which planning accepts only in the form the built-in pools build, qualifies with phase 0. The vacuum rule is a sufficient condition, so a generator that fails it may still keep \(|0\rangle\), and planning refuses it all the same. Every generator of the four built-in pools passes the vacuum rule with n_spatial_orbitals of 2, 3 and 4 and, for the pools that need an occupation-number reference, a half-filled one. With a failing generator, planning raises ApplicabilityError naming the first such generator by its pool index. "prefix" and "controlled" control the generators and need no vacuum phase. The rule concerns the generators after \(B\), and \(B\) can be any reference preparation. When \(B\) is a supplied circuit, planning counts the synthesis of its controlled form against max_products and max_bytes, as for "controlled". Classical execution and a Plan that measures nothing prepare no joint state, so for them both values are inactive, the vacuum rule is not required and the request is kept. The vacuum phase itself is read from each generator circuit when a query builds it. Before the query runs, the library checks that the angles and global phases of its generator circuits are finite, that half of each rotation angle of a term of only I and Z letters, which enters the vacuum phase, is exactly a binary64 number and, under "vacuum_restore", that the adjoints can be formed from the existing circuit definitions. A query that fails these checks raises ApplicabilityError before it is measured, and the data already obtained stay in the Run. FixedGCIM and QFD have no generator chain from which planning reads a vacuum phase. A FixedGCIM basis state is one supplied preparation, and QFD records and computes no vacuum phase for its time-evolution step. Both therefore raise ApplicabilityError for both values in quantum execution, exact or with shots, unless the Plan measures nothing, which for these two methods happens when \(A\) is a multiple of the identity. The upper bound on appended circuit instructions that ADAPT checks before any query (ADAPT queries) includes the adjoint circuits and phase gates of these constructions. It counts instructions, not the native two-qubit depth that max_native_two_qubit_depth limits for each submitted circuit.
• Sampled, real-only. With shots, the default acquisition=PencilAcquisition(quadratures="auto") of FixedGCIM and ADAPT leaves out the ancilla Y settings and the imaginary Hadamard tests of the off-diagonal pairs when the library proves from the form of the supplied inputs, without forming any state, that \(A\) and every state the method can prepare are real in the computational basis, so that \(H\) and \(S\) are real. The unmeasured imaginary parts are then recorded as exact zeros. The proof concerns the ideal states, and on a noisy device real-only measurement does not show that the imaginary parts are zero. The proof needs every nonzero Pauli coefficient of \(A\) to be real and its Pauli word to contain an even number of Y letters. Every FixedGCIM basis state, and the ADAPT reference, must be an occupation-number state, a tensor product of one-qubit states or a vector, and the supplied amplitudes of the last two must all have zero imaginary part. ADAPT also needs every pool generator to have purely imaginary Pauli coefficients on words with an odd number of Y letters, so that each exponential \(\exp(\theta A_j)\) at a real angle \(\theta\) is a real matrix. A supplied circuit never qualifies. ADAPT measures \(A\) after its pauli_coefficient_cutoff has removed the Pauli terms whose coefficients have magnitude at or below the cutoff, so its proof says nothing about imaginary parts that those terms contribute. quadratures="complex" measures both parts, and quadratures="real" makes planning raise ApplicabilityError when the proof fails. plan.reconstruction.acquisition records the choice in quadratures and the first input that failed the proof in reason. Exact and classical execution compute both parts whatever quadratures requests. A Plan that measures nothing makes no proof, keeps the request unused and raises no error for "real". FixedGCIM, QFD and ADAPT's default growth measure nothing when \(A\) is a multiple of the identity. Bounded Ritz growth on such an \(A\) still measures overlaps whenever its basis can hold more than the reference, so the request applies to it as to any other measured Plan.
• Classical. execution="classical" obtains \(V^\dagger(A-cI)V\) and \(V^\dagger V\) directly from the input, with \(c=\operatorname{tr}(A)/d\), or the identity coefficient of a Pauli input. It removes \(c\) from the Pauli table or from the stored diagonal before any matrix–vector product, and it performs no reference eigensolve or unnecessary Pauli conversion.

Exact quantum FixedGCIM and ADAPT therefore need a target that runs these pair computations, such as Aer's statevector target or NWQ-Sim's CPU/SV target. On other targets, use shots=N or execution="classical". All paths share the projected analysis. Missing data for any real or imaginary part that the Plan measures leaves the pencil incomplete.

For \(A=cI+A_0\), where \(c\) is the identity coefficient, \(I\) the identity and \(A_0\) the remaining operator, the same overlap observations supply \(H=cS+H_0\), with no separately measured identity term and no coordinate-identity substitution. The solve uses \(H_0\) and adds \(c\) to the Ritz values afterwards, and result.pencil records \(H\) (Proposition 14). Inside the solve, \(cS\) would enter the product \(X^\dagger HX\) of the canonical orthogonalization, whose roundoff grows like \(u\lvert c\rvert/s_{\min}\) for the smallest kept overlap eigenvalue \(s_{\min}\), so a large offset on a nearly dependent basis would shift the Ritz values by far more than roundoff. Here \(X\) is the canonical orthogonalizer and \(u\) is unit roundoff. Like classical FixedGCIM, classical ADAPT removes the identity component \(c=\operatorname{tr}(A)/d\) from the Pauli table or from the stored diagonal of a dense or sparse matrix before it projects, where \(d\) is the original dimension of \(A\). The solve adds \(c\) back, so an identity offset of \(A\) does not enter the rounding of the projected entries.

Number of circuits and shots

A complete pencil on \(b\) basis states has \(b(b+1)/2\) pairs \(i\le j\), of which \(b\) are diagonal. Exact quantum execution uses one preparation per pair. For \(G\) nonempty qubit-wise commuting groups, grouped measurement uses \(G\) settings on each diagonal pair and \(2G\) on each off-diagonal pair, \(b^2G\) in total, each with shots shots (fixed_basis.py::sampled_settings). Real-only measurement halves only the off-diagonal settings, to \(b(b+1)G/2\) in total, so three basis states with \(G=2\) need 12 settings under real-only measurement and 18 otherwise. ADAPT measures in each round only the pairs of basis states that no earlier round measured. pencil.matrix_element_count counts the scalar matrix elements before any pair or setting expansion, one termwise Hadamard test each. For \(M\) pairs, \(D\) of them diagonal, and \(L\) nonidentity Pauli terms, it is \(2(M-D)(1+L)+DL\) with both parts measured, because an off-diagonal pair needs the real and imaginary parts of the overlap and of every term, while a diagonal pair needs only the real part of each term. Real-only measurement gives \((M-D)(1+L)+DL\). For a complete pencil the two counts are \(b^2L+b(b-1)\) and \(bL+b(b-1)(L+1)/2\) (pencil.py::matrix_element_count).

Pilot and main shot allocation

FixedGCIM(sampling=PencilSampling(total_shots=12000, pilot_fraction=.2)) replaces shots=N with a budget of total_shots shots that two stages share, and QFD(sampling=...) does the same for its basis of evolved states. Both total_shots and pilot_fraction are required. The basis, the grouped settings and their circuits, and the projected pencil and its solve are those of shots=, and only the number of shots of each setting changes. A pilot stage measures every setting with uniform shots, the solve of the pilot pencil decides the main-stage shots of every setting once, and the main stage then measures every setting again with fresh shots. Only the main-stage measurements form the final \(H\) and \(S\), and with them result.eigenvalue and result.pencil. The pilot measurements stay in the Run, and the pilot shots count toward the Run's shot total and cost like the main ones.

For a pencil with \(m>0\) settings and \(B\) = total_shots, the pilot stage takes \(B_p=\min(B-2m,\max(2m,\lfloor fB\rfloor))\) shots and the main stage the remaining \(B-B_p\), where \(f\) is the exact rational value of the binary64 number pilot_fraction. Each stage can therefore give every setting at least two shots. For FixedGCIM and QFD, \(m\) equals len(plan.experiments), and planning gives \(m\) when it refuses a budget below \(4m\). pilot_fraction sets the pilot budget \(B_p\) through this formula, and the main stage gets \(B-B_p\), whose measurements alone form the final \(H\) and \(S\). The pilot gives every setting \(\lfloor B_p/m\rfloor\) shots and the first \(B_p-m\lfloor B_p/m\rfloor\) settings in the pencil's setting order one more. For FixedGCIM and QFD this is Plan order. ADAPT uses the grouped settings of its final basis, as described below. From the pilot pencil, solved at overlap_cutoff, the method estimates for each setting \(s\) the per-shot variance \(v_s\) of the first-order change that one shot of \(s\) makes in the lowest Ritz value, including the covariance of the overlap and Hamiltonian parts of the same shots. For main-stage shot numbers \(N_s\) that sum to \(B-B_p\), a linear model at the pilot solution predicts the energy variance \(\sum_sv_s/N_s\), and by Cauchy–Schwarz

\[ \sum_s\frac{v_s}{N_s}\ge\frac{\left(\sum_s\sqrt{v_s}\right)^2}{B-B_p}, \]

with equality when \(N_s\) is proportional to \(\sqrt{v_s}\). With at least two shots per setting, the minimum over real-valued \(N_s\) is at \(N_s=\max(2,\sqrt{v_s/\lambda})\), with \(\lambda>0\) chosen so that the \(N_s\) sum to \(B-B_p\). The main-stage shot numbers round this solution deterministically to integers of at least two that spend \(B-B_p\) exactly. Settings at the two-shot floor get two shots, the others get their real-valued \(N_s\) rounded down, and the remaining shots go one each to the settings above the floor where one more shot lowers the local variance \(\sum_sv_s/N_s\) the most, with ties in setting order. These integers need not give the smallest local variance over integer shot numbers. The allocation counts every shot at the same cost, whatever its circuit.

The allocation uses the coefficients and the offset-free Ritz values \(e_0\le e_1\) of the pilot solve (ProjectedPencil.offset_free_eigenvalues), the shifted Ritz values of Read the result, which the solve computes from \(H_0\) and \(S\) before it adds \(c\). For the same pilot measurements, in the same order and with the same budgets, changing only the identity coefficient \(c\) of \(A\) therefore changes no main-stage shot number, provided every physical quantity that the pilot analysis forms stays valid and finite. These include the physical \(H\) and Ritz values of the pilot pencil and, formed before the pilot solve, the variance of each entry of the physical Hamiltonian \(H=H_0+cS\), computed as the physical_h field of a SampledPairVariance record that the allocation does not keep. That variance contains \(c^2\) times the overlap variance, so it can stop being finite while \(H\) and the Ritz values are still finite. If one of these quantities fails this condition, the allocation raises an error, which can come when the variances are formed, before a pilot pencil exists. The module docstring of pencil_sampling.py derives \(v_s\) and states the integer rounding.

The main stage takes the uniform split of \(B-B_p\) instead, without another pilot, in the following cases, and fallback_reason records the first that applies:

• "pilot_solve_failed": the pilot solve fails.
• "pilot_solution_nonfinite": a pilot Ritz value or coefficient is not finite.
• "pilot_lowest_ritz_degenerate": the two lowest offset-free pilot Ritz values are equal.
• "pilot_variance_unusable": a required \(v_s\) is unavailable, not finite or negative.
• "pilot_variances_all_zero": every \(v_s\) is zero, which does not mean that the uncertainty is zero.

A pilot Ritz value outside the sampled enclosure of Read the result is recorded with the pilot pencil and does not change the allocation. A failed measurement follows the Run's failure and resume rules instead.

result.shot_allocation (ShotAllocation) records the two stage budgets, the shots of every setting in the setting order of the pencil, the estimates \(v_s\), the fallback reason, the predicted local variance \(\sum_sv_s/N_s\), the rounded gap \(e_1-e_0\) of the two lowest offset-free pilot Ritz values and the pilot pencil, and print(result) of a FixedGCIM or QFD Result summarizes it in one line. The decision is saved with the Run before the first main-stage measurement. Resuming or loading the Run and calling result.analyze(overlap_cutoff=...) use the saved decision, so a new cutoff changes only the main pencil. An interrupted allocation or final analysis is repeated with run.resume(reanalyze=True), which submits no measurement again.

The local variance describes the linear model at the pilot solution, the gap describes the pilot solve, and neither is a confidence bound. The local variance leaves out the noise of the pilot estimates, a change of kept rank at the cutoff, near-degenerate Ritz values and device bias. No tolerance is applied to the gap, and a nonzero gap does not show that the lowest Ritz value is simple. Given the pilot, the main-stage shot numbers are fixed, so the fresh main-stage measurements support the usual fixed-sample analysis of their moments. The lowest Ritz value of the main pencil is still a nonlinear estimate that is not shown to be unbiased, and a nominal Gaussian interval around it is not shown to have exact coverage. To compare with shots=, match total shots including the pilot. A pilot can lower the main-stage budget enough to change the kept rank, so the allocation alone does not show a more accurate eigenvalue at a fixed total budget.

For FixedGCIM and QFD, planning refuses a budget below \(4m\), two pilot and two main shots per setting. Before any shot it also checks the grouped analysis of the measurements of one stage against max_analysis_work and max_bytes. Because the main-stage shot numbers are not known at planning, the check uses the larger stage budget \(\max(B_p,B-B_p)\) for either stage, so it can refuse a budget whose actual allocation would fit. The message gives the required value. Raise max_analysis_work or max_bytes to that value, or reduce total_shots, the number of basis states or the number of grouped Pauli terms. Before the first pilot circuit is prepared, the Run's ExecutionLimits must also leave room for \(2m\) more circuit preparations and \(2m\) more circuit runs, each counted against max_total_circuits, and for \(B\) more shots within max_total_shots.

Because every setting is measured in a pilot stage and again in a main stage whose shots the pilot decides, nwqlib.estimate(plan) reports the totals of the whole Plan, such as its settings and shots, as unavailable with the reason that plan.method.query_catalog_reason() returns, whereas the shot budget is total_shots and one prepared setting, whose stage and shots are fixed, is estimated normally from its preparation record as for ADAPT queries. print(plan) shows shots=2 for such a FixedGCIM or QFD Plan. This value marks a sampled Plan with the minimum of two shots per setting and stage, and it does not state the shots of the settings. total_shots is the budget of the two stages, and result.shot_allocation holds the shots of every setting. For FixedGCIM and QFD, when \(A\) is a multiple of the identity, no setting is measured, both stages have zero shots, effective_fraction is "not_applicable" and the eigenvalue is the identity coefficient. FixedGCIM and QFD planning refuses sampling together with shots= or with execution="classical". For these two methods, prepare(plan, settings="all") is refused too, because the main-stage shots depend on the pilot, and prepare(plan) prepares the next setting.

ADAPT(sampling=PencilSampling(total_shots=B, pilot_fraction=f)) keeps shots=N for every circuit that grows the basis and adds the pilot and main stages described above once growth has stopped. They measure the basis that growth ended with again, with the separate budget \(B\). For ADAPT, \(m\) is the number of grouped settings of that basis, counted as in Number of circuits and shots. It depends on the basis, not only on \(A\), because bounded growth on a multiple of the identity measures overlap settings once its basis has two or more states. The settings are ordered by the pairs \(i\le j\) of basis states, row by row, then by quadrature and group, which is also the order in which result.group_moments lists the main-stage records. When \(m>0\) and both stages complete, the Run spends \(B\) shots in addition to those of the growth. A Run that ends during the stages counts the shots already submitted, up to \(B\) stage shots. A basis with \(m=0\) adds no shots. result.stop_reason remains the reason that ended growth. After the main stage, only its measurements form result.eigenvalue, result.pencil and result.group_moments, while result.history keeps the analyses of the growth. result.shot_allocation records the decision, and the measurements of the growth and of the pilot stay in the Run. The main-stage measurements are new, but the growth chose the basis from the earlier measurements, so the estimate is conditional on that basis. Without sampling, a run with shots takes the eigenvalue from the same measurements that chose the basis. print(plan) of an ADAPT Plan shows the shots=N of its growth, print(result) shows no allocation line, and prepare(plan) prepares the first query of the growth, if there is one, and never a stage setting.

The stages follow every scientific stop of the growth, invalid_sampled_evidence included, and the main stage can still end without an eigenvalue, for example when its solve fails. A circuit that max_native_two_qubit_depth refuses, or a cancellation with run.cancel(), ends the Run without the remaining stages. The Result then reports the last pencil that the growth solved, if there is one, together with the decision if it was made. The Run counts every shot already submitted, those of an unfinished stage included. If no projected analysis has started, as when a circuit needed for the first pencil is refused, result.eigenvalue is None and result.analysis_attempts is 0. A cancellation during the stages returns at once a Result built from the measurements that the Run had collected, without waiting for a submitted stage job. For such a job, data that the Run has retrieved and saved can be read with run.completed_observation(attempt), including data obtained by a later explicit refresh. These data are not added to that Result. Once the Result exists, run.resume() returns without asking the backend about its jobs again, so a job that completes later enters neither the Run nor the Result by itself.

For ADAPT, planning refuses sampling with execution="classical". When the Plan measures anything, ADAPT supports sampling with a positive shots= and the default grouped measurement, and planning also refuses shots=None and matrix_measurement="termwise". A quantum Plan that measures nothing, such as the default growth for a multiple of the identity, accepts sampling with shots=None or termwise measurement, runs no stage and records zero pilot and main shots. The final basis is known only after the growth, so the checks that planning makes before the first circuit use the largest basis the Plan can reach, 2*min(max_iterations, pool size) states, or min(max_basis_size, max_iterations + 1) states with bounded Ritz growth. Planning requires \(B\ge4m^\ast\) for the \(m^\ast\) settings of that largest basis and checks the grouped analysis of one stage of at most \(B\) shots on it against max_bytes and max_products. A run that stops with a smaller basis may have needed fewer shots, and the stages use the \(m\le m^\ast\) settings of the basis it ends with. A refusal names total_shots, any exceeded limit and the parameters that set the largest basis. Before the first pilot circuit, the Run's ExecutionLimits must leave room for \(2m\) more circuit preparations, \(2m\) more circuit runs and \(B\) more shots. The Run checks this room once growth has stopped, so a Run without it is refused before the first stage shot, with the shots of the growth spent, and continues after run.extend_limits(...) and run.resume().

When the overlap matrix is refused

A valid rank-deficient trial space is distinguished from an invalid deterministic Gram matrix. The solve refuses a deterministic Gram matrix only when an overlap eigenvalue lies below minus the sum of the derived error of its entries and the eigensolver roundoff. For a sampled overlap, it instead solves on the span of overlap eigenvectors whose eigenvalues are strictly greater than the positive overlap_cutoff, without rejecting the matrix because other overlap eigenvalues are negative (canonical orthogonalization). The entry error is as follows, with unit roundoff \(u=2^{-53}\) and \(\gamma_n=nu/(1-nu)\) (Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., doi:10.1137/1.9780898718027, Lemma 3.1).

• Classical entries. The error comes from the inner products of length \(d\), \(\sqrt2\,\gamma_{2d}\operatorname{tr}(S)\) to first order.
• Exact quantum FixedGCIM and ADAPT. The error is the largest row sum of the off-diagonal overlap entry bounds, each from the saved-state error (saved_state_error) of its own pair's preparation record (Engineering constants). Only off-diagonal overlap bounds enter this Gram error budget. An unresolved preparation record leaves the bound unavailable, and the solve then allows no input error beyond its own roundoff.

Exact FixedGCIM pairs use X gates in the circuit for the ancilla flips. Their preparation records supply a first-order state-error budget when the executed circuit operations and parameters satisfy the backend's roundoff model, and the pair computation accounts for the classical contractions. Supplied matrix gates are synthesized when controlled for off-diagonal pairs, whose finite bounds assume the library's exact-to-rounding synthesis convention. An uncontrolled matrix in a diagonal pair can leave that pair's bound unavailable, and unresolved roundoff exclusions of gate parameters, such as an excessive U-gate phase sum, can also leave off-diagonal bounds unavailable.

Choose the measurement and the basis

The first seven rows below change only how the matrix elements of a fixed basis are measured. They keep the ideal \(H\) and \(S\), so with exact readout they give the same Ritz values up to rounding, while the circuits, the settings, the shots and the finite-shot covariance differ. PencilSampling keeps the basis and the settings and chooses the shots of each setting from a pilot stage. Bounded Ritz growth builds a different basis and therefore solves a different projected problem.

Input or construction	Recommended or default	Evidence
Real \(A\), occupation, product or real-vector states and, for ADAPT, generators whose exponentials are real at real angles, as for chemistry with the built-in pools	quadratures="auto", the default, which resolves to real-only measurement (Sampled, real-only)	E1
Complex amplitudes, a supplied circuit or a QFD basis	quadratures="auto", which resolves to complex measurement	E2
Basis states whose preparations begin with the same blocks, such as nested ADAPT chains and their common reference	pair_preparation="prefix", the default	E3
Supplied states with no common leading block	The same "prefix" request, which then controls both preparations of the pair	E4
A Pauli Hamiltonian with fewer qubit-wise commuting groups than terms	matrix_measurement="grouped", the default	E5
ADAPT with an occupation reference and a pool that passes the vacuum rule, when controlled generators are costly	Try pair_preparation="vacuum_reference"	E3
ADAPT with shallow chains and a pool that passes the vacuum rule	Compare pair_preparation="vacuum_restore" on your circuits	E3
Settings whose contributions to the energy variance differ widely, with enough shots for a pilot stage	Try sampling=PencilSampling(...)	E6
ADAPT whose product chain would exceed the circuit depth you can run	basis_growth=BoundedRitzGrowth(...), which changes the basis	E7

Each row rests on one comparison on one input of one or two system qubits, run with NWQLib 1.0.3.post2, Python 3.12.14, Qiskit 2.6.0 and Aer 0.18.0 on macOS arm64. A sampled comparison is one run with seed 7, so its energies carry a sampling error that one run does not measure, and they rank no option. No comparison used a noisy device or a noise model, so none shows a smaller device error. Gate counts and depths are those of the circuits compiled for Aer's statevector target, and qubit routing and the gate set of another target can change them.

• E1. FixedGCIM on \(A=0.5\,ZZ+ZI+0.7\,IZ+0.3\,XX+0.2\,XI\) of Bounded Ritz growth, whose terms the library places in \(G=3\) groups, with the basis ingest_occupation(bits, num_qubits=2) for bits (0, 0), (0, 1) and (1, 1) and 4096 shots per setting on Aer's counts simulator. "auto" resolves to real-only measurement with 18 settings and 73,728 shots, and quadratures="complex" measures 27 settings with 110,592 shots, as \(b(b+1)G/2\) and \(b^2G\) give for \(b=3\). The two sampled values, −1.2320 and −1.2231, lie on either side of the value −1.2270 that exact readout gives for the same basis. Real-only measurement does not show that the imaginary parts are zero on a noisy device.
• E2. FixedGCIM on \(A=\operatorname{diag}(1,-1)\) with the basis \((1,1)/\sqrt2\), \((1,i)/\sqrt2\) of the first example and 1000 shots. Planning resolves complex measurement and names, in reason, the vector with a complex amplitude. Here \(H_{01}=(1-i)/2\) has the imaginary part \(-1/2\), which real-only measurement would record as zero. QFD always measures both parts (QFD).

• E3. ADAPT with exact readout on \(A=0.5\,ZI+0.25\,XX+0.375\,IY\) from the occupation reference with qubit 0 occupied, with the pool \(i(0.5\,II+0.5\,XX+0.5\,YY+0.25\,ZZ)\) and \(0.5\,i(XY-YX)\), both of which pass the vacuum rule. With \(g_1\) and \(g_2\) generator 0 at the angles 0.3 and 0.9 and \(h\) generator 1 at 0.6, the table gives the two-qubit gates and the depth of the compiled circuit of four pairs, its final save instruction included. The two-qubit depth equals the two-qubit gate count in every cell. The counts include the phase gates on the ancilla, and they show no inequality for other circuits or devices.

  Pair of chains	Construction	Two-qubit gates and depth
  \(()\), \((g_1)\)	"prefix"	9, depth 16
  \(()\), \((g_1)\)	"controlled"	11, depth 20
  \(()\), \((g_1)\)	"vacuum_reference"	5, depth 9
  \(()\), \((g_1)\)	"vacuum_restore"	5, depth 9
  \(()\), \((g_1,h)\)	"prefix"	15, depth 29
  \(()\), \((g_1,h)\)	"controlled"	17, depth 33
  \(()\), \((g_1,h)\)	"vacuum_reference"	9, depth 21
  \(()\), \((g_1,h)\)	"vacuum_restore"	9, depth 21
  \((g_1)\), \((g_2,h)\)	"prefix"	24, depth 44
  \((g_1)\), \((g_2,h)\)	"controlled"	26, depth 47
  \((g_1)\), \((g_2,h)\)	"vacuum_reference"	24, depth 44, as "prefix"
  \((g_1)\), \((g_2,h)\)	"vacuum_restore"	15, depth 27
  \((g_1,h)\), \((g_2)\)	"prefix"	24, depth 44
  \((g_1,h)\), \((g_2)\)	"controlled"	26, depth 46
  \((g_1,h)\), \((g_2)\)	"vacuum_reference"	24, depth 44, as "prefix"
  \((g_1,h)\), \((g_2)\)	"vacuum_restore"	19, depth 38

• E4. FixedGCIM on the \(A\) of E1 with the basis ingest_occupation((1, 0), num_qubits=2), the vector \((0.6,0,0,0.8)\) and the same occupation state again, with 64 shots and the default "prefix". plan.reconstruction.acquisition.pairs records the "prefix" route with one shared block for the pair of the two identical occupation states and the "controlled" route for the other two pairs. Each supplied state is one block, so only identical preparations share a prefix.

• E5. ADAPT on \(A=ZI+IZ+0.5\,XX\), with \(L=3\) terms in \(G=2\) groups, from \(|00\rangle\) with the pool \(iXY\), \(iYX\), \(\theta=\pi/8\), one selection and 256 shots per setting on Aer's counts simulator. Both measurements select the same generator and end with a basis of two states. Grouped measurement used 9 settings and 2,304 shots in all, and termwise measurement 13 settings and 3,328 shots. On the final basis alone, real-only grouped measurement needs \(b(b+1)G/2=6\) settings and termwise measurement \(bL+b(b-1)(L+1)/2=10\) Hadamard tests. The sampled values, −1.117 with grouped and −2.468 with termwise measurement, are both far from the exact-readout value −2.0616 of the same method, and the termwise value lies below the lowest eigenvalue of \(A\), as a sampled pencil can. One run at this shot count compares no accuracy. Grouped measurement estimates the terms of one group from the same shots, so their covariances enter the energy variance, and fewer settings do not show a smaller error under every noise model.
• E6. FixedGCIM on the input of E1 with PencilSampling(total_shots=12000, pilot_fraction=0.2). The pilot stage gives its 2,400 shots to the 18 settings as 133 or 134 each, and the main stage gives each setting between 2 and 2,316 of its 9,600 shots, with no fallback. The two stages ran 36 circuits, one per setting and stage. Uniform measurement of the same 18 settings with 666 shots each, 11,988 shots in all, ran 18 circuits. The two values, −1.2329 and −1.2348, against −1.2270 from exact readout of the same basis, come from one run each and show no difference in accuracy. The pilot shots count toward the total, and the local variance that the allocation predicts, 2.7e-4 in squared energy units, is not a confidence bound.
• E7. The example of Bounded Ritz growth, with exact readout on Aer, against the default growth with the same pool, \(\theta=\pi/8\) and max_iterations=3. Bounded growth ends with four chains whose CX count bounds are at most 2, and gives −1.2769679189111813. The default growth selects all three generators and ends with six basis states, and gives −1.2769679189111725. Both runs prepared 21 queries. Both values equal the lowest eigenvalue −1.2769679 of \(A\) up to roundoff, because both bases span the two-qubit space. The default product chain here also needs at most two CX gates, so the example shows the two bases and their cost, not a depth saving. On a larger problem the two methods solve different projected problems and can give different energies.

Adaptive generator coordinates

from nwqlib import Eigenproblem, solve
from nwqlib.algorithms.gcim import ADAPT
from nwqlib.operators import ingest_pauli

problem = Eigenproblem(A=[[1, 0], [0, -1]])
pool = (ingest_pauli((("Y", 1j),), num_qubits=1),)
method = ADAPT(initial_state=[1, 1], pool=pool)
result = solve(problem, method=method, seed=7)
print(result.selected, result.stop_reason, result.eigenvalue)

(0,) pool_exhausted -0.999999999999999

ADAPT screens the current product state with \(\langle[A,A_j]\rangle\), chooses the largest absolute gradient, then expands the generator-coordinate basis. Here the full operator is \(A=Z\) and the sole pool generator is \(A_j=iY\), where X, Y and Z are the one-qubit Pauli matrices, giving \([Z,iY]=2X\), so the initial \(|+\rangle\) gradient is 2, and the expanded basis reaches the eigenvalue \(-1\) before the one-generator pool is exhausted. The gradients are evaluated on the product state, not on the projected Ritz state. For \(k\ge2\) selections, write \(U_\ell=\exp(\theta_\ell A_{j_\ell})\) using the selected indices and angles defined in the opening paragraph. The ordered basis is the reference \(|\phi_{\rm ref}\rangle\), the \(k\) individual states \(U_\ell|\phi_{\rm ref}\rangle\) in selection order, the full product \(U_k\cdots U_1|\phi_{\rm ref}\rangle\), then the prefixes \(U_r\cdots U_1|\phi_{\rm ref}\rangle\) for \(r=2,\ldots,k-1\) in increasing order. This gives \(1+k+1+(k-2)=2k\) entries. At \(k=1\) the basis contains the reference and \(U_1|\phi_{\rm ref}\rangle\), and at \(k=0\) it contains only the reference. The count is the number of listed states, not their rank. ADAPT obtains every number it needs through queries. A query is one evaluation, on the backend or by the classical routine, at fixed generators and angles, for example the screening gradients of one product state, matrix elements of one pair of basis states or one product-state energy. With shots, a matrix-element query covers one grouped measurement setting, for one pair of basis states, one group of Hamiltonian terms and one real or imaginary part. Under termwise measurement it covers one Hadamard test of a pair, for one Pauli term or the overlap and one real or imaginary part. With shots, a screening or energy query covers one measurement group of its labels.

Pool order decides only between absolute gradients that are equal as binary64 numbers, and the lowest index wins. Gradients that are equal in exact arithmetic, for example because a symmetry of \(A\), the reference and the pool maps one generator to another, usually differ in binary64 by roundoff, so the summation order, such as the order of the Hamiltonian's Pauli terms, can decide which of them wins. When a symmetry relates the tied generators, each choice gives the same energies in exact arithmetic. For example, take six qubits on a ring, \(A=\sum_j(h_1Z_j+h_2X_jX_{j+1}+h_3X_j)\) with \(h_1=\sqrt2/3\), \(h_2=\pi/7\) and \(h_3=e/5\) on every site, the product reference \(\cos(0.3)|0\rangle+\sin(0.3)|1\rangle\) on every qubit and the pool \(iY_j\). The six gradients agree to 50 digits. With the Pauli terms passed through ingest_pauli in the order of this sum for \(j=0,\ldots,5\), and qubit 0 as the rightmost character of each label, classical screening gives five different binary64 values and selects \(iY_5\). Reversing the order of the terms gives four values and selects \(iY_3\) (Measured evidence).

Defaults and stopping

Defaults are theta=pi/4, gradient_norm_floor=1e-8, energy_change_tolerance=1e-6, overlap_cutoff=1e-12, stop_criterion="flat_counter" and max_iterations=8. With execution="classical", stop_criterion="residual_norm" instead tests whether the normalized Ritz state \(\psi\) of the current basis satisfies \(\lVert(A-EI)\psi\rVert_2\le\texttt{residual_norm_tolerance}\) after at least one selection. The default threshold is 1.6e-3 in the Hamiltonian's energy unit. These are operational settings, not guarantees of convergence or of ground-state coverage. With a fixed angle, the selected generators and the energy after a given iteration depend on the sign of each generator. Replacing a generator \(A_j\) by \(-A_j\) replaces the factor \(\exp(\theta A_j)\) by \(\exp(-\theta A_j)\) in the product state defined in the opening paragraph, where \(\theta\) is the fixed angle, so the spanned subspace and the later screening can differ.

The default flat-counter stopping rule follows Zheng et al., arXiv:2312.07691v3, Methods. It compares consecutive projected-energy changes with energy_change_tolerance and stops after min(t_auto_fraction * unselected_count, t_user) flat rounds. The defaults are t_auto_fraction=.2 and t_user=10. The threshold \(T\) is recomputed after each selection from the current number of unselected generators, and because the count is an integer, the rule stops as soon as the count is at least \(\lceil T\rceil\). Counting starts with the reference-to-first-basis energy change and resets after a large or unavailable change. This heuristic does not prove a bound on the physical energy error. Replacing it with better-supported stopping evidence is planned research work.

With both default limits, max_iterations=8 and t_user=10, the flat counter can fire only if the original pool has at most 42 generators. After k selections a pool of P generators has P-k unselected members and a flat count of at most k, so the stop needs k >= .2(P-k). The iteration cap is checked first and ends the run at step 8, which makes step 7 the last step for this stop and gives P <= 42. Larger pools reach the iteration cap before this particular stop can fire. This does not affect the other stopping conditions.

The paper's accompanying QuGCM implementation, adapt_gcim.py, lines 303–398, selects each pool index at most once, and NWQLib follows this ADAPT-GCiM selection convention. Reusing a fixed-angle generator would repeat its single-generator reference basis vector, although its cumulative product vector could be new, so allowing reuse would change the trial-space construction and the pool-exhaustion rule. Grimsley's ADAPT-VQE step 7 (arXiv:1812.11173v2) addresses a different ansatz with optimized parameters. The same QuGCM source uses max(10, floor(.1*unselected)) and an initial warmup for stopping, which differs from Zheng's published Methods rule (arXiv:2312.07691v3). NWQLib's flat counter follows the published rule above.

Generator pools

Pass a named pool to pool as one of the strings "spin_adapted_sd", "uccsd_sd", "qeb_sd" or "ceo_ovp" (the OVP-CEO pool). A named pool needs n_spatial_orbitals, and pools other than "spin_adapted_sd" need an occupation-number reference state. The named pools keep their established ordering and preparation restrictions. The Pauli coefficient cutoffs and any anti-Hermitian correction of a pool generator within input_symmetry_tolerance are recorded with the size of what they removed and the identifier of the original input.

The default spin-adapted pool uses the signed operators of Zheng et al., arXiv:2312.07691v3, Appendix E.1, Eqs. (E1)–(E3). It enumerates spatial pairs with p <= q, r <= s and the first pair no later than the second in lexicographic order. The double-excitation builders write interleaved ladder operators, so their signs follow from anticommutation, adjoint subtraction and positive coefficient normalization together. NWQLib's tests compare every operator in Table V (Table 5 in the npj version) with the pool, including signs and repeated indices. Eq. (E4)'s q < p convention describes elementary circuit decomposition and does not prescribe the spin-adapted pool's enumeration.

Bounded Ritz growth

ADAPT(basis_growth=BoundedRitzGrowth(max_chain_length=L)) replaces the gradient screening and the basis described above with a basis of short generator products that grows by one product per round, chosen by a score computed from the current Ritz pair. Every basis state then needs at most \(L\) generator exponentials after the reference, basis_depth_limit can also bound their CX count, and the Result records what each round examined and what it obtained. Here \(A=0.5\,ZZ+ZI+0.7\,IZ+0.3\,XX+0.2\,XI\), with qubit 0 as the rightmost character of each label, the reference is \(|00\rangle\) and the pool holds \(iYI\), \(iIY\) and \(iXY\) with indices 0, 1 and 2:

from math import pi

from nwqlib import Eigenproblem, solve
from nwqlib.algorithms.gcim import ADAPT, BoundedRitzGrowth
from nwqlib.operators import ingest_pauli
from nwqlib.problems import ingest_occupation

terms = (("ZZ", .5), ("ZI", 1.0), ("IZ", .7), ("XX", .3), ("XI", .2))
pool = tuple(ingest_pauli(((p, 1j),), num_qubits=2) for p in ("YI", "IY", "XY"))
growth = BoundedRitzGrowth(max_chain_length=3, basis_depth_limit=2)
method = ADAPT(initial_state=ingest_occupation((0, 0), num_qubits=2), pool=pool, theta=pi / 8,
               max_iterations=3, basis_growth=growth)
result = solve(Eigenproblem(A=ingest_pauli(terms, num_qubits=2)), method=method, seed=7)
print([[index for index, _ in chain] for chain in result.basis])
for row in result.history[:-1]:
    print([([index for index, _ in c.chain], c.reason, c.depth_bound) for c in row.candidates])
print(result.stop_reason, result.eigenvalue)

[[], [2], [2, 1], [2, 1, 0]]
[([0], 'ranked', 0), ([1], 'ranked', 0), ([2], 'accepted', 2)]
[([2, 0], 'ranked', 2), ([2, 1], 'accepted', 2), ([2, 2], 'above_depth_limit', 4)]
[([2, 1, 0], 'accepted', 2), ([2, 1, 1], 'below_independence_cutoff', 2), ([2, 1, 2], 'above_depth_limit', 4)]
max_iterations -1.2769679189111813

Each basis state is a chain \(\exp(\theta A_{g_k})\cdots\exp(\theta A_{g_1})|\phi_{\rm ref}\rangle\) of \(k\le L\) pool generators at the fixed angle \(\theta\), and the reference is the empty chain. result.basis lists the chains in basis order, each as a tuple of (pool index, theta) items, the generator that acts first listed first, and the example prints only their pool indices. A round extends one chain of the basis, its parent, by one more generator, which acts last, and the resulting chain is an extension of that parent. The parents of a round are the reference and the parent_window - 1 chains added last (default parent_window=4), except chains that already have \(L\) generators. The round visits the parents in turn, repeatedly, in order of basis position, and each visit extends the parent by its next pool generator, so each parent continues through the pool in pool order from where it stopped. The turn also continues across rounds, because a round starts with the first parent after the one whose extension the previous round examined last, and wraps to the reference when no later parent remains. The examination stops after candidates_per_round (default 4) extensions that are new to the basis and within the depth limit, after max_candidate_scans (default 256) examined extensions or when no parent has a generator left. The round then evaluates each of these extensions \(v\). It obtains the matrix elements of \(v\) with the \(b\) current basis states and its diagonal element and computes the score below, then adds the extension with the largest score, the earliest examined one among equal scores, and the enlarged basis is solved with ADAPT's overlap_cutoff.

A generator can occur several times in a chain, also twice in a row, because \(\exp(\theta A_g)\exp(\theta A_g)=\exp(2\theta A_g)\) is in general a different state from \(\exp(\theta A_g)\), and each occurrence counts toward \(L\) and toward the CX count. For a generator with \(A_g^2=-I\), such as \(i\) times one Pauli word, like \(iXY\) in the example, \(\exp(\theta A_g)^2=2\cos\theta\,\exp(\theta A_g)-I\). So if a basis chain \(\beta\) extends a basis chain \(\alpha\) by \(g\), repeating \(g\) gives the state \(2\cos\theta\,|\beta\rangle-|\alpha\rangle\), which adds no direction to the span of the basis, and its part outside the kept overlap subspace is zero whenever that subspace contains both states. With candidates_per_round=1 and without shots, a round whose only evaluated extension does not pass independence_cutoff with a finite score ends the growth with no_admissible_candidate. With shots, a round whose only evaluated extension fails the sampled eligibility rule ends the growth with noise_inconclusive, and a round whose only evaluated extension passes the rule without a finite score ends the growth with no_admissible_candidate. A new parent starts at pool generator 0, so after a first added chain of generator 0 the second round examines this repeat first only when that round is reached, max_chain_length is at least 2 and parent_window is at least 2, because a chain that already has max_chain_length generators, or one outside the window, is not a parent, and it evaluates the repeat only when the repeat is within the depth limit. A larger candidates_per_round lets a round evaluate other extensions as well, but does not guarantee that one of them passes independence_cutoff, or with shots the sampled eligibility rule, with a finite score. An extension equal to a chain of the basis is not examined further. The default basis holds the reference, the single-generator states and the prefixes of one product, while a bounded basis can hold products that start from different parents, so its states need not be prefixes of one another.

In the example, \(\exp(\theta\,iYI)\) and \(\exp(\theta\,iIY)\) are one-qubit rotations without CX gates and \(\exp(\theta\,iXY)\) has a decomposition into one-qubit gates and two CX gates, so basis_depth_limit=2 excludes every chain with two \(iXY\) factors, whose CX count bound is 4. The first round examines the three single generators, all from the reference. The reference then has no generator left, so the second round extends only [2] and the third only [2, 1]. The adjacent repeat [2, 1, 1] is the state \(\exp(2\theta\,iIY)\exp(\theta\,iXY)|00\rangle\). It lies in the span of the three basis states in exact arithmetic, because the rotations \(\exp(\varphi\,iIY)\) of \(\exp(\theta\,iXY)|00\rangle\) stay in one plane, so its part outside that span is at rounding level and it is not scored. The four chains span the two-qubit space, and the final value is the lowest eigenvalue of \(A\) up to roundoff.

The score uses the current lowest Ritz value \(E\) and Ritz state \(\psi\). Let \(q\) be the normalized part of \(v\) outside the kept overlap subspace, the span of the overlap directions that the current solve keeps above its cutoff, and \(\eta\) the squared norm of that part before normalization. With \(\delta=\langle q|A|q\rangle-E\) and \(\rho=\lvert\langle q|A-E|\psi\rangle\rvert\), the score \((\sqrt{\delta^2+4\rho^2}-\delta)/2\) is the amount by which the lowest eigenvalue of \(A\) projected on \(\psi\) and \(q\) lies below \(E\). The kept subspace together with \(q\) contains both states, so by Rayleigh–Ritz its lowest Ritz value is at or below that eigenvalue, and the score is a lower bound on the decrease that an exact solve of the current kept subspace together with \(q\) would give. It does not rank the extensions by the energy that a solve of each enlarged basis would give. The solve that follows can keep a different overlap subspace, so the decrease it gives, recorded as energy_decrease in the next entry of result.history, can differ from the score and is negative when the energy rose. A score of zero does not show that an extension cannot lower the energy, because an extension coupled only to an excited Ritz state has \(\rho=0\) while adding it can lower the lowest Ritz value. The score is formed from \(A_0=A-cI\) and the Ritz value before \(c\) is added, as in the solve (How matrix elements are obtained), so a large identity coefficient cannot round it away. An extension is scored only when \(\eta\) exceeds independence_cutoff, whose default is ADAPT's overlap_cutoff. This floor does not guard against sampling noise, so with shots the rule below decides which extensions are scored.

With shots, the matrix elements of the extension and of the current basis are sample means, and the division by \(\sqrt\eta\) and \(\eta\) in the score can enlarge their noise without bound, also when the sampled \(\eta\) exceeds independence_cutoff. Each evaluated extension is then scored only when it passes a sampled eligibility rule, which uses only overlap settings already measured, among the current basis states and between them and the extension. Let \(T\) map coordinates of the kept overlap subspace to the basis, with \(T^\dagger ST=I\), let \(s_i=\langle\phi_i|v\rangle\) be the sampled overlaps of the extension \(v\) with the basis states, \(t=T^\dagger s\) and \(a=Tt\), and let \(\hat\eta=1-t^\dagger t\) be the sampled squared norm of its part outside the kept subspace, which can be negative and is not clipped. Each measured real or imaginary part of an overlap, with sample-mean variance \(V\) from the \(N\) shots of its setting, enters as \(\tilde V=\max\{V,(N+1)^{-2}\}\). Here \(1/(N+1)\) is a chosen scale for the resolution of one count among outcomes in \([-1,1]\), and its square keeps a part whose \(N\) outcomes were all equal from counting as noise-free. It is not a lower bound on the true variance. The diagonal entries \(S_{ii}=1\) and the imaginary parts that real-only measurement leaves out are exact and contribute nothing. With the fixed multiplier \(k=3\), the rule forms the following quantities:

• The first-order spread \(\sigma\) of \(\hat\eta\), from the first-order change \(a^\dagger(\delta S)a-2\operatorname{Re}(a^\dagger\delta s)\) of \(\eta\) under errors \(\delta S\) of the current overlap matrix and \(\delta s\) of the extension's overlaps. \(\sigma^2\) sums, over the measured parts, the squared derivative of this change with respect to the part times its \(\tilde V\).
• The ratio \(\epsilon=k\sqrt{\nu_S}/\lambda_*\), where \(\nu_S\) is the largest, over the rows of \(S\), of the sum of \(\tilde V\) over the measured parts of the row's off-diagonal entries, and \(\lambda_*\) the smallest overlap eigenvalue that the current solve keeps.
• \(d=k\sqrt{b_s}\), where \(b_s=\sum_i\lVert T_{i,:}\rVert^2(\tilde V_{i,\mathrm{re}}+\tilde V_{i,\mathrm{im}})\) collects the extension's overlap terms, and \(x=\lVert t\rVert\).

An extension passes the rule when \(\epsilon<1/2\) and the margin \(\hat\eta-R\) exceeds independence_cutoff, with the resolution \(R=k\sigma+(d+\epsilon x)^2/(1-\epsilon)\). For \(\epsilon<1\), the second term of \(R\) bounds the change of \(\eta\) beyond first order whenever \(\lVert T^\dagger(\delta S)T\rVert\le\epsilon\) and \(\lVert T^\dagger\delta s\rVert\le d\), and it accounts for the noise at \(t=0\), where the first-order spread vanishes. The rule computes \(\epsilon\) and \(d\) from variance estimates rather than from such norm bounds and assumes independently measured settings, so it is an empirical rule, not a confidence interval or a proof that the extension is linearly independent. It does not cover the rotation of the kept subspace by the noise of \(S\), the choice of the extension from earlier data, or device bias and drift. An extension that fails the rule has the reason screening_noise_unknown when a setting has fewer than two shots or no finite nonnegative variance, or when one set of counts serves two different overlap entries, screening_metric_unresolved when \(\epsilon\ge1/2\) or a needed quantity is not finite, and screening_independence_unresolved when the margin is at most independence_cutoff, which does not show that the extension depends on the basis. The rule adds no shot or repeated measurement and no full-state computation, and its work is of order \(b^2\) for each extension. It does not control the noise of the Hamiltonian elements or of the current Ritz pair, so the score of an extension that passes is a raw estimate, not a validated energy decrease, and the round ranks these extensions by it as with exact elements.

More basis states add terms to the row sums \(\nu_S\), and a small \(\lambda_*\) enlarges \(\epsilon\). As a rough planning estimate, not a recommended shot count, real-only overlaps with \(\tilde V\approx1/N\) give \(\epsilon<1/2\) only when each setting has \(N>36(b-1)/\lambda_*^2\) shots, which for \(\lambda_*\approx1\) is \(N>540\) at \(b=16\) and \(N>1116\) at \(b=32\), and measuring both parts doubles these numbers. Each extension's margin still decides whether it is scored. When no evaluated extension of a round passes the rule, growth stops with noise_inconclusive on the current basis. Without sampling the run then reports the energy of that basis's last solve and adds no shots. With ADAPT.sampling the pilot and main stages of Pilot and main shot allocation measure that basis again, and a completed main stage gives the reported energy.

Each evaluated extension needs \(b\) off-diagonal pairs and one diagonal pair, so with \(b\) basis states a round obtains the matrix elements of at most candidates_per_round * (b + 1) new pairs of states. In exact quantum execution each pair is one query, and with shots each extension takes \(G(b+1)\) grouped settings for \(G\) measurement groups under real-only measurement and \(G(2b+1)\) otherwise (Number of circuits and shots). The matrix stage of the next round finds the elements of the added chain already obtained. The other evaluated extensions are not used again, while their work stays counted in the Run and their observations stay in its data. In the example the three rounds obtain 7, 6 and 8 pairs, the 7 including the reference's diagonal element. In quantum execution, before a round obtains the elements of its evaluated extensions, the Run checks the number of their new queries against its preparation and circuit limits, with shots also shots times that number against its shot limit, and in exact execution the work of their pair reductions against max_products, as it does for a matrix stage. A round above those limits is rejected before its first query, and after run.extend_limits(...), run.resume() continues it. run.extend_limits(...) raises the Run's limits, while ADAPT.max_products is a field of the Method that this call does not change, and the resumed round checks each of its exact pair reductions against max_products as described next. A resumed round may already be partly run, so, as for a resumed matrix stage, the Run then skips the check of the whole round and checks each query when it runs (Continue and reanalyze). run.resume() without raising a limit therefore runs the round's queries until one of them exceeds a Run limit, and each exact pair reduction, not the round's total, is checked against max_products. Without basis_depth_limit, examining an extension only compares chains. With it, the uncontrolled circuit of each examined generator is built to read its CX count, and the count is kept for later rounds. While the round obtains the matrix elements of its evaluated extensions without interruption, their queries reuse those uncontrolled circuits and build the controlled and inverse forms separately. Bounded growth screens no gradients and forms no commutators \([A,A_j]\), so gradient_norm_floor has no effect.

basis_depth_limit=D excludes every extension whose CX count bound exceeds \(D\), and a bound equal to \(D\) is allowed. The bound adds the CX counts of the reference preparation and of every generator occurrence of the chain, each decomposed into one-qubit gates and CX before routing, with no cancellation between them, so it bounds the CX depth of the chain's preparation. An occupation or product reference counts 0, and every direct vector on \(n\) qubits counts \(2(2^n-2)\), a computational basis state given as a vector included, because the amplitudes are not read. Planning raises ValueError for a supplied or matrix-product-state reference, whose count is unknown, for a reference whose count exceeds \(D\), and for classical execution. The counts are read from the generator circuits that quantum execution builds, so the limit is supported for quantum execution only, and classical bounded growth without a limit is unaffected. An extension whose generator circuit contains an instruction without a known CX count is excluded like one above \(D\). The limit changes the basis that ADAPT can reach. It is not the native two-qubit depth of the executed circuits, whose joint-state preparation and routing add depth, and which ADAPT.max_native_two_qubit_depth limits.

The following operational limits end bounded growth. None of them shows that the energy has converged or reached a requested accuracy, and a small score among the examined extensions says nothing about the extensions that no round examined.

• max_basis_size or max_iterations: the basis reached max_basis_size states or max_iterations rounds each added a chain. A round that adds no chain ends the run with one of the reasons below, so the basis holds at most min(max_basis_size, max_iterations + 1) states.
• flat_counter: the energy changes stayed below energy_change_tolerance as in Defaults and stopping, with the pool size as a fixed scale in place of the number of unselected generators.
• basis_depth_limit: the round examined at least one extension that is not already in the basis, and every such extension exceeds basis_depth_limit or has no CX count.
• parent_window_exhausted: no parent of the round has a generator left, and every extension the round examined was already in the basis.
• candidate_scans_exhausted: the round examined max_candidate_scans extensions, all already in the basis.
• no_admissible_candidate: without shots, no evaluated extension has \(\eta\) above independence_cutoff and a finite score. With shots, at least one evaluated extension passed the sampled eligibility rule and none of these has a finite score.
• noise_inconclusive: with shots, no evaluated extension passed the sampled eligibility rule.

When a round finds no extension to evaluate, basis_depth_limit takes precedence over parent_window_exhausted, and that over candidate_scans_exhausted, and each of these three reasons describes only the extensions the round examined. Without shots, a solve that keeps no overlap direction also stops the run with no_usable_overlap_subspace. With shots, a solve that the eigensolver rejects, including one that keeps no overlap direction, stops it with invalid_sampled_evidence. Execution states of result.stop_reason such as pending_acquisition, cancelled and native_depth_refused remain possible.

result.history[k].candidates lists the extensions that the round after analysis \(k\) examined, in examination order, each with its chain, reason, depth_bound, eta (\(\eta\), or the unclipped \(\hat\eta\) with shots), coupling (\(\rho\)), gap (\(\delta\)) and score, and with shots also the quantities of the eligibility rule, such as eta_sigma (\(\sigma\)), metric_epsilon (\(\epsilon\)), eta_resolution (\(R\)) and eta_margin. The added extension has reason "accepted". history[k].energy_decrease and history[k].kept_rank give the decrease and the kept overlap rank of each solve, without clipping a rise. result.frontier_cursor holds, for each parent of the last round, its basis position and the next pool index it would take, ordered so that the parent whose extension the round examined last comes last, and the next round starts after it. Saving, loading, run.resume() and reanalysis at another cutoff keep the chains, the cursor and these records, examine no completed round again and obtain no matrix element twice. result.verify raises ValueError for residual and sector checks, which rebuild the default basis, while supplied-reference-energy checks and ProjectedVerificationOptions apply to a bounded basis.

ADAPT raises ValueError for BoundedRitzGrowth together with product-state optimization or stop_criterion="residual_norm", because bounded growth keeps every angle at \(\theta\) and computes no full-space residual. With shots, bounded growth is implemented for the default grouped measurement, whose first group's settings supply the overlap variances of the eligibility rule, and planning raises ValueError when termwise measurement, acquisition=PencilAcquisition(matrix_measurement="termwise", pair_preparation="controlled"), is requested. Use the grouped measurement, or the default growth basis_growth=None for a termwise comparison.

Optional product-state optimization

Optional optimization requires a period, an iteration count and a finite limit on total energy evaluations, which counts distinct evaluated parameter points. Quantum optimization uses the bounded generalized or Pauli-insertion shift rules (Proposition 15). Classical product optimization evaluates energy and gradient together with a forward and backward sweep, sharing one Hamiltonian action per parameter point. The derivative is taken of the same generator action and normalization that build the product state. BFGS minimizes the surrogate product energy \(\langle\chi_k|A|\chi_k\rangle\) over the angles of the product state defined in the opening paragraph. The method saved in result.plan.method records optimize_every_m, the optimization period in selected generators, optimize_rounds_n, the largest number of iterations of one BFGS run, and optimize_max_evaluations, the largest total number of energy evaluations. result.optimizer_rounds and result.energy_queries record the cumulative completed BFGS iterations and counted energy queries, while result.history records those counters at each projected analysis. The projected GCiM energy is then obtained by solving the pencil for the basis at the optimized angles. The energy at a parameter point already evaluated is reused without new seeds. SciPy's BFGS keeps no internal state that a resumed Run (one execution of a Plan, saved so that it can be reopened) could restore, so after a resume it restarts from the best completed point, and the result records the restart because the search path may then differ.

Continue and reanalyze

result.history keeps the settings, choices and observations that drove each decision, and post-hoc cutoff analysis cannot rewrite that history. run.resume() retrieves the jobs already submitted before it starts new work, and run.resume(reanalyze=True) explicitly recovers an interrupted projected analysis from complete saved data. A job whose state is unknown is not submitted again.

Each ADAPT round begins with its matrix stage, which obtains the matrix elements of the current basis. Before the matrix stage of quantum ADAPT starts its first query, the Run checks a lower bound on the stage's new queries against its preparation, circuit and shot limits, counting one preparation and one circuit for each query, and with shots shots shots for each query. In exact execution a query is one pair reduction (pencil.pair_count), and the classical work of reducing these pairs is also checked against max_products. With shots a query is one grouped setting (fixed_basis.sampled_setting_count) under the default grouped measurement, and one scalar Hadamard test (pencil.matrix_element_count) only under PencilAcquisition(matrix_measurement="termwise"). The bound leaves out the pairs of basis states that earlier stages obtained. When \(s\) of the \(b\) current basis states belong to earlier bases, it is the count for a complete pencil on \(b\) states minus the count on \(s\) states, both as given in Number of circuits and shots. Without product-state optimization each basis contains the one before it, so the stage that grows the basis from four to six states has a bound of \(20G\) grouped settings for \(G\) measurement groups, or \(11G\) under real-only measurement. A stage whose bound exceeds those limits is rejected before any of its work runs, and the error message gives the counted, requested and needed totals. Each query is also checked against the limits when it runs. After run.extend_limits(...), run.resume() runs the stage with the same seeds as a Run that was never rejected. run.extend_limits(...) raises the Run's limits, while ADAPT.max_products is a field of the Method that this call does not change. The stage in progress when a Run resumes may be partly run, so it skips the stage check and relies on the check of each query, which in exact execution includes the check of each pair reduction against max_products. result.save(...) and loading in the current format keep the pool, the compiler data of the selected generators (the angle-independent data from which their circuits are built) and enough observations to analyze the Result again, without planning again.

Verify an ADAPT result

Passing AdaptVerificationOptions to result.verify(checks=...) selects residual, sector or supplied-reference-energy checks. ProjectedVerificationOptions checks only the saved projected data. EnergyShiftOptions.for_result also supports FixedGCIM's original dense, CSR and CSC classical inputs. Its explicit matrix-entry evidence uses the same energy-shift check as Lanczos, keeps complex off-diagonals and implicit zeros, and performs no Pauli conversion or solve. A zero residual can belong to an excited eigenstate (Result 17). Sector expectations are not membership proofs, and signed spin contamination stays signed. A sampled result carries no joint statistical coverage or ground-state claim that was not supplied with it.

Circuit-level verification builds one Qiskit circuit for each selected generator, applied after the reference preparation, and records that construction work. Classical verification reuses the prefix states \(U_r\cdots U_1|\phi_{\rm ref}\rangle\) saved with the Result, also after the Result is loaded from a file. Both accumulate the Ritz state in basis order and free each basis or prefix state once it has contributed to the sum and the next state has been evolved from it. Their check against max_bytes covers these states, for circuit-level verification the largest gate matrix applied at once, and the arrays the saved Result already holds. Circuit-level residual and sector checks rebuild each basis state as a full statevector of \(2^n\) amplitudes on \(n\) qubits.

Before any state is computed, verification builds and decomposes, once, each circuit it will simulate. Each dense unitary in a supplied reference circuit is first replaced by its exact synthesis, whose work (_dense_synthesis.dense_synthesis_size) is counted against max_products and max_bytes before it starts and adds to the total. The limit check then counts each decomposed k-qubit gate as \(1+4^k+2^n(2^k+2)\) operations, for its instruction copy, its matrix, its application and two state copies, and adds \(2^n\) for the reference's zero state and for each nonzero global phase. The decomposition itself and Qiskit's workspace are outside this count. With an occupation reference, the count gives these sizes:

Check	Qubits	Operations
Circuit-level check on the spin-adapted pool, any selection of up to eight generators	8	at most 91.1 million, plus the shared Pauli action of the Hamiltonian for a residual
One shared-index double	12	up to 104 million
One four-distinct double	12	up to 56.7 million
One pair-split double	12	up to 10.3 million
A direct-vector reference alone	12	about 335 million

These counts depend on Qiskit's decomposition of the circuits. The Qiskit version and the circuits of each row, including the generator selection that gives the first row's maximum, were not recorded. The counts indicate the size of a check, and another Qiskit version can decompose the circuits into other gates and give other counts.

The default AdaptVerificationOptions.max_products=1_000_000_000 limits the counted state and operator work of verification, including the exact synthesis of dense unitaries in a supplied reference circuit. A check above it needs an explicitly raised max_products.

Chemistry inputs and reference energies

build_gcim_chemistry_problem(...) returns the molecular Hamiltonian, the occupation ordering and an explicit reference panel, which holds the RHF energy and the energy or status of each requested MP2, CCSD and CASCI reference. Its eigenproblem() keeps Hartree units and the recorded constant offset, and adapt_method(...) uses its occupation preparation. MP2, CCSD and CASCI run only when requested. Molecular-orbital signs from the eigensolver are arbitrary and can differ between LAPACK builds, and each orbital's sign fixes the sign of the generators that contain it. The builder therefore makes the first AO coefficient of each orbital above ORBITAL_PHASE_RELATIVE_THRESHOLD of its largest magnitude positive. This removes the dependence on the eigensolver's sign choice and the generator sign flips that choice could cause. It fixes signs only. Degenerate orbitals can still be returned in another rotation of their subspace, and binary64 rounding can still decide between generators whose gradients tie in exact arithmetic (Adaptive generator coordinates).

To show the stored reference panel next to an obtained energy, build a Chemistry References ReportSection from the chemistry input's metadata and the Result's eigenvalue:

from nwqlib import solve
from nwqlib.algorithms.gcim import (
    build_gcim_chemistry_problem,
    chemistry_reference_diagnostic,
    chemistry_reference_report_section,
)

data = build_gcim_chemistry_problem("H 0 0 0; H 0 0 0.74", basis="sto-3g",
                                    reference_methods=("mp2", "ccsd"))
result = solve(data.eigenproblem(), method=data.adapt_method(),
               execution="classical", seed=7)
print(result.eigenvalue)
diagnostic = chemistry_reference_diagnostic(data.metadata,
                                            energy=result.eigenvalue)
section = chemistry_reference_report_section(diagnostic)
print("\n".join(section.lines))

-1.1372838344885012
label: application context, not a validation gate
E_HF: -1.1167593074
E_MP2: -1.12989738099
E_CCSD: -1.13728399861
full-space CCSD correlation fraction: 0.999992003683

The first line is the ADAPT energy estimate result.eigenvalue in Hartree. The diagnostic reads only the stored reference panel and the supplied energy E, and the section only formats that diagnostic. Neither starts an RHF, MP2, CCSD or CASCI calculation or a simulation. The section lists E_HF, E_MP2 and E_CCSD, plus E_CASCI for an active space. A missing value shows its recorded status and reason, such as not_requested, failed or non_converged, or unavailable when no status is stored. When the needed energies are present, separate lines give the full-space CCSD correlation fraction \((E-E_{HF})/(E_{CCSD}-E_{HF})\) and the active-space CASCI fraction \((E-E_{HF})/(E_{CASCI}-E_{HF})\). A zero denominator is reported as undefined: E_CCSD == E_HF, or undefined: E_CASCI == E_HF for the CASCI fraction. section.data holds the same diagnostic as a mapping, and correlation_fraction(E, E_HF, E_CCSD) evaluates one ratio directly, raising ValueError for a zero denominator. These references give application context. They neither prove a bound on the total error nor identify the ground state. result.report() does not include this section, so print or store it next to the report.

Measured evidence

Example and section	Settings and environment
Six-qubit ring tie (Adaptive generator coordinates)	Classical ADAPT with max_iterations=1 and seed 7, the Hamiltonian passed through ingest_pauli as stated, read from the first round's gradients in result.history. Python 3.12.14, NumPy 2.5.2, SciPy 1.18.1, Qiskit 2.5.2 and Aer 0.17.2 on macOS arm64.
Chemistry reference panel (Chemistry inputs and reference energies)	The example as shown, with PySCF 2.14.0 and the versions of the row above.

Limits

• The eigenvalue is a projected value and does not establish the full-space ground energy (Read the result).
• The default total-error assessment remains INCONCLUSIVE, because six error sources have no automatic bound.
• The sampled enclosure diagnostic is not a sampling or accuracy bound.
• ADAPT's defaults and flat-counter stopping are operational settings, not convergence guarantees, and fixed-angle results depend on the sign of each generator (Defaults and stopping).
• The bounded Ritz growth score is a lower bound on the decrease of an exact solve of the current kept subspace together with the extension's part outside it, not a ranking by the energy of the next solve, and its stop reasons are operational. With shots its eligibility rule is empirical and assumes independently measured settings, not a confidence interval or an independence proof, and the sampled score is not a validated energy decrease (Bounded Ritz growth).
• A zero residual can belong to an excited state, and sector expectations are not membership proofs (Verify an ADAPT result).
• Exact quantum execution needs a target that runs the pair computations, such as Aer's statevector target or NWQ-Sim's CPU/SV target.
• Chemistry reference energies give context and no error bound.

Limitations and open work lists the open items for this method.

Numerical scope and planning limits

FixedGCIM(max_admission_steps=10_000_000), the default, limits the planning work for each circuit description the method builds and does not change the quantum operations. This default allows a sampled FixedGCIM plan with three basis states and a 20-qubit, 8,083-term molecular Hamiltonian in 1,285 measurement groups, or 11,565 circuits with both the real and imaginary parts measured. When the planning work exceeds this limit, the error message names the field, the planning stage that exceeded it and its count. Raise the field to the reported count, and further when the count is marked as a lower bound, because later planning stages can need more (planning work limit).

Classical GCiM fills normalized trial columns \(V\) and shifted actions \(W=(A-cI)V\) in two preallocated complex arrays, which together hold at most \(32bd\) bytes for \(b\) trial columns of length \(d\), the original operator dimension, with \(A\), \(c\) and \(I\) as in How matrix elements are obtained. A state object that appears more than once in the basis shares one stored column and action. Conjugating inner products form \(H_0=V^\dagger W\) and \(S=V^\dagger V\) without an additional state block. The projection's byte count adds the \(b\)-by-\(b\) matrices and the preparation, operator and action workspace. The stored physical Hamiltonian is \(H=H_0+cS\). Pauli input acts on all columns with one grouped Pauli action. Dense input is shifted in row blocks, each copied into the same buffer. CSR and CSC input is shifted once per Run into one read-only sparse matrix, which holds a full copy of the stored values and shares the index pattern of the input only when every diagonal entry is stored. Its bytes are given in Engineering constants. The ordered columns and their complete pencil \((V^\dagger(A-cI)V,\,V^\dagger V)\) are kept, also when the columns are linearly dependent.

ADAPT planning checks the size of the growing basis before any query. After \(k\ge1\) selections the basis holds \(2k\) listed states, while at \(k=0\) it contains only the reference. Thus 2*min(pool size, max_iterations) bounds the projected matrix dimension and must not exceed max_basis_size (default 64). For an explicit pool this bound is checked before conversion and commutator formation. A named pool is checked after it is enumerated, under that enumeration's own byte and product limits. A bounded Ritz growth basis holds at most min(max_basis_size, max_iterations + 1) states instead. Planning raises ApplicabilityError for a Problem that names a fixed subspace, because the basis of ADAPT changes from round to round.

ADAPT queries

Classical ADAPT screening reuses the action of \(A-cI\) on the current normalized product state when the matrix-element computation has already produced it. In exact quantum execution, one query on the full product state \(|\chi_k\rangle\) supplies both its diagonal Hamiltonian entry and the commutator expectations of the unselected generators. The matrix stage runs this query before the projected analysis, and the screening that follows reuses its values. The diagonal overlap is one under the normalized-preparation convention. Both kinds of value come from one simulation, whose labels, preparation record and data content hash are kept through continuation and loading. A round that stops after its matrix stage has therefore already paid for the screening values of that query.

Quantum queries use the compiler data of the generators, the pool indices and angles of the product states, the ancilla control and phase, and the measured register. Before any query runs, ADAPT checks against the limits an upper bound on the number of circuit instructions that building the generators can append. That bound uses the system width, the Pauli term count and the supported circuit family of the generators. A generator's compiler data is constructed when that generator is selected and reused across its later angles. The count of appended instructions and the later cost of synthesizing them into backend gates are reported separately. Classical queries use the same processed input and iteration state, with explicit state and action accounting. Only a quantum Plan forms the commutators \([A,A_j]\) whose Pauli labels its screening queries measure. Only a quantum Plan with shots also groups those labels and the Hamiltonian's into qubit-wise commuting measurement groups, whose comparisons are counted against max_products.

nwqlib.estimate(plan) treats the experiments of an ADAPT Plan in plan.experiments, such as screen, pair and energy, as query templates, because the iteration sets the generators, angles, Pauli term or measurement group and register widths of each query, and the number of queries depends on the outcomes. A count is reported as zero only when the steps and measurement batches of the templates show that no query can contribute to it. These zeros are the exact evaluations with shots, the shots in exact execution, and in classical execution every count of circuit operations and measurements, including settings, shots and exact evaluations. Every other count, such as the settings, gate counts, measurements and resets of quantum execution, the classical work and the construction work, is unavailable with the reason that plan.method.query_catalog_reason() returns. In quantum execution the qubit and bit widths are formulas in the ancilla and readout widths that the iteration sets for each query, and each width and memory peak describes the templates at their declared parameter values rather than the largest query of a Run (Adaptive Methods). The resources of one query are estimated from its own construction, which its preparation record gives, as here for the first query that prepare builds:

from nwqlib import Eigenproblem, plan, prepare
from nwqlib.algorithms.gcim import ADAPT
from nwqlib.operators import ingest_pauli
from nwqlib.resources import estimate as estimate_construction

problem = Eigenproblem(A=ingest_pauli((("Z", 1.0),), num_qubits=1))
method = ADAPT(initial_state=[1, 1], pool=(ingest_pauli((("Y", 1j),), num_qubits=1),))
adapt_plan = plan(problem, method=method, shots=100, seed=7)
with prepare(adapt_plan).run as run:
    record = run.data.receipts[0]
query_estimate = estimate_construction(record.realization.selected_construction(adapt_plan))
print(query_estimate.quantity("settings"))
print(query_estimate.quantity("shots"))

settings: 1 count [exact]
  basis=selected_logical; lifecycle=planned; population=dynamic terminal MeasurementBatch setting visits
shots: 100 count [exact]
  basis=selected_logical; lifecycle=planned; population=terminal MeasurementBatch sampled shots

run.data.receipts holds the preparation record of each query that the Run has prepared, so the same call estimates any later query. For a query of an exact Plan the estimate gives one exact evaluation instead of the shots, and gate counts without a counting formula stay unavailable (Estimate resources).

Reanalysis and saving check the saved parts of each query's observations against the obtained data and do not run the quantum circuits again (saved run mechanics).

Source and code map

The projected Hill–Wheeler equation is from Zheng et al., Physical Review Research 5, 023200 (2023), Eq. (13) in the numbering of arXiv:2212.09205v1, DOI. The adaptive method is from Zheng et al., npj Quantum Information 10, 127 (2024), arXiv:2312.07691v3. Equation, section and page numbers refer to the arXiv version named in each row. Rows marked NWQLib have no paper source, and the docstring of the code in the last column states the derivation or rule. A bare file name in the Code column is in src/nwqlib/algorithms/gcim/, except _projected_eigensolver.py, which is in src/nwqlib/, and fermionic_pool.py and fermionic_circuits.py, which are in src/nwqlib/subroutines/.

Step	Source	Location	Code
Projected Hill–Wheeler problem H f = E S f	Zheng et al. (2023), arXiv:2212.09205v1	Eq. (13), p. 2	_projected_eigensolver.py
Matrix elements and their Pauli-sum form	Zheng et al. (2023), arXiv:2212.09205v1	Eqs. (14)–(15), p. 3 Eqs. (33)–(34), p. 5	pencil.py: matrix_elements, assemble_pencil
Hadamard-test amplitudes and per-shot variance	Zheng et al. (2024), arXiv:2312.07691v3	App. G, Eqs. (G4)–(G5) and (G8), pp. 14–15	adapt_acquisition.py: construct_adapt
Grouped measurement in sampled FixedGCIM and ADAPT: joint-state ancilla X/Y quadratures, qubit-wise-commuting system groups, weighted group moments with their pooled variances and covariance, and the byte and work formula of one grouped analysis pass	NWQLib	Derivation in the docstrings of sampled_settings, reduce_group, pooled_statistics and grouped_pass_requirements, which both methods call, and of _joint_state_blocks, _sampled_construction and read_groups, which build the FixedGCIM circuits and reduce their counts	fixed_basis.py: sampled_settings, reduce_group, pooled_statistics, grouped_pass_requirements, FixedGCIM.plan, _analyze_groups
Grouped measurement in sampled ADAPT: the joint-state circuit of each grouped setting and the orientation and pooling of its counts	NWQLib	Derivation in the docstrings of construct_adapt, _grouped_moments and _grouped_statistics, which pool through the FixedGCIM functions of the row above	adapt_acquisition.py: construct_adapt, _grouped_moments, _grouped_statistics, _matrix_entries
Structural proof of real \(H\) and \(S\) for real-only measurement	NWQLib	Conditions on the Pauli coefficients, the stored state amplitudes and the generator terms, in the docstrings of the functions in the last column	fixed_basis.py: real_hamiltonian_reason, real_state_reason, real_generator_reason, resolve_acquisition
Pilot and main shot allocation (PencilSampling)	NWQLib	Derivation of \(v_s\), the floored proportional rule and its integer rounding in the module docstring of pencil_sampling.py, and the two stages in the docstrings of pilot_decision and advance_stages	pencil_sampling.py: influence_variances, pilot_split, floored_allocation, pilot_fallback_reason pencil_workflow.py: pilot_decision, advance_stages adapt.py: _admit_final_stages
Joint-state pair identity and weighted-pencil reduction	NWQLib	One joint preparation per pair, with the ancilla on qubit 0, its zero branch holding phi_i and its one branch phi_j. Entry bounds and pooling are in the pair_reducer docstrings	pair_reducer.py: entry_bounds, pooled_bound, pair_work, pair_bytes fixed_basis.py: _exact_pair_construction, _exact_overlap_allowance
ADAPT joint states that control only the reference (pair_preparation="vacuum_reference" and "vacuum_restore")	Yoshioka et al., arXiv:2407.14431v3, for the controlled reference preparation NWQLib for the vacuum rule, the vacuum phase and its correction, the restoration with the adjoint and the instruction bounds	Fig. 1b and main text, pp. 2–3 Derivations in the docstrings of _admit_generator_routes, _vacuum_phase, _realized_adjoint and _append_vacuum_pair and the comment above the bound in _construction	adapt.py: _admit_generator_routes, _construction adapt_acquisition.py: pair_route, _vacuum_phase, _realized_adjoint, _append_vacuum_pair
Positive-overlap truncation (canonical orthogonalization)	Epperly, Lin and Nakatsukasa, arXiv:2110.07492v2 Zheng et al. (2024), arXiv:2312.07691v3	Algorithm 1.1, Sec. 1.2, p. 5, and Theorem 2.7, p. 14 App. F, Eqs. (F1)–(F3), p. 14	_projected_eigensolver.py
Identity offset removed before the projected solve	NWQLib (Proposition 14)	Exact shift covariance of the pencil	fixed_basis.py: solve_pencil
Sampled Ritz enclosure diagnostic	NWQLib	Pauli L1 bound and Rayleigh–Ritz interlacing	fixed_basis.py: processed_pauli_enclosure
Gradient screening with \(\langle[A,A_i]\rangle\)	Grimsley et al. (2019), arXiv:1812.11173v2 Anastasiou et al., arXiv:2306.03227v3 NWQLib screens only unselected generators, while Grimsley et al. step 7 keeps selected operators in the pool.	Sec. II.B, steps 5–7, pp. 3–4 Eq. (2), p. 3	adapt_inputs.py: prepare_adapt_inputs adapt_acquisition.py: drive_adapt, _screen_gradients, _offset_free_column
Surrogate product state, fixed angle and two new basis states per selection	Zheng et al. (2024), arXiv:2312.07691v3	Main text, p. 4 METHODS, p. 9	adapt_acquisition.py: basis_chains adapt.py: ADAPT
Bounded Ritz growth: examined extensions, CX count bound, score (sqrt(delta**2 + 4 rho**2) - delta)/2 and sampled eligibility	NWQLib	Round-robin examination, whole-chain comparison, the CX table, the two-state Rayleigh–Ritz bound and the first- and second-order change of eta in the docstrings of frontier, generator_cx_weight, reference_cx_bound, candidate_score, ritz_gain, overlap_proxy, eta_gradients, sampled_metric and sampled_eligibility	bounded_growth.py: frontier, generator_cx_weight, reference_cx_bound, current_projection, candidate_score, ritz_gain, overlap_proxy, eta_gradients, sampled_metric, sampled_eligibility adapt.py: BoundedRitzGrowth adapt_acquisition.py: drive_adapt
Flat-counter stop	Zheng et al. (2024), arXiv:2312.07691v3	METHODS, p. 9, T = min(T_auto, T_usr), with T_auto equal to 20% of the unselected operators	adapt.py: ADAPT adapt_acquisition.py: drive_adapt
Intermittent, truncated BFGS optimization	Zheng et al. (2024), arXiv:2312.07691v3	App. H, p. 17 Table III caption, p. 8	optimization.py: optimize_product
Classical product energy and its Taylor-model gradient	NWQLib	Normalized adjoint of the fixed-step Taylor factors, derivation in the normalized_taylor_energy_gradient docstring	optimization.py: normalized_taylor_energy_gradient, taylor_action_derivative adapt_acquisition.py: _classical_energy_gradient
Generalized parameter-shift rule	Wierichs et al., arXiv:2107.12390v3	Sec. 2.1, Eqs. (1)–(3) Sec. 3.2, Eqs. (16)–(18) Sec. 3.4, Eq. (24) App. B.1	optimization.py: _generator_shift_rule
Pauli insertion derivative	Schuld et al., arXiv:1811.11184v1, for the two-eigenvalue identity NWQLib for inserting after the complete exponential (Proposition 15)	Sec. III.A, Theorem 1, Eq. (8), and Eqs. (13)–(14), pp. 3–4 Sec. III.B, Eqs. (16)–(17), p. 5	optimization.py: _generator_shift_rule adapt_acquisition.py: construct_adapt
Spin-orbital order and default spin-adapted GSD pool	Zheng et al. (2024), arXiv:2312.07691v3 Magoulas and Evangelista, arXiv:2511.13485v2	App. E, Eqs. (E1)–(E3), p. 13 Sec. III, Eqs. (3), (6), (7)	fermionic_pool.py: enumerate_spin_adapted_gsd_pool
UCCSD-SD, QEB-SD and OVP-CEO pools	Romero et al., arXiv:1701.02691v2 Zheng et al. (2024), arXiv:2312.07691v3 Yordanov et al., arXiv:2011.10540v2 Ramôa et al., arXiv:2407.08696v3	Eqs. (8)–(9), p. 3, and Eq. (17), p. 4 App. E, Eqs. (E4)–(E5), p. 13 Eqs. (17)–(18), p. 6 Sec. II.B.3, Eqs. (23)–(24), p. 8	fermionic_pool.py: enumerate_uccsd_sd_pool, enumerate_qeb_sd_pool, enumerate_ceo_ovp_pool
Orientation, sqrt(2) normalization, factor order and signs of the compact circuits	Magoulas and Evangelista, arXiv:2511.13485v2 Zheng et al. (2024), arXiv:2312.07691v3	Sec. III, Eqs. (2), (5)–(7), p. 4, and Sec. V, text after Eq. (15), p. 6 App. E, Eqs. (E2), (E3) and (E8), p. 13	fermionic_circuits.py: module docstring, _append_excitation
Pair/split double circuit	Magoulas and Evangelista, arXiv:2511.13485v2	Sec. V, Eqs. (14)–(15), (21), (22), (25)–(27), pp. 6–7	fermionic_circuits.py: build_generator_circuit
Four-distinct double circuits	Magoulas and Evangelista, arXiv:2511.13485v2	Sec. VI, Table I, p. 13, and Table II, p. 14 Supplemental Tables SI–SII, pp. S16–S18	fermionic_circuits.py: build_generator_circuit
Shared-index double circuits	NWQLib (Proposition 16)	Exact block synthesis with fermionic swaps and Givens rotations	fermionic_circuits.py: _append_shared_index
Classical generator exponential exp(theta A)	NWQLib	Scaled degree-18 Taylor steps with ||step A|| <= 1/2, derivation beside GENERATOR_TAYLOR_DEGREE	fermionic_pool.py: apply_generator_exponential, _taylor_steps adapt_actions.py: packed_exponential, the same steps on the classical Plan's packed tables
Residual nearness interval and Kato–Temple heuristic	NWQLib, from standard Hermitian residual bounds (Result 17)	Derivation in the residual-metrics docstring	adapt_verification.py: _residual_metrics, verify
Sector expectations with S^2 = S_- S_+ + S_z (S_z + 1)	NWQLib	Interleaved spin-orbital order of the default pool	sector.py: sector_expectations fermionic_pool.py: apply_spin_squared
Byte and work formulas	NWQLib	Units in Byte and work budgets. Each function maps every term of its formula to the arrays or passes it counts and names any untuned limit term	pencil.py: _projected_solve_bytes, _projected_solve_work fixed_basis.py: FixedGCIM.plan, _plan_classical, grouped_pass_requirements, group_requirements, read_groups adapt.py: _construction, selected_kernels, generator_pool_slot_bound adapt_records.py: _exponential_action_work, _energy_gradient_work adapt_actions.py: action_requirements, build_action_tables, action_sizes, query_work adapt_inputs.py: _pool_enumeration_requirements, array_record_bytes adapt_acquisition.py: bind_classical adapt_verification.py: verify, _state_frontier sector.py: _sector_work fermionic_pool.py: _reference_pool_size
Molecular Hamiltonian and RHF orbital signs	NWQLib, with PySCF RHF and integrals	ORBITAL_PHASE_RELATIVE_THRESHOLD convention	chemistry.py: build_gcim_chemistry_problem