Skip to content

QSP and QSVT

Compute symmetric quantum signal processing (QSP) phase factors for a real Chebyshev target, and build QSVT circuits that apply the polynomial to a block-encoded matrix: a real Chebyshev transform, the Hamiltonian evolution exp(-i t A), the 1/x polynomial of the quantum linear-system solver (QLS) and Dalzell's kernel reflection. The conventions define the normalization and phase signs used here. Import the functions and records from nwqlib.subroutines.qsp. Five entries import only from their submodule: chebyshev_norming_sup_bound from nwqlib.subroutines.qsp.phases, QSPPreparedEvolution, prepare_qsp_evolution and qsp_evolution_error_terms from nwqlib.subroutines.qsp.evolution, and plan_kernel_reflection from nwqlib.subroutines.qsp.shortcut.

import numpy as np
from qiskit.quantum_info import Operator
from nwqlib.subroutines.block_encoding import build_block_encoding
from nwqlib.subroutines.qsp import (
    build_real_chebyshev_encoding,
    solve_symmetric_qsp_phases,
)

solution = solve_symmetric_qsp_phases([0.0, 0.0, 0.0, 0.5])  # f = T_3 / 2
encoding = build_block_encoding(np.diag([0.6, -0.2]))  # alpha = 0.6
circuit = build_real_chebyshev_encoding(encoding, solution.phases)
stride = 2 ** (circuit.num_qubits - encoding.system_qubits)
block = Operator(circuit).data[::stride, ::stride]
print(np.round(block.real, 6))

[[0.5      0.      ]
 [0.       0.425926]]

A / alpha is diag(1, -1/3), and T_3(x)/2 is 0.5 at 1 and 23/54 = 0.425926... at -1/3. The ancilla qubits come first, so the block is read at the stride of their dimension.

Phase fitting

solve_symmetric_qsp_phases accepts finite, real, definite-parity Chebyshev targets with max |f| <= 1. The module text below states the Wx convention, the solver and their sources (Martyn, Rossi, Tan and Chuang, arXiv:2105.02859v5, Dong, Meng, Whaley and Lin, arXiv:2002.11649v2, and Dong, Lin, Ni and Wang, arXiv:2307.12468v1).

Wx-convention quantum signal processing: evaluation and phase factors.

[MRTC] is Martyn, Rossi, Tan, and Chuang, arXiv:2105.02859v5, [DMWL] is Dong, Meng, Whaley, and Lin, arXiv:2002.11649v2, and [DLNW] is Dong, Lin, Ni, and Wang, arXiv:2307.12468v1. Equation and section numbers refer to these arXiv versions, those of [DLNW] to the PDF of that version.

Convention (project standard, [MRTC] arXiv:2105.02859v5, Sec. II.A, Eqs. (1)-(3) and Theorem 1, and App. A.1, Eq. (A2)): the signal operator is W(x) = [[x, i sqrt(1-x^2)], [i sqrt(1-x^2), x]], phases apply as e^{i phi Z}, and a phase vector (phi_0, ..., phi_d) realizes

P(x) = <0| e^{i phi_0 Z} prod_{j=1..d} [ W(x) e^{i phi_j Z} ] |0>.

Conversion identities implemented here (each guarded by a unit test):

  • Reflection convention: with R(x) = [[x, sqrt(1-x^2)], [sqrt(1-x^2), -x]] ([MRTC] arXiv:2105.02859v5, Eq. (A4)) and W(x) = i e^{-i pi/4 Z} R(x) e^{-i pi/4 Z} (the inverse of [MRTC] arXiv:2105.02859v5, Eq. (14)), the Wx product equals i^d times the reflection product with phases phi'_0 = phi_0 - pi/4, phi'_d = phi_d - pi/4, and interior phi'_j = phi_j - pi/2. wx_phases_to_reflection returns those shifted phases plus the compensating global phase d * pi/2. [MRTC] arXiv:2105.02859v5, App. A.2, Eq. (A5), uses the same interior and final shifts but sets phi'_0 = phi_0 + (2d-1) pi/4, folding the i^d factor into the first phase. That form reproduces the <0|.|0> entry. Keeping the factor as a global phase makes the identity hold for the whole 2 x 2 product.
  • Qiskit sign: RZGate(theta) = exp(-i theta Z / 2) exactly, so applying e^{i phi Z} is RZGate(-2 phi) with no residual global phase. The QSVT circuit builders apply projector phases e^{i phi (2 Pi - 1)} as RZGate(+2 phi) on a signal qubit that is flipped onto the projector subspace, which is the same identity.

Phase factors follow the symmetric-QSP optimization of [DMWL] arXiv:2002.11649v2: symmetric phase vectors phi_j = phi_{d-j} (Sec. III.1, Theorem 2, and Eq. (24)), the mean-square residual of Re P(x) against the target polynomial on positive Chebyshev nodes (Sec. III.2, Eq. (23)), the initialization (pi/4, 0, ..., 0, pi/4) of Eq. (27), which realizes Re P = 0 and whose free phases (pi/4, 0, ..., 0) are Eq. (28) (Sec. III.4), and L-BFGS descent (Sec. III.5, Algorithm 1). This module departs from the paper in four places. The objective uses four times as many nodes as free phases (QSP_SOLVER_GRID_MULTIPLIER). A damped Gauss-Newton polish supplies the terminal digits after L-BFGS. If L-BFGS from the Eq. (27) point leaves an objective-node residual at or above the tolerance, the same damped iteration runs from that point, which is the Newton method of [DLNW] arXiv:2307.12468v1, Sec. 3, Eq. (3.1) and Algorithm 3.1. Their zero reduced phases are the [DMWL] start written in an imaginary-part convention ([DLNW] Sec. 2.2 and Eq. (2.3)). Acceptance is decided on a separate dense verification grid with the Chebyshev norming bound rather than on the objective.

The Newton start exists for targets whose sup norm lies near 1. The larger parity target of the Hamiltonian evolution of build_qsp_evolution_encoding has a sup norm near 1/(1 + margin) with a margin of 1e-3, and every QLS phase target has the norming bound 1/(1 + QLS_TARGET_MARGIN), with the same margin. [DMWL] arXiv:2002.11649v2, Sec. IV.1, divides the Jacobi-Anger series by 2, so their Hamiltonian-simulation tests of the fixed start stay at max|f| <= 1/2, and Sec. IV.5 reports a Hessian condition number that grows as max|f| approaches 1. [DLNW] arXiv:2307.12468v1, Sec. 2.2, states that optimization methods such as L-BFGS lose their convergence guarantee near max|f| = 1, and Sec. 4.3 (Fig. 8) reports Newton's method from the same start converging for 1 - max|f| down to 1e-9. The engineering constants record the NWQLib sweeps behind this order, for the evolution targets and for the QLS targets in the QLS_TARGET_MARGIN row. In both families L-BFGS missed the tolerance on some sampled targets, and the Newton start converged on each of them.

