Recurrent neural networks (RNNs) offer linear-time scaling with sequence length while requiring only constant memory, yet they struggle to capture long-range dependencies due to vanishing gradients and limited receptive fields. To address these limitations, we introduce a second-order recurrent model in which the standard neuron-to-neuron communication is replaced by a spatially evolving field governed by (discretized) partial differential equations. Drawing inspiration from the role of cortical waves in brain computation, this mechanism allows structured spatiotemporal patterns to serve as an implicit, high-capacity memory. We show that the resulting model is equivalent to a structured infinite-order RNN in which the current state depends explicitly on its entire history of past states, yielding an effectively unbounded receptive field with a fixed number of parameters. We further derive constructive conditions to ensure marginal stability, constraining the gradient spectrum on the unit circle and thereby eliminating vanishing and exploding gradients. Empirically, the proposed architecture outperforms other recurrent models on long-horizon benchmarks while using substantially fewer parameters, demonstrating that spatial dynamics can effectively bridge the gap between efficient inference and long-term memory.
Figures & tables
Figure 1: Schematic of the model architectures. ( a ) SpatialRNN: The neuron states ht (represented by white nodes) are coupled through a shared spatiotemporal communication medium ψt (color coded), where information propagates locally over space and time. Each node writes to and reads from the medium, inducing wave-like interactions that encode past activity. ( b ) HORNN: Nodes correspond to neuron states ht , and interactions are explicitly represented by edges connecting the current state to a finite set of past states ht−k , for k=1,…,d . For d=1 , the model reduces to a standard RNN. Block diagrams of the corresponding discretized computations are shown on the right of each panel.
Figure 2: Gradient analysis of SpatialRNNs and HORNNs. ( a ) Eigenvalue spectrum of the gradient dynamics for standard RNNs (top) and SpatialRNNs (bottom). Red and green markers indicate the eigenvalues computed for D=I and D∈D , respectively. The spectrum of standard RNNs is sensitive to variations in D , whereas SpatialRNNs satisfying Theorem 3.3 allow for arbitrary eigenvalue placement with invariant magnitude, including on the unit circle. ( b ) Gradient norm versus time horizon. ( c ) Gradient norm (at lag m=20 ) versus the number of parameters (with 8 neurons). The results in (b,c) are averaged over 10 trajectories with uniformly sampled inputs.
Figure 3: Performance of SpatialRNN and HORNNs on the copy task problem, shown as a function of sequence length T .
Model (# units)
Acc. (%)
Params
unitary RNN (512) ( Arjovsky et al., 2016 )
91.4 (max)
9k
LSTM (256) ( Helfrich et al., 2018 )
92.9 (max)
270k
GRU (256) ( Chang et al., 2017 )
94.1 (max)
200k
NWM (256) ( Keller and Welling, 2023 )
94.8±1.1
50k
coRNN (256) ( Keller and Welling, 2023 )
95.0±2.4
134k
SpatialRNN (88)
95.5±0.2
13k
Table 1: psMNIST performance. Mean ± std over 3 runs when available; otherwise, maximum accuracy. Only single-layer architectures are accounted for.
Model (# units)
Acc. (%)
Params
LSTM (128) ( Kag et al., 2020 )
11.6
116k
GRU (128) ( Keller et al., 2024 )
43.8
88k
iRNN (529) ( Keller et al., 2024 )
51.3
336k
wRNN (256) ( Keller et al., 2024 )
55.0
435k
SpatialRNN (88)
57.2
21k
LRNN (128) ( Keller et al., 2024 )
57.4
46k
Table 2: npCIFAR10 performance for T=1000 . Maximum test accuracy over 3 trials.
Figure 4: SpatialRNN validation accuracy on the npCIFAR10 task for varying T .
Appendix figures & tables8 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 5: Evolution of the parameters {αi}i=1k during training on the copy task with sequence length T=128 . The dashed lines show the learned values of αi for a SpatialRNN architecture with block sizes [8, 8, 8, 8], while the solid line indicates training accuracy. Performance remains at chance-level until α4→1 , after which learning rapidly progresses.
Figure 6: Gradient norm ∥∂hT−t∂L∥ as a function of the time horizon, shown for the psMNIST classification task. Matrix norms are set such that ∣∣Bi∣∣2<0.4,∣∣Wψ,i∣∣2<0.4 and ∣∣Wi∣∣2<1.8 .
Figure 7: Frobenius norm of the exact state-to-state Jacobians Gmτ=∂hτ+m/∂hτ , propagated via Eqs. ( 13 )–( 14 ) along 100 psMNIST validation sequences ( τ=0 , T=784 ), for a trained SpatialRNN (left) and the same architecture at initialization (right). Thin lines show individual sequences, while the black line and the shaded band show the geometric mean and ±1 standard deviation (in log scale).
Model
Log-linear slope (per step)
Change over 783 steps
CV of ∥Gmτ∥F
Trained SpatialRNN
(7.98±7.75)⋅10−5
≈6%
0.071±0.010
Untrained SpatialRNN
(2.60±1.51)⋅10−4
≈23%
0.072±0.032
Appendix
Table 3: Empirical LTV gradient statistics over 100 psMNIST validation sequences (mean ± standard deviation across sequences). Slopes are fitted for m≥50 , and the coefficient of variation (CV) measures the dispersion of ∥Gmτ∥F along each sequence.
Figure 8: Relative Frobenius error ∥Gmexact−Gmrec∥F/∥Gmexact∥F between the exact Jacobians Gm=∂hT−1/∂hT−1−m , computed from the stored forward states, and those computed from the backward-reconstructed states, as a function of the time lag m , in single (float32) and double (float64) precision. (a) Non-dissipative SpatialRNN (single block, n=40 , α1=1 ): the error remains at machine precision over the entire horizon, with no exponential accumulation. (b) SpatialRNN with dissipative blocks (block sizes {22,6,6,6} , α=[1,0.5,0.5,0.5] ): round-off is amplified exponentially, at a measured rate of 1.42 per backward step, matching the predicted αmin−1/2=1.414 . Vertical lines mark the lags at which the median error first exceeds 10−3 , yielding m=44 in single and m=101 in double precision. Since the reconstruction starts from (ψT−1,ψT−2) and proceeds backward, the state at time t is recovered after T−1−t backward steps, so m is also the depth of the backward recursion. Curves show the median over 20 input sequences of length T=784 and shaded bands the 10 – 90 percentile range; the same effective matrices are used in both precisions, so the two curves differ only in the arithmetic precision of the recursions. Both models are randomly initialized, with the diagonal blocks rescaled so that condition 2 of Theorem 3.3 holds.
Model
Neurons
Parameters
Mean Accuracy(%) ± std.
SpatialRNN [32, 8, 8, 8]
56
6,815
99.04 ± 0.07
SpatialRNN [8, 8, 8, 8]
32
2,915
96.90 ± 4.34
SpatialRNN [4, 4, 4, 4]
16
1,179
69.33 ± 16.54
HORNN (T=[1,5,10,15])
64
17,744
13.39 ± 3.07
HORNN (T=[1,5,10])
64
13,648
14.37 ± 2.72
HORNN (T=[1,5])
64
9,552
11.05 ± 1.20
Appendix
Table 4: Copy task architecture configurations and performance. Mean accuracy is evaluated at the maximum sequence length T=1024 . Average and standard deviation over 10 runs.
Model
Neurons
Parameters
Mean Accuracy
SpatialRNN [128, 8, 8, 8]
152
36,497
96.34±0.09
SpatialRNN [64, 8, 8, 8]
88
12,849
95.49±0.22
SpatialRNN [32, 8, 8, 8]
56
5,633
94.37±0.13
Appendix
Table 5: psMNIST task architecture configurations and performance. Mean accuracy is evaluated over 3 runs, ordered by mean accuracy.
Model
Neurons
Parameters
Max Accuracy
SpatialRNN [128, 8, 8, 8]
152
50,937
57.8
SpatialRNN [64, 8, 8, 8]
88
21,209
57.2
SpatialRNN [32, 8, 8, 8]
56
10,953
54.1
SpatialRNN [16, 8, 8, 8]
40
6,977
50.8
SpatialRNN [8, 8, 8, 8]
32
5,277
45.6
Appendix
Table 6: npCIFAR10 task architecture configurations and performance. Maximum accuracy is evaluated at the sequence length T=1000 over 3 runs.
College of Pharmacy Chungnam National University Daejeon, South Korea · Department of Computer Science & Engineering Chungnam National University Daejeon, South Korea