🧠 Reservoir Neural Networks with Fourier Layer

By ACS Lab, ITMO University · 2025

Official repository for the paper:

A. Kovantsev, R. Vysotskiy (2025). Advanced Reservoir Neural Network Techniques for Chaotic Time Series Prediction. SSRN 5481760.

We present ESN‑F — the Echo State Network (ESN) enhanced with Fourier features and polynomial expansion for forecasting chaotic / nonlinear time series. The approach keeps the reservoir untrained and learns only a ridge readout, while enriching inputs with periodic (sin/cos) and nonlinear bases that improve long‑horizon stability.

Use cases include finance/economics, risk modeling, and other non‑stationary domains.

• Why this project?

• Install

• Quickstart

• Generative (multi-step) forecasting

• Multivariate example

• Background & Theory

  • Echo State Network (leaky ESN)
  • Fourier & Polynomial feature blocks (FAN)
  • Ridge readout & objective
  • Forecasting strategies
  • Why it helps with chaos & non‑stationarity

• Predictability metrics (optional)

• 📐 Predictability features (formulas)

• Hyperparameters

• Complexity & scaling

• Reproducibility

• 📊 Experimental Results (from the paper)

• Project layout

• Contributing

• License & citation

• Troubleshooting / FAQ

Long‑horizon forecasting on chaotic / weakly stationary signals is tricky: standard RNNs tend to drift, while purely statistical baselines miss nonlinear structure. EnhancedESN_FAN keeps a random, untrained reservoir for rich dynamics and augments the readout with deterministic Fourier harmonics and polynomial features. The linear readout is trained via ridge regression, keeping training fast, convex, and robust.

git clone https://github.com/CapitalistGeorge/chaotic_library.git
cd chaotic_library

python -m venv .venv

# macOS/Linux
source .venv/bin/activate

# Windows (PowerShell)
# .venv\Scripts\Activate.ps1

python -m pip install -U pip
pip install -e .

Planned: after the API is frozen we’ll publish to PyPI, so you can pip install esn-fan (package name TBD). The import path in the examples below assumes a module enhanced_esn_fan.py at the project root; adjust if packaged differently.

import numpy as np
from chaotic_library import EnhancedESN_FAN  # adjust import path if packaged differently

# 1) Build training data (shape: [n_timesteps, input_dim])
# Univariate example: input_dim=1 → X is 2D with one column, y is 1D
T = 1200
noise = 0.1 * np.random.randn(T)
signal = np.sin(2*np.pi*np.arange(T)/50) + 0.25*np.sin(2*np.pi*np.arange(T)/7) + noise

X = signal[:-1].reshape(-1, 1)   # features are previous value(s)
y = signal[1:]                   # next-step target

# 2) Initialize Enhanced ESN + FAN features (Fourier + polynomial)
esn = EnhancedESN_FAN(
    input_dim=1,          # number of input features per timestep (columns of X)
    reservoir_size=800,
    spectral_radius=0.95,
    sparsity=0.1,
    ridge_alpha=1e-2,
    leaking_rate=0.3,
    poly_order=2,
    fan_terms=8,
    random_state=42,
    clip_value=3.0,
)

# 3) Fit and one-step-ahead predictions (teacher forcing / open loop)
esn.fit(X, y)
y_hat = esn.predict(X[:100])        # shape: (100,) for univariate target
print("y_hat shape:", np.asarray(y_hat).shape)

# Seed with the last observed input row (shape: [1, input_dim])
seed = X[-1:].copy()

# Produce next 300 steps autoregressively
future = esn.predict(seed, generative_steps=300)   # shape: (300,) for univariate

# Concatenate history + forecast for plotting
full = np.concatenate([signal, future.ravel()])

import numpy as np
from chaotic_library import EnhancedESN_FAN

# Suppose you have 3 exogenous drivers + the main signal → input_dim=4
n = 2000
main = np.sin(2*np.pi*np.arange(n)/30) + 0.05*np.random.randn(n)
x1 = np.cos(2*np.pi*np.arange(n)/100)
x2 = np.sin(2*np.pi*np.arange(n)/7)
x3 = 0.01*np.arange(n)  # slow trend proxy

# Features: use current drivers + lagged main as input; predict next main
X = np.column_stack([main[:-1], x1[:-1], x2[:-1], x3[:-1]])  # shape (n-1, 4)
y = main[1:]                                               # shape (n-1,)