L-BFGS runs in scipy.optimize and each Newton step is one numpy.linalg.lstsq solve. Phase-factor computation is classical preprocessing of poly(d) cost, independent of any Hilbert-space dimension.

SymmetricQSPPhases 

SymmetricQSPPhases(*, phases: tuple[float, ...], degree: int, parity: int, target_chebyshev_coefficients: tuple[float, ...], max_residual: float, residual_sup_bound: float, evaluations: int)

Solved symmetric QSP phase factors with their residual and its sup-norm bound.

solve_symmetric_qsp_phases returns it. Its phases field holds the result. The fields below are read-only.

Attributes:

  • phases (tuple[float, ...]) –

    Full symmetric Wx phase vector (phi_0, ..., phi_d).

  • degree (int) –

    Polynomial degree d.

  • parity (int) –

    Target parity (0 even, 1 odd).

  • target_chebyshev_coefficients (tuple[float, ...]) –

    Chebyshev coefficients (c_0, ..., c_d) of the target f that the phases fit.

  • max_residual (float) –

    Largest |Re P - f| on the dense verification grid.

  • residual_sup_bound (float) –

    Chebyshev-grid norming bound on sup_{x in [-1,1]} |Re P(x) - f(x)|.

  • evaluations (int) –

    Number of objective and Jacobian evaluations and final verification calls made, across both starts, polish steps and line searches.

to_dict 

to_dict() -> dict[str, Any]

Return a JSON-like phase-solution record.

chebyshev_grid 

chebyshev_grid(num_points: int) -> ndarray

Return the Chebyshev grid x_j = cos(pi (j + 1/2) / M), j = 0, ..., M - 1, with M = num_points.

Parameters:

  • num_points (int) –

    Positive number M of nodes.

Returns:

  • grid ( ndarray ) –

    The M nodes in index order.

Raises:

  • ValueError –

    If num_points is not positive.

chebyshev_norming_sup_bound 

chebyshev_norming_sup_bound(grid_max: float, *, degree: int, num_points: int) -> float

Bound the sup norm of a polynomial of degree degree from its maximum on a Chebyshev grid.

Returns grid_max / cos(pi d / (2 N)), an upper bound on sup_{x in [-1, 1]} |p(x)| when grid_max is the largest |p| on the N nodes of chebyshev_grid(N). This is the Chebyshev-node norming bound of Ehlich and Zeller (1964), doi:10.1007/BF01111276, Satz 2, Eqs. (12)-(14), pp. 42-43, which use exactly these Chebyshev-zero nodes. The extension to complex coefficients follows by rotating a maximum value to the real axis and applying the real-polynomial inequality. Sunderhauf et al., arXiv:2507.15537v1 Eq. (25), state this factor for equidistant x points, where it is not valid. Equispaced angles give the required nodes x_j=cos(theta_j), not an equispaced x grid. Ehlich and Zeller, doi:10.1007/BF01111276, Satz 1, treat equidistant x nodes separately. The QLS guide gives an exact-rational counterexample to the Eq. (25) extension.

num_points is the number N of nodes on the whole interval [-1, 1], as opposed to the positive-half grid of the phase solver. The inequality requires strictly more nodes than the polynomial degree. The result is an exact-arithmetic bound evaluated in binary64, not an interval-arithmetic enclosure.

chebyshev_polynomial_sup_bound 

chebyshev_polynomial_sup_bound(chebyshev_coefficients: Any, *, max_degree: int = 256, max_bytes: int = DEFAULT_INPUT_BYTES) -> float

Return an upper bound on sup_{x in [-1, 1]} |f(x)| for f = sum_k c_k T_k.

The series is sampled on N = 64 (d + 1) Chebyshev zeros, and the sampled maximum is divided by cos(pi d / (2N)) as in chebyshev_norming_sup_bound. The QLS rescale, the evolution target scale and the input check of solve_symmetric_qsp_phases all use this bound. max_degree (default 256) and max_bytes (default 10 GB) are checked before the grid is formed.

Parameters:

  • chebyshev_coefficients (array_like) –

    Chebyshev coefficients (c_0, ..., c_d) of f.

  • max_degree (int, default: 256 ) –

    Default 256. Largest accepted degree d.

  • max_bytes (int, default: DEFAULT_INPUT_BYTES ) –

    Default 10 GB (decimal, 10_000_000_000 bytes). Limit on the grid arrays.

Returns:

  • bound ( float ) –

    The upper bound on sup_{x in [-1, 1]} |f(x)|.

evaluate_qsp_polynomial 

evaluate_qsp_polynomial(phases: Any, x_values: Any) -> ndarray

Evaluate the Wx-convention QSP polynomial P(x) on a grid.

Parameters:

  • phases (Any) –

    Phase vector (phi_0, ..., phi_d).

  • x_values (Any) –

    Signal values in [-1, 1].

Returns:

  • ndarray –

    Complex P(x) values, one per grid point.

wx_phases_to_reflection 

wx_phases_to_reflection(phases: Any) -> tuple[ndarray, float]

Convert Wx phases to reflection-convention phases ([MRTC] arXiv:2105.02859v5, Eq. (14)).

The shifts are those of [MRTC] arXiv:2105.02859v5, App. A.2, Eq. (A5), except that the i^d factor is returned as a global phase (see the conversion identities in the module description).

Parameters:

  • phases (array_like) –

    Wx phase vector (phi_0, ..., phi_d) with d >= 1.

Returns:

  • reflection_phases ( ndarray ) –

    The reflection-convention phases.

  • global_phase ( float ) –

    global_phase = d * pi / 2 compensates the i^d factor, so that the reflection-form product times e^{i global_phase} equals the Wx-form product.

Raises:

  • ValueError –

    If fewer than two phases are given.

solve_symmetric_qsp_phases 

solve_symmetric_qsp_phases(chebyshev_coefficients: Any, *, residual_tolerance: float = QSP_SOLVER_RESIDUAL_TOLERANCE, max_degree: int = 256, max_evaluations: int = 20000, max_bytes: int = DEFAULT_INPUT_BYTES, limit_name: str = 'max_evaluations') -> SymmetricQSPPhases

Solve symmetric Wx phase factors whose Re P(x) matches a real definite-parity Chebyshev target.

Implements the [DMWL] optimization (arXiv:2002.11649v2, Sec. III): symmetric phases, positive Chebyshev node objective, (pi/4, 0, ..., 0, pi/4) initialization and L-BFGS, followed by a damped Gauss-Newton polish. If the objective-node residual from that start is at or above residual_tolerance or nonfinite, the damped Newton iteration of [DLNW] arXiv:2307.12468v1 (Sec. 3, Eq. (3.1) and Algorithm 3.1) runs from the same initialization, and the better of the two results is verified. Near max|f| = 1, where the evolution and QLS targets lie, L-BFGS from that point has no convergence guarantee ([DLNW] Sec. 2.2). The returned phases satisfy Re P(x) ~= f(x) in the Wx convention. residual_sup_bound bounds sup |Re P - f| over [-1, 1] by the Chebyshev norming inequality. The degree and the known simultaneous numerical arrays are checked before optimization.

