QPE¶

The QPE methods estimate an eigenvalue of a Hamiltonian \(H\), or an eigenphase of a unitary \(U\), from the signal \(z_p=\langle\psi|U^p|\psi\rangle\) of a prepared state \(\psi\), with \(U=e^{-i\tau H}\) for a Hamiltonian. Real and imaginary Hadamard tests give \(z_p\) (Ding and Lin, arXiv:2211.11973v2, Fig. 1 and Eqs. (4)–(5)). Each method estimates a different quantity.

QCELS fits one complex exponential \(z_p\approx ae^{-ip\theta}\) at integer powers \(p\) and returns the energy \(\theta/\tau\) (Ding and Lin, arXiv:2211.11973v2, Eqs. (2) and (9)–(12), Result 3).
SPE, statistical phase estimation, returns the lowest energy in the prepared spectral support, at the first point where a filtered estimate of the spectral cumulative distribution crosses half a declared lower bound \(\eta\) on the target component's probability (Wan, Berta and Campbell, arXiv:2110.12071v2, Eqs. (2) and (6) and Algorithm 1, Result 4).
RFE, randomized Fourier estimation, draws powers uniformly, forms sampled Fourier coefficients of the signal at \(K\) frequencies \(j=0,\dots,K-1\) and returns the grid angle \(2\pi j/K\) of the frequency whose sampled coefficient has the largest magnitude (Kshirsagar, Katabarwa and Johnson, arXiv:2209.11322v3, Eqs. (5)–(7) and Algorithm 1, Result 4).
RWPE, random-walk phase estimation, updates a Gaussian estimate of the phase from one measured bit per step and reports its mean with a Gaussian interval (Granade and Wiebe, arXiv:2208.04526v1, Eq. (7) and Algorithm 1, Result 5).

The preparation determines which part of the spectrum is observed. A single-mode estimate establishes neither that the value is the ground energy nor that it covers the dominant component. The four methods share the controlled-power construction and the execution steps. For the lowest eigenvalue in a trial space built from one initial state, use Chebyshev Lanczos, and for supplied or adaptively grown trial states, use GCiM.

from nwqlib import Eigenproblem, solve
from nwqlib.algorithms.qpe import QCELS, SPE, RFE, RWPE

problem = Eigenproblem(A=[[0.2, 0.0], [0.0, 0.7]])
result = solve(problem, method=QCELS(initial_state=[1, 0]), seed=7)
print(result.eigenvalue, result.phase)

0.20000000000000012 0.8714285714285713

The prepared state \((1,0)\) is the eigenvector of eigenvalue 0.2, and result.phase is the same eigenvalue as an eigenphase of \(U=e^{-i\tau H}\) in turns (Read the result). This default QCELS example runs a two-qubit circuit on the Aer simulator. Other choices include SPE(initial_state=[1, 0], overlap_lower_bound=1.), RFE(initial_state=[1, 0]) and RWPE(initial_state=[1, 0]). Replace A and initial_state to treat your own problem. The eigenvalue notebook (examples/gcim_lanczos_qpe_eigenvalue_intro.ipynb) applies QCELS to a stretched H4 chain together with ADAPT-GCiM and Lanczos.

Choose an estimator¶

Method	Returns	What its accuracy statement assumes	Interval	Default measurement work
`QCELS`	Energy of a one-frequency least-squares fit to \(z_p\)	Ding–Lin's accuracy results, Theorems 1–2 (arXiv:2211.11973v2, p. 14), assume \(p_0>0.71\), the squared overlap of the preparation with the target eigenvector, and their data set uses consecutive powers (Eq. (6)). NWQLib does not check \(p_0\) (When the accuracy theorems apply).	None	22 elementary settings, two quadratures at each of 11 powers
`SPE`	Lowest energy in the prepared support	A declared lower bound \(\eta\) on the target component's probability, which planning does not verify. The search guarantee of Wan, Berta and Campbell (arXiv:2110.12071v2, Theorem 1 with Eq. (11)) does not transfer to NWQLib's grid scan (SPE).	None	32 paired draws, 64 quadrature draws
`RFE`	Grid angle \(2\pi j/K\) of the frequency, among \(K\) frequencies \(j=0,\dots,K-1\), with the largest-magnitude sampled Fourier coefficient	Theorem 2.1 of arXiv:2209.11322v3 (p. 6), for an eigenstate, needs more than 1100 single-shot draws for any phase tolerance \(\epsilon<\pi/2\) and failure probability \(\delta_{\mathrm{RFE}}<1\). The default number of random power draws, \(M=97\), is a finite workload without that guarantee (RFE).	None	97 paired draws, 194 quadrature draws
`RWPE`	Mean of a Gaussian phase estimate	The basic Gaussian walk of Granade and Wiebe (arXiv:2208.04526v1, Eq. (7)), whose mean can move at most about 2.96 prior widths (RWPE).	Nominal 95-percent interval of the moment-matched Gaussian model	14 one-shot steps

None of these intervals is frequentist coverage or a total physical-error bound (Intervals and sampling error).

Inputs¶

The input is either an Eigenproblem together with initial_state on the method, or a SpectralEstimation, which holds its Hamiltonian or unitary and its initial state. With an Eigenproblem, each method requires initial_state. A SpectralEstimation already holds its initial state, so the method then rejects an initial_state of its own. Changing the method's initial_state does not change the target of the Eigenproblem. An initial state that is an excited eigenstate can return that excited energy while the problem still asks for the smallest eigenvalue, and the Result states this limitation. RWPE needs a Hamiltonian, because its evolution times are not integer powers of a unitary (RWPE).

A small dense Hamiltonian \(A\) whose dimension \(d\) is not a power of two is padded to the next power of two with a Hermitian dummy block with diagonal entries \(\operatorname{tr}(A)/d\), and the initial state is padded with zeros. The Result still refers to the original problem and initial state. A unitary is padded with an identity block. Scalar targets use one control qubit and no system register or state-preparation block.

Dense input uses controlled_power_backend="dense_exact" by default. A Hamiltonian power is the exact spectral power, formed from one eigendecomposition shared by all powers, so the algorithmic error recorded for its controlled block is 0. For a dense unitary input \(A\), the powers are built from its unitary polar factor \(V\). Every requested power, and every gap between consecutive powers of the controlled trajectory (How the signal is measured), is a power of this same \(V\), so the target operator at power \(p\) is \(V^p\). A supplied unitary is accepted when every entry of \(U^\dagger U-I\) is at most 1e-8 in magnitude. The polar decomposition, the floating-point matrix powers and the synthesis of the controlled gates contribute numerical error. Each controlled matrix is synthesized as its unitary polar factor to binary64 rounding (Controlled dense unitaries). Under the default max_work, each dense block fits on up to seven system qubits (Numerical scope and planning limits).

Read the result¶

result.phase is in turns, in \([0,1)\), with \(Uv=e^{2\pi i\,\mathrm{phase}}v\). For a Hamiltonian, \(U=e^{-i\tau H}\), so result.eigenvalue is in the original matrix units and includes the minus sign of the conversion. Unitary-only input has no Hamiltonian eigenvalue without a chosen logarithm branch, and accessing that field raises AttributeError. result.value is phase for an Eigenphase output and eigenvalue otherwise.

