Unconditional stability of a recurrent neural circuit implementing divisive normalization

Authors: shivang  rawat
david j heeger
stefano  martiniani
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 02:24:45

Abstract: Stability in recurrent neural models poses a significant challenge, particularly in developing biologically plausible neurodynamical models that can be seamlessly trained. Traditional cortical circuit models are notoriously difficult to train due to expansive nonlinearities in the dynamical system, leading to an optimization problem with nonlinear stability constraints that are difficult to impose. Conversely, recurrent neural networks (RNNs) excel in tasks involving sequential data but lack biological plausibility and interpretability. In this work, we address these challenges by linking dynamic divisive normalization (DN) to the stability of "oscillatory recurrent gated neural integrator circuits'' (ORGaNICs), a biologically plausible recurrent cortical circuit model that dynamically achieves DN and that has been shown to simulate a wide range of neurophysiological phenomena. By using the indirect method of Lyapunov, we prove the remarkable property of unconditional local stability for an arbitrary-dimensional ORGaNICs circuit when the recurrent weight matrix is the identity. We thus connect ORGaNICs to a system of coupled damped harmonic oscillators, which enables us to derive the circuit's energy function, providing a normative principle of what the circuit, and individual neurons, aim to accomplish. Further, for a generic recurrent weight matrix, we prove the stability of the 2D model and demonstrate empirically that stability holds in higher dimensions. Finally, we show that ORGaNICs can be trained by backpropagation through time without gradient clipping/scaling, thanks to its intrinsic stability property and adaptive time constants, which address the problems of exploding, vanishing, and oscillating gradients. By evaluating the model's performance on RNN benchmarks, we find that ORGaNICs outperform alternative neurodynamical models on static image classification tasks and perform comparably to LSTMs on sequential tasks.

Main paper: 
 introduction: 
 Deep neural networks (DNNs) have found widespread use in modeling tasks from experimental systems neuroscience. The allure of DNN-based models lies in their ease of training and the flexibility they offer in architecting systems with desired properties [1] [2] [3] . In contrast, neurodynamical models like the Wilson-Cowan [4] or the Stabilized Supralinear Network (SSN) [5] are more biologically plausible than DNNs, but these models confront considerable training challenges due to the lack of stability guarantees for high-dimensional problems. Training recurrent neural networks (RNNs), by comparison, is more straightforward thanks to ad hoc regularization techniques like layer nor-malization, batch normalization, and gradient clipping/scaling, which help stabilize training without imposing strict stability constraints. Conversely, neurodynamical models require enforcing hard stability constraints while maintaining biological plausibility. In lower dimensions, it is relatively straightforward to derive constraints on model parameters that ensure a dynamically stable system [6, 7] . However, for high-dimensional systems, this becomes significantly more challenging, as integrating these hard constraints into the optimization problem is more complex [8, 9] . Stability is generally advantageous in DNNs, as it is linked to improved generalization, mitigation of exploding gradient problems, increased robustness to input noise, and simplified training techniques [10] .The divisive normalization (DN) model was developed to explain the responses of neurons in the primary visual cortex (V1) [11] [12] [13] [14] , and has since been applied to diverse cognitive processes and neural systems [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] . Therefore, DN has been proposed as a canonical neural computation [25] that is linked to many well-documented physiological [26, 27] and psychophysical [28, 29] phenomena. DN models various neural processes: adaptation [30, 31] , attention [32] , automatic gain control [33] , decorrelation, and statistical whitening [34] . The defining characteristic of DN is that each neuron's response is divided by a weighted sum of the activity of a pool of neurons (Eq. 2, below) like when normalizing the length of a vector. Due to its wide applicability and ability to explain a variety of neurophysiological phenomena, we argue that this characteristic should be central to any neurodynamical model. Both the Wilson-Cowan and SSN models have been shown to approximate DN responses [5, 35] , but only approximately in certain parameter regimes.Normalization techniques have been extensively adopted for training DNNs, demonstrating their ability to stabilize, accelerate training, and enhance generalization [36] [37] [38] . Divisive normalization can be viewed as a comprehensive normalization strategy, with batch and layer normalization being specific instances [39] . Models implementing DN have shown superior performance compared to common normalization methods (Batch [36] , Layer [37] , Group [40] ) in tasks such as image recognition with convolutional neural networks (CNNs) [41] and language modeling with RNNs [39, 42] . Despite the foundational role of these techniques in deep learning algorithms, their implementation is ad hoc, limiting their conceptual relevance. They serve as practical solutions addressing the limitations of current machine learning frameworks rather than offering principled insights derived from understanding cortical circuits.It has been proposed that DN is achieved via a recurrent circuit [11, 13, [43] [44] [45] [46] [47] . Oscillatory recurrent gated neural integrator circuits (ORGaNICs) are rate-based recurrent neural circuit models that implement DN dynamically via recurrent amplification [47, 48] . Since ORGaNICs' response follows the DN equation at steady-state, its steady-state response captures the full range of aforementioned neural phenomena explained by DN . ORGaNICs have further been shown to simulate key time-dependent neurophysiological and cognitive/perceptual phenomena under realistic biophysical constraints [47, 48] . Additional phenomena not explained by DN [49] can in principle be integrated into the model. In this paper, however, we focus on the effects of DN on the dynamical stability of ORGaNICs. Despite some empirical evidence that ORGaNICs are highly robust, the question of whether the model is stable for arbitrary parameter choices, and thus whether it can be robustly trained on ML tasks by backpropagation-through-time (BPTT), remains open.Here, we establish the unconditional stability -applicable across all parameters and inputs -of a multidimensional two-neuron-types ORGaNICs model when the recurrent weight matrix is the identity. We prove this result, detailed in Section 4, by the indirect method of Lyapunov: we perform linear stability analysis around the model's analytically-known normalization fixed point and reduce the stability problem to that of a high-dimensional mechanical system, whose stability is defined in terms of a tractable quadratic eigenvalue problem. We then address the stability of the model with an arbitrary recurrent weight matrix in Section 5. While the indirect method of Lyapunov becomes intractable for such a system, we provide proof of unconditional stability for a two-dimensional circuit with an arbitrary recurrent weight and offer empirical evidence supporting the claim of stability for high-dimensional systems.ORGaNICs can be viewed as biophysically plausible extensions of Long Short Term Memory units (LSTMs) [3] and Gated Recurrent Units (GRUs) [50] , RNN architectures that have been widely used in ML applications [3, [51] [52] [53] [54] . The main differences are that ORGaNICs operate in continuous time and have built-in dynamic normalization (via recurrent gain modulation) and built-in attention (via input gain modulation). Thus, we expect that ORGaNICs should be able to solve relatively sophisticated tasks [47] . Here, we demonstrate (Section 6) that by virtue of their intrinsic stability, ORGaNICs can be trained on sequence modeling tasks by BPTT, in the same manner as traditional RNNs (unlike SSN that instead requires costly specialized training strategies [55] ), despite implementing power-law activations [5] . Moreover, we show that ORGaNICs trained by naive BPTT (i.e., without gradient clipping/scaling or other ad hoc strategies) achieve performance comparable to LSTMs on the tasks that we consider, despite no systematic hyperparameter tuning. 
related work: 
 Trainable biologically plausible neurodynamical models: There have been several attempts to develop neurodynamical models that mimic the function of biological circuits and that can be trained on cognitive tasks. Song et al. [56] incorporated Dale's law into the vanilla RNN architecture, which was successfully trained across a variety of cognitive tasks. Building on this, Soo et al. [57] developed a technique for such RNNs to learn long-term dependencies by using skip connections through time. ORGaNICs is a model that is already built on biological principles and can learn long-term dependencies intrinsically by tuning the (intrinsic or effective) time constants, therefore it does not require the method used in [57] . Soo et al. [55] introduced a novel training methodology (dynamics-neural growth) for SSNs and demonstrated its utility for tasks involving static (timeindependent) stimuli. However, this training approach is costly and difficult to scale (because SSNs, unlike ORGaNICs, are not unconditionally stable), and its applicability on tasks with dynamically changing inputs remains unclear.Dynamical systems view of RNNs: The stability of continuous-time RNNs has been extensively studied and discussed in a comprehensive review by Zhang et al. [58] . Recent advancements have focused on designing architectures that address the issues of vanishing and exploding gradients, thereby enhancing trainability and performance. A central idea in these designs is to achieve better trainability and generalization by ensuring the dynamical stability of the network. Moreover to avoid the problem of vanishing gradients the key idea is to constrain the real part of the eigenvalues of the linearized dynamical system to be close to zero, which facilitates the propagation and retention of information over long durations of time. Chang et al. [59] and Erichson et al. [60] achieve this by imposing an antisymmetric constraint on the recurrent weight matrix. Meanwhile, Rusch et al. [61, 62] propose an architecture based on coupled damped harmonic oscillators, resulting in a second-order system of ordinary differential equations that behaves similarly to how ORGaNICs behave in the vicinity of the normalization fixed point, as we show in Section 4. Despite their impressive performance on various sequential data benchmarks, these models lack biological plausibility due to their use of saturating nonlinearities (instead of normalization) and unrealistic weight parameterizations. 
model description: 
 In its simplest form, the two-neuron-types ORGaNICs model [47, 48] with n neurons of each type can be written as,EQUATIONwhere y ∈ R n and a ∈ R n are the membrane potentials (relative to an arbitrary threshold potential that we take to be 0) of the excitatory (y) and inhibitory (a) neurons, evolving according to the dynamical equations defined above with ẏ and ȧ denoting the time derivatives. The notation ⊙ denotes element-wise multiplication of vectors, and squaring, rectification, square-root, and division are also performed element-wise. 1 is an n-dimensional vector with all entries equal to 1. z ∈ R n is the input drive to the circuit and is a weighted sum of the input, x ∈ R m , i.e., z = W zx x. The firing rates, y ± = ⌊±y⌋ 2 and a + = ⌊a⌋ are rectified (⌊.⌋) power functions of the underlying membrane potentials. For the derivation of a general model with arbitrary power-law exponents, including the Eq. 1, see Appendix A. Note that the term y + -y -serves the purpose of defining a mechanism for reconstructing the membrane potential (which can be negative, depending on the sign of the input) from the firing rates y ± that are strictly nonnegative. y + and y -are the firing rates of neurons with complementary receptive fields such that they encode inputs with positive and negative signs, respectively. Note that only one of these neurons fires at a given time. In ORGaNICs, these neurons have a single dynamical equation for their membrane potentials, where the sign of y indicates which neuron is active. Neurons with such complementary (anti-phase) receptive fields are found adjacent to each other in the visual cortex [63] , and we hypothesize that such complementary neurons are ubiquitous throughout the neocortex. b ∈ R + * n and b 0 ∈ R + * n are the input gains for the external inputs z and σ fed to neurons y and a, respectively. R + * is the set of positive real numbers, {x ∈ R | x > 0}. σ ∈ R + * n determines the semisaturation of the responses of neurons y by contributing to the depolarization of neurons a. τ y ∈ R + * n and τ a ∈ R + * n represent the time constants of y and a neurons.In addition to receiving external inputs, both y and a neurons receive recurrent inputs, represented by the last term in both of the equations. W r ∈ R n×n is the recurrent weight matrix that captures lateral connections between the y neurons. This recurrent input is gated by the a neurons, via the term (1 -a + ). Similarly, the nonnegative normalization weight matrix, W ∈ R * n×n , encapsulates the recurrent inputs received by the a neurons. The differential equations are designed in such a way that when W r = I and b = b 0 = b 0 1 (i.e., with all elements equal to a constant b 0 ), the principal neurons follow the normalization equation exactly (and approximately when W r ̸ = I) at steady-state,EQUATION)⌊z⌋ 2 and ⌊-z⌋ 2 represent the contribution of neurons with complementary receptive fields to the normalization pool, and ⌊z⌋ 2 + ⌊-z⌋ 2 = z 2 is the contrast energy of the input. Note that the recurrent gain, (1 -a + ), is a particular nonlinear function of the output responses/activation designed to achieve DN, while the input gain, b + , is an input gate that can implement an attention mechanism. 
