QHD¶

Reference for the QHD Method, its schedules and initial states, its readout and its fault-tolerant resource estimates. Import every name on this page from nwqlib.algorithms.qhd, for example from nwqlib.algorithms.qhd import QHD, circuit_resources.

The QHD guide explains the model, the readout and the measured evidence, and its source map gives the paper, version and equation behind each step. The augmented-Lagrangian layer for constrained problems and box refinement are on QHD constrained problems and box refinement. The QHD entry below has a complete example.

Task	Entry
Configure and run QHD on an `Optimization`	`QHD` with `solve` or `plan`
Choose the binary circuit's synthesis	`BinarySynthesis`
Choose the schedule a(t), b(t)	`QuadraticSchedule`, `CubicSchedule`, `ShiftedCubicSchedule`
Choose the initial state	`KineticGroundState`, `UniformState`, `GaussianState`
Read the best observed and most probable points	`QHDAnalysis`
Compare the result with a grid minimum or a reference evolution	`QHDVerification` with `result.verify(checks=...)`
Count rotations and T gates and bound the product-formula error	`circuit_resources`, `evolution_bound`, `run_resources`

Configure QHD¶

QHD ¶

Bases: Method

Quantum Hamiltonian Descent (QHD) method for the minimum of an objective over a box.

Build it with keyword arguments and pass it as method= with an Optimization, for example solve(problem, method=QHD(num_grid_points=4)). Every argument is optional. The result is a QHDAnalysis. Its candidate is the best observed grid point and objective the objective there, in the Problem's coordinates and unit, and most_probable_coordinates is the grid point where the evolution concentrated probability. QHD is not a global optimizer. Both points belong to a finite grid of the box, and no setting implies a continuous global-optimization guarantee.

The Method evolves H(t) = a(t) K + b(t) V, the Hamiltonian of Leng et al., arXiv:2303.01471v1, Eq. (1), with a = exp(phi_t) and b = exp(chi_t), K the kinetic operator -Delta/2 on the grid and V the objective values, one step per time step, and reads out grid points. Each step is a product formula on the circuit routes and in the ir_product and split_step flavors, and a numerical evaluation of the exponential of the step's Hamiltonian in the default schrodinger flavor. The default quadratic schedule is the example of QHDOPT (Kushnir et al., arXiv:2409.03121v1, Sec. 2.1), and the one-hot encoding and the shifted-cubic schedule come from Wu et al., arXiv:2605.12066v1. The rows and notes below say where a setting departs from Leng et al.'s Algorithm 1. Its step 6 measures one grid point, while QHD reports the best observed and the most probable point of the observed population. The QHD guide explains the grids, schedules, initial states and readout, and its source map gives the paper location of each step or marks it as standard or NWQLib's own.

max_work and max_bytes limit the work and bytes that planning counts. They do not control undocumented SymPy, SciPy or Qiskit costs.

Attributes:

num_grid_points (int) –

Default 2. Grid points K per variable, at least 2. The one-hot periodic grid needs an even K of at least 4, and the binary encoding needs K = 2**b. K is also the one-hot register width per variable, and the binary encoding uses b qubits per variable.
encoding (Literal['one_hot', 'binary']) –

Default "one_hot", the established circuit. How each variable's grid index is stored in qubits: "one_hot", K qubits with one excitation, or "binary", b qubits for K = 2**b, which requires boundary="periodic". The Encodings note below describes both registers.
binary_synthesis (BinarySynthesis) –

Default BinarySynthesis(). Circuit choices of the binary encoding: the potential and kinetic phase diagonals, the QFT bit reversal and the AQFT cutoff (BinarySynthesis). The one-hot encoding requires the default record.
num_steps (PositiveInt) –

Default 1. Positive number of evolution steps.
total_time (Real) –

Default 1.0. Positive total evolution time. Each step lasts total_time/num_steps.
schedule (Schedule) –

Default QuadraticSchedule(), with gamma 0.3. Kinetic weight a(t) and potential weight b(t) of H(t) = a(t) K + b(t) V, one of QuadraticSchedule, CubicSchedule or ShiftedCubicSchedule. Its kind names the formula. A schedule parameter stays fixed when num_steps changes, and no schedule carries a convergence guarantee.
coefficient_rule (Literal['midpoint', 'integrated']) –

Default "midpoint". The weights each step uses: "midpoint", the point values a(t) and b(t) at the step midpoint, or "integrated", the step averages of their interval integrals. The Product formula note below explains both.
trotter_order (Literal[1, 2]) –

Default 2. Product-formula order, 1 or 2, shared by circuit selection and the numerical ir_product reference. The split_step flavor is always symmetric. The Product formula note below gives the factor order.
boundary (Literal['dirichlet', 'periodic']) –

Default "dirichlet", which accepts the default K = 2 with the one-hot encoding. Boundary condition of the finite-difference kinetic term, "dirichlet" or "periodic". The periodic grid takes no boundary points, and with the one-hot encoding it needs an even K of at least 4. The Grids note below gives the grid points.
include_boundary_points (StrictBool) –

Default False. Place the Dirichlet grid on both box endpoints instead of the interior. Not available with the periodic boundary.
kinetic_model (Literal['finite_difference', 'spectral']) –

Default "finite_difference", which every route applies. Kinetic operator K of each variable: "finite_difference", the three-point stencil of -Delta/2, or "spectral", the Fourier operator of -Delta/2. The spectral model needs the periodic boundary and a route that applies Fourier energies. The Kinetic models note below gives their eigenvalues and routes.
rotation_threshold (Nonnegative) –

Default 0.0, which keeps every computed nonzero angle. Angle in radians, nonnegative, below which compiled rotations are omitted. The initial-state preparation is never pruned. The Rotation threshold note below says which rotations each encoding omits and how the Plan records them.
initial_state (InitialState) –

Default KineticGroundState(). State that every route starts from, a product of nonnegative per-variable amplitude vectors: KineticGroundState(), UniformState() or GaussianState(center=..., widths=...). Leng et al.'s Algorithm 1, step 3, names the uniform and a Gaussian state. The row "QHD Method defaults" of the engineering constants gives the evidence for the default and its limits.
initial_state_preparation (Literal['structured', 'qiskit_state_preparation', 'none']) –

Default "structured", the library's circuit for the encoding and state, whose CX count the Plan records. How quantum execution prepares the initial state. "qiskit_state_preparation" appends a Qiskit StatePreparation, and "none" selects a resource-only construction, which classical execution rejects. The State preparation note below describes each.
theory_flavor (Literal['schrodinger', 'ir_product', 'split_step']) –

Default "schrodinger". Classical evolution of the K**d grid amplitudes under execution="classical": "schrodinger", the exponential of each step's Hamiltonian, "ir_product", a numerical reference of the compiled product, or "split_step", a symmetric split-step product. The Classical flavors note below compares them.
keep_state (StrictBool) –

Default False. Keep the final state array when the output and execution support it. It requires exact readout, without shots.
max_bytes (PositiveInt) –

Default 10 GB (decimal, 10_000_000_000 bytes). Positive limit on the known bytes of state, grid and numerical workspace that planning counts. It does not bound the process memory.
max_work (PositiveInt) –

Default 1_000_000_000. Positive limit on the counted work of objective tables, products and classical evolution, in planning work units, the scalar and array operations that planning declares. A work unit is not a measured CPU operation or a time. A refusal names the stage and the amount to set (cost).

Encodings

"one_hot" places variable j on the K qubits j K to j K + K - 1 with one excitation, the one-hot encoding of Wu et al., arXiv:2605.12066v1, Sec. IV.A, with the occupation operators of their Eqs. (9)-(10). The one-hot encoding is the default as the established circuit.

"binary" places it on the b qubits j b to j b + b - 1, grid index n_j = sum_l 2**l n_(j,l), so every outcome is a grid point. This is the register of Leng et al.'s Algorithm 1, and it applies the kinetic term as F^dagger exp(-i alpha diag(E)) F with a QFT on each register, the form of their kinetic step, Eq. (E.4). It requires boundary="periodic", whose kinetic operator the QFT diagonalizes, and K = 2**b, the register's state count.

Grids

