Like E/I balance, synaptic sparsity has also been linked to disease states, and synaptic pruning has even been hypothesized as a major underlying mechanism in Schizophrenia and ASD, with excessive pruning believed to occur in Schizophrenia and - inversely - insufficient pruning in ASD [Salter and Stevens, 2017, Sakai, 2020]. While there is likely a great deal of nuance in how this manifests due to the massive amount of heterogeneity between brain regions and even between different neuron types, that is all the more reason for closer study of sparsity in theory and particularly in a multi-area RNN context, which will of course be intertwined with study of stability properties. Pruning-related hypotheses are so compelling because they can connect observations from across genomics, immunology, developmental neuroscience, connectomics, and cellular neurobiology with ideas from cognitive science, machine learning, and theoretical neuroscience, thereby creating an opportunity for rich interdisciplinary crosstalk and model iteration.

The relevance of sparsity has indeed been a discussion topic in neuroscience theory for many decades, including results on the energy efficiency of sparse coding, which reduces the total amount of neural activity needed to represent information [Levy and Baxter, 1996], as well as results on learning rules that can lead to sparse representations [Foldiak, 1990]. Of course when taken to an extreme sparse coding is not without its drawbacks, and it is not the case that concepts are each represented by activity of a single neuron in the brain. However, neuroscientific evidence does suggest that the brain leans more towards a sparse code than a dense one, and a number of computational neuroscience works have demonstrated that observed neural responses can be emergent properties of a nonnegative sparse coding scheme [Beyeler et al., 2019]. [Beyeler et al., 2019] additionally point out that the degree to which neural activity matches a sparse representation model can vary by brain region, yet another argument for greater study of "networks of networks".

