QHD

Quantum Hamiltonian descent (QHD) is a quantum algorithm for minimizing an objective over a box, from Leng, Hickman, Li and Wu, Quantum Hamiltonian Descent, arXiv:2303.01471v1. NWQLib evolves a state on a finite grid of the box under a grid form of the Hamiltonian of their Eq. (1), written here as H(t) = a(t) K + b(t) V, with K the kinetic operator -Delta/2, V the objective and a(t), b(t) the weights of a schedule, and reads out grid points.

QHD reports two points, which can differ. The most probable point is the grid point of largest readout probability, up to a numerical tie window. The candidate is the observed valid point with positive weight that has the least evaluated objective. valid_probability is the unconditional probability of the outcomes that encode a grid point. An observed candidate does not establish a continuous or global optimum.

Use QHD for an objective written as a SymPy expression with bounds on every variable. For d variables with K grid points each, classical evaluation stores K**d amplitudes and a one-hot circuit uses d*K qubits, or d*log2(K) with the binary encoding. Equality and inequality constraints go through the augmented-Lagrangian layer, and box refinement repeats QHD on shrinking boxes. Terms defines the words this guide uses in a narrow sense.

import sympy as sp
from nwqlib import solve
from nwqlib.problems import Optimization
from nwqlib.algorithms.qhd import QHD

x = sp.Symbol("x", real=True)
problem = Optimization(objective=(x - 0.2)**2,
                       variables=(x,), bounds=((-1.0, 1.0),))
result = solve(problem, method=QHD(), seed=7)
analysis = result.analyze()
print(analysis.most_probable_coordinates, analysis.most_probable_objective)
print(analysis.candidate, analysis.objective, analysis.valid_probability)

(0.33333333333333326,) 0.017777777777777767
(0.33333333333333326,) 0.017777777777777767 1.0

The default grid of [-1, 1] has the two interior points -1/3 and 1/3. The most probable point and the candidate are both 1/3, the grid point nearest the minimizer 0.2, and valid_probability is 1.0.

For this example the defaults run a two-qubit Aer circuit with two grid points per variable, one second-order step, total time 1, the quadratic schedule with gamma 0.3, midpoint coefficients and the kinetic ground state as initial state, prepared by the structured one-hot chain. Use QHD(num_grid_points=3, num_steps=4) to change the finite scientific model. A positive shots selects full-register counts. Absent shots uses exact probabilities. keep_state=True explicitly saves exact amplitudes for a later fidelity comparison. Counts cannot supply those amplitudes.

Using QHD

• Read the most probable grid point, where QHD concentrated probability, from most_probable_coordinates and most_probable_objective of result.analyze(). The constrained and refinement layers read this point by default. A numerical tie window can select a point just below the largest computed probability.
• Read the best observed point, the grid point with the least evaluated objective among those with positive observed weight, from candidate and objective. With exact readout it has the least evaluated objective over the whole grid when every grid point has positive weight, and with counts it is the best sampled grid point.
• Choose the grid type with boundary and include_boundary_points before setting num_grid_points, the number K of points per variable. One-hot encoding uses d * K qubits for d variables, and encoding="binary" uses d * log2(K) and requires a periodic grid with K a power of two. Classical joint arrays have K**d entries, while quantum simulation depends on the encoded qubit count.
• Use total_time to change evolution duration, include_boundary_points=True for the endpoint grid and boundary="periodic" for the periodic grid (the default is the interior Dirichlet grid), schedule with coefficient_rule to select the schedule weights, and initial_state to select the start. These change the finite scientific model, and so does kinetic_model="spectral" on the periodic grid. encoding="binary" and binary_synthesis change the circuit that approximates the model. The QHD API documents each configuration field and its work and byte limits. No option implies a continuous global-optimization guarantee.
• For a classical run, pass execution="classical" to solve and choose QHD(theory_flavor=...). The default "schrodinger" evolves the finite-difference grid model directly, and its work grows with the schedule weights and the objective's range. "split_step" has work per step that depends on neither but adds a splitting error, so check how its result changes as num_steps grows. "ir_product" is the numerical reference for the product formula that the circuit applies.
• After a work or byte refusal, raise each exceeded max_work or max_bytes to its reported amount. A complete amount clears that check. An amount reported as “at least” is a lower bound, because counting may have stopped early, and a limit set to it can refuse again, so raise that limit further, for example to twice the amount. Instead of raising a limit, or after another kind of refusal, follow the remedy that the message gives, and recheck the model and its accuracy when the remedy changes the grid, schedule or objective scale. Explicit verification has separate limits in QHDVerification.
• Add constraints by passing a ConstrainedOptimization to solve_augmented_lagrangian, and repeat QHD on shrinking boxes with refine_box or the constrained solver's refinement option. Read feasibility, complementarity and termination alongside the objective. Neither procedure establishes a continuous global optimum.
• Save a Result with result.save(path) and reopen it with load_result, or with load_augmented_lagrangian for a constrained result and load_box_refinement for a standalone refinement result. To continue constrained or refinement work after an interruption, run it with directory=path and resume it with resume_augmented_lagrangian or resume_box_refinement. Resuming depends on the recorded inner Run being recoverable.

Three notebooks apply QHD. examples/qhd_optimization_intro.ipynb applies it to a nonconvex objective of two variables, and its Section 7 adds an inequality constraint on the unit disk and solves it with the augmented-Lagrangian layer. examples/qhd_scientific.ipynb runs box refinement on the shifted Ackley function and the augmented-Lagrangian layer with refinement on a constrained Rastrigin function of Wu et al., arXiv:2605.12066v1, on binary periodic grids, and reports the solution quality, the shots needed to find a good point and the cost. Section 4 of examples/resource_estimation_at_scale.ipynb, "QHD: one-hot versus binary encoding", compares the CX count, the arbitrary rotations and the evolution bound of one step of a three-variable quadratic in both encodings, up to a 96-qubit one-hot register at K = 32 (fault-tolerant resources).

Each measured comparison, default, timing and bound figure in this guide gives its settings where it appears, and Measured evidence gives the environment and the remaining settings and limits.

plan(problem, method=QHD(...)) computes the objective's support tables and the schedule's step weights once, without importing Qiskit, building a circuit or running a reference evolution. The returned Plan holds the chosen construction and its costs, computed before any circuit exists, together with its original Problem, Method, output and readout. prepare(plan).circuits holds the built Qiskit circuit. result.analyze() reuses the same observations, and saving and loading use the stored tables and phases without evaluating the objective or replaying the schedule.

A Plan fixes its Method, output, accuracy, execution, shots and seed, and passing a non-None value for any of these to solve(plan, ...) raises. Call plan again to change them. Backend, progress and Run limits remain execution controls.

The QHD circuit uses one measurement setting named qhd. A counts measurement with S repetitions counts as S shots, one setting and zero exact evaluations. Exact probabilities or amplitudes count as zero sampled shots, one setting and one exact evaluation. A kept-amplitude readout is read from the experiment's amplitude declaration.

Per-circuit resources describe one prepared circuit. Shot-weighted resource totals also depend on the requested repetitions. Zero sampled shots for an exact evaluation does not specify a hardware repetition budget, so a hardware cost that needs that budget can remain unknown with a reason.

Scientific model and ordering

The Hamiltonian follows Leng, Hickman, Li and Wu, Quantum Hamiltonian Descent, arXiv:2303.01471v1, Eq. (1), with the finite-difference stencil and diagonal potential of Eqs. (F.7)/(F.9). The one-hot map is described by Wu et al., arXiv:2605.12066v1, Sec. IV.A, Eqs. (9)–(10). NWQLib combines it with the schedule and coefficient rule described in Schedules and coefficient rule. Leng's Algorithm 1 uses radix-2 registers, Fourier kinetic evolution and left-endpoint weights, while Eq. (F.36) uses symmetric Hamming-weight states. Those encoding-specific resource bounds do not describe the one-hot dK-qubit register or NWQLib's binary circuit, whose costs Binary encoding derives. The method module docstring lists each departure from Leng et al.'s Algorithm 1 with its reason.

boundary and include_boundary_points select one of three grids. Each has K points per variable, with the local index i = 0..K-1 increasing with the coordinate, and every part of QHD reads the same coordinates, spacing and links from grid.OneHotGrid.

Grid	Points x_i	Spacing h	Links per variable
Dirichlet interior, the default	lower + (i+1) h	(upper-lower)/(K+1)	K-1, a chain
Dirichlet endpoints, include_boundary_points=True	lower + i h, so x_0 = lower and x_(K-1) = upper	(upper-lower)/(K-1)	K-1, a chain
Periodic, boundary="periodic"	lower + i h, with upper excluded	(upper-lower)/K	K, a cycle

The table gives the coordinates in exact arithmetic. grid.OneHotGrid.grid_value evaluates them in binary64, so on the endpoint grid the last coordinate lower + (K-1) h can differ from upper by rounding in either direction, for example 0.9999999999999998 on [-1, 1] with K = 50 and 0.7000000000000001 on [0.1, 0.7] with K = 38. Planning requires the computed coordinates of every grid to be finite and strictly increasing, the first interior coordinate to lie above lower, the first periodic or endpoint coordinate to equal lower, and the last periodic or interior coordinate to lie below upper. It checks only that the last endpoint coordinate is finite, because rounding moves it off upper for many ordinary boxes (grid.OneHotGrid).

Each variable contributes kinetic diagonal 1/h² and hopping -1/(2h²) on each of its links. Variable j occupies qubits j*K through (j+1)*K-1. In Qiskit bit strings qubit 0 is the rightmost character, and restricted arrays use lexicographic ordered grid-index tuples. The caller's variables and bounds determine original coordinate order.

Interior points represent the free amplitudes of a homogeneous Dirichlet problem, since the two endpoint amplitudes are fixed at zero. The centered stencil approximates -psi''/2 with local error O(h²) for smooth data. This spatial choice is independent of the time discretization of the schedule. The endpoint grid reproduces the endpoint mesh j/r, j=0..r, of arXiv:2303.01471v1, Eqs. (F.5) and (F.9), with K=r+1. The paper states the vanishing boundary condition with Eq. (F.4) but keeps the endpoint amplitudes in the (r+1)-square matrices of Eq. (F.9). On the endpoint grid the endpoint amplitudes evolve as part of the truncated chain and the missing neighbors lie outside the declared box, so only the interior grid imposes the vanishing boundary condition, by keeping the zero endpoint values outside the register. Use the interior convention when zero amplitude at the box endpoints is the intended boundary condition.

The periodic grid treats the box as one period of length L = upper - lower. The point upper is the same point as lower, so the grid excludes it, and a wrap link joins x_(K-1) to x_0. Each variable's kinetic matrix is then circulant, with 1/h² on the diagonal and -1/(2h²) at (i, i±1 mod K), and the circuit applies the wrap link's XX+YY hopping with the same coefficient -1/(4h²) as the other links. For variable j that gate acts on qubits j*K and j*K+K-1, which are neighbors only on hardware with ring connectivity, and the Plan's CX count includes it before routing. This is the one-hot kinetic term of Liu et al., arXiv:2607.16996v1, Eq. (12). The stencil approximates -psi''/2 with local error O(h²) at every point for a smooth L-periodic function. For an objective that is not periodic on the box, the wrap link still moves probability directly between x_0 = lower and x_(K-1) = upper - h. Every row of the circulant matrix sums to zero, so the uniform initial state is its eigenvector with kinetic energy zero.

The periodic grid rejects include_boundary_points=True, and with the one-hot encoding it requires an even K of at least 4. The second-order one-hot product splits each variable's links into two layers of disjoint links. The wrap link has index K-1, which falls in the odd layer without meeting another odd link only when K is even. This is a limit of the current link compiler, not of the periodic Laplacian. K = 2 is excluded as a degenerate cycle, whose wrap link would repeat the link between points 0 and 1. The rule also holds for first order and for classical execution, so a one-hot Method accepts the same periodic grids for its quantum circuit and for its classical reference. The binary encoding uses no links and accepts K = 2 (Binary encoding).

Planning expands the objective in monomials and groups them by the variables they contain, one support table per group plus a constant. It expands about the point of each box nearest the origin, which is the origin itself whenever the box contains it. Expanded about the origin, (x - 1e7)**2 on a box of width 2 near 1e7 would store a table and a constant near 1e14 each, and their binary64 sum would keep the objective only to about 0.016. About the box point, the terms of that example stay below the squared box width. A box away from the origin can therefore produce support groups that the plain expansion does not have, such as single-variable tables from a product term.

Planning evaluates each support group into a finite float64 array. Those stored values define the grid objective. Grid readout and refinement scores read the corresponding entries. The augmented-Lagrangian layer also checks the original supplied functions at every point it compares, and explicit grid verification evaluates the original objective independently. These checks use the original functions and can differ from the tables by numerical evaluation error.

Step k of num_steps covers [k dt, (k+1) dt] with dt = total_time/num_steps and applies the kinetic and potential weights that the coefficient rule selects. First order uses potential then forward hopping. Second order uses potential halves around odd-link half, even-link full, odd-link half kinetic steps. Disjoint links in a parity commute. The fused XX+YY kernel keeps the coefficient sign. Compact number projectors implement exp[-i angle*(product(n)-I/2**support)]. The omitted physical identity phase is restored once. The step block declares no support for arbitrary external control or for an adjoint. On an 8-by-8 periodic grid of a two-variable quadratic with a cosine term, with the quadratic schedule, midpoint weights and T = 1 (measured evidence), the second-order product had the higher fidelity to theory_flavor="schrodinger" with the same weights at every step count from 8 to 64 in both encodings, for example 0.99996 against 0.911 for the one-hot encoding at 64 steps. The one-hot first-order fidelity did not rise monotonically over that range (0.72, 0.32, 0.59 and 0.91 at 8, 16, 32 and 64 steps), so check step refinement on the chosen model. A comparison at equal step count does not compare equal gate costs.

rotation_threshold defaults to zero, preserving every computed nonzero angle. An explicit positive value drops rotations strictly below that angle in radians, every compiled hopping or projector block on the one-hot encoding and only the Walsh-string Rz rotations on the binary encoding, whose QFT controlled phases and dense diagonals stay, while aqft_cutoff truncates the QFT instead (circuit synthesis). The reconstruction stores dropped count, raw absolute-angle sum and physical phase sources, and pruning_error_bound, an operator-norm bound on how far pruning moves the emitted product. A dropped one-hot block differs from the identity by at most its angle, so for one-hot the bound is the exact angle sum, rounded upward once. The bound also includes the contributions that planning omits below the normal binary64 range (normal binary64 range). None of these is an objective-gap bound. Nonfinite and non-real objective values reject before pruning, including small nonzero imaginary parts. Zero objectives remain legal.

execution="classical" explicitly evaluates the restricted Schrodinger model. QHD(theory_flavor="ir_product") evaluates the numerical product reference, and QHD(theory_flavor="split_step") the symmetric split-step product of the split-step evaluation. These routes use Kd amplitudes, whereas exact circuit readout has 2(dK) one-hot probabilities, or the same K*d for the binary encoding. A classical evaluation of the model is not quantum execution, and by itself it is not an independent baseline. theory_flavor="schrodinger" numerically evaluates the exponential of each step's Hamiltonian with the selected step weights, which for the midpoint rule is the exponential midpoint rule. It has no splitting error. For a time-dependent schedule it keeps the time-discretization error of the coefficient rule and is not the exact time-ordered evolution. For QuadraticSchedule(gamma=0.0) the Hamiltonian is constant, so there is no time-ordering or coefficient-quadrature error, but the numerical evaluation keeps its floating-point and solver error.

theory_flavor="schrodinger" applies SciPy's expm_multiply, whose roundoff lets the total probability of the unitary evolution differ from one by an amount that grows with each step's generator norm. On the one-hot encoding, theory_flavor="ir_product" applies each stored block's analytic action at its stored angle, a phase on the projector's grid slice or a rotation of the two slices that a hopping link joins (theory._run_ir_product), and on the binary encoding it has the state-operation budget on computed phase arrays and selected stored QFT angles. The kernel therefore records a probability_window in its KernelApplication, derived from those norms or that budget as the engineering constants describe, and the reported masses must lie within it. The window is a first-order bound under the assumptions that its docstring states about SciPy's parameter choice, not a measured error. The split-step kernel derives its window from its own state budget. Numerical guarantees states the binary and split-step budgets and the budget that each route uses.

theory_flavor="schrodinger" and theory_flavor="split_step" evolve the potential of the support tables alone. The objective constant c adds b(t) c I to the Hamiltonian, which commutes with every operator and so contributes only the global phase exp(-i c sum_k dt b_k) over the step weights. The kernel reads the probabilities before it multiplies the kept state by that phase, so a constant such as 1e16, near which adjacent binary64 numbers are 2 apart, cannot round the objective's variation away. With keep_state=True planning rejects a constant whose phase c sum_k dt b_k exceeds the binary64 range, such as c = 1e308 over two unit steps, or whose phase-error budget (Numerical guarantees), evaluated upward, reaches pi, because the budget then gives no guarantee on the kept state's phase (method._admit_constant_phase). Without keep_state the same objective is accepted. The quantum circuit and the IR-product reference apply the compiled identity phase, the one-hot sum of the kinetic diagonal, the projector identities and the constant, or the constant alone for the binary encoding, before every readout, so there planning rejects a phase beyond the binary64 range for every readout. The probabilities are accepted for any finite phase. The candidate, the objectives and the records keep the constant. Numerical guarantees gives the budgets behind these phase rules.

Schedules and coefficient rule

The Hamiltonian at time t is H(t) = a(t) K + b(t) V, with K the kinetic operator -Delta/2 and V the objective. QHD.schedule selects the weights a and b as one of three records. Each record holds only its own parameter, and its kind field names the formula and is part of the Method's content hash and of the saved Method.

Record (kind)	a(t)	b(t)	Parameter	Source
QuadraticSchedule ("quadratic"), the default	1/(1+gamma*t²)	1+gamma*t²	gamma >= 0, default 0.3	QHDOPT, arXiv:2409.03121v1, Sec. 2.1
CubicSchedule ("cubic")	2/(s+t³)	2t³	s > 0, no default	Leng et al., arXiv:2303.01471v1, Eq. (C.4). Liu et al., arXiv:2607.16996v1, Eq. (92), write this form with s = 1 as their QHD-C schedule, but their code runs the shifted form of the next row under that name (circuit synthesis)
ShiftedCubicSchedule ("shifted_cubic")	(2/(s+t))³	2t³	s > 0, no default	Wu et al., arXiv:2605.12066v1, Eq. (15), a replication option. Liu et al.'s code runs it as qhd-c, and with s = 1 only this form reproduces the Rz counts of their Table II (circuit synthesis)

For gamma > 0 the quadratic ratio a/b = (1+gamma*t²)**(-2) tends to zero as t increases, satisfying the ratio condition after Eq. (1) of Leng et al. arXiv:2303.01471v1. The value 0.3 is an illustrative working point also used in the optimization notebook, and QuadraticSchedule(gamma=0.0) selects the undamped, time-independent model. The cubic schedule follows the ODE model of Nesterov's accelerated gradient method, with s removing the singularity of 2/t³ at t = 0 (Leng et al.). Eliminating the momentum from the classical Hamiltonian a|p|²/2 + b f(x) gives x'' - (a'/a) x' + a b ∇f = 0. The cubic schedule has damping 3t²/(s+t³) and gradient factor 4t³/(s+t³), which tends to 4. The shifted cubic schedule has damping 3/(s+t) and gradient factor 16t³/(s+t)³, which tends to 16. Wu et al. present their Eq. (15) as the QHD-C schedule of Leng et al., but it defines different dynamics from Leng's Eq. (C.4), so a result for one schedule does not carry over to the other.

The parameter s belongs to the model, not to the time step. Refining num_steps keeps s fixed, because changing s changes the Hamiltonian. Leng et al. set s to their time step 0.001. The code behind the published results of Wu et al. counts N_t time points, so its time step is T/(N_t - 1) = 10/49999, and it sets s = T/N_t = 10/50000 = 2e-4, close to that step. A reproduction therefore passes that value of s explicitly. No schedule name carries a convergence guarantee for a finite grid or evolution time.

coefficient_rule selects each step's weights. "midpoint", the default, uses the point values a(t_k) and b(t_k) at the step midpoint t_k = (k+1/2) dt. "integrated" uses the step averages A_k/dt and B_k/dt of the interval integrals A_k of a and B_k of b over the step, so the step's kinetic and potential exponents equal A_k and B_k, the first Magnus term of H(t). Second order then applies B_k/2 in each potential half and A_k/2, A_k and A_k/2 to the odd, even and odd link layers. Exact coefficient integrals remove the quadrature error of the coefficients but not the time-ordering error, which is nonzero when K and V do not commute and a/b varies in time. The error model therefore lists midpoint_time for the midpoint rule and time_ordering for the integrated rule. The two rules differ most where a varies quickly. On the first interval [0, s] of the shifted cubic schedule the kinetic integral is 3/s², while the midpoint value s a(s/2) is 64/(27 s²), about 21% lower. This compares coefficients, not final states. The left-endpoint weights of Leng et al.'s Algorithm 1 would leave every route first order in the step, even with symmetric placement, so no left-endpoint rule is offered (split-step evaluation gives the Magnus term that such weights miss).

The time-ordering error budget takes the per-step minimum of the available derivative and integral bounds.

In the convergence check of the split-step evaluation, with the quadratic schedule and both cubic schedules at s = 1 (measured evidence), both rules were second order and their errors differed by less than 1%, so the midpoint rule, which needs no interval integrals, remains the default. The check does not cover a small s, where the coefficients differ most and where the midpoint sum of a depends strongly on the step count (initial states). The interval integrals avoid the cancellation of their closed forms when a step is short compared with its start time, and the docstrings of the schedules module derive their relative error bounds, at most 48 units of roundoff on the documented parameter range.

The classical schrodinger flavor makes one expm_multiply call per step, and planning counts the work of each call from the norm of its generator, so the work of one evolution grows with sum dt a and sum dt b. Both cubic forms raise it, because their kinetic weight is largest near t = 0 when s is small and b = 2t³ grows with the total time. Planning checks the total against max_work before any evolution. theory_flavor="split_step" counts a fixed amount per step instead, independent of a, b and the objective's range. Planning limits gives both counts.

The Schrodinger limit check depends on the generator norm, including the kinetic scale and the relative-potential range. A smaller time step can reduce a finite step norm, but it cannot repair an earlier overflow while forming the weighted kinetic or potential terms. A coarser grid or wider box can reduce the kinetic scale. For the shifted cubic schedule, increasing s lowers the kinetic weight. These changes affect the model.

At fixed grid, tables, start and step count, the declared split-step kernel work does not increase with the potential's amplitude or the schedule weights. Its numerical error budget can still change. Switching theory_flavor therefore changes both the computational procedure and its approximation.

Split-step classical evaluation

QHD(theory_flavor="split_step") with execution="classical" evolves the K**d grid amplitudes with one symmetric (Strang) product per step, exp(-i dt b_k V/2) exp(-i dt a_k K) exp(-i dt b_k V/2), where (a_k, b_k) are the step weights of the coefficient rule. Under "integrated" the kinetic exponent is the step integral A_k and the two potential halves are equal, B_k/2 each, which needs only the whole-step integrals. The kinetic factor is applied exactly in the eigenbasis of the kinetic operator, one transform pair per variable along that variable's axis of the lexicographic array. This is the symmetric second-order form of the pseudo-spectral step of Leng et al., arXiv:2303.01471v1, Eq. (C.3), which applies the potential and then the kinetic factor with left-endpoint weights and is first order. It starts from the Method's initial_state, as the other theory_flavor values do.

kinetic_model selects the kinetic operator K. "finite_difference", the default, is the stencil that every route applies. On both Dirichlet grids its eigenvectors are the sine vectors sqrt(2/(K+1)) sin(pi (j+1) r/(K+1)) of orthonormal DST-I, with eigenvalues (2/h²) sin²(pi r/(2(K+1))) for r = 1..K, where each grid uses its own spacing h and the endpoint grid keeps its endpoint amplitudes dynamical. On the periodic grid the Fourier modes diagonalize the circulant stencil with eigenvalues (2/h²) sin²(pi k/K) (Result 37). "spectral", available on the periodic grid only, is the Fourier operator of -Delta/2 with eigenvalue 2 pi² q²/L² for period L = K h and the signed index q = k for k < ceil(K/2) and k - K otherwise, the order of SciPy's fftfreq. For K = 8 the modes carry q = 0, 1, 2, 3, -4, -3, -2, -1. The kernel calls scipy.fft.dst(type=1, norm="ortho"), which is its own inverse, and scipy.fft.fft and ifft with norm="ortho". The Dirichlet boundary rejects the spectral model, because its Fourier modes do not vanish at the box ends. The one-hot circuit, the one-hot theory_flavor="ir_product" and theory_flavor="schrodinger" reject it, because they apply the finite-difference stencil, while the binary circuit and the binary theory_flavor="ir_product" apply the same eigenvalues between a QFT and its inverse (Binary encoding).

The two kinetic models are different operators on the same grid. For |q| <= K/2 their energies satisfy 0 <= E_sp - E_FD <= 2 pi⁴ q⁴ h²/(3 L⁴) (Proposition 43), so at a fixed physical momentum the two energies agree to O(h²), while at the Nyquist mode their ratio is pi²/4 and the norm of their difference grows as 1/h². For the same initial state, Duhamel's formula gives ||psi_FD(T) - psi_sp(T)|| <= ∫ |a(t)| ||(K_FD - K_sp) psi_sp(t)|| dt, so the difference between the two models depends on the momentum content along the whole trajectory, not only on the initial or final positions. The error model of a spectral Plan lists the source kinetic_model.

The time error has three separate parts, for a fixed grid, objective and schedule parameter and weights a and b with four bounded derivatives in time. First, the midpoint rule's frozen point values miss the first Magnus correction -i (dt³/24)(a'' K + b'' V) at the step midpoint, which the integrated rule's exact integrals remove. A left-endpoint rule would miss -i (dt²/2)(a' K + b' V) and be first order even with symmetric placement, so it is not offered. Second, the Magnus term Omega_2 = (dt³/12)(a b' - a' b)[K, V] + O(dt⁵) remains under both rules, so exact coefficient integrals do not make a noncommuting problem exact in time. Third, the symmetric product adds the Baker-Campbell-Hausdorff terms (i A² B/12)[K, [K, V]] + (i A B²/24)[V, [K, V]] of the splitting, which theory_flavor="schrodinger" does not have. Each step therefore has local error O(dt³), and the global error at fixed total time is O(dt²) under both rules.

The constants depend on commutators and schedule derivatives, so they are not uniform as s tends to zero, the grid is refined or a penalty grows. A nonsmooth objective, such as a PHR penalty or Ackley's cusp, does not change the order in dt of this finite model, but error estimates that assume a smooth potential do not apply to it. Refining num_steps keeps the schedule parameter s fixed, which a convergence check of one Hamiltonian requires. For [K, V] = 0, which on a grid means an objective that is constant on the grid points, the integrated rule is exact at any step count. On the one-dimensional two-mode cosine of Liu et al., arXiv:2607.16996v1, on an eight-point periodic grid with the kinetic ground state, T = 1 and s = 1 for both cubic forms (measured evidence), the split-step error against an independent adaptive ODE solution fell by a factor of about 4.0 from 128 to 256 steps for all three schedules under both rules, and the errors of the two rules differed by less than 1% at each step count. The time-ordering and splitting terms that the integrated rule keeps have the same local order dt³ as the quadrature error it removes.

The split-step work per step does not depend on a, b or the objective's values. Planning limits gives its work and byte counts.

A larger penalty can increase the range of the augmented objective and, on the Schrodinger path, the planning work counted from its generator norm. The range need not grow monotonically because the objective, linear multiplier term and quadratic penalty can cancel. For the two-dimensional Ackley function on [-5, 5]² mapped to the unit square, with K = 32 on the default Dirichlet interior grid, 1000 midpoint steps and total time 10 (measured evidence), planning counts 19 times as many units of work for theory_flavor="schrodinger" as for split-step under the quadratic schedule and 1.9 million times as many under ShiftedCubicSchedule(s=2e-4), while the split-step count stays at 4.7e7. These are ratios of nominal units of planning work, which the Plan records as estimates with SciPy's internal work unknown, not measured runtimes or operation counts.

Compare theory_flavor="schrodinger" and theory_flavor="split_step" at a common final-state error, because they need different step counts. theory_flavor="schrodinger" has no splitting error, while theory_flavor="split_step" applies its kinetic factor through transforms whose work per step does not depend on the penalty or the schedule weights, but its commutator terms can require more steps as a penalty grows.

One comparison used the two-dimensional Dirichlet interior grid with K = 32 or 64, the potential (x-1)² + (y-1)² + (rho/2) max(0, x² + y² - 1)² on [-1, 1]² with rho from 1 to 512, T = 1, the kinetic ground state, integrated coefficients of the quadratic schedule (gamma 0.3) or of either cubic schedule with s = 1, step counts from 16 to 2048, and max_work and max_bytes raised to 1e10 (measured evidence). The targets were phase-aligned state errors of 1e-3 and 1e-4. theory_flavor="schrodinger" often reached a target with fewer steps, but split-step reached it with less counted planning work and less wall time in every case. At 1e-3 the ratio of Schrodinger to split-step wall time was 2.0 to 9.0 at K = 32 and 4.3 to 34 at K = 64. It was at least ten in 10 of the 12 K = 64 cases at 1e-3, in 5 of 12 at 1e-4 and in none at K = 32, and a larger penalty could lower the ratio, because split-step then needed more steps. These were single timings on one machine.