Parameters:

  • chebyshev_coefficients (array_like) –

    Finite real Chebyshev coefficients (c_0, ..., c_d) of the target f. Entries of the wrong parity must be zero, and max |f| must not exceed 1.

  • residual_tolerance (float, default: QSP_SOLVER_RESIDUAL_TOLERANCE ) –

    Default 1e-12. Largest acceptable |Re P - f| on the dense verification grid.

  • max_degree (int, default: 256 ) –

    Default 256. Largest accepted target degree.

  • max_evaluations (int, default: 20000 ) –

    Default 20000. Limit on the combined objective and Jacobian evaluations and final verification calls across both starts, polish steps and line searches.

  • max_bytes (int, default: DEFAULT_INPUT_BYTES ) –

    Default 10 GB (decimal, 10_000_000_000 bytes). Bound on known simultaneous numerical arrays. It does not measure process memory, the internal workspace of external libraries or elapsed time.

  • limit_name (str, default: 'max_evaluations' ) –

    Default "max_evaluations". Name of the caller's option that sets max_evaluations, used in the error messages.

Returns:

  • solution ( SymmetricQSPPhases ) –

    The full phase vector in solution.phases, with the residual and its sup-norm bound.

Raises:

  • ValueError –

    For empty, mixed-parity or out-of-range targets, or when the solver cannot reach residual_tolerance (with diagnostics).

Examples:

For f(x) = x / 2, a degree-one target, both symmetric phases are -pi/6 = -0.5235987756..., and Re P(0.3) = 0.15:

>>> import numpy as np
>>> from nwqlib.subroutines.qsp import (
...     evaluate_qsp_polynomial, solve_symmetric_qsp_phases)
>>> solution = solve_symmetric_qsp_phases([0.0, 0.5])
>>> print([round(phase, 10) for phase in solution.phases])
[-0.5235987756, -0.5235987756]
>>> values = evaluate_qsp_polynomial(solution.phases, [0.3, -0.8])
>>> print(np.round(values.real, 10))
[ 0.15 -0.4 ]

Hamiltonian evolution and circuit construction

build_qsvt_circuit and build_real_chebyshev_encoding apply given Wx phases to a block encoding. build_qsp_evolution_encoding builds a block encoding of exp(-i t A) from one of a Hermitian A. jacobi_anger_expansion computes its polynomial without Qiskit, for finite real tau>0 and 0<epsilon<1, with degree at most max_degree. Within the search range H=min(ceil(tau),max_degree)+200, it chooses the smallest degree whose complete tail bound meets epsilon, and breaks ties by the smaller tail and then by the earlier last order of the finite Bessel sum. The real-argument remainder bound and an analytic exclusion leave at most 202 candidates for that last order, and the Bessel values are computed once through order H. The tail bound includes the omitted finite orders and the infinite suffix, with 2R_T per parity. If the analytic remainders of all passing candidates underflow, the first passing polynomial is returned with None tail bounds and slack, and its degree is not claimed to be the smallest certified one. The search raises an error when no degree within its finite range meets epsilon. The limit of 96*(H+1) bytes on the known arrays and NumPy's check of addressable array size are applied before allocation. The bound is evaluated in ordinary binary64 arithmetic. prepare_qsp_evolution solves the phases within one evaluation limit, and qsp_evolution_error_terms gives the error terms without a circuit. The assumptions and limits below apply to the builders.

QSVT circuits and Jacobi-Anger Hamiltonian-evolution synthesis.

[GSLW] is Gilyen, Su, Low, and Wiebe, arXiv:1806.01838v1, and [MRTC] is Martyn, Rossi, Tan, and Chuang, arXiv:2105.02859v5. Equation, lemma and theorem numbers refer to these arXiv versions.

