chc.residual

Stability tier: stable.

Learned residual backends (Strategy): MLP, RBF Kolmogorov-Arnold, graph, control-affine, and two structured backbones that carry a guarantee -- port-Hamiltonian (passive, Lyapunov-stable) and Lipschitz (certified bounded gain). GP is future.

PortHamiltonianEnergy = Literal['mlp', 'icnn'] module-attribute

Which energy the port-Hamiltonian residual carries: a plain MLP, or an input-convex one.

ZeroResidual

Bases: Module

A residual that contributes nothing — recovers the pure known dynamics.

MLPResidual

Bases: Module

Learned residual r_θ(x, u) backed by an MLP over the concatenated [x, u].

One interchangeable backend; the point of the abstraction is that KAN / RBF / GP slot in here without touching dynamics or control.

ControlAffineResidual

Bases: Module

r_θ(x, u) = a_θ(x) + B_θ(x) u, both coefficients linear in control_affine_features.

This is the plant class the rest of the library already speaks about: chc.plan.certify_safety reads the control channel off the Jacobian at u = 0 and chc.spine requires a control-affine plant, so identifying a residual in this class puts the identification layer and the safety layer on the same object rather than on two different linearisations of it.

It is also the class in which the control channel is a parameter rather than a derivative of a black box, which is what makes it estimable by a partialling-out moment: see chc.dynamics_id.fit_causal_residual. A general r_θ(x, u) fitted by prediction error has no such guarantee -- under a confounded logging policy it learns the observational response.

The channel's basis has its own degree, channel_degree, which defaults to degree. With one basis for both, degree = 1 makes the channel affine in the state, B_θ(x) = C_0 + sum_l C_l x_l, which contains the constant channel the orthogonality results §18/§19 cover without being restricted to it. channel_degree = 0 is that constant channel beside a drift of any degree: n m channel coefficients where the affine one has n m (n + 1). Higher degrees make either polynomial in the state.

control_channel(x)

B_θ(x) -- the (out_dim, control_dim) response of the state rate to the action.

drift_jacobian(x)

∂a_θ/∂x at x -- the (out_dim, n) local linearisation of the drift.

The companion to control_channel, and it exists because an MPC plans on both halves of a_θ(x) + B_θ(x) u while only the second is identified causally by chc.dynamics_id.fit_causal_residual.

This is the drift at u = 0 and nothing else: at channel_degree >= 1 the channel depends on the state, so the vector field the MPC integrates has Jacobian ∂a_θ/∂x + ∂(B_θ(x) u)/∂x and this method returns only the first term. Reading stability off it is then a statement about where the actuator's coordinates put their zero rather than about the plant: under u = alpha v + beta the fitted class is closed and the drift absorbs beta times the state-dependent part of the channel, so a plant that decays everywhere its actuator can reach can still show a positive drift eigenvalue. On BOPTEST, a setpoint-actuated zone reporting its action in [15, 25] °C came back at +6.42 here while closed_loop_jacobian at the setpoint it actually held read -1.40.

For a stability check use closed_loop_jacobian over the admissible action set.

closed_loop_jacobian(x, u)

∂(a_θ + B_θ u)/∂x at (x, u) -- the linearisation the horizon actually follows.

Equals drift_jacobian exactly at channel_degree = 0, where the channel is constant and the two questions coincide. Everywhere else the difference is the whole of the coordinate dependence described above, and it is the difference that decides whether a plan is being integrated against a decaying model or an extrapolation.

RBFKANLayer

Bases: Module

One Kolmogorov-Arnold layer with RBF edge functions (FastKAN-style).

Each input-output edge is a learnable 1D map phi(x) = sum_g c_g rbf_g(x) + w silu(x) over a fixed radial-basis grid; the output is bias_j + sum_i phi_{ji}(z_i). Each edge is an extractable 1D curve (interpretable) and cheap to evaluate.

chc.symbolic performs that extraction, and carries the two caveats it needs: an edge's intercept is a gauge (a constant moves freely between an edge and the bias, so only the total is identified), and this layer represents only additively separable functions -- for an interaction it silently returns the best additive fit, whose error is bounded BELOW by r^2 on [-r, r]^2.

KANResidual

Bases: Module

Learned residual r_θ(x, u) backed by a Kolmogorov-Arnold layer (RBF edges).

A drop-in ResidualModel alternative to MLPResidual; interpretable, and the same training / adjoint machinery applies unchanged.

SpectralResidual

Bases: Module

Residual that IS a circulant operator on a periodic grid, parameterised by its kernel.