SPE returns the lowest value in the prepared spectral support. For a Hamiltonian that is the lowest energy. For a unitary it is the lowest principal eigenphase angle \(\theta\), with \(Uv=e^{i\theta}v\) and \(\theta\) scanned on \([-\pi/2,\pi/2)\). When \(\tau\lvert E\rvert<\pi/2\) for every prepared energy \(E\) of a unitary \(U=e^{-i\tau H}\), \(\theta=-\tau E\), so SPE on \(U\) returns the phase of the highest energy of the prepared support. Supply \(H\) itself, or \(U^\dagger=e^{i\tau H}\), to target the lowest energy.

estimator_value is the value the estimator computes, an energy for a Hamiltonian and a principal eigenphase in radians for a unitary. The phase and the interval are converted from it once. For an Eigenphase output from a Hamiltonian, the energy endpoints of the interval are multiplied by \(-\tau/(2\pi)\), which reverses their order, and the interval is recentered around the reported phase. For an Eigenvalue output, the interval stays in energy units. The interval is not wrapped, because wrapping each endpoint separately would change its meaning. It is the estimator's own interval and neither bounds the total error nor identifies the ground state.

Exact readout, the default of QCELS, SPE and RFE when shots is not given, computes each ancilla expectation without shot sampling, on the simulator or by classical execution, subject to numerical rounding. For quantum execution with exact readout, each saved ancilla expectation of the executed trajectory must have magnitude at most 1 up to a binary64 roundoff window, or the observation is rejected. The window bounds roundoff in the saved expectation and is not an accuracy bound for the estimate (Readout tolerance branches). It does not cover the operations or conditions that the circuit's preparation record, the saved record of the circuit as it was built with its operation counts, lists as roundoff exclusions, because the roundoff derivation does not treat them.

Field	Meaning
`eigenvalue`, `phase`, `value`	The energy in the original units, the eigenphase of \(U\) in turns, and the requested output, as described above.
`estimator_value`	Energy in operator units for a Hamiltonian, or principal eigenphase in radians for a unitary. None when no estimate was identified.
`interval`	RWPE's nominal model interval, None for the other estimators (Intervals and sampling error).
`fit`	QCELS search record: effective grid, objective evaluations, residual and interval method.
`gaussian`	RWPE's Gaussian phase mean and standard deviation.
`exposure`	One `QPEExposure` per sampled SPE or RFE query, that is, per measured (power, quadrature) setting, with its draw `multiplicity` and its `requested` and `received` shots.
`requested_exposure_complete`	Whether every sampled query returned its requested shots.
`minimum_effective_shots_per_draw`	The smallest `received/multiplicity` over the queries, as an exact (numerator, denominator) pair.
`rfe_coordinate_variance_upper`	An upper bound, in units of squared ancilla expectation, on the received-count variance proxy of each real and imaginary coordinate of the sampled Fourier coefficients. It concerns the conditional shot contribution under independent stationary shots, not the random-draw term, and gives no energy confidence interval.
`stop_reason`, `missing`, `complete`	Why measurement stopped or no estimate is identified (for QCELS, also that the fit is aliased), the planned queries without data, and whether all planned data and an estimate are present (for RWPE, whether every planned step was processed).

Each sampled SPE or RFE query contributes its mean over the shots returned, weighted by its draw multiplicity, the number of times the random schedule drew its power, and every sampled query needs at least one valid shot. Under outcome-independent sampling with stationary independent shots, the expected weighted contribution is the multiplicity times the query signal, and its shot variance scales as the square of the multiplicity divided by the received shots. Sampling bounds use these shot numbers and keep the uncertainty from drawing Fourier powers separate. When every query receives its requested shots, the pooled data have the distribution of separate draws measured with equal shots, while a query that returns fewer shots has a different shot-noise distribution. The finite cumulative-distribution crossing and the Fourier peak do not provide an energy confidence interval. The recorded measurements and seeds describe the pooled shots of each query, and a record for each individual draw would require separate measurements. NWQLib's two-qubit tests count the circuit submissions to Aer and their shots, verify multiplicities and check dense synthesis reuse without computing reference spectral weights. They do not establish mixed-spectrum accuracy or interval coverage.

Intervals and sampling error¶

QCELS returns interval=None, because the finite complex fit has no matched uncertainty model for the measurement. result.fit keeps its residual, effective grid and objective-evaluation count, and a flat objective has no identified estimate. SPE and RFE also return interval=None, because a cumulative-distribution crossing or a Fourier peak does not supply a confidence interval. RWPE reports the nominal 95-percent interval of its moment-matched Gaussian model, with its phase mean and standard deviation in result.gaussian. The Gaussian approximation, the prepared spectral weights and controlled-evolution error limit its interpretation. It is not a frequentist coverage or total physical-error guarantee.

The error model of the Plan, the construction and its costs computed before any circuit exists, lists every QPE error source as unknown (How a PASS is reached). With exact readout, result.facts states a zero sampling error for QCELS, whose schedule is fixed, and for an RWPE Result that reports its prior without an update. SPE and RFE average over randomly drawn frequencies or powers, and each RWPE update uses one random bit, which classical execution emulates from the exact signal. Their estimates can therefore change with the seed even with exact readout. The sample-count conditions of Wan, Berta and Campbell, arXiv:2110.12071v2, Eq. (11), and of Kshirsagar, Katabarwa and Johnson, arXiv:2209.11322v3, Theorem 2.1, control this error. NWQLib attaches no bound derived from them, so result.assess(component="sampling", ...) stays INCONCLUSIVE with sampling among the remaining sources.

Settings and measurement work¶

Method	Estimator controls and defaults	Meaning
`QCELS`	`num_times=16`, `max_time=None`, `grid_size=4096`	Integer powers plus zero and paired quadratures. `None` selects maximum power 10. A supplied `max_time` selects `floor(max_time/tau)` (unitary: `floor(max_time)`). Complex least squares uses a finite grid with bounded bracket refinement.
`SPE`	Required `overlap_lower_bound`, `num_samples=32`, `fourier_degree=11`, `filter_beta=6`, `grid_size=4096`	Sample positive odd Fourier frequencies proportional to their coefficient magnitudes, measure paired quadratures, and locate the first crossing of \(\eta/2\) by the approximate cumulative distribution.
`RFE`	`num_samples=97`, `num_frequencies=49`	Draw M independent uniform powers from 0 through K-1, measure both quadratures for each draw, and select the largest sampled Fourier coefficient.
`RWPE`	`max_steps=14`, `prior_mean=0`, `prior_std=pi`	Basic Gaussian random walk with one datum per continuous-time experiment. Zero steps reports the prior without measurement.

In the RFE row, M is num_samples and K is num_frequencies. Each method also accepts the shared controls tau (default None, chosen as Time step and energy range describes), controlled_power_backend (default "auto", which picks "trotter_error_budgeted" for Pauli input and "dense_exact" otherwise), controlled_power_error_budget (default 1e-3, the operator-norm error allowed for each power) and pauli_pruning_rtol (default 1e-12). SPE's maximum Fourier frequency is 2*fourier_degree+1. The common max_power separately caps the absolute evolution time in units of tau, including RWPE's noninteger times.

