Skip to content

Work · Lyapunov time

Lyapunov time

A forecast horizon of "347 steps" means nothing on its own: it depends on the integration step, the system, and the sampling rate. resdag.metrics reports horizons in Lyapunov times instead — the natural clock of a chaotic attractor — so a number is comparable across systems, integration steps, and papers.

Why Lyapunov times

On a chaotic system, two trajectories that start infinitesimally close separate on average like \(e^{\lambda_\mathrm{max} t}\), where \(\lambda_\mathrm{max}\) is the largest Lyapunov exponent. One Lyapunov time

\[ t_\lambda = \frac{1}{\lambda_\mathrm{max}} \]

is the timescale over which that separation grows by a factor of \(e\). It is the system's own unit of "how far ahead is hard", so it makes forecast horizons directly comparable: a valid prediction of \(6\,t_\lambda\) on the Lorenz system and \(6\,t_\lambda\) on Rössler represent the same predictive skill, even though the raw step counts differ by an order of magnitude.

A horizon in steps converts to Lyapunov times by the headline formula

\[ \text{VPT}_{\Lambda t} = \text{steps} \cdot dt \cdot \lambda_\mathrm{max}, \]

where \(dt\) is the integration step. This is steps_to_lyapunov; its inverse is lyapunov_to_steps, and steps_per_lyapunov_time gives \(1/(dt\,\lambda_\mathrm{max})\) — how many samples span one \(t_\lambda\).

Canonical exponents

For resdag's built-in chaotic generators the literature \(\lambda_\mathrm{max}\) values are tabulated in CANONICAL_LAMBDA_MAX (and, keyed by the Mackey-Glass delay \(\tau\), in MACKEY_GLASS_LAMBDA):

Generator \(\lambda_\mathrm{max}\) (per unit \(t\)) \(dt\)
lorenz (σ, ρ, β = 10, 28, 8/3) 0.9056 0.02
rossler (a, b, c = 0.2, 0.2, 5.7) 0.0714 0.05
mackey_glass (τ = 17) 0.0086 1.0

The exponent is in inverse natural-time units, so multiply by \(dt\) to turn a step count into Lyapunov times. Passing a dataset= name to the metrics below looks these up for you.

Measuring forecast quality

The Valid Prediction Time (VPT) is the standard forecast-skill metric on chaotic systems: the number of contiguous steps a closed-loop forecast tracks the truth before its normalised error (NRMSE) first crosses a threshold. valid_prediction_time returns it in raw steps by default, or in Lyapunov times with units="lyapunov" — pass dt and lambda_max explicitly, or a known dataset= name to infer both.

import torch

import resdag as rd
from resdag.datasets import lorenz
from resdag.metrics import reservoir_lyapunov, steps_to_lyapunov, valid_prediction_time
from resdag.utils import prepare_esn_data

torch.manual_seed(0)

# A canonical Lorenz-63 series at dt = 0.02.
series = lorenz(6000)  # (1, 6000, 3)
warmup, train, target, f_warmup, val = prepare_esn_data(
    series, warmup_steps=500, train_steps=4000, val_steps=1000
)

# Train an Ott-augmented ESN in one algebraic pass.
model = rd.ott_esn(reservoir_size=400, feedback_size=3, output_size=3, seed=0)
rd.ESNTrainer(model).fit((warmup,), (train,), {"output": target})

# Forecast aligns one-to-one with the validation split.
preds = model.forecast(f_warmup, horizon=val.shape[1])

# VPT in raw steps and in Lyapunov times (dt / lambda_max inferred from the name).
vpt_steps = valid_prediction_time(preds, val)
vpt_lyap = valid_prediction_time(preds, val, units="lyapunov", dataset="lorenz")
print(f"VPT = {vpt_steps:.0f} steps  =  {vpt_lyap:.2f} Lyapunov times")

# The same conversion done by hand from the canonical exponent.
print("check:", round(steps_to_lyapunov(vpt_steps, dt=0.02, lambda_max=0.9056), 2))

# The reservoir's own maximal Lyapunov exponent (edge-of-chaos diagnostic).
lam = reservoir_lyapunov(model, f_warmup, horizon=500)
print(f"reservoir lambda_max = {lam:.4f} per step  (< 0 => ESP-satisfying)")

Running it prints (exact numbers depend on the BLAS backend):

VPT = 347 steps  =  6.28 Lyapunov times
check: 6.28
reservoir lambda_max = -0.5955 per step  (< 0 => ESP-satisfying)

The horizon of ~6 Lyapunov times is the meaningful figure — you can quote it alongside any other chaotic-system result and it stays comparable regardless of the integration step used.

The reservoir's own exponent

reservoir_lyapunov estimates the maximal Lyapunov exponent of the trained reservoir's own driven dynamics — a distinct quantity from the attractor's \(\lambda_\mathrm{max}\) above. It runs Benettin's algorithm with QR reorthonormalisation of the analytic tangent map along an input-driven trajectory, in native torch with no optional dependency.

It is the standard edge-of-chaos diagnostic. The Echo State Property requires the driven dynamics to contract — trajectories from different initial states must converge — so a healthy ESN has a maximal exponent that is negative but close to zero (it sits just below the edge). A positive value signals a reservoir whose state does not forget its initial condition, which breaks the ESP. The returned exponent is in inverse-step units (per reservoir update); multiply by \(1/dt\) to compare it with a natural-time attractor exponent.

lam = reservoir_lyapunov(model, f_warmup, horizon=500)   # scalar, k=1
spectrum = reservoir_lyapunov(model, f_warmup, horizon=500, k=4)   # top-4 as (4,) array

Pass a single ESNLayer or a model containing exactly one ESN reservoir; the reservoir must use the tanh activation (the only one whose analytic Jacobian is wired in).

Arbitrary systems — the dynamics extra

For systems that are not in the canonical tables, largest_lyapunov dispatches on its argument:

  • a generator name ("lorenz", "rossler", "mackey_glass") returns the canonical constant — no dependency;
  • a measured series gives a noisy, data-driven diagnostic estimate (opt in with from_="data", and pass dt);
  • a tsdynamics System gives a high-accuracy Benettin spectrum estimate.
from resdag.metrics import largest_lyapunov

largest_lyapunov("lorenz")                 # 0.9056  — canonical, no dependency
largest_lyapunov("mackey_glass", tau=17)   # 0.0086  — canonical, no dependency

Optional dependency: resdag[dynamics]

The series- and System-based paths of largest_lyapunov require the optional tsdynamics package, which needs Python ≥ 3.12:

pip install "resdag[dynamics]"   # or: uv sync --extra dynamics

Every tsdynamics import is lazy and guarded, so this extra is only needed to estimate an exponent for an arbitrary system. The Lyapunov-time conversions, the canonical constants, valid_prediction_time, and reservoir_lyapunov all work with no extra installed.

Next