With ShiftedCubicSchedule(s=2e-4) and the same raised limits (measured evidence), the limit check refused theory_flavor="schrodinger" before evolution at every tested step count from 32 to 2048, while split-step at K = 64 and rho 512 needed about 16,384 steps to reach an error near 1e-4. These measurements give no general crossover rule. Check step refinement on the chosen problem, particularly before reusing a Schrodinger step count with split-step.

For example, the augmented-Lagrangian run on the constrained Rastrigin problem of Wu et al., arXiv:2605.12066v1, with K = 32 on the Dirichlet interior grid of [-5, 5]², T = 10, the quadratic schedule with midpoint coefficients, the kinetic ground state and three levels of refinement in each round (measured evidence), stopped after one round at objective 33.2 with 200 split steps. With 3,200 split steps it reached objective 7.97, the point that theory_flavor="schrodinger" reaches with 200 steps, and 6,400 steps changed neither the point nor the number of rounds and levels.

The split-step kernel's mass window and the upper limit of its tie window (tie_window_ceiling) come from the Plan's state budget (split_step.state_error), and the kernel selects its most probable point with the smaller budget that split_step.evolve observes on the state's own trajectory. Numerical guarantees derives both.

On a Dirichlet grid the uniform state has probability in excited kinetic modes, which pick up phases of about E_k ∫ a dt while a(t) is large, as Initial states and preparation explains. The midpoint rule approximates that integral by sum dt a(t_k), whose value changes strongly with the step count under a cubic schedule with small s, whereas the integrated rule uses the exact integral. On the periodic grid the uniform state is the zero-momentum eigenstate, whose kinetic energy is zero in both kinetic models.

Initial states and preparation

QHD.initial_state selects the state that every route starts from, and QHD.initial_state_preparation selects how quantum execution prepares it. Each state is a product of nonnegative per-variable amplitude vectors on the grid points of that variable. Planning evaluates the vectors and the error bound of the classical start vector once and stores both in the Plan (QHDReconstruction.initial_amplitudes and initial_state_error), so solving and analysis do not evaluate them again. Quantum execution prepares each stored vector in its variable's register. The classical routes start from the tensor product of the stored vectors in lexicographic grid order for a general product state, and fill 1/sqrt(K**d) directly for the uniform state and the kinetic ground state on the periodic grid, whose construction error of 2u is then its own start term.

Record (kind)	Amplitude at grid point i of variable j, before normalization	Source
UniformState() ("uniform")	1	Leng et al., arXiv:2303.01471v1, Algorithm 1 step 3
KineticGroundState() ("kinetic_ground"), the default	sin(pi (i+1)/(K+1)) on the Dirichlet grids, 1 on the periodic grid	Ground state of the finite-difference kinetic operator, derived in the record's docstring
GaussianState(center=c, widths=sigma) ("gaussian")	exp(-(x_ji - c_j)²/(2 sigma_j²)) at the grid coordinate x_ji	Warm start. Wu et al., arXiv:2605.12066v1, Sec. VI, state a uniform start, while the code behind their published refinement results starts every layer after the first from this state, centered at the best point found so far, which BoxRefinement(level_initial_state="best_point_gaussian") offers

On both Dirichlet grids the kinetic operator of one variable is T/h², where T is the K-by-K tridiagonal matrix with diagonal 1 and off-diagonals -1/2. The interior grid and the grid with boundary points therefore share the ground-state vector, whose energy is 2 sin²(pi/(2(K+1)))/h², and on the grid with boundary points it does not vanish at the box endpoints. On the periodic grid the one-variable operator is the circulant C/h² with C = I - (S + S^T)/2 for the cyclic shift S. Its Fourier eigenvalues 2 sin²(pi r/K)/h² vanish only for r = 0, whose eigenvector is constant, so there KineticGroundState() is the uniform state with kinetic energy zero. For several variables the ground state is the product of the per-variable vectors. Result 37 gives these spectra.

The center of a Gaussian is a point in the coordinates of the solved problem, which are the original coordinates for solve and the unit coordinates of each level box in a search-model box refinement, and it may lie outside the box. Its widths are in the same coordinates. The amplitude, not the probability, has width sigma_j, so the continuous probability density, proportional to exp(-(x_j - c_j)²/sigma_j²), has standard deviation sigma_j/sqrt(2) along variable j. Any finite center and any finite positive widths are accepted. On the periodic grid the Gaussian uses the chart distance |x_ji - c_j| between coordinates of the solved problem's box [lower, upper), not the distance on the circle, so for a center near one end of the period the points near the other end, its neighbors through the wrap link, get the small amplitudes of their chart distance (Limitations and open work).

Each variable's exponents (x_ji - c_j)²/(2 sigma_j²) are shifted by their minimum before exponentiation, so the largest amplitude is exactly 1 before normalization and the state never underflows to zero, however far the center lies from the box. If every exponent of a variable exceeds the binary64 range, the grid point at the smallest distance from the center gets the whole amplitude of that variable, and points at equal rounded distances share it. The construction error bound of the start vector is finite for every input. It grows with the exponents, about 10u times the exponent per entry. For a center far outside the box with steeply decaying amplitudes, the error is measured against the largest entry, which stays exact, and the bound stays near 7u per variable. When the rounding of the exponents cannot be bounded, or when rounded distances tie, the bound becomes the distance of two nonnegative vectors of norm about one, about sqrt(2), and the tie window of the classical kernel then covers every point.

The choice matters most when the kinetic weight a(t) is large at early times, as for the cubic schedules with small s. On a Dirichlet grid the uniform state is not an eigenstate of the kinetic operator. Its probability outside the kinetic ground state is 1 - 2 cot²(pi/(2(K+1)))/(K(K+1)) for one variable, 0.165 for K = 32 and 0.177 for K = 64, with limit 1 - 8/pi² ≈ 0.189, and for two variables it is 0.303 and 0.323. While a(t) is large, each kinetic eigencomponent of energy E picks up a phase of about E ∫a dt. Under the midpoint rule the discrete value of this integral depends strongly on the step count. For the shifted cubic schedule with s = 2e-4 and total time 10, the midpoint sum sum dt a is about 6.0e5 with 1000 steps and 2.6e7 with 10000 steps, while the integral is about 1.0e8. A result from the uniform state then depends on the time discretization through this component. The kinetic ground state is an eigenstate of the kinetic operator, so on a Dirichlet grid it has no such component. On the periodic grid the uniform state is itself the kinetic ground state, so either record starts there without this dependence. A Gaussian contains several kinetic eigencomponents and does not share this property. It concentrates the initial probability near a chosen point, which suits a warm start such as a later layer of a box refinement.

KineticGroundState() is the default. It is the ground state of the kinetic term, whose weight relative to the potential, a(t)/b(t), is largest at t = 0, and it has no excited kinetic component whose phase depends on the step count. Under the quadratic schedule both terms are present at t = 0, so it is not an eigenstate of H(0), and the choice carries no adiabatic or global-optimization guarantee. On the periodic grid it is the uniform state, and on the Dirichlet grids its structured one-hot preparation has the same 3(K - 1) CX per variable as the uniform state's.

In a comparison with the uniform state on two-variable Dirichlet grids, with the classical Schrodinger model, the quadratic schedule and exact readout (measured evidence), the kinetic ground state gave the smallest best-point distance in all six standalone refinement runs with the most probable point, on a double well, an anisotropic quadratic and the Ackley function, and a much better refined constrained Rastrigin point, distance 0.0013 against 0.034 with the most probable point. It is not better everywhere. The refined run with the most probable point on a mixed-scale problem took 5 rounds instead of 3, and the unrefined mode_or_mean run on the unit disk ended farther from the optimum, distance 0.023 against 0.00079. theory_flavor="split_step" with 3,200 steps gave the same refined constrained Rastrigin distances, 0.0013 against 0.034, so theory_flavor="schrodinger" and theory_flavor="split_step" share this default. UniformState() remains available for the uniform start of Leng et al.

initial_state_preparation="structured", the default, prepares the initial state with the library's structured circuit for the Method's encoding and initial state, and the Plan records its CX count. Each encoding supplies its own construction. On the one-hot encoding it is a linear chain in each register. With amplitudes alpha_i and remaining norms r_m = sqrt(sum_(i >= m) alpha_i²), an X gate excites the register's first qubit, and link m applies a controlled RY with cos(theta_m/2) = alpha_(m-1)/r_(m-1) followed by a CX that moves the excitation to qubit m. For the uniform state this is the standard linear W-state construction. The chain stops after the register's last positive amplitude, because every later link would act as the identity. It also stops before the first link whose rotation parameter theta_m/2 would fall below the normal binary64 range, and the omitted links are counted in the state_preparation error source (initial_state.chain_selection, fault-tolerant resources). It emits three CX per link, at most 3d(K-1) in total, and the Plan's CX count covers exactly the emitted links (Proposition 53). "qiskit_state_preparation" appends one Qiskit StatePreparation per register with the one-hot embedding of that register's amplitudes, a vector of 2**K entries, and the Plan then records no CX count. "none" selects a quantum resource-only construction without preparation.

Binary encoding

QHD(encoding="binary", boundary="periodic", num_grid_points=2**b) stores grid index n_j = sum_l 2**l n_(j,l) of variable j in the b qubits j*b to j*b+b-1, bit l on qubit j*b+l. The register has 2**(d*b) = K**d states and each of them is a grid point, so invalid_mass is zero and valid_mass is one up to roundoff. The encoding requires the periodic grid, whose kinetic operator the quantum Fourier transform (QFT) diagonalizes, and K a power of two, the number of states of a b-qubit register. Any b ≥ 1 is accepted. At K = 2 each point's two neighbors are the other point, so the finite-difference stencil is (I - X)/h**2 with energies 0 and 2/h**2, and one H conjugation applies it. A Dirichlet kinetic term would need a sine-transform circuit, which the binary encoding does not have. Qiskit's statevector index of a grid tuple is z = sum_j n_j K**j, with variable 0 in the low-order qubits, while the classical kernels use the lexicographic index i = sum_j n_j K**(d-1-j). The permutation P|i> = |z> reverses the order of the variable axes and keeps the bits within each variable, and analysis, the classical product and verification convert between the two orders with it.

The binary encoding changes the circuit and not the finite model. Each step applies the potential and kinetic factors in the order of the one-hot compiler, first order U_K exp(-i B V) and second order exp(-i B V/2) U_K exp(-i B V/2), and under the integrated coefficient rule both potential halves apply B_k/2. The kinetic factor of variable j is F^dagger diag(exp(-i alpha E_k)) F on its register, where F is Qiskit's QFT |n> -> K**(-1/2) sum_k exp(2 pi i k n/K) |k>. Qiskit's sign is opposite to that of the forward discrete Fourier transform, but both kinetic models below depend on k only through min(k, K - k), so E_(-k mod K) = E_k and this conjugation is the same operator as the classical Fourier-space evolution (binary.kinetic_table). With exact QFTs and no pruned rotation the factor applies each variable's periodic kinetic operator exactly. The one-hot product instead splits the kinetic term over overlapping links, so the two encodings give different circuits at a finite step count, and each is compared with its own product formula.

For d variables with K grid points each, one-hot encoding uses d*K qubits. Binary encoding uses d*log2(K) qubits and requires K=2**b with b>=1 on a periodic grid. Both have K**d valid grid states. The full one-hot Hilbert space also contains invalid states, whereas binary encoding uses every basis state for a grid point. Changing boundary conditions or grid size changes the finite problem.

A local quantum simulation checks its encoded qubit count against max_simulation_qubits before allocating the readout state space.

Kinetic models

The kinetic factor applies the periodic eigenvalues E_k of kinetic_model in Fourier-index order, the same values as the split-step evaluation: E_k = 2 sin**2(pi k/K)/h**2 for "finite_difference", the eigenvalues of the circulant stencil, so the binary circuit and the one-hot circuit approximate the same model, and E_k = 2 pi**2 q_k**2/L**2 of the signed index q_k for "spectral". That section bounds the difference between the two models. The spectral model runs on the binary circuit, the binary theory_flavor="ir_product" and theory_flavor="split_step".

Liu et al., arXiv:2607.16996v1, Sec. IV E, introduce the quadratic kinetic phase to replace the dense kinetic diagonal by one- and two-bit terms, and study it numerically in Sec. V C. Their Eq. (89) approximates the finite-difference eigenvalue by the continuum form, and their Eqs. (90)–(91) expand k**2 with the unsigned index k = sum_l 2**l k_l. The unsigned square assigns (K - 1)**2 to the index K - 1, which is momentum -1, for example 49 instead of 1 at K = 8, and it breaks the symmetry E_(-k) = E_k that makes the QFT conjugation equal to the Fourier-space evolution. The low-momentum form of Eq. (89) is the signed one, so NWQLib implements only the signed index. Its weights w_l = 2**l for l < b - 1 and w_(b-1) = -2**(b-1) give q**2 = (1 + sum_l w_l**2)/4 + (1/2) sum_l w_l Z_l + (1/2) sum_(l<m) w_l w_m Z_l Z_m, still b single-Z and b(b-1)/2 two-Z strings (Proposition 43).

Circuit synthesis

BinarySynthesis holds the circuit choices. qft_bit_reversal="relabel", the default, omits the swap layers of both QFTs. The swap-free QFT is W = R F, with R the bit reversal of the register, and for any diagonal D, W^dagger (R D R) W = F^dagger D F. Applying the phase table v[rev_b(k)] between the swap-free QFT and its exact inverse therefore gives the same kinetic factor with 6 floor(b/2) fewer CX, and it leaves no permutation on the position register (Proposition 44). "swap" keeps Qiskit's swap layers.

aqft_cutoff=m keeps the controlled phases between wires at distance r <= m, whose angles are pi/2**r, and drops the rest, the approximate QFT of Coppersmith (arXiv:quant-ph/0201067v1). This is Qiskit's approximation_degree = max(0, b - 1 - m), and None, the default, keeps the exact QFT. One truncated QFT differs from the exact one by at most e_F = min(2, sum_(r=m+1..b-1) (b - r) 2 sin(pi/2**(r+1))) in operator norm, a matched pair around an exact diagonal by at most 2 e_F, and the product of N_s steps with one conjugation per variable by at most min(2, 2 N_s d e_F) (Proposition 44). The reconstruction records that total as aqft_error_bound, and the error model lists aqft_truncation. At b = 6 the bound 2 e_F of one conjugation is 0.196 for m = 4, 0.980 for m = 3 and the cap 2 for m <= 2.

Each support table S becomes a real diagonal on |S| b qubits, addressed by z_S = sum_t n_(j_t) K**t over its variables in increasing order. potential and kinetic_phase select how each diagonal is synthesized.

• "dense_diagonal" passes the phases -x v[z] of exponent x to the exact phase-diagonal construction of Shende, Bullock and Markov (quant-ph/0406176v5, Theorem 7, p. 10, whose multiplexed Rz with k select bits costs 2**k CX by the count after their Theorem 8, p. 11), 2**n - 2 CX on n qubits whatever the values, and none for an all-zero table.
• "walsh_rotations" expands the table in Pauli-Z strings with the Walsh coefficients c_m = 2**(-n) sum_z (-1)**popcount(m & z) v[z], as Liu et al. do in their Sec. IV A (Eqs. (53)–(54)), and applies one parity-ladder Rz(2 x c_m) per string, 2 (w - 1) CX for weight w, the rotation of their Eqs. (57)–(58) and the CX count of the sentence after them. The identity coefficient becomes a global phase. The kinetic tables use their exact Pauli expansions, (b - 1) 2**(b-1) CX for finite difference, whose nonzero strings all contain the top bit, and b(b - 1) for the signed spectral square. A support table's coefficients come from a fast Walsh–Hadamard transform of its rounded values, so a coefficient that is zero for the exact polynomial can come out as a rounding residue of order u times the table's scale. Such a coefficient is not a proved zero and is emitted unless rotation_threshold removes it.
• "min_cx", the default, compares the two counts for each table at the exponent of each block and takes the smaller, Walsh rotations on a tie. For the kinetic phase this selects Walsh rotations for the spectral model at every b and for finite difference at b = 2, and the dense diagonal for finite difference from b = 3. It minimizes the CX count of each table separately. It does not minimize Rz or T counts, depth, routed CX on limited connectivity, or the CX of a jointly synthesized potential whose strings share parity networks.

b	QFT, relabel	QFT, swap	Dense kinetic phase	Finite-difference Walsh	Spectral Walsh
1	0	0	0	0	0
2	2	5	2	2	2
3	6	9	6	8	6
4	12	18	14	24	12
5	20	26	30	64	20
6	30	39	62	160	30

The table gives CX counts at Qiskit optimization level 0. One kinetic factor costs two QFTs plus its phase diagonal, and an AQFT cutoff m replaces the controlled-phase part of a QFT by 2 sum_(r<=m) (b - r). rotation_threshold drops the Walsh rotations whose angle magnitude is strictly below it, records their number and absolute angle sum, and records pruning_error_bound, half that sum, because ||Rz(theta) - I|| <= |theta|/2. The bound sums the dropped binary64 angles exactly, halves the sum exactly, adds the exact error bounds of the rotations and dense phase entries omitted below the normal binary64 range and rounds the result upward once, so it is never below the exact value and is exactly zero when nothing is removed. The threshold never prunes the dense diagonal.

For a two-variable objective whose one coupled table has all 1023 Walsh strings, such as the normalized Ackley function of Liu et al. at K = 32, Walsh rotations need 8194 CX per potential factor and the dense diagonal 1022, and at K = 64, with all 4095 strings, 40962 and 4094. Because min_cx takes the cheaper of the two counts for each table, its CX count never exceeds that of Walsh rotations with the same computed coefficients and pruning. Fewer CX gates do not imply fewer T gates or a lower physical cost, which also depend on the rotation angles and the allocation of the synthesis error (fault-tolerant resources). Kinetic factors and preparation add their own costs, and a second-order step applies the potential twice. A polynomial table has few nonzero strings in exact arithmetic but rounding residues on the others, so Walsh rotations pay off there only with a threshold that removes the residues.

Liu et al.'s Table II counts use Walsh-rotation potentials, dense kinetic phases, swap layers and first-order steps. These options with rotation_threshold=1e-12, ShiftedCubicSchedule(s=1), total_time=10 and 100 steps on the periodic box [0, 1]² (measured evidence) reproduce the per-step binary CNOT and Rz counts of their Table II for the four two-dimensional problems at N = 32 and 64. The Rz count is the number of emitted rotations farther than 1e-12 from a multiple of pi/2, averaged over the 100 steps. Table II prints the integer part of that average, for example 1,189 for their Ackley function at N = 32, where the average is 1189.72, the authors' 100-step total of 118972 in their analysis notebook (experiments/res_ana/20260613_nisq_res_ana.ipynb of their repository) divided by 100. The paper's Sec. VI A and the caption of Table II name the QHD-C schedule, which their Eq. (92) defines as the cubic schedule, while the authors' code applies the shifted cubic form under the name qhd-c. With CubicSchedule(s=1), the form of Eq. (92), the CX counts are the same, but the Rz average is 4.02 per step lower for each N = 32 row and 7.52 lower for each N = 64 row, because the kinetic weight decides which dense kinetic-phase angles lie within 1e-12 of a multiple of pi/2.

With the same settings the defaults relabel and min_cx give the same Rz counts as these options under either schedule and need fewer CX per step for all four problems (measured evidence), for example 1162 instead of 8358 for their Ackley function and 230 instead of 254 for their coupled quadratic at N = 32. Without the threshold, rounding residues in the two-variable table of the coupled quadratic make min_cx choose the dense diagonal for that table, 1222 CX per step at N = 32 (measured evidence), so a comparison with a published count must state its pruning rule.

The CX count of a binary Plan is the sum of its recorded block counts, the count of the circuit transpiled to CX and single-qubit gates at optimization level 0 before routing. Its rotation count includes three phase gates per controlled phase, one Rz per Walsh string and the 2**n - 1 multiplexor angles of a nonzero dense diagonal. Both counts are upper bounds for later optimization, and angle-specific Clifford or T replacement can only lower the rotation count. The structured preparation, initial_state_preparation="structured", is one H per qubit, which prepares the uniform grid state exactly and with no CX. On the periodic grid the default kinetic ground state is the same state. A Gaussian initial state needs initial_state_preparation="qiskit_state_preparation", one Qiskit StatePreparation of K amplitudes per register, and the Plan then records no CX count.

The classical evaluations theory_flavor="schrodinger" and theory_flavor="split_step" evolve the same finite model for either encoding.

Binary ir_product applies the selected QFT operations and computed phase arrays as a numerical reference for the emitted product. For a Walsh diagonal, reconstructing the phase array from stored rotation angles introduces floating-point error. The first-order floating-point budget bounds the subsequent state operations relative to the computed arrays and stored QFT angles, under its arithmetic assumptions. It does not include every error from forming those arrays or the Walsh reconstruction discrepancy from the exact emitted rotations.

IR verification compares against this numerical reference. The current tie window and verification fidelity decision do not propagate the Walsh reconstruction bound R_W into an uncertainty interval. A small reported discrepancy therefore does not, by itself, prove fidelity to the exact emitted circuit.

The binary reference acts on the K**d lexicographic state, reading qubit j*b+l as bit l of the variable-j axis. The schrodinger reference evolves the finite-difference stencil, so schrodinger_fidelity of a spectral Plan compares two kinetic models, and its relation says so.

Readout

For probabilities or amplitudes, the candidate minimizes evaluated values over the positive-weight valid outcomes present in the readout. When every grid point has positive observed weight, it is a full-grid evaluated minimum. Otherwise the comparison guarantees minimization only over the observed positive-weight subset. Exact readout removes shot sampling, but does not make objective evaluation or evolution exact.

For counts, the candidate is the best sampled valid point and the mode is the most frequent valid point, with lexicographic tie handling. Pooling counts compares integer frequencies. It does not provide a bound on the mode of the underlying distribution. valid_count records the number of returned outcomes that decode to a grid point and returned_shots the number of shots returned, both as integers summed over all observation chunks (ObservationChunk, the stored result of one measurement). When returned_shots is positive, valid_mass is their quotient rounded once. Both are None under exact readout.

An attached exact-readout summary leads with the computed probability maximizer (the probability_maximizer_* fields), followed by the tie representative (the most_probable_* fields) with its deficit and window and by the mode status (mode_status), and then gives the candidate. An attached counts summary leads with the best sampled valid candidate and follows it with the same three lines. Its status line reads Empirical mode status and states that counts give no proof of the population mode. A Result without its Plan leads its summary with the candidate, described as the least evaluated objective among observed valid points, followed by the three lines. The status line appears even when the deficit is zero. If no positive valid outcome is available, one line gives the unavailable reason in place of both points.

Read the mode through most_probable_indices, most_probable_coordinates, most_probable_probability and most_probable_objective, with the numerical tie rule below.

Ties are decided among the valid points with positive observed probability, the same set from which the candidate is chosen. Analysis first finds the computed maximum probability M and then selects the lexicographically smallest grid-index tuple whose probability p satisfies M - p <= most_probable_tie_window. Every point is compared with M, because near-equality between neighbors is not transitive. most_probable_deficit records the computed M - p of the selected point, so a positive value means that a point within the window was chosen over the computed maximum. Other points can lie within the window even when the deficit is zero, so the selection does not identify a unique mode. The window bounds how far a computed difference of two probabilities can lie from the difference in the model that the path evaluates. Numerical guarantees derives the window from a state budget, and State budgets gives the budget of each route.

probability_maximizer_indices, probability_maximizer_coordinates, probability_maximizer_probability and probability_maximizer_objective describe the computed probability maximizer, with the lexicographically first index on exact ties. The most_probable_* fields describe the representative selected by the numerical tie rule. mode_status is resolved when that rule accepts only one positive observed valid point, unresolved when it accepts another, and unavailable when there is no positive valid outcome or no derived window. A zero deficit can still have unresolved status.

For exact readout, resolved status gives a unique computed maximum on the compared set of points. If the window bounds pairwise probability errors for the stated numerical reference, it also establishes that reference's unique maximum on that set. Omitted probabilities are covered only when the observation establishes their values and the same error bound separates them. For counts, resolved status means a unique largest pooled integer count. It supplies no confidence statement about the sampling population's mode.

Under the assumption that every pairwise probability difference has error at most the window W, a selected representative has model-probability deficit at most 2W when its accepted subtraction is exact. For the other binary64 subtractions the bound is W + W/(1-2**-53), provided the subtraction is finite, and either bound can be capped at 1 for normalized probabilities (Proposition 45). A resolved status uses strict rejection of every competitor and remains valid under the monotonicity of rounded subtraction. These statements concern probability ordering in the specified readout model, not objective optimality.

The window of theory_flavor="schrodinger"'s expm_multiply calls grows with the number of calls and steeply with each call's generator norm, which grows with the step duration, the inverse squared grid spacing and the objective's range. Evaluated from its formula (measured evidence), it is about 3e-13 for three grid points and two steps, 1e-11 for a 3-by-3 grid with 80 steps and 3e-8 for 39 grid points with two steps. The split-step window grows with the number of steps and with the angles that the state's own modes and points pick up. With a small shifted-cubic parameter a window can exceed the gap between the most probable points of two wells, as the Plan's window in the K = 16 example in Numerical guarantees does, so inspect mode_status beside most_probable_deficit and most_probable_tie_window. The summary prints all three. Exact Aer probabilities give 2.8e-11, 7.6e-11 and 1.5e-10 for 36, 98 and 190 circuit instructions (engineering constants). Numerical guarantees compares the split-step windows with high-precision evaluations.

When a preparation record lies outside the derivation of the tie window, for example with an excluded unitary instruction from initial_state_preparation="qiskit_state_preparation" or an unknown instruction count, most_probable_tie_window is None, most_probable_tie_window_unavailable gives the reason, which the summary line repeats, and the computed probabilities are compared exactly. Roundoff then decides between points that the model makes equal, with no guarantee that the selected point is the model's most probable one. Counts are pooled as integers per grid point over all chunks, compared exactly with window and deficit zero, and divided by the total returned shots only afterwards. Equal counts are equal observed frequencies, not evidence that the underlying probabilities are equal or statistically indistinguishable. Two readout cases are not yet checked: the tie window with the NWQ-Sim per-instruction constant, and the readout roundoff of a QHD circuit that has ancilla qubits beyond its register.

position_mean and position_standard_deviation give, for each variable j, the mean <x_j> = sum_i x_j(i) m_j(i) / valid_mass and the standard deviation sigma_j = sqrt(sum_i (x_j(i) - <x_j>)**2 m_j(i) / valid_mass) conditional on a valid outcome, where x_j(i) are the grid coordinates and m_j the stored unconditional marginal, row j of marginals.array, the read-only (d, K) float64 array of the FrozenArray field marginals. marginals[j] gives the same read-only row, and iteration over marginals gives the rows in the order of the variables. Liu et al., arXiv:2607.16996v1, Eq. (94), define these moments as expectations in the final state. They are computed on demand from the marginals and the Plan's grid, are not stored or computed by reports, and are None when no valid mass was observed. They describe a result and are no criterion for the agreement of two states or for optimization success. The distributions {-1, +1} with probability 1/2 each and {-2, 0, +2} with probabilities {1/8, 3/4, 1/8} have the same mean 0 and variance 1 at total-variation distance 1, and they put probability 0 and 3/4 on the point 0. On the periodic grid both moments use the coordinates of [lower, upper) and therefore depend on where the period is cut. Mass split between x_0 and x_(K-1), which are neighbors through the wrap link, gives a mean near the middle of the box.

Each marginal entry is the accurately summed unconditional weight of its grid slice. Position moments and refinement rules use those entries. A change in reduction order can change the final bits and can change a decision at a numerical tie or threshold. Candidate objective comparisons continue to use the stored support entries combined by objective_at. The observed and valid probability masses are fsum sums over their outcomes, and counts are summed as integers before one division.

Analysis preserves full-register valid and invalid bins. To decode the raw readout directly, read a probabilities or counts chunk through chunk.histogram(): index_list() gives the outcome indices, whose bit 0 is the first observed qubit for probabilities and classical bit 0 for counts, and weights gives their probabilities or counts. No valid mass yields no candidate and no most probable point. The candidate is the observed valid point with positive weight that has the least evaluated objective, and ties go to the lexicographically smallest grid-index tuple. Marginals, candidate probability and most-probable probability are unconditional, and the expected objective is conditional on observed valid mass. A positive-weight finite mean must lie between its observed minimum and maximum. An endpoint-difference evaluation handles overflow/cancellation without clipping an invalid result into the domain. Register maps, settings, the content hash of the Plan and the associations of completed contributions are checked before decoding. Attaching or loading a Result also checks the coordinates and objectives of the candidate and the most probable point against the Plan's grid and the existing support tables. Negative objectives remain legal when the objective produces them. These checks neither evolve a state nor search the grid. The most-probable selection adds a second pass over the observed points without storing anything per point.

Cost and work limits

The structured one-hot preparation uses O(dK) gates and forms no 2K vector, and the structured binary preparation, one H per qubit, uses d b gates and prepares the uniform state only. The explicit initial_state_preparation="qiskit_state_preparation" alternative constructs 2K one-hot or K binary amplitudes per register. initial_state_preparation="none" is a quantum resource-only construction.

max_work and max_bytes refuse symbolic expansion, grid-table, circuit preparation or explicit reference work before it occurs. Their defaults are 1 billion known scalar/action units and 10 GB known arrays. They do not bound process RSS or undocumented SDK workspace. The shared ExecutionLimits separately bound measurements and stored data.