The Dirichlet interior grid, the default, has the vanishing boundary values of Leng et al.'s Eq. (F.4), and the endpoint grid of include_boundary_points=True keeps the endpoint amplitudes, as their Eq. (F.9) does. The periodic grid, the mesh of their Eq. (E.1), is x_i = lower + i (upper - lower)/K, i = 0..K-1, with upper identified with lower, and adds the wrap link between the points K-1 and 0. It takes no boundary points, and with the one-hot encoding it needs an even num_grid_points of at least 4. Dirichlet is the default because it accepts the default K = 2 with the one-hot encoding.

Product formula

coefficient_rule="midpoint" uses the point values a(t) and b(t) at each step midpoint. "integrated" uses the step averages of the interval integrals, so each step's kinetic and potential exponents equal A and B over the step, the first Magnus term. Leng et al.'s Eq. (E.3) and Algorithm 1 take the left endpoint, which leaves the product first order in the step even with symmetric placement, so no left-endpoint rule is offered (Split-step classical evaluation).

trotter_order=1 uses the potential-then-kinetic ordering of Leng et al.'s Algorithm 1. Order 2 uses the symmetric composition of Strang, SIAM J. Numer. Anal. 5 (1968) 506-517, doi:10.1137/0705041, with a potential half on either side of the kinetic factor. The one-hot kinetic factor is itself approximated by a product over links, in increasing link order for order 1 and by an odd-half, even-full, odd-half split for order 2. Binary kinetic factors use Fourier conjugation of their phase diagonal.

Kinetic models

"finite_difference" is the three-point stencil of -Delta/2 from the central difference of Leng et al.'s Eqs. (F.7) and (F.9), with diagonal 1/h**2 and -1/(2 h**2) on each link. Its eigenvalues are (2/h**2) sin(pi r/(2 (K + 1)))**2, r = 1..K, on both Dirichlet grids (each with its own h), and (2/h**2) sin(pi k/K)**2, k = 0..K-1, on the periodic grid.

"spectral" is the Fourier operator of -Delta/2 that Leng et al.'s Eq. (E.4) and Algorithm 1 apply, on the periodic grid of period L = K h. Its eigenvalue is 2 pi**2 q**2/L**2 for the Fourier mode k whose signed index q is k below ceil(K/2) and k - K from there on. The two models agree at low momentum, 0 <= E_sp - E_FD <= 2 pi**4 q**4 h**2/(3 L**4), and differ by the factor pi**2/4 at the Nyquist mode.

The spectral model requires the periodic boundary and a route that applies the Fourier energies: the classical split_step flavor, or the binary encoding's quantum circuit and ir_product flavor. The one-hot circuit, the one-hot ir_product flavor and the schrodinger flavor apply the finite-difference stencil, and planning refuses the spectral model there. The finite-difference model is the default because every route applies it.

Rotation threshold

On the one-hot encoding every compiled hopping or number-projector block whose nonzero angle magnitude is below rotation_threshold is omitted. On the binary encoding it removes only Walsh-string Rz rotations. The QFT controlled phases and the angles of a dense phase diagonal stay, and binary_synthesis.aqft_cutoff truncates the QFT instead. The initial-state preparation is never pruned. Each omission is an approximation, which the Plan records with the dropped count and an operator-norm bound on its effect, pruning_error_bound.

State preparation

initial_state_preparation="structured" emits the library's structured circuit for the encoding and initial state, and the Plan records its CX count as a ResourceLaw. On the one-hot encoding it is the nonnegative amplitude chain of O(dK) gates, which for the uniform state is the linear W-state chain. The chain uses computed rotation parameters and can stop at its lower-range cutoff. The error bound of the omitted links is recorded per register in range_omissions.chains and included in the state_preparation error source. On the binary encoding it is one H per qubit, with no CX, which prepares the uniform state and the periodic kinetic ground state, the same state. Other states need "qiskit_state_preparation" there.

"qiskit_state_preparation" appends one Qiskit StatePreparation per register, of 2**K one-hot or K binary amplitudes. "none" selects a quantum resource-only construction without preparation. Classical execution prepares no circuit and rejects "none".

Classical flavors

"schrodinger" numerically evaluates each step's Hamiltonian exponential with SciPy's expm_multiply, whose work grows with the step's kinetic and potential weights and the objective's range. "ir_product" evaluates a numerical reference of the compiled product. It applies each one-hot block's analytic restricted action at its stored angle, and for the binary encoding the QFT operations with computed diagonal phase arrays. "split_step" applies the symmetric (Strang) product of half potential phases and the kinetic factor, diagonal in the kinetic eigenbasis and evaluated numerically by orthonormal DST-I on the Dirichlet grids and the FFT on the periodic grid. Its work per step depends on neither the schedule weights nor the objective's range, and it adds a splitting error that "schrodinger" does not have.

Raises:

ValueError –

If boundary="periodic" is combined with include_boundary_points=True, or with the one-hot encoding and an odd num_grid_points or one below 4. If kinetic_model="spectral" is used with a Dirichlet boundary. If encoding="binary" is used without the periodic boundary or with a num_grid_points that is not a power of two. If binary_synthesis differs from its default with the one-hot encoding.

Examples:

The minimum of (x - 0.2)**2 + (y + 0.2)**2 on [-1, 1]**2 is 0 at (0.2, -0.2), which is a point of the 4-point Dirichlet interior grid (-0.6, -0.2, 0.2, 0.6) of each axis. With exact readout on the default Aer simulator the candidate is the point of least evaluated objective among the grid points with positive probability, here that grid minimum. The evolution also puts the most probability there.

>>> import sympy as sp
>>> from nwqlib import solve
>>> from nwqlib.problems import Optimization
>>> from nwqlib.algorithms.qhd import QHD
>>> x, y = sp.symbols("x y", real=True)
>>> problem = Optimization(objective=(x - 0.2)**2 + (y + 0.2)**2,
...                        variables=(x, y),
...                        bounds=((-1.0, 1.0), (-1.0, 1.0)))
>>> method = QHD(num_grid_points=4, num_steps=8, total_time=4.0)
>>> result = solve(problem, method=method, seed=7)
>>> print([round(v, 6) for v in result.candidate],
...       round(result.objective, 12))
[0.2, -0.2] 0.0
>>> print([round(v, 6) for v in result.most_probable_coordinates],
...       result.mode_status)
[0.2, -0.2] resolved

BinarySynthesis ¶

Bases: Record

Circuit choices of the binary encoding: phase diagonals, QFT bit reversal and AQFT cutoff.

Build it with keyword arguments and pass it as QHD(encoding="binary", binary_synthesis=BinarySynthesis(...)). Every argument is optional, and the defaults build exact circuits. Every choice builds the same exact operator unless aqft_cutoff or a positive QHD.rotation_threshold selects an approximation, whose operator-norm bound the Plan records as aqft_error_bound and pruning_error_bound. The circuit synthesis section of the guide gives the CX counts of each choice.

Attributes:

potential (SynthesisChoice) –

Default "min_cx". Synthesis of each support table's phase diagonal: "dense_diagonal", "walsh_rotations", or "min_cx", the one of the two with fewer CX for each table, Walsh rotations on a tie. The Synthesis choices note below gives their CX counts and what "min_cx" does not minimize.
kinetic_phase (SynthesisChoice) –

Default "min_cx". The same choice for each variable's kinetic phase diagonal in the Fourier basis.
qft_bit_reversal (Literal['relabel', 'swap']) –

Default "relabel", which omits the swap layers of both QFTs and applies the phase table permuted by bit reversal between the swap-free QFT and its exact inverse, the same operator. "swap" keeps Qiskit's swap layers, 3 floor(b/2) CX per QFT.
aqft_cutoff (StrictInt | None) –

Default None, which keeps every controlled phase of the QFT. A nonnegative integer m keeps the controlled phases between wires at distance r <= m, whose angles are pi/2**r, and drops the rest, Qiskit's approximation_degree = max(0, b - 1 - m). m = 0 drops every controlled phase.

Synthesis choices

"dense_diagonal" is the exact phase diagonal of Shende, Bullock and Markov, quant-ph/0406176v5, Theorem 7, 2**n - 2 CX on n qubits for a nonzero table and none for an all-zero one. "walsh_rotations" applies one parity-ladder Z rotation per nonzero Walsh coefficient, 2 (w - 1) CX for a weight-w string. "min_cx" compares the two emitted CX counts of each table at its step exponent and takes the smaller, preferring "walsh_rotations" on a tie. It minimizes that CX count table by table, not the Rz or T count, the depth, routed CX on limited connectivity, or the CX of a jointly synthesized potential.

