LCHS

Use LinearDynamics(A=..., initial_state=..., time=..., source=...) with LCHS to approximate the solution of du/dt = -A u + b at the final time. Choose the physical solution, a state vector, a norm, an observable or samples through the output quantities. The LCHS guide explains finite quadrature, the evolution backends, scale recovery and the error components.

from nwqlib.algorithms.lchs import LCHS, LCHSRefinement, LCHSVerification

Solve a linear ODE

LCHS

Bases: Method

LCHS method for the solution of a linear ODE given by LinearDynamics.

Build it with keyword arguments and pass it as method=, for example solve(LinearDynamics(A=A, initial_state=u0, time=T), method=LCHS()). Every argument is optional. The result is an LCHSAnalysis whose solution, for the default Solution output, approximates u(T) = exp(-A T) u0 + (integral_0^T exp(-A s) ds) b, the solution of du/dt = -A u + b at elapsed time T, with u0 the initial_state and b the source (An, Childs and Lin, arXiv:2312.03916v2, here ACL, Eq. (2)).

Linear combination of Hamiltonian simulation (LCHS) writes exp(-A T) as an integral over real k of g(k) exp(-i T (k L + H)), with A = L + iH, L = (A + A^dagger)/2 and H = (A - A^dagger)/(2i) (ACL Eqs. (3)-(4), (6) and (60)), and replaces the integral by a finite sum over quadrature nodes k_j. L must be positive semidefinite, or the default make_l_psd=True shifts it. The default kernel is ACL's near-optimal kernel, Eq. (7), with beta = 0.75, on composite Gauss quadrature. lchs_kernel also offers the Cauchy kernel (An, Liu and Lin, doi:10.1103/PhysRevLett.131.150603, restated as ACL Eqs. (5) and (181)) and the Low–Somma kernel (Low and Somma, arXiv:2508.19238v2, Eq. (6)), whose integral representation is approximate (their Eq. (8)).

approximation_tolerance bounds only the kernel-integral and k-quadrature components of the error. It does not bound the total error of the physical output, which also depends on input scaling, PSD recovery, state preparation, evolution and floating-point error. The default dense_exact backend computes each branch matrix exp(-i t (k_j L + H)) classically and builds one dense SELECT circuit, which applies branch j when its address register holds j. The product-formula and qsp_block_encoding backends build more structured circuits, each with its own limit checks.

Vector outputs (Solution, StateVector) need exact amplitude readout, so they take no shots. Samples needs execution="quantum" and positive shots. execution="classical" takes no shots and does not run the qsp_block_encoding backend. result.analyze() takes no settings, because the Plan fixes the recovery scale and coordinates. result.verify(checks=...) takes one LCHSVerification or LCHSRefinement, extra computation that plan, solve and analysis never run. The LCHS guide explains kernel and quadrature selection, the evolution backends and each error component.

Examples:

The exact solution of this two-coordinate problem at T = 0.1 is exp(-0.1 A) [1, 0], approximately [0.96082565, -0.00484022]. The default run on local Aer differs from it by about 0.000821 in the L2 norm.

>>> import numpy as np
>>> from nwqlib import LinearDynamics, solve
>>> from nwqlib.algorithms.lchs import LCHS
>>> problem = LinearDynamics(A=[[.4, .15], [.05, .25]],
...                          initial_state=[1, 0], time=.1)
>>> result = solve(problem, method=LCHS())
>>> print(np.round(result.solution, 4))
[ 0.9601-0.j -0.0051+0.j]

Attributes:

• approximation_tolerance (Real) –

  Default 0.01, strictly between 0 and 1. Construction tolerance of the kernel-quadrature pair, in the operator norm for a unit input before PSD growth. The default pair gives half to the cutoff tail and half to the k quadrature. QSP synthesis, and each node of trotter_error_budgeted, separately receives 0.1 times this value. Planning refuses a coefficient table whose binary64 rounding u*alpha exceeds 0.1 times this value, with u = 2**-53 and alpha the coefficient 1-norm. It is not a tolerance on the total physical-output error.

• lchs_kernel (Any) –

  Default "near_optimal_eq7" with beta=0.75. Kernel g(k) of the integral, given as a name, as a dict {"implementation": name, "parameters": {...}} or as a ProviderConfig. "near_optimal_eq7" is ACL Eq. (7) with beta strictly between 0 and 1. "cauchy_density" is the Cauchy kernel and takes no parameters. "low_somma_f2" is the Low–Somma kernel with shift c > 0, default 1, whose coefficient 1-norm grows with c. It requires approximation_tolerance at most 4/5 and gives the kernel-integral and quadrature components a third of it each (Low and Somma, Theorems 2-4).