The justifying plant is chc.transport.advection_diffusion_field -- a translation- invariant vector field on a periodic domain. Every such linear operator is circulant, every circulant is diagonalised exactly by the DFT, and its eigenvalues are the DFT of its first column (validation/spectral_circulant.mac STEP 1). So this backbone's hypothesis class is exactly the class of translation-invariant linear maps -- n parameters, not n^2.

Two properties no MLP backbone has, both structural rather than matters of fit quality:

• Translation equivariance to machine precision. r(roll(x, s)) == roll(r(x), s) for every shift, by construction. An MLP can be trained towards it and never attains it.
• An operator norm that is ATTAINED. operator_norm returns max_k |lambda_k|, the exact spectral norm, with the maximising Fourier mode as an explicit witness. Compare LipschitzResidual, whose constant comes from the Schur bound sigma_max(W) <= sqrt(||W||_1 ||W||_inf) -- valid, but with no witness and, measured, slack by more than an order of magnitude. Because gains multiply exactly under composition (Brahmagupta-Fibonacci, proofs/spectral_circulant.v gain_multiplies), the Result 28/30 rollout tube built on this constant is tight rather than merely valid.

The parameterisation is the kernel, not the symbol, deliberately: the map kernel -> circulant is a bijection, whereas a free half-spectrum carries two unidentified imaginary parts (at DC and, on an even grid, at Nyquist) that irfft silently discards -- the same gauge problem chc.symbolic prices for KAN edges.

A circulant is square, so out_dim is not a free parameter here; and the control channel is a circulant too, which is what lets it represent a NONLOCAL actuator. control_dim must therefore be 0 or state_dim.

symbol()

The state operator's eigenvalues as a half spectrum -- complex, length n//2 + 1.

operator_norm()

Exact spectral norm of the state operator: max_k |lambda_k|, attained not bounded.

GraphResidual

Bases: Module

Message-passing residual over a fixed graph -- a GNN backend for networked/spatial dynamics.

The state is n_nodes blocks of node_dim. Each node's update is an MLP of its features, the mean of its neighbours' encoded features (one message-passing round), and the control. It is permutation-equivariant and parameter-shared, learning a coupling a pointwise MLP re-learns per node. The adjacency is frozen (stop_gradient).

PortHamiltonianResidual

Bases: Module

Port-Hamiltonian residual x' = (J - R) grad H(x) + g(x) u -- passive, Lyapunov-stable.

Encodes energy + dissipation + a control port: J = A - A^T skew (lossless interconnection), R = L L^T >= 0 dissipation, H_theta a scalar energy network, g(x) the input matrix. With no input the energy obeys H' = dH . (J - R) dH = -dH . R dH <= 0 (the skew term dH . J dH = 0, dH := grad H). Cleanest when the known part is itself port-Hamiltonian (as DampedOscillator is), so known + residual stays port-Hamiltonian; as a pure additive correction the bound is on the residual's own contribution. Pure autograd (one scalar gradient), no inner solve, so it composes with the discrete/diffrax adjoints. Prefer over plain HNN/LNN (no port or dissipation).

energy decides what H' <= 0 is worth, and the difference is not cosmetic (Result 68, validation/convex_port_hamiltonian.mac). The inequality confines the state to { H <= H(x0) }, which bounds the state only if that set is bounded:

• "mlp" (the default, and what every model built before this option is) reads a tanh MLP out linearly, so H is BOUNDED in x. Every sublevel set above its supremum is the whole space, invariant_radius is inf, and grad H decays in all directions -- so the unforced flow stalls wherever it started rather than returning, and the energy has many near-equilibria. H' <= 0 is still true; it just certifies nothing about the state.
• "icnn" builds an input-convex energy with a quadratic floor, H(x) >= (eps/2)|x|^2, which makes { H <= c } a ball of radius sqrt(2c/eps) -- forward invariant, and a number the caller can print -- and makes the critical point unique, since the gradient of an eps-strongly convex function is injective.

The default stays "mlp" because changing it would silently move the semantics of every model already trained against this class; convex_energy_certificate measures the gap.

invariant_radius(level)

Radius of the ball the sublevel set { H <= level } sits inside.

sqrt(2 level/eps) for an input-convex energy, because H(x) >= (eps/2)|x|^2; inf for a plain MLP, whose bounded range leaves every high sublevel set equal to the whole space. The infinity is the honest answer, not a missing field.

energy_gradient(x)

grad_x H(x) -- the port-Hamiltonian co-energy vector.

structure_matrices()

The interconnection J (skew) and dissipation R (PSD) matrices.

LipschitzResidual

Bases: Module

Residual with a CERTIFIED Lipschitz constant L in [x, u] -- bounded off-support gain.