Schedules¶

The guide's schedules section compares the three schedules. QHD.coefficient_rule selects whether a step uses the point values at its midpoint or the step averages of the interval integrals. Mathematics derives the point values, step integrals and derivative bounds that planning and evolution_bound use, with their rounding.

QuadraticSchedule ¶

Bases: Record

Quadratic schedule a(t) = 1/(1 + gamma t**2), b(t) = 1 + gamma t**2, the QHD default.

Build it as QuadraticSchedule(gamma=0.3) and pass it as QHD(schedule=...). gamma is the only argument and is optional. Kushnir, Leng, Peng, Fan and Wu give phi_t = -log(1 + gamma t**2) and chi_t = log(1 + gamma t**2) as an example in Sec. 2.1 of QHDOPT, arXiv:2409.03121v1, and report that schedules of this form work well for many test problems. For gamma > 0 the ratio a/b = (1 + gamma t**2)**-2 tends to zero, the condition after Eq. (1) of Leng et al., arXiv:2303.01471v1. gamma = 0 gives the time-independent model a = b = 1. A decaying ratio does not by itself establish convergence on a finite grid or in a finite time. The guide's schedules section compares the three schedules.

The methods give the point values, interval integrals and derivative bounds that planning and evolution_bound use, with their rounding in units of u = 2**-53.

Attributes:

kind (Literal['quadratic']) –

Default "quadratic", the only accepted value. Formula identifier.
gamma (Nonnegative) –

Default 0.3, an illustrative working point and not an accuracy choice. Nonnegative finite rate, in inverse squared time units.

CubicSchedule ¶

Bases: Record

Cubic schedule a(t) = 2/(s + t**3), b(t) = 2 t**3 of Leng et al., Eq. (C.4).