• k_quadrature (Any) –

  Default "composite_gauss" with truncation_multiplier=1. Quadrature rule in k, given like lchs_kernel. "composite_gauss" chooses the cutoff and the Gauss–Legendre panels from the tolerance, and truncation_multiplier, at least 1, enlarges the cutoff. "signed_binary_uniform" is a uniform grid of 2**num_qubits nodes with spacing 2**lsb_position, both required, and carries no error bound. "symmetric_uniform_trapezoid" is the Low–Somma trapezoid rule and takes no parameters. The allowed pairs are near_optimal_eq7 with composite_gauss (both component bounds hold) or with signed_binary_uniform, cauchy_density with signed_binary_uniform, and low_somma_f2 with symmetric_uniform_trapezoid (both component bounds hold, and the synthesis stages are not budgeted with them). Building the Method refuses any other pair.

• hamiltonian_evolution_backend (Literal['dense_exact', 'trotter', 'trotter_error_budgeted', 'qsp_block_encoding']) –

  Default "dense_exact". How each branch exp(-i t (k_j L + H)) is built. "dense_exact" computes the branch matrix from the node's Hermitian eigensystem and controls its circuit. "trotter" uses a sorted-Pauli product formula of order trotter_order with trotter_steps steps. "trotter_error_budgeted" gives each node the smallest step count whose product-formula bound plus the Pauli pruning bound fits 0.1 times approximation_tolerance (Childs et al., doi:10.1103/PhysRevX.11.011020, Propositions 9-10 and Sec. V B). "qsp_block_encoding" evolves under the joint generator with quantum signal processing on block encodings of L and H (Pocrnic et al., arXiv:2506.20760v2, Section IV) and needs quantum execution. A PeriodicStencil A needs "trotter" with trotter_order=2 and quantum execution.

• lcu_select_implementation (Literal['auto', 'multiplexor', 'structured', 'branch_controlled']) –

  Default "auto". How the SELECT circuit applies the branches. "dense_exact" accepts "auto" or "branch_controlled", and "qsp_block_encoding" accepts only "auto". With "trotter_error_budgeted", "auto" gives "multiplexor", and "structured" is refused, because per-node step counts break the affine address structure. With "trotter", "auto" gives "structured" only when an affine address structure is eligible, as checked from the k-node addresses and the real Pauli coefficients of the generator, and its CX count is lower than the multiplexor's, and gives "multiplexor" otherwise. An explicit "structured" without that eligibility is refused. A PeriodicStencil A accepts "auto" or "structured".

• dense_control_route (Literal['gatewise', 'whole_matrix', 'auto']) –

  Default "auto". How a construction controls a dense unitary on two or more qubits. It applies to a dense_exact branch on its address bits and to a dense-dilation QSP child on the combine qubit of the joint generator. "gatewise" synthesizes the unitary and lets Qiskit control each synthesized gate, "whole_matrix" synthesizes the controlled matrix, and "auto" takes the whole-matrix route for one control and the gate-wise route for more. On the whole-matrix route a constant-source branch is one unitary, its evolution times the matrix of its input preparation. The QSP pass that holds the generator is controlled gate-wise on every route.

• trotter_steps (PositiveInt | None) –

  Default 1. Fixed step count of the "trotter" backend, a positive integer at most max_trotter_steps, or None when trotter_synthesis_tolerance chooses it. Exactly one of the two is set. "trotter_error_budgeted" chooses its own counts and requires the default 1. "dense_exact" and "qsp_block_encoding" apply no product formula, so they treat any value as 1, and planning warns about a value other than 1.

• trotter_synthesis_tolerance (Real | None) –

  Default None. Positive tolerance on the Strang synthesis bound of a PeriodicStencil A, in the same unit-input operator norm as approximation_tolerance. Planning chooses the smallest step count r with B/r**2 at most this value, where B = T**3 sum_j |c_j| (2 |k_j| diffusion + |potential|)**3 / 3 over the nodes k_j and coefficients c_j (Childs et al., doi:10.1103/PhysRevX.11.011020, Proposition 10, Eq. (121), relaxed by NWQLib). It requires hamiltonian_evolution_backend="trotter" and trotter_steps=None, and other inputs refuse it.

• trotter_order (Literal[1, 2]) –

  Default 2. Order of the product formula of the "trotter" and "trotter_error_budgeted" backends, 1 for the Lie formula or 2 for the symmetric Suzuki (Strang) formula. "dense_exact" and "qsp_block_encoding" accept only 2.

• max_trotter_steps (PositiveInt) –

  Default 100_000. Upper limit on the product-formula step count, fixed or chosen automatically.

• duhamel_nodes (PositiveInt) –

  Default 8. Number of Gauss–Legendre nodes of the time integral of a constant source, one panel on [0, T] (ACL Eq. (72)). It is never reduced automatically to fit a limit.

• make_l_psd (StrictBool) –

  Default True. Rule for an eigenvalue of L below -psd_tolerance*||L||_2. True shifts L by s = -lambda_min + psd_tolerance*||L||_2 and multiplies the result by the growth factor exp(s*T). False refuses the input.