QCELS, SPE and RFE, the static methods, fix their measurement schedule when the Plan is made and by default use exact readout of one controlled trajectory (How the signal is measured). A positive shots sets the shots per elementary QCELS setting and per random SPE or RFE draw, so an SPE or RFE setting drawn m times receives m times shots, and backend limits on shots per circuit apply to that combined number. RWPE uses one shot per feedback step and defaults to 14 steps. It rejects a different requested shot count, because its update is defined for one Bernoulli datum.

The default shot workload is 22 elementary settings for QCELS and 14 one-shot measurements for RWPE. SPE and RFE draw 32 and 97 paired frequencies respectively, so their 64 and 194 quadrature draws request 64 and 194 times shots coherent shots. SPE and RFE make one query per distinct drawn (power, quadrature), so U distinct drawn powers give 2U queries, and at RFE's defaults the expected number of distinct settings is about 84.7. With quantum execution and exact readout, every query is observed once, on one controlled trajectory. A repeated sampled query receives the combined number of fresh shots.

With quantum execution, RWPE plans max_steps circuit preparations, measurements and one-shot observations, and it checks all three totals before preparing the first feedback query. The caller can cancel the Run, the execution of the Plan that records its submissions and observations and can be saved and reopened, before all steps are done. The feedback angle depends on the preceding outcomes, so only the next circuit is prepared at a time.

When a trajectory repeats controlled applications of one base unitary, the number of applications is the largest requested power. A trajectory of dense gap powers has its own matrix-construction and gate costs. Resource reports name the construction and count each executed increment once. A quantum QCELS, SPE or RFE Plan with exact readout is therefore submitted to Aer once, and the operation count in its preparation record counts the executed circuit operations through the last measurement point once.

Classical execution¶

execution="classical" evaluates the estimator signal from an exact eigendecomposition and the prepared spectral weights. The signals of the static methods take one eigendecomposition, one projection of the initial state and \(O(D)\) scalar work per distinct power. Classical RWPE then emulates one Bernoulli datum from each exact signal with the seed saved for that step. The exact signal remains in the observation, and the emulated bit does not count as a quantum shot. Within a Run, RWPE and the dense powers built for circuits reuse an eigendecomposition already computed for the same operator. The classical routine (a host kernel in the code) records how many eigendecompositions and projections it performed, separately from the conservative work estimate declared before it runs. Classical execution expands neither a supplied preparation circuit nor Pauli input into dense form, and it rejects both.

Inspect the circuits and their resources¶

prepare(plan).circuits exposes the circuit of the setting that the method submits first. prepare(plan, settings="all") exposes the circuits of every setting of a QCELS, SPE or RFE Plan on the local backends that Run on a backend names. RWPE rejects it, because each feedback angle depends on the outcomes before it. Printing or estimating does not build another circuit.

Each Suzuki step records its operation and CX counts for the Pauli terms kept after pruning and for its step_time. On the common grid of an exact trajectory, each power position refers to that shared step and records its cumulative step count. In the selected_logical basis a step counts 2L operations for L terms, one controlled Pauli rotation per term in the forward sweep of the symmetric formula and one in the reverse sweep. In the cx basis it costs \(4\sum_jw_j\) CX, where \(w_j\) is the number of nonidentity symbols of term \(j\) (powers.suzuki_step_cx, Proposition 7). A controlled rotation of support w becomes a parity chain and its inverse, 2(w-1) CX, around a controlled RZ of 2 CX. The CX count is exact for Qiskit 2.6.0 decomposition into cx and one-qubit gates at optimization level 0 without a coupling map, before cancellation, rotation merging or routing. The identity phase P, the feedback RZ and the Hadamard gates are one operation each and cost no CX. With a preparation that has no CX, such as an occupation state, a Hamiltonian made only of the identity term therefore has an exact zero CX total. The counts of a repeated step are multiplied by the number of times it is executed. A common-grid exact trajectory counts its increments through the final position once, while sampled experiments count their own preparations and shots. Thus estimate(plan) returns the whole-workload operation count and estimate(plan, context=ResourceContext(basis="cx")) its CX count. The whole-workload counts stay unavailable when the initial-state preparation has no counting formula. A dense power supplies neither formula, so its operation and CX counts are unavailable.

Time step and energy range¶

For static Hamiltonian queries, integer powers identify energy only modulo \(2\pi/\tau\). The automatic choice of \(\tau\) uses a bound on the spectrum that needs no eigensolve. Without a supplied tau, planning chooses \(\tau=0.9\,r/R\), where \(r\) is the principal radius of the estimator, stated below, and \(R\) is the dense row-sum bound or the Pauli L1 norm. A unitary has no \(\tau\). QCELS and RFE use a principal radius of \(\pi\). SPE also needs a filter transition margin \(\delta\), since its scan evaluates \(F(x-\tau E)\) near the periodic wrap boundary. A sufficient condition is \(\tau R+\delta\le\pi/2\), where \(R\) bounds the spectrum, and SPE's automatic radius is \(\min(\pi/3,\pi/2-\delta)\), followed by the common 0.9 interior margin. Theorem 3 of Wan, Berta and Campbell, arXiv:2110.12071v2, supplies \(\delta\) from beta and the declared overlap, and its degree assumption is checked separately. An unavailable assumption or a supplied time outside that domain appears in the Result. A crossing at the scan boundary is unidentified and is not reported as an endpoint energy. For raw-unitary input, SPE has no established interval that encloses the prepared spectrum. RWPE uses continuous times and a Gaussian prior, whose finite exploration range (RWPE) is checked against the same spectral bound \(R\).

How the signal is measured¶

With quantum execution and exact readout, the static methods prepare the control and system once and then apply one controlled trajectory, a sequence of controlled increments of \(U\) whose cumulative powers reach each planned power in turn. The ancilla X and Y expectations at each power position give the real and imaginary parts of \(z_p\), with \(\langle Y_a\rangle=\operatorname{Im}z_p\). Sampled X and Y settings use separate preparations, and RWPE follows its feedback-dependent measurement schedule. SPE and RFE keep the draw multiplicity of each power. Each query records the experiment that measures it and, on the trajectory, its measurement point, which the real and imaginary queries of one power share. When a point has no data, its queries remain missing, and neither analysis nor loading reruns part of the trajectory to replace them.

For a Pauli-sum Hamiltonian the controlled powers are product formulas. Exact static QPE trajectories use one ordered common Suzuki step on an integer time grid whose bounds meet every requested power's error budget, and the Plan records cumulative step counts and per-power errors. Sampled QCELS, SPE and RFE settings prepare separate circuits and choose a step count for each power independently, and the Plan records each power's count and block, with no common step. Both routes reuse one count of Pauli commutators. Every power must fit its controlled-evolution error budget (controlled_power_error_budget) after adding pruning, product-formula error at the represented step time, represented-time displacement, rotation-angle formation and identity-phase discrepancy. The exact trajectory counts its step increments once, while resource estimates for sampled Plans multiply each setting's own count by its requested shots.

