Cross-Scale Reservoir Computing for large spatio-temporal forecasting and modeling

TL;DR

This paper tackles forecasting large-scale, high-resolution spatiotemporal data by capturing dynamics across multiple spatial scales. It introduces a cross-scale reservoir computing framework with hierarchical layers spanning coarse-to-fine grids, where information flows top-down from coarser layers and bottom-up via local inputs; each layer is trained sequentially, reducing computational load. On sea surface temperature data, the method yields improved long-horizon forecasts relative to a single-layer RC, with gains driven by slow quasi-periodic modes that propagate across layers and by an emergent linear regime in coarser layers described by $\tau \dot{\mathbf r} \approx (\mathbf{W} + \mathbf{W}_{\text{in}}\mathbf{W}_{\text{out}})\mathbf{r} - \mathbf{r}$. This linearization enables a modal decomposition via eigenpairs of the effective connectivity $\tilde{\mathbf W}=\mathbf{Z}\boldsymbol{\Lambda}\mathbf{Z}^{-1}$, offering interpretable insights and potential computational advantages for multi-resolution forecasting and other domains such as neuroscience.

Abstract

We propose a new reservoir computing method for forecasting high-resolution spatiotemporal datasets. By combining multi-resolution inputs from coarser to finer layers, our architecture better captures both local and global dynamics. Applied to Sea Surface Temperature data, it outperforms standard parallel reservoir models in long-term forecasting, demonstrating the effectiveness of cross-layers coupling in improving predictive accuracy. Finally, we show that the optimal network dynamics in each layer become increasingly linear, revealing the slow modes propagated to subsequent layers.

Figures (3)

• Figure 1: a, A dataset is organized into different resolution levels, forming a hierarchical architecture. b, Sketch of a three-layer hierarchical model architecture. Each resolution level compose a two-dimensional spatial domain represented as a grid of coordinates, each corresponding to an independent time series of a dynamical variable. The domain is partitioned into local subregions each serving as the learning target for a distinct reservoir network. Information from distant regions is passed through the previous layer at a lower resolution. In particular, looking at panel c, In each layer, the highlighted reservoir $\mathcal{R}$ (superscript indicates layer hierarchy) receives additional input from the output of the corresponding reservoir in the layer above (red contour), excluding contributions from the region that overlaps directly with its own coordinates (white area). This ensures that only non-local coarse-grained dynamics are passed downward, promoting hierarchical context without redundancy. The orange area highlights the neighboring coordinates which are included as additional input features to incorporate local spatial context as in pathak2018model.
• Figure 2: a, Average total error ($\mathrm{RMSE}_{t\le T}$) as a function of time $T$ for different numbers of layers, $N_L$. b, Ratio of the $\mathrm{RMSE}_{t\le T}$ for various $N_L$ to the reference case with $N_L = 3$, shown as a function of time. c, Absolute autocorrelation of the dataset versus time lag, averaged over all spatial locations, for different numbers of removed principal components, $P$. d, Ratio of the $\mathrm{RMSE}_{t\le T}$ at a fixed time of 10 weeks, plotted as a function of $P$. e, Difference in $\mathrm{RMSE}_{t \leq T}$ between $N_L = 3$ and $N_L = 1$ for individual geographical locations at $T = 50$ weeks, averaged over 25 runs. In all panels, solid lines denote the mean and shaded areas indicate the SEM, computed over $30$ independent runs (unless otherwise stated).
• Figure 3: a, Maximum network activity across all nodes for three different layers (L1,L2,L3). The solid line represents the mean and the shaded area represents the standard deviation across 25 independent system realizations. b, The eigenvalue spectrum of the effective reservoir that reproduces the dynamics of Layer 1 (L1). The color of each point indicates its decomposition weight. The data shown is from a single system realization. c, Decomposition weights distribution for the effective reservoir shown in (b). d, The real (black) and imaginary (red, blue) components of the complex coniugate linear mode pair associated with the highest decomposition weight.