Lurie Networks with Robust Convergent Dynamics
Authors: Carl R Richardson
Matthew C. Turner
Steve R. Gunn
Venue: Accepted by TMLR ()

Abstract: The Lurie network is a novel and unifying time-invariant neural ODE. Many existing continuous-time models, including recurrent neural networks and neural oscillators, are special cases of the Lurie network in this context. Mild constraints on the weights and biases of the Lurie network are derived to ensure a generalised concept of stability is guaranteed. This generalised stability measure is that of k-contraction which permits global convergence to a point, line or plane in the neural state-space. This includes global convergence to one of multiple equilibrium points or limit cycles as observed in many dynamical systems including associative and working memory. Weights and biases of the Lurie network, which satisfy the k-contraction constraints, are encoded through unconstrained parametrisations. The novel stability results and parametrisations provide a toolset for training over the space of k-contracting Lurie network's using standard optimisation algorithms. These results are also leveraged to construct and train a graph Lurie network satisfying the same convergence properties. Empirical results show the improvement in prediction accuracy, generalisation and robustness on a range of simulated dynamical systems, when the graph structure and k-contraction conditions are introduced. These results also compare favourably against other well known stability-constrained models and an unconstrained neural ODE.

Published in Transactions on Machine Learning Research (04/2025)
Lurie Networks with Robust Convergent Dynamics
Carl R. Richardson 2
cr2g 1022t0ooso..ncc.
Matthew C. Turner
m. w.t.tmmer@sstoon.coo ac. uk
Steve R. Gunn!
srgBecs.sooon.ce
School of Electronics and Computer Science, University of Southampton, Southampton, UK The Alan Turing Institute, London, UK
Reviewed on OpenReview: https: openreiiew. net/ forum?? id- 3Jm4dbbKGZ
Abstract
The Lurie network is novel and unifying time-invariant neural ODE. Many existing continuous-tine models, including recurrent neural networks and neural oscillators, are special cases of the Lurie network in this context. Mild constraints on the weights and biases of the Lurie network are derived to ensure generalised concept of stability is guaranteed. This generalised stability measure is that of k-contraction which permits global convergence to a point, line or plane in the neural state-space. This includes global convergence to one of multiple equilibrium points or limit cycles as observed in many dynamical systems including associative and working memory. Weights and biases of the Lurie network, which satisfy encoded through unconstrained paramotrisations. The
the k-contraction constraints, novel stability results and parametrisations provide a toolset for Iraining over the space of k-contracting Lurie network's using standard optimisation algorithms. These results are also leveraged to construct and train a graph Luric network satisfying the same convergence properties. Empirical results show the improvement in prediction accuracy, generalisation and robustness on a range of simulated dynamical systems, when the graph structure and k-contraction conditions are introduced. These results also compare favourably against other well known stability-constranne models and unconstrained neural ODE.
1 Introduction
A Lurie system is class of nonlinear ordinary differential equations (ODE) comprising linear time-invariant (LTI) component interconnected with potentially time-varying, nonlinearity. Such systems are ubiquitous throughout the sciences and engineering, including machine learning (ML) and neuroscience Pauli et al. (2021); Lessard et al. (2016); Lanthaler al. (2024); Wilson & Cowan (1972). When modelling time-invariant. dynamical systems, many ML models including linear state space models (LSSM), recurrent neural networks (RNN) and some graph neural networks (GNN) are special cases of Lurie system ($3.1, $B.2). Such models wide range of tasks have proven to be highly expressive demonstrated by their successful application such sequential processing Kozachkov et al. (2022a); Erichson et al. (2020); Chang et al. (2019); Gu et al. (2021), computer vision Erichson et al. (2020); Chang et al. (2019); Gu et al. (2021), language modelling Gu et al. (2021) and computational chemistry Rusch et al. (2022).
Consider dynamical systems neuroscience. The brain organises its representations of the world and carries out complex functions through collective interactions of simpler modules Kandel al. (2000). More abstractly, this be viewed graph structured dynamical system. Convergent dynamics in widespread regions of the central nervous system are thought to play a crucial role in: forming some of these representations Khona
Named after Anatolii Isakovich Larrie and sometimes spelt Lur'e Lurye.

---

Published in Transactions on Machine Learning Research (04/2025)
& Ficte (2022), processing information over extended periods Vogt et al. (2022), learning Kozachkov et al. (2020); Centorrino et. al. (2022), memory storage Hopfield (1982; 1984); Kozachkov et al. (2022b); Pals et. al. (2024) and enhancing the robustness of each of these functions Khona & Fiete (2022). As summarised in Khona & Ficte (2022), convergent dynamics in the brain take several forms. For example, ncural circuits with nultiple equilibrium points (bistable and multi-stable) have been observed in the anterlor lateral motor cortex of a rat Pict et al. (2017) and are conjectured to appear in the mammalian hippocampus and auditory cortex. Some theories of associative memory also believe reconstruction of learned pattern is obtained by flow to equilibrium points Krotov & Hopfield (2020); Sharma et al. (2022); Kozachkov et al. (2023). Limit cycles are another form of convergent dynamics for which there are numerous examples within the central nervous system. These include working memory Kozachkov et al. (2022b), which is thought to arise from the sustained spiking of neurons Ashwin et al. (2016), and sleep cycle generation Adamantidis et al. (2019). The graph structure and convergent dynamics are also shared with many other dynamical systems such as chemical processes Ofir et al. (2023), opinion dynamics and power systems Ofir ct al. (2024). As a result, encoding such structural and dynamical properties as an inductive bias is molivated by neuroscience, for a general ML framework, and for learning models of dynamical systems.
The convergence and stability analysis of dynamical systems has been well-studied in the control theory literature. A pertinent example is the absolute stability problem where the nonlinearity of the Lurie system is unknown, but assumed to be sector-bounded or slope-restricted. The goal is to find conditions on the model parameters which ensure the trajectories of all Lurie systems, with nonlinearities in the assumed class, uphold a chosen definition of convergence. Approaches to this problem can be classified as Lyapunov analysis Khalil (2002); Park (2002; 1997). Zames-Falb multipliers Zames & Falb (1968); Turner & Drummond (2019); Carrasco et al. (2016); Drummond et al. (2024) or k-contraction analysis Zhang & Cui (2013); Ofir et al. (2023; 2024). Lyapunov and Zames-Falb multiplier methods are primarily designed to analyse the convergence to an equilibrium point, whereas k-contraction methods analyse a variety of global convergence behaviours including convergence to points, lines and planes. As many ML models are examples of Lurie systems and many activation functions are sector-bourded or slope-restricted (Drummond et al., 2022, Table 1), much of the literature on the absolute stability problem is applicable to systems involving neural net works Pauli et al. (2021); Richardson et al. (2023; 2024); Fazlyab et al. (2020).
Although designing networks with convergent dynamics is well motivated, ensuring such a property requires constraints on the network parameters, which can be detrimental if too restrictive. With this in mind, this work focuses on: (i) using k-contraction analysis to derive mild constraints on the weights of the Lurie network which ensure global convergence to a point, line or plane in the neural state space, for all Lurie networks with slope-restricted activation functions; (ii) establishing unconstrained parametrisations of these conditions which allows the Lurie network to be trained using gradient based optimisation algorithms whilst limiting the search space to weights which satisfy the k-contraction conditions; (iii) constructing a graph Lurie network (GLN) from individual Lurie networks and deriving constraints the graph coupling matrix which ensure the k-contraction property is maintained. Similar unconstrained parametrisations are also derived. We empirically test the difference in prediction accuracy, generalisation and robustness of the unconstrained Lurie network, k-contracting Lurie network and GLN on a range of simulated dynamical systems dat datasels. We also compare against other well known stability architectures, which are special cases of the Lurie network, and an unconstrained neural ODE.
2 Preliminaries
2.1 Notation
For two integers i < j, we define [i,]] + The set of non-negative real numbers is denoted by R4. Symmetric matrices of dimension n are denoted by S" with the positive definite subset denoted by ST. All other positive definite subsets are denoted by a + subscript. Square diagonal matrices are denoted by D and n x m diagonal matrices are symbolised by D A positive definite (semi-definite) matrix P is sometimes indicated by P> 0 (P ž 0). Negative definite (semi-definite) matrices are indicated analogously. The set of n x n orthogonal and skew-symmetrie matrices are respectively denoted by O(n) and Skew(n). For W E Raxm, the ordered singular values are represented by a (i( 2 ... 2 o miin(n.m) (W) 2 0
2

---

