A Few Months in Reservoir Computing

I spent the last few months contributing to ReservoirComputing.jl, the reservoir computing package inside SciML. Four models merged: SVESM, LIFESN, ResESN, and RMNResESN. This post covers what reservoir computing actually is, how the package works under the hood, and what each of those models does.

If you've never heard of reservoir computing: it's a way to train recurrent neural networks without backpropagating through time. The recurrent layer is random and frozen. You only train the readout. Training is a single least-squares solve — no epochs, no learning rate, no gradient descent. The whole thing runs on CPU, fits in L2 cache, and works on chaotic time series. It's weird and it's elegant and it's been around since 2001. Here's how it works.

Echo state networks

The canonical model. Input \(\mathbf{u}(t) \in \mathbb{R}^D\), reservoir state \(\mathbf{r}(t) \in \mathbb{R}^N\) where \(N\) is large — typically a few hundred to a few thousand — and output \(\hat{\mathbf{y}}(t) \in \mathbb{R}^O\).

$$\begin{aligned} \tilde{\mathbf{r}}(t) &= \tanh(W_{\text{in}} \mathbf{u}(t) + W_r \mathbf{r}(t-1) + \mathbf{b}) \\ \mathbf{r}(t) &= (1 - \alpha)\mathbf{r}(t-1) + \alpha \tilde{\mathbf{r}}(t) \\ \hat{\mathbf{y}}(t) &= W_{\text{out}} \mathbf{r}(t) \end{aligned}$$

\(W_{\text{in}} \in \mathbb{R}^{N \times D}\) is the input weight matrix. \(W_r \in \mathbb{R}^{N \times N}\) is the recurrent weight matrix — sparse, random, fixed forever. \(\mathbf{b}\) is the bias. \(\alpha \in (0, 1]\) is the leak rate: \(\alpha = 1\) means no leak (the original Jaeger formulation), \(\alpha < 1\) introduces exponential smoothing with timescale \(1/\alpha\). \(W_{\text{out}} \in \mathbb{R}^{O \times N}\) is the readout — the only thing that gets trained.

Randomize the hard part. Train the easy part. The reservoir does the nonlinear projection. You just learn to read it.

The echo state property

The ESP is what makes this well-posed. The reservoir state must asymptotically forget its initial condition. For any two initial states driven by the same input:

$$\lim_{t \to \infty} \|\mathbf{r}_1(t) - \mathbf{r}_2(t)\| = 0$$

Without input and bias, the reservoir map reduces to:

$$\mathbf{r}(t) = (1 - \alpha)\mathbf{r}(t-1) + \alpha \tanh(W_r \mathbf{r}(t-1))$$

Since \(\|\tanh(\mathbf{x})\| \le \|\mathbf{x}\|\), if \(\|W_r\|_2 \cdot \alpha < 1\) the map is a contraction and the ESP holds by Banach's fixed-point theorem. That's the operator-norm sufficient condition. In practice people use the spectral radius: \(\rho(W_r) < 1\) is necessary (proved by Jaeger), and \(\rho(W_r) \in [0.8, 1.2]\) works fine because the input drive and bias can stabilize reservoirs above 1.0. The package enforces this at init:

$$W_r := W_r \cdot \frac{\rho_{\text{target}}}{\rho(W_r)}, \quad \rho(W) = \max\{|\lambda_i| : \lambda_i \in \text{eig}(W)\}$$

For continuous-time reservoirs the condition tightens. You need \(\|W_r\|_2 < \alpha\), not \(\rho(W_r) < \alpha\). Spectral radius is purely a discrete-time artefact. Operator 2-norm is the actual continuous bound.

Training

Training an ESN is almost disappointingly simple. Drive the reservoir with your input sequence \(\{\mathbf{u}(1), \ldots, \mathbf{u}(T)\}\). Collect the reservoir states into a matrix \(R = [\mathbf{r}(1) \cdots \mathbf{r}(T)] \in \mathbb{R}^{N \times T}\). Throw away the first \(W\) columns as washout — the reservoir needs a few steps to forget its initial state. Then solve:

