Paper deep dive
Divisive Normalization Shapes Low-Rank Slow Manifolds for Continuous Working Memory
Zhaotian Gu, Jie Su, Weiwei Wang, Chang Liu, Tianyi Qian, Dahui Wang
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The ability to robustly maintain and update continuous variables is a hallmark of working memory. While classical continuous attractor networks suffer from severe fine-tuning fragility, standard artificial recurrent neural networks (RNNs) like GRUs and LSTMs typically fail to stably learn continuous manifolds, instead shattering the state space into discretized point attractors. To bridge this gap, we draw inspiration from divisive normalization, a canonical neural computation widely observed across cortical circuits, and propose the Recurrent Divisive Normalization Network (RDNN), a minimal and algebraically isolated model of dynamic division. Through dynamical systems analysis on canonical working memory tasks, we demonstrate that this biophysical constraint allows the network to converge to robust, high-fidelity slow manifolds. Furthermore, we analyze the gradient dynamics of divisive normalization during Backpropagation Through Time (BPTT), showing that it introduces an activity-dependent local gradient scaling. This scaling dampens parameter updates in highly active regimes, which empirically aligns with a significant self-compression of the network's effective rank, confining the recurrent dynamics to a tight, low-dimensional subspace while avoiding the optimization pathologies associated with explicit low-rank factorization. Finally, ablations demonstrate that while subtractive inhibition can maintain static memories, divisive normalization is mathematically essential to prevent manifold shattering under time-varying inputs. Our findings identify divisive normalization not merely as a biological artifact, but as a critical computational mechanism for learning high-fidelity continuous representations.
Tags
Links
- Source: https://arxiv.org/abs/2608.01947v1
- Canonical: https://arxiv.org/abs/2608.01947v1
Trouble viewing inline? Open PDF directly →
Full Text
134,321 characters extracted from source content.
Expand or collapse full text
Divisive Normalization Shapes Low-Rank Slow Manifolds for Continuous Working Memory Zhaotian Gu1† *, Jie Su3† , Weiwei Wang2,3, Chang Liu1, Tianyi Qian3, Dahui Wang1,2* 1 School of System Science, Beijing Normal University 2 State Key Laboratory of Cognitive Neuroscience and Learning, Beijing Normal University 3 Qiyuan Laboratory zhaotiangu, weiwei.wang, liu_chang@mail.bnu.edu.cn wangdh@bnu.edu.cn sujie, qiantianyi@qiyuanlab.com † Equal Contribution * Corresponding author Abstract The ability to robustly maintain and update continuous variables is a hallmark of working memory. While classical continuous attractor networks suffer from severe fine-tuning fragility, standard artificial recurrent neural networks (RNNs) like GRUs and LSTMs typically fail to stably learn continuous manifolds, instead shattering the state space into discretized point attractors. To bridge this gap, we draw inspiration from divisive normalization, a canonical neural computation widely observed across cortical circuits, and propose the Recurrent Divisive Normalization Network (RDNN), a minimal and algebraically isolated model of dynamic division. Through dynamical systems analysis on canonical working memory tasks, we demonstrate that this biophysical constraint allows the network to converge to robust, high-fidelity slow manifolds. Furthermore, we analyze the gradient dynamics of divisive normalization during Backpropagation Through Time (BPTT), showing that it introduces an activity-dependent local gradient scaling. This scaling dampens parameter updates in highly active regimes, which empirically aligns with a significant self-compression of the network’s effective rank, confining the recurrent dynamics to a tight, low-dimensional subspace while avoiding the optimization pathologies associated with explicit low-rank factorization. Finally, ablations demonstrate that while subtractive inhibition can maintain static memories, divisive normalization is mathematically essential to prevent manifold shattering under time-varying inputs. Our findings identify divisive normalization not merely as a biological artifact, but as a critical computational mechanism for learning high-fidelity continuous representations. 1 Introduction The ability to maintain and manipulate continuous variables over time, such as spatial location, head direction, or accumulated evidence, is a hallmark of biological working memory. In computational neuroscience, this cognitive function is classically modeled using Continuous Attractor Networks (CANs) (Wimmer et al., 2014; Khona and Fiete, 2022). However, CANs are notoriously brittle: they suffer from the fine-tuning problem, where infinitesimal perturbations in recurrent dynamics can destroy the continuous manifold of equilibria, causing the stored memory to rapidly drift or collapse (Seung, 1996; Koulakov et al., 2002). In parallel, while standard artificial Recurrent Neural Networks (RNNs) like LSTMs and GRUs excel at sequential tasks, recent dynamical systems analyses reveal that they struggle to stably learn true continuous attractors. Instead, they often discretize the continuous state space, shattering the manifold into localized point attractors separated by saddles (Jordan et al., 2021). This raises a fundamental question: What computational mechanisms allow biological neural networks to robustly represent continuous working memory, and how can we integrate them into learnable computational models? To bridge this gap, we draw inspiration from divisive normalization (DN), a canonical neural computation widely observed across cortical circuits for dynamic gain control and contrast invariance (Carandini and Heeger, 2012). While previous works have incorporated DN into recurrent neural circuits to study biophysical properties and Lyapunov stability, most notably in ORGaNICs (Heeger and Mackey, 2019; Rawat et al., 2024), these models often feature complex multi-component dynamics that obscure the specific, isolated role of division under optimization. Furthermore, how DN shapes the gradient flow during training, the topological landscape of learned slow manifolds, and the representation dimensionality remains unexplored. In this work, we propose the Recurrent Divisive Normalization Network (RDNN), a minimal, algebraically isolated model of dynamic division designed to address these questions. Rather than optimizing for raw task performance, our objective is to systematically investigate how this biophysical mechanism structurally shapes learned continuous working memory representations. By training the RDNN on canonical continuous working memory tasks, including autonomous maintenance (memory-guided saccade) and input-driven updating (angular velocity integration), we investigate how this structural constraint shapes the learned neural representations. Through dynamical systems analysis, we uncover that divisive normalization is not merely a biological artifact, but a critical computational mechanism for high-fidelity continuous memory. Furthermore, this mechanism induces a significant self-compression of the network’s effective rank, forcing the recurrent dynamics into a tight, low-dimensional subspace without requiring explicit low-rank parameterization. Our main contributions are summarized as follows: • Biologically-inspired Architecture for Continuous Memory: We introduce the RDNN and demonstrate that it learns robust slow manifolds for continuous working memory, overcoming the topological discretization inherent in standard GRU and LSTM architectures. • Emergent Low-Rank Dynamics: We analyze the gradient dynamics of divisive normalization under BPTT, showing that the dynamic divisor attenuates gradient updates in high-activity regimes. This local scaling restricts major parameter updates to a subset of states, and we show empirically that RDNNs converge to representations with significantly lower effective ranks than baselines, and avoids the non-convex optimization anomalies associated with explicit low-rank factorization. • Mechanistic Ablation of Inhibition: Through systematic ablation, we delineate the functional niches of divisive versus subtractive inhibition. We demonstrate that while subtractive inhibition can maintain static memories, multiplicative gain control is essential to prevent manifold shattering under time-varying external inputs. 2 Related Works Divisive Normalization In neuroscience, divisive normalization is a canonical cortical operation that implements gain control and robust encoding. Classical studies (Heeger, 1992; Carandini and Heeger, 2012) demonstrated that V1 neurons divide their inputs by a pooled activity, explaining nonlinear contrast effects. Recent work has quantitatively fit DN models to neural data (Sawada and Petrov, 2017) and even learned orientation-specific normalization kernels from visual cortex (Burg et al., 2021). In parallel, machine learning uses related normalization layers (Ioffe and Szegedy, 2015; Xu et al., 2019), but the neuroscience view emphasizes explicit dynamic divisive feedback. These findings underscore that normalization circuits filter out noise and shape responses before further processing. Several recent models incorporate DN into recurrent memory circuits. Heeger and Mackey (2019) introduced ORGaNICs (Oscillatory Recurrent Gated Neural Integrator Circuits), a biologically plausible RNN framework with dynamic DN that unifies sensory gain control, gated integration, and working-memory maintenance. Subsequently, Rawat et al. (2024) analyzed these circuits to prove their unconditional Lyapunov stability under broad conditions, showing they can be trained with backpropagation without special tricks and perform comparably to LSTMs on sequential tasks. However, both ORGaNICs and its subsequent analyses focus primarily on hand-crafted biophysical dynamics and analytical stability under constant inputs, without addressing the optimization dynamics or isolating the algebraic impact of the divisive operation itself. While ORGaNICs integrates multiple biophysical components, the complex interactions between these mechanisms make it conceptually difficult to isolate the individual contribution of division to learned representations. In contrast, our proposed RDNN serves as a minimal, structurally isolated abstraction of DN, allowing us to demonstrate how pure dynamic division intrinsically guides the network toward low-dimensional continuous representations. Continuous Attractor Networks and Slow Manifolds A classic model for analog memory is the continuous attractor network (Khona and Fiete, 2022; Wimmer et al., 2014), which can store graded values indefinitely. Such models have been used to explain spatial working memory in prefrontal cortex, for example by linking the precision of recall to the width of a bump attractor (Wimmer et al., 2014). However, continuous attractors suffer from a severe structural instability known as the fine-tuning problem (Seung, 1996; Koulakov et al., 2002): even infinitesimal perturbations in recurrent dynamics can destroy the continuous manifold of equilibria, causing the stored representation to rapidly drift or collapse. To address this fragility, rather than enforcing perfect fine-tuning, recent work has shifted toward approximate continuous attractors, i.e. persistent, low-dimensional slow manifolds that survive bifurcations and bound memory drift over behaviorally relevant timescales (Ghazizadeh and Ching, 2021; Ságodi et al., 2024). In particular, Ságodi et al. (2024) established a theoretical foundation showing that physical systems only require a slow manifold close to a continuous attractor to remain functionally robust, effectively resolving the ontological dilemma of ideal CANs. We build directly upon this approximate attractor framework. While Ságodi et al. (2024) analyzed the survival of slow manifolds in standard architectures under perturbations, we investigate how specific biologically motivated mechanisms naturally shape these slow manifolds. Meanwhile, dynamical analyses of standard gated architectures, such as Gated Recurrent Units (GRUs) viewed through continuous-time systems, have highlighted their structural inability to stably learn or represent continuous attractors (Jordan et al., 2021), further underscoring the need for alternative biologically inspired architectures with multiplicative feedback. Low-rank and low-dimensional RNNs A related line of work studies how structured connectivity produces low-dimensional dynamics. Mastrogiuseppe and Ostojic (2018) introduced the low-rank RNN framework: if recurrent weights have a low-dimensional component plus noise, the network’s dynamics are confined to a low-dimensional subspace determined by that structure. They show one can design minimal connectivity to implement particular computations (e.g. integration or oscillation). More recently, Mastrogiuseppe et al. (2025) extended this to stochastic inputs: they prove that low-dimensional inputs keep the activity low-rank, but high-dimensional noise can inflate the network’s activity dimension beyond the nominal rank. Bridging local constraints with global behavior, Shao et al. (2025) demonstrated that specific local synaptic patterns in excitatory-inhibitory networks mathematically induce low-rank structures that dominate the overall recurrent dynamics. In general, tools like fixed-point analysis (Sussillo and Barak, 2013) reveal that trained RNNs operate via a few fixed points and slow trajectories. These insights explain how E-I networks with normalization or other constraints could generate the simple low-dimensional manifolds observed in memory tasks. 3 Proposed Model To investigate how biological constraints shape the neural representations of continuous working memory, we draw inspiration from divisive normalization, a common mechanism for dynamic gain control in the brain. As illustrated in Figure 1A, the proposed Recurrent Divisive Normalization Network (RDNN) consists of two interacting populations: a principal excitatory population ∈ℝ+HR ^H_+ (representing the working memory state) and an auxiliary inhibitory pool ∈ℝ+HG ^H_+ (providing dynamic gain control), where H is the hidden dimension. Figure 1: Overview of the proposed RDNN and its performance on continuous working memory tasks. A The RDNN architecture, featuring an excitatory population R and an inhibitory pool G for divisive gain control (square-headed lines). B Schematics of the continuous working memory tasks: angular velocity integration (B1) and memory-guided saccade (B2), with colors indicating the encoded angular variable. C Test NMSE (dB) across different architectures, activations, and hidden sizes. Error bars denote the standard deviation across different random seeds. The continuous-time dynamics of the RDNN are governed by a system of coupled stochastic differential equations (SDEs). To emulate the intrinsic stochasticity and synaptic noise characteristic of biological neural systems (Faisal et al., 2008), we inject state noise into the dynamics. The dynamics of the network is defined as: τd _RdR =(−+f()+(t)+)dt+σd, = (-R+ Jf(R)+I(t) η+G )dt+ _RdW_R, (1) τd _GdG =(−+f())dt+σd. = (-G+Wf(R) )dt+ _GdW_G. Here, (t)∈ℝHI(t) ^H denotes the external input stimulus at time t; f(⋅)f(·) is an elementwise nonlinear activation function (e.g., ReLU or Tanh) that maps the state variables to firing rates; ∈ℝ+H×HJ _+^H× H is the recurrent excitatory weight matrix (see the arrows between RiR_i in Figure 1A), and ∈ℝ+H×HW _+^H× H is the projection weight from the excitatory to the inhibitory population (the arrows from RiR_i to GiG_i in Figure 1A). The division operation (⋅)(⋅) (·)(·) is performed elementwise and represents the shunting inhibition exerted by G on R (square-headed connections in Figure 1A); ∈ℝ+H η _+^H is a semi-saturation constant vector that prevents division by zero and sets the neurons’ baseline responsiveness; τ,τ∈ℝ+H _R, _G _+^H are the membrane time constants of the two populations; and ddW_R and ddW_G are independent standard Wiener processes scaled by the noise intensities σ _R and σ _G, respectively, which model the state noise. In our model, parameters ,,,τ,τJ,W, η, _R, _G are all learnable and constrained to be positive, ensuring that the network operates as a positive dynamical system consistent with biological principles such as Dale’s law (Fitzpatrick et al., 1987) (see Appendix B). Unlike standard RNNs that rely on gating mechanisms, the RDNN employs a divisive bottleneck. The inhibitory state G dynamically tracks the overall network activity via f()Wf(R). When the recurrent drive f()Jf(R) or external input (t)I(t) fluctuates, G proportionally scales the denominator, thereby stabilizing the network against explosive positive feedback while preserving the relative activation patterns. As we will demonstrate, this structural inductive bias naturally induces low-rank, high-fidelity slow manifolds. 4 Slow Manifold Dynamics in Continuous Working Memory Tasks To understand the computational mechanisms underlying continuous working memory, we trained the proposed RDNN and standard gated models (GRU, LSTM) across various hidden sizes (64, 128, 256) and activation functions (ReLU, Tanh) on two canonical tasks (Figure 1B): Angular Velocity Integration (input-driven updating) and Memory-guided Saccade (autonomous maintenance). While the RDNN achieves comparable or superior task performance to the baselines (Figure 1C), we seek to uncover the fundamental differences in their learned representations. To this end, we systematically analyze their slow manifold dynamics, focusing on the topological distribution of fixed points, eigenvalue spectra, uniform flow norms, and manifold reliability (see Appendix D for detailed analytical methods). Notably, as detailed in Appendix C, both GRU and LSTM failed to converge on the memory-guided saccade task, and the ReLU-LSTM variant was untrainable due to gradient instability. This failure reflects their structural bias toward discrete attractors, which may cause catastrophic state drift during prolonged zero-input delays (Jordan et al., 2021). Consequently, our analysis for the memory-guided saccade task focuses exclusively on the RDNN. Furthermore, to assess the generalizability of our findings, we extend these dynamical and performance evaluations to alternative DN-equipped networks (ORGaNICs, see Appendix K), Neural ODEs (see Appendix L), and a higher-dimensional coupled task (Double Angular Velocity Integration, see Appendix O). These evaluations indicate that divisive normalization consistently facilitates the formation of stable slow manifolds across different model implementations and higher-dimensional toroidal tasks, whereas continuous-time baselines lacking this mechanism fail to maintain normal hyperbolicity. Figure 2: Dynamical topology and spectral properties of the learned slow manifolds. A Angular Velocity Integration. Left: Representative fixed-point distributions on the output ring. Middle: Statistical counts of topological states across random seeds. Right: Eigenvalue modulus density. RDNN shows a structured bimodal distribution indicating normal hyperbolicity, whereas baselines exhibit diffuse spectra without a clear gap. B Memory-guided Saccade (RDNN only). In this autonomous maintenance task, the RDNN manifold is overwhelmingly dominated by marginal fixed points (left, middle), approximating a theoretical continuous attractor. The corresponding eigenvalue spectrum (right) displays a strict spectral gap, indicating near-perfect time-scale separation. 4.1 Input-Driven Integration: Rotational Dynamics vs. Discretized States The angular velocity integration task requires the network to continuously rotate its internal state in response to external inputs. Here, the topological strategies diverge significantly between architectures (Figure 2A, middle). GRU and LSTM networks predominantly form alternating discrete stable fixed points and saddle points. This suggests that gated RNNs effectively discretize the continuous variable, forcing the state to jump between localized basins of attraction (Ghazizadeh and Ching, 2021). In contrast, the RDNN exhibits a more fluid topological structure. Its state space is frequently dominated by saddle points outnumbering stable fixed points, supplemented by marginal fixed points and limit cycles, indicating a flatter energy landscape operating near a bifurcation, e.g., saddle-node on an invariant circle (Strogatz, 2018; Seung, 1996; Maheswaranathan et al., 2019). Rather than acting as strict barriers, these saddles and marginal points function as slow channels (Strogatz, 2018; Ermentrout and Kopell, 1986), or stable heteroclinic channels (Rabinovich et al., 2008). This fluid topology allows external inputs to smoothly drive the state along the manifold without being trapped in deep localized attractors (Khona and Fiete, 2022). This distinction is further illuminated by the complex eigenvalue spectrum of the Jacobians (Figure A6). The RDNN displays a highly structured, tripartite spectrum: the majority of eigenvalues cluster strictly along the real axis, while the remaining modes form off-axis conjugate pairs spanning the interior, alongside a distinct set of boundary conjugate pairs near the unit circle (|z|≈1,Im(z)≠0|z|≈ 1,Im(z)≠ 0). These boundary conjugate pairs provide the marginally stable rotational dynamics required for continuous integration (Mastrogiuseppe and Ostojic, 2018; Sussillo and Barak, 2013), while the interior conjugate pairs govern damped transient oscillations that stabilize state transitions along the manifold (Hennequin et al., 2014). Conversely, the spectra of GRU and LSTM resemble a diffuse cloud, lacking both a clear spectral gap and structured rotational modes (Rajan and Abbott, 2006; Farrell et al., 2022). Consequently, the RDNN maintains a significantly lower uniform norm than the baselines (Figure A7, left), reflecting a smoother manifold that is less susceptible to input-induced drift. 4.2 Autonomous Maintenance: Emergence of Approximate Continuous Attractors In the autonomous memory-guided saccade task, the RDNN exhibits dynamics that closely approximate a theoretical continuous attractor. Topological classification of the identified slow manifolds (Figure 2B, middle) reveals that the RDNN’s state space is overwhelmingly dominated by marginal fixed points. This neutral stability along the manifold is a defining characteristic of continuous attractors (Seung, 1996; Khona and Fiete, 2022), allowing the network to maintain arbitrary continuous values without drifting toward discrete attractors (Compte, 2000; Koulakov et al., 2002). This near-critical dynamic regime is quantitatively supported by the uniform norm of the vector field restricted to the manifold (Figure A7, right). The RDNN achieves an exceptionally low uniform norm (on the order of 10−310^-3 to 10−410^-4), indicating an extremely flat energy landscape. Instance-level analysis (Figure A9) further confirms that this continuous attractor-like topology is consistently learned across different random seeds and hidden sizes, demonstrating that divisive normalization provides a strong inductive bias for high-fidelity autonomous memory. 4.3 Spectral Properties, Activations, and Manifold Reliability The eigenvalue modulus density (Figure 2A and 2B, right) highlights the normal hyperbolicity of the RDNN (Fenichel, 1979; Chaudhuri et al., 2019), though its strictness depends inherently on the task regime. In the autonomous memory-guided saccade task, the RDNN exhibits a sharp bimodal distribution with a strict spectral gap between peaks near 0 (fast dissipation) and 1 (slow maintenance). As derived in Appendix E.3, this separation is actively sculpted by BPTT and noise-induced regularization, catalyzed by the dynamic gain control of divisive scaling. However, in the input-driven angular velocity integration task, this gap is partially filled with intermediate eigenvalues. This non-zero density is a functional necessity: it reflects the intermediate time-scale modes required to continuously process external velocity inputs and couple them to the slow manifold (Bondanelli and Ostojic, 2020; Hennequin et al., 2014). Interestingly, the choice of activation function subtly modulates this spectrum during active integration. Specifically in the angular velocity integration task, the ReLU activation in the RDNN produces a significantly higher peak at 0 compared to Tanh. This suggests that ReLU’s hard-thresholding, when stabilized by divisive normalization, induces a sparser representation during dynamic updating (Glorot et al., 2011). Strictly silent neurons provide zero gradients, thereby instantly annihilating orthogonal noise and creating a harder manifold during continuous integration (Pennington et al., 2017; Engelken et al., 2023; Laurent and von Brecht, 2017; Ahmad and Hawkins, 2016). Finally, this structured dynamical bias translates directly to training reliability. As shown in Figure A5, the RDNN successfully forms reliable, topologically correct ring manifolds in approximately 90% of the trained instances. In contrast, GRU and LSTM models yield reliable manifolds in less than 50% of cases, further underscoring the necessity of divisive normalization for robust continuous working memory. 5 Low-Rank Analysis The algebraic rank of the recurrent connectivity matrix fundamentally bounds and shapes the dimensionality of the learned state-space trajectories in dynamical systems theory (Mastrogiuseppe and Ostojic, 2018). To understand whether the low-dimensional representations discovered before are a fundamental property of divisive normalization or merely a coincidental feature of training, we conduct a comprehensive low-rank analysis. We first examine the unconstrained effective rank of different architectures. Motivated by the naturally emergent low-rank structure of the RDNN, we then investigate the effects of forcing an explicit low-rank bottleneck by training an H=128H=128 RDNN with explicitly factorized recurrent weights, =12J=J_1J_2 and =12W=W_1W_2, across varying bottleneck ranks r∈8,16,24,32r∈\8,16,24,32\. Figure 3: Effective rank scaling and explicit low-rank training dynamics. A The 99% energy effective rank as a function of physical hidden size (H). Across both tasks, the RDNN maintains an exceptionally compact, sub-linear scaling, whereas standard gated baselines (GRU, LSTM) scale linearly and exploit nearly the entire physical state space (GRU and LSTM are absent in the memory-guided saccade task due to non-convergence). B Training convergence curves (Loss vs. Epoch) of the explicitly factorized low-rank RDNN (H=128H=128) across bottleneck ranks r∈8,16,24,32r∈\8,16,24,32\. The memory-guided saccade task (right) illustrates the higher-rank convergence anomaly, where wider explicit bottlenecks (e.g., r=32r=32) suffer from wider optimization variance and slower descent compared to tighter bottlenecks (e.g., r=8r=8). Error bands and error bars represent the standard deviation across random seeds. 5.1 Effective Rank Scaling and the Self-Compression Phenomenon We first compare the 99% energy effective rank of the unconstrained models across different hidden sizes (Figure 3A and Table 1), and a complementary component-wise analysis is provided in Appendix J to verify the robustness of this joint metric. While the effective ranks of the GRU and LSTM baselines scale linearly with the physical hidden size H (essentially exploiting all available degrees of freedom), the RDNN exhibits a highly sub-linear, flat scaling. For instance, at H=256H=256 in the integration task, the RDNN maintains an effective rank of approximately 35.4∼38.035.4 38.0, whereas the baselines require a rank of over 240240. Mathematically, the inverse gradient scaling of the divisive term acts as an activity-dependent gradient attenuator during backpropagation, dampening the magnitude of parameter updates in highly active directions (see Appendix E.2 and E.4). While this local scaling does not formally minimize a global rank metric, it restricts the effective directions of parameter updates, which empirically correlates with the observed low-rank representations. Table 1: The 99% energy effective rank (mean ± std over different random seeds) across different network architectures, activation functions, and hidden sizes (H). Model Activation Angular Velocity Integration Memory-guided Saccade H=64H=64 H=128H=128 H=256H=256 H=64H=64 H=128H=128 H=256H=256 RDNN (Ours) relu 16.8±2.416.8± 2.4 24.0±2.524.0± 2.5 38.0±6.738.0± 6.7 19.2±2.519.2± 2.5 27.0±7.227.0± 7.2 31.0±5.031.0± 5.0 tanh 17.2±1.917.2± 1.9 22.6±2.622.6± 2.6 35.4±5.335.4± 5.3 16.0±2.016.0± 2.0 22.6±3.522.6± 3.5 28.3±5.628.3± 5.6 GRU relu 61.0±0.061.0± 0.0 121.0±0.0121.0± 0.0 240.8±0.4240.8± 0.4 — — — tanh 60.0±0.060.0± 0.0 120.6±0.5120.6± 0.5 240.2±0.4240.2± 0.4 — — — LSTM tanh 61.0±0.061.0± 0.0 122.0±0.0122.0± 0.0 244.0±0.0244.0± 0.0 — — — When we enforce an explicit bottleneck of rank r on the H=128H=128 RDNN, a self-compression phenomenon emerges during training. As illustrated in Figure A11, regardless of how wide we set the explicit bottleneck r∈8,16,24,32r∈\8,16,24,32\, the trained networks consistently compress their effective ranks back to a much lower, intrinsic dimensionality: approximately 3.5∼4.03.5 4.0 for the integration task, and ∼2.0 2.0 for the saccade task. This extreme self-compression indicates that the intrinsic dimensionality of the slow manifold for these continuous tasks is fundamentally low. Intriguingly, the effective rank of these explicit low-rank networks is significantly lower than that of the unconstrained RDNN, and this hyper-compression is observed to be an outcome of the coupled BPTT gradient updates in the factorized parameter space, which contribute to an accelerated singular value decay. This process appears to amplify the local gradient-scaling effects introduced by divisive normalization (see Appendix E.5.2). 5.2 The Higher-Rank Convergence Anomaly A key feature of low-rank parameterization is its impact on the optimization landscape. Counter-intuitively, our training curves (Figure 3B) reveal a higher-rank convergence anomaly: in many instances, networks with higher explicit ranks (e.g., r=32r=32) are significantly more difficult to optimize and converge more slowly than those with lower explicit ranks (e.g., r=8r=8). This is particularly evident in the memory-guided saccade task (Figure 3B, right), where the r=32r=32 configurations exhibit a massive variance in training loss across seeds and suffer from prolonged convergence plateaus, whereas the r=8r=8 networks converge more rapidly and stably. As derived in Appendix E.5.3, this anomaly can be understood through the lens of non-convex matrix factorization. A wider bottleneck (r=32r=32) increases the dimensionality of the search space, introducing a vast number of redundant scaling symmetries, flat saddle points, and degenerate directions in the product space 12J_1J_2 (Li et al., 2019; Valavi et al., 2020). Conversely, a tight bottleneck (r=8r=8) restricts the gradient search strictly within the low-dimensional subspace relevant to the continuous attractor, effectively smoothing the optimization landscape (Arora et al., 2019). In contrast, the unconstrained RDNN is parameterized directly without product-form weights. This flat parameterization avoids these non-convex scaling symmetries; as a result, increasing the hidden size H of the unconstrained RDNN leads to standard over-parameterization, which mathematically smoothes the loss landscape and facilitates optimization (Du et al., 2019; Cooper, 2018). Consequently, the unconstrained RDNN scales gracefully with H, while the factorized network suffers from the bottleneck anomaly (see Appendix E.5.4). 5.3 Spectral Integrity and the Degradation of Normal Hyperbolicity Although the explicit low-rank networks achieve highly competitive test performance (Figure A10), forcing an explicit factorization has noticeable consequences on their spectral properties. As shown in the eigenvalue modulus density (Figure A12) and complex plane scattering (Figure A13), the bimodal spectrum, comprising a sharp peak near 0 (fast modes) and a peak near 1 (slow modes), is largely preserved, demonstrating the robustness of the divisor-gated architecture. However, in the memory-guided saccade task (Figure A12, right), we observe a subtle degradation of the spectral gap in some explicit low-rank networks. Unlike the unconstrained RDNN, which exhibits a strict zero density in the middle of the spectrum, some explicit low-rank variants show a non-zero leakage of eigenvalues in the intermediate region between the two peaks. Dynamically, this spectral leakage indicates that the hard factorization limits the network’s ability to perfectly coordinate the feedback loops between the excitatory state R and the inhibitory state G (Murphy and Miller, 2009; Mastrogiuseppe and Ostojic, 2018). This restriction introduces weak, intermediate-time-scale modes that do not decay rapidly enough, slightly compromising the normal hyperbolicity of the slow manifold and, to some extent, explaining the minor performance gap compared to the unconstrained RDNN (Chaudhuri et al., 2019; Ságodi et al., 2024) (Figure A10). 6 Ablation Study: Divisive vs. Subtractive Inhibition In biological neural circuits, inhibition primarily manifests in two distinct forms: subtractive inhibition, which shifts the membrane potential linearly, and divisive inhibition, which multiplicatively scales the neural response (Carandini and Heeger, 2012; Holt and Koch, 1997). To isolate the specific computational benefits of divisive normalization, we conduct an ablation study by replacing the divisive gate in the RDNN with a subtractive inhibitory term (referred to as the Subtractive Network). Specifically, the continuous-time excitatory dynamics are modified to: τd=(−+f()+(t)−)dt+σd,←max(,), _RdR= (-R+Jf(R)+I(t)-G )dt+ _RdW_R, ← (0,R), (2) where the state is continuously rectified to enforce non-negative representations. While both networks achieve competitive task performance (Figure A14), we investigate whether the multiplicative nature of divisive gain control is strictly necessary for the underlying slow manifold dynamics. Figure 4: Comparison of dynamical topology between divisive and subtractive networks. Statistical counts of topological types across random seeds and configurations. A Angular Velocity Integration. Under continuous velocity inputs, the Subtractive Network’s manifold fragments into a portfolio heavily populated by stable fixed points (blue). In contrast, the RDNN exhibits a more fluid topology governed by saddles (acting as slow channels). B Memory-guided Saccade. In this autonomous task ((t)=I(t)=0), both networks successfully sustain marginal fixed points (red), demonstrating that subtractive inhibition is sufficient for static memory maintenance when dynamic gain control is not required. Figure 5: The 99% energy effective rank of the merged recurrent weight matrices is plotted against hidden sizes (H∈64,128,256H∈\64,128,256\). While both the RDNN (blue) and the Subtractive Network (orange) maintain a sub-linear scaling with H, the RDNN’s effective rank is consistently and significantly lower across both tasks. This quantitative gap indicates that multiplicative divisive normalization more effectively confines the recurrent dynamics to a lower-dimensional subspace during training compared to additive subtractive inhibition, potentially due to its activity-dependent gradient-attenuation properties. 6.1 Attractor Shattering in Input-Driven Integration The main divergence between the two inhibitory mechanisms emerges in the input-driven angular velocity integration task. Notably, the Subtractive Network’s slow manifold shatters into a discretized state space heavily populated by stable fixed points (Figure 4A), which approximate or outnumber saddle points across most configurations. This topological fragmentation indicates that the continuous manifold has broken into deep, localized point-attractor basins that easily trap the network state during integration. In contrast, the RDNN’s topology is predominantly governed by saddle points. In driven dynamical systems, rather than acting as strict barriers, these saddles function as slow channels or ghost attractors (Strogatz, 2018), facilitating smooth, continuous state transitions along the manifold without trapping the trajectory. While both networks exhibit similar eigenvalue modulus density, with a prominent eigenvalue modulus peak near |λ|≈1|λ|≈ 1, to satisfy the BPTT memory constraint during integration (Figure A15, left), their underlying dimensionalities differ significantly. As quantitatively corroborated by the elevated effective rank of the Subtractive Network (Figure 5 and Table A1), the ablated network recruited additional physical dimensions to support its dynamics. As derived in Appendix E.5.1, this dimensional compensation stems from a gradient-level decoupling: unlike divisive normalization, additive inhibition fails to multiplicatively scale the BPTT gradients. While the hard-thresholding of ReLU still suppresses strictly silent neurons, the subtractive network lacks the multiplicative gain control of divisive normalization that, in the RDNN, contributes to dynamically attenuating gradients along active normal dimensions (|λnormal|→0| _normal|→ 0). Without this dynamic gradient-scaling capacity to confine its dynamics, the optimizer in the subtractive network tends to construct redundant slow pathways across a fragmented, higher-dimensional landscape, hindering the formation of a strict one-dimensional ring for continuous integration under time-varying inputs. 6.2 Autonomous Maintenance and Effective Rank Compression Conversely, in the autonomous memory-guided saccade task, where external input is absent during the delay period, the Subtractive Network successfully forms marginal fixed points and maintains a strict spectral gap, closely mirroring the RDNN (Figure 4B and A15, right). Indeed, both networks exhibit nearly identical fixed-point portfolios (dominated by marginal fixed points) and eigenvalue modulus density spectra. This similarity in the autonomous regime ((t)=0I(t)=0) is mathematically expected: under zero-input conditions, the primary computational requirement is to balance recurrent excitation with inhibitory feedback to achieve marginal stability (λmax≈1 _ ≈ 1). Since the input is absent during the delay, the unique scale-invariance and gain-control properties of divisive normalization are not put to the test, allowing both linear subtractive and non-linear divisive inhibition to sustain a functionally equivalent ring attractor (Ben-Yishai et al., 1995; Zhang, 1996; Ayaz and Chance, 2009). As derived in Appendix E.5.1, while time-varying inputs cause discontinuous state-switching in the Heaviside gating that shatters the ring (Hahnloser et al., 2000), this gating remains static under zero-input conditions, preserving the manifold. However, quantitative analysis of the 99% energy effective rank (Figure 5 and Table A1) reveals a fundamental difference between the networks. While the Subtractive Network successfully avoids the linear rank explosion seen in standard gated RNNs (exhibiting sub-linear scaling with hidden size H), its effective rank remains consistently higher than that of the RDNN across all configurations. This suggests that the multiplicative scaling of divisive normalization, by its activity-dependent gradient attenuation, contributes to effectively guiding the recurrent dynamics into a tighter, lower-dimensional subspace than additive subtraction. 6.3 Functional Implications in Neuroscience These findings provide a compelling computational perspective on the coexistence of divisive and subtractive inhibition in biological systems, delineating their distinct functional niches. Subtractive inhibition is highly suited for static thresholding, baseline noise filtering, and sparse coding (Carandini and Heeger, 2012; Silver, 2010). By shifting the neural activation curve to the right, subtractive mechanisms act as an effective noise gate that filters out weak spontaneous background activity and sharpens sensory tuning curves. This makes subtractive inhibition computationally efficient for static memory maintenance in autonomous regimes (as demonstrated in our memory-guided saccade task) and binary state switching. However, when a circuit must continuously integrate time-varying external inputs, additive subtractive mechanisms fail to scale dynamically with the input magnitude, leading to attractor shattering. Divisive normalization, by providing contrast invariance and dynamic gain control (Heeger, 1992), preserves the structural integrity of the low-dimensional manifold under external drive. This suggests that divisive normalization is not merely a biological artifact, but a computational necessity for high-fidelity, input-driven continuous working memory. Nonetheless, our long-horizon analysis (see Appendix N) reveals that the flat manifolds of divisive normalization are susceptible to asymptotic random-walk drift over extended durations, whereas subtractive discretization bounds long-term error via an error-correcting “lock-in” effect (Koulakov et al., 2002; Brody et al., 2003). This suggests a functional synergy in cortical circuits, where divisive normalization is utilized for contrast-invariant dynamic tracking, while subtractive mechanisms provide noise-resistant stability over long timescales. 7 Discussion and Conclusion Divisive Normalization as a Canonical Computation Our findings provide a novel dynamical systems perspective on why divisive normalization (DN) is widely regarded as a canonical neural computation (Carandini and Heeger, 2012). While traditionally understood as a mechanism for sensory gain control and contrast invariance, our results demonstrate its role in shaping the topological landscape of recurrent memory circuits. By its dynamic gain control and the resulting activity-dependent gradient attenuation that appears to contribute to a lower effective rank, DN allows neural populations to robustly maintain continuous attractors under time-varying inputs, a computational feat where subtractive inhibition structurally fails. This suggests that DN is computationally essential not just for sensory encoding, but as a fundamental building block for the stable maintenance and manipulation of continuous cognitive variables. Local Gradient Scaling and Optimization Dynamics Our low-rank analysis sheds light on how different dimensionality reduction strategies influence the optimization landscape of recurrent networks. While the unconstrained RDNN, through its activity-dependent divisive bottleneck, shows an empirical tendency to maintain a low effective rank, explicitly factorizing the recurrent weights introduces severe optimization challenges. Specifically, we observed a higher-rank convergence anomaly where wider explicit bottlenecks paradoxically hinder training. This aligns with non-convex matrix factorization theory, which shows that over-parameterized product spaces introduce redundant scaling symmetries and degenerate saddle points (Li et al., 2019). Conversely, the unconstrained network avoids these non-convexities and benefits from standard over-parameterization, which mathematically smoothes the loss landscape (Du et al., 2019). These findings suggest that biologically inspired multiplicative gain control offers a highly optimization-friendly route to discovering low-dimensional cognitive representations without relying on hard architectural bottlenecks. Representations of Continuous Manifolds From a machine learning perspective, our study suggests how neural architectures can be designed to better support continuous state spaces. While standard gated RNNs (e.g., GRUs, LSTMs) tend to discretize continuous manifolds into fragmented point attractors, the RDNN indicates that multiplicative divisive feedback can serve as a beneficial inductive bias. This mechanism allows the network to maintain low-rank, continuous representations without the optimization difficulties typical of explicit low-rank factorization. Although our current validation is restricted to canonical cognitive tasks, this architectural feature suggests a potential foundation for continuous-time models requiring both stable memory maintenance and continuous integration. Further investigation is necessary to determine how these mechanisms scale to more complex, high-dimensional engineering applications such as spatial navigation or robotics. Limitations and Future Directions Despite these theoretical insights, our study has several limitations that pave the way for future research. First, while Backpropagation Through Time (BPTT) provides a mathematically rigorous framework to understand the emergence of low-rank dynamics and spectral compression, it lacks biological plausibility. Future work should investigate whether more physiological synaptic plasticity rules, such as Hebbian learning or predictive coding approximations (Lillicrap et al., 2020), can similarly induce these topological properties in DN-equipped networks. Second, the current model focuses exclusively on continuous working memory. Given that biological systems seamlessly integrate both discrete and continuous information (Constantinidis et al., 2018), exploring how divisive normalization might support mixed discrete-continuous memory representations is a compelling next step. Additionally, because we model divisive normalization in isolation, our framework does not yet address the potential synergy of multiple co-existing inhibitory and homeostatic mechanisms, such as subtractive thresholding and synaptic scaling (Silver, 2010; Renart et al., 2003; Carandini and Heeger, 2012). Investigating the optimization and dynamical properties of such multi-component inhibitory networks remains an important avenue for future study. Finally, as a theoretical and computational model, our predictions regarding the eigenvalue spectra and effective rank of recurrent circuits require further empirical validation using large-scale neural population recordings from behaving animals (Chaudhuri et al., 2019). Conclusion In conclusion, we introduced the RDNN to bridge the gap between biological gain control and artificial recurrent memory. Through dynamical and optimization analyses, we established that divisive normalization naturally induces low-rank, normally hyperbolic slow manifolds. By overcoming the attractor shattering and fine-tuning problems inherent in standard architectures, our work provides a unified theoretical framework for understanding and modeling robust continuous working memory in both biological and artificial systems. References S. Ahmad and J. Hawkins (2016) How do neurons operate on sparse distributed representations? A mathematical theory of sparsity, neurons and active dendrites. arXiv. External Links: Document Cited by: §4.3. S. Arora, N. Cohen, W. Hu, and Y. Luo (2019) Implicit regularization in deep matrix factorization. In Proceedings of the 33rd International Conference on Neural Information Processing Systems, Cited by: 3rd item, §E.5.1, §E.5.2, §5.2. A. Ayaz and F. S. Chance (2009) Gain Modulation of Neuronal Responses by Subtractive and Divisive Mechanisms of Inhibition. Journal of Neurophysiology 101 (2), p. 958–968. External Links: ISSN 0022-3077, 1522-1598, Document Cited by: §6.2. R. Ben-Yishai, R. L. Bar-Or, and H. Sompolinsky (1995) Theory of orientation tuning in visual cortex.. Proceedings of the National Academy of Sciences 92 (9), p. 3844–3848. External Links: ISSN 0027-8424, 1091-6490, Document Cited by: §6.2. G. Bondanelli and S. Ostojic (2020) Coding with transient trajectories in recurrent neural networks. PLOS Computational Biology 16 (2), p. e1007655. External Links: ISSN 1553-7358, Document Cited by: §4.3. C. D. Brody, R. Romo, and A. Kepecs (2003) Basic mechanisms for graded persistent activity: discrete attractors, continuous attractors, and dynamic representations. Current Opinion in Neurobiology 13 (2), p. 204–211. External Links: ISSN 09594388, Document Cited by: Appendix N, §6.3. M. F. Burg, S. A. Cadena, G. H. Denfield, E. Y. Walker, A. S. Tolias, M. Bethge, and A. S. Ecker (2021) Learning divisive normalization in primary visual cortex. PLOS Computational Biology 17 (6), p. e1009028. External Links: ISSN 1553-7358, Document Cited by: §2. M. Carandini and D. J. Heeger (2012) Normalization as a canonical neural computation. Nature Reviews Neuroscience 13 (1), p. 51–62. External Links: ISSN 1471-003X, 1471-0048, Document Cited by: §1, §2, §6.3, §6, §7, §7. R. Chaudhuri, B. Gerçek, B. Pandey, A. Peyrache, and I. Fiete (2019) The intrinsic attractor manifold and population dynamics of a canonical cognitive circuit across waking and sleep. Nature Neuroscience 22 (9), p. 1512–1520. External Links: ISSN 1097-6256, 1546-1726, Document Cited by: §4.3, §5.3, §7. R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud (2018) Neural ordinary differential equations. In Advances in Neural Information Processing Systems, S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (Eds.), Vol. 31. Cited by: Appendix L. A. Compte (2000) Synaptic Mechanisms and Network Dynamics Underlying Spatial Working Memory in a Cortical Network Model. Cerebral Cortex 10 (9), p. 910–923. External Links: ISSN 14602199, Document Cited by: §4.2. C. Constantinidis, S. Funahashi, D. Lee, J. D. Murray, X. Qi, M. Wang, and A. F.T. Arnsten (2018) Persistent Spiking Activity Underlies Working Memory. The Journal of Neuroscience 38 (32), p. 7020–7028. External Links: ISSN 0270-6474, 1529-2401, Document Cited by: §7. Y. Cooper (2018) The loss landscape of overparameterized neural networks. arXiv. External Links: Document Cited by: item 2, §5.2. L. Dinh, R. Pascanu, S. Bengio, and Y. Bengio (2017) Sharp minima can generalize for deep nets. In Proceedings of the 34th International Conference on Machine Learning, D. Precup and Y. W. Teh (Eds.), Proceedings of Machine Learning Research, Vol. 70, p. 1019–1028. Cited by: item 1. S. Du, J. Lee, H. Li, L. Wang, and X. Zhai (2019) Gradient descent finds global minima of deep neural networks. In Proceedings of the 36th International Conference on Machine Learning, K. Chaudhuri and R. Salakhutdinov (Eds.), Proceedings of Machine Learning Research, Vol. 97, p. 1675–1685. Cited by: item 2, §5.2, §7. R. Engelken, F. Wolf, and L. F. Abbott (2023) Lyapunov spectra of chaotic recurrent neural networks. Physical Review Research 5 (4), p. 043044. External Links: ISSN 2643-1564, Document Cited by: §4.3. G. B. Ermentrout and N. Kopell (1986) Parabolic Bursting in an Excitable System Coupled with a Slow Oscillation. SIAM Journal on Applied Mathematics 46 (2), p. 233–253. External Links: ISSN 0036-1399, 1095-712X, Document Cited by: §4.1. A. A. Faisal, L. P. J. Selen, and D. M. Wolpert (2008) Noise in the nervous system. Nature Reviews Neuroscience 9 (4), p. 292–303. External Links: ISSN 1471-003X, 1471-0048, Document Cited by: §3. M. Farrell, S. Recanatesi, T. Moore, G. Lajoie, and E. Shea-Brown (2022) Gradient-based learning drives robust representations in recurrent neural networks by balancing compression and expansion. Nature Machine Intelligence 4 (6), p. 564–573. External Links: ISSN 2522-5839, Document Cited by: §4.1. N. Fenichel (1979) Geometric singular perturbation theory for ordinary differential equations. Journal of Differential Equations 31 (1), p. 53–98. External Links: ISSN 00220396, Document Cited by: §4.3. D. Fitzpatrick, J. S. Lund, D. E. Schmechel, and A. C. Towles (1987) Distribution of GABAergic neurons and axon terminals in the macaque striate cortex. Journal of Comparative Neurology 264 (1), p. 73–91. External Links: ISSN 0021-9967, 1096-9861, Document Cited by: §3. R. Ge, C. Jin, and Y. Zheng (2017) No spurious local minima in nonconvex low rank problems: a unified geometric analysis. In Proceedings of the 34th International Conference on Machine Learning, D. Precup and Y. W. Teh (Eds.), Proceedings of Machine Learning Research, Vol. 70, p. 1233–1242. Cited by: 2nd item. E. Ghazizadeh and S. Ching (2021) Slow manifolds within network dynamics encode working memory efficiently and robustly. PLOS Computational Biology 17 (9), p. e1009366. External Links: ISSN 1553-7358, Document Cited by: §2, §4.1. X. Glorot, A. Bordes, and Y. Bengio (2011) Deep sparse rectifier neural networks. In Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, G. Gordon, D. Dunson, and M. Dudík (Eds.), Proceedings of Machine Learning Research, Vol. 15, Fort Lauderdale, FL, USA, p. 315–323. Cited by: §E.5.1, §4.3. K. Greff, R. K. Srivastava, J. Koutnik, B. R. Steunebrink, and J. Schmidhuber (2017) LSTM: A Search Space Odyssey. IEEE Transactions on Neural Networks and Learning Systems 28 (10), p. 2222–2232. External Links: ISSN 2162-237X, 2162-2388, Document Cited by: 1st item. S. Gunasekar, B. Woodworth, S. Bhojanapalli, B. Neyshabur, and N. Srebro (2017) Implicit regularization in matrix factorization. In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, Red Hook, NY, USA, p. 6152–6160. External Links: ISBN 978-1-5108-6096-4 Cited by: §E.4, §E.5.1, §E.5.2. R. H. R. Hahnloser, R. Sarpeshkar, M. A. Mahowald, R. J. Douglas, and H. S. Seung (2000) Digital selection and analogue amplification coexist in a cortex-inspired silicon circuit. Nature 405 (6789), p. 947–951. External Links: ISSN 0028-0836, 1476-4687, Document Cited by: §6.2. D. J. Heeger and W. E. Mackey (2019) Oscillatory recurrent gated neural integrator circuits (ORGaNICs), a unifying theoretical framework for neural dynamics. Proceedings of the National Academy of Sciences 116 (45), p. 22783–22794. External Links: ISSN 0027-8424, 1091-6490, Document Cited by: §1, §2. D. J. Heeger (1992) Normalization of cell responses in cat striate cortex. Visual Neuroscience 9 (2), p. 181–197. External Links: ISSN 0952-5238, 1469-8714, Document Cited by: §2, §6.3. G. Hennequin, T. P. Vogels, and W. Gerstner (2014) Optimal Control of Transient Dynamics in Balanced Networks Supports Generation of Complex Movements. Neuron 82 (6), p. 1394–1406. External Links: ISSN 08966273, Document Cited by: §4.1, §4.3. G. R. Holt and C. Koch (1997) Shunting Inhibition Does Not Have a Divisive Effect on Firing Rates. Neural Computation 9 (5), p. 1001–1013. External Links: ISSN 0899-7667, 1530-888X, Document Cited by: §6. S. Ioffe and C. Szegedy (2015) Batch normalization: Accelerating deep network training by reducing internal covariate shift. In Proceedings of the 32nd International Conference on Machine Learning, F. Bach and D. Blei (Eds.), Proceedings of Machine Learning Research, Vol. 37, Lille, France, p. 448–456. Cited by: §2. I. D. Jordan, P. A. Sokół, and I. M. Park (2021) Gated Recurrent Units Viewed Through the Lens of Continuous Time Dynamical Systems. Frontiers in Computational Neuroscience 15, p. 678158. External Links: ISSN 1662-5188, Document Cited by: 2nd item, §1, §2, §4. R. Jozefowicz, W. Zaremba, and I. Sutskever (2015) An empirical exploration of recurrent network architectures. In Proceedings of the 32nd International Conference on Machine Learning, F. Bach and D. Blei (Eds.), Proceedings of Machine Learning Research, Vol. 37, Lille, France, p. 2342–2350. Cited by: 1st item. M. Khona and I. R. Fiete (2022) Attractor and integrator networks in the brain. Nature Reviews Neuroscience 23 (12), p. 744–766. External Links: ISSN 1471-003X, 1471-0048, Document Cited by: §1, §2, §4.1, §4.2. A. A. Koulakov, S. Raghavachari, A. Kepecs, and J. E. Lisman (2002) Model for a robust neural integrator. Nature Neuroscience 5 (8), p. 775–782. External Links: ISSN 1097-6256, 1546-1726, Document Cited by: Appendix N, §1, §2, §4.2, §6.3. T. Laurent and J. von Brecht (2017) A recurrent neural network without chaos. In International Conference on Learning Representations, Cited by: §4.3. Q. V. Le, N. Jaitly, and G. E. Hinton (2015) A Simple Way to Initialize Recurrent Networks of Rectified Linear Units. arXiv. External Links: Document Cited by: Appendix B, 1st item. X. Li, J. Lu, R. Arora, J. Haupt, H. Liu, Z. Wang, and T. Zhao (2019) Symmetry, Saddle Points, and Global Optimization Landscape of Nonconvex Matrix Factorization. IEEE Transactions on Information Theory 65 (6), p. 3489–3514. External Links: ISSN 0018-9448, 1557-9654, Document Cited by: 2nd item, item 1, §5.2, §7. Y. Li, T. Ma, and H. Zhang (2018) Algorithmic regularization in over-parameterized matrix sensing and neural networks with quadratic activations. In Proceedings of the 31st Conference on Learning Theory, S. Bubeck, V. Perchet, and P. Rigollet (Eds.), Proceedings of Machine Learning Research, Vol. 75, p. 2–47. Cited by: 1st item. T. P. Lillicrap, A. Santoro, L. Marris, C. J. Akerman, and G. Hinton (2020) Backpropagation and the brain. Nature Reviews Neuroscience 21 (6), p. 335–346. External Links: ISSN 1471-003X, 1471-0048, Document Cited by: §7. N. Maheswaranathan, A. Williams, M. Golub, S. Ganguli, and D. Sussillo (2019) Universality and individuality in neural dynamics across large populations of recurrent networks. In Advances in Neural Information Processing Systems, H. Wallach, H. Larochelle, A. Beygelzimer, F. dAlché-Buc, E. Fox, and R. Garnett (Eds.), Vol. 32. Cited by: §4.1. F. Mastrogiuseppe, J. Carmona, and C. K. Machens (2025) Stochastic activity in low-rank recurrent neural networks. PLOS Computational Biology 21 (8), p. e1013371. External Links: ISSN 1553-7358, Document Cited by: §2. F. Mastrogiuseppe and S. Ostojic (2018) Linking Connectivity, Dynamics, and Computations in Low-Rank Recurrent Neural Networks. Neuron 99 (3), p. 609–623.e29. External Links: ISSN 08966273, Document Cited by: §2, §4.1, §5.3, §5. B. K. Murphy and K. D. Miller (2009) Balanced Amplification: A New Mechanism of Selective Amplification of Neural Activity Patterns. Neuron 61 (4), p. 635–648. External Links: ISSN 08966273, Document Cited by: §5.3. J. Pennington, S. Schoenholz, and S. Ganguli (2017) Resurrecting the sigmoid in deep learning through dynamical isometry: theory and practice. In Advances in Neural Information Processing Systems, I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Eds.), Vol. 30. Cited by: §4.3. M. I. Rabinovich, R. Huerta, P. Varona, and V. S. Afraimovich (2008) Transient Cognitive Dynamics, Metastability, and Decision Making. PLoS Computational Biology 4 (5), p. e1000072. External Links: ISSN 1553-7358, Document Cited by: §4.1. K. Rajan and L. F. Abbott (2006) Eigenvalue Spectra of Random Matrices for Neural Networks. Physical Review Letters 97 (18), p. 188104. External Links: ISSN 0031-9007, 1079-7114, Document Cited by: §4.1. S. Rawat, D. J. Heeger, and S. Martiniani (2024) Unconditional stability of a recurrent neural circuit implementing divisive normalization. In Advances in Neural Information Processing Systems, A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Paquet, J. Tomczak, and C. Zhang (Eds.), Vol. 37, p. 14712–14750. External Links: Document Cited by: Appendix K, Appendix K, §1, §2. A. Renart, P. Song, and X. Wang (2003) Robust Spatial Working Memory through Homeostatic Synaptic Scaling in Heterogeneous Cortical Networks. Neuron 38 (3), p. 473–485. External Links: ISSN 08966273, Document Cited by: §7. Á. Ságodi, G. Martín-Sánchez, P. Sokół, and I. M. Park (2024) Back to the continuous attractor. In Advances in Neural Information Processing Systems, A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Paquet, J. Tomczak, and C. Zhang (Eds.), Vol. 37, p. 66856–66906. External Links: Document Cited by: Appendix N, §D.2, §2, §5.3. T. Sawada and A. A. Petrov (2017) The divisive normalization model of V1 neurons: a comprehensive comparison of physiological data and model predictions. Journal of Neurophysiology 118 (6), p. 3051–3091. External Links: ISSN 0022-3077, 1522-1598, Document Cited by: §2. H. S. Seung (1996) How the brain keeps the eyes still. Proceedings of the National Academy of Sciences 93 (23), p. 13339–13344. External Links: ISSN 0027-8424, 1091-6490, Document Cited by: §1, §2, §4.1, §4.2. Y. Shao, D. Dahmen, S. Recanatesi, E. Shea-Brown, and S. Ostojic (2025) Impact of Local Connectivity Patterns on Excitatory-Inhibitory Network Dynamics. PRX Life 3 (2), p. 023008. External Links: ISSN 2835-8279, Document Cited by: §2. R. A. Silver (2010) Neuronal arithmetic. Nature Reviews Neuroscience 11 (7), p. 474–489. External Links: ISSN 1471-003X, 1471-0048, Document Cited by: §6.3, §7. S. H. Strogatz (2018) Nonlinear Dynamics and Chaos. 0 edition, CRC Press. External Links: Document, ISBN 978-0-429-96111-3 Cited by: §4.1, §6.1. D. Sussillo and O. Barak (2013) Opening the Black Box: Low-Dimensional Dynamics in High-Dimensional Recurrent Neural Networks. Neural Computation 25 (3), p. 626–649. External Links: ISSN 0899-7667, 1530-888X, Document Cited by: §E.5.1, §2, §4.1. H. Valavi, S. Liu, and P. Ramadge (2020) Revisiting the landscape of matrix factorization. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, S. Chiappa and R. Calandra (Eds.), Proceedings of Machine Learning Research, Vol. 108, p. 1629–1638. Cited by: §5.2. K. Wimmer, D. Q. Nykamp, C. Constantinidis, and A. Compte (2014) Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory. Nature Neuroscience 17 (3), p. 431–439. External Links: ISSN 1097-6256, 1546-1726, Document Cited by: §1, §2. B. Woodworth, S. Gunasekar, J. D. Lee, E. Moroshko, P. Savarese, I. Golan, D. Soudry, and N. Srebro (2020) Kernel and rich regimes in overparametrized models. In Proceedings of Thirty Third Conference on Learning Theory, J. Abernethy and S. Agarwal (Eds.), Proceedings of Machine Learning Research, Vol. 125, p. 3635–3673. Cited by: §E.4. H. Wu, J. Xu, J. Wang, and M. Long (2021) Autoformer: Decomposition transformers with auto-correlation for long-term series forecasting. In Advances in Neural Information Processing Systems, M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. W. Vaughan (Eds.), Vol. 34, p. 22419–22430. Cited by: Appendix M. N. Xiong, L. Ding, and S. S. Du (2024) How over-parameterization slows down gradient descent in matrix sensing: The curses of symmetry and initialization. In The Twelfth International Conference on Learning Representations, Cited by: 2nd item. J. Xu, X. Sun, Z. Zhang, G. Zhao, and J. Lin (2019) Understanding and improving layer normalization. In Advances in Neural Information Processing Systems, H. Wallach, H. Larochelle, A. Beygelzimer, F. dAlché-Buc, E. Fox, and R. Garnett (Eds.), Vol. 32. Cited by: §2. K. Zhang (1996) Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory. The Journal of Neuroscience 16 (6), p. 2112–2126. External Links: ISSN 0270-6474, 1529-2401, Document Cited by: §6.2. B. Zhao, R. Walters, and R. Yu (2026) Symmetry in neural network parameter spaces. Transactions on Machine Learning Research. External Links: ISSN 2835-8856 Cited by: 1st item. Appendix A Tasks We employ three canonical continuous working memory tasks to evaluate the capacity of recurrent networks to form and maintain continuous working memory. These tasks require the network to encode one or more circular variables θ∈[0,2π)θ∈[0,2π) and output their corresponding 2D or 4D Cartesian coordinates. Angular Velocity Integration This task tests the network’s ability to continuously update its internal representation based on a dynamic input. The total length of a trial is 256 time steps. The input is a one-dimensional continuous angular velocity signal v(t)v(t) generated by superimposing three sinusoidal waves. For each wave, the frequency, phase, and amplitude magnitude are sampled from uniform distributions (0.2,2.0)U(0.2,2.0) Hz, (0,2π)U(0,2π), and (0.5,2.0)U(0.5,2.0), respectively, with a random sign. A constant bias is added to ensure the total rotation is uniformly sampled from (1.2π,4.0π)U(1.2π,4.0π). The network is provided with the initial position [sin(θ0),cos(θ0)][ ( _0), ( _0)] (θ0∼(0,2π) _0 (0,2π)) at t=0t=0 and must output the integrated position at each subsequent step (Figure A1). An additional linear mapping is applied in this task to initialize the hidden state onto the initial position along the ring from which the network needed to integrate from. Figure A1: Visualization of the angular velocity integration task. A The 2D output trajectory from RDNN, with the color gradient indicating the instantaneous encoded angle. B Time-course of the true and predicted [sin(θ),cos(θ)][ (θ), (θ)] values for a single trial performed by RDNN. Memory-guided Saccade This task evaluates the network’s ability to maintain a static memory over a prolonged, variable delay period without external input. The total length of a trial is 512 time steps. The input is a 3-dimensional sequence: the first two dimensions provide the target angle [sin(θ),cos(θ)][ (θ), (θ)] (θ∼(0,2π)θ (0,2π)) for a brief stimulus period (15 steps), followed by a variable zero-input delay period sampled from a discrete uniform distribution 50,400U\50,400\, and the fixation signal as the third dimension. A Go Cue is then presented in the third input dimension for 5 steps. Crucially, the network is trained to output [0,0][0,0] (fixation) during the stimulus and delay periods, and must rapidly decode the remembered angle to the output layer only after receiving the Go Cue, with the loss masked out during the 5-step cue transition (Figure A2). Figure A2: Visualization of the memory-guided saccade task. A RDNN’s 2D output trajectories across multiple target angles. The network maintains a zero-output fixation state at the origin before executing rapid radial saccades to the target angles. B Time-course of a single trial. Double Angular Velocity Integration This task is a multi-dimensional, coupled extension of the single-dimensional angular velocity integration task, designed to evaluate how the network generalizes to higher-dimensional continuous manifolds (Figure A3). Instead of tracking a single angle, the network must simultaneously and independently integrate two separate 1D angular velocity signals, v1(t)v_1(t) and v2(t)v_2(t), which are generated independently following the same stochastic sinusoidal process as the single integration task. The input is a 2D velocity vector, and the target output is a 4D vector tracking both integrated angles simultaneously: [sin(θ1,t),cos(θ1,t),sin(θ2,t),cos(θ2,t)][ ( _1,t), ( _1,t), ( _2,t), ( _2,t)]. Mathematically, this continuous integration task maps the recurrent state space onto a two-dimensional torus attractor (T2=S1×S1T^2=S^1× S^1) embedded in a four-dimensional output space. The network is initialized with a 4D initial position at t=0t=0. Figure A3: Visualization of the double angular velocity integration task. Left: RDNN’s time-course of the true (dashed lines) and predicted (solid lines) angles for the two independent circular variables (θ1 _1 and θ2 _2). Right: The 4D predicted state trajectory projected onto a 3D parameteric torus manifold (T2=S1×S1T^2=S^1× S^1, grey mesh). Appendix B Discretization and Implementation of RDNN Dynamics To train the RDNN using backpropagation through time (BPTT) in standard deep learning frameworks, we discretize the continuous-time stochastic differential equations (Eq. 1) using the Euler-Maruyama method. Input Projection and Output Decoding To interface the recurrent core with external task variables, the raw input sequence t∈ℝDinx_t ^D_in at each time step is mapped to the network’s internal input current t∈ℝHI_t ^H via a linear projection layer t=int+inI_t=W_inx_t+b_in, where in∈ℝH×DinW_in ^H× D_in and in∈ℝHb_in ^H are trainable weights and biases. Similarly, the task-specific predictions t∈ℝDouty_t ^D_out are decoded linearly from the principal excitatory state tR_t: t=outt+outy_t=W_outR_t+b_out. Discrete-Time Recurrent Updates Let Δt t be the integration time step. We define the discrete update rates (or leak factors) as R=ΔtτR α_R= t _R and G=ΔtτG α_G= t _G. The discrete-time iterative updates for the hidden states tR_t and tG_t at step t are formulated as: ~t R_t =f(t−1)+t+t−1+R,t, = Jf(R_t-1)+I_t η+G_t-1+ ξ_R,t, (A-1) t _t =(1−R)⊙t−1+R⊙~t, =(1- α_R) _t-1+ α_R R_t, ~t G_t =f(t−1)+G,t, =Wf(R_t-1)+ ξ_G,t, t _t =(1−G)⊙t−1+G⊙~t, =(1- α_G) _t-1+ α_G G_t, where ⊙ denotes element-wise multiplication, and R,t,G,t∼(0,σ2) ξ_R,t, ξ_G,t (0,σ^2I) are the injected Gaussian state noises at each time step. Parameterization and Biological Constraints To ensure the system remains a positive dynamical system (consistent with Dale’s principle and the biological reality of excitatory and inhibitory identities), we enforce strict positivity on specific parameters during training. In our implementation, this is achieved by defining unconstrained raw parameters and applying the softplus function during the forward pass: =Softplus(raw),=Softplus(raw),=Softplus(raw)+ϵ,J=Softplus(J_raw), =Softplus(W_raw), η=Softplus( η_raw)+ε, (A-2) where ϵ=10−5ε=10^-5 is a small constant for numerical stability. Similarly, to ensure that the update rates R α_R and G α_G strictly bound the leaky integration between 0 and 11, we parameterize them using the sigmoid function: R=σ(R,raw),G=σ(G,raw). α_R=σ( α_R,raw), α_G=σ( α_G,raw). (A-3) Initialization Strategy Proper initialization is crucial for training recurrent networks with positive feedback loops (Le et al., 2015). To prevent the product of positive matrices from causing numerical explosion in early training phases, we initialize rawJ_raw and rawW_raw from a normal distribution (−ln(H),1/H)N(- (H),1/ H). This ensures that the initial spectral radius of the recurrent drive is well-behaved. Furthermore, R,raw α_R,raw and G,raw α_G,raw are initialized with a negative mean (e.g., −2.2-2.2), yielding an initial update rate of α≈0.1α≈ 0.1. This explicitly biases the network towards functioning as a slow integrator at the beginning of training, which is essential for learning long-term dependencies in continuous working memory tasks. Appendix C Training Details General Training Setup All models were trained using the Adam optimizer (β1=0.9 _1=0.9, β2=0.999 _2=0.999) with a batch size of 64 for a total of 100 epochs with training data generated online. During training, we applied gradient clipping with a maximum norm of 1.01.0 to prevent gradient explosion. We also injected Gaussian state noise (σ=0.1σ=0.1) into the hidden dynamics at each time step during the training phase. While training low-rank RDNN, we extended the training epochs to 200 to ensure convergence, as the reduced parameterization can lead to slower learning dynamics. For the single and double angular velocity integration task, the models were optimized using the standard Mean Squared Error (MSE) over the entire sequence of length T: ℒint=1T∑t=1T‖t−^t‖22,L_int= 1T _t=1^T\|y_t- y_t\|_2^2, (A-4) where ^t y_t is the network’s prediction at time t. For the one-dimensional task, the target is t=[sin(θt),cos(θt)]y_t=[ ( _t), ( _t)], whereas for the double integration task, the target is the 4D concatenated vector tracking both independent angles: t=[sin(θ1,t),cos(θ1,t),sin(θ2,t),cos(θ2,t)]y_t=[ ( _1,t), ( _1,t), ( _2,t), ( _2,t)]. For the memory-guided saccade task, we utilized a masked MSE loss defined as: ℒsac=∑t=1Tmt‖t−^t‖22D∑t=1Tmt,L_sac= _t=1^Tm_t\|y_t- y_t\|_2^2D _t=1^Tm_t, (A-5) where mt∈0,1m_t∈\0,1\ is a binary mask (mt=0m_t=0 during the Go Cue period, and mt=1m_t=1 otherwise), and D=2D=2 is the output dimension. This formulation ensures that the gradient is normalized correctly by the exact number of active output elements. Hyperparameter Selection and Model Configurations To ensure a fair comparison, we conducted a pilot hyperparameter search for the learning rate. For each combination of network architecture, hidden size (H∈64,128,256H∈\64,128,256\), and activation function, we evaluated learning rates from the set 10−2,5×10−3,10−3\10^-2,5× 10^-3,10^-3\. Each configuration was trained for 10 epochs, and the learning rate 5×10−35× 10^-3 that yielded the lowest training loss was selected for the full 100-epoch training. All reported results and dynamical analyses are based on 5 independent runs with different random seeds for each configuration. Training a single network took about 10∼2010 20 minutes on CPU, and ∼800 800MB of memory. Baseline Initialization For the baseline GRU and LSTM models, we employed specific initialization strategies to improve their ability to capture long-term dependencies. The hidden-to-hidden projection weights were initialized from a normal distribution (0,σ2)N(0,σ^2) with σ=1/Hσ=1/ H. Furthermore, following standard practices for training gated RNNs on memory tasks, the biases corresponding to the update gate in the GRU and the forget gate in the LSTM were explicitly initialized to 1.01.0. This biases the gates towards retaining past information during the early stages of training. Explicit Low-Rank Initialization For the explicitly factorized Low-Rank RDNN, proper initialization is critical to prevent numerical instability caused by the product of positive matrices. We initialize the raw factor matrices (1,2,1,2J_1,J_2,W_1,W_2) from a normal distribution with a negative mean of −1.0-1.0 and standard deviations scaled by 1/r1/ r and 1/H1/ H, respectively. Since these parameters are passed through a Softplus function during the forward pass, the negative mean ensures that the initial effective weights are sufficiently small (Softplus(−1.0)≈0.31Softplus(-1.0)≈ 0.31). This prevents the explosive positive feedback that would otherwise occur when multiplying two strictly positive matrices in the early stages of training. Furthermore, the raw update rates R,raw α_R,raw and G,raw α_G,raw are initialized with a mean of −2.2-2.2. After applying the sigmoid function, this yields an initial update rate of α≈0.1α≈ 0.1, and explicitly biases the network to operate as a slow integrator at the onset of training, facilitating the learning of long-term dependencies in continuous working memory tasks. Exceptions and Convergence Issues During our experiments, we observed two notable training exceptions: • LSTM with ReLU: We excluded the LSTM architecture paired with the ReLU activation function from our analysis. Empirical observations indicated that this specific combination was highly susceptible to gradient explosion and remained largely untrainable across all tested learning rates (Le et al., 2015; Greff et al., 2017; Jozefowicz et al., 2015). • Baseline Failures on the Memory-guided Saccade Task: While the proposed RDNN successfully and stably converged on the memory-guided saccade task, both the GRU and LSTM baselines completely failed to learn the task across all tested learning rates. As illustrated by the training loss curves in Figure A4, the RDNN demonstrates a steady decrease in loss, whereas the GRU and LSTM remain stuck at high error levels. From a dynamical systems perspective, maintaining an arbitrary continuous variable over a prolonged zero-input delay requires synthesizing a perfectly flat energy landscape (i.e., a continuous attractor). Standard gated RNNs are structurally biased toward discrete point attractors; consequently, their hidden states inevitably drift into localized basins during the long delay, resulting in large, uncorrectable readout errors (Jordan et al., 2021). While they might converge with task-specific tricks, the RDNN naturally solves this via its divisive gain control. Consequently, the analysis for the memory-guided saccade task was exclusively performed on the RDNN configurations. Figure A4: Training loss curves across epochs. Error bands indicate standard deviations across different hidden sizes, activation functions, and random seeds. Appendix D Analysis Methods D.1 Evaluation Metric On trained networks, we report the normalized mean squared error (NMSE) of the network prediction y compared to the ground truth target y across 1024 test trials with state noise: NMSE=[‖^−‖22][‖22].NMSE= E [\| y-y\|_2^2 ]E [\|y\|_2^2 ]. (A-6) This metric quantifies the proportion of variance in the target that is not captured by the network’s predictions, with lower values indicating better performance. D.2 Dynamical System Analysis To rigorously characterize the underlying computational mechanisms of the trained networks, we analyzed their autonomous dynamics (i.e., network evolution without external inputs, t=I_t=0) initialized from task-relevant states. The analysis pipeline extracts the slow manifold, identifies topological structures (fixed points and limit cycles), and computes the spectral properties of the system. Identification of the Slow Manifold To identify the slow manifold (approximate continuous attractor) embedded in the high-dimensional state space, we simulated the autonomous dynamics for an extended period (32×32× task length) of 1024 trajectories without state noise starting from states collected at the end of task trials. We computed the instantaneous drift speed in the output space as vt=‖t+1−t‖2v_t=\|y_t+1-y_t\|_2. Points with drift speeds below a progressive relaxation threshold (starting from 10−310^-3 of the maximum trajectory speed) were selected as candidate slow points. To ensure a uniform representation of the ring manifold, we discretized the angular output space [−π,π)[-π,π) into 256 equal bins. Within each angular bin, the state exhibiting the minimum drift speed was selected. This procedure yields an ordered set of states that empirically traces the one-dimensional slow manifold. For the Double Angular Velocity Integration task, this procedure is extended to a two-dimensional grid of 32×3232× 32 equal bins across the joint angular space [−π,π)×[−π,π)[-π,π)×[-π,π), where the state with the minimum drift speed in each bin is selected to reconstruct the T2T^2 torus manifold. Fixed Point Detection and Topological Classification Fixed points were initially located by analyzing the one-dimensional angular flow on the identified slow manifold. We computed the angular displacement Δθ=θt+1−θt θ= _t+1- _t for each state. Candidate fixed points were identified at zero-crossings where the sign of Δθ θ flipped between adjacent states. To classify the topology of these candidates, we computed the exact Jacobian matrix =∂F()∂J= ∂ F(h) at each point using exact automatic differentiation, where F(⋅)F(·) denotes the discrete-time recurrent update function. Let |λ|max|λ|_ and |λ|min|λ|_ be the maximum and minimum moduli of the eigenvalues of J, respectively. Given a small tolerance ϵ=10−3ε=10^-3, the fixed points were classified as follows: • Stable Fixed Point: |λ|max<1−ϵ|λ|_ <1-ε. • Unstable Fixed Point: |λ|min>1+ϵ|λ|_ >1+ε. • Saddle Point: The spectrum contains both |λ|>1+ϵ|λ|>1+ε and |λ|<1−ϵ|λ|<1-ε. • Marginal Fixed Point: |λ|max≈1|λ|_ ≈ 1 (within [1−ϵ,1+ϵ][1-ε,1+ε]), indicating neutral stability along the manifold, a hallmark of perfect continuous attractors. Limit Cycle Detection and Floquet Analysis For tasks involving continuous integration (e.g., angular velocity), the network may form stable limit cycles rather than fixed points. We detected limit cycles by analyzing the tail of the autonomous trajectories. We searched for periodic orbits by comparing adjacent trajectory chunks of varying periods p. A cycle of period p was identified if the relative L2L_2 distance between consecutive chunks fell below a strict tolerance (5×10−35× 10^-3) and the oscillation amplitude was non-trivial. To assess the stability of the detected limit cycles, we performed Floquet analysis. We computed the Monodromy matrix M as the product of the Jacobians along the periodic orbit: =∏k=1pkM= _k=1^pJ_k. The limit cycle is classified as stable if the maximum modulus of the eigenvalues of M (Floquet multipliers) is strictly less than 1, and uncertain otherwise. Eigenvalue Spectrum To evaluate the time-scale separation and normal hyperbolicity of the dynamics, we extracted the full complex spectrum of the Jacobians. For all identified slow points on the manifold, we calculated the eigenvalues of the Jacobian matrix J, analyzed the kernel density estimation of their moduli and derived the spectrum on the complex plane. Uniform Norm of the Manifold Flow To theoretically bound the short-term memory error (Ságodi et al., 2024), we calculated the uniform norm (L∞L_∞ norm) of the vector field restricted to the slow manifold. Specifically, for all states ∈ℳh on the identified manifold, we computed the maximum 1-step drift magnitude in the output space: ‖Δ‖∞=max∈ℳ‖outF()−out‖2.\| \|_∞= _h \|W_outF(h)-W_outh\|_2. (A-7) Manifold Reliability Assessment For tasks on 1D ring manifolds, to automatically filter out degenerate solutions or networks that failed to form a coherent ring topology, we implemented a reliability metric. A manifold was deemed reliable if it successfully covered at least 25% of the angular bins (coverage≥0.25coverage≥ 0.25) and contained a minimum absolute number of valid slow points (e.g., ≥max(16,Nbins/8)≥ (16,N_bins/8)). Networks failing this criterion were considered to have collapsed into trivial point attractors. Figure A5 summarizes the reliability of the learned manifolds across different architectures. Figure A5: Reliability of the learned slow manifolds across architectures. A manifold is classified as reliable if it successfully forms a coherent ring topology covering a sufficient proportion of the angular space. D.3 Effective Rank of Recurrent Weight Matrices To quantify the intrinsic dimensionality of the learned recurrent dynamics and evaluate the structural inductive biases of different architectures, we computed the effective rank of their recurrent weight matrices through constructing a unified recurrent weight matrix mergedW_merged for each model by concatenating its constituent hidden-to-hidden projection matrices. Specifically, for the proposed RDNN, the merged matrix was formed by concatenating the effective excitatory and inhibitory recurrent weights after applying the positivity constraints: merged=[Softplus(raw)∥Softplus(raw)]W_merged=[Softplus(J_raw) (W_raw)]. For the GRU baseline, mergedW_merged consisted of the concatenated hidden-state weights for the reset, update, and new gates. Similarly, for the LSTM, it comprised the hidden-state weights for the input, forget, candidate, and output gates. We performed Singular Value Decomposition (SVD) on mergedW_merged to extract its singular values σ1≥σ2≥⋯≥σn≥0 _1≥ _2≥…≥ _n≥ 0. To robustly measure the dimensionality of the subspace that dominates the recurrent computations, we defined the effective rank as the 99% energy rank. This is calculated as the minimum number of singular values k required to capture 99% of the total spectral energy (i.e., the sum of squared singular values): k=minr|∑i=1rσi2∑j=1nσj2≥0.99.k= \r\; |\; _i=1^r _i^2 _j=1^n _j^2≥ 0.99 \. (A-8) Appendix E Training and Optimization Dynamics under BPTT In this section, we present the brief mathematical formulations of the Recurrent Divisive Normalization Network (RDNN), its subtractive variant, and its explicitly factorized low-rank counterpart. We derive their gradient dynamics under Backpropagation Through Time (BPTT) and provide a theoretical treatment of the emergent low-rank connectivity and optimization landscapes observed in our experiments. E.1 Discrete-Time State Formulations Let t∈ℝ+HR_t _+^H and t∈ℝ+HG_t _+^H be the excitatory and inhibitory state vectors at time step t. We denote the element-wise division by ⊘ and the element-wise product by ⊙ . The activation function is denoted by f(⋅)f(·), and t∈ℝHI_t ^H is the input current. The discrete-time update rule of the RDNN is formulated as: t _t =(−R)⊙t−1+R⊙[(f(t−1)+t)⊘(+t−1)], =(1- α_R) _t-1+ α_R [ (Jf(R_t-1)+I_t ) ( η+G_t-1) ], (A-9) t _t =(−G)⊙t−1+G⊙[f(t−1)]. =(1- α_G) _t-1+ α_G [Wf(R_t-1) ]. For the Subtractive Network, the divisive gate is replaced by an additive inhibitory subtraction, and a rectification threshold is applied to the updated state to enforce non-negative firing rates: t=max(,(−R)⊙t−1+R⊙(f(t−1)+t−t−1)).R_t= (0,\,(1- α_R) _t-1+ α_R (Jf(R_t-1)+I_t-G_t-1 ) ). (A-10) For the Low-Rank RDNN with explicit rank r<Hr<H, the recurrent weight matrices J and W are explicitly factorized into product forms: =12,=12.J=J_1J_2, =W_1W_2. (A-11) where 1,1∈ℝH×rJ_1,W_1 ^H× r and 2,2∈ℝr×HJ_2,W_2 ^r× H. E.2 Gradient Dynamics and BPTT Modulation Let ℒ=∑t=1Tℒt(t)L= _t=1^TL_t(R_t) be the total scalar loss over the sequence. Under BPTT, the gradient of the loss with respect to the recurrent excitatory matrix J is given by: ∂ℒ∂=∑t=1T∑i=1H∂ℒ∂Rt,i∂Rt,i∂. = _t=1^T _i=1^H ∂ R_t,i ∂ R_t,i . (A-12) To isolate the direct influence of the divisive normalization at step t, we compute the partial derivative of Rt,iR_t,i with respect to JjkJ_jk (where j,k∈1,…,Hj,k∈\1,…,H\): ∂Rt,i∂Jjk=(1−αR,i)∂Rt−1,i∂Jjk+αR,i[δijf(Rt−1,k)ηi+Gt−1,i−Ji,⋅f(t−1)+It,i(ηi+Gt−1,i)2∂Gt−1,i∂Jjk], ∂ R_t,i∂ J_jk=(1- _R,i) ∂ R_t-1,i∂ J_jk+ _R,i [ _ijf(R_t-1,k) _i+G_t-1,i- J_i,·f(R_t-1)+I_t,i( _i+G_t-1,i)^2 ∂ G_t-1,i∂ J_jk ], (A-13) where δij _ij is the Kronecker delta. This formulation reveals a fundamental mathematical property of the RDNN: both the direct feedforward gradient update and the indirect recurrent feedback update are scaled inversely by the dynamic divisor. Specifically, the direct gradient contribution at step t is scaled by: ∂ℒt∂Jjk∝αR,jηj+Gt−1,j. _t∂ J_jk _R,j _j+G_t-1,j. (A-14) E.3 Gradient-Driven Time-Scale Separation and Normal Hyperbolicity We now mathematically analyze how the optimization of RDNN under noisy working memory tasks naturally shapes the bimodal eigenvalue modulus spectrum (|λ|≈1|λ|≈ 1 and |λ|≈0|λ|≈ 0) and normal hyperbolicity observed in Figures 2. Let R(k)=∂k∂k−1J_R(k)= _k _k-1 be the localized Jacobian of the excitatory state update at step k, which can be expressed from the tangent linear flow of the update equation (ignoring indirect cross-coupling through G for clarity of the primary recurrent flow) as: R(k)=diag(−R)+diag(R)diag(+k−1)diag(f′(k−1)).J_R(k)=diag(1- α_R)+diag( α_R)diag ( 1 η+G_k-1 )Jdiag (f (R_k-1) ). (A-15) The gradient of the loss ℒTL_T at the end of the trial with respect to the hidden state at an earlier step t is propagated backward through time via the chain product of these Jacobians: ∂ℒT∂t=∂ℒT∂T∏k=t+1TR(k). _T _t= _T _T _k=t+1^TJ_R(k). (A-16) Convergence of Tangent Modes (|λtangent|≈1| _tangent|≈ 1) For continuous working memory tasks where information must be preserved over a long horizon T−t≫1T-t 1, the loss function ℒTL_T penalizes representation drift along the task-relevant manifold ℳM. • If the spectral radius ρ of the Jacobian product along the tangent direction tangentv_tangent of ℳM satisfies ρ(∏R(k))<1ρ ( _R(k) )<1, the gradient vanishes exponentially (∂ℒT∂t→ _T _t 0), and the network fails to learn long-term dependencies. • If ρ>1ρ>1, the gradient explodes, disrupting training stability. Consequently, gradient-based optimization mathematically drives the parameters toward a stable critical regime along the manifold, forcing the tangent modes to satisfy: R(k)tangent≈λtangenttangent,with |λtangent|≈1.J_R(k)v_tangent≈ _tangentv_tangent, | _tangent|≈ 1. (A-17) This explains the sharp, robust peak in the eigenvalue modulus density near 1.01.0 (Figures 2). Convergence of Normal Modes (|λnormal|→0| _normal|→ 0) During training, we inject Gaussian state noise t∼(,σ2) ξ_t (0,σ^2I) into the hidden dynamics. In the linear approximation, the forward propagation of the off-manifold perturbation covariance t _t along the normal (orthogonal) directions normalv_normal to the manifold ℳM is governed by the discrete-time Lyapunov-like update: t=R(t)t−1R(t)T+σ2. _t=J_R(t) _t-1J_R(t)^T+σ^2I. (A-18) To minimize the expected loss [ℒ]E[L], the optimizer must suppress the cumulative propagation of this state noise off the manifold, preventing representation diffusion. This imposes a selective gradient bias that drives the eigenvalues corresponding to the normal directions toward zero: R(k)normal≈λnormalnormal,with |λnormal|→0.J_R(k)v_normal≈ _normalv_normal, | _normal|→ 0. (A-19) This leads to the massive accumulation of eigenvalues near 0 (Figures 2). For configurations utilizing the ReLU activation, there exists f′(Rt−1,i)=0f (R_t-1,i)=0 for any neuron i below the threshold. This structurally forces the diagonal derivative matrix diag(f′(t−1))diag(f (R_t-1)) to be sparse, setting the corresponding columns of the Jacobian R(k)J_R(k) to exactly zero and contributing to the massive eigenvalue peak at 0. The Role of Divisive Scaling in Facilitating Spectral Separation While standard gating RNNs (like GRU/LSTM) struggle to perfectly separate these eigenvalues due to rigid linear coupling, the RDNN’s divisive normalization provides a unique mathematical template. When a noisy perturbation ξ pushes the state off the manifold, it temporarily inflates the overall network activity. The inhibitory population instantly tracks this expansion, increasing Gk−1,iG_k-1,i. Through the term diag(+k−1)diag ( 1 η+G_k-1 ), this dynamic increase in the divisor mathematically scales down the localized Jacobian R(k)J_R(k) for those off-manifold modes, dragging their eigenvalues towards 0. Conversely, when the state is unperturbed and resides on the manifold, the divisor remains balanced, preserving |λtangent|≈1| _tangent|≈ 1. This state-dependent gain control provides the exact mathematical mechanism that allows the gradient descent optimizer to easily sculpt the strict spectral gap and robust normal hyperbolicity observed in our RDNN models. E.4 Local Gradient Scaling and Its Empirical Relation to Low-Rank Connectivity The inverse scaling of the gradients by the dynamic denominator (+t−1)( η+G_t-1) acts as an activity-dependent gradient attenuator during gradient descent (Gunasekar et al., 2017; Woodworth et al., 2020). • Gradient Attenuation in High-Activity Regimes: When the network is highly active along certain directions in the state space, the inhibitory pool t−1G_t-1 is driven strongly via f(t−1)Wf(R_t-1), resulting in large values of Gt−1,jG_t-1,j. Consequently, the gradient updates for the corresponding rows Jj,⋅J_j,· are heavily suppressed by the factor 1/Gt−1,j1/G_t-1,j. • Selective Parameter Updates: Weight updates are concentrated almost exclusively on neurons that are marginally active (where Gt−1,jG_t-1,j is small but Rt−1,kR_t-1,k is non-zero). This creates a competitive winner-take-all environment in the parameter space. • Connection to Low-Rank Biases: While local gradient scaling does not mathematically equate to minimizing a formal rank metric or singular-value entropy, restricting the directions in parameter space that receive significant gradient updates can guide the network toward a lower-dimensional subspace. In deep learning theory, certain multiplicative or gated update rules have been associated with implicit biases toward sparsity or low-rank structures (Arora et al., 2019). Here, the dynamic scaling suppresses updates along highly active directions, which heuristically encourages the self-compression observed in our empirical rank measurements (as verified empirically in Table 1 and Figure 3A), without requiring explicit structural constraints. E.5 Optimization and Dynamical Landscapes of Ablated Networks E.5.1 Subtractive Normalization and Attractor Shattering We now mathematically analyze why the Subtractive Network sustains a stable ring attractor in the autonomous memory-guided saccade task but shatters into discrete point attractors under the input-driven angular velocity integration task, and why it consistently exhibits an elevated effective rank. In the discrete-time formulation, the Jacobian of the Subtractive Network at step t is given by: sub(t)=t[diag(−R)+diag(R)diag(f′(t−1))],J_sub(t)= _t [diag(1- α_R)+diag( α_R)Jdiag (f (R_t-1) ) ], (A-20) where t=diag(Θ((−R)⊙t−1+R⊙(f(t−1)+t−t−1)))∈0,1H×H _t=diag ( ((1- α_R) _t-1+ α_R (Jf(R_t-1)+I_t-G_t-1) ) )∈\0,1\^H× H is a diagonal gating matrix defined by the Heaviside step function Θ(⋅) (·) originating from the state rectification boundary. A continuous ring attractor requires a smooth 1-parameter family of fixed points satisfying the marginal stability condition, where the maximum eigenvalue modulus of the Jacobian is strictly bounded at unity (|λ|max≈1|λ|_ ≈ 1). 1. Autonomous Regime (t=I_t=0): The gating matrix t _t remains static and spatially uniform over the persistent representation. Under these zero-input conditions, the network only needs to satisfy the marginal stability condition at a single unperturbed operating point. This is easily achieved through standard linear balancing of recurrent excitation and subtraction, preserving the ring topology. 2. Input-Driven Regime (t≠I_t 0): As the external velocity input tI_t continuously fluctuates over time, the term within the Heaviside function is directly modulated. This forces the diagonal gating matrix t _t to undergo discrete, discontinuous state-switching. Consequently, as external inputs toggle the binary states in t _t, the effective recurrent connectivity undergoes abrupt structural changes. Without a dynamic divisor to smoothly normalize the Jacobian, its eigenvalues experience uncompensated discrete jumps. This mathematically violates the marginal stability constraint, triggering a sequence of saddle-node bifurcations that shatter the smooth ring manifold into discrete, highly stable point-attractor basins (sinks) separated by saddles (Figure 4A). To understand why the Subtractive Network maintains a consistently higher effective rank than the RDNN (as illustrated in Figure 5 and Table A1), we analyze the BPTT gradient updates of the subtractive recurrent weights. Denoting the local loss at step t by ℒtL_t, the direct gradient of ℒtL_t with respect to the subtractive recurrent weights JjkJ_jk is given by: ∂ℒt∂Jjk=et,j⋅αR,j⋅Θt,j⋅f(Rt−1,k), _t∂ J_jk=e_t,j· _R,j· _t,j· f(R_t-1,k), (A-21) where et,j=∂ℒt∂Rt,je_t,j= _t∂ R_t,j is the backpropagated error, and Θt,j=Θ((1−αR,j)Rt−1,j+αR,j(j,⋅f(t−1)+It,j−Gt−1,j))∈0,1 _t,j= ((1- _R,j)R_t-1,j+ _R,j (J_j,·f(R_t-1)+I_t,j-G_t-1,j ) )∈\0,1\ (A-22) is the binary activation state of neuron j after rectification. Contrast this with the gradient updates in the RDNN, where the updates are scaled continuously by the multiplicative divisor: ∂ℒt∂Jjk∝αR,jηj+Gt−1,j. _t∂ J_jk _R,j _j+G_t-1,j. (A-23) In the RDNN, the dynamic denominator (ηj+Gt−1,j)( _j+G_t-1,j) acts as an activity-dependent spectral compressor. When the network is active, Gt−1,jG_t-1,j is large, attenuating the weight updates along those dimensions and restricting parameter changes to a sparse, low-rank subspace. In the Subtractive Network, however, the additive inhibition Gt−1,jG_t-1,j only modulates the gradient via the binary gate Θt,j _t,j. For all active units along the persistent representation (Θt,j=1 _t,j=1), the magnitude of the gradient update is completely decoupled from the scale of the inhibition Gt−1,jG_t-1,j. Without this continuous multiplicative scaling, the optimization path lacks the activity-dependent gradient attenuation mechanism provided by divisive normalization. Consequently, the weight updates tend to propagate across a wider, less constrained subspace of ℝH×HR^H× H, leading to the elevated effective rank observed in our quantitative scaling analysis. BPTT Optimization and Dimensional Compensation This gradient-level decoupling also elucidates the discrepancy between the topological fragmentation and the seemingly similar eigenvalue modulus spectra in the angular velocity integration task (Figure 4A and A15). To minimize the integration loss over a long horizon T, BPTT must prevent gradient vanishing along the active trajectory by forcing the task-relevant eigenvalues to satisfy |λ|≈1|λ|≈ 1. In the RDNN, the dynamic denominator + η+G acts as a spectral compressor, driving orthogonal (normal) directions to |λnormal|→0| _normal|→ 0 (normal hyperbolicity, see Appendix E.3). This confines the slow modes strictly to the low-dimensional tangent space of the ring, allowing saddles to act as efficient slow channels for continuous integration (Sussillo and Barak, 2013). In the Subtractive Network, the continuous manifold shatters into a high density of discrete stable point-attractors (sinks). To bridge these discrete wells and prevent the state representation from decaying during input-driven transitions, BPTT is forced to preserve information across a much wider subspace. While the post-update hard-thresholding of the ReLU function still enforces structural zero eigenvalues for strictly silent neurons, as t _t zeros out their entire corresponding rows in the Jacobian, accounting for the visible density peak near |λ|≈0|λ|≈ 0 (Glorot et al., 2011), the subtractive gradient lacks the multiplicative scaling to dynamically contract the active normal dimensions. Therefore, the optimizer cannot efficiently suppress non-tangent directions among the active subpopulation (Arora et al., 2019; Gunasekar et al., 2017). Consequently, while both networks exhibit a prominent density peak near |λ|≈1|λ|≈ 1 to satisfy the BPTT memory constraint, the Subtractive Network achieves this by recruiting a higher number of physical dimensions to construct redundant slow pathways across the fragmented landscape. This dimensional compensation is mathematically manifested as the significantly elevated effective rank observed in Figure 5. E.5.2 Spectral Alignment and Hyper-Compression in Low-Rank RDNN Here we investigate the dynamics underlying the observation that the trained Explicit Low-Rank RDNN exhibits an effective rank that is significantly lower than both the permitted bottleneck rank r and the unconstrained RDNN’s effective rank. When the recurrent weights are factorized as =12J=J_1J_2 and =12W=W_1W_2, the gradients of the BPTT loss ℒL with respect to the factorized matrices are: ∇1ℒ=(∇ℒ)2T,∇2ℒ=1T(∇ℒ). _J_1L=( _JL)J_2^T, _J_2L=J_1^T( _JL). (A-24) Under gradient descent, this coupled product formulation induces strong spectral alignment (or co-adaptation) between the column space of 1J_1 and the row space of 2J_2. In deep matrix factorization theory (Arora et al., 2019; Gunasekar et al., 2017), this multiplicative gradient flow acts as an accelerator for singular value decay. Specifically, updates to the dominant singular vectors are mutually reinforced by the transpose factors: Δσi()∝σi(1)⋅σi(2). _i(J) _i(J_1)· _i(J_2). (A-25) In the context of the RDNN, these factorized gradient dynamics appear to interact with and amplify the local gradient-scaling effects of divisive normalization. This combined effect correlates with the network exhibiting a pronounced hyper-compression, where it effectively discards many available dimensions of the bottleneck r and concentrates 99% of its spectral energy into an exceptionally tight, low-dimensional subspace, empirically observed to be of rank ≈2∼4≈ 2 4. E.5.3 High-Rank Convergence Anomaly in Explicit Low-Rank Factorization For the Explicit Low-Rank RDNN, the gradient of the loss with respect to the factorized matrices (e.g., 1J_1 and 2J_2) is given by: ∂ℒ∂1=∂ℒ∂2T,∂ℒ∂2=1T∂ℒ∂. _1= J_2^T, _2=J_1^T . (A-26) While this factorization guarantees rank()≤rrank(J)≤ r, it introduces severe non-convexities into the optimization landscape: • Scaling Symmetries and Flat Valleys: For any invertible matrix ∈ℝr×rP ^r× r, the transformation (1,−12)(J_1P,P^-1J_2) leaves the product J invariant. This continuous symmetry creates infinitely many flat directions (valleys) in the loss landscape. As dictated by the cross-dependent gradients in Eq. A-26, drifting along these valleys can induce severe scale imbalances between the factors (e.g., ‖2‖→0\|J_2\|→ 0). When this occurs, the gradient for the other factor (∂ℒ∂1∝2T _1 _2^T) vanishes and stalls the optimization process (Li et al., 2018; Zhao et al., 2026). • The High-Rank Convergence Anomaly: When the permitted explicit bottleneck rank r is unnecessarily large (e.g., r=32r=32) relative to the task’s true low-dimensional manifold (which we showed requires only rank≈2∼4rank≈ 2 4), the number of redundant symmetries and saddle points scales quadratically with r (Li et al., 2019; Ge et al., 2017). This over-parameterization of the bottleneck introduces a vast number of degenerate saddle points and flat plateaus, severely hindering convergence (Xiong et al., 2024) (as observed in the high variance of r=32r=32 in Figure 3B). A tight bottleneck (r=8r=8) minimizes these redundant degrees of freedom, smoothing out the optimization landscape and enabling rapid convergence. E.5.4 Linear Parameterization and the Absence of Anomaly in Unconstrained RDNN We analyze why the unconstrained (full-rank) RDNN does not suffer from the convergence anomaly as its hidden size H increases, unlike the explicitly factorized network. In the unconstrained RDNN, the parameters are represented directly by the flat recurrent matrices ,∈ℝH×HJ,W ^H× H, avoiding any product-form parameterization. 1. Absence of Scaling Symmetries: The parameter-to-weight mapping is a trivial identity. Thus, the unconstrained formulation does not introduce the continuous non-convex scaling symmetries (1,−12)(J_1P,\,P^-1J_2) that create flat valleys and zero-gradient directions (Dinh et al., 2017; Li et al., 2019). 2. Landscape Smoothing via Over-Parameterization: In deep learning theory, increasing the dimensionality H of a directly parameterized weight matrix J leads to standard over-parameterization. This over-parameterization mathematically smoothes the loss landscape by creating highly connected, high-dimensional paths that eliminate bad local minima and degenerate saddle points (Du et al., 2019; Cooper, 2018). Consequently, larger unconstrained RDNNs (e.g., H=256H=256) train more rapidly and stably than smaller ones, whereas larger explicit bottlenecks r in factorized networks exacerbate the non-convex factorization pathology. Appendix F Extended Slow Manifold Analysis Figure A6: Complex eigenvalue spectra of the Jacobians evaluated on the slow manifolds. In the angular velocity integration task, the RDNN exhibits a highly structured, tripartite spectrum consisting of a dense cluster along the real axis (governing fast dissipation), off-axis conjugate pairs in the interior (stabilizing transient transitions), and boundary conjugate pairs near the unit circle (facilitating continuous rotational dynamics). In the memory-guided saccade task, the RDNN suppresses these imaginary components, aligning its slow modes along the real axis to maintain static memory. In contrast, GRU and LSTM exhibit diffuse spectral clouds lacking a clear spectral gap or structured rotational modes. Figure A7: The uniform norm quantifies the maximum drift speed along the manifold; a lower value indicates a flatter energy landscape and a higher-fidelity continuous attractor. The RDNN maintains a substantially lower uniform norm than the baselines in the integration task (left), and achieves near-zero drift (on the order of 10−310^-3 to 10−410^-4) in the autonomous saccade task (right). Markers indicate the mean across reliable instances, with error bars denoting the standard deviation. Figure A8: Instance-level equilibrium classification for the angular velocity integration task. The heatmap displays the absolute count of different dynamical states (fixed points and limit cycles) for each trained instance (seed) across various configurations. While GRU and LSTM instances are strictly dominated by alternating stable fixed points and saddles (indicating discretized state spaces), RDNN instances exhibit a more fluid topology, frequently developing marginal fixed points and stable limit cycles to support continuous integration. Figure A9: Instance-level equilibrium classification for the memory-guided saccade task. The heatmap details the dynamical states of the RDNN across all trained seeds and configurations. The network consistently converges to solutions overwhelmingly dominated by marginal fixed points. This instance-level consistency confirms that divisive normalization provides a strong inductive bias for forming approximate continuous attractors for autonomous memory maintenance. Appendix G Extended Analysis of Explicit Low-Rank Parameterization Figure A10: Test NMSE of explicitly factorized RDNNs. Normalized Mean Squared Error (NMSE, in dB) on the test set of 1024 trajectories for the H=128H=128 RDNN with explicitly factorized recurrent weights (J and W) across bottleneck ranks r∈8,16,24,32r∈\8,16,24,32\, compared with the unconstrained Full-Rank RDNN. Left: Angular Velocity Integration task. Right: Memory-guided Saccade task. Error bars denote the standard deviation across different random seeds. Even under an extremely tight bottleneck (e.g., r=8r=8), the RDNN maintains highly competitive task performance, demonstrating the intrinsic low-dimensional nature of the continuous working memory representation. Figure A11: Self-compressed effective rank of explicitly factorized RDNNs. The post-training 99% energy effective rank of the H=128H=128 RDNN with explicitly factorized weights across activation functions (relu and tanh). Left: Angular Velocity Integration task. Right: Memory-guided Saccade task. Despite having an explicit parameter bottleneck of r∈8,16,24,32r∈\8,16,24,32\, the trained models consistently compress their effective dimensionality to approximately 3.5∼4.03.5 4.0 (Angular Velocity Integration) and ∼2.0 2.0 (Memory-guided Saccade), indicating that the network actively discards redundant dimensions to align with the task’s low-dimensional manifold. Figure A12: Jacobian eigenvalue modulus density of explicit low-rank RDNNs. Kernel density estimate (KDE) of the eigenvalue modulus (|λ||λ|) evaluated on the slow manifold for the H=128H=128 explicitly factorized RDNNs across ranks r∈8,16,24,32r∈\8,16,24,32\ and activations (relu as dashed, tanh as solid lines). Left: Angular Velocity Integration task. Right: Memory-guided Saccade task. While the bimodal structure (representing slow-fast separation) is largely preserved, some configurations in the saccade task exhibit minor eigenvalue leakage in the intermediate region, reflecting a slight degradation of the strict spectral gap due to explicit low-rank constraints. Figure A13: Complex plane eigenvalue spectra of explicit low-rank RDNNs. Scatter plots of the complex Jacobian eigenvalues evaluated on the slow manifold across different bottleneck ranks r∈8,16,24,32r∈\8,16,24,32\ (rows represent different explicit ranks). Left: Angular Velocity Integration task. Right: Memory-guided Saccade task. The tripartite spectral structure in the integration task (real-axis clustering, interior and boundary conjugate pairs) and the strict real-axis alignment in the saccade task are largely preserved with explicit factorization, though higher ranks introduce minor off-axis scattering on the saccade manifold. Appendix H Extended Ablation Analysis Table A1: The 99% energy effective rank (mean±stdmean over multiple seeds) across RDNN and Subtractive Network. Model Activation Angular Velocity Integration Memory-guided Saccade H=64H=64 H=128H=128 H=256H=256 H=64H=64 H=128H=128 H=256H=256 RDNN relu 16.8±2.416.8± 2.4 24.0±2.524.0± 2.5 38.0±6.738.0± 6.7 19.2±2.519.2± 2.5 27.0±7.227.0± 7.2 31.0±5.031.0± 5.0 tanh 17.2±1.917.2± 1.9 22.6±2.622.6± 2.6 35.4±5.335.4± 5.3 16.0±2.016.0± 2.0 22.6±3.522.6± 3.5 28.3±5.628.3± 5.6 Subtractive Network relu 24.0±1.024.0± 1.0 37.6±1.537.6± 1.5 51.0±2.751.0± 2.7 24.0±2.524.0± 2.5 29.4±4.129.4± 4.1 40.25±3.940.25± 3.9 tanh 25.2±2.325.2± 2.3 40.2±3.640.2± 3.6 52.2±3.352.2± 3.3 21.0±2.221.0± 2.2 29.6±2.229.6± 2.2 37.0±5.737.0± 5.7 Figure A14: Performance comparison (Test NMSE) of divisive (RDNN) and subtractive networks. Test set Normalized Mean Squared Error (NMSE, in dB) of 1024 trajectories on both tasks across different hidden sizes and activation functions. While the Subtractive Network achieves highly competitive or marginally lower test NMSE compared to the RDNN, its underlying dynamical representations with external inputs are highly discretized (as shown in Figures 4 and A17), contrasting with the fluid, low-rank manifolds learned by the RDNN. Error bars denote the standard deviation across different random seeds. Figure A15: Eigenvalue modulus density spectra of divisive and subtractive networks. Kernel density estimates (KDE) of the Jacobian eigenvalue moduli (|λ||λ|) evaluated on the slow manifolds. Left (Angular Velocity Integration): Both networks exhibit a prominent peak near |λ|≈1|λ|≈ 1 to satisfy the BPTT memory constraint for continuous integration. Right (Memory-guided Saccade): Both networks display nearly identical, sharply separated bimodal spectra with a strict spectral gap, confirming that both mechanisms can achieve marginal stability in the absence of external input fluctuations. Figure A16: Complex plane eigenvalue spectra of the Subtractive Network. Scatter plots of the complex Jacobian eigenvalues evaluated on the slow manifold of the Subtractive Network. Left: Angular Velocity Integration task. Right: Memory-guided Saccade task. Unlike the RDNN (cf. Figure A6), the Subtractive Network’s eigenvalues remain heavily concentrated along the real axis even during continuous integration (left), lacking the boundary-hugging complex conjugate pairs required for smooth, marginally stable rotational dynamics. Figure A17: The heatmaps display the absolute count of different dynamical states (fixed points and limit cycles) of the Subtractive Network across all trained random seeds and configurations for the angular velocity integration task (left) and the memory-guided saccade task (right). In the integration task, the Subtractive Network consistently converges to a discretized state space dominated by stable fixed points and saddles across every single seed, demonstrating that attractor shattering is a structural consequence of additive subtraction. In the autonomous saccade task (right), it robustly maintains marginal fixed points, matching the qualitative behavior of the RDNN under zero-input conditions. Appendix I Robustness Analysis To evaluate the stability of the learned representations, we test the RDNN’s performance under varying intensities of state noise σ∈[0.0,0.3]σ∈[0.0,0.3] and analyze its sensitivity to the training learning rate. As illustrated in Figure A18, the test NMSE (dB) exhibits a gradual and stable degradation as the noise level increases across all evaluated hidden sizes (H∈64,128,256H∈\64,128,256\), with larger networks consistently maintaining lower error rates on both tasks. This stable performance degradation under noise is consistent with the normal hyperbolicity of the learned slow manifolds, which contributes to suppressing off-manifold perturbations. Additionally, our learning rate sensitivity analysis for the H=128H=128 configuration (Table A2) shows that the network is sensitive to optimization parameters; while a learning rate of LR=5×10−3LR=5× 10^-3 yields the best test performance on both tasks, lowering the learning rate to LR=10−3LR=10^-3 leads to a notable decrease in performance. Figure A18: Test NMSE (dB) of the RDNN across different hidden dimensions (H) as a function of state noise intensity σ. Table A2: Test NMSE (dB) of the H=128H=128 RDNN trained with different learning rates. Learning Rate Angular Velocity Integration Memory-guided Saccade LR=10−2LR=10^-2 -14.5 ± 0.8 -27.3 ± 1.4 LR=5×10−3LR=5× 10^-3 -15.6 ± 0.7 -29.8 ± 1.0 LR=10−3LR=10^-3 -10.2 ± 0.5 -21.0 ± 1.1 Appendix J Rank of Component Matrices We computed the 99%99\% energy effective rank of each individual component matrix in isolation (e.g., the excitatory matrix J and inhibitory matrix W individually for the RDNN and the Subtractive Network, and the individual gating matrices for the GRU and LSTM baselines). As shown in Figure A19, the effective ranks of these individual matrices exhibit the same scaling behavior as their merged counterparts. Individual RDNN weights maintain a compressed, sub-linear scaling with H, whereas individual gating weights of the GRU and LSTM baselines scale linearly, nearly matching the physical hidden size. Similarly, individual recurrent weights of the Subtractive Network scale sub-linearly with H, though their ranks remain slightly higher than those of the RDNN across configurations. Furthermore, the effective rank of any individual component matrix is similar in magnitude to the merged rank of the concatenated matrix reported in Table 1. This consistency indicates that our findings regarding the low-rank properties of the RDNN are robust and independent of whether the recurrent weights are analyzed in a merged or component-wise manner. Figure A19: The 99%99\% energy effective rank of individual gate/component weights vs. hidden sizes (H=64,128,256H=64,128,256) for the Angular Velocity Integration task (A) and the Memory-guided Saccade task (B) across different architectures. Appendix K Supplementary Experiments on ORGaNICs To investigate whether the robust slow manifold dynamics observed in the RDNN are generalizable properties of divisive normalization (DN) as a canonical computational principle, we conduct supplementary experiments on ORGaNICs (Rawat et al., 2024), a biologically plausible recurrent circuit model that dynamically implements DN. As shown in Figures A20 and A21, ORGaNICs successfully converges and achieves highly competitive test performance. Crucially, the topological distributions and eigenvalue modulus spectra (Figure A22) demonstrate that ORGaNICs forms slow manifolds with dynamical landscapes highly similar to the RDNN, characterized by the dominance of marginal fixed points in the autonomous regime and saddles under external drive. This robust behavior is theoretically aligned with the proof in Rawat et al. (2024) that the DN-equipped recurrent loop is mathematically equivalent to coupled damped harmonic oscillators with guaranteed local asymptotic stability, thereby empirically reinforcing the canonical role of divisive normalization in stabilizing continuous working memory. However, we also identify key structural discrepancies: the effective rank of ORGaNICs’ recurrent weights scales linearly and remains substantially higher than that of the RDNN (Figure A23), and its spectrum lacks a distinct dissipative peak near zero. Unlike our RDNN, which represents a minimal model designed to strictly isolate the algebraic role of division, ORGaNICs incorporates distinct biophysical mechanisms, such as implicit multiplier-based gating for normalization and sluggish inhibitory time constants. A precise mathematical characterization of how these diverse physiological designs modulate the representation dimensionality and spectral structure remains an important open question for future study. Figure A20: Training loss curves of the ORGaNICs model. Error bands represent the standard deviation across different random seeds. Figure A21: Performance comparison (Test NMSE) between ORGaNICs and the RDNN. Figure A22: Steady-state topology (left) and eigenvalue modulus density (right) of the ORGaNICs model. Figure A23: Effective rank of individual recurrent components in the ORGaNICs model. Appendix L Supplementary Experiments on Neural ODEs To address continuous-time baselines, we evaluate a Neural Ordinary Differential Equation (Neural ODE) model (Chen et al., 2018) implemented via the torchdiffeq framework. At each step t with input tx_t, the hidden state (τ)h(τ) evolves over the interval τ∈[0,0.1]τ∈[0,0.1] according to: d(τ)dτ=2σ(1[(τ)∥t]+1)+2 dh(τ)dτ=W_2σ (W_1[h(τ) _t]+b_1 )+b_2 (A-27) where [⋅∥⋅][· ·] denotes concatenation, and σ(⋅)σ(·) represents the Tanh or ReLU activation. The state is integrated using a 4th-order Runge-Kutta (RK4) solver, and decoded linearly as t=outt+outy_t=W_outh_t+b_out. While the Neural ODE stably converges to a low loss during training (Figure A24), its test generalization remains limited, exhibiting a higher test NMSE compared to the proposed RDNN (Figure A25). This is explained by its underlying dynamical properties: the effective rank of its recurrent components scales linearly with the hidden size (Figure A27), indicating a lack of the mechanisms that contribute to the lower effective rank in our proposed model. Furthermore, its steady states are heavily dominated by unstable or marginally stable limit cycles (Uncertain LCs), and its eigenvalue spectrum lacks a peak near 0 (Figure A26), reflecting a collapse of normal hyperbolicity. Furthermore, we observe severe physical anomalies during testing (Figure A28), such as rapid state decay during the reporting phase (left) or continuous state drift during the zero-input delay period (right). This suggests that continuous-time modeling alone is insufficient to stabilize continuous representations, highlighting the importance of specific architectural inductive biases, such as divisive normalization. Figure A24: Training loss curves of the Neural ODE. Error bands represent the standard deviation across different random seeds. Figure A25: Performance comparison (Test NMSE) between the Neural ODE and the RDNN. Figure A26: Steady-state topology and eigenvalue modulus density of the Neural ODE. Top: Angular Velocity Integration. Down: Memory-guided Saccade. Left: Statistical counts of state/cycle types across configurations. The steady states are heavily dominated by unstable or marginally stable limit cycles (Uncertain LCs). Right: Kernel density estimates (KDE) of the Jacobian eigenvalue modulus (|λ||λ|). Figure A27: Effective rank of individual recurrent components in Neural ODE and RDNN. For the Neural ODE (left), the effective ranks of both raw parameter matrices (1W_1 and 2W_2) scale linearly with H, while RDNN (right) maintains a flat, sub-linear scaling. Figure A28: Trajectory anomalies of the Neural ODE on the memory-guided saccade task. Left: The model exhibits rapid state decay and instability immediately after the Go Cue (black dashed line). Right: The network fails to maintain the memory during the long delay period, exhibiting continuous state drift prior to the Go Cue. Appendix M Supplemantary Experiments on Transformer Models To evaluate attention-based baselines, we implement task-specific wrappers for the Autoformer architecture (Wu et al., 2021). The core model consists of an Encoder with elayers=2e_layers=2 and a Decoder with dlayers=1d_layers=1. We set the model dimension dmodeld_model equal to the hidden size H, and utilize nheads=4n_heads=4 attention heads, a series decomposition moving average window of 2525, and a dropout rate of 0.10.1 with GeLU activation. Except for the learning rate, which we sweep over 5×10−3,10−3,5×10−4,10−4\5× 10^-3,10^-3,5× 10^-4,10^-4\ to evaluate convergence, all other training hyperparameters are kept identical to the RDNN training setup. As illustrated in Figure A29, the Autoformer exhibits severe convergence stagnation, with its training loss plateauing at highly elevated levels across all tested learning rates. This optimization failure might be attributed to Autoformer’s core design assumptions: its progressive series decomposition blocks, which utilize average pooling, tend to smooth out and smear the non-stationary step transitions and cumulative updates required for path integration, while its auto-correlation attention is structurally biased toward discovering global periodicity rather than local, causal temporal integration. Furthermore, we note that the continuous-time dynamical systems analysis is mathematically inapplicable to Transformer-based models. Because these feedforward attention architectures lack a recurrent state equation t=F(t−1,t)h_t=F(h_t-1,x_t), they do not possess a recursive hidden state space in which localized fixed points and tangent flows can be defined. Figure A29: Training loss curves of the Autoformer under different learning rates. Error bands represent the standard deviation across different random seeds. Appendix N Long-Horizon Stability We evaluate the long-horizon stability of the learned manifolds by testing the trained networks far beyond their training horizons. Specifically, we evaluate the models on the angular velocity integration task (Figure A30, left) and the memory-guided saccade task (Figure A30, right) up to 2048 time steps. The long-term error dynamics demonstrate a canonical trade-off in attractor networks, closely aligning with the theoretical framework of Ságodi et al. (2024). While the uniform norm of the flow field bounds the memory error in the short behavioral timescale, the long-term asymptotic error is dominated by the topology of the dynamics. In the integration task (Figure A30, left), the RDNN’s error grows continuously over 2000 steps due to the diffusion along its flat, continuous manifold. Conversely, the Subtractive Network and gating RNNs shatter into discrete basins; although this discretization introduces a small rounding error, these discrete basins act as an error-correcting mechanism (the lock-in effect) that bounds the long-term error (Koulakov et al., 2002; Brody et al., 2003). In the Saccade task (Figure A30, right), both the RDNN and the Subtractive Network display nearly identical error degradation up to 2000 steps. This is consistent with our findings in Section 6.2, confirming that under zero-input conditions, both models form functionally equivalent continuous ring attractors and thus suffer from the same long-term diffusion. Figure A30: Long-horizon prediction error and attractor diffusion across task regimes. Error bands represent the standard deviation across different runs. Left: Angular Velocity Integration task. Right: Memory-guided Saccade task. Appendix O Double Angular Velocity Integration Task To investigate whether our findings generalize to higher-dimensional, coupled continuous systems, we evaluate the networks on the Double Angular Velocity Integration task, which requires the simultaneous integration of two independent circular variables on a two-dimensional torus manifold. As shown in Figure A31, the proposed RDNN maintains a superior performance compared to standard baselines, demonstrating its robustness in high-dimensional continuous tracking. This performance advantage is explained by the underlying topological and spectral properties of the learned manifolds. Our analysis (Figure A32, top) reveals that the RDNN is predominantly governed by stable and uncertain limit cycles, reflecting a fluid, continuous flow on the torus surface. In contrast, the GRU and LSTM are overwhelmingly dominated by stable fixed points, indicating that they solve the double integration by discretizing the 2D torus into a rigid grid of point-attractor basins, while the Neural ODE is dominated by uncertain limit cycles on a structurally unstable, thick manifold. This is further corroborated by their complex eigenvalue spectra (Figure A32, bottom) and effective ranks (Figure A33): the RDNN exhibits a highly structured, low-rank spectrum, while the other networks display diffuse, high-rank spectral clouds. These results indicate that the biological divisive normalization bottleneck is critical for stabilizing higher-dimensional continuous attractors. Figure A31: Test set NMSE on the Double Angular Velocity Integration task. Figure A32: Steady-state topology and complex eigenvalue spectra for the Double Angular Velocity Integration task. Top: Statistical counts of state/cycle types across configurations. Bottom: Jacobian eigenvalues plotted on the complex plane. Figure A33: Effective rank of component matrices on the Double Angular Velocity Integration task. Appendix P Parameter Count and Computational Cost To evaluate the computational and memory overhead of the networks, we analyze their analytical parameter counts, asymptotic complexities, and leading-order operations. Table A3 details the parameter scaling of the RDNN, which behaves quadratically as (H2)O(H^2). A comparative evaluation across architectures under various hidden sizes H∈64,128,256H∈\64,128,256\ (assuming Din=1D_in=1 and Dout=2D_out=2) is summarized in Table A4. Because all evaluated architectures are state-space networks, they share an identical single-step asymptotic time complexity of (H2)O(H^2) and a transient active memory space complexity of (H)O(H) per step during inference. However, their exact leading-order operations differ: the RDNN requires 4H24H^2 recurrent floating-point operations (FLOPs) per step to compute its excitatory and inhibitory updates, which is lower than the TaskGRU (6H26H^2), TaskLSTM (8H28H^2), and ORGaNICs (12H212H^2, due to its multiple gating and modulation matrices). Furthermore, because the NeuralODE requires four evaluations of its derivative network per step under the 4th-order Runge-Kutta (RK4) integration, its step-wise computation scales as 16H216H^2 FLOPs. Table A3: Analytical parameter breakdown of the RDNN. Parameter Component Parameter Count Input Projection (in,inW_in,b_in) H⋅Din+H· D_in+H Output Projection (out,outW_out,b_out) Dout⋅H+DoutD_out· H+D_out Recurrent Excitatory (J) H2H^2 Recurrent Inhibitory (W) H2H^2 Semi-saturation Vector ( η) H Update Rates (R,G α_R, α_G) 2H2H Total Parameter Count 2H2+H(Din+Dout+4)+Dout2H^2+H(D_in+D_out+4)+D_out Table A4: Computational complexity and parameter comparison across architectures. Architecture Hidden Size (H) Parameter Count Time Complexity Leading-Order FLOPs (per step) RDNN H=64H=64 8.6 K (H2)O(H^2) 4H24H^2 H=128H=128 33.7 K H=256H=256 132.9 K GRU H=64H=64 12.8 K (H2)O(H^2) 6H26H^2 H=128H=128 50.2 K H=256H=256 198.7 K LSTM H=64H=64 17.0 K (H2)O(H^2) 8H28H^2 H=128H=128 66.8 K H=256H=256 264.7 K NeuralODE (RK4) H=64H=64 8.5 K (H2)O(H^2) 16H216H^2 H=128H=128 33.4 K H=256H=256 132.4 K ORGaNICs H=64H=64 29.3 K (H2)O(H^2) 12H212H^2 H=128H=128 115.8 K H=256H=256 461.1 K