Despite the fact that sparse coding has been a popular topic in the computational neuroscience literature, there has only recently been a reemergence of works focused on sparsity in the weights of trainable RNN models [Khona et al., 2022]. This is likely a result of the practical success of artificial neural networks over the last decade, which have in large part focused on fully connected models. While it is known that sparsity can improve generalization ability (for reasons akin to Occam's razor), and regularization techniques like $L1$ loss are directly related to sparsity, most deep learning architectures remain ultimately very densely connected [Goodfellow et al., 2016]. A subfield of deep learning research that has recently increased in popularity focuses on achieving strong (albeit often slightly worse) performance on benchmark tasks using networks that are much smaller than the typical state of the art architectures employed. This might involve building up networks with fewer units and fewer weights between them via an evolutionary inspired algorithm [Gaier and Ha, 2019], or it might involve "pruning" the majority of weights based on early training results [Frankle and Carbin, 2019]. It also can involve architecture design that draws from theoretical and/or neuroscience inspired principles, such as the work presented in this chapter.

Ultimately, there are many practical reasons for a deep learning researcher to care about sparsity of network weights, including efficiency of training, ability to easily deploy a trained network on weak/small hardware, and the interpretability of network behavior. It is plausible that artificial neural networks could one day be a common tool not only for predictions, but also to assist data scientists in forming mechanistic insights. However it will first be necessary to substantially further the mathematics and the science of deep learning; sparse network topologies will be an important component to understand for this pursuit, and given the clear connections to neurobiology this is an obvious area for computational neuroscientists and machine learning experts to join forces.

Understanding neural synchronization.
Besides stability analysis, control theory also has deep ties to study of system synchrony – a recurring topic in neurobiology and likely a critical topic to elucidate if we hope to understand sleep, itself a critical component of healthy (and diseased) human functioning. Because contraction analysis defines stability in terms of convergence of different system trajectories to each other as opposed to convergence of system trajectories to equilibrium states, it is possible to pose questions about whether connected systems will synchronize instead as questions about whether some "auxiliary system" representing a difference between their states will contract (see section <ref> for more on contraction).

Sleep broadly synchronizes neural activity, with certain stages of sleep driving slow oscillations at frequencies that are otherwise not usually observed in a healthy, alert brain [Timofeev and Chauvette, 2013]. It is well known that sleep plays a key role in learning and memory, with experience replay observed via neural activity in sleep, and clear cognitive deficits measured with even a moderate lack of sleep [Alhola and Polo-Kantola, 2007]. Further, it is likely that sleep plays some role in maintaining neural stability, possibly by "cleaning up" neural state while closing off inputs/outputs; interestingly, sleep deprivation has been a well-documented seizure trigger for Epilepsy patients and animal studies have demonstrated that REM sleep has a protective effect against seizures [Malow, 2004]. It can be difficult to probe the causal direction of sleep perturbations, so it would be immensely helpful to gain a better understanding of how synchronization can theoretically affect both network state and weight updates (and vice versa) through nonlinear dynamical systems modeling.

Unsurprisingly, an array of different sleep abnormalities have been found across psychiatric disorders. While it is of course possible that altered sleep dynamics are an outcome of psychiatric disorders rather than a cause, it is more plausible that causation runs in both directions, generating an unwanted positive feedback loop in patients. For example, poor sleep can increase the severity of hallucinations in Schizophrenia, and psychotic episodes can even occur in otherwise healthy adults due to severe sleep deprivation [Waters et al., 2018]. Interestingly, sleep abnormalities seen prior to a seizure (e.g. abnormally long sleep spindles in slow wave sleep [Picard et al., 2007]) are the opposite of those typically seen in Schizophrenia patients (i.e. shortened sleep spindles [Ferrarelli et al., 2007]). This is consistent with the more general "mirror image" phenomena of Autism and related disorders having many inverse properties of Schizophrenia and related disorders, most notably results related to hyper- versus hypo- connectivity and excitatory versus inhibitory imbalance. It perhaps hints that architectural differences may lead to the differing sleep dynamics, which may in turn be a causal factor for associated symptoms.

Besides sleep, there are additional emerging reasons to be interested in periodic, synchronizing inputs to neural systems. The use of flashing lights or tones of a particular frequency can entrain brain activity, and may have effects such as clearing protein aggregates associated with Alzheimer’s Disease [Iaccarino et al., 2016]. A theoretical understanding of how such inputs may alter brain activity and even plasticity would be beneficial in guiding future research on behavioral and medical applications of this type of protocol. Along the same lines, invasive procedures to alter brain activity - most notably deep brain stimulation (DBS) - often involve periodic inputs in practice as well. More broadly, DBS can be thought of as a control system, and thus control theory tools have obvious potential relevance. Recall that the DBS system trialed in chapter <ref> aimed to disrupt hypothesized corticostriatal hypersynchrony in OCD by driving ventral striatal and motor cortical regions at slightly different frequencies, but inadvertently increased coherence empirically observed between the regions in the pilot subject.

Design of DBS systems is of course highly complex in practice, and even the settings of established paradigms require experimental fine tuning for each patient. At the same time, DBS often involves chronic stimulation at a particular frequency, including the psychiatric treatment applications with the most evidence to date in OCD and depression. As such, theoretical study of model system synchronization behavior in response to a recurring input could indeed prove useful for future understanding of why these DBS paradigms work, possibly even assisting directly in design of new paradigms. By contrast, while stability properties are highly relevant for the biological brain and studying stability of RNNs provides valuable insights into network topologies of interest, provably stable systems can still demonstrate a wide variety of initial responses to transient inputs, and thus convergent behavior is not always of the utmost practical relevance. One could also imagine leveraging both of these perspectives by asking questions on which network configurations are more or less likely to encourage "resting state" hyper- or hypo- synchrony between regions in a multi-area RNN model, to further our theoretical understanding of connection patterns that might lead to psychiatric symptoms. This could then be coupled with analysis of convergent behavior of such a system in response to extended DBS-like input, in order to form new hypotheses for specific DBS system designs. By working closely with neurologists, network models could be built iteratively based on their ability to predict properties of neural activity like corticostriatal coherence.

§.§.§ Dynamical models of behavior

It is highly plausible that network configurations resulting in abnormal connectivity levels between brain regions may also require different inputs to modulate synchrony than would be effective in a more typical topology. Additionally, it is highly plausible that starkly different topology aberrations could produce observable animal behavior that is only subtly different – especially in terms of our current classification systems. There is a great deal of nuance in the extent to which a disease state is actually pathological in a particular environment as well, and although DBS is presently reserved for incredibly severe treatment resistant conditions, as technology advances and less invasive alternatives become available there will be many advantages to a closed loop system that can utilize information about system outputs and their environmental context, not just the system state information available through neural recording. A better understanding of behavior is therefore a crucial (and often underappreciated) component to the future of not only DBS (as discussed in chapter <ref>), but also neuroscientific modeling more broadly.

In the introductory chapter, I presented arguments that behavior in psychiatry disease is particularly underrated amongst computational neuroscientists. Study of aberrant system output can be a highly effective reverse engineering tool across technical disciplines, yet careful characterization of human "output" in disease states is frequently labeled as not neuroscientific, and the vast majority of computational modeling work focuses on normal learning (which incidentally, normal is not the best either). Whether one cares about disease for medical purposes or not, naturally occurring variation is a prime opportunity to drive hypothesis generation in entirely different directions than typical incremental question answering would. In my experience with non-psychiatrist neurobiology researchers, it is an open secret that "disease analogies" are mostly used disingenuously and somewhat begrudgingly for the purposes of grantsmanship. Both psychiatry and neuroscience could substantially benefit if models of disease were more often seriously engaged with for the inherent opportunities that they hold.

On the other side, psychiatry has largely failed to formalize behavioral observations in terms of dynamical systems. As pointed out by [Salvi et al., 2021], while psychiatrists often utilize information about behavioral fluctuations qualitatively when evaluating a patient, the vast majority of psychiatry research focuses on cross-sectional studies that treat disease properties as static traits. In order to build dynamical models of behavior that can be analyzed using control theory and that can be integrated with established work in computational neuroscience, it will first be necessary to collect richer and - most critically - more longitudinal behavioral data, in both humans and animal models and in both health and disease. With more precise and much more temporally dense ways of measuring system output, we can begin the iterative process of forming dynamical systems models of psychiatric symptomatology and updating them based on their ability to predict real world readouts. Collection of longitudinal naturalistic behavioral signals from psychiatric patients was indeed a theme throughout the rest of the thesis.

Note that although this framework could directly involve neural network models (and perhaps should in the case of animal models), it is also possible to define models of psychiatric symptom dynamics and employ control theory to better understand them whilst entirely abstracting away any neurological components. Such work could have immediate relevance to clinical aims, for example evaluation of treatment planning as a control system – what sort of external inputs might be able to inhibit a "bad" positive feedback loop, and what sort might trigger one? Studies often discuss correlation versus causation, but complicated causal webs are not as often considered, at least not beyond the mere mention they might exist. Dynamical models of behavior could go a long way in forming intuition about how such systems operate, bringing much needed nuance to study interpretation. Further, dynamical models of behavior - even when not explicitly neuroscientific - would provide fundamental groundwork for eventually connecting complicated neural network models to the behavior that they produce.

It is often the case in theory that a proof involves "translating" a system to some other equivalent representation in which it is much easier to prove the desired claim. Mathematics can serve as a rigorous and interdisciplinary language for describing ideas and forming downstream predictions based on them. The rest of this chapter will focus on neural network theory and deep learning science, but it is worth emphasizing that for phenomena as complex and dynamic as human (or any animal) behavior, we should really be putting a lot more focus on establishing such a language. This applies both to the neuroscientists that think of behavior as an afterthought and to the psychiatrists that focus so heavily on static diagnostic labels.

§.§ Moving forward

Neuro-inspired machine learning has profoundly altered many fields such as computer vision, natural language processing, and computational neuroscience [Bengio et al., 2017, Hassabis et al., 2017]. While deep learning models can be remarkably powerful, they are for the most part `black boxes'. This opaqueness may be dangerous in safety-critical applications like autonomous driving or human-centered robotics, and it limits conceptual progress, progress that would have the potential to in turn inform neuroscience. Such drawbacks are of particular relevance for machine learning applications in psychiatry, which encounter roadblocks with both practical clinical concerns and limitations in the foundational knowledge collectively available at present. Therefore, not only do many open questions remain to be addressed in understanding artificial neural networks, but neuroscientists ought to be triply motivated to address them.

Stability carries a number of theoretical benefits, especially for the robustness and predictability of system behavior, and is thus a crucial component in the future of deep learning science. As mentioned, stability also has implications for the biological brain, and important neuroscience concepts like sparsity go hand-in-hand with stability; it is likely that some notion of stability will be necessary for fully understanding brain function. There are many reasons to study stability in RNN models, and many open questions that remain.

While there is a rich history of stability-related results in neuroscience theory [Hopfield, 1982, Seung, 1996, Murphy and Miller, 2009, Wiener, 1948], proven stability conditions for continuous-time nonlinear RNNs remain quite restrictive – whether in the weights they allow or in other prior assumptions such as restricting the allowed nonlinearities to exclude some popular activation functions like ReLU. At the same time, incorrect RNN stability theorems have been found in the machine learning literature, and the field as a whole has struggled to prove certain seemingly simple conjectures.

After introducing the mathematical background on the RNN dynamics we study in section <ref>, I will discuss existing stability conditions for this RNN model (<ref>). Within our results in section <ref>, I will present counterexamples to a handful of problematic previously published conditions, in addition to dissecting the technical challenges that have made proving conjectured stability conditions in this model so difficult (<ref>). Suffice to say there are still a large number of open theoretical questions about continuous-time nonlinear RNN stability conditions, which is one area our work makes progress in.

Of course, perhaps the largest contribution of our work is not in the study of individual nonlinear RNN stability properties, but rather in the study of stability conditions for combinations of RNNs, as well as the application of such conditions in machine learning practice. There is a shockingly small amount of prior literature on "RNNs of RNNs" across the board, whether considering neuroscience theory or deep learning science. In the next subsection I will thus overview why modular systems are so salient and some of the major research questions about networks of networks the community should hope to address.

§.§.§ Study of network combinations: a forefront for future interdisciplinary breakthroughs

The combination and reuse of primitive "modules" has enabled a great deal of progress in computer science, engineering, and biology. Modularity is particularly apparent in the structure of the brain, as different parts are specialized for different functions [Kandel et al., 2000], and these areas must perform in concert with each other to achieve the goals of the organism. Furthermore, connectivity patterns between such regions likely play an important role in human intelligence, as models accounting for global representations have predicted cognitive ability beyond what could be explained through local representations alone [Anderson and Barbey, 2023]. Recent work in modeling macaque neural recording data has also shown that modular neural networks were much more capable of predicting  parietofrontal cortex activity during a grasping task than traditional neural networks were [Michaels et al., 2020].

Most experimental studies throughout the history of neuroscience have focused on a single brain area in association with a single behavior [Abbott and Svoboda, 2020]. Similarly, RNN models of brain function have mostly been limited to a single RNN modeling a single area. However, neuroscience is entering an age where recording from many different brain areas simultaneously during complex behaviors is possible. As experimental neuroscience has shifted towards multi-area recordings, computational techniques for analyzing, modeling, and interpreting these multi-area recordings have begun to emerge [Mashour et al., 2020, Abbott and Svoboda, 2020, Perich et al., 2021, Semedo et al., 2019, Yang and Molano-Mazón, 2021]. Despite this, RNN theory has lagged behind.

In the field of deep learning, there has been a similar paucity of study on so called "networks of networks", despite the clear neuroscience analogies and the small conceptual distance from the idea of "deep" multi-layer architectures that also take advantage of compositionality. There have been successful architectures with hints of modularity – such as the Generative Adversarial Network (GAN), which typically utilizes a two network structure where one network is rewarded for "tricking" the other [Goodfellow et al., 2016]. The groundbreaking success of AlphaGo was similarly enabled by a multi-network structure [Silver et al., 2016]. However, there has been minimal work to date on leveraging modularity in highly recurrent systems, as well as minimal mathematical or scientific work on how modularity impacts network learning.

At the same time, there is increasing demand for machine learning systems that can accomplish multimodal tasks, as well as for deep learning agents that can multi-task, continually learn, accurately assess uncertainty, behave in an interpretable way, etc. Artificial neural networks have lead to paradigm shifting advances in a number of domain-specific areas like image recognition and language processing, but applications that require "putting it all together" have not been quite so earth shattering to date. In the light of the many open questions that remain about modular deep learning systems, this is not surprising. For example, while one might intuitively expect that training a particular module to perform better (e.g increasing visual acuity) would improve the performance of the overall system, this is not necessarily the case. In the context of self-driving cars, it has been shown recently that additional training for one module while holding the other modules' parameters fixed can actually make overall performance worse [Wu et al., 2021].

Establishing sufficient conditions for preservation of functionality in the construction and training of modular systems is therefore an interesting and open question. Our research makes several advances in this domain by providing conditions under which training one (or several) modules while holding the others fixed provably preserves the stability of the overall system – a minimal requirement for robust and predictable behavior. To do so, we utilize a special notion of nonlinear stability called contractive stability.

A conceptual introduction to contractive stability.
In section <ref> below, I will provide a mathematical primer on contraction analysis, along with additional context from nonlinear stability theory. Here I focus on contraction analysis from the perspective of how it enabled our work, as well as its ties to biological sciences. Most notably contractive stability has potential applications in evolutionary biology and - as we show throughout this work - neurobiology. However, contraction analysis can also be used in many other domains, such as in modeling the dynamics of bacterial quorum sensing or the spread of infectious disease.

A major benefit of contractive stability is that there are a number of ways in which contracting systems can be combined whilst preserving stability of the overall system. Indeed, themes of modularity are apparent in biological evolution and in the structure of the brain. It is thought that the majority of traits that have developed over the last 400+ million years are the result of evolutionary forces acting on regulatory elements that combine core components, rather than mutations to the core components themselves. This mechanism of action makes meaningful variation in population phenotypes much more feasible to achieve, and is appropriately titled "facilitated variation" [Gerhart and Kirschner, 2007]. In addition to the biological evidence for facilitated variation, computational models have demonstrated that this approach produces populations that are better able to generalize to new environments [Parter et al., 2008], an ability that will be critical to further develop in deep learning systems.

[Combinations of stable systems are not necessarily stable.]Combinations of stable systems are not necessarily stable. When simulating neural dynamics in MATLAB, it can be easily seen that maintaining stability in the final network when combining two stable networks requires great care. Here we see how each of the neuron voltages varies over an input to two different stable RNNs (top), with values staying in check. When the two (black and green) RNNs are connected together, we can see that random conjoining connections will cause the voltages in the new combined RNN to rapidly rise in response to the same input and ultimately "explode" (bottom left), whereas using negative feedback conjoining connections derived from the contraction combination properties will keep the new combined RNN's behavior stable (bottom right).

Although the importance of building modular systems is clear, in general there is no guarantee that a combination of stable systems will itself be stable (Figure <ref>). Thus the tractability of these proposed evolutionary processes hinges on some mechanism for ensuring stability of combinations. As contraction analysis tools allow complicated contracting systems to be built up recursively from simpler elements, this form of stability is well suited for neural and other biological systems [Slotine and Liu, 2012]. Moreover, contracting combinations can be made between systems with very different dynamics, as long as those dynamics are each contracting – thereby facilitating a great deal of flexibility in constructing contracting modular systems [Slotine and Lohmiller, 2001]. In our particular case, this also means that the contracting RNNs discussed herein could be combined with physical systems with quantified contraction properties.

§.§.§ Chapter aims

As discussed, cognitive models are moving increasingly towards the study of multi-area dynamics, but many questions remain on how to best train and evaluate such networks [Yang et al., 2018, Yang and Molano-Mazon, 2021]. Understanding how different brain regions interact harmoniously to produce a unified percept/action will require new ideas and analysis tools.

In the presented work, we leveraged the hierarchical and especially negative feedback combination properties of contracting systems in order to develop a framework for constructing complex nonlinear "RNNs of RNNs" that preserve stability features of their component subnetworks, something not generally guaranteed with system combinations (Figure <ref>). In this process, we also pushed forward the literature on stability conditions for the individual nonlinear RNN dynamics in continuous time, and subsequently provided empirical proof of concept results for a deep learning implementation of our architecture design on challenging sequence learning tasks.

Ultimately, our primary contributions are four-fold:

* Novel contraction conditions for continuous-time nonlinear RNNs, to use in conjunction with the combination condition (section <ref>).

* A novel parameterization for linear negative feedback combination of contracting nonlinear RNNs that enables direct optimization using standard deep learning libraries (section <ref>).

* Identification of flaws in stability proofs from prior literature on nonlinear continuous-time RNNs, along with technical discussion of some of the challenges inherent in proving popular conjectured conditions for RNN stability (section <ref>).

* Experiments demonstrating that our 'network of networks', implemented based on our theoretical results, sets a new state of the art (SOTA) for stability-constrained RNNs on benchmark sequential processing tasks – also performing well against overall SOTA benchmarks with a much smaller number of trainable parameters (section <ref>).

As part of my thesis chapter, I additionally provide extended commentary on the numerous future directions that our work can enable (section <ref>).

§ MATHEMATICAL PRELIMINARIES

Before presenting our mathematical results (section <ref>), I will first establish the necessary mathematical background to understand them. Recall that we focus on proving the stability of a given fixed RNN (or RNN of RNNs), not on the stability of the weight update process. Though at times we will discuss how to guarantee that a network remains provably stable throughout the training process per some of our conditions, the primary perspective remains on stability certificates for networks with static connections. The individual nonlinear RNNs we study, characterized mainly by their weight matrix $\mathbf{W}$, will be defined along with further model explanation in section <ref>. I will then introduce details of the specific nonlinear control tool we use - contraction analysis - and present background on the combination properties enjoyed by contracting systems (<ref>), which we leverage in conjunction with the individual RNN model to create our modular "RNNs of RNNs". Finally, section <ref> overviews prior literature relevant to the chosen RNN model and to applications of contraction analysis in neuroscience, in greater technical specifics enabled by the rest of the mathematical preliminaries.

§.§ Nonlinear model definitions

This chapter focuses on analysis of rate-based neural networks. Unlike spiking neural networks, these models are continuous and smooth. We consider the following RNN with update dynamics introduced by [Wilson and Cowan, 1972], which may be viewed as an approximation to a more biophysically-detailed spiking network:

\begin{equation}\label{eq:RNN}\tag{nonlinear RNN dynamics}
\begin{split}
\tau \dot{\mathbf{x}} = f(\mathbf{x}) = -\mathbf{x} + \mathbf{W}\phi(\mathbf{x}) + \mathbf{u}(t)
\end{split}
\end{equation}

Here $\tau > 0$ is the time constant of the network [Dayan and Abbott, 2005], and the vector $\mathbf{x} \in \mathbb{R}^n$ contains the "voltages" of all $n$ neurons in the network. The voltages are converted into firing-rates through a static nonlinearity $\phi$. We only consider monotonic activation functions with bounded slope: in other words, $0 \leq \phi' \leq g$. We do not restrain the sign of the nonlinearity. The matrix $\mathbf{W} \in \mathbb{R}^{n \times n}$ contains the synaptic weights of the RNN. It is this matrix that ultimately determines the stability of (<ref>), and will be a main target for our analysis. Finally, $\mathbf{u}(t)$ is the potentially time-varying external input into the network, capturing both explicit input into the RNN from the external world, as well as unmodeled dynamics from other brain areas.

Common example nonlinearities $\phi(x)$ that satisfy the given constraints include hyperbolic tangent (tanh) and rectified linear unit (ReLU), both of which have $g=1$. Though most of our theorems can handle more general $g$ (as will be included within the proof details), I will assume $g=1$ in much of the discussion in this chapter, for more simple interpretability. Note that some prior literature on the (<ref>) network instead restricts $\phi' > 0$, for example the stability analyses by [Matsuoka, 1992]. Such results will not necessarily apply to networks using activation functions like ReLU, which has 0 slope for negative input values.

§.§.§ Network Jacobian

Contraction analysis, to be reviewed next (section <ref>), requires the Jacobian $\mathbf{J}$ of the system in question. The Jacobian is the matrix of partial derivatives of the dynamics with respect to the system state – in a 2 unit system, we could express (<ref>) in terms of $f(x_{1})$ and $f(x_{2})$ separately, and then see the Jacobian would be $\mathbf{J} =$

\begin{align*}
\begin{bmatrix}
\frac{\partial f(x_{1})}{\partial x_{1}} & \frac{\partial f(x_{1})}{\partial x_{2}} \\
\frac{\partial f(x_{2})}{\partial x_{1}} & \frac{\partial f(x_{2})}{\partial x_{2}}
\end{bmatrix}
\end{align*}

We can see that for our system (<ref>), the leak term $-\mathbf{x}$ will contribute only diagonal entries to the Jacobian, as each unit's leak only directly affects itself. Similarly, the input $\mathbf{u}(t)$ is not state-dependent and thus will have 0 impact on the Jacobian. However the connection weights $\mathbf{W}$ allow different units to influence each other, so that this term accounts for most of $\mathbf{J}$ here. Note that when we take $\phi(x) = x$ i.e. a linear system, the partial derivative matrix of the term would be $\mathbf{W}$. The nonlinearity of course complicates things.

Ultimately, we find through the chain rule that the system (<ref>) has Jacobian:

\begin{equation*}
\begin{split}
\mathbf{J} = \mathbf{WD} - \mathbf{I}
\end{split}
\end{equation*}

where $\mathbf{D}$ is a diagonal matrix with entry $D_{ii} = \phi'(x_{i})$. $\mathbf{D}$ thus represents the derivative of the activation function evaluated at the current state. Note that we can allow different activation functions for different units without loss of generality, as long as each function satisfies our constraint $0 \leq \phi' \leq g$. $\mathbf{D}$ is therefore a matrix that can vary over time due to inputs or internal dynamics, but it will always be nonnegative and each entry will always be upper bounded by $g$.

Equivalent representations.
(<ref>) is equivalent to another commonly used class of RNNs where the term $\mathbf{W}\mathbf{x} + \mathbf{u}$ appears inside the nonlinearity $\phi$. See the proof by [Miller and Fumarola, 2012] for details. Our mathematical results therefore apply equally well to both types of RNNs. Because this other class of RNNs has Jacobian $\mathbf{DW} - \mathbf{I}$, we can use either ordering when plugging $\mathbf{J}$ into contraction analysis, and any resulting conditions will apply to both dynamics formats.

§.§ Contraction analysis

Unlike linear stability analysis, there are many different notions of stability for a nonlinear system, and it can sometimes be difficult to say anything at all about the stability of a particular nonlinear dynamics. Lyapunov analysis for example relies on finding an appropriate Lyapunov function to prove that a system is Lyapunov stable, but failing to find a Lyapunov function does not prove that a general system is not Lyapunov stable [Slotine and Li, 1991]. As discussed within the background section <ref>, another complication of nonlinear stability analysis is that combinations of stable systems need not be stable in general, and even when combined in a seemingly "stable" way, when nonlinearities are involved it can be extremely difficult to verify that stability was actually maintained.

Here we focus on a special type of nonlinear stability, known as contractive stability, as introduced by [Lohmiller and Slotine, 1998]. Contractive stability is a strong form of dynamical stability which implies many other forms of stability, such as input-to-state stability [Sontag, 2010]. Loosely speaking, a dynamical system is said to be globally contracting if any two of its trajectories converge to each other exponentially, regardless of initial conditions (Figure <ref>). Contraction analysis is convenient because it allows us to analyze convergence behavior without worrying about the identity of equilibrium points like in e.g. Lyapunov analysis. This also enables the combination properties of contractive stability, which we make heavy use of throughout the present work.

[Contraction theory is based on a differential analysis of convergence.]Contraction theory is based on a differential analysis of convergence. By focusing on exponential convergence of any two neighboring trajectories, contraction analysis allows convergence and limit behavior to be treated separately. While the traditional Lyapunov analysis for stability of nonlinear systems can be seen as "virtual classical mechanics", contraction analysis can be seen as "virtual fluid mechanics". This explanatory cartoon has been reproduced from [Lohmiller and Slotine, 1998], the paper introducing contraction. It is fair use, and the senior author is of course also the senior author on our paper [Kozachkov et al., 2022]

A primary advantage of contraction analysis is that it is directly applicable to non-autonomous systems, which the vast majority of recurrent models are [Slotine and Lohmiller, 2001, Slotine, 2003]. Contraction analysis has found wide application in nonlinear control theory [Manchester and Slotine, 2017], synchronization [Pham and Slotine, 2007], and robotics [Chung and Slotine, 2009], but has only recently begun to find application in neuroscience and machine learning [Boffi et al., 2020, Wensing and Slotine, 2020, Kozachkov et al., 2020, Revay and Manchester, 2020, Jafarpour et al., 2021, Kozachkov et al., 2022, Centorrino et al., 2022, Burghi et al., 2022, Kozachkov et al., 2022]. Contraction analysis is particularly useful for neuroscience because it is directly applicable to systems with external inputs. It also has the potential to be closely linked with evolutionary biology through its combination properties [Gerhart and Kirschner, 2007, Slotine and Liu, 2012].

I will next provide a summary of the math behind contraction analysis as originally reported by [Lohmiller and Slotine, 1998], along with the specific contraction condition that will be explored in this chapter for our (<ref>) network (subsection <ref>). I will then review two of the core contraction combination properties - hierarchical combos and negative feedback combos - in terms of technical details, to set up how we will later make our 'networks of networks' (subsection <ref>).

§.§.§ Contraction math

We will begin by considering contractive stability in Euclidean space. The core idea behind contraction analysis is that stability can be analyzed differentially by evaluating if nearby trajectories will converge to each other – so we consider an infinitesimal displacement at a fixed time and then "observe" the subsequent behavior. Specifically, we will measure the squared distance between the displaced and original trajectories over time, and aim to demonstrate exponential convergence (Figure <ref>). If they do exponentially converge, and we can say this for similar displacements of any trajectory of the system, then the system is globally contracting in Euclidean space.

[Lohmiller and Slotine, 1998] use this to derive the condition that if $\mathbf{J} + \mathbf{J^{T}} \preceq -\beta$ uniformly, then the system with Jacobian $\mathbf{J}$ is globally contracting at rate $\beta$ ($> 0$). Note that throughout we take $\prec$ to correspond to an inequality on the maximum eigenvalue of the matrix on the left hand side. Recall also that the symmetric part of a matrix $\mathbf{A} = \frac{\mathbf{A} + \mathbf{A^T}}{2}$, and all matrices can be broken down into the sum of their symmetric and skew-symmetric parts. Symmetric matrices have all real eigenvalues, and a matrix is referred to as 'positive definite' if it is symmetric with only positive eigenvalues.

The core idea of contraction analysis does not require such a specific distance metric as Euclidean space however. Contractive stability can involve a wide variety of distance functions, and these distance functions can even vary over the course of the "observation". It is of course necessary that the distance metric always evaluates to 0 when given identical points, and evaluates to a positive number when not given identical points. Basically, contraction is guaranteeing that even in Euclidean space the trajectories will converge, but it is possible the path will look squiggly, while in some other valid space the trajectories are in fact exponentially converging to each other with no "backtrack" possible – like what we saw in Euclidean space for the original simple case.

Specifically, [Lohmiller and Slotine, 1998] show that the non-autonomous system

\[\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x},t) \]

is contracting if there exists a metric $\mathbf{M}(\mathbf{x},t)  = \mathbf{\Theta}(\mathbf{x},t)^T\mathbf{\Theta}(\mathbf{x},t) \succ 0$ such that uniformly

\[\dot{\mathbf{M}} + \mathbf{M}\mathbf{J} + \mathbf{J}^T\mathbf{M} \preceq -\beta \mathbf{M} \]

where $\mathbf{J} = \frac{\partial \mathbf{f}}{\partial \mathbf{x}}$ and $\beta > 0$.

Essentially, we want to find some uniformly positive definite $\mathbf{M}$ that will make our system Jacobian $\mathbf{J}$ uniformly satisfy $\dot{\mathbf{M}} + \mathbf{M}\mathbf{J} + \mathbf{J}^T\mathbf{M} \prec 0$ in order to prove global contraction. This chapter focuses entirely on using a constant metric to show global contraction of our network (for reasons that will be discussed in section <ref>), so we deal exclusively here with a static $\mathbf{M}$ and thus $\dot{\mathbf{M}} = 0$.

This leads us to the following contraction condition for the (<ref>) system, which we will work to prove is true for particular restrictions on $\mathbf{W}$ and with particular valid metrics $\mathbf{M}$:

\begin{equation}\label{eq:JacobCond}
\begin{split}
\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} \prec 0
\end{split}
\end{equation}

In other words, we want to show that there exists a positive definite $\mathbf{M}$ such that $\mathbf{MWD}_{sym} - \mathbf{M}$ has all negative eigenvalues for a given network, for all possible values that $\mathbf{D}$ can take. As $\mathbf{MWD}_{sym} - \mathbf{M}$ is a symmetric matrix, an equivalent statement is that if we can show this is negative definite uniformly, we have proven that (<ref>) with weights $\mathbf{W}$ is globally contracting in metric $\mathbf{M}$.

§.§.§ Modularity principles

As mentioned, a number of combination properties for contracting systems have been previously proven [Slotine and Lohmiller, 2001, Slotine, 2003]. These properties allow for different types of connections between different stable systems to be formed while preserving contractivity in the combined system. It is possible to combine systems with different underlying dynamics, systems that are contracting in different metrics, etc. through these properties. Here we focus on negative feedback and hierarchical system combinations, which have straightforward neuroscientific interpretations and which can allow for highly complex combined networks to be built up through recursive construction (Figure <ref>).

[Contractive stability implies a modularity principle.]Contractive stability implies a modularity principle. Because contraction analysis tools allow complicated contracting systems to be built up recursively from simpler elements, this form of stability is well suited for understanding biological systems. Contracting combinations can be made between systems with very different dynamics, as long as those dynamics are contracting.

A major goal of our results (section <ref>) was to take individual RNNs described by (<ref>) that are provably contracting, and then apply these combination properties to build up "RNNs of RNNs" that are themselves contracting. I will therefore introduce the prior more general results on contracting combinations of relevance here.

Feedback combination.
Consider two systems, independently contracting in constant metrics $\mathbf{M}_1$ and $\mathbf{M}_2$, which are combined in feedback:

\[\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x},t) + \mathbf{B}\mathbf{y} \]

\[\dot{\mathbf{y}} = \mathbf{g}(\mathbf{y},t) + \mathbf{G}\mathbf{x} \]

\begin{equation*}
\mathbf{B} = -\mathbf{M}^{-1}_1\mathbf{G}^T\mathbf{M}_2
\end{equation*}

then the combined system is contracting as well. This may be seen as a special case of the feedback combination derived by [Tabareau and Slotine, 2006]. In particular, this describes a linear negative feedback connection between the systems.

Hierarchical combination.
Consider again two systems, independently contracting in some metrics, which are combined in hierarchy:

\[\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x},t)\]

\[\dot{\mathbf{y}} = \mathbf{g}(\mathbf{y},t) + \mathbf{h}(\mathbf{x},t) \]

where $\mathbf{h}(\mathbf{x},t)$ is a function with bounded Jacobian. Then this combined system is contracting in a diagonal metric, as shown by [Lohmiller and Slotine, 1998]. By recursion, this extends to hierarchies of arbitrary depth. Note also this permits for nonlinear connection between the systems.

§.§ Prior research on RNN model stability

While there has been minimal theory to date on stability conditions for multi-area nonlinear RNNs, there of course has been a good deal of prior literature on the individual RNN model described by (<ref>), as well as prior application of contraction analysis to other recurrent models. This section serves to establish more technical background on such prior results, to better contextualize our work on (<ref>).

We first note that the 'Echo-State Condition' introduced by [Jaeger, 2001] is equivalent to discrete-time contraction in the identity metric. A later generalization of this condition included a diagonal metric [Buehner and Young, 2006]. In the context of neuroscience, contraction analysis has been applied to analyzing the dynamics of winner-take-all networks [Rutishauser et al., 2011, Rutishauser et al., 2015] as well as networks with synaptic plasticity [Kozachkov et al., 2020]. In the machine learning context, Miller and Hardt recently rederived an `echo-state property' for discrete recurrent models, and went on to prove that these contracting recurrent models could be well-approximated by feedforward networks in certain cases [Miller and Hardt, 2018].
More recently still, in a series of papers Revay, Wang, and Manchester applied contraction analysis to discrete-time recurrent networks [Revay and Manchester, 2020, Revay et al., 2021, Revay et al., 2020], expanding the class of models considered by [Miller and Hardt, 2018].

In addition to contraction (which amounts to a strong form of non-autonomous exponential stability) there has been a considerable amount of work attempting to enforce weaker forms of stability in RNNs, such as autonomous stability [Erichson et al., 2021, Chang et al., 2019]. Unfortunately, as we discuss in section <ref>, autonomous stability does not in general imply non-autonomous stability, so it is not clear what stability properties these RNNs possess when driven with external input.

There is also a line of work which uses orthogonal weight matrices to avoid the vanishing/exploding gradient problem during training. Orthogonality is typically ensured during training via a parameterization [Arjovsky et al., 2016, Lezcano-Casado and Martınez-Rubio, 2019] – for example, by exploiting the fact that the matrix exponential of a skew-symmetric matrix is orthogonal, as is done in [Lezcano-Casado and Martınez-Rubio, 2019]. These works focus on parameterizations of individual RNNs which ensure stability is preserved throughout training, and our use of Theorem <ref> does draw from these works. By contrast, while we also introduce novel stability conditions for individual RNNs, our work focuses on 'combined' RNNs (i.e network of networks) and provides a parameterization on the connections between stable subnetworks such that training preserves the overall network stability (Figure <ref>).

§ MATHEMATICAL RESULTS

The mathematical results of this chapter fall under 4 major categories, to be covered in the following subsections:

* New stability conditions on the weights of the general network equation (<ref>), each guaranteeing global contraction in a constant metric when satisfied (subsection <ref>).

* Application of prior contraction analysis results to develop a deep learning applicable methodology for combining contracting networks described by (<ref>) into a larger recurrent "RNN of RNNs" that is itself provably stable (subsection <ref>).

* A discussion of conditions under which pruning might disrupt stability, to facilitate future work on directly sparsifying multi-area RNNs, a topic of potential relevance both for more organic extensions on our pilot experimental results to be discussed, as well as for better understanding neuroscientific principles (subsection <ref>).

* Counterexamples to a few previously published stability "conditions" for (<ref>), alongside a deeper explanation for why certain conjectured stability conditions on these dynamics have been so hard to prove thus far (subsection <ref>).

§.§ Novel contraction conditions for individual, continuous-time RNNs

This subsection focuses on presenting, explaining, and proving our novel stability conditions for the nonlinear RNN described by (<ref>). As previously proven conditions have been quite restrictive, these theoretical results are valuable in their own right. However, they are also important to our broader aims, because the new conditions can guarantee that (<ref>) is globally contracting in a constant metric, and moreover can provide such a metric – the properties that allow us to apply prior results in contraction analysis to construct stable "networks of networks" out of the component RNNs satisfying the conditions here. The two major contraction conditions on (<ref>) that we derive have quite different pros and cons, and were employed in different ways experimentally.

As will be covered in section <ref>, Theorem <ref> was used for the "Sparse Combo Net" stable multi-area architecture, which demonstrated the best performance in our hands. The condition is very simple to check because of its deep link to linear stability, and it is satisfied by a number of matrices with interesting sparse structures. Theorem <ref> will thus be the major focus of exposition in this subsection, covered in depth in <ref> next. The other condition I will review here (in <ref>) was used for our "SVD Combo Net" architecture. Theorem <ref> has the advantage that matrices satisfying it can be written in terms of straightforward operations between other matrices, the latter of which are all either expressions that do not require special constraints or contraction properties that can be preset. The condition can therefore be enforced throughout training of a network by performing optimization on components of said parameterization. By contrast, we could find no clear way to directly parameterize Theorem <ref>.

§.§.§ Linear stability of weight magnitudes as a nonlinear contraction condition

As will be discussed in more depth in section <ref>, interpretation of negative feedback is the major complication in proving stability of a nonlinear system described by (<ref>). We can think of this intuitively using the rectified linear unit (ReLU) activation function as an example. ReLU is a valid $\phi$ for this system and it is a very common activation function to use in modern deep learning systems. From a neuroscience perspective, ReLU was originally designed to move on from the classic sigmoid function interpretation of unit outputs as an "all or nothing" signal, instead framing them to be interpreted as neuron firing rates. As such, ReLU outputs 0 for inputs $\leq 0$, but otherwise just outputs the inputs, making the derivative instead the binary function.

When we compare ReLU activation to an otherwise identical linear network, it becomes obvious that if everything were constrained to be positive, the networks would be functionally equivalent. In reality we do not assume control of initial conditions or of inputs, but if the network can "survive" the most adversarially large all-positive inputs, it makes intuitive sense that negative inputs would not "break" global contraction. Linear stability is already well-understood and easy to check, so taken together this is what inspired the question of whether we could find stability conditions on $\mathbf{W}$ in terms of the linearized version of the system: first beginning with the narrow class of excitatory-only RNNs, but then extending the idea to a more general condition on magnitudes of weights in the network.

We ultimately ended up with the following theorem, which shows that a linearly stable $\mathbf{|W|} - \mathbf{I}$ implies the original nonlinear system with weight matrix $\mathbf{W}$ is also stable, and in fact will be stable in the same constant metric, a property we will leverage later.

Let $\mathbf{|W|}$ denote the matrix formed by taking the element-wise absolute value of $\mathbf{|W|}$.  If there exists a positive, diagonal $\mathbf{M}$ such that:

\[\mathbf{M}(\mathbf{|W|}-\mathbf{I}) + (\mathbf{|W|}-\mathbf{I})^T\mathbf{M} \prec 0 \]

then (<ref>) is contracting in metric $\mathbf{M}$. If $W_{ii} \leq 0$, then $|W|_{ii}$ may be set to zero to reduce conservatism.

This condition is particularly straightforward in the common special case where the network does not have any self weights, with the leak term driving stability. While it can be applied to a more general $\mathbf{W}$, the condition will of course not be met if the network was relying on highly negative values on the diagonal of $\mathbf{W}$ for linear stability. As demonstrated by counterexample in the proof of Theorem <ref>, it can be impossible to use the same metric $\mathbf{M}$ for the nonlinear RNN in such cases, which breaks the proof method. However as the issue arises when only some of the self weights are highly negative, one could consider instead allowing a gain on the leak term (which would scale the identity matrix $\mathbf{I}$) if designing a system.

Theorem <ref> allows many weight matrices with low magnitudes or a generally sparse structure to be verified as contracting in the nonlinear system (<ref>), by simply checking a linear stability condition.  This is particularly useful because sparsity is a recurring theme in neuroscience, both enabling stability in theory [Kozachkov et al., 2020], and occurring in practice, with only $\sim 0.00001\%$ of potential connections forming synapses in the human brain [Slotine and Liu, 2012]. This could explain why there has been little issue approximating nonlinear systems with linear ones in many computational neuroscience papers.

Beyond verifying contraction, Theorem <ref> actually provides a metric, with little need for additional computation. Not only is it of inherent interest that the same metric can be shared across systems in this case, it is also of use in machine learning, where stability certificates are becoming increasingly necessary for some applications. More generally, there are a variety of systems of interest that meet this condition but do not meet the well-known maximum singular value condition, including those with a hierarchical structure.

One downside of the condition is that it is not straightforward to enforce throughout training in a deep learning framework. However, there are a number of ways that one could encourage the condition to be maintained throughout training, and it is easy to confirm the final network is stable under these criteria before using it in applications. A number of regularization techniques that work to keep network weights small and/or sparse would likely help with maintaining the condition, as would strategic initialization methods or training paradigms that restrict which weights can be updated. This is even more relevant in light of the recent impacts of the Lottery Ticket Hypothesis, which states that given a neural network with a particular initialization, aggressive pruning after training can identify much smaller networks that are similarly capable of learning the task from scratch [Frankle and Carbin, 2019]. One could imagine using pruning to enforce Theorem <ref> post-hoc, or in an inverse scenario, aiming to build up small stable RNNs capable of strong task performance via mechanisms like the Weight Agnostic Neural Networks introduced by [Gaier and Ha, 2019].

Of course in our multi-area framework (to be detailed in section <ref>), it is possible to train the linear negative feedback connections between nonlinear RNNs meeting the Theorem <ref> condition without updating the nonlinear RNNs themselves, like a more complex version of the historical echo state network idea [Jaeger, 2001] and reminiscent of the evolutionary concept of facilitated variation [Gerhart and Kirschner, 2007]. This is the approach we took with our pilot experiments, which demonstrated great early success particularly in the context of the simplicity of this system (<ref>).

I will now restate the Theorem <ref> in different but equivalent terms, before proceeding to provide a proof of the claims. The primary principle underlying the proof is the fact that linear stability is equivalent to diagonal stability for Metzler (i.e. nonnegative off-diagonal) matrices [Narendra and Shorten, 2010]. This principle is used in conjunction with well-established linear algebra inequalities. To skip these details, please see instead the next Theorem's introduction below, in subsection <ref>.

Proof of Theorem <ref>.
Stated another way:

Given an anti-diagonal (i.e. no self weights) $\mathbf{W}$ in the nonlinear RNN, if its elementwise absolute value $\mathbf{|W|}$ would produce a linearly stable system (i.e. $\mathbf{|W|} - \mathbf{I}$ is Hurwitz), then the nonlinear system with $\mathbf{W}$ is guaranteed to be globally contracting. Furthermore, it is globally contracting in the same diagonal metric(s) in which the $\mathbf{|W|}$ system is linearly stable.

Recall that if a positive definite $\mathbf{M}$ exists such that $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} \prec 0$ uniformly for a given $\mathbf{W}$, then the RNN system is contracting in that metric $\mathbf{M}$. As noted, '$\prec$' refers throughout this chapter to a bound on the maximum eigenvalue of the matrix on the left hand side. Recall also that $\mathbf{D}$ is a diagonal matrix containing the activation function derivatives evaluated at the current state for each unit (which must be nonnegative and at most 1 for the functions considered here). See section <ref> for further introduction on mathematical notation, contraction analysis, and the Jacobian of the system in question.