$$W_{\text{out}} = \arg\min_W \|Y_{\text{target}} - WR\|_F^2 + \beta\|W\|_F^2$$ $$W_{\text{out}} = Y_{\text{target}} R^\top (R R^\top + \beta I)^{-1}$$

Tikhonov regularization. \(\beta\) is typically in \([10^{-8}, 10^{-2}]\). You need it — \(R R^\top\) is \(N \times N\) and usually rank-deficient because reservoir states are correlated across time. Without \(\beta\), the pseudoinverse amplifies noise in degenerate directions.

The solve is a Cholesky factorization of the \(N \times N\) matrix \(R R^\top + \beta I\). That's \(O(N^3)\). For \(N = 300\), \(T = 5000\), total training time is sub-millisecond on a single core. Not gradient descent. One linear solve.

After training, prediction is autoregressive: feed the model's own output back as input. For chaotic systems the quality decays — errors grow exponentially at the Lyapunov rate — but for short-to-medium horizons an ESN is competitive with models that cost orders of magnitude more to train.

Memory capacity

There's a formal metric for how far back a reservoir can remember. Drive it with i.i.d. input \(u(t) \sim \mathcal{U}(-1, 1)\). For a given delay \(k\), train a separate readout to reconstruct \(u(t-k)\) from the current state \(\mathbf{r}(t)\):

$$\hat{y}_k(t) = W_{\text{out},k} \cdot \mathbf{r}(t), \quad \text{MC}_k = 1 - \frac{\text{MSE}(\hat{y}_k, u(t-k))}{\text{Var}(u)} \in [0,1]$$

Total memory capacity \(\text{MC} = \sum_{k=1}^{\infty} \text{MC}_k \le N\). That's a hard bound, proved by Jaeger. Typical numbers: \(\rho = 1, \alpha = 0.5, N = 100\) gives \(\text{MC} \approx 10{-}15\).

There's a tradeoff. Higher spectral radius → more memory. But higher spectral radius also → more tanh saturation → less capacity for nonlinear separation. Smaller \(\alpha\) shifts the memory profile toward longer delays. Every model variant below is basically tinkering with this tradeoff in a different way.

The package: how it's built

ReservoirComputing.jl lives inside the SciML ecosystem — the same ecosystem that gives you 200+ ODE solvers, physics-informed neural networks, and adjoint sensitivity methods. The package itself uses Lux.jl, a functional deep learning framework where layers are parameter containers and the forward pass is a pure function.

Julia's type system does heavy lifting here. Models are concrete structs with compile-time-known field tuples:

abstract type AbstractReservoirComputer <: Lux.AbstractLuxLayer end

@concrete struct ESN <: AbstractEchoStateNetwork{(
    :reservoir, :states_modifiers, :readout)}
    reservoir::ESNCell
    states_modifiers
    readout::LinearReadout
end

The @concrete macro auto-generates constructors and property accessors. The field tuple (:reservoir, :states_modifiers, :readout) is a compile-time constant — the training pipeline unpacks models without reflection. The cell type holds the actual weights:

@concrete struct ESNCell <: AbstractEchoStateNetworkCell
    activation
    in_dims  ::Int
    out_dims ::Int
    init_bias
    init_reservoir
    init_input
    init_state
    use_bias ::StaticBool       # True() or False() at type level
end

function (cell::ESNCell)(x, ps, st)
    new_r = tanh.(ps.input_layer(x) .+ ps.reservoir_layer(st.r) .+ ps.bias)
    st_new = (r = (1 - cell.leak_rate) .* st.r .+ cell.leak_rate .* new_r,)
    return st_new.r, st_new
end

StaticBool is one of those Julia idioms that feels like cheating. True() and False() are types, not values. If use_bias = False(), the compiler eliminates the bias branch entirely at compile time. No runtime check. No branch prediction miss. The bias addition literally doesn't exist in the generated machine code.