stability analysis of high-dimensional model with identity recurrent weights: 
 We consider the stability of the general high-dimensional ORGaNICs (Eq. 1) when the recurrent weight matrix is identity, W r = I. We first simplify the dynamical system by noting that y +y -= y and y + + y -= y 2 yielding the following equations,EQUATIONFor identity recurrent weights, we have a unique fixed point, given by,EQUATIONSince the normalization weights in the matrix W are nonnegative, at steady-state we have a s > 0, so that ⌊a s ⌋ = √ a s , and the corresponding firing rates at steady-state are,EQUATIONNote that we recover the normalization equation, Eq. 2, if b = b 0 = b 0 1. Since the fixed points of y and a neurons are known analytically, to prove that this fixed point is locally asymptotically stable (i.e., the responses converge asymptotically to the fixed point), we apply the indirect method of Lyapunov at this fixed point [64] . This method allows us to analyze the stability of the nonlinear system in the vicinity of the fixed point by studying the corresponding linearized system. The Jacobian matrix J ∈ R 2n×2n about (y s , a s ), defining the linearized system, is given by,EQUATIONwhere D(x) is a diagonal matrix of appropriate size with the elements of the vector x on the diagonal. A necessary and sufficient condition for local stability is that the real parts of all eigenvalues of this matrix are negative. We thus proceed by computing the characteristic polynomial for the Jacobian, p J (λ) ≡ det(J -λI). The roots of this polynomial, found by setting p J (λ) = 0, are the eigenvalues of the system. Consider the block matrix,EQUATIONNotice that A 11 and A 12 are diagonal and therefore they commute, i.e., A 11 A 12 = A 12 A 11 , so we have that det(J -λI) = det(A 22 A 11 -A 21 A 12 ) which is a property of the determinant of block matrices with commuting blocks [65] . Therefore, the characteristic polynomial of the linearized system after expansion of the terms and simplification is given by,EQUATIONFinding the roots of this polynomial is thus a quadratic eigenvalue problem of the form L(λ) ≡ det(λ 2 I + λB + K) = 0, which has been studied extensively [66] [67] [68] [69] . L(λ) can be interpreted as the characteristic polynomial associated with a system of linear second-order differential equations with constant coefficients of the form Iẍ + B ẋ + Kx = 0. Therefore, proving the stability of our system (i.e., Re(λ) < 0 for {λ : L(λ) = 0}), is equivalent to proving the asymptotic stability ofIẍ + B ẋ + Kx = 0.Tisseur et al. [67] and Kirillov et al. [69] list a set of constraints on the damping matrix, B, and stiffness matrix, K, that yield a stable system, but they are not directly applicable to our system.In the context of a high-dimensional mechanical system, our system falls under the category of gyroscopically stabilized systems with indefinite damping. Few results are known about the conditions leading to the stability of such systems. By constructing a Lyapunov function, we prove (Appendix B) the following stability theorem that is directly applicable to our system, following an approach similar to Kliem et al. [70] . Theorem 4.1. For a system of linear differential equations with constant coefficients of the form,EQUATIONwhere B ∈ R n×n and K ∈ R n×n is a positive diagonal matrix (hence K ≻ 0), the dynamical system is globally asymptotically stable if B is Lyapunov diagonally stable.Since the stiffness matrix,EQUATION)is a positive diagonal matrix, a sufficient condition for stability of the system is that the damping matrix, B, given by, We now show that, indeed, B has a convergent regular splitting for all combinations of the circuit parameters and for all inputs. We have already shown that B is a Z-matrix, therefore, the first condition of the theorem is satisfied. Next, we consider the following splitting B = M -N with M = B 1 + B 2 and N = B 3 . Since B 1 and B 2 are positive diagonal matrices, M -1 is nonnegative, while N is also nonnegative because B 3 has all positive entries. Therefore, the only condition left to satisfy is that the spectral radius of M -1 N is smaller than 1, or that the matrix is convergent.EQUATIONThe matrix S = M -1 N = (B 1 + B 2 ) -1 B 3 can be written as,EQUATIONWe prove the following theorem (Appendix D) which directly applies to S, Theorem 4.3. A matrix A of the form A = D(t) W D (u/ (v + Wu)) is convergent (i.e., its spectral radius is less than 1), if W ∈ R n×n and t, u, v ∈ R n satisfy 0 < t i < 1, u i ≥ 0, v i > 0 and w ij ≥ 0 for all i, j.Defining t → 1/(1 + (τ a /τ y ) ⊙ b 2 0 ⊙ σ 2 + W (b 2 ⊙ z 2 )), u → b 2 ⊙ z 2 and v → b 2 0 ⊙ σ 2, it can be seen that they satisfy the constraints of the theorem, and thus S is convergent. This implies that B has a convergent regular splitting and, as a result, the linearized dynamical system is unconditionally globally asymptotically stable for all the values of parameters and inputs. Further, the global asymptotic stability of linearization implies the local asymptotic stability of the normalization fixed point for ORGaNICs. This result holds even when the neurons have different time constants, regardless of their type, as no assumptions were made about the time constants. This finding is significant for machine learning, particularly for designing architectures based on ORGaNICs. It allows neurons/units to integrate information at varying time scales while maintaining a stable circuit that performs normalization dynamically. Moreover, analytical expressions for eigenvalues can be obtained in the following case, Theorem 4.4. Let W r = I, the normalization matrix be given by W = αE, where E is the all-ones matrix, and the parameters are scalars, i.e., τ y = τ y 1, τ a = τ a 1, b 0 = b 0 1, and σ = σ1. Under these conditions, the eigenvalues of the system admit closed form solutions (detailed in Appendix C).This result is particularly useful for neuroscience as it elucidates the connection between ORGaNICs parameters and the strength and frequency of oscillatory activity. Since we followed a direct Lyapunov approach to prove Theorem 4.1 as shown in Appendix B, we can derive an energy (viz., Lyapunov function) for ORGaNICs as shown in Appendix H. Theorem 4.5. When W r = I, the energy (Lyapunov function) minimized by ORGaNICs in the vicinity of the normalization fixed point, is given by,EQUATION)Where t i are the diagonal entries of T and y si (a si ) are the steady-state values of neurons y i (a i ).Specifically, for a two-dimensional model (one y neuron and one a neuron) this expression simplifies to reveal that ORGaNICs behave like a damped harmonic oscillator with energy,EQUATIONThis result demonstrates that ORGaNICs minimize the residual of the instantaneously reconstructed gated input drive ( √ ay -bz), while also ensuring that the principal neuron's response, y, achieves DN. The balance between these objectives is governed by the parameters and the external input strength. With fixed parameters, weaker inputs, z, cause the model to prioritize input matching over normalization, whereas stronger inputs increasingly engage the normalization objective. 
stability analysis for arbitrary recurrent weights: 
 Now, we relax the constraint that the recurrent weight matrix must be identity, allowing W r ̸ = I, and see how the stability result changes. This leads to the following set of equations, The linear stability analysis becomes intractable for a general W r because we no longer have a closed-form analytical expression for the steady states of y and a. Additionally, the characteristic polynomial cannot be expressed in a way similar to Eq.8. Nevertheless, for a two-dimensional system,EQUATIONEQUATIONwe can prove the following, with a detailed analysis provided in Appendix E. Theorem 5.1. Given that the recurrence is contracting, i.e., 0 < w r ≤ 1, when z > 0 (z < 0) there exists a unique fixed point with y s > 0 (y s < 0) and a s > 0, and it is asymptotically stable. Theorem 5.2. Given that the recurrence is expansive, i.e., w r > 1, there are either 1 or 3 fixed points of which at least one is asymptotically stable. When z > 0 (z < 0) there exists exactly 1 fixed point with y s > 0 (y s < 0) and a s > 0, and it is asymptotically stable. If b 0 σ > 1 -1/w r , there are no additional fixed points. If b 0 σ < 1 -1/w r , there exist either 0 or 2 additional fixed points with y s < 0 (y s > 0) and a s > 0 whose stability cannot be guaranteed.We plot the phase portraits for these different cases in Fig. 1 . The key takeaway is that there is always a fixed point (y s , a s ) with a s > 0 and y s having the same sign as z. This fixed point is asymptotically stable regardless of the value of w r . Based on these results and the proven stability of arbitrary dimensional ORGaNICs when W r = I (as shown in Section 4), we conjecture that Conjecture 5.3. Consider high-dimensional ORGaNICs with an arbitrary recurrent weight matrix W r and no constraints on the remaining parameters. If the norm of the input drive satisfies ||z|| ≤ 1, and the maximum singular value of W r is constrained to be 1, then the system possesses at least one asymptotically stable fixed point.This conjecture is supported by empirical evidence showing consistent stability, as ORGaNICs initialized with random parameters and inputs under these constraints have exhibited stability in 100%Algorithm 1 Iterative scheme for the fixed point when the maximum singular value of W r is 11: Input: ORGaNICs parameters, input (z), tolerance (ϵ), maximum iterations (N ) 2: Output: Approximation to the fixed point (y s , a s )3: a ← σ 2 ⊙ b 2 0 + W (W r (b ⊙ z)) 2 // initial approximation for a neurons 4: y ← (W r (b ⊙ z)) / √ a // initial approximation for y neurons 5: k ← 0 6: while ||y -b ⊙ z -1 - √ a ⊙ (W r y) || > ϵ and k < N do 7: y ← I -W r + D √ a W r -1 (b ⊙ z) // y neurons update 8: a ← b 2 0 ⊙ σ 2 + W y 2 * a // a neurons update 9:k ← k + 1 10: end while 11: return (y, a) of trials, see Fig. 4 . We further speculate that ORGaNICs may be typically stable beyond this regime as we find that 100% of trials yield a stable circuit when the constraint on the maximum singular value of W r is increased to 2, but it becomes unstable when it is increased to 3. 
experiments: 
 We provide further empirical evidence in support of Conjecture 5.3 that ORGaNICs is asymptotically stable by showing that stability is preserved when training ORGaNICs using naïve BPTT on two different tasks: 1) static classification of MNIST, 2) sequential classification of pixel-by-pixel MNIST. Because these ML tasks have no relevance for neurobiological or cognitive processes, we relax one aspect of the biological plausibility of ORGaNICs, specifically, allowing arbitrary (learned) nonnegative values for the intrinsic time constants. 1 We first show that we can train ORGaNICs on the MNIST handwritten digit dataset [72] presented to the circuit as a static input. This setting corresponds to evolving the responses of the neurons dynamically until they reach a fixed point solution and using the steady-state firing rates of the principal neurons to predict the labels, akin to deep equilibrium models [73] . While the fixed point of the circuit is known when W r = I (given by Eq. 89), we allow W r to be learnable and parameterized it to have a maximum singular value of 1. This constraint allows us to find the fixed point responses of all the neurons without simulation, using a fixed point iteration scheme (Algorithm 1) that converges with great accuracy in a few (less than 5) steps, see Fig. 4 & 5. We provide an intuition for why this algorithm works with empirical evidence of fast convergence in Appendix G.We trained ORGaNICs on this task (details provided in Appendix I.1) and compared its performance to SSN [5] trained by dynamics-neutral growth [55] . We found that ORGaNICs perform better than SSN with the same model size, and on par with an MLP (Table 1 ). We analyzed the eigenvalues of the Jacobian matrix of the trained model and consistently found the largest real part to be negative (Fig. 5 ), indicating stability. Moreover, we found that stability was maintained during training (Fig. 6 ).We trained unconstrained ORGaNICs by naïve BPTT on a classification task of sequential MNIST (sMNIST), proposed by Le et al. [74] . This is a challenging task because it involves long-term dependencies and requires the architecture to maintain and integrate information over long timescales. Briefly, the task involves the presentation of pixels of MNIST images sequentially (one pixel for each timestep) in scanline order, and at the end of the input the model has to predict the digit that was presented. There is a more complicated version of this task, permuted sequential MNIST, in which the pixels of all images are permuted in some random order before being presented sequentially. We train ORGaNICs with different hidden layer sizes (number of y neurons) on these two tasks by discretizing the rectified ORGaNICs with arbitrary recurrence, Eq. 87, which has all the properties that we have derived for the main model. Since an unstable fixed point is undesirable in such a task, as it may lead to diverging trajectories, we prefer the rectified model (Appendix F) over the main model. We proved that the 2D rectified ORGaNICs (Eq. 102) does not exhibit an unstable fixed point for positive inputs, as it can also be seen in Fig 1 . The hidden states of the neurons are initialized with a uniform random distribution (for more details, see Appendix I.2). Additionally, we make the input gains b and b 0 dynamical with their ODEs given by,EQUATIONWe achieved slightly better performance than LSTMs on sMNIST with a smaller model size and comparable performance on permuted sMNIST, without hyperparameter optimization and without gradient clipping/scaling (Table 2 ). We found that the trajectories of y are bounded when it is trained on the sequential task (Fig. 7 ), indicating stability. We also show that the training of ORGaNICs is stable and does not require gradient clipping when the intrinsic time constants of the neurons are fixed (Table 2 ). 
