Performance: continuous-time reservoirs (PR4 of #397)

[Saswatsusmoy]

Summary

Tracking issue for performance work on continuous-time reservoirs (ContinuousESN / SciMLProblemReservoir / RCODEReservoirExt), corresponding to PR4 of the SciML fellowship plan on #397.

PR1–PR3 landed the correct continuous core (#446, #450, #456). This issue tracks data-driven optimizations so each change can land as a small, reviewable PR with a before/after measurement.

Profiling-only harness (no package code changes): benchmarks/continuous_profile/ (local for now; can be proposed separately if useful upstream).

Motivation

On Apple M4 / Julia 1.12.5 / single-thread OpenBLAS, with reltol=1e-6, abstol=1e-8, Float64, use_bias=true, seed 17:

Workload	Setup	Median
W3	ContinuousESN collectstates, N=500, T=1000, Tsit5	~13.5 s
W4	same, Euler dt=1	~42 ms
W5 / W5t	Rodas5 vs Tsit5, N=100, T=200	~1.57 s / ~0.12 s
W7	AR predict, N=500, 250 steps, Tsit5	~0.76 s
W10b	discrete ESN collectstates, N=500, T=1000	~42 ms

Key ratios:

• Tsit5 / Euler (matched N=500, T=1000): ~320×
• Continuous Tsit5 / discrete ESN (N=500): ~350×
• Rodas5 / Tsit5 (matched N=100, T=200): ~13×

Stage timers show solve ≈ 100% of wall time on every continuous path. Grid construction, ZOH, NamedTuple merges, and modifiers are negligible once warm.

Profiling findings (granular)

1. Adaptive stepping dominates, not the RHS wrapper

• W3 (Tsit5, N=500, T=1000): 330 537 RHS evaluations (nf), 21 262 accepted steps, 33 827 rejected.
• W4 (Euler, same problem): 1 001 RHS evaluations.
• Default RHS _continuous_esn_rhs! is already 0 allocations / call (mul! + in-place broadcast). Re-optimizing matvec alone will not hit a 3× target if NFE counts stay the same.

2. Teacher-forced path is one long solve

collectstates builds a ZOH input, one ODEProblem, then solve(...; saveat=sample_ts). Flamegraphs for W3 land in:

• OrdinaryDiffEqTsit5 perform_step!
• residual / error-control helpers
• tanh / FastBroadcast
• mul! (matvec) at larger N

3. Autoregressive path rebuilds the problem every step

W7: 250 independent short solves. Per-step costs (alloc audit):

Stage	Bytes / step
solve	~76 KB
ODEProblem build	~3 KB
ConstantInterpolation input	~0.7 KB

Profile.Allocs is dominated by fresh ODEProblem / ODEIntegrator / ConstantInterpolation types each step — classic integrator reuse surface (init once + reinit!).

4. Implicit solvers without sparsity are expensive

W5 Rodas5 reports 2978 Jacobians (njacs) and ~334k RHS calls at N=100, T=200, with dense finite-difference Jacobian cost. Sparse jac_prototype from rand_sparse is the natural fix.

5. Allocation audit (secondary to wall time)

Site	Notes
reduce(hcat, sol.u) in _sample	~half of oneshot bytes on W3; wall time tiny when warm
solve solution / integrator buffers	other half on W3
AR per-step rebuild	main AR allocation driver

ZOH returns a true SubArray view (no vector copy of the input column).

6. Train is collectstates-bound

W8 train!: 99.8% in collectstates; ridge ~18 ms.

Non-goals (this issue)

• Continuous predict warmup / cold-start (deferred separately on feat(models): ContinuousESN #456 discussion).
• LSM / DDE stretch (PR5–PR6 of Continuous-time reservoir computing #397).
• Merging the local profiling harness into the package tree (optional follow-up; keep root Project.toml free of BenchmarkTools etc.).
• Changing the continuous ODE form (eq. (5) / Lukoševičius 2012 is settled).

References

• Parent: Continuous-time reservoir computing #397
• Core types: feat: SciMLProblemReservoir core types + collectstates dispatch helper #446
• Extension: feat: continuous-time reservoirs via RCODEReservoirExt #450
• ContinuousESN model: feat(models): ContinuousESN #456
• Local profile report (reproducible via harness): benchmarks/continuous_profile/results/REPORT.md
• Canonical ODE: Lukoševičius 2012 §3.2.6 eq. (5)

Environment used for numbers above

• Julia 1.12.5, -O3 --check-bounds=no
• BLAS: OpenBLAS ILP64, 1 thread
• CPU: Apple M4
• Package git SHA at profile time: 2f47eeb5a4b3876c2f72b032217434fc84aa007c

cc @MartinuzziFrancesco

[hiSandog]

Validation would be most useful if it reports the first incompatible dimensions together with the relevant configuration names, rather than allowing a later matrix operation to fail generically. Tests could cover each invalid relationship independently so users can tell which input needs correction.

[Saswatsusmoy]

Hey @MartinuzziFrancesco — short update before any performance PRs.

Profiling (no package speed changes yet): on ContinuousESN collectstates N=500 / T=1000, wall time is ~100% in solve. Adaptive Tsit5 is ~13s with ~3e5 RHS evals; grid Euler (dt=1) is ~40ms (~1e3 RHS), same ballpark as discrete ESN. RHS is already 0-alloc. AR predict still does continuous free-run but rebuilds the problem every step. Train is collectstates-bound.

Constraint: keep real SciML/ODE continuous RC (eq. 5, continuous (u(t)), terminal sampling) — no silent discrete bypass.

Need your call:

1. Keep adaptive as default and add an opt-in fixed-step solve path for uniform grids, or change defaults?
2. Next: AR integrator reuse (init/reinit!) vs warmup (separate PR as you preferred) vs docs/NFE dig only for now?

Happy to follow whichever order you prefer.