Each linear layer's weight is divided by a rigorous spectral-norm upper bound sigma_max(W) <= sqrt(||W||_1 ||W||_inf) (Schur), so every layer is <= 1-Lipschitz; with the 1-Lipschitz tanh the composition is <= 1-Lipschitz, and one learnable output scale sets the overall constant L = softplus(log_scale). This turns the soft Lipschitz penalty in chc.uncertainty into a by-construction invariant (invariants over guards): non-explosive rollouts and an honest L for the pessimism / regret bounds. The bound is on the residual's own contribution; total-system contraction still needs the known Jacobian. Pure matmuls, adjoint-friendly.

lipschitz_constant()

The certified Lipschitz constant L of r_theta with respect to [x, u].

ContractiveResidual

Bases: Module

Residual with a CERTIFIED negative one-sided Lipschitz (log-norm) mu < 0 -- contracting.

r(x,u) = -(rate + softplus(eta)) * x + rate * g(x,u) with g a Schur-normalized tanh-MLP (||dg/dx||_2 <= 1, as in LipschitzResidual). By subadditivity of the log-norm and mu_2(-D) = -min_i D_ii, the state-Jacobian obeys mu_2(J_r) <= -(rate + min softplus(eta)) + rate = -min softplus(eta) < 0 -- a certified contraction rate, not a penalty. Unlike the non-negative ||.||-Lipschitz L of LipschitzResidual (rollout bound e^{L*T}), a negative log-norm mu = -c gives a UNIFORMLY BOUNDED CONTINUOUS-time radius eps/c (no e^{L*T}). For EXPLICIT EULER the discrete contraction factor is q = sqrt(1 + 2*mu*dt + L^2*dt^2) (NOT 1+mu*dt; L=lipschitz_constant() the full Lipschitz constant), so contraction needs the step dt < 2c/L^2 (sufficient); under it the rollout radius stays bounded and tends to eps/c as dt -> 0 (see chc.uncertainty.contractive_rollout_bound; Rocq contractive_euler.v). mu bounds only the residual's contribution; total-system contraction needs mu_2(J_known)+margin < 0 (a gate to check offline). Pure matmuls; the L2/eigen bound is exact, general-n verification would use the L-inf (Gershgorin) variant.

contraction_rate()

The certified contraction rate c = |mu| = min_i softplus(eta_i) > 0.

lipschitz_constant()

Full Lipschitz constant L of r in x (needed for the explicit-Euler CFL).

||J_r|| <= ||-(rate + softplus(eta)) I|| + rate*||dg/dx|| <= (rate + max softplus(eta)) + rate since g is 1-Lipschitz -- so L = 2*rate + max softplus(eta). Sets the Euler contraction step dt < 2c/L^2 in chc.uncertainty.contractive_rollout_bound.

PortHamiltonianCertificate dataclass

Numeric evidence that a PortHamiltonianResidual is passive (Lyapunov-stable).

ConvexEnergyCertificate dataclass

Result 68: a DECREASING energy is not a BOUNDED one, and only a coercive one gives a radius.

H' <= 0 holds for both energies -- it is an identity of the port-Hamiltonian structure, not a property of the network -- so the inequality cannot tell them apart. What separates them is whether the sublevel set it confines the state to is bounded.

LipschitzCertificate dataclass

Numeric evidence that a LipschitzResidual respects its certified constant L.

DampingInjectionCertificate dataclass

Evidence that damping injection makes the port-Hamiltonian residual dissipate energy.

SpectralResidualCurve dataclass

Evidence that a circulant backbone earns its place on a translation-invariant plant.

control_affine_features(x, degree)

Monomials of a single state up to total degree, bias first.

Owned here rather than shared with chc.causal's batch version because it indexes the parameters of ControlAffineResidual: the estimator that fits those parameters (chc.dynamics_id.fit_causal_residual) has to use this exact basis, and a second implementation is a silent way for the two to disagree about what a coefficient means.

port_hamiltonian_certificate(seed=0, state_dim=3, control_dim=1, n=64)

Sample states and check passivity: skew J adds no energy, R >= 0, and H' <= 0.

convex_energy_certificate(seed=0, state_dim=2, control_dim=1, convexity=0.25, n_start=49, spread=3.0, dt=0.05, n_step=4000)

RESULT 68 (A20) -- the port-Hamiltonian residual needs a COERCIVE energy, not a falling one.

Derived in validation/convex_port_hamiltonian.mac, proved in proofs/convex_port_hamiltonian.v. port_hamiltonian_certificate already checks the three things that make the structure passive: the skew form vanishes, R >= 0, and H' <= 0. All three hold for ANY energy network, because they are identities of (J - R) grad H -- which is exactly why they cannot distinguish a useful energy from a useless one.