Implements Jacobi-Anger evolution using the shared block-encoding records and the block-encoding and QSP conventions:

  • build_qsvt_circuit: projector-controlled phase QSVT ([MRTC] arXiv:2105.02859v5, Sec. II.C-II.D, Eq. (27), Fig. 3, Theorems 3-4, and [GSLW] arXiv:1806.01838v1, Theorem 17 and Lemma 19): alternate U_BE / U_BE^dagger with e^{i phi (2 Pi - 1)} phases, where Pi projects onto the all-zero block-encoding ancillas. The phase operator is implemented by conjugation with X. A signal qubit is flipped onto the Pi subspace (X controlled on the ancillas being all zero), rotated by RZGate(+2 phi) (= e^{-i phi Z}, the Qiskit sign RZGate(theta) = exp(-i theta Z / 2)), and flipped back. Wx phases, in the convention dictionary of [MRTC] arXiv:2105.02859v5, App. A, are converted with wx_phases_to_reflection and the i^d factor becomes circuit global phase, so the encoded block is exactly the Wx polynomial P(A / alpha).
  • build_real_chebyshev_encoding: block-encodes Re P(A / alpha) by averaging the +Phi and -Phi passes over one pair ancilla ([GSLW] arXiv:1806.01838v1, Corollary 18, Eq. (33)). Both branches share the same ladder of alternating U_BE and U_BE^dagger applications, because a CNOT from the pair qubit onto the signal qubit conjugates every RZ and so negates the reflection phases on the pair-one branch (X RZ(theta) X = RZ(-theta)). A Z on the pair qubit for odd degree fixes the (-1)^d of i^d <0|U_R(-Phi')|0> = (-1)^d conj(P). The shared ladder is NWQLib's circuit for Corollary 18's controlled U_Phi / U_-Phi pair.
  • build_qsp_evolution_encoding: Jacobi-Anger split ([GSLW] arXiv:1806.01838v1, Lemma 57, with the Bessel tails of Eqs. (53)-(56)) cos(tau x) = J_0(tau) + 2 sum (-1)^k J_2k(tau) T_2k(x) (even) and sin(tau x) = 2 sum (-1)^k J_2k+1(tau) T_2k+1(x) (odd), truncation degree from the Bessel tail bound 2 sum_{k>d} |J_k(tau)| <= eps with the recorded slack, one parity ancilla combining the passes with coefficients (1/2, -i/2) (the construction in the proof of [GSLW] arXiv:1806.01838v1, Theorem 58), and 3-step oblivious amplitude amplification A R A^dagger R A ([GSLW] arXiv:1806.01838v1, Theorem 28 with n = 3), which realizes -(3B - 4 B B^dagger B), where B is the all-zero-ancilla block before amplification. The circuit adds a compensating pi global phase. At block amplitude a = 1/(2s) the OAA output (3a - 4a^3) is quadratically insensitive to the recorded target rescale s, so the margin costs only the recorded amplitude deficit.
  • build_control_diagonal_generator_encoding: block encoding of the effective Hamiltonian of the linear combination of Hamiltonian simulation (LCHS) method, adapted from Pocrnic, Johnson, Katabarwa, and Wiebe, arXiv:2506.20760v2, Lemma 7 (p. 15).

The qubitized walk of Low and Chuang, arXiv:1610.06546v3, is an alternative route to Hamiltonian evolution that this module does not implement.

JacobiAngerExpansion 

JacobiAngerExpansion(*, tau: float, epsilon: float, degree: int, cos_degree: int, sin_degree: int, cos_coefficients: tuple[float, ...], sin_coefficients: tuple[float, ...], tail_bound: float | None, tail_slack: float | None, cos_tail_bound: float | None, sin_tail_bound: float | None, analysis_terminal: int, analytic_infinite_tail_bound: float | None)

Truncated Jacobi-Anger expansion of exp(-i tau x) with its Bessel tail bounds.

jacobi_anger_expansion returns it. The polynomial is in cos_coefficients and sin_coefficients, and its truncation error bound is tail_bound. The fields below are read-only.

Attributes:

  • tau (float) –

    Effective evolution time alpha * t.

  • epsilon (float) –

    Requested truncation tolerance for the combined expansion.

  • degree (int) –

    Smallest degree at or above the supplied min_degree whose complete tail bound, with a positive analytic remainder, meets epsilon within the finite search range of jacobi_anger_expansion. Ties choose the smaller computed complete tail, then the earlier analysis_terminal. When the analytic remainder underflows, the returned polynomial has None tail bounds, and its degree is not claimed to be the smallest certified one.

  • cos_degree (int) –

    Even-part polynomial degree.

  • sin_degree (int) –

    Odd-part polynomial degree.

  • cos_coefficients (tuple[float, ...]) –

    Chebyshev coefficients of the truncated cosine.

  • sin_coefficients (tuple[float, ...]) –

    Chebyshev coefficients of the truncated sine.

  • tail_bound (float | None) –

    Upper bound on the combined Bessel tail 2 sum_{k>degree} |J_k|, computed as the finite sum through analysis_terminal plus analytic_infinite_tail_bound.

  • tail_slack (float | None) –

    epsilon / tail_bound, the margin to the truncation boundary. A degree selected with slack near one can change under a different Bessel-function implementation, so recorded integer degrees are meaningful only together with this slack.

  • cos_tail_bound (float | None) –

    Complete even-part finite-plus-infinite tail bound.

  • sin_tail_bound (float | None) –

    Complete odd-part finite-plus-infinite tail bound.

  • analysis_terminal (int) –

    Last Bessel order included in the finite sums used for the reported tail bounds.

  • analytic_infinite_tail_bound (float | None) –

    Combined analytic infinite-tail term added across both parity bounds, or None when that mathematically nonzero term is not representable in binary64. On that path, all tail bounds and tail_slack are also None.

to_dict 

to_dict() -> dict[str, Any]

Return a JSON-like expansion record.

QSPPreparedEvolution 

QSPPreparedEvolution(*, expansion: JacobiAngerExpansion, cos_solution: SymmetricQSPPhases, sin_solution: SymmetricQSPPhases, scale: float, target_norming_bounds: tuple[float, float], margin: float, attempted_margins: tuple[float, ...])

The solved phases of both parity parts of a QSP evolution, with the common target scale.

prepare_qsp_evolution returns it, and qsp_evolution_error_terms reads it. It holds numbers only, no circuit. The coefficient and phase tuples are immutable. The fields below are read-only.

Attributes:

  • expansion (JacobiAngerExpansion) –

    The Jacobi-Anger truncation whose parity targets were solved.

  • cos_solution (SymmetricQSPPhases) –

    Phases realizing Re P = cos-series / scale.

  • sin_solution (SymmetricQSPPhases) –

    Phases realizing Re P = sin-series / scale.

  • scale (float) –

    Common target divisor s = max(1, cos bound, sin bound) * (1 + margin) >= 1.

  • target_norming_bounds (tuple[float, float]) –

    Norming bounds on the cos and sin series before division by scale.

  • margin (float) –

    The target margin, from QSP_EVOLUTION_TARGET_MARGINS = (1e-3, 3e-3), for which both phase solves converged.

  • attempted_margins (tuple[float, ...]) –

    Every margin tried, in order, including failures.

recovery_scale 

recovery_scale

Return 1 / (3a - 4a^3) with a = 1/(2 scale).

Multiplying recovered postselected amplitudes by this value removes the known OAA amplitude deficit.

jacobi_anger_expansion 

jacobi_anger_expansion(tau: float, epsilon: float, *, min_degree: int = 1, max_degree: int = 256, max_bytes: int = DEFAULT_INPUT_BYTES) -> JacobiAngerExpansion

Truncate the Jacobi-Anger expansion by the Bessel tail bound.

exp(-i tau x) = cos(tau x) - i sin(tau x) with cos(tau x) = J_0 + 2 sum_{k>=1} (-1)^k J_2k T_2k and sin(tau x) = 2 sum_{k>=0} (-1)^k J_2k+1 T_2k+1 ([GSLW] Lemma 57, arXiv:1806.01838v1, pp. 49-50). Truncating both parities at k <= d leaves error at most 2 sum_{k>d} |J_k(tau)| because |T_k| <= 1 on [-1, 1] (Eqs. (53)-(54)). The bound adds the Bessel magnitudes through the selected terminal T, the last order included in its finite sum, and bounds the rest by the real-argument analytic remainder R_T=(tau/2)^(T+1)/((T+1)!*(1-tau/(2(T+2)))). Its derivation from [GSLW] Eq. (55), DLMF 10.9.4 and a geometric suffix is in _analytic_jacobi_anger_remainder. For eligible degrees min_degree <= d <= min(max_degree, T - 1), the complete combined bound is C(d,T)=2 sum_{k=d+1..T} |J_k(tau)| + 4R_T. Each parity tail contains its omitted finite orders and 2R_T for the whole analytic suffix, so their sum conservatively counts the suffix twice.

Within the search range H=min(ceil(tau),max_degree)+200, the candidates are the pairs (d,T) with a positive remainder and finite C(d,T)<=epsilon, ordered by (d,C(d,T),T), that is, smallest degree first, then smallest complete tail, then earlier terminal. For T>=3, T<=tau implies R_T>1/4 whenever q<1, so such a terminal cannot meet epsilon<1. Only T=2 needs a separate candidate when min_degree=1 and tau<3. The proof, given in the source of this function, leaves at most 202 terminals. The minimum is taken over these candidates within H, evaluated with ordinary binary64 sums and special functions, and is not a global or outward-rounded minimum.

If every passing remainder underflows to zero, the function returns the polynomial of the first passing terminal and its first passing degree, with all tail bounds and tail_slack set to None. A candidate with a positive remainder takes priority over this fallback. If neither exists, the search raises ValueError. The scaling of [GSLW] arXiv:1806.01838v1, Cor. 60, d = Theta(tau + log(1/eps)/log(e + log(1/eps)/tau)), is the reference form of the resulting degree. min_degree is a lower limit on the returned degree. The phase solver needs min_degree=2 so the even (cos) part is never a constant. max_degree bounds the returned polynomial degree. The Bessel magnitude table is evaluated once through H, at most max_degree + 201 orders, and the fixed limit on candidate terminals keeps the work of the suffix sums at O(202*(H+1)). Scalar checks of the remainder reject a search that cannot succeed before the Bessel values are allocated. max_bytes bounds the known numerical arrays, 96*(H+1) bytes, separately from NumPy's limit on addressable array size. It excludes bounded Python bookkeeping, the workspace of the special functions and process memory.

Parameters:

  • tau (float) –

    Finite positive real effective evolution time alpha * t.

  • epsilon (float) –

    Combined truncation tolerance, strictly between 0 and 1.

  • min_degree (int, default: 1 ) –

    Default 1. Smallest returned degree.

  • max_degree (int, default: 256 ) –

    Default 256. Largest returned degree.

  • max_bytes (int, default: DEFAULT_INPUT_BYTES ) –

    Default 10 GB (decimal, 10_000_000_000 bytes). Limit on the known numerical arrays.

Returns:

  • expansion ( JacobiAngerExpansion ) –

    Cosine and sine Chebyshev coefficients with complete tail bounds and tail_slack, or a polynomial with None tail evidence when the suffix underflows.

Raises:

  • ValueError –

    If tau, epsilon or a degree or byte limit is invalid, the arrays exceed the byte or size limit, or no eligible degree meets epsilon inside the bounded search.

build_qsvt_circuit 

build_qsvt_circuit(encoding: BlockEncoding, wx_phases: Any) -> QuantumCircuit

Build the projector-controlled-phase QSVT circuit for Wx phases.

The all-zero block (signal + block-encoding ancillas) equals the Wx polynomial P(A / alpha) exactly (for Hermitian encoded blocks and definite-parity-compatible phases), as in [MRTC] arXiv:2105.02859v5, Theorem 3. Trivial (all-zero) Wx phases realize the Chebyshev polynomial T_d(A / alpha) ([MRTC] arXiv:2105.02859v5, Sec. II.A, p. 3). That known answer checks the phase, sign and ordering conventions of this module at circuit level.

Parameters:

  • encoding (BlockEncoding) –

    Block encoding of the operator to transform, including its normalization and ancilla layout.

  • wx_phases (Any) –

    Wx-convention phase vector (phi_0, ..., phi_d).

Returns:

  • circuit ( QuantumCircuit ) –

    Circuit on 1 + num_ancillas + system_qubits qubits with the signal qubit first (register order: signal, block ancillas, system).

build_real_chebyshev_encoding 

build_real_chebyshev_encoding(encoding: BlockEncoding, wx_phases: Any) -> QuantumCircuit

Block-encode Re P(A / alpha) with the +Phi and -Phi passes on one shared ladder.

This is [GSLW] arXiv:1806.01838v1, Corollary 18, Eq. (33): Re P = (P + P^*)/2, and the phases -Phi realize P^*. One pair ancilla in H .. H averages the +Phi and -Phi QSVT passes. Both branches share the U_BE ladder because a CNOT from the pair qubit negates every reflection phase on the pair-one branch, and a Z on the pair qubit for odd degree absorbs the (-1)^d of the conjugated pass. Register order: pair, signal, block ancillas, system.

Parameters:

  • encoding (BlockEncoding) –

    Block encoding of the operator to transform, as for build_qsvt_circuit.

  • wx_phases (array_like) –

    Wx-convention phase vector (phi_0, ..., phi_d).

Returns:

  • circuit ( QuantumCircuit ) –

    The block encoding of Re P(A / alpha), with register order pair, signal, block ancillas, system.

build_control_diagonal_generator_encoding 

build_control_diagonal_generator_encoding(encoding_l: BlockEncoding, encoding_h: BlockEncoding | None, *, l_diagonal: Any, h_diagonal: Any | None = None, max_work: int = 1000000000, max_bytes: int = DEFAULT_INPUT_BYTES, dense_control_route: str = 'auto') -> BlockEncoding

Block-encode the joint generator D_L (x) L + D_H (x) H.

Following the effective-Hamiltonian SELECT route of Pocrnic et al., arXiv:2506.20760v2, Lemma 7 (p. 15), adapted to the library's separate BE(L)/BE(H) children, the LCHS node coefficients (k, 1) become real diagonals D_L, D_H on a control register that joins the encoded operator's system space. Length-one diagonals give the scalar generator k L + H. Each child branch is a one-ancilla multiplexed-RY encoding of its rescaled diagonal tensored with the child encoding, and one LCU combine qubit sums the branches, so

alpha = max|D_L| * alpha_L + max|D_H| * alpha_H.

A child with fewer ancillas uses a prefix of the shared ancilla register, and the idle ancillas keep its all-zero block. Encoded system order: control register first (low-order), then the child system, so the encoded matrix is kron(L, D_L) + kron(H, D_H) in numpy index convention. Zero diagonal entries make their control state evolve trivially under the subsequent Hamiltonian simulation (identity padding).

Parameters:

  • encoding_l (BlockEncoding) –

    Block encoding of the Hermitian part L.

  • encoding_h (BlockEncoding | None) –

    Block encoding of H, or None when H = 0.

  • l_diagonal (Any) –

    Real diagonal for the L branch, power-of-two length 2**control_qubits indexed by the control-register basis.

  • h_diagonal (Any | None, default: None ) –

    Real diagonal for the H branch (same length), required exactly when encoding_h is supplied.

  • max_work (int, default: 1000000000 ) –

    Default 1_000_000_000. Limit on the work of the exact dense syntheses that controlling the two branches makes, and of Qiskit's control of the synthesized gates, in the work units of Exact dense synthesis. A dense_dilation child is synthesized and controlled once. A single active branch is not controlled here.

  • max_bytes (int, default: DEFAULT_INPUT_BYTES ) –

    Default 10 GB (decimal, 10_000_000_000 bytes). Limit on the working and kept bytes of those syntheses and control steps.

  • dense_control_route (str, default: 'auto' ) –

    Default "auto". "gatewise", "whole_matrix" or "auto", the route by which the combine qubit controls a dense_dilation child (dense control route). "auto" takes the whole-matrix route for this one control. The child's dilation is then synthesized with the control, and the branch's diagonal rotation is controlled gate-wise. A child whose circuit is not one dense unitary is controlled gate-wise on either route.

Returns:

  • encoding ( BlockEncoding ) –

    BlockEncoding of the joint generator, with implementation "control_diagonal_generator_lcu". An all-zero diagonal on one branch degenerates to the other branch alone.

prepare_qsp_evolution 

prepare_qsp_evolution(*, tau, epsilon, expansion=None, max_degree=256, max_evaluations=20000, max_bytes=DEFAULT_INPUT_BYTES, limit_name='max_evaluations')

Solve the QSP phases of both parity parts of a Jacobi-Anger expansion within one evaluation limit.

Both parity targets are divided by one common scale s = max(1, cos bound, sin bound) * (1 + margin), so the parity passes keep the relative weights of the Jacobi-Anger expansion. Margins are tried in the order of QSP_EVOLUTION_TARGET_MARGINS, 1e-3 then 3e-3, and the first one for which both phase solves converge is kept. Only a phase-convergence failure moves to the next margin. max_evaluations covers both parities and every attempted target margin, including failed solves and both starts of each solve. Sup-norm and Bessel preprocessing are bounded by degree and known numerical-array bytes. A supplied expansion is used directly, without another Bessel selection. This function and jacobi_anger_expansion do not import Qiskit.

Parameters:

  • tau (float) –

    Positive effective evolution time alpha * t.

  • epsilon (float) –

    Truncation tolerance, strictly between 0 and 1.

  • expansion (JacobiAngerExpansion | None, default: None ) –

    Default None, which computes jacobi_anger_expansion(tau, epsilon, min_degree=2). A supplied expansion must have the same tau and epsilon.

  • max_degree (int, default: 256 ) –

    Default 256. Largest accepted degree.

  • max_evaluations (int, default: 20000 ) –

    Default 20000. Limit on the phase-solver evaluations summed over both parities, both starts of each solve and every attempted margin.

  • max_bytes (int, default: DEFAULT_INPUT_BYTES ) –

    Default 10 GB (decimal, 10_000_000_000 bytes). Limit on known numerical arrays.

  • limit_name (str, default: 'max_evaluations' ) –

    Default "max_evaluations". Name of the caller's option that sets max_evaluations, used in the error messages.

Returns:

Raises:

  • ValueError –

    If tau <= 0, epsilon is not in (0, 1), the degree exceeds max_degree, the evaluations are exhausted, or the phases do not converge at any registered margin.

qsp_evolution_error_terms 

qsp_evolution_error_terms(prepared, *, evolution_time, child_error_bound)

Return the error terms of a prepared QSP evolution without building a circuit.

NWQLib's derivation. Before amplification, the parity qubit combines the cosine and sine passes with coefficients 1/2 and -i/2. The resulting half block, named for its amplitude a = 1/(2s) <= 1/2, is B = a (V + E) with V = exp(-i t A), the target at effective time tau = alpha t. The polynomial part obeys ||E|| <= tail_c + tail_s + s (resid_c + resid_s), where tail_c and tail_s are the parity Bessel tail bounds and resid_c and resid_s the phase-solver residual bounds. The child encoding represents A only within child_error_bound. For Hermitian A and A', [GSLW] arXiv:1806.01838v1, Lemma 61, gives ||exp(-i t A) - exp(-i t A')|| <= abs(t) ||A - A'||. That physical-time term is added to e_B = a ||E|| without the factor a, which is conservative because a <= 1/2. OAA maps B to 3B - 4 B B^dagger B. Expanding around aV gives ||block - (3a - 4a^3) V|| <= 3 e_B (1 + 4a^2) + 12 a e_B^2 + 4 e_B^3, which is at most 6 e_B + 6 e_B^2 + 4 e_B^3 for a <= 1/2.

error_bound is the raw BlockEncoding bound against exp(-i t A): the amplitude deficit 1 - (3a - 4a^3) plus that residual. Code that divides recovered amplitudes by the known factor 3a - 4a^3 uses compensated_recovery_error_bound, the residual bound times recovery_scale, instead. A missing child bound or an unrepresentable Bessel tail makes the dependent terms None.

Parameters:

  • prepared (QSPPreparedEvolution) –

    Output of prepare_qsp_evolution.

  • evolution_time (float) –

    Physical evolution time t.

  • child_error_bound (float | None) –

    Operator-norm error bound of the child encoding of A, or None when unknown.

Returns:

  • terms ( dict ) –

    amplitude (a), oaa_amplitude_factor (3a - 4a^3), deficit, polynomial_part_error, child_error, oaa_residual_error_bound, recovery_scale, compensated_recovery_error_bound and error_bound.

build_qsp_evolution_encoding 

build_qsp_evolution_encoding(encoding: BlockEncoding, *, evolution_time: float, epsilon: float, require_hermitian_evidence: bool = False, max_work: int = 1000000000, max_bytes: int = DEFAULT_INPUT_BYTES, _prepared: QSPPreparedEvolution | None = None) -> BlockEncoding

Synthesize a block encoding of exp(-i t A) from BE(A).

Follows the proof of [GSLW] arXiv:1806.01838v1, Theorem 58: Jacobi-Anger parity split (Lemma 57) at effective time tau = alpha * t, one real-Chebyshev pass per parity (each a shared ladder over the +-Phi pair), a parity ancilla combining them with coefficients (1/2, -i/2), and 3-step OAA (Theorem 28, n = 3) raising the amplitude from a = 1/(2s) to 3a - 4a^3, close to one because s is close to one. The error terms are derived in qsp_evolution_error_terms.

GSLW, arXiv:1806.01838v1, Lemma 61 (p. 53), bounds the change in evolution by abs(t) * ||A - A_encoded|| when both generators are Hermitian. The QSVT polynomial calculus also needs a Hermitian encoded block. For an externally constructed encoding, the caller can record these separately as target_operator_is_hermitian and encoded_operator_is_hermitian after establishing them. A dense numerical check is left to the caller, in the representation of the source and within its resource budget, because a circuit matrix needs memory exponential in the circuit width.

Dense Hermitian evidence uses exact entry equality, so even an asymmetry from rounding during assembly is rejected. If the intended target is Hermitian, the caller can explicitly choose (A + A.conj().T)/2 before encoding. That changes the supplied target and belongs to the caller's input model. The check never performs this projection implicitly.

Parameters:

  • encoding (BlockEncoding) –

    Block encoding of a Hermitian operator A (the QSVT eigenvalue calculus assumes Hermitian encoded blocks).

  • evolution_time (float) –

    Evolution time t > 0.

  • epsilon (float) –

    Jacobi-Anger truncation tolerance for exp(-i tau x).

  • require_hermitian_evidence (bool, default: False ) –

    Require construction metadata declaring both A and the encoded block Hermitian. By default, missing evidence is treated as an assumption that the caller has checked, recorded in the output with a warning. Explicitly non-Hermitian inputs are rejected in either mode. This option reads metadata only and never forms a circuit matrix or applies an operator to a state.

  • max_work (int, default: 1000000000 ) –

    Default 1_000_000_000, the default max_work of the block-encoding builders. Limit on the work of the exact dense syntheses that controlling the two passes makes, and of Qiskit's control of the synthesized gates, in the work units of Exact dense synthesis. A child that holds a dense UnitaryGate, such as a dense_dilation encoding, has that unitary synthesized once for each pass, and its adjoint once for each pass of degree two or more, and the parity control unrolls the synthesized gates at every query of the pass. The parity control takes the gate-wise route whatever route the joint generator uses, because its controlled object is the whole pass. A gate that an earlier control already synthesized, such as a branch of a two-child joint generator, is controlled again outside this count, and LCHS planning counts it.

  • max_bytes (int, default: DEFAULT_INPUT_BYTES ) –

    Default 10 GB (decimal, 10_000_000_000 bytes). Limit on the working and kept bytes of those syntheses and control steps.

  • _prepared (QSPPreparedEvolution | None, default: None ) –

    Default None. Phase solution that LCHS planning computes with prepare_qsp_evolution and supplies, solved for tau = encoding.alpha * evolution_time and the same epsilon. A caller leaves it at its default.

Returns:

  • encoding ( BlockEncoding ) –

    BlockEncoding of exp(-i t A) with subnormalization 1, ancillas num_ancillas + 3 (parity, pair, signal + block ancillas), the documented error_bound, and full synthesis metadata (degrees, recorded tail slack, rescale, solver residuals, amplitude deficit, raw OAA residual bound, compensated-recovery scale and bound, and query counts).

Raises:

  • ValueError –

    If evolution_time is not positive, epsilon is not in the open interval (0, 1), the metadata declares the target or encoded operator non-Hermitian, require_hermitian_evidence is True and either declaration is missing, _prepared was solved for another tau or epsilon, or the dense syntheses would exceed max_work or max_bytes. With _prepared left at None, also when prepare_qsp_evolution fails for tau = encoding.alpha * evolution_time: a degree above 256, exhausted evaluations, or no target margin that converges.

  • TypeError –

    If require_hermitian_evidence is not a bool.

Warns:

  • UserWarning –

    When either Hermitian declaration is missing and require_hermitian_evidence is False.

Inverse polynomials

Certified 1/x polynomial construction for QSVT matrix inversion.

When QLS chooses its polynomial, it fits an odd polynomial P approximating 1/(kappa x) on the domain D(a) = {x : a <= |x| <= 1} with a = 1/kappa. The fit is an odd-Chebyshev least-squares fit on Chebyshev nodes restricted to the domain, and its degree is found by doubling and then bisection until an affine Chebyshev-grid norming bound establishes max_D |kappa * x * P(x) - 1| <= epsilon_inv. The fit procedure is NWQLib's own construction.

The domain parameter kappa must exceed 1. QLS passes its polynomial_kappa, the encoding's kappa_be raised to at least QLS_POLYNOMIAL_KAPPA_FLOOR = 1.01 so that a perfectly conditioned A still gives a domain with a < 1. kappa_be is max(1, alpha / sigma_min(A)) for kappa="auto" and otherwise the supplied kappa, which must cover that value. With kappa="auto", the analytic periodic encoding passes alpha / min(sigma_min, encoded gap), rounded outward (QLS periodic encoding). Every kappa below is that domain parameter.

Two choices make the polynomial usable by QSVT and by physical recovery. The target 1/(kappa x) has magnitude at most one on D, so P needs only a small rescale into the QSP amplitude domain. The criterion is a relative residual. For each eigenvalue (Hermitian A) or singular value (general A, through its odd singular-value transform) s of A/alpha in D, |P(s) - 1/(kappa s)| <= epsilon_inv / (kappa |s|). Summing over components, the transformed vector y of b satisfies ||y - (kappa A/alpha)^{-1} b|| <= epsilon_inv ||(kappa A/alpha)^{-1} b||. Childs, Kothari, and Somma [CKS], arXiv:1511.02306v2, and Sunderhauf, Nemeth, Walayat, Patterson, and Berntson [SNWPB], arXiv:2507.15537v1, instead bound the absolute error |p(x) - 1/x|. Equation, lemma and section numbers refer to these arXiv versions.

InverseChebyshevFit 

InverseChebyshevFit(*, coefficients: tuple[float, ...], degree: int, kappa: float, epsilon_inv: float, certificate: float, certificate_grid_points: int, lsq_node_count: int)

Certified odd-Chebyshev fit of 1/(kappa * x) on D(1/kappa).

The QLS Method's planning produces it, and no public function of this module returns a fitted instance. Its coefficients define the polynomial P. The fields below are read-only.

Attributes:

  • coefficients (tuple[float, ...]) –

    Full Chebyshev coefficient vector (c_0..c_d) of the fitted P. Even entries are exactly zero.

  • degree (int) –

    Odd degree d chosen by the bounded search.

  • kappa (float) –

    Domain parameter, greater than 1, defining the domain D = {1/kappa <= |x| <= 1}. QLS passes its polynomial_kappa.

  • epsilon_inv (float) –

    Requested certificate tolerance.

  • certificate (float) –

    Chebyshev-grid norming upper bound on max_D |kappa * x * P(x) - 1|, evaluated in binary64, not interval arithmetic.

  • certificate_grid_points (int) –

    Affine Chebyshev-grid size behind the bound.

  • lsq_node_count (int) –

    Least-squares node count of the accepted fit.

Dalzell kernel reflection

Kernel-reflection polynomial of Dalzell's shortcut quantum linear solver.

Dalzell, "A shortcut to an optimal quantum linear system solver", arXiv:2406.12086v2. Equation and page numbers refer to that version. Algorithm 1 (p. 5) applies kernel reflection (KR) to |e_n> through QSVT on the augmented matrix G_t of Eq. (11). KR uses the even degree-2 ell polynomial K of App. B.3, Eq. (62), built from the filter F of App. B.2, Eq. (52), with Delta = 1/kappa and ell from Eq. (6). Lemma 3 gives K(0) = 1 and |K| <= 1 on [-1, 1], and p. 4 states -1 <= K(x) <= -1 + 4 eta/(1 + eta) for 1/kappa <= x <= 1. Lemma 3, item 2, prints this upper bound for |K(x)|, which cannot hold for eta < 1/3 because the bound is then negative. The proof, through Eq. (63), bounds K(x) itself. That signed bound is the statement of p. 4 and the one behind Eqs. (5) and (17). The QLS Method's circuit and its classical numerical model use this polynomial.

KernelReflectionPolynomial 

KernelReflectionPolynomial(*, coefficients: tuple[float, ...], kr_ell: int, kappa_be: float, kr_eta: float)

Kernel-reflection polynomial K of Dalzell arXiv:2406.12086v2, App. B.3, Eq. (62).

plan_kernel_reflection returns a coefficient-free plan, and the QLS Method's planning fills in the coefficients. The domain gap is Delta = 1/kappa_be. The fields below are read-only.

Attributes:

  • coefficients (tuple[float, ...]) –

    Chebyshev coefficients of the even kernel-reflection polynomial, or empty for a coefficient-free sizing plan.

  • kr_ell (int) –

    Integer half-degree from Dalzell arXiv:2406.12086v2, Eq. 6.

  • kappa_be (float) –

    Reciprocal of the gap Delta = 1/kappa_be, a lower bound on the nonzero singular values of the normalized input. It is not a matrix condition number. QLS passes its polynomial domain parameter, which is at least its encoded inverse-gap parameter alpha/sigma_min(A) or the caller's bound on it (the kappa setting of QLS). That parameter can exceed sigma_max(A)/sigma_min(A) (QLS guide).

  • kr_eta (float) –

    Kernel-reflection approximation parameter in (0, 1].

degree 

degree: int

Even KR degree, structurally 2 * kr_ell (Dalzell arXiv:2406.12086v2, Eq. 6).

dalzell_eta_from_precision 

dalzell_eta_from_precision(epsilon_inv: float) -> float

Return the known-norm shortcut choice eta = eps / sqrt(2).

Dalzell (arXiv:2406.12086v2, pp. 5-6, before Eq. (22)) chooses this kernel-reflection parameter for the norm estimate equal to the solution norm, t = ||x||. Eq. (20) bounds the output trace distance by eta / cos(theta_t), and cos(theta_t) = 1/sqrt(2) at that t, so the bound equals eps. The shortcut solvers use this eta for every numeric t. For t != ||x|| the same bound is eta / cos(theta_t) with theta_t = arctan(||x||/t) (Eq. (7)).

plan_kernel_reflection 

plan_kernel_reflection(kappa_be: float, eta: float) -> KernelReflectionPolynomial

Return the coefficient-free sizing plan with the half-degree of Dalzell arXiv:2406.12086v2, Eq. 6.

kr_ell = ceil(kappa_be * ln(2 / eta) / 2). This is the upper bound in Dalzell arXiv:2406.12086v2, Eq. (51), with Delta = 1/kappa_be, so it satisfies the degree condition under which Lemma 3 holds with the filter value F(Delta) <= eta.

Parameters:

  • kappa_be (float) –

    Reciprocal of the gap Delta, greater than 1.

  • eta (float) –

    Kernel-reflection approximation parameter, in (0, 1].

Returns:

  • plan ( KernelReflectionPolynomial ) –

    Empty coefficients, the half-degree kr_ell above, and the arguments as kappa_be and kr_eta. The polynomial degree is 2*kr_ell.

Raises:

  • ValueError –

    If kappa_be <= 1, or if eta is not in (0, 1].

Assumptions and limits of the evolution builders

Hermitian assumption of build_qsp_evolution_encoding:

  • It requires Hermitian target and encoded generators. These two assumptions support the QSVT polynomial calculus and the perturbation bound of Gilyén, Su, Low and Wiebe (GSLW), arXiv:1806.01838v1, Lemma 61, p. 53.
  • The check reads existing construction metadata. Real Pauli coefficients and conjugate-paired periodic shifts establish the assumption from their terms.
  • For an external encoding without this metadata, the default mode assumes that the caller has checked both operators, issues a warning and records hermitian_premise="caller_assumption". The numerical error bound is then conditional on that assumption. require_hermitian_evidence=True instead requires both metadata entries. An explicit non-Hermitian entry is rejected in either mode.
  • The check does not expand a circuit or apply an operator. A dense complex128 matrix needs 16 * 4**n bytes, which is 16 TiB at 20 system qubits and 16 PiB at 25, before ancillas or temporary arrays, so an optional numerical check of an external operator is left to the caller, in the caller's representation and within the caller's resource budget. A declaration for the target alone does not establish that an approximate encoded block is Hermitian.

Synthesis limits of the controlled passes and branches:

  • Controlling a parity pass of build_qsp_evolution_encoding, or a branch of build_control_diagonal_generator_encoding, synthesizes each dense UnitaryGate that it holds exactly, such as the unitary of a dense_dilation child. On the gate-wise route, which the parity control of a pass always takes, Qiskit then controls the synthesized gates at every occurrence.
  • Both builders take max_work (default 1,000,000,000) and max_bytes (default 10,000,000,000), and they check the work and bytes of those syntheses and of Qiskit's control of them before the first synthesis (checks of the exact synthesis).
  • build_control_diagonal_generator_encoding takes dense_control_route, whose default "auto" synthesizes a dense_dilation child together with the combine qubit's control and controls the branch's diagonal rotation gate-wise (dense control route). The parity control of a pass is gate-wise on every route.
  • A gate that an earlier control already synthesized or unrolled, such as a branch of a two-child joint generator, is controlled again by each pass outside this count, and LCHS planning counts that control.

Missing error bounds:

  • When a mathematically nonzero analytic tail is unrepresentable in binary64, the expansion's tail bounds and slack are None. Circuit construction remains available. Its BlockEncoding.error_bound and dependent OAA error claims are None, while the known normalization, recovery scale and query counts remain available.
  • A missing child error bound propagates only to claims that depend on that child.
  • NaN, infinity and negative physical error inputs are invalid, not missing bounds.

Source map

The evolution steps follow Gilyén, Su, Low and Wiebe, arXiv:1806.01838v1. The QLS source map collects the phase-convention, phase-solving, inverse-polynomial and kernel-reflection rows. The Newton step of the phase solver is derived in the docstring of phases._damped_newton, and the engineering constants record the solver sweeps of the evolution targets and the QLS targets.

Step Location in GSLW Code
Jacobi-Anger expansion of cos(tau x) and sin(tau x), parity tails Lemma 57, Eqs. (53)-(54) jacobi_anger_expansion
Bessel remainder beyond the analysis terminal Real-argument inequality in Eq. (55), derived from DLMF 10.9.4 and summed as a geometric suffix in Proposition 24 jacobi_anger_expansion
Parity combination with coefficients (1/2, -i/2) and 3-step OAA Proof of Theorem 58, and Theorem 28 with n = 3 build_qsp_evolution_encoding
Child encoding error in physical time Lemma 61 qsp_evolution_error_terms
OAA residual and amplitude-deficit bound NWQLib derivation in the docstring qsp_evolution_error_terms
Reference query scaling Corollary 60 jacobi_anger_expansion