model = EnhancedESN_FAN(
    input_dim=4,
    reservoir_size=500,
    spectral_radius=0.9,
    sparsity=0.1,
    ridge_alpha=0.1,
    leaking_rate=0.4,
    poly_order=2,
    fan_terms=6,
    random_state=7,
)

model.fit(X, y)
pred = model.predict(X[:128])              # one-step predictions (teacher forcing)
gen  = model.predict(X[-1:], generative_steps=200)  # recursive forecast

Reservoir state $\mathbf{x}_t \in \mathbb{R}^{N_r}$ evolves under a fixed random dynamical system:

$$ \begin{aligned} \tilde{\mathbf{x}}_t &= \tanh\left( \mathbf{W} \mathbf{x}_{t-1} + \mathbf{W}_{\text{in}} [1; u_t] \right), \\ \mathbf{x}_t &= (1-\alpha) \mathbf{x}_{t-1} + \alpha \tilde{\mathbf{x}}_t, \end{aligned} $$

where:

• $u_t$ is the input (e.g., components of $X_t$)
• $[1; u_t]$ denotes a bias‑augmented input
• $\alpha$ is the leaking rate (leaking_rate parameter)

To satisfy the echo state property (state forgets initial conditions), scale the reservoir so that its spectral radius $\rho(\mathbf{W})$ is near 1 (practically 0.7–1.2 with leakage; spectral_radius parameter).

We collect features at time $t$ by concatenating the reservoir state with deterministic blocks:

$$ \mathbf{z}_t = \big[ \mathbf{x}_t \mid \phi_{\text{poly}}(u_t) \mid \phi_{\text{Fourier}}(u_t) \big] \in \mathbb{R}^{D} $$

Notation:

• $\mathbf{W}$: reservoir weight matrix
• $\mathbf{W}_{\text{in}}$: input weight matrix
• $\phi_{\text{poly}}$: polynomial features
• $\phi_{\text{Fourier}}$: Fourier features
• $D$: total feature dimension

• PolynomialFeatures of degree $d$ (poly_order parameter): $[u_t, u_t^2, \dots, u_t^d]$ per input dimension (no extra bias term; bias provided separately).
• Fourier (FAN) features with harmonics $k=1..K$ (fan_terms parameter): for each input dimension, compute $\sin(2\pi k X)$ and $\cos(2\pi k X)$. These inject periodic structure explicitly, so the reservoir does not have to "discover" it from scratch.

Only the final linear readout $\mathbf{W}_{\text{out}} \in \mathbb{R}^{D\times m}$ is trained via ridge regression:

$$ \min_{\mathbf{W}_{\text{out}}} \left|\mathbf{Y} - \mathbf{Z}\mathbf{W}_{\text{out}}\right|_2^2 + \lambda \left|\mathbf{W}_{\text{out}}\right|_2^2 $$

$$ \Rightarrow\quad \mathbf{W}_{\text{out}} = (\mathbf{Z}^\top\mathbf{Z} + \lambda \mathbf{I})^{-1}\mathbf{Z}^\top\mathbf{Y}. $$

Columns of $\mathbf{Z}$ should be standardized for numerical stability (the implementation uses StandardScaler).

• Teacher forcing / open loop (default in predict(X)): one‑step predictions using the provided inputs.
• Generative / recursive (predict(seed, generative_steps=m)): feed back model outputs as inputs to generate future steps.
• Hybrid (future option): recursive core with direct corrections for selected horizons.

• Reservoir provides a rich, fading memory of nonlinear histories.
• Fourier layer anchors periodic structure → less burden on the reservoir.
• Polynomial bias stabilizes local trends and offsets.
• Ridge readout tames variance and keeps training convex & fast.

You can compute these to cluster series by predictability and adapt hyperparameters:

• Hurst exponent (H) — persistence (>0.5) vs anti‑persistence (<0.5) vs (=0.5) random walk
• Correlation dimension (D₂) — attractor dimension (Grassberger–Procaccia)
• Max Lyapunov exponent (λₘₐₓ) — sensitivity to initial conditions
• Kolmogorov–Sinai entropy (KSE) — information production rate
• # Prevailing harmonics — count strong spectral peaks (e.g., via periodogram)

Use the cluster to pick reservoir_size, spectral_radius, and fan_terms. For highly chaotic signals (large λₘₐₓ), prefer slightly lower spectral_radius and stronger regularization (ridge_alpha).

NOTATION:

• $\bar{x}_\tau$ - sample mean on window of length $\tau$
• $\theta(\cdot)$ - Heaviside step function
• $\rho(i,j)$ - distance in reconstructed phase space (delay embedding optional)
• $x_i' = x_i - x_{i-1}$ (first difference)