The broadcast operations (.+ and tanh.) are fused into a single loop with zero intermediate allocations. On a Sapphire Rapids Xeon with AVX-512, this runs at ~2.5 TOPS/core on Float32.

The training pipeline has three stages: collectstates (drive reservoir, gather \(R\)), train! (solve for \(W_{\text{out}}\)), predict (autoregressive rollout or teacher-forced). Any type that satisfies AbstractReservoirComputer inherits all three. Add a new cell type, implement the forward pass, and everything downstream just works.

SVESM — SVM readout

Shi & Han, IEEE TNN 2006. The simplest change you can make: swap the linear readout for an SVM.

Linear readouts work well for regression but classification sometimes benefits from a maximum-margin decision boundary. The reservoir states for class \(k\) form a set \(R_k = \{\mathbf{r}(t) : y(t) = k\} \subset \mathbb{R}^N\). An SVM finds the hyperplane that maximizes distance to the nearest point in each class.

Primal:

$$\min_{\mathbf{w}, b, \xi} \frac{1}{2}\|\mathbf{w}\|^2 + C\sum_i \xi_i \quad \text{s.t.} \quad y_i(\mathbf{w}^\top \mathbf{r}_i + b) \ge 1 - \xi_i,\; \xi_i \ge 0$$

Dual, solved via SMO (Platt 1998):

$$\max_{\alpha} \sum_i \alpha_i - \frac{1}{2}\sum_{i,j} \alpha_i\alpha_j y_i y_j K(\mathbf{r}_i, \mathbf{r}_j) \quad \text{s.t.} \quad 0 \le \alpha_i \le C,\; \sum_i \alpha_i y_i = 0$$

The reservoir already projects inputs into nonlinear features, so a linear kernel \(K(\mathbf{r}_i, \mathbf{r}_j) = \mathbf{r}_i^\top \mathbf{r}_j\) is usually sufficient. That was the original paper's insight. RBF and polynomial kernels are exposed as kwargs but default to linear. LIBSVM.jl does the actual solve; the model wrapper is just SVESM(in_dims, res_dims, out_dims) with an ESNCell + SVMReadout. Training calls svmtrain instead of ridge regression. The reservoir doesn't know the difference.

LIFESN — information-flow initialization

Nonlinear Dynamics 2025. Standard ESN initialization is lazy: every neuron connects to ~5% of other neurons with random weights drawn from a uniform distribution. This treats all neurons as interchangeable. LIFESN says: measure the actual information flow between neurons, then wire the ones that talk to each other.

Transfer entropy from process \(X\) to \(Y\):

$$T_{X \to Y} = \sum p(y_{t+1}, y_t^{(k)}, x_t^{(l)}) \cdot \log_2 \frac{p(y_{t+1} \mid y_t^{(k)}, x_t^{(l)})}{p(y_{t+1} \mid y_t^{(k)})}$$

\(y_t^{(k)}\) is the \(k\)-length history of \(Y\). \(T_{X \to Y}\) measures how much knowing \(X\)'s past reduces uncertainty about \(Y\)'s future, beyond what \(Y\)'s own past already tells you. It's directed (\(T_{X \to Y} \neq T_{Y \to X}\)) and zero iff \(X\) and \(Y\) are conditionally independent.

The algorithm:

1. Drive a temporary ESN with white noise for a few hundred steps. Collect states.
2. For every ordered neuron pair \((i, j)\), estimate \(T_{r_i \to r_j}\) using kernel density estimation on time-delayed state triplets.
3. Form a directed adjacency matrix: \(A_{ij} = 1\) if \(T_{r_i \to r_j} > \theta\), else 0. \(\theta\) is set at the 95th percentile of the null distribution from surrogate data — randomly shuffled time indices that preserve marginals while breaking actual information flow. Standard Schreiber surrogate test.
4. For connected pairs, set \((W_r)_{ij} \propto T_{r_i \to r_j}\). Disconnected pairs get zero.
5. Rescale to target spectral radius.