A QHD work or byte refusal reports the stage, its count and the exceeded limits. Raise a limit as Using QHD describes for a complete count and for a lower bound. A raised limit clears only its own check, so other planning or execution checks can still refuse, and raising a limit beyond a lower bound, for example to twice the displayed amount, gives progress without guaranteeing that the next attempt succeeds. Symbolic expansion and support-table limit checks can report lower bounds. An explicit verification request has its own QHDVerification work and byte limits, separate from the QHD Method's limits. Its work count can also be truncated, while its displayed byte amount is already the complete declared amount.

Choose a remedy for the named stage. Fewer variables or grid points reduce array sizes. Fewer product steps reduce repeated binary-compilation work at fixed support structure. Changing an objective scale, schedule or grid changes the finite evolution as well as its cost, so recheck the resulting model and accuracy.

The one-hot potential applies one number projector per nonzero support-table entry. For a projector on s qubits the structured provider takes the cheaper of a multi-controlled phase and the phase diagonal of 2**s - 2 CX, so it never needs more CX than the parity expansion of the same projector, (s - 2) 2**s + 2 CX. A tie selects the multi-controlled phase. The diagonal is cheaper for s = 4 to 7 (14 against 20 CX at s = 4), and the multi-controlled phase is cheaper again from s = 8 (220 against 254 CX), because its count grows quadratically in s while the diagonal's grows as 2**s. On the shifted Ackley function of Wu et al., arXiv:2605.12066v1, on a 32-by-32 periodic grid of [-5, 5)², every projector acts on two qubits, one in each register, and both counts are 2,048 CX per potential application. Projectors on more qubits can save CX, for example 1,890 against 3,942 for the product (1 + x_1)(1 + x_2)(1 + x_3)(1 + x_4) on a three-point Dirichlet grid. These counts exclude kinetic evolution and preparation, a second-order step applies the potential twice, and the per-projector choice is not optimal over a joint synthesis of the whole potential.

Planning limits gives the work and byte counts behind each check.

Fault-tolerant resources

Planning a one-hot quantum circuit with the structured or no initial-state preparation also records the number of arbitrary rotations of that circuit, a ResourceLaw with metric arbitrary_rotations in the selected_logical basis next to the CX count (resources.rotation_law). The count reads the stored blocks and preparation amplitudes and expands each emitted gate by its Qiskit 2.5.2 definition, or by NWQLib's phase-diagonal multiplexor, into Clifford gates and single-axis rotations, before any cross-gate simplification or routing. A rotation whose binary64 angle equals the stored value of 0, ±pi/2, ±pi, ±3pi/2 or ±2pi is a Clifford gate, one at ±pi/4 is an exact T gate, and every other nonzero angle is an arbitrary rotation, the classification that NWQLib's resource totals use (resources).

A fused hopping link has two rotations at its stored angle. A number projector on s qubits has 1, 3 and 7 rotations for s = 1, 2 and 3, and 2**s - 1 rotations with the phase-diagonal provider at 4 ≤ s ≤ 7. From s = 8 the multi-controlled phase has 4 s - 5 angle-dependent rotations plus the fixed rotations of Qiskit's multi-controlled X gates. At s = 8 and a generic angle these are 147 arbitrary rotations and 116 exact T gates, where the 2**s - 1 formula would give 255 rotations and no T gate. The amplitude chain adds two rotations per emitted link, and for the uniform state the two of its last link are exact T gates. Rotations that rotation_threshold removed are absent from the count, and their possible effect is the pruning entry of the error sources below, not a proved zero. Qiskit's StatePreparation recipe has no rotation count, as it has no CX count.

A binary Plan records its own rotation count (circuit synthesis) in the same selected_logical basis, so the default resource estimate reports the rotation counts of both encodings. That count includes the rotation gates before angle classification and is therefore an upper bound on the arbitrary rotations of its circuit.

circuit_resources(plan, synthesis_epsilon=E) returns the per-circuit record QHDCircuitResources without building, transpiling or compiling a circuit. It holds the arbitrary rotations and exact T gates of the emitted circuit, the budgeted Clifford replacement with its T estimate (synthesis), the evolution bound (evolution) and the error sources (error_sources). \(E\) is the total operator-norm error budget of one circuit execution for replacing rotations by Clifford gates and for synthesizing the others. It has no default value, because the budget is the user's accuracy choice, and without it synthesis is None, with no Clifford replacement or T estimate. For a binary Plan the call first checks the work and bytes of classifying its rotations against the Plan's max_work and max_bytes and raises ValueError when they are exceeded, and run_resources checks each inner binary Plan the same way with the run's kept inner Results held (rotation-count checks).

With \(N>0\) arbitrary rotations, every rotation whose distance \(d_C=2\sin(|\delta|/4)\) from the nearest Clifford rotation is at most \(q=E/N\) becomes that Clifford, which costs no T gate, and adds its distance bound to the exact accumulated error \(C\). Here \(\delta\) is the rotation angle minus the nearest multiple of \(\pi/2\), and \(d_C\) is the phase-minimized operator distance. The distance is evaluated outward, in exact rational arithmetic with the binary64 neighbors of \(\pi\) and the alternating Taylor bound of the sine. The simpler bound \(|\delta|/2\) would be at most 0.65 percent looser and about twice as fast (measured evidence), and resources._distance_bound records that comparison with its timings. The enclosure of \(\delta\) widens with the multiple of \(\pi/2\), so a rotation by a huge angle is replaced only when the budget covers that width.

Let \(E\) be the positive finite binary64 synthesis budget and \(N\) the initial arbitrary-rotation count. If \(N=0\), the replacement threshold and per-rotation tolerance are None and the replacement error is zero. For \(N>0\), candidate Clifford replacements are tested against the exact rational share \(q=E/N\). If \(B\) replacements have exact accumulated distance bound \(C\) and \(M=N-B\) arbitrary rotations remain, the record stores

\[ E_C=\operatorname{up64}(C),\qquad R=\operatorname{down64}(E-E_C). \]

For \(M>0\), the per-rotation synthesis tolerance is

\[ \epsilon_{\rm rot}=\operatorname{down64}(R/M). \]

The subtraction and division are formed from exact rational values before their directed conversion. The synthesis entry of the error sources is the same stored \(R\). Consequently the stored fields satisfy

\[ E_C+M\epsilon_{\rm rot}\le E_C+R\le E. \]

With no remaining arbitrary rotations, no per-rotation tolerance is needed. If \(M>0\) but the downward-rounded tolerance is zero, circuit_resources and run_resources refuse that budget, and the Plan is unaffected. The stored per-rotation tolerance can be below the initial exact share \(q\) because both directed roundings preserve the total budget.

The replacement describes an approximated construction, while the emitted circuit keeps every rotation, so the record reports both the emitted count and the remaining rotation count \(M\) (rotations) with the replacement error \(E_C\) (replacement_error). The real logarithmic allocation and the stored rounding are distinguished in the allocation derivation.

The T estimate is \(T_{\rm exact}+3M\log_2(1/\epsilon)\), with \(\epsilon\) denoting the stored \(\epsilon_{\rm rot}\) above, the leading term of the typical Ross–Selinger count (arXiv:1403.2975v3) with intercept zero, where \(T_{\rm exact}\) counts the exact T gates. It is a ResourceLaw with metric t, basis clifford_t and interpretation estimate, whose assumptions state that the circuit's angles are treated as typical Ross–Selinger instances. QHD's angles repeat and are structured, so a large rotation count gives the estimate no statistical support, and it carries no error bar. The logarithmic model applies to per-rotation budgets below 1, and for a larger budget the record gives the reason instead of an estimate.

For a binary Plan the record classifies the computed angles of every block the same way: three phase gates at pi/2**(r+1) for each controlled phase of each QFT, exact T gates at r = 1, one Rz per Walsh string, and the nonzero multiplexor angles of a dense diagonal, recomputed with the helper's own arithmetic. This costs one synthesis of each stored block, the table transforms of planning, plus the Gray-code angles of every dense-diagonal occurrence, O(n 2**n) on n qubits, and builds no circuit.

