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 QLS and Dalzell's kernel reflection. The conventions define the normalization and phase signs used here. Import the functions from nwqlib.subroutines.qsp, or from the submodule (phases, evolution, inverse, shortcut) named in each entry's path.

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; 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 (classical_preprocessing_scope: "qsp_phase_factor_optimization").

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. The answer is phases. 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, ...]) –

The solved-for coefficient vector.
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).

chebyshev_norming_sup_bound ¶

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

Bound a degree-degree polynomial's sup norm from 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--Zeller, doi:10.1007/BF01111276, Satz 1, treats 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.

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 module docstring).

Returns:

ndarray –

(reflection_phases, global_phase) such that the reflection-form
float –

product times e^{i global_phase} equals the Wx-form product;
tuple[ndarray, float] –

global_phase = d * pi / 2 compensates the i^d factor.

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, vendor workspaces 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. Without Qiskit, jacobi_anger_expansion computes its polynomial, with the degree and the Bessel search bounded and the analytic infinite-tail term kept, and prepare_qsp_evolution solves its phases within one evaluation limit. qsp_evolution_error_terms gives its error terms without a circuit. The assumptions and limits below apply to the evolution 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; [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 realized by the X-conjugation trick: flip a signal qubit onto the Pi subspace (X controlled on the ancillas being all zero), rotate RZGate(+2 phi) (= e^{-i phi Z}, the Qiskit sign RZGate(theta) = exp(-i theta Z / 2)), and flip 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 U_BE ladder: a CNOT from the pair qubit onto the signal qubit conjugates every RZ, negating 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 - 4BB^dagB). 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 LCHS effective Hamiltonian, 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 feasible degree at the first certified tail terminal. This is not a global minimum over all possible terminals.
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 integer degree anchors 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 evaluated directly.
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 registered target margin that succeeded.
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 first admissible terminal T directly and bounds the rest by an analytic remainder (NWQLib's power-series bound, for which [GSLW] Eq. (55) is the sharper real-argument form). Each parity tail carries the whole analytic suffix, so the combined tail_bound counts it twice, which is conservative. The scaling of [GSLW] arXiv:1806.01838v1, Cor. 60, d = Theta(tau + log(1/eps)/log(e + log(1/eps)/tau)), is the documented reference form for the resulting degree. min_degree floors the selected degree; the phase solver needs min_degree=2 so the even (cos) part is never a constant. max_degree bounds the returned polynomial and the finite search: at most max_degree + 200 Bessel orders are explored. max_bytes bounds known numerical arrays, not special-function workspaces or process memory.

Parameters:

tau (float) –

Effective evolution time alpha * t.
epsilon (float) –

Truncation tolerance for the combined expansion.
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 their tail bounds and tail_slack.

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:

QuantumCircuit –

Circuit on 1 + num_ancillas + system_qubits qubits with the
QuantumCircuit –

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) via the shared-ladder +-Phi pair.

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.

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 downstream 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:

BlockEncoding –

BlockEncoding of the joint generator (implementation
BlockEncoding –

"control_diagonal_generator_lcu"; an all-zero diagonal on one
BlockEncoding –

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 registered order and the first one for which both phase solves converge is kept. Only a phase-convergence failure moves to the next margin. The allowance 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:

prepared ( QSPPreparedEvolution ) –

Both phase solutions, the scale and the margins tried.

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 selected QSP evolution without building a circuit.

NWQLib's derivation. Before amplification the combined half block is B = a (V + E) with a = 1/(2s) and 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.

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 a caller premise, recorded in the output with a warning. Explicitly non-Hermitian inputs are rejected in either mode. This option reads metadata only and never materializes 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 on every route of the joint generator, 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.

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 belongs to the source representation and its resource budget, since a circuit matrix requires space exponential in its width.

Dense Hermitian evidence uses exact entry equality, so even an assembly roundoff asymmetry 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.

Returns:

BlockEncoding –

BlockEncoding of exp(-i t A) with subnormalization 1, ancillas
BlockEncoding –

num_ancillas + 3 (parity, pair, signal + block ancillas), the
BlockEncoding –

documented error_bound, and full synthesis metadata (degrees,
BlockEncoding –

recorded tail slack, rescale, solver residuals, amplitude deficit,
BlockEncoding –

raw OAA residual bound, compensated-recovery scale and bound, and
BlockEncoding –

query counts).

Inverse polynomials¶

Certified 1/x polynomial construction for QSVT matrix inversion.

QLS polynomial selection builds an odd polynomial P approximating 1/(kappa x) on the domain D(a) = {x : a <= |x| <= 1} with a = 1/kappa: odd-Chebyshev least squares on domain-restricted Chebyshev nodes, degree found by double-then-bisect against an affine Chebyshev-grid norming bound on 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).

QLS polynomial selection produces it, and 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) –

Selected odd degree d from 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 circuit and host-model users are the QLS Method's circuit and its numerical model.

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 QLS 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; 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 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 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 stored terms.
For an external encoding with missing evidence, the default 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 belongs to the caller's representation and 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 arXiv:1806.01838v1	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	NWQLib power-series bound. Eq. (55) is the sharper real-argument form.	`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`