Stabilize generic integrated process covariance over long gaps

[seldak]

For IntegratedMatern52(scale=1, sigma=1), the generic process covariance at dt=300 is finite but inaccurate: the integral variance is 114020.673294, while independent quadrature of the physical impulse response gives 711.341753. This patch brings the result into agreement with quadrature. A finiteness check alone does not detect the discrepancy.

The following reproduces the result on smolgp 0.3.2 at 22bc8c19232140aca4f7903a37b687e073a70e5e, with JAX x64 on CPU:

import jax
jax.config.update("jax_enable_x64", True)
import smolgp

kernel = smolgp.kernels.IntegratedMatern52(scale=1.0, sigma=1.0)
print(kernel.process_noise(0.0, 300.0)[-1, -1])
# Base revision: approximately 114020.673294
# Physical impulse-response quadrature: approximately 711.341753

Case	Before	After / independent reference
IntegratedMatern52, dt=300, integral variance	114020.673294 (finite)	711.341753
IntegratedExp, dt=700, integral variance	1397 (finite)	1397
IntegratedExp, dt=1000, integral variance	nonfinite	1997
Two-instrument IntegratedExp GP, long-gap log likelihood	NaN	-4.44626673600823

The new regression integrates exp(F*t) L Qc L.T exp(F.T*t) directly with SciPy quadrature, so its reference calculation does not reuse a Van Loan block.

Change

The generic integral covariance previously extracted F3/H2/K1 from a full-gap auxiliary exponential containing both F and -F. Those auxiliary blocks can overflow or lose accuracy even when the physical covariance remains representable.

The replacement changes both discretization and assembly: append one integral state to the base system, compute its short-step transition and covariance, then compose the physical pairs with the classical identities A(2h)=A(h)@A(h) and Q(2h)=Q(h)+A(h)@Q(h)@A(h).T. Extract the cross-covariance and integral variance from this augmented covariance. The existing process_noise method replicates those blocks across instruments.

The appended state integrates base coordinate 0, matching the existing integrated model's design, observation and transition conventions. The public VanLoan helper, Phibar_from_VanLoan, and the analytic SHO path are unchanged.

Validation and limits

python -m pytest tests/test_long_gap_noise.py tests/test_proofs.py tests/test_integrated.py -q passes 23 tests. The one warning is the existing instrument-count reinitialization test. New coverage includes the finite-but-inaccurate Matérn case, independent physical quadrature for Exp/Matérn32/Matérn52, analytic OU values and parameter gradients, JIT/vmap, zero-time behavior and derivative, semigroup composition, multiple integral states, and a full two-instrument likelihood.

These regressions use x64 reference calculations. The CI test workflow at the base revision also sets JAX_ENABLE_X64=1.

Test environment: Python 3.12.14, JAX/jaxlib 0.11.1, NumPy 2.5.3, SciPy 1.18.1, Equinox 0.13.8; tinygp pinned at bea3c43e47c339cb20bc14b91cbd0e55a592f7e6.

The helper permits at most 32 doublings by default; negative/nonfinite elapsed time or a required scaling beyond that budget returns NaN matrices. A static loop preserves reverse-mode differentiation, and the norm-based scaling decision is stopped from differentiation. This is an explicit numerical policy, not a guarantee for arbitrary ill-conditioned or physically overflowing systems.