• psd_tolerance (Real) –

  Default 1e-12, positive. Relative window of the PSD check. An eigenvalue of L at or above -psd_tolerance*||L||_2 is accepted without a shift. The window is the scale of the eigensolver's rounding (LAPACK Users' Guide, 3rd ed., Sec. 4.7), so the decision does not depend on the time unit. It does not permit a negative eigenvalue below that window.

• initial_state_preparation (Literal['direct', 'mps_circuit']) –

  Default "direct", which prepares the initial state exactly. "mps_circuit" prepares it approximately with a layered circuit from a matrix-product-state (MPS) compression. Quantum execution with a constant source, or with a PeriodicStencil A, requires "direct".

• initial_state_mps_max_bond_dim (PositiveInt | None) –

  Default None, no limit. Upper limit on the kept bond dimension of the initial-state compression.

• initial_state_mps_threshold (Real) –

  Default 1e-14, positive. Singular-value truncation threshold of the initial-state compression.

• initial_state_mps_num_layers (PositiveInt) –

  Default 2, positive. Number of layers of the initial-state MPS circuit.

• lcu_state_preparation (Literal['direct', 'mps_circuit']) –

  Default "direct". Preparation of the coefficient state sqrt(|c_j|/alpha) that weights the branches, with the same choices and restrictions as initial_state_preparation. The circuit error of "mps_circuit" stays unevaluated until an explicit validation, even when the compression discards no weight.

• lcu_mps_max_bond_dim (PositiveInt | None) –

  Default None, no limit. Upper limit on the kept bond dimension of the coefficient-state compression.

• lcu_mps_threshold (Real) –

  Default 1e-14, positive. Singular-value truncation threshold of the coefficient-state compression.

• mps_num_layers (PositiveInt) –

  Default 2, positive. Number of layers of the coefficient-state MPS circuit.

• max_dense_select_slots (PositiveInt) –

  Default 256. Upper limit on the padded address count, a power of two, of the dense SELECT circuit of "dense_exact". Planning refuses a larger construction and never changes the grid, Duhamel nodes or tolerance to fit.

• max_bytes (PositiveInt) –

  Default 10 GB (decimal, 10_000_000_000 bytes). Upper limit on the known bytes of input and numerical workspace.

• max_svd_work (PositiveInt) –

  Default 1e8 (100_000_000). Upper limit on the counted work of state and coefficient compression, including the construction of an "mps_circuit" layered circuit.

• max_spectral_work (PositiveInt) –

  Default 1e8. Upper limit on the counted work of the eigenvalue check of L for dense A, d**3 for physical dimension d, so the default accepts a dense A of dimension at most 464 when L is nonzero. A periodic stencil has an analytically PSD L and no eigenvalue check.

• max_quadrature_work (PositiveInt) –

  Default 1e8. Upper limit on the counted work of building the k nodes and weights. For "composite_gauss" it covers the completed panel search plus the construction of the chosen rule, not every scalar probe.

• max_select_work (PositiveInt) –

  Default 1e8. Upper limit on the counted work of building the SELECT circuit, including the branch matrices and their synthesis for "dense_exact" and the Pauli decomposition of the product-formula backends. The LCHS guide gives the sizes that fit this default.

• max_readout_work (PositiveInt) –

  Default 2_000_000_000. Upper limit, inclusive, on the comparisons of grouping sampled Pauli terms and the input visits of exact scalar readout (shots=None). Planning checks and records the grouping comparisons. Before the circuits run, each exact readout is checked against what remains, and completed exact readouts count against it.

• max_qsp_degree (PositiveInt) –

  Default 256. Upper limit on the QSP polynomial degree of "qsp_block_encoding".

• max_qsp_evaluations (PositiveInt) –

  Default 20_000. Upper limit on the residual evaluations of the QSP phase solver.