discussion: 
 Summary: While extensive research has been aimed at identifying highly expressive RNN architectures that can model complex data, there has been little advancement in developing robust, biologically plausible recurrent neural circuits that are easy to train and perform comparably to their artificial counterparts. Regularization techniques such as batch, group, and layer normalization have been developed and are implemented as ad hoc add-ons making them biologically implausible. In this work, we bridge these gaps by leveraging the recently proposed ORGaNICs model which implements divisive normalization (DN) dynamically in a recurrent circuit. We establish the unconditional stability of an arbitrary dimensional ORGaNICs circuit with an identity recurrent weight matrix (W r ), with all of the other parameters and inputs unconstrained, and provide empirical evidence of stability for ORGaNICs with arbitrary W r . Since ORGaNICs remain stable for all parameter values and inputs, we do not need to resort to techniques that are restrictive in parameter space, or that require designing unrealistic structures for weight matrices. ORGaNICs' intrinsic stability mitigates the issues of exploding and oscillating gradients, enabling the use of "vanilla" BPTT without the need for gradient clipping, which is instead required when training LSTMs. Moreover, ORGaNICs effectively address the vanishing gradient problem often encountered when training RNNs. This is achieved by processing information across various timescales, resulting in a blend of lossy and non-lossy neurons, while preserving stability. The model's effectiveness in overcoming vanishing gradients is further evidenced by its competitive performance against architectures specifically designed to address this issue, such as LSTMs.Dynamic normalization: Normalization techniques, such as batch and layer normalization, are fundamental in modern ML architectures significantly enhancing the training and performance of CNNs. However, a principled approach to incorporating normalization into RNNs has remained elusive. While layer normalization is commonly applied to RNNs to stabilize training, it does not influence the underlying circuit dynamics since it is applied a-posteriori to the output activations, leaving the stability of RNNs unaffected. Furthermore, DN has been shown to generalize batch and layer normalization [39] , leading to improved performance [39, 41, 42] . ORGaNICs, unlike RNNs with layer normalization, implement DN dynamically within the circuit, marking the first instance of this concept being applied and analyzed in ML. Our work demonstrates that embedding DN within a circuit naturally leads to stability, which is greatly advantageous for trainability. This stability, a consequence of dynamic DN, sets ORGaNICs apart from other RNNs by providing both output normalization and model robustness. As a result, ORGaNICs can be trained using BPTT, achieving performance on par with LSTMs. The key insight is that the dynamic application of DN not only enhances training efficiency but also improves model robustness. This illustrates how the incorporation of neurobiological principles can drive advances in ML.Interpretability: In the proof of stability, we establish a direct connection between ORGaNICs and systems of coupled damped harmonic oscillators, which have long been studied in mechanics and control theory. This analogy not only enables us to derive an interpretable energy function for ORGaNICs (Eq. 13), providing a normative principle of what the circuit aims to accomplish, but also sheds light on the link between normalization and dynamical stability of neural circuits. For a relevant ML task, having an analytical expression for the energy function allows us to quantify the relative contributions of the individual neurons in the trained model, offering more interpretability than other RNN architectures. For instance, Eq. 13 shows that the ratio of time constants (τ y /τ a ) for E-I neuron pairs determines how much weight a neuron assigns to divisive normalization relative to aligning its responses with the input drive z. This insight provides a clear functional role for each neuron in the trained model. Moreover, since ORGaNICs are biologically plausible, we can understand how the various components of the dynamical system might be computed within a neural circuit [48] , bridging the gap between theoretical models and biological implementation, and offering a means to generate and test hypotheses about neural computation in real biological systems (which we will be reporting elsewhere).Limitations: Although the stability property pertains to a continuous-time system of nonlinear differential equations, typical implementations for tasks with sequential data involve an Euler discretization of these equations for training purposes. This might lead to a stiff dynamical system, potentially causing numerical instabilities and explosive dynamics, highlighting the importance of carefully parameterizing time constants and choosing a small enough time step to maintain stable dynamics. The proof of unconditional stability is only tractable for the two-dimensional circuit and the high-dimensional circuit with W r = I. Therefore, we can only conjecture the stability of ORGaNICs for arbitrary W r , based on these two limiting cases and on empirical evidence. In the current form, the weight matrices of the input gain modulators, W by , W ba , W b0y , and W b0a , are each n × n. As a result, the number of parameters grows more rapidly with the hidden state size compared to other RNNs. To mitigate this, we plan to explore using compact and/or convolutional weights to prevent a significant increase in the number of parameters as the hidden state size expands.Attention mechanisms in ORGaNICs: ORGaNICs have a built-in mechanism for attention: modulating the input gain b (e.g., Eq. 17), coupled with DN. This attention mechanism aligns with experimental data on both increases in the gain of neural responses and improvements in behavioral performance [19, 20, 32, [76] [77] [78] [79] [80] [81] [82] [83] [84] [85] . Moreover, this mechanism performs a computation that is analogous to that of an attention head in ML systems (including transformers [2] ) as both operate by changing the gain over time. In ORGaNICs, DN replaces the softmax operation typically used in an attention head.Future work: This study has explored only a single layer of ORGaNICs for the sequential tasks. Future work will examine how stacked layers with feedback connections, similar to those in the cortex, perform on benchmarks for sequential modeling and also on cognitive tasks with long-term dependencies. We have thus far shown that ORGaNICs can address the problem of long-term dependencies by learning intrinsic time constants. Future investigations will assess the performance of ORGaNICs for tasks with long-term dependencies by learning to modulate the responses of the a and b neurons to control the effective time constant of the recurrent circuit (without changing the intrinsic time constants) [47] , i.e., implementing a working memory circuit capable of learning to maintain and manipulate information across various timescales.Here, we derive a generalized 2-neuron types (excitatory and inhibitory) ORGaNICs model for a high-dimensional input. The system presented in Eq. 1 is a special case of this generalized model where p = 2 and a + = ⌊a⌋. Assuming W is the normalization weight matrix, and z is the input drive, we can write the normalization equations for principal neurons with complementary receptive fields as,EQUATIONNote that typically the exponent of the input p ∼ 2 for cortical neurons. ⌊z⌋ p and ⌊-z⌋ p represent the contribution of neurons with complementary receptive fields to the normalization pool. Mathematically, we have, ⌊z⌋ p + ⌊-z⌋ p = |z| p .Here, we derive, for a general p, the dynamical equations that have the fixed point defined by the normalization equation above. First, it is important to distinguish between the membrane potentials and the firing rates of neurons. The coarse (low-pass filtered) membrane potential of a given type of neuron is denoted by the vector, y, a, with the corresponding firing rates of y + (y -), a + . The instantaneous firing rates of the neurons are obtained from the corresponding coarse membrane potentials by applying rectification, denoted by ⌊.⌋, and a power law (sub/supra-linear) activation for different types of neurons [86] [87] [88] . Therefore, for a set of membrane potentials x the instantaneous firing rates are x + = ⌊x⌋ α . Specifically for principal neurons, we have, y + = ⌊y⌋ p and y -= ⌊-y⌋ p . Combining the firing of principal neurons with the complementary receptive fields, y + and y -, Eq. 18 can be alternatively written as,EQUATIONNow for the principal neuron y j receiving an input drive z j , we can rewrite the normalization equation for each neuron as,EQUATIONIn the ORGaNICs paradigm [47] , the steady-state activity of the principal neurons (single equation for complementary receptive fields) is a weighted sum of input drive and recurrent drive,τ yj dy j dt = -y j + b j z j Weighted input drive + 1 -a + j k w r jk (y + k ) 1/p -(y - k ) 1/pWeighted recurrent drive (21) Here w r are the weights of the recurrent weight matrix W r encoding the recurrent/lateral connections between the principal neurons y; y + k is the firing rate of the principal neuron k, given by y + k = ⌊y k ⌋ p , and y - k = ⌊-y k ⌋ p is the firing rate of the complementary principal neuron k. b j is the gain to the input drive z j which simulates attention (realized via gain modulation). (1 -a + j ) is the gain to the recurrent drive which controls the recurrent amplification. Now, we find the dynamics of the inhibitory neurons a j with firing rates a + j , that yield stable dynamics with the fixed point given by Eq. 18. First, we assume that the recurrent weight matrix W r = I. Also, note that the following identity holds:(y + k ) 1/p -(y - k ) 1/p = ⌊y k ⌋ -⌊-y k ⌋ = y k .Therefore, Eq.21 can be simplified to,EQUATIONAt steady-state, the fixed-points (y s j , a s j ), satisfy the following relationship,EQUATIONTaking modulus and raising both sides to power p, we get,EQUATIONor in vector form,EQUATIONFrom Eq. 19, we know thatEQUATIONSince we can write the element-wise product between two vectors x 1 and x 2 ,x 1 ⊙ x 2 = x 2 ⊙ x 1 = D(x 1 )x 2 .Where D(x 1 ) is a diagonal matrix with elements of the vector x 1 on the diagonal. Therefore, using this fact we can rewrite the equation above as,EQUATIONSubstituting the expression for |z| p into Eq. 25, we get,EQUATIONWe simplify this equation as,EQUATIONNow we assume that all of the entries of b are equal to a constant b 0 , therefore, we have,EQUATIONElement-wise multiplying both sides by a s / (a + s ) p on the left, we get,EQUATIONThis equation is true at the steady-state, but it is also in a form similar to that of y, such that we have a weighted input drive and a weighted recurrent drive, therefore, we propose the following dynamical equation for a,EQUATIONThis equation naturally follows Eq. 31 at steady-state and thus also follows the normalization equation Eq. 19. Note that the equation above is true for any choice of a + , but cortical neurons have been experimentally observed to have an exponent close to p = 2, i.e., y + = ⌊y⌋ 2 and y -= ⌊-y⌋ 2 . Therefore, we get a particularly simple form of the equations when p = 2 and a + = ⌊a⌋. In this case, the equation for a is given by,EQUATIONWe reintroduce the recurrent weight matrix W r and to simulate the effect of attention by gain modulation, we define different gains for y and a neurons. Additionally, replacing y 2 → y + + y - and y → y + -y -yields the dynamical system that we analyze,EQUATIONFinally, in the most general form, i.e., any choice of p and nonlinearity for a + , the dynamical system of ORGaNICs is given by,EQUATIONTheorem B.1. For a system of linear differential equations with constant coefficients of the form,EQUATION)where B ∈ R n×n and K ∈ R n×n is a positive diagonal matrix (hence K ≻ 0), the dynamical system is globally asymptotically stable if B is Lyapunov diagonally stable.Proof. The stability of the system is defined by solving the following associated quadratic eigenvalue problem,EQUATIONThe spectrum of L(λ), i.e., {λ ∈ C : det(L(λ)) = 0} are also known as the eigenvalues of the system. The system defined by Eq. 36 is globally asymptotically stable if all of the eigenvalues have negative real parts. We take the direct Lyapunov approach to prove the stability of the linear system. We write Eq. 36 in the matrix form as ż = Jz, where,EQUATIONWe will first prove the Lyapunov stability (Re(λ J ) ≤ 0) of this system to find the appropriate block diagonal matrices and then we will prove global asymptotic stability (Re(λ J ) < 0). To prove Lyapunov stability (Appendix B.2), we propose a Lyapunov function V (z) = z ⊤ Pz, where P is a block positive definite matrix defined as follows,EQUATIONwhere T ∈ R n×n is a positive diagonal matrix such that TB + B ⊤ T ≻ 0 (notation for positive definite matrix) and A ∈ R n×n a flexible symmetric positive definite matrix that we will find using the second Lyapunov criteria. Note that such a matrix T exits since B is defined to be Lyapunov diagonally stable. It can be easily seen that P ≻ 0 using the first criteria in Section B.1 since A ≻ 0 and T ≻ 0. Additionally, since T ≻ 0, it is invertible. Now, for Lyapunov stability, we need V (z) = z ⊤ PJ + J ⊤ P z ≤ 0. Therefore, we find A such that Q = -PJ + J ⊤ P is positive semi-definite.Q = -PJ + J ⊤ P = - A 0 0 T 0 I -K -B - 0 -K I -B ⊤ A 0 0 T = 0 KT -A TK -A TB + B ⊤ T (40)We want to define A ≻ 0, such that Q ⪰ 0. Using the second criteria from Section B.1, we need TB + B ⊤ T ≻ 0 and -(KT -A)(TB + B ⊤ T) -1 (TK -A) ⪰ 0. The first condition is satisfied by the definition of T. For the second condition to be satisfied, an obvious candidate for A is TK. Note that both T and K are positive definite and diagonal, therefore they commute (KT = TK) and A is symmetric and positive definite. When A = TK, the LHS of the second condition becomes 0 ⪰ 0. Therefore, the system is Lyapunov stable. Now, we prove the global asymptotic stability of the system by again using the direct Lyapunov approach. We propose the Lyapunov function of the same form as before, i.e., V (z) = z ⊤ Pz, where P is a positive definite matrix. But for global asymptotic stability, we need a more stringent condition on the Lyapunov function, V (z) = z ⊤ PJ + J ⊤ P z < 0. Drawing inspiration from the previous exercise, we consider the following form of the matrix P, P = TK ϵI ϵI THere, ϵ > 0 is a scalar whose magnitude is to be determined based on the Lyapunov criteria for asymptotic stability. First, we need P ≻ 0. Applying the first criteria from Section B.1, we want ϵ to satisfy, TK -ϵ 2 T -1 ≻ 0. Second, for V (z) < 0, we want Q = -PJ + J ⊤ P to be positive definite. Q is given by,EQUATIONAgain, we apply the first criteria from Section B.1.Q ≻ 0 if and only if TB + B ⊤ T -2ϵI ≻ 0 and 2ϵK -ϵ 2 B TB + B ⊤ T -2ϵI -1 B ⊤ ≻ 0.To simplify notation we replace TB + B ⊤ T with a positive definite matrix M. Therefore, we have to prove that there exists an ϵ > 0 which satisfies the following conditions,• TK -ϵ 2 T -1 ≻ 0 • M -2ϵI ≻ 0 • 2K -ϵB (M -2ϵI) -1 B ⊤ ≻ 0Assuming t i and k i to be the diagonal values of the positive diagonal matrices T and K, we get the following two conditions, ϵ < min (t i √ k i ) and ϵ < α/2, where α is the smallest eigenvalue of M, or α = min (λ M ). Now, for the third condition, we will use a number of facts about positive definite matrices which are all listed in [89] . We first consider the following matrix inequality, M ⪰ αI. This notation is equivalent to saying that M -αI is positive semi-definite or M -αI ⪰ 0. Therefore, we have,EQUATIONAssuming, ϵ is small enough such that the matrices on LHS and RHS are positive definite, we have,EQUATIONSince B is nonsingular, it is full rank, therefore, we have,EQUATIONMultiplying both sides by ϵ, we get,EQUATIONNotice that BB ⊤ is a positive definite matrix. Let β be the maximum eigenvalue of BB ⊤ , or β = max (λ BB ⊤ ). Therefore, we can add an upper-bound matrix to the inequality as follows,EQUATIONNow, we find ϵ such that,EQUATIONThe range of values for which the above inequality is true is,EQUATIONIf ϵ satisfies the inequality above, we have,EQUATIONTherefore, we have, 2K ≻ ϵB (M -2ϵI) -1 B ⊤ and the third condition required for ϵ is satisfied.This implies that there exists a range of ϵ, given by,EQUATIONfor whichEQUATIONis a valid Lyapunov function for asymptotic stability. Therefore, the dynamical system is globally asymptotically stable.For a symmetric block matrix of the form,EQUATIONwith A = A ⊤ and C = C ⊤ . If C is invertible the following two properties hold, [90] ,• P ≻ 0 if and only if C ≻ 0 andA -BC -1 B ⊤ ≻ 0. • If C ≻ 0, then P ⪰ 0 if and only if A -BC -1 B ⊤ ⪰ 0.Consider a non-linear autonomous dynamical system defined as ẋ = f (x), with a point of equilibrium at x = 0. Where x ∈ D ⊆ R n is the system state vector and f (x) : D → R n is a continuous vector field on D (contains origin). The dynamical system is called Lyapunov stable if there exists a real scalar function V (x) : R n → R, also known as the Lyapunov function, such that it satisfies the following conditions,• V (x) = 0, if and only if x = 0.• V (x) > 0, if and only if x ̸ = 0.• V (x) ≤ 0, ∀ x ̸ = 0. Note that for asymptotic stability, we require the strict inequality V (x) < 0.For a linear dynamical system of the form ẋ = Jx, where J ∈ R n×n , with a point of equilibrium at x = 0. Consider a Lyapunov function V (x) of the form x ⊤ Px, such that P ≻ 0. By the definition of a positive definite matrix, it satisfies the first two conditions, namely, V (x) = x ⊤ Px = 0 when x = 0 and V (x) = x ⊤ Px > 0 when x ̸ = 0. For the third condition, consider V (x),EQUATIONFor stability, We need V (x) = x ⊤ PJ + J ⊤ P x ≤ 0. This is satisfied when the matrix PJ + J ⊤ P ⪯ 0. In summary, a linear dynamical system ẋ = Jx is Lyapunov stable if there exists a positive definite matrix P, such that PJ + J ⊤ P is negative semi-definite.Here we show that when all of the normalization weights in the system are equal, to value α, and the various parameters are scalars, i.e., τ y = τ y 1, τ a = τ a 1, b 0 = b 0 1 and σ = σ1, we can derive a closed-form analytical expression for all of the eigenvalues. Considering these assumptions, we can break the determinant in Eq. 8 into a diagonal and non-diagonal part as follows,EQUATIONConsider the non-diagonal part of the matrix in the determinant, (λ/τ a )WD y 2 s . Since W is a matrix with all entries equal to a positive constant, α, it is rank 1. Therefore, it can be written as the following outer product,EQUATIONTherefore, the non-diagonal part of the matrix can be written as,EQUATIONwhereu = [1, 1, ..., 1] ⊤ and v = b 2 ⊙ z 2 / σ 2 b 2 0 1 + W b 2 ⊙ z 2 .We use the matrix determinant lemma which states that,EQUATIONThe matrix A is given by,EQUATIONHere, ||x||, represents the Euclidean norm of x. Therefore, A = δI, where δ is a quadratic scalar polynomial in λ. Now using Eq. 58, we can writeEQUATIONSolving for the eigenvalues, det(J -λI) = 0, we get 2(n -1) repeated solutions by equation δ n-1 = 0, where each solution is found by solving δ = 0,EQUATIONThis gives us the following strictly negative eigenvalues,EQUATIONThe potentially complex eigenvalues are given by solving for zeroes of the first part of the factorized determinant,λ 2 + λ σ 2 b 2 0 τ a (σ 2 b 2 0 + α||b ⊙ z|| 2 ) + σ 2 b 2 0 + α||b ⊙ z|| 2 τ y + σ 2 b 2 0 + α||b ⊙ z|| 2 τ y τ a = 0 (63)The eigenvalues found by solving this equation have negative real parts (as expected) for all choices of parameters since the coefficient of λ and the constant term of the quadratic equation are positive for all choices of parameters and inputs. Therefore, this solution does admit complex roots (thereby oscillations) for certain choices of the parameters which can be found by solving this quadratic equation.Theorem D.1. A matrix A of the form,EQUATION)is convergent, i.e., its spectral radius is less than one, if W ∈ R n×n and t, u, v ∈ R n with the additional constraints 0 < t i < 1, u i ≥ 0, v i > 0 and w ij ≥ 0 for all i, j.Proof. Assuming the constraints mentioned in the theorem, we first notice that,EQUATION)is a nonnegative matrix. Further, because D(t) is a positive diagonal matrix with all entries less than 1, we notice that the following is true element-wise,EQUATION)We denote the spectral radius of A, max{|λ| : λ = σ(A)}, by ρ(A), where σ(A) is the spectrum of A and make use of the following inequality for the spectral radius of nonnegative matrices., Theorem D.2. (Theorem 8.1.18 from [91] ) Let X and Y be nonnegative matrices with spectral radius ρ(X) and ρ(Y), respectively. If X ≤ Y (x ij ≤ y ij , ∀ i, j), then ρ(X) ≤ ρ(Y).Therefore, we have,EQUATION)Now the elements of the matrix are given by,EQUATIONLet s ij be the i, j element of this matrix. Upon multiplication of the matrices above, we find that,EQUATION)Now we use the following theorem that provides an upper bound for the spectral radius based on the entries of the matrix.Theorem D.3. (Theorem 8.1.26 from [91] ) Let S be a nonnegative matrix. Then for any positive vector x ∈ R n with entries x j , we have,EQUATIONWe pick x j = v j k w jk u k , which is a positive vector, and apply the theorem to the matrix. Substituting s ij and x j , we get the following bound for the spectral radius,ρ W D u v + Wu ≤ max 1≤i≤n 1 v i + k w ik u k n j=1 w ij u j v j + k w jk u k v j + k w jk u k (71)Upon simplification, we get,EQUATION)The inequality is true, because, for all i, the denominator is larger than the numerator because of the extra positive term v i . Therefore,EQUATION)and A is a convergent matrix.Here we consider the dynamical stability of the 2D model containing one neuron each of y and a when the recurrent scalar, w r can take any positive value. The dynamical system to be analyzed is given by,EQUATIONwith a positive real constraint on the following set of parameters, τ y , τ a , b, b 0 , σ, w. We first notice the symmetry in the dynamical system about y = 0. Replacing, ŷ → -y, we get the mirrored dynamical system,EQUATION)These equations are the same as Eq. 74, up to a sign change in z. Therefore, we can derive analogous conditions for stability for z < 0 once the conditions for z > 0 are known.The steady-state of Eq. 74, (y s , a s ) satisfies the following equations,EQUATION)EQUATIONSince the RHS of the second equation is always positive, if a root exists, we have a s > 0. Therefore, we can remove the rectification around a s and look for positive solutions for a s . Upon rearranging the terms, we get,EQUATIONSubstituting y s from Eq. 78 in Eq. 79, we get the following quartic equation in m = √ a s ,EQUATIONFor a valid fixed point (y s , a s ), i.e., y s ∈ R and a s ∈ R + * , a necessary condition is that a s must be positive and real, which in turn implies m must be positive and real. Theorem E. Proof. =⇒ Given m > b 0 σ, we have m > 0 because b 0 > 0 and σ > 0. ⇐= Since m is a root of the quartic, it also satisfies Eq. 79. Therefore, we have m = b 0 σ/ 1 -wy 2 s . Since we are given m > 0, this implies that 0 < 1 -wy 2 s < 1, or m > b 0 σ. Therefore, m > 0 is equivalent to m > b 0 σ. This also implies that there exists no fixed point that satisfies 0 < m < b 0 σ.Further, the Jacobian matrix at the fixed point of the dynamical system, in terms of the parameters and the fixed point is given by,EQUATIONFrom the linear stability theory, we know that a fixed point (y s , a s ) is asymptotically stable when the real part of the eigenvalues of J are less than 0, i.e., Re(λ J ) < 0. For a 2D system, this is equivalent to the conditions: Tr(J) < 0 and det(J) > 0. The trace of the Jacobian matrix is given by,EQUATION)The determinant of the Jacobian matrix is given by,EQUATIONNow we consider the different conditions on w r and the input drive z for stability and state the various cases as theorems along with their proofs, E.1 Contractive constraint on recurrence (0 < w r ≤ 1)• z > 0 : There exists a unique fixed point with y s > 0 and a s > 0 and it is asymptotically stable.Proof. The existence and uniqueness of the fixed point, along with its stability properties, are established by Lemma E.4. Additionally, the positivity of a s (i.e., a s > 0) ensures that the expression 1 -w r + w r √ a s > 0. Consequently, given that z > 0, it follows from Eq. 78 that y s > 0.• z < 0 : There exists a unique fixed point with y s < 0 and a s > and it is asymptotically stable.Proof. Since this condition becomes equivalent to z > 0, up to a sign change in y (Eq. 75), it is straightforward to see why this is true. Lemma E.4. Given 0 < w r ≤ 1 and z ∈ R, there exists a unique fixed point and it is asymptotically stable.Proof. We observe that given 0 < w r ≤ 1, if there exists a valid fixed point, which satisfies a s > 0, Tr(J) in Eq. 82 is less than 0 (since 1 -w r + w r √ a s > 0 for any √ a s > 0), and det(J) in Eq. 83 is greater than 0 for all combinations of parameters and z ∈ R. Therefore, we need to find the constraints on the parameters that allow for at least one fixed point. The fixed point, m = √ a s , satisfies the quartic polynomial in Eq 80. Due to Theorem E.1, for a valid fixed point, we need to find positive real roots that satisfy m > 0.The sequence of signs of the coefficients of the polynomial when 0 ≤ w r < 1 is given by (+, +, ±, -, -). We use Descartes' Rule of Signs, which states the following: The number of positive real roots of a polynomial p(x) is either equal to the number of sign changes (omitting the zero coefficients) between consecutive non-zero coefficients of p(x), or it is less than this by a multiple of 2. Since there is exactly one sign change from left to right in the sequence, regardless of the sign of the coefficient of m 2 , we know that the equation above has exactly one real positive root for m. This further implies that a s has exactly one positive root and the corresponding fixed point (y s , a s ), is locally dynamically stable for all of the combinations of parameters and z ∈ R.Note that the result also holds for w r = 1. There is exactly one positive fixed point (y s , a s ) and it is given by, i.e., y s = bz/ b 2 0 σ 2 + wb 2 z 2 and a s = b 2 0 σ 2 + wb 2 z 2 .E.2 Expansive constraint on recurrence (w r > 1)Now we consider the case when w r > 1. We first find an alternative form of the determinant of the Jacobian (Eq.83) that is more amenable to the case, w r > 1,EQUATION)Rewriting the trace,EQUATIONThe properties of stability can be summarized in the following two cases,• z > 0 : There exists exactly one fixed point with y s > 0 and a s > 0 and it is asymptotically stable. Further, if b 0 σ > 1 -1/w r , there exist no additional fixed points. But if b 0 σ < 1 -1/w r , then there exist either two or no fixed points with y s < 0 and a s > 0 and they may or may not be stable.Proof. Since the fixed point satisfies Eq. 78, there are two possibilities, either 1 -w r + w r √ a s > 0 or 1 -w r + w r √ a s < 0. We consider them separately, ▶ 1-w r +w r √ a s > 0 : If z > 0, we have y s > 0 because y s satisfies Eq. 78. Further due to Lemma E.5, this fixed point is unique and asymptotically stable for all combinations of parameters. ▶ 1 -w r + w r √ a s < 0 : If z > 0, we have y s < 0 because y s satisfies Eq. 78. Further, if we are given b 0 σ > 1 -1/w r , no root exists. We prove this by contradiction, assuming a root exists that satisfies √ a s < 1 -1/w r , since b 0 σ > 1 -1/w r , this root satisfies 0 < √ a s < b 0 σ, but due to Lemma E.3 no such root exists, therefore we have a contradiction. Now if we have b 0 σ < 1 -1/w r , the root for m = √ a s must satisfy b 0 σ < m < 1 -1/w r . Since either one or three roots satisfy m > 0 (equivalently m > b 0 σ) in Eq. 80 and exactly one root satisfies m > 1 -1/w r (Lemma E.5), we conclude that number of roots for m in the interval (b 0 σ, 1 -1/w r ) are either two or none.• z < 0 : There exists exactly one fixed point with y s < 0 and a s > 0 and it is asymptotically stable. Further, if b 0 σ > 1 -1/w r , there exist no additional fixed points. But if b 0 σ < 1 -1/w r , then there exist either two or no fixed points with y s > 0 and a s > 0 and they may or may not be stable.Proof. Since this condition becomes equivalent to z > 0, up to a sign change in y (Eq. 75), it is straightforward to see why this is true.Lemma E.5. Given w r > 1, there exists a unique fixed point, (y s , a s ), that satisfies 1-w r +w r √ a s > 0 and it is asymptotically stable.Proof. The fixed point, m = √ a s , satisfies the quartic polynomial in Eq 80. Due to Theorem E.1, for a valid fixed point, we need to find positive real roots that satisfy m > 0.The sequence of signs of the coefficients of the polynomial is given by (+, -, ±, +, -). Regardless of the sign of the coefficient of m 2 , using Descartes' Rule of Signs, we find that since there are 3 sign changes, there are either 1 or 3 positive real roots for m. Since, we are given w r > 1 and 1 -w r + w r √ a s > 0, these roots are valid only if they satisfy 1 -w r + w r m > 0, or, m > 1 -1/w r .We make a transformation m → m + (1 -1/w r ) which gives us,EQUATIONThe number of valid roots is given by the number of positive real roots of this polynomial. The signs of the coefficients are given by, (+, +, ±, -, -). There is exactly one sign change, hence there is a unique fixed point that satisfies 1 -w r + w r √ a s > 0.Now we prove that this fixed point is asymptotically stable. It can be easily seen that Tr(J) < 0 in Eq. 85 and det(J) > 0 in Eq. 84 when 1 -w r + w r √ a s > 0 and w r > 1. Therefore, this unique fixed point is asymptotically stable.Here, we present the stability results for the model with only the positive part of the complementary receptive fields present. The model is given by the following dynamical equations,EQUATIONWe will follow the same procedure for stability analysis as we did for the main model and state the key steps in the various derivations for stability.We analyze the stability of the system when W r = I. Upon substituting y + → ⌊y⌋ 2 and a + → ⌊a⌋, we get the following dynamical system,EQUATIONThe fixed point is given by,EQUATIONAlso, ⌊y s ⌋ is given by,EQUATIONThe Jacobian matrix, J, about this fixed point is given by,EQUATIONHere 1 z≥0 is an indicator function that returns a vector the size of z with 1 at locations where z i ≥ 0 and 0 elsewhere. Similarly, 1 z<0 returns a vector the size of z with 1 at locations where z i < 0 and 0 elsewhere. Further J -λI is given by,EQUATIONSince A 11 and A 12 commute, we can write, det(J -λI) = det(A 22 A 11 -A 21 A 12 ). Upon simplification, we get,EQUATIONNote that, here we used the fact that y s has the same sign as z, as seen from Eq. 89. The dynamical system is stable if all eigenvalues of this characteristic polynomial have negative real parts. Just like the main model, we map this to a quadratic eigenvalue problem of the form L(λ) = det(λ 2 I+λB+K) = 0. Now, we see if the conditions of Theorem 4.1 are met. The stiffness matrix is given by,EQUATIONClearly, K is a positive diagonal matrix for all choices of parameters and input. Now we check if the B is a Lyapunov diagonally stable matrix, which is implied when B admits a regular convergent splitting, i.e., it has a representation of the form B = M-N, where M -1 and N have all nonnegative entries and M -1 N has a spectral radius smaller than 1.EQUATION)Since M is a positive diagonal matrix, all the entries of M -1 are nonnegative. Also, since all the matrices involved in the definition of N are nonnegative, therefore their product, N, is nonnegative. We are only left to prove that M -1 N = S has a spectral radius smaller than 1. Consider the matrix S,EQUATIONTheorem D puts a bound of 1 on the spectral radius of this matrix. Define t → 1/ 1 + (τ a /τ y ) ⊙ √ a s ⊙ 1 z≥0 + 1 z<0 , u → ⌊b ⊙ z⌋ 2 and v → b 2 0 ⊙ σ 2 , we notice that they follow the constraints of the theorem, therefore, S is convergent. This implies that B has a convergent regular splitting, therefore, the linearized dynamical system is unconditionally globally asymptotically stable (and nonlinear dynamical system is locally asymptotically stable) across all the values of parameters and inputs.Following a procedure similar to that in Appendix C, when all of the normalization weights in the system are equal, to value α, and the parameters are scalars, i.e., τ y = τ y 1, τ a = τ a 1, b 0 = b 0 1 and σ = σ1, we can write the analytical expressions for eigenvalues. The characteristic polynomial is given by,EQUATIONwhere, u = [1, 1, ..., 1] ⊤ , v = ⌊b ⊙ z⌋ 2 / σ 2 b 2 0 1 + W⌊b ⊙ z⌋ 2 ; n 1 and n 2 are the number of nonnegative and negative values, respectively, in the input drive z; δ 1 and δ 2 are given by,δ 1 = λ 2 + λ 1 τ a + σ 2 b 2 0 + α||⌊b ⊙ z⌋|| 2 τ y + σ 2 b 2 0 + α||⌊b ⊙ z⌋|| 2 τ y τ a δ 2 = λ 2 + λ 1 τ a + 1 τ y + 1 τ y τ a (98)Simplification of Eq. 97 gives us,det(J -λI) = δ 1 - λ τ a α||⌊b ⊙ z⌋|| 2 σ 2 b 2 0 + α||⌊b ⊙ z⌋|| 2 δ n1-1 1 δ n2 2 (99)Since the characteristic polynomial is a product of quadratic polynomials, we can solve them analytically. The strictly negative eigenvalues are given by,EQUATIONThe potentially complex eigenvalues are given by the solution to the following quadratic equation,EQUATIONF.3 Linear stability analysis of the two-dimensional modelThe dynamical system to consider is,EQUATIONSince a valid fixed point must have a s > 0, the fixed point (y s , a s ) satisfies,EQUATIONWe divide this into two cases,• y s > 0 : The fixed point is given by the equations,EQUATIONA fixed point, (y s , a s ), is valid only if it satisfies y s ∈ R + * and a s ∈ R + * . The Jacobian matrix is given by, Note that this is equivalent to the main model, with the additional constraint of y s > 0, whose stability analysis is presented in Appendix E.• y s < 0 : Eq. 103 & 104 yield us a unique fixed point y s = bz and a s = b 2 0 σ 2 . Since y s 0 and y s = bz, this is only possible when z < 0. The Jacobian matrix about (y s , a s ) is given by,J = -1 τy 0 0 -1 τa (108)Since the eigenvalues of J are λ J = -1/τ y , -1/τ a are real and both negative, this fixed point is always stable.Combining the cases above and results already established in Appendix E, we characterize the stability of the system as follows, F.3.1 Contractive constraint on recurrence (0 < w r ≤ 1)• z > 0 : There exists a unique fixed point with y s > 0 and a s > 0 and it is asymptotically stable.• z < 0 : There exists a unique fixed point with y s < 0 and a s > 0, given by, (bz, b 2 0 σ 2 ), and it is asymptotically stable. F.3.2 Expansive constraint on recurrence (w r > 1)• z > 0 : There exists a unique fixed point with y s > 0 and a s > 0 and it is asymptotically stable.• z < 0 : There exists exactly one fixed point with y s < 0 and a s > 0 given by, (bz, b 2 0 σ 2 ), and it is asymptotically stable. Further, if b 0 σ > 1 -1/w r , there exist no additional fixed points. But if b 0 σ < 1 -1/w r , then there exist either two or no fixed points with y s > 0 and a s > 0 and they may or may not be stable.In this section, we present an iterative approach to finding the fixed point for ORGaNICs with an arbitrary recurrent weight matrix. We show that this algorithm converges in a few steps (2-10) with great accuracy. We consider the system given by Eq. 15,EQUATIONThe fixed point of this system (y s and a s ) is found by solving the following simultaneous equations,EQUATIONEQUATIONThese equations do not admit a closed-form analytical solution when W r ̸ = I. We first find a good approximation for the initialization of y s and a s and then define the iterative algorithm. The equation for y s can be written in terms of a s as,EQUATIONNow applying the Woodbury matrix identity, which states thatEQUATIONto the inverse in Eq. 112 with A = I, U = I, C = D √ a s -I and V = W r , we get,EQUATIONWe approximate the above equation by assuming that W r is a symmetric matrix with the eigendecomposition given by QΛQ ⊤ , with Q ⊤ Q = QQ ⊤ = I and Λ = D (λ) is a diagonal matrix with the eigenvalues as its diagonal entries. This gives us the following approximation,EQUATIONWe approximate the eigenvalues, λ, by the maximum eigenvalue of the W r and assume the entries of √ a s are identical. This gives us the following initial guess for y s .EQUATIONFor the initial guess of a 0 s , we use Eq. 111 and plug in the following on the RHS a s → b 2 0 ⊙ σ 2 and the corresponding y s found by using Eq. 116. This gives us the following,EQUATIONNext we update the y s and a s by performing the following iterations derived using Eq. 110 & 111. For instance, (y 1 s , a 1 s ) are given by,EQUATION)This procedure is summarized in Algorithm 2. Substituting λ m = 1 in Eq. 116 & 117 yields simpler initial conditions and gives us Algorithm 1. Even though we had assumed that W r should be symmetric, in practice we find that this algorithm leads to fast convergence even for non-symmetric matrices. The fast convergence is owed to the fact that we have a good initial approximation of the solution. We also find that this iteration scheme works only for recurrent weight matrices with a maximum singular value of 1.Here, we find the energy (Lyapunov function) that is minimized by the dynamics of the ORGaNICs in the vicinity of the normalization fixed point. We consider the dynamical system with W r = I, which is given by Eq. 3. Upon linearizing about the fixed point we get the following linear dynamical system,EQUATIONAlgorithm 2 Iterative scheme for finding the fixed point 1: Input: ORGaNICs parameters and input (z), Tolerance ϵ, maximum iterations N 2: Output: Approximation to the fixed point (y s , a s ) 3: a ← σ 2 ⊙ b 2 0 + W λm-λ 2 m +λ 2 m √ b 2 0 ⊙σ 2 2 ⊙ b 2 0 ⊙ σ 2 // initialy ← I -W r + D √ a W r -1 (b ⊙ z) // y update 8: a ← b 2 0 ⊙ σ 2 + W y 2 * a // a update 9:k ← k + 1 10: end while 11: return (y, a)This system is dynamically equivalent (admits the same eigenvalues) to a system of coupled harmonic oscillators with the following equations,EQUATIONWe can rewrite the linear system in terms of the position, x, and the velocity, v,EQUATIONSince this system is of the form Iẍ + B ẋ + Kx = 0, Eq. 36, the energy of this dynamical system is given by V(z) = z ⊤ Pz, or, V (x, v) = x ⊤ v ⊤ TK ϵI ϵI T x v = x ⊤ (TK) x + v ⊤ Tv + 2ϵ x ⊤ v (122)Here, T is any positive diagonal matrix such that TB + B ⊤ T ≻ 0 and K is also a positive diagonal matrix given by D √ a s / (τ y ⊙ τ a ) . Now, for a valid Lyapunov function, we can take ϵ to be arbitrarily small, Eq. 51. Therefore, the energy minimized by the dynamical system is given by V (x, v) = x ⊤ (TK) x + v ⊤ Tv. This is a high-dimensional version of the energy of a damped harmonic oscillator. For a single oscillator, we have V (x, v) = t(kx 2 + v 2 ) which is proportional to the total energy (kinetic + potential) of the oscillator.We now express this energy in terms of the variables relevant to ORGaNICs, i.e., we find V (y, a). First, we denote the Jacobian matrices in RHS of Eq. 119 & 121 by A & B, respectively. We note the simple fact that A & B are related by a similarity transformation (a change of basis). This means that there exists an invertible matrix U, such that A = U -1 BU and the corresponding transform is given by [x v] ⊤ = U[y -y s a -a s ] ⊤ .Assuming that U is invertible, we can write this equation as UA = BU. To solve this, we consider a block matrix representation of U and find the following solution,EQUATIONThis change of basis gives us the transformation,EQUATIONor,EQUATIONSubstituting these expressions into V (x, v) = x ⊤ (TK) x + v ⊤ Tv, and assuming the diagonal entries of the matrix T to be t i and substituting the diagonal entries of K, k i → √ a si /(τ y i τ ai ), we get the energy in of y and a,EQUATION(126) We notice that Taylor expanding the term √ a i y i about √ a si y si and ignoring the second order terms, we get,EQUATIONTherefore the energy function, V (y, a), is given by,EQUATIONNotice that √ a si y si = b i z i . This gives us the following expression for the energy function,EQUATIONFurther, for an ORGaNICs model containing one y and one a neuron, after removing the proportionality constants, the energy function is given by,EQUATIONAfter plugging in the steady-state values, we get,EQUATIONFor this system, it is easy to verify that this is a valid Lyapunov function and is minimized by the dynamics of the circuit. We now demonstrate that it has the properties of a Lyapunov function. First, V (y s , a s ) = 0 and V (y s , a s ) > 0 ∀ y ̸ = y s & a ̸ = a s , this can be easily seen from Eq. 130. Second, we need to show that, V (y, a) < 0 ∀ y ̸ = y s & a ̸ = a s . Using Eq. 119 & 131 V (y, a), we can write the total time derivative of the energy to be, EQUATIONfor all the choices of parameters and ∀ y ̸ = y s , a ̸ = a s , we have V (y, a) < 0, therefore, it is a valid Lyapunov function and can be interpreted as the energy that decreases with time via the dynamics of ORGaNICs.The code (written in PyTorch [92] ) to produce all the results can be found at https://github.com/martiniani-lab/dynamic-divisive-norm. For both the static input and the sequential input, we train ORGaNICs on the MNIST handwritten digit dataset [72] , and to the best of our knowledge, it does not pose any privacy concern and has been used widely by the ML community freely. The simulations were performed on an HPC cluster. All of the models were trained on a single A100 (80GB) GPU. We use Adam optimizer [93] with default parameters for minimizing the loss function.We performed a random split of 57,000 training samples and 3,000 validation samples and picked the model with the largest validation accuracy for testing. To make a direct comparison to SSN [55] , we use the same architecture structure as theirs. First, we train an autoencoder (Table 3 ) to reduce the dimensionality of MNIST images to 40 by using a mean-squared loss function. Then, we use this 40-dimensional vector as to ORGaNICs and train it using the cross-entropy loss function. We additionally make the input gain b dependent on the input x, b = f (W bx x), where f is sigmoid. A layer of ORGaNICs is given by Eq. 15,EQUATIONThe "output" of a layer is the steady-state firing rate of the neuron with the positive receptive field, i.e., y + s = ⌊y s ⌋ 2 . We parameterize W r to have a maximum singular value of 1 and instead of simulating the dynamical system to find the fixed point, we use the iterative Algorithm 1 with a maximum number of steps = 10. More details about the parameters are given in Table 4 ; kaiming uniform initialization is used from [94] . Additional hyperparameters are given in Table 7 . We train ORGaNICs in a single-layer setting with the number of y neurons encoding the input, N 1 = 50, 80. We also train two-layer ORGaNICs (Table 5 ) with N 1 = 120 and N 2 = 60 neurons in each layer. The model is trained using backpropagation and takes approximately 10 min to fully train. We performed a random split of 57,000 training samples and 3,000 validation samples and picked the model with the largest validation accuracy for testing. The unpermuted sequential MNIST task is defined as follows: for a given 28×28 image, we flatten it to get a one-dimensional, 784 timestep-long input. Then these pixels are presented as an input (x i , one pixel at each time-step i) to the discretized ORGaNICs model given by the equation,y i+1 = y i + ∆t τ y ⊙ -y i + b i ⊙ (W zx x i ) + 1 -a + i ⊙ W r y + i a i+1 = a i + ∆t τ a ⊙ -a i + b 2 0,i ⊙ σ 2 + W y + i ⊙ a +2 i b i+1 = b i + ∆t τ b ⊙ i + f (W bx x i + W by y i + W ba a i ))b 0,i+1 = b 0,i + ∆t τ b0 ⊙ (-b 0,i + f (W b0x x i + W b0y y i + W b0a a i ))(135)When we are done presenting the pixels we use y + 784 = ⌊y 784 ⌋ 2 to make the predictions. To make this more challenging we also train ORGaNICs on permuted sMNIST where we first permute the pixels of all the images in some random order and the rest of the task is the same. Instead of parametrizing τ , we parametrize ∆t/τ y = 0.05 * f (p y ), ∆t/τ a = 0.01 * f (p a ) and ∆t/τ b = 0.1 * f (p b ) and ∆t/τ b0 = 0.1 * f (p b0 ), so we can control the dimensionless relative time constants. In practice, we find it is better to make the a neurons sluggish compared to y. This is based on the intuition given by the two-dimensional phase portrait for different relative time constants Fig. 3 . All the parameters (including W r ) are unconstrained for this task with initialization specified in Table 6 . Since ORGaNICs are stable, we did not need to use gradient clipping for training, which is commonly used for LSTMs. Additionally, we train the model using a StepLR learning rate scheduler with parameters given in Table 7 . The model is trained using backpropagation through time (BPTT) and takes approximately 30 hours to fully train. Note that since this is implemented in a feedforward manner, this is a cascading system with no feedback, hence we can perform the stability analysis of the two layers independently. 38th Conference on Neural Information Processing Systems (NeurIPS 2024).Python code for this study is available at https://github.com/martiniani-lab/dynamic-divisive-norm. 
 
 Appendix:  