Liu et al., arXiv:2607.16996v1, Sec. VI A, replace every rotation within an angle threshold of 0, ±pi/2 or pi by a Clifford gate for their rotation counts of each encoding, and their analysis script (experiments/res_ana/analyze_periodic_qasm_stats.py of the authors' repository) sets that threshold to 1e-12. They divide a total synthesis budget of 1e-4 evenly over the remaining rotations of each benchmark circuit (Sec. VI B). The fixed 1e-12 window suits a rotation count. Its error per replacement is about 5e-13, but its total grows with the number of replacements and is counted in no budget. NWQLib derives the replacement threshold from the user's budget instead, so every replacement is counted in that same budget, and the emitted count compares angles exactly. Their even split is the rule here too, and their numerical budget is a benchmark setting rather than a library default.

For small circuits nwqlib.backends.nwqec.compile_logical gives the exact Clifford+T inventory of the compiled circuit within its default cap of 100,000 operations (Estimate fault-tolerant resources). On (x - 1/5)**2 over [-1, 1] with three interior points, one first-order step of duration 0.17, QuadraticSchedule(gamma=0.0) and the uniform start, the circuit has 8 arbitrary rotations and 2 exact T gates. At E = 1e-4 no rotation is replaced and the estimate is 2 + 24 log2(8e4) = 392.9 T. NWQEC 0.1.2 with rz_err="total" and epsilon=1e-4 (measured evidence) compiled the same circuit to 348 T and T-inverse gates, at the synthesis stage before fusion and in the final compiled circuit alike. The estimate is 13% above that count on this one instance, a model discrepancy rather than an error bar. Section 4 of the resource estimation notebook applies these counts and the evolution bound to one step of a three-variable quadratic in the one-hot and binary encodings at K = 8, 16 and 32, up to 96 one-hot qubits, without building a circuit.

A requested budget is not an achieved error. NWQEC 0.1.2 removes Rz angles within a fixed 1e-4 of trivial values, groups angles to four significant digits and passes them to its GridSynth as decimal strings (std::to_string), all outside the requested tolerance. Compiling H followed by Rz(pi/2 + 5e-5) with rz_err="total" and epsilon=1e-8 (measured evidence) gave H followed by S, no T gate and a phase-aligned operator error of 2.5e-5, 2500 times the request. LogicalCompilation.requested_epsilon is therefore labeled as requested, and the synthesis entry of the error sources is the requested budget, with the achieved synthesis and compiler errors unavailable.

QHDCircuitResources.error_sources names every error source of the circuit against its finite model separately, each with the status bound, estimate, requested, unavailable or not_applicable, in the order of the stages they compare. The finite model is the Hamiltonian a(t) T + b(t) V on the valid one-hot subspace or the whole binary register, with the exact schedule functions, the stored support tables and constant and the grid's binary64 spacings, over the nominal steps [k dt, (k + 1) dt] that end at the exact product N dt. The first entry, kinetic_model, instead compares that finite model with the finite-difference one. It is not applicable to the finite-difference stencil, which every one-hot circuit applies, and unavailable for the spectral model of a binary circuit, whose difference from the stencil depends on the momentum content of the whole trajectory (split-step classical evaluation). No entry concerns grid discretization or the optimization gap.

evolution_bound (QHDEvolutionBound) bounds the first four stages in the operator norm, with the physical phase included: time ordering, midpoint quadrature, the coefficient residual and splitting. The splitting bound applies the commutator bounds of Childs, Su, Tran, Wiebe and Zhu (arXiv:1912.08854v3, Sec. 5.1, Propositions 15–16) to the emitted groups (Proposition 40). First order pays |alpha beta| C/2 + alpha**2 Gamma/2 per step for the potential against the links and for the links among themselves, with Gamma = (K - 2) sum_j w_j**2 on the Dirichlet chain, (K - 1) sum_j w_j**2 on the periodic cycle and w_j = 1/(2 h_j**2). Second order pays alpha**2 beta D_T/12 + alpha beta**2 D_V/24 + alpha**3 (J_E/12 + J_O/24), whose last term is the split between the odd and even link layers and is nonzero even for a constant potential. The binary product has no link term, because its exact QFT conjugation applies each variable's kinetic factor whole, so first order pays |alpha beta| C/2 and second order alpha**2 beta D_T/12 + alpha beta**2 D_V/24, K = 2 included. Here alpha and beta are the exact products of dt and the stored kinetic and potential weights of the step.

The objective constant enters the error sources as a global phase, so a large phase-sensitive bound can include error that leaves the ideal readout probabilities unchanged and does not by itself establish poor readout accuracy.

For each time step the time-ordering error budget uses the smaller of an available derivative bound and an available integral bound (Proposition 39). For \(H(t)=a(t)T+b(t)V\) with finite Hermitian \(T,V\), nonnegative \(a,b\), step width \(\Delta\) and a commutator bound \(C\ge\|[T,V]\|\), write \(A=\int a\) and \(B=\int b\) over that step. The integral budget is

\[ w_{\rm int}=\min\{2,CAB/2\}. \]

If \(a_*\) and \(b_*\) bound the coefficients and \(A_1\) and \(B_1\) bound their first derivatives on the step, the derivative budget is

\[ w_{\rm der}=\Delta^3(a_*B_1+b_*A_1)C/12. \]

The code takes the per-step minimum of the available bounds, accumulates it with outward rounding and applies the unitary-distance cap of 2. Exact step integrals are used for the shifted cubic schedule and for the quadratic schedule with \(\gamma=0\). Other schedules use the existing derivative bound. An unavailable commutator bound remains unavailable.

Under the midpoint rule the quadrature bound \(\Delta^3(A_2\mu_T+B_2\mu_V)/24\) uses second-derivative bounds \(A_2,B_2\) and raw operator norms \(\mu_T=\|T\|,\mu_V=\|V\|\). The coefficient residual compares those values with the stored exponents. It is an exact rational discrepancy under the midpoint rule and for every rational interval integral, and a first-order estimate from the integral routine's documented error for the quadratic (\(\gamma>0\)) and cubic kinetic integrals, which makes evolution conditional on it. The commutator bounds \(C\), \(D_T\) and \(D_V\) come from the table ranges and, for the finite-difference stencil of either encoding, from the largest potential change across each link, and the layer constants from the \(K\)-point graph. The spectral model takes the range bounds alone, with the kinetic spread \(\pi^2/(4h^2)\) evaluated at the binary64 number above \(\pi\).

Computing evolution_bound costs one pass over the stored tables and steps, O(sum_S |S| K**|S| + d K + N), with no state or matrix. Every value is evaluated outward from exact rationals of the binary64 inputs. The one-variable six-point comparison used four steps and total time 0.5. The Dirichlet box was [-3.5, 3.5] and the periodic box was [-2.5, 3.5), giving grid spacing 1 in both cases. Across the recorded schedules, rules, orders and boundaries, the measured bound-to-distance ratios were 1.55–3.56. These are ratios from this comparison (measured evidence), not a bound on ratios for arbitrary inputs. For two binary variables at K = 4 and four steps to time 0.25 it exceeded the distance 2.10 to 5.06 times for finite difference and 5.38 to 11.92 times for the spectral model, over both coefficient rules and both orders.

The later stages compare the split product with the quantum circuit (circuit_errors). angle_formation bounds the stored angles of the kept blocks against their intended exponents at the nominal durations dt or dt/2, exactly. Planning accepts dt/2 as a normal number, so a stored half step is exact. rotation_pruning bounds the blocks the circuit omits, pruned by rotation_threshold or omitted below the normal range, at their exact intended exponents, since a hopping block differs from the identity by at most its exponent and a projector by at most (1 - 2**-s) times it. gate_parameters compares the fused hopping gate's own parameter with twice the stored angle, which is zero because the builder doubles the stored angle's own product, and relies on the projector and chain rotations being exact power-of-two scalings of their stored angles, which planning accepts in the normal range. identity_phase bounds the circuit's global-phase bookkeeping against the exact identity phase, from the phase-error budget of the running phase total (Numerical guarantees) and the budget of the circuit's phase assignments, about 2u per radian plus 20u per assignment, with 21u more for each phase-diagonal projector, where u = 2**-53, and the identity phases of omitted projectors and constant contributions.

state_preparation is an estimate added once. Its direction term is zero for the uniform state, about 11u per register for the kinetic ground state on a Dirichlet grid, and for a Gaussian the Plan's aggregate initial-state error u * initial_state_error plus about 3u per register, capped at sqrt(2). The chain's rounding adds a first-order (K (K - 1)/2 + pi (K - 1)) u per register, and a chain that stops before a link whose parameter would fall below the normal range adds the exact error of the omitted links. The phase-diagonal provider wraps the phase of every projector on 4 to 7 qubits through NumPy's complex exponential and argument, which document no uniform error constant, so diagonal_wrap is unavailable when such a projector is emitted. Potential compilation and the rotation count compute this wrapped phase with the same 2**s phase array and NumPy operations as the circuit construction, phases=(0,...,0,-angle) followed by numpy.angle(numpy.exp(1j*phases)), so the final entry is reproduced bit for bit. Range selection accepts an exactly zero wrapped phase and requires a nonzero phase to remain normal after its power-of-two scaling. The Clifford replacement error \(E_C\) is an outward bound that describes the approximated construction, and the synthesis entry is the same stored \(R=\operatorname{down64}(E-E_C)\), the only entry not capped at 2. The one-hot circuit applies no Fourier transform, so its AQFT entry is not applicable.

For a normalized input, entries that bound adjacent stages add by telescoping to a bound on the final state's 2-norm error, and a common measurement then changes by at most that amount in total variation. An estimate makes such a sum conditional, an unavailable entry leaves it unavailable, and the record forms no total. Sampling and device noise lie outside it.

A binary circuit has its own later stages, in this order. angle_formation bounds the blocks built from their computed angles against the exact stored-weight split product, before pruning (circuit_errors.binary_angle_formation). It bounds every Walsh string, from the butterfly error gamma_n a of each coefficient of a support table, plus the coefficient itself when its normalization falls below the normal range, or an enclosure of each analytic kinetic coefficient, the exponent's rounding and the angle products, then the dense phases, with enclosures of the kinetic energies, and the QFT angles formed from binary64 pi, 2 d N_s eps_pi (b - 2 + 2**(1-b)) over the kinetic conjugations. The Walsh identity phases stay with identity_phase. The table reductions return upper bounds on the exact absolute sums of the stored binary64 entries. A correctly rounded fsum is moved upward before division, and every later positive operation is evaluated outward. The kinetic enclosures reuse one analytic enclosure per distinct trigonometric factor and propagate interval products outward. The table sums can exceed their exact values by a few units in the last place and the propagated kinetic radii by more, so the resulting budget may be larger than an exact-rational evaluation but cannot be smaller under these arithmetic assumptions (Proposition 44). On seven small binary cases (measured evidence) the bound exceeded an independent exact evaluation of these errors 2.6 to 10.5 times, and 4.5 to 8.4 times on the four with K = 4.

aqft then compares the full stored-angle QFTs with their truncation, the reconstruction's aqft_error_bound, when a cutoff truncates them and is not applicable otherwise, and rotation_pruning is the reconstruction's pruning_error_bound, half the dropped Rz angles plus the exact errors of the rotations and dense phase entries omitted below the normal range. gate_parameters is zero, because the Walsh angles and controlled phases reach their gates unchanged or halved exactly. diagonal_wrap is unavailable when a dense diagonal is emitted, for the same NumPy wrap and for the recursive means and Gray-code angles of its circuit construction. identity_phase adds the phase-error budget of the constant's running phase total, the formation budget of the Walsh identity phases and the circuit's assignments of those phases and of the physical phase (method._native_phase_allowance). state_preparation is zero, because the H layer prepares the uniform state exactly.

Up to the quantum circuit, before Clifford replacement and synthesis, the entries after kinetic_model of a circuit without a dense phase diagonal add to a total against its finite model, conditional where an entry is an estimate. A one-hot circuit with a number projector on 4 to 7 qubits and a binary circuit with a dense diagonal have no such circuit total, because their diagonal_wrap entry is unavailable (Limitations and open work). A binary circuit has a dense diagonal where dense_diagonal is chosen, or where min_cx picks it because it needs fewer CX, as for the finite-difference kinetic phase from b = 3. Beyond the quantum circuit, the achieved synthesis and compiler errors are unavailable for every circuit.

The bounds rest on platform assumptions that each entry states. These are round-to-nearest binary64 arithmetic with gradual underflow, a correctly rounded math.fsum for the binary phase and Walsh sums, and at most one ulp of error in math.sin for the AQFT bound and in math.sin, math.cos and NumPy's sin for the finite-difference enclosures of the binary angle formation, whose spectral enclosures need no trigonometric assumption. The one-ulp assumption is conditional. The bounds were tested against independent evaluations on macOS arm64 and Linux aarch64 (measured evidence), and NumPy's SIMD sin on Linux x86-64 hosts has not been checked (dependency issues).

Planning omits a contribution whose exact value would be nonzero and below the normal binary64 range, and adds its exact action to one entry of the error sources (normal binary64 range).

Each augmented-Lagrangian round and refinement level uses its own tables, spacing and exponents. ALResources.arbitrary_rotations and RefinementResources.arbitrary_rotations sum the per-circuit count over the rounds or levels that prepared their circuit, one circuit each and not multiplied by shots, as their cx fields do. A round or level that prepared no circuit, such as one that stopped before its preparation, whose Run's creation raised before the Run recorded its run-log header, or that ran classically, contributes zero. One whose preparation raised after that header, or whose preparation count is unknown, makes the sum unknown, since its record does not establish that the circuit was built. For a quantum run these fields therefore equal the body total arbitrary_rotations of run_resources.

run_resources(result, synthesis_epsilon=E) reads each round, or each level of a round's refinement, and reports two totals: the compiled circuits counted once (arbitrary_rotations, t_estimate), and every circuit weighted by the raw shots set aside for its attempts (shot_arbitrary_rotations, shot_t_estimate), an upper bound on the executed shots, each of which runs the whole circuit with its preparation. The rotation counts come from the recorded resources and the T estimates from the Plan objects of the live inner results. A binary run records its upper-bound rotation count, so its rotation_interpretation is upper_bound, while its T estimates use the exact classification of each circuit. A round or level that adds a positive recorded count without a live inner result, such as a level that stopped the run, also makes it upper_bound, because its record does not keep that count's interpretation. A resumed run records the same counts as an uninterrupted one. A value that is unknown for a round or level with a circuit, such as the T estimate of a level whose Run raised and left no inner result, makes the affected totals unavailable with the reason, while the entries keep the per-circuit values of the others. Exact readout uses no shots, while hardware would need some, so its shot totals are unavailable, with a reason that says to multiply each entry's per-circuit values by the intended shots.

Explicit verification

result.verify(checks=QHDVerification(comparisons=("grid_minimum",))) returns (receipt, facts). facts holds the checked fact of each comparison. grid_minimum reports the candidate's gap within a freshly evaluated binary64 objective grid, in the Problem's objective unit. The explicit grid reference evaluates the original objective and required constraints on the declared original-variable grid. It chooses the lexicographically first point attaining the least computed feasible objective. Feasibility uses the normalized residual test at the stored tolerance. All reported reference values come from the same evaluated table. The result describes this finite numerical grid and supplies neither a continuous optimum nor a proven expression-evaluation error bound. The table is one array evaluation of the original objective at the original coordinates, so its entries can differ from a scalar evaluation of the same expression by an ulp or more, which can change a minimum or a tolerance decision near a threshold. It is an evaluated-value diagnostic. Converting it to a bound for the mathematical objective requires an objective-evaluation error bound. The comparison can also report minimum_success_mass, the observed weight of grid points whose evaluated objective differs from the evaluated minimum by at most minimum_tolerance. With zero tolerance, this counts ties in the evaluated values. The mass is computed from a kept state (keep_state=True), exact circuit probabilities or counts, and a classical Result without a kept state reports it as unavailable.

receipt also reports most_probable_gap, the gap of the most probable point, which the augmented-Lagrangian and refinement layers read by default, from the same evaluated grid. Certificate.with_verification compares them with objective_gap_tolerance and the infidelity tolerances (Check accuracy and verify a result). Fidelity choices schrodinger_fidelity and ir_product_fidelity require an existing state and run only the named restricted model. The latter requires a product already selected in the Plan. Raw infidelity and its reduction-derived roundoff window stay visible. Only a negative value within that window can become zero. Neither a grid comparison nor conditional fidelity supplies missing discretization, leakage or global optimization guarantees.

receipt and each fact record the computation that produced the Result and the reference it is compared with. Repeating the original theory_flavor is a fresh replay against the stored state. It checks consistency and can detect changed stored data, but does not provide independent accuracy evidence. A different theory_flavor compares models on a shared grid/schedule. For a split-step Result, schrodinger_fidelity compares against the unsplit finite-difference model, so it measures the splitting error and, with kinetic_model="spectral", also the difference between the two kinetic models. Quantum-versus-one-hot-IR checks construction consistency, and quantum-versus-Schrodinger compares the circuit to the restricted evolution. Quantum-versus-binary-IR compares the readout with a numerical reference that applies selected QFT operations and computed phase arrays. Walsh phase reconstruction can differ from the exact product of the stored rotations, and its R_W bound is not propagated into this comparison's fidelity uncertainty. Passing the tolerance therefore has only the scope of this numerical-reference comparison. Every requested evolution is checked against its work and byte limits. Reporting or loading verification records performs no replay.

Compare a circuit with a reference that uses the same kinetic model before reading a difference as an implementation error. The signed spectral kinetic, the finite-difference kinetic and an AQFT cutoff define separate comparisons. On the two-dimensional Ackley function of Liu et al., arXiv:2607.16996v1, with K = 8 per axis and 512 second-order steps (measured evidence), the full-QFT binary circuit had fidelity 0.99999996 to its own finite-difference model, while the finite-difference and spectral models differed by fidelity 0.977 and aqft_cutoff=1 lowered the circuit's fidelity to 0.975. Report fidelity and total-variation distance for dynamics, and for optimization an event with a declared tolerance, such as a feasible objective within a stated distance of a reference value. A sampled best point, the most probable point, an evaluated grid minimum under full positive-support exact readout and a continuous optimum are different quantities. State the sampling budget, the grid and time resolution, the preparation, every repeated circuit and any augmented-Lagrangian or refinement work with the result. Reproducing a published table also needs the configuration the authors executed and their data, which agreement on a different model does not supply.

Constrained problems

solve_augmented_lagrangian solves a ConstrainedOptimization, the box problem with equalities h_i(x) = 0 and inequalities g_j(x) <= 0, by an augmented-Lagrangian sequence of QHD box problems. Each round selects an inner objective, obtains a point from its QHD solve or box refinement, projects that point to the original variables and updates multipliers and a penalty. Automatic inequality selection can plan additional candidates before the selected inner solve. The returned ConstrainedQHDResult holds the portable AugmentedLagrangianRecord and the inner QHD Results. The rounds read points of the finite grid of the box, except that the mode_or_mean rule can take an off-grid mean position, and no status of the run is a statement about the continuous problem. Li, Fan and Han, arXiv:2508.02969v1, also combine an augmented Lagrangian with QHD, and simulate the QHD dynamics on classical hardware with the simulated bifurcation algorithm.

import sympy as sp
from nwqlib.problems import ConstrainedOptimization
from nwqlib.algorithms.qhd import (
    AugmentedLagrangian,
    QHD,
    constrained_grid_minimum,
    solve_augmented_lagrangian,
)

x, y = sp.symbols("x y", real=True)
problem = ConstrainedOptimization(
    objective=(x - 1)**2 + (y - 1)**2,
    variables=(x, y),
    bounds=((0.0, 1.0), (0.0, 1.0)),
    inequalities=(x**2 + y**2 - 1,),
)
result = solve_augmented_lagrangian(
    problem,
    qhd=QHD(num_grid_points=4, num_steps=80, total_time=10.0),
    options=AugmentedLagrangian(stationarity=True),
    execution="classical",
    seed=7,
)
print(result.termination, result.candidate, result.objective)
print(result.best.evaluation.infeasibility)
print(result.multipliers())
print(constrained_grid_minimum(result).gap)
if result.last is not None:
    evaluation = result.last.evaluation
    print(evaluation.stationarity, evaluation.stationarity_unavailable)

feasible_complementary (0.6000000000000001, 0.8) 0.1999999999999999
2.220446049250313e-16
((), (0.5600000000000009,))
0.0
0.49600000000000166 None

The first two rounds read the infeasible point (0.8, 0.8), where g = 0.28, and raise the multiplier to 0.56 and the penalty to 2. The third round reads (0.6, 0.8) on the circle, and the run stops with status feasible_complementary at objective 0.2. Its infeasibility, max(||h||_inf, ||g_+||_inf) on normalized residuals, is the rounding residual g = 2.2e-16 of that point, which the default feasibility tolerance 1e-9 accepts. The explicit reference finds the same point as the least evaluated objective among feasible points of the stated grid, with a signed difference of 0. A stopping status does not assess optimality. With QHD(num_grid_points=4) and the other Method defaults, this example uses one step and total time 1. The first round's most probable point (0.6, 0.6) is feasible with a zero multiplier, so the run stops there with the same status, 0.12 above the least evaluated objective among feasible points of the stated grid. The optimization notebook, Section 7, runs this problem and compares the point with the continuous solution.

For this example, the comparison uses the least evaluated objective among feasible points of the stated grid. Its signed difference is an evaluated-value diagnostic. A bound for the mathematical objective requires bounds for objective-evaluation error and an appropriate feasible comparison set.

result.termination names why the run ended.

termination	The run ended because
feasible_complementary	the default feasibility-and-complementarity test, described below, passed
feasible	termination="feasibility" was set and the point passed the feasibility test alone
iteration_limit	max_iterations rounds completed
no_valid_point	an inner result had no valid grid point
budget_exhausted	the remaining cumulative limits cannot cover another round
inner_failed	an inner planning, preparation or execution raised in a later round, for example because a larger penalty makes planning exceed max_work or a Run refuses a remainder that passed the round check, or f, h or g was not finite and real at a later round's selected grid point or a point checked by its refinement

In the first round no round has completed, so the original exception propagates, and a configuration that QHD rejects surfaces at once.

The statuses no_valid_point, budget_exhausted and inner_failed update nothing, and the completed rounds stay in the result. best is the feasible round with the least objective, or the round with the least violation when none is feasible. The multipliers and the penalty belong to last, the last round that chose a point, and do not form a KKT pair with best.

Round k keeps the normalized PHR effective objective L_k for the outer update. With F=f/s_f, H_i=h_i/s_{h_i}, G_j=g_j/s_{g_j}, entering multipliers lambda_bar and mu_bar, and penalty rho, it is

\[ L_k(x)=F(x)+\sum_i\bar\lambda_iH_i(x)+\frac\rho2\sum_iH_i(x)^2 +\sum_j\frac{[\bar\mu_j+\rho G_j(x)]_+^2-\bar\mu_j^2}{2\rho}. \]

The equality part is Eq. (6) of Wu et al., arXiv:2605.12066v1. The inequality part is the Powell–Hestenes–Rockafellar (PHR) term (Rockafellar, doi:10.1007/BF01580138). Under inequality_form="phr", it is used directly and adds no slack register. A converted inequality instead adds one explicit normalized slack coordinate to the inner problem, with the cap, finite-grid error and projected outer update of Proposition 54. The PHR constant -mu_bar_j**2/(2 rho) contributes only a global phase, and L_k equals F at an exactly feasible point complementary to the entering inequality multipliers. Term by term, L_k is Eq. (10.3) of Birgin and Martinez, doi:10.1137/1.9781611973365 (Sec. 10.1, p. 114). Their Eq. (4.3) completes the squares and exceeds L_k by the sums of lambda_bar_i**2/(2 rho) and mu_bar_j**2/(2 rho), which change neither its minimizer nor its multiplier update.

After the round, the tentative multipliers are lambda+ = lambda + rho h(x) and mu+ = max(0, mu + rho g(x)) (Birgin and Martinez, Practical Augmented Lagrangian Methods for Constrained Optimization, doi:10.1137/1.9781611973365, Eqs. (4.7)–(4.8), and Wu et al. Eq. (8)). With multiplier_bounds the next round uses them clipped to the bounds (the book's Algorithm 4.1, Step 4). Without bounds it uses them unchanged, and the result states that no safeguard was used. Algorithm 4.1 also starts the multipliers inside the bounds, so when the equality interval of multiplier_bounds excludes zero, a problem with equalities needs explicit equality_multipliers, and the default None, which starts every equality multiplier at 0, is refused. The penalty stays after the first round and whenever max(||h||_inf, ||V||_inf), with V_j = min(-g_j, mu_j/rho), falls to at most reduction_ratio times its previous value, and otherwise grows by penalty_growth up to max_penalty (Algorithm 4.1, Step 3, Eq. (4.9)). Wu et al. keep rho fixed or increase it by the update rules of Nocedal and Wright, Numerical Optimization, doi:10.1007/978-0-387-40065-5, and their code doubles it after every round, which penalty_update="every_iteration" reproduces. With the same growth factor, the test of Eq. (4.9) never gives a round a larger penalty than growth after every round does.

Every test acts on the normalized quantities f/s_f, h_i/s_{h_i} and g_j/s_{g_j} with the scales objective_scale, equality_scales and inequality_scales. Their default 1 assumes a dimensionless problem. Writing a constraint as c g <= 0 with scale c s_g, c > 0, leaves every normalized quantity, and so every decision of the run, unchanged in exact arithmetic. In binary64 the tables and residuals of the two runs agree only up to rounding, even when c is a power of two, and rounding can change a decision. The layer builds L_k symbolically and divides each term whose scale is not 1 by that scale as a binary64 number before QHD planning expands L_k, so SymPy rounds the term's coefficients to 53 bits and rounds their products, while a term of scale 1 keeps exact numbers such as 1/5. The tables of L_k can then differ in their last bits and change a point chosen among nearly tied values. For another c the normalized residuals also differ by the rounding of c g and of the division by c s_g. In either case a test that compares a value exactly at its threshold can go either way.

An unscaled factor of 1000 instead multiplies the quadratic penalty by 10**6. The record keeps original-unit residuals next to the normalized ones, and result.multipliers() converts multipliers to original units with lambda_i = (s_f/s_{h_i}) lambda~_i and mu_j = (s_f/s_{g_j}) mu~_j. It forms each product exactly and rounds it once, and raises for a value outside the normal binary64 range, which the summary then reports instead of a number. Every decision uses the normalized multipliers, so an initial multiplier whose original-unit value lies outside that range does not stop the run, and result.multipliers(tentative=False) raises for it in the same way when round 0 is the last round.

inner_point selects the point each round reads. most_probable, the default, is the most probable valid grid point, where QHD concentrated probability, read from exact probabilities or from observed counts alike with the tie rule of Readout. best_observed is the QHD candidate. Within a level, best_observed chooses the least observed table value of the objective that level solved. A search-scaled level solves its normalized search objective, while a physical level solves the round's inner objective, L_k under the PHR policy and L(x, s) with an inequality representation. The added-constant example of box refinement shows how these values can tie.

mode_or_mean compares the most probable joint point with the full mean position conditional on a valid outcome, using the round's inner objective, L_k under the PHR policy and L(x, s) with an inequality representation. In an unrefined round the mean's value is computed from the original f, h and g at its projected x and, when present, its mean slack coordinates. The mean replaces the grid point only when its value is strictly smaller, so an exact tie keeps the grid point. An invalid mean evaluation keeps the grid point and records why. Wu et al.'s code instead takes the mean unless the grid point's value is strictly smaller and evaluates both points with one function. The unrefined layer compares a grid table value with a separately evaluated mean value, so roundoff can affect their order (_outer.mean_point). Refinement compares the two recorded relative values of its inner objective. On a periodic axis the coordinate mean depends on the box's cut, as the position moments of Readout do. No point rule proves inner minimization.

A refined round compares the recorded relative values of its inner objective at the completed levels and chooses the least one, with the earlier level winning an exact tie. The inner objective is L_k under the PHR policy and L(x, s) with an inequality representation. Under every point rule the refinement boxes depend on the distribution. most_probable is the default because it keeps a point that the computed distribution defines. With the default kinetic initial state and theory_flavor="schrodinger" (measured evidence) it matched best_observed on exact refined constrained Rastrigin and on the standalone double-well and Ackley comparisons, and on sampled refined constrained Rastrigin it tied twice and gave a better point once in three seeds. With theory_flavor="split_step" and 3,200 steps (measured evidence) it again matched best_observed on exact refined constrained Rastrigin, and in three sampled seeds each rule gave the better point once and they tied once. With the uniform initial state (measured evidence) best_observed gave better points on refined constrained Rastrigin and on sampled standalone Ackley, so it remains an explicit option. The point rule, the initial state and the refinement defaults are the same for theory_flavor="schrodinger" and theory_flavor="split_step". Feasibility and complementarity describe the selected point and its multipliers, not the accuracy of the evolution.

The default stopping test, termination="feasibility_and_complementarity", ends the run with status feasible_complementary when max(||h||_inf, max_j |min(-g_j, mu+_j)|) <= epsilon_c and max(||h||_inf, ||g_+||_inf) <= epsilon_f, the normalized Eqs. (10.7)–(10.8) of the book's Algencan. The projected stationarity of its Eq. (10.6) is not a stopping condition, so the status is not convergence in Algencan's sense and does not assess optimality. Even with exact inner minimization it can stop away from the feasible grid minimum. On the grid {0, 1} with f = (x - 1)**2/10, g = x - 1 <= 0, mu = 1/2 and rho = 1, L(0) = -1/40 < L(1) = 0, so the run stops at x = 0 with mu+ = 0, while the feasible grid minimum x = 1 has an objective 0.1 lower (Proposition 47).

The book's penalty measure of Eq. (4.9) is not offered as a stop. With the default kinetic initial state and the most probable or best observed point (measured evidence), continuing past the default stop to a stop on that measure found the same best feasible points on five exact refined problems and six sampled runs. Where the two stops differed, continuing took more rounds and a larger final penalty, for example 11 rounds and rho 512 against 10 rounds and rho 256 on constrained Rastrigin. With the uniform initial state or mode_or_mean, some runs did find better points by continuing. Neither test proves stationarity or optimality. termination="feasibility" stops on feasibility alone with status feasible, and penalty_update="every_iteration" grows the penalty after every round. These are the two rules of Wu et al.'s code, which stops when the violation is strictly below the tolerance, while NWQLib stops when it is at most the tolerance.

Before the first round, the layer checks the expanded support tables of f, h and g on the preprocessed grid. When expansion leaves a supplied restriction without unconditional structural coverage, it also scans the term's original additive summands on their own support grids and checks a rounding-aware bound on their additions. An inconclusive represented range uses the same check, with a limit-checked whole-term scan if the addition bound is inconclusive. This is a numerical finiteness check. Terms checked only through their expanded tables can still fail in their original arithmetic.

The original expressions are evaluated at the selected round point and any attempted mean. A failed mean falls back to the grid point and records the reason. Refinement also checks every point compared at each level. Finite rounded values do not prove the exact real domain, equality of the original and expanded numerical evaluations, or finiteness throughout the continuous box. A mathematically singular expression can have a finite rounded value. An inactive Piecewise branch may produce a warning or a nonfinite intermediate during array evaluation without invalidating a finite real selected table value. A nonfinite or nonreal selected array value is refused, and an exception in a required original-expression scan or point evaluation can still cause refusal.

Request the diagnostic with options=AugmentedLagrangian(stationarity=True). Read it from result.last.evaluation.stationarity when result.last is present. If the value is unavailable, result.last.evaluation.stationarity_unavailable gives the reason. Stationarity is optional and does not enter the default feasibility-and-complementarity stopping rule.

The summary labels the selected result as “Best feasible point” or “Least-violation point, no feasible round”, as appropriate. Read the best point separately from the last round's evaluation and tentative multiplier update. Unknown resource counts are reported with their reason.

The explicit grid reference evaluates the original objective and required constraints on the declared original-variable grid. It chooses the lexicographically first point attaining the least computed feasible objective. Feasibility uses the normalized residual test at the stored tolerance. All reported reference values come from the same evaluated table. The result describes this finite numerical grid and supplies neither a continuous optimum nor a proven expression-evaluation error bound. constrained_grid_minimum(result) checks its work and bytes against its limits, evaluates f, h and g on all K**d grid points in C-order slabs sized from its max_bytes, and keeps no table of all points. Its signed difference subtracts this freshly evaluated minimum from the result's stored best objective. That difference can be negative, including when the stored point belongs to the evaluated feasible grid, because the two values can have different rounding errors. An off-grid or infeasible stored point gives no grid-optimality conclusion. If no grid point passes computed feasibility, indices, point, objective and gap are None. The gap is also None when the result has no best point. The helper does not support a result that uses box refinement. The solver never calls it. See the evaluated-gap derivation and the stationarity equations.

A grid-point rule can satisfy an equality tolerance only if its grid contains a point whose computed normalized equality residual passes that tolerance. Increasing the penalty or updating multipliers does not add grid points. Failure of this condition prevents feasibility at that resolution, but does not imply that selected points oscillate or that the penalty grows without bound. Safeguards, stopping tests, work limits or an inner error can end the run. A problem with equalities must set feasibility_tolerance explicitly.

With inequalities only, the default is 1e-9 on normalized residuals. It classifies the residual of a returned point and does not claim nine-digit accuracy of its position. A grid point strictly inside the feasible set has zero violation however far the boundary lies from the grid, and a grid point on the boundary can have a rounding residual, such as g = 2.2e-16 at (0.6, 0.8) in the example above, which a positive tolerance accepts. For x² + y² - 1 <= 0 on the interior grids of [-1, 1]² with K = 16, 32 and 64, the smallest positive residuals, 1/289, 1/1089 and 9/4225, lie far above both 1e-9 and 1e-6, so both values classify these grids alike, and with the default point rule the tested refined runs on the disk and on constrained Rastrigin (measured evidence) stopped at the same points and rounds with either value under theory_flavor="schrodinger" and theory_flavor="split_step".

A finer grid or a changed box can make a passing point available. A larger feasibility tolerance can permit a nonzero residual, which changes the acceptance criterion. A mean point rule can return an off-grid point, so the fixed-grid obstruction does not apply to it in the same way. The choice between such values matters where a point can have a small positive violation well above rounding, as the off-grid mean of mode_or_mean can. One such refined disk run (measured evidence) stopped after 11 rounds with violation 2.9e-7 at tolerance 1e-6, and reached the limit of 15 rounds with violation 1.04e-9 at 1e-9, so a caller who accepts a larger violation sets feasibility_tolerance explicitly.

For a feasible inequality, let s_j=-g_j(x)/s_{g,j}>=0 be its normalized slack and mu_j^+>=0 its tentative normalized multiplier. This complementarity slack is computed from the original residual at the projected point and is independent of the inner slack coordinate that a converted round sampled. The complementarity test compares min(s_j,mu_j^+) with the tolerance. A strictly positive slack alone does not force a small multiplier. The test forces mu_j^+<=epsilon_comp only when s_j>epsilon_comp. When the slack itself is within tolerance, a larger multiplier can still pass that component. Finite-grid multipliers are dual iterates of the grid problem, not the continuous KKT multipliers.

The reported refusal thresholds belong to the stated example and configuration. Under max_work = 100_000_000 and with zero multipliers (measured evidence), classical planning of the example above (4 grid points, 80 steps, total time 10) passes up to rho = 214 and refuses rho = 215, and with one step and total time 1 it passes up to rho = 2**21. With the settings of the QHD optimization notebook (6 grid points, 80 steps, total time 10) on its double well restricted to the disk x² + y² <= 1/2, it passes up to rho = 512 and refuses rho = 1024, where the largest projector angle of quantum planning is 923 rad. A refused round ends the run with inner_failed.

absorb_bounds=True turns a one-variable linear inequality a x_j + b <= 0 into the box bound x_j <= -b/a for a > 0 or x_j >= -b/a for a < 0. On the Dirichlet grids the absorbed boundary becomes a Dirichlet boundary that the interior grid never reaches, so only include_boundary_points=True can represent an optimum on it. With boundary="periodic" the tightened box becomes one period. Its lower edge is the grid point x_0 and its upper edge is excluded, so an optimum on an absorbed lower bound is representable and one on an absorbed upper bound is not. The wrap link joins x_0 to x_(K-1) = upper - h across the absorbed bound. The layer compares the exact cuts with the box before rounding anything. An empty intersection is infeasible and raises, and a single point fixes the variable, which the grid cannot represent, so it raises with the remedy of substituting the value. A cut whose exact value is not a binary64 number, such as x <= 1/3, stays a constraint of the rounds and leaves the box unchanged, because a rounded bound would either admit infeasible points or remove feasible ones.

The summary and report() name each absorbed inequality with its bound, the box the rounds searched and whether the grid reaches an absorbed bound, and they list any initial multiplier given for an absorbed inequality, which the rounds do not use. They list the multiplier estimates by the problem position of their constraint, for example inequalities {1: 0.0} when inequality 0 was absorbed, while result.multipliers() returns them in the order of record.preprocessing.inequalities.

Each round records its measurement work (circuit preparations and attempts, shots, data bytes) apart from its classical work (table evaluations, classical evolution, construction, synthesis and the layer's own evaluations), and a total is unknown when any round's count is unknown. table_evaluations includes automatic-selection candidates that were discarded. An unrefined round counts its reused Plan once. A refined round adds all selection candidates to its level planning counts and each search-model level's unscaled table stage. If an attempt fails after table acceptance and its partial count is unknown, the total is unavailable with a reason.

A round or level whose preparation raised records the counts that its Run recorded before the error, attempts of every status included, read from the run log of its closed Run folder in a run with a directory, and its CX bound and rotation count are unknown. When the Run's creation raised before it committed its run-log header, nothing was prepared or measured, so the round or level records zero for every count, its CX bound and rotation count included, except its data bytes. These are the stored size of the files in its Run folder apart from the run log (run.sqlite, its lock and its SQLite rollback journal), or unknown with the error when the folder cannot be read. When the run log itself cannot be read, the folder may hold a committed header and counted work, so every count that the round or level reads from its Run (circuit preparations and attempts, shots, data bytes, classical evolution, construction and synthesis) is unknown, with the read error as the reason, and so are its CX bound and rotation count, while its table evaluations and the layer's own evaluations stay known. Without a run directory the counts of a round or level whose preparation raised are unknown.

limits caps the whole run, and each round's Run receives the remaining circuits, shots, data bytes and synthesis work. result.save(path) writes constrained.json, problem.pickle and each inner Result under iterations/<k>/result/, and load_augmented_lagrangian(path) reads them back without planning, evaluating or measuring. Planning limits gives the sums of limits over a whole run.

solve_augmented_lagrangian(..., refinement=BoxRefinement(...)) runs box refinement on each round's inner problem. Its initial box is the preprocessed x box followed by the round's slack boxes. It chooses the completed level with the least recorded relative inner-objective value, with earlier levels winning exact ties, and projects that level's selected point to x. Within a level, best_observed minimizes the observed table of the objective that level solved, including its search normalization when selected. The original f, h and g are evaluated at the reported projected point, and the outer multiplier and penalty update is unchanged. Refinement thus runs inside each multiplier round, the order of the reproduction scripts of Wu et al., arXiv:2605.12066v1. Without refinement each round is one QHD solve. The best and last points are chosen among rounds by the rules above, not among levels.

The refinement's point_rule reads every level and must equal inner_point, and a run with two different rules is refused before any work.

With exact evaluations and full observation of the first product grid, best_observed gives an inner-objective value no larger than that grid's inner minimum. Under the cap, mesh and arithmetic assumptions of Proposition 54, its projection has PHR value at most the first original grid's PHR minimum plus the initial slack-grid error bound. With table error at most epsilon_T in inner-objective units and relative-comparison error at most epsilon_R, the additional budget is 2 epsilon_T + 2 epsilon_R. The comparison term is unnecessary with one level. These numerical error bounds are assumptions that the API does not prove. Partial observation, another point rule, or a numerically flat level does not supply the full-grid assumption. See also selection across refinement levels.

For the equality problem f = (x - 11/16)**2 + (y - 7/16)**2, 4 (x + y - 1/8) = 0 on [-1, 1]**2 with scales 2 and 4, Proposition 51 bounds the distance of each round's point from the KKT point (3/16, -1/16) through the multiplier and penalty entering the round, under its assumption that the round's first level selects a minimum of the analytic quadratic L_k over its grid and that the round selects its best level by L_k. The library's tests confirm that assumption by checking the analytic L_k minima at the reported level points. In a classical run (K = 7, 20 steps over T = 6, the uniform initial state, the search model with gain 1, three levels and threshold 0.5, seed 11) the stopping round has the bound 0.0978, and its best level returns (1/4, 0) at distance 0.0884. The last level and the level with the least f of that round both return (29/128, 5/128) at distance 0.109, which also passes the stopping test and would have ended the run there. The bound and these distances are properties of that example. General refined rounds select by recorded relative values of their inner objective.

limits stays cumulative over all rounds and levels, and each level's Run receives what the earlier rounds and levels left. Level z of round k plans from the z-th child of the round's random stream, with spawn key (k, z - 1). A refinement that completed a level always gives the round its point and update. When it stops with budget_exhausted, inner_failed or no_valid_point, the run stops after this round, because the remainder that refused a level also refuses the next round's first level, and an inner failure or a result without a valid point ends a run without refinement too. The run then ends with feasible_complementary or feasible when the round's point meets the stopping test, and with the refinement's reason otherwise, which the round's refinement record keeps in either case.

A refinement that completed no level leaves the round without a point or an update, and the run ends with the refinement's reason. Besides no_valid_point and inner_failed, this can be flat_objective or unresolved_objective, when the inner objective's tables on the round's initial box show no variation or none the refinement can resolve, or width_floor, when an initial side is at the refinement's width floor. Each round checks its initial inner box, consisting of the original preprocessed x box and that round's slack boxes. An error in the first level of round 0 propagates. Any later error of a level ends the run after keeping the completed rounds and levels, with inner_failed, or with the success status when its round kept a point that meets the stopping test. No round or level is retried. A stall split, when the refinement options ask for one, acts within each round's refinement with the options' budget, and a round whose refinement stopped with split_limit goes on like one that stopped with box_unchanged.

A run with refinement executes at most R L inner-level Plan objects, with R = max_iterations and L = max_levels, and automatic inequality selection can plan additional candidates. A level that stops or fails counts among the max_levels attempts. Planning limits gives the sums of limits over such a run.

Each round's resources combine the refinement's level totals, the representation-selection table evaluations and the layer's checks of the original functions. table_evaluations includes each search-model level's unscaled table stage and solved Plan, and every automatic-selection candidate. A total is unknown when a contributing count is unknown. Each round's record nests its refinement record (ALIteration.refinement), the summary lines and report() give the number of levels and the stopping reason of every round, and save writes each level's Result under iterations/<k>/levels/<z>/result/. constrained_grid_minimum refuses a run with refinement, whose levels search nested grids that its enumeration of the preprocessed box's grid does not cover.

Slack variables for inequality constraints

AugmentedLagrangian(inequality_form="phr") is the default. It keeps the Powell–Hestenes–Rockafellar (PHR) term for every kept inequality. inequality_form="slack" requests explicit slack coordinates after the eligible quadratic and constant branch tests. inequality_form="auto" selects a representation by the execution rule below. Each converted inequality appends one normalized slack coordinate to the inner problem, after the original variables and in constraint order.

An unsimplified PHR term involving r variables needs a support table with K**r entries. If the normalized inequality is additively separable after grouping by variable, expanding its slack term produces supports of at most two variables. For a general inequality, products have supports contained in unions of the original supports, and a slack cross term can involve all r variables plus the slack. QHD checks the expanded supports against its limits before evaluating their tables. The cap construction and the slack-elimination identity are in Proposition 54.

Planning through plan_augmented_lagrangian (measured evidence) gave 131,148 CX per circuit for PHR and 666 for forced slack with eight original variables. With fourteen original variables, PHR was refused before its table was formed, with at least 4,026,532,346 units of planning work counted against max_work=100_000_000, while forced slack required 1,830 CX per circuit. Automatic selection chose slack in both cases. Each was a round-0 quantum Plan for f = sum_i (x_i - 1/2)**2, g = sum_i x_i - 3 <= 0 on [0, 1]**n, binary periodic K=4, one step, total time 0.001, seed 7, unit scales, zero entering multipliers and penalty 1, with the other QHD and layer defaults. The slack cap was 3 and its periodic box was [0, 4]. The check used Python 3.12.14, Qiskit 2.5.2 and Aer 0.17.2. These are planning and resource-count results, not executed trajectories or runtime measurements.

The slack changes the inner problem. The outer state remains that of the PHR algorithm. For G_j = g_j/s_{g,j}, entering multiplier mu_bar_j >= 0 and penalty rho > 0, define the PHR term P_j. The continuous slack minimizer and its value satisfy

\[ P_j(t)=\frac{[\bar\mu_j+\rho t]_+^2-\bar\mu_j^2}{2\rho}, \qquad s_j^*(x)=[-G_j(x)-\bar\mu_j/\rho]_+, \qquad \min_{s_j\ge0}\left\{\bar\mu_j(G_j(x)+s_j)+\frac\rho2(G_j(x)+s_j)^2\right\}=P_j(G_j(x)). \]

The identity also holds with an upper cap containing this minimizer. The round selects a joint inner point, projects it to x, and evaluates the original f, h and kept g there. Its tentative multipliers are lambda+ = lambda_bar + rho H(x) and mu+ = max(0, mu_bar + rho G(x)). The safeguard, penalty rule, stopping tests and stationarity diagnostic use those projected original residuals. The sampled slack and G_j(x) + s_j do not replace them. The slack coordinate has the units of normalized G_j, so the corresponding original-unit expression is g_j(x) + s_{g,j} s_j.

Proposition 54 proves the partial-minimization identity and distinguishes it from the dynamics of a joint QHD evolution. A finite slack grid, its additional kinetic operator and the selected point can change the projected distribution and multiplier trajectory. They do not establish the inner-minimization assumptions of the convergence or objective-gap statements in Result 46 and Proposition 47.

The preprocessing tables supply a lower bound ell_j and upper bound M_j for each kept normalized inequality. If G_j = c_0 + sum_A T_A(x_A) is an exact support decomposition on the stated grid, the bounds are

\[ \ell_j=c_0+\sum_A\min T_A, \qquad M_j=c_0+\sum_A\max T_A, \qquad U_{0,j}=[-\ell_j]_+. \]

The cap U_0 contains every required continuous slack on that grid for every nonnegative entering multiplier and positive penalty. Grouping all terms of each variable before taking extrema gives the exact grid extrema for an additively separable inequality in exact arithmetic. Overlapping supports can make the bounds loose. The implementation sums extrema of the stored unnormalized tables and divides by the positive scale exactly as rational numbers, then rounds ell_j downward and M_j upward. These bounds apply to the tabulated representation. They do not include errors from table evaluation or symbolic normalization relative to the mathematical original function. cap_scope="grid" means the grid of the original preprocessed box. It proves nothing for the continuous box, an off-grid mean or a later refinement grid. The cap conditions and the arithmetic qualification are in Proposition 54.

Each slack uses the Method's common num_grid_points, boundary and include_boundary_points. There is no separate slack-grid setting. A Dirichlet slack box has upper end U = U_0. A periodic slack box uses the upward-rounded margin U >= K/(K-1) U_0, so the ideal last node U(1-1/K) reaches U_0. The margin improves endpoint coverage. The periodic kinetic still couples the first and last slack nodes, so it defines a periodic search model, with no hard wall at zero.

Before a conversion, eligible rounds test ell_j + mu_bar_j/rho >= 0 and use the direct quadratic term mu_bar_j G_j + (rho/2) G_j**2 when it holds. Otherwise, M_j + mu_bar_j/rho <= 0 permits the constant term -mu_bar_j**2/(2 rho). The comparisons use exact rational arithmetic on the recorded binary64 values. These tests run only with inequality_form="slack" or eligible quantum auto, without refinement and with a grid-point rule. They are skipped under phr, under classical auto, with a user-supplied GaussianState under auto, with refinement, and with mode_or_mean.

A zero cap creates no slack axis. It gives the quadratic branch when that branch is eligible and otherwise leaves the PHR form. A zero cap means the bounded representation is nonnegative on the stated grid, not that the inequality is inactive. The inequality stays in the outer residual and multiplier lists. Forced slack refuses a positive-cap conversion when its range, slack grid or recorded error bound cannot be represented within the implementation's numerical requirements. Automatic selection records the unavailable conversion and keeps the current form.

For each converted inequality, SlackAxis.error_bound stores an upward-rounded evaluation of the slack-grid excess formula, using this round's entering multiplier and penalty. Its interpretation as a bound on the excess from finite-grid rather than continuous slack minimization requires the cap and mesh assumptions below. With h = SlackAxis.spacing and the recorded upper bound M_j, the formulas are

Grid convention	Ideal spacing	Recorded per-axis formula
dirichlet_interior	h = U/(K+1)	h max(0, mu_bar_j + rho M_j) + rho h**2/2
dirichlet_endpoints	h = U/(K-1)	rho h**2/8
periodic, with the cap margin	h = U/K	rho h**2/8

The implementation forms each formula exactly from the stored binary64 spacing, multiplier, penalty and range bound and rounds its result upward. InnerRepresentation.error_bound sums the stored per-axis bounds and rounds upward again. It is zero when no slack axis was added. No accuracy tolerance is applied to this number. A per-axis bound outside the finite binary64 range makes that conversion unavailable, and an unrepresentable aggregate bound raises an error under either slack or auto.

For a valid cap, sufficient mesh assumptions are that the smallest interior node is at most h and every required continuous slack is within h of a node, or, for the endpoint and periodic formulas, that zero is present and every required slack is within h/2 of a node. The ideal grids with the exact spacing relations in the table satisfy those conditions, with the margin on the periodic grid. The implementation uses rounded spacings and coordinates and supplies no separate budget for their effect on endpoint coverage or distances. Rounding the final formula upward therefore does not by itself bound the excess on the actual binary64 grid. Table-evaluation and normalization errors are also outside this number. It proves no bound on the quality of a QHD-selected point, a probability-distribution error or the continuous constrained optimum. The initial-grid formula also does not automatically apply after refinement changes the original or slack grid. Proposition 54 states the exact excess and the additional assumptions needed for an inner-minimizer or objective-gap statement.

All inner point rules act in the augmented coordinates. most_probable selects a joint grid mode. Its projection need not be the mode of the marginal distribution on x. best_observed selects the least observed table value of the objective that the level solved. A physical level solves L(x, s), while a search-model level solves its positive normalization. mode_or_mean compares the joint mode and the full conditional mean (E[x], E[s]) using the inner objective, with an exact tie going to the grid point. An unavailable mean falls back to the grid point and records mean_unavailable.

The outer ALEvaluation.point and indices contain only the original coordinates and indices. slack_point and slack_indices record the selected slack coordinates and indices when present. For any round with an InnerRepresentation, including one with no added slack, effective_value is the PHR value L_k evaluated from the original f, h and kept g at the projected point, and effective_value_source is evaluated. The selected inner value and its table or evaluated source are stored separately as inner_value and inner_value_source. Under the phr policy there is no representation record, those inner fields are None, and the existing table-versus-evaluated meaning of effective_value applies.

The recorded probability, tie deficit and tie window belong to the selected joint inner grid point. They are not marginal probabilities or tie diagnostics after summing over the slack coordinates. A selected mean has no grid indices or point probability. In exact arithmetic, the slack objective is at least the PHR objective at its projection for nonnegative slacks, under the stated branch assumptions. The two recorded values can have different evaluation and coordinate-transformation errors, so that mathematical inequality proves no comparison of the stored numbers. See Proposition 54.

Automatic conversion is enabled only for execution="quantum", including Plan objects that will run on Aer and quantum planning without execution. It makes one deterministic pass through kept inequalities in original order. Each trial adds one slack to the representation accepted so far. When neither Plan is refused, it accepts the trial only if planned classical work and CX per circuit are both no larger and at least one is strictly smaller. Equal costs or a tradeoff between the two costs keep the current form. An accepted trial with a known CX count can replace a current representation whose planning was refused. A trial without a recorded CX count is not accepted. An oversized PHR table can therefore be refused structurally before allocation while a slack Plan within the limits is selected.

The classical-work comparison is the work returned by QHD's symbolic/table limit check plus the circuit construction work. It includes that candidate's monomial, initial-state, step-row, table and compilation work. It is a planning metric, not a runtime measurement or a sum of the work of all trials. Shared layer preprocessing has its own limit check. FormTrial records the original inequality position, current_work, current_cx, trial_work, trial_cx, accepted and a reason. Unavailable counts are None. After an acceptance the next trial compares against the new current representation, so the pass does not claim to find the cheapest subset of conversions.

Under execution="classical", automatic selection keeps every inequality in PHR form for both theory_flavor="schrodinger" and theory_flavor="split_step". Their restricted state has K**d amplitudes and becomes K**(d+m_s) with m_s added slack axes. This is a policy against increasing the state dimension, not a theorem that slack always increases every planning cost. Forced slack remains available on either route. Quantum planning uses this circuit-cost rule even when the eventual backend is a statevector simulator, so a lower CX count alone does not promise a faster Aer run. The resolved representation is fixed for the round and all its refinement levels, then reconsidered at the next round. The change of inner model is subject to the scope in Proposition 54.

With m_s slack variables, box refinement operates on all d + m_s inner coordinates. The inner variable order, objective and representation stay fixed for the round, while later levels can shrink both original and slack boxes. The round projects its selected level point only after the inner comparison. Every new round starts from the original preprocessed x box and that round's slack boxes. An initial-grid cap need not cover a new x grid, and a shrunken slack interval need not contain the conditional minimizer. The recorded slack-grid error therefore concerns the initial preprocessed grid under its stated mesh and arithmetic assumptions. Proposition 54 gives the conditional comparison with that grid.

UniformState and KineticGroundState define factors for every axis and extend to the augmented problem. A user-supplied GaussianState has no configured extension to slack axes. Automatic selection keeps PHR with that state. Forced slack is refused during setup whenever any inequality remains after absorption, even if a later branch test could have avoided an axis. A best_point_gaussian constructed by refinement uses the augmented best point and remains subject to the refinement and preparation options' existing requirements.

Read round.representation.forms together with record.preprocessing.inequalities to map forms to original constraint positions. representation.slacks gives each converted position, variable name, the cap before any periodic margin, the box end upper, spacing and error_bound. The representation also stores the inner variable order, grid convention and size, cap_scope, aggregate error, trial comparisons, selection reason and the content hash of the inner objective. round.effective_objective continues to identify L_k. With a representation, round.inner_range_bound describes the inner Plan's table-range sum and effective_range_bound is None. With refinement, those round-level range fields are None and the levels hold their own range information. ConstraintPreprocessing.bounds and constraint lists remain in original coordinates.

report() exposes inequality_forms by original constraint position, the aggregate slack-grid bound and the selection reason. The complete representation is in the report's record. The summary names each round's forms and shows the bound when slacks are present. Under the phr policy the representation and the report's inequality_forms are None.

Saving a constrained result writes constrained.json, problem.pickle and the inner Results, and a durable constrained run keeps its outer record in controller.json. Method record formats names the format of each file. The loaders require their supported format. They do not convert older formats. The loaders and resume read the SymPy objective, variables and constraints from problem.pickle, and each inner Plan's objective from its saved Result or Run, the same way that load_result reads the objective of a QHD Result. They rebuild each expression node from its saved class and arguments, first with SymPy's automatic evaluation and, when evaluation changes the node, again with evaluation disabled, and they keep that second build when it reproduces the saved node. An objective built with evaluate=False as x + (x + 11/16)/(x + 11/16) therefore reopens with its quotient, which is undefined at x = -11/16 and which evaluation would replace by 1, and the checks of the reopened problem and Plan objects against their records compare the expressions as supplied. The outer record of a run directory holds a round's representation from before its first inner Run. Resume uses that saved representation and checks the reconstructed augmented objective and saved Plan/Result associations. It does not rerun automatic selection for a partially completed round. Each refined Result remains associated with its saved level Plan, including a search-model level's transformed objective. Use load_augmented_lagrangian for a result archive and resume_augmented_lagrangian for a run directory.

plan_augmented_lagrangian(problem, *, qhd, options, execution, shots, seed) returns (preprocessing, representation, plan) for an unrefined round 0. The representation is None under phr. The call performs preprocessing, representation selection and QHD planning, with no Run, preparation, measurement, evolution or grid-reference solve. nwqlib.estimate(plan) and circuit_resources(plan) describe the selected inner Plan, including its slack axes. Resource-only planning with initial_state_preparation="none" is supported. The call has no refinement argument and does not predict a refined first level, later adaptive representations, an executed trajectory or exact future costs. The whole-run sums of limits belong to the solve record and are not fields of this three-item tuple.

constrained_grid_minimum continues to enumerate only the original K**d grid and applies the original feasibility test. It imposes no slack-grid equality. It is an explicit reference, is never called by planning or solving, and refuses a result with box refinement.

A comparison of the three forms (measured evidence) exercised quantum auto on three cases of f = sum_i (x_i - 1/2)**2 with the affine inequality g = sum_i x_i - n/4 <= 0 on [0, 1]**n at n = 3 and K = 4, using KineticGroundState, joint-mode readout (most_probable), no refinement and exact noiseless Aer readout. The three cases used the one-hot Dirichlet interior grid, the one-hot Dirichlet endpoint grid and the binary periodic grid, with finite-difference kinetic, QuadraticSchedule(gamma=0.3), midpoint coefficients, total time 6 and 64 steps. All three quantum auto runs matched PHR's best feasible objective and the independently enumerated original-grid minimum under the outer settings that its measured-evidence row states. Separate round-0 planning cases reached fourteen original variables, where slack was accepted and PHR was refused under that row's planning limits, without executing the large quantum trajectories. The binary periodic planning values with eight and fourteen variables agree with the planning check earlier in this section.

For the one-hot interior case, total circuit CX decreased from 301,986 under PHR to 121,533 under auto, while observed Aer elapsed time increased from 7.986 s to 21.122 s. These sequential, unrandomized runs used Python 3.12.14, Qiskit 2.5.2 and Aer 0.17.2 on macOS 26.7 arm64, and their elapsed times included evidence extraction and the grid-reference call. They provide runtime observations, not a controlled speed benchmark. The quantum process reached about 833.5 MiB of resident memory at its cumulative high-water mark, which is not an isolated peak for either policy or case. The measured-evidence row's 256 MiB simulator-state limit applied to Aer's quantum-state storage, and its max_bytes of the same size applied to each QHD check. Neither limits total process memory.

The public default stays phr, with auto and slack as explicit choices, because the automatic rule has no solution-quality criterion. In the classical case with the same objective and the squared-coordinate inequality g = sum_i x_i**2 - 3/8 <= 0 at n = 3 on the binary periodic grid with K = 4, the best feasible normalized objective was 3/16 for forced slack and 1/8 for PHR, a loss of 1/16. The final forced-slack round selected s = 1/8 while the conditional minimizer was s* = 0. As Proposition 54 explains, the slack-grid bound concerns minimization over slack and does not by itself bound the inner suboptimality of a selected joint mode. Classical auto stayed with PHR in that case, and the quantum round-0 cost rule also rejected the slack trial.

If you choose auto or slack, compare the best feasible objective with constrained_grid_minimum when the unrefined original grid can be enumerated within its limits. A comparison with PHR alone can miss a gap shared by both formulations. In the classical case of the affine inequality at n = 3 on the one-hot Dirichlet endpoint grid with K = 8, PHR and forced slack both returned a best feasible objective of 0.301020408163, about 0.0816 above the original-grid minimum 43/196.

Continue an interrupted constrained run

solve_augmented_lagrangian(..., directory=path) saves a run in a directory, which must not exist yet, so that recoverable inner Runs can continue after an interruption. Each inner Run is a saved Run with its own run log (Continue an interrupted run). The rounds are the same as without a directory, except that each saved Run also stores its run log and inputs, which count against max_data_bytes, so a run whose data limit binds can stop at an earlier round. The directory's own files do not count against max_data_bytes, which caps the data that the inner Runs record, and no cumulative disk limit applies to the directory. Saved run mechanics describes the files in the directory and how they are written.

resume_augmented_lagrangian(path, backend=...) continues such a run. Its backend is that of the original call, and None stands for AerBackend() under quantum execution, as it did there. A backend with another configuration raises ValueError before anything is reopened, planned or committed, because the Runs that resume creates would otherwise run on another backend than the saved Runs, which the records could not show. A configuration carries no noise model, and a new AerBackend.from_noise_model(...) binding has another noise_model_id. For a noisy Aer run in a new process, pass AerBackend(noise_model_id=...) with the noise_model_id of the original binding, which that error names. Resume then uses the noise model saved in the directory, with its errors and basis gates, for the Runs it creates, so every round runs on the Aer target of the original Runs. A resume given the original bound backend, in the process that holds it, uses that backend's model for the Runs it creates. Treat the directory as read-only, like a saved Run folder (Saved folders are read-only).

Resume reopens the recorded inner Run with its recorded Plan and run log. For a round with an inequality representation, the outer record holds that representation from before the first inner Run. Resume reuses its forms, variable order, slack ranges and the content hash of the inner objective and validates the saved Plan association before continuing. Completed work stays in the outer record and is counted once. A recoverable Run can continue, including a remote result that becomes available after interruption. Some interrupted local preparation, measurement or classical evolution cannot be completed from the saved Run.

After successful recovery, the records equal those of an uninterrupted run with a directory and the same seed, apart from the fields that differ between any two executions, namely the Run and Result identifiers, the stored-data byte counts, whose text of wall times and measured timings varies in length, and the record identifiers that contain these fields. A run whose max_data_bytes is almost spent can therefore stop at another round, as two uninterrupted runs with a directory can. Every Run that the resume call creates or continues reports to its progress callback. For a run that has ended, resume_augmented_lagrangian returns its result without planning, evaluating or measuring, and result.save(path) then writes the archive that load_augmented_lagrangian reads.

The directory covers process interruption and does not synchronize files with fsync for power-loss recovery. The Run keeps the interrupted step counted against its limits and records an interrupted attempt as uncertain. See Continue an interrupted run for its recovery rules.

An unfinishable Run normally raises with a recovery note. With end_at_unfinishable=True, a handled unfinishable case after a completed round or level can return termination="inner_failed" with the completed results and counted work. If no round or level has completed, the error still propagates. A folder without a committed run-log header fails during reopen, and this option does not convert that failure to a terminal result.

For a headerless Run folder, follow the error's instruction about removing that named inner folder before retrying. Resume does not remove it or recreate it automatically. When an existing run log cannot be read, keep the folder and retry when it is readable. Its resource counts can be unknown.

A constrained-result archive and a run directory serve different purposes. Load a saved constrained-result archive with load_augmented_lagrangian. To continue a constrained run from its directory, use resume_augmented_lagrangian, then save its returned constrained result to an archive if needed.

Open a saved standalone refinement archive with load_box_refinement. Continue a refinement from its run directory with resume_box_refinement, then call save on its returned result to create a separate result archive. A headerless inner Run folder fails during reopen, including with end_at_unfinishable=True. Follow the recovery note for that specific folder.

Only one process at a time can continue a directory, and a second one raises at once.

Box refinement

refine_box repeats QHD on shrinking boxes, following Wu et al., arXiv:2605.12066v1, Sec. V. Each level is one ordinary QHD Plan, run with prepare and submit as solve runs it, so each level result keeps the Plan, observations and preparation records of a single solve. execution, shots, backend, seed and progress have their solve meaning, and each level plans with its own child of SeedSequence(seed). limits caps the whole refinement, and each level Run receives what the earlier levels left of the circuits, shots, stored data and synthesis work. The scientific notebook refines the shifted Ackley function of that paper and compares the result with the known minimum.

import sympy as sp

from nwqlib.algorithms.qhd import BoxRefinement, QHD, refine_box
from nwqlib.problems import Optimization

x = sp.Symbol("x", real=True)
problem = Optimization(
    objective=(x - sp.Rational(1, 5))**2,
    variables=(x,),
    bounds=((-1.0, 1.0),),
)
result = refine_box(
    problem,
    qhd=QHD(
        num_grid_points=4,
        num_steps=8,
        total_time=1.0,
        theory_flavor="split_step",
    ),
    options=BoxRefinement(),
    execution="classical",
    seed=7,
)
print(result.candidate, result.objective, result.termination)

(0.20000000000000018,) 0.0 box_unchanged

The refinement stops with box_unchanged, when the next box equals the current one, at a point within rounding of the minimizer 0.2.

print(result) summarizes a standalone BoxRefinementResult from its stored levels, selected point, termination and resources. result.report() returns the corresponding structured report. result.save(path) stores the original problem, refinement record and completed level Results in a new archive, and load_box_refinement(path) loads that archive without planning or execution. To preserve an in-progress refinement, pass directory=path to refine_box and continue it from that directory with resume_box_refinement(path, backend=...).

Each level divides every marginal by the valid mass (Eq. (12)) and keeps one index interval per variable. Following Eq. (13), the interval starts at the most probable index and adds the more probable neighbor until its conditional mass m_j reaches mass_threshold or it covers the axis. The paper breaks no ties. NWQLib starts at the smallest index among equal largest values and adds the lower neighbor on equal values, which makes the interval a deterministic function of the marginal (refinement._axis_interval). Each kept grid point represents its centered cell, whose faces are the midpoints between its coordinate and those of its neighbors, x_i - h/2 and x_i + h/2 on a uniform grid of spacing h. An interval that reaches an end index keeps that face of the box. The next box is the product of the kept cells. The paper's Eq. (14) lets each point represent the cell to its right instead, which moves each bound by half a cell.

On a periodic grid the interval does not wrap from the last index to the first, so the next box is a sub-box that does not cross the identified faces of the level box. When the distribution has mass on both sides of those faces, the interval grows across the axis and the side keeps its full width, so the result depends on where the period is cut. The next level keeps the periodic kinetic on its sub-box as a choice of the search, although the objective need not be periodic there.

With counts, print(result) and result.report() automatically give a lower bound on each selected region's mass in that level's backend-sampled distribution conditioned on valid decoding. The default confidence is 95%, simultaneous over the configured run. For one-hot readout, the distribution is conditional on one excitation in every variable register. A bound of one for the full valid grid does not mean that every raw outcome decoded validly. Binary readout treats every complete register outcome as valid.

report = result.report()                     # total failure probability 0.05
confidence = report["confidence"]
result.report(failure_probability=0.01)      # total failure probability 0.01
result.report(failure_probability=None)      # omit the statistical report

Fix the failure probability alpha before inspecting the counts to interpret the simultaneous confidence statement. The configured horizon H is max_levels for standalone refinement and max_iterations * max_levels for a constrained run, even after early termination. Each level uses its actual inner dimension, including slack coordinates. Exact readout, no completed counts levels, a constrained run without refinement, or explicit None gives confidence=None.

The report allocates half of alpha to Hoeffding bounds over every candidate axis interval and half to one-sided Clopper–Pearson bounds over every candidate joint box, then selects the larger available lower bound. This correction covers the box chosen from the same samples. At a stall split, the selected region is the chosen side including the valley, with every other axis whole. Proposition 49 gives the assumptions and proof of this coverage result (Wu et al., arXiv:2605.12066). The claim requires independent draws from the level's fixed conditional sampling distribution. It does not cover correlated shots, changes of sampling population within a batch or stopping shots according to the positions already observed.

Each RefinementLevel stores the integer number of valid draws S as valid_count, the raw returned count as returned_shots, the selected region's axis hit counts as region_axis_counts, and its joint hit count M as region_count. The report uses M/S for the selected region's empirical mass. Its confidence["levels"] entries give the region's indices, counts, both lower bounds, the selected bound and method, and any reason CP was unavailable. An augmented-Lagrangian entry also names the round. Computing the report reads these saved scalars without scanning observations or repeating a solve.

The existing axis_masses, joint_mass_bound and joint_mass describe the level's rounded observed distribution. joint_mass_bound evaluates the marginal union formula max(0, 1 - sum_j (1 - m_j)), while joint_mass uses the joint observations and can be larger because it includes their correlations. For counts, both describe the empirical distribution conditioned on valid decoding. At a split they still describe the ordinary interval box, which is the whole level box, and split.region_masses gives the two side masses. The coverage report instead uses the selected side's integer counts.

The displayed statistical bounds are binary64 approximations. CP can be unavailable outside its stated numerical range, in which case the report uses Hoeffding with its original half of the failure budget. Full-support mass is one by the definition of the conditional distribution. Neither a high empirical mass nor its population lower bound establishes that the global minimizer lies in the region, and bounds from different levels cannot be multiplied into coverage of the final box under the first level's distribution.

The sufficient shot budgets help interpret a requested mass. For example, with K=64, d=2, H=60 and alpha=.05, 2295 valid draws suffice for a CP lower bound of .99 when every valid draw hits the selected region. This is conditional on the observed all-hit event. It neither guarantees that a future batch will hit only that region nor gives a raw-shot budget for noisy one-hot readout. The report uses the counts obtained and does not change the requested shots.

point_rule selects the point each level reports: most_probable (the default, the candidate rule of Sec. V), best_observed (the level result's candidate) or mode_or_mean, which keeps whichever of the most probable point and the conditional mean position has the smaller recorded objective, compared as relative_objective, the most probable point on ties. The mean is generally not a grid point, so its objective is evaluated off the grid, and when that value is not finite and real the level keeps the grid point and records why in mean_unavailable. The rules are shared with the augmented-Lagrangian layer. For refinement within a constrained round, the inner objective is L_k under the PHR policy and L(x, s) with an inequality representation. A search-model level solves its positive normalization, while a physical level solves the inner objective itself. Comparisons between levels use recorded relative values of the inner objective.

Each level copies its inner readout's mode_status. This describes the joint grid-point tie rule independently of the marginal intervals that choose the next box. An unresolved mode leaves the ordinary refinement and stall-split rules unchanged. In an augmented-Lagrangian round, a selected mode representative carries the status of its inner Result, or of the selected best level with refinement. With slack variables, that status and the copied tie deficit describe the joint point before projection onto the original variables. They do not describe the mode or probability deficit of the marginalized original-variable distribution.

Every reported point and box face comes from the coordinates of the level's own grid (OneHotGrid.grid_value with the Method's boundary). A physical level reports and evaluates the objective at its grid coordinate, the point whose probability QHD reports. A search-model level evaluates F(a + D u) at its unit coordinate u, in the expanded form that its tables evaluate, which for exact coefficients is F at the exact image a + D u of that point up to the evaluation's rounding, and records u as unit_point. Its reported point is the image rounded once to binary64, a display value. F at the rounded point can differ from the recorded objective, or be undefined. For -1/(x - c) with c the binary64 number 1.2, on [1, 2] with K = 4, the first image is c + 2**-54, which rounds onto the pole c.

At its grid point each level reads the stored entries of its unscaled support tables and sums them and the constant with fsum as the classical kernel sums the tables, so it records exactly the tabulated value. At an off-grid mean it evaluates each expanded support term of the original objective separately and sums the terms and the constant with fsum. It does so at the grid coordinate or mean for the physical model and at the unit point for the search model, and the best point is the one whose level has the least recorded relative_objective, the objective minus one constant for the whole refinement (see the invariances below), the earlier level on ties. With exact readout, best_observed is the level's grid point with the least binary64 table value of the level objective, the smallest index on ties, whenever every grid point has positive probability, whatever the evolution did. The added-constant example below shows how these values can tie. The most probable point shows where the evolution concentrated probability, and the mode rule does not establish a mathematical minimum of the level objective.

scaling selects what each level solves, "search_model" by default. "physical" solves the original objective on the level box. Its kinetic weight 1/h**2 grows as the box shrinks, so the distribution can spread over the smaller box until an interval covers its whole axis and the box stops changing. "search_model" solves, on the unit box u in [0, 1]^d, the normalized objective V(u) = kappa (F(a + D u) - c)/E, where a is the lower corner and D the diagonal matrix of side lengths of the level box, and kappa is the gain BoxRefinement.gain, which is potential_gain or 8 when that is left unset. The unit grid has the same spacing at every level, so QHD applies the same dimensionless kinetic at every level (Proposition 48).

E is the sum of the ranges of the support tables of F(a + D u), rounded up, an upper bound on its range over the grid, and c is its constant term plus the table minima, so (F - c)/E lies in [0, 1] on the grid for these tables. Finding E and c takes one planning of the unscaled level objective before the planning of V, and table_evaluations in the level resources counts that extra table work. V is formed from the support expressions of F(a + D u) without its constant term, with the exact table minima and the exact kappa/E, so a large constant, as in 1e16 + x, leaves V unchanged whether SymPy holds it as an integer or as a binary64 number. The planning of V evaluates its tables anew. For tables rounded once, each value of V then lies within about 4 u kappa times conditioning of the value that the unscaled tables give, with u = 2**-53 and conditioning = table_magnitude/energy_scale, the sum of the largest absolute table values divided by E. This estimate illustrates simple expressions. It is not a bound for an arbitrary expression, and the classical kernel rounds each potential value once more when it sums the terms. A level's conditioning shows how much of its normalized potential can be rounding, and no stop bounds it.

A level whose table ranges sum to zero has E = 0 and nothing to normalize, and it stops the refinement as flat_objective. A level whose table ranges sum to a positive value no larger than the units in the last place of the tables' largest absolute values stops as unresolved_objective. This is a resolution rule. The tables cannot tell such a variation from evaluation error, and dividing by E would turn that error into a potential of order one. It does not show that the objective is constant. exp(2**-52 x) on the endpoint grid {0, 1} is strictly increasing and stops there, as does sin(x)**2 + cos(x)**2, which equals 1, on [-1, 1] with K = 6. Unresolved levels can also appear at deep levels of an objective such as cos(x), whose table values keep their size while their variation shrinks with the box. A polynomial objective keeps its variation, because the substitution x = a + D u is exact and its value at the corner a goes into the constant term rather than the tables.

potential_gain multiplies the normalized potential after E is formed, which gives the level Hamiltonian a(t) T_u + kappa b(t) (F - c)/E with the unit-grid kinetic T_u. Multiplying F itself would cancel in the normalization. A larger kappa narrows the ground state of a local well, but a finite-time evolution need not follow the ground state, and a larger potential difference also limits the transfer out of a basin in which the state starts, for constant coefficients, so the gain has no general monotone effect. On the double well (2 x**2 - 1)**2 + 3 x/5 + 2 (y - 3/10)**2 + 6 x y/5 on [-1.2, 1.2]**2, with K = 6, T = 10, 80 steps, the quadratic schedule with gamma = 0.3, the uniform initial state, threshold 0.99, the most probable point, classical execution, max_levels=10 and max_no_improve=10 (measured evidence), kappa = 1 stopped after five levels with best-point error 0.0826 in the infinity norm. At its last level both end cells of each axis held more than 1% of the probability, so every interval of mass 0.99 covered the whole axis. kappa = 8 kept shrinking for ten levels to error 0.00325, and kappa = 64 stopped after seven levels at 0.0791, or after six with 160 or 320 steps, where kappa = 1 and 8 kept their results. With the default max_no_improve=2 all three gains stop after three levels at 0.0826. No value is best for every problem.

The default 8 comes from a comparison of gains 1 and 8 on two-variable Dirichlet grids with the classical Schrodinger model, T = 10, 200 steps, the quadratic schedule, the kinetic initial state, threshold 0.99, the most probable point, max_levels=10, max_no_improve=10 and max_work = max_bytes = 10**10 (measured evidence), where gain 8 gave the smaller infinity-norm distance of the best point after ten levels on all three standalone problems. The distances were 5.0e-6 against 2.2e-4 on the same double well with K = 12, 3.0e-6 against 3.9e-4 on the anisotropic quadratic (x - 3/10)**2 + (y - 2/5)**2 + x y/2 on [0, 1] x [-1/2, 3/2] with K = 12, and 3.0e-8 against 1.4e-5 on the two-variable Ackley function of Wu et al., arXiv:2605.12066v1, Eq. (16), on [-5, 5]**2 with K = 32. Its minimizer is the binary64 point (0.961528396769231, 0.3696594771697266) of examples/generators/qhd_scientific.py. That paper's Sec. VI.A reports this seed-123 shift as x* ≈ (0.962, 0.370), and the point is the fourth of the successive uniform(-1, 1, size=2) draws from numpy.random.RandomState(123) that the authors' experiment script, which the paper does not link, makes one per test function, the Ackley draw after those of the quadratic, Rosenbrock and Rastrigin functions.

Search refinement uses a fixed unit-grid kinetic operator and a level-dependent normalized potential. Its Hamiltonian has the form a(t) T_u + kappa b(t) (F-c)/E. The levels share this form, but their sampled potentials can differ. Physical refinement keeps physical coordinates, so shrinking a box can increase the kinetic scale through the inverse square of the grid spacing.

The step count must resolve both the kinetic and potential evolution. A larger gain can increase potential phases and product-formula error, so the default gain does not establish that a chosen step count is adequate. The gain-8 convergence evidence uses classical Schrodinger evolution. For a circuit-product comparison, also check the product approximation and the numerical-reference limitations of ir_product.

With best_observed and the uniform initial state, gain 1 did better on the exact-readout double well and on one of three sampled Ackley seeds, whose 256 outcomes per level were drawn from each computed distribution with max_no_improve=2 (measured evidence), so potential_gain stays configurable. The physical model solves the original objective and has no gain, so it rejects an explicit potential_gain.

The search model is invariant under a translation of the box and objective, a positive scaling of each coordinate, a positive multiple of the objective and a constant added to it, for every fixed kappa (Proposition 48). Two such runs agree bit for bit when SymPy's exact arithmetic reduces their level objectives to the same expression, the multiple is a power of two, and every level's coordinates and faces map exactly, as with dyadic data. Otherwise they agree to rounding, which can change the choice at a tie or at the threshold. It is a different Hamiltonian at each level, not a change of variables of the original one. In u the original kinetic is -1/2 sum_j L_j**-2 d²/du_j², so a box with sides (1, 2) has physical kinetic weights (1, 1/4), where the model uses (1, 1). Search-model results therefore make no claim about the dynamics of the original QHD Hamiltonian.

These invariances concern the problem each level solves. The refinement decides between levels (the best level, the no-improvement count, the mode_or_mean choice and the split scores) by the recorded objective minus one exact constant C for the whole refinement, the constant term of the first level's expansion, and records these values as relative_objective and relative_scores. An added constant enters every level's constant term and C alike and cancels exactly when the constant and every coefficient of the objective are exact SymPy numbers, such as integers and rationals, so for the same level results it changes none of these decisions, in either model.

The level results themselves can change in the physical model with point_rule="best_observed", whose level point is the QHD candidate. The candidate is the observed valid point with positive weight that has the least evaluated objective (Readout), and QHD compares binary64 objective values that include the constant, so next to a large constant these values can tie at every grid point, and the candidate is then the lexicographically smallest grid-index tuple. Adding 2**60 + 1/7 to (x - 2/7)**2 + 3 (y + 1/9)**2 + x y/5 on [-0.9, 1.3] x [-1.1, 0.7], with K = 4, 10 steps over T = 4, mass threshold 0.7 and classical execution (measured evidence), made every physical level with best_observed report index (0, 0). In the same example the other two point rules of the physical model and all three rules of the search model, whose V holds no constant, gave the same levels as without the constant.

A binary64 Float anywhere in the objective makes SymPy form each level's constant term in 53-bit arithmetic at the size of the constant, so next to a large constant the compared values can lose the objective's variation as the reported values do. Writing such coefficients with sympy.Rational, which keeps their exact binary64 values, avoids this. A positive multiple of the objective multiplies the compared values, exactly when it is a power of two. The reported objective stays in original units and can lose the objective's variation next to a large constant. For example, (x - 1/3)**2 + 2**54 rounds to 2**54 at every point of [0, 1], where binary64 numbers are 4 apart.

With the default level_initial_state="configured" every level plans the level problem with the qhd argument, its initial state and preparation recipe included, so a GaussianState given as QHD.initial_state is read in the coordinates of the level problem. In the search model these are the unit coordinates u, so its center and widths are unit coordinates, not original ones, and the Gaussian keeps its place and widths relative to every level box. In the physical model they are the original coordinates, so the Gaussian stays at one point with fixed widths while the box shrinks, and a later box can exclude its center.

level_initial_state="best_point_gaussian" starts every level after the first from a Gaussian centered at the refinement's best point so far, as the code behind the published refinement results of Wu et al. does apart from its clipping of the center. Its width is level_gaussian_width times each side of the level box, 1/6 when left unset, so in the search model the Gaussian has width 1/6 and center (p - a)/L in the unit coordinates of the new box [a, a + L], with p the best point. The center is not clipped. A best point outside the new box, for example in a region that a stall split gave up, gives a Gaussian centered outside it, which the state still normalizes on the grid. The first level starts from QHD.initial_state, and each later level records its Gaussian as initial_state.

The option helps on some problems and hurts on others. On (x - 3/10)**2 over [0, 1] with K = 16, theory_flavor="split_step", T = 1, the quadratic schedule with gamma = 0.3, gain 8, max_levels=6 and max_no_improve=6 (measured evidence), a kinetic first level followed by best-point Gaussians ended at best-point error 0.00668, against 0.0534 with the kinetic ground state at every level. Both made six solves, and both results were the same at 128, 512 and 2048 steps. With the default max_no_improve=2 the Gaussians stopped after three levels at 0.112. On the double well above at gain 1, with K = 12, T = 10, 200 steps, a uniform first level, max_levels=10 and max_no_improve=10 (measured evidence), the Gaussians stopped after six levels at error 0.116, against 0.0062 with the uniform state and 0.00022 with the kinetic ground state at every level, and after a kinetic first level they reached 0.00137. In that two-variable comparison, which also covered the anisotropic quadratic and the Ackley function at gains 1 and 8, the kinetic ground state at every level gave the smallest error in all six cases, whether the Gaussians followed a uniform or a kinetic first level, so "configured" stays the default. Only the width 1/6 has been compared.

stall_split decides what a level does when its next box equals its box. When both end cells of an axis hold more than 1 - eta of its conditional marginal, every contiguous interval of mass at least eta contains both ends (Proposition 50). No interval rule can then shrink that axis, and when every axis of a level is kept whole the default, stall_split="none", stops the refinement with box_unchanged. With stall_split="best_region" such a level is split instead, at most max_splits times (default 1), unless it is the last of max_levels, since no later level would solve a region. It is not split either when the next level could not run. When the no-improvement count reaches max_no_improve or the remaining limits cannot cover a level, the refinement stops with no_improvement or budget_exhausted before any valley search, and when neither region a split can keep passes the width floor it stops with width_floor, in each case without spending an evaluation on the split.

The split looks for an index v of an axis marginal p that lies below a larger value on each side, and among these takes the one with the largest ratio min(L, R)/p(v), where L and R are the largest values of p below and above v. The ratio measures the rise from the valley to the lower of its two peaks, does not change when the marginal is rescaled, and bounds neither the discarded mass nor the objective in either region.

Only resolved valleys take part. With exact readout, the gap between the smaller peak and the valley must exceed the level's point-difference window (most_probable_tie_window) divided by the valid mass, plus the rounding of the marginal sums. With shots, the smaller peak's count must exceed the valley's count by more than sqrt(2 N log(2 H d K/alpha)), with N the valid shots, H max_levels, d variables, K grid points per variable and alpha = 0.01, which keeps the probability of accepting any valley that sampling cannot resolve at most 1% over the whole refinement when each level's shots are independent draws from one population. The threshold follows from Hoeffding's inequality for each cell's empirical mass (doi:10.1080/01621459.1963.10500830, Theorem 1, Eq. (2.3), p. 15, with the lower tail of Eq. (1.4), p. 13) and a union bound over the cells and levels (refinement._valley_admission). It concerns the population that the backend samples, so correlated shots, drift between pooled batches or a systematic device error against the ideal model need a different claim. The augmented-Lagrangian layer runs one refinement per round, so over R rounds the bound is R times 1%. A level without a resolved valley stops with box_unchanged, as without the option.

Each strict side of the valley is scored by the original objective at its most probable joint grid point, evaluated as the level's reported grid point is, which takes two objective evaluations. The next level solves the side with the lower score (compared as relative_scores) together with the valley cell, so the cut is the face of the valley's centered cell toward the other side, and the split discards only that side's probability. Both regions are unions of centered cells of the level grid and keep every other side of the box. The level records the split in split, a StallSplit with the axis, the valley, the cut coordinate, both regions with their scored points, scores, relative scores and probabilities, and the discarded probability.

A classical level keeps its joint distribution only in a kept state, so a classical refinement with splits needs keep_state=True. A periodic grid with a positive split budget raises ValueError, because one cut opens a periodic axis but does not divide it into two regions. An unchanged box after the last split stops the refinement with split_limit, whether or not that level has a valley. A stalled level that does not split states why in split_declined, namely that the level is the last of max_levels, that the split budget is spent, that the next level cannot run, that the readout has no tie window, or that no valley passes the resolution screen. The field explains the stop reported in termination and stays None without the option.

The split gives up the probability of the discarded region and with it the ordinary mass threshold. In exact arithmetic, each strict side of the selected valley holds at least the mass of the end cell that the greedy interval added last, so either region a split can keep holds less than eta whenever exact greedy growth reached eta only with that cell (Proposition 50). The greedy rule compares a running binary64 sum with eta, so a kept region can hold slightly more than eta. The split replaces only the box_unchanged stop, so for the same observations the levels before it are those of the run without it, and the best relative objective, taken over all levels, is never larger than that run's. The risk lies in the continuation. Each side's score is one objective value at the point where the evolution put the most probability, which is weak evidence about the side's minimum whenever the distribution is not concentrated near the minima. A narrow, deep minimum away from that point, for example in the discarded region, can make the discarded region the better one without changing any value the rule reads. A broad basin can do the same when the distribution is spread out, as on a coarse grid or at a small gain, because the most probable point of each side can then lie far from that side's minimum.

On the double well with K = 4, 20 steps, T = 10, gain 1, the uniform initial state and exact quantum readout (measured evidence), the first level's sides scored 1.441 and 1.487, and the split kept y <= 0.48, whose best grid value is -0.631, while the discarded region held the grid point (-0.72, 0.72) of value -0.700 and the continuous minimizer. A peak score therefore cannot justify discarding a region in general, which is why "none" is the default.

On the double well above, with the search model at kappa = 1, keep_state=True, classical execution, at most ten levels and max_no_improve=10 (measured evidence), the run without a split stops after five levels at best-point error 0.0826. With stall_split="best_region" and max_splits=1 the fifth level splits x between its peaks near x = -0.66 and 0.42, discards 2.4% of that level's distribution, and continues in the lower region to error 0.0108 after ten levels, against 0.00325 with kappa = 8 and no split. The intervals, the split and the best point were the same with 160 and 320 steps. A different trigger, which refine_box does not offer, would split as soon as one axis is kept whole. On the double well it would split x at level 1, and the same cut with the ordinary y interval reached error 0.000857 in ten solves after discarding 5.7% of that level's distribution. That trigger changes the run before any stall, so nothing bounds its best relative objective by the default run's.

Before each level, refinement stops, in this order, when

1. max_levels levels are complete (level_limit),
2. the next box equals the last one (box_unchanged, or split_limit after max_splits stall splits),
3. a side of the box is at its width floor (width_floor),
4. max_no_improve consecutive levels passed without a strict decrease of the best relative objective (no_improvement), or
5. the remaining limits cannot cover a level (budget_exhausted). A quantum level needs one circuit and the requested shots, and every other cumulative limit needs a positive remainder.

A level also stops it when its tables show no variation (flat_objective) or no resolvable variation (unresolved_objective), when its planning or Run raises (inner_failed, with the exception recorded), for example because the data the Run requests to store exceed what remains of the data limit, or when its result has no valid point (no_valid_point). Completed levels stay in the result, and the work of the stopping level's tables and Run is counted in its resources. When the planning or Run of the first level raises, nothing has completed, so the original exception propagates, as in the augmented-Lagrangian layer, and a configuration that QHD rejects surfaces at once.

The width floor is a resolution condition. A side is resolved while the rounded midpoint of each pair of adjacent grid coordinates lies strictly between them and the first and last midpoints lie strictly inside the faces. Every grid point of a resolved side then has its own binary64 coordinate, and the next box is not empty. A side of n grid intervals reaches the floor at a width of the order of n ulp(M), with M the largest magnitude of its coordinates. Distinct coordinates need not give distinct objective values, which near a smooth minimum stop differing well before the floor. The grid itself accepts a box only when every kinetic coefficient from 1/h**2 to 1/(4 h**2) is a normal binary64 number with a finite denominator and its binary64 coordinates are strictly increasing, so the physical model raises for a box such as [-1e308, 1e308], whose width overflows, or [0, 1e-200], whose h**2 underflows, and stops with width_floor when a later level's box fails this check. The search model solves every level on the unit grid, where this check always passes, and maps the unit points into the box, for example into [0, 1e-200].

The result's resources keep measurement work (circuit preparations and attempts, shots, stored Run data) apart from classical work (table evaluations, the classical kernel, circuit construction, exact synthesis and the refinement's own objective evaluations). max_plannings is the largest number of QHD plannings the options allow, two per search-model level and one per physical level. Each operation of a planning is checked against max_work and max_bytes of the QHD configuration, so the table, support, kernel and construction work each stay within max_plannings times max_work, a bound known before the first level. The best point is a finite-grid search result, not a continuous or global optimum.

Continue an interrupted refinement

refine_box(..., directory=path) saves a refinement in a run directory with each level's Run under levels/<z>/run/, and resume_box_refinement(path, backend=...) reopens its recorded inner Run with its recorded Plan and run log. Each level also stores its table-stage data, which a resumed level reads instead of evaluating the tables again (saved run mechanics). The outer record, its checks, the continuation from the last completed level, the progress of the resume call and the interruption scope with end_at_unfinishable are those of the augmented-Lagrangian layer (constrained problems). The level that follows the last completed one starts from the box, best level, no-improvement count, split count, spent limits and random stream of the uninterrupted refinement, and from a stop that the last completed level already decided.

Limits

• QHD reads out points of a finite grid model. A candidate, a most probable point or an evaluated grid minimum is not a continuous or global optimum, and no option implies a continuous global-optimization guarantee (Readout, Explicit verification).
• The binary encoding needs the periodic grid and K a power of two. The one-hot periodic grid needs an even K of at least 4, a limit of the current link compiler (Binary encoding, Scientific model and ordering).
• The spectral kinetic model needs the periodic grid and runs on the binary circuit, the binary theory_flavor="ir_product" and theory_flavor="split_step" only (Kinetic models).
• A classical evaluation of the model is not quantum execution, and by itself it is not an independent baseline (Scientific model and ordering).
• With initial_state_preparation="qiskit_state_preparation", which a Gaussian start needs on the binary encoding, the Plan records no CX count, and an excluded unitary instruction in the preparation record leaves the readout without a tie window (Circuit synthesis, Readout).
• The T estimate carries no error bar, and a requested synthesis budget is not an achieved error. Achieved synthesis and compiler errors are unavailable for every circuit, and a one-hot circuit with a number projector on 4 to 7 qubits or a binary circuit with a dense diagonal has no circuit error total (Fault-tolerant resources).
• The AQFT bound and the binary angle-formation bound assume at most one ulp of error in the sine and cosine functions they use, which has not been checked for NumPy's SIMD sin on Linux x86-64 hosts (Fault-tolerant resources).
• A stopping status of the augmented-Lagrangian layer does not assess optimality, and the default stopping test does not include projected stationarity (Constrained problems).
• Refinement masses describe each level's own distribution and do not show that the global minimizer lies in the next box, and search-model results make no claim about the dynamics of the original QHD Hamiltonian (Box refinement).
• A run directory covers process interruption, not power loss (Continue an interrupted constrained run).
• Two readout cases are not yet checked, the tie window with the NWQ-Sim per-instruction constant and the readout roundoff of a QHD circuit that has ancilla qubits beyond its register (Readout).

Limitations and open work lists the current limitations and the open work of QHD.

Terms

This guide, the augmented-Lagrangian layer, box refinement and their docstrings use these words in narrow senses.

Term	Meaning in NWQLib	Code
layer, round, level	Two senses of layer. A layer of QHD is a function that runs a sequence of QHD Plan objects, solve_augmented_lagrangian or refine_box. In the one-hot kinetic compiler a layer is a set of disjoint links applied as one commuting factor. A round is one step of the augmented-Lagrangian layer, its QHD solve or refinement with the multiplier and penalty update that follows. A level is one QHD Plan of a box refinement on one box, while the transform level count L of theory_flavor="split_step", ceil(log2 K) for the FFT (split_step.transform_levels), is an unrelated quantity.	algorithms/qhd/_outer.py, algorithms/qhd/kinetic.py KineticCompiler
layer work	The augmented-Lagrangian layer's preprocessing and original-function evaluations, counted against its limits. Preprocessing counts expansion monomials, support-table work K**r (N+r) and required original summand scans on their own supports, with a limit-checked whole-term scan when needed. A point evaluation of a term with N tree nodes in d original coordinates costs N+d units. The layer checks the preprocessing and all the round checks it sets aside against its limits once for the whole run. QHD planning and representation trials have separate checks.	algorithms/qhd/constrained.py _setup, _admit_layer_work
effective objective L_k	The normalized PHR objective of round k, with the multiplier and penalty terms of Eq. (10.3) of Birgin and Martinez, doi:10.1137/1.9781611973365. It governs the projected outer update. A round with an inequality representation evaluates it from the original functions at its projected point for ALEvaluation.effective_value.	algorithms/qhd/constrained.py _effective_objective, _effective_value
inner objective L(x, s), slack coordinate	The objective the round's inner problem uses, with the chosen PHR, quadratic, constant or slack term for each kept inequality. Converted inequalities add normalized nonnegative slack coordinates. Point selection and refinement compare this inner objective before projection to the original variables. Proposition 54 states its relation to L_k.	algorithms/qhd/constrained.py _inner_objective, _round_problem, algorithms/qhd/constrained_records.py InnerRepresentation
inner representation, cap scope	InnerRepresentation records each kept inequality's form, the ordered slack axes, grid convention, trial costs and the content hash of the inner objective. Its cap_scope="grid" refers to the preprocessed original grid and the tabulated representation, with the cap and arithmetic assumptions of Proposition 54. Register and refinement dimensions include its slack axes.	algorithms/qhd/constrained_records.py InnerRepresentation, SlackAxis, FormTrial
slack cap, slack box	For a kept normalized inequality whose preprocessing tables give the lower bound ell_j, the cap is U_{0,j} = [-ell_j]_+. It contains every required continuous slack on the grid of the preprocessed box for every nonnegative entering multiplier and positive penalty. A zero cap creates no slack axis. A Dirichlet slack box has upper end U = U_0. A periodic slack box uses the upward-rounded margin U >= K/(K-1) U_0, so the ideal last node U(1-1/K) reaches U_0. The cap conditions and the arithmetic qualification are in Proposition 54.	algorithms/qhd/constrained.py _inequality_range, _slack_axis, algorithms/qhd/constrained_records.py SlackAxis
inequality policy, automatic selection	AugmentedLagrangian.inequality_form sets the policy, which InnerRepresentation.policy records. phr, the default, keeps the PHR term for every kept inequality. slack requests explicit slack coordinates after the eligible quadratic and constant branch tests. auto converts only under execution="quantum", in one deterministic pass through kept inequalities in original order. It accepts a trial only if planned classical work and CX per circuit are both no larger and at least one is strictly smaller, or when a trial whose planning passed, with a known CX count, replaces a current representation whose planning was refused. Under execution="classical" or with a user-supplied GaussianState, automatic selection keeps every inequality in PHR form.	algorithms/qhd/constrained_records.py AugmentedLagrangian, algorithms/qhd/constrained.py _round_problem, _planned_costs, _accepts
planning-only entry	plan_augmented_lagrangian(problem, *, qhd, options, execution, shots, seed) returns (preprocessing, representation, plan) for an unrefined round 0. The representation is None under phr. The call performs preprocessing, representation selection and QHD planning, with no Run, preparation, measurement, evolution or grid-reference solve. It does not predict a refined first level, later adaptive representations, an executed trajectory or exact future costs.	algorithms/qhd/constrained.py plan_augmented_lagrangian
tentative and safeguarded multipliers	The tentative multipliers are lambda+ = lambda + rho h and mu+ = max(0, mu + rho g) at the point of a round. The next round uses them clipped to MultiplierBounds, the safeguarded multipliers, or unchanged without bounds.	algorithms/qhd/constrained.py
normalized residual, normalized infeasibility	A residual divided by its scale, h_i/s_{h_i} or g_j/s_{g_j}. The normalized infeasibility is max(||h||_inf, ||g_+||_inf) over the normalized residuals.	algorithms/qhd/constrained.py
complementarity measure, penalty measure	Both are computed from normalized residuals. The complementarity measure max(||h||_inf, max_j |min(-g_j, mu+_j)|) with the tentative multipliers is the left side of the stopping test of Birgin and Martinez, Eq. (10.7). The penalty measure max(||h||_inf, ||V||_inf) with V_j = min(-g_j, mu_j/rho) decides whether the penalty grows (Eq. (4.9)).	algorithms/qhd/constrained.py _update
point rule	The rule by which a layer reads one point from a QHD result, most_probable, best_observed or mode_or_mean.	algorithms/qhd/_outer.py
joint point, projected point	A converted QHD inner problem reads coordinates (x, s). The outer point is its projection x. The stored point probability and mode tie diagnostics describe the joint point, while f, h, g and the PHR value are evaluated at x.	algorithms/qhd/constrained.py _projected, algorithms/qhd/constrained_records.py ALEvaluation
candidate, most probable point	QHD evolves a state on a finite grid and reports its valid readout outcomes. The most probable point describes the readout's mode representative. The candidate is the observed valid point with positive weight that has the least evaluated objective. These can be different points. The most probable point can be below the computed maximum within the numerical tie window. Wu et al., arXiv:2605.12066v1, Sec. V, call the grid point of largest probability the candidate solution, which NWQLib calls the most probable point.	algorithms/qhd/method.py
tie window, tie-window ceiling	The bound probability_difference_window on the error of a computed difference of two probabilities. The mode representative is the first observed valid grid point with positive weight in lexicographic order whose computed probability deficit from the maximum over that same set of points passes the numerical tie window. Its computed probability can therefore be below the maximum. A classical kernel records the window of its Plan's state budget as tie_window_ceiling, and the window it observes must not exceed that ceiling.	_validation.py probability_difference_window, algorithms/qhd/method.py
observed state budget	The state budget of theory_flavor="split_step" that weights each phase-angle error by the state's own mode and point probabilities. It is at most the Plan's budget.	algorithms/qhd/split_step.py evolve
absorption, preprocessed box	Absorption turns a one-variable linear inequality into a box bound (absorb_bounds=True). The preprocessed box is the original-coordinate box after absorption. Each round starts from that x box and appends its selected slack boxes to form the inner box.	algorithms/qhd/constrained.py
run directory, outer record	A run directory is the directory given to a layer as directory=. It holds the inner Runs, the problem and the outer record, controller.json, which keeps the arguments, the configuration of the backend of every inner Run and the completed rounds or levels, together with a constrained round's committed inequality representation, and is committed by write-then-rename.	algorithms/qhd/_durable.py
unfinishable Run	An inner Run that cannot finish without new work counted against its limits, because its preparation failed or was interrupted after its work was counted, or because its measurement or classical evolution cannot deliver its outcome, such as an interrupted local evolution or a lost acknowledgement that the backend cannot reconcile. An unfinishable Run normally raises with a recovery note. With end_at_unfinishable=True, a handled unfinishable case after a completed round or level can return termination="inner_failed" with the completed results and counted work. If no round or level has completed, the error still propagates. A folder without a committed run-log header fails during reopen, and this option does not convert that failure to a terminal result.	algorithms/qhd/_durable.py unrecoverable
search model, potential gain	A search-model level solves the dimensionless objective kappa (F(a + D u) - c)/E on the unit box, and kappa is its potential gain.	algorithms/qhd/refinement.py
relative objective	F - C, the objective minus the constant term C of the first level's expansion, whose recorded values refinement compares between levels and split regions. Within a constrained round, F is the round's inner objective, L_k under the PHR policy and L(x, s) with an inequality representation. The least recorded relative value wins across completed levels, with the earlier level winning an exact tie. Within a level, best_observed selects the least observed table value of that level's solved objective, including its search normalization when selected.	algorithms/qhd/refinement_records.py RefinementLevel
level conditioning	The sum of the largest absolute table values of a level divided by E, the sum of their ranges. For a search-model level it shows how much of the normalized potential can be rounding, and for a physical level how much of the potential variation can be rounding.	algorithms/qhd/refinement_records.py RefinementLevel
width floor	The side length at which the rounded midpoints of adjacent grid coordinates no longer all lie strictly between them, with the first and last strictly inside the ends of the side. A level at its width floor stops the refinement.	algorithms/qhd/refinement.py
stall split, valley, resolution screen	A stall split divides a level whose box stops changing between two peaks of one axis marginal, at the valley between them, and continues on one side. The resolution screen accepts only valleys that the readout resolves, by the tie window or, for counts, a Hoeffding threshold.	algorithms/qhd/refinement.py
support expression, support table	A support expression holds the terms of the expanded objective whose variables form one support set S. Its support table holds its values on the grid tuples of S as one float64 array in C order, first variable of S most significant, which planning evaluates once and every grid computation reads by index, together with the minimum, maximum and largest magnitude taken from it.	algorithms/qhd/objective.py ObjectiveDecomposer, algorithms/qhd/records.py SupportValues
step weights, step averages	The step weights (t_k, a_k, b_k) are the schedule values that step k applies under the coefficient rule. The step averages A_k/dt and B_k/dt of the interval integrals are the weights of the integrated rule.	algorithms/qhd/schedules.py step_weights
kept state	The final state stored with keep_state=True, whose global phase the phase rules of planning govern.	algorithms/qhd/method.py QHD
running phase total, identity event	The running phase total is the physical phase of a compiled QHD product, the sum of the identity parts that its blocks omit, with its rounding allowance. On the binary encoding it holds the objective constant's phase alone, since the diagonal blocks keep their identity phases. An identity event is one of those contributions, such as the identity part of one number projector. The objective constant enters the running phase total as a global phase, so a large phase-sensitive bound can include error that leaves the ideal readout probabilities unchanged and does not by itself establish poor readout accuracy.	algorithms/qhd/compiler.py, algorithms/qhd/method.py
error sources	The entries of QHDCircuitResources.error_sources, one per error source of a circuit against its finite model, each with its status.	algorithms/qhd/resources.py
normal range, range rule, range check, lower-range omission	The normal range holds the binary64 magnitudes from the smallest normal number, 2 to the power -1022 (sys.float_info.min), to the largest binary64 number (sys.float_info.max). The range rule requires a formed product, quotient or power-of-two scaling to be exactly zero or in the normal range, and the range check is planning's check of that rule. A lower-range omission drops a contribution formed from objective data or the initial state whose exact value would be nonzero and below the normal range, and adds its exact action to the error sources.	algorithms/qhd/validation.py, algorithms/qhd/records.py QHDRangeOmissions
rotation classification	The classification of a circuit's rotations into Clifford gates, exact T gates and arbitrary rotations by their binary64 angles.	algorithms/qhd/resources.py rotation_population
Walsh block	One occurrence of a Walsh-rotation diagonal in a binary step, a support table or a kinetic phase.	algorithms/qhd/binary.py
link, wrap link	A link is a pair of neighboring grid points of one variable, which one hopping term of the one-hot kinetic operator couples and one step of the amplitude chain connects. The wrap link joins x_(K-1) and x_0 of the periodic grid.	algorithms/qhd/grid.py
structured preparation, amplitude chain	The structured preparation is the library's exact initial-state circuit of each encoding. On the one-hot encoding it is the amplitude chain, one controlled RY and one CX per link.	algorithms/qhd/initial_state.py append_amplitude_chain

Numerical guarantees

This section collects the roundoff budgets behind the phase rules for kept states (Scientific model and ordering), the mass window of each classical kernel and the tie window of the most probable point (Readout). The docstrings of the named functions give the full derivations.

Objective constant and running phase total

On theory_flavor="schrodinger" and theory_flavor="split_step" the objective constant c contributes the global phase exp(-i c sum_k dt b_k), which the kernel applies to the kept state after it reads the probabilities. With finite intermediates, normal relative-error arithmetic and correctly rounded math.fsum, forming the constant's phase angle has error at most (3u + 3u**2 + u**3) |c| sum_k |dt b_k| rad, about 3u |c| sum_k |dt b_k|. If the final nonzero phase product is below the accepted normal range, the kernel replaces it by zero and adds its exact omitted magnitude to the phase-error budget. Evaluating the exponential and multiplying the kept state have their own roundoff allowances (method._constant_phase). The probabilities are read before this global phase is applied. As Scientific model and ordering states, with keep_state=True planning rejects a constant whose phase-error budget, evaluated upward, reaches pi, because the budget then gives no guarantee on the kept state's phase (method._admit_constant_phase).

The kept states of the quantum circuit and the IR-product reference follow the same pi rule with the phase-error budget of their running phase total, the compiled identity phase, which is the one-hot sum of the kinetic diagonal, the projector identities and the constant, or the constant alone for the binary encoding. The budget counts contribution formation and accumulation, using Neumaier's compensated sum for the one-hot running phase total and math.fsum for the binary one. It also adds the magnitudes of identity contributions omitted below the normal range. On the circuit route it additionally counts the circuit's phase assignments and their reduction modulo the binary64 tau (method._compiled_phase_allowance, method._native_phase_allowance, method._admit_compiled_phase). Proposition 42 derives this budget.

A binary diagonal keeps its identity phase inside its block, so the binary routes add two terms from the Walsh diagonals. Identity-phase formation is counted on both routes. Let X be the exact intended exponent, x its stored value, c0* the exact intended identity coefficient, c0 its accepted computed value, and phi the emitted identity phase. With eps_x >= |x - X| and rho0 >= |c0 - c0*|, the budget is |X| rho0 + eps_x |c0| + |phi + x c0|. For a potential table on n qubits, use rho0 = gamma_n a + d0, where gamma_n = n u/(1 - n u), u = 2**-53, n u < 1, a is the mean absolute table value, and d0 is the error of omitting the identity coefficient at normalization. Kinetic tables use the radius from binary.kinetic_identity_enclosure. For a potential table without normalization or phase omission, normal finite arithmetic gives the first-order estimate (n + 2) u |x| a. The general bound also covers cancellation and omitted identities. The circuit route also counts each Walsh block's global-phase assignment, and the IR-product route the reconstruction of each diagonal's phases from its kept rotation angles (binary.WalshPhaseTerms).

The budget is evaluated upward at every operation. For fixed contribution counts, the first-order formation-and-accumulation budget is 3u (S_C + S_P) + 5u S_K, where the three S terms are the absolute sums of the constant, projector and kinetic contributions. The implemented finite bound also contains gamma_(M-1)**2 for M phase contributions and the corresponding term for summing the kinetic coefficient over the variables. These terms depend on the counts and need not be negligible for long schedules. Omitted identities, binary block phases and circuit phase assignments carry their additional budgets. For the two-dimensional Ackley replication with K = 32, 100 integrated steps of ShiftedCubicSchedule(s=2e-4) (measured evidence), whose phase of about -2.2e11 rad is almost all kinetic diagonal, it is 1.2e-4 rad, while c = 1e17 over two unit steps gives 67 rad on the IR-product route and 111 rad on the circuit route, so that kept state is rejected.

Tie window

The tie window of Readout bounds how far a computed difference of two probabilities can lie from the difference in the model that the path evaluates. It is formed from a state budget delta and the relative error e of each computed probability. If the computed state differs from that model's unit state by at most delta in 2-norm, the window is 2*delta + delta**2 + e*(1 + delta)**2 (Proposition 45). The selection compares the computed M - p with the window (method._summarize). The subtraction is exact when p >= M/2 (Sterbenz's theorem, Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., doi:10.1137/1.9780898718027, Theorem 2.5). For smaller p its rounding is at most u (M - p), and the comparison adds no budget for it. Rounding is monotone and the window is a binary64 number, so an exact deficit at most the window always rounds to at most the window, while a deficit just above the window can round down to it and be accepted. Readout states the resulting bound on the probability deficit of the selected point.

State budgets

Every classical kernel evaluates its probabilities as np.abs(state)**2, so e = 5u. The Schrodinger kernel takes delta from the same first-order budget as its mass window, described in the next paragraph. A one-hot IR-product Plan takes delta from the direct block budget below (theory.onehot_product_state_error), and the binary IR-product kernel from the state-operation budget on computed phase arrays and selected stored QFT angles (theory.binary_product_state_error), in each case applied to the exact initial state.

The Schrodinger state budget compares with the ordered exponentials of dt*(a*Khat+b*diag(Vstar)), where Khat is the stored kinetic matrix and Vstar is the exact sum of stored support entries. Its constant stored kinetic diagonal is reproduced by an ordered scalar sum. The budget adds broadcast summation error and generator-formation error to the SciPy numerical allowance. Generator formation includes the three-rounding relative term and an absolute budget for subnormal products under round-to-nearest binary64 arithmetic with gradual underflow and finite intermediates. The max_work check uses the same shifted generator norm. Unitary propagation adds the perturbations in state 2-norm. The result remains conditional and first order under the assumptions of expm_multiply_state_error, including its computed-trace-shift assumption.

If the intermediate abs(total_time/num_steps)*(abs(a)*J+abs(b)*W_vec) overflows before the gamma_3 factor, lower total_time or increase num_steps at fixed weights, scale the objective down, or widen the box at fixed num_grid_points to reduce the stored kinetic bound J, rechecking the potential magnitude bound W_vec and all recomputed schedule weights so the unscaled products and sum are finite and the intermediate is below the binary64 maximum with room for outward rounding. The table/assembly or outward centering refusal reports abs(dt), the step weights a and b, J, and W_vec at the failing step.

The binary ir_product flavor's first-order roundoff budget, relative to the computed phase arrays and selected stored QFT angles, is u (s + sum of block terms + 5), where s is the construction bound of the start vector (2 for the uniform state, initial_state.restricted_state_error), 3 per H gate, 5 per controlled phase and per diagonal, 2 (3 b + 5 C_p) + 5 for a kinetic block with C_p controlled phases per QFT, 5 for a potential block and 5 for the final phase. It gives the kernel's mass window and most-probable tie window, as the state budget of each other theory_flavor value does for its kernel.

The one-hot product reference is the ordered product of analytic projector and hopping blocks at their stored binary64 angles, applied to the exact initial state. Under normal round-to-nearest arithmetic and the stated one-ulp elementary-function assumption, its first-order state-operation budget is u*(start+10*B_P+6*B_H+5) (theory.onehot_product_state_error). The counts include every executed occurrence and the final term covers the physical-phase application. Parameter formation, circuit construction, time splitting and grid errors are outside this budget. The mass and mode windows use the same state-to-probability conversions as the other classical kernels.

The binary IR floating-point budget describes state operations on computed phase arrays and selected stored QFT angles. Its first-order arithmetic model excludes phase-array formation and the discrepancy between a reconstructed Walsh array and the exact action of the stored rotations. The kept-state phase-error budget includes R_W, but R_W is not propagated into IR readout accuracy, the mode tie window or verification uncertainty.

The split-step kernel's mass window and the upper limit of its tie window come from the Plan's state budget, split_step.state_error (Proposition 41), relative to the exact split-step product with the stored weights, the exact eigenvalues and the exact sum of the support tables, applied to the exact initial state. Like the budgets of the other theory_flavor values, it starts from the construction error of the initial vector (initial_state.restricted_state_error). Each step counts 5 u L for each of its 2d transforms and 5u for each of its d + 2 phase multiplications, with u = 2**-53. Rounding its phase angles adds 15 u |dt a_k| max E per variable and (M + 1) u |dt b_k| W for the two potential halves, with W = sum_S max |v_S| over the M support tables v_S. The transform term is a stated assumption about SciPy's transform backend, ducc0.fft in SciPy 1.18.1, ||fl(T x) - T x|| <= 5 u L ||x||, which a test checks against 80-digit sums on sample lengths without proving it for every length. The budget is first order in u. Its relative-roundoff terms assume round-to-nearest binary64 arithmetic without overflow, with normal nonzero quantities wherever a relative error is used, so an evolution whose intermediate values underflow would need absolute terms that the budget does not include.

The angle terms dominate when the kinetic integral is large near t = 0. With ShiftedCubicSchedule(s=1e-3) and the integrated rule on a five-point grid of spacing 1/2 (measured evidence), the first kinetic angles reach 3e7 rad, the kernel's state from the kinetic ground state differs from a 30-digit evaluation of the same product by 4.3e-10, and the Plan's budget is 5.0e-8. For the s=1e-3 five-point example, the trajectory-dependent budget observed by split_step.evolve is 3.6e-9. It weights angle errors by the state's mode and point probabilities. The Plan's 5.0e-8 budget instead uses an upper bound available before evolution.

state_mass_window and probability_difference_window convert the budget into the mass window and the upper limit of the tie window exactly as for the other theory_flavor values. A budget above about 1.3e154 makes both windows overflow although the evolution can finish, so planning refuses every classical Plan whose budget or windows are not finite binary64 numbers, with a ValueError that names the window and the budget (method._admit_host_windows). A grid whose K**d exceeds the binary64 range is also refused by a ValueError. It names the limit that refuses the Plan first, max_work or max_bytes of the kernel's check, and otherwise the grid size and num_grid_points.

For f(x)=(x-1/3)**2 on eight periodic spectral points of [-1,1), with GaussianState(center=(0.2,), widths=(0.3,)) and three integrated split steps (measured evidence), ShiftedCubicSchedule(s=1e-150) to T=1 gives a Plan budget of approximately 5.3e287 and refuses because the required windows are nonfinite. With the same objective, grid, initial state and step rule, QuadraticSchedule(gamma=0) to T=1e150 gives a budget of approximately 1.3e137 and solves. A finite budget this large does not give a useful accuracy guarantee.

The split-step kernel takes delta from the budget that it observes on its own trajectory (split_step.evolve), relative to the exact split-step product applied to the exact initial state. It counts the same transform and phase-multiplication roundoff as the Plan's budget, but it weights each phase-angle error by the state that the phase acts on, 15u |dt a_k| ||E_j z|| for the kinetic phase of axis j, with z the normalized state after that axis's forward transform, and u |dt b_k/2| (2 ||V z|| + (M - 1) W) for each potential half, with W the sum of the largest magnitudes of the M support tables. The potential phase of one step's second half keeps every magnitude of its input, so the first half of the next step reuses that input's moment ||V z||, and N steps need N + 1 potential moments. Each weighted term is a root mean square of the angle errors over the state's own probabilities (Proposition 41). This budget is at most the Plan's budget (split_step.state_error), which counts the largest eigenvalue and table magnitude against the whole state and still sets the mass window and the upper limit of the tie window.

Every classical kernel records its delta as the scalar readout_state_error and records the Plan's window as tie_window_ceiling in its KernelApplication, and analysis and loading convert the scalar again and require the window to be at most that limit, without evolving anything. The budget starts from the construction error of the initial vector, 2u for the uniform state and for the kinetic ground state on the periodic grid, about (9d - 1)u for the kinetic ground state on the Dirichlet grids and, for a Gaussian, a finite bound computed from its exponents that grows with their size (initial_state.restricted_state_error).

Exact circuit probabilities take delta from the preparation record, (G*c/2 + 5)*u for the record's G circuit operations (on Aer, without save instructions) with the simulator's per-instruction roundoff constant c and u = 2**-53, and use e = 2u, plus a pooling term when several chunks are averaged. Stored circuit amplitudes use the same delta with e = 5u for np.abs(z)**2. These windows are first-order worst-case bounds.

When several chunks are pooled into QHD's masses, analysis divides each of the M stored bin values by the chunk count (one rounding), adds it into running sums (at most M additions per value) and forms the valid mass with one fsum (one more rounding), so each mass moves by at most (M + 2)*u, with u = 2**-53. Each chunk's own total lies within the window of its preparation record, which bounds the simulator and readout roundoff. This allowance widens the mass check only. The tie window of the most probable point has its own pooling term gamma_C = C u/(1 - C u) for the C roundings of each point's pooled probability. M counts the stored entries of each probabilities or counts chunk, which include stored zero values, so M is not the number of nonzero values. Another chunk adds its number of stored values.

A preparation record with an excluded unitary instruction has no tie window (Readout). The instruction count does not supply a window, because the per-instruction roundoff allowance of the derivation does not cover an arbitrary dense unitary. On a register of at most 6 qubits (b <= 64 amplitudes), the normwise rounding bound of applying the matrix, about 2 sqrt(2) b**1.5 u, stays below Aer's roundoff allowance of about 3479u per instruction, with coefficients 181, 512 and 1448 at b = 16, 32 and 64. There the obstacle is the error of the matrix that Qiskit supplies, which the record cannot bound. From 7 qubits on, the rounding bound alone exceeds the allowance, with coefficient 4096 at b = 128 (Limitations and open work).

Window sizes

The split-step window is 4.1e-14, against the Plan's 5.7e-14, for the objective x³/8 - x²/2 + x/4 on the five-point interior grid of [-1.5, 1.5] with three steps over total time 0.75 under the quadratic schedule (measured evidence). On the periodic spectral grid of [-1, 1) with K = 16, the objective (x² - 1/2)² + 5e-8 x, ShiftedCubicSchedule(s=2e-4), 1000 integrated steps and total time 10 (measured evidence), the first kinetic integral is about 1.0e8. The Plan's window, 1.05e-4, exceeds the gap of 6.65e-5 between the most probable points of the two wells and would let the lexicographic tie rule select the lower well. The observed window is 1.3e-11, a 40-digit evaluation of the same product puts the largest error of a computed probability at 1.5e-14, and the selection is the reference maximum. Most of that window is the transform and phase-multiplication term, whose constants are qualification assumptions. The angle terms can also stay far above the actual error, because an angle error on a mode that holds nearly all of the probability acts almost as a global phase, which the budget still counts.

On a 7-by-7 Dirichlet interior grid of the disk potential of the comparison of the classical evaluations at rho 1, with ShiftedCubicSchedule(s=2e-4), 200 integrated steps and T = 1 (measured evidence), the observed window is 2.9e-6 from the uniform state and 8.1e-7 from the kinetic ground state, against the Plan's 2.1e-5, while a 50-digit evaluation of the same product put the largest error of a probability difference at 4.8e-9 and 7e-16.

Normal binary64 range

Planning establishes the range assumptions of the error sources (validation._normal_range). It checks that durations, weights, coefficients, angles and identity products are zero or normal binary64 numbers, where each is formed and before pruning, so half steps and projector scalings are exact. It also checks that the grid coordinates are distinct and increasing and that a Walsh diagonal's reconstruction cannot overflow. The stored tables, the objective constant and the initial amplitudes need only be finite. A contribution formed from them whose exact value would be nonzero and below 2**-1022 is omitted, and its exact action is counted in one entry of the error sources. A one-hot projector is omitted together with its identity phase, counted in rotation_pruning and identity_phase, and a constant contribution in identity_phase. A binary nonidentity Walsh coefficient is counted in angle_formation, a Walsh rotation or dense phase entry in rotation_pruning, and a Walsh identity coefficient or identity phase in identity_phase. A structured one-hot chain stops before the first link whose parameter theta/2 would be that small, counted in state_preparation. QHDReconstruction.range_omissions reports the counts and their errors (its *_charge fields) with the policy name lower_range_omission/1. For example, -exp(-100 x**2) on [-5, 5] with K = 64 has grid values near 1.6e-315 at the box edges, and its quantum Plan omits those projectors in both potential halves of a step, at a total error near 3.2e-315 per step.

Kinetic coefficients, energies and one-hot kinetic phases are not omitted, while the rotations and phase entries of a binary kinetic block follow the same product rule as a potential block's. A grid spacing h and kinetic weight a that put a/h**2 outside the range raise with that remedy, for example K = 64 on a periodic binary box of width 6.4e154. Under the integrated rule, a potential interval integral or step average outside the range names its remedy too, fewer steps or a longer total_time below the range and a shorter total_time above it, because b does not decrease in time and more steps reduce a late step's integral but not its average. For example, one step of either cubic schedule with s = 1 at total_time=2e103 is refused with a shorter total_time as the remedy, and total_time=1e70 plans. A kinetic integral or step average outside the range names the quantity, the step and the schedule but no remedy, because a short first step, late times and a large s can each put it below the range, and a small s above it. Every other range failure raises and names the quantity, for example a step total_time/num_steps below 2**-1021, a table 1e308 x whose Walsh butterfly overflows or a box [1e16, 1e16 + 4) whose coordinates coincide. Error bounds themselves may be subnormal and are rounded upward.

Planning limits

QHD checks the work and bytes of each stage against max_work and max_bytes before it runs the stage. This section gives those counts. The work counts are units of planning work, not measured CPU instructions or elapsed time, and the byte counts are declared amounts, not measured peak memory (engineering constants).

Tables, schedule rows and initial state

Planning checks the evaluation of the initial amplitudes and their stored payload, dK units of planning work and the byte count below (method._initial_state_bytes), and runs it before any table work. With it, before the compiler forms any schedule row, planning checks the rows, 2 units per midpoint step or the schedule's interval_work per integrated step, and 768 bytes per step for each row and its share of the records (method._STEP_BYTES), and it counts 4096 bytes for each compiled block occurrence (method._BLOCK_BYTES). Planning and the Schrodinger kernel read the generator norms of the steps one at a time (method._generator_norms), so no per-step norm is held.

Support tables are evaluated in bounded C-order slabs. Their stored float64 entries define the grid objective. The planning byte count includes the supported expression evaluator's workspace as well as the final table payload.

QHD planning limits both support-table construction and the temporary storage used to compute the content hashes of records. The limit check includes the stored arrays, support metadata, evaluation workspace and array serialization. Construction and serialization use separate peak calculations because their temporary buffers have different lifetimes.

The initial-state byte count was calibrated with 64-bit CPython 3.12.14 and NumPy 2.5.2 on macOS arm64. For \(d\) variables and \(K\) points per variable, the byte count is

\[ B_{\rm initial}=336dK+444d+104+4096(d+1)\ {\rm bytes}. \]

The coefficient multiplies the \(dK\) one-variable grid entries, not the \(K^d\) joint grid size. It is a declared count of Python objects and arrays, calibrated on that environment, not a platform-independent peak-memory theorem.

Classical kernels

Planning counts each expm_multiply call of theory_flavor="schrodinger" by the shifted norm of its assembled generator, dt (a H + b R) plus twice the step's table and assembly perturbation, where H is the stored kinetic column bound and R bounds the sum of the support-table ranges of the objective, so the work of one evolution grows with sum dt a and sum dt b. The integrated rule adds the evaluation of its interval integrals, at most 129 units of planning work per step for the cubic schedule and 2 otherwise.

The Schrodinger stencil construction counts W_stencil=28S+4M+(37d+12)D+40(d+1) logical visits (theory.restricted_kinetic_work), where D=K**d, S counts raw COO entries and M counts final CSR entries. The count covers the raw fill, index and mask passes, qualified sorted COO-to-CSR conversion, duplicate reduction and possible compaction. The 40(d+1) bundle covers scalar construction bookkeeping and five conversion endpoints. Planning uses the boundary-free upper counts S=3dD and M=(1+2d)D, giving (129d+16)D+40(d+1). These units count named element and record visits, not elapsed time. The docstring of theory.restricted_kinetic_work and the engineering constants give the itemized count and its dependency assumptions.

The classical kernel's limit check adds 1024 bytes for each of its d K + 3 d + 13 returned scalar records (method._SCALAR_BYTES), and 262,144 bytes once per call for Python objects of bounded total size that no size formula counts, such as the fixed objects of the readout and of the returned records (method._KERNEL_CALL_BYTES). Each allowance is set above sizes measured with tracemalloc, about 400 to 500 bytes per step, 2,400 to 3,500 per block, 600 per scalar and at most 156,775 bytes per kernel call beyond the arrays that the kernel's counts include (measured evidence), which depend on the Python runtime (engineering constants). Each classical kernel's work count adds the work of forming the start vector from the stored amplitudes once per evolution, D units for the uniform fill and d K + sum_(j=2)^d K**j + D for a product of the stored factors (initial_state.start_vector_work).

theory_flavor="split_step" counts W_start + D (d + M + 7) + d K + N_s (2 d L D + (2 d + 13) D + 7 d K + d + 4) units for D = K**d, the start vector's construction W_start (initial_state.start_vector_work), M support tables, N_s steps and the transform level count L = ceil(log2 K) for the FFT or ceil(log2(2 (K + 1))) for DST-I, and 64 D + 32 (d + 1) K + 8 T_max + 8 N_s + 1024 K + 4096 (d + 1) bytes, where T_max is the largest table and the last two terms are measured allowances for SciPy's transform scratch and Python objects (split_step.sizes). The transform term 2 d L D is a nominal proxy, the level count times the entries transformed, not a count of the backend's operations. The terms 7 D once, (d + 7) D + 4 d K + d + 4 per step and 8 N_s bytes count the reductions and per-step error terms of the observed state budget (Numerical guarantees), whose N_s steps take N_s + 1 potential moments. Apart from the transform term, each per-step term counts one unit per element of an array pass that the kernel runs, for example four passes of K per axis to weight a kinetic moment and two for a kinetic phase, while the once-per-evolution terms are nominal counts. Proposition 41 derives this work count.

Binary circuit construction

Binary circuit construction counts its declared scalar and array-element arithmetic, including phase wrapping and omitted-rotation error arithmetic. The engineering constants give the complete count and define its unit. These are units of planning work, not measured CPU instructions or elapsed time.

Compact binary planning checks its compilation work against max_work before BinaryModel or compile_binary_steps runs. The count includes kinetic setup, identity enclosures, each potential table's absolute mean and identity radius, Walsh normalization losses, requested synthesis trials, omission-error arithmetic and final numerical reductions. A requested dense block with N entries costs at most 2N+7 units. A Walsh-capable request costs at most 7N+80 units, including a min_cx trial that ultimately selects dense synthesis. Quantum binary planning, except with initial_state_preparation="none", also checks the bytes of the whole circuit-construction phase, the source records, the binary model and the circuit graph with one block's construction workspace, as construct_qhd checks it before allocating (native._binary_native_bytes, QHD._select_binary_native).

Classical binary ir_product planning also checks the run's model phase against QHD.max_bytes. It sets aside bytes for the Plan's source records, 48 bytes per grid state for the live state buffers, 16(sum_potential E_t + K + E_max) bytes for phase buffers, and the binary model baseline, where E_t is a potential table's size and E_max=max(K,E_1,...,E_T). Other planning phases can require a larger limit. These amounts use the qualified object sizes in the engineering constants.

Rotation-count checks

For a binary Plan, circuit_resources first checks the classification of its rotations against the Plan's QHD.max_work and QHD.max_bytes and raises ValueError when the classification exceeds them. The classification synthesizes each distinct stored block once and groups its angle magnitudes with np.unique. Its work and bytes cover the grouping of the stored blocks, the angle formation of each dense diagonal by its full phase-table length even when no rotation survives, the grouping arrays, the Python magnitude and count objects, the final sort and the validation of the classified rotations (resources._census_sizes). The classification shares the byte limit with the binary model, its construction cache and the Plan's source records. Quantum binary planning checks the same classification against the Plan's own limits after compilation releases its model, with the source tables counted once, so a quantum binary Plan that passed planning is inspected under the same limits. A one-hot Plan whose projector blocks use the diagonal provider checks the bytes of its classification beside its source records, the classified rotations and one wrapped-phase cache with its largest miss (resources._admitted_onehot_population), at planning and at inspection. run_resources checks each inner binary Plan's classification the same way, with the kept states of the run's inner Results and the tables of its other inner Plan objects held, so it can refuse a Plan that passed planning when those kept Results do not fit beside the classification.

Constrained and refinement runs

For an unrefined run with R = max_iterations, W = max_work and B = max_bytes, the recorded sums of limits are (4 R + t R + 1) W units of planning work and (4 R + t R) B array bytes. Here t = 6(m+1)-4 for quantum auto without a user GaussianState, with m the number of kept inequalities, and t = 0 otherwise. The layer's preprocessing and point evaluations share the final W limit, and the extra t term covers bounded selection trials. These sums of per-category limits are not runtime or peak-memory bounds.

For a run with refinement, with R = max_iterations and L = max_levels, W = max_work, B = max_bytes, a = 6 for search-model levels and a = 4 for physical levels, the sums of limits are (a R L + t R + 1) W units of planning work and (a R L + t R) B array bytes. For quantum auto without a user GaussianState, t = 6(m+1) with m kept inequalities. Otherwise t = 0. These are sums of per-category limits over attempts, not counts of limit checks or a peak-memory bound. The refinement's geometry, marginal and joint-mass processing, its own decompositions and evaluations of the inner objective at its level points, and the JSON from which the records' content hashes are computed lie outside these bounds. The number of JSON values of the run record grows in proportion to the accepted work, at most 14 d + 77 per completed level of d inner variables, including slack variables, and 8 d + 13 more for a level that made a stall split, while its bytes have no separate budget (engineering constants).

A kept state's dense probability grid is stored for a refinement level only when the level's QHD.max_bytes and QHD.max_work allow it. Otherwise each use streams the kept state in chunks of 4096 grid points and forms the same probabilities without the grid.

Measured evidence

The measurements that this guide quotes ran with Python 3.12.14, NumPy 2.5.2, SciPy 1.18.1, SymPy 1.14.0, Qiskit 2.5.2 and Aer 0.17.2 on macOS arm64, unless the table states otherwise.

Measurement and section	Settings that the text does not state
First- and second-order products on an 8-by-8 periodic grid (Scientific model and ordering)	Box [0, 1.6) x [-0.7, 1.1), objective 3(x - 0.2)² + 2(y + 0.1)² + 0.7 x y + 0.2 cos(3x), the uniform initial state, which is the kinetic ground state of the periodic grid, the exact QFT and no pruning. The fidelity is |<psi_ref|psi>|² of normalized states, with theory_flavor="schrodinger" at the same step weights as the reference.
Convergence of the two coefficient rules (Schedules and coefficient rule, Split-step classical evaluation)	Quadratic schedule with gamma 0.3, step counts 16 to 256. The error is the phase-aligned 2-norm distance from a DOP853 solution with relative tolerance 2e-12.
Planning work of theory_flavor="schrodinger" and theory_flavor="split_step" (Split-step classical evaluation)	Planning only, with max_work = 1e20 and max_bytes = 1e12 so that planning reports the count instead of refusing it, and the quadratic schedule with gamma 0.3 or ShiftedCubicSchedule(s=2e-4) where the text names it.
theory_flavor="schrodinger" and theory_flavor="split_step" at a common state error (Split-step classical evaluation)	Rho 1, 8, 64 and 512. The reference is the split-step product with 16,384 steps, within 2.1e-7 of 8,192 steps, and within 5.4e-8 of a DOP853 solution on four checked cases. A step count qualifies when its error plus that reference difference is below the target, and each ratio compares the least wall times among the qualifying step counts, from one timing per run that includes planning. The references of ShiftedCubicSchedule(s=2e-4) use 49,152 steps at K = 64 and 65,536 at K = 32.
Augmented-Lagrangian runs on constrained Rastrigin with theory_flavor="split_step" (Split-step classical evaluation)	Constrained Rastrigin of Wu et al., arXiv:2605.12066v1, Eqs. (17)–(18), the most_probable point, refinement of at most three levels per round with gain 8, max_no_improve=2 and mass threshold 0.99, penalty from 1, at most 15 rounds, normalized tolerance 1e-6 and max_work = max_bytes = 10**10.
Defaults of the initial state, the point rule and the stopping test (Initial states and preparation, Constrained problems)	theory_flavor="schrodinger" with T = 10 and 200 steps, and theory_flavor="split_step" with 3,200 steps for constrained Rastrigin and 320 for the equality problem of Proposition 51, both with the quadratic schedule (gamma 0.3), midpoint weights and max_work = max_bytes = 10**10 per solve. Standalone refinement ran at most ten levels, with max_no_improve=10 for the initial-state comparison and 2 for the point-rule comparison. Each augmented-Lagrangian round refined at most three levels with gain 8, max_no_improve=2 and mass threshold 0.99, from penalty 1, for at most 15 rounds at normalized tolerance 1e-6. Sampled runs drew 256 outcomes per level from each computed distribution with seeds 7, 19 and 43, by composing one-level classical solves. The stopping comparison continued each run to a stop on the penalty measure of Eq. (4.9) of Birgin and Martinez. Its five exact problems were constrained Rastrigin, the unit disk (K = 16), the off-grid equality x + y = 0.3 (K = 7), the mixed-scale problem (K = 7) and that equality problem, and its six sampled runs that equality problem and constrained Rastrigin with the three seeds.
Feasibility tolerance (Constrained problems)	Classical execution on the Dirichlet interior grid with KineticGroundState and QuadraticSchedule(gamma=0.3), refinement in every round with the search model, gain 8, at most three levels and mass threshold 0.99, T = 10, theory_flavor="schrodinger" at 200 steps and theory_flavor="split_step" at 3,200. The unit disk used K = 16 on [-1, 1]² and stopped at round 7, and constrained Rastrigin K = 32 on [-5, 5]² and round 10. The mode_or_mean disk run used that rule as both the inner point and the refinement point rule, midpoint coefficients, seed 11, max_iterations=15 and equal feasibility and complementarity tolerances, under both theory_flavor values. Its stop under 1e-6 was read from the rounds of the run under 1e-9, because the tolerances enter only the stopping test, and a direct Schrodinger run under 1e-6 stopped after the same 11 rounds with the same violation.
Penalty limits of classical planning (Constrained problems)	The largest projector angle is that of quantum planning at rho = 1024.
Slack and PHR planning with eight and fourteen variables (Slack variables for inequality constraints)	Planning only, through plan_augmented_lagrangian with execution="quantum", once with each of inequality_form="phr", "slack" and "auto", under max_work = 100_000_000 and max_bytes = 10_000_000_000.
Comparison of the three inequality forms (Slack variables for inequality constraints)	Each case named here minimizes f = sum_i (x_i - 1/2)**2 on [0, 1]**n under one inequality g <= 0, and each ran with PHR, forced slack and auto. Twelve of the twenty classical cases use the affine inequality g = sum_i x_i - n/4 at n = 3, 4 and 6 with K = 4 and at n = 3 with K = 8, each on the one-hot Dirichlet interior grid, the one-hot Dirichlet endpoint grid and the binary periodic grid. Two cases at n = 3 on the binary periodic grid with K = 4 use the squared-coordinate inequality g = sum_i x_i**2 - 3/8 and the coupled inequality g = x_0 x_1 + x_1 x_2 + x_0/4 - 1/4. The other six use the affine inequality at n = 3 and K = 4 on the endpoint and periodic grids, with inner_point="best_observed", with inner_point="mode_or_mean" or with three-level refinement. The three quantum cases are the affine cases at n = 3 and K = 4 on the three grids. Classical cases used theory_flavor="split_step" and quantum cases the circuit, both with the finite-difference kinetic, KineticGroundState, QuadraticSchedule(gamma=0.3), midpoint coefficients, total time 6, 64 steps, exact readout, seed 7, and max_work = 200_000_000 and max_bytes = 268_435_456 per QHD check. The outer layer ran at most eight rounds from zero multipliers and penalty 1, with growth 2, reduction ratio 0.25, penalty_update="on_insufficient_decrease", feasibility and complementarity tolerances 1e-6, the default max_penalty of 1e9, absorb_bounds=False, stationarity=True, unit scales and no multiplier safeguard. Refinement used the search model, gain 8, three levels, max_no_improve=2, mass threshold 0.99, the most_probable point and no stall split. The quantum cases used the noiseless AerBackend with shots=None, the simulator-state limit ExecutionLimits(simulator_memory_mb=256) and the default 20-qubit simulation cap. The round-0 planning cases with eight and fourteen variables used g = sum_i x_i - 3 on the three grids with K = 4, one step, total time 0.001, max_work = 100_000_000 and the same max_bytes, and gave the same binary periodic values under the limits of the row above. The trajectories, CX totals, elapsed times, objective gaps and memory figures come from executed runs.
Table II counts of Liu et al. (Circuit synthesis)	Planning only, with first-order steps and midpoint weights.
Windows evaluated from their formulas (Readout)	The 3-by-3 grid with 80 steps gives 9.0e-12 at T = 1 and 1.3e-11 at T = 10. The Aer windows are those of the formula for 36, 98 and 190 circuit instructions. Each is a first-order bound, not a measured error.
Split-step windows and budgets (State budgets, Window sizes)	Seed 7, and seed 1 for the K = 16 example. The comparisons with 30-, 40- and 50-digit evaluations kept the state (keep_state=True) and used mpmath 1.3.0. The five-point example uses f(x)=x**3/8-x**2/2+x/4+3/8 on the Dirichlet interior grid of [-1.5,1.5], spacing 0.5, three integrated steps to T=0.75, ShiftedCubicSchedule(s=1e-3), the kinetic ground state, keep_state=True and seed 7. The two eight-point periodic spectral examples use (x-1/3)**2, [-1,1), three integrated steps and GaussianState(center=(0.2,), widths=(0.3,)). The s=1e-150, T=1 case refuses and the gamma=0, T=1e150 case solves. These are Plan budget calculations, distinct from the five-point trajectory budget and measured state discrepancy.
Phase-error budgets of the running phase total (Numerical guarantees)	The Ackley replication is the two-dimensional Ackley function on [-5, 5]² mapped to the unit square, on the one-hot Dirichlet interior grid with T = 10, theory_flavor="ir_product" with keep_state=True, seed 7 and max_work = 1e18. The c = 1e17 case is x² + c with K = 3, QuadraticSchedule(gamma=0.0), two steps and T = 2.
Planning byte counts (Planning limits)	Initial-state calibration used 64-bit CPython 3.12.14 and NumPy 2.5.2 on macOS arm64, with Gaussian states at d = 1 to 3 and K = 2 to 4096 and the kinetic ground state. The 336dK+444d+104+4096(d+1) byte count covers one-variable grid entries. Peaks traced with tracemalloc after one warm-up call, with Pydantic 2.13.5. The per-call size was traced with the kernel that reads its generator norms and compiled blocks one at a time, and a second trace of the same cases without the one-hot IR product at K >= 32 gave at most 155,938 bytes. The cases are listed in engineering constants.
Cost of the Clifford distance bound (Fault-tolerant resources)	Apple M3 Max, best of three runs, on 20,000 distances and on the 155,063 rotations of the one-hot Ackley circuit of Liu et al. at K = 64 under CubicSchedule(s=1), with E = 1e-4.
NWQEC compilations (Fault-tolerant resources)	NWQEC 0.1.2.
Tightness of the evolution and angle-formation bounds (Fault-tolerant resources)	The six-point one-variable comparison used f(x)=x**2/4+x**3/16+7/8, four steps, total time 0.5, Dirichlet [-3.5,3.5] or periodic [-2.5,3.5), and a DOP853 reference with relative tolerance 2e-13. The two-variable binary comparison used four points per axis on [-1,1)**2, four steps, total time 0.25, f(x,y)=x**2/4+y**2/8+x*y/3+7/8, gamma=0.3 and DOP853 relative tolerance 1e-12. Its recorded ratios were 2.10–5.06 for finite differences and 5.38–11.92 for spectral kinetic energy. The seven angle-formation cases plan x**2/3 - x*y/5 + y/7 + 1/9 on the binary periodic grid of [0, 1.7) x [-1, 1.2) with three steps to total time 0.45, QuadraticSchedule(gamma=0.2) and the uniform initial state, at (K, kinetic model, synthesis of the potential and the kinetic phase, order) = (4, finite difference, Walsh, 2), (4, spectral, Walsh, 1), (4, finite difference, dense, 2), (4, spectral, dense, 1), (2, finite difference, Walsh, 2), (8, finite difference, Walsh, 1) and (8, finite difference, dense, 2), and compare with exact rational and 60-digit evaluations.
Platform assumptions of the error sources (Fault-tolerant resources)	Full test suites on macOS arm64 and Linux aarch64. The Linux host had the same Python and package versions. The suites include the angle-formation and AQFT tests against exact evaluations.
Spatial models of the binary circuit (Explicit verification)	T = 1, the quadratic schedule with gamma 0.3 and midpoint weights. The references are DOP853 solutions of the same models with relative tolerance 2e-12.
Refinement examples (Box refinement)	Seed 7, seed 3 for the example with exact quantum readout and seed 11 with at most seven levels for the added-constant example. The gain comparison of the default used max_work = max_bytes = 10**10, and the sampled Ackley runs drew 256 outcomes per level from each computed distribution with seeds 7, 19 and 43.

The library's tests check QHD against independent references:

• 3- and 4-qubit amplitudes against independent finite differences, with the meaning of coordinates and phases.
• On the periodic grid, the classical model against a 30-digit eigendecomposition of an independently built circulant Hamiltonian, the emitted hopping blocks against the circulant stencil, exact 4-qubit circuit probabilities against the classical model as the step count doubles, and the grid reference of an augmented-Lagrangian run against an enumeration of the periodic points.
• The split-step kernel against a 30-digit Strang product whose kinetic exponentials come from eigendecompositions of independently built stencils, the DST-I and Fourier eigenbases and the spectral energies, second-order time convergence against an adaptive ODE solution, exactness of the integrated rule for commuting terms, convergence to theory_flavor="schrodinger", and SciPy's transforms against 80-digit sums.
• The full 4-qubit binary evolution operators, global phase included, against independently written products of both kinetic models for every synthesis choice, bit-reversal choice and product order, and the 6-qubit operators of both models with the default choices. Further binary checks cover the signed square, the CX counts against transpiled circuits, the AQFT bound against the measured truncation error at b = 6, and circuit amplitudes, circuit probabilities and the classical binary product against the independent state.
• Time-step refinement and sampled probabilities on 3- and 5-qubit circuits, and compact provider circuits of at most eight qubits.

Source map

Rows marked standard or NWQLib have no paper location. The docstring of the named code writes out the relation.

Scientific step	Source and location	Code
Time-dependent Hamiltonian exp(phi_t)(-Delta/2) + exp(chi_t) f(x)	Leng et al., arXiv:2303.01471v1, Eq. (1), p. 3	compiler.QHDCompiler
Quadratic schedule a = 1/(1+gamma*t²), b = 1+gamma*t²	Kushnir, Leng, Peng, Fan and Wu, QHDOPT, arXiv:2409.03121v1, Sec. 2.1, the form it recommends for general nonconvex problems. The default schedule, with the illustrative gamma 0.3	schedules.QuadraticSchedule
Cubic schedule a = 2/(s+t³), b = 2t³	Leng et al., arXiv:2303.01471v1, Eq. (C.4), p. 32. Liu et al., arXiv:2607.16996v1, Eq. (92), write it with s = 1 as QHD-C	schedules.CubicSchedule
Shifted cubic schedule a = (2/(s+t))³, b = 2t³	Wu et al., arXiv:2605.12066v1, Eq. (15). Offered for replication with its own name, because its dynamics differ from those of Leng et al.'s Eq. (C.4), which Wu et al. cite as its source. Liu et al.'s code runs this form as qhd-c for the QHD-C schedule of their Eq. (92), arXiv:2607.16996v1, and with s = 1 only this form reproduces the Rz counts of their Table II	schedules.ShiftedCubicSchedule
Midpoint and integrated coefficient rules	NWQLib choice. The integrated step exponents are the first Magnus term of H(t). Leng et al., arXiv:2303.01471v1, Eq. (E.3) and Algorithm 1 use the left endpoint	schedules.step_weights, compiler.QHDCompiler
Interval integrals without cancellation, 16-point Gauss–Legendre panels for cubic A	NWQLib derivation in the code's docstrings	schedules.QuadraticSchedule.kinetic_integral, schedules.CubicSchedule.kinetic_integral
Central-difference kinetic stencil and diagonal grid potential	Leng et al., arXiv:2303.01471v1, Eqs. (F.7) and (F.9), Remark 6, pp. 46–47, and Eq. (F.15), p. 47, for the potential of d variables	kinetic.KineticCompiler, theory.restricted_kinetic_sparse
Vanishing boundary values and the endpoint mesh	Leng et al., arXiv:2303.01471v1, text after Eq. (F.4) and Eq. (F.5), p. 46	grid.OneHotGrid
Periodic grid, wrap link and circulant kinetic matrix	Liu et al., arXiv:2607.16996v1, Sec. III, Eq. (12), whose one-hot hopping includes the wrap term X_(N-1) X_0 + Y_(N-1) Y_0, and the grouping of its hopping strings into two commuting layers in the same section. The endpoint-exclusive coordinates are the mesh of Leng et al., arXiv:2303.01471v1, Eq. (E.1), p. 43, and Algorithm 1 step 1, p. 44, which assumes a periodic boundary condition on the box (p. 44). The even-K condition is NWQLib's, derived in the code's docstrings	grid.OneHotGrid, kinetic.KineticCompiler, theory.restricted_kinetic_sparse, method.QHD
Sum over variables and the constant stencil diagonal as a phase	Leng et al., arXiv:2303.01471v1, sentence introducing Eq. (F.14), p. 47	compiler.QHDCompiler
One-hot occupation operators and the single-variable potential	Wu et al., arXiv:2605.12066v1, Sec. IV.A, Eqs. (9)–(10), p. 5	potential.PotentialCompiler
Products of single-variable operators for product terms	Wu et al., arXiv:2605.12066v1, Sec. IV.A, text after Eq. (10), p. 5	potential.PotentialCompiler
General support term tabulated over its grid tuples	NWQLib extension of the row above	objective.ObjectiveDecomposer, potential.PotentialCompiler
Expansion about the box point nearest the origin	NWQLib choice against cancellation, explained in the expansion_centers docstring	objective.expansion_centers, compiler.QHDCompiler._objective_grid_values
XX+YY hopping between neighboring one-hot states	Standard two-qubit identity	kinetic.KineticCompiler
First-order step with potential applied before kinetic	Leng et al., arXiv:2303.01471v1, Eq. (E.3), p. 43, and Algorithm 1 step 5, p. 44	compiler.QHDCompiler.build_step_pauli_ir
Default second-order step and odd/even link layers	Standard symmetric (Strang) splitting. Leng et al., arXiv:2303.01471v1, Algorithm 1, use a first-order step. Liu et al., arXiv:2607.16996v1, Sec. V A, build first- or second-order steps and split the one-hot hopping into even and odd link layers, as in their Eq. (13)	compiler.QHDCompiler.build_step_pauli_ir, kinetic.KineticCompiler
Uniform initial superposition	Leng et al., arXiv:2303.01471v1, Algorithm 1 step 3, p. 44	initial_state.UniformState
Kinetic ground state sin(pi (i+1)/(K+1)) on both Dirichlet grids	Standard eigenvectors of the tridiagonal stencil, derived in the code's docstring	initial_state.KineticGroundState
Gaussian warm start, its shifted exponents and its error bound	Leng et al., arXiv:2303.01471v1, Algorithm 1 step 3, p. 44, name a Gaussian state as one possible start. NWQLib record for the Gaussian start of the Wu et al. code (arXiv:2605.12066v1), with the derivation in the code's docstrings	initial_state.GaussianState, initial_state._gaussian_beta
Structured one-hot preparation: the amplitude chain, its zero-tail stop, its lower-range cutoff and its 3L CX count	Standard linear W-state construction, generalized to nonnegative amplitudes in the code's docstring. The cutoff and its error are an NWQLib derivation in the code's docstring	initial_state.append_amplitude_chain, initial_state.chain_selection, initial_state.chain_links
Classical start vector and its construction error	NWQLib derivation in the code's docstrings	initial_state.restricted_state, initial_state.restricted_state_error
Number-projector phase and its CX count	Phase diagonal: Shende, Bullock and Markov, quant-ph/0406176v5, Theorem 7, p. 10, with 2**k CX for a multiplexed Rz of k select bits by the count after their Theorem 8, p. 11. Multi-controlled phase: Qiskit's MCPhaseGate definition, no paper location	nwqlib.subroutines.hamiltonian_evolution.pauli_evolution.structured_number_projector_provider
Arbitrary rotations and exact T gates of the emitted circuit	Qiskit 2.5.2 gate definitions, including synth_mcx_n_dirty_i15, and the phase-diagonal multiplexor, counted in the code's docstrings. No paper location	resources.rotation_population, resources.rotation_law
Clifford distance 2 sin(abs(delta)/4), budgeted replacement and even split of the synthesis budget	Standard identities, derived in the code's docstrings	resources.clifford_distance, resources.synthesis_projection
Leading-order T count 3 log2(1/eps) per synthesized rotation	Ross and Selinger, arXiv:1403.2975v3, typical T count of one Rz approximation	resources.T_PER_PRECISION_BIT, resources.synthesis_projection
Splitting bounds of the one-hot and binary first- and second-order products	Childs, Su, Tran, Wiebe and Zhu, arXiv:1912.08854v3, Sec. 5.1, Propositions 15–16, Eqs. (145) and (152), applied to the emitted groups in the code's docstring	evolution_bounds.evolution_bound
Time-ordering and midpoint-quadrature bounds of a step, with the schedules' derivative bounds	Duhamel's identity and the midpoint-rule remainder (Proposition 39), with the schedules' bounds in QHD schedule values, integrals and derivative bounds	evolution_bounds._schedule_terms, schedules.QuadraticSchedule.derivative_bounds
Commutator bounds from table ranges, link differences and the layer graph	Standard norm inequalities, derived in the code's docstrings	evolution_bounds._norm_inputs, evolution_bounds._graph_constants
Coefficient residual of the stored step exponents	Exponential perturbation inequality with exact midpoint and rational-integral discrepancies and the schedules' documented integral error, derived in the code's docstring	evolution_bounds._coefficient_term
Block-angle formation, omitted blocks, hopping parameters, circuit phase bookkeeping and structured preparation	NWQLib derivations in the code's docstrings, with Qiskit 2.5.2's rem_euclid phase setter	circuit_errors.block_errors, circuit_errors.phase_allowance, circuit_errors.preparation_error
Binary angle formation: Walsh coefficients, kinetic coefficients and energies, dense phases and QFT angles	Butterfly rounding: Higham, SIAM J. Sci. Comput. 14 (1993) 783–799, doi:10.1137/0914050, Sec. 3, Eq. (3.6), p. 788, and the pairwise summation bound of Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., doi:10.1137/1.9780898718027, Sec. 4.2, Eq. (4.6), p. 83. The enclosures and the stage chain are NWQLib derivations in the code's docstring	circuit_errors.binary_angle_formation, circuit_errors.binary_phase_allowance
Position measurement and valid/invalid register decoding	Leng et al., arXiv:2303.01471v1, Algorithm 1 step 6, p. 44. Invalid-mass accounting is NWQLib's	decoding.decode_onehot_basis_index, decoding.decode_basis_index, method._summarize
Binary register, its decoding and the variable-axis permutation to lexicographic order	NWQLib, derived in the code's docstrings	binary.decode_register_index, binary.lexicographic_register_indices, binary.lexicographic_register_array, binary.address_table
QFT circuit and its bit reversal	Nielsen and Chuang, Quantum Computation and Quantum Information (2000, ISBN 978-0-521-63503-5), Sec. 5.1 and Fig. 5.1, in the gate order of Qiskit's synth_qft_full	binary.qft_gates
Relabeled phase table between the swap-free QFT and its inverse	NWQLib derivation, W^dagger R D R W = F^dagger D F	binary.kinetic_table
Approximate QFT	Coppersmith, arXiv:quant-ph/0201067v1 (IBM report RC 19642). The operator-norm bound per QFT and its accumulation over conjugations are NWQLib's	binary.qft_error_bound, binary.compile_binary_steps
QFT sign of the binary kinetic factor	NWQLib derivation, E_(-k mod K) = E_k for both periodic energies	binary.kinetic_table
Quadratic low-momentum kinetic phase and its Pauli expansion	Liu et al., arXiv:2607.16996v1, Sec. IV E, Eqs. (89)–(91), where Eq. (90) writes the unsigned index. NWQLib uses the signed index, the low-momentum form of Eq. (89)	binary.signed_square_walsh
Pauli expansion of the finite-difference kinetic phase	NWQLib derivation	binary.kinetic_walsh
Walsh expansion of a support-local diagonal, one parity-ladder rotation per Z string	Liu et al., arXiv:2607.16996v1, Sec. IV A, Eqs. (53)–(54), with the rotation of Eqs. (57)–(58) and its CX count in the sentence after them. The CX count, the pruning rule and the per-table min_cx choice are NWQLib's	binary.walsh_coefficients, binary.PhaseTable.synthesize
Exact phase diagonal of a support table or kinetic phase, 2**n - 2 CX	Shende, Bullock and Markov, quant-ph/0406176v5, Theorem 7, p. 10, with 2**k CX for a multiplexed Rz of k select bits by the count after their Theorem 8, p. 11	binary.PhaseTable.dense_cx, nwqlib.subroutines._multiplexors.append_control_diagonal_phases
Classical binary product and its roundoff budget	NWQLib derivation in the code's docstrings	theory.run_binary_product, theory.binary_product_state_error
Most probable valid point and its tie window	NWQLib derivation in the docstrings of the named code	method._summarize, method._readout_window, split_step.evolve (observed split-step budget), nwqlib._validation.probability_difference_window, nwqlib._validation.native_state_error, nwqlib.execution.PreparedArtifact.state_error
Conditional position mean and standard deviation	Liu et al., arXiv:2607.16996v1, Eq. (94), as expectations in the final state. The conditioning on a valid outcome is NWQLib's. Computed on demand and reported, not used as an agreement or success criterion, because distributions with equal moments can differ in total variation by 1	records.QHDAnalysis.position_mean, records.QHDAnalysis.position_standard_deviation
Restricted Schrodinger reference	Standard exponential midpoint rule under the midpoint coefficient rule, and the first Magnus exponent of each step under the integrated rule	method._evolve_restricted
Objective constant as one global phase of theory_flavor="schrodinger" and theory_flavor="split_step"	NWQLib derivation in the code's docstring. c I commutes with every step operator, so it factors out exactly	method._constant_phase, method._admit_constant_phase, split_step.potential_diagonal
Phase-error budget of the compiled physical phase and the pi rule for kept circuit-route states	NWQLib derivation in the code's docstrings: the formation of each identity contribution, Neumaier's compensated summation (Neumaier, Z. Angew. Math. Mech. 54 (1974) 39–51, doi:10.1002/zamm.19740540106), whose recurrences are those of the cascaded summation of Ogita, Rump and Oishi (SIAM J. Sci. Comput. 26 (2005) 1955–1988, doi:10.1137/030601818, Algorithm 4.1). Their Proposition 4.5 bounds the equivalent Algorithm 4.4 by u |X| + gamma_(M-1)**2 sum_j |x_j|, and the budget uses the slightly larger residual bound u |X| + (1 + u) gamma_(M-1)**2 sum_j |x_j| that compiler.QHDCompiler._record_global_phase writes out. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., doi:10.1137/1.9780898718027, Sec. 4.3, Eq. (4.10), gives a backward error bound for the same class of summation. For the quantum circuit the derivation also counts each phase assignment's addition, reduction modulo the binary64 tau by Qiskit 2.5.2 and remainder rounding, 2u (Y + |physical_phase|) + 20u (B + 1). For the binary encoding, the formation budget F_W of the Walsh identity phases and the reconstruction budget R_W of the emitted phases	compiler.QHDCompiler._record_global_phase, method._compiled_phase_allowance, method._native_phase_allowance, method._admit_compiled_phase, binary.compile_binary_steps, binary.kinetic_identity_enclosure
Range rule of planning and upward evaluation of budgets	Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., doi:10.1137/1.9780898718027, Theorem 2.2 and Eq. (2.4) for the relative rounding of a result in the normal range, and Eq. (2.8) for gradual underflow. Evaluating each budget upward, one binary64 number above each rounded sum, product or quotient, is NWQLib's rule, stated in the code's docstrings	validation._normal_range, method._add_up
Lower-range omission of contributions formed from objective data and the initial state, and their errors	NWQLib derivation in the code's docstrings, from ||exp(-i a G) - I|| <= |a| ||G|| for Hermitian G, the projector split P = (P - 2**-s I) + 2**-s I and the distance of a shortened amplitude chain	validation._underflows, potential.PotentialCompiler.select_occurrences, binary.walsh_admission, binary.PhaseTable.synthesize, initial_state.chain_selection, records.QHDRangeOmissions
Classical-kernel state budget, mass window and tie window	Al-Mohy and Higham, SIAM J. Sci. Comput. 33 (2011) 488–511, doi:10.1137/100788860, Algorithm 3.2, Eqs. (3.5)–(3.15), Table 3.1, Lemma 4.1 and Eq. (4.7), with the NWQLib first-order derivation in the expm_multiply_call_roundoff and expm_multiply_state_error docstrings. The stored-table Schrodinger budget is the NWQLib derivation in the method._schrodinger_step_bounds and method._schrodinger_state_error docstrings, the split-step budget the one in the split_step.state_error docstring, and the direct one-hot product budget the one in the theory.onehot_product_state_error docstring	method._host_state_error, method._schrodinger_state_error, method._schrodinger_step_bounds, method._schrodinger_kinetic_bounds, theory.onehot_product_state_error, method._host_generator_norms, method._host_tie_window, initial_state.restricted_state_error, nwqlib._validation.expm_multiply_state_error, nwqlib._validation.state_mass_window, split_step.state_error
Split-step product with the kinetic factor in its eigenbasis	Leng et al., arXiv:2303.01471v1, Eq. (C.3), p. 32, the first-order pseudo-spectral step with left-endpoint weights. The symmetric placement is Strang's splitting, SIAM J. Numer. Anal. 5 (1968) 506–517, doi:10.1137/0705041, chosen because it is second order with the same number of transforms per step. The time-error terms and the equal potential halves of the integrated rule are NWQLib's derivation in the module docstring	split_step.evolve, method._evolve_restricted
DST-I and Fourier eigenbases of the stencils and the signed spectral index	Standard eigenvector identities, derived in the code's docstring. The orthonormal conventions are those of SciPy's scipy.fft.dst, fft and fftfreq	split_step.kinetic_eigenvalues, split_step._transform
Transform accuracy assumption and split-step work and bytes	NWQLib qualification constant for SciPy 1.18.1's ducc0.fft backend, checked against 80-digit sums, a nominal work proxy and a byte count with a scratch allowance measured as high-water-mark growth with one transform worker	split_step.TRANSFORM_ROUNDOFF, split_step.sizes
Fidelity reduction roundoff window	NWQLib derivation in the standard rounding model of Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., doi:10.1137/1.9780898718027, Secs. 2.2, 3.1, 3.4 and 3.6	verification._infidelity, nwqlib._numerics.normalized_fidelity_with_window
Hamming encoding, not used here	Leng et al., arXiv:2303.01471v1, Definition 7 and Eq. (F.36), p. 55	none
Augmented Lagrangian of equality constraints and the multiplier update	Wu et al., arXiv:2605.12066v1, Sec. III.B, Eqs. (6)–(8), applied to normalized residuals, so that rescaling a constraint together with its scale leaves every decision unchanged in exact arithmetic, and in binary64 rounding can change a decision near a tie or a threshold (Constrained problems)	constrained._effective_objective, constrained._update
PHR term of inequality constraints	Rockafellar, doi:10.1007/BF01580138, and Birgin and Martinez, doi:10.1137/1.9781611973365, Eq. (10.3), Sec. 10.1, p. 114. The default phr policy uses this term directly. Its constant -mu_bar_j**2/(2 rho) contributes a global phase and makes L_k equal the normalized objective at an exactly feasible point complementary to the entering multipliers. Eq. (4.3) differs by the equality and inequality square-completion constants. Converted inner slack terms have the partial-minimization relation in Proposition 54, and the outer update stays PHR at the projected point	constrained._effective_objective, constrained._effective_value, constrained._update
Explicit inner slack variables, caps and finite-grid excess	Proposition 54, including its support-table arithmetic qualification and periodic cap margin. The identity concerns partial minimization and does not equate joint QHD dynamics with PHR dynamics	constrained._inequality_range, constrained._slack_axis, constrained._branch, constrained._inner_objective
Inequality representation selection and projected readout	NWQLib policy. Eligible quantum automatic selection requires no larger planned classical work and CX count with one strict improvement, or an accepted trial with a CX count when current planning is refused. Inner point rules act on the joint problem and the outer update uses original functions at the projection	constrained._planned_costs, constrained._accepts, constrained._round_problem, constrained._choose, constrained._refined_choice
Initial-grid numerical finiteness check	NWQLib numerical rule. Finite support tables and a represented-range check, supplemented when required by complete supplied-summand scans and an addition bound, with a limit-checked whole-term fallback. Finite rounded values do not prove the exact real domain	constrained._setup, constrained._support_tables, constrained._domain_check, constrained._addition_certificate, constrained._admit_layer_work
Tentative multipliers, penalty test and multiplier safeguard	Birgin and Martinez, doi:10.1137/1.9781611973365, Algorithm 4.1 and Eqs. (4.7)–(4.9), pp. 33–34. The default penalty test replaces the growth after every round of Wu et al.'s code, whose text cites the update rules of Nocedal and Wright, doi:10.1007/978-0-387-40065-5. That growth stays available as every_iteration. With the same growth factor the test never gives a round a larger penalty than growth after every round. A larger penalty can increase the range of the augmented objective and, on the Schrodinger path, the planning work counted from its generator norm. The range need not grow monotonically because the objective, linear multiplier term and quadratic penalty can cancel. Clipping keeps the tentative multipliers whenever they lie inside the bounds, the choice of the book's Assumption 7.6 (Sec. 7.5, p. 64), and keeping rho when the test holds is the choice that its Theorem 5.2 (p. 42) and Theorem 7.2 with Assumption 7.9 (Sec. 7.7, p. 70) assume. Those convergence results also need a Step 1 point that satisfies the book's Assumption 5.1 (p. 41) or 6.1 (p. 48), which no point rule of the layer proves, and max_penalty and the default multiplier_bounds=None depart from Algorithm 4.1 as they assume it, so the layer claims none of them	constrained._update, constrained._next_penalty, constrained._safeguard
Stopping test and projected-stationarity diagnostic	Birgin and Martinez, doi:10.1137/1.9781611973365, Sec. 10.2.2, Eqs. (10.6)–(10.8), pp. 116–117, applied to normalized residuals and, for Eq. (10.6), in box coordinates. Eq. (10.6) is a diagnostic and not a stop, because a finite-grid minimizer generally has a positive projected stationarity	constrained.solve_augmented_lagrangian, constrained._stationarity
Multipliers in original units	KKT condition of the scaled problem, with the scaling of Birgin and Martinez, doi:10.1137/1.9781611973365, Eq. (10.4)	constrained.ConstrainedQHDResult.multipliers
Normalized marginals and the per-axis interval of mass eta	Wu et al., arXiv:2605.12066v1, Sec. V, Eqs. (12)–(13), p. 7. The paper breaks no ties, and the tie rules are NWQLib's	refinement._axis_interval
Next box from centered cells of the kept grid points, with faces at the midpoints of adjacent grid coordinates	NWQLib choice. Eq. (14) of the same section uses left-endpoint cells	refinement._next_box
Joint box mass lower bound max(0, 1 - sum_j (1 - m_j))	Standard union bound	refinement._joint_mass_bound
Simultaneous selected region coverage from Hoeffding and one-sided Clopper–Pearson bounds, each using half of alpha across all configured levels	Wu et al., arXiv:2605.12066, with the full derivation from Hoeffding, doi:10.1080/01621459.1963.10500830, and Clopper–Pearson, doi:10.1093/biomet/26.4.404, in Proposition 49	_coverage.interval_radius, _coverage.cp_lower, _coverage.bounds_from_counts, _coverage.coverage
Search model kappa (F(a + D u) - c)/E on the unit box, its rounding budget and the potential gain kappa	NWQLib definition, derived in the code's docstring	refinement._level_problem, _outer.range_bound
Flat level (zero table ranges) and unresolved level (table ranges within the units in the last place of the table values), a resolution rule	NWQLib rule, stated in the code's docstring	refinement._resolution_stop
Box width floor, the binary64 resolution of adjacent grid coordinates, and the grid's check of its spacing and kinetic coefficient	NWQLib derivations in the code's docstrings	refinement._resolved, grid.OneHotGrid.__post_init__
Stall condition (both end cells above 1 - eta) and the kept mass of a split, the valley ratio min(L, R)/p(v), its resolution screen (readout window, or Hoeffding's inequality for counts), the cut face and the choice of the lower-scoring side	NWQLib derivations in the code's docstrings	refinement._stall_split, refinement._clearest_valley, refinement._valley_admission
Point rules of both layers: the candidate, the most probable point, and the valid-mass mean when its objective is strictly smaller	Wu et al., arXiv:2605.12066v1, Sec. V, states the most probable point, the default because it keeps a point that the computed distribution defines. The other two rules are NWQLib definitions, derived in the code's docstrings	_outer.grid_point, _outer.mean_point

Code map

• method.py: the QHD Method fields, planning (support tables, compiled blocks and operation-size checks), the classical restricted evolution and the readout analysis.
• records.py: the stored planning data (support tables, compiled blocks, running phase total, pruning and AQFT bounds), the result record with its consistency checks, and the explicit verification options.
• compiler.py, objective.py, potential.py, kinetic.py, schedules.py: grouping of the objective by variable support, the grid tables evaluated once, the schedule records with their interval integrals and step weights, and the compiled product-formula blocks.
• binary.py: the binary register order, the QFT gates and their AQFT bound, the Fourier energies of both kinetic models, the Walsh and dense syntheses of the phase diagonals with their CX counts, and the binary step compiler.
• initial_state.py, native.py: the initial-state records with their classical start vector and its construction error, the circuit preparation recipes of both encodings, the Qiskit circuit of either encoding and the restored identity phase.
• grid.py, decoding.py, theory.py: coordinates, spacing, kinetic links and one-hot qubit order, decoding of full-register outcomes of either encoding, the restricted matrices of the classical references and the classical binary product.
• split_step.py: the split-step classical evaluation (theory_flavor="split_step"), its kinetic eigenvalues and transforms, its state budget and its work and byte counts.
• archive.py, verification.py: saving and loading, and the explicit grid and fidelity comparisons.
• constrained.py, constrained_records.py: the augmented-Lagrangian layer, support-table and original-summand preprocessing, inequality representation selection, slack ranges and error records, joint readout and projection, PHR outer updates, refinement of the inner problem, the original-grid reference, and constrained records and archives.
• constrained.plan_augmented_lagrangian: preprocessing, representation selection and an unrefined round-0 Plan without a Run, for resource planning.
• refinement.py, refinement_records.py: box refinement over repeated QHD Plan objects, its search model, box rule, stall split and stopping, and its options, level records, split record and result.
• resources.py, evolution_bounds.py, circuit_errors.py: the rotation count of the emitted circuit, the budgeted Clifford replacement and T estimate, the per-circuit error sources and the run totals, the splitting, schedule and coefficient bounds of the compiled product, and the block-angle, omitted-block, phase and preparation entries of the quantum circuit.
• _outer.py: the rules that the augmented-Lagrangian layer and box refinement share, namely their common arguments, the random streams of each inner Plan, the remainder of the cumulative limits, the counts read from an inner Run, the point that a point rule reads, and how a round's refinement ends the augmented-Lagrangian run.
• _durable.py: the run directory of both layers, namely its layout, the outer record, its pre-measurement inequality representation and its commit by write-then-rename, the directory lock, the backend of the inner Runs with its saved noise model, the reopening of an unfinished round's or level's Run on resume, and the counts of a Run whose preparation raised.