• max_admission_steps (PositiveInt) –

  Default 1_000_000, ten times the shared default. Upper limit on the planning work of checking each Program of the Plan (NWQLib's description of a circuit as named steps), counting its stored fields and the structural and lifecycle work of one check. Summing a Program's resource counts may use up to 24 times this value. Raising it changes no quantum operation. See the planning work limit.

Raises:

• ValueError –

  If the kernel-quadrature pair is not allowed or a kernel or quadrature parameter is invalid, if both or neither of trotter_steps and trotter_synthesis_tolerance are set, if trotter_synthesis_tolerance is set without the "trotter" backend, if "trotter_error_budgeted" gets trotter_steps other than 1, if trotter_steps exceeds max_trotter_steps, or if trotter_order=1 is set without a product-formula backend.

• TypeError –

  If lchs_kernel or k_quadrature is not a name, a dict or a ProviderConfig.

ProviderConfig

ProviderConfig(*, implementation: str, parameters: Mapping[str, ProviderParameter] = dict())

Choice of an LCHS kernel or k-quadrature rule, with its explicit parameters.

Build it with keyword arguments and pass it as lchs_kernel= or k_quadrature= of LCHS or of resolve_lchs_coefficient_plan, for example LCHS(lchs_kernel=ProviderConfig(implementation="near_optimal_eq7", parameters={"beta": 0.7})). implementation is the only required argument. LCHS also accepts the name alone or a dict with these two keys. The lchs_kernel and k_quadrature rows of LCHS list the names and their parameters.

Attributes:

• implementation (str) –

  Required. Kernel or quadrature name. Building this record does not check the name, and LCHS refuses an unknown one.

• parameters (Mapping[str, ProviderParameter]) –

  Default empty. Parameter names mapped to JSON scalars or None, stored sorted and read-only. The named rule supplies the defaults and range checks when LCHS resolves it.

ProviderParameter

ProviderParameter: TypeAlias = str | int | float | bool | None

Type of one kernel or quadrature parameter value, a JSON scalar or None.

Read the result

LCHSAnalysis

Bases: Result

Solution or requested output of du/dt = -A u + b from an LCHS run.

solve returns it for an LCHS method, and load_result reopens a saved one. For the default Solution output the answer is solution, the physical solution u(T) in the original coordinates, with its norm and phase, in the problem's unit. A StateVector output is in state_vector. NormSquared, QuadraticForm and NormalizedExpectation give the scalar value, and Samples gives samples. Counts do not recover a complex vector. When the output cannot be recovered, unavailable gives the reason, for example a normalized output at zero physical mass. print(result) shows the value with the construction tolerance and the kernel bounds, each labeled as not a total error bound. result.analyze() takes no settings. The fields below are read-only. The fields of Result are present too.

The kernel-integral and k-quadrature bounds are facts of the Plan (result.plan.facts), and a classical Result also carries them in result.facts. They are physical L2 bounds of error components, not a bound on the total error of the output. No reference solution is computed unless result.verify is called.

Attributes:

• value (Real | None) –

  Requested scalar of a NormSquared, QuadraticForm or NormalizedExpectation output, physical or unit-normalized as that output defines it. None for vector and Samples outputs and when unavailable.

• unavailable (Text | None) –

  Reason the requested output cannot be recovered, or None.

• norm_squared (Nonnegative | None) –

  Recovered physical norm squared ||u(T)||**2, when available.

• physical_scale (PhysicalScale | None) –

  Physical scale of the obtained vector, when representable.

• numerator (Real | None) –

  Observable value, the normalized ratio for a NormalizedExpectation and the recovered physical quadratic form for a QuadraticForm.

• numerator_frame (Literal['physical', 'unit'] | None) –

  Normalization of numerator, "unit" for the normalized ratio or "physical" for the quadratic form.

• artifact (ArtifactManifest | None) –

  Description of the stored solution or state array. Printing the result does not load the array, and solution and state_vector read it.

• applications (tuple[KernelApplication, ...]) –

  Numerical operator applications of the run, each with its error-component facts.

• references (Literal['not_run']) –

  Always "not_run". An independent reference needs result.verify(checks=LCHSVerification(...)).

• submitted_shots (Count | None) –

  Shots requested for a sampled output, if applicable.

• returned_shots (Count | None) –

  Shots returned, before selection of the physical coordinates.

• selected_shots (Count | None) –

  Returned shots kept for the physical output. For Samples, the sum of the sample counts.

• samples (LCHSSamples | None) –

  For a Samples output, the distinct original-coordinate indices (samples.indices, increasing, padding coordinates excluded) and their positive counts (samples.counts), both int64 arrays. They do not cover all submitted shots.

• reduction (LCHSProjectedMoments | None) –

  Saved statistics of the exact scalar readout (shots=None): the complete, success and physical-slice masses and the projected moment, as mantissa-exponent pairs.

• groups (tuple[LCHSGroupMoments, ...]) –

  Weighted moments of each sampled measurement setting, with its Pauli labels, basis and returned and selected shots.

solution

solution

The physical solution u(T) of a Solution output, a complex array.

It is in the original coordinates and keeps the physical norm and phase. With no evolution (time equal to initial_time, or a zero input) it is the supplied initial vector.

Raises:

• AttributeError –

  If the requested output was not Solution.

• ValueError –

  If the array is unavailable, with the reason.

state_vector

state_vector

The state vector of a StateVector output, a complex array.

It is in the original coordinates, normalized and phased as the StateVector output requests.

Raises:

• AttributeError –

  If the requested output was not StateVector.

• ValueError –

  If the array is unavailable, with the reason.

Check a result

LCHSVerification compares the physical solution with an independent reference, and LCHSRefinement evaluates error components that planning leaves unknown. Each runs only when passed to result.verify(checks=...).

LCHSVerification

Bases: Record

Reference comparison for the physical solution of an LCHS Result.

Build it with keyword arguments and pass it to result.verify(checks=...), for example result.verify(checks=LCHSVerification(reference="expm")). reference is the only required argument. The call returns (receipt, facts). facts[0].fact.value.value is the discrepancy between the Result's physical solution u and the reference u_ref, and the receipt records every computed value and call count. The check is named name, and "ivp_closed_form" adds name + ".reference_consistency". Certificate.with_verification turns facts into PASS when the value is at most threshold, FAIL when it is larger, and INCONCLUSIVE when the value or the threshold is missing. Planning, solving and ordinary analysis never compute a reference.

The Result must hold a physical solution vector with its phase, from a Solution output or a physical StateVector, and A must be a dense matrix, never converted from a compact input. The work limits cover the whole call, including both references of "ivp_closed_form" and every node of "selected_grid". Dense work is a size-based count, not a floating-point operation count. Reference arrays are temporary, and only discrepancies and call counts are kept. A discrepancy is numerical evidence for this input, not a proven error bound. The LCHS guide shows a complete check.

Attributes:

• reference (Literal['expm', 'ivp', 'closed_form', 'ivp_closed_form', 'selected_grid']) –

  Required. "expm" compares with exp(-A T) u0 for dynamics without a source (An, Childs and Lin, arXiv:2312.03916v2, here ACL, Eq. (2) with b = 0). "ivp" solves du/dt = -A u + b with SciPy's RK45 under explicit rtol, atol and max_rhs_calls. "closed_form" evaluates ACL Eq. (2) for a constant source with one exponential of [[-A*T, b*T], [0, 0]]. "ivp_closed_form" runs both and also reports their signed consistency. These three need a source. A matrix with 1-norm above 2**37, or with an exponential that is not finite in binary64, gives an unknown discrepancy. "selected_grid" reproduces the saved finite sum, with its nodes and product-formula steps, so its discrepancy excludes the kernel-integral, k- and Duhamel-quadrature and product-formula errors.

• name (Text) –

  Default "reference_error". Name of the reported discrepancy and prefix of its companion values.

• metric (Literal['absolute_l2', 'relative_l2']) –

  Default "absolute_l2", the norm ||u - u_ref||_2 without a global-phase fit. "relative_l2" divides it by ||u_ref||_2.

• threshold (Nonnegative | None) –

  Default None. Optional nonnegative pass threshold of the discrepancy check.

• consistency_threshold (Nonnegative | None) –

  Default None. Optional nonnegative threshold of the consistency check between the IVP and closed-form references, only for "ivp_closed_form".

• rtol (Real | None) –

  Default None. Positive RK45 relative tolerance, at least SciPy's floor of 100 times the binary64 machine epsilon. Required for "ivp" and "ivp_closed_form" and refused otherwise.

• atol (Real | None) –

  Default None. Positive RK45 absolute tolerance. Required for "ivp" and "ivp_closed_form" and refused otherwise.

• max_rhs_calls (PositiveInt | None) –

  Default None. Upper limit on RK45 right-hand-side evaluations. Required for "ivp" and "ivp_closed_form" and refused otherwise.

• max_bytes (PositiveInt) –

  Default 10 GB (decimal, 10_000_000_000 bytes). Upper limit on the known arrays of this call.

• max_dense_work (PositiveInt) –

  Default 1e8 (100_000_000). Upper limit on the counted dense work of this call. For "selected_grid" it also covers reading the stored product-formula node table and the node actions.

• max_node_evaluations (PositiveInt) –

  Default 4096. Upper limit on the number of k nodes times operator applications evaluated.

• max_pf_operations (PositiveInt) –

  Default 100_000. Upper limit on the product-formula rotations or elementary gates that "selected_grid" evaluates.

Raises:

• ValueError –

  If an IVP reference lacks rtol, atol or max_rhs_calls, if rtol is below SciPy's floor, if these are given for another reference, or if consistency_threshold is given without "ivp_closed_form". result.verify also raises for "expm" with a source, for the source references without one, and for a non-dense A.

LCHSRefinement

Bases: Record

Error-component bounds of an LCHS Result at its saved grid and step counts.

Build it with keyword arguments and pass it to result.verify(checks=...), for example result.verify(checks=LCHSRefinement(components=("duhamel",))). components is the only required argument. The call returns (receipt, facts). facts holds the refined physical L2 components and their sum carried into the requested output, algorithmic_approximation, which result.assess can use. A request for "spectral_norms" alone leaves facts empty. The receipt holds every computed value. Refinement defines no pass or fail check.

Planning leaves two components unknown for a dense A, the Duhamel remainder of a constant source and the synthesis error of fixed-step "trotter". Refinement evaluates them at the saved quadrature, Duhamel nodes and step counts and replaces the planned values of the same name. It never computes a solution vector, chooses a new grid or step count, or converts a compact input to a dense matrix. algorithmic_approximation stays None while any component is unknown, because a missing term is not zero. It covers the algorithmic components only. Floating-point error of the backend run and sampling error remain separate terms.

Attributes:

• components (tuple[Literal['spectral_norms', 'duhamel', 'fixed_pf'], ...]) –

  Required. Distinct analyses, at least one. "spectral_norms" evaluates ||A||, ||A + shift*I|| and ||H|| and needs a dense A. "duhamel" evaluates the Gauss–Legendre remainder of the source time integral (DLMF Eqs. 3.5.19 and 3.5.21). "fixed_pf" evaluates the product-formula commutator bound at the saved step counts (Childs et al., doi:10.1103/PhysRevX.11.011020, Propositions 9-10) and needs the "trotter" or "trotter_error_budgeted" backend.

• name (Text) –

  Default "refinement". Prefix of the reported values that are not error components, such as the spectral norms.

• dense_validation (StrictBool) –

  Default False. Also evaluate dense commutator norms for "fixed_pf", as a check of the Pauli-triangle bound that does not replace it. Requires "fixed_pf", and a periodic stencil refuses it.

• max_bytes (PositiveInt) –

  Default 10 GB (decimal, 10_000_000_000 bytes). Upper limit on the known arrays of this call.

• max_dense_work (PositiveInt) –

  Default 1e8 (100_000_000). Upper limit on the counted dense work of this call.

• max_structural_work (PositiveInt) –

  Default 1e8. Upper limit on the counted work of reading the stored node table, preparing the count of Pauli commutator terms and evaluating them.

• max_node_evaluations (PositiveInt) –

  Default 4096. Upper limit on the number of k nodes times operator applications evaluated.

• max_steps (PositiveInt) –

  Default 100_000. Largest saved product-formula step count this call accepts.

Raises:

• ValueError –

  If components is empty or repeats a name, or if dense_validation is set without "fixed_pf".

Inspect the coefficient table

resolve_lchs_coefficient_plan

resolve_lchs_coefficient_plan(problem_or_plan, *, lchs_kernel=None, k_quadrature=None, approximation_tolerance=None, max_bytes=DEFAULT_INPUT_BYTES)

Return the LCHS coefficient table of a problem without a source, or of an LCHS Plan.

For a LinearDynamics without a source, it computes the k nodes, the weights and the coefficients c_j of the finite sum sum_j c_j exp(-i T (k_j L + H)) u0, with L after any PSD shift, from the given kernel, quadrature and tolerance. It builds no SELECT circuit or QSP phases, computes no output vector and runs nothing on a backend. For a Plan made by plan(problem, method=LCHS(...)), it returns the table that the Plan executes, so an analysis of the coefficient tensor always describes those coefficients. Call decompose_mps on the result to study the MPS compression of the coefficient state, as in the coefficient MPS analysis of the LCHS guide.

Parameters:

• problem_or_plan (LinearDynamics | Plan) –

  A LinearDynamics without a source and with positive elapsed time, or an LCHS Plan with a coefficient table for an input without a source.

• lchs_kernel (ProviderConfig | None, default: None ) –

  Kernel choice. None selects "near_optimal_eq7", as in LCHS. Must be None for a Plan.

• k_quadrature (ProviderConfig | None, default: None ) –

  Quadrature choice. None selects "composite_gauss", as in LCHS. Must be None for a Plan.

• approximation_tolerance (float | None, default: None ) –

  Construction tolerance, strictly between 0 and 1. None selects 0.01, as in LCHS. Must be None for a Plan.

• max_bytes (int, default: DEFAULT_INPUT_BYTES ) –

  Default 10 GB (decimal, 10_000_000_000 bytes). Upper limit on the arrays of the Cartesian parts and the PSD check, tested before they are formed. A PeriodicStencil A forms no matrix and uses the bound ||L|| <= mass + 4*diffusion.

Returns:

• coefficients ( LCHSCoefficientPlan ) –

  The nodes, weights, coefficients and component bounds. For a dense A the coefficients include the PSD growth factor exp(shift*T), divided by 2**e when the factor exceeds binary64, as recorded in derived_parameters["homogeneous_coefficient_compensation"].

Raises:

• ValueError –

  If the problem has a source or zero elapsed time, if approximation_tolerance is outside (0, 1), or if a Plan is given with any of the three choices or holds no coefficient table for an input without a source.

• TypeError –

  If problem_or_plan is neither a LinearDynamics nor a Plan.

LCHSCoefficientPlan

LCHSCoefficientPlan(*, requested_lchs_kernel: ProviderConfig, resolved_lchs_kernel: ResolvedProviderConfig, requested_k_quadrature: ProviderConfig, resolved_k_quadrature: ResolvedProviderConfig, lchs_kernel_formula_id: str, k_quadrature_formula_id: str, lchs_kernel_fixed_parameters: Mapping[str, ProviderParameter], kernel_parameter_selection: Mapping[str, str], derived_parameters: Mapping[str, ProviderParameter], nodes: tuple[float, ...], base_weights: tuple[float, ...], coefficients: tuple[complex, ...], coefficient_l1_norm: float, quadrature: LCHSQuadraturePlan, compatibility: LCHSPairCompatibility, mps_decomposition: object | None = None)

Nodes, weights and coefficients of the finite LCHS sum, shared by every backend.

resolve_lchs_coefficient_plan returns it. The finite sum is sum_j c_j exp(-i T (k_j L + H)) u0 over the nodes k_j and coefficients c_j. The fields below are read-only.

Attributes:

• requested_lchs_kernel (ProviderConfig) –

  Kernel choice as given, showing which parameters the caller supplied.

• resolved_lchs_kernel (ResolvedProviderConfig) –

  Kernel used, with its defaults filled in.

• requested_k_quadrature (ProviderConfig) –

  Quadrature choice as given.

• resolved_k_quadrature (ResolvedProviderConfig) –

  Quadrature used, or the algebraic reduction for an exactly zero L.

• lchs_kernel_formula_id (str) –

  Versioned identifier of the kernel formula.

• k_quadrature_formula_id (str) –

  Versioned identifier of the node and weight formula.

• lchs_kernel_fixed_parameters (Mapping[str, ProviderParameter]) –

  Constants of the kernel formula that are not user parameters, such as Low–Somma j=2 and y=1.

• kernel_parameter_selection (Mapping[str, str]) –

  Where each parameter value came from, normally "user" or "provider_default", with a profile identifier where one applies.

• derived_parameters (Mapping[str, ProviderParameter]) –

  Quantities derived for the pair, such as the Low–Somma gamma, R, h and J. A problem without a source also records its PSD-shift compensation here.

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

  Ordered real k values, as in quadrature.

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

  Integration weights before multiplication by the kernel.

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

  Complex coefficients c_j in node order. For a problem without a source they include the PSD growth factor exp(shift*T), divided by 2**e when the factor exceeds binary64, and the physical recovery carries 2**e. Their phases are applied in the SELECT circuit.

• coefficient_l1_norm (float) –

  sum_j |c_j|, the normalization of the coefficient state. It is neither the physical output norm nor an error estimate.

• quadrature (LCHSQuadraturePlan) –

  Node and address data of the rule, with its component bounds approximate_lchs_error_bound (kernel integral, the cutoff tail for the default kernel) and quadrature_error_bound, each for a unit input before PSD growth. None means unavailable, not zero.

• compatibility (LCHSPairCompatibility) –

  Pair status, certificate_status "certified", "uncertified", "incomplete" or "exact_algebraic" for the zero-L reduction, with its reason and unbudgeted_error_stages. None of these alone certifies a complete physical solution.

• mps_decomposition (object | None) –

  Attached TT-SVD (MPS) decomposition of the normalized, padded coefficient-state amplitudes, or None. decompose_mps reuses it only for the same tensor ordering and truncation settings.

prep_amplitudes

prep_amplitudes(*, max_bytes=DEFAULT_INPUT_BYTES)

Return the coefficient-state amplitudes sqrt(|c_j|/alpha), zero-padded.

alpha is coefficient_l1_norm, and the array length is the next power of two at or above len(coefficients). The complex phases of c_j are not in this state. The SELECT circuit applies them (An, Childs and Lin, arXiv:2312.03916v2, Lemma 24, Eq. (178), with the phases moved from the two preparation oracles into SELECT).

Parameters:

• max_bytes (int, default: DEFAULT_INPUT_BYTES ) –

  Default 10 GB (decimal, 10_000_000_000 bytes). Upper limit on the arrays, 96 bytes per padded address, checked first.

Returns:

• amplitudes ( ndarray ) –

  Read-only amplitudes, one per padded address.

decompose_mps

decompose_mps(max_bond_dim=None, threshold=1e-14, *, reuse=None, max_bytes=DEFAULT_INPUT_BYTES, max_svd_work=100000000)

Return the TT-SVD (MPS) decomposition of the prep_amplitudes tensor.

A stored or supplied decomposition is reused only when it was made from this exact normalized tensor with the same max_bond_dim and threshold. Otherwise a new TT-SVD runs under max_bytes and max_svd_work. No LCHS solve or circuit is run.

Parameters:

• max_bond_dim (int | None, default: None ) –

  Upper limit on the kept bond dimension. None means no limit.

• threshold (float, default: 1e-14 ) –

  Singular-value truncation threshold.

• reuse (MPSDecomposition | None, default: None ) –

  Decomposition to reuse instead of the attached mps_decomposition.

• max_bytes (int, default: DEFAULT_INPUT_BYTES ) –

  Default 10 GB (decimal, 10_000_000_000 bytes). Upper limit on the known arrays.

• max_svd_work (int, default: 100000000 ) –

  Upper limit on the counted TT-SVD work.

Returns:

• decomposition ( MPSDecomposition ) –

  The TT cores, bond dimensions and discarded weight, the sum of the squared singular values discarded.

Raises:

• TypeError –

  If reuse is not an MPSDecomposition.

• ValueError –

  If reuse does not match the tensor or the settings, instead of being replaced silently.

record

record() -> dict[str, Any]

Return the kernel, quadrature and pair choices with their bounds as a plain dict.

The keys name the requested and resolved kernel and quadrature, their parameters and formula identifiers, the node order, panel width and address structure, the pair's certificate status, the component bounds approximate_lchs_error_bound and quadrature_error_bound, and the unbudgeted_error_stages.

ResolvedProviderConfig

ResolvedProviderConfig(*, implementation: str, parameters: Mapping[str, ProviderParameter])

Kernel or quadrature choice after its defaults and checks are applied.

It appears in the resolved_lchs_kernel and resolved_k_quadrature fields of an LCHSCoefficientPlan. The fields below are read-only.

Attributes:

• implementation (str) –

  Name of the rule used. For an exactly zero L the algebraic reduction "unitary_reduction" replaces the requested quadrature name.

• parameters (Mapping[str, ProviderParameter]) –

  Every parameter after defaults, stored sorted and read-only. Together with implementation, they identify the configuration.

Lower-level functions

cartesian_decomposition

cartesian_decomposition(matrix: Any, *, eigcheck: bool = True, psd_tolerance: float = DEFAULT_PSD_TOLERANCE, max_bytes=DEFAULT_MAX_BYTES, max_spectral_work=100000000) -> tuple[ndarray, ndarray]

Return the Hermitian parts (L, H) of A = L + iH used by LCHS.

L = (A + A^dagger)/2 and H = (A - A^dagger)/(2i), as in An, Childs and Lin, arXiv:2312.03916v2, Eqs. (3)-(4). LCHS needs L positive semidefinite (PSD), which keeps ||exp(-A t)|| <= 1 (same paper, Lemma 21, Eq. (162)). With eigcheck=True the function also checks that condition. eigcheck=False only decomposes. It does not establish the PSD condition or replace the check of LCHS, which tests its own eigenvalue endpoints.

Parameters:

• matrix (array_like) –

  Finite square matrix A.

• eigcheck (bool, default: True ) –

  Check that L is PSD within psd_tolerance.

• psd_tolerance (float, default: DEFAULT_PSD_TOLERANCE ) –

  Default 1e-12. Relative window of the check. An eigenvalue of L at or above -psd_tolerance*||L||_2 passes.

• max_bytes (int, default: DEFAULT_MAX_BYTES ) –

  Default 10 GB (decimal, 10_000_000_000 bytes). Upper limit on the arrays of the check, 64 d**2 + 16 d bytes for dimension d, tested before they are formed. The count covers the known NumPy arrays, not the eigensolver's own LAPACK workspace.

• max_spectral_work (int, default: 100000000 ) –

  Upper limit on the eigenvalue check of a nonzero L, counted as d**3 work units, which are not timings.

Returns:

• parts ( tuple[ndarray, ndarray] ) –

  L and H as complex arrays.

Raises:

• ValueError –

  If A is not a finite square matrix, or, with eigcheck, if L has an eigenvalue below -psd_tolerance*||L||_2 or the check exceeds max_spectral_work or max_bytes.

duhamel_quadrature

duhamel_quadrature(final_time: float, node_count: int) -> tuple[ndarray, ndarray]

Return Gauss–Legendre nodes and weights on [0, T] for the source time integral.

One Gauss–Legendre panel on [0, T] discretizes the Duhamel integral of An, Childs and Lin, arXiv:2312.03916v2, Eq. (2), the single-panel case of their Eq. (72). LCHS uses duhamel_nodes of these nodes for a constant source. The rule is exact for polynomials of degree 2*node_count - 1. Its remainder bound follows DLMF Eqs. 3.5.19 and 3.5.21, and LCHSRefinement(components=("duhamel",)) evaluates it for a Result.

Parameters:

• final_time (float) –

  Positive, finite elapsed time T.

• node_count (int) –

  Number of nodes, at least 1.

Returns:

• rule ( tuple[ndarray, ndarray] ) –

  The nodes T*(x_i + 1)/2 and weights T*w_i/2, from the Gauss–Legendre rule (x_i, w_i) on [-1, 1].

Raises:

• ValueError –

  If final_time is not positive and finite or node_count is below 1.

Limits

• A and the source are constant in time. General time-dependent A or source is not implemented.
• The Hermitian part L of A must be positive semidefinite, or the default make_l_psd=True shifts it and restores the growth factor.
• approximation_tolerance and the kernel-integral and k-quadrature bounds are component bounds. They do not bound the total error of the physical output, and floating-point, backend and model errors remain separate.
• The default dense_exact backend computes its branch matrices classically and accepts at most max_dense_select_slots=256 padded addresses.
• Reference checks need a dense A and a physical solution vector.

Limitations and open work lists the open items. The LCHS source map gives the paper, equation and implementing function of each step. To add a method of your own, see Extending NWQLib.