For proof of the main claim, we will first prove that given a nonnegative (i.e. all elements nonnegative) Hurwitz $\mathbf{W}$, the system is contracting. Note that if a nonnegative $\mathbf{W}$ is linearly stable at all, it must be diagonally stable [Narendra and Shorten, 2010] – it then follows that the system with linear activation (meaning $\mathbf{D}=\mathbf{I}$ for all time) must be contracting in some diagonal metric. So, there exists a positive definite diagonal $\mathbf{M}$ such that $\frac{\mathbf{MW} + \mathbf{W^{T}M}}{2} - \mathbf{M} \prec 0$. To extend to the nonlinear activation case, We will show that given any diagonal $\mathbf{D}$ with all elements between 0 and 1 (inclusive), we can guarantee that $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} \prec 0$, meaning the nonlinear system is contracting in the same diagonal metric $\mathbf{M}$. Note that the proof applies equally well regardless of the $\mathbf{M}$ used, as long as it is a diagonal contraction metric for the linear system; therefore the nonlinear system will share all diagonal contraction metrics with the linear counterpart.

To prove the claim for the nonnegative Hurwitz $\mathbf{W}$, we will make use of the spectral radius function, $\rho$. Per Perron-Frobenius, for a real nonnegative symmetric matrix $\mathbf{A}$, $\rho(\mathbf{A})$ will always be equivalent to $\lambda_{max}(\mathbf{A})$ (i.e. the maximum eigenvalue of $\mathbf{A}$), which will be real. Furthermore, the spectral radius is monotone for nonnegative matrices $\mathbf{A}$ and $\mathbf{B}$, meaning $\rho(\mathbf{A}) \leq \rho(\mathbf{B})$ if $\mathbf{A} \leq \mathbf{B}$ elementwise. More broadly, we will use the fact that $\rho(\mathbf{A}) \leq \rho(\mathbf{|A|}) \leq \rho(\mathbf{B})$ if $\mathbf{|A|} \leq \mathbf{B}$ elementwise at multiple times throughout this proof. A proof of that inequality can be found in the Linear Algebra text by [Meyer, 2000], in Example 7.10.2.