Decision: This paper proves the unconditional stability of the ORGaNICs model, a biologically plausible RNN that implements divisive normalization. The authors provide rigorous mathematical analysis, particularly under specific conditions like the identity recurrent weight matrix, and demonstrate empirical stability across various tasks, including static and sequential classification benchmarks. Reviewers raised concerns about the broader applicability of the stability proof and the reliance on empirical evidence for more general cases. However, the authors addressed these concerns with additional experiments. While all reviewers voted for acceptance, the paper could benefit from clearer exposition and more comprehensive comparisons.

Reviews: review: paper_summary: paper_summary: Recurrent neural networks are widely applied in both solving machine learning tasks and modeling neural recordings.  However, they can be difficult to train due to exploding/vanishing gradients and other sources of instability, often requiring hyperparameter optimization.  In this paper, the authors argue that the ORGaNICs model is capable of solving sequential tasks without the need for hyperprameter optimization due to the conjectured property of ``unconditional stability,'' which guarantees the existence of a stable fixed point in the dynamical system for any model parameters.  The biological plausibility of the ORGaNICs dynamics also suggests that a similar dynamical system could guarantee stability in biological neural circuits for processing sensory information and generating motor sequences.

The stability property is proven rigorously in the case where the recurrent weights are equal to the identity and in the case where only a single recurrent and inhibitory unit are included in the model.  The authors provide new insights into the behavior of ORGaNICs by providing a Lyapunov function for the dynamics near the fixed point in the case where $\mathbf{W}_r = \mathbf{I}$.  An algorithm for quickly finding a stable fixed point of the ORGaNICs dynamical system in general is provided.  This algorithm is used to probe empirically for the fixed points of the ORGaNICs dynamics with general recurrent weights $\mathbf{W}_r$.

