Capabilities and limitations

Capabilities and limitations#

A candid summary of what lyapax does and does not do, so you can decide quickly whether it fits your problem.

What lyapax can do#

  • Full or partial Lyapunov spectra for ODEs. Track all \(d\) exponents or only the leading \(k \le d\), for a hand-written JAX right-hand side or a ModelSpec/build_jax_dfun-compiled symbolic one; cost scales with \(k\), not \(d\) (see How lyapax is implemented).

  • Coupled networks. Arbitrary topology (a weight matrix), any coupling rule - built-in linear / sigmoidal / Kuramoto, or a user-written callable - with fixed-step Euler / Heun / RK4 / RK6 integration.

  • Fixed-delay DDE networks. Both a single uniform delay shared by every edge (network_dde_problem(..., tau=...)) and a genuine per-edge heterogeneous delay matrix via the lower-level simulator path - the latter currently only with the built-in linear coupling (see the gaps below).

  • Grid-snapped or Hermite-interpolated DDE history - fast-and-simple versus exact-\(\tau\), higher-order-accurate; see How lyapax is implemented for the tradeoff.

  • Discrete chaotic maps (logistic, Hénon, or any user map) - same engine, dt=1.0 per iterate, no integrator involved.

  • Adaptive-step ODE integration (optional adaptive extra, lyapax.adaptive.diffrax_adaptive_step) - drops into ode_problem alongside the fixed-step builtins with no other code changes; ode_problem only (not network_problem, not DDEs), see the gaps below.

  • Batched parameter / initial-condition sweeps via jax.vmap (sweep_lyapunov_spectrum): a whole parameter grid as one XLA call.

  • Kaplan-Yorke (Lyapunov) dimension (lyapax.core.kaplan_yorke_dimension)

    • a pure post-processing readout of an existing LyapunovResult.exponents, no extra simulation. Raises rather than silently understating the answer if given a partial spectrum (d_total=) whose cumulative sum never crosses zero within the tracked exponents.

  • Differentiating an exponent w.r.t. a system parameter (jax.grad/jax.jacfwd(lyapunov_spectrum(...).exponents[i])) - reliable for non-chaotic or short-horizon systems (e.g. gradient-based tuning of a parameter toward a target exponent); see the caveat below for chaotic trajectories.

  • Transparent GPU execution - JAX picks the backend; no code changes. Whether the GPU is faster depends on problem size.

  • A validation suite anchored to independent sources - exact eigenvalues, structural invariants, published literature values; see Validation: how lyapax’s results are checked.

What lyapax cannot do#

  • No stiff (implicit) ODE integration. The built-in integrators are fixed-step (Euler / Heun / RK4 / RK6); adaptive-step explicit Runge-Kutta is available for single, uncoupled ODEs (ode_problem) via the optional adaptive extra (lyapax.adaptive.diffrax_adaptive_step, backed by diffrax) - not for network_problem (rejected outright, even for a non-delayed network) or DDEs, and there is still no implicit solver, so genuinely stiff systems are out of scope regardless. The computed exponents are those of the numerical time-dt map (or, for the adaptive integrator, of the accepted-step sequence under a given rtol/atol), not the exact continuous flow - checking convergence in dt (or tolerance) is the caller’s responsibility (the gallery’s speed-and-accuracy and adaptive-ODE examples show how). The adaptive integrator’s internal step-size control is a dynamic-trip-count while_loop, so differentiating through it needs jax.jacfwd (forward-mode); jax.grad/jax.jacrev (reverse-mode) do not work through it.

  • DDE delays must be known and fixed. No state-dependent or distributed delays. Grid-snapped mode further rounds \(\tau\) to an integer multiple of dt; only interpolate=True removes that restriction, and only on the uniform-delay path.

  • Per-edge delays cannot be combined with a custom coupling rule. Heterogeneous per-edge delay matrices currently work only with the built-in linear coupling; a delayed per-edge Kuramoto network, for example, would need an edge-aware coupling signature that does not exist yet.

  • No stochastic / noise-driven Lyapunov exponents. Noise injection is a deliberate non-goal: it is non-smooth and not meaningfully compatible with the deterministic tangent-space spectrum computed here.

  • No PDEs or spatiotemporal chaos. State is always a finite-dimensional vector (or an (n_state_vars, n_nodes) array for networks); there is no spatial discretization support.

  • history column ordering is not stable for near-degenerate exponents. Columns are ordered once, by the final row; exponents that nearly cross can swap order mid-run. See LyapunovResult.history’s docstring.

  • dfun_str model specs are compiled with exec(). ModelSpec string right-hand sides are code-generated without sanitization - fine for specs you write yourself, not safe for specs built from untrusted input.

  • Gradients through a chaotic trajectory are numerically useless, even though they don’t raise an error. jax.grad/jax.jacfwd of an exponent w.r.t. a parameter differentiate through the whole unrolled state trajectory, so on a chaotic system the result inherits that trajectory’s own exponential sensitivity - it grows roughly like exp(lambda_max * horizon) and is unreliable well before it overflows (a known chaotic-sensitivity-analysis phenomenon, not a lyapax bug - see lyapax.core’s module docstring). Safe for non-chaotic/short- horizon systems; always sanity-check a chaotic-system gradient against finite differences before trusting it.

Adaptive integration is not a speed optimization#

Measured, not assumed: at matched accuracy (not matched dt/rtol, which isn’t a fair comparison across integrators with different meanings of “step size”), diffrax_adaptive_step was 2-4x slower than fixed-step rk4/rk6 on a small chaotic system (Lorenz), up to ~4.5x slower on GPU for that same small system, and still ~4x slower even on a relaxation oscillator (Van der Pol) with sharply time-varying stiffness - the scenario adaptive stepping exists for. The overhead ratio does shrink at large state dimension (d=1500-3000), approaching rough parity, but never produced a clear win in any case tested.

Root cause: the internal accept/reject jax.lax.while_loop pays real per-step overhead (PID controller update, solver-step bookkeeping) that a fixed-step raw step doesn’t, and avoiding extra fixed-step steps only pays off in wall-clock terms if each one was expensive to begin with - which isn’t the case for any system small enough to fit the exact/published references this page and Validation: how lyapax’s results are checked rely on.

Use the adaptive integrator for tolerance-driven accuracy control (no manual dt-convergence sweep) or because it’s the only forward-mode- differentiable adaptive option here - not for throughput.