Learning infinite context windows in recurrent architectures via spatial neural computing

# Learning infinite context windows in recurrent architectures via spatial neural computing
Source: [https://arxiv.org/html/2610.10690](https://arxiv.org/html/2610.10690)
Arthur N\. Montanariarthur\.montanari@northwestern\.eduAffiliation:Northwestern UniversityEarl K\. Millerekmiller@mit\.eduAffiliation:Massachusetts Institute of TechnologyAdilson E\. Mottermotter@northwestern\.eduAffiliation:Northwestern UniversityJorge Cortésjcortes@ucsd\.eduAffiliation:UC San Diego

###### Abstract

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\.

## 1Introduction

For decades, recurrent neural networks \(RNNs\) have served as the cornerstone of sequence modeling\. These models rely on hidden states that are recurrently updated to process sequential inputs, acting as a dynamical memory of past information\. In contrast to stateless architectures, RNNs offer the theoretical promise of compressing an unbounded context into a fixed\-size state, resulting inO⁡\(1\)O\(1\)memory cost with respect to the sequence lengthTTandO⁡\(T\)O\(T\)computational complexity\. However, the practical efficacy of classical RNNs is severely constrained by the vanishing and exploding gradient problem\([Bengio et al\., 1994](https://arxiv.org/html/2610.10690#bib.bib5)\)\. From a dynamical systems perspective, this issue stems from the spectrum of the recurrence matrix\. When its eigenvalues are stable \(\|λ\|<1\|\\lambda\|<1\) or unstable \(\|λ\|\>1\|\\lambda\|\>1\), perturbations respectively decay or grow exponentially, leading to a loss of long\-term dependencies during training\. While unitary and orthogonal RNNs\([Arjovsky et al\., 2016](https://arxiv.org/html/2610.10690#bib.bib25);[Ceni, 2025](https://arxiv.org/html/2610.10690#bib.bib24);[Henaff et al\., 2016](https://arxiv.org/html/2610.10690#bib.bib2);[Kerg et al\., 2019](https://arxiv.org/html/2610.10690#bib.bib4);[Le et al\., 2015](https://arxiv.org/html/2610.10690#bib.bib39)\)can mitigate this problem by constraining the eigenvalues to bemarginally stable\(\|λ\|=1\|\\lambda\|=1\), the nonlinearities of the activation functions perturb the effective dynamics away from the unit circle, reintroducing vanishing/exploding gradients\. Moreover, such constraints limit the model expressivity and convergence speed\([Kerg et al\., 2019](https://arxiv.org/html/2610.10690#bib.bib4);[Ribeiro et al\., 2020](https://arxiv.org/html/2610.10690#bib.bib23)\)\. To bypass these spectral issues, architectures such as the long short\-term memory \(LSTM\)\([Hochreiter and Schmidhuber, 1997](https://arxiv.org/html/2610.10690#bib.bib9)\)and gated recurrent unit \(GRU\)\([Cho et al\., 2014](https://arxiv.org/html/2610.10690#bib.bib10)\)introduce gating mechanisms with linear gradient updates\. These additive gates prevent gradients from exponentially decaying due to repetitive multiplicative operations, but they also come at the cost of increased architecture complexity and parameter count, leading to higher computational overhead\([Chung et al\., 2014](https://arxiv.org/html/2610.10690#bib.bib26)\), limited parallelization\([Li et al\., 2018](https://arxiv.org/html/2610.10690#bib.bib7)\), and implicit timescale biases\([Neil et al\., 2016](https://arxiv.org/html/2610.10690#bib.bib27)\)\.

Recently, the transformer architecture has largely supplanted RNNs due to its ability to model long\-range dependencies via self\-attention\. However, the transformer’s supremacy comes at a steep cost: the key\-value memory grows linearly,𝒪⁡\(T\)\\mathcal\{O\}\(T\), during inference and the attention mechanism scales quadratically,𝒪⁡\(T2\)\\mathcal\{O\}\(T^\{2\}\)\([Kwon et al\., 2023](https://arxiv.org/html/2610.10690#bib.bib29);[Vaswani et al\., 2017](https://arxiv.org/html/2610.10690#bib.bib6)\)\. These scaling properties effectively reintroduce the “von Neumann gap” in modern artificial intelligence \(AI\) as a significant portion of inference time and energy is spent on memory access rather than computation\([Dao et al\., 2022](https://arxiv.org/html/2610.10690#bib.bib28)\)\. This computational bottleneck has renewed interest in dynamical\-systems\-based architectures\([Sieber et al\., 2024](https://arxiv.org/html/2610.10690#bib.bib37)\), such as state\-space models \(SSMs\)\([Gu and Dao, 2024](https://arxiv.org/html/2610.10690#bib.bib12);[Gu et al\., 2022](https://arxiv.org/html/2610.10690#bib.bib11);[Somvanshi et al\., 2025](https://arxiv.org/html/2610.10690#bib.bib32)\)and oscillator networks\([Effenberger et al\., 2025](https://arxiv.org/html/2610.10690#bib.bib36);[Keller and Welling, 2023](https://arxiv.org/html/2610.10690#bib.bib14);[Keller et al\., 2024](https://arxiv.org/html/2610.10690#bib.bib18);[Rusch and Rus, 2025](https://arxiv.org/html/2610.10690#bib.bib22);[Rusch et al\., 2022](https://arxiv.org/html/2610.10690#bib.bib13)\)\. Although these approaches improve memory\-computation trade\-offs, they typically rely on dissipative dynamics\. As a result, they do not provide mechanisms to maintain nonvanishing gradients overarbitrarily long horizons\([Gu and Dao, 2024](https://arxiv.org/html/2610.10690#bib.bib12);[Rusch and Mishra, 2021a](https://arxiv.org/html/2610.10690#bib.bib3)\)\. Moreover, controlling memory\-retention and forgetting mechanisms in these models remains challenging, often requiring increased parameterization or complex hybrid architectures to achieve sufficient expressivity\([Karuvally et al\., 2026](https://arxiv.org/html/2610.10690#bib.bib34);[Wang et al\., 2025](https://arxiv.org/html/2610.10690#bib.bib35)\)\.

Taken together, these limitations raise a fundamental question:What structural features allow recurrent architectures to maintain nonvanishing gradient modes over long horizons despite input\-dependent nonlinear perturbations?Through a spectral analysis of the gradient dynamics, we show that the answer lies inhowthe nonlinearity enters the gradient equations, rather than in the spectrum of the recurrence matrices\. In conventional first\-order recurrent architectures \(e\.g\., RNNs\), these perturbations directly affect the coefficients that control modal decay, making robust marginal stability impossible regardless of the parameterization\. In contrast, a second\-order recurrent architecture can have additional, input\-independent coefficients that selectively fix the magnitude of gradient modes to enforce marginal stability while allowing inputs to perturb the remaining dynamics\. Such a separation arises naturally, for example, in wave propagation dynamics, since its inertial term is determined by the medium itself rather than by the input signal driving it\. This connection suggests that wave dynamics may be potentially used as a natural computational substrate for memory retention—an idea further motivated by increasing evidence on the role of cortical traveling waves in encoding, maintaining, and transforming information through distributed spatiotemporal activity patterns\([Benigno et al\., 2023](https://arxiv.org/html/2610.10690#bib.bib20);[Davis et al\., 2021](https://arxiv.org/html/2610.10690#bib.bib19);[Ermentrout and Kleinfeld, 2001](https://arxiv.org/html/2610.10690#bib.bib16);[Singer, 2021](https://arxiv.org/html/2610.10690#bib.bib17)\)\.

In this paper, we propose aspatial computing model, called SpatialRNN, in which recurrent neuron\-to\-neuron communication is mediated through a spatially evolving medium governed by \(discretized\) partial differential equations \(PDEs\)\. This mechanism allows structured spatiotemporal patterns to be encoded and processed directly in the state of the medium, providing a dynamical substrate for retaining long\-term temporal dependencies\. In fact, we rigorously prove that the SpatialRNN is mathematically equivalent to an RNN in which the current statehth\_\{t\}depends explicitly on all its past states\{ht−1,ht−2,ht−3,…\}\\\{h\_\{t\-1\},h\_\{t\-2\},h\_\{t\-3\},\\ldots\\\}\. The SpatialRNN therefore provides an implicit and parameter\-efficient realization of an unbounded receptive field\. Crucially, the second\-order gradient recurrence induced by the spatial medium can be designed to support nonvanishing gradient modes over arbitrarily long time horizons, ensuring robust marginal stability\. To achieve this, we derive constructive conditions that constrain the \(frozen\-time\) gradient modes on circles of prescribed radii for every admissible nonlinear perturbation, providing explicit control over memory retention \(unit radius\) and forgetting \(smaller radii\) through tunable decay rates\. Our empirical results show that the gradient norms of trained SpatialRNNs remain nearly constant over long input trajectories and that the SpatialRNN achieves competitive performance on long\-horizon benchmarks while often requiring considerably fewer parameters, highlighting the role of spatial computing as an efficient computational paradigm\.

## 2Spatial Computing Model

We propose the SpatialRNN, a recurrent architecture designed to store long\-term dependencies in sequential data by fundamentally changing how neurons communicate among themselves\. In standard RNNs, information is propagated through a learned recurrence matrixWhW\_\{h\}, resulting in instantaneous, all\-to\-all interactions between hidden \(neuron\) statesht∈ℝnh\_\{t\}\\in\\mathbb\{R\}^\{n\}\. Instead, the SpatialRNN replaces this direct neuron\-to\-neuron communication with a dynamic communication mediumψ⁡\(x,t\)\\psi\(x,t\)through which information propagates over space and time\. Each neuron interacts locally with this medium, which evolves according to a PDE\. This induces structured, wave\-like information propagation, allowing information to persist within the evolving medium rather than being repeatedly compressed into a finite\-dimensional hidden state\.

Formally, the medium evolves according to the second\-order linear PDE:

ρ⁡\(x\)​∂2ψ∂t2\+γ⁡\(x\)​∂ψ∂t=ℒ⁡\(ψ\)\+β⁡\(h⁡\(t\)\),\\rho\(x\)\\frac\{\\partial^\{2\}\\psi\}\{\\partial t^\{2\}\}\+\\gamma\(x\)\\frac\{\\partial\\psi\}\{\\partial t\}=\\mathcal\{L\}\(\\psi\)\+\\beta\(h\(t\)\),\(1\)wherex∈ℝ2x\\in\\mathbb\{R\}^\{2\}denotes the 2\-dimensional spatial coordinate andψ⁡\(x,t\)\\psi\(x,t\)is the state of the communication medium\. Here,ℒ\\mathcal\{L\}is a spatial linear operator modeling local interactions in the medium \(e\.g\.,ℒ\\mathcal\{L\}could be the LaplacianΔ=∂2ψ∂x12\+∂2ψ∂x22\\Delta=\\frac\{\\partial^\{2\}\\psi\}\{\\partial x\_\{1\}^\{2\}\}\+\\frac\{\\partial^\{2\}\\psi\}\{\\partial x\_\{2\}^\{2\}\}\),ρ⁡\(x\)\\rho\(x\)andγ⁡\(x\)\\gamma\(x\)are spatially\-dependent functions controlling the local propagation and damping dynamics, respectively, andβ\\betamaps the neural activityh⁡\(t\)h\(t\)onto the spatial domain, acting as a forcing term on the mediumψ⁡\(x,t\)\\psi\(x,t\)\.

To implement this model in discrete time, we discretize the PDE \([1](https://arxiv.org/html/2610.10690#S2.E1)\), yielding the following second\-order linear difference equation \(LDE\):

ψt=W¯​ψt−1\+W¯¯​ψt−2\+B​ht−1,\\psi\_\{t\}=\\overline\{W\}\\psi\_\{t\-1\}\+\\overline\{\\overline\{W\}\}\\psi\_\{t\-2\}\+Bh\_\{t\-1\},\(2\)whereW¯,W¯¯∈ℝn×n\{\\overline\{W\}\},\\overline\{\\overline\{W\}\}\\in\\mathbb\{R\}^\{n\\times n\}define the medium dynamics, andB∈ℝn×nB\\in\\mathbb\{R\}^\{n\\times n\}acts as a forcing term driven by the neural stateht∈ℝnh\_\{t\}\\in\\mathbb\{R\}^\{n\}\. The SpatialRNN integrates this medium into the recurrent computation structure of an RNN, leading to our proposed model architecture:

SpatialRNN​\{ht=Φ⁡\(Wx​xt\+Wψ​ψt\+bh\),ψt=W¯​ψt−1\+W¯¯​ψt−2\+B​ht−1,yt=softmax​\(Wyh​ht\+Wyψ,1​ψt\+Wyψ,2​ψt−1\+by\)\.\\displaystyle\\text\{SpatialRNN\}\\begin\{cases\}h\_\{t\}=\\Phi\(W\_\{x\}x\_\{t\}\+W\_\{\\psi\}\\psi\_\{t\}\+b\_\{h\}\),\\\\ \\psi\_\{t\}=\\overline\{W\}\\psi\_\{t\-1\}\+\\overline\{\\overline\{W\}\}\\psi\_\{t\-2\}\+Bh\_\{t\-1\},\\\\ y\_\{t\}=\\text\{softmax\}\(W\_\{y\}^\{h\}h\_\{t\}\+W^\{\\psi,1\}\_\{y\}\\psi\_\{t\}\+W^\{\\psi,2\}\_\{y\}\\psi\_\{t\-1\}\+b\_\{y\}\)\.\\end\{cases\}\(3\)Here,xt∈ℝmx\_\{t\}\\in\\mathbb\{R\}^\{m\}is the time\-dependent input signal,ht∈ℝnh\_\{t\}\\in\\mathbb\{R\}^\{n\}is the hidden neuron state, andyt∈ℝpy\_\{t\}\\in\\mathbb\{R\}^\{p\}is the output state\. The model is parameterized by the \(trainable\) input\-to\-state weightsWx∈ℝn×mW\_\{x\}\\in\\mathbb\{R\}^\{n\\times m\}, medium\-to\-state weightsWψ∈ℝn×nW\_\{\\psi\}\\in\\mathbb\{R\}^\{n\\times n\}, medium weightsW¯,W¯¯∈ℝn×n\{\\overline\{W\}\},\\overline\{\\overline\{W\}\}\\in\\mathbb\{R\}^\{n\\times n\}, output weightsWyh,Wyψ,1,Wyψ,2∈ℝp×nW\_\{y\}^\{h\},W^\{\\psi,1\}\_\{y\},W^\{\\psi,2\}\_\{y\}\\in\\mathbb\{R\}^\{p\\times n\}, and biasesbh∈ℝnb\_\{h\}\\in\\mathbb\{R\}^\{n\},by∈ℝpb\_\{y\}\\in\\mathbb\{R\}^\{p\}\. The activation functionΦ:ℝn↦ℝn\\Phi\\mathrel\{\\mathop\{\\mathchar 58\\relax\}\}\\mathbb\{R\}^\{n\}\\mapsto\\mathbb\{R\}^\{n\}is applied elementwise, withΦ⁡\(z\)=\(Φ⁡\(z1\),…,Φ⁡\(zn\)\)⊤\\Phi\(z\)=\(\\Phi\(z\_\{1\}\),\\ldots,\\Phi\(z\_\{n\}\)\)^\{\\top\}, and can represent any standard nonlinearity \(e\.g\., tanh or ReLU\)\. For classification tasks, a softmax output is taken without loss of generality\. For initialization, we seth0=ψ0=ψ−1=0h\_\{0\}=\\psi\_\{0\}=\\psi\_\{\-1\}=0\. The mediumψt∈ℝn\\psi\_\{t\}\\in\\mathbb\{R\}^\{n\}serves as aspatiotemporal memory, allowing information to propagate and interact over time without repeated compression\. Fig\.[1](https://arxiv.org/html/2610.10690#S2.F1)a illustrates the spatiotemporal dynamics of the SpatialRNN with a 2\-dimensional medium\.

![Refer to caption](https://arxiv.org/html/2610.10690v1/fig_diagram.png)Figure 1:Schematic of the model architectures\. \(a\)SpatialRNN:The neuron stateshth\_\{t\}\(represented by white nodes\) are coupled through a shared spatiotemporal communication mediumψt\\psi\_\{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 stateshth\_\{t\}, and interactions are explicitly represented by edges connecting the current state to a finite set of past statesht−kh\_\{t\-k\}, fork=1,…,dk=1,\\ldots,d\. Ford=1d=1, the model reduces to a standard RNN\. Block diagrams of the corresponding discretized computations are shown on the right of each panel\.### 2\.1Higher\-Order and Standard RNNs

A straightforward approach to extend the long\-term memory of an RNN is to introduce direct dependencies of the current statehth\_\{t\}on multiple time\-delayed states\{ht−1,ht−2,…​ht−d\}\\\{h\_\{t\-1\},h\_\{t\-2\},\\ldots h\_\{t\-d\}\\\}\. This concept leads tohigher\-order RNNs\(HORNNs\)\([Soltani and Jiang, 2016](https://arxiv.org/html/2610.10690#bib.bib1)\), which generalize standard RNNs by explicitly incorporating a finite history of past hidden states \(Fig\.[1](https://arxiv.org/html/2610.10690#S2.F1)b\)\. In this formulation, the update rule is written as a discrete convolution over past states, parameterized by a set of \(trainable\) recurrence matrices𝒲=\{Wk\}k=1d\\mathcal\{W\}=\\\{W\_\{k\}\\\}\_\{k=1\}^\{d\}, withWk∈ℝn×nW\_\{k\}\\in\\mathbb\{R\}^\{n\\times n\}:

HORNN​\{ht=Φ⁡\(Wx​xt\+∑k=1dWk​ht−k\+bh\),yt=softmax​\(Wyh​ht\+by\)\.\\displaystyle\\text\{HORNN\}\\begin\{cases\}h\_\{t\}=\\Phi\(W\_\{x\}x\_\{t\}\+\\sum\_\{k=1\}^\{d\}W\_\{k\}h\_\{t\-k\}\+b\_\{h\}\),\\\\ y\_\{t\}=\\text\{softmax\}\(W\_\{y\}^\{h\}h\_\{t\}\+b\_\{y\}\)\.\\end\{cases\}\(4\)In practice, HORNNs exhibit a clear trade\-off between memory capacity and parameter efficiency\. The orderd∈ℕd\\in\\mathbb\{N\}defines the temporal receptive field of the model, but also increases the number of parameters linearly, since each additional lagkkintroduces a new matrixWkW\_\{k\}\. Moreover, for any finitedd, dependencies beyond lagddare not modeled directly and must propagate through intermediate states, leading again to exponential decay of long\-range memories for all lags beyonddd\(see results in Sec\.[3\.1](https://arxiv.org/html/2610.10690#S3.SS1)and Fig\.[2](https://arxiv.org/html/2610.10690#S3.F2)for detailed theoretical analysis\)\. Ford=1d=1, the HORNN model reduces to thestandard RNN:

Standard RNN​\{ht=Φ⁡\(Wx​xt\+Wh​ht−1\+bh\),yt=softmax​\(Wy​ht\+by\)\.\\displaystyle\\text\{Standard RNN\}\\begin\{cases\}h\_\{t\}=\\Phi\(W\_\{x\}x\_\{t\}\+W\_\{h\}h\_\{t\-1\}\+b\_\{h\}\),\\\\ y\_\{t\}=\\text\{softmax\}\(W\_\{y\}h\_\{t\}\+b\_\{y\}\)\.\\end\{cases\}\(5\)

### 2\.2Equivalence between SpatialRNNs and∞\\infty\-order RNNs

The effective temporal context window of a finite\-order HORNN \([4](https://arxiv.org/html/2610.10690#S2.E4)\) is limited by its receptive fielddd\. Extending it to an infinite horizond→∞d\\to\\inftytheoretically solves this limitation but becomes computationally intractable for long\-sequence tasks due to the prohibitively large parameter and memory requirements\. Our key insight is that the SpatialRNN provides an implicit parameterization of this infinite history\. Thus, instead of storing past states or learning a separate matrixWkW\_\{k\}for every lagkk, the SpatialRNN generates the effective convolution kernel\{Wk\}k=1∞\\\{W\_\{k\}\\\}\_\{k=1\}^\{\\infty\}via the dynamics of the spatial medium, as formalized below\.

###### Theorem 2\.3\(Equivalence between medium state and infinite history\)\.

Let\{ht\}t≥1\\\{h\_\{t\}\\\}\_\{t\\geq 1\}be a sequence of hidden states and letψt\\psi\_\{t\}be governed by a second\-order LDE \([2](https://arxiv.org/html/2610.10690#S2.E2)\) with initial conditionsψ0=ψ−1=0\\psi\_\{0\}=\\psi\_\{\-1\}=0\. Then, the spatial stateψt\\psi\_\{t\}is equivalent to a convolution over the entire history of hidden states:

ψt=∑k=1t−1Wk​ht−k,\\psi\_\{t\}=\\sum\_\{k=1\}^\{t\-1\}W\_\{k\}h\_\{t\-k\},\(6\)where the sequence of weight matrices\{Wk\}k≥1\\\{W\_\{k\}\\\}\_\{k\\geq 1\}are defined recursively by

W1=B,W2=W¯B,andWk=W¯Wk−1\+W¯¯Wk−2,k≥3\.W\_\{1\}=B,\\quad W\_\{2\}=\\overline\{W\}B,\\quad\\text\{and\}\\quad W\_\{k\}=\\overline\{W\}W\_\{k\-1\}\+\\overline\{\\overline\{W\}\}W\_\{k\-2\},\\,\\,k\\geq 3\.\(7\)

###### Proof\.

See Appendix[A\.1](https://arxiv.org/html/2610.10690#A1.SS1)\. The result follows by recursively expanding \([2](https://arxiv.org/html/2610.10690#S2.E2)\) and matching the coefficients for each laght−kh\_\{t\-k\}\. In particular, the sequence of matrices \([7](https://arxiv.org/html/2610.10690#S2.E7)\) corresponds to the impulse response\([Ogata and others, 1995](https://arxiv.org/html/2610.10690#bib.bib50)\)of the LDE \([2](https://arxiv.org/html/2610.10690#S2.E2)\) \(i\.e\., the system response to a unit input at each lag\)\. ∎

Theorem[2\.3](https://arxiv.org/html/2610.10690#S2.Thmtheorem3)establishes a direct connection between the SpatialRNN and the∞\\infty\-order RNN, highlighting two fundamental advantages of our approach ast→∞t\\rightarrow\\infty\. Onmemory efficiency, while the∞\\infty\-order RNN requires full history\{hk\}k=1t−1\\\{h\_\{k\}\\\}\_\{k=1\}^\{t\-1\}to compute the next statehth\_\{t\}, the SpatialRNN compresses this entire history into the medium statesψt\\psi\_\{t\}andψt−1\\psi\_\{t\-1\}, restoring𝒪⁡\(1\)\\mathcal\{O\}\(1\)memory complexity during inference\. Onparameter efficiency, an∞\\infty\-order RNN requires infinitely many unique parameters\{Wk\}k=1∞\\\{W\_\{k\}\\\}\_\{k=1\}^\{\\infty\}, while the SpatialRNN uses only four fixed\-size matrices\{W¯,W¯¯,B,Wψ\}\\\{\\overline\{W\},\\overline\{\\overline\{W\}\},B,W\_\{\\psi\}\\\}to implicitly encode this infinite sequence of matrices\. The SpatialRNN can thus be viewed as a feasible realization of astructured∞\\infty\-order RNN that is both memory and parameter efficient\. By virtue of this equivalence, the SpatialRNN has an unbounded receptive field, allowing hidden states at arbitrary lags to influence the current state; whether this influence is retained during training is determined by the gradient dynamics, as we discuss next\.

## 3Gradient Dynamics and Conditions for Long\-Term Memory

RNNs are typically trained via backpropagation through time \(BPTT\)\. For a standard first\-order RNN architecture, the gradient of a sequence\-to\-output lossℒ\\mathcal\{L\}with respect to the weightsWWcan be expressed via the chain rule as

∂ℒ∂W=∂ℒ∂yT​∂yT∂hT​∑k=1T∂hT∂hk​∂hk∂W,\\frac\{\\partial\\mathcal\{L\}\}\{\\partial W\}=\\frac\{\\partial\\mathcal\{L\}\}\{\\partial y\_\{T\}\}\\frac\{\\partial y\_\{T\}\}\{\\partial h\_\{T\}\}\\sum\_\{k=1\}^\{T\}\\frac\{\\partial h\_\{T\}\}\{\\partial h\_\{k\}\}\\frac\{\\partial h\_\{k\}\}\{\\partial W\},\(8\)whereTTis the sequence length\. This decomposition makes explicit that the loss gradients are governed by the state\-to\-state Jacobian matrices∂hT∂hk\\partialderivative\{h\_T\}\{h\_k\}, which quantify how perturbations propagate across time\. We can recursively define themm\-step state transition Jacobian matrix as follows:

Gmt\\displaystyle G\_\{m\}^\{t\}:=∂ht\+m∂ht=∏k=0m−1∂ht\+k\+1∂ht\+k⟹Gmt=At\+m−1Gm−1t,\\displaystyle\\mathrel\{\\mathop\{\\mathchar 58\\relax\}\}=\\frac\{\\partial h\_\{t\+m\}\}\{\\partial h\_\{t\}\}=\\prod\_\{k=0\}^\{m\-1\}\\frac\{\\partial h\_\{t\+k\+1\}\}\{\\partial h\_\{t\+k\}\}\\quad\\implies\\quad G\_\{m\}^\{t\}=A\_\{t\+m\-1\}G\_\{m\-1\}^\{t\},\(9\)whereAt:=∂ht\+1∂ht∈ℝn×nA\_\{t\}\\mathrel\{\\mathop\{\\mathchar 58\\relax\}\}=\\frac\{\\partial h\_\{t\+1\}\}\{\\partial h\_\{t\}\}\\in\\mathbb\{R\}^\{n\\times n\}andG1t=∂ht\+1∂htG\_\{1\}^\{t\}=\\frac\{\\partial h\_\{t\+1\}\}\{\\partial h\_\{t\}\}\. Eq\. \([9](https://arxiv.org/html/2610.10690#S3.E9)\) constitutes adiscrete\-time linear time\-varying \(LTV\) dynamical systembecause the matricesAtA\_\{t\}depend on the time\-varying input sequencextx\_\{t\}\. As the temporal distancemmgrows, the asymptotic behavior of‖Gmt‖\\\|G\_\{m\}^\{t\}\\\|determines whether the model can effectively learn long\-term dependencies\. If‖Gmt‖\\\|G\_\{m\}^\{t\}\\\|decays to zero \(vanishing gradients\) or grows unboundedly \(exploding gradients\), the training fails to capture long\-horizon correlations\([Bengio et al\., 1994](https://arxiv.org/html/2610.10690#bib.bib5)\)\. Hereafter, we refer to the evolution ofGmtG\_\{m\}^\{t\}across time as the gradient dynamics\.

To build intuition, we recall that, for a discrete\-time linear time\-invariant \(LTI\) systemzm=A​zm−1z\_\{m\}=Az\_\{m\-1\}with spectral radiusρ⁡\(A\)=maxλ∈σ⁡\(A\)⁡\|λ\|\\rho\(A\)=\\max\_\{\\lambda\\in\\sigma\(A\)\}\|\\lambda\|, the state decays exponentially ifρ⁡\(A\)<1\\rho\(A\)<1\(vanishing gradients\), grows unboundedly ifρ⁡\(A\)\>1\\rho\(A\)\>1or if nonsimple eigenvalues lie on the unit circle \(exploding gradients\), and remains bounded and nondecaying ifρ⁡\(A\)=1\\rho\(A\)=1with all unit\-modulus eigenvalues semisimple \(marginal stability\)\. Our goal is to identify the structural conditions under which the gradient dynamics are marginally stable for every admissible input\-dependent perturbation, allowing for arbitrarily long memory retention\. However, establishing marginal stability conditions for arbitrary LTV systems, such as Eq\. \([9](https://arxiv.org/html/2610.10690#S3.E9)\), is generally challenging\([Sun, 2008](https://arxiv.org/html/2610.10690#bib.bib52)\)\. To circumvent this challenge, we adopt instead afrozen\-time approximation, in which we hold the instantaneous transition matrixAtA\_\{t\}fixed at each reference timettand analyze the spectrum of the resulting LTI system\. In contrast, the full LTV gradient dynamics evolve through the productGtm=At\+m−1⋯At\+1AtG^\{t\}\_\{m\}=A\_\{t\+m\-1\}\\cdot\\cdot\\cdot A\_\{t\+1\}A\_\{t\}\. The frozen\-time approximation thus serves as a tractable proxy for practical stability, since the marginal stability of every frozen\-time system is a necessary, but not sufficient, condition for the marginal stability of the full LTV system under arbitrary switching of the perturbations\.

In what follows, we show that finite\-order HORNNs, including standard RNNs, cannot robustly maintain marginal stability, and therefore inherently fall into regimes of vanishing or exploding gradients \(Sec\.[3\.1](https://arxiv.org/html/2610.10690#S3.SS1)\)\. In contrast, we prove that the SpatialRNN, under suitable conditions, allows for robust eigenvalue placement on the unit circle, thereby achieving robust marginal stability of the frozen\-time gradient dynamics \(Sec\.[3\.2](https://arxiv.org/html/2610.10690#S3.SS2)\)\.

### 3\.1Nonrobustness of RNN and HORNN Gradients

We now analyze the gradient dynamics of HORNNs\. The key issue of these models is the lack of robustness: even when parameters are carefully initialized to have marginally stable modes, small variations induced by inputs and nonlinearities inevitably push the system into vanishing or exploding regimes\. Adopting a frozen\-time approximation, we fix a reference timettand examine the gradient dynamics ofGmtG\_\{m\}^\{t\}\. As in Eq\. \([9](https://arxiv.org/html/2610.10690#S3.E9)\), we characterize the gradient dynamics of HORNNs through the state\-to\-state JacobianGmt=∂ht\+m∂htG^\{t\}\_\{m\}=\\partialderivative\{h\_\{t\+m\}\}\{h\_t\}, yielding the following recurrence equation:

Gmt=Dmt​∑k=1min⁡\{m,d\}Wk​Gm−kt,m≥1,G\_\{m\}^\{t\}=D\_\{m\}^\{t\}\\sum\\nolimits\_\{k=1\}^\{\\min\\\{m,d\\\}\}W\_\{k\}G\_\{m\-k\}^\{t\},\\quad m\\geq 1,\(10\)initialized withG0t=𝕀G\_\{0\}^\{t\}=\\mathbb\{I\}andGmt=0G\_\{m\}^\{t\}=0form<0m<0\. Here,DmtD\_\{m\}^\{t\}is a diagonal matrix defined as

Dmt=diag​\(Φ′​\(Wx​xt\+m\+∑k=1dWk​ht\+m−k\+bh\)\),D\_\{m\}^\{t\}=\\text\{diag\}\\Big\(\\Phi^\{\\prime\}\\Big\(W\_\{x\}x\_\{t\+m\}\+\\sum\\nolimits\_\{k=1\}^\{d\}W\_\{k\}h\_\{t\+m\-k\}\+b\_\{h\}\\Big\)\\Big\),\(11\)whereΦ′\\Phi^\{\\prime\}denotes the derivative of the activation functionΦ\\Phi\. In the special case of a standard RNN \(d=1d=1\), Eq\. \([10](https://arxiv.org/html/2610.10690#S3.E10)\) reduces toGmt=Dmt​Wh​Gm−1tG\_\{m\}^\{t\}=D\_\{m\}^\{t\}W\_\{h\}G\_\{m\-1\}^\{t\}\.

Crucially, for both HORNNs and standard RNNs, each matrixDmtD\_\{m\}^\{t\}depends on the time\-varying inputxt\+mx\_\{t\+m\}and can therefore vary over time in a manner that is effectively uncontrolled\. However, their entries are constrained by the range of the activation function derivative\. This implies that, at each time step,DmtD\_\{m\}^\{t\}lies in the set

𝒟=\{D∈ℝn×n\|D=diag\(d1,…,dn\),di∈\[Φmin′,Φmax′\]\},\\mathcal\{D\}=\\\{D\\in\\mathbb\{R\}^\{n\\times n\}\\,\|\\,D=\\operatorname\{diag\}\(d\_\{1\},\\ldots,d\_\{n\}\),\\,\\,d\_\{i\}\\in\[\\Phi^\{\\prime\}\_\{\\min\},\\Phi^\{\\prime\}\_\{\\max\}\]\\\},whereΦmin′=infz∈ℝΦ′​\(z\)\\Phi^\{\\prime\}\_\{\\min\}=\\inf\_\{z\\in\\mathbb\{R\}\}\\Phi^\{\\prime\}\(z\)andΦmax′=supz∈ℝΦ′​\(z\)\\Phi^\{\\prime\}\_\{\\max\}=\\sup\_\{z\\in\\mathbb\{R\}\}\\Phi^\{\\prime\}\(z\)denote the bounds of the functionΦ′\\Phi^\{\\prime\}\. To assess robustness, we consider the family of systems obtained by fixingDmt=DD\_\{m\}^\{t\}=D,∀m\\forall m, with some constant matrixD∈𝒟D\\in\\mathcal\{D\}\. Each such choice corresponds to a frozen\-time \(LTI\) systemGmt=D​∑kWk​Gm−ktG\_\{m\}^\{t\}=D\\sum\_\{k\}W\_\{k\}\\,G\_\{m\-k\}^\{t\}\. A necessary condition for robust marginal stability\([Liberzon, 2003](https://arxiv.org/html/2610.10690#bib.bib8)\)is that the frozen\-time gradient dynamics remain marginally stable for everyD∈𝒟D\\in\\mathcal\{D\}\.

###### Proposition 3\.1\(Spectrum of frozen\-time HORNN gradients\)\.

For any fixedD∈𝒟D\\in\\mathcal\{D\}, a valueλ∈ℂ\\lambda\\in\\mathbb\{C\}is an eigenvalue of the HORNN frozen\-time gradient dynamics \([10](https://arxiv.org/html/2610.10690#S3.E10)\) if and only if

det⁡\(𝕀−D​W^​\(λ\)\)=0,whereW^​\(λ\)=∑k=1dWk​λ−k\.\\det\\big\(\\mathbb I \- D\\hat W\(\\lambda\)\\big\)=0,\\quad\\text\{where\}\\quad\\hat\{W\}\(\\lambda\)=\\sum\\nolimits\_\{k=1\}^\{d\}W\_\{k\}\\lambda^\{\-k\}\.\(12\)

###### Proof\.

The proof exploits thez\{z\}\-transform of the linear recurrence system and is given in Appendix[A\.1](https://arxiv.org/html/2610.10690#A1.SS1)∎

The characteristic equation \([12](https://arxiv.org/html/2610.10690#S3.E12)\) reveals that the eigenvaluesλ\\lambdaare roots of a polynomial whose coefficients depend linearly on the diagonal entries ofDD\. Consequently, variations inDDdirectly perturb the spectrum, preventing robust marginal stability, as we formalize next\.

###### Proposition 3\.2\(Non\-existence of marginally stable HORNNs\)\.

Suppose the derivative of the activation functionΦ\\Phitakes values arbitrarily close to 0\. For any HORNN \([4](https://arxiv.org/html/2610.10690#S2.E4)\), defined by a set of weight matrices\{Wk\}k=1d\\\{W\_\{k\}\\\}\_\{k=1\}^\{d\}with finited∈ℕd\\in\\mathbb\{N\}, there exists a matrixD∈𝒟D\\in\\mathcal\{D\}such that the gradient dynamics are not marginally stable \(i\.e\., no rootλ\\lambdaof Eq\. \([12](https://arxiv.org/html/2610.10690#S3.E12)\) satisfies\|λ\|=1\|\\lambda\|=1\)\.

###### Proof\.

See proof in Appendix[A\.1](https://arxiv.org/html/2610.10690#A1.SS1)∎

Proposition[3\.2](https://arxiv.org/html/2610.10690#S3.Thmtheorem2)identifies a structural limitation of finite\-order recurrent architectures, since its proof relies only on the fact thatDDmultiplies the entire recurrence operator\. The result is therefore independent of the specific parameterization of the recurrence matrices and extends to a broad class of architectures, including orthogonal\([Helfrich et al\., 2018](https://arxiv.org/html/2610.10690#bib.bib30);[Henaff et al\., 2016](https://arxiv.org/html/2610.10690#bib.bib2)\), nonnormal\([Kerg et al\., 2019](https://arxiv.org/html/2610.10690#bib.bib4)\), and diagonal\([Li et al\., 2018](https://arxiv.org/html/2610.10690#bib.bib7)\)RNNs\. As a result, such parameterization constraints may mitigate or delay gradient decay, but they cannot make the gradient modes robust to nonlinear perturbations because the same input\-dependent factorDDdirectly perturbs the coefficients governing modal decay\. For instance, even orthogonal RNNs, which are initialized with eigenvalues on the unit circle\|λ⁡\(Wh\)\|=1\|\\lambda\(W\_\{h\}\)\|=1, are not robustly marginally stable once nonlinearities are introduced: the input\-dependent factorDDperturbs the transition matrix toD​WhDW\_\{h\}, generically shifting its spectrum away from the unit circle \(leading to contractive or expansive gradient modes depending on the realization ofDD\)\. Fig\.[2](https://arxiv.org/html/2610.10690#S3.F2)demonstrates the lack of robustness in HORNNs for finitedd, showing the vanishing of gradients over long sequences\.

![Refer to caption](https://arxiv.org/html/2610.10690v1/fig_gradients.png)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 forD=𝕀D=\\mathbb\{I\}andD∈𝒟D\\in\\mathcal\{D\}, respectively\. The spectrum of standard RNNs is sensitive to variations inDD, whereas SpatialRNNs satisfying Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)allow for arbitrary eigenvalue placement with invariant magnitude, including on the unit circle\. \(b\) Gradient norm versus time horizon\. \(c\) Gradient norm \(at lagm=20m=20\) versus the number of parameters \(with 8 neurons\)\. The results in \(b,c\) are averaged over 10 trajectories with uniformly sampled inputs\.
### 3\.2Robustness of SpatialRNN Gradients

The gradient dynamics in SpatialRNNs is given by \(see Appendix[A\.2](https://arxiv.org/html/2610.10690#A1.SS2)for the derivation\):

Gmt=Dmt​Wψ​Ψmt,m≥1,G\_\{m\}^\{t\}=D^\{t\}\_\{m\}W\_\{\\psi\}\\Psi\_\{m\}^\{t\},\\quad m\\geq 1,\(13\)whereDmt=diag​\(Φ′​\(Wx​xt\+m\+Wψ​ψt\+m\+bh\)\)D\_\{m\}^\{t\}=\\text\{diag\}\(\\Phi^\{\\prime\}\(W\_\{x\}x\_\{t\+m\}\+W\_\{\\psi\}\\psi\_\{t\+m\}\+b\_\{h\}\)\), andΨmt\\Psi\_\{m\}^\{t\}follows the recurrence

\(ΨmtΨm−1t\)=\(W¯\+B​Dm−1t​WψW¯¯𝕀0\)​\(Ψm−1tΨm−2t\),m≥2,\\displaystyle\\begin\{pmatrix\}\\Psi\_\{m\}^\{t\}\\\\ \\Psi\_\{m\-1\}^\{t\}\\end\{pmatrix\}=\\begin\{pmatrix\}\\overline\{W\}\+BD\_\{m\-1\}^\{t\}W\_\{\\psi\}&\\overline\{\\overline\{W\}\}\\\\ \\mathbb\{I\}&0\\end\{pmatrix\}\\begin\{pmatrix\}\\Psi\_\{m\-1\}^\{t\}\\\\ \\Psi\_\{m\-2\}^\{t\}\\end\{pmatrix\},\\quad m\\geq 2,\(14\)initialized withΨ1t=B\\Psi\_\{1\}^\{t\}=BandΨ0t=0\\Psi\_\{0\}^\{t\}=0\. In contrast to RNNs and HORNNs, the SpatialRNNs gradient dynamics are governed by thesecond\-orderrecurrence system \([14](https://arxiv.org/html/2610.10690#S3.E14)\)\. In this formulation, the nonlinearity termDmtD^\{t\}\_\{m\}acts as a structural perturbation to the gradient dynamics rather than a scaling factor applied at each step\. This distinction is crucial as the asymptotic behavior of the gradients is determined by the spectral properties of the transition matrix associated withΨmt\\Psi\_\{m\}^\{t\}, rather than by the instantaneous diagonal entries ofDmtD^\{t\}\_\{m\}\. Under a frozen\-time approximation, we can derive the sufficient conditions forrobustgradient behavior of the SpatialRNN\. In particular, the following result shows that the eigenvaluesλ∈ℂ\\lambda\\in\\mathbb\{C\}of the resulting companion matrix \([14](https://arxiv.org/html/2610.10690#S3.E14)\) can be constrained to lie on circles offixed radiusαi\\alpha\_\{i\}which are invariant to the specific realization ofDm−1t∈𝒟D\_\{m\-1\}^\{t\}\\in\\mathcal\{D\}\.

###### Theorem 3\.3\(Robustness of SpatialRNN gradients\)\.

Let the matricesW¯\\overline\{W\}andW¯¯\\overline\{\\overline\{W\}\}be defined as block upper\-triangular matrices with blocks sizesnin\_\{i\},i=1,…,ki=1,\\ldots,k, such that

W¯=\(W¯1⋆⋆0⋱⋆00W¯k\)andW¯¯=\(−α1​𝕀n1⋆⋆0⋱⋆00−αk​𝕀nk\),\\overline\{W\}=\\begin\{pmatrix\}\\overline\{W\}\_\{1\}&\\star&\\star\\\\ 0&\\ddots&\\star\\\\ 0&0&\\overline\{W\}\_\{k\}\\end\{pmatrix\}\\quad\\text\{and\}\\quad\\overline\{\\overline\{W\}\}=\\begin\{pmatrix\}\-\\alpha\_\{1\}\\mathbb\{I\}\_\{n\_\{1\}\}&\\star&\\star\\\\ 0&\\ddots&\\star\\\\ 0&0&\-\\alpha\_\{k\}\\mathbb\{I\}\_\{n\_\{k\}\}\\end\{pmatrix\},\(15\)whereW¯i∈ℝni×ni\\overline\{W\}\_\{i\}\\in\\mathbb\{R\}^\{n\_\{i\}\\times n\_\{i\}\},αi∈ℝ\>0\\alpha\_\{i\}\\in\\mathbb\{R\}\_\{\>0\},∑i=1kni=n\\sum\_\{i=1\}^\{k\}n\_\{i\}=nis the model size, and⋆\\stardenotes arbitrary entries\. Suppose the following conditions are satisfied:

1. 1\.Wi¯=Wi¯⊤\\overline\{W\_\{i\}\}=\\overline\{W\_\{i\}\}^\{\\top\},B=blkdiag⁡\{B1,…,Bk\}B=\\operatorname\{blkdiag\}\\\{B\_\{1\},\\ldots,B\_\{k\}\\\}, for block matricesBi∈ℝni×niB\_\{i\}\\in\\mathbb\{R\}^\{n\_\{i\}\\times n\_\{i\}\}, andWψ=Σ​B⊤=blkdiag⁡\{Wψ,1,…,Wψ,k\}W\_\{\\psi\}=\\Sigma B^\{\\top\}=\\operatorname\{blkdiag\}\\\{W\_\{\\psi,1\},\\ldots,W\_\{\\psi,k\}\\\}for some diagonal matrixΣ∈diag⁡\(ℝn\)\\Sigma\\in\\operatorname\{diag\}\(\\mathbb\{R\}^\{n\}\)\.
2. 2\.‖W¯i‖2\+Φ′¯​‖Bi‖2​‖Wψ,i‖2<2​αi​∀i\\\|\\overline\{W\}\_\{i\}\\\|\_\{2\}\+\\overline\{\\Phi^\{\\prime\}\}\\\|B\_\{i\}\\\|\_\{2\}\\\|W\_\{\\psi,i\}\\\|\_\{2\}<2\\sqrt\{\\alpha\_\{i\}\}\\ \\forall i, whereΦ′¯:=max⁡\{\|Φmin′\|,\|Φmax′\|\}\\overline\{\\Phi^\{\\prime\}\}\\mathrel\{\\mathop\{\\mathchar 58\\relax\}\}=\\max\\\{\|\\Phi^\{\\prime\}\_\{\\min\}\|,\|\\Phi^\{\\prime\}\_\{\\max\}\|\\\}\.

Then, for any matrixDm−1t∈𝒟D\_\{m\-1\}^\{t\}\\in\\mathcal\{D\}, the companion matrix of the SpatialRNN gradient dynamics \([14](https://arxiv.org/html/2610.10690#S3.E14)\) has exactly2​ni2n\_\{i\}simple eigenvalues of moduli\|λ\|=αi\|\\lambda\|=\\sqrt\{\\alpha\_\{i\}\}for each blocki=1,…,ki=1,\\dots,k\.

###### Proof\.

See Appendix[A\.1](https://arxiv.org/html/2610.10690#A1.SS1)\. The proof exploits the block\-triangular structure ofW¯\\overline\{W\}andW¯¯\\overline\{\\overline\{W\}\}to factor the characteristic eigenvalue equation into a product ofkkindependent polynomials\. Condition 1 ensures that the eigenvalues ofW¯i\+Bi​Dm−1t​Wψ,i\\overline\{W\}\_\{i\}\+B\_\{i\}D^\{t\}\_\{m\-1\}W\_\{\\psi,i\}are real, while condition 2 ensures that the discriminant of the resulting quadratic equations remains negative, leading to underdamped \(oscillatory\) modes\. Consequently, the eigenvalues form complex conjugate pairs whose magnitude is determined by the constantαi\\sqrt\{\\alpha\_\{i\}\}\. ∎

The block\-triangular structure of the weight matrices enforces a directed information flow among thekkblocks in the PDE dynamics\. Within each block, neurons interact fully, while information flows strictly upstream among blocks\. This structure allows the parameterαi\\alpha\_\{i\}to directly control the gradient dynamics at each blockii\. By choosingαi=1\\alpha\_\{i\}=1, the SpatialRNN achieves marginal stability across all admissible nonlinear perturbations, enabling long\-term gradient propagation and memory retention\. Settingαi<1\\alpha\_\{i\}<1, in turn, leads to exponential decay at the controlled rate, providing an explicit forgetting mechanism at a desired rate\. For a SpatialRNN satisfying the conditions of Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)\(withα1=1\\alpha\_\{1\}=1and a single blockk=1k=1\), Fig\.[2](https://arxiv.org/html/2610.10690#S3.F2)b shows that gradient norms remain stable over long sequences of time\-varying matricesDmt∈𝒟D\_\{m\}^\{t\}\\in\\mathcal\{D\}\. Importantly, the parametersαi\\alpha\_\{i\}can either be fixed a priori or learned via BPTT\.

SpatialRNNs, similarly to SSMs, can improve the training memory efficiency of recurrent models\. In standard RNNs, BPTT demands𝒪⁡\(T\)\\mathcal\{O\}\(T\)memory to cache intermediate hidden states for gradient computation\. For SpatialRNNs, the medium dynamics are invertible, so that intermediate states can be reconstructed backward in time from the final states, reducing this requirement to𝒪⁡\(1\)\\mathcal\{O\}\(1\)in exact arithmetic\. In finite precision, however, the numerical reliability of this reconstruction depends on the dynamical regime\. In the non\-dissipative case \(αi=1,∀i\\alpha\_\{i\}=1,\\,\\forall i\), the medium dynamics preserve a symplectic form\([Marsden and West, 2001](https://arxiv.org/html/2610.10690#bib.bib51)\), enabling backward reconstruction to remain accurate to machine precision\. By contrast, dissipative blocks \(αi<1\\alpha\_\{i\}<1\) exponentially amplify round\-off errors under backward reconstruction and therefore require caching or checkpointing\. Further details on this analysis are presented in Appendix[B](https://arxiv.org/html/2610.10690#A2)\.

## 4Experiments

We evaluate the memory capabilities of the SpatialRNN on benchmark datasets with long temporal dependencies, comparing it against other recurrent models\. Experimental details \(datasets, model architectures, and training procedures\) are included in Appendix[C](https://arxiv.org/html/2610.10690#A3); no regularizers \(e\.g\., dropout or normalization\) have been used, highlighting the intrinsic capabilities of the SpatialRNN without auxiliary optimization\. All experiments use standard BPTT\. Code to reproduce all experiments can be found in[anonymous\.4open\.science/r/SpatialRNN\-code\-748C](https://anonymous.4open.science/r/SpatialRNN-code-748C)\.

Copy task\.As a proof\-of\-concept, we compare SpatialRNNs and HORNNs on the copy task\([Hochreiter and Schmidhuber, 1997](https://arxiv.org/html/2610.10690#bib.bib9);[Arjovsky et al\., 2016](https://arxiv.org/html/2610.10690#bib.bib25)\), a benchmark designed to test the ability to store information over long sequences and accurately retrieve it\. The input sequenceX=\[Xdata,Xdist,Xtrig\]∈\{0,±1\}D×TX=\[X\_\{\\text\{data\}\},X\_\{\\text\{dist\}\},X\_\{\\text\{trig\}\}\]\\in\\\{0,\\pm 1\\\}^\{D\\times T\}consists ofL=5L=5binary tokensXdataX\_\{\\text\{data\}\}with dimensionD=4D=4\(representing the data to be stored\), followed byT−2​LT\-2LdistractorsXdistX\_\{\\text\{dist\}\}\(padded with zeros\), and finally a sequence ofLLtriggersXtrigX\_\{\\text\{trig\}\}\(padded with ones\)\. The models are trained to reconstruct theXdataX\_\{\\text\{data\}\}at the finalLLsteps\. Here, we evaluate the model performance for sequences of increasing lengthT∈\{2p\}p=110T\\in\\\{2^\{p\}\\\}\_\{p=1\}^\{10\}\.

Fig\.[3](https://arxiv.org/html/2610.10690#S4.F3)compares the SpatialRNN against HORNNs \(64 neurons\) with increasing receptive fields, fromd=1d=1\(standard RNN, 5\.5k parameters\) tod=15d=15\(17\.7k parameters\)\. For the SpatialRNN, we impose the conditions of Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)withk=4k=4block matrices of sizesni∈\{32,8,8,8\}n\_\{i\}\\in\\\{32,8,8,8\\\}\(6\.8k parameters\) andα1=1\\alpha\_\{1\}=1to guarantee marginal stability\. Consistent with the numerical

Figure 3:Performance of SpatialRNN and HORNNs on the copy task problem, shown as a function of sequence lengthTT\.gradient analysis in Fig\.[2](https://arxiv.org/html/2610.10690#S3.F2)b, HORNN performance degrades with increasingTT, with larger receptive fieldsdddelaying, but not preventing, this deterioration\. In contrast, the SpatialRNN maintains near\-perfect \(≈99%\{\\approx\}99\\%\) validation accuracy across all sequence lengths, despite having fewer parameters\. This behavior is consistent with the robust marginal stability established in Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)\(withα1=1\\alpha\_\{1\}=1\), which enables long\-term memory retention\. The HORNN comparison also serves as an ablation of state dimension versus PDE\-mediated structure: Add\-order HORNN has a lifted state of dimensionn​dnd, so a second\-order HORNN already matches the order and state dimension of the SpatialRNN \(2​n2n, corresponding toψt,ψt−1\\psi\_\{t\},\\psi\_\{t\-1\}\), while models withd≥5d\\geq 5greatly exceed it\. Nevertheless, all finite\-order HORNNs remain subject to Proposition[3\.2](https://arxiv.org/html/2610.10690#S3.Thmtheorem2)and eventually lose accuracy in the copy task for largeTT, indicating that the observed advantage of the SpatialRNN cannot be attributed simply to higher recurrence order or increased state dimension\. Interestingly, even whenαi\\alpha\_\{i\}is left unconstrained for all blocks, training exhibits a sharp transition in performance when the model learns a mix of blocks withαi<1\\alpha\_\{i\}<1andαj≈1\\alpha\_\{j\}\\approx 1, allowing it to simultaneously retain information and forget distractors \(see Appendix[A\.4](https://arxiv.org/html/2610.10690#A1.SS4)and Fig\.[5](https://arxiv.org/html/2610.10690#A1.F5)\)\.

Permuted sequential MNIST \(psMNIST\)\.We evaluate long\-range dependency modeling on the psMNIST classification task\([Le et al\., 2015](https://arxiv.org/html/2610.10690#bib.bib39);[Arjovsky et al\., 2016](https://arxiv.org/html/2610.10690#bib.bib25)\)\. Here, MNIST images are first flattened into sequences of lengthT=784T=784with feature dimensionD=1D=1, such that pixels are fed to the model one at a time\. The sequence is then normalized and randomly permuted, destroying the local spatial structure of the images and forcing the model to rely entirely on long\-term temporal memory for classification\.

Table 1:psMNIST performance\. Mean ± std over 3 runs when available; otherwise, maximum accuracy\. Only single\-layer architectures are accounted for\.We evaluate two SpatialRNN configurations: an 88\-neuron model \(block sizes\{64,8,8,8\}\\\{64,8,8,8\\\}, 12\.8k parameters\) and a 152\-neuron model \(block sizes\{128,8,8,8\}\\\{128,8,8,8\\\}, 36\.5k parameters\)\. In both cases, we enforceα1=1\\alpha\_\{1\}=1on the larger block to guarantee marginal stability\. The SpatialRNNs achieve the highest accuracy among the single\-layer recurrent baselines,96\.3±0\.196\.3\\pm 0\.1% \(Table[1](https://arxiv.org/html/2610.10690#S4.T1)\), outperforming standard gated architectures such as LSTM and GRU, as well as recent oscillatory models including NWM & coRNN\. Notably, this performance is achieved with roughly four to ten times fewer parameters\.

Noise padded CIFAR\-10 \(npCIFAR10\)\.In this benchmark\([Chang et al\., 2019](https://arxiv.org/html/2610.10690#bib.bib38)\), CIFAR\-10 images \(of size32×32×332\\times 32\\times 3\) are normalized and reshaped into sequences of lengthL=32L=32with feature dimensionD=96D=96, so that each timestep contains one image row across all color channels\. Standard normal noise is then appended to reach a total sequence length ofT=1000T=1000, acting as distractors\.

We evaluate two SpatialRNN configurations: a 88\-neuron model \(block sizes\{64,8,8,8\}\\\{64,8,8,8\\\}, 21\.2k parameters\) and an 152\-neuron model \(block sizes\{128,8,8,8\}\\\{128,8,8,8\\\}, 50\.9k parameters\), where we once again enforceα1=1\\alpha\_\{1\}=1for marginal stability\. The 88\-neuron SpatialRNN achieves competitive accuracy while operating with a fraction of the parameters required by standard gated and continuous\-time baselines \(Table[2](https://arxiv.org/html/2610.10690#S4.T2)\)\. Most importantly, Fig\.[4](https://arxiv.org/html/2610.10690#S4.F4)shows that performance is invariant to the amount of appended noise, indicating that accuracy is mostly dictated by the spatial feature extraction of the initial sequence rather than by temporal degradation over the sequence lengthTT\.

Table 2:npCIFAR10 performance forT=1000T=1000\. Maximum test accuracy over 3 trials\.Figure 4:SpatialRNN validation accuracy on the npCIFAR10 task for varyingTT\.
## 5Conclusion

In this work, we introduced a spatial computing model in which neuron states communicate through a dynamically evolving spatial medium rather than through direct recurrent connections\. We showed that the medium state implicitly encodes the entire history of past hidden states, giving the model an unbounded receptive field with a fixed number of parameters\. The second\-order dynamics induced by the medium also provide a structural mechanism for controlling gradient propagation: under suitable constraints, the gradient modes can be confined to circles of prescribed radii despite input\-dependent nonlinear perturbations, enabling long\-term memory retention\. Together, these results establish spatial dynamics as a parameter\-efficient mechanism for combining long\-range temporal dependencies with robust gradient propagation\.

From a theoretical perspective, our analysis identifies a structural principle underlying long\-term memory in recurrent architectures\. We show that robustness to input\-dependent nonlinear perturbations is governed not simply by the spectrum or parameterization of the recurrence matrices, but by how the nonlinearity enters the gradient dynamics \(Proposition[3\.2](https://arxiv.org/html/2610.10690#S3.Thmtheorem2)and Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)\)\. In finite\-order RNNs, the nonlinear term directly perturbs the coefficients governing modal decay, precluding robust marginal stability\. In the SpatialRNN, the second\-order recurrence separates the input\-dependent perturbation from the coefficients that determine the magnitudes of the gradient modes, allowing those magnitudes to remain invariant\. This distinction provides a general design principle for recurrent architectures capable of sustaining long\-lived gradient modes\.

Related work\.Neural ordinary differential equations\([Chen et al\., 2018](https://arxiv.org/html/2610.10690#bib.bib21)\)generalize RNNs by extending discrete\-time recurrences to continuous\-time dynamics, enabling richer state evolution and improved control over gradient behavior\. A closely related line of work considers oscillatory\([Rusch and Mishra, 2021a](https://arxiv.org/html/2610.10690#bib.bib3);[Rusch et al\., 2022](https://arxiv.org/html/2610.10690#bib.bib13)\)and second\-order\([Effenberger et al\., 2025](https://arxiv.org/html/2610.10690#bib.bib36);[Rusch and Mishra, 2021b](https://arxiv.org/html/2610.10690#bib.bib43)\)dynamical systems, in which neurons evolve as coupled oscillators\. When constrained to local interactions, these models give rise to wave\-like dynamics, as in neural wave machines\([Keller and Welling, 2023](https://arxiv.org/html/2610.10690#bib.bib14)\), which have demonstrated strong performance on spatially structured data\([Jacobs et al\., 2025](https://arxiv.org/html/2610.10690#bib.bib15)\)and long\-sequence tasks\([Keller and Welling, 2023](https://arxiv.org/html/2610.10690#bib.bib14)\)\. These models often exploit traveling waves to encode and propagate temporal information\([Benigno et al\., 2023](https://arxiv.org/html/2610.10690#bib.bib20);[Jacobs et al\., 2025](https://arxiv.org/html/2610.10690#bib.bib15);[Shervani\-Tabar et al\., 2026](https://arxiv.org/html/2610.10690#bib.bib46)\)\. However, their stability depends directly on discretization parameters, and thus the gradient dynamics exhibit residual exponential scaling with the time stepΔ​t\\Delta t\([Rusch et al\., 2022](https://arxiv.org/html/2610.10690#bib.bib13)\)\. Related ideas have also been explored in other state\-based architectures, including reservoir computing\([Benigno et al\., 2023](https://arxiv.org/html/2610.10690#bib.bib20)\)and SSMs\([Rusch and Rus, 2025](https://arxiv.org/html/2610.10690#bib.bib22)\)\. Our approach is closely related to these frameworks but differs fundamentally since traveling waves are encoded in a shared medium, rather than emerging from interactions between local oscillatory units; these approaches coincide only under certain restrictive conditions \(Appendix[A\.5](https://arxiv.org/html/2610.10690#A1.SS5)\)\.

Limitation\.A key limitation of our analysis is the gap between the frozen\-time and the full LTV gradient dynamics\. The frozen\-time guarantees of Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)are necessary, but not sufficient, for marginal stability under arbitrary switching, since products of matrices with unit\-modulus spectra may still grow or decay\. Gradients may thus grow over long horizons, although we observe empirically that this drift is substantially attenuated and does not appear exponential\. In particular, the exact Jacobians of a trained SpatialRNN change by only around6%6\\%over783783steps along psMNIST trajectories \(see Appendix[A\.7](https://arxiv.org/html/2610.10690#A1.SS7)for further numerical analysis\)\. Asymptotic stability of the SpatialRNN gradients can be enforced for the full LTV system under additional structural constraints \(Theorem[A\.3](https://arxiv.org/html/2610.10690#A1.Thmtheorem3)\), yet establishing necessary and sufficient conditions for LTV marginal stability remains an open problem\.

More broadly, our framework provides an explicit mechanism through which spatiotemporal dynamics endow recurrent models with long\-lived memory—a longstanding goal in the AI community\([Keller et al\., 2026](https://arxiv.org/html/2610.10690#bib.bib47)\)\. However, memory alone does not determine overall performance\. The information capacity of the medium is effectively bounded by its dimension and numerical precision \(Remark[2\.4](https://arxiv.org/html/2610.10690#S2.Thmtheorem4)\), while feature extraction and expressivity remain architecture\-dependent\. This result is evident in the npCIFAR10 benchmark \(Fig\.[4](https://arxiv.org/html/2610.10690#S4.F4)\), where the SpatialRNN robustly maintains performance across varying sequence lengths yet remains limited in classification accuracy\. Computational efficiency presents another trade\-off: Each SpatialRNN step involves a few more matrix\-vector products than standard RNNs, and its nonlinear recurrence is evaluated sequentially \(unlike the associative scans in linear SSMs\)\.

Future work\.A natural direction for future work is thus to integrate spatial computing mechanisms into more expressive state\-based architectures such as SSMs, which can scale efficiently on modern hardware\([Gonzalez et al\., 2026](https://arxiv.org/html/2610.10690#bib.bib42)\); for instance, by lifting their recurrence equations to second order\. We expect this paradigm to be best suited to streaming problems requiring a fixed\-size state, long retention, and selective forgetting\. Relevant applications include sensor monitoring, long\-horizon control, and spatiotemporal signal prediction, where the locality and geometry of the medium can provide useful inductive biases\. We also expect this formulation to find applications in physical computing, where wave\-based systems, such as those exhibiting solitons, serve as computational substrates\([Hughes et al\., 2019](https://arxiv.org/html/2610.10690#bib.bib45);[Marcucci et al\., 2020](https://arxiv.org/html/2610.10690#bib.bib44);[Montanari et al\., 2026](https://arxiv.org/html/2610.10690#bib.bib40)\)\. In this context, the symplectic structure induced by non\-dissipative PDE dynamics can provide an additional advantage\([Chen et al\., 2020](https://arxiv.org/html/2610.10690#bib.bib41)\), enabling exact backward state reconstruction and reducing the memory cost of backpropagation toO⁡\(1\)O\(1\)with respect to the sequence lengthTT\.

## Broader Impact Statement

This work studies how a spatially evolving medium, governed by discretized PDEs, can mediate memory retention and gradient propagation in recurrent models\. Its contributions are theoretical results and small\-scale benchmarks, and we do not foresee direct negative societal impacts beyond those common to general\-purpose sequence models\. By retaining long\-range dependencies in a fixed\-size state, rather than through gating mechanisms or a key\-value memory that grows with sequence length, the SpatialRNN may help reduce the memory cost of long\-horizon inference\. We also expect the formulation to be relevant to analog and wave\-based physical computing substrates, where the medium dynamics could be realized directly in hardware\. Whether these advantages translate into energy savings in practice remains to be established\.

#### Acknowledgments

This work was supported by the ARO Grant No\. W911NF\-24\-1\-0228\.

## References

- Arjovskyet al\.\(2016\)M\. Arjovsky, A\. Shah, and Y\. BengioUnitary evolution recurrent neural networks\.InInternational Conference on Machine Learning,pp\. 1120–1128\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1),[Table 1](https://arxiv.org/html/2610.10690#S4.T1.2.2.1),[§4](https://arxiv.org/html/2610.10690#S4.p2.1),[§4](https://arxiv.org/html/2610.10690#S4.p5.1)\.
- Bengioet al\.\(1994\)Y\. Bengio, P\. Simard, and P\. FrasconiLearning long\-term dependencies with gradient descent is difficult\.IEEE Transactions on Neural Networks5\(2\),pp\. 157–166\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1),[§3](https://arxiv.org/html/2610.10690#S3.p1.4)\.
- Benignoet al\.\(2023\)G\. B\. Benigno, R\. C\. Budzinski, Z\. W\. Davis, J\. H\. Reynolds, and L\. MullerWaves traveling over a map of visual space can ignite short\-term predictions of sensory input\.Nature Communications14,pp\. 3409\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p3.1),[§5](https://arxiv.org/html/2610.10690#S5.p3.1)\.
- Ceni \(2025\)A\. CeniRandom orthogonal additive filters: a solution to the vanishing/exploding gradient of deep neural networks\.IEEE Transactions on Neural Networks and Learning Systems36\(6\),pp\. 10794–10807\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1)\.
- Changet al\.\(2019\)B\. Chang, M\. Chen, E\. Haber, and E\. H\. ChiAntisymmetricRNN: a dynamical system view on recurrent neural networks\.InInternational Conference on Learning Representations,Cited by:[§4](https://arxiv.org/html/2610.10690#S4.p7.1)\.
- Changet al\.\(2017\)S\. Chang, Y\. Zhang, W\. Han, M\. Yu, X\. Guo, W\. Tan, X\. Cui, M\. Witbrock, M\. A\. Hasegawa\-Johnson, and T\. S\. HuangDilated recurrent neural networks\.InAdvances in Neural Information Processing Systems,Vol\.30\.Cited by:[Table 1](https://arxiv.org/html/2610.10690#S4.T1.2.4.1)\.
- Chenet al\.\(2018\)R\. T\. Chen, Y\. Rubanova, J\. Bettencourt, and D\. K\. DuvenaudNeural ordinary differential equations\.InAdvances in Neural Information Processing Systems,Vol\.31\.Cited by:[§5](https://arxiv.org/html/2610.10690#S5.p3.1)\.
- Chenet al\.\(2020\)Z\. Chen, J\. Zhang, M\. Arjovsky, and L\. BottouSymplectic recurrent neural networks\.InInternational Conference on Learning Representations,Cited by:[§5](https://arxiv.org/html/2610.10690#S5.p6.1)\.
- Choet al\.\(2014\)K\. Cho, B\. Van Merriënboer, Ç\. Gulçehre, D\. Bahdanau, F\. Bougares, H\. Schwenk, and Y\. BengioLearning phrase representations using RNN encoder–decoder for statistical machine translation\.InConference on Empirical Methods in Natural Language Processing,pp\. 1724–1734\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1)\.
- Chunget al\.\(2014\)J\. Chung, C\. Gulcehre, K\. Cho, and Y\. BengioEmpirical evaluation of gated recurrent neural networks on sequence modeling\.arXiv:1412\.3555\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1)\.
- Daoet al\.\(2022\)T\. Dao, D\. Fu, S\. Ermon, A\. Rudra, and C\. RéFlashattention: fast and memory\-efficient exact attention with io\-awareness\.InAdvances in Neural Information Processing Systems,Vol\.35\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1)\.
- Daviset al\.\(2021\)Z\. W\. Davis, G\. B\. Benigno, C\. Fletterman, T\. Desbordes, C\. Steward, T\. J\. Sejnowski, J\. H\. Reynolds, and L\. MullerSpontaneous traveling waves naturally emerge from horizontal fiber time delays and travel through locally asynchronous\-irregular states\.Nature Communications12,pp\. 6057\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p3.1)\.
- Effenbergeret al\.\(2025\)F\. Effenberger, P\. Carvalho, I\. Dubinin, and W\. SingerThe functional role of oscillatory dynamics in neocortical circuits: a computational perspective\.Proceedings of the National Academy of Sciences of the U\.S\.A\.122\(4\),pp\. e2412830122\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1),[§5](https://arxiv.org/html/2610.10690#S5.p3.1)\.
- Ermentrout and Kleinfeld \(2001\)G\. B\. Ermentrout and D\. KleinfeldTraveling electrical waves in cortex: insights from phase dynamics and speculation on a computational role\.Neuron29\(1\),pp\. 33–44\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p3.1)\.
- Gonzalezet al\.\(2026\)X\. Gonzalez, L\. Kozachkov, D\. M\. Zoltowski, K\. L\. Clarkson, and S\. LindermanPredictability enables parallelization of nonlinear state space models\.InAdvances in Neural Information Processing Systems,Vol\.39\.Cited by:[§5](https://arxiv.org/html/2610.10690#S5.p6.1)\.
- Gu and Dao \(2024\)A\. Gu and T\. DaoMamba: linear\-time sequence modeling with selective state spaces\.InFirst Conference on Language Modeling,Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1)\.
- Guet al\.\(2022\)A\. Gu, K\. Goel, and C\. ReEfficiently modeling long sequences with structured state spaces\.InInternational Conference on Learning Representations,Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1)\.
- Heet al\.\(2016\)K\. He, X\. Zhang, S\. Ren, and J\. SunDeep residual learning for image recognition\.InIEEE Conference on Computer Vision and Pattern Recognition,pp\. 770–778\.Cited by:[§C\.4](https://arxiv.org/html/2610.10690#A3.SS4.p4.1)\.
- Helfrichet al\.\(2018\)K\. Helfrich, D\. Willmott, and Q\. YeOrthogonal recurrent neural networks with scaled Cayley transform\.InInternational Conference on Machine Learning,pp\. 1969–1978\.Cited by:[§3\.1](https://arxiv.org/html/2610.10690#S3.SS1.p6.1),[Table 1](https://arxiv.org/html/2610.10690#S4.T1.2.3.1)\.
- Henaffet al\.\(2016\)M\. Henaff, A\. Szlam, and Y\. LeCunRecurrent orthogonal networks and long\-memory tasks\.InInternational Conference on Machine Learning,pp\. 2034–2042\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1),[§3\.1](https://arxiv.org/html/2610.10690#S3.SS1.p6.1)\.
- Hochreiter and Schmidhuber \(1997\)S\. Hochreiter and J\. SchmidhuberLong short\-term memory\.Neural Computation9\(8\),pp\. 1735–1780\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1),[§4](https://arxiv.org/html/2610.10690#S4.p2.1)\.
- Hugheset al\.\(2019\)T\. W\. Hughes, I\. A\. Williamson, M\. Minkov, and S\. FanWave physics as an analog recurrent neural network\.Science Advances5\(12\),pp\. eaay6946\.Cited by:[§5](https://arxiv.org/html/2610.10690#S5.p6.1)\.
- Ioffe and Szegedy \(2015\)S\. Ioffe and C\. SzegedyBatch normalization: accelerating deep network training by reducing internal covariate shift\.InInternational Conference on Machine Learning,pp\. 448–456\.Cited by:[§C\.4](https://arxiv.org/html/2610.10690#A3.SS4.p4.1)\.
- Jacobset al\.\(2025\)M\. Jacobs, R\. C\. Budzinski, L\. Muller, D\. E\. Ba, and T\. A\. KellerTraveling waves integrate spatial information through time\.In8th Annual Conference on Cognitive Computational Neuroscience,Cited by:[§5](https://arxiv.org/html/2610.10690#S5.p3.1)\.
- Kaget al\.\(2020\)A\. Kag, Z\. Zhang, and V\. SaligramaRNNs incrementally evolving on an equilibrium manifold: a panacea for vanishing and exploding gradients?\.InInternational Conference on Learning Representations,Cited by:[Table 2](https://arxiv.org/html/2610.10690#S4.T2.2.2.1)\.
- Karuvallyet al\.\(2026\)A\. Karuvally, F\. Nowak, T\. A\. Keller, C\. A\. Alonso, T\. Sejnowski, and H\. T\. SiegelmannBridging expressivity and scalability with adaptive unitary SSMs\.InAdvances in Neural Information Processing Systems,Vol\.39\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1)\.
- Kelleret al\.\(2026\)T\. A\. Keller, L\. Muller, T\. J\. Sejnowski, and M\. WellingA spatiotemporal perspective on dynamical computation in neural information processing systems\.arXiv:2409\.13669v2\.Cited by:[§5](https://arxiv.org/html/2610.10690#S5.p5.1)\.
- Kelleret al\.\(2024\)T\. A\. Keller, L\. Muller, T\. Sejnowski, and M\. WellingTraveling waves encode the recent past and enhance sequence learning\.InInternational Conference on Learning Representations,Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1),[Table 2](https://arxiv.org/html/2610.10690#S4.T2.2.10.1),[Table 2](https://arxiv.org/html/2610.10690#S4.T2.2.3.1),[Table 2](https://arxiv.org/html/2610.10690#S4.T2.2.4.1),[Table 2](https://arxiv.org/html/2610.10690#S4.T2.2.5.1),[Table 2](https://arxiv.org/html/2610.10690#S4.T2.2.7.1)\.
- Keller and Welling \(2023\)T\. A\. Keller and M\. WellingNeural wave machines: learning spatiotemporally structured representations with locally coupled oscillatory recurrent neural networks\.InInternational Conference on Machine Learning,pp\. 16168–16189\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1),[Table 1](https://arxiv.org/html/2610.10690#S4.T1.2.5.1),[Table 1](https://arxiv.org/html/2610.10690#S4.T1.2.6.1),[§5](https://arxiv.org/html/2610.10690#S5.p3.1)\.
- Kerget al\.\(2019\)G\. Kerg, K\. Goyette, M\. Puelma Touzel, G\. Gidel, E\. Vorontsov, Y\. Bengio, and G\. LajoieNon\-normal recurrent neural network \(NNRNN\): learning long time dependencies while improving expressivity with transient dynamics\.InAdvances in Neural Information Processing Systems,Vol\.32\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1),[§3\.1](https://arxiv.org/html/2610.10690#S3.SS1.p6.1)\.
- Kwonet al\.\(2023\)W\. Kwon, Z\. Li, S\. Zhuang, Y\. Sheng, L\. Zheng, C\. H\. Yu, J\. Gonzalez, H\. Zhang, and I\. StoicaEfficient memory management for large language model serving with pagedattention\.InProceedings of the 29th Symposium on Operating Systems Principles,pp\. 611–626\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1)\.
- Leet al\.\(2015\)Q\. V\. Le, N\. Jaitly, and G\. E\. HintonA simple way to initialize recurrent networks of rectified linear units\.arXiv:1504\.00941\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1),[§4](https://arxiv.org/html/2610.10690#S4.p5.1)\.
- Liet al\.\(2018\)S\. Li, W\. Li, C\. Cook, C\. Zhu, and Y\. GaoIndependently recurrent neural network \(indRNN\): building a longer and deeper rnn\.InIEEE Conference on Computer Vision and Pattern Recognition,pp\. 5457–5466\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1),[§3\.1](https://arxiv.org/html/2610.10690#S3.SS1.p6.1)\.
- Liberzon \(2003\)D\. LiberzonSwitching in systems and control\.Vol\.190,Springer\.Cited by:[§A\.6](https://arxiv.org/html/2610.10690#A1.SS6.p4.1.1),[§3\.1](https://arxiv.org/html/2610.10690#S3.SS1.p2.2)\.
- Marcucciet al\.\(2020\)G\. Marcucci, D\. Pierangeli, and C\. ContiTheory of neuromorphic computing by waves: machine learning by rogue waves, dispersive shocks, and solitons\.Physical Review Letters125\(9\),pp\. 093901\.Cited by:[§5](https://arxiv.org/html/2610.10690#S5.p6.1)\.
- Marsden and West \(2001\)J\. E\. Marsden and M\. WestDiscrete mechanics and variational integrators\.Acta Numerica10,pp\. 357–514\.Cited by:[Appendix B](https://arxiv.org/html/2610.10690#A2.p1.1),[§3\.2](https://arxiv.org/html/2610.10690#S3.SS2.p4.1)\.
- Montanariet al\.\(2026\)A\. N\. Montanari, F\. Bullo, D\. Krotov, and A\. E\. MotterEnergy\-based dynamical models for neurocomputation, learning, and optimization\.arXiv:2604\.05042\.Cited by:[§5](https://arxiv.org/html/2610.10690#S5.p6.1)\.
- Neilet al\.\(2016\)D\. Neil, M\. Pfeiffer, and S\. LiuPhased lstm: accelerating recurrent network training for long or event\-based sequences\.InAdvances in Neural Information Processing Systems,Vol\.29\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1)\.
- Ogataet al\.\(1995\)K\. Ogataet al\.Discrete\-time control systems\.Vol\.2,Prentice Hall Englewood Cliffs\.Cited by:[§A\.1](https://arxiv.org/html/2610.10690#A1.SS1.p3.3.1),[§2\.2](https://arxiv.org/html/2610.10690#S2.SS2.p2.1.1)\.
- Ribeiroet al\.\(2020\)A\. H\. Ribeiro, K\. Tiels, L\. A\. Aguirre, and T\. SchönBeyond exploding and vanishing gradients: analysing RNN training using attractors and smoothness\.InInternational Conference on Artificial Intelligence and Statistics,pp\. 2370–2380\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p1.1)\.
- Ruschet al\.\(2022\)T\. K\. Rusch, B\. Chamberlain, J\. Rowbottom, S\. Mishra, and M\. BronsteinGraph\-coupled oscillator networks\.InInternational Conference on Machine Learning,pp\. 18888–18909\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1),[§5](https://arxiv.org/html/2610.10690#S5.p3.1)\.
- Rusch and Mishra \(2021a\)T\. K\. Rusch and S\. MishraCoupled oscillatory recurrent neural network \(coRNN\): an accurate and \(gradient\) stable architecture for learning long time dependencies\.InInternational Conference on Learning Representations,Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1),[Table 2](https://arxiv.org/html/2610.10690#S4.T2.2.9.1),[§5](https://arxiv.org/html/2610.10690#S5.p3.1)\.
- Rusch and Mishra \(2021b\)T\. K\. Rusch and S\. MishraUnicoRNN: a recurrent model for learning very long time dependencies\.InInternational Conference on Machine Learning,pp\. 9168–9178\.Cited by:[§5](https://arxiv.org/html/2610.10690#S5.p3.1)\.
- Rusch and Rus \(2025\)T\. K\. Rusch and D\. RusOscillatory state\-space models\.InInternational Conference on Learning Representations,Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1),[§5](https://arxiv.org/html/2610.10690#S5.p3.1)\.
- Shervani\-Tabaret al\.\(2026\)N\. Shervani\-Tabar, S\. L\. Brincat, M\. Lundqvist, and E\. K\. MillerEmergent traveling waves in neural circuits\.bioRxiv\.Cited by:[§5](https://arxiv.org/html/2610.10690#S5.p3.1)\.
- Sieberet al\.\(2024\)J\. Sieber, C\. A\. Alonso, A\. Didier, M\. N\. Zeilinger, and A\. OrvietoUnderstanding the differences in foundation models: attention, state space models, and recurrent neural networks\.InAdvances in Neural Information Processing Systems,Vol\.37\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1)\.
- Singer \(2021\)W\. SingerRecurrent dynamics in the cerebral cortex: integration of sensory evidence with stored knowledge\.Proceedings of the National Academy of Sciences of the U\.S\.A\.118\(33\),pp\. e2101043118\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p3.1)\.
- Soltani and Jiang \(2016\)R\. Soltani and H\. JiangHigher order recurrent neural networks\.arXiv:1605\.00064\.Cited by:[§2\.1](https://arxiv.org/html/2610.10690#S2.SS1.p1.1)\.
- Somvanshiet al\.\(2025\)S\. Somvanshi, M\. M\. Islam, M\. S\. Mimi, S\. B\. B\. Polock, G\. Chhetri, and S\. DasFrom S4 to Mamba: a comprehensive survey on structured state space models\.arXiv:2503\.18970\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1)\.
- Srivastavaet al\.\(2014\)N\. Srivastava, G\. Hinton, A\. Krizhevsky, I\. Sutskever, and R\. SalakhutdinovDropout: a simple way to prevent neural networks from overfitting\.Journal of Machine Learning Research15\(56\),pp\. 1929–1958\.Cited by:[§C\.4](https://arxiv.org/html/2610.10690#A3.SS4.p4.1)\.
- Sun \(2008\)Z\. SunA note on marginal stability of switched systems\.IEEE Transactions on Automatic Control53\(2\),pp\. 625–631\.Cited by:[§3](https://arxiv.org/html/2610.10690#S3.p2.1)\.
- Vaswaniet al\.\(2017\)A\. Vaswani, N\. Shazeer, N\. Parmar, J\. Uszkoreit, L\. Jones, A\. N\. Gomez, Ł\. Kaiser, and I\. PolosukhinAttention is all you need\.InAdvances in Neural Information Processing Systems,Vol\.30\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1)\.
- Wanget al\.\(2025\)F\. Wang, Z\. Zhang, X\. Zhang, Z\. Wu, T\. Mo, Q\. Lu, W\. Wang, R\. Li, J\. Xu, X\. Tang,et al\.A comprehensive survey of small language models in the era of large language models: techniques, enhancements, applications, collaboration with LLMs, and trustworthiness\.ACM Transactions on Intelligent Systems and Technology16\(6\),pp\. 1–87\.Cited by:[§1](https://arxiv.org/html/2610.10690#S1.p2.1)\.

## Appendix ASupplementary Theoretical Analysis

### A\.1Proofs of Main Theoretical Results

We provide the proofs for the main theoretical statements presented in the main text\.

###### Proof of Theorem[2\.3](https://arxiv.org/html/2610.10690#S2.Thmtheorem3)\.

Consider the medium stateψ\\psigiven by Eq\. \([6](https://arxiv.org/html/2610.10690#S2.E6)\) and substitute it into both sides of the LDE \([2](https://arxiv.org/html/2610.10690#S2.E2)\):

∑k=1t−1Wk​ht−k=W¯​∑k=1t−2Wk​ht−1−k\+W¯¯​∑k=1t−3Wk​ht−2−k\+B​ht−1\\displaystyle\\sum\_\{k=1\}^\{t\-1\}W\_\{k\}h\_\{t\-k\}=\\overline\{W\}\\sum\_\{k=1\}^\{t\-2\}W\_\{k\}h\_\{t\-1\-k\}\+\\overline\{\\overline\{W\}\}\\sum\_\{k=1\}^\{t\-3\}W\_\{k\}h\_\{t\-2\-k\}\+Bh\_\{t\-1\}=W¯​∑k=2t−1Wk−1​ht−k\+W¯¯​∑k=3t−1Wk−2​ht−k\+B​ht−1\.\\displaystyle=\\overline\{W\}\\sum\_\{k=2\}^\{t\-1\}W\_\{k\-1\}h\_\{t\-k\}\+\\overline\{\\overline\{W\}\}\\sum\_\{k=3\}^\{t\-1\}W\_\{k\-2\}h\_\{t\-k\}\+Bh\_\{t\-1\}\.Equating the coefficients in both sides for eachht−kh\_\{t\-k\}, we obtainW1=BW\_\{1\}=Bfork=1k=1,W2=W¯​W1W\_\{2\}=\\overline\{W\}W\_\{1\}fork=2k=2, andWk=W¯​Wk−1\+W¯¯​Wk−2W\_\{k\}=\\overline\{W\}W\_\{k\-1\}\+\\overline\{\\overline\{W\}\}W\_\{k\-2\}fork≥3k\\geq 3\. ∎

###### Proof of Proposition[3\.1](https://arxiv.org/html/2610.10690#S3.Thmtheorem1)\.

Consider the frozen\-time recurrence systemGmt=D​∑kWk​Gm−ktG\_\{m\}^\{t\}=D\\sum\_\{k\}W\_\{k\}\\,G\_\{m\-k\}^\{t\}, form≥1m\\geq 1, andG0t=𝕀G\_\{0\}^\{t\}=\\mathbb\{I\}\. Define thezz\-transformG^​\(z\):=∑m=0∞Gmt​z−m\\hat\{G\}\(z\)\\mathrel\{\\mathop\{\\mathchar 58\\relax\}\}=\\sum\_\{m=0\}^\{\\infty\}G^\{t\}\_\{m\}z^\{\-m\}\. Then, it follows from the recurrence system that

G^​\(z\)\\displaystyle\\hat\{G\}\(z\)=G0\+∑m=1∞D⁡\(∑k=1mWk​Gm−kt\)​z−m\\displaystyle=G\_\{0\}\+\\sum\_\{m=1\}^\{\\infty\}D\\left\(\\sum\_\{k=1\}^\{m\}W\_\{k\}G^\{t\}\_\{m\-k\}\\right\)z^\{\-m\}\(16\)=𝕀\+D​∑k=1∞Wk​z−k​∑m=k∞Gm−kt​zm−k\\displaystyle=\\mathbb\{I\}\+D\\sum\_\{k=1\}^\{\\infty\}W\_\{k\}z^\{\-k\}\\sum\_\{m=k\}^\{\\infty\}G^\{t\}\_\{m\-k\}z^\{m\-k\}\(17\)=𝕀\+D⁡\(∑k=1∞Wk​z−k\)​\(∑j=0∞Gjt​z−j\)\\displaystyle=\\mathbb\{I\}\+D\\left\(\\sum\_\{k=1\}^\{\\infty\}W\_\{k\}z^\{\-k\}\\right\)\\left\(\\sum\_\{j=0\}^\{\\infty\}G^\{t\}\_\{j\}z^\{\-j\}\\right\)\(18\)=𝕀\+D​W​\(z\)​G​\(z\)\.\\displaystyle=\\mathbb\{I\}\+DW\(z\)G\(z\)\.\(19\)Re\-arranging to make the function explicit, we get

G⁡\(z\)=\(𝕀−D​W^​\(z\)\)−1\.G\(z\)=\(\\mathbb\{I\}\-D\\hat\{W\}\(z\)\)^\{\-1\}\.\(20\)It is well\-known that the poles ofG⁡\(z\)G\(z\)are precisely the valuesz∈ℂz\\in\\mathbb\{C\}that satisfy the characteristic equationdet⁡\(𝕀−D​W^​\(z\)\)=0\\det\(\\mathbb I \- D \\hat W\(z\)\)=0\([Ogata and others, 1995](https://arxiv.org/html/2610.10690#bib.bib50)\)\. These values characterize the modes of the linear recurrence system, and therefore correspond to the eigenvaluesλ\\lambdagoverning the asymptotic behavior of the gradient dynamics\. ∎

###### Proof of Proposition[3\.2](https://arxiv.org/html/2610.10690#S3.Thmtheorem2)\.

We proceed by contradiction\. Suppose the system is robustly marginally stable, i\.e\., all rootsλ\\lambdaof the characteristic equation \([12](https://arxiv.org/html/2610.10690#S3.E12)\) satisfy\|λ\|=1\|\\lambda\|=1for everyD∈𝒟D\\in\\mathcal\{D\}\.

Consider the subsetΩ=\{D∈𝒟∣D=γ𝕀,γ∈ℝ≠0\}\\Omega=\\\{D\\in\\mathcal\{D\}\\mid D=\\gamma\\mathbb\{I\},\\,\\gamma\\,\\in\\,\\mathbb\{R\}\_\{\\neq 0\}\\\}and the auxiliary domainΩ¯=\{D−1∣D∈Ω\}\\overline\{\\Omega\}=\\\{D^\{\-1\}\\mid D\\in\\Omega\\\}\. SinceΦ′\\Phi^\{\\prime\}takes values arbitrarily close to 0, it follows thatΩ¯\\overline\{\\Omega\}contains scalar matrices of the formΔ=δ​𝕀\\Delta=\\delta\\mathbb\{I\}for arbitrarily large \(in norm\)\|δ\|\>0\|\\delta\|\>0\. Note that all matricesD∈ΩD\\in\\Omegaare invertible\. TakingΔ=D−1\\Delta=D^\{\-1\}we can re\-write the characteristic equation \([12](https://arxiv.org/html/2610.10690#S3.E12)\) for such matrices

det⁡\(I−D​W^​\(λ\)\)=det⁡\(D​\(Δ−W^​\(λ\)\)\)=det⁡\(D\)​det⁡\(Δ−W^​\(λ\)\)=0\.\\det\(I \- D\\hat\{W\}\(\\lambda\)\)=\\det\(D\(\\Delta\- \\hat\{W\}\(\\lambda\)\)\)=\\det\(D\)\\det\(\\Delta\- \\hat\{W\}\(\\lambda\)\)=0\.\(21\)SinceDDis nonsingular, this equation is equivalent todet⁡\(Δ−W^​\(λ\)\)=0\\det\(\\Delta\- \\hat\{W\}\(\\lambda\)\)=0\.

Recall thatW^​\(λ\)=∑k=1dWk​λ−k\\hat\{W\}\(\\lambda\)=\\sum\_\{k=1\}^\{d\}W\_\{k\}\\lambda^\{\-k\}\. Now, consider the unit circle𝕊1=\{z∈ℂ:\|z\|=1\}\\mathbb\{S\}^\{1\}=\\\{z\\in\\mathbb\{C\}\\mathrel\{\\mathop\{\\mathchar 58\\relax\}\}\|z\|=1\\\}, and let us defineρW^​\(λ\)=maxi⁡\|μi​\(W^​\(λ\)\)\|\\rho\_\{\\hat\{W\}\}\(\\lambda\)=\\max\_\{i\}\|\\mu\_\{i\}\(\\hat\{W\}\(\\lambda\)\)\|as the spectral radius of the polynomial matrixW^​\(λ\)\\hat\{W\}\(\\lambda\), whereμi\\mu\_\{i\}denotes its eigenvalues\. SinceW^​\(λ\)\\hat\{W\}\(\\lambda\)is a finite polynomial inλ−1\\lambda^\{\-1\}, it is continuous on𝕊1\\mathbb\{S\}^\{1\}\. Because𝕊1\\mathbb\{S\}^\{1\}is compact,ρW^​\(λ\)\\rho\_\{\\hat\{W\}\}\(\\lambda\)achieves a maximum valueM:=supλ∈𝕊1ρW^​\(λ\)<∞M\\mathrel\{\\mathop\{\\mathchar 58\\relax\}\}=\\sup\_\{\\lambda\\in\\mathbb\{S\}^\{1\}\}\\rho\_\{\\hat\{W\}\}\(\\lambda\)<\\infty\. We can thus choose a matrixΔ∗=δ​𝕀∈Ω\\Delta^\{\*\}=\\delta\\mathbb\{I\}\\in\\Omegasuch thatδ\>M\\delta\>M, and takeD∗=\(Δ∗\)−1∈𝒟D^\{\*\}=\(\\Delta^\{\*\}\)^\{\-1\}\\in\\mathcal\{D\}\. The characteristic equation becomesdet⁡\(δ​𝕀−W^​\(λ\)\)=0\\det\(\\delta\\mathbb\{I\} \- \\hat\{W\}\(\\lambda\)\)=0, implying thatλ\\lambdais a root if and only ifδ\\deltais an eigenvalue ofW^​\(λ\)\\hat\{W\}\(\\lambda\)\.

Under the assumption of robust marginal stability, there must exist a rootλ0\\lambda\_\{0\}on the unit circle \(i\.e\.,\|λ0\|=1\|\\lambda\_\{0\}\|=1\)\. Then,δ\\deltais an eigenvalue ofW^​\(λ0\)\\hat\{W\}\(\\lambda\_\{0\}\), which implies that

δ≤ρW^​\(λ0\)≤M\.\\delta\\leq\\rho\_\{\\hat\{W\}\}\(\\lambda\_\{0\}\)\\leq M\.\(22\)This contradicts the choiceδ\>M\\delta\>M\. Therefore, there exists a matrixD=\(Δ∗\)−1∈𝒟D=\(\\Delta^\{\*\}\)^\{\-1\}\\in\\mathcal\{D\}such that the characteristic equation has no roots on the unit circle\. Consequently, the system is not robustly marginally stable\. ∎

###### Proof of Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)\.

The second\-order system \([14](https://arxiv.org/html/2610.10690#S3.E14)\) has characteristic equation

det⁡\(z2​𝕀−z⁡\(W¯\+B​Dm−1t​Wψ\)−W¯¯\)=0\.\\det\(z^2\\mathbb\{I\} \- z\(\\overline\{W\} \+ B D\_\{m\-1\}^t W\_\{\\psi\}\) \- \\overline\{\\overline\{W\}\}\)=0\.\(23\)By assumption, matrices \(W¯,W¯¯,B,Wψ\\overline\{W\},\\overline\{\\overline\{W\}\},B,W\_\{\\psi\}\) have block upper\-triangular structure with the same block partition\(n1,…,nk\)\(n\_\{1\},\\ldots,n\_\{k\}\)\. It thus follows that the matrixM⁡\(z\):=z2​𝕀−z⁡\(W¯\+B​Dm−1t​Wψ\)−W¯¯M\(z\)\\mathrel\{\\mathop\{\\mathchar 58\\relax\}\}=z^\{2\}\\mathbb\{I\}\-z\(\\overline\{W\}\+BD\_\{m\-1\}^\{t\}W\_\{\\psi\}\)\-\\overline\{\\overline\{W\}\}is also block upper\-triangular\. Therefore,

det⁡\(M⁡\(z\)\)=∏i=1kdet⁡\(Mi​\(z\)\),\\det\(M\(z\)\)=\\prod\_\{i=1\}^\{k\}\\det\(M\_i\(z\)\),\(24\)whereMi​\(z\)M\_\{i\}\(z\)denotes theiith diagonal block ofM⁡\(z\)M\(z\)\. Using the assumptionW¯¯i=−αi​𝕀ni\\overline\{\\overline\{W\}\}\_\{i\}=\-\\alpha\_\{i\}\\mathbb\{I\}\_\{n\_\{i\}\}, Eq\. \([24](https://arxiv.org/html/2610.10690#A1.E24)\) implies thatdet⁡\(M⁡\(z\)\)=0\\det\(M\(z\)\)=0if and only if

det⁡\(z2​𝕀−z⁡\(W¯i\+Bi​Dm−1t​Wψ,i\)\+αi​𝕀ni\)=0,for some​i∈\{1,…,k\}\.\\det\(z^2\\mathbb I \- z\(\\overline W\_i \+ B\_i D\_\{m\-1\}^t W\_\{\\psi, i\}\) \+ \\alpha\_i \\mathbb I\_\{n\_i\}\)=0,\\,\\,\\text\{for some \}i\\in\\\{1,\\ldots,k\\\}\.\(25\)Therefore, the set of eigenvalues of the full systemM⁡\(z\)M\(z\)is equal to the union of the eigenvalues of each subsystemMi​\(z\)M\_\{i\}\(z\)\.

For each blockii, letNi=W¯i\+Bi​Dm−1t​Wψ,iN\_\{i\}=\\overline\{W\}\_\{i\}\+B\_\{i\}D\_\{m\-1\}^\{t\}W\_\{\\psi,i\}\. SinceW¯i\\overline\{W\}\_\{i\}is symmetric,Dm−1tD\_\{m\-1\}^\{t\}is diagonal, andWψ,i=Σi​Bi⊤W\_\{\\psi,i\}=\\Sigma\_\{i\}B\_\{i\}^\{\\top\}\(withΣi\\Sigma\_\{i\}theiith diagonal block ofΣ\\Sigma\),NiN\_\{i\}is a real symmetric matrix\. It thus holds thatNiN\_\{i\}is diagonalizable and has real eigenvaluesμj\\mu\_\{j\},j=1,…,nij=1,\\dots,n\_\{i\}\. We can thus decouple the characteristic equation \([25](https://arxiv.org/html/2610.10690#A1.E25)\) intonin\_\{i\}independent quadratic equations:

z2−μjz\+αi=0,forj=1,…,ni\.z^\{2\}\-\\mu\_\{j\}z\+\\alpha\_\{i\}=0,\\quad\\text\{for \}j=1,\\dots,n\_\{i\}\.\(26\)The roots of these quadratic equations are given by

zj,±=μj±μj2−4​αi2\.z\_\{j,\\pm\}=\\frac\{\\mu\_\{j\}\\pm\\sqrt\{\\mu\_\{j\}^\{2\}\-4\\alpha\_\{i\}\}\}\{2\}\.\(27\)Given thatDm−1t∈𝒟D\_\{m\-1\}^\{t\}\\in\\mathcal\{D\}, we have‖Dm−1t‖2≤Φ′¯\\\|D\_\{m\-1\}^\{t\}\\\|\_\{2\}\\leq\\overline\{\\Phi^\{\\prime\}\}\. Therefore,

‖W¯i‖2\+Φ′¯​‖Bi‖2​‖Wψ,i‖2<2​αi\\displaystyle\\\|\\overline\{W\}\_\{i\}\\\|\_\{2\}\+\\overline\{\\Phi^\{\\prime\}\}\\\|B\_\{i\}\\\|\_\{2\}\\\|W\_\{\\psi,i\}\\\|\_\{2\}<2\\sqrt\{\\alpha\_\{i\}\}\(28\)‖W¯i\+Bi​Dm−1t​Wψ,i‖2<2​αi\\displaystyle\\\|\\overline\{W\}\_\{i\}\+B\_\{i\}D\_\{m\-1\}^\{t\}W\_\{\\psi,i\}\\\|\_\{2\}<2\\sqrt\{\\alpha\_\{i\}\}\(29\)\|μj\|<2​αi\.\\displaystyle\|\\mu\_\{j\}\|<2\\sqrt\{\\alpha\_\{i\}\}\.\(30\)
This condition ensures that the discriminant of Eq\. \([26](https://arxiv.org/html/2610.10690#A1.E26)\) is negative \(i\.e\.,μj2−4​αi<0\\mu\_\{j\}^\{2\}\-4\\alpha\_\{i\}<0\), which ensures that the rootszj,±z\_\{j,\\pm\}are simple and have a nonzero imaginary part\. Moreover, the absolute values of the roots can be seen to lie on the complex circle of radiusαi\\sqrt\{\\alpha\_\{i\}\},

\|zj,±\|2=\(μj2\)2\+\(4​αi−μj22\)2=αi,\|z\_\{j,\\pm\}\|^\{2\}=\\left\(\\frac\{\\mu\_\{j\}\}\{2\}\\right\)^\{2\}\+\\left\(\\frac\{\\sqrt\{4\\alpha\_\{i\}\-\\mu\_\{j\}^\{2\}\}\}\{2\}\\right\)^\{2\}=\{\\alpha\_\{i\}\},\(31\)for any matrixDm−1t∈𝒟D\_\{m\-1\}^\{t\}\\in\\mathcal\{D\}\. ∎

### A\.2Derivation of the SpatialRNN Gradient Dynamics

We fix the reference timettand define the Jacobian matricesGmt=∂ht\+m∂htG\_\{m\}^\{t\}=\\frac\{\\partial h\_\{t\+m\}\}\{\\partial h\_\{t\}\}andΨmt=∂ψt\+m∂ht\\Psi^\{t\}\_\{m\}=\\frac\{\\partial\\psi\_\{t\+m\}\}\{\\partial h\_\{t\}\}\. Using the chain rule, it follows that the gradients of the SpatialRNN \([3](https://arxiv.org/html/2610.10690#S2.E3)\) are given by the coupled LTV recurrence equations:

\{Gmt=Dmt​Wψ​Ψmt,Ψtm=W¯Ψtm−1\+W¯¯Ψtm−2\+BGm−1,m≥2,\\begin\{cases\}G\_\{m\}^\{t\}=D\_\{m\}^\{t\}W\_\{\\psi\}\\Psi^\{t\}\_\{m\},\\\\ \\Psi^\{t\}\_\{m\}=\\overline\{W\}\\Psi^\{t\}\_\{m\-1\}\+\\overline\{\\overline\{W\}\}\\Psi^\{t\}\_\{m\-2\}\+BG\_\{m\-1\},\\quad m\\geq 2,\\end\{cases\}\(32\)whereΨ1t=B\\Psi^\{t\}\_\{1\}=B,Ψ0t=0\\Psi^\{t\}\_\{0\}=0, andDmt=diag​\(Φ′​\(Wx​xt\+m\+Wψ​ψt\+m\+bh\)\)D\_\{m\}^\{t\}=\\text\{diag\}\\left\(\\Phi^\{\\prime\}\\left\(W\_\{x\}x\_\{t\+m\}\+W\_\{\\psi\}\\psi\_\{t\+m\}\+b\_\{h\}\\right\)\\right\)\. Substituting the first equation into the second one, we obtain

Ψmt\\displaystyle\\Psi^\{t\}\_\{m\}=W¯​Ψm−1t\+W¯¯​Ψm−2t\+B⁡\(Dm−1t​Wψ​Ψm−1t\)\\displaystyle=\\overline\{W\}\\Psi^\{t\}\_\{m\-1\}\+\\overline\{\\overline\{W\}\}\\Psi^\{t\}\_\{m\-2\}\+B\(D\_\{m\-1\}^\{t\}W\_\{\\psi\}\\Psi^\{t\}\_\{m\-1\}\)=\(W¯\+B​Dm−1t​Wψ\)​Ψm−1t\+W¯¯​Ψm−2t,\\displaystyle=\(\\overline\{W\}\+BD\_\{m\-1\}^\{t\}W\_\{\\psi\}\)\\Psi^\{t\}\_\{m\-1\}\+\\overline\{\\overline\{W\}\}\\Psi^\{t\}\_\{m\-2\},which can be written in matrix form as Eq\. \([14](https://arxiv.org/html/2610.10690#S3.E14)\)\.

### A\.3Robustness of the Gradient Dynamics for∞\\infty\-order RNNs

The following result obtains the gradient dynamics for an∞\\infty\-order RNN whose recurrence matrices are determined by Eq\. \([7](https://arxiv.org/html/2610.10690#S2.E7)\)\. Following Theorem[2\.3](https://arxiv.org/html/2610.10690#S2.Thmtheorem3), such∞\\infty\-order RNN is equivalent to a SpatialRNN, and hence we can directly prove that its gradient dynamics are equivalent to those of the SpatialRNN\.

###### Proposition A\.1\.

Let\{Gmt\}m≥0\\\{G\_\{m\}^\{t\}\\\}\_\{m\\geq 0\}be the gradient sequence defined by the HORNN recurrence \([10](https://arxiv.org/html/2610.10690#S3.E10)\), with delay transition matricesWk=Wψ​W~kW\_\{k\}=W\_\{\\psi\}\\tilde\{W\}\_\{k\}such that\{W~k\}k=1∞\\\{\\tilde\{W\}\_\{k\}\\\}\_\{k=1\}^\{\\infty\}satisfy the relation \([7](https://arxiv.org/html/2610.10690#S2.E7)\) andWψW\_\{\\psi\}is invertible\. For any frozen\-time matrixD∈𝒟D\\in\\mathcal\{D\},λ∈ℂ\\lambda\\in\\mathbb\{C\}is an eigenvalue of the recurrence system \([10](https://arxiv.org/html/2610.10690#S3.E10)\) if and only if

det⁡\(λ2​𝕀−λ⁡\(W¯\+B​D​Wψ\)−W¯¯\)=0\.\\det\(\\lambda^2\\mathbb I \- \\lambda\(\\overline W \+ B D W\_\\psi\) \- \\overline\{\\overline W\}\)=0\.\(33\)

###### Proof\.

Recall thatW^​\(λ\)=∑k=1∞Wk​λ−k=∑k=1∞Wψ​W~k​λ−k\\hat\{W\}\(\\lambda\)=\\sum\_\{k=1\}^\{\\infty\}W\_\{k\}\\lambda^\{\-k\}=\\sum\_\{k=1\}^\{\\infty\}W\_\{\\psi\}\\tilde\{W\}\_\{k\}\\lambda^\{\-k\}and let us definez=λ−1z=\\lambda^\{\-1\}\. SinceWψW\_\{\\psi\}is invertible, we can factor it out ofW^​\(z\)=∑k=1∞Wψ​W~k​zk\\hat\{W\}\(z\)=\\sum\_\{k=1\}^\{\\infty\}W\_\{\\psi\}\\tilde\{W\}\_\{k\}z^\{k\}as follows:

Wψ−1​W^​\(z\)\\displaystyle W\_\{\\psi\}^\{\-1\}\\hat\{W\}\(z\)=W~1​z\+W~2​z2\+∑k=3∞W~k​zk\\displaystyle=\\tilde\{W\}\_\{1\}z\+\\tilde\{W\}\_\{2\}z^\{2\}\+\\sum\_\{k=3\}^\{\\infty\}\\tilde\{W\}\_\{k\}z^\{k\}\(34\)Wψ−1​W^​\(z\)\\displaystyle W\_\{\\psi\}^\{\-1\}\\hat\{W\}\(z\)=W~1​z\+W~2​z2\+∑k=3∞\(W¯​W~k−1\+W¯¯​W~k−2\)​zk\\displaystyle=\\tilde\{W\}\_\{1\}z\+\\tilde\{W\}\_\{2\}z^\{2\}\+\\sum\_\{k=3\}^\{\\infty\}\(\\overline\{W\}\\tilde\{W\}\_\{k\-1\}\+\\overline\{\\overline\{W\}\}\\tilde\{W\}\_\{k\-2\}\)z^\{k\}\(35\)Wψ−1​W^​\(z\)\\displaystyle W\_\{\\psi\}^\{\-1\}\\hat\{W\}\(z\)=W~1​z\+W~2​z2\+W¯​z​∑k=3∞W~k−1​zk−1\+W¯¯​z2​∑k=3∞W~k−2​zk−2\\displaystyle=\\tilde\{W\}\_\{1\}z\+\\tilde\{W\}\_\{2\}z^\{2\}\+\\overline\{W\}z\\sum\_\{k=3\}^\{\\infty\}\\tilde\{W\}\_\{k\-1\}z^\{k\-1\}\+\\overline\{\\overline\{W\}\}z^\{2\}\\sum\_\{k=3\}^\{\\infty\}\\tilde\{W\}\_\{k\-2\}z^\{k\-2\}\(36\)Wψ−1​W^​\(z\)\\displaystyle W\_\{\\psi\}^\{\-1\}\\hat\{W\}\(z\)=W~1​z\+W~2​z2\+z​W¯​W^​\(z\)−z2​W¯​W~1\+z2​W¯¯​W^​\(z\)\\displaystyle=\\tilde\{W\}\_\{1\}z\+\\tilde\{W\}\_\{2\}z^\{2\}\+z\\overline\{W\}\\hat\{W\}\(z\)\-z^\{2\}\\overline\{W\}\\tilde\{W\}\_\{1\}\+z^\{2\}\\overline\{\\overline\{W\}\}\\hat\{W\}\(z\)\(37\)Wψ−1​W^​\(z\)​\(𝕀−z​W¯−z2​W¯¯\)\\displaystyle W\_\{\\psi\}^\{\-1\}\\hat\{W\}\(z\)\\left\(\\mathbb\{I\}\-z\\overline\{W\}\-z^\{2\}\\overline\{\\overline\{W\}\}\\right\)=W~1​z\+\(W~2−W¯​W~1\)​z2\\displaystyle=\\tilde\{W\}\_\{1\}z\+\(\\tilde\{W\}\_\{2\}\-\\overline\{W\}\\tilde\{W\}\_\{1\}\)z^\{2\}\(38\)Wψ−1​W^​\(z\)​\(𝕀−z​W¯−z2​W¯¯\)\\displaystyle W\_\{\\psi\}^\{\-1\}\\hat\{W\}\(z\)\\left\(\\mathbb\{I\}\-z\\overline\{W\}\-z^\{2\}\\overline\{\\overline\{W\}\}\\right\)=B​z\+\(W¯​B−W¯​B\)​z2\\displaystyle=Bz\+\\cancel\{\(\\overline\{W\}B\-\\overline\{W\}B\)\}z^\{2\}\(39\)W^​\(z\)\\displaystyle\\hat\{W\}\(z\)=Wψ​\(𝕀−z​W¯−z2​W¯¯\)−1​B​z\.\\displaystyle=W\_\{\\psi\}\(\\mathbb\{I\}\-z\\overline\{W\}\-z^\{2\}\\overline\{\\overline\{W\}\}\)^\{\-1\}Bz\.\(40\)
From Eq\. \([12](https://arxiv.org/html/2610.10690#S3.E12)\), we have that the poles ofG⁡\(z\)G\(z\)are the solutions todet⁡\(𝕀−D​W^​\(z\)\)=0\\det\(\\mathbb I \- D\\hat W\(z\)\)=0\. Substituting Eq\. \([40](https://arxiv.org/html/2610.10690#A1.E40)\) into it, we have that

det⁡\(𝕀−D​Wψ​\(𝕀−z​W¯−z2​W¯¯\)−1​B​z\)\\displaystyle\\det\(\\mathbb\{I\} \- D W\_\{\\psi\} \(\\mathbb\{I\} \- z\\overline\{W\} \- z^2\\overline\{\\overline\{W\}\}\)^\{\-1\} B z\)=0\\displaystyle=0det⁡\(\(𝕀−z​W¯−z2​W¯¯\)−B​D​Wψ​z\)\\displaystyle\\det\(\(\\mathbb\{I\} \- z\\overline\{W\} \- z^2\\overline\{\\overline\{W\}\}\) \- B D W\_\{\\psi\} z\)=0\\displaystyle=0det⁡\(𝕀−z⁡\(W¯\+B​D​Wψ\)−z2​W¯¯\)\\displaystyle\\det\(\\mathbb\{I\} \- z\(\\overline\{W\} \+ B D W\_\{\\psi\}\) \- z^2\\overline\{\\overline\{W\}\}\)=0,\\displaystyle=0,for everyD∈𝒟D\\in\\mathcal\{D\}\. ∎

### A\.4Emergence of Marginal Stability During SpatialRNN Training

Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)shows that settingαi=1\\alpha\_\{i\}=1is a required condition to guarantee marginal stability within theii\-th block\. Importantly,αi\\alpha\_\{i\}can either be fixed a priori or treated as a trainable parameter\. In the latter case, we show empirically that the SpatialRNN model can spontaneously learn to setα≈1\\alpha\\approx 1during BPTT training\. Fig\.[5](https://arxiv.org/html/2610.10690#A1.F5)illustrates the evolution of the decay parametersαi\\alpha\_\{i\}alongside the model’s accuracy over the training period\. Early in the training, the model fails to capture long\-range dependencies, remaining stuck at chance\-level performance\. The learning dynamics only exhibit a sharp transition when gradient descent pushes at least one parameter \(in this case,α4\\alpha\_\{4\}\) close to11\.

Figure 5:Evolution of the parameters\{αi\}i=1k\\\{\\alpha\_\{i\}\\\}\_\{i=1\}^\{k\}during training on the copy task with sequence lengthT=128T=128\. The dashed lines show the learned values ofαi\\alpha\_\{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\\alpha\_\{4\}\\to 1, after which learning rapidly progresses\.While the SpatialRNN is expressive enough to learn the marginally stable regime \(α→1\\alpha\\to 1\), training can be slow when starting from a random initialization\. To bypass this issue, we explicitly fixαi=1\\alpha\_\{i\}=1for at least one block in all experimental evaluations\.

### A\.5Connection between SpatialRNN and Oscillator Networks

Enforcing a block\-triangular structure in the weight matrices allows us to decompose the SpatialRNN’s gradient eigenvalue problem \([14](https://arxiv.org/html/2610.10690#S3.E14)\) intokkindependent blocks, offering robust control over the moduli of the eigenvalues of each block \(Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)\)\. This construction allows the modes \(eigenvectors\) to vary at each timestep, promoting information mixing and higher model expressivity\.

An alternative method to ensure robust control of the gradient spectrum is to ensure simultaneous diagonalizability of the medium update matricesW¯\\overline\{W\}andW¯¯\\overline\{\\overline\{W\}\}\. As we show next,this condition essentially decouples the wave equation intonnindependent harmonic oscillators\. In this case, the eigenvectors remain fixed over time, preventing information mixing among modes and thereby limiting expressivity\. Thus, we do not employ this approach in our experiments\.

###### Theorem A\.2\(Robust SpatialRNN gradients through simultaneous diagonalization\)\.

Let the medium update matricesW¯\\overline\{W\}andW¯¯\\overline\{\\overline\{W\}\}be simultaneously diagonalizable such thatW¯=S​Λ¯​S−1\\overline\{W\}=S\\overline\{\\Lambda\}S^\{\-1\}andW¯¯=S​Λ¯¯​S−1\\overline\{\\overline\{W\}\}=S\\overline\{\\overline\{\\Lambda\}\}S^\{\-1\}for an invertible matrixS∈ℝn×nS\\in\\mathbb\{R\}^\{n\\times n\}, withΛ¯=diag⁡\(μ1,…,μn\)\\overline\{\\Lambda\}=\\operatorname\{diag\}\(\\mu\_\{1\},\\dots,\\mu\_\{n\}\)forμi∈ℝ\\mu\_\{i\}\\in\\mathbb\{R\}andΛ¯¯=diag⁡\(−α1,…,−αn\)\\overline\{\\overline\{\\Lambda\}\}=\\operatorname\{diag\}\(\-\\alpha\_\{1\},\\dots,\-\\alpha\_\{n\}\)forαi∈ℝ\>0\\alpha\_\{i\}\\in\\mathbb\{R\}\_\{\>0\}\. Suppose the following conditions are satisfied:

1. 1\.B=S​ΛBB=S\\Lambda\_\{B\}andWψ=Λψ​S−1W\_\{\\psi\}=\\Lambda\_\{\\psi\}S^\{\-1\}, whereΛB=diag⁡\(β1,…,βn\)\\Lambda\_\{B\}=\\operatorname\{diag\}\(\\beta\_\{1\},\\dots,\\beta\_\{n\}\)andΛψ=diag⁡\(γ1,…,γn\)\\Lambda\_\{\\psi\}=\\operatorname\{diag\}\(\\gamma\_\{1\},\\dots,\\gamma\_\{n\}\)\.
2. 2\.\|μi\|\+Φ′¯​\|βi​γi\|<2​αi\|\\mu\_\{i\}\|\+\\overline\{\\Phi^\{\\prime\}\}\|\\beta\_\{i\}\\gamma\_\{i\}\|<2\\sqrt\{\\alpha\_\{i\}\}for everyiith eigenvalue, whereΦ′¯:=max⁡\{\|Φmin′\|,\|Φmax′\|\}\\overline\{\\Phi^\{\\prime\}\}\\mathrel\{\\mathop\{\\mathchar 58\\relax\}\}=\\max\\\{\|\\Phi^\{\\prime\}\_\{\\min\}\|,\|\\Phi^\{\\prime\}\_\{\\max\}\|\\\}\.

Then, for any linear matrixDm−1t∈𝒟D\_\{m\-1\}^\{t\}\\in\\mathcal\{D\}, the SpatialRNN gradient dynamics decouple intonnindependent second\-order modes\. For each modeii, the associated characteristic equation yields exactly two simple complex conjugate eigenvalues with a time\-invariant modulus\|λ\|=αi\|\\lambda\|=\\sqrt\{\\alpha\_\{i\}\}\.

###### Proof\.

Recall that the characteristic equation of the second\-order gradient dynamics of the SpatialRNN \([14](https://arxiv.org/html/2610.10690#S3.E14)\) is given by

det⁡\(z2​𝕀−z⁡\(W¯\+B​Dm−1t​Wψ\)−W¯¯\)=0\.\\det\(z^2 \\mathbb\{I\} \- z\(\\overline\{W\} \+ B D\_\{m\-1\}^t W\_\\psi\) \- \\overline\{\\overline\{W\}\}\)=0\.\(41\)Substituting the conditionsW¯=S​Λ¯​S−1\\overline\{W\}=S\\overline\{\\Lambda\}S^\{\-1\},W¯¯=S​Λ¯¯​S−1\\overline\{\\overline\{W\}\}=S\\overline\{\\overline\{\\Lambda\}\}S^\{\-1\},B=S​ΛBB=S\\Lambda\_\{B\}, andWψ=Λψ​S−1W\_\{\\psi\}=\\Lambda\_\{\\psi\}S^\{\-1\}yields

det⁡\(S⁡\[z2​𝕀−z⁡\(Λ¯\+ΛB​Dm−1t​Λψ\)−Λ¯¯\]​S−1\)=0\.\\det\(S \[z^2 \\mathbb\{I\} \- z\(\\overline\{\\Lambda\} \+ \\Lambda\_B D\_\{m\-1\}^t \\Lambda\_\\psi\) \- \\overline\{\\overline\{\\Lambda\}\}\] S^\{\-1\}\)=0\.\(42\)SinceDm−1t∈𝒟D\_\{m\-1\}^\{t\}\\in\\mathcal\{D\}is diagonal, it commutes with the diagonal matricesΛB\\Lambda\_\{B\}andΛψ\\Lambda\_\{\\psi\}\. LetD~=ΛB​Dm−1t​Λψ\\tilde\{D\}=\\Lambda\_\{B\}D\_\{m\-1\}^\{t\}\\Lambda\_\{\\psi\}, which is a diagonal matrix with entriesd~i=βi​γi​di\\tilde\{d\}\_\{i\}=\\beta\_\{i\}\\gamma\_\{i\}d\_\{i\}, wheredi=\[Dm−1t\]i​id\_\{i\}=\[D\_\{m\-1\}^\{t\}\]\_\{ii\}and\|di\|≤Φ′¯\|d\_\{i\}\|\\leq\\overline\{\\Phi^\{\\prime\}\}\. By the similarity invariance of the determinant, the system fully decouples intonnindependent quadratic equations:

z2−\(μi\+βiγidi\)z\+αi=0,fori=1,…,n\.z^\{2\}\-\(\\mu\_\{i\}\+\\beta\_\{i\}\\gamma\_\{i\}d\_\{i\}\)z\+\\alpha\_\{i\}=0,\\quad\\text\{for \}i=1,\\dots,n\.\(43\)The roots for each modeiiare:

zi,±=\(μi\+βi​γi​di\)±\(μi\+βi​γi​di\)2−4​αi2\.z\_\{i,\\pm\}=\\frac\{\(\\mu\_\{i\}\+\\beta\_\{i\}\\gamma\_\{i\}d\_\{i\}\)\\pm\\sqrt\{\(\\mu\_\{i\}\+\\beta\_\{i\}\\gamma\_\{i\}d\_\{i\}\)^\{2\}\-4\\alpha\_\{i\}\}\}\{2\}\.\(44\)To guarantee simple, complex conjugate roots \(i\.e\., underdamped oscillatory modes\), the discriminant must be strictly negative:\(μi\+βi​γi​di\)2<4​αi\(\\mu\_\{i\}\+\\beta\_\{i\}\\gamma\_\{i\}d\_\{i\}\)^\{2\}<4\\alpha\_\{i\}\. Applying the triangle inequality and the activation derivative bound, we have that

\|μi\+βi​γi​di\|≤\|μi\|\+\|βi​γi\|\|di\|≤\|μi\|\+Φ′¯​\|βi​γi\|\.\|\\mu\_\{i\}\+\\beta\_\{i\}\\gamma\_\{i\}d\_\{i\}\|\\leq\|\\mu\_\{i\}\|\+\|\\beta\_\{i\}\\gamma\_\{i\}\|\|d\_\{i\}\|\\leq\|\\mu\_\{i\}\|\+\\overline\{\\Phi^\{\\prime\}\}\|\\beta\_\{i\}\\gamma\_\{i\}\|\.\(45\)By assumption,\|μi\|\+Φ′¯​\|βi​γi\|<2​αi\|\\mu\_\{i\}\|\+\\overline\{\\Phi^\{\\prime\}\}\|\\beta\_\{i\}\\gamma\_\{i\}\|<2\\sqrt\{\\alpha\_\{i\}\}\. Thus, the discriminant is strictly negative for any LTV sequence matrixDm−1t∈𝒟D\_\{m\-1\}^\{t\}\\in\\mathcal\{D\}\. Consequently, the roots are complex conjugates, and their squared magnitude evaluates to

\|zi,±\|2=\(μi\+βi​γi​di2\)2\+4​αi−\(μi\+βi​γi​di\)24=αi\.\|z\_\{i,\\pm\}\|^\{2\}=\\left\(\\frac\{\\mu\_\{i\}\+\\beta\_\{i\}\\gamma\_\{i\}d\_\{i\}\}\{2\}\\right\)^\{2\}\+\\frac\{4\\alpha\_\{i\}\-\(\\mu\_\{i\}\+\\beta\_\{i\}\\gamma\_\{i\}d\_\{i\}\)^\{2\}\}\{4\}=\\alpha\_\{i\}\.\(46\)Therefore,\|zi,±\|=αi\|z\_\{i,\\pm\}\|=\\sqrt\{\\alpha\_\{i\}\}, establishing robust eigenvalue placement for all modes invariant to the time\-varying sequence\. ∎

### A\.6LTV Stability Conditions for the SpatialRNN Gradient Dynamics

The conditions derived in Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)ensure marginal stability of the gradient dynamics for any frozen\-time matrixD∈𝒟D\\in\\mathcal\{D\}\. However, given that the input sequences\{xt\}t≥0\\\{x\_\{t\}\\\}\_\{t\\geq 0\}are unknown and therefore\{Dmt\}m≥0\\\{D\_\{m\}^\{t\}\\\}\_\{m\\geq 0\}can switch arbitrarily, marginal stability of the frozen\-time systems is only a necessary, but not sufficient, condition for the marginal stability of the switching LTV gradient equations \([13](https://arxiv.org/html/2610.10690#S3.E13)\)\. Empirically, we observe that gradient norms can exhibit mild drift over long horizons, although deviations are relatively small, as shown in Fig\.[6](https://arxiv.org/html/2610.10690#A1.F6)\.

Figure 6:Gradient norm‖∂ℒ∂hT−t‖\\\|\\frac\{\\partial\\mathcal\{L\}\}\{\\partial h\_\{T\-t\}\}\\\|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\|\|B\_\{i\}\|\|\_\{2\}<0\.4,\\ \|\|W\_\{\\psi,i\}\|\|\_\{2\}<0\.4and‖W¯i‖2<1\.8\|\|\\overline\{W\}\_\{i\}\|\|\_\{2\}<1\.8\.In order to prevent exploding gradients, we provide further constraints on the SpatialRNN’s weight matrices to ensureasymptotic stabilityof the actual LTV gradient dynamics\. These results are weaker than Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)in the sense that we cannot control the decay rate of the system, nor can we guarantee marginal stability \(αi=1\\alpha\_\{i\}=1\), yet they provide concrete bounds on the gradient dynamics\.

###### Theorem A\.3\(LTV stability of the SpatialRNN gradients\)\.

Let the matricesW¯\\overline\{W\}andW¯¯\\overline\{\\overline\{W\}\}be defined as block diagonal matrices with blocks sizesnin\_\{i\},i=1,…,ki=1,\\ldots,k, such that

W¯=\(W¯1000⋱000W¯k\)andW¯¯=\(−α1​𝕀n1000⋱000−αk​𝕀nk\),\\overline\{W\}=\\begin\{pmatrix\}\\overline\{W\}\_\{1\}&0&0\\\\ 0&\\ddots&0\\\\ 0&0&\\overline\{W\}\_\{k\}\\end\{pmatrix\}\\quad\\text\{and\}\\quad\\overline\{\\overline\{W\}\}=\\begin\{pmatrix\}\-\\alpha\_\{1\}\\mathbb\{I\}\_\{n\_\{1\}\}&0&0\\\\ 0&\\ddots&0\\\\ 0&0&\-\\alpha\_\{k\}\\mathbb\{I\}\_\{n\_\{k\}\}\\end\{pmatrix\},\(47\)whereW¯i∈ℝni×ni\\overline\{W\}\_\{i\}\\in\\mathbb\{R\}^\{n\_\{i\}\\times n\_\{i\}\},αi∈\(0,1\)\\alpha\_\{i\}\\in\(0,1\),∑i=1kni=n\\sum\_\{i=1\}^\{k\}n\_\{i\}=nis the model size\. Suppose the following conditions are satisfied:

1. 1\.Wi¯=Wi¯⊤\\overline\{W\_\{i\}\}=\\overline\{W\_\{i\}\}^\{\\top\},B=blkdiag⁡\{B1,…,Bk\}B=\\operatorname\{blkdiag\}\\\{B\_\{1\},\\ldots,B\_\{k\}\\\}, for block matricesBi∈ℝni×niB\_\{i\}\\in\\mathbb\{R\}^\{n\_\{i\}\\times n\_\{i\}\}, andWψ=Σ​B⊤=blkdiag⁡\{Wψ,1,…,Wψ,k\}W\_\{\\psi\}=\\Sigma B^\{\\top\}=\\operatorname\{blkdiag\}\\\{W\_\{\\psi,1\},\\ldots,W\_\{\\psi,k\}\\\}for some diagonal matrixΣ∈diag⁡\(ℝn\)\\Sigma\\in\\operatorname\{diag\}\(\\mathbb\{R\}^\{n\}\)\.
2. 2\.‖W¯i‖2<1\+αiαi\|\|\\overline\{W\}\_\{i\}\|\|\_\{2\}<\\frac\{1\+\\alpha\_\{i\}\}\{\\sqrt\{\\alpha\_\{i\}\}\}andΦ′¯​‖Bi‖2​‖Wψ,i‖2<\(1−αi\)​\(1−αi\(1\+αi\)2​‖W¯i‖22−‖W¯i‖21\+αi\)​∀i\\overline\{\\Phi^\{\\prime\}\}\|\|B\_\{i\}\|\|\_\{2\}\|\|W\_\{\\psi,i\}\|\|\_\{2\}<\(1\-\\alpha\_\{i\}\)\\left\(\\sqrt\{1\-\\frac\{\\alpha\_\{i\}\}\{\(1\+\\alpha\_\{i\}\)^\{2\}\}\|\|\\overline\{W\}\_\{i\}\|\|\_\{2\}^\{2\}\}\-\\frac\{\|\|\\overline\{W\}\_\{i\}\|\|\_\{2\}\}\{1\+\\alpha\_\{i\}\}\\right\)\\ \\forall i, whereΦ′¯:=max⁡\{\|Φmin′\|,\|Φmax′\|\}\\overline\{\\Phi^\{\\prime\}\}\\mathrel\{\\mathop\{\\mathchar 58\\relax\}\}=\\max\\\{\|\\Phi^\{\\prime\}\_\{\\min\}\|,\|\\Phi^\{\\prime\}\_\{\\max\}\|\\\}\.

Then, the LTV gradient dynamicsΨmt\\Psi\_\{m\}^\{t\}are uniformly asymptotically stable for any matrix sequence\{Dmt\}m≥0,Dmt∈𝒟\\\{D\_\{m\}^\{t\}\\\}\_\{m\\geq 0\},\\ D\_\{m\}^\{t\}\\in\\mathcal\{D\}\.

###### Proof\.

The block\-diagonal structure of the parameter matrices naturally decouples the gradient dynamics intokkindependent LTV subsystems\. For theii\-th block, the recurrence is given byxm\(i\)=Am,i​xm−1\(i\)x^\{\(i\)\}\_\{m\}=A\_\{m,i\}x^\{\(i\)\}\_\{m\-1\}, wherexm\(i\)=\[Ψm,itΨm−1,it\]x^\{\(i\)\}\_\{m\}=\\begin\{bmatrix\}\\Psi\_\{m,i\}^\{t\}\\\\ \\Psi\_\{m\-1,i\}^\{t\}\\end\{bmatrix\},Am,i=\[Wm,i−αi​II0\]A\_\{m,i\}=\\begin\{bmatrix\}W\_\{m,i\}&\-\\alpha\_\{i\}I\\\\ I&0\\end\{bmatrix\},Wm,i=W¯i\+W~m,iW\_\{m,i\}=\\overline\{W\}\_\{i\}\+\\tilde\{W\}\_\{m,i\}, andW~m,i=Bi​Dm−1t​Wψ,i\\tilde\{W\}\_\{m,i\}=B\_\{i\}D\_\{m\-1\}^\{t\}W\_\{\\psi,i\}\. By definition, the spectral norm of the perturbation is bounded:‖W~m,i‖2≤w~i\\\|\\tilde\{W\}\_\{m,i\}\\\|\_\{2\}\\leq\\tilde\{w\}\_\{i\}for anyDm−1t∈𝒟D\_\{m\-1\}^\{t\}\\in\\mathcal\{D\}\.

To guarantee stability under arbitrary sequential variations, we construct a time\-invariant quadratic Lyapunov function for each LTV subsystem:Vi​\(x\)=x⊤​Pi​xV\_\{i\}\(x\)=x^\{\\top\}P\_\{i\}x\. If a fixed positive\-definite matrixPiP\_\{i\}exists such that the energy differenceΔ​Vi,m=Am,i⊤​Pi​Am,i−Pi\\Delta V\_\{i,m\}=A\_\{m,i\}^\{\\top\}P\_\{i\}A\_\{m,i\}\-P\_\{i\}is strictly negative definite at every step, uniform stability is guaranteed\([Liberzon, 2003](https://arxiv.org/html/2610.10690#bib.bib8)\)\. We definePiP\_\{i\}as

Pi=\[I−γi​W¯i−γi​W¯iαi​I\],whereγi=αi1\+αi\.P\_\{i\}=\\begin\{bmatrix\}I&\-\\gamma\_\{i\}\\overline\{W\}\_\{i\}\\\\ \-\\gamma\_\{i\}\\overline\{W\}\_\{i\}&\\alpha\_\{i\}I\\end\{bmatrix\},\\quad\\text\{where\}\\quad\\gamma\_\{i\}=\\frac\{\\alpha\_\{i\}\}\{1\+\\alpha\_\{i\}\}\.\(48\)By the Schur complement,Pi≻0P\_\{i\}\\succ 0if and only ifαi​I−γi2​W¯i2≻0\\alpha\_\{i\}I\-\\gamma\_\{i\}^\{2\}\\overline\{W\}\_\{i\}^\{2\}\\succ 0\. This holds precisely when‖W¯i‖2<1\+αiαi\\\|\\overline\{W\}\_\{i\}\\\|\_\{2\}<\\frac\{1\+\\alpha\_\{i\}\}\{\\sqrt\{\\alpha\_\{i\}\}\}, which is imposed by condition 2\.

Next, we evaluate the differenceΔ​Vi,m\\Delta V\_\{i,m\}\. ExpandingAm,i⊤​Pi​Am,iA\_\{m,i\}^\{\\top\}P\_\{i\}A\_\{m,i\}and substitutingγi\\gamma\_\{i\}forces the linearW¯i\\overline\{W\}\_\{i\}terms in the off\-diagonal blocks to perfectly cancel out, yielding:

Δ​Vi,m=\[S11−αi​W~m,i−αi​W~m,i−αi​\(1−αi\)​I\],\\Delta V\_\{i,m\}=\\begin\{bmatrix\}S\_\{11\}&\-\\alpha\_\{i\}\\tilde\{W\}\_\{m,i\}\\\\ \-\\alpha\_\{i\}\\tilde\{W\}\_\{m,i\}&\-\\alpha\_\{i\}\(1\-\\alpha\_\{i\}\)I\\end\{bmatrix\},\(49\)whereS11=W~m,i2\+11\+αi​\(W¯i​W~m,i\+W~m,i​W¯i\)\+1−αi1\+αi​W¯i2−\(1−αi\)​IS\_\{11\}=\\tilde\{W\}\_\{m,i\}^\{2\}\+\\frac\{1\}\{1\+\\alpha\_\{i\}\}\(\\overline\{W\}\_\{i\}\\tilde\{W\}\_\{m,i\}\+\\tilde\{W\}\_\{m,i\}\\overline\{W\}\_\{i\}\)\+\\frac\{1\-\\alpha\_\{i\}\}\{1\+\\alpha\_\{i\}\}\\overline\{W\}\_\{i\}^\{2\}\-\(1\-\\alpha\_\{i\}\)I\. Sinceαi∈\(0,1\)\\alpha\_\{i\}\\in\(0,1\), the bottom\-right blockS22=−αi​\(1−αi\)​IS\_\{22\}=\-\\alpha\_\{i\}\(1\-\\alpha\_\{i\}\)Iis strictly negative definite\. Applying the Schur complement toΔ​Vi,m\\Delta V\_\{i,m\}, we requireS11−S12​S22−1​S12⊤≺0S\_\{11\}\-S\_\{12\}S\_\{22\}^\{\-1\}S\_\{12\}^\{\\top\}\\prec 0\. Thus, we have that

11−αi​W~m,i2\+11\+αi​\(W¯i​W~m,i\+W~m,i​W¯i\)\+1−αi1\+αi​W¯i2−\(1−αi\)​I≺0\.\\frac\{1\}\{1\-\\alpha\_\{i\}\}\\tilde\{W\}\_\{m,i\}^\{2\}\+\\frac\{1\}\{1\+\\alpha\_\{i\}\}\(\\overline\{W\}\_\{i\}\\tilde\{W\}\_\{m,i\}\+\\tilde\{W\}\_\{m,i\}\\overline\{W\}\_\{i\}\)\+\\frac\{1\-\\alpha\_\{i\}\}\{1\+\\alpha\_\{i\}\}\\overline\{W\}\_\{i\}^\{2\}\-\(1\-\\alpha\_\{i\}\)I\\prec 0\.\(50\)To enforce strict negative definiteness, we bound the spectral norm\. Using the triangle inequality and the property‖W¯i​W~m,i\+W~m,i​W¯i‖2≤2​‖W¯i‖2​‖W~m,i‖2\\\|\\overline\{W\}\_\{i\}\\tilde\{W\}\_\{m,i\}\+\\tilde\{W\}\_\{m,i\}\\overline\{W\}\_\{i\}\\\|\_\{2\}\\leq 2\\\|\\overline\{W\}\_\{i\}\\\|\_\{2\}\\\|\\tilde\{W\}\_\{m,i\}\\\|\_\{2\}, we obtain the sufficient scalar inequality

11−αi​‖W~m,i‖22\+2​‖W¯i‖21\+αi​‖W~m,i‖2\+1−αi1\+αi​‖W¯i‖22−\(1−αi\)<0\.\\frac\{1\}\{1\-\\alpha\_\{i\}\}\\\|\\tilde\{W\}\_\{m,i\}\\\|\_\{2\}^\{2\}\+\\frac\{2\\\|\\overline\{W\}\_\{i\}\\\|\_\{2\}\}\{1\+\\alpha\_\{i\}\}\\\|\\tilde\{W\}\_\{m,i\}\\\|\_\{2\}\+\\frac\{1\-\\alpha\_\{i\}\}\{1\+\\alpha\_\{i\}\}\\\|\\overline\{W\}\_\{i\}\\\|\_\{2\}^\{2\}\-\(1\-\\alpha\_\{i\}\)<0\.\(51\)Multiplying by\(1−αi\)\(1\-\\alpha\_\{i\}\), we can solve this quadratic inequality for‖W~m,i‖2\\\|\\tilde\{W\}\_\{m,i\}\\\|\_\{2\}\. The following strict upper bound, corresponding to the positive root of Eq\. \([51](https://arxiv.org/html/2610.10690#A1.E51)\), guarantees stability:

‖W~m,i‖2<\(1−αi\)​\(1−αi\(1\+αi\)2​‖W¯i‖22−‖W¯i‖21\+αi\)\.\\\|\\tilde\{W\}\_\{m,i\}\\\|\_\{2\}<\(1\-\\alpha\_\{i\}\)\\left\(\\sqrt\{1\-\\frac\{\\alpha\_\{i\}\}\{\(1\+\\alpha\_\{i\}\)^\{2\}\}\\\|\\overline\{W\}\_\{i\}\\\|\_\{2\}^\{2\}\}\-\\frac\{\\\|\\overline\{W\}\_\{i\}\\\|\_\{2\}\}\{1\+\\alpha\_\{i\}\}\\right\)\.\(52\)Condition 2 guarantees this inequality holds, ensuringΔ​Vi,m≺0\\Delta V\_\{i,m\}\\prec 0uniformly over time\. Summing the block\-wise Lyapunov functions,V⁡\(x\)=∑i=1kVi​\(x\(i\)\)V\(x\)=\\sum\_\{i=1\}^\{k\}V\_\{i\}\(x^\{\(i\)\}\)forms a global time\-invariant Lyapunov function, establishing uniform asymptotic stability for the full LTV system\. ∎

While condition 1 can be relaxed to the block\-triangular structure considered in Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3), the most significant change is the bound on the matrix norms\. From Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3), to ensure marginal stability we requireαi→1\\alpha\_\{i\}\\to 1\. However, due to condition 2 of Theorem[A\.3](https://arxiv.org/html/2610.10690#A1.Thmtheorem3), we would need to impose the constraint‖Bi‖2​‖Wψ,i‖2→0\\\|B\_\{i\}\\\|\_\{2\}\\\|W\_\{\\psi,i\}\\\|\_\{2\}\\to 0for the state\-to\-medium and medium\-to\-state matrices\. Such constraint would render the SpatialRNN untrainable\. Thus, instead we enforce these matrices to have small \(but nonvanishing\) norms, which completely prevents exploding gradients at the cost of introducing some degree of gradient decay\. Consequently, under Theorem[A\.3](https://arxiv.org/html/2610.10690#A1.Thmtheorem3), the LTV gradient dynamics of the SpatialRNN are asymptotically stable, but not necessarily marginally stable\.

Input\-independent volume contraction\.Although Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)is a frozen\-time result, part of the structure that it exploits carries over exactly to the LTV gradient dynamics\. Since the input\-dependent termB​Dm−1t​WψBD\_\{m\-1\}^\{t\}W\_\{\\psi\}only enters the first block of the companion matrix in Eq\. \([14](https://arxiv.org/html/2610.10690#S3.E14)\), the determinant of the gradient transition matrix does not depend on the input sequence\.

###### Proposition A\.4\(Input\-independent volume contraction\)\.

LetW¯¯\\overline\{\\overline\{W\}\}be defined as in Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3),𝒜m\\mathcal\{A\}\_\{m\}denote the companion matrix of the gradient dynamics \([14](https://arxiv.org/html/2610.10690#S3.E14)\) at stepmm, andΦm:=𝒜m𝒜m−1⋯𝒜2\\Phi\_\{m\}\\mathrel\{\\mathop\{\\mathchar 58\\relax\}\}=\\mathcal\{A\}\_\{m\}\\mathcal\{A\}\_\{m\-1\}\\cdots\\mathcal\{A\}\_\{2\}\. Then, for any sequence\{Dmt\}m≥1\\\{D\_\{m\}^\{t\}\\\}\_\{m\\geq 1\}, withDmt∈𝒟D\_\{m\}^\{t\}\\in\\mathcal\{D\}, it holds that

\|det⁡Φm\|=\(∏i=1kαini\)m−1,m≥2\.\|\\det\\Phi\_\{m\}\|=\\Big\(\\prod\\nolimits\_\{i=1\}^\{k\}\\alpha\_\{i\}^\{n\_\{i\}\}\\Big\)^\{m\-1\},\\quad m\\geq 2\.\(53\)In particular, ifαi=1\\alpha\_\{i\}=1for alli=1,…,ki=1,\\ldots,k, then‖Φm‖2≥1\\\|\\Phi\_\{m\}\\\|\_\{2\}\\geq 1for allm≥2m\\geq 2\.

###### Proof\.

Exchanging the two block rows of𝒜m\\mathcal\{A\}\_\{m\}, each of sizenn, yields a block lower\-triangular matrix with diagonal blocks𝕀\\mathbb\{I\}andW¯¯\\overline\{\\overline\{W\}\}\. Since this exchange is a permutation with sign\(−1\)n\(\-1\)^\{n\}, andW¯¯\\overline\{\\overline\{W\}\}is block upper\-triangular with diagonal blocks−αi​𝕀ni\-\\alpha\_\{i\}\\mathbb\{I\}\_\{n\_\{i\}\}, it follows thatdet⁡𝒜m=\(−1\)n​det⁡W¯¯=∏i=1kαini\\det\\mathcal\{A\}\_\{m\}=\(\-1\)^\{n\}\\det\\overline\{\\overline\{W\}\}=\\prod\_\{i=1\}^\{k\}\\alpha\_\{i\}^\{n\_\{i\}\}, independently ofDm−1tD\_\{m\-1\}^\{t\}\. The first claim then follows from the multiplicativity of the determinant\. For the second claim, note that\|det⁡Φm\|\|\\det\\Phi\_\{m\}\|equals the product of the singular values ofΦm\\Phi\_\{m\}; if\|det⁡Φm\|=1\|\\det\\Phi\_\{m\}\|=1, then its largest singular value satisfies‖Φm‖2≥1\\\|\\Phi\_\{m\}\\\|\_\{2\}\\geq 1\. ∎

Consequently, forαi=1\\alpha\_\{i\}=1, the LTV gradient dynamics of the SpatialRNN are volume preserving for any input sequence, and the medium gradients cannot vanish along all directions\. In contrast, the companion matrix of the HORNN gradient dynamics \([10](https://arxiv.org/html/2610.10690#S3.E10)\) has determinant±det⁡\(Dmt\)​det⁡\(Wd\)\\pm\\det\(D\_m^t\)\\det\(W\_d\), which depends on the input and vanishes as the activations saturate\. Volume preservation is, however, necessary but not sufficient for marginal stability, as it does not preclude growth along some directions compensated by decay along others\.

### A\.7Empirical LTV Gradient Dynamics along Input Trajectories

To quantify the gap between the frozen\-time analysis and the full LTV gradient dynamics, we compute the exact state\-to\-state JacobiansGmτ=∂hτ\+m/∂hτG^\{\\tau\}\_\{m\}=\\partial h\_\{\\tau\+m\}/\\partial h\_\{\\tau\}of a SpatialRNN along real input trajectories\. To this end, we propagate the recurrence equations \([13](https://arxiv.org/html/2610.10690#S3.E13)\)–\([14](https://arxiv.org/html/2610.10690#S3.E14)\) using the activation derivativesDmτD^\{\\tau\}\_\{m\}realized by each input sequence\. In contrast with the loss gradients in Fig\.[6](https://arxiv.org/html/2610.10690#A1.F6), this measurement does not depend on the loss, the readout, or the label, and thus directly probes the LTV dynamics whose frozen\-time spectrum is controlled by Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)\. We consider a SpatialRNN with block sizes\{128,8,8,8,8\}\\\{128,8,8,8,8\\\}trained on psMNIST, withα1=1\\alpha\_\{1\}=1fixed and learned decay parametersαi∈\[0\.31,0\.46\]\\alpha\_\{i\}\\in\[0\.31,0\.46\]fori≥2i\\geq 2\(which reached96\.2%96\.2\\%validation accuracy during training\), and the same architecture at initialization \(withαi=0\.2\\alpha\_\{i\}=0\.2fori≥2i\\geq 2\)\. For each model, we evaluate‖Gmτ‖F\\\|G^\{\\tau\}\_\{m\}\\\|\_\{F\}forτ=0\\tau=0andm≤783m\\leq 783along100100psMNIST validation sequences in double precision, rescaling the recursion at every step to prevent overflow and underflow\. For each sequence, we fit a log\-linear regression to‖Gmτ‖F\\\|G^\{\\tau\}\_\{m\}\\\|\_\{F\}form≥50m\\geq 50, whose slope provides an empirical estimate of the top Lyapunov exponent of the LTV gradient dynamics\.

Fig\.[7](https://arxiv.org/html/2610.10690#A1.F7)and Table[3](https://arxiv.org/html/2610.10690#A1.T3)show that the Jacobian norms neither vanish nor explode over the full horizon\. For the trained model, the average slope is7\.98⋅10−57\.98\\cdot 10^\{\-5\}per step, corresponding to a change of only≈6%\{\\approx\}6\\%in the Jacobian norm across783783steps, compared with≈23%\{\\approx\}23\\%at initialization\. Training thus moves the model closer to the marginal regime predicted by the frozen\-time analysis and reduces the residual drift due to the time\-varying dynamics\.

Figure 7:Frobenius norm of the exact state\-to\-state JacobiansGmτ=∂hτ\+m/∂hτG^\{\\tau\}\_\{m\}=\\partial h\_\{\\tau\+m\}/\\partial h\_\{\\tau\}, propagated via Eqs\. \([13](https://arxiv.org/html/2610.10690#S3.E13)\)–\([14](https://arxiv.org/html/2610.10690#S3.E14)\) along100100psMNIST validation sequences \(τ=0\\tau=0,T=784T=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\\pm 1standard deviation \(in log scale\)\.Table 3:Empirical LTV gradient statistics over100100psMNIST validation sequences \(mean±\\pmstandard deviation across sequences\)\. Slopes are fitted form≥50m\\geq 50, and the coefficient of variation \(CV\) measures the dispersion of‖Gmτ‖F\\\|G^\{\\tau\}\_\{m\}\\\|\_\{F\}along each sequence\.

## Appendix BSymplectic Form Preservation of the SpatialRNN

When the weight matrices are constrained to ensure marginal stability \(αi=1\\alpha\_\{i\}=1\), the corresponding gradient dynamics can maintain a symplectic form\([Marsden and West, 2001](https://arxiv.org/html/2610.10690#bib.bib51)\)\. A map is symplectic if it preserves the canonical geometry of phase space \(i\.e\., it conserves a bilinear form that encodes position\-momentum interactions\) and therefore preserves the phase\-space volume\. Unlike dissipative dynamics, symplectic dynamics preserve phase\-space volume and are invertible, so information is not irreversibly lost through volume contraction\. Consequently, this property leads to numerically stable and time\-reversible evolution, allowing past states to be recovered exactly from future ones \(as shown below in Theorem[B\.3](https://arxiv.org/html/2610.10690#A2.Thmtheorem3)\)\. Here, we first formalize the condition for symplectic form preservation:

###### Proposition B\.1\(Symplectic preservation of second\-order LDEs\)\.

Consider the unforced medium dynamics

ψt\+1=W¯​ψt\+W¯¯​ψt−1,\\psi\_\{t\+1\}=\\overline\{W\}\\psi\_\{t\}\+\\overline\{\\overline\{W\}\}\\psi\_\{t\-1\},\(54\)and define the canonical position\-momentum coordinateszt=\(qt,pt\)Tz\_\{t\}=\(q\_\{t\},p\_\{t\}\)^\{T\}, whereqt=ψtq\_\{t\}=\\psi\_\{t\}andpt=1Δ​t​\(ψt−ψt−1\)p\_\{t\}=\\frac\{1\}\{\\Delta t\}\(\\psi\_\{t\}\-\\psi\_\{t\-1\}\)\. In these coordinates, the evolution defines a one\-step linear mapzt\+1=M​ztz\_\{t\+1\}=Mz\_\{t\}\. This map is symplectic \(i\.e\.,M⊤​J​M=JM^\{\\top\}JM=Jfor the canonical bilinear formJJ\) if and only ifW¯¯=−𝕀\\overline\{\\overline\{W\}\}=\-\\mathbb\{I\}andW¯⊤=W¯\\overline\{W\}^\{\\top\}=\\overline\{W\}\.

###### Proof\.

The one\-step mapMMis derived from the canonical coordinate definitions:

qt\+1\\displaystyle q\_\{t\+1\}=ψt\+1=W¯​qt\+W¯¯​\(qt−Δ​t​pt\)=\(W¯\+W¯¯\)​qt−Δ​t​W¯¯​pt,\\displaystyle=\\psi\_\{t\+1\}=\\overline\{W\}q\_\{t\}\+\\overline\{\\overline\{W\}\}\(q\_\{t\}\-\{\\Delta t\}p\_\{t\}\)=\(\\overline\{W\}\+\\overline\{\\overline\{W\}\}\)q\_\{t\}\-\{\\Delta t\}\\overline\{\\overline\{W\}\}p\_\{t\},pt\+1\\displaystyle p\_\{t\+1\}=Δ​t−1​\(qt\+1−qt\)=1Δ​t​\(W¯\+W¯¯−𝕀\)​qt−W¯¯​pt\.\\displaystyle=\{\\Delta t\}^\{\-1\}\(q\_\{t\+1\}\-q\_\{t\}\)=\\frac\{1\}\{\\Delta t\}\(\\overline\{W\}\+\\overline\{\\overline\{W\}\}\-\\mathbb\{I\}\)q\_\{t\}\-\\overline\{\\overline\{W\}\}p\_\{t\}\.Therefore, the dynamics can be written as a linear mapzt\+1=M​ztz\_\{t\+1\}=Mz\_\{t\}, with the propagator

M=\(W¯\+W¯¯−Δ​t​W¯¯1Δ​t​\(W¯\+W¯¯−𝕀\)−W¯¯\)≡\(ABCD\)\.M=\\begin\{pmatrix\}\\overline\{W\}\+\\overline\{\\overline\{W\}\}&\-\{\\Delta t\}\\overline\{\\overline\{W\}\}\\\\ \\frac\{1\}\{\\Delta t\}\(\\overline\{W\}\+\\overline\{\\overline\{W\}\}\-\\mathbb\{I\}\)&\-\\overline\{\\overline\{W\}\}\\end\{pmatrix\}\\equiv\\begin\{pmatrix\}A&B\\\\ C&D\\end\{pmatrix\}\.\(55\)The map is symplectic ifM⊤​J​M=J=\(0𝕀−𝕀0\)M^\{\\top\}JM=J=\\begin\{pmatrix\}0&\\mathbb\{I\}\\\\ \-\\mathbb\{I\}&0\\end\{pmatrix\}\. This requires four block\-matrix conditions: 1\)A⊤​D−C⊤​B=𝕀A^\{\\top\}D\-C^\{\\top\}B=\\mathbb\{I\}, 2\)A⊤​C−C⊤​A=0A^\{\\top\}C\-C^\{\\top\}A=0, 3\)B⊤​D−D⊤​B=0B^\{\\top\}D\-D^\{\\top\}B=0, andB⊤​C−D⊤​A=−𝕀B^\{\\top\}C\-D^\{\\top\}A=\-\\mathbb\{I\}\.

Condition 1\. Substituting the expressions for\(A,B,C,D\)\(A,B,C,D\), we have that

\(W¯\+W¯¯\)⊤​\(−W¯¯\)−\(1Δ​t​\(W¯\+W¯¯−𝕀\)\)⊤​\(−Δ​t​W¯¯\)\\displaystyle\(\\overline\{W\}\+\\overline\{\\overline\{W\}\}\)^\{\\top\}\(\-\\overline\{\\overline\{W\}\}\)\-\(\\frac\{1\}\{\\Delta t\}\(\\overline\{W\}\+\\overline\{\\overline\{W\}\}\-\\mathbb\{I\}\)\)^\{\\top\}\(\-\{\\Delta t\}\\overline\{\\overline\{W\}\}\)=−\(W¯⊤\+W¯¯⊤\)​W¯¯\+\(W¯⊤\+W¯¯⊤−𝕀\)​W¯¯\\displaystyle=\-\(\\overline\{W\}^\{\\top\}\+\\overline\{\\overline\{W\}\}^\{\\top\}\)\\overline\{\\overline\{W\}\}\+\(\\overline\{W\}^\{\\top\}\+\\overline\{\\overline\{W\}\}^\{\\top\}\-\\mathbb\{I\}\)\\overline\{\\overline\{W\}\}=\(−W¯⊤−W¯¯⊤\+W¯⊤\+W¯¯⊤−𝕀\)​W¯¯\\displaystyle=\(\-\\overline\{W\}^\{\\top\}\-\\overline\{\\overline\{W\}\}^\{\\top\}\+\\overline\{W\}^\{\\top\}\+\\overline\{\\overline\{W\}\}^\{\\top\}\-\\mathbb\{I\}\)\\overline\{\\overline\{W\}\}=−W¯¯,\\displaystyle=\-\\overline\{\\overline\{W\}\},implyingW¯¯=−𝕀\\overline\{\\overline\{W\}\}=\-\\mathbb\{I\}\.

Condition 2\. UsingW¯¯=−𝕀\\overline\{\\overline\{W\}\}=\-\\mathbb\{I\}from condition 1, the block matrices simplify toA=W¯−𝕀A=\\overline\{W\}\-\\mathbb\{I\}andC=1Δ​t​\(W¯−2​𝕀\)C=\\frac\{1\}\{\\Delta t\}\(\\overline\{W\}\-2\\mathbb\{I\}\), yielding the following expression for this condition:

\(W¯−𝕀\)⊤​1Δ​t​\(W¯−2​𝕀\)\\displaystyle\(\\overline\{W\}\-\\mathbb\{I\}\)^\{\\top\}\\frac\{1\}\{\\Delta t\}\(\\overline\{W\}\-2\\mathbb\{I\}\)−1Δ​t​\(W¯−2​𝕀\)⊤​\(W¯−𝕀\)\\displaystyle\-\\frac\{1\}\{\\Delta t\}\(\\overline\{W\}\-2\\mathbb\{I\}\)^\{\\top\}\(\\overline\{W\}\-\\mathbb\{I\}\)=1Δ​t​\[\(W¯⊤−𝕀\)​\(W¯−2​𝕀\)−\(W¯⊤−2​𝕀\)​\(W¯−𝕀\)\]\\displaystyle=\\frac\{1\}\{\\Delta t\}\\left\[\(\\overline\{W\}^\{\\top\}\-\\mathbb\{I\}\)\(\\overline\{W\}\-2\\mathbb\{I\}\)\-\(\\overline\{W\}^\{\\top\}\-2\\mathbb\{I\}\)\(\\overline\{W\}\-\\mathbb\{I\}\)\\right\]=1Δ​t​\[−W¯⊤\+W¯\]\\displaystyle=\\frac\{1\}\{\\Delta t\}\[\-\\overline\{W\}^\{\\top\}\+\\overline\{W\}\]=0,\\displaystyle=0,implyingW¯⊤=W¯\\overline\{W\}^\{\\top\}=\\overline\{W\}\.

Finally, conditions 3 and 4 are automatically satisfied by these constraints \(with condition 4 being the transpose counterpart of condition 1\)\. ∎

### B\.1Exact Backward Reconstruction of the SpatialRNN States

We can now formalize the structural conditions under which the SpatialRNN admitsO⁡\(1\)O\(1\)memory cost during training\. In contrast to standard RNNs, where the intermediate states must be stored explicitly for backpropagation, the SpatialRNN enables exact reconstruction of its full trajectory from the final state\. This property follows from the second\-order dynamics of the medium, which induce a reversible state evolution\. When the medium is non\-dissipative \(W¯¯=−𝕀\\overline\{\\overline\{W\}\}=\-\\mathbb\{I\}andW¯=W¯⊤\\overline\{W\}=\\overline\{W\}^\{\\top\}\), the forward dynamics preserve a symplectic structure \(Proposition[B\.1](https://arxiv.org/html/2610.10690#A2.Thmtheorem1)\), implying non\-dissipative information flow and time\-reversibility\. In particular, for blocks withαi=1\\alpha\_\{i\}=1, no information is lost over time, while forαi<1\\alpha\_\{i\}<1the invertibility of the medium dynamics still permits, in exact arithmetic, the recovery of past states despite dissipative behavior\. As a result, intermediate states can be reconstructed by evolving the system backward in time, eliminating the need to cache activations during the forward pass\. The following theorem explicitly constructs this backward mapping\.

###### Theorem B\.3\(Exact history reconstruction\)\.

Suppose the weight matrices satisfy the conditions of Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)\. Then, given the input sequence\{xt\}t=1T\\\{x\_\{t\}\\\}\_\{t=1\}^\{T\}, the complete trajectory\{ψt,ht\}t=1T\\\{\\psi\_\{t\},h\_\{t\}\\\}\_\{t=1\}^\{T\}can be exactly reconstructed by iteratively recovering past states starting from the final boundary states\(ψT,ψT−1\)\(\\psi\_\{T\},\\psi\_\{T\-1\}\)\.

###### Proof\.

Assume that during the backward pass at steptt, the statesψt\\psi\_\{t\}andψt−1\\psi\_\{t\-1\}are available in memory, along with the given input sequence\{xt\}t=1T\\\{x\_\{t\}\\\}\_\{t=1\}^\{T\}\. We show that this information suffices to uniquely reconstruct the preceding states\. First, we note that the hidden state at timet−1t\-1is perfectly reconstructed from the forward update equation:

ht−1=Φ⁡\(Wx​xt−1\+Wψ​ψt−1\+bh\)h\_\{t\-1\}=\\Phi\(W\_\{x\}x\_\{t\-1\}\+W\_\{\\psi\}\\psi\_\{t\-1\}\+b\_\{h\}\)\(56\)
Next, isolatingψt−2\\psi\_\{t\-2\}in the medium dynamics yields

W¯¯​ψt−2=ψt−W¯​ψt−1−B​ht−1\\overline\{\\overline\{W\}\}\\psi\_\{t\-2\}=\\psi\_\{t\}\-\\overline\{W\}\\psi\_\{t\-1\}\-Bh\_\{t\-1\}\(57\)Under the assumption of Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3), we have that∃W¯¯−1\\exists\\ \\overline\{\\overline\{W\}\}^\{\-1\}\. Applying its inverse gives a unique solution for the prior medium state:

ψt−2=W¯¯−1​\(ψt−W¯​ψt−1−B​ht−1\)\.\\psi\_\{t\-2\}=\\overline\{\\overline\{W\}\}^\{\-1\}\\left\(\\psi\_\{t\}\-\\overline\{W\}\\psi\_\{t\-1\}\-Bh\_\{t\-1\}\\right\)\.\(58\)This establishes that the pair\(ψt,ψt−1\)\(\\psi\_\{t\},\\psi\_\{t\-1\}\)uniquely determines\(ψt−1,ψt−2\)\(\\psi\_\{t\-1\},\\psi\_\{t\-2\}\), and hence the state at timet−2t\-2\. By recursively applying this procedure fromt=Tt=Tdown tot=2t=2, we recover the complete forward sequence\{ψt,ht\}t=1T\\\{\\psi\_\{t\},h\_\{t\}\\\}\_\{t=1\}^\{T\}, without requiring storage of intermediate states\. ∎

### B\.2Numerical Precision of the Reconstruction

Theorem[B\.3](https://arxiv.org/html/2610.10690#A2.Thmtheorem3)holds in exact arithmetic\. In finite precision, round\-off is propagated by the backward recursion, whose conditioning is governed by the inverse medium dynamics\. Under the conditions of Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)the modes of the medium have modulusαi\\sqrt\{\\alpha\_\{i\}\}, so the backward recursion, which inverts them, amplifies round\-off byαi−1/2\\alpha\_\{i\}^\{\-1/2\}per step\. For a toleranceεtol\\varepsilon\_\{\\mathrm\{tol\}\}and machine epsilonεmach\\varepsilon\_\{\\mathrm\{mach\}\}, the useful horizon is therefore

mmax≈log⁡\(εtol/εmach\)log\(αmin−1/2\),αmin=mini⁡αi\.m\_\{\\max\}\\;\\approx\\;\\frac\{\\log\\\!\\big\(\\varepsilon\_\{\\mathrm\{tol\}\}/\\varepsilon\_\{\\mathrm\{mach\}\}\\big\)\}\{\\log\\\!\\big\(\\alpha\_\{\\min\}^\{\-1/2\}\\big\)\},\\qquad\\alpha\_\{\\min\}=\\min\_\{i\}\\alpha\_\{i\}\.\(59\)This is a property of the recursion rather than of the data, so we isolate it in a controlled setting: we compare the JacobiansGm=∂hT−1/∂hT−1−mG\_\{m\}=\\partial h\_\{T\-1\}/\\partial h\_\{T\-1\-m\}obtained from the stored forward states with those obtained from the backward\-reconstructed states, for randomly initialized models driven by2020synthetic scalar input sequences of lengthT=784T=784\(Fig\.[8](https://arxiv.org/html/2610.10690#A2.F8)\)\. In the non\-dissipative case \(αi=1​∀i\\alpha\_\{i\}=1\\;\\forall i\) the recursion has all its modes on the unit circle, no amplification occurs, and the relative error stays below10−510^\{\-5\}in single and10−1310^\{\-13\}in double precision over the entire horizon\. With dissipative blocks \(α=\[1,0\.5,0\.5,0\.5\]\\alpha=\[1,0\.5,0\.5,0\.5\]\) the error grows exponentially at a measured rate of1\.421\.42per step, matching the1/0\.5=1\.4141/\\sqrt\{0\.5\}=1\.414predicted by the modal analysis, and exceeds10−310^\{\-3\}after4444steps in single and101101steps in double precision, in agreement with \([59](https://arxiv.org/html/2610.10690#A2.E59)\)\. BecauseW¯\\overline\{W\}andW¯¯\\overline\{\\overline\{W\}\}are block*upper*\-triangular, the marginally stable block is not shielded: it reads from the dissipative blocks and inherits their divergence\. Consequently,𝒪⁡\(1\)\\mathcal\{O\}\(1\)\-memory BPTT is practical only when all blocks are non\-dissipative; models withαi<1\\alpha\_\{i\}<1require caching, or checkpointing everyK<mmaxK<m\_\{\\max\}steps, which brings the memory use to𝒪⁡\(T/K\)\\mathcal\{O\}\(T/K\), and could be used for parallelization\.

\(a\) Single block,α1=1\\alpha\_\{1\}=1 ![Refer to caption](https://arxiv.org/html/2610.10690v1/figures/reversibility_single.jpeg)

\(b\) Four blocks,α=\[1,0\.5,0\.5,0\.5\]\\alpha=\[1,0\.5,0\.5,0\.5\] ![Refer to caption](https://arxiv.org/html/2610.10690v1/figures/reversibility_multi.jpeg)

Figure 8:Relative Frobenius error‖Gmexact−Gmrec‖F/‖Gmexact‖F\\\|G\_\{m\}^\{\\mathrm\{exact\}\}\-G\_\{m\}^\{\\mathrm\{rec\}\}\\\|\_\{F\}/\\\|G\_\{m\}^\{\\mathrm\{exact\}\}\\\|\_\{F\}between the exact JacobiansGm=∂hT−1/∂hT−1−mG\_\{m\}=\\partial h\_\{T\-1\}/\\partial h\_\{T\-1\-m\}, computed from the stored forward states, and those computed from the backward\-reconstructed states, as a function of the time lagmm, in single \(float32\) and double \(float64\) precision\.\(a\)Non\-dissipative SpatialRNN \(single block,n=40n=40,α1=1\\alpha\_\{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\}\\\{22,6,6,6\\\},α=\[1,0\.5,0\.5,0\.5\]\\alpha=\[1,0\.5,0\.5,0\.5\]\): round\-off is amplified exponentially, at a measured rate of1\.421\.42per backward step, matching the predictedαmin−1/2=1\.414\\alpha\_\{\\min\}^\{\-1/2\}=1\.414\. Vertical lines mark the lags at which the median error first exceeds10−310^\{\-3\}, yieldingm=44m=44in single andm=101m=101in double precision\. Since the reconstruction starts from\(ψT−1,ψT−2\)\(\\psi\_\{T\-1\},\\psi\_\{T\-2\}\)and proceeds backward, the state at timettis recovered afterT−1−tT\-1\-tbackward steps, sommis also the depth of the backward recursion\. Curves show the median over2020input sequences of lengthT=784T=784and shaded bands the1010–9090percentile 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](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)holds\.

## Appendix CExperimental Details

In this section, we detail the datasets, training procedures, and benchmarking used in Section[4](https://arxiv.org/html/2610.10690#S4), along with representative data samples\. All code for reproduction is available at: \[[anonymous\.4open\.science/r/SpatialRNN\-code\-748C](https://anonymous.4open.science/r/SpatialRNN-code-748C)\]

### C\.1SpatialRNN implementation details

Activation function and decay parameters\.To implement the SpatialRNN in practice, we set the activation function toΦ=tanh\\Phi=\\tanh\. In all experiments, we enforce the structural and spectral conditions outlined in Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3)to guarantee the retention of long\-term dependencies\. We thus setα1=1\\alpha\_\{1\}=1\(typically corresponding to the largest block\), while treating the remaining decay parameters\{α2,α3,…\}\\\{\\alpha\_\{2\},\\alpha\_\{3\},\\ldots\\\}as trainable\.

Masking and initialization\.To satisfy the block upper\-triangular constraints of Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3), we apply static binary masks during each forward pass\. We use block\-diagonal masks forBB, strictly diagonal mask forΣ\\Sigma, setWψ=Σ​B⊤W\_\{\\psi\}=\\Sigma B^\{\\top\}, and use block triangular masks for bothW¯\\bar\{W\}andW¯¯\\bar\{\\bar\{W\}\}\.

The input and output projection matrices \(Wx,WyW\_\{x\},W\_\{y\}\) are initialized using Xavier uniform distribution\. For the medium dynamics, the block\-diagonal components ofW¯\\bar\{W\}are initialized using orthogonal initialization \(gain of0\.50\.5\), whileBBandΣ\\Sigmaare initialized as the identity matrix\. The upper\-triangular interaction blocks are initialized with Xavier uniform distribution with a reduced gain of0\.10\.1to not dominate the intra\-block dynamics\.

Spectral Bounding\.To strictly satisfy the stability bounds at every iteration, we dynamically project the weight matrices using their spectral norms\. Following Theorem[3\.3](https://arxiv.org/html/2610.10690#S3.Thmtheorem3), we apply a multiplicative scaling factorγi\\gamma\_\{i\}to the blocks to constrain the gradient spectrum:

γi=min⁡\(1,2​αi−ϵ‖W¯i‖2\+‖Bi‖2​‖Wψ,i‖2\),\\gamma\_\{i\}=\\min\\left\(1,\\frac\{2\\sqrt\{\\alpha\_\{i\}\}\-\\epsilon\}\{\\\|\\bar\{W\}\_\{i\}\\\|\_\{2\}\+\\\|B\_\{i\}\\\|\_\{2\}\\\|W\_\{\\psi,i\}\\\|\_\{2\}\}\\right\),\(60\)whereϵ=10−4\\epsilon=10^\{\-4\}ensures the inequality remains strict\. Alternatively, we can also use a decoupled bounding scheme where the norms are independently constrained via separate scaling factors \(e\.g\.,‖W¯i‖2≤cW\\\|\\bar\{W\}\_\{i\}\\\|\_\{2\}\\leq c\_\{W\}and‖Bi‖2,‖Wψ,i‖2≤cB\\\|B\_\{i\}\\\|\_\{2\},\\\|W\_\{\\psi,i\}\\\|\_\{2\}\\leq c\_\{B\}\)\. Such constraints tend to yield improved performance, as it more directly enforces—or heuristically approximates—the LTV stability conditions derived in Theorem[A\.3](https://arxiv.org/html/2610.10690#A1.Thmtheorem3)\.

### C\.2HORNN implementation details

Thedd\-order HORNN with delay matrices\{Wk\}k=1d\\\{W\_\{k\}\\\}\_\{k=1\}^\{d\}\(as presented in Eq\. \([4](https://arxiv.org/html/2610.10690#S2.E4)\)\) can be generalized to use a sparse set of delays\{Wk\}𝒯\\\{W\_\{k\}\\\}\_\{\\mathcal\{T\}\}for𝒯⊂\{1,…,d\}\\mathcal\{T\}\\subset\\\{1,\\ldots,d\\\}\. This generalization leads to the update equation:

ht=Φ⁡\(Wx​xt\+∑k∈𝒯Wk​ht−k\+bh\)\.h\_\{t\}=\\Phi\\left\(W\_\{x\}x\_\{t\}\+\\sum\_\{k\\in\\mathcal\{T\}\}W\_\{k\}h\_\{t\-k\}\+b\_\{h\}\\right\)\.\(61\)Such model can more efficiently evaluate long\-range memory retention without the parameter overhead of dense delays\. The effective temporal receptive field of the model is bounded bydmax=max⁡\(𝒯\)d\_\{\\max\}=\\max\(\\mathcal\{T\}\)\.

Implementation and initialization\.Consistent with the SpatialRNN architecture, we takeΦ=tanh\\Phi=\\tanh\. To ensure stable variance scaling across the sparse temporal connections, the input projectionWxW\_\{x\}, output projectionWyW\_\{y\}, and all sparse recurrence matrices\{Wk\}k∈𝒯\\\{W\_\{k\}\\\}\_\{k\\in\\mathcal\{T\}\}are initialized using a Xavier uniform distribution\. The hidden state biasbhb\_\{h\}and output biasbyb\_\{y\}are strictly initialized to zero\.

### C\.3Copy Task

Dataset generation\.The copy task evaluates the model’s ability to perfectly retain information over extended temporal horizons\. We define a token dimensionD=4D=4and a vocabulary ofV=2D=16V=2^\{D\}=16tokens\. The data tokens are embedded as vectors in\{−1,1\}D\\\{\-1,1\\\}^\{D\}\.

The input sequenceXXof total lengthTTis constructed as a concatenation of three phases:X=\[Xdata,Xdist,Xtrig\]X=\[X\_\{\\text\{data\}\},X\_\{\\text\{dist\}\},X\_\{\\text\{trig\}\}\], whereXdata∈\{−1,1\}D×LX\_\{\\text\{data\}\}\\in\\\{\-1,1\\\}^\{D\\times L\}consists ofL=5L=5tokens sampled uniformly at random from the data tokens,Xdist∈\{0\}D×TdistX\_\{\\text\{dist\}\}\\in\\\{0\\\}^\{D\\times T\_\{\\text\{dist\}\}\}consists ofTdistT\_\{\\text\{dist\}\}distractor tokens padded with zeros, andXtrig∈\{1\}D×LX\_\{\\text\{trig\}\}\\in\\\{1\\\}^\{D\\times L\}consists ofL=5L=5trigger tokens padded with ones\.

Training and evaluation\.Models are trained to minimize the cross\-entropy loss over the finalLLretrieval steps, comparing the generated outputs \(after the whole distractor sequence has been fed into the model\) to the initialLLdata tokensXdataX\_\{\\text\{data\}\}\. The outputs of the models for the firstL\+TdistL\+T\_\{\\text\{dist\}\}are discarded\. At each epoch, we generateNtrain=3000N\_\{\\text\{train\}\}=3000training sequences,Ntest=800N\_\{\\text\{test\}\}=800testing andNval=800N\_\{\\text\{val\}\}=800validation sequences\. We take a batch size of6464, and use the Adam optimizer with a learning rate of0\.0020\.002\. Training proceeds for a maximum of30003000epochs, using an early stopping patience of300300epochs\.

Model details\.For the SpatialRNN models, the temporal decay parameters are initialized toαi=0\.8\\alpha\_\{i\}=0\.8for all blocksii, except for the first block which is strictly fixed toα1=1\.0\\alpha\_\{1\}=1\.0throughout training to ensure at least one non\-dissipative, marginally stable memory stream\. We performed a hyperparameter search over various block size configurations\. Table[4](https://arxiv.org/html/2610.10690#A3.T4)details the specific sizes evaluated, their corresponding parameter counts, and the resulting accuracy on the fullT=1024T=1024sequence length\. We bound the matrix norms by‖W¯i‖2≤1\.7\\\|\\bar\{W\}\_\{i\}\\\|\_\{2\}\\leq 1\.7and‖Bi‖2,‖Wψ,i‖2≤0\.2\\\|B\_\{i\}\\\|\_\{2\},\\\|W\_\{\\psi,i\}\\\|\_\{2\}\\leq 0\.2\.

Table 4:Copy task architecture configurations and performance\. Mean accuracy is evaluated at the maximum sequence lengthT=1024T=1024\. Average and standard deviation over 10 runs\.
### C\.4Permuted Sequential MNIST \(psMNIST\)

Dataset and preprocessing\.The psMNIST task takes MNIST images and flattens them into sequences of lengthT=784T=784with feature dimensionDin=1D\_\{\\text\{in\}\}=1\. To ablate local spatial structure, a fixed random permutationπ∈S784\\pi\\in S\_\{784\}\(seed 42\) is applied to the time steps of all sequences, isolating the task’s reliance on long\-range temporal dependencies\. The vector values are then normalized viax↦x−μσx\\mapsto\\frac\{x\-\\mu\}\{\\sigma\}\(μ=0\.1307,σ=0\.3081\\mu=0\.1307,\\sigma=0\.3081\)\.

Training and hyperparameters\.The images are fed to the model sequentially, and after all the sequence has been processed, a single integer \(encoded as the maximum one\-hot vector\)Dout=10D\_\{\\text\{out\}\}=10is generated\. Models are trained to minimize cross\-entropy loss overNruns=3N\_\{\\text\{runs\}\}=3independent initializations\. Optimization utilizes Adam with a learning rate of10−410^\{\-4\}and batch size of128128\. Training proceeds for a maximum of10001000epochs with an early stopping patience of300300epochs\. We use the standardNtrain=60,000N\_\{\\text\{train\}\}=60,000training images andNtest=10,000N\_\{\\text\{test\}\}=10,000testing images\.

Temporal decay parameters are initialized toαi=0\.8\\alpha\_\{i\}=0\.8for all blocksii, except for the first block which is strictly fixed toα1=1\.0\\alpha\_\{1\}=1\.0\. We bound the matrix norms by‖W¯i‖2≤1\.7\\\|\\bar\{W\}\_\{i\}\\\|\_\{2\}\\leq 1\.7and‖Bi‖2,‖Wψ,i‖2≤0\.2\\\|B\_\{i\}\\\|\_\{2\},\\\|W\_\{\\psi,i\}\\\|\_\{2\}\\leq 0\.2and clip the global gradient norm to 1\.0\. To produce the results for Table[1](https://arxiv.org/html/2610.10690#S4.T1), we performed a simple hyperparameter search over various block size configurations\. Table[5](https://arxiv.org/html/2610.10690#A3.T5)details the specific sizes evaluated, their corresponding parameter counts, and the resulting mean accuracy across 3 independent runs, ordered by performance\.

Table 5:psMNIST task architecture configurations and performance\. Mean accuracy is evaluated over 3 runs, ordered by mean accuracy\.Compared models\.We evaluate our model against state\-of\-the\-art RNN architectures designed for long\-range dependencies\. When previous methods \(e\.g\., NWM, coRNN\) have been evaluated multiple times in the literature, we prioritize results derived from multi\-run aggregations to ensure statistical robustness\. Additionally, to guarantee a fair comparison, we restrict all baselines to single\-layer models trained under standard optimization protocols\. By omitting complex learning rate schedules and regularizers, such as dropout\([Srivastava et al\., 2014](https://arxiv.org/html/2610.10690#bib.bib48)\), batch normalization\([Ioffe and Szegedy, 2015](https://arxiv.org/html/2610.10690#bib.bib49)\)or step decay\([He et al\., 2016](https://arxiv.org/html/2610.10690#bib.bib53)\), we isolate the intrinsic memory capacity of the recurrent mechanism from the performance gains typically afforded by increased network depth or training enhancements\.

### C\.5Noise Padded CIFAR\-10 \(npCIFAR10\)

Dataset and preprocessing\.The npCIFAR10 task evaluates the model’s ability to extract spatial features and retain them over prolonged noisy horizons\. We use the standard CIFAR\-10 dataset \(50,00050,000training and10,00010,000testing images on 10 classes\)\. Each image is normalized using the channel\-wise means\(0\.4914,0\.4822,0\.4465\)\(0\.4914,0\.4822,0\.4465\)and deviations\(0\.2023,0\.1994,0\.2010\)\(0\.2023,0\.1994,0\.2010\), so the entries are scaled approximately to\[−1,1\]\[\-1,1\]\.

The normalized3×32×323\\times 32\\times 32images are then permuted and reshaped into sequences of lengthL=32L=32with feature dimensionDi​n=32×3=96D\_\{in\}=32\\times 3=96\. In this representation, each timestept∈\{1,…,32\}t\\in\\\{1,\\dots,32\\\}corresponds to an image row across all color channels\. To introduce long\-term memory requirement, we append random vectors following a standard normal distribution𝒩⁡\(0,I96\)\\mathcal\{N\}\(0,I\_\{96\}\)for the remainingT−32T\-32timesteps, yielding a final sequence lengthTT\.

Training and hyperparameters\.The models’ firstT−1T\-1outputs are ignored, and the final outputyT∈\{0,1\}10y\_\{T\}\\in\\\{0,1\\\}^\{10\}\(after processing the complete input sequence of lengthTT\) is compared to the image class using the standard cross\-entropy loss\. While the primary benchmark results \(Table[2](https://arxiv.org/html/2610.10690#S4.T2)\) evaluate the standardT=1000T=1000regime, we also perform an ablation across varying sequence lengthsT∈\{50,100,200,500,1000\}T\\in\\\{50,100,200,500,1000\\\}to show that the amount of noise distractors does not affect the task performance, but the limitation comes from the classification performance of the model\. Models are trained using the Adam optimizer with a learning rate of10−410^\{\-4\}, global gradient clipping of 1\.0 and a batch size of256256\. Training is executed for a maximum of30003000epochs acrossNr​u​n​s=3N\_\{runs\}=3independent initializations\. We employ an early stopping patience of300300epochs based on validation performance\.

Model details\.For the SpatialRNN models, the temporal decay parameters are initialized toαi=0\.8\\alpha\_\{i\}=0\.8for all blocksii, except for the first block which is strictly fixed toα1=1\.0\\alpha\_\{1\}=1\.0throughout training to ensure at least one non\-dissipative, marginally stable memory stream\. The specific model configurations used for Table[2](https://arxiv.org/html/2610.10690#S4.T2)consist of block sizes\{128,8,8,8\}\\\{128,8,8,8\\\}for the SpatialRNN of 152 neurons \(51k parameters\), and block sizes\{64,8,8,8\}\\\{64,8,8,8\\\}for the SpatialRNN of 88 neurons, resulting in 21k parameters\. We performed a hyperparameter search over various block size configurations\. Table[6](https://arxiv.org/html/2610.10690#A3.T6)details the specific sizes evaluated, their corresponding parameter counts, and the resulting maximum accuracy \(over 3 runs\) on the fullT=1000T=1000sequence length\. We bound the matrix norms by‖W¯i‖2≤1\.7\\\|\\bar\{W\}\_\{i\}\\\|\_\{2\}\\leq 1\.7and‖Bi‖2,‖Wψ,i‖2≤0\.2\\\|B\_\{i\}\\\|\_\{2\},\\\|W\_\{\\psi,i\}\\\|\_\{2\}\\leq 0\.2\.

Table 6:npCIFAR10 task architecture configurations and performance\. Maximum accuracy is evaluated at the sequence lengthT=1000T=1000over 3 runs\.Compared models\.To appropriately benchmark the SpatialRNN, we select a strong suite of competitors previously evaluated on the npCIFAR10 task\. These include traditional gated architectures \(LSTM, GRU\) to demonstrate the breakdown of standard gating mechanisms over long sequences, as well as state\-of\-the\-art continuous\-time and oscillatory models \(e\.g\., coRNN, LEM, Lipschitz RNN\) that specifically target the vanishing gradient problem\. To ensure a rigorous and equivalent comparison of intrinsic memory capacity, all evaluated baseline models are restricted to single\-layer architectures\.

### C\.6Hardware and Compute

All experiments were implemented in PyTorch and executed on a private workstation equipped with two AMD EPYC 9534 64\-core processors \(total of 128 physical CPU cores and 256 threads\) and 755 GB RAM\. Model training and evaluation were performed on CPU resources only; no GPU acceleration or multi\-GPU model parallelism was used throughout\. The reported experiments are computationally lightweight relative to large\-scale sequence\-modeling workloads\. The copy\-task experiments used synthetically generated sequences and small recurrent models, with the largest reported SpatialRNN architecture containing 6,815 parameters and the largest HORNN architecture containing 17,744 parameters\. The psMNIST \(resp\. npCIFAR\) experiments used single\-layer SpatialRNNs with up to 36,497 \(resp\. 50,937\) parameters, trained for at most 1000 \(resp\. 3000\) epochs with early stopping over three independent initializations\. The copy\-task experiments accounted for approximately 1,600 CPU\-thread\-hours, including sweeps over sequence length and recurrent baselines\. The psMNIST experiment required approximately 1,600 CPU\-thread\-hours across the three reported initializations and model\-size configurations\. The npCIFAR10 experiments required approximately 3,200 CPU\-thread\-hours across the three reported initializations and sequence\-length ablations\. Overall, using up to 20 hardware threads for each experiment type, we estimate that the total compute required for the final reported experiments was approximately two weeks on the reported CPU workstation\.