Published in Transactions on Machine Learning Research (04/2025)
1-Contracting Opinion Dynamics
2-Contracting Hopfield Network
3-Contracting Simple Attractor
1.0
1.0
0.5
0.5
X3
X3
X3
0.0
005
-0.5
-1.0
-1.0
1.0
1.0
2.5
-1.0
-0.5
0.0
-0.5
0.0
X2
-0.5
0.0
-0.5
0.0
X2
-3
-5.0
-2.3
0.0
10-10
100
Figure Trajectories from three dynamical systems satisfying the k-contraction property. Crosses denote the initial condition and stars denote equilibirum points.
and for W E Rnx", the ordered eigenvalues are denoted by \ 2 (W). The k-multiplicative and k-additive compound matrices of W (see $A.1 for definitions) are respectively denoted by W and W[k], The Jacobian of a function [(1, x) is denoted by z)), The scaled 2-norm of a vector x E R with respect to (w..t)) an invertible scaling matrix O E Rnxx is defined by |x/2. == 10x|2, and the matrix measure inducod by the scaled 2-norm
(2,0 == )2(OWW-1)= A1 =1 ( +
2
2.2 k-contraction Analysis
In this work, we leverage k-contraction analysis Wu et al. (2022); Muldowney (1990), the geometrical generalisation of contraction analysis Lohmiller Slotine (1998), as a tool for controlling convergence the neural state space. Intuitively, k-contraction implies the volume of k-dimensional bodies exponentially converges to zero when governed by the system dynamics. Alternatively, this could be thought of as exponential convergence to (k - 1)-dimensional subspace. When k 1, this reduces to standard contraction Lohmiller & Slotine (1998), which implies that all trajectories exponentially converge to single trajectory. For a general time-varying dynamical system, satisfying the k-contraction property does not guarantee stability. However, for Sime-invariant dynamical systems, it has been shown that for every bounded solution: 1-contraction implies global convergence to a unique equilibrium point Tohmiller & Slotine (1998), 2-contraction implies global convergence to an equilibrium point, which not nccessarily unique but must be connocted along a line Muldowney (1990), and 3-contraction, under certain assumptions, implies convergence to non-unique attractor in 2d subspace Cecilia al. (2023). Three examples of k-contracting dynamics are presented in Figure 1.
Time-invariant dynamical systems which satisfy the k-contraction property for A E (1,2,3) have several desirable properties for ML models. They can exhibit a wide range of complex convergent behaviours such as multi-stable and orbitally stable systems Zoboli al. (2024). This suggests that model can be expressive whilst satisfying the k-contraction conditions, particularly for higher values of & where the constraints are less restrictive, as highlighted in $3.2. The k-contraction property also implies inherent. robustness as the trajectories only converge to a finite number of long term behaviours. Next, present the fundamental k-contraction result from Wu et al. (2022), see SA.3 for the formal definition of k-contraction.
Theorem Fix E [1, n] and consider the nontinear system à = f(t, x) with f R4 x R + R continuously differentiable. If there erists ni > and an invertible matriz E xa such that
22 (J (t. S V E R and t R4
(1)
then the nonlinear system is k-contracting in the 2-norm w.r.t the metric p==0 PT..
3

---

Published in Transactions on Machine Learning Research (04/2025)
This result has two features: (i) it requires the existence of an invertible matrix P. In the simplest case, can expect solution P -T,, to exist. For other systems, such simple solutions will not exist and more general matrices such as E S4 will be required, making the proofs more difficult; (ii) it requires the use of compound matrices ($A.1). For matrix W E R* the matrix W with k E [1, min m)] will have the size (2) (*) which is typically much larger and more computationally difficult to work with. A technical introduction to compound matrices, k-contraction analysis and how they relate is presented in $A. In $3.2 we derive results which verify (1) for the special case of nonlinear systems of the form (2).
3 Lurie Network
A Lurie network defined by (2) with weights A E R=xn B E Rmxm, E Rmxx and biases by E R",
by E Rm
x(t) Az(t) + (v(()) b.
vt) Cx(() + by
x00) 2o
(2)
The model has a biased linear component interconnected with a nonlinearity of the form +(v)
(0)())) ... [)))) where (b:) assumed to be slope-restricted with an upper bound g 0, such that 0 < (v) 3 gIm- This separation of the linear and nonlinear components is useful for analysis. Activation functions which satisfy this slope-restricted assumption include the hyperbolic tangent (tanh) and the rectified linear unit (ReLU). For simplicity, we assume the same scalar nonlinearity is applied element-wise and drop the subscript. The proposed Lurie network is a very general model as highlighted by the special cases in $3.1 and its relationship to deep foodforward neural networks as outlined in $B.1. Finally, it is imporlant to observe that the model is time-invariant, this implies that if Theorem is satisfied, the model will inherit the appealing convergence and robustness properties stated in $2.2. This may be viewed as limitation for modelling neural dynamics since the brain is subject to various external inputs; however, it is common to assume that, at least on the time scale of interest, the dynamics evolve in a lime-invariant manner Khona & Fiete (2022).
3.1 Example Lurie Networks
When applied to time- invariant dynamical systems, many models from the ML literature become special cases of (2). In this setting, the timo-varying external inputs are replacod with trainable biases. A subset of examples are presented next with further examples included in $B.2. As the results in $3.2 apply to any model of the form (2), they can also be applied to these special cases.
A stability constrained RNN Erichson et al. (2020) where the paramcters A, C are expressed a weighted sum of symmetric and skew-symmetrie terms in order to control the eigenvalues of the Jacobian. The remaining components are B I,, b, 0 and o(-) tanh(-).
Lipschitz RNN:
A,C E ((1 - (W + WT) + 3(W W In W E R 0.5 1, 7 1,,
Antisymmetric RNN: A constrained RNN Chang al. (2019) with the same motivations as above. Related to (2) by A 0, B In, Ce (W- 744x W Skew(n) y > 0}, by and o(-) tanh(-).
SVD Combo: A 1-contracting graph coupled RNN Kozachkov et al. (2022a) where denotes the number of individual RNNs and denotes the state dimension of each RNN. A special case of (2) with A = T, - aT == and by = O. The graph coupling matrix is denoted by L and B is block diagonal where each block contains the synaptic weights
the individual RNNs. Both of these matrices are expressed by special parametrised forms to ensure the individual RNNs and graph coupled RNN are 1-contracting. As in (2), the nonlincarity is required to be slope-restricted.
3.2 k-contraction Analysis of Lurie Networks
Two sufficient results which satisfy Theorem 1 and guarantoe (2) is k-contracting are presented next. Conditions were derived in (Ofir et al., 2024, Theorem 2) which verify Theorem 1 for Lurie network with D and by = 0. Theorem 2 extends them to account for A Ra and by * O. Refer to $C.1 for the proof.
4

---

Published in Transactions on Machine Learning Research (04/2025)
Theorem Consider the Lurie network (2) with o(y) == [0i 6m (ym) being slope-restricted
such that 0 3 J+(y) 3 glm. Pix k E [1, and define == (2k)- ):(A AT). If ag < 0 and
of (C) < aat
(3)
then (2) k-contructing in the 2-norm w.r.t the metric P - In.
The additional freedom permitted by k-contraction over standard contraction is highlighted by the summation of the eigenvalues and singular values. In I-contraction, Theorem requires the largest eigenvalue of the symmetric component of be negative whereas for & n|, this condition on becomes incrementally more relaxed as k is increased. Equation (3) illustrates a similar relaxation of the constraints on B and C. Theorem 2 has several appealing features: (i) it does not require the computation of the troublesome compound matrices; (ii) it provides a way of embedding the k-contraction properly into the structure of a Lurie network based on fairly simple unconstrained parametrisations of the weights, as shown in $4; (iii) the biases are not present in the condition, so are naturally unconstrained; (iv) the result does not rely on symmetries of the parameters. In many Hopfield-basd based models of associative memory. symmetry in the parameters is needed to make the existence of a global energy function mathematically tractable; however, this simplification is biologically unrealistic and limits the model's expressive power. The limitation of the result is that only Lurie networks which are k-contracting in scalar metric be verified.
We now present second result which addresses the scalar metric drawback; however, it comes at the cost of strong constraints on the weights B and C. See $C.2 for the proof.
Theorem 3 Consider the Lurie network (2) with ((y) == shn (y1)
@a(ua)] being slope-restricted such
that 0 < (y) 3 gln. Fir E nI. If B E D", C B-1 and
p(k) A + (A p(k) + 2kg <
(4)
then (2) k-contracling in the 2-norm w.r.l the metric PET D
Theorem 3 improves upon Theorem 2 in the sense that Lurie networks which are k-contracting in diagonal metric can now be verified; however, only when C = B- Due to this constraint, Theorem 3 has very limited practical use; however, may prove to be thalretically insightful for addressing the scalar metric drawback. Finally, it is important highlight that Theorem 2 and Theorem 3 apply the class of slope-restricted nonlincarities, so these results address the absolute stability problem for the k-contraction property.
3.3 Graph Lurie Networks
Many larger scale dynamical systems such as molccular, social, biological, and financial networks Hamilton et al. (2017) naturally have graph structure. To make the Lurie network more applicable to these problems, graph coupling term is introduced to model set of g interacting Lurie networks. To illustrate this, q independent Luric networks can be modelled by
a(1) = (Ix + (x + bR
v(1) = Cax(() + bb
x<0) = xo
(5)
where the weights are AG == blockdiag(Ap ( Ag) E Rmxxqn Bg := blockdiag(B B E Rmxam
CG == (Cadiiaag(C....... Cg) E Ramxxn and the biases are by E Rmm by E Ramh The graph Lurie network (GIN) is then defined by
#(t) (Ao L)x(t)+ Bob(v(t)) + b
y(t Cgx(t) + by
x(0) 20
(6)
where = [,,]] is a block matrix with block Ly E x connecting Lurie network Lurie network j. The state and nonlinearity of the GLN are defined by x E R0= and Φ Rem - R2 where the states of the independent Lurie networks have been stacked into single state. Interestingly, the GLN (6) is actually
5

---

Published in Transactions on Machine Learning Research (04/2025)
special case of a Lurie network; however, the networks get larger, does the search space. Thus, imposing any assumptions which respect the structure of the problem can reduce the search space and lead to more robust and generalisable models. Any further prior knowledge about the graph can be encoded through constraints on the graph coupling matrix.
3.4 k-contraction Analysis of Graph Lurie Networks
In this section, we assumee that q independent Lurie networks (2) are k-contracting in the 2-norm w.r.t metrics P1 P.....,PP. p.
This is equivalent to (5) k-contracting in the 2-norm w.r.t the metric P = blockdiag(P1..... ,,,....,PP)))) P.)). Theorem 2 and Theorem 3 provide two results which can, respectively, verify this w.r.t a scalar metric Pj = pn and a diagonal metric P3 E Dt. for E [1. Other results may be used providing they apply to systems of the form (2). Under this assumption, Theorem provides constraint, on the graph coupling term which ensures the GLN is k-contracting when constructed from q independently k-contracting Lurie networks. The proof is detailed in SC.3.
Fir E [1,n]- Consider a GLN (6) where the q independent Lurie networks are col- lectively defined by the weights Ag == (. blockdiag (CnniiiggCC......... and biases by RRn by E Ram and are k-contracting in the 2-norm w.r.t the metric P = blockdiag( b... Pg)- the graph coupling matrix I. Ranxon satisfies
Theorem
Bg == blockdiag(BI,... BB), B,), Cg ==
p(k) I + p(t) 30
(7)
then (6) is also k- k-contracting in the 2-norm w.r.t the metric P.
Remark 1 It should be noted that Theorem 4 could be applied more generally for constructing k-contracting graph coupled systems. Providing the q subsystems are k-contracting in the 2-norm w.r.t some metric P= blockdiag( (1: ....., Pq); then Theorem is applicable when coupling them through the graph coupling term. It is not necessary for the q subsystems be Lurie networks.
4 Parametrisation of kscoonraatiing Lurie Networks
To train a k-contracting Lurie network using gradient based optimisers, parametrisations which express the constrained weights in terms of unconstrained variables must be found. To formalise this idea, define the sets a2(g, k) and (k: P)- As the biases do not appear in these sets, they are naturally unconstrained. The set (2 k) contains all the weights of the Lurie network which satisfy Theorem 2.
k2lg,k) { (A, ak
1 M(++AT) O, == g² Σ (B)) < <aik)
Σ
(8a)
2k
1=1
(4(k, = [ Pmmxqn
p(k) T + ( =0]
(8b)
The set (k, P) contains the graph coupling matrices which satisfy Theorem 4. A parametrisation associated with Theorem was not established due to its limitations mentioned earlier. The next results present two different parametrisations of the set (2(g, k). See SC.4 for the proofs which leverage the eigenvalue and singular value decompositions.
Thcorem Given g > 0, E [1,m], UA. UB.UB Ve e O(n), VB, e O(m), EB E Dam, Ee E D
YA E Skew(n), GA E D and define
UUAAUU
1 Y
4 V
Σ사
In GA
(9a)
2
k.
B Vn
C. UcENV
(9b)
then (A, B, C)C) kg

---

Published in Transactions on Machine Learning Research (04/2025)
Theorem 6 Given g 0, & [1, UA.UB. Ve e O(n), VB,UC E O(m), EB E D Ee E Df". YA E (k)w(n), ΣAT E DD-11, GA2 > 0, Gas E D and define
1
A UALAUX + YA
B ::= BBBUEEVV
C:=UoICVe
(10a)
AA == blockdiag AA1 A2 A3
AA E D
(10b)
k-1
A2 == -Akkz Σ (SAI)) - -A
As == min(2A1. MA2) - -GA8
(10c)
then (A, (.B,C))
Remark 2 In both resulls, a mapping belween the sets Skew(.) and O(-) was exploited to eapress the orthogonal matrices in terms of unconstroined skew- bhew-yymmetric matrices Lezcano- & Kartanez-Rubio (2019). The remaining variables are unconstruined or simply require positive elements, which can be obtained by laking the absolule value of unconstruined variables.
In Theorem 5 and Theorem 6, the B and C matrices unconstrained since they are simply expressed by their singular value decomposition. The variables representing the singular values of B and C directly upper bound ag. The only source of conservatism in the parametrisations is introduced through the definition of ΣA- In Theorem 5, the constraint on og is enforced by a uniform negative constraint on the diagonal elements of ΣA. this assumption is true, this significantly speeds up the learning process; however, it can be prohibitive if not. Theorem 6 allows the largest (k - 1) eigenvalues to be unconstrained (E^1), meaning non-Hurwitz A matrices are encapsulated in the parametrisation. The constraint on aj is implemented by the A2 block and the Es block ensures the remaining eigenvalues are less than the other k. As l is hyperparameter of the model, it is useful to note that the set (2 k - 1) intersects with N2(g, k) when the following holds
+ 1 +
k k-1
k-1
< 41
k (
h O
(11)
We expect this to hold for many realistic systems since (11) just requires Ax(A + AT) ) to be sulficiently negative to counteract the product of the growth of B and C in the 4th direction. Hence, if the best value of k is unknown for a given application, then setting k = 3 allows the model to search over many of the stable l-cotrrcciin systems parametrised by Theorem 5 or Theorem 6, where I E [1, k). See $C.5 for the proof.
The next result provides an unconstrained parametrisation of the set (k P) (8b). Coupling this with the earlier parametrisation allows the k-contracting GLN to be trained with gradient based optimisers. See $C.6 for the proof.
Theorem Given E [I, GLI, G12 E Ronxon = blockdiag(@)...... θ₂) where O, E 8 for j E [1, P= T and define
L:= L=-ui-P-1C P
(12)
where (G12 + ( = 0, then L P),
Theorem 7 provides flexibility in the choice of GLI, GL2 define different graph structures, where the constraint on G12 is equally constraint the sum of it's k-largest eigenvalues (Fact 7). This constraint could be parametrised in simlar manner to the symmetric component of A in Theorem 5 and Theorem 6; in which case, Theorem 7 is satisfied and the graph has a directed all-to-all structure, including self-loops.
Remark 3 If G12 0 along with the main diagonal blocks of GLI, then (12) is boundary condition of (7). Thus, Theorem 7 is satisfied and the self-loops are removed from the graph structure.
One consideration when constructing k-contracting Lurie networks is the number of additional parameters required for the parametrisation. A Lurie network, defined in (2), has total parameter count. Nr,
7

---

Published in Transactions on Machine Learning Research (04/2025)
n² + 2nm + + m; whereas k-contracting Lurie network, constructed according to Theorem 5 or Theorem 6, has total parameter count NK I 2n²² + m² The number of parameters only increases significantly when m is large compared to n; however, throughout the literature, many special cases of a Lurie network set n= m (sce $3.1); in which case, Nk becomes marginally smaller than Nz. Furthermore, the all-to-all graph coupling term has Ne = 2(qa)² paramelers or Nc = qn²(q - 1) if the self-loops are removed according to Remark 3. This analysis highlights that ensuring the Lurie network and GLN are k-contracting comes at a minimal computational expense.
5 Empirical Evaluation
In $1, it highlighted that range of convergent dynamics are prevalent in many dynamical systems and neural processes. In the following section we test the impact of the additional expressivity of the Lurie network and the importance of encoding this general convergence property for modelling a range of dynamical systems. The prediction accuracy, generalisation and robustness of the proposed models were used to compare performance. Refer to Richardson et al. (2025) for additional experiments regarding the Lurie network's ability to form representations.
5.1 Datasets and Training
We first consider three time-invariant dynamical systems: (i) an opinion dynamics model of social network where all opinions agree and thus converge unique equilibrium point; (ii) a Hopfield network of associative memory with two stable equilibrium points and one unstable; (iii) generic simple attractor which could be used to model the stored patterns in working memory. Each system has n 3 states allowing for visualisation of the ground truth and predictions. For each dynamical system, we have a dataset including 1, 000 trajectories, sampled every 0.01s over a 20s interval. The test sets were formed by holding out 100 trajectories. The imput to each model was the initial condition sampled from a uniform distribution with the domain (11 for the opinion/Hopficld HHopficld datasets and (-3, +3)3 for the simple attractor. The full trajectory was then used the target to train the model. An illustration of these datasets can be seen in Figure 1 and further details about the data generation can be found in Appendix D.
To test the out of distribution generalisation and robustness, two additional datasets were generated for cach of the opinion, Hopfield and attractor Lasks. These differ from the training datasets in the following ways: (i) they include 100 trajectories over a. 30s interval; (ii) the initial conditions were sampled from a uniform distribution over the intervals 1 < [0)( < 4 for the opinion/ Hopfield datasets and < [0)( < for the opinion/Hopfield
simple attractor, where i E (1,2, These trajectories are 10s longer than the training data with initial conditions also sampled outside the training distribution. To test the robustness, noise sampled from the standard normal distribution was added to the initial conditions of one of these datasets, before generating the trajectories.
For the second set of experiments, we consider two 30 dimensional dynamical systems. The first is a graph- coupled (GC) Hopfield network, formed by connecting 10 previously described Hopfield networks through a graph coupling matrix. To ensure the convergence property was preserved, the matrix was expressed by (12) with GLI sampled from a uniform distribution and G12 = 0. The second system was a graph-coupled attractor, constructed in the same way. The datasets were generated as described in the previous two paragraphs; however, they each included 30, 000 trajectories.
The training settings used are explicitly detailed in Appendix D.. For all models and all datasets, the mean squared error (MSE) loss was used alongside the Adam optimiser. All code was implemented in PyTorc and can be found at attps://githir www./C/RRiicarrsooonnniiiioeeeee waw.iaarrdd//rrreeeeettesss
5.2 k-contracting Lurie Network
In this section, we compare the k-contracting Luric network against five other continuous-time models: (i) the unconstrained Lurie network, for testing the importance of the k-contraction constraints; (ii) an unconstrained neural ODE with two hidden layers, each comprised of 20 neurons and ReLU activations; (iii) the three
8

---

Published in Transactions on Machine Learning Research (04/2025)
constrained continuous-time models detailed in $3.1 where the SVD combo network is comprised of single node. The numerical integration was performed using the Euler method with step size d = 1 x 10-2. Besides the neural ODE, each model had tanh activations, which is slope-restricted with g = 1.
For the opinion and Hopfield datasets, the k-contracting Lurie network was constructed according to Theorem whereas Theorem was used for the simple attractor. We set k 1 for the opinion dynamics, 1 = 2 for the mulit-stable Hopfield network and k 3 for the simple attractor. Similar results were obtained when setting k = 3 for all examples.
Table MSE on a test sot of 100 trajectorics. The lowest MSE presented alongside the mcan and standard deviation caleulated after training each model N = times on single T4 GPU (Google Colab).
Model
MSE (mean std, best)
Opinion
Hopfield
Attractor
L-Luuie Network
3.00
x 10
(1.5 0 0.26) x 10
x 10
Lurie Network
4.0. x 10
(3.6 0.39)
x 10-1
(5.1 x 10
Neural ODE
(2.0 0.433 x 10-4
12.5
x 10-2
1..0 x 10-2
Lipschitz RNN
(3.9+33.0,.88) x 10-2
22.9+2000030)
x 10-1
1.48 2.0, x 10-2
SVD Combo
(306
x 10-4
(2.9
x 10
3.12 2.74
Antisym. RNN
(30.0=0..,,8.0))
10-2
(4.3
x 10-1
+
Table compares the MSE on the test set of each task. The k-contracting Lurie network achieved the best, MSE on two out of three examples. The importance of the k-contraction conditions is particularly clear when comparing the mean and standard deviation with that of the unconstrained Lurie network. These conditions clearly reduce the search space to tractable region to optimise over as the MSE of the unconstrained Lurie network is at least order of magnitude worse than its k-contracting counterpart. The other models perform one would expect: (i) the neural ODE demonstrates strong accuracy across all tasks; (ii) SVD combo performs well on the opinion dataset, where the 1-contraction assumption is valid, but struggles on the others; (iii) the Lipschitz RNN struggles on the attractor dataset whilst the antisymmetric RNN struggles across the board due to the A matrix being fixed at zero and the cigenvalues of the C matrix being fixed to almost purely imaginary values. Figures 3, 6, show a random sample of trajectories from each test set, along with the predictions of each model.
Table 2: MSE on a new test set of 100 trajectories where: (i) the initial conditions were sampled outside the range used for training and the trajectories were 10s longer than the training set; (ii) additionally, noise sampled from the standard normal distribution was added to the initial conditions.
Model
MSE
MSE (noisy inputs)
Opinion
Hopfield
Attractor
Opinion
Hopfield
Attractor
k-Lurie Network
2.9 x 10
5.6 10-2
2.3 10
2.0 10
3.2 x 10-1
1.28
Lurie Network
6.4 x 10-2
1.8 x 10-1
5.96
2.7 x 10-1
4.4 x 10-1
6.79
Neural ODE
1.2 x 10-2
1.09
2.31
3.9 x 10-2
1.63
4.76
Lipschitz RNN
2.3 x 10-1
7.3 x 10-1
6.2 x 10-1
2.9 x 10-1
9.7 x 10-1
1.71
SVD Combo
1.0 x 10-2
2.38
20.9
3.3 x 10-2
7.97
30.30
Antisym. RNN
6.43
5.25
52.1
7.29
6.33
52.9
Table 2 compares the generalisation and robustness of the models each task. No new models were trained, instead the best models from Table 1 were directly applied to these out of distribution and noisy datasets ($5.1). The k-contracting Luric network performs the best on all of these datascts and for some, it still demonstrates a MSE of an order of maguilude lower than the next best model. The MSE of the unconstrained Lurie network and neural ODE tended to drop off for these datasets whereas the MSE of the constrained models, excluding the antisymmetric RNN, tended to stay fairly consistent when their assumptions were valid. Figures 4, 7, 10 show a random sample of foise-free trajectories and predictions for each dataset, whilst. Figures 5, 8, 11 repeat. the same for the noisy inputs. The k-contracting Lurie network predicted the
9

---

Published in Transactions on Machine Learning Research (04/2025)
correct long-term behaviour most accurately under all conditions; even when noise was added, the error predominantly present during the initial transient.
5.3 k-contracting Graph Lurie Network
This section repeats the same experiments as the previous section, but for two 30-dimensional graph-coupled (GC) dynamical systems: the GC Hopfield network and the GC simple attractor ($5.1). The state of these datasets is significantly larger than those used in other dynamical systems datascts such as: (i) the LASA dataset Lemme et al. (2015) where the 2-d trajectories typically stacked to form 4-d or 8-d trajectories; (ii) simulated datasets of the 2, 4 or 8 link pendulums which, respectively, have 4, 8 or 16 dimension trajectories.
For both datasets, the GLN was constructed according to Theorem 7 and Remark 3 where n = m = 3 and q = 10. The individual Lurie networks were constructed according to Theorem for the GC Hopfield network, with k = 2 and Theorem for the GC attractor, with 3. The neural ODE was formed using two layers with 100 neurons and ReLU activations. The SVD combo leveraged similar graph structure to the one used in this paper. The other models were the same in the previous section but with (gn)-dinensional state.
Table 3: MSE test set of 100 trajectories. The lowest MSE is presented alongside the mcan and standard deviation calculated after training each model N = 3 times on single A100 GPU (Google Colab).
Model
MSE (mean + std, best)
GC Hopfield
GC Attractor
GLN
0.016 = 0.0003, 0.016
0.293 0. 091999, 0.015
k-Lurie Network
0.238 + 0.0011, 0.237
0.737 + 0.0108, 0.723
Lurie Network
2.537 + 1.3327, 1.157
291.8 191.84, 21.83
Neural ODE
0.138 = 0.0291, 0.114
3.445 = 0.3626, 2.942
Lipschitz RNN
0.124 = 0.0161, 0.105
0.658 = 0.004, 0.433
SVD Combo
0.339 + 0.1363, 0.229
3.024 + 1.1240, 1.435
Antisym. RNN
0.448 + 0.0023, 0.444
6.386 + 0.0066, 6.380
Table 3 shows the GLN had lower MSE than all other models by factor of 10. Comparing the GLN, k-contracting Lurie network and the unconstrained Lurie network highlights the improvements due to the k-contraction conditions and the graph structure. Table also indicates that the GLN generalised and remained robust to noise, even in these nigh-dimensional systems. The same can be said for the k-contracting Lurie network which achieved the second lowest MSE the generalisation and robustness tests. The inherent. graph structure of the SVD combo may be the reason behind its improved ranking in the GC Hopfield tasks whereas the neural ODE particularly struggled to generalise for the GC attractor. Possible explanations behind the poor performance of the constrained benchmark models are suggested in the next section.
Table 4: MSE on new test set of 100 trajectories where: (i) the initial conditions were sampled outside the range used for training and the trajectories were 10s longer than the training set; (ii) additionally, noise sampled from the standard normal distribution was added to the initial conditions.
Model
MSE
MSE (noisy imputs)
GC Hopfield
GC Attractor
GC Hopfield
GC Attractor
GLN
0.08
1.67
0.26
2.85
k-Lurie Network
0.77
6.10
1.05
6.90
Lurie Network
20350
5638
25481
6368
Neural ODE
2.12
24.93
2.76
25.11
Lipschitz RNN
3.58
6.17
4.25
7.84
SVD Combo
0.84
11.40
1.07
12.30
Antisym. RNN
5.32
50.68
6.21
52.10
10

---

Published in Transactions on Machine Learning Research (04/2025)
6
Related Work
Several constrained continuous-time RNN models exist in the ML. literature. The model structure of three notable examples were presented in $3.1 as they happen to be special cases of a Lurie network, when modelling lime-invariant dynamical systems. The anlisymmetric RNN Chang el al. (2019) and the Lipschitz RNN Erichson et al. (2020) were designed to address the exploding and vanishing gradient problem Pascanu et al. (2013). The antisymmetric RNN did so by parametrising the RNN such that the real cigenvalues of the Jacobian were zero. It achieved this by setting A = 0 and restricting C' to being skew-symmetric. Whilst this does prevent the gradients from exploding and vanishing, it restricts the dynamics which the model can learn to purely oscillatory behaviour. The Lipschitz RNN has more relaxed parametrisation. This model constructs the A and C matrices such that they are both convex combinations of symmetric and skew-symmetric matrices. However, the weight, of the symmetrie matrix can only vary between 0 and 0.5, whereas the weight of the skew-symmetrie matrix can vary between 0.5 and 1. Again, this addresses the vanishing and exploding gradient problem, but the model cannot encode dynamics which are predominately decaying or growing. Finally, the RNN proposed in Kozachkov et al. (2022a) has biological motivations and encodes 1-contracting dynamics. This implies that all possible trajectories will exponentially converge, making the model robust to input disturbances.
Whilst the models above address variety of problems, each model is quite limited in the range of dynamics they can learn, which is a problem when trying to design generalised model for learning time-invariant dynamical systems. With respect to (2), the Lurie network has more flexibility than all of these models. Firstly, it includes all three weight, matrices (A, B, C)) and biases (ba, by); whereas the models mentioned above fix at. least one of these parameters. Secondly, the constraints imposed, and the corresponding parametrisations, allow the model to learn a variety of dynamics including, but not limited to, those mentioned above. The only limitation is that the dynamics must converge in some way. This includes certain types of chaotic systems, for example Thomas' cyclically symmetric attractor Thomas (1999). At the edge of the chaotic regime, the trajectories converge to strange attractor, which could be modelled by 3-cottracting Lurie network.
The Lurie network is also related to a class of feed-forward models named implicil or equilibrium networks. These models use an implicit equation to express the relationship between the model output, layer outputs and model input in a compact vectorised form El Ghaoui et al. (2021). Like the Lurie network, these models can be represented by the interconnection of a linear time-invariant system and a nonlinearity. This makes analysis tools from control theory, such as Lipschitz bounds, applicable to these models Fazlyab et al. (2019). An additional connection is that the solution to the implicit equations correspond to equilibrium points of a Lurie system Revay et al. (2020).
As mentioned in the introduction, the k-contraction constraints used in this paper have an interesting connection to some propertics observed in biological learning systems. A 2-contracting model can replicate the behaviour of associative memory, where every stored paltern corresponds to an equilibrium Kozachkov et al. (2023). Furthermore, a 3-contracting model can replicate the dynamics of working memory, where patterns are retained as attractor states Kozachkov et al. (2022b). Due to the 2-contraction (3-contraction) constraints, the equilibrium points (attractor states) must lie on a line (plane). As a result, the conditions developed in this paper could be of interest. to neuroscientists and MT. rescarchers interested in memory storage and retrieval Ramsauer et al. (2020); Hopfield (1984); Krotov & Hopfield (2020). The formation of the GLN also has connections to biology. For example, constructing larger systems from smaller modules can be motivated by evolutionary biology Simon (1962) where the name facilitated variation Gerhart & Kirschner (2007) is usod to describe the development of traits in response to adaptation of the regulatory clements that connect modules, rather than the core components themselves. This was investigated in Kozachkov et al. (2022a) where the weights of the individual RNNs satisfied 1-contraction condition but were fixed. Only the graph coupling weights were learnt during training.
The relationship between properties guaranteed by k-contraction analysis and those observed in associative and working memory suggest the k-contracting Lurie network possesses a number of appealing properties for an ML model; hence, may be suitable for wider class of ML problems. This proposition supported by the successful application of the coastaty RNN models from $3.1 (special cases of Lurie worw on wide range of ML. tasks. Since the Lurie network is more structured, time-invariant example
11

---

Published in Transactions on Machine Learning Research (04/2025)
of a neural ODE Chen et al. (2018), it will also be applicable to similar array of tasks. Beyond modelling time-invariant dynamical systems, this includes image classification and continuous normalising flows Kidger (2022). The only limiting requirement is that the input must be passed in through the initial condition which in some cases, such as image classification, may require a pro-processing layer.
7 Conclusion
The Lurie network was presented as novel and unifying time-invariant neural ODE. Stability results were presented which verified a generalised form of stability of the Lurie network, namely k-contraction. These results were incorporated in a principled approach for constructing k-contracting Lurie networks and graph Lurie networks which, ullike existing orbitally stable systems, such as those observed associative and working memory. Empirical results showed the benefit of the graph structure and the k-contraction constraints through improved prediction accuracy, out of distribution generalisation and robustness models also compared favourably against well known stability-constrained models, which have been shown to be special cases of the Lurie network, and an unconstrained neural ODE. Future thcoretical work will try to expand the class of systems the k-contracting Lurie network can be optimised over, by obtaining similar results for systems which are k-contracting in a diagonal metric. It would also be interesting to find conditions on the inputs of time-varying Luric network for which the convergence properties of the k-contracting Luric network are also upheld, making the Lurie network applicable to sequential processing tasks. Finally, due to the connections with biological neural processes, would like to empirically investigate the performance of the model wider range of machine learning tasks. The recently proposed working memory benchmark Sikarwar & Zhang (2024) is of particular interest due to its relationship with 3-contracting dynamics.
models, could be trained to model multi-stable and
range of simulated dynamical systems. These same
Acknowledgments
This work was supported in part by the Defence Science and Technology Laboratory (DSTL) and the UK Research and Innovation (UKRI) Centre of Machine Intelligence for Nano-electronie Devices and Systems [FP/S024298/1]-
References
Antoine R Adamantidis, Carolina Gutierrez Herrera, and Thomas Gent. Oscillating circuitries in the sleeping brain. Nature Reviews Neuroscience, 20(12):746-702 2019.
Peter Ashwin, Stephen Coombes, and Rachel Nicks. Mathematical frameworks for oscillatory network dynamics in neuroscience. The Journal of Mathemalical Neuroscience, 6:1-92, 2016.
Eyal Bar-Shalom, Omri Dalin, and Michael Margaliot. Compound matrices in systems and control theory: a tutorial. Mathematics of Control, Signals, and Systems, pp. 1 15, 2023.
Joaquin Carrasco, Matthew C Turner, and William P Heath. Fames-Falb multipliers for absolute stability:
From O'Shea's contribution to convex searches. European Journal of Control, 28:1 19, 2016.
Andreu Cecilia, Samuele Zoboli, Daniele Astolfi, Ulysse Serres, and Vincent Andrieu. Generalized Lyapunov conditions for k-contraction: analysis and feedback design. 2023.
Veronica Centorrino, Francesco Bullo, and Giovanni Russo. Contraction analysis of Hopfield neural networks with Hebbian learning. In 2022 IEFE 61st Conference Decision and Control (CDC), pp. 622 0277 IEFE, 2022.
Bo Chang, Minmin Chen, Eldad Haber, and Ed Н Chi. Antisymmetricrm:: A dynamical system view on recurrent neural networks. arXiv preprint arxiv: Xk::1902.0968. 2019.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
12

---

Published in Transactions on Machine Learning Research (04/2025)
Ross Drummond, Matthew C Turner, and Stephen R Duncan. Reduced-order neural network synthesis with robustness guarantees. IEEF Transactions on Newral Netwoorks and Learning Systems, 2022.
Ross Drummond, Chris Guiver, and Matthew Turner. Exponential input-to-state stability for Lur'e systems via integral quadratic constraints and Zames-Fale multipliers. IMA Journal of Mathematical Control and Information, pp. dnae003, 2024.
Laurent Ghaoui, Fangda Gu, Bertrand Travacca, Armin Askari, and Alicia Tsai. Implicit deep learning.
SIAM Journal on Mathematics of Data Science, 3(3):930 958, 2021.
N Benjamin Erichson, Omri Azencot, Alejandro Queiruga, Liam Hodgkinson, and Michael W Mahoncy. Lipschitz recurrent neural networks. arXiv preprint ar xrx:22006.1 12070, 2020.
Mahyar Fazlyab, Alexander Robey, Hamed Hassani, Manfred Morari, and George Pappas. Efficient and accurate estimation of lipschitz constants for deep neural networks. Advances in neural information processing systems, 32, 2019.
Mahyar Fazlyab, Manfred Morari, and George J Pappas. Safety verification and robustness analysis of neural networks via quadratic constraints and semidefinite programming. IEEE Transactions on Automatic Control, 67(1):1 15, 2020.
John Gerhart and Mare Kirschner. The theory of facilitated variation. Proceedings of the National Academy of Sciences, 104(suppl 1):8582 8589, 2007.
Albert Gu, Tri Dao, Stefano Ermon, Atri Rudra, and Christopher Ré. Hippo: Recurrent memory with optimal polynomial projections. Advances in neural information processing systems, 33:1474- 1487, 2020.
Albert Gu, Karan Goel, and Christopher Ré. Efficiently modelling long sequences with structured state spaces. arxie preprint orXiv:2111 11.00996, 2021.
William L Hamilton, Rex Ying, and Jure Leskovee. Representation learning on graphs: Methods and applications. arXiv preprint orXiv: 701.17990000 709.05584. 2017.
John J Hopfield. Neural networks and physical systems with emergent collective computational abilities. Proceedings of the national academy of sciences, 79(3):23544255x 1982.
John J Hopfield. Neurons with graded response have collective computational properties like those of two-state neurons. Proceedings of the national academy of sciences, (11)0)33088-300 1984.
Roger Horn and Charles R Johnson. Topics in matrix analysis. Cambridge university press, 1994.
Eric R Kandel, James H Schwartz, Thomas M Jessell, Steven Siegelbaum, A James Hudspeth, Sarah Mack, et al. Principles of neural science, volume 4. McGraw-hill New York, 2000.
Hassan K Khalil. Nonlinear systems. Patsence Hall, 115, 2002.
Mikail Khona and Ila R Fiele. Attractor and integrator networks the brain. Nature Reviews Neuroscience, 23( 23(1))7744-766 2022.
Patrick Kidger. On neural differential equations. arXiv preprint arXiv:2202.024 02435, 2022.
Leo Kozachkov, Mikael Lundqvist, Jean-Jacques Slotine, and Earl K Miller, Achieving stable dynamics in neural circuits. PLoS computational biology, 16(8):e1007659, 2020.
Leo Kozachkov, Michaela Ennis, and Jean-Jacques Slotine. RNNs of RNNs: Recursive construction of stable assemblies of recurrent neural networks. Advances in Neural Information Processing Systems, 35: 30512-30527, 2022a.
Leo Kozachkov, John Tauber, Mikacl Lundqvist, Scott Brincat, Jean-Jacques Slotine, and Earl K Miller. Robust and brain-like working memory through short-term synaptic plasticity. PLOS Computational Biology, 18(12):e1010776, 2022b.
13

---

Published in Transactions on Machine Learning Research (04/2025)
Leo Kozachkov, Jean-Jacques Slotine, and Dmitry Krotov. Neuron-astrocyte associative memory. arXiv preprinl 0rXii:2311.088 1.01335, 2023.
Dmitry Krotov and John Hopfield. Large associative memory problem in neurobiology and machine learning. arXiv preprint Xrx:2:2008, 06996, 2020.
Samuel Lanthaler, T Konstantin Rusch, and Siddhartha Mishra. Neural oscillators are universal. Advances in Neural Information Processing Systems, 36, 2024.
Andre Lemme, Yaron Meirovitch, M Khansari- Zadeh, Tamar Flash, Aude Billard, and Jochen Steil. Open- source benchmarking for learned reaching motion generation in robotics. Paladyn, Journal of Behavioral Robotics, 611:0000001155520150000 2015.
Laurent Lessard, Benjamin Recht, and Andrew Packard. Analysis and design of optimization algorithms via integral quadratic constraints. SIAM Journal on Oplimization, 26(1):57 95, 2016.
Mario Lezcano-Casado and David Martmez-Rubio. Cheap orthogonal constraints in neural networks: A simple parametrization of the orthogonal and unitary group. In International Conference on Machine Learning, pp. 3794 -803. PMT.R, 2019.
Winfried Lohmiller and Jean-Jacques E Slotine. On contraction analysis for non-linear systems. Aulomalicn, 34(6):683-00 (66666663-69,,, 1998.
James S Muldowney. Compound matrices and ordinary differential equations. The Rocky Mountain Journal of Mathematics, 357-872, 1990.
Ron Ofir, Jean-Jacques Slotine, and Michael Margaliot. k-contraction in generalized Lurie system. arXiv preprint arXiv:2309.0754, 2023.
Ron Ofir, Alexander Ovseevich, and Michael Margaliot. Contraction and k-contraction in Lurie systems with applications to networked systems. Automatica, 159:111341, 2024.
Matthijs Pals, Jakob H Macke, and Omri Barak. Trained recurrent neural networks develop phaso-locked
limit cycles in working memory task. PLOS Computational Biology, 20(2):e1011852, 2024.
Poogyeon Park. A revisited Popov criterion for nonlinear Lur'e systems with sector-restriction.. International Journal of Control, 68(3):461 470, 1997.
PooGyeon Park. Stability criteria of sector-and slope-restricted Lur'e systems. IEEE Transactions on Automatic Control, (722)308-3333 2002.
Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. On the difficulty of training recurrent neural networks. In International conference on machine learning, pp. 1310 1318. Pmlr, 2013.
Patricia Pauli, Dennis Gramlich, Julian Berberich, and Frank Allgöwer. Linear systems with neural network nonlinearities: Improved stability analysis via acausal Zames-Fall Falb multipliers. In 2021 60th IEEE Conference on Decision and Control (CDC), pp. 3611-3618. IEEE, 2021.
Kaare Brandt Petersen, Michael Syskind Pedersen, et al. The matrix cookbook. Technácal University of Denmark, 7(15):510, 2008.
Alex T Piet, Jeffrey C Erlich, Charles D Kopec, and Carlos D Brody. Rat prefrontal cortex inactivations during decision making explained by bistable attractor dynamics. Neural computation, 29(11))2861-266 2017.
Hubert Ramsauer, Bernhard Schäfl, Johannes Lehner, Philipp Seidl, Michael Widrich, Thomas Adler, Lukas Gruber, Markus Holzleitner, Milena Pavlović, Geir Kjetil Sandve, et al. Hopfield networks is all you need. arXiv preprint arXtv:2008.0220 2020.
Max Revay, Ruigang Wang, and Ian Manchester. Lipschitz bounded equilibrium networks, arXiv preprint arXiv:2010.0173: 2020.
14

---

Published in Transactions on Machine Learning Research (04/2025)
Carl Richardson, Matthew Turner, Steve Gunn, and Ross Drummond. Strengthened stability analysis of discrete-time Lurie systems involving ReLU neural networks. In 6th Annual Learning for Dynamics & Control Conference, pp. 209 221. PMLR, 2024.
Carl R Richardson, Matthew C Turner, and Steve R Gunn. Strengthened Circle and Popov Criteria for the stability analysis of feedback systems with ReLU neural networks. IEEE Control Systems Letters, 2023.
Carl R Richardson, Matthew C. Turner, and Steve Gunn. Lurie networks with k-contracting dynamics. wtt.pporrrreeeee..oooorommr New Frontiers in Associative Memories Workshop at ICLR, 2025. URL https://openrevieet net/forum?
id=RaAYeCxj1u.
T Konstantin Rusch, Ben Chamberlain, James Rowbottom, Siddhartha Mishra, and Michael Bronstein. Graph-coupled oscillator networks. In International Conference on Machine Learning, pp. 18888-18909. PMLR, 2022.
Sugandha Sharma, Sarthak Chandra, and la Fiete. Content addressable memory without catastrophic forgetting by heteroassociation with fixed scaffold. In International Conference on Machine Learning, pp- 19658-196822 PMLR, 2022.
Ankur Sikarwar and Mengmi Zhang. Decoding the enigma:
benchmarking humans and Als on the many facets of working memory. Advances in Neural Information Processing Systems, 36, 2024.
Herbert A Simon. The architecture of complexity. Proceedings of the American philosophical society, 106(6): 467-482, 1962.
René Thomas, Deterministic chaos seen in terms of feedback circuits: Analysis, synthesis, labyrinth chaos. International Journal of Bifureation and Chaos, 9(10):1889 1905, 1999.
Matthew C Turner and Ross Drummond. Analysis of MIMO Lurie systems with slope restricted nonlincarities using concepts of external positivity. In 2019 IEEE 58th Conference on Decision and Control (CDC), pp. 163-168. IEEE, 2019.
Mathukumalli Vidyasagar. Nonlinear systems analysis. SIAM, 2002.
Ryan Vogt, Maximilian Puelma Touzel, Shlizerman, and Guillaume Lajoie. On Lyapunov exponents for RNNs: Understanding information propagation using dynamical systems tools. Frontiers in Applied Mathematics and Statistics, 8:818799, 2022.
Hugh R Wilson and Jack D Cowan. Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical journal, 12(1):1 24, 1972.
Chengshuai WI,, Tlya Kanevskiy, and Michael Margaliot. k-contraction: Theory and applications. Automatico, 136:110048, 2022.
George Zames and PL Falb. Stability conditions for systems with monotone and slope-restricted nonlinearities. SIAM Journal on Control, 6(1):89-108, 1968.
Xiaojiao Zhang and Baotong Cui. Synchronization of Lurie system based on contraction analysis. Applied Muthemalics and Compulation, 223:180 190, 2013.
Samuele Zoboli, Andreu Cecilia, and Sophie Tarbouriech. Quadratic abstractions for k-contraction. 2024.
15

---

Published in Transactions on Machine Learning Research (04/2025)
A Extended Preliminaries
As many of the tools used in this paper are not well-known in the machine learning community, a condensed background is presented here, based on results from Bar-Shalom et al. (2023); Wu et al. (2022). The first section is on compound matrices, an important algebraic tool needed for generalising contraction analysis k-contraction analysis. Following that, geometric interpretation of the k-compound matrix is stated through its relationship to the volume of parametrised set. Finally, the results in these two sections are utilised to define and provide intuition for the convergence of k-contracting dynamics.
A.1 Compound Matrices
In this section, we document several known definitions and algebraic results related to compound matrices. The results are included without proof; the interested reader should refer to Bar-Shalom et al. (2023) for a more detailed tutorial on the topic.
Let n be a positive integer and fix & E [1, n]. The ordered set of incrensing sequences of k integers from [1, is denoted by Q(k, n). For example: Q(3, 4) = [11, (1, 4), (1, 4), (2)
Now consider matrix sub-matrix obtained by taking the entries of W along the indexed by and columns indexed by As example, if and (3(1, 4)] would then be given by
R For E Q(k, n) and B E Q(k, m), the matrix | denotes the k
4, then Q(2, 4) {(1,2), (1, (1, 4), (2, 3), (2, 4), (3, 4)}. The sub-matrix
W13
W14
W[( (3, 4)]
1223
W24
The k-minors of the matrix W are defined as (( = det(W(a)))).
Definition 1 (I-multiplicative compound) Let W E Raxm and fúr k E mi,min m)). The k- multiplicative compound of W, denoted W(k), the (x) x (7) matria containing all the k-minors of W ordered lezicographieally
For example, if we have n = m and then a, 3 E Q(2,3) = (1,3, and
(11,))
W (11,))
(((L,
(
(111,3) (113)) (11)
W (12.))
((2233) ((2.3)
Some important special cases include
( W
((n) = det(W)
(k) = =
W E D" + E D
((k)
(13)
with s = (w). Next, we present a series of algebraic results concerned with the k-multiplicative compound.
Fact 1 (Cauchy-Binet Formula) U E Raxm, V E pmxx and E [1, min(n, , m, p)), then (UV)(M) (())
Fact. 2 Fix h E [1, mi). As a consequence of Definition 1, if W e paxx then (WT)(*) (W(K))T
Fact 3 Fix ls E [I. If W e Paxn is non-singular, then by Theorem
(W (()) = -W(
Fact 4 Fix l E [1, min(n, (,,nnn,,,,p) p)]. If W e R U = Roxa and V E Raxp, then by Theorem 1 (UWV)( (14)
16

---

Published in Transactions on Machine Learning Research (04/2025)
Fact 5 Fix k E [1,m]. An implication of Theorem 1 is that if W R with eigenvalues A An then the eigenvalues of W(k) are the (x) products
II Ar 1:. < < mn
We now introduce the definition of second compound matrix, the k-additive compound, and set of algebraic results related to it.
Definition 2 (k-additive compound) Let W E Aaxx and k E [I. The k-additive compound of W is the matrix defined by
d (h2 + (eW)
W
de
Special cases include
W W
Whn = tr(W)
(pT,) kpI,
We WED" D - wll E
(14)
with 8 = (2)) Like before, we now present some useful algebraic results related to the k-additive compound.
Fact 6 If W e gaax and l E [1, then as a consequence of Definition 2
(Wil))T =(WIl))
Fact 7 Pir k E [], For W E xaxx with eigenvalues A Xn, the eigenvalues of Wil are the (
sums
k
Au 15 < ... < mn]
An important consequence of Fact 7 is that if positive definite (semi-definite), then this property is upheld by WIM), Opposite conclusions can be drawn if W is negative definite (scmi-defiitte.
Fact 8 Fix k E [1.m]. If U.U. V. E Rnxn then
(U + VV = (/(k] + V
Fact 9 Pix ls E [1, miniinn, p)). W E Raxn, U Roxn, V Raax and UV Ip, then
(UWV) wik)
A.2 Volume of kssets
This section aims to provide a clear geometric interpretation of the k-multiplicative compound. The section begins by defining a k-set (the codomain of function dependent variables) before presenting Theorem 8, the key result. which exposes the relationship between the volume of k-set and the k-multipliattive compound of the Jacobian. A k-parallclotope is then shown results, are presented without proof. Refer to Wu el al. (2022) for more information.
example. Like before, these cxisting
Definition 3 (k-sets) Consider compact set D c. R and a continuous differentiable map D R" with n]. The codomain of V is given by the porametrised set
V(D) =
Rd
(15)
Since D is compact and V()) is continuous, (DD) is closed set.
17

---

Published in Transactions on Machine Learning Research (04/2025)
4(D;
4:D 3p
X3
Figure 2: The 3-parallelotope with vertices 21, 22, my E R parametrised by the unit cube, D.
Theorem (Volume of k-sets) Fir & E [], Consider a compact set D c R and a continuously differentiable map & D + R². The volume of the parametrised set (15) is given by
vol (I( S 1 (r)
where (r)
(((r)
@v(r) Dre
is the Jacobian of V. Note that as Jo D + Rnnk implies J (r)
8r1
D + R with s (+)
Theorem 8 states that the volume of k-set is governed by the k-multiplicative compound of the Jacobian. An important geometrical feature, is that for & E (1, the volume of the k-set is equivalent the standard notions of length, area and volume. We now consider the k-parallelotope as an example of a k-set.
Definition 4 (k-parallelotopp) Fix l E n] and let vectors 2,,..., Ia E R². The parallelotope generated by these veclors (and the zero verlex) is the set given by
P(a
Based on Definition 4, the k-parallelotope is k-set with the following compact domain D and continuous differentiable function ((r; ,.., TR), as illustrated for k = 3 in Figure 2.
D D: fre R r 10 [1,
(r
Σ naa
The Jacobian of the k-parallelotope is (r(r)= X and from Theorem 8, the volume of the k-parallelotope is given by [
A.3 k-contraction Analysis
Consider the timo-varying nonlincar system
5= s(1,))
(16)
where f R4 R n It is assumed throughout that continuously differentiable w.r.t x. Fix [1,n] and let & denote the unit simplex.
Sh
(r R ra 2 0 and ri + m < 1}
The convex combination of set of initial conditions .i,..., 2b+1 E R is defined by h S + Rd
h(1 Tk+1)
raa; + (1 Σ
Ik+1
i=1
18

---

Published in Transactions on Machine Learning Research (04/2025)
The set h(S*) is k-set and can be thought of as k-dimensional body of states representing initial conditions of (16). We now define w;((,r) r) measure of the sensitivity of solution to (16) at time i, to a change in the initial condition h(r), caused by a change in Ti.
h(r)) where w;(0,r) 0h(rr 2; -k41 for all i E [1,
wi (;((,r) = Ör;
Ör:
We are now ready to present the definition of k-contraction, followed by its geometric interpretation.
Definition 5 (k-contraction) Fix & E [L.n]. The nonlinear system (16) is k-contracting if there exists an 1) > O and a vector norm such that for any 21.-.., E R and any E S the mapping W R4 S Raxk defined by (t, = [w](1,r)
wx(1,r) satisfies
(t,r) S (H VEER
To explain the geometric meaning of this definition, pick domain D C S* and recall that h(D) is k-set representing k-dimensional bodies of initial conditions for (16); thus, 2 æt,,, ((D)) == (z(t, h(r)) :r E D) is a k-set describing how k-dimensional bodies evolve over time. We now leverage Theorem 8, to show how the volume of these bodies evolves over time when governed by (16).
vol (=(t,h(D)) / Jz (t, ( dr
/ (a(())(r) Dr
aa(e.Mor)) (k)
dr
Ore
When (16) is k-contracting, the volume of these bodies is upper bounded by the initial volume scaled by an exponentially decaying term.
vol (e((,((D)) s sxp(-ntt) / [w((k) r))rr
- = exp(-nt) [(2] xx+1) (#x- xx+1) (K) So dr
Therefore, k-contraction of (16) implies the volume of k-dimensional bodies xtt, ((D) converges to zero at an exponential rate. This can also be interpreted as the volume of k-dimensional bodies is contracting or converging to a (k: - 1)-dimensionall subspace. Figure 1 in the main paper gives an illustration of (1, contracting systems where the 1-contracting system converges to point, the 2-contracting system converges to line and the 3-contracting system converges to plane.
Many existing k-contraction results, including Theorem 1, are expressed in terms of matrix measures. An overview of their definitions and properties can be found in (Vidyasagar, 2002, Section 2.2). Theorem provides a sufficient condition for verifying k-contraction in the 2-norm w.r.t a metric P. The 2-norm was chosen for this work due to its relationship with the eigenvalues of its argument, but other norms could be chosen. Furthermore, one may apply an invertible linear transformation θ to w;(t, and the k-contraction analysis may be performed in this new domain whilst implying the same property holds in the original domain. This idea is made clear in Lohmiller & Slotine (1998) and translates analogously to the k-contraction case. When using such an invertible transformation, the system is said to be k-contracting w.r.t the metric P = OT.
19

---

Published in Transactions on Machine Learning Research (04/2025)
B Lurie Network
B.1 Relationship Between Lurie Networks and Deep Feedforward Models
Consider the vector field å = f(z) defined by the following L-layer feedforward network
1 = WL + yy
wr-1 W Wil-1011 ((r,,2 + bz,-1
ur-2 W 2 + br-
(17)
U2 = W20(u1) + b2 d1 W + bj
To illustrate the superior expressivity of the Lurie network, we wish to show that a special case of (2) can approximate the deep feodforward network (17). An alternative expression for (17) is
0z + (uLLL by
eüg- - + 1-1 ( bb-1
eŭr-2 - + W 2-20 - bb-2
(18)
eig -12 W2 b2
W bj
where - O. Defining new state [z WL-1 UL-2 Uy U2] and output vector y
JuL-1 UL-2 U2 U1] it is clear that (18) is a special case of the Lurie network (2) with state-space matricos and biases defined by the sparse structures below.
0
0
WI
0
0
I
0
0
WL-1
0
A
0
0
0
B
0
0
0
W2.
0
0
0
b1
0
0
100
1
0
C
0
0
0
0
0
b =
bb-2
by
W
0 0
This is just one realisation of a Lurie network which approximates a feedforward network. Other permutations of the state would result in different realisations of the Luric networks weights and biases. Finally, due the division by e, it would not be possible train Lurie network to have the exact form of feedforward werw however, this analysis shows that it is possible approximate the strucuture of deep feedforward network with Lurie network.
20

---

Published in Transactions on Machine Learning Research (04/2025)
B.2 Further Examples
This example is from the graph ML literature Lanthaler et al. (2024). The state of the general neural oscillator is governed by a second order ODE; however, it's equivalent first order representation takes the form (2) with onc possible realisation given by C21 E R by = 0 and
Neural Oscillators:
0 0 B = 0 1
C 0 C21
=
by = by2
The solution then passed through an affine readout layer.
Graph Coupled Oscillators: Another example of second order ODE from the graph ML literature matrix, but this can simply be recast in vectorised form which relates to (2) if a linear coupling function is chosen. This requires the weights to have block matrix form, where one possible realisation is defined by C21 E Rnxn, d, = by 0 and
Rusch et al. (2022). The state is defined
0
I -al
0 0 T
0
0 C21
0
=
B
C=
-77
When the external input is replaced by nonlinear output feedback (i.e., a(t) (v)) and D = 0, the linear state-space layer used in S4 Gu et al. (2021) and Hippo Gu et al. (2020) is related to (2) with A being a lower triangular Hippo matrix, B E Raxm C E R and by = by = 0.
LSSM:
21

---

Published in Transactions on Machine Learning Research (04/2025)
C Proofs
C.1 Proof of Theorem 2
We aim to verify Theorem 1 for the particular case where the nonlinear system is described by the Lurie
network (2). Our proof begins with Theorem 9, which restates (Ofir et al., 2024, Theorem 1). This result is sufficient to satisfy Theorem for systems of the form (19).
Theorem Fin [1, n] and consider the system below.
a = Ax(t)) (t,
Ca
(19)
If there exists mh 12 0 and an invertible O E Rn such that
p(k) A + (A p(t) + ( ((B (00T)))))00C000101000000 ((k) - P(k)
(20)
and
[+] I- V R4 and R.
(21)
where = = then (19) is k-contracting in the 2-norm w.r.t the metric P == θ O.
We first need to express the Lurie network in the form (19). By (3) there exists y < 0 satisfying
k
0< <2 < of and g² Σ (B)) < ~
(22)
Using 7, we can express (2) in the form (19) through the definitions below, where the dependence on has been dropped from V
A A
B = Inn
: = In
(xx) -B=---BBBCaar + by) - --= bz
(23)
The next step is to verify (20). Subbing (23) into the left hand side of (20) and assuming θ T results in the first equality. Setting P: := pln with p > 0 results in the second. Now must leverage some of the facts presented in $A.1. Using the relevant special cases from (13) and (14) leads to equality three and consequently applying Fact 6 and Fact 8 results in equality four. Re-applying (14) and Fact results in the final equality.
p(()) A + ( p(k) + (() ((7*P) (P-1))k) ((--)))))) )(x)
(R),) A + (AAM)
I,) ((²2Tn)) + (1)) (p² (a)
=p* (AlK) + (A +??
=pA((A+A + +k(?"pp-"))
= = (1"p+r [K]
If the matrix above negative definite, then (20) satisfied for some suitably chosen m 0. This is true when
(A+A + AT + ( x
By Fact 7, the inequality above can equivalently expressed as a condition on the sum the largest eigenvalues the matrix inside the k-compound operator. Leveraging (Potersen al., 2008, Eq. 285) allows separate p from the eigenvalues the symmetric component of A, resulting in the equality below.
A + -02p+p-111 ))))) -
Σ
X + AT) <0
22

---

Published in Transactions on Machine Learning Research (04/2025)
By the definition of ag in Theorem this simplifies to
r²p² + 2app + <<0
For salisfying (22), the quadratic inequality always emits al least one solution p -0k -1
The final step is to verify (21). The Jacobian of v as defined in (23), is
Jo(x) = - -B--
For θ = =p In and the definitions from (23), the left hand side of (21) reduces to
= (p J Jy Jv -oo-pp1110 I,) [+]
If the matrix above is negative definite, then (21) is satisfied for some suitably chosen 7/2 0. Repeating the same steps above, this negative definite requirement reduces to the incquality below.
1 J Ju) - kp 1=1
oJ( kp <0
1=1
Subbing in the definition of Jy and applying the well-known property of singular values (Horn & Johnson, 1994, Theorem 3.3.14), then (21) true if
Σ [d((B))((J))))(10
By the assumption made the slope of o, this inequality will always be satisfied if (3) holds.
C.2 Proof of Theorem 3
For this proof, aim to directly verify Theorem 1 for θ E D where f represents the Lurie network (2). We start by substituting the Jacobian of the Lurie network into the left hand side of (1), followed by the application of Fact 8 to obtain the second equality. The subadditivity property of the matrix measure μ2 was then leveraged to split the terms (Vidyasagar, 2002, Section 2.2). As the second term difficult to manipulate, we rely on the simplifying assumption C= B-1 in order to apply Fact 9. As O, B, Jo E D", so are both of their k-compound counterparts (13) (14), which means the θ, terms cancel out. We then apply the property (²(+) S || lla (Vidyasagar, 2002, Theorem 16) and Fact 7, which allows us to leverage the slope restricted assumption on 4. By the relevant special case of the k-additive compound (14), can directly calculate the 2-norm. Finally, kg can be incorporated in #22 as shown in the final inequality.
22 ()
(t,x)) (2,00 ((A+ (J
= 12,000 (A + (BJ+()))))
12,00 (A 12,0 (BB)
12,00 (Alk)) (2(( BB(k) J!h B @ -(1))
12,000 (A (x)
S 12,000 (A (1( 1M11
12,000 (Alk)) + kg
(k(((()) 1 (klle(()) lg
If the final inequality is negative, then Theorem 1 will be satisfied for some suitably chosen n. This is equivalent to the matrix incquality below.
(0(0)) - (A) 22gln)
2
Multiplying on the left. by (k) and the right by (*) results in (4).
23

---

Published in Transactions on Machine Learning Research (04/2025)
C.3 Proof of Theorem 4
The aim of this proof is lo verily Theorem 1 when S' is the GLN (6). We begin by expressing the Jacobian sum of the Jacobian of the independent Lurie networks, Jnddep: and Jacobian of the coupling term Jcouple
Jy(x) (1ddep()) + (x)
where Jindap (x) = AG + Ca and Jenupte(x) = L. Subbing y into the left. hand side of (1) and applying the subadditivity property μ2 results in
12,0000 (J S (2,0 (Iidepp) 12, (Jcouple)
(J
As we assume the q independent Lurie networks are k-contracting in the 2-norm w.r.t the metric P blockdiag(P1,) Pa), know that (2,000 (J indep)) < O. Under this assumption, Theorem is satisfied if
12,000 (Jcoupie) 0
This is equivalent to the matrix incquality below, where Jeouple has been subbed in.
0 (k) r, 1. 0 O (+) (x)1))))))) 30
2
Multiplying on the left by (+(() and on the right by (()) results in (7).
C.4 Proof of Theorem 5 and Theorem 6
The aim of this proof is express the weights of the Lurie network (2) such that Theorem 2 is al always satisfied. More formally, this requires A, B, C (g k) to always hold. As both theorems share the same proof until the final step, where Σa is defined, they have been combined into one proof.
To expose the singular values of B,C, we leverage the singular value decomposition, as in (9b) and (10a). This requires the matrices UB, Uc, VB, Ve to be orthogonal. We can immodiately use the unconstrained parameurisation of the orthogonal class from dezcano-Casade matrices unconstrained symmetric matrices. The matrices ΣB, Ec contain the singular values of the respective matrix, hence ΣB Dmm and Eg DDn" We also treat these unconstrained sets since any element can be obtained by taking the absolute value of an unconstrained diagonal matrix with the same shape.
Martmez--ubio (2019) to express these
To verify Theorem 2, we can combine 0g < 0 and (3) into one inequality reprosenting the intersection of the two sets.
MAA 1 (T) -2k
2 zk where 2₁ : g² Σ (B))
(24)
Thanks to the definition of B and C, the right hand side is a function of the hyperparameters g, and elements of the parameters En Ee, so can be easily computed using sort and sum functions.
To impose this constraint directly on the eigenvalues of the symmetric component of A, we express A as a sum of symmetric and skew-symmetric matrices. The skew-symmetric matrix is unconstrained, this can be left as it is. Finally, we apply the eigenvalue decomposition ion of symmetric matrix to obtain (9a) (10a). This expression allows us to directly place the constraint above on the diagonal matrix ΣA-
This is where Theorem 5 and Theorem differ. Defining ΣA as in (9a) ensures
for all [1,m] A(A+AT) -2 2k k
which guarantees both conditions Theorem 2 will be satisfied; however, conservatism is introduced as all the eigenvalues must be negative. Theorem 6 was established to address this issue. Defining EA in (10b)
24

---

Published in Transactions on Machine Learning Research (04/2025)
guarantees both conditions of Theorem 2 will be satisfied via (24). The definition of ΣA is split into racondt block for the first (k - 1)-eigenvalues, a block for (A + AT) which is defined to ensure (24) holds, and finally a block for the remaining eigenvalues which must be defined to ensure the k eigenvalues involved in (24) are the largest.
C.5 Intersection Proof
The set k2 is defined by (8a). For h: E n], the set of (k - 1)-contracting systems is denoted by
D2(g, k - 1)= (A. 04-1 <0, zk-1 a 1)}
Theorem 10 For E (2,3), the sets l2zgg k - 1) and l22(g, k) intersect when
1
k-
1 k
g (B) <
\ +
k ax
ak
(25)
4k
Proof:
We first note the following useful relationships, where A ((A+AA))
1 k k
ag
A
a
= o') + #k-1
(26)
2k
From (26) it is clear that a < 0 is true when
A < -2k - 1)ax-1
A subset these solutions occurs when the first condition l2(g. k - 1) holds. That is ag < 0.
From (26) it is also clear that 2μ < afkk true when
2k-1 < a (o)
If ak -1) 1) aakk- (C) then a subset of these solutions occurs when the second condition of O2(g, k: - 1) holds. After rearranging, this is equivalent to
g² (o aat 1
The RHS can equivalently be expressed by a quadratic equation in d₁
k
2
a - 2 1
11
ds
[kk-
4:
02 (k-
2h:
b:
1
X2
+
1
Okk 124 +
(k 11)
a
2
k 16-00 2 a (k - 1)
4122
k2
k2
1
k
(k 1k-1) 1)2
X
+
a
aq
1d
ak-
2
_ (k- 1)
2
4k
k
k
2
+
k --1
1dx
h k-
Ok-1
2
4k
k
k
QE ak-))
Therefore, l2(g. & 1) intersects with (g k) if (25) holds.
25

---

Published in Transactions on Machine Learning Research (04/2025)
Finding the roots of Q& (Ax, ak-1) will allow to check what intervals of d the Qx(RR ak-1) 0. This is required since the LHS of (25) must be positive.
(k: 1
h
b = 2(1 - k)ck-1 1 + 2k
1:2
k-1
1:2
0
k-1
= 2(1 k)ax-1 + 2/ k(h - 1)ax-1
2 - k + likk - 11) cs-1
The roots for and are given by
X = (-2
X/ (-4 2v6)a²
Given 0x-1 < O, will test the intervals A2 (-2,, and A (-4t to see if the co-domain of Qa (AR, 0k-1) is posilive. These intervals ensure d < 0 (a requirement for or < and 0r-1 < 0).
Using x2 = 01 and À = a2 as the test points, we see that
1
1 2
Q2(a1;a))
goi
+
01
2
01
(27)
2
2
8 >
0
(28)
2
2 3 02
Q2(02;02)
02
+
02
(29)
12
3
02
2
0
(30)
As a result, for k E [2,3), 3], the sets 02(g, k - 1) and [2(g, k) intersect if 0 E 2(- - [- /kkk 1k-kk-11 1))ax-1 1
and (²o < (A 11- If (A, B, C) E 02(g, 2) then X2 must be in this interval since A2 S d By the same logic, if (, E l2 3) then A must be in this interval since λ I X S A
Hence, the requirements for intersection can be simplified just (25).
C.6 Proof of Theorem 7
This proof aims show that Theorem is satisfied when L defined as in (12). First, multiply (7) on the left and right by O- (k), Sequentially applying Fact 6, Fact and Fact results in
Subbing in the definition of L from (12) and recalling that P = θ θ leads to
(GL2 + TTa) 30
which is assumed to hold.
26

---

Published in Transactions on Machine Learning Research (04/2025)
D Extended Empirical Evaluation
D.1
Data
The data was synthctically generated by numerically integrating over the analytical models of cach dynamical system. The integration was performed using the Euler method with step size = x 10-2
D.1.1
Opinion Dynamics
This model was presented in Ofir et al. (2024). It has the following state space equations
a 1.51, + (.5)(Cz) +
where
+0.2
C
0.2
and the tanh function is applied clement-wise. at b.
The model is 1-contracting and has unique equilibrium point
D.1.2 Hopfield Network
This model is variation of the Hopfield network presented in Ofir et al. (2024). It has the following state space equations
i =
Bb(x)
where
1
1
17
B
1
1
1
1
1
and the tanh function is also applied elementtwise. The model is 2-contracting and has two stable equilibrium
points: e,
[0.79 0.79 0.79 e2 -ei; and an unstable equilibrium point e3 0.
D.1.3 Simple Attractor
This model was presented Cecilia et al. (2023). It has the following state space equations
t Ax + (C((x))
where
0
-2
0
0
0
0
A =
-1
0
-1
B
0
0
0
C
0.5
0
0.5
0.5
0
0
and (0=) 33 is the nonlinearity applied element-wise This function is not slope-restricted, so the simple attractor does not satisfy the assumptions of the Lurie network. The model is 3-contracting and has several attractor states.
27

---

Published in Transactions on Machine Learning Research (04/2025)
D.2
Training
The default training setlings common to the opinion, Hopfield and attractor dalasets are presented in Table 5. The only parameters which varied between models were the learning rate and epoch which it was cut. Deviations from the default settings are detailed in Table 6. These values were chosen based on observations during training. No hyperparameter sweep was performed. The same details are also presented for the graph-coupled Hopfield network and graph-coupled attractor datasets Table and Table When training the models for the opinion, Hopfield and attractor datasets, a single T4 GPU (accessed through Google Colab) used. A single A100 GPU (also accessed through Google Colab) used for training the models the graph-coupled Hopfield and graph-coupled attractor datasets.
Table 5: Default training set settings for opinion, Hopfield and attractor datasets.
Parameter
Valuc
Batches
10
Batch size
100
Test split
0.1
Epochs
100
Criterion
Mean squared error
Optimiser
Adam (default scttings)
Learning rate (LR.)
1 10-2 (no docay)
Table 6: Deviations from the default settings for the opinion, Hopfield and attractor datasets. LR decay inchudes the size of the docay followed by the epoch the decay was made.
Dataset
Model
LR
LR Decay
Opinion Dynamics
Lurie Network
5 X 10-3
Opinion Dynamics
Antisymmetric RNN
5 x 10-3
Hopfield Network
k-Lurie Network
x 10 10
Hopfield Network
Neural ODE
10
0.1, 75
Hopfield Network
Antisymmetric RNN
10
Simple Attractor
Neural ODE
10-3
Simple Attractor
Antisymmetric RNN
5 10-3
Table 7: Default training settings for the graph-coupled Hopfield and attractor datasets.
Parameter
Value
Batches
15
Batch size
2000
Test split
Epochs
100
Criterion
Mean squared error
Optimiser
Adam (default settings)
Learning rate (LR.)
5 x 10-3 (no cuts)
28

---

Published in Transactions on Machine Learning Research (04/2025)
Table 8: Deviations from the default settings for the graph-coupled Hopfilld and attractor datasets. LR decay includes the size of the decay followed by the epoch the decay made.
Dataset
Model
LR
LR decay
GC Hopfield Network
Lipschitz RNN
5 5 10
0.2, 60
GC Hopfield Network
GLN
3 10
3 60
GC Simple Attractor
GLN
1 10
GC Simple Attractor
k-Lariee Network
1 10-2
GC Simple Attractor
Neural ODE
5 5 10
0.5, 40
D.3
Further Results
Ground Truth
k-Lurie Net.
Lurie Net.
Neural ODE
0.5
0.5
0.5
0.5 X3
0.0
0.0
0.0
0.0
-0.5
0.5
005
-0.5
1.0
0.5
0.5
0.5
0.5
-0.5
0.0
-0.5
X2
-0.5
X2
-0.5
a.5
x1
0.5
x1
0.5
Lipschitz RNN
SVD Combo Net.
Antisym. RNN
0.0
0.0
0.5
0.0
X3
-055
-05
1.
1.0
1.0
0.5
-0.5
0.0
-0.5
0.0
0.5
Figure 3: Random sample of trajectories from the opinion dynamics test set. Predictions made by the model associated with the best MSE in Table 1. Crosses denote initial conditions, stars denote equilibrium points.
29

---

Published in Transactions on Machine Learning Research (04/2025)
Ground Truth
k-Lurie Net.
Lure Net.
Neural ODE
Lipschitz RNN
SVD Combo Net.
Antisym. RNN
200
2
Figure 4: Random sample of 30s trajectories from the noise-free, out of distribution opinion dynamics dataset. Predictions are made by the model associated with the best MSE Table 1. Crosses denote initial conditions, stars denote cquilibrium points.
Noise free Ground Truth
k-Lurie Net.
Lurie Net.
Neural ODE
2
2
Lipschitz RNN
SVD Comba Net.
Antisym. RNN
4yo 0
20247
0 2
0
Figure 5: Random sample 30s trajectorics from the noisy, out of distribution opinion dynamics dataset.
Predictions are made by the model associated with the best MSE Table Crosses denote initial conditions, stars denote equilibrium points.
30

---

Published in Transactions on Machine Learning Research (04/2025)
Ground Truth
k-Lurie Net.
Lurie Net.
Neural ODE
0.5
0.5
0.5
0.0
0.0
0.0
0.0
-1.5
-1.5_0.5
-1.0
0.5
-1.0_0.5
-1.0
0.0
0.0
0.0
0.0
0.5
0.5
0.5
0.5
Lipschitz
SVD Combo Net.
Antisym. RNN
0.5
0.0
-1
-1
0
X2
11
X2
X1
X1
Figure 6: Random sample of trajectories from the Hopfield network test set. Predictions are made by the model associated with the best MSE in Table Crosses denote initial conditions, stars denote equilibrium points.
Ground Truth
k-uurie Net.
Lurie Net.
Neural ODE
Lipschitz RNN
SVD Comba Net.
Antisym. RNN
5.0 25 05
00
255
2,5 50
Figure 7: Random sample of 30s trajectories from the noisc-free, out of distribution Hopfield network dataset. Predictions are made by the model associated with the best MSE Table Crosses denote initial conditions, stars denote equilibrium points.
31

---

Published in Transactions on Machine Learning Research (04/2025)
Noise Free Ground Truth
k Lurie Net.
Lurie le Net.
Neural ODE
Lipachitz RNN
SVD Combo Net.
Antisym. RNN
2. aD
25
0.0
-50
Figure 8: Random sample of 30s trajectories from the noisy, out of distribution Hopfield network dataset.
Predictions are made by the model associated with the best MSE in Table 1. Crosses denote initial conditions, stars denote equilibrium points.
Ground Truth
k-Lurie Net.
Lurie Net.
Neural ODE
22
2.5
-5
0.0
5
0
x1
-2.5
55
Lipschitz RNN
SVD Combo Net.
Antisym. RNN
X3
0
-2
2
2
-5
-2
0
0
2
Figure 9: Random sample of trajectories from the simple attractor test set. Predictions made by the model associated with the best MSE in Table 1. Crosses denote initial conditions.
32

---

Published in Transactions on Machine Learning Research (04/2025)
Ground Truth
k-Lurie Net.
Lure Net.
Neural ODE
25
00
25
10
10
-10
-10
10
-0
Lipschitz RNN
SVD Combo Net.
Antisym. RNN
0.0
10
3.5
-5,0
-15
5
Figure 10: Random sample of 30s trajectories from the noise-free, out of distribution simple attractor dataset. Predictions are made by the model associated with the best MSE Table 1. Crosses denote initial conditions.
Noise free Ground Truth
Lurie Net.
Lurie Net.
Neural ODE
50
25
25
00
00
-30
-5.0
10
10
-10
-10
-20-10
-10
16
Lipschitz RNN
SVD Combo Net.
Antisym. RNN
5.0
50
00
00
-25
-2
-10
Figure 11: Random sample of 30s trajectorics from the noisy, out of distribution simple attractor datasct. Predictions are made by the model associated with the best MSE Table 1. Crosses denote initial conditions.
33


Reviews:
[{"review_id": "7N7yEf3nty", "reviewer": "Reviewer_cQ12", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n- Introduction the Lurie Network as a unifying architecture for modeling nonlinear dynamical systems, covering existing models such as RNNs and Neural Oscillators.\n- Introduction of k-contraction as a general stability measure, enabling global convergence to a point, line, or plane in the neural state space.\n- Rigorous sufficient condition for Lurie Network and Graph Lurier Network (GLN) to be k-contracting.\n- Introduction of unconstrained parametrizations that allow the k-contracting Lurie Network and GLN to be trained with standard gradient-based optimization algorithms.\n\n### Weaknesses\n- Lack of discussion on chaotic systems: the paper does not explicitly discuss chaotic systems, which are crucial in many nonlinear dynamical applications, for example, the Lorenz attractor is 3-contracting, but there is no discussion on whether the proposed framework can model or analyze such systems.\n- Experimental validation limited to toy datasets: While neuroscientific applications and real-world systems are mentioned in Introduction, the experiments are conducted only on toy datasets.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces the Lurie Network as a unifying architecture for modelling time-invariant nonlinear dynamical systems.  Several typical architectures, such as recurrent neural networks (RNNs) and neural oscillators, can be described using the Lurie Networks as special cases.  \nThis paper rigorously derives sufficient conditions under which a Lurier Network has k-contraction properties.  This result provides a theoretical guarantee of global convergence to a point, line, or plane in the state-space, and extends the applicability of stability analysis to broader classes of dynamical systems.\nThis paper also derives a special parametrization of the Lurie  Network, ensuring the Lurie Network possesses the k-contraction property. Also this parametrization is unconstrained and enables simple gradient-based training of the Lurier Network with k-contraction property.\n\nThis paper also propose the Graph Lurier Network (GLN) to capture interactions between interconnected subsystems possesing the k-contraction property.  This paper provides the corresponding theoretical results for GLN.\n\nThis paper presents empirical results that imply the effectiveness of the Lurie Network with the special parametrization and k-contraction property.  This paper provides the numerical evidence of the high-performance in learning problem of dynamical systems in comparison with other existing architectures, for example, Neural ODE.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Since chaotic systems are inherently unpredictable over long horizons,  the long-term prediction conducted in this paper will not work, however, qualitative discussions (e.g., how k-contraction might capture local stability properties or the limits of predictability) would strengthen the contribution.\n\n- Although real-data experiments might be beyond the scope of this paper, a discussion on how the Lurie Network and GLN could be applied to real data in neuroscience would improve the paper.\n\n### Minor error\n- In Section 2.1, the authors use $\\mathcal{SO}(n)$ and $so(n)$ to describe the orthogonal and skew-symmetric matrices, but they are usually denote the special orthogonal group (orthogonal matrices with determinant 1) and its Lie algebra.  The author should follow typical notation of the theory of Lie group.\n- p.5, l.5, there is double \"it\".\n- p.18, in the formulas in center, use roman $\\exp$ not italic $exp$.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "93ghO7cwuH", "reviewer": "Reviewer_AryQ", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is very well-written, with clear motivation. The mathematical proofs were easy to follow, and the final theorems are intuitive. I do not see any weaknesses, given the target audience of this paper.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper presents an thorough analysis of the k-contraction properties of Lurie Networks, an important class of nonlinear models which appear across many different domains. The paper shows that Lurie Networks which are parameterized to be k-contracting have favorable (and sometimes superior) performance properties compared to 1-contracting networks trained on the same task.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "None.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Ej605c08yc", "reviewer": "Reviewer_8uBL", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. Authors introduce and explain techniques from the control theory that are not well-known in the ML community. I believe these techniques can be beneficial for other researchers.\n2. The motivation is well-articulated and the organisation of the paper is clear enough to understand the essential ideas.\n3. Obtained theoretical results are well-suited for practical realisation.\n4. Provided code and comments in the appendix allows for reproduction of results.\n5. Appendix contains many technical details and definitions that can be useful for readers unfamiliar with k-contraction analysis.\n\nWeaknesses (explained in details in the section below):\n1. k-contraction property seems too restrictive especially for Hopfield networks and for some cases it is not selective enough\n2. Numerical part of the article is not convincing\n3. There are some problems with organisation of the text", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Authors propose a Neural ODEs called Lurie Network. A parametric form of Lurie Networks encompasses several classes of continuous RNNs, Hopfield networks, neural oscillators, etc. To ensure stability authors use k-contraction analysis and derive several convenient parametrisations that guarantee k-contraction property for Lurie Network. It is demonstrated on the set of synthetic benchmarks that constrained Lurie Networks are competitive with other architectures.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. **Clarity and organisation**\n   1. Authors gather theoretical facts and definitions in the appendix. Unfortunately, the main text of the article does not reference appendix too often. Below are several examples where a reference to the appendix would be appropriate:\n      1. [Section 2.1] $k-$multiplicative and $k$-additive compound of matrix are not widely known, consider adding references to Definition 1 and Definition 2 in Appendix A\n      2. [Theorem 1] It is hard to understand the content of Theorem 1 without explicit definition of what $k$-contraction means.\n   2. Examples of Lurie networks in Section 3.1. are all RNNs. Can the authors please provide some other examples, e.g., Hopfield networks and Linear State Space Models (mentioned in the introduction), Neural Oscillators (mentioned in the abstract), etc?\n   \n2. **k-contraction is too restrictive and not selective enough**\n   1. It seems to me that $k$-contraction is convenient but too restrictive. I think this should be discussed in the article. I will provide several examples:\n      1. $2$-contraction illustrated in Figure 2 implies that all steady states of Hopfield network with this property are confined to the line. It seems that this condition is too restrictive for associative memory to be practically useful. Given that it is unfair to claim in the Section 6 that \"A 2-contracting model can replicate the behaviour of associative memory, where every stored pattern corresponds to an equilibrium\". \n      2. If one applies Theorem 2 to the Lurie network with $A=-I$, the condition for $1$-contraction is $\\\\left\\\\|B\\right\\\\|_{F} < \\frac{1}{L}$ where $L$ is Lipschitz constant of nonlinearity used. It seems that under a much milder sufficient condition $\\\\left\\\\|B\\\\right\\\\|_2<\\frac{1}{L}$ the resulting ODE still has a single global steady state.\n   2. I also think that $k$-contraction is not selective enough. I am not entirely sure, but from Theorem 2 it looks like $k+1$-contraction implies $k$-contraction. If this is the case, $3$-contraction from Figure 1 can also be used to learn Hopfield network (with memories restricted to the line) and contraction with single globally stable steady state.\n   \n3. **Numerical part**\n   1. I think that training time provided in Table 1 and Table 2 is misleading. I performed my own experiments with Hopfield data provided by authors. The results can be found in https://colab.research.google.com/drive/1arYUr0RkWGR-NH8mT_o-aUHtzdQPmUfO?usp=sharing. What I observed for NeuralODE:\n      1. For batch size $100$, $100$ epoch will take about $12$ minutes on GPU and $8$ minutes on CPU using jax, diffrax, optax and equinox.\n      2. Training time will depend on the batch size, solver, number of checkpoints, adjoint method, optimiser, precise details of architecture, hardware used, etc. It is virtually impossible to claim that one method is superior to the other one: it will require very expensive gridsearch on the space of hyperparameters.\n\n\t Beside these details, there are also different training strategies that are more suitable for NeuralODEs. For example, one can train with flow matching avoiding direct solution of ODE on the training phase. In my view, all that makes it hard to perform a fair comparison of different methods, so I suggest extending discussion on this and removing training time from the tables.\n\t \n   2. What is the reason Neural ODE has a much larger number of parameters? Why are there no experiments on the influence of network size? This information is of interest also for the out of distribution generalisation discussed in Table 2.\n   \n   4. There is an implicit assumption that we know which Lurie network to use for a particular dataset. What happens if we use a $1$-contracting or $3$-contracting network for Hopfield dataset? What happens if we use a $2$-contracting network for an attractor dataset? What happens if Hopfield dataset has three attractors not confined to a single line but we still use the $2$-contracting Lurie network?\n   \n   5. One claim of the authors is that constraints of $k$-contractive Lurie networks are beneficial for accuracy, robustness, and numerical efficiency. However, from the empirical data it is not clear whether the constraints help in learning and that they are beneficial for out-of-distribution generalisation. If one compares Table 1 and Table 2, contractive Lurie network is two orders less accurate on Option and Attractor datasets and $5\\times$ less accurate for Hopfield dataset. Unconstrained Lurie network lost one order of accuracy for both Option and Attractor and performed better on Hopfield dataset. Given that some may argue that unconstrained Lurie network is more robust. It is also not clear why the unconstrained Lurie network has fewer parameters than the $k$-contractive one. The same goes for all other benchmarks: unconstrained Lurie network is smaller and clearly undertrained.\n   \n   7. Another interesting part in comparison of unconstrained Lurie and constrained Lurie networks is that they are implemented differently. In the code authors provided they use all kinds of matrices with special structures to define parameters of $k$-contractive Lurie network and just linear layers to implement unconstrained Lurie networks. The difference in the parametrisation leads to distinct training dynamics. I suggest authors use the same parametrisation (and the same number of parameters) for both constrained and unconstrained Lurie networks to remove the difference between training dynamics. For example since in constrained Lurie network authors represent matrix $B$ as $U_B\\Sigma_B V_B^\\top$, the same expression should be used to parametrise $B$ in the unconstrained Lurie networks since SVD is unique (up to a complex sign) and exist for arbitrary matrix.", "broader_impact_concerns": "na", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: I recommend acceptance without revision. The paper introduces an elegant and general framework for learning stable neural dynamics via Lurie Networks, characterized by k-contracting properties. The authors provide strong theoretical underpinnings, including new sufficient conditions and an unconstrained parametrization that facilitates efficient training. Despite concerns about limitations in empirical breadth and application to real-world data, the response to reviewers was detailed and constructive, and the revised manuscript adequately addresses key points. The clear exposition and reproducibility further enhance its suitability for publication. While the practical significance may benefit from future real-world applications, the methodological contributions are substantial and timely.
Decision: Accept (poster)