Build it as CubicSchedule(s=...) and pass it as QHD(schedule=...). s is required and has no default, because it is a model parameter that the caller chooses. This is Eq. (C.4) of Leng et al., arXiv:2303.01471v1, p. 32, whose two weights follow the ODE model of Nesterov's accelerated gradient method, with s removing the singularity of 2/t**3 at t = 0. Liu et al., arXiv:2607.16996v1, Eq. (92), write it with s = 1 and call it QHD-C, but their code runs the shifted form (ShiftedCubicSchedule) under that name, and with s = 1 only that form reproduces the Rz counts of their Table II (circuit synthesis). s is a model parameter. Leng et al. set it to their time step, but here it stays fixed when the step count changes, because changing s changes the Hamiltonian. For the classical Hamiltonian a |p|**2/2 + b f(x), Hamilton's equations x' = a p and p' = -b grad f give x'' = (a'/a) x' - a b grad f, so the damping is -a'/a and the gradient factor a b. Here they are 3 t**2/(s + t**3) and 4 t**3/(s + t**3), which tends to 4, and they differ from those of ShiftedCubicSchedule.

The methods give the point values, interval integrals and derivative bounds that planning and evolution_bound use, with their rounding in units of u = 2**-53.

Attributes:

kind (Literal['cubic']) –

Default "cubic", the only accepted value. Formula identifier.
s (Positive) –

Required. Positive finite regularization, in cubed time units.

ShiftedCubicSchedule ¶

Bases: Record

Shifted cubic schedule a(t) = (2/(s + t))**3, b(t) = 2 t**3 of Wu et al., Eq. (15).

Build it as ShiftedCubicSchedule(s=...) and pass it as QHD(schedule=...). s is required and has no default, because it is a model parameter that the caller chooses. This is Eq. (15) of Wu et al., arXiv:2605.12066v1, Sec. VI, the schedule of their published results, offered for replication. The paper calls s a small regularization parameter, and the code behind its results sets s = T/N_t = 2e-4 for T = 10 and the paper's N_t = 50,000 time steps. That code uses N_t as the number of time points, so its step T/(N_t - 1) is close to s but not equal. The paper presents the schedule as the QHD-C schedule of Leng et al., arXiv:2303.01471v1, but its dynamics differ from Leng's Eq. (C.4), which is CubicSchedule. The damping -a'/a and gradient factor a b (CubicSchedule derives them) are 3/(s + t) and 16 t**3/(s + t)**3, which tends to 16 instead of 4. Liu et al., arXiv:2607.16996v1, run this form under the name QHD-C in their code, and with s = 1 only this form reproduces the Rz counts of their Table II (circuit synthesis). As for the cubic form, s stays fixed when the step count changes.

The methods give the point values, interval integrals and derivative bounds that planning and evolution_bound use, with their rounding in units of u = 2**-53.

Attributes:

kind (Literal['shifted_cubic']) –

Default "shifted_cubic", the only accepted value. Formula identifier.
s (Positive) –

Required. Positive finite time shift, in time units.

Initial states¶

Each initial state is a product of nonnegative per-variable amplitude vectors, and QHD.initial_state_preparation selects how quantum execution prepares it. The guide's initial states section compares them. Mathematics gives the construction error of each vector.

KineticGroundState ¶

Bases: Record

Ground state of each variable's kinetic operator, the default initial state of QHD.

Pass it as QHD(initial_state=KineticGroundState()), which is the default. It takes no arguments. The row "QHD Method defaults" of the engineering constants gives the reasons for this default and its limits.

On both Dirichlet grids the kinetic operator of one variable, restricted to the grid, is T/h**2 with the K-by-K tridiagonal matrix T of diagonal 1 and off-diagonals -1/2. The interior grid and the grid with boundary points differ only in the spacing h, so they share the eigenvectors of T. For v_i = sin(r pi (i + 1)/(K + 1)), i = 0..K-1 and r = 1..K, the identity sin(a - b) + sin(a + b) = 2 sin(a) cos(b) gives (T v)_i = v_i - (v_(i-1) + v_(i+1))/2 = (1 - cos(r pi/(K + 1))) v_i, where the missing neighbors v_(-1) = sin(0) and v_K = sin(r pi) vanish as the stencil requires. The eigenvalues 2 sin(r pi/(2 (K + 1)))**2/h**2 increase with r, so r = 1 is the ground state, sin(pi (i + 1)/(K + 1)) with every entry positive and energy 2 sin(pi/(2 (K + 1)))**2/h**2. The kinetic operator of d variables is the sum of the one-variable operators on separate tensor factors, so its ground state is the product of these vectors and its energy the sum of theirs. The relation is standard linear algebra, derived here in full. On the grid with boundary points the vector does not vanish at the box endpoints, because that grid's missing neighbors lie outside the box.

On the periodic grid the operator is C/h**2 with the circulant C = I - (S + S^T)/2 and the cyclic shift S. The Fourier vectors f_j = exp(2 pi i r j/K) satisfy (C f)_j = (1 - cos(2 pi r/K)) f_j, so the eigenvalues 2 sin(pi r/K)**2/h**2, r = 0..K-1, are zero only for r = 0, whose eigenvector is the constant vector. The ground state is therefore the uniform state with energy zero, and this record then prepares exactly what UniformState prepares. The spectral kinetic model (QHD.kinetic_model="spectral"), which the Method accepts on the periodic grid only, has the same Fourier eigenvectors with eigenvalues 2 pi**2 q**2/L**2 for the signed frequency index q and period L = K h. They also vanish only for q = 0, so the uniform state is the ground state under both periodic kinetic models.

Attributes:

kind (Literal['kinetic_ground']) –

Default "kinetic_ground", the only accepted value. Initial-state identifier.

UniformState ¶

Bases: Record

Uniform initial state: equal amplitude on every valid grid point.

Pass it as QHD(initial_state=UniformState()). It takes no arguments. Each variable has amplitude 1/sqrt(K) on each of its K grid points, so the state has 1/sqrt(K**d) on every grid tuple. This is the uniform superposition of Leng et al., arXiv:2303.01471v1, Algorithm 1 step 3, and the initial state that Wu et al., arXiv:2605.12066v1, Sec. VI, state for their simulations. On the Dirichlet grids it has a component outside the kinetic ground state (KineticGroundState), which the guide's initial states section discusses.

Attributes:

kind (Literal['uniform']) –

Default "uniform", the only accepted value. Initial-state identifier.

GaussianState ¶

Bases: Record

Gaussian warm start prod_j exp(-(x_j - c_j)**2/(2 sigma_j**2)), normalized on the grid points.

Build it with keyword arguments, for example GaussianState(center=(0.5, -0.2), widths=(0.3, 0.3)), and pass it as QHD(initial_state=...). center and widths are both required, with one entry per Problem variable. The center c is a point in the coordinates of the Problem that the Plan solves, and the widths sigma_j are in the units of each of its variables. For solve these are the original coordinates. Box refinement plans each level's own problem with the QHD configuration's state, and its default search model solves in the unit coordinates u = (x - a)/L of each level box [a, a + L], so there c and sigma are unit coordinates. The center may lie outside the box. The amplitude of a grid tuple is the product of the per-variable factors at its grid coordinates, and the state is that product normalized over the valid grid points.

Leng et al., arXiv:2303.01471v1, Algorithm 1 step 3, name a Gaussian state as one choice of initial state, and the code behind the refinement results of Wu et al., arXiv:2605.12066v1, starts every refinement level after the first from this amplitude, centered at the best point found so far (BoxRefinement(level_initial_state="best_point_gaussian")). The amplitude, not the probability, has width sigma, so the continuous probability density, proportional to exp(-(x_j - c_j)**2/sigma_j**2), has standard deviation sigma_j/sqrt(2) along variable j. A Gaussian contains several momentum components of the kinetic operator, so it does not avoid the time-discretization dependence that the kinetic ground state avoids (initial states). On the periodic grid the distance is the chart difference |x_ji - c_j| of the coordinates in [lower, upper), not the distance on the circle, so the amplitudes do not continue across the wrap link. For a center near one end of the period, the points near the other end, which are its neighbors through that link, get the amplitudes of their chart distance. A minimum-image variant is open work (Limitations and open work).

Every finite center and every finite positive width is accepted. Each variable's exponents z_ji = (x_ji - c_j)**2/(2 sigma_j**2) are shifted by their minimum before exponentiation, so the largest amplitude is exactly 1 before normalization and the state never underflows to zero. The construction error bound of the classical start vector is finite for every input. It grows with the exponents, about 10u times the exponent per entry, with u = 2**-53, and it is capped by the distance of two nonnegative vectors of norm about one, about sqrt(2). The per-variable bound of variable_errors is infinite where the exponents give no bound.

Attributes:

kind (Literal['gaussian']) –

Default "gaussian", the only accepted value. Initial-state identifier.
center (tuple[Real, ...]) –

Required. Finite center coordinate of each variable, in the Problem's variable order and the coordinates of the solved Problem, which are unit coordinates for a search-model refinement level.
widths (tuple[Positive, ...]) –

Required. Positive finite width sigma of each variable, in the same coordinates.

Raises:

ValueError –

If center is empty or center and widths differ in length. Planning also raises when they do not have one entry per Problem variable.

Read the result¶

QHDAnalysis ¶

Bases: Result

Readout of a QHD run: best observed point, most probable point and probability masses.

solve returns it for a QHD method, and load_result reopens a saved one. The answer is candidate, the valid grid point with positive observed weight that has the least evaluated objective, and objective, the objective there from the stored support tables, in the Problem's objective unit. most_probable_coordinates is the grid point where the evolution concentrated probability, which the augmented-Lagrangian and refinement layers read by default. The two points can differ. Both are None when no valid outcome was observed, and neither establishes a continuous or global optimum. print(result) shows both points, the mode status and the valid probability. The fields below are read-only. The fields of Result are present too. The guide's readout section explains both points and the tie window.

The candidate minimizes the evaluated binary64 objective among valid points with positive observed weight, with the smallest grid index on objective ties. For exact readout these are positive-probability points, and for counts they are sampled points. Exact probability readout does not make objective evaluation exact or guarantee that every grid point has positive probability. The most probable point describes concentration of the observed distribution under the tie rule of the Mode status note below. Its fields are unavailable exactly when the candidate is.

Attributes:

value (Real | None) –

Objective at the candidate, read from the stored support tables, or None when no valid outcome was observed. The same as objective.
candidate_indices (tuple[Count, ...] | None) –

Grid index of the candidate for each variable.
candidate_coordinates (tuple[Real, ...] | None) –

Candidate grid point in the original coordinates. The same as candidate.
candidate_probability (Nonnegative | None) –

Unconditional probability of the candidate.
most_probable_indices (tuple[Count, ...] | None) –

Grid index for each variable of the most probable valid point. Among the valid points with positive observed probability, it is the lexicographically smallest index tuple whose probability lies within most_probable_tie_window of the computed maximum. Counts compare pooled integer counts exactly.
most_probable_coordinates (tuple[Real, ...] | None) –

That point in the original coordinates.
most_probable_probability (Nonnegative | None) –

Its unconditional probability.
most_probable_objective (Real | None) –

Objective at that point, read from the stored support tables.
most_probable_deficit (Nonnegative | None) –

Computed maximum probability minus most_probable_probability, at most the tie window. A positive value means a point within the window was chosen over the computed maximum. Other points can lie within the window even when it is zero, so the selection does not identify a unique mode.
most_probable_tie_window (Nonnegative | None) –

Derived bound on the error of a computed difference of two probabilities on this readout path, 0 for counts, or None when no bound was derived. The classical split_step kernel derives it from the state budget it observed on its trajectory, at most the Plan's window.
most_probable_tie_window_unavailable (Text | None) –

Why no tie window was derived, for example a preparation record outside the roundoff derivation. The selection then compared the computed probabilities exactly, so roundoff can decide between points that the model makes equal.
probability_maximizer_indices (tuple[Count, ...] | None) –

Lexicographically first valid grid-index tuple attaining the largest positive observed weight. Exact readout compares the computed probabilities. Counts compare pooled integers before division by returned_shots. None without a positive valid outcome.
probability_maximizer_coordinates (tuple[Real, ...] | None) –

Coordinates of the computed probability maximizer on the Plan's grid, or None with its indices.
probability_maximizer_probability (Nonnegative | None) –

Unconditional computed probability, or empirical frequency rounded once from its pooled count and returned_shots, at the computed probability maximizer. None with its indices.
probability_maximizer_objective (Real | None) –

Evaluated objective from the Plan's support tables at the computed probability maximizer. None with its indices.
mode_status (ModeStatus) –

Resolution of the most-probable selection among positive observed valid points, resolved, unresolved or unavailable. With an available window, resolved means only one point passes the maximum-to-point tie test, and unresolved means another point passes it. Counts use equality of pooled integer counts. Unavailable means no positive valid outcome or no derived numerical window. This status supplies no sampling-error or optimization guarantee.
expected_objective (Real | None) –

Objective mean conditional on a valid outcome.
marginals (FrozenArray) –

Unconditional probability of each grid index per variable, a FrozenArray whose .array is the read-only float64 array of shape (d, K) in variable order and grid order. Each row sums to valid_mass. Analysis accumulates it from the observed points, and for classical execution reads it from the marginals that the kernel returns.
valid_mass (Nonnegative) –

Probability of outcomes that encode a grid point, one excitation per register for one-hot and every outcome for binary. The same as valid_probability.
invalid_mass (Nonnegative) –

Probability of all other register outcomes, zero for binary.
observed_mass (Nonnegative) –

Total observed probability.
valid_count (Count | None) –

For counts, the number of returned outcomes that encode a grid point, summed as integers over all chunks. None for exact readout. With positive returned_shots, valid_mass is valid_count/returned_shots rounded once.
returned_shots (Count | None) –

For counts, the shots that the chunks returned, summed as integers, which can be fewer than the Plan requested. None for exact readout.
missing (tuple[Text, ...]) –

Reasons the observation is incomplete or yields no valid candidate.
applications (tuple[KernelApplication, ...]) –

Records of the classical kernel's run, for classical execution.
artifact (ArtifactManifest | None) –

Saved final state, when keep_state requested it.

Mode status

Write O for the valid points with positive observed weight, q_i for the stored probability of point i (for counts the exact empirical frequency N_i/R of the pooled integer count N_i and the integer returned_shots R, compared through N_i), M for the largest q_i over O, m for the computed probability maximizer (the probability_maximizer_* fields) and s for the tie representative (the most_probable_* fields). Under an available window W the representative is the lexicographically first point passing the maximum-to-point tie test fl(M - q_i) <= W, with fl the binary64 subtraction, and for counts the lexicographically first point with the largest pooled count, W = 0. Without a window the stored probabilities are compared exactly, so s = m. mode_status is unavailable when O is empty or W is unavailable, unresolved when W is available and a point of O other than s passes the same test, and resolved otherwise, that is when exactly one point passes it. A zero deficit occurs both at a well-separated maximum and at a tie, so the status is evaluated from the observed weights together with the selection and cannot be reconstructed from the deficit.

Guarantee

Suppose W bounds the pairwise errors of the chosen reference population, abs((q_i - q_j) - (p_i - p_j)) <= W for every compared i, j. If exactly one point passes the test, it is m, hence s = m, and fl(q_m - q_t) > W for every other t. Rounding is monotone and fixes the stored binary64 W, so an exact difference at most W would round to at most W, and the rejection implies the exact real difference q_m - q_t > W, so p_m - p_t >= q_m - q_t - W > 0. resolved therefore certifies a unique maximum of the stored q and, under the stated pairwise error model, a unique maximum of p on O. Exact subtraction, including Sterbenz's condition q_t >= q_m/2, suffices for this conclusion but is not necessary for it. Rounding can make the status conservatively unresolved, since a rounded subtraction can accept an exact gap slightly above W, and it cannot make it falsely resolved under the pairwise assumption. For counts, a single largest integer N_i gives a unique maximum of the exact empirical frequencies N_i/R, and a count tie is unresolved. Neither result is a confidence statement about an underlying distribution. unresolved says that this window does not separate the representative from every competitor. It does not prove equal probabilities, multiple exact modes, a large numerical error or an inferior objective. resolved says nothing about optimization quality, a continuous minimizer, physical-model discrepancy, hardware bias or finite-shot uncertainty. The population is O. If an exact readout completely represents a valid grid distribution, an omitted known-zero q_i cannot exceed the positive maximum, so m is also a computed maximizer on the full valid grid, but a certificate about p on O extends to omitted points only when their pairwise error bounds apply and fl(M - 0) > W, namely M > W. An incomplete observation does not make omitted probabilities known zero.

candidate ¶

candidate

Candidate point in the original coordinates, the same as candidate_coordinates.

objective ¶

objective

Original objective at the candidate, the same as value.

valid_probability ¶

valid_probability

Unconditional probability of the outcomes that encode a grid point, the same as valid_mass.

position_mean ¶

position_mean

Mean grid coordinate of each variable conditional on a valid outcome, or None without valid mass.

<x_j> = sum_i x_j(i) m_j(i) / valid_mass, with m_j the stored unconditional marginal of variable j and x_j(i) the Plan's grid coordinates. Liu et al. arXiv:2607.16996v1 Eq. (94) define the mean over the whole final state. Dividing by valid_mass conditions it on a valid outcome, NWQLib's choice, because an invalid one-hot outcome encodes no position. It is computed on demand from the marginals and the attached Plan's grid only, with no evolution, decoding, objective evaluation or measurement, and is not stored. The mean and position_standard_deviation describe the observed distribution. They are no criterion for agreement of two states or for optimization success, because two distributions can share mean and variance at total-variation distance 1. For example, {-1, +1} with probability 1/2 each and {-2, 0, +2} with {1/8, 3/4, 1/8} both have mean 0 and variance 1, and put probability 0 and 3/4 on the point 0. On the periodic grid the coordinates are those of the chart [lower, upper), so both moments 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, where it may have no mass.

position_standard_deviation ¶

position_standard_deviation

Standard deviation of each variable's grid coordinate conditional on a valid outcome, or None.

sigma_j = sqrt(sum_i (x_j(i) - <x_j>)**2 m_j(i) / valid_mass) with <x_j> from position_mean, conditioned on a valid outcome as the mean is. Liu et al. arXiv:2607.16996v1 Eq. (94) define the standard deviation over the whole final state. The deviations are scaled by their largest magnitude before squaring, so coordinates far from the origin do not overflow. None when valid_mass is zero.

Check a result¶

QHDVerification ¶

Bases: Record

Options of an explicit QHD check: grid minimum, or fidelity against a reference evolution.

Build it with keyword arguments, for example QHDVerification(comparisons=("grid_minimum",)), and pass it to result.verify(checks=...), which returns (receipt, facts). comparisons is the only required argument. solve never runs these checks. Replaying the original classical flavor checks consistency, not independent accuracy, and the receipt names what produced the result and the reference it compares. The guide's explicit verification section describes what each comparison checks and what it does not establish. The default tolerances are untuned, and the engineering constants register them.

Attributes:

comparisons (tuple[Literal['grid_minimum', 'schrodinger_fidelity', 'ir_product_fidelity'], ...]) –

Required. Distinct checks to run, a nonempty tuple of "grid_minimum", "schrodinger_fidelity" and "ir_product_fidelity". grid_minimum evaluates the original objective on every grid point. Each fidelity choice runs one restricted evolution against the stored state.
objective_gap_tolerance (Nonnegative) –

Default 1e-12. Pass threshold, in objective units, for the candidate's evaluated gap to the least reevaluated binary64 value.
schrodinger_infidelity_tolerance (Real) –

Default 0.1, between 0 and 1. Pass threshold for the infidelity against the restricted Schrodinger evolution.
ir_product_infidelity_tolerance (Real) –

Default 1e-9, between 0 and 1. Pass threshold for the computed infidelity against the numerical reference of the compiled product. Binary phase reconstruction and reference-evolution error are not propagated into this tolerance decision.
minimum_tolerance (Nonnegative) –

Default 0.0. Window around the least reevaluated binary64 objective value that minimum_success_mass uses. Zero counts ties in the evaluated table. It counts exactly the mathematical minimizers when the evaluations are exact and this tolerance is zero.
max_bytes (PositiveInt) –

Default 10 GB (decimal, 10_000_000_000 bytes). Positive limit on the known array bytes of this check, separate from the QHD Method's limits.
max_work (PositiveInt) –

Default 1_000_000_000. Positive limit on the known work of this check, separate from the QHD Method's limits.

Raises:

ValueError –

If comparisons is empty or repeats a choice.

Estimate fault-tolerant resources¶

circuit_resources(plan, synthesis_epsilon=E) reads a quantum QHD Plan and builds no circuit. The guide's fault-tolerant resources section states what the rotation count includes, how the synthesis budget is split and which error sources are bounded.

circuit_resources ¶

circuit_resources(plan, *, synthesis_epsilon=None)

Count the rotations of one QHD circuit, estimate its T gates and list its error sources.

circuit_resources(plan, synthesis_epsilon=E) reads a QHD Plan with execution="quantum" and the structured or no initial-state preparation, in either encoding, and builds, transpiles and compiles no circuit. The record holds the arbitrary rotations and exact T gates of the emitted circuit, the evolution_bound of its product formula and every error source against the finite model, each with its status. With a synthesis budget E it also replaces by Clifford gates the rotations that the budget pays for and estimates the T count of the rest, the leading term T_exact + 3 M log2(1/epsilon_rot) of the typical Ross-Selinger count (arXiv:1403.2975v3). The guide's fault-tolerant resources section gives the replacement rule, the budget split and a comparison with a compiled count.

Error sources, in the order of the stages they compare. The finite model uses the exact schedule functions with the stored tables, constant and binary64 spacings on the Plan's grid, and the compiled product uses the stored step weights.

kinetic_model: the replacement of the finite-difference stencil by the finite model's kinetic operator. Not applicable to the finite-difference model, which every one-hot circuit applies. Unavailable for the spectral model of a binary circuit, against which every later entry is measured. The two kinetic operators differ by a norm that grows as 1/h**2, and Duhamel's formula bounds the change of the state only through the momentum content of the whole trajectory (Split-step classical evaluation), which the record does not have.
time_ordering, midpoint_quadrature (not applicable under the integrated coefficient rule), coefficient_rounding and product_formula: the stages of evolution_bound. The coefficient residual is exact except for the quadratic (gamma > 0) and cubic kinetic integrals, whose first-order estimate makes it an estimate.
angle_formation: one-hot, the stored angles of the kept blocks against their intended exponents, exact. Binary, every computed Walsh angle before pruning, the dense phases and the QFT angles formed from binary64 pi against the exact stored-weight split product, under the one-ulp assumption for math.sin, math.cos and NumPy's sin, Walsh coefficients omitted below the normal range included.
aqft: not applicable to the one-hot circuit, which applies no Fourier transform, and to a binary circuit with exact QFTs. With a cutoff it is the Plan's aqft_error_bound, which bounds the full stored-angle QFTs against their truncation.
rotation_pruning: one-hot, the blocks the circuit omits, pruned by rotation_threshold or omitted below the normal binary64 range, counted at their exact intended exponents. Binary, the Plan's pruning_error_bound, half the dropped Rz angles plus the exact contributions of the rotations and dense phase entries omitted below the normal range.
gate_parameters: one-hot, the fused hopping gate's own parameter against twice the stored angle, zero since the builder doubles the stored angle's own product, with the projector and chain rotations exact power-of-two scalings of their stored angles, which planning checks to lie in the normal binary64 range. Binary, zero, since every angle reaches its gate unchanged or halved exactly.
diagonal_wrap: the phase wrap angle(exp(-i a)) of a phase diagonal, the one-hot projectors on 4 to 7 qubits or the binary dense diagonals, whose circuit construction then forms recursive means and Gray-code angles. NumPy documents no uniform error constant for its complex exponential and argument, so it is unavailable when the circuit has such a block and not applicable otherwise.
identity_phase: the circuit's global-phase bookkeeping against the exact identity phase, the identity actions omitted below the normal range included.
state_preparation: one-hot, the prepared state against the exact target, a first-order estimate whose chain-cutoff contribution is exact. Binary, zero, since the H layer prepares the uniform state exactly. Not applicable without preparation.
clifford_replacement: E_C of the synthesis projection, an outward bound of the replaced rotations' distances from their Clifford gates, describing the approximated construction, not the emitted circuit.
synthesis: the requested allowance E - E_C, rounded downward, of the remaining rotations. NWQEC 0.1.2 applies a fixed 1e-4 angle cleanup, groups angles to four significant digits and passes them as decimal strings, all outside the requested GridSynth tolerance, so an achieved synthesis error is unavailable (fault-tolerant resources).
compiler: other compiler transformations, unavailable.

For a normalized input these add by telescoping to a bound on the final state's 2-norm error when every entry is a bound on adjacent stages, against the finite model without kinetic_model and against the finite-difference model with it. Each entry compares two adjacent whole evolutions, so the entries add by the triangle inequality through the intermediate evolutions, with no independence assumption. Within an entry, for ideal factors U_j and their implementations U~_j, the identity U~_L ... U~_1 - U_L ... U_1 = sum_j U~_L ... U~_(j+1) (U~_j - U_j) U_(j-1) ... U_1 bounds the product's error by the sum of the factors' errors, because the unitaries around each difference have norm one. The preparation error adds once, because later unitaries keep its norm. Any common measurement then changes by at most that amount in total variation, capped at 1, because measurement cannot increase the trace distance sqrt(1 - |<psi|psi~>|**2) of two pure states, which is at most their phase-aligned 2-norm distance. An entry built from an estimate makes such a sum conditional on it, and an unavailable entry leaves it unavailable, so the record forms no total. Sampling, device noise and the optimization gap are outside these error sources.

Cost. The one-hot work is one pass over the stored blocks with a few exact rational operations per block, and the bound's pass over tables and steps. The binary work reads the stored Walsh phase terms, synthesizes each distinct stored block once for its angles, with the table transforms of planning, adds the O(n 2**n) Gray-code angles of each distinct dense-diagonal block on n qubits and two upward reductions per table for the angle formation, and builds no circuit. The binary rotation count is first checked against the Plan's QHD.max_work and QHD.max_bytes.

Parameters:

plan (Plan) –

A QHD Plan from nwqlib.plan(problem, method=QHD(...)) with execution="quantum", the default, and initial_state_preparation "structured" or "none".
synthesis_epsilon (float | None, default: None ) –

Total operator-norm allowance E of one circuit execution for Clifford replacement plus rotation synthesis, positive and finite. It has no default, because a budget is the user's accuracy choice. Without it the record has no synthesis projection or T estimate.

Returns:

record ( QHDCircuitResources ) –

Per circuit execution, the rotation counts in arbitrary_rotations and exact_t, the T estimate in synthesis.t, the evolution bound in evolution and the error sources in error_sources.

Raises:

ValueError –

If the Plan is classical, if its state preparation has no rotation count ("qiskit_state_preparation"), if the binary rotation count exceeds the Plan's QHD work or byte limit with the model's baseline, or if a budget leaves a per-rotation tolerance that rounds down to zero.

Examples:

One first-order step of (x - 1/5)**2 on three interior grid points of [-1, 1] has 8 arbitrary rotations and 2 exact T gates. At E = 1e-4 no rotation is replaced, so epsilon_rot = 1e-4/8 and the estimate is 2 + 24 log2(8e4) = 392.9 T. The evolution bound is the splitting term alone, because QuadraticSchedule(gamma=0.0) makes the Hamiltonian constant.

>>> import sympy as sp
>>> from nwqlib import plan
>>> from nwqlib.problems import Optimization
>>> from nwqlib.algorithms.qhd import (
...     QHD, QuadraticSchedule, UniformState, circuit_resources)
>>> x = sp.Symbol("x", real=True)
>>> problem = Optimization(objective=(x - sp.Rational(1, 5))**2,
...                        variables=(x,), bounds=((-1.0, 1.0),))
>>> method = QHD(num_grid_points=3, total_time=0.17, trotter_order=1,
...              schedule=QuadraticSchedule(gamma=0.0),
...              initial_state=UniformState())
>>> record = circuit_resources(plan(problem, method=method),
...                            synthesis_epsilon=1e-4)
>>> print(record.arbitrary_rotations, record.exact_t)
8 2
>>> print(round(record.synthesis.t.value.value, 1))
392.9
>>> print(record.evolution.evolution, record.evolution.evolution_status)
0.07225000000000002 bound

evolution_bound ¶

evolution_bound(plan)

Bound the error of a QHD Plan's compiled product against its finite model, in the operator norm.

evolution_bound(plan) returns a QHDEvolutionBound for a Plan whose step blocks were compiled, which quantum execution and the ir_product flavor do, in either encoding. circuit_resources includes the same record. Its evolution field bounds the compiled product with the stored exponents against the exact time-ordered evolution, with the physical phase included, conditionally where the coefficient residual is an estimate. It adds four stages: time ordering, midpoint quadrature, the coefficient residual and splitting. The work is one pass over the stored tables and steps, O(sum_S |S| K**|S| + d K + N), with no state or matrix.

Splitting (Childs et al., arXiv:1912.08854v3, Eqs. (145) and (152)). For dimensionless Hermitian generators G_j in application order and R_j = sum_(q > j) G_q, the first-order product errs by at most (1/2) sum_j ||[R_j, G_j]|| and the symmetric second-order product by at most (1/12) sum_j ||[R_j, [R_j, G_j]]|| + (1/24) sum_j ||[G_j, [G_j, R_j]]||. Unitarity makes these valid with no exponential growth factor.

First order applies beta_k V and then every link generator alpha_k L_(j,e) in emitted order, with L_(j,e) = -w_j (|p><q| + |q><p|) of norm w_j = 1/(2 h_j**2), the restriction of -(XX + YY)/(4 h_j**2). XX and YY on one link commute, so the fused hopping gate adds no split, and all projector blocks commute. The potential factor against all links gives |alpha_k beta_k| C/2, and the links among themselves alpha_k**2 Gamma/2 with Gamma = sum_j sum_e ||[sum_(f > e) L_(j,f), L_(j,e)]||. On a chain each link but the last has one later adjacent link, whose commutator has norm w_j**2. On an even cycle the first link has two later adjacent links, which give two disjoint skew-symmetric 2-by-2 blocks of norm w_j**2 together, and every other nonfinal link has one. Disjoint links and different variables commute, so Gamma = (K - 2) sum_j w_j**2 on the Dirichlet chain and (K - 1) sum_j w_j**2 on the periodic cycle (zero for K = 2), and

s_k = |alpha_k beta_k| C/2 + alpha_k**2 Gamma/2.

Second order applies exp(-i beta_k V/2) exp(-i alpha_k O/2) exp(-i alpha_k E) exp(-i alpha_k O/2) exp(-i beta_k V/2) after regrouping the per-variable layers, which commute across variables, for the odd and even layer sums O and E. The symmetric bound for the three generators beta_k V, alpha_k O and alpha_k E gives

where the first two terms hold every kinetic/potential commutator, since O + E is T up to an identity, and the last the internal link split, which a bound from [T, V] alone would miss (it is nonzero for a constant V). C, D_T and D_V come from the table ranges and, for the finite-difference stencil, from the largest potential change across each link, and the layer constants J_E and J_O from the K-point graph.

Binary encoding. An exact Fourier conjugation applies each variable's kinetic factor F^dagger exp(-i alpha_k diag(E_j)) F whole, and the bit-reversal relabeling conjugates the phase table as well, so it leaves T_j unchanged. Kinetic factors on different variables commute, and so do the potential tables. First order therefore has the two generators beta_k V and alpha_k T, and second order the symmetric beta_k V/2, alpha_k T, beta_k V/2, which give

with no link term, K = 2 included. The dense and unpruned Walsh diagonals implement the same factors. An approximate QFT, pruning and the rounding of the computed angles are later entries of the circuit's error sources (circuit_resources).

Schedule and coefficients. For step k of width Delta, w_k is the time-ordering budget, the smaller of the available integral bound min(2, C A B/2) and derivative bound Delta**3 (a* B1 + b* A1) C/12, q_k the midpoint-quadrature bound Delta**3 (A2 mu_T + B2 mu_V)/24 under the midpoint rule, and e_coeff,k the residual between the exact coefficients and the stored exponents. The guide's fault-tolerant resources section states when each bound is available. The totals are splitting = min(2, sum_k s_k), schedule = min(2, sum_k (w_k + q_k)), coefficient_residual = min(2, sum_k e_coeff,k) and evolution = min(2, splitting + schedule + coefficient_residual), which is conditional on the coefficient estimate where one is used. For the fixed smooth schedules the schedule bound falls as N**-2 at fixed final time, first-order splitting as N**-1 and second-order splitting as N**-2. A small schedule parameter s or a fine grid can make the bounds large, and the cap 2 is then valid but uninformative.

Parameters:

plan (Plan) –

A QHD Plan with compiled step blocks, from nwqlib.plan(problem, method=QHD(...)) with quantum execution or theory_flavor="ir_product".

Returns:

bound ( QHDEvolutionBound ) –

The norm inputs, the per-stage bounds and their total evolution with its evolution_status.

Raises:

ValueError –

If the Plan has no compiled step blocks, so there is no product to bound.

run_resources ¶

run_resources(result, *, synthesis_epsilon=None)

Total the rotations and T estimates of a quantum augmented-Lagrangian run or box refinement.

run_resources(result, synthesis_epsilon=E) reads a ConstrainedQHDResult, with or without box refinement in each round, or a BoxRefinementResult, run with execution="quantum". Each started round, or each started level of a round's refinement, contributes the circuit preparations, the shots set aside and the per-circuit rotation count that its resources recorded (ALResources, RefinementResources), and the T estimate of its live inner result's Plan. A round or level whose Run raised, and a refinement level that stopped the run, keep their recorded counts but have no inner result, so their T estimate is unknown when they prepared a circuit. A Plan with Qiskit's state preparation has no rotation count, which its resources record as unknown. Nothing is planned, compiled or executed. The work is one pass over each inner Plan's stored blocks, which for the binary encoding synthesizes each distinct block's diagonal once more to read its angles, after the rotation count is checked against that Plan's QHD limits.

Parameters:

result (ConstrainedQHDResult | BoxRefinementResult) –

The run with its inner results, as solve_augmented_lagrangian, refine_box or their load and resume functions return it.
synthesis_epsilon (float | None, default: None ) –

The per-circuit budget E of circuit_resources, or None for rotation totals without T estimates.

Returns:

record ( QHDRunResources ) –

The body totals arbitrary_rotations and t_estimate, the shot-weighted totals shot_arbitrary_rotations and shot_t_estimate, and one entry per started round or level.

Raises:

ValueError –

If the run used classical execution, which selects no circuit, or if an inner binary Plan's rotation count exceeds that Plan's QHD work or byte limit.

QHDCircuitResources ¶

Bases: Record

Fault-tolerant resources and error sources of one QHD circuit, from circuit_resources.

circuit_resources returns it for a quantum QHD Plan. All counts are per circuit execution. Shots and outer rounds multiply them, which run_resources totals. The fields below are read-only.

Attributes:

plan_id (ContentID) –

Content hash of the Plan.
arbitrary_rotations (Count) –

Arbitrary rotations of the emitted circuit, classified exactly by angle, from the stored blocks for one-hot and from the computed block angles for binary. For one-hot this is the value of rotation_law.
exact_t (Count) –

Exact T and T-inverse occurrences of the emitted circuit.
rotation_law (ResourceLaw) –

The Plan's arbitrary_rotations count, a ResourceLaw. The binary encoding's count includes the rotation gates before angle classification, exact T gates and the omitted zero angles of dense diagonals, so it is an upper bound on arbitrary_rotations.
synthesis (QHDSynthesisProjection | None) –

The budgeted Clifford replacement and T estimate, a QHDSynthesisProjection, or None without a synthesis budget.
evolution (QHDEvolutionBound) –

Splitting and schedule bounds, a QHDEvolutionBound.
error_sources (tuple[QHDErrorSource, ...]) –

The replacement of the finite-difference kinetic, then every error source of the circuit against the finite model, in stage order, each a QHDErrorSource.

QHDSynthesisProjection ¶

Bases: Record

Budgeted Clifford replacement and a leading-order T estimate for one QHD circuit.

circuit_resources(plan, synthesis_epsilon=E) returns it as synthesis. The fields below are read-only. The projection describes an approximated construction. The emitted circuit is unchanged, and replaced rotations are Cliffords only in the construction this record prices. With N arbitrary rotations, every rotation whose distance from the nearest Clifford rotation is at most q = E/N becomes that Clifford, and the remaining budget is split evenly over the M rotations left. The T estimate is T_exact + 3 M log2(1/epsilon_rot), the leading term of the typical Ross-Selinger count (arXiv:1403.2975v3) with intercept zero. It treats the circuit's angles 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 guide's fault-tolerant resources section derives the rule.

Attributes:

epsilon (Nonnegative) –

E, the total operator-norm allowance of one circuit execution for Clifford replacement plus rotation synthesis.
candidates (Count) –

N, the arbitrary rotations of the emitted circuit.
threshold (Nonnegative | None) –

q = E/N, the initial per-rotation allowance, rounded downward, None when N is zero.
replaced (Count) –

Candidates whose distance bound 2 sin(|delta|/4), with delta the angle minus the nearest multiple of pi/2, is at most q, replaced by Cliffords.
replacement_error (Nonnegative) –

E_C, the exact sum of their distance bounds, rounded upward.
rotations (Count) –

M = N - replaced, the rotations left to synthesize.
rotation_epsilon (Nonnegative | None) –

Per-rotation allowance obtained by rounding E minus the stored replacement_error downward, dividing that remaining allowance by M exactly, and rounding downward again. None when M is zero. The stored replacement error plus M times this allowance is at most E.
exact_t (Count) –

T and T-inverse occurrences of the emitted circuit, which need no synthesis and spend no budget.
t (ResourceLaw | None) –

The T estimate as a ResourceLaw (metric t, basis clifford_t, interpretation estimate, precision E), or None.
t_unavailable (Text | None) –

Why there is no estimate, for example a per-rotation budget of 1 or more, where the logarithmic model does not apply, or None.

QHDErrorSource ¶

Bases: Record

One source of error of a QHD circuit against its finite model, with its own status.

circuit_resources returns one per source in error_sources. The fields below are read-only. The kinetic_model source instead compares that finite model with the finite-difference one. A bound is an upper bound on the operator-norm change that the source causes, or on the state 2-norm change for the preparation. An estimate is a value without that guarantee, which description explains. A requested value is a budget the user asked for, never an achieved error. An unavailable source has no value and says why, and not_applicable names a source that this construction does not have. The sources add by telescoping only when each compares adjacent stages of one chain, and the record forms no combined total (circuit_resources lists the stages).

Attributes:

name (ErrorSourceName) –

The source, one of the stages that circuit_resources lists.
status (Literal['bound', 'estimate', 'requested', 'unavailable', 'not_applicable']) –

bound, estimate, requested, unavailable or not_applicable, as above.
value (Nonnegative | None) –

The value, None exactly for unavailable and not_applicable.
description (Text) –

What the value compares and on what it rests.

QHDEvolutionBound ¶

Bases: Record

Operator-norm bounds on the compiled product of one QHD Plan against its finite model.

evolution_bound returns it, and circuit_resources includes it as evolution. The fields below are read-only. The answer is evolution, a bound on the compiled product with the stored exponents against the exact time-ordered evolution of the finite model, with evolution_status saying whether it is a bound or conditional. The finite model is H(t) = a(t) T + b(t) V on the K**d grid states, with the exact schedule functions, T the Plan's kinetic operator on its grid and V the stored support tables plus the constant. The bounds say nothing about the symbolic objective between or at the grid points, the continuum problem or an optimization gap.

The norm inputs and the splitting, time-ordering and midpoint-quadrature values are upper bounds in the spectral norm on domain, evaluated outward from exact rationals of the binary64 inputs, so they are never below the real-arithmetic formulas. The coefficient residual is a bound or a first-order estimate, as coefficient_status says. A norm input is None when it exceeds the binary64 range, and the bounds it enters are then the cap 2.

Attributes:

domain (Literal['valid_one_hot_subspace', 'full_binary_register']) –

"valid_one_hot_subspace", the states with one excitation per variable register, or "full_binary_register", every state of the d b binary qubits after the permutation between circuit and lexicographic variable order, with the physical identity phase included.
reference (Text) –

The finite model: encoding, kinetic model, boundary, grid points, spacings, step count and time step, schedule, coefficient rule and order.
formula (Literal['one_hot_first_order', 'one_hot_second_order', 'binary_first_order', 'binary_second_order']) –

"one_hot_first_order" for the potential factor followed by every link in emitted order, "one_hot_second_order" for the symmetric potential, odd-link, even-link product, and "binary_first_order", "binary_second_order" for the potential and whole-kinetic factors of the binary encoding with exact QFTs.
commutator (Nonnegative | None) –

C, a bound on ||[T, V]||.
kinetic_nested (Nonnegative | None) –

D_T, a bound on ||[T, [T, V]]||.
potential_nested (Nonnegative | None) –

D_V, a bound on ||[V, [V, T]]||.
kinetic_norm (Nonnegative | None) –

mu_T, a bound on ||T||.
potential_norm (Nonnegative | None) –

mu_V, a bound on ||V|| including the constant.
hopping (Nonnegative | None) –

Gamma, the sequential link commutator sum of one-hot first order, None for second order and for the binary encoding.
even_nested (Nonnegative | None) –

J_E, a bound on ||[E, [E, O]]|| of the even and odd link layers of one-hot second order, None otherwise.
odd_nested (Nonnegative | None) –

J_O, a bound on ||[O, [O, E]]||, None as J_E is.
norm_methods (tuple[tuple[Text, Text], ...]) –

(input, method) pairs naming which bound attained each minimum: range for the table-range bounds C_0 and D_(T,0), neighbor for the neighbor-difference bounds C_edge and D_(V,edge), and commutator for 2 tau C (D_T) and 2 nu C (D_V). mu_T and mu_V come from the kinetic diagonal and the table extremes, Gamma from its closed form and J_E, J_O from the K-point graph.
splitting (Nonnegative) –

min(2, sum_k s_k), the product-formula bound.
time_ordering (Nonnegative) –

min(2, sum_k w_k), the time-ordering bound.
midpoint_quadrature (Nonnegative | None) –

min(2, sum_k q_k) under the midpoint rule, None under the integrated rule.
schedule (Nonnegative) –

min(2, sum_k (w_k + q_k)).
coefficient_residual (Nonnegative | None) –

min(2, sum_k e_coeff,k), or None when an assumption of its estimate is missing.
coefficient_status (Literal['bound', 'estimate', 'unavailable']) –

"bound" when every step's residual is an exact rational discrepancy, "estimate" when some step uses the first-order integral estimate, "unavailable" when the estimate's assumption fails.
coefficient_unavailable (Text | None) –

Why the coefficient residual is None, or None.
evolution (Nonnegative | None) –

min(2, splitting + schedule + coefficient_residual), for the compiled product with the stored exponents against the exact time-ordered evolution of the finite model: a bound when evolution_status is "bound", conditional on the coefficient estimate when it is "conditional", and None without the coefficient residual.
evolution_status (Literal['bound', 'conditional', 'unavailable']) –

"bound", "conditional" when it rests on the coefficient estimate, or "unavailable".
work (Count) –

Scalar work: table entries reduced, support-axis comparisons, graph rows and steps.

QHDRunResources ¶

Bases: Record

Rotation and T totals of an augmented-Lagrangian run or a box refinement over its circuits and shots.

run_resources returns it. The fields below are read-only. Two totals are kept apart. The body totals count every compiled circuit once (circuits of each entry), the static inventory of what a compiler synthesizes. The shot totals weight each circuit by the raw shots set aside for its attempts of every status, sum_e shots_e R_e, because each executed shot runs the whole circuit, state preparation included, and the shots set aside bound the executed ones from above. Exact readout uses no shots, while hardware would need some, so its shot totals are None and unavailable says to multiply each entry's per-circuit values by the intended shots. Per circuit the counts are those of circuit_resources. Each round and level uses its own tables, spacing and exponents, so its own circuit is counted rather than the first one multiplied. A total is None when an entry with a positive or unknown multiplicity has no per-circuit value, and unavailable names it. An entry with zero multiplicity contributes zero. A sum of T estimates is an estimate.

The rotation totals add the recorded rotation counts. For the binary encoding the count includes the rotation gates before angle classification, exact T gates and the omitted zero angles of dense diagonals included, so those totals are upper bounds, while the T estimates use the exact angle classification of circuit_resources.

Attributes:

synthesis_epsilon (Nonnegative | None) –

The per-circuit budget E of the T estimates, or None.
rotation_interpretation (Literal['exact', 'upper_bound']) –

"upper_bound" when the run uses the binary encoding, a live Plan's rotation count is an upper bound, or a round or level without a live result adds a positive recorded count with a positive or unknown multiplicity, because its record does not keep the interpretation of that count. "exact" otherwise.
entries (tuple[QHDRunEntry, ...]) –

One QHDRunEntry per started round, or per started level of a refinement, in order.
arbitrary_rotations (Count | None) –

sum_e circuits_e R_e, with R_e the entry's arbitrary rotations.
shot_arbitrary_rotations (Count | None) –

sum_e shots_e R_e, None under exact readout.
t_estimate (Nonnegative | None) –

sum_e circuits_e T_e, with T_e the entry's T estimate.
shot_t_estimate (Nonnegative | None) –

sum_e shots_e T_e, None under exact readout.
unavailable (tuple[tuple[Text, Text], ...]) –

(total, reason) for every total that is None.

QHDRunEntry ¶

Bases: Record

The circuit of one augmented-Lagrangian round or refinement level and its multiplicities.

run_resources returns one per started round or level in entries. The fields below are read-only.

Attributes:

label (Text) –

"round k", "level z" or, for a round with box refinement, "round k level z".
plan_id (ContentID | None) –

Content hash of the inner Plan, or None when no Plan was chosen or the record does not name one.
circuits (Count | None) –

Circuit preparations of its Run, the compiled bodies.
shots (Count | None) –

Raw shots set aside by its attempts of every status, which bound the shots it executed from above.
arbitrary_rotations (Count | None) –

Arbitrary rotations of its circuit, the value that its resources recorded, which for the binary encoding counts rotation gates before angle classification.
t_estimate (Nonnegative | None) –

T estimate of its circuit, as circuit_resources forms it.
unavailable (tuple[tuple[Text, Text], ...]) –

(field, reason) for every per-circuit value that is None.

Limits¶

QHD is not a global optimizer. The candidate and the most probable point are points of a finite grid of the box, and the candidate is the least evaluated objective only among the observed points with positive weight.
max_work and max_bytes limit the work and bytes that planning counts. They do not bound the process memory or undocumented SymPy, SciPy or Qiskit costs.
The error sources of circuit_resources compare the circuit with its finite model. None concerns grid discretization or the optimization gap, and the record forms no total.
The T estimate is the leading term of the typical Ross-Selinger count, without an error bar, and a requested synthesis budget is not an achieved error.
Limitations and open work lists the open work on the QHD models.

Entries on other pages¶

The augmented-Lagrangian layer and box refinement are documented on QHD constrained problems and box refinement: