🧠 memory-esn: Long-Memory Echo State Networks

Reservoir Computing for Time-Series Forecasting with Long-Range Dependence 📈✨

Long-Memory Echo State Networks for Time-Series Forecasting

A reservoir-computing library that augments the classic Echo State Network with a dedicated memory reservoir for long-range dependence. Two memory mechanisms are provided — a fractional-difference filter (fESN) and a wavelet-smoothed fractional filter (wESN) — while keeping all internal weights fixed after random initialization and training only a linear ridge readout. BLAS-backed Cython kernels make it fast; pure-NumPy fallbacks keep it portable.

The vanilla reservoir captures short-term nonlinear dynamics (exponential forgetting); the memory reservoir sustains long-term dynamics (polynomially decaying weights). Their states are concatenated and mapped to the output.

Class hierarchy

PersistenceMixin              save() / load()  (joblib)
│
BaseESN                       single reservoir + ridge readout
│
MultiESN                      N reservoirs, one input each (composes BaseESN)
└── DoubleReservoirESN        fixed to 2 reservoirs;  fit([X1, X2], y)
    ├── fESN                  memory reservoir <- ((1-B)^d - 1) u(t)
    └── wESN                  memory reservoir <- ((1-B)^d - 1) MODWT_smooth(u(t))

Reservoir 1 (vanilla) always sees the raw series; reservoir 2 (memory) sees a pure-past, long-memory filter of it. fESN/wESN take a single series and build the memory input internally; DoubleReservoirESN and MultiESN take the inputs explicitly.

Plus TimeSeriesDataset — sliding-window train/val/test splitting with leakage-free scaling (none / minmax / standard / log).

Defaults follow the paper

Out of the box the models match the paper's construction: uniform weight initialization, reservoirs sparsified to 90% zeros then rescaled to a spectral radius < 1, and a memory filter f(t) = ((1-B)^d - 1) u(t) = Σ_{k≥1} ω_k u(t-k) (standard fractional differencing with the present term removed, so the memory reservoir is driven purely by the past). Optional isotropic state noise η(t) ~ N(0, σ² I) is available via noise=σ.

Install

pip install -e .                       # installs deps; builds Cython kernels if a
                                       # C compiler is available
# or build the fast kernels in place:
python setup.py build_ext --inplace

The compiled kernels are optional — pure-NumPy fallbacks are used automatically when they are not present (with a one-time warning). Check which is active:

import memory_esn as m
m.USING_CYTHON_RESERVOIR, m.USING_CYTHON_FRACDIFF, m.USING_CYTHON_MODWT

Dependencies: numpy, scipy, scikit-learn, joblib, PyWavelets.

Quickstart

import numpy as np
from memory_esn import BaseESN, fESN

X = np.cumsum(np.random.randn(1000, 3), axis=0) * 0.1
y = np.sin(X[:, :1])
Xtr, Xte, ytr, yte = X[:800], X[800:], y[:800], y[800:]

# Single reservoir
esn = BaseESN(n_reservoir=200, spectral_radius=0.95, random_state=0)
esn.fit(Xtr, ytr, washout=100)
y_pred = esn.predict(Xte)

# Fractional ESN (raw + fractionally-differenced branch)
fesn = fESN(n_reservoir=(200, 150), d=0.5, K=100, random_state=0)
fesn.fit(Xtr, ytr, washout=100)
y_pred = fesn.predict(Xte, continuation=True)   # continuation keeps output length

See examples/quickstart.py for one example per class.

Wavelet options (wESN)

wESN decomposes each input feature with a MODWT (undecimated, shift-invariant) and feeds the selected component(s) to reservoir 2. The wavelet is validated at construction:

• Discrete wavelets (usable): haar, db1..db38, sym2..sym20, coif1..coif17, bior*, rbio*, dmey — 106 in total (pywt.wavelist(kind='discrete')).
• Continuous wavelets (mexh, morl, gaus*, cmor, …) raise NotImplementedError — support is planned for a future release.

Pick which decomposition components form reservoir 2's input via wavelet_components:

from memory_esn import wESN

# default: level-J smooth (approximation) trend — one component
wESN(wavelet='db4', wavelet_level=2)                       # components = [A2]

# a single detail level
wESN(wavelet='db4', wavelet_level=3, wavelet_components=2) # components = [D2]

# MULTIVARIATE: several levels stacked -> reservoir 2 gets multiple channels
wESN(wavelet='db4', wavelet_level=3,
           wavelet_components=[1, 2, 'smooth'])                  # components = [D1, D2, A3]

# shortcuts
wESN(wavelet='haar', wavelet_level=2, wavelet_components='all')      # D1, D2, A2
wESN(wavelet='haar', wavelet_level=2, wavelet_components='details')  # D1, D2

wavelet_level is the decomposition depth J; detail levels are 1..J and 'smooth' is the level-J approximation. Selecting n components on an F-feature series gives reservoir 2 an F × n-channel (multivariate) input.

MODWT normalization — wavelet_norm='modwt' (default) uses the standard Percival–Walden rescaling (Haar smooth filter [1/2, 1/2]); wavelet_norm='classic' uses the DWT orthonormal filters ([1/√2, 1/√2]). The two differ only by a per-level constant absorbed by the random memory-input weights, so they are model-equivalent.

Multivariate fESN (multiple differencing orders)

fESN has the analogous knob: d accepts a single order (default, univariate) or a list of orders that are stacked as channels.

from memory_esn import fESN

fESN(d=0.5)                 # univariate (default)
fESN(d=[0.3, 0.5, 0.8])     # multivariate: 3 orders stacked

Selecting n orders on an F-feature series gives reservoir 2 an F × n-channel input. Defaults for both classes are univariate: fESN d=0.5, wESN wavelet_components='smooth'.

wESN also accepts a list d, which cross-products with the wavelet components: n_features × n_components × n_d channels. For example wESN(wavelet_components=[1, 2, 'smooth'], d=[0.3, 0.5]) on a 2-feature series → 2 × 3 × 2 = 12 channels into reservoir 2.

Weight initialization & reservoir options

Every reservoir weight is drawn from a configurable distribution (memory_esn.weights):

parameter	controls	options	default
input_init	input weights W_in	uniform, gaussian, bernoulli, laplace	uniform
reservoir_init	reservoir weights W	same	uniform
bias_init	neuron bias	same	uniform
input_scaling	scale of W_in (= U[-ψ, ψ] for uniform)	float	0.5
spectral_radius	W rescaled to this ρ	float <1	0.9
sparsity	fraction of W set to zero (φ)	[0,1)	0.9
bias_scaling	scale of bias (0 = none)	float	0.0
noise	std σ of state noise η(t)~N(0,σ²I)	float ≥ 0	0.0

In MultiESN each is broadcastable — a scalar applies to all reservoirs, a list of length n_reservoirs sets them per reservoir. In the two-reservoir models a scalar or a (vanilla, memory) pair is accepted. random_state follows the same rule: a single int derives distinct per-reservoir seeds (reproducible but not identical); a list sets them explicitly.

Paper ↔ code notation

paper	meaning	code
p, q	vanilla / memory reservoir sizes	n_reservoir=(p, q)
ρ_x, ρ_m	spectral radii	spectral_radius=(ρ_x, ρ_m)
φ_x, φ_m	sparsification proportion	sparsity=(φ_x, φ_m)
ψ_x, ψ_m	input scalings	input_scaling=(ψ_x, ψ_m)
σ_x, σ_m	state-noise std	noise=(σ_x, σ_m)
d	fractional-differencing order	d
K	filter truncation (lags)	K
ζ	washout ratio T₀/T	washout=int(ζ*T)
H	forecast horizon	number of columns of y

The per-horizon readouts W_out^(h) are fitted jointly as a multi-output ridge on a single design matrix Φ = [x; m; u; f] ∈ ℝ^{p+q+2}, with RidgeCV selecting λ by efficient leave-one-out — matching the paper's shared-Cholesky LOOCV.

Notes on alignment (fESN / wESN)

The fractional-difference filter needs K past samples. With skip_initial=True (default) and continuation=False, predictions are K steps shorter than the input. With continuation=True, the model prepends stored history so the output matches the input length — this is the normal mode for streaming / iterative forecasting. wESN.predict uses continuation=True by default (it also re-smooths using training history to avoid wavelet boundary effects).

Speed-up kernels

memory_esn/_speedups/ holds three Cython extensions (with NumPy fallbacks):

kernel	what it does
esn_reservoir_optimized	reservoir state update via BLAS dgemv, nogil, 12 activations
fractional_diff	binomial fractional-difference weights + windowed filtering
modwt_fast	maximal-overlap DWT (smoothing) via circular convolution

Tests

pip install pytest
python -m pytest tests/          # run from the project root

The suite includes Cython-vs-NumPy numerical-equivalence checks, so the fast and fallback paths are guaranteed to agree.

Repository layout

memory_esn/            the package (import as `memory_esn`)
  base.py              BaseESN, PersistenceMixin
  multi.py             MultiESN
  double.py            DoubleReservoirESN
  fractional.py        fESN  (alias FractionalESN)
  wavelet.py           wESN  (alias WaveletESN)
  dataset.py           TimeSeriesDataset
  _speedups/           Cython kernels + NumPy fallbacks
examples/quickstart.py
tests/test_package.py
setup.py               builds the Cython extensions
garbage/               the previous, pre-refactor files (archived)

Authors

• Rahul Goswami | Shinjini Paul | Palash Ghosh | Tanujit Chakraborty

Citation

If you use memory-esn in your research, please cite it (see CITATION.cff).

License

MIT — see LICENSE.