Now observe that $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M}$ is the sum of two symmetric matrices, and therefore symmetric. It is Metzler, but the diagonals can (and will) contain negative values due to the $-\mathbf{M}$ contributed by the leak term (otherwise it would not be linearly stable in the first place). To work with a fully nonnegative matrix instead, we will choose some positive integer $c$ such that $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I}$ has all nonnegative diagonals. Given the system assumptions, it will always be possible to use $c = max(\mathbf{M})$ i.e. the maximum element in the diagonal contraction metric used for the linear system. Adding $c\mathbf{I}$ will increase the eigenvalue of the resulting matrix by exactly $c$ [Meyer, 2000], so that without loss of generality we can instead prove $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I} \prec c$ uniformly to show that the nonlinear RNN system is contracting in metric $\mathbf{M}$.

$\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I}$ meets all our desired conditions for utilizing the spectral radius. Furthermore, we know that $\frac{\mathbf{MW} + \mathbf{W^{T}M}}{2} - \mathbf{M} + c\mathbf{I} \prec c$ by the same logic, and that $\frac{\mathbf{MW} + \mathbf{W^{T}M}}{2} - \mathbf{M} + c\mathbf{I}$ is also a real symmetric nonnegative matrix. We will therefore aim to show that $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I} \preceq \frac{\mathbf{MW} + \mathbf{W^{T}M}}{2} - \mathbf{M} + c\mathbf{I}$ uniformly via the spectral radius inequality, thereby demonstrating that $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I} \prec c$, as desired.

To do so, note that for any nonnegative $\mathbf{W}$, the diagonal $\mathbf{D}$ with values between 0 and 1 can only possibly decrease the elements of $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2}$ (recall $\mathbf{M}$ in this case is a positive diagonal matrix). Because the other terms on the left hand side of the above inequality are constant regardless of $\mathbf{D}$, $\mathbf{D}=\mathbf{I}$ is guaranteed to produce the largest possible matrix elementwise. Due to the monotonicity of spectral radius in nonnegative matrices, the max eigenvalue of $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I}$ must then always be $\leq$ to that when $\mathbf{D}=\mathbf{I}$ – which is exactly the linear RNN condition given. Thus the maximum eigenvalue of $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I}$ can never go above $c$ regardless of how $\mathbf{D}$ varies, completing our lemma for Hurwitz nonnegative networks.

To return to the main claim, suppose instead the given was that $\mathbf{|W|}$ is Hurwitz, where $\mathbf{|W|}$ represents the elementwise absolute value of the weight matrix $\mathbf{W}$. Analogous to the simpler case, we know based on this information that there is a diagonal $\mathbf{M}$ such that $\frac{\mathbf{M|W|} + \mathbf{|W|^{T}M}}{2} - \mathbf{M} \prec 0$. If we plug in the preceding spectral radius inequality here, we have that $\rho(\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I}) \leq \rho(|\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I}|)$. Our aim is to use these known properties to show that the original $\mathbf{W}$ is contracting in the nonlinear system with metric $\mathbf{M}$, i.e. $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} \prec 0$ uniformly.

We will begin with an assumption that $\mathbf{W}$ has nonnegative diagonals. In that case, we can pull $c\mathbf{I} - \mathbf{M}$ (which will be a positive diagonal matrix) outside of the absolute value sign. Then observe that taking the absolute value after taking the symmetric part of a matrix will always produce a result that is elementwise less than or equal to the matrix that results from taking the absolute value beforehand (because anti-symmetric elements cannot cancel out when the absolute value is taken first). $\mathbf{M}$ and $\mathbf{D}$ are nonnegative diagonal matrices, so that $\mathbf{M|W|D} + \mathbf{D|W|^{T}M} = \mathbf{|MWD|} + \mathbf{|DW^{T}M|}$. This gives us the final relationship $|\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I}| \leq \frac{\mathbf{M|W|D} + \mathbf{D|W|^{T}M}}{2} - \mathbf{M} + c\mathbf{I}$ elementwise.

These matrices can now be plugged into the above spectral radius inequality to get $\rho(\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I}) \leq \rho(|\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} + c\mathbf{I}|) \leq \rho(\frac{\mathbf{M|W|D} + \mathbf{D|W|^{T}M}}{2} - \mathbf{M} + c\mathbf{I})$. Thus the maximum magnitude eigenvalue of the original system will always be less than or equal to that of the system with $\mathbf{|W|}$, where we can already show stability across possible $\mathbf{D}$ per our first lemma.

The above assumed that $\mathbf{W}$ had no negative diagonal entries. To still apply the technique even with a $\mathbf{W}$ that has negative self weights, we can simply 0 those out – as this only makes the system less likely to be stable. More formally, we will express $\mathbf{W}$ as $\mathbf{C} + \mathbf{S}$, where $\mathbf{S}$ is a diagonal matrix containing all negative self weights (0 on the corresponding diagonals that do not have a negative value), and $\mathbf{C} = \mathbf{W} - \mathbf{S}$ is all the connections and any remaining positive self weights (note these must be $< 1$ for hope at stability anyway).

Plugging the new expression for $\mathbf{W}$ into the symmetric part of the Jacobian and rearranging terms, we have $\mathbf{J} = (\mathbf{M}(\mathbf{C} + \mathbf{S})\mathbf{D})_{sym} - \mathbf{M} = (\mathbf{MCD})_{sym} + (\mathbf{MSD})_{sym} - \mathbf{M}$. Applying Weyl's inequality, $\mu_{2}(\mathbf{MCD} + \mathbf{MSD} - \mathbf{M}) \leq \mu_{2}(\mathbf{MCD} - \mathbf{M}) + \mu_{2}(\mathbf{MSD})$. With $\mathbf{M}$ and $\mathbf{D}$ as nonnegative diagonal matrices and $\mathbf{S}$ a nonpositive one, $\mu_{2}(\mathbf{MSD}) \leq 0$, so that we have the final inequality $\lambda_{max}(\mathbf{J}_{sym}) \leq \lambda_{max}((\mathbf{MCD})_{sym} - \mathbf{M})$. Therefore any stability condition that holds for the connection matrix $\mathbf{C}$ must also hold for the larger $\mathbf{W}$ matrix.

It is worth noting that highly negative diagonal values in $\mathbf{W}$ will prevent the same metric $\mathbf{M}$ from being used for the nonlinear system. Therefore the method used in this proof cannot feasibly be adapted to further relax the treatment of the diagonal part of $\mathbf{W}$. The intuitive reason behind this is that in the symmetric part of the Jacobian, $\frac{\mathbf{M}\mathbf{W}\mathbf{D} + \mathbf{D}\mathbf{W}^{T}\mathbf{M}}{2} - \mathbf{M}$, the diagonal self weights will also be scaled down by small $\mathbf{D}$, while the leak portion $-\mathbf{M}$ remains untouched by $\mathbf{D}$.

Now we actually demonstrate a counterexample, presenting a $2 \times 2$ symmetric Metzler matrix $\mathbf{W}$ that is contracting in the identity in the linear system, but cannot be contracting in the identity in the nonlinear system (<ref>):

\begin{align*}
\mathbf{W} = \begin{bmatrix}
-9 & 2.5 \\
2.5 & 0
\end{bmatrix}
\end{align*}

To see that it is not possible for the more general nonlinear system with these weights to be contracting in the identity, take
$\mathbf{D} = \begin{bmatrix}
0 & 0 \\
0 & 1
\end{bmatrix}$. Now

\begin{align*}
(\mathbf{W}\mathbf{D})_{sym} - \mathbf{I} = \begin{bmatrix}
-1 & 1.25 \\
1.25 & -1
\end{bmatrix}
\end{align*}

which has a positive eigenvalue of $\frac{1}{4}$.

Note that this proof has focused on the case $g=1$ (i.e. activation function derivative upper bounded at 1), which most of the commonly used activation functions at present satisfy. However, a version of the condition can be derived for more general $g$ as well (supplemental section <ref>).

§.§.§ A contraction condition for guaranteed stability throughout training

While there are many interesting network structures that satisfy Theorem <ref>, and a number of potential methods to push networks towards satisfying Theorem <ref> during training, strictly enforcing that condition throughout optimization is not straightforward. This motivates our next theorem, which derives a different continuous-time contraction condition that does have a direct parameterization that can easily be plugged into modern optimization libraries.

The theorem states:

If there exists a positive diagonal matrix $\mathbf{P}$ such that:

\[g^2\mathbf{W}^T\mathbf{P}\mathbf{W} - \mathbf{P} \prec 0 \]

then (<ref>) is contracting in metric $\mathbf{P}$.

This condition is an extension of checking the maximum singular value of $\mathbf{W}$, but allowing a constant metric that can expand the set of matrices it can detect as stable – demonstrating another advantage of contraction analysis (Figure <ref>). It reduces to the former condition when $\mathbf{M}=\mathbf{I}$. A discrete time version of this contraction condition had been previously proven by [Revay and Manchester, 2020], but this is the first proof in continuous time to our knowledge.

[Our novel single-RNN conditions capture a larger set of stable nonlinear network weight matrices.]Considering different contraction metrics can capture a larger set of stable network weight matrices. Linear projection (using principal components analysis) of randomly generated matrices satisfying the stability condition in Theorem <ref> with $\mathbf{M} = \mathbf{I}$ and $\mathbf{M} = \mathbf{P}$.

For a direct parameterization of the theorem, note the above condition on $\mathbf{W}$ is satisfied if $\mathbf{W}$ can be written as:

\[\mathbf{W} = g^{-2}\mathbf{P}^{-1/2}\mathbf{U}\mathbf{S}\mathbf{V}^T\mathbf{P}^{1/2}\]

where $\mathbf{U},\mathbf{V}$ are orthogonal matrices and $\mathbf{S}$ is a non-negative diagonal matrix with $ 0 \leq S_{ii} < 1$. The matrix $e^{\mathbf{N}}$ with $\mathbf{N} = -\mathbf{N}^T$ is an orthogonal matrix, so this can be used to parameterize $\mathbf{U}$ and $\mathbf{V}$.

The proof of Theorem <ref> is brief and presented next. Our work on combination properties will be reported on following the proof, in section <ref>.

Proof of Theorem <ref>.

Consider the generalized Jacobian:

\[\mathbf{F} = \mathbf{P}^{1/2}\mathbf{J} \mathbf{P}^{-1/2} = -\mathbf{I} + \mathbf{P}^{1/2}\mathbf{W}\mathbf{P}^{-1/2}\mathbf{D}  \]

where $\mathbf{D}$ is a diagonal matrix with $\mathbf{D}_{ii} = \frac{d\phi_i}{dx_i} \geq 0$. Using the subadditivity of the matrix measure $\mu_2$ of the generalized Jacobian we get:

\[\mu_2(\mathbf{F}) \leq -1 + \mu_2(\mathbf{P}^{1/2}\mathbf{W}\mathbf{P}^{-1/2}\mathbf{D}) \]

Now using the fact that $\mu_2(\cdot) \leq ||\cdot||_2$ we have:

\[\mu_2(\mathbf{F}) \leq -1 + ||\mathbf{P}^{1/2}\mathbf{W}\mathbf{P}^{-1/2}\mathbf{D})||_2 \leq -1 + g||\mathbf{P}^{1/2}\mathbf{W}\mathbf{P}^{-1/2}||_2 \]

Using the definition of the 2-norm, imposing the condition $\mu_2(\mathbf{F}) \leq 0$ may be written:

\[g^2\mathbf{W}^T\mathbf{P}\mathbf{W} - \mathbf{P} \prec 0 \]

which completes the proof.

§.§ Generalizing to the network of networks: "RNNs of RNNs"

The RNN in (<ref>) is a subnetwork of the overall model we introduce here. Our goal is to combine these subnetworks into a 'network of networks' in a manner that preserves stability of the overall network as long as the underlying modules are stable. We will make use of both of the theorems from section <ref> to enforce contraction in subnetworks for our experiments, but other contraction conditions could also be used. Through the rest of this section on network combinations, we will thus assume we are given a set of subnetworks that are provably contracting in a given set of corresponding metrics, and we will focus only on combining them in a provably stable way.

Unsurprisingly, we will utilize the contraction combination properties discussed within section <ref> to do so. We focus on negative feedback combinations in particular due to their inherently recurrent structure. For reasons hinted at in discussing Theorem <ref>, which will be covered in substantially more detail with Theorem <ref> below, it is difficult to utilize negative feedback connections in a provably stable way when those connections are nonlinear. However, stability of linear negative feedback connections between the subnetworks follows directly from the existing contraction combination properties.

Note that the connections within each subnetwork will remain nonlinear as described by (<ref>), it is only the connections between subnetworks that take the form of linear negative feedback. Therefore the system on the whole remains nonlinear. The linearity assumption can also be motivated by the fact that RNNs have been found to be well-approximated by linear systems in many neuroscience contexts [Sussillo and Barak, 2013, Langdon and Engel, 2022].

Here, we formally apply the negative feedback contraction combination property to our model, and prove that it can be parameterized in a way that allows the strengths of different negative feedback connections within the network to be learned using modern optimization tools. We thus make extensive use of the following theorem:

Consider a collection of $p$ subnetwork RNNs governed by (<ref>). Assume that these RNNs each have hidden-to-hidden weight matrices $\{ \mathbf{W}_1,\dots,\mathbf{W}_p \}$ and are independently contracting in metrics $\{ \mathbf{M}_1,\dots,\mathbf{M}_p \}$. Define the block matrices $\Tilde{\mathbf{W}} \equiv \text{BlockDiag}(\mathbf{W}_1,\dots,\mathbf{W}_p )$ and $\Tilde{\mathbf{M}} \equiv \text{BlockDiag}(\mathbf{M}_1,\dots,\mathbf{M}_p )$, as well as the overall state vector $\Tilde{\mathbf{x}}^T \equiv (\mathbf{x}^T_1 \cdots \mathbf{x}^T_2)$. Then the following 'network of networks' is globally contracting in metric $\Tilde{\mathbf{M}} $:

\begin{equation}\label{eq:comboRNN}
\begin{split}
\tau \dot{\Tilde{\mathbf{x}}} = -\Tilde{\mathbf{x}} + \Tilde{\mathbf{W}}\phi(\Tilde{\mathbf{x}}) + \mathbf{u}(t) + \mathbf{L}\Tilde{\mathbf{x}} \\
\mathbf{L} \equiv \mathbf{B} - \Tilde{\mathbf{M}}^{-1}\mathbf{B}^T\Tilde{\mathbf{M}}
\end{split}
\end{equation}

Where $\mathbf{B}$ is an arbitrary square matrix, and can be treated as a trainable tensor in a deep learning framework.

In the experiments of section <ref>, we freeze the block diagonal elements of $\mathbf{B}$ at 0, to prevent linear negative feedback connections from forming within the subnetworks, thereby preserving the interpretation of our model as presented conceptually.

In theory, the feedback combination property, taken together with hierarchical combinations, may be used as combination primitives for recursively constructing complicated networks of networks while automatically maintaining stability. The recursion comes from the fact that once a modular system is shown to be contracting it may be treated as a single contracting system, which may in turn be combined with other contracting systems, ad infinitum. Note that while the negative feedback currently requires linear interareal connections, hierarchical interareal connections can be nonlinear [Lohmiller and Slotine, 1998]. Although we do not formally explore this in the paper, it is another potentially fruitful direction in which our architecture could be taken.

To prove Theorem <ref>, we show that in the appropriate metric the connections in $\mathbf{L}$ boil down to pure antisymmetric negative feedback, i.e. $\mathbf{L}_{ij} = -\mathbf{L}_{ji}^T$. Recall that the metric refers to the "space" in which we know the system has all trajectories exponentially converge. When $\mathbf{M} = \mathbf{I}$, this is simply Euclidean space and we would indeed look for skew-symmetry directly. It is of course possible for a system to be contracting in multiple different metrics, though we only need one to use the Theorem and it is considered a given input there. The intuition behind the negative feedback combination property itself is that because contraction analysis relies on analyzing the symmetric part of Jacobian matrices, the skew-symmetry of $\mathbf{L}_{ij}$ 'cancels out' when  computing the symmetric part, and leaves the stability of the subnetwork RNNs untouched.

§.§.§ Proof of Theorem <ref>

Consider the differential Lyapunov function:

\[V = \frac{1}{2}\delta\mathbf{x}^T\Tilde{\mathbf{M}}\delta\mathbf{x} \]

The time-derivative of this function is:

\[\dot{V} = \delta\mathbf{x}^T\Tilde{\mathbf{M}}\delta\dot{\mathbf{x}} = \delta\mathbf{x}^T (\underbrace{\Tilde{\mathbf{M}}\Tilde{\mathbf{J}} + \Tilde{\mathbf{J}}^T\Tilde{\mathbf{M}}}_{\text{Jacobian of RNNs before interconnection}} + \underbrace{\Tilde{\mathbf{M}}\mathbf{L} + \mathbf{L}^T\Tilde{\mathbf{M}}}_{\text{Interconnection Jacobian}} )\delta\mathbf{x}\]

Since we assume that the RNNs are contracting in isolation, the first term in this sum is less that the slowest contracting rate of the individual RNNs, which we call $\lambda > 0$. The second term in the sum is the zero matrix, by construction. Thus the time-derivative of $V$ is upper-bounded by:

\[\dot{V} \leq -2\lambda V\]

This implies that the network of networks is contracting with rate $\lambda$.

§.§ Relation of stability conditions with neuron pruning

Throughout this chapter, the word "subnetwork" is generally used to refer to an individual nonlinear RNN as defined in (<ref>), a set of which will then be combined using linear connections that meet the combination property constraints detailed above in <ref>. However, for this subsection only, "subnetwork" will refer to a network defined by taking only a subset of the units in a given RNN with (<ref>), and considering an RNN defined by the same dynamics but with the unchosen neurons entirely excluded (i.e. considering $\mathbf{W}$ to be the principal submatrix corresponding to only the selected units, and keeping the states and inputs of only the selected units).

While we have mainly focused on sufficient conditions for the system (<ref>) to be contracting regardless of (eligible) activation function, identifying necessary conditions is also of great potential use. Interestingly, one such necessary condition, presented in Theorem <ref>, has clear interpretations in neuroscience and machine learning, mirroring some of the implications of the combination properties described above, but now in the inverse case of pruning a system rather than building one up. Furthermore, the necessary condition again stems from the presence of both positive and negative connections between neurons.

In order for a linearly stable RNN to be contracting in all of the nonlinear activation functions permitted by (<ref>), every subnetwork must also be linearly stable. Put another way, it is necessary for all principal submatrices of the matrix $\mathbf{W}$ to have only eigenvalues with negative real part.

Additionally, we show that if the nonlinear system (<ref>) is contracting in a diagonal metric, all of its subnetworks are also contracting in a diagonal metric.

Note that as linear activation is a possible function for $\phi$ in (<ref>), it is definitionally true that any $\mathbf{W}$ must be Hurwitz to possibly be stable in the more general case we consider. Thus the necessary condition here could also be formulated as "$\mathbf{W}$ and all of its principal submatrices must be linearly stable". Obviously it may be possible for the network to still be stable in a specific activation function of interest, but that is not the goal of our investigation, and furthermore the popular ReLU function (which we also use in our experiments) is strictly more complicated than the linear $\phi$, and for nonlinear stability with ReLU we do require this necessary condition, as can be seen in its proof.

Theorem <ref> also connects back to Theorem <ref>. All diagonally stable linear networks will meet the necessary condition for potential nonlinear stability, because diagonal stability guarantees that all principal submatrices are also diagonally stable [Narendra and Shorten, 2010]. This will of course include all Hurwitz Metzler matrices, as we know it must due to the sufficient condition of Theorem <ref>. Furthermore, Theorem <ref> shows that the same property is true in the (<ref>) system – if it is contracting in a diagonal metric, all of its subnetworks are also contracting in a diagonal metric. Because Theorem <ref> guarantees a diagonal metric for the nonlinear system, we know that removing units from an RNN that satisfies Theorem <ref> will not impact stability. This is not generally true, as one can imagine a system that is relying on a highly inhibitory unit to enforce stability.

Enforcing stability of subnetworks has a number of functional benefits. Robustness to loss of units has obvious safety implications, and furthermore the ability to remove units may also play a role in learning and development. Dropout, a powerful regularization technique for deep neural networks, involves randomly removing a fraction of units for each training round [Bengio et al., 2017]. Many Neural Architecture Search algorithms include removing units as a potential step. Even in biology, a large number of neurons undergo programmed cell death during human development, and this pattern is conserved across many animals [Hutchins and Barger, 1998]. It is likely important not only to be able to combine stable blocks without concern for stability of the overall system, but also to be able to take blocks away without concern. Additionally, work on conditions related to the pruning of neurons could pave the way for more complex conditions dealing with pruned weights, a topic of great recent interest in neurobiology and machine learning [Frankle and Carbin, 2019].

Recall that our system (<ref>) allows for different activation functions to be used with different units. The core idea of the Theorem <ref> proof exploits that fact – although a similar effect could be replicated in a network with all ReLU activation, through targeting highly negative inputs that will put the unit below the activation threshold only at those units we want to "deactivate". Further intuition and then the mathematical details of the proof follow. To skip to the next topic on incorrect prior "stability" theorems in machine learning, see section <ref>.

§.§.§ Proof of Theorem <ref>

From an adversarial perspective, unstable subnetworks could be exploited by the activation functions purely via linear gains, so this necessary condition does not even require a dynamic $\mathbf{D}$. If there is an unstable subnetwork within our network, we can simply assign linear activation to those units in that subnetwork, and constant activation to all other units. This will produce a $\mathbf{D}$ in the Jacobian that zeros outputs from any units not in the unstable subnetwork, essentially generating a network that is the hierarchical combination of the unstable subnetwork with each other individual unit. Therefore the unstable portion of the network will have its behavior unchecked, and as it provides inputs to the other units can drive the entire network to "blow up".

Suppose that some principal submatrix $\mathbf{W}_{S}$ of an $n \times n$ weight matrix $\mathbf{W}$ has an eigenvalue with real part in the right half plane. Without loss of generality (as units can be renumbered), assume that indices 0 through $m-1$ of $\mathbf{W}$ are also in $\mathbf{W}_{S}$, and indices $m$ through $n-1$ are not in $\mathbf{W}_{S}$. Then we can set $\phi_{i}$ to be the identity function for $i < m$, and set $\phi_{i}$ to be a constant function $f(x) = 1$ for $i \geq m$. Now $\mathbf{D}$ is a static diagonal matrix with 1 on diagonals up to index $m-1$, and 0 on diagonals after. $\mathbf{W}\mathbf{D}$, which can also be viewed as a linear system and therefore evaluated by linear stability criteria, is equal to:

\begin{align*}
\begin{bmatrix}
\mathbf{W}_{S} & \mathbf{0} \\
\mathbf{W}_{S_{out}} & \mathbf{0}
\end{bmatrix}
\end{align*}

where $\mathbf{W}_{S_{out}}$ is a matrix containing the connections from the first $m$ neurons to the last $n-m$ neurons.

Observe that this is a hierarchical combination, with $\mathbf{W}_{S}$ receiving no inputs. Therefore it is impossible for the overall system $\mathbf{W}\mathbf{D}$ to be stable, as it contains an unstable system at the top of its hierarchy. This completes the proof of the necessary condition.

As mentioned, we will additionally prove that diagonal stability of the nonlinear system (<ref>) guarantees all subnetworks of that system will also be stable in a diagonal metric. This extends a known result for linear diagonal stability [Narendra and Shorten, 2010] to the nonlinear RNN, and is of particular relevance here due to both results such as Theorem <ref> that guarantee diagonal stability of the nonlinear system, and also the theoretical benefits of subnetwork stability described in section <ref>.

To begin, suppose we are given an RNN with weights $\mathbf{W}$, contracting in known diagonal metric $\mathbf{M}$. Then we have:

\begin{align*}
(\mathbf{M}\mathbf{W}\mathbf{D})_{sym} - \mathbf{M} \prec 0
\end{align*}

Now suppose we want to keep an arbitrary subset, $S$, of the neurons in $\mathbf{W}$. We will call the principal submatrix of $\mathbf{W}$ containing only the indices in $S$ $\mathbf{W}_{S}$, and we will prove that this network is contracting in the analogously defined $\mathbf{M}_{S}$. As of course only the activation functions that act on the kept neurons will remain relevant, we will need to show:

\begin{align*}
(\mathbf{M}_{S}\mathbf{W}_{S}\mathbf{D}_{S})_{sym} - \mathbf{M}_{S} \prec 0
\end{align*}

Observe that with diagonal $\mathbf{M}$, taking the principal submatrix of the overall Jacobian is equivalent to taking the principal submatrix of each component independently, i.e.

\begin{align*}
(\mathbf{M}_{S}\mathbf{W}_{S}\mathbf{D}_{S})_{sym} - \mathbf{M}_{S} = ((\mathbf{M}\mathbf{W}\mathbf{D})_{sym} - \mathbf{M})_{S}
\end{align*}

This can be seen in steps, first focusing on the $\mathbf{M}_{S}\mathbf{W}_{S}\mathbf{D}_{S}$ term, as this equals $(\mathbf{M}\mathbf{W}\mathbf{D})_{S}$ due to $\mathbf{M}$ and $\mathbf{D}$ both being diagonal. $\mathbf{M}$ simply weights each row by the corresponding diagonal element, so it is irrelevant whether a row is discarded before or after multiplication - and analogously for $\mathbf{D}$ on the columns.

Next, note that the symmetric part can also be taken before or after subsetting, as it is simply the sum of a matrix and its transpose multiplied by a scalar. The sum operates elementwise, so again there is no reason for discard time to matter. This can also be applied to subtraction of the $\mathbf{M}_{S}$ term, completing the transformation.

To complete the proof, we use Cauchy's Interlace Theorem as described by [Hwang, 2004] to show that:

\begin{align*}
((\mathbf{M}\mathbf{W}\mathbf{D})_{sym} - \mathbf{M})_{S} \prec (\mathbf{M}\mathbf{W}\mathbf{D})_{sym} - \mathbf{M}
\end{align*}

with the latter term already known to be uniformly negative definite.

Because we are concerned with the maximum eigenvalues of symmetric matrices, we can guarantee we will only be working with real values. Thus Cauchy's Interlace Theorem can be applied directly, as it states that given a square Hermitian matrix $\mathbf{A}$ and a principal submatrix of $\mathbf{A}$, $\mathbf{B}$, then $\lambda_{max}(B) < \lambda_{max}(A)$.

§.§ Identifying counterexamples to previously published RNN stability theorems

Several recent papers in ML, e.g  [Haber and Ruthotto, 2017, Chang et al., 2019], claim that a sufficient condition for stability of the nonlinear system:

\[\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x},t)\]

is that the associated Jacobian matrix $\mathbf{J}(\mathbf{x},t) = \frac{\partial \mathbf{f}}{\partial \mathbf{x}}$ has eigenvalues whose real parts are strictly negative, i.e:

\[\max_i \text{Re}(\lambda_i(\mathbf{J}(\mathbf{x},t)) \leq -\alpha\]

with $\alpha>0$. However, this claim is generally false - see Section 4.4.2 in [Slotine and Li, 1991]. Recall that contraction analysis evaluates only the symmetric part of the Jacobian. For any matrix $\mathbf{A}$, the largest eigenvalue of the symmetric part will always be $\geq Re(\lambda_{max}(\mathbf{A}))$.

It is of course not necessarily the case that certain claimed conditions are inherently false, as we do not have counterexamples yet that directly demonstrate that e.g. a diagonally stable $\mathbf{W}$ was actively unstable in the (<ref>) system. We simply point out that some existing proofs make demonstrably false claims, and thus we have no guarantee that (<ref>) will be globally contracting just because $\mathbf{W}$ is diagonally stable. In the rest of this section, we provide more detail on why certain conditions have been difficult to prove, and what additional constraints might be imposed in order for such properties to guarantee stability.

In the specific case of the RNN (<ref>), it appears that the eigenvalues of the symmetric part of $\mathbf{W}$ do provide information on global stability in a number of applications. For example, [Matsuoka, 1992] showed that if $\mathbf{W}_{sym} = \frac{\mathbf{W}+\mathbf{W^T}}{2}$ has all its eigenvalues less than unity, and the input $\mathbf{u}$ is constant, then (<ref>) has a unique, globally asymptotically stable fixed point. This condition also implies that the real parts of the eigenvalues of the Jacobian are uniformly negative. Moreover, it has been shown that setting the symmetric part of $\mathbf{W}$ almost equal to zero (yet slightly negative) led to rotational, yet stable dynamics in practice [Chang et al., 2019]. This leads us to the following theorem, which shows that if the slopes of the activation functions change sufficiently slowly as a function of time, then the condition in [Matsuoka, 1992] in fact implies global contraction of (<ref>).

Let $\mathbf{D}$ be a positive, diagonal matrix with $D_{ii} = \frac{d\phi_i}{dx_i}$, and let $\mathbf{P}$ be an arbitrary, positive diagonal matrix. If:

\[ (g\mathbf{W}-\mathbf{I})\mathbf{P} + \mathbf{P}(g\mathbf{W}^T-\mathbf{I}) \preceq -c\mathbf{P} \hspace{.5cm} \text{and} \hspace{.5cm} \dot{\mathbf{D}} - cg^{-1}\mathbf{D} \preceq -\beta\mathbf{D}\]

for $c,\beta > 0$, then (<ref>) is contracting in metric $\mathbf{D}$ with rate $\beta$.

We stress however, that it is an open question whether or not diagonal stability of $\mathbf{W}$ implies that (<ref>) is contracting. Note that the ReLU activation function has slope only $=$ 0 or 1, such that targeted inputs can theoretically cause rapidly changing $\mathbf{D}$, thus making it much less likely to satisfy Theorem <ref>.

It has in fact repeatedly been conjectured that diagonal stability of $g\mathbf{W}-\mathbf{I}$ is a sufficient condition for global contraction of the (<ref>) system, for example by [Revay et al., 2020]. However this has been difficult to prove, with multiple mistaken proof attempts appearing and later disappearing in various preprints, and many months ourselves struggling with what seems like a relatively simple claim. To better characterize this conjecture, we present Theorem <ref>, which shows by way of counterexample that diagonal stability of $g\mathbf{W}-\mathbf{I}$ can not imply global contraction in a constant metric for the general (<ref>) system.

Satisfaction of the condition

\begin{align*}
g\mathbf{W}_{sym} - \mathbf{I} \prec 0
\end{align*}

is not sufficient to show global contraction of the general nonlinear RNN (<ref>) in any constant metric. High levels of antisymmetry in $\mathbf{W}$ can make it impossible to find such a metric, which we demonstrate via a $2 \times 2$ counterexample of the form

\begin{align*}
\mathbf{W} = \begin{bmatrix}
0 & -c \\
c & 0
\end{bmatrix}
\end{align*}

with $c \geq 2$ when $g=1$.

Note that $g\mathbf{W}_{sym} - \mathbf{I} = g\frac{\mathbf{W} + \mathbf{W}^{T}}{2} - \mathbf{I} \prec 0$ is equivalent to the condition for contraction of the system (<ref>) with linear activation in the identity metric. It is also possible to show that if $g\mathbf{W}_{sym} - \mathbf{I} \prec 0$ were a sufficient condition, it would follow that contraction of the linear system in any diagonal metric is a sufficient condition – hence the discussion of diagonal stability of $\mathbf{W}$ as a conjectured condition. Recall that linear diagonal stability of $\mathbf{W}$ does at least meet the necessary condition of Theorem <ref>, whereas for Hurwitz $\mathbf{W}$ without diagonal stability, it may still be impossible to achieve the more general nonlinear stability. It would be necessary to check for linear stability of all the principal submatrices independently in that case, as described in section <ref> (and that would only be a necessary condition, not a sufficient one).

Theorem <ref> means that attempts to prove the conjectured condition with a constant contraction metric are in actuality futile. It is possible for a system to be contracting in a varying metric, and it is certainly possible that the conjecture is true because of this. However, it will be extremely difficult to prove in the general case, because with a metric that varies, the derivative of that metric also appears in the contraction condition inequality. For the present system, that will involve a $\dot{\mathbf{D}}$ term; but to bound $\dot{\mathbf{D}}$ is complicated, because it is implicitly a function of the input (as $D_{ii} = \phi_i''\dot{x}_i$).

Theorem <ref> thus suggests that future work might want to either focus more on local notions of contraction or introduce additional constraints on the (<ref>) system to make proving '$\mathbf{W}_{sym} - \mathbf{I} \prec 0 \implies$ global contraction' more approachable. One might also consider searching for other forms of conditions that can serve new purposes despite being less permissive than $\mathbf{W}_{sym} - \mathbf{I} \prec 0$ – akin to those we presented in section <ref>.

It is worth highlighting that Theorem <ref> has deep ties to our methods for combining RNNs as well. If we allowed the negative feedback connections to be nonlinear in the same way that we allow nonlinearities in (<ref>), the exact same problem would be encountered. Because the combined system would be fully expressible in the terms of (<ref>) in that case, and Theorem <ref> has just shown we cannot easily prove global contraction with high levels of skew-symmetry, we would therefore lose our stability guarantee. This also underscores the relevance of continued theoretical research along these lines.

The proof for Theorem <ref> follows in subsection <ref>. I will then present the proof for Theorem <ref>, along with intuition on the dynamics that make this Theorem true, in subsection <ref>. That will conclude the mathematical results; to proceed to the experimental work, where the theory was applied for proof of concept success in sequence learning benchmarks, see section <ref>.

§.§.§ Proof of Theorem <ref>

Consider the differential, quadratic Lyapunov function:

\[V = \delta\mathbf{x}^T\mathbf{P}\mathbf{D} \delta\mathbf{x}\]

where $\mathbf{D} \succ 0 $ is as defined above. The time derivative of $V$ is:

\[ \begin{split} \dot{V} = \delta\mathbf{x}^T \mathbf{P}\dot{\mathbf{D}}\delta\mathbf{x} +\delta\mathbf{x}^T (-2\mathbf{P}\mathbf{D} + \mathbf{P}\mathbf{D}\mathbf{W}\mathbf{D} + \mathbf{D}\mathbf{W}^T\mathbf{D}\mathbf{P} )\delta\mathbf{x}  \end{split} \]

The second term on the right can be factored as:

\[\begin{split}\delta\mathbf{x}^T (-2\mathbf{P}\mathbf{D} + \mathbf{P}\mathbf{D}\mathbf{W}\mathbf{D} + \mathbf{D}\mathbf{W}^T\mathbf{D}\mathbf{P} )\delta\mathbf{x} = \\
\delta\mathbf{x}^T\mathbf{D} (-2\mathbf{P}\mathbf{D}^{-1} + \mathbf{P}\mathbf{W} + \mathbf{W}^T\mathbf{P} )\mathbf{D} \delta\mathbf{x} \leq \\
\delta\mathbf{x}^T\mathbf{D} (-2\mathbf{P}g^{-1} + \mathbf{P}\mathbf{W} + \mathbf{W}^T\mathbf{P} )\mathbf{D} \delta\mathbf{x} = \\
\delta\mathbf{x}^T\mathbf{D} [\mathbf{P}(\mathbf{W} - g^{-1}\mathbf{I}) + (\mathbf{W}^T - g^{-1}\mathbf{I})\mathbf{P} ]\mathbf{D} \delta\mathbf{x} \leq \\
-cg^{-1}\delta\mathbf{x}^T\mathbf{P}\mathbf{D}^2 \delta\mathbf{x}
\end{split} \]

where the last inequality was obtained by substituting in the first assumption above. Combining this with the expression for $\dot{V}$, we have:

\[ \begin{split} \dot{V} \leq \delta\mathbf{x}^T \mathbf{P}\dot{\mathbf{D}}\delta\mathbf{x} -cg^{-1}\delta\mathbf{x}^T \mathbf{P}\mathbf{D}^2\delta\mathbf{x}  \end{split} \]

Substituting in the second assumption, we have:

\[ \begin{split} \dot{V} \leq \delta\mathbf{x}^T \mathbf{P}(\dot{\mathbf{D}} -cg^{-1} \mathbf{D}^2)\delta\mathbf{x}  \leq -\beta \delta\mathbf{x}^T\mathbf{P}\mathbf{D}\delta\mathbf{x} = -\beta V \end{split} \]

and thus $V$ converges exponentially to $0$ with rate $\beta$.

§.§.§ Proof of Theorem <ref>

The main intuition behind the counterexample of this theorem is that high levels of antisymmetry can prevent a constant metric from being found in the nonlinear system. This is because (in the simplest and most ubiquitous $g=1$ case) $\mathbf{D}$ is a diagonal matrix with values between 0 and 1, so the primary functionality it can have in the symmetric part of the Jacobian is to downweight the outputs of certain neurons selectively. In the extreme case of all 0 or 1 values (which is in fact the possible values that can be found in $\mathbf{D}$ when $\phi = $ ReLU), we can think of this as selecting a subnetwork of the original network, and taking each of the remaining neurons to be single unit systems receiving input from that subnetwork. For a given static configuration of $\mathbf{D}$ (think linear gains), this is a hierarchical system that will be stable if the subnetwork is stable (which it is for diagonally stable $\mathbf{W}$). But as $\mathbf{D}$ can evolve over time with many of the considered nonlinearities – including ReLU, where targeted inputs could theoretically switch an entry of $\mathbf{D}$ to 0 or 1 at will – we would need to find a constant metric that can serve completely distinct hierarchical structures simultaneously in order to accommodate those cases where negative feedback enabled high weight magnitudes. As it turns out, such a constant metric does not exist even in a simple 2 neuron negative feedback example.

To sketch the intuition of the proof in terms of matrix algebra, observe that in the constant $\mathbf{M}$ contraction condition for the system (i.e. $\frac{\mathbf{MWD} + \mathbf{DW^{T}M}}{2} - \mathbf{M} \prec 0$), the matrix $\mathbf{D}$ can zero out columns of $\mathbf{W}$ but not their corresponding rows. So for a given weight pair $w_{ij}, w_{ji}$, which has entry in $\mathbf{W}_{sym} = \frac{w_{ij} + w_{ji}}{2}$, if $D_{i} = 0$ and $D_{j} = 1$, the $i,j$ entry in $(\mathbf{W}\mathbf{D})_{sym}$ will be guaranteed to have lower magnitude if the signs of $w_{ij}$ and $w_{ji}$ are the same, but guaranteed to have higher magnitude if the signs are different (because we are removing one and not the other). Thus if the linear system would be stable based on magnitudes alone $\mathbf{D}$ poses no real threat (as we already saw with Theorem <ref>), but if the linear system requires antisymmetry to be stable, $\mathbf{D}$ can make proving contraction quite complicated (if possible at all).

The nonlinear system is globally contracting in a constant metric if there exists a symmetric, positive definite $\mathbf{M}$ such that the symmetric part of the Jacobian for the system, $(\mathbf{M}\mathbf{W}\mathbf{D})_{sym} - \mathbf{M}$ is negative definite uniformly. Therefore $(\mathbf{M}\mathbf{W}\mathbf{D})_{sym} - \mathbf{M} \prec 0$ must hold for all possible $\mathbf{D}$ if $\mathbf{M}$ is a constant metric the system globally contracts in with any allowed activation function, as some combination of settings and inputs to obtain a particular $\mathbf{D}$ can always be found, with ReLU in particular allowing for $\mathbf{D}$ to take the form of various permutations of 0s and 1s over time as desired by using targeted inputs. We will move forward assuming that $\phi = ReLU$ then, and demonstrate by counterexample that there exists a diagonally stable $\mathbf{W}$ that cannot meet the contraction condition using the same $\mathbf{M}$ when $\mathbf{D}$ takes different configurations of 0s and 1s, values that could naturally occur for the system as defined.

Thus to prove the main claim, we present here a simple 2-neuron system that is contracting in the identity metric with linear activation function, but which can be shown to have no $\mathbf{M}$ that simultaneously satisfies the $(\mathbf{M}\mathbf{W}\mathbf{D})_{sym} - \mathbf{M} \prec 0$ condition for two different possible $\mathbf{D}$ matrices.

To begin, take

\begin{align*}
\mathbf{W} = \begin{bmatrix}
0 & -2 \\
2 & 0
\end{bmatrix}
\end{align*}

Note that any off-diagonal magnitude $\geq 2$ would work, as this is the point at which $\frac{1}{2}$ of one of the weights (found in $\mathbf{W}_{sym}$ when the other is zeroed) will have magnitude too large for $(\mathbf{W}\mathbf{D})_{sym} - \mathbf{I}$ to be stable.

Looking at the linear system, we can see it is contracting in the identity because

\begin{align*}
\mathbf{W}_{sym} - \mathbf{I} = \begin{bmatrix}
-1 & 0 \\
0 & -1
\end{bmatrix} \prec 0
\end{align*}

Now consider $(\mathbf{M}\mathbf{W}\mathbf{D})_{sym} - \mathbf{M}$ with $\mathbf{D}$ taking two possible values of

\begin{align*}
\mathbf{D}_{1} = \begin{bmatrix}
1 & 0 \\
0 & 0
\end{bmatrix}
\hspace{5mm} and \hspace{5mm}
\mathbf{D}_{2} = \begin{bmatrix}
0 & 0 \\
0 & 1
\end{bmatrix}
\end{align*}

We want to find some symmetric, positive definite
$\mathbf{M} = \begin{bmatrix}
a & m \\
m & b
\end{bmatrix}$
such that $(\mathbf{M}\mathbf{W}\mathbf{D}_{1})_{sym} - \mathbf{M}$ and $(\mathbf{M}\mathbf{W}\mathbf{D}_{2})_{sym} - \mathbf{M}$ are both negative definite.

Working out the matrix multiplication, we get

\begin{align*}
(\mathbf{M}\mathbf{W}\mathbf{D}_{1})_{sym} - \mathbf{M} = \begin{bmatrix}
2m - a & b - m \\
b - m & -b
\end{bmatrix}
\end{align*}

\begin{align*}
(\mathbf{M}\mathbf{W}\mathbf{D}_{2})_{sym} - \mathbf{M} = \begin{bmatrix}
-a & -(a + m) \\
-(a + m) & -2m - b
\end{bmatrix}
\end{align*}

We can now check necessary conditions for negative definiteness on these two matrices, as well as for positive definiteness on $\mathbf{M}$, to try to find an $\mathbf{M}$ that will satisfy all these conditions simultaneously. In this process we will reach a contradiction, showing that no such $\mathbf{M}$ can exist.

A necessary condition for positive definiteness in a real, symmetric $n \times n$ matrix $\mathbf{X}$ is $x_{ii} > 0$, and for negative definiteness $x_{ii} < 0$. Another well known necessary condition for definiteness of a real symmetric matrix is $|x_{ii} + x_{jj}| > |x_{ij} + x_{ji}| = 2|x_{ij}| \hspace{.2cm} \forall i \neq j$. See [373] for more info on these conditions.

Thus we will require $a$ and $b$ to be positive, and can identify the following conditions as necessary for our 3 matrices to all meet the requisite definiteness conditions:

\begin{align}
2m < a \label{D1diag}
\end{align}

\begin{align}
-2m < b \label{D2diag}
\end{align}

\begin{align}
|2m - (a + b)| > 2|b - m| \label{D1offdiag}
\end{align}

\begin{align}
|-2m - (a + b)| > 2|a + m| \label{D2offdiag}
\end{align}

Note that the necessary condition for $\mathbf{M}$ to be PD, $a + b > 2|m|$, is not listed, as it is automatically satisfied if (<ref>) and (<ref>) are.

It is easy to see that if $m = 0$, conditions (<ref>) and (<ref>) will result in the contradictory conditions $a > b$ and $b > a$ respectively, so we will require a metric with off-diagonal elements. To make the absolute values easier to deal with, we will check $m > 0$ and $m < 0$ cases independently.

First we take $m > 0$. By condition (<ref>) we must have $a > 2m$, so between that and knowing the signs of all unknowns are positive, we can reduce many of the absolute values. Condition (<ref>) becomes $a + b - 2m > |2b - 2m|$, and condition (<ref>) becomes $a + b + 2m > 2a + 2m$, which is equivalent to $b > a$. If $b > a$ we must also have $b > m$, so condition (<ref>) further reduces to $a + b - 2m > 2b - 2m$, which is equivalent to $a > b$. Therefore we have again reached contradictory conditions.

A very similar approach can be applied when $m < 0$. Using condition (<ref>) and the known signs we reduce condition (<ref>) to $2|m| + a + b > 2b + 2|m|$, i.e. $a > b$. Meanwhile condition (<ref>) works out to $a + b - 2|m| > 2a - 2|m|$, i.e. $b > a$.

Therefore it is impossible for a single constant $\mathbf{M}$ to accommodate both $\mathbf{D}_{1}$ and $\mathbf{D}_{2}$, so that no constant metric can exist for $\mathbf{W}$ to be contracting in when a nonlinearity is introduced that can possibly have derivative reaching both of these configurations. One real world example of such a nonlinearity is ReLU. Given a sufficiently high negative input to one of the units and a sufficiently high positive input to the other, $\mathbf{D}$ can reach one of these configurations. The targeted inputs could then flip at any time to reach the other configuration.

An additional condition we could impose on the activation function is to require it to be a strictly increasing function, so that the activation function derivative can never actually reach 0. We will now show that a very similar counterexample applies in this case, by taking

\begin{align*}
\mathbf{D}_{1*} = \begin{bmatrix}
1 & 0 \\
0 & \epsilon
\end{bmatrix}
\hspace{5mm} and \hspace{5mm}
\mathbf{D}_{2*} = \begin{bmatrix}
\epsilon & 0 \\
0 & 1
\end{bmatrix}
\end{align*}

Note here that the $\mathbf{W}$ used above
produced a $(\mathbf{W}\mathbf{D})_{sym} - \mathbf{I}$ that just barely avoided being negative definite with the original $\mathbf{D}_{1}$ and $\mathbf{D}_{2}$, so we will have to increase the values on the off-diagonals a bit for this next example. In fact anything with magnitude larger than 2 will have some $\epsilon > 0$ that will cause a constant metric to be impossible, but for simplicity we will now take

\begin{align*}
\mathbf{W}_{*} = \begin{bmatrix}
0 & -4 \\
4 & 0
\end{bmatrix}
\end{align*}

Note that with $\mathbf{W}_{*}$, even just halving one of the off-diagonals while keeping the other intact will produce a $(\mathbf{W}\mathbf{D})_{sym} - \mathbf{I}$ that is not negative definite. Anything less than halving however will keep the identity metric valid. Therefore, we expect that taking $\epsilon$ in $\mathbf{D}_{1*}$ and $\mathbf{D}_{2*}$ to be in the range $0.5 \geq \epsilon > 0$ will also cause issues when trying to obtain a constant metric.

We will now actually show via a similar proof to the above that $\mathbf{M}$ is impossible to find for $\mathbf{W}_{*}$ when $\epsilon \leq 0.5$. This result is compelling because it not only shows that $\epsilon$ does not need to be a particularly small value, but it also drives home the point about antisymmetry - the larger in magnitude the antisymmetric weights are, the larger the $\epsilon$ where we will begin to encounter problems.

Working out the matrix multiplication again, we now get

\begin{align*}
(\mathbf{M}\mathbf{W}_{*}\mathbf{D}_{1*})_{sym} - \mathbf{M} = \begin{bmatrix}
4m - a & 2b - m - 2a\epsilon \\
b - m - 2a\epsilon & -4m\epsilon - b
\end{bmatrix}
\end{align*}

\begin{align*}
(\mathbf{M}\mathbf{W}_{*}\mathbf{D}_{2*})_{sym} - \mathbf{M} = \begin{bmatrix}
4m\epsilon - a & -(2a + m - 2b\epsilon) \\
-(2a + m - 2b\epsilon) & -4m - b
\end{bmatrix}
\end{align*}

Resulting in two new main necessary conditions:

\begin{align}
|4m - a - b - 4m\epsilon| > 2|2b - m - 2a\epsilon| \label{epsD1offdiag}
\end{align}

\begin{align}
|4m\epsilon - a - b - 4m| > 2|2a + m - 2b\epsilon| \label{epsD2offdiag}
\end{align}

As well as new conditions on the diagonal elements:

\begin{align}
4m - a < 0 \label{epsD1diag}
\end{align}

\begin{align}
-4m - b < 0 \label{epsD2diag}
\end{align}

We will now proceed with trying to find $a,b,m$ that can simultaneously meet all conditions, setting $\epsilon = 0.5$ for simplicity.

Looking at $m=0$, we can see again that $\mathbf{M}$ will require off-diagonal elements, as condition (<ref>) is now equivalent to the condition $a + b > |4b - 2a|$ and condition (<ref>) is similarly now equivalent to $a + b > |4a - 2b|$.

Evaluating these conditions in more detail, if we assume $4b > 2a$ and $4a > 2b$, we can remove the absolute value and the conditions work out to the contradicting $3a > 3b$ and $3b > 3a$ respectively. As an aside, if $\epsilon > 0.5$, this would no longer be the case, whereas with $\epsilon < 0.5$, the conditions would be pushed even further in opposite directions.

If we instead assume $2a > 4b$, this means $4a > 2b$, so the latter condition would still lead to $b > a$, contradicting the original assumption of $2a > 4b$. $2b > 4a$ causes a contradiction analogously. Trying $4b = 2a$ will lead to the other condition becoming $b > 2a$, once again a contradiction. Thus a diagonal $\mathbf{M}$ is impossible

So now we again break down the conditions into $m > 0$ and $m < 0$ cases, first looking at $m > 0$. Using condition (<ref>) and knowing all unknowns have positive sign, condition (<ref>) reduces to $a + b - 2m > |4b - 2(a + m)|$ and condition (<ref>) reduces to $a + b + 2m > |4a - 2(b - m)|$. This looks remarkably similar to the $m = 0$ case, except now condition (<ref>) has $-2m$ added to both sides (inside the absolute value), and condition (<ref>) has $2m$ added to both sides in the same manner. If $4b > 2(a + m)$ the $-2m$ term on each side will simply cancel, and similarly if $4a > 2(b - m)$ the $+2m$ terms will cancel, leaving us with the same contradictory conditions as before.

Therefore we check $2(a + m) > 4b$. This rearranges to $2a > 2(2b - m) > 2(b - m)$, so that from condition (<ref>) we get $b > a$. Subbing condition (<ref>) in to $2(a + m) > 4b$ gives $8b < 4a + 4m < 5a$ i.e. $b < \frac{5}{8}a$, a contradiction. The analogous issue arises if trying $2(b - m) > 4a$. Trying $2(a + m) = 4b$ gives $m = 2b - a$, which in condition (<ref>) results in $5b - a > |6a - 6b|$, while in condition (<ref>) leads to $5a > 8b$, so (<ref>) can further reduce to $5b - a > 6a - 6b$ i.e. $11b > 7a$. But $b > \frac{7}{11}a$ and $b < \frac{5}{8}a$ is a contradiction. Thus there is no way for $m > 0$ to work.

Finally, trying $m < 0$, we now use condition (<ref>) and the signs of the unknowns to reduce condition (<ref>) to $a + b + 2|m| > |4b - 2(a - |m|)|$ and condition (<ref>) to $a + b - 2|m| > |4a - 2(b + |m|)|$. These two conditions are clearly directly analogous to in the $m > 0$ case, where $b$ now acts as $a$ with condition (<ref>) being $b > 4|m|$. Therefore the proof is complete.

It is worth mentioning that an additional condition we could impose on the activation function is to require it to be a strictly increasing function, so that the activation function derivative can never actually reach 0. This is sometimes an assumption made on the (<ref>) network in the literature for simplicity. However, even with such an assumption, the counterexample of Theorem <ref> still holds, as I prove in supplemental section <ref>. An example activation function that would fit this definition in real life is softplus, a smooth approximation to ReLU.

§ SEQUENCE LEARNING EXPERIMENTS

As an empirical proof of concept, we used our mathematical results to train provably stable "network of networks" architectures on sequence learning benchmark tasks. We focused on three sequential image classification tasks often used to evaluate RNN performance: sequential MNIST, permuted sequential MNIST, and sequential CIFAR10. These tasks measure information processing ability over long sequences [Le et al., 2015], and are particularly common to use in benchmarking provably stable RNNs, as stability in theory has a tradeoff with the ability to integrate information over long timescales – though in practice stable RNNs have largely performed well with sequence processing [Miller and Hardt, 2018].

In these tasks, images are presented pixel-by-pixel and the network makes a prediction at the end of the sequence. For permuted seqMNIST, pixels are input in a fixed but random order, to increase the distance between relevant pixels and thus make the sequence learning harder. Note that while MNIST (black and white handwritten single digit recognition) and CIFAR10 (RGB photos belonging to one of 10 classes - 6 different animals and 4 different vehicle types) are extremely easy problems for modern computer vision, we are turning them into quite difficult sequence learning tasks. Even the simpler MNIST involves 784 pixels when flattened, and the network receives just 1 pixel per timestep, only making a prediction at the very end of the sequence. Because our network is a vanilla RNN, meaning it does not include any engineered mechanism for memory and can only process sequential information based on the state evolution described by (<ref>), a high level of performance on sequence classification would be especially impressive.

In fact, the long short-term memory (LSTM) recurrent architecture with engineered gating mechanisms for retaining information was designed because vanilla RNNs were not capable of learning even moderately difficult sequence tasks in practice, including the original sequential MNIST [Le et al., 2015]. The work by [Le et al., 2015] developed a weight initialization method that allowed vanilla RNNs to compete with LSTMs on sequential MNIST; in that paper, they also introduced the harder permuted version of the task, as standard seqMNIST performance had exceeded $90\%$. Their performance on permuted seqMNIST saturated at $\sim 80\%$, while our architecture with a similar number of units (but many fewer trainable weights) achieved $\sim 90\%$ on the same task – supplemental section <ref> contains an exhaustive account of our experimental results across all trials.

Of course, models in the years since the work of [Le et al., 2015] have gotten increasingly large and in many cases have added further engineering solutions such as imposing attention mechanisms. Though permuted seqMNIST remains a good task for experimentation on smaller vanilla RNNs, it is not a difficult problem for modern machine learning at large. However, seqCIFAR10 remains quite challenging, with exceedingly few methods (none of which are vanilla RNNs) exceeding $75\%$ best test accuracy performance to date, despite often using network sizes $10x$ (or more) larger than the largest network we include here. As such, our experiments show not only state of the art (SOTA) performance amongst provably stable RNNs on these tasks, but also demonstrate highly competitive sequence learning performance against a wider range of methods.

Thus beyond the benefits of having a stability certificate, our architecture stands on its own due to the sequence learning performance it achieves for its size and for its straightforward vanilla RNN formulation, not to mention the potential for interpretability afforded by the multi-area structure. While the tasks are ultimately toy tasks in the sense that their difficulty is entirely contrived, it is important to highlight that they are indeed challenging benchmarks, and our pilot experimental results bode well for future applications of the network to real world sequence learning problems. More broadly, the experiments highlight the relevance of our major theoretical contributions, and motivate a number of other ways they could be applied (to be covered in detail within the discussion section <ref>).

In this section, I will first detail how we formalized our theoretical results within a modern deep learning framework, to train stable "RNNs of RNNs" in practice (<ref>). I will then report on our best performing architecture, including comparison to prior benchmark scores and results on the consistency of its performance over repeated trials (<ref>). Finally, I will present results on how architecture settings such as network sparsity levels (<ref>) and network size (<ref>) modulated performance, alongside discussion of how these results might inform future applications of "RNNs of RNNs".

§.§ Network architecture and methods summary

Recall that a major theme of our work is establishing foundational results for RNNs with multi-area structure. In the mathematical results of section <ref>, we focused on linear negative feedback connections between (given) stable nonlinear subnetworks, because such connections were shown to guarantee stability of the larger network, and can be easily trained using the derived parameterization with modern deep learning libraries. Here, we realized that work empirically. In particular, the matrix $\mathbf{B}$ in Theorem <ref> was made a trainable tensor in PyTorch, used in conjunction with contraction metrics of the individual subnetworks to calculate the network's linear all-to-all negative feedback weight matrix $\mathbf{L}$ (as defined in the Theorem).

In order to enforce the linear negative feedback weights $\mathbf{L}$ to form connections only between different subnetworks (but not within a subnetwork), we masked the diagonal blocks of $\mathbf{B}$ to freeze those self-connections at 0; thereby preserving an interpretation of network structure that separates the longer range feedback that different "regions" provide to each other from the weights internal to each "region". Note that because each weight has dependency on its opposite direction counterpart, we can train only half of the matrix $\mathbf{B}$ without losing degrees of freedom in the network weights $\mathbf{L}$, so to improve the efficiency of the training process we ultimately only trained the lower-triangular blocks of $\mathbf{B}$ while masking the other entries.

That was one major piece of our "RNN of RNNs" architecture design, and was used consistently throughout these pilot experiments. The other major piece would of course be the subnetworks making up the combined network, taken as a given when defining the feedback connection scheme. A benefit of our overarching framework is the level of flexibility allowed in the makeup of these subnetworks, flexibility which we leveraged to only a small extent here. It is required that every subnetwork included be definable in terms of the nonlinear RNN dynamics (<ref>), and that every subnetwork included be provably globally contracting in an obtainable metric $\mathbf{M}$. But subnetworks could span a range of sizes, activation functions, time constants, and even contraction metrics. Furthermore, for any subnetwork defined in a way that allows contraction to remain guaranteed throughout training, it is possible for those internal nonlinear subnetwork weights to be updated jointly with the negative feedback training.

In these experiments, we used our novel contraction conditions for individual nonlinear RNNs - proved in section <ref> - to enforce the needed constraints on the subnetworks. For all models, the subnetworks' dynamics were each defined according to (a discrete-time version of) equation (<ref>), using ReLU activation function across the board. Because the two considered stability conditions have quite different pros and cons however, they led to two distinct "network of networks" approaches that we explored empirically.

For Theorem <ref>, we used a bespoke initialization scheme for $\mathbf{W}$ that would guarantee contraction of each subnetwork (and provide a corresponding contraction metric), but we did not train the internal subnetwork weights at all (Figure <ref>B). For Theorem <ref>, we used a more standard initialization and jointly trained all parameters of the model (Figure <ref>C).

[Summary of proof-of-concept task structures and implemented "network of networks" architectures.]Summary of task structure and network architectures. Images from MNIST (or CIFAR10) were flattened into an array of pixels and fed sequentially one by one into the modular 'network of networks', with classification based on the output at the last timestep. For MNIST, each image was also permuted in a fixed manner (A). The subnetwork 'modules' of our architecture were constrained to meet either Theorem <ref> via sparse initialization (B) or Theorem <ref> via direct parameterization (C). Linear negative feedback connections were trained between the subnetworks according to Theorem <ref>.

We refer to networks with subnetworks constrained by Theorem <ref> as 'Sparse Combo Nets' (because sparse hidden-to-hidden weights more readily satisfy the condition), and to networks with subnetworks constrained by Theorem <ref> as 'SVD Combo Nets'. Throughout the experimental results we use the notation '$p \times n$ network' - such a network consists of $p$ distinct subnetwork RNNs, with each such subnetwork RNN containing $n$ units.

Note that task input at each timestep was first passed through a (trainable) linear input layer, represented by a weight matrix of size 1 (or 3 for RGB CIFAR pixels) by $N$, where $N = p*n$ is the total number of units in the Combo RNN. On the flip side, the outputs of the Combo RNN neurons at the final timestep in the sequence were used for prediction via passing them through a (trainable) linear output layer. The output layer was represented by a matrix of size $N$ by 10, with 10 being the number of different classes contained in both the MNIST and CIFAR10 datasets. This allowed for computation of cross entropy loss during training based on the true class of each input. At test time, the output index with the maximum value was selected as the predicted label (Figure <ref>A).

Unless specified otherwise, all networks were trained for 150 epochs, using an Adam optimizer with initial learning rate 1e-3 and weight decay 1e-5. The learning rate was cut to 1e-4 after 90 epochs and to 1e-5 after 140 epochs. Once we identified the most promising settings through our experiments, we ran final repetition trials training the best architecture for 200 epochs with learning rate cuts after epochs 140 and 190.

Because Sparse Combo Net reliably demonstrated better performance on the sequence learning tasks, and this was also the architecture that I primarily worked with, I will only briefly report on the SVD Combo Net results in this chapter. For details on the direct parameterization that allowed those subnetwork weights to be trained while maintaining contraction guarantees, see [Kozachkov et al., 2022]. On the other hand, I will next elaborate on the implementation of Sparse Combo Net in greater detail, before moving on to results (<ref>).

§.§.§ More on sparsity

As mentioned, we were not able to find a parameterization to continuously update the internal nonlinear subnetwork weights during training in a way that preserved a contraction guarantee through Theorem <ref>, and thus the Sparse Combo Net did not have the weights within its individual modules trained – relying only on training the linear connections between modules (Figure <ref>), as well as the linear input and output layers. This is reminiscent of the facilitated variation theory of evolution discussed within the background section <ref>, in that we train this combination network architecture only through training connections between contracting subnetworks, without changing the component modules at all.

[Example Sparse Combo Net weights before and after training.]Example $3$x$16$ Sparse Combo Net weight matrices, before and after training. Nonlinear intra-subnetwork weights are initialized using a set sparsity, and do not change in training (A). Linear inter-subnetwork connections are constrained to be antisymmetric with respect to the overall network metric, and are updated in training (B).

Fortunately, it is easy to randomly generate matrices with a particular likelihood of meeting the Theorem <ref> condition by selecting an appropriate sparsity level and limit on entry magnitude. Sparsity in particular is of interest due to its relevance in neurobiology and machine learning, so it is convenient that the condition makes it easy to verify stability of many different sparse RNNs. As $g=1$ for ReLU activation, we check potential subnetwork matrices $\mathbf{W}$ by simply verifying linear stability of $|\mathbf{W}| - \mathbf{I}$, keeping only those subnetworks that satisfy.

Because every RNN meeting the condition has a corresponding well-defined stable linear time-invariant (LTI) system contracting in the same metric, it is also easy to find a metric to use in our training algorithm, another advantage of Theorem <ref>. Solving for $\mathbf{M}$ in $-\mathbf{I} = \mathbf{M}\mathbf{A} + \mathbf{A}^{T}\mathbf{M}$ will produce a valid metric for any stable LTI system $\mathbf{A}$ [Slotine and Li, 1991]. Thus the $\mathbf{M}$ calculated when $\mathbf{A} = |\mathbf{W}| - \mathbf{I}$ will be a contraction metric for our nonlinear subnetwork, and can be used as needed in training the negative feedback weights per Theorem <ref>.

To solve for this $\mathbf{M}$ in practice, I integrated $e^{\mathbf{A}^{T} t} \mathbf{Q} e^{\mathbf{A} t} dt$ from 0 to $\infty$ (with scipy.integrate.quad). For efficiency reasons, and due to the guaranteed existence of a diagonal metric in the case of Theorem <ref> (as Hurwitz Metzler matrices are diagonally stable), integration was only performed to solve for the diagonal elements of $\mathbf{M}$. For safety, a check was added prior to training to confirm that the initialized network indeed satisfied Theorem <ref> with metric $\mathbf{M}$. However, it was never triggered by our initialization method.

To randomly generate the fixed nonlinear subnetworks to use for each Sparse Combo Net instantiation, I required each individual subnet weight matrix $\mathbf{W}$ to have all diagonal entries $= 0$ and only $s\%$ of off-diagonal entries non-zero. Each non-zero entry was then sampled uniformly between $-z$ and $z$. Throughout this chapter, I will therefore refer to Sparse Combo Nets containing subnetworks with just $s\%$ of their weights non-zero as having an "$s\%$ sparsity level". Both $s$ and $z$ were settings of the architecture explored during our experiments, though they were always the same across subnetworks in a given Sparse Combo Net.

Selection of these values should depend to some degree on the size of the subnetworks being used, but there are a wide variety of options that can be considered. While $s$ and $z$ need to be reasonable enough values for a subnetwork of a given size to possibly meet the Theorem <ref> condition, the random sample and the verification of the Theorem are both computationally cheap, so that the proportion of random $\mathbf{W}$ satisfying Theorem <ref> can be small if desired, as the code will continue to resample until the intended number of contracting subnetworks have been found. Note that with a smaller $s$ a larger $z$ (i.e. entry magnitude) can still be used while maintaining stability – and vice versa, though there are other theoretical reasons to prefer sparser weights with higher magnitudes to a significant extent.