Statement of need/field

Statement of need/field#

Why lyapax exists, and how it relates to the other tools that compute Lyapunov spectra. This mirrors the “Statement of need” and “State of the field” sections of the JOSS paper (joss/paper.md), kept here so the motivation is visible from the documentation site as well.

Statement of need#

Computing a Lyapunov spectrum reliably requires more than the textbook formula: one must integrate the state and a set of orthonormalized tangent vectors together, re-orthonormalize periodically (QR), and average the logarithmic growth of each direction over a long horizon in a numerically stable way. Implementing this correctly for a new model - especially a high-dimensional network or a delay system - is error-prone, and existing tools each cover only part of the space.

lyapax fills a specific gap: a Lyapunov-spectrum engine that

  1. works on any differentiable JAX map, so the model is ordinary Python/jax.numpy code rather than a symbolic expression that must be built and compiled ahead of time;

  2. is matrix-free, propagating only the k tangent directions actually requested via jax.jvp instead of forming a dense \(d \times d\) Jacobian, so partial spectra of high-dimensional systems are cheap - and, unlike a symbolic Jacobian, this cost is paid at run time, not compile time (see the note on network size below);

  3. treats parameter sweeps as a first-class, vectorized operation through jax.vmap; and

  4. runs unchanged on CPU or GPU, with no code change to move a model from one to the other.

It is aimed at researchers in nonlinear dynamics, computational neuroscience, and physics who need Lyapunov spectra of custom ODE, map, or DDE models - particularly coupled oscillator networks and delay-coupled systems - without implementing or hand-linearizing the tangent-propagation machinery themselves.

State of the field#

  • jitcode and jitcdde compile a symbolic right-hand side to C and are the established Python references for ODE and DDE Lyapunov spectra respectively. They are fast and mature once built, but the model must be expressed as a symbolic expression rather than arbitrary code, neither has a GPU backend, and their integrator choice is fixed to SciPy’s adaptive solvers (dopri5, RK45, dop853, RK23, BDF, LSODA, Radau, vode) - there is no fixed-step classical RK4/Heun option, so a step-for-step same-algorithm comparison against lyapax’s own fixed-step integrators isn’t possible without patching the library (see Benchmarks’s performance section). The symbolic-then-compile step also does not scale to larger, densely coupled networks: building and compiling jitcode for a fully-connected 50-node Kuramoto network alone took 54 seconds in a direct test, before a single integration step ran - the symbolic differentiation needed for the tangent/Jacobian machinery grows with the number of nonzero coupling terms (~2,450 of them at just 50 dense nodes), so this cost only gets worse at the network sizes (hundreds to thousands of nodes) where lyapax’s GPU/vmap path is actually meant to be used. lyapax pays no such compile-time symbolic-differentiation cost at any network size - jax.jvp computes the tangent direction at run time from ordinary Python code.

  • ChaosTools.jl, part of DynamicalSystems.jl, is a mature, broad-scope Julia toolbox for nonlinear dynamics, including Lyapunov spectra, but it lives outside the Python/NumPy ecosystem and is not designed around autodiff or GPU execution. Its OrdinaryDiffEq.jl backend does expose fixed-step classical methods (RK4(), Vern6()), which does let a genuine same-algorithm CPU comparison be constructed (see Benchmarks) - but lyapunovspectrum itself has no GPU path at all: Julia’s GPU story for ODEs (DiffEqGPU.jl) is built around parallelizing an ensemble of independent trajectories on the device, not accelerating one small system’s tangent propagation, which doesn’t match how a single Lyapunov-spectrum run is structured. Getting a ChaosTools.jl-based computation onto a GPU at all would mean bypassing its own API and hand-writing the Benettin/QR loop directly against CUDA kernels.

  • Ad hoc SciPy-based approaches (custom scripts using finite differences, hand-derived Jacobians, or manually coded variational equations) are common in practice but are one-off, unvalidated, and not reusable across projects.

  • Diffrax provides differentiable ODE/SDE/CDE solvers natively in JAX, including adaptive-step and stiff methods, but it is a general differential-equation solver library, not a Lyapunov-spectrum package; it does not itself provide tangent-space QR renormalization, partial-spectrum tracking, or DDE ring-buffer handling.

lyapax is, to the author’s knowledge, the first JAX-native package purpose-built for Lyapunov spectra that accepts ordinary user-written JAX functions (rather than symbolic expressions) and targets the autodiff/JIT/GPU workflow directly, while remaining a narrowly scoped complement to - not a replacement for - these tools. A concrete, reproducible cross-check of lyapax against jitcode, jitcdde, and ChaosTools.jl on shared benchmark systems is in Benchmarks.