The authors test their prediction empirically in numerical experiments where they train an ORGaNICs network on two tasks: Static and sequential MNIST classification.  In both cases, they show that ORGaNICs obtains competitive performance without hyperparameter optimization.
strength_weakness: strength_weakness: The main theorem makes the very restrictive assumption that the recurrent weights $\mathbf{W}_r = \mathbf{I}$, such that they exactly cancel the intrinsic decay rate of the recurrent neurons.  While there are neural circuits where all recurrent interactions take place through inhibitory interneurons (e.g. cerebellar granule cells interacting through inhibitory Golgi cells), the intrinsic decay must also be taken into account.  If the assumption can be relaxed to say that $\mathbf{W}_r = \alpha \mathbf{I}$ for some $\alpha \in [0, 1]$ then I would be inclined to raise my score.

Also, the argument for unconditional stability of the ORGaNICs model is not well-supported in the paper's current form.  Beyond the case $\mathbf{W}_r = \mathbf{I}$ and the two-dimensional case, there this is tested only for the MNIST task with constant inputs.  In these tests, it also appears that the fixed-point-finding algorithm is used to quickly discover the stable fixed point.  However, it would be stronger to show that the dynamics naturally converge to the fixed point in a reasonable amount of time when started from randomly generated initial configuration.  In fact, it is possible that the dynamics naturally lead the system to a stable limit cycle and never converge to the stable fixed point, even if it is guaranteed to exist.  We ask the authors to clarify the methods used to locate the fixed point of the ORGaNICs network in their experiments and show that the system's dynamics tend toward this fixed point and not another limit cycle in a reasonable amount of time.  The numerical training also enforces a maximum singular value of $\mathbf{W}_r$, which is a strong condition required for stability.  ``Towards Unconditional Stability...'' would be a more appropriate title given the strong constraints required for stability guarantees in both the analytical proofs and numerical tests.