$$ H = \frac{\ln\left( R(\tau) / S(\tau) \right)}{\ln(\alpha \cdot \tau)} $$

where:

$$ \begin{aligned} R(\tau) &amp;= \max_{1 \leq t \leq \tau} \left[ \sum_{i=1}^{t} (x_i - \bar{x}_\tau) \right] - \min_{1 \leq t \leq \tau} \left[ \sum_{i=1}^{t} (x_i - \bar{x}_\tau) \right] \\ S(\tau) &amp;= \sqrt{ \frac{1}{\tau} \sum_{t=1}^{\tau} (x_t - \bar{x}_\tau)^2 } \end{aligned} $$

Note: In classical R/S analysis, $H$ is the slope of $\ln(R/S)$ vs $\ln\tau$ (i.e., $\alpha=1$). Including a constant $\alpha$ is equivalent up to offset.

Definition via entropy-rate upper bound:

$$ \begin{aligned} h_\mu(T,\xi) &= - \lim_{n\to\infty} \frac{1}{n} \times \sum_{i_1,\dots,i_n} \mu( T^{-1}C_{i_1} \cap \dots \cap T^{-n}C_{i_n} ) \cdot \ln \mu(\dots) \\ h_\mu^{KS}(T) &= \sup_{\xi} h_\mu(T,\xi) \end{aligned} $$

$$ d_2 = \lim_{r\to 0} \lim_{m\to\infty} \frac{\ln C(r)}{\ln r} $$

where:

$$ C(r) = \frac{1}{m(m-1)} \times \sum_{i=1}^{m} \sum_{j=i+1}^{m} \theta\left( r - \rho(i,j) \right) $$

• State update: $O(T \cdot N_r \cdot s)$ with sparsity fraction $s$ (dense → $O(T \cdot N_r^2)$)
• Readout training: build $\mathbf{Z} \in \mathbb{R}^{T \times D}$; solve ridge via Cholesky/QR: $\sim O(D^3)$ (usually $D \ll T$)
• Memory: $O(T \cdot D)$ if keeping all features; use chunked/online solvers for very long series

• Fix random_state for weights and reservoirs.
• Standardize inputs and feature matrix consistently across train/forecast.
• Log: hyperparameters, seeds, and package versions.
• Provide notebooks that mirror experiments and regenerate figures.

In the Bad cluster, ESN‑F beats LGBM by ≥1 pp in 27% of series (LGBM better in 15%; remainder negligible).

Chaotic traits of the real‑estate series (for interpretation): Hurst 0.65, Noise 0.99, Corr. dimension 1.33, max Lyapunov 0.01, KSE 1.84, (N_{Fh}=30).

.
├── src/
│   └── chaotic_library/
│       ├── __init__.py            # public API (EnhancedESN_FAN, chaos measures, version, etc.)
│       ├── enhanced_esn_fan.py    # ESN-FAN implementation
│       └── chaotic_measures.py    # Hurst, Lyapunov, entropy, dimensionality, etc.
├── tests/                         # unit tests
├── .github/workflows/             # CI (linting, tests)
├── requirements.txt               # runtime/dev dependencies
├── pyproject.toml                 # packaging metadata (build system, project info)
├── README.md
└── LICENSE

• Run linters/formatters before committing (e.g., ruff check . / ruff format .).
• Add/extend tests in tests/.
• For new feature blocks, include a minimal notebook demo.
• Keep figures reproducible from notebooks where possible.

License: MIT — see LICENSE.

Citation (placeholder): If you use this repository, please cite the corresponding preprint/paper.

@misc{esn_fan_2025,
  title   = {Enhanced Echo State Network with Fourier Analysis Network (FAN) Features},
  author  = {Kovantsev, A. and Vysotskiy, R.},
  year    = {2025},
  note    = {preprint},
  howpublished = {URL: add when available}
}

Q: My recursive (generative) forecast saturates or explodes. A: Increase ridge_alpha, decrease spectral_radius, and consider a slightly larger clip_value (2–5). Also try lowering leaking_rate for longer memory.

Q: Shapes? A: X must be 2D: (n_timesteps, input_dim). For univariate, reshape with reshape(-1, 1). y can be 1D (univariate) or 2D (multi‑output).

Q: Scaling consistency between train and predict? A: The model uses internal scalers. Ensure that polynomial and Fourier features at prediction time are computed in a way consistent with training. If you modify the code, apply the same input scaling before feature generation in all paths (teacher forcing and generative).