On the exact trajectory, the nonnegative integer power p targets the exact time \(t_p=p\cdot\mathrm{val}(\tau)\), where \(\mathrm{val}(\tau)\) is the real number that the binary64 time unit \(\tau\) represents. A displayed evolution time may round this product. Each power records its cumulative step count and a bound that includes pruning, the product formula at the represented step time, represented-time displacement, and applicable phase and rotation-angle formation errors. Planning checks that subtotal at every queried power. The identity term acts on the control through phase increments targeting \(-c_I(p_k-p_{k-1})\,\mathrm{val}(\tau)\), because removing it would change the measured interference signal (Proposition 7). The circuit description contains the shared controlled evolution and its X and Y measurement points. Each sampled Hadamard setting ends with its own phase gate and final Hadamard before measurement, and circuit construction, inspection and resource totals use these same gates.

RWPE independently chooses a step count for each continuous power and checks the complete bound at the exact signed target time \(\mathrm{val}(p)\,\mathrm{val}(\tau)\). Its total_error is the upward rounding of the exact sum of pruning, product-formula error at the represented step time, represented-time displacement, rotation-angle formation and the discrepancy of the identity phase that the circuit uses. The step block's epsilon includes its own product-formula and rotation-angle terms. The complete subtotal must fit the controlled-evolution error budget, and planning fails without retry when this recheck fails. Circuit synthesis and execution roundoff remain outside this structural bound.

QCELS¶

How the fit works¶

QCELS fits \(z_p\approx ae^{-ip\theta}\) with \(\theta=E\tau\), and the raw-unitary phase is \(-\theta\). Eliminating \(a=\operatorname{mean}_p(z_pe^{ip\theta})\) (Ding–Lin arXiv:2211.11973v2, Eqs. (11)–(12)) reduces the complex least-squares objective in Ding–Lin Eq. (2) to one variable (Result 3). grid_size is a requested minimum. The effective grid G=max(grid_size, 8*(max(p)-min(p))+1) covers one phase period with at least eight points per period of the fastest objective term. Every discrete local minimum receives QCELS_BRACKET_EVALUATIONS=64 objective evaluations of golden-section search inside its adjacent periodic bracket, and the estimate is the best point among all grid points and all refined points. This finite search does not guarantee the global minimum.

When the accuracy theorems apply¶

Ding–Lin's accuracy results, Theorems 1 and 2 (arXiv:2211.11973v2, p. 14), assume \(p_0>0.71\), where \(p_0\) is the squared overlap of the preparation with the target eigenvector. NWQLib does not check \(p_0\). Smaller overlaps need their Fourier-filtered variant (Sec. IV, p. 15), which this method does not implement. NWQLib fits one schedule rather than the multi-level schedule of their Algorithm 1, whose times double as in their Eq. (35). Their data set, Eq. (6), uses the consecutive powers \(0,\ldots,N-1\). When the maximum power exceeds num_times, NWQLib spreads num_times integer powers over 1 through that maximum, a schedule their theorems do not describe.

Choose num_times and max_time¶

With max_time=None the maximum power is 10. A supplied max_time \(T\) gives the maximum power \(P=\lfloor T/\tau\rfloor\), or \(\lfloor T\rfloor\) for a unitary, and a max_time below one positive power is rejected. When num_times is at least \(P\), the schedule holds the consecutive powers 0 through \(P\). Otherwise it keeps power 0 and spreads num_times integer powers over 1 through \(P\).

The point estimator is dimensionless, so changing energy and time units together preserves its physical meaning for the same powers. The default integer schedule is independent of units. With a supplied max_time, the schedule is preserved when the floored dimensionless ratios agree. The strict floor has discontinuities, and binary64 unit conversion at an integer boundary can change the final power by one. Inputs immediately below that boundary must remain below it, since no near-integer snapping or decimal reinterpretation is applied.

A spread schedule can fit a replica of the energy

A spread schedule can make the objective nearly periodic, and the Result does not report that. For maximum power \(P\), its positive powers are \(p_j=\lfloor1+j\bar g\rfloor\) with \(\bar g=(P-1)/(N_t-1)\), where \(N_t\) is num_times. At \(\theta+2\pi/\bar g\) each term \(e^{ip_j\theta}\) turns by \(2\pi(1-f_j)/\bar g\), where \(f_j\) is the fractional part of \(1+j\bar g\), and the power-0 term does not turn. Every turn lies in \([0,2\pi/\bar g]\), about one radian for \(\bar g\) near 6, so the rotated terms stay nearly aligned and the objective there keeps most of the depth of the minimum at \(\theta\). With a mixed spectrum, finite shots or evolution error, such a replica minimum can win the fit and move the energy by \(2\pi/(\bar g\tau)\). On the H4 problem of the eigenvalue example notebook (Hartree–Fock reference, \(\tau=0.2\)/Ha, num_times=32), max_time=38/Ha gives \(P=190\) and \(\bar g=6.10\), and the fit lands 5.156 Ha above the ground energy, next to the replica spacing \(2\pi/(\bar g\tau)=5.153\) Ha, while 36 and 40/Ha stay within 7 mHa (Measured evidence). Consecutive powers 0 through \(P\), which the schedule uses when num_times \(\ge P\), have no such near period.

Aliasing report¶

If the powers with a nonzero summed signal share a common spacing \(g>1\), the objective has period \(2\pi/g\) and the Result reports the aliasing. The planned schedule always contains powers 0 and 1, so this happens only when the signal vanishes at some powers.

When a signal counts as zero¶

A power's summed signal counts as zero when it lies within the error that its exact means and the summation allow, sqrt(2)*mean_window*W plus the rounding of the sum for \(W\) summed samples. mean_window is the largest probability_window of the preparation records for exact probabilities and NUMERICAL_RELATION_RTOL for classical execution, and count data use an exact-zero test (Engineering constants). Exact data whose signal vanishes in exact arithmetic therefore report aliasing or a flat objective despite their roundoff, and RFE applies the same test before it calls a sampled Fourier spectrum flat.

Reanalyze the fit¶

import numpy as np
from nwqlib import Eigenproblem, solve
from nwqlib.algorithms.qpe import QCELS

result = solve(
    Eigenproblem(A=np.diag([.4, -.7])),
    method=QCELS(initial_state=[0, 1], tau=.2, num_times=16,
                 max_time=100),
    execution="classical", seed=7,
)
print(result.eigenvalue, result.fit.effective_grid, result.fit.evaluations)
refined = result.analyze(grid_size=8192)  # same observations
print(refined.eigenvalue)

-0.7000000000000006 4096 31168
-0.7000000000000006

The prepared state \((0,1)\) is the eigenvector of eigenvalue \(-0.7\). Post-hoc result.analyze(grid_size=...) changes the QCELS or SPE analysis grid using the existing observations, and no reanalysis measures new samples. Increasing grid_size is checked against the original method's limits. If the intended later analysis exceeds them, construct a new Plan with the needed limits before measurement. result.analyze(grid_size=8192) reuses the current observations within this example's max_work. RFE's K defines both the power distribution and the Fourier frequencies, so changing it requires a new Plan, as do sample counts, the initial state and every RWPE setting. seed records separate PCG64 streams for the method and the backend, and reanalysis does not draw from them.

SPE¶

SPE uses the Fourier filter of Wan, Berta and Campbell, arXiv:2110.12071v2, PDF Eqs. (A1)–(A3), corresponding to HTML Eqs. (16)–(18). Write \(F_j=-ia_j\) for its Fourier coefficient at each positive odd frequency \(j\), where \(a_j=|F_j|\), and let \(S=\sum_ja_j\). Drawing \(M\) frequencies \(j_1,\ldots,j_M\) with probabilities \(a_j/S\) and measuring their complex overlaps gives the estimate \(\widehat C(x)=1/2+(2S/M)\sum_{m=1}^M\operatorname{Im}(e^{ij_mx}z_{j_m})\) of the approximate cumulative distribution at the dimensionless threshold \(x\). For a Hamiltonian threshold energy \(E\), \(x=\tau E\) and \(z_j\) estimates \(\langle\psi|e^{-ij\tau H}|\psi\rangle\), with the Hamiltonian \(H\), prepared state \(\psi\) and time step \(\tau\) defined above. For raw-unitary input, \(x\) is the phase threshold and the formula uses the complex conjugate of the measured overlap. Negative frequencies follow by conjugation. Grouping repeated draws multiplies each exact or pooled contribution by its draw count. The terminal coefficient at frequency \(2d+1\) contains \(I_d(\beta)\) alone, where \(d\) is fourier_degree, \(\beta\) is filter_beta and \(I_d\) is the modified Bessel function of the first kind of order \(d\), as specified by PDF Eq. (A2), HTML Eq. (17). Algorithm 1 instead samples \(j\) over all of \(\{0,\pm1,\pm3,\ldots\}\) in proportion to \(\lvert F_j\rvert\). NWQLib keeps \(F_0=1/2\) exactly and pairs each positive draw with its conjugate, which also gives an unbiased estimate of the approximate cumulative distribution under unbiased overlap measurements.

The declared overlap lower bound \(\eta\) is an assumption supplied by the caller. Planning does not diagonalize an operator or expand a preparation to verify it. For the filter domain, Theorem 3 of Wan, Berta and Campbell, arXiv:2110.12071v2, uses \(\epsilon_1=\epsilon_2=\eta/8\) and \(\epsilon_3=\eta/2\) in PDF Eqs. (A6)–(A7) and (A12), HTML Eqs. (21)–(22) and (27). This gives a conditional filter error of at most \(3\eta/8\) away from a transition half-width \(\delta\), leaving \(\eta/8\) below the \(\eta/2\) decision threshold for sampling error. The Result reports \(\delta\) separately from the grid spacing. These analytic bounds are evaluated in binary64 and do not include coefficient roundoff or controlled-evolution error. Degree and beta are tested against Theorem 3, and a failed assumption appears in the Result.

NWQLib scans a fixed grid using one shared sample set, which permits reanalysis without measurement. The paper also reuses one sample set, but it makes only \(s=O(\log(1/\delta))\) decisions chosen as in binary search (Sec. III, p. 4), where \(\delta\) is the transition half-width defined above, and gives each decision failure probability \(\xi/s\), where \(\xi\) is the paper's total failure probability, so a union bound limits the total failure probability of Theorem 1 to \(\xi\). A scan over \(G\) grid points would need failure probability \(\xi/G\) per decision in Eq. (11). The default sample count is not derived from Eq. (11), so the search guarantee does not transfer. The paper's randomized LCU compiler and gate-complexity claim are also outside this controlled-evolution implementation.

RFE¶

RFE uses the randomized Fourier estimator of Kshirsagar, Katabarwa and Johnson, Sec. 2, Eqs. (5)–(7). A draw \(k\) contributes \(z_ke^{-2\pi ikj/K}/M\) to frequency \(j\). Grouping the identical draws of one power into one exact or pooled observation replaces their terms by the draw multiplicity times that observation's term. Uniform draws make the result an unbiased estimate of the finite Fourier coefficient. For an eigenstate with \(U|\psi\rangle=e^{i\phi}|\psi\rangle\), where \(\phi\) is the eigenphase angle in radians, the measured quadrature convention fixes the ideal signal \(z_k=e^{ik\phi}\). Postprocessing returns the estimate \(\widehat\phi\), the angle \(2\pi j/K\) of the frequency \(j\) of the largest-magnitude sampled coefficient reduced to \([-\pi,\pi)\). The reported result.phase is \(j/K\) turns, so \(\widehat\phi\) is \(2\pi\,\)result.phase reduced to \([-\pi,\pi)\), and the Hamiltonian energy estimate is \(-\widehat\phi/\tau\). The arXiv:2209.11322v3 note on p. 1 records that Eq. (4) and the output line of Algorithm 1 were corrected for a missing sign. With the paper's S gate before the final Hadamard, outcome 0 has probability \((1-\sin(k\phi))/2\), so the raw imaginary outcome has mean \(-\sin(k\phi)\) and the printed output is \(2\pi-2\pi j/K\). NWQLib's imaginary quadrature uses ancilla phase \(-\pi/2\) and measures \(+\operatorname{Im}z_k\) directly, so its output needs no reflection. Work is \(O(UK)\) for the U distinct drawn powers, with \(O(K)\) workspace. These sizes depend on the estimator budget, not the Hilbert-space dimension.

The default K gives a grid spacing of \(2\pi/49\) in \(\widehat\phi\), or \(1/49\) turn in result.phase, and maximum power 48. M sets the random-draw count. Its 194 quadrature draws become one query per distinct (power, quadrature) with both exact and sampled readout, and sampled readout requests 194 times shots coherent shots. The paper's precision theorem is Theorem 2.1 of arXiv:2209.11322v3, p. 6. For an eigenstate it gives an \(\epsilon\)-accurate phase estimate with probability above \(1-\delta_{\mathrm{RFE}}\), where \(\epsilon>0\) is the phase tolerance and \(0<\delta_{\mathrm{RFE}}<1\) is the theorem's failure probability, when \(K\ge\lceil2\pi/\epsilon\rceil\) and \(M\ge\lceil(81\pi^2/2)\ln(8\pi/(\delta_{\mathrm{RFE}}\epsilon))\rceil\) single-shot draws. For any \(\epsilon<\pi/2\) and \(\delta_{\mathrm{RFE}}<1\) that bound exceeds 1100 draws, so the default \(M=97\) is a finite workload without the theorem's guarantee. A successful finite Fourier fit establishes neither the eigenstate assumption nor the sample budget and supplies no RFE uncertainty interval. A dominant peak in a mixed spectrum need not represent its ground energy.

RWPE¶

RWPE implements the basic Gaussian update of Granade and Wiebe, arXiv:2208.04526v1, Eq. (7) and Algorithm 1. Its coordinate is the unwrapped phase \(\phi=-\tau E\). At step \(k\), \(\sigma_k=\sigma_0((e-1)/e)^{k/2}\), where \(\sigma_0\) is prior_std, the relative evolution time is \(1/\sigma_k\), and one bit \(d\) changes the mean by \((1-2d)\sigma_k/\sqrt e\). The physical evolution is \(e^{-i\tau H/\sigma_k}\), so a bare unitary with integer-power access is insufficient. All times are known from the prior and the step count, while the feedback depends on previous bits. The update uses constant scalar workspace at every Hilbert-space dimension.

The Hadamard circuit uses the feedback \(-\mathrm{mean}/\sigma-\pi/2\), giving \(P(0)=(1+\sin((\phi-\mathrm{mean})/\sigma))/2\). Integrating this likelihood against the Gaussian prior gives the positive \(d=0\) mean shift in Eq. (7a) of arXiv:2208.04526v1, because \(\mathbb E[x\sin x]=e^{-1/2}\) for a standard normal \(x\). In the paper's likelihood, Eq. (2), \(P(d|\omega)=\cos^2(t(\omega-\omega_{\mathrm{inv}})/2+d\pi/2)\), this feedback is \(\omega_{\mathrm{inv}}=\mathrm{mean}+\pi\sigma/2\). Algorithm 1 prints \(\omega_{\mathrm{inv}}=\mathrm{mean}-\pi\sigma/2\) alongside the same update, which gives \(P(0)=(1-\sin((\omega-\mathrm{mean})/\sigma))/2\) and would lower the mean after \(d=0\). At \(t=1/\sigma\) the paper's own exact update, Eq. (6a), p. 4, reproduces Eq. (7a) with \(\omega_{\mathrm{inv}}=\mathrm{mean}+\pi\sigma/2\) and gives the opposite sign with \(\omega_{\mathrm{inv}}=\mathrm{mean}-\pi\sigma/2\). The sentence that introduces Eq. (7) on p. 5 states \(\omega_{\mathrm{inv}}=\mathrm{mean}\), for which Eq. (6a) leaves the mean unchanged. An independent Gaussian-integral check fixes the sign from the circuit probability.

This implementation omits Algorithm 2's consistency checks and unwinding. As discussed after Eq. (8), the basic walk has a finite exploration range. The mean can move at most \(\sigma_0/(\sqrt e-\sqrt{e-1})\), about 2.96 prior widths (the paper prints 2.95), so a shrinking Gaussian width alone does not establish convergence. The Result reports when the scaled spectral radius and the prior mean reach that range.

RWPE's planned times do not determine its future feedback angles. A device-profile estimate (estimate(plan, profile=...), see Check device fit and run time) and a compare(...) row therefore give no forecast for the circuits with those angles, and a scan over Plans with different shot counts reports an objective that sums the forecasts of the individual circuits as unavailable, with the reason. These forecasts neither choose a feedback angle nor run the update.

Verify a result¶

result.verify(checks=QPEVerification(tolerance=...)) explicitly computes a nominal comparison of the spectrum and its eigenvalue clusters with a QR projector. SPE compares the lowest cluster and uses its declared overlap bound unless minimum_overlap is supplied. The other estimators compare the largest prepared cluster, with default minimum overlap 0.9. The comparison reports discrepancy, overlap, numerical adjustment, resolution and how many reference computations it performed or reused. It returns (receipt, facts), a verification record and the two checked facts, component_error against tolerance in the output unit and overlap_deficit = max(0, minimum_overlap - weight) against zero, for Certificate.with_verification (Check accuracy and verify a result). group_atol clusters in the unit of the compared eigenvalues, the Problem's energy unit for a Hamiltonian and radians of eigenphase angle for a unitary, whatever the output. reuse_only=True requires spectral and reference data already computed for the same target and initial state, including data reloaded after save and load, instead of computing them. The saved arrays themselves are not rechecked (Saved folders are read-only). A reference for a supplied circuit is computed only with materialize_preparation=True. These dense comparisons run only on request, within their own max_work and max_bytes limits, and never expand Pauli input automatically.

For the static estimators, method_resolution is the Result's grid spacing in output units. It includes a later supported analyze(grid_size=...) call and is not a fit-error bound. RWPE's method resolution is its current Gaussian standard deviation in output units. Neither quantity is a proven total-error bound, and verification leaves the original Result and Plan unchanged.

Continue and save¶

For RWPE, prepare builds only the circuit of the next query. Before each new step, RWPE saves in the Run its iteration state, which consists of the Gaussian mean, the number of processed steps, the next evolution time, feedback angle and backend seed, and the sequence number of the Run checkpoint of the pending step. submit and Run.wait advance this state on the backend of the Run. Each completed measurement updates the Gaussian exactly once. A submitted measurement whose outcome is not yet known is not submitted again, and after cancellation its recorded cost remains visible.

Before submission, the Run rejects a measurement whose fixed seed would repeat the sampling stream of a known measurement, whether completed, uncollected, uncertain or in the same batch. After completion, the job identifiers that the provider returned are checked as well. During recovery the same completed measurement can be read again, but it gives no second posterior update.

result.save(path) and load_result(path) keep the concrete method, the evolutions it chose (including the polar factor of a unitary input), the observations, the Gaussian moments and the existing numerical caches. Loading rejects a conflicting method before it loads state data. Saved formats lists the archive format names and versions. run.save(path) and load_run(path, backend=...) keep the Run's progress, the random-number state and the circuit caches. Loading performs no schedule search, eigensolve, estimator replay or circuit construction, and reanalysis reuses the unchanged data.

result.analyze(...) measures nothing new. For RWPE it reads the Gaussian mean and the number of processed steps from the saved run state, and its width follows Granade and Wiebe, arXiv:2208.04526v1, Eq. (7b), without replaying observations. Without saved run state, the RWPE estimate is reported as missing. Its stop reason is the one the run recorded when it ended.

Measured evidence¶

Measurement and section Settings and environment

QCELS replica on H4 (Choose num_times and max_time) The H4 chain of the eigenvalue notebook: four hydrogen atoms 2.0 Å apart, STO-3G, all eight spin orbitals active, its dense 256-by-256 Hamiltonian and the Hartree–Fock basis state. QCELS(tau=0.2, num_times=32, max_time=36, 38 or 40) with the default grid_size, execution="classical" and seed 7. The errors are relative to the CASCI (full CI) energy, −1.8977806 Ha: −6.4 mHa at 36/Ha, 5.156 Ha at 38/Ha and −0.073 mHa at 40/Ha. Python 3.12.14, NumPy 2.5.2, SciPy 1.18.1, PySCF 2.14.0 and Qiskit 2.5.2 on macOS arm64.

Limits¶

A single-mode estimate establishes neither the ground energy nor coverage of the dominant component, and the preparation determines which part of the spectrum is observed.
QCELS, SPE and RFE return no interval, and RWPE's interval is a nominal Gaussian model interval (Intervals and sampling error).
The papers' accuracy theorems need conditions that NWQLib does not check or does not meet by default (Choose an estimator). QCELS's theorems assume \(p_0>0.71\), which is not checked, and consecutive powers. SPE's search guarantee assumes the paper's binary-search decisions, not a grid scan. RFE's theorem needs more than 1100 draws, against a default of 97.
A spread QCELS schedule can fit a replica of the energy, and the Result does not report it (Choose num_times and max_time).
The basic RWPE walk has a finite exploration range and omits Algorithm 2's consistency checks and unwinding (RWPE).
Dense powers are explicit finite constructions, not scalable oracles. Under the default max_work, each dense block fits on up to seven system qubits.
The default sampling assessment stays INCONCLUSIVE, because NWQLib attaches no bound from the papers' sample-count conditions.

Limitations and open work lists the open items for these methods.

Numerical scope and planning limits¶