The TE estimation costs \(O(N^2 \cdot T_{\text{cal}} \cdot k^2)\). For \(N = 300, T_{\text{cal}} = 1000, k = 5\): about 225 million operations, 2–5 seconds on one core. This only runs once at initialization. Inference cost is identical to a standard ESN — the weight matrix is just structured differently.

LIFESNCell adds calibration parameters to the cell struct. The init_reservoir field is a closure that runs the TE estimation during Lux.setup. After setup, the cell is indistinguishable from a standard ESNCell. Same update equation, same training pipeline. Only the connectivity pattern is different.

ResESN — decoupled timescales

Neurocomputing 2024. The standard ESN has one dial for temporal dynamics: the leak rate \(\alpha\). It controls both how much old state you keep (\(1-\alpha\)) and how much new activation you admit (\(\alpha\)). These are antagonistic. If you want long memory, you need small \(\alpha\). But small \(\alpha\) also makes the reservoir slow to respond to new inputs. You can't have both.

ResESN splits \(\alpha\) into two independent parameters:

$$\mathbf{r}(t) = (1 - \alpha)\mathbf{r}(t-1) + \beta \tanh(W_{\text{in}}\mathbf{u}(t) + W_r\mathbf{r}(t-1) + \mathbf{b})$$

\(\alpha\): state retention (skip connection). \(\beta\): input sensitivity (transform scale). The ESP condition becomes \(\|W_r\|_2 \cdot \beta < 1\).

In the frequency domain, the reservoir state is a low-pass filter. The cutoff frequency is \(f_c = \alpha/(2\pi(1-\alpha))\) (discrete-time units). \(\beta\) controls the passband gain. You can set \(\alpha = 0.1\) (memory ~55 steps back, \(f_c \approx 0.018\)) and \(\beta = 1.5\) (strong response to new inputs) simultaneously. That's impossible in the standard formulation where \(\alpha = \beta\).

Linearizing around the origin (where \(\tanh(\mathbf{x}) \approx \mathbf{x}\)):

$$\mathbf{r}(t) \approx ((1 - \alpha)I + \beta W_r)\mathbf{r}(t-1) + \beta W_{\text{in}}\mathbf{u}(t) + \beta\mathbf{b}$$

The effective recurrent matrix is \((1-\alpha)I + \beta W_r\) with eigenvalues \((1-\alpha) + \beta\lambda_i\). ESP requires \(|(1-\alpha) + \beta\lambda_i| < 1\) for all \(i\). This is tighter than \(\rho(W_r) < 1\) — large \(\beta\) can push sub-spectral-radius eigenvalues past the stability boundary. The constructor checks this at init.

ResESNCell stores alpha::Float64 and beta::Float64 instead of a single leak_rate. Default \(\alpha = \beta = 0.7\) recovers the standard ESN exactly. Tests verify identical output to floating-point precision at those defaults, and confirm that \(\alpha = 0.1, \beta = 1.5\) produces different output. The decoupling is real.

RMNResESN — gated memory + residual skip

ESANN 2024. The most architecturally involved of the four. Combines two ideas: a Reservoir Memory Network (explicit memory state with per-dimension gating) wrapped in a residual ESN (decoupled skip connections).

A standard ESN's memory is distributed — every neuron carries a decaying trace of past inputs, all with the same decay rate \(\alpha\). There's no way to say "remember this specific thing for a long time but forget that other thing quickly." The RMN adds a dedicated memory vector \(\mathbf{m}(t) \in \mathbb{R}^M\) (typically \(M = N\)) with per-dimension forget gates:

$$\begin{aligned} \mathbf{g}(t) &= \sigma(W_g \mathbf{r}(t-1) + \mathbf{b}_g) \quad &&\mathbf{g}(t) \in (0,1)^M \\ \mathbf{m}(t) &= \mathbf{g}(t) \odot \mathbf{m}(t-1) + (1 - \mathbf{g}(t)) \odot \mathbf{r}(t-1) \end{aligned}$$

Hadamard product. \(\sigma\) is logistic sigmoid. Each dimension \(j\) independently decides its retention rate \(g_j(t)\) based on the current reservoir state. \(g_j \to 1\): preserve old memory. \(g_j \to 0\): overwrite with current reservoir state.