For the sequential MNIST task, the authors claim that they are able to achieve LSTM-like performance without hyperparameter optimization or gradient clipping/scaling.  During these experiments, the neural time constants are treated as parameters that can be learned.  This result is different than showing that for any fixed values of the time constants, ORGaNICs can be trained without issue.  The authors should verify that fixing the time constants and then training also leads stable training and no need for gradient clipping.  This is important to maintain a claim of biological plausibility where time constants are intrinsic properties of the neurons and not trained.

The static and sequential MNIST tasks also represent a rather limited scope in which to test for unconditional stability in practice.  It's possible (though I agree, unlikely) that for other more unwieldy tasks unconditional stability will not hold.  Either a proof of unconditional stability for general recurrent weight matrices or a wider range of numerical tests are necessary to make the case for unconditional stability of the ORGaNICs dynamics.
questions: questions: What hyperparameters are present in the ORGaNICs model as trained on the static and sequential MNIST tasks?  While the time constants are learned, how is the learning rate for naive BPTT selected?  Would this hyperparameter plausibly need to be optimized to achieve performance competitive with a well-tuned LSTM?

In the verification of the iterative fixed-point finding algorithm, how is the ``true'' fixed point determined?  By running the dynamics from a random starting point to the fixed point, or by running the dynamics from an initial estimate obtained using the iterative fixed point finding algorithm?
limitations: limitations: Even for static inputs, there are scenarios in which a biological circuit may want to produce a non-static output.  For example, in central pattern generators of motor sequences. If indeed ORGaNICs always converges to a stable fixed point, this may represent too constrained of a circuit to explain many interesting neural phenomena.
paper_summary: paper_summary: The paper discusses the development and analysis of "Oscillatory Recurrent Gated Neural Integrator Circuits" (ORGaNICs), a biologically plausible model of recurrent neural networks that implements divisive normalization. The authors prove the unconditional local stability of ORGaNICs with an identity recurrent weight matrix using the indirect method of Lyapunov function. They also demonstrate empirical stability for higher-dimensional circuits. The model's performance is evaluated on static and dynamic classification tasks, showing comparable results to LSTMs without specialized training strategies, thanks to its intrinsic stability properties.
strength_weakness: strength_weakness: 1. Generalization: The paper does not extensively discuss how well the results might generalize to more complex or different types of tasks beyond the tested benchmarks.
2. Performance on Benchmarks: The proposed model performance still has a gap compared to machine learning models.
3. Scalability: parameter size scales too fast with model size and there is no obvious way for impovement.
questions: questions: 1. Can this model provide more interpretability in relevant machine learning tasks? 
2. Besides normalization and stability, can this E-I model offer more insights into neuroscience compared to machine learning models?
3. Many existing machine learning models have already implemented various forms of divisive like normalization such as layer normalization. What are advantages or innovations of the divisive normalization implemented in this model.
4. In comparison to standard machine learning models, biophysically plausible models typically offer superior interpretability. However, the paper offers limited discussion on this topic.
limitations: limitations: 1. In the sequence modeling benchmarks, only 2 tasks are evaluated. It is recommended to include a broader range of sequence modeling tasks to allow for more comprehensive comparisons, such as those where SSM excels in long sequence modeling or where LSTM performs well in formal language tasks like the parity task.
paper_summary: paper_summary: This paper studies the stability properties of a model of cortical circuits which was introduced in 2019 (and I wasn't yet aware of): the ORGaNIC model by Heeger & Mackey. This LSTM-like model uses a simple set of differential equations that unify several phenomena observed in cortex, including normalization (Carandini & Heeger, 2012), and in some sense map LSTMs onto plausible cortical circuitry. This paper presents two main theoretical results: (i) that the model is locally stable around its unique fixed point when the recurrent matrix is the identity, and that (ii) existence and uniqueness of a stable fixed point can be shown for a general 2D model (really just one "principal neuron", and a gate variable) under certain conditions on the parameters. The paper wraps up with a few training experiments on (regular/sequential/permuted) MNIST classification, showing that ORGaNICs perform well.
strength_weakness: strength_weakness: Overall this is a highly esoteric paper for which I have difficulty assessing potential impact (hence my high uncertainty rating). The theoretical results are impressively detailed but seem fairly limited in scope (how are the W_r = I and 2D cases relevant for either ML or neuroscience?). I am also not super convinced by the utility of the broader conjecture on stability at the end of section 5. The experiments are fairly limited in scope and breadth, too (e.g. no mention of hyperparameter tuning for the SSN comparison).

The abstract claims that it is “thanks to its intrinsic stability property” that ORGaNICs perform comparably to LSTMs; while this seems like a sensible hypothesis I don't think the paper really shows that, yet I believe the paper would be much stronger if that was shown to be true. In this respect, I wonder in what way the comparison to SSNs (e.g. Soo & Lengyel) helps make this case; is it because SSNs also perform a form of normalization at steady-state, much like ORGaNICs, but have otherwise no stability guarantees? Is it even true that SSN training is brittle / often fails because the network loses stability, thus giving rise to a stronger-than-usual tradeoff between training stability and learning rate magnitude? How could the author rule out that there might be another fundamental difference between these 2 models that has little to do with stability and yet underlies the performance difference? After all, no formal stability guarantees exist for standard LSTMs as far as I am aware, and that doesn't prevent them from training very well on most tasks (was gradient clipping even necessary for the LSTMs in those particular experiments?). In summary, I suppose the paper lacks a couple of convincing ablations to make the point that stability guarantees increase training performance.

On clarity of exposition: this is a pretty hard paper to follow, it's very dense and doesn't really offer the hierarchy of exposition that a reader would want to see -- the theoretical results are presented in a very flat way, with a permanent back and forth between details and bigger picture that I found hard to navigate.
questions: questions: - The abstract (and main text) says you trained your models using "simple backpropagation through time" but appendix I says you trained all models using Adam -- in principle there is no contradiction between these two claims, as backprop is also involved in Adam, but I wonder why you wrote "simple backpropagation" (which could easily be misunderstood as "simple gradient descent", i.e. not Adam). What's non-simple backpropagation?

- The formatting of references is annoying :D please use parentheses to clearly separate references from the main text.

- typo on l104: "dynamics-neural growth" → "dynamics-neutral growth"
limitations: limitations: Some technical limitations are stated in the Discussion section.

*EDIT*: following the rebuttal, I am raising my score to a 5.
paper_summary: paper_summary: This work is based on ORGaNICs, a particular type of RNN architecture that implements divisive normalization (from neuroscience). The authors explore whether ORGaNICs are stable enough to be meaningfully (and stably) trained by gradient descent. Due to their stability, which the authors have proven theoretically for specific cases and empirically in general, the authors claim that ORGaNICs can learn long-term dependencies and solve vanishing and exploding gradients found in other neurodynamical models.
strength_weakness: strength_weakness: Subjectively, I believe that the scope of this work is too limited for the audience at NeurIPS, but I will not consider this in my final score. 

The general theme of my objective issues with this work is the slight overselling of this work:

1. Fundamentally, ORGaNIC is an architecture. Song et al. [2016] implemented Dale's law by manipulating the recurrent weight matrix, which is a **method** applicable to any architecture. Soo et al. [2024] trained vanilla RNNs on long-term dependencies using skip connections through time, which once again is a **method** applicable to any architecture. So the claim that these "**models** do not accurately reflect cortical dynamics" in lines 105-106 do not make sense to me. To illustrate my point, can the authors comment on whether there are any issues if Soo et al. [2024] is applied to ORGaNICs (beyond the fact that it is not needed)? 

2. This work seems to be strongly driving the point that divisive normalization is a key criteria for biological plausibility. There are a variety of biological properties that the brain exhibits. For example, there are 22 different properties of V1 identified in Marques et al. [2020], of which specific versions of DN make up some of them. It is not convincing to select a model based on one particular effect found in V1. 

Marques et al. [2020]: https://www.biorxiv.org/content/10.1101/2021.03.01.433495v2 

3. The authors claim that ORGaNICs is stable and therefore solves exploding and vanishing gradients, which does not make sense. The way the authors used the word "stable" here is in the language of dynamical systems. "Stability" in the context of exploding and vanishing gradients has a different meaning. The fact that the model can learn sequential MNIST (gradient stability) does not mean that it is due to its (dynamical) stability. A completely (dynamically) stable architecture can still have vanishing gradients. Lines 115-116 suggest to me that the authors are somewhat aware of this but not completely.
questions: questions: ORGaNIC begins with the word "oscillatory". Can the authors comment on whether the network must always be oscillatory? I see that the analysis looks into stable and unstable fixed points, which means that the answer is no, which makes the name confusing to me.

Likewise, can the authors elaborate on the "gate" in ORGaNIC and how it is the "biophysically plausible extension" (line 86) of gating in machine learning? Specifically, why it is biophysically plausible, and why it is an extension.
limitations: limitations: In line with my issue of overselling this work, there really is not any honest discussion of genuine limitations in the discussion section, which the authors claim they have done so in the checklist. Discretization problems (lines 350-354) and using single-layers (line 355) are not specific to this work at all. I would encourage the authors to talk about genuine limitations. For example, I suspect but cannot prove, that there would be some increased computational complexity compared to vanilla RNNs, since having multiple differential equation for a single model means that more things need to be done in a single time step.

This work, focusing on a biologically plausible model, shows the model performing machine learning experiments (in contrast, in Heeger and Mackey [2019] the model is performing neuroscience experiments). 

Also, while this model implements divisive normalization, as I mentioned above there are many biological phenomena that the brain exhibits. It would be unreasonable to ask the authors to state what ORGaNICs cannot do, but at least I would like to ask the authors if ORGaNICs can do everything that SSNs (which they cited) can do? 

These are the limitations that I believe readers would like to see.
paper_summary: paper_summary: Stability is a critical notion in the understanding of dynamical systems as well as for learning them. This paper studies the stability of a biologically plausible recurrent cortical model, ORGaNICs, that implements divisive normalization. More precisely:
- it demonstrates the local stability of two specific subclasses of the model.
- it shows that the model performs competitively on standard machine learning tasks.
strength_weakness: strength_weakness: The main weakness of the paper is that it might not be fully suited to the NeurIPS community in its current form:
- the introduction of the dynamics lacks some intuition on what the dynamics achieves. Equation 2 partially does so, but is likely not sufficient to fully convey the intuition. If possible, connecting to existing RNN architectures would help. Introducing the Lypanuov energy earlier could also help.
- many notions used in Section 4 might deserve a better introduction as NeurIPS is not a dynamical systems conference. This holds for "Lyapunov diagonally stable", "Z-matrix", "M-matrix" and, to some extent, to "stability" and "indirect method of Lyapunov". To make their results more accessible, the authors may want to provide intuitive definitions of these terms and why they are important.

Additionally, the connections made to the machine learning literature are sometimes imprecise:
- in L43, it is written that stability helps mitigate the vanishing gradient problem, which is inaccurate (e.g. a system with very fast time constants is stable but suffers from vanishing gradients). On top of that, stable recurrent networks having slow time constants can still suffer from some form of exploding gradients [Zucchet and Orvieto, 2024](https://arxiv.org/abs/2405.21064).
- in L61, the authors write that divisive normalization generalizes batch and layer normalization. While there are definitely some links, this is not true: divisive normalization affects the recurrent dynamics whereas both batch and layer normalizations are usually applied after the recurrent layer. In particular, those normalization schemes do not affect the stability of the system.
- in L349, divisive normalization is directly compared to the form of attention used in Transformers. Transformers normalize over the time axis, which is highly implausible and different from the type of attention mentioned here. This point deserves to be made more precise.

Finally, there are some additional links to existing ML literature that would nicely fit within the paper:
- the kind of neural network studied in Section 6.1 is known as [Deep equilibrium models, DEQ](https://arxiv.org/abs/1909.01377). One of the main difficulties that come with training such models is to [keep the dynamics stable](https://arxiv.org/abs/2106.14342) throughout learning. The studied architecture may be less prone to these behaviors: adding the DEQ baseline to the experiments and monitoring the stability of the dynamics during training would be interesting additions to Section 6.1.
- in Section 6.2, the baseline architectures are rather "old". The state-of-the-art networks on tasks requiring modeling long-range dependencies are state-space models (e.g. [S4](https://arxiv.org/abs/2111.00396) or [LRU](https://arxiv.org/abs/2303.06349)). Additionally, sequential models are known to be extremely sensitive to initializations on tasks such as sMNIST, see [Amos et al. 2023](https://arxiv.org/abs/2310.02980). It is therefore important to describe the initialization scheme of the recurrence in the main text, as it is likely a main driver of the performance. 

I hope these remarks can help the authors to improve their paper.
questions: questions: - The model considered has the same number of E and I neurons, which is, to the best of my knowledge, far from being the case in the brain. Is the system still stable when the number of I neurons is much smaller than the number of E neurons?
limitations: limitations: Limitations are properly addressed.

Comments: title: Thank you again
Updating statements and summary of key results
Thank you
Good response
Comment
Discussion of limitations
Raising Score
Early response
comment: I thank the authors for all the time and effort put into addressing my remarks, and I appreciate them for their diligence.
As we noted in our original response, the model’s limitations had already been stated in the original manuscript. We will summarize them in the Discussion.

The paper investigates a neurobiologically plausible RNN model that achieves divisive normalization (DN) through recurrent excitation. This model has been shown to recapitulate a broad range of neurophysiological observations, but the application of this theory to neuroscience is not the focus of the current paper. The core contributions can be summarized at a lay level as follows:
- **Unconditional Stability:** ORGaNICs is shown to be unconditionally stable under mild constraints (this is highly nontrivial). We develop mathematical machinery to prove this in a couple of limiting cases, by expressing ORGaNICs as a mechanical (specifically, gyroscopically stabilized) system – to do this we even need to prove new results on mechanical systems.
- **Implications for ML:** This stability allows ORGaNICs to be trained using naïve BPTT without gradient clipping/scaling, performing comparably to LSTMs. The point is not simply that ORGaNICs is at least as good as LSTMs, but that imposing DN dynamically makes training more efficient and robust.  Moreover, ORGaNICs trained by naïve BPTT fares well when compared to SSN trained by a DNG (the comparison is to SSN, not to DNG).
- **Dynamic Normalization:** Unlike standard RNNs, where normalization is imposed a-posteriori to the output layer, ORGaNICs implements normalization dynamically: this is the first time this idea is tested in ML and analyzed theoretically. We believe that this is an important contribution. In contrast layer/batch normalization does not affect the stability of the system.
- **Theoretical Insights:** The connection to mechanical systems enables deriving a Lyapunov function, offering a normative principle for understanding the model's dynamics. 

Much of the reviewer’s concerns are focused on the language used in the Related Work section which we will rephrase as follows (adapting the reviewer’s first suggestion):

(Song et al., 2016) incorporated Dale’s law into the vanilla RNN architecture, which was successfully trained across a variety of cognitive tasks. Building on this, (Soo et al., 2024) developed a technique for such RNNs to learn long-term dependencies by using skip connections through time. ORGaNICs is a model that is already built on biological principles, and can learn long-term dependencies intrinsically, therefore it does not require the method used by (Soo et al., 2024).
 
Regarding the remaining concerns. The reviewer raises concerns about modeling the full range of V1 phenomena, but V1 is not the focus of this paper. There is no mention of V1 or the visual system anywhere (except when discussing the history of DN). Many of the phenomena mentioned by the reviewer have been or can be incorporated into ORGaNICs (see paragraph below). Application of the ORGaNICs theory to model V1 and several other neural systems will be done in future work. Regardless, we welcome the reviewer’s suggestion and we will add the following text:

DN has been proposed as a canonical neural computation (Carandini & Heeger, 2012) and is linked to many well-documented physiological (Brouwer & Heeger, 2011; Cavanaughet et al., 2002) and psychophysical (Xing & Heeger, 2000; Petrov et al., 2005) phenomena. DN models diverse neural processes: adaptation (Wainwright et al., 2002; Westrick et al., 2016), attention (Reynolds & Heeger, 2009)), automatic gain control (Heeger, Simoncelli & Movshon, 1996), decorrelation and statistical whitening (Lyu, & Simoncelli, 2009). Since ORGaNICs’ response follows the DN equation at steady-state, it already incorporates this wide variety of neural phenomena. ORGaNICs have been shown to capture some of the dynamics of neural activity (Heeger & Zemlianova 2020; Rawat et. al., Cosyne 2024). Additional phenomena not explained by DN (Duong et. al., NeurIPS 2024) can in principle be integrated into the model. Regardless, more work needs to be done, of course, to explain the full range of neurophysiological phenomena. In this paper, however, we focus on the effects of DN on the dynamical stability of ORGaNICs.
 
When we state that “we have no evidence (theoretical or experimental) of additional limitations” we are referring specifically to ORGaNICs as an instance of RNNs for ML (because that is the focus of this paper), not to limitations of the theory as a neurobiological circuit model (which we intend to publish separately). We have no evidence that there is something that ORGaNICs couldn’t do when compared to an LSTM.
  
Finally, we assure that the paper has been written with the most genuine intentions, purely driven by curiosity. In fact, the core contributions are conceptual and technical (in the form of theorems), rather than ad-hoc benchmarks. We are grateful for the reviewer’s constructive comments which we believe will help clarify the key contributions of the paper.
As the discussion period is coming to an end, I realize that this back and forth might have required too much time and effort from the authors, so I will post my concluding remarks now, giving the authors full benefit of the doubt in outstanding issues.

Overall, the main rebuttal has addressed many doubts that I have, and I will raise my score by 1. This work represents progress for previously handcrafted models in neuroscience to be trained and actually perform tasks, which is a meaningful contribution to neuroscience.
We thank the reviewer for their feedback. Regarding the exposition, we will enhance the paper by providing clearer, more intuitive explanations of terms like “Lyapunov diagonally stable,” “Z-matrix,” “M-matrix,” and “indirect method of Lyapunov,” as suggested by reviewer ma5B. However, keeping the technical details is crucial. One of the core innovations of the paper is the technical approach that we developed to prove the stability of the recurrent circuit model (ORGaNICs), connecting the linear stability of a neural circuit model to the dynamic analysis of mechanical systems. As such, we are not merely going off on a tangent, but achieving a rigorous proof of stability for a neural circuit model by linking two disparate fields. This has important implications for the derivation of a Lyapunov function and thus a normative principle that makes ORGaNICs interpretable. We need to be technical enough to make this link precise and for the description of the approach to be correct and reproducible. We note that there are many mathematical papers published in NeurIPS that are richer in theorems and technical details/jargon than this work. We tried to strike a balance to reach a broader audience, so we welcome the reviewer’s comment and will clarify the writing further.

Regarding the comparison to SSN, the results by Song 2022 were not hyper-tuned either. We used the same learning rate, optimizer configuration, and equal number of layers and number of units in each layer to make the comparison fair. We will state this clearly in the manuscript. Moreover, we will conduct hyperparameter tuning to facilitate future comparisons and ensure our results are robust.
 
Nevertheless, we would like to underscore once again that the numerical comparison of ORGaNICs with SSN or LSTMs are not the main point of the paper. Specifically, the core contributions of the paper can be summarized at a lay level as follows:
 
- We take an existing neurobiologically plausible RNN model designed to achieve divisive normalization (DN) exactly via recurrent excitation. This model has been shown to recapitulate a broad range of neurophysiological observations, not by design but as a result of imposing a specific circuit function (i.e., DN). Given the resemblance to LSTMs, we hypothesize that this model should also do well on ML tasks and ask what is the impact of imposing DN on the dynamics and trainability of an RNN. 

- We discover empirically that ORGaNICs is unconditionally stable under very mild constraints (this is highly nontrivial). We thus develop mathematical machinery to prove this unconditional stability in a couple of limiting cases (d=2 for all parameter values, and $\mathbf{W}_r = \mathbf{I}$ for generic d and all other parameter values) by expressing ORGaNICs as a mechanical (specifically, gyroscopically stabilized) system – to do this we even need to prove new results on mechanical systems.

- But what are the implications of this unconditional stability, if any? Stability is crucial to ensure the trainability of RNNs that typically require gradient clipping/scaling. We discover that by virtue of its stability, ORGaNICs can be trained on a standard ML task (sMNIST and psMNIST) by naïve BPTT, without gradient clipping/scaling. In fact, naïve BPTT  works so well that on a first attempt, it does about as well as LSTMs (trained with gradient clipping/scaling and hyper-tuned) – it is not luck, the result is robust with respect to re-initialization of the training. The point then is not simply that ORGaNICs is at least as good as LSTMs, but that imposing DN dynamically makes training effortless. Note that unlike standard RNNs, where normalization is imposed a-posteriori to the output layer, ORGaNICs implements normalization dynamically: this is the first time this idea is tested in ML and analyzed theoretically (as pointed out enthusiastically by reviewer ma5B). Moreover, ORGaNICs trained by naïve BPTT fares well when compared to SSN trained by a specialized technique. The reason we choose SSN is that this model is also neurobiologically inspired and it has also been shown to approximate normalization in certain parameter ranges.

- Finally, the connection to mechanical systems allows us to derive a Lyapunov function that provides a normative principle for the dynamics of the models (i.e., a means of interpreting why and how the model tends to the stable fixed point).
 
We will better summarize the paper’s key contributions in the paper’s introduction and discussion.

We kindly ask that the reviewer assess our paper on its own merit, and not “reflecting a comparative rating taking into account the strength of the field across all my review assignments”. This request is in line with the NeurIPS FAQ that states “Q: Can I accept or reject all the papers in my stack?” “A: Please accept and reject papers based on their own merits. You do not have to match the conference acceptance rate.”
Thank you for your rebuttal. I can see the potential impact of this work a bit more clearly now, and will raise my score to a 5, also reflecting a comparative rating taking into account the strength of the field across all my review assignments. Regarding clarity of exposition, I would encourage the authors to do more than just “revis[ing] the introduction of the paper to provide a clearer roadmap of the main results of the paper” -- my concern is that the main text itself is hard going, with the authors often going off on a tangent, and including a plethora of very technical details where the reader would instead like to be given a clear synthesis. Also:

> There was no hyperparameter tuning for any of the tasks – ORGaNICs simply outperforms SSNs on the first attempt at the task we considered.

Surely the comparison to SSN is meaningless unless you (at least) hyper-tune the SSN results; otherwise, how do we know if you weren't just lucky here? In general though, I think it's a good idea to hyper-optimize even one's own method, as this saves a lot of time downstream to people who would like to compare to your method. For example, if your answer to my concern above is that you, in fact, did not run your own SSN simulations but directly imported the (hyper-tuned) results from a previous SSN paper, then you are relying on SSN authors having done their job, but you are not doing yours.
I thank the authors for the follow up response. I am glad that authors finally understand the position I am coming from, so I will continue with additional comments.

My main point, since the very first review I posted here, is that the tone of the work is a little overblown, citing the following examples:
1. criticizing method papers that are not specific to any architecture to promote your own model (point 1 of my weaknesses that remains unaddressed)
2. absolutely judging every model based on whether it can or cannot do DN (point 2 of my weaknesses that the authors doubled down on in the reply, highlighting the importance of DN)

which leads me to feel like the paper is not coming from a genuine standpoint of trying to contribute to academia, but rather trying to sell itself. Also, it is appalling that the authors can claim "we have no evidence (theoretical or experimental) of additional limitation" when V1 has so many interesting phenomena (see bolded below of biological traits that ORGaNIC does not have). Again, from my purely sincere point of view for this work to read like a good paper, might I suggest a more contributing narrative that fosters collaboration and improvement as a field instead of relying on criticism:

- the vanilla RNN architecture was adapted to be biologically realistic by [Song et al. 2016] by incoporating Dale's Law. [Soo et al. 2024] developed a method for such RNNs to learn long-term dependencies. ORGaNIC is a model that is already built on biological principles, and can learn long-term dependencies intrinsically, therefore not needing any of those results. 
(those papers provided methods, not models to be compared with ORGaNIC)
- DN is a phenomenon found in the brain, and ORGaNIC is able to express this effect. The biological visual system also expresses other phenomena, such as **surround suppression, adaptation, attention-based mechanisms, foveal and peripheral vision, retinotopic mapping, columnar structure, binocular and monocular effects, unique processing of color and plasticity**. Other models have been built to study those traits, but here we focus on DN as it is one of the more interesting effects that is believed to stem purely from dynamical effects.

(Edit: To be clear, these are examples of how to write a collaborative narrative. I am in no way forcing the authors to include anything in the paper.)

Once again, I am glad that the authors are now giving genuine limitations, and I hope that they can continue to accept my constructive criticism on the points above. I will reply again right before the deadline to make my final decision.
We thank the reviewer for their feedback and for increasing the score.
In addition to the general limitations (discretization problem and using single-layers), we also discussed other limitations throughout the text. For clarity, we will summarize them again in the Discussion in a new version of the manuscript. These are,

- **Limitation #1:** We pointed out very clearly that we cannot rigorously prove stability for the case of a general recurrent weight matrix. So our most general results rely on empirical evidence obtained through extensive numerical tests presented in the original manuscript, as well as additional results presented in response to the reviewers (see Fig.1 of the attached pdf for a demonstration of stability from random initial conditions).

- **Limitation #2:** We are not yet taking advantage of the modulators to control the effective time constant and instead are learning the intrinsic time constants. This potential limitation was pointed out by reviewer NGE1 who asked for further validation of the model with fixed time constants. We welcomed the suggestion and performed further numerical tests in which we found that ORGaNICs achieves good performance even when we fix the intrinsic time constants (as shown in Table 1 of the attached pdf).

- **Limitation #3:** In the current work, the weight matrices $ \mathbf{W_{by}}$, $ \mathbf{W_{ba}}$, $ \mathbf{W_{b_0 y}}$ and $\mathbf{W_{b_0 a}}$ are $n \times n$. This leads to the number of parameters increasing faster than other RNNs with hidden state size, as is evident from Table 2. This was pointed out by reviewer 2CUH. We will stress this limitation further in the updated version of the work. In future work, we plan to test the performance with compact and/or convolutional weights so that the # of parameters does not increase markedly with the hidden state size.

To the extent that we tested the model for the purpose of demonstrating stability, and trainability by SGD without gradient clipping/scaling in a standard ML task, these are the only limitations that we encountered. We will add a section to the Discussion to summarize them clearly. 

Based on our current understanding and results, we have no evidence (theoretical or experimental) of additional limitations. In future work, we will address the challenges noted above, and benchmark the model extensively on additional ML tasks. If the reviewer has specific questions that we have not thought of, we would be happy to investigate them and perform additional numerical tests as we did in response to other reviewers.
The authors have addressed most of my concerns with added experiments.  I believe that this is a significant contribution even without a proof of unconditional stability in the most general case.  I've raised my score to a 6.
I thank the authors for the thorough response. I will reply again in a few days after carefully reading them. But I need to post this early response first to give the authors a chance to reply.

In my original review, I mentioned the point about the lack of any honest discussion of limitations, and then suggested three possible points for the authors to discuss as a limitation. Perhaps the authors felt the need to defend them as if they were criticisms, which is why I am writing this post to clarify: those points are not criticisms and I am truly encouraging them to give limitations of their work just like in any complete piece of research. 

Basically I am saying that in the "limitations" part of the response, I still do not see any limitations.
I acknowledge the rebuttal and keep my score as it is.