What H' <= 0 buys is confinement to { H <= H(x0) }, and that is a bound on the state only when the set is bounded. A tanh MLP read out linearly is bounded in x, so its high sublevel sets are the whole space and the confinement is vacuous; and because its gradient decays in every direction it has many critical points, so the "Lyapunov function" certifies convergence to something without saying where. An input-convex energy with a quadratic floor gives H(x) >= (eps/2)|x|^2, hence the forward-invariant ball |x(t)| <= sqrt(2 H(x0)/eps) -- a number the caller can print -- and a unique critical point, because grad H of an eps-strongly convex function is injective.

The arms, each able to fail: the passivity identity by autodiff against the algebra; the range of each energy on a box and on a box ten times larger, which saturates for one and grows quadratically for the other; the smallest Hessian eigenvalue over states spanning that whole sweep; the quadratic floor on the outermost shell; the realised excursion of the unforced flow against the predicted radius; and the number of distinct near-equilibria reached from a grid of starts.

lipschitz_certificate(seed=0, state_dim=3, control_dim=1, out_dim=3, n=200)

Sample pairs and confirm ||r(a) - r(b)|| <= L * ||a - b|| for the certified L.

damping_injection_certificate(seed=0, state_dim=3, kappa=1.0, horizon=40, dt=0.02)

Apply u = -kappa*g^T dH to a PortHamiltonianResidual and confirm the energy decays.

The shipped port_hamiltonian_certificate checks only AUTONOMOUS (u=0) passivity. Here the control adds dissipation: the closed-loop rate Hdot = dH^T (J-R) dH - kappa*(g^T dH)^2 = -dH^T R dH - kappa*(g^T dH)^2 <= 0 (skew term exactly 0), so H strictly decays along the closed loop. Machine-checked in proofs/port_hamiltonian_lyapunov.v (scalar/2x2 case). HONEST SCOPE: this is the RESIDUAL field's closed loop, and H (an MLP) is not enforced positive-definite, so it decays to argmin H, not necessarily 0 -- convergence to the origin needs a PD energy; for the full hybrid f_known + r the whole-system claim needs f_known matching or energy-orthogonality (see the proof header).

fit_spectral_residual(model, xs, us, ys)

Fit a SpectralResidual by CLOSED-FORM least squares, per Fourier mode.

A circulant is linear in its kernel, so gradient descent is the wrong tool: in the Fourier basis the design decouples completely and each bin k is its own two-parameter complex regression y_k = lambda_k x_k + mu_k u_k. One 2x2 normal system per bin, solved exactly -- no learning rate, no step count, no seed. Same argument as the ordinary least squares behind chc.symbolic.extract_symbolic, and it stops being available the moment a nonlinearity is stacked on top.

That this matters is a MEASUREMENT, not a preference. The circulant kernel of a spectral diffusion operator has entries of order nu*n^2/L^2 -- 134.8 at the certificate's settings -- and a first-order optimiser starting near zero needs O(nu*n^2/lr) steps to travel that far. Fit by Adam on the same budget as the MLP, this backbone loses; fit properly, it wins by four orders of magnitude. The optimiser, not the hypothesis class, decided the naive comparison.

irfft projects away any imaginary residue at the DC and Nyquist bins. That is the correct projection rather than a loss: for real data those bins of x, u and y are all real, so the exact solution there is real too, and what is discarded is float noise.

spectral_residual_certificate(*, seed=0, n=64, length=1.0, speed=0.8, viscosity=0.01, n_train=256, n_test=128, steps=400, lr=0.02, rollout=40, dt=0.002)

Does a circulant backbone beat an MLP on the plant that justifies it -- and how, exactly.

This is the backbone's kill-criterion: if the MLP matches it on its own home ground, the backbone is deleted and the finding recorded.

The plant is chc.transport's periodic advection-diffusion field plus a Gaussian-smoothed control channel -- translation-invariant by construction, so the truth lies IN the circulant hypothesis class. That makes the fit-quality arms a test of inductive bias rather than of approximation power (the MLP carries ~130x more parameters and is not handicapped), and it makes the EQUIVARIANCE arm the decisive one: that gap cannot be closed by more training.

Each backbone is fit the way it should be: the circulant by the closed-form per-mode least squares of fit_spectral_residual, the MLP by Adam. spectral_adam_test_mse records what happens when the circulant is instead fit by Adam on the MLP's own budget -- it loses, and the reason is arithmetic rather than structural (see that function's docstring).