This is a stripped-down LSTM: forget gate only. No input gate — the "new content" is just the reservoir state. No output gate — the memory is exposed directly. The gate parameters \(W_g\) and \(\mathbf{b}_g\) are trained via logistic regression, not backprop. The training signal comes from the mismatch between the current reservoir state and the delayed reconstruction target.

Combined with ResESN, the full system:

$$\begin{aligned} \tilde{\mathbf{r}}(t) &= \tanh(W_{\text{in}}\mathbf{u}(t) + W_r\mathbf{r}(t-1) + W_m\mathbf{m}(t-1) + \mathbf{b}) \\ \mathbf{r}(t) &= (1 - \alpha)\mathbf{r}(t-1) + \beta\tilde{\mathbf{r}}(t) \\ \mathbf{g}(t) &= \sigma(W_g\mathbf{r}(t-1) + \mathbf{b}_g) \\ \mathbf{m}(t) &= \mathbf{g}(t) \odot \mathbf{m}(t-1) + (1 - \mathbf{g}(t)) \odot \mathbf{r}(t-1) \\ \hat{\mathbf{y}}(t) &= W_{\text{out}}[\mathbf{r}(t);\; \mathbf{m}(t)] \end{aligned}$$

The readout sees the concatenated state \([\mathbf{r}(t); \mathbf{m}(t)] \in \mathbb{R}^{2N}\). Reservoir = short-timescale nonlinear dynamics. Memory = long-timescale linear trends, per-dimension adaptive timescales. The \(W_m \mathbf{m}(t-1)\) term means memory can modulate the reservoir's nonlinear processing. It's functionally top-down attention even though the paper doesn't use that vocabulary.

The combined state-space:

$$\begin{bmatrix} \mathbf{r}(t) \\ \mathbf{m}(t) \end{bmatrix} = \begin{bmatrix} (1-\alpha)I + \beta J_r(t)W_r & \beta J_r(t)W_m \\ I - \operatorname{diag}(\mathbf{g}(t)) & \operatorname{diag}(\mathbf{g}(t)) \end{bmatrix} \begin{bmatrix} \mathbf{r}(t-1) \\ \mathbf{m}(t-1) \end{bmatrix} + \begin{bmatrix} \beta J_r(t)(W_{\text{in}}\mathbf{u}(t) + \mathbf{b}) \\ \mathbf{0} \end{bmatrix}$$