max_power=100000 and max_trotter_steps=1000000 cap finite evolution work. max_work=1000000000 bounds the work units counted for numerical routines, schedules and the RWPE iteration state, where a work unit is one scalar operation in a count that planning derives from the problem sizes by formula before the computation runs. max_bytes=10_000_000_000 bounds the arrays whose sizes are known in advance. Neither bounds the CPU time or memory of LAPACK or Qiskit. RWPE's update counts 2 units when it starts and 32 for each choice and each update, and its last update follows the last measurement, so planning checks this whole count, 2+64*max_steps units, against max_work before the first query. On its product-formula route, RWPE also checks the independent exact recheck of all max(1, max_steps) relative times during planning. That recheck, the commutator count and the RWPE update each have their own comparison against max_work, so the update total alone does not cover every planning stage. On resume, an interrupted choice or update is counted again. The ExecutionLimits shared by all methods cover circuit preparation, measurement and data. Symbolic evolutions and repeated steps remain compact until the Qiskit circuits are built explicitly.

The QCELS analysis limit check counts at most (1 + QCELS_BRACKET_EVALUATIONS)*G = 65*G objective evaluations, the \(G\) grid points plus one bracket for each possible local minimum. Each evaluation counts 32 scalar size units per complex point, so the QCELS formula is \(32\cdot65\cdot GN\) units plus 64 per scalar quadrature sample for the reduction. Here \(N\) counts distinct complex points after quadrature reduction. These are conservative operation-size units, not measured CPU instructions. Grid and scores use \(O(G)\) storage and row temporaries \(O(N)\), without a \(G\times N\) matrix. Limits are checked before the schedule is built or measured and again before reanalysis. The default has \(N=11\), \(G=4096\) and 22 quadrature queries, giving \(32\cdot65\cdot4096\cdot11+64\cdot22=93{,}717{,}888\) units. A limit of 90 million rejects this analysis count before measurement, while 100 million and the default max_work=1000000000 accept this check. Other planning stages have separate limit checks.

When a limit check fails, the error reports the required work and bytes and the supplied max_work and max_bytes. NUMERICAL_RELATION_RTOL is a dimensionless tolerance of the library's input checks, defined in nwqlib._validation, and not an accuracy setting of the method.

Commutator count and common step¶

Before counting commutators, exact-trajectory product-formula planning sets aside an upper bound on common-step work and logical bytes from the powers, term count and input exponents. It uses the full Pauli-triangle expression when that combined work and every stage's bytes fit. Otherwise it uses the second-order degree-capped suffix relaxation, with pair tests and weighted degrees only. The Plan names the expression it used (Proposition 36). For \(L\) terms on \(q\le64\) qubits, let \(N_{\rm pair}=L(L-1)/2\) and \(N_{\rm nested}=L(L-1)(2L-1)/6\). Packed tests and magnitude visits give initial work bounds \(2Lq+2(N_{\rm pair}+N_{\rm nested})\) for full and \(2Lq+5N_{\rm pair}+2L\) for relaxed. At the explicit budget max_work=100000000 and q=20, the count alone allows at most 531 full or 6316 relaxed terms. Sampled settings and RWPE set aside no common-step bound. Their exact recheck has its own work limit before the arithmetic runs, using the formula in Engineering constants.

On the product-formula route, default sampled schedules contain 11 distinct powers for QCELS, at most 12 for SPE, and at most 49 for RFE. For \(K\) powers, including zero, and \(L\) kept nonidentity terms, the independent recheck accepts 32L(K+1)C_T + (320K+128)C_S + 128K_a C_I <= max_work, where \(K_a\) counts nonzero powers when \(L\) is positive and \(C_T\), \(C_S\) and \(C_I\) are the term, scalar and candidate-inversion width costs defined in Engineering constants. With costs bounded by \((1,27,100)\), a max_work of 100000000, a tenth of the default, accepts this recheck throughout the default commutator-count term range, including schedules through 64 powers. At costs exactly \((8,27,100)\) and that limit, an RFE schedule containing all 49 allowed powers accepts at most 7,731 terms at this stage. Wider arithmetic can reduce that limit further. The coefficients, times and pruning mass determine the widths, so \(K\) and \(L\) alone do not guarantee acceptance. The commutator count and the other planning requirements remain separate checks.

An exact-trajectory Plan must also fit the common-step bound determined by its power schedule, coefficient exponents, pruning loss, identity coefficient, time scale, per-power error budgets and step cap, so these commutator-count figures are not its term limits. Above the full-path threshold, an accepted relaxed initial bound can still lead to a full count after the pair structure is known. Other planning work, storage and step limits also apply.

The exact trajectory's common-step selection checks its exact arithmetic against the limits before constructing trial step counts. Its work count uses the exponent ranges of the inputs and includes trial counts that max_trotter_steps later rejects. It computes the kept coefficient mass and the one-step rotation-angle discrepancy once, then adds each identity phase that the circuit uses once to the exact running prefix. The arithmetic reductions take \(O(L+K)\) visits for \(L\) kept terms and \(K\) requested powers. The limit check also includes sorting, a conservative \(K^2\) count for integer-key collisions, and an integer-width factor. Both commutator-count checks include the same common-step bound in max_work, whose default is 1_000_000_000, so a commutator-count variant that passes its check also fits the later exact-arithmetic work limit. The full expression is preferred when the combined count fits, and otherwise planning uses the relaxation when it passes. Each prefix is accepted only when its complete exact subtotal fits its error budget. max_bytes separately limits the logical size of the objects held in memory at once during the exact stage.

Dense powers¶

Planning forms the polar factor V once per Plan from one SVD, counted as \(9D^3+8D^2\) units of planning work, and the classical route and verification diagonalize V, not A. The exact trajectory constructs one block per distinct positive gap between consecutive powers, and a sampled Plan one block per distinct nonzero power. A gap of k powers is constructed by repeated squaring of V and built as a controlled dense block. Each gap's limit check covers its matrix work, its synthesis and the cached matrices held in memory. Planning counts each block's matrix arithmetic in \(D^3\) units per product of two \(D\)-square matrices, where \(D\) is the dimension of the system register. A unitary gap or power k counts the \(\mathrm{bit\_length}(k)+\mathrm{popcount}(k)-2\) products of an independent repeated squaring, although the squares of V are shared within the Plan and a later block forms only the squares it lacks. Its bytes include every cached square, \(16JD^2\) for the largest square index \(J\) of the Plan's exponents. A Hamiltonian power takes one product and also counts the eigendecomposition, \(8D^3+32D^2\) units, because the first power that is constructed computes it and the others reuse it. Planning adds \(11\cdot8^m+(m^2+5m+256)\cdot4^m\) units for the synthesis on \(m=n+1\) qubits with \(n\) system qubits. Each block is checked against max_work on its own, so the default max_work accepts dense blocks on up to seven system qubits. The Plan's total over its blocks is not compared with max_work.

Source and code map¶

The code is in src/nwqlib/algorithms/qpe/:

method.py: method settings, original inputs, the circuit description, readout and static reconstruction.
records.py: queries, powers and results, with their input, schedule and error-allowance associations.
powers.py: the signed-power numerical and circuit constructors.
numerical.py: the QCELS, SPE, RFE and RWPE scalar routines and their source conventions.
controller.py: RWPE's iteration state between steps, its Gaussian mean and pending measurement point, saved in the Run.
archive.py, verification.py: saving and loading of the data and caches, and the nominal reference comparison of result.verify.
subroutines/trotterization/error_budget.py: commutator bounds and the exact choice of an integer step count for a given evolution time. On the common grid, QPE uses those bounds to choose sufficient cumulative counts of the shared step and to check the represented steps against them.
subroutines/qpe/coherent.py: the separate coherent phase-register circuit.

The QPE methods build their controlled product-formula steps in powers.py from Qiskit Pauli evolution gates. They use only the bounds of the Trotterization subroutine, not its circuit builders.

Each row below names the scientific step, paper version, equation and code in algorithms/qpe unless a subroutine path is given. Here z_p is the Hadamard-test signal at evolution multiplier p, and sigma is RWPE's Gaussian phase width. The docstrings of the listed code derive implementation-specific formulas and state their assumptions.

Step	Source	Location	Code
Real and imaginary Hadamard tests (W = I or S†) give `z_p`. The exact trajectory reads the same quadratures as ancilla X and Y at each power position	Ding–Lin, arXiv:2211.11973v2	Fig. 1, p. 3. Eqs. (4)–(5), p. 8	`method._quantum_program`, `method._trajectory_samples`, `method._static_estimate`
QCELS objective and eliminated amplitude	Ding–Lin, arXiv:2211.11973v2	Eq. (2), p. 4. Eqs. (9)–(12), p. 9	`numerical.qcels`
QCELS search	Ding–Lin, arXiv:2211.11973v2, with NWQLib grid	Sec. II.2, p. 10, suggests up to ceil(T_max) grid points with gradient ascent. NWQLib uses eight cells per unit of power span and golden-section brackets of `QCELS_BRACKET_EVALUATIONS` evaluations	`numerical.qcels`, `numerical.qcels_grid_size`
QCELS schedule and validity	Ding–Lin, arXiv:2211.11973v2	Eq. (6), p. 8. Algorithm 1, p. 13, and its doubling times, Eq. (35), p. 14, not implemented. Theorems 1–2, p. 14, need `p0 > 0.71`. Sec. IV, p. 15, not implemented	`numerical.planned_power_schedule`, `method.QCELS`
SPE random schedule and filter	Wan–Berta–Campbell, arXiv:2110.12071v2	PDF Eqs. (A1)–(A2), p. 7, HTML Eqs. (16)–(17), Algorithm 1. Positive odd frequencies sampled by coefficient magnitude, with paired quadratures and conjugate negative-frequency contributions	`numerical.spe_fourier_weights`, `numerical.planned_power_schedule`
SPE estimator	Wan–Berta–Campbell, arXiv:2110.12071v2	Eq. (2), p. 2. Eq. (6), p. 3. Algorithm 1, p. 4. First CDF crossing of eta/2 over a finite grid of decisions, in place of the O(log(1/delta)) binary-search decisions of Sec. III, p. 4. Theorem 3, p. 10, gives the filter domain. Controlled-evolution synthesis error is accounted for separately	`numerical.spe_cdf`, `numerical.spe`, `numerical.spe_filter_domain`
RFE signal	Kshirsagar–Katabarwa–Johnson, arXiv:2209.11322v3	Eqs. (1)–(6), p. 4	`numerical.planned_power_schedule`, `method._static_estimate`
RFE estimator	Kshirsagar–Katabarwa–Johnson, arXiv:2209.11322v3	Algorithm 1 and Eq. (7), p. 5. Uniform powers in [0,K), sampled Fourier coefficients, largest magnitude. Sign note, p. 1. Sample bound of Theorem 2.1, p. 6	`numerical.rfe`
RWPE likelihood	Granade–Wiebe, arXiv:2208.04526v1	Eq. (2), p. 2	`numerical.rwpe_choose`, `numerical.rwpe_update`
RWPE experiment choice	Granade–Wiebe, arXiv:2208.04526v1	Algorithm 1, p. 3. Time 1/sigma and one Bernoulli datum. Feedback sign follows Eq. (2) and the positive d=0 shift in Eq. (7a)	`numerical.rwpe_choose`, `controller.run_rwpe`
RWPE posterior and estimate	Granade–Wiebe, arXiv:2208.04526v1	Eqs. (5)–(7), pp. 4–5. Gaussian moment updates and mean estimate. Algorithm 2 unwinding is absent	`numerical.rwpe_scale`, `numerical.rwpe_update`, `numerical.rwpe_estimate`
Second-order Suzuki step and term order	Childs et al., Phys. Rev. X 11, 011020 (2021), doi:10.1103/PhysRevX.11.011020	Prop. 10, p. 011020-22 (arXiv:1912.08854v3 Prop. 16, p. 39)	`powers.rotation_schedule`
Commutator bounds	Childs et al., doi:10.1103/PhysRevX.11.011020 (Phys. Rev. X version)	Prop. 9, Eq. (120). Prop. 10, Eq. (121). p. 011020-22 (arXiv:1912.08854v3 Sec. 5.1, Prop. 15, Eq. (145), p. 38, and Prop. 16, Eq. (152), p. 39)	`subroutines/trotterization/error_budget.py`
Step-count selection	Childs et al., doi:10.1103/PhysRevX.11.011020, with NWQLib's common-grid specialization	Sec. V B for the independent minimum-count rule. The common-grid derivation gives sufficient prefix counts, each checked against the represented-step bounds	`powers.select_powers`, `subroutines/trotterization/error_budget.py`
CX of a controlled Suzuki step	NWQLib, from Qiskit 2.6.0 `qiskit/circuit/library/pauli_evolution.py::PauliEvolutionGate.control` and `crates/synthesis/src/pauli_evolution.rs` (Proposition 7)	Basis change, parity chain of w-1 CX and its inverse around a controlled RZ of 2 CX, derivation in the docstring	`powers.suzuki_step_cx`, `powers._make`
Pruning budget, controlled identity phase, automatic tau, non-power-of-two padding	NWQLib	Docstrings of the listed code	`method._split_pauli_terms`, `powers.select_powers`, `method._quantum_program`, `method._safe_time`, `method._QPEMethod.plan` (Pauli input), `method._QPEMethod._inputs`
Classical nominal signal	NWQLib	One dense eigendecomposition and the prepared spectral weights give the noiseless ancilla means	`method._bind_host`, `method._host_construction`, `numerical.overlap`
Coherent phase register and inverse QFT	Cleve–Ekert–Macchiavello–Mosca, arXiv:quant-ph/9708016v1	Sec. 5. Fig. 6 and Eq. (5.1), p. 10. QFT of Eq. (4.1), p. 8. Readout and its 4/pi**2 bound, Eqs. (5.2)–(5.4), p. 11. The builder docstring maps the paper's bit order to Qiskit's	`subroutines/qpe/coherent.py`

Independent checks cover two-level signed quadratures, the controlled global phase, physical-unit conversion, finite likelihoods, non-power-of-two coordinates, rejection of a reused sampling stream and recovery of an interrupted RWPE run. Circuit-level checks of the QPE methods use at most four qubits. The separate coherent QPE subroutine check uses seven, six phase qubits and one system qubit.