Carleman Linearization¶

Category C2 · Module oracq.algorithms.qnlss.carleman · Stage V2

Overview¶

Lifts a polynomial-nonlinear ODE into a finite-dimensional linear system. Write

\[ u'=\sum_{p=0}^{D}F_p\,u^{\otimes p},\qquad u^{\otimes 0}=1, \]

where \(F_0\) is a constant forcing vector, \(F_1\) a linear operator, and \(F_p\) (\(p\ge 2\)) a contraction that is linear in the whole tensor input. Letting \(y_k=u^{\otimes k}\), the product rule gives

\[ y_k'=\sum_{p=0}^{D}\sum_{j=0}^{k-1}\bigl(I^{\otimes j}\otimes F_p\otimes I^{\otimes(k-j-1)}\bigr)\,y_{k+p-1}. \]

The truncation order \(K\) keeps \(y_0,\dots,y_K\) and drops terms sourced from orders greater than \(K\), yielding the linear generator \(G_K\); the finite matrix construction is deterministic, and how well it approximates the original nonlinear equation depends on the truncation and the applicability conditions (for the quantum convergence/efficiency basis see the Carleman paper arXiv:2011.03185v4; see also the differential equations chapter §6).

Interface and input model¶

PolynomialODE(width, coefficients, initial, initial_norm=1.0)
carleman_lift(problem, *, cutoff=2)
carleman_initial(problem, *, cutoff=2)
carleman_qode(problem, time, linear_solver, *, cutoff=2)

API entries: PolynomialODE, carleman_lift, carleman_initial, carleman_qode

width: the physical register width \(n\) (dimension \(d=2^n\), capped at 64).
coefficients: ((p, the BE of F_p), ...), with non-negative, non-repeating degrees; the BE width of \(F_p\) must be \(\max(1,p)\cdot n\) (the zero-signal block contract padded to a square: the first \(d\) rows are valid, the remaining rows are declared zero).
initial / initial_norm: the StatePreparation (SP) of the physical initial state and the raw vector norm; the norm determines the relative amplitudes of the orders in the lifted initial state — it is not a decorative attribute.
linear_solver: the three-argument protocol (generator, initial, time) -> StateOracle, usually injected as linear_qode(...).

The three entries return, respectively: a BlockEncoding of the lifted generator \(G_K\) (attributes algorithm="carleman_lift", cutoff, coefficient_assumption), a StatePreparation of the lifted initial state, and the final StateOracle (attributes algorithm="carleman_qode", cutoff, truncation_assumption="Carleman tail pending").

Implementation notes¶

The register layout is the padded layout: target = data region (K·n bits, low) + level register (K.bit_length() bits, high); the valid data of order \(k\) occupies only the first \(k\cdot n\) bits, the rest being zero — not the compact \(1+d+\dots+d^K\) addressing; the resource difference is kept in reports.

For each \((k,p,\text{position})\) (\(1\le k\le K\), source order \(k+p-1\le K\), position \(0..k-1\)), carleman_lift generates one placement BE: the level register routes from the source order to the output order (XOR flips the differing bits), \(F_p\) acts on the \(p\) consecutive \(n\)-bit groups starting at position, and when \(p>1\) the output stays in the lowest \(n\) bits of the group with the zero rows moved to the high end of the data region; summing all terms with an LCU gives \(G_K\).

carleman_initial prepares the normalized lifted vector

\[ \frac{(1,\,u_0,\,u_0^{\otimes 2},\,\dots,\,u_0^{\otimes K})}{\sqrt{\sum_{k=0}^{K}\lVert u_0\rVert^{2k}}}, \]

calling the initial-state oracle multiple times, controlled by level (work reserves \(K\times\) the initial state’s work width). carleman_qode hands \(G_K\) and the lifted initial state to the injected linear solver, then uses select_subspace to select the level 1 channel (the selection value 1 << ((cutoff-1)*width) lands exactly on the lowest bit of the level register, while also requiring the high data bits to be zero), with the selection condition folded into the signal.

Applicability boundary: \(G_K\) is not guaranteed to be dissipative, even if the original PDE or \(F_1\) is. Connecting Schrödingerization can pass \(G_K\) directly (judge the recovery region separately); connecting LCHS must confirm \(\mathrm{Hermitian}(-G_K)\succeq 0\), or apply a global shift \(\mu\ge\alpha_{G_K}\) acting on the whole lifted system (including level 0), and the shift cost should be visible in the application configuration. The truncation error of the tensor tail terms is not quantified (truncation_assumption is marked pending).

Validation approach¶

Category C2 (approximately continuous semantics; the acceptance criteria are in ../development/validation-plan.md §2). Three layers of evidence:

Structure: the Carleman branch of tests/core/test_differential.py:DifferentialStructureTests.test_four_methods_keep_input_oracles — with PolynomialODE(1, ((1, a), (2, f2)), initial) taking abstract BEs as coefficients and CBMD as the linear solver, the open slots input_A/input_b/input_F2 are all kept, a placement module with carleman_row_level == 1 and carleman_column_level == 2 exists, and a serialization round trip holds.
Numerical: tests/core/test_differential.py is positioned as a structure test; numerical cross-checks of truncation convergence belong to the unified ODE convergence benchmark (see known gaps).
Binding: same-type registration and contract checks — the coefficient BE widths, degree uniqueness, and initial-state layout are validated at generation time in PolynomialODE.__post_init__, with violations raising ValidationError.

Known gaps and planned stages¶

A unified “analytically solvable ODE family” convergence benchmark is missing (the same problem through all solvers, including a Carleman truncation-order sweep); it is assigned to the stage V2 convergence-sweep framework, consistent with the carleman.py row of validation-coverage.md.

Numerical validation¶

Paper-grade numerical experiments are in tests/verification/verify_ode.py (the ode group; this file covers the parts of carleman.py). All classical references are independent: the lifted generator \(G_K\) assembled in numpy basis vector by basis vector on the padded layout per the product rule, the lifted initial state constructed analytically from tensor powers, the truncated linear system via scipy.linalg.expm, and the exact nonlinear solution via scipy.integrate.solve_ivp (Riccati, rtol=1e-12); the quantum programs run on the real backends reference, rir-pysparq, adapter-pysparq, and OriginIR-ext (UniQC full amplitude / to_matrix).

Experiment design: (a) lifted block — the Riccati problem \(u'=-u+u^2\) (\(F_1=-I\), \(F_2\) the contraction matrix \(C e_{i_0+2i_1}=\delta_{i_0,i_1}e_{i_0}\), \(d=2\)) with truncation \(K=2\); the zero-signal block of \(G_K\) is extracted (unitary to_matrix + effective_block, cross-checked with a reference zero-state sweep) and compared against the independent assembly; (b) lifted initial state — \((1,u_0,u_0^{\otimes2})/Z\) (\(u_0=(0.6,0.8)\), \(r=0.5\)) exhaustively enumerated at full amplitude on four paths; (c) end to end — carleman_qode with an injected three-argument protocol solver (a minimal truncated-Taylor \(e^{Gt}\) block-encoding assembly inside the verification script, a real quantum program, degree 3), \(t=0.2\), the physical channel against a full-stack numpy simulation; (d) truncation trend — the truncation errors of the \(K=1,2\) quantum runs against the exact solve_ivp Riccati solution, with \(K=3\) filled in classically.

Key metrics:

Case	Scale	Paths	Metric	Value
carleman-lift-block	\(K\)=2, 10 qubits	originir+to_matrix, reference	maximum element-wise deviation of the \(G_K\) block	3.4e-16
carleman-initial-state	4 qubits	four paths	maximum full-amplitude deviation	0.0
carleman-riccati-endtoend	18 qubits	reference, rir	implementation error (vs independent simulation)	2.8e-17
ditto			Taylor remainder (degree 3, informational)	1.4e-4
ditto			Carleman truncation error (vs the exact nonlinear solution)	1.1e-3
ditto			quantum direction vs the exact truncated-linear solution direction	4.8e-5
carleman-cutoff-trend	\(t=0.2\)	quantum (K=1,2) + classical (K=3)	truncation error K=1 / K=2 / K=3	9.3e-3 / 1.1e-3 / 1.05e-4

The truncation error drops by about one order of magnitude per added order — direct numerical evidence of Carleman truncation convergence. The success-subspace probability (all signals zero) also matches the classical value pointwise in the cross-check.

Reproduction:

PYTHONPATH=src /home/agony/projects/qcfd-dev/quantum-cfd-software/.venv/bin/python tests/verification/verify_ode.py

Artifacts: out/verification/ode.json (5 carleman-* cases in total).