\(J_r(t) = \operatorname{diag}(\tanh'(\cdots))\). Near the origin (\(J_r \approx I\)), the recurrent matrix is block-triangular. The ESP depends on \(\alpha, \beta, \rho(W_r)\), and gate saturation. Spectral radius heuristic holds empirically; the formal proof is messier because the gate makes the linearization state-dependent.

Implementation: four types, two composable cells.

MemoryResESNCell  <: AbstractEchoStateNetworkCell  # memory-augmented reservoir update
RMNCell           <: AbstractReservoirComputer      # standalone gating, wraps any reservoir cell
RMNESN            <: AbstractReservoirComputer      # RMNCell + ESNCell
RMNResESN         <: AbstractReservoirComputer      # RMNCell + ResESNCell

The key design decision was making RMNCell a standalone wrapper rather than baking memory into the reservoir cell. This means you can later do RMN + DeepESN or RMN + EuSN without rewriting the gating logic. Composability over convenience.

Test coverage: \(\mathbf{g}(t) \in (0,1)\) bounds check, \(\mathbf{m}(t)\) shape contracts, determinism under fixed seed, autoregressive rollout on a synthetic delayed-XOR task (output = XOR of inputs at \(t\) and \(t-50\), explicitly impossible without explicit memory), and equivalence collapse to ResESN when \(W_m = 0\) and the gate is frozen at \(\mathbf{g}(t) = \mathbf{1}\).

How they relate

ESN
├─ SVESM           readout: ridge → SVM
├─ LIFESN          init: random → TE-driven
├─ ResESN          dynamics: α=β → α,β decoupled
└─ RMNResESN       architecture: + memory, + gate, inherits ResESN
   └─ RMNESN       memory on vanilla ESN, subsumed at α=β

The modifications are orthogonal axes: readout, initialization, temporal dynamics, architectural depth. Any combination works — LIFESN + SVESN, ResESN + SVESN, RMN + LIFESN — because each model touches a different component of the pipeline. No new code needed, just constructor chaining. The package's test infrastructure is designed for this: any type implementing AbstractReservoirComputer automatically gets state collection, training, prediction, and basic metrics.

What it actually costs

Per-step inference:

$$\begin{aligned} W_r\mathbf{r}(t) &: O(N^2),\; \text{sparse at density } \rho_s \\ W_{\text{in}}\mathbf{u}(t) &: O(ND) \\ \tanh + \text{leak} &: O(N) \\ \text{Total for } N = 300, D = 10, \rho_s = 5\% &: \sim\!5 \times 10^4 \text{ FLOPs/step} \end{aligned}$$

Sapphire Rapids, 2.5 GHz, AVX-512 (32 FP32 ops/cycle): ~0.04 μs/step theoretical. In practice with Julia + LoopVectorization: ~0.1–0.3 μs/step. 10,000-step prediction: ~2–3 ms.

Training:

$$\begin{aligned} \text{State collection} &: T \cdot O(N^2) \;\approx\; 3 \times 10^7 \text{ FLOPs} \quad (N = 300, T = 5000) \\ \text{Readout solve} &: O(N^3) \;\approx\; 9 \times 10^6 \text{ FLOPs} \\ \text{Total} &: \sim\!4 \times 10^7 \text{ FLOPs} \;\approx\; 15{-}20\;\mu\text{s} \end{aligned}$$

Memory: \(W_r\) (sparse) ~100 KB, \(W_{\text{in}}\) ~10 KB, \(W_{\text{out}}\) ~10 KB. Total under 200 KB for \(N = 300\), entirely in L2 cache. \(N = 10{,}000\) at 1% sparsity: ~800 KB, fits in L3. This thing runs on microcontrollers. The bottleneck is the readout Cholesky at \(N \approx 5000\). Past that you switch to conjugate gradient via LinearSolve.jl.

Where this fits in SciML

SciML is a differential equation ecosystem that happens to contain ML. OrdinaryDiffEq.jl has 200+ solvers. NeuralPDE embeds neural networks inside PDEs. DiffEqFlux does adjoint sensitivity training. ReservoirComputing.jl provides the recurrent models that sit at the junction of data-driven and equation-driven modeling.

Concrete cross-package examples:

• ESN as ODE right-hand side: \(\dot{\mathbf{x}} = f_{\text{ESN}}(\mathbf{x}, \mathbf{u}, t)\). Train readout on observations, solve with Tsit5 for interpolation.
• Universal differential equations: known physics + ESN learned residual. Joint training via DiffEqFlux.
• LIFESN transfer entropy matrix as structural prior for a Neural ODE's connectivity graph.

The SciMLProblemReservoir type formalizes this. Any ODEProblem becomes a reservoir. Any reservoir becomes a SciML component. The interface is the abstraction.

References

Jaeger, H. (2001). The "echo state" approach to analysing and training recurrent neural networks. GMD Report 148.
Lukoševičius, M. (2012). A practical guide to applying echo state networks. Neural Networks: Tricks of the Trade, 2nd ed.
Martinuzzi, F., Rackauckas, C., et al. (2022). ReservoirComputing.jl: An efficient and modular library for reservoir computing models. JMLR 23(288).
Shi, Z. & Han, M. (2006). Support vector echo-state machine for chaotic time-series prediction. IEEE TNN 17(3).
[LIFESN] (2025). Local information flow echo state network. Nonlinear Dynamics.
[ResESN] (2024). Residual echo state networks with decoupled skip connections. Neurocomputing.
[RMNResESN] (2024). Reservoir memory networks with residual connections. ESANN.