Gradient Flossing: Improving Gradient Descent through Dynamic Control of Jacobians
Authors: Rainer Engelken
Venue: NeurIPS 2023 poster (2023)
Keywords: exploding/vanishing gradients
Lyapunov exponents
Lyapunov spectrum
chaos
RNN
condition number
Jacobian
TL;DR: We propose a novel method to analyze and ameliorate the vanishing/exploding gradient problem based on concepts form dynamical systems theory (Lyapunov exponents)

Abstract: Training recurrent neural networks (RNNs) remains a challenge due to the instability of gradients across long time horizons, which can lead to exploding and vanishing gradients. Recent research has linked these problems to the values of Lyapunov exponents for the forward-dynamics, which describe the growth or shrinkage of infinitesimal perturbations. Here, we propose gradient flossing, a novel approach to tackling gradient instability by pushing Lyapunov exponents of the forward dynamics toward zero during learning. We achieve this by regularizing Lyapunov exponents through backpropagation using differentiable linear algebra. This enables us to "floss" the gradients, stabilizing them and thus improving network training. We show that gradient flossing controls not only the gradient norm but also the condition number of the long-term Jacobian, facilitating multidimensional error feedback propagation. We find that applying gradient flossing before training enhances both the success rate and convergence speed for tasks involving long time horizons.
For challenging tasks, we show that gradient flossing during training can further increase the time horizon that can be bridged by backpropagation through time. Moreover, we demonstrate the effectiveness of our approach on various RNN architectures and tasks of variable temporal complexity. Additionally, we provide a simple implementation of our gradient flossing algorithm that can be used in practice. Our results indicate that gradient flossing via regularizing Lyapunov exponents can significantly enhance the effectiveness of RNN training and mitigate the exploding and vanishing gradients problem.

Gradient Flossing: Improving Gradient Descent through Dynamic Control of Jacobians
Rainer Engelken
Zuckerman Mind Brain Behavior Institute
Columbia University
New York, USA
re2365@columbia.co edu
Abstract
Training recurrent neural networks (RNNs) remains challenge due to the insta- bility of gradients across long time horizons, which can lead to exploding and vanishing gradients. Recent research has linked these problems the values of Lyapunov exponents for the forward-dynamies, which describe the growth or shrinkage of infinitesimal perturbations. Here, we propose gradient flossing, a novel approach to tackling gradient instability by pushing Lyapunov exponents of the forward dynamics toward zero during learning. We achieve this by regu- larizing Lyapunov exponents through backpropagation using differentiable linear algebra. This enables us to "floss" the gradients, stabilizing them and thus improv- ing network training. We demonstrate that gradient flossing controls not only the gradient norm but also the condition number the long-term Jacobian, facilitating multidimensional error feedback propagation. We find that applying gradient flossing prior to training enhances both the success rate and convergence speed for tasks involving long time horizons. For challenging tasks, we show that gradient flossing during training can further increase the time horizon that can be bridged by backpropagation through time. Moreover, we demonstrate the effectiveness of our approach on various RNN architectures and tasks of variable temporal com- plexity. Additionally, we provide simple implementation of our gradient flossing algorithm that can be used in practice. Our results indicate that gradient flossing via regularizing Lyapunov exponents can significantly enhance the effectiveness of RNN training and mitigate the exploding and vanishing gradients problem.
1 Introduction
Recurrent neural networks are commonly used both in machine learning and computational neu- roscience for tasks that involve input-to-output mappings over sequences and dynamic trajectories. Training is often achieved through gradient descent by the backpropagation of error information across time steps I 4. This amounts unrolling the network dynamics in time and recursively applying the chain rule to calculate the gradient of the loss with respect to the network parameters. Mathematically, evaluating the product of Jacobians of the recurrent state update describes how error signals travel across time steps. When trained on tasks that have Long-range temporal dependencies, recurrent neural networks are prone to exploding and vanishing gradients the exponential amplification or attenuation of recursive derivatives of recurrent network states over many time steps. Intuitively, evaluate how an output error depends on small parameter change at much earlier point in time, the error information has to be propagated through the recurrent network states iteratively. Mathematically, this corresponds to a product Jacobians that describe how changes in recurrent network state depend changes in the previous network state. Together, this product forms the long-term Jacobian. The singular value spectrum of the long-term Jacobian reg-
g These arise from
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

ulates how well error signals can propagate backwards along multiple time steps, allowing temporal credit assignment. A close mathematical correspondence of these singular values and the Lyapunov exponents of the forward dynamics was established recently [9] I0
12]. Lyapunov exponents
characterize the asymptotic average rate of exponential divergence or convergence of nearby initial
conditions and are cornerstone of dynamical systems theory 13 14. We will use this link to improve the trainability of RNNS.
Previous approaches that tackled the problem of exploding or vanishing gradients have suggested solutions at different levels. First, specialized units such as LSTM and GRU were introduced, which have additional latent variables that can be decoupled from the recurrent network states via multiplicative (gating) interactions. The gating interactions shield the latent memory state, which can therefore transport information across multiple time steps [3] 6
[5]. Second, exploding gradients can be avoided by gradient clipping. which re-scales the gradient norm [6] or their individual elements [17] if they become too large [8]. Third, normalization schemes like batch normalization, layer norm and group norm prevent saturated nonlinearities that contribute to vanishing gradients [19 20 21
Fourth, it was suggested that the problem of exploding/vanishiin gradients can be ameliorated by specialized network architectures, for example, antisymmetric networks [22], orthogonal/unitary
initializations (23 24 25 coupled oscillatory RNNs (26), Lipschitz RNNs (27,, structured state space models 288 29
30 31,, echo state networks 32 33 (recurrent) highway networks [3] 35], and
stable limit cycle neural networks [T 36, 37.. Fifth, for large networks, a suitable choice of weights can guarantee a well-conditioned Jacobian at initialization [23] 38, 39, 40 42, 45
44 45] 46 47).
These initializations are based on mean- field methods, which become exact only in the netwokk limit. Such initialization schemes have also been suggested for gated networks [45). However, even when initializing the network with well-behaved gradients, gradients will typically not retain their stability during training once the network parameters have changed.
Here, we propose novel approach to tackling this challenge by introducing gradient flossing, technique that keeps gradients well-behaved throughout training. Gradient flossing is based on recently described link between the gradients of backpropagation through time and Lyapunov exponents, which are the time-averaged logarithms of the singular values of the long-term Jacobian 19]11 37). Gradient flossing regularizes or several Lyapunov exponents to keep them close to zero during training. This improves not only the error gradient norm but also the condition number of the long-term Jacobian. As result, error signals can be propagated back over longer time horizons. We first demonstrate that the Lyapunov exponents can be controlled during training by including an additional loss term. We then demonstrate that gradient flossing improves the gradient norm and effective dimension of the gradient signal. We find empirically that gradient flossing improves test accuracy and convergence speed on synthetic tasks over a range of temporal complexities. Finally, we find that gradient flossing during training further helps to bridge long-time horizons and show that it combines well with other approaches to ameliorate exploding and vamishing gradients, such as dynamic mean-field theory for initialization, orthogonal initialization and gated units. Our contributions include:
Gradient flossing, novel approach to the problem of exploding and vanishing gradients in recurrent neural networks based on regularization of Lyapunov exponents.
Analytical estimates of the condition number of the long-term Jacobian based on Lyapunov exponents.
Empirical evidence that gradient flossing improves training tasks that involve bridging long time horizons.
2 RNN Gradients and Lyapunov Exponents
We begin by revisiting the established mathematical relationship between the gradients of the loss function, computed via backpropagation through time, and Lyapunov exponents 9 12] and how it relates to the problem of vanishing and exploding gradients. In backpropagation through time, network parameters 0 are iteratively updated by stochastic gradient descent such that a loss L, is locally reduced [ 2 3 4. For RNN dynamics h,+1 = fy (h, x.+1), with recurrent network state h, external input x, and parameters 0, the gradient of the loss L, with respect to 0 is evaluated by
'Code in Julia using Flux 77 78 available http://ghhhbccooo/iiiieeeeeeells

---

unrolling the network dynamics in time. The resulting expression for the gradient is given by:
ac.
ac, r=t-1
Σ r=
0h,'+1 Th
Th,
ac.
Σ Ti(h,)
Ih, 00
(1)
00
0h,
00
Oh;
where T(( is composed of product of one-step Jacobians D, = Sh,
11,
f-1
ohr't1 Th
[-1
T, = II
D, D.
T'=T
Due the chain of matrix multiplications in Tr, the gradients tend vanish explode exponentially with time. This complicates training particularly when the task loss at time dependents on inputs x states h from many time steps prior which creates long temporal dependencies [5] well error signals can propagate back in time is constrained by the tangent space dynamics along trajectory h2, which dictate how local perturbations around each point on the trajectory stretch, rotate, shear, or compress as the system evolves.
6 [ How
The singular values of the Jacobian's 's product Tr, which determine how quickly gradients vanish or explode during backpropagation through time, are directly related to the Lyapunov exponents of the forward dynamics (| [2] Lyapunov exponents Å₁ 2 A2.... 2 Ax defined as the asymptotic time-averaged logarithms the singular values of the long-term Jacohian 13 48
A lim log(or,)
where Titt denotes the ith singular value of Thb,,) with 01.1 2 2 ... N. (See Appendix [ for details). This means that positive Lyapunov exponents in the forward dynamics correspond to exponentially exploding gradient modes, while negative Lyapunov exponents in the forward dynamics correspond to exponentially vanishing gradient modes.
In summary, the Lyapunov exponents give the average asymptotic exponential growth rates of infinitesimal perturbations in the tangent space of the forward dynamics, which also constrain the signal propagation in backpropagation for long time horizons. Lyapunov exponents close zero in the forward dynamics correspond to tangent space directions along which error signals are neither drastically attenuated amplified in backpropagation through time. Such close-to-neutral modes in the tangent dynamics can propagate information reliably across many time steps.
3 Gradient Flossing: Idea and Algorithm
We now leverage the mathematical connection established between Lyapunov exponents and the prevalent issue of exploding and vanishing gradients for regularizing the singular values of the long-term Jacobian. We term this procedure gradient flossing. To prevent exploding and vanishing gradients, we constrain Lyapunov exponents to be close to zero. This ensures that the corresponding directions in tangent space grow and shrink on average only slowly. This leads to a better-conditioned long-term Jacobian Th,)) We achieve this by using the sum of the squares of the first k largest Lyapunov exponent '₁ d2 ... ÅE as a loss function:
Bossing
X?
(4)
and evaluate the gradient obtained from backpropagation through time:
i Lasssing
axa
(5)
00
ae
This might seem like an ill-fated enterprise, as the gradient expression in Eq suffers from its own problem of exploding and vanishing gradients. However, instead of calculating the Lyapunov exponents by directly evaluating the long-term Jacobian T, (Fq(2), we use an established iterative reorthonormalization method involving QR decomposition that avoids directly evaluating the ill-
conditioned long-term Jacobian 22 50]
3

---

First, we evolve an initially orthonormal system Q. = = ... in the tangent space along the trajectory using the Jacobian D, This means to calculate
0h,+1
Sh,
Q++ D.Q.
(6)
at every time-step. Second, we extract the exponential growth rates using the QR decomposition,
Qs+1 Qs+1 R*
which decomposes Qs+1 uniquely into the product of an orthonormal matrix Qs+1 of size N x so Q Qs+1 1 kx and an upper triangular matrix R of size k l with positive diagonal elements. Note that the QR decomposition does not have to be applied at every step. just sufficiently often, ie., once every tons such that Q does not become ill ill-conditioned.
The Lyapunov exponents given by time-averaged logarithms of the diagonal entries of R
49
50
A lim log
I
R
lim log R
This way, the Lyapunov exponent can be expressed in terms of a temporal average over the diagonal elements of the R'-matris of a QR decomposition of the iterated Jacobian. To propagate the gradient of the square of the Lyapunov exponents backward through time in gradient flossing, we used an analytical expression for the pullback of the QR decomposition [5T]:
The backward pass of the QR
decomposition is given by [51 52 53
54
Q
Q copyltu(M)) R R-T
(8)
where M R.R' Q Q and the copyltu function generates symmetric matrix by copying the lower triangle of the input matrix to its upper triangle, with the element
Mmax(i.)).minf( [T] 52, 53 54]. We denote here adjoint variable as T = åc/ØT. A simple
implementation of this algorithm in pseudocode is:
Algorithm 1 Algorithm for gradient flossing of ls tangent space directions
initialize h, Q for 1 E do
for f = + do h + (h,
D dh.
Q D Q
if = 0 (mod toNs) then
Q. R t qr Q)
+= log R
end if end for
A = 7
2ct1 t O.
86
end for
For clarity, we described gradient flossing in terms of stochastic gradient descent, but we actually implemented with the ADAM optimizer using standard hyperparameters 77, 31 and B2. An example implementation in Julia (77 using Flux (78) is available here). Note that this algorithm also works for different recurrent network architectures. In this case, the Jacobians has size
n, where the number of dynamic variables of the recurrent network model. For example, in case of a single recurrent network of LSTM units, the Jacobian has size 2N 2N 12 46) The Jacobian matrix D can either be calculated analytically it can be obtained via automatic differentiation.
4 Gradient Flossing: Control of Lyapunov Exponents
In Fig we demonstrate that gradient flossing can set one or several Lyapunov exponents target value via gradient descent with the ADAM optimizer in random Vanilla RNNs initialized with different weight variances. The N units of the recurrent neural network follow the dynamics
h f(h,, x++1) Wo(h,) Vx
(9)

---

Th th
Th
II
number of flossed
Thy
32
16
10-3-
1
-Wi
By:
10-5
dbc
ahu-1
dbu
dbn
abu-1
h
hu
Pu-.)
W
20
40
60
100
Epochs
0.0
0
-0.5
2
2,0m
A.1m
-1.0
-4
number of flossed A)
k= 32
-1.5
target
6
16
actual
-2.0
60
Epochs
UN
Figure 1: Gradient flossing controls Lyapunov exponents and gradient signal propagation A) Exploding and vanishing gradients in backpropagation through time arise from amplifica- tion/attenuation of product of Jacobians that form the long-term Jacobian T, -I==~ 0h,'+1 ah, function of training epochs. Minimizing the
B) First Lyapunov exponent of Vanilla RNN mean squared error between estimated first Lyapunov exponent and target Lyapunov exponent A - 0.5, O by gradient descent. 10 Vanilla RNNs were initialized with Gaussian recurrent weights Way ~ No0, g² //N where values of were drawn g ~ Uni 1). C) Gradient flossing minimizes the square of Lyapunov exponents over epochs. D) Full Lyapunov spectrum of Vanilla RNN after a different number of Lyapunov exponents are pushed to zero via gradient flossing. Note, the variability of the Lyapunov exponents that were flossed. Parameters: network size 32 with 10 network realizations. Error bars in C indicate the 25% and 75% percentiles and solid line shows median.
The initial entries of W are drawn independently from Gaussian distribution with variance g²/ //N, where is a gain parameter that controls the heterogeneity of weights. We here use the transfer function d(x) = tanh(x:) (See appendix B for gradient flossing with ReLU and LSTM units). X, is sequence of inputs and V is the input weight. X, is a stream of i.i.d. Gaussian input T, NO0 1) and the input weights V are (0, 1). Both W and V are trained during gradient flossing.
and
In BB show that for randomly initialized RNNs, the Lyapunov exponent can be modified by gradient flossing to match a desired target value. The networks were initialized with 10 different values of initial weight strength g chosen uniformly between 0 and 1. During gradient flossing, they quickly approached three different target values of the first Lyapunov exponents X {-1, -0.5, 0} within less than 100 training epochs with batch size B = 1. We note that gradient flossing with positive target X Jarget seems not to arrive a positive Lyapunov exponent λς-
Fig [ shows gradient flossing for different numbers of Lyapunov exponents k. Here, during gradient- descent, the sum the squares of 1, 16, or Lyapunov exponents used loss in gradient flossing (see Fig [ Fig ID shows the Lyapunov spectrum after flossing, which now has 1, 16, or 32 Lyapunov exponents close zero. We conclude that gradient flossing can selectively manipulate one, several, or all Lyapunov exponents before or during network training. Gradient flossing also works for RNNs ReLU and LSTM units (See appendix [B] Further, we find that the computational bottleneck of gradient flossing is the QR decomposition, which has a computational complexity of O (N k-²). both the forward pass and in the backward pass. Thus, gradient flossing of the entire Lyapunov spectrum is computationally expensive. However, as we will show, not all Lyapunov exponents need be flossed and only short episodes of gradient flossing are sufficient for significantly improving the training performance.
5

---

5
Gradient Flossing: Condition Number of the Kong-Terma Jacobian
condition number k2(Te())
B
1034
10
initial after flossing theory
10
1025
(T(T(())
10
simulátions
K2 (theory)
10.0
10t4
10a2
Initial after flossing theory simulations
1010
10¹
10
15
10a
100
200
300
10c
1010
10
1030
time horizon
tangent space dimension m
Xy (numerical)
Figure 2: Gradient flossing reduces condition number of the long-term Jacobian A) Condition number N2 long-term Jacobian (hz) as a function of time horizon - T at initialization (blue) and after gradient flossing (orange). Direct numerical simulations are done with arbitrary precision floating point arithmetic (transparent lines) with 256 bits per float, asymptotic theory based on Lyapunov exponents (dashed lines) (O) B) Condition number for different number of tangent space dimensions m. Simulations (dots) and Lyapunov exponent based theory (dashed lines) at initialization (blue) and after gradient flossing (orange). Gradient flossing increases the number of tangent space dimensions available for backpropagation for a given condition number (Grey dotted line guide for eye for K2 10 First 15 Lyapunov exponents condition number obtained via direct numerical simulations vs. Lyapunov exponent-based. Colors denote the number of flossed Lyapunov exponents k. Parameters: g 1, batch size N 80, epochs 500, T = 500, gradient flossing for Ey = 500 epochs. Input x, identical to delayed XOR task in Fi.
flossed. C) Comparison of
Gradient flossing improves the condition number of the long-term Jacobian which constrains the error signal propagation across long time horizons in backpropagation (Fig 2). The condition number R2 of linear map A measures how close the map is being singular and is given by the ratio of the largest
A well-conditioned Jacobian is essential for efficient and fast learning 23
55
56
(A) According to the singular value max and the smallest singular values min. so K2 R2(A) (A) = aam(AA
rule of thumb given in [57], if R2(AA) 10P, can anticipate losing least p digits of precision when solving the equation Ax = b. Note that the long-term Jacobian T; is composed of a product of Jacobians, which generically makes ill-conditioned. To nevertheless quantify the condition number numerically, we use arbitrary-precision arithmetic with 256 bits per float. We find numerically that the condition number of T, exponentially diverges with the number of time steps (Fig ( We compare the numerically measured condition number X2 with an asymptotic approximation of the condition number based on Lyapunov exponents that are calculated in the forward pass and find good match (Fig AA).
Our theoretical estimate of the condition number K2 of an orthonormal system Q of size N x that is temporally evolved by the long-term Jacobian T, is:
01 (T,())) so exp ((A1 (m)(1-T))
K2 (.+r) = =2 (T(((h))).)))) O. T (T,(T)))))
(10)
where (h,)) and om (T, (hr)) are the first and math singular value of the long-term Jacobian. We note that this theoretical estimate of the condition number follows from the asymptotic definition of Lyapunov exponents and should be exact in the limit of long times. We find that gradient flossing roduces the condition number by a factor whose magnitude increases exponentially with time (orange expect that gradient flossing has stronger effect problems with a long
in Fig(A) Thus, time horizon to bridge. We will later confirm this numerically.
Moreover, Lyapunov exponents enable the estimation of the number of gradient dimensions available for the backpropagation of error signals. Generally, the long- term Jacobian ill-conditioned, however, the Lyapunov spectrum provides for given number of tangent space dimensions an estimate of the condition number. This indicates how close to singular the gradient signal for a given number of tangent space dimensions is. Given fixed acceptable condition by noise level or doating-point precision-ww we observe that gradient flossing increases the number of usable tangent space dimensions for backpropagation (Fig 23).
for example,

---

Finally, we show that the asymptotic estimate of the condition number based on Lyapunov exponents can even predict differences in condition number that originate from finite network size N (Fig[ZC). We emphasize that this goes beyond mean-field methods, which become exact only in the large- network limit N- + oo and usually do not capture finite-size effects [58]
(see appendix G
Initial Gradient Flossing Improves Trainability
delayed copy task
delayed XOR task
10-4
10-2
test test loss
gradient flossing gr gradient flassing
gradient flossing
test
10-3
gradient flossing
10 3
10
2000
4000
6000
8000
10000
0
2000
4000
6000
8000
10000
Epochs
Epochs
delayed copy task
delayed XOR task
10
10 2
10-3
test loss
10
test
10-1
10
20
40
60
80
10
20
30
40
50
60
70
80
task difficulty (delay d)
task difficulty (delay d)
Figure 3: Gradient Gladientflossing improves trainability on tasks that involve long time horizons A) Test error for Vanilla RNNs trained on delayed copy task 3/2 = Zd-d for d = 40 with and without gradiens flossing flossing. Solid lines are medians across 5 network realizations. B) Same A for delayed XOR task with y/- = xt-d//2 dt=d// C) Mean final test loss as a function of task difficulty (delay d) for delayed copy task. D) Mean final test loss as function of task difficulty (delay d) for delayed XOR task. Parameters: g = 1, batch size = 16, N = 80, epochs = 10a T 300, gradient flossing for E 500 epochs on k 75 before training. Shaded regions in C and D indicate the 20% and 80% percentiles and solid line shows mean. Dots are individual runs Task loss: MSE(y, 9).
We next present numerical results on two tasks with variable spatial and temporal complexity. demonstrating that gradient flossing before training improves the trainability of Vanilla RNNs. We call gradient flossing before training in the following preflossing. For preflossing, we first initialize the network randomly, then minimize Lasssing = Σ X using the ADAM optimizer and subsequently train on the tasks. We deliberately do not use sequential MNIST or similar toy tasks commonly used to probe exploding/vanishing gradients, because we want a task where the structure of long-range dependencies in the data is transparent and can be varied as desired.
First, we consider the delayed copy task. where a scalar stream of random input numbers must be reproduced by the output y delayed by d time steps, i.c. = It-d. Although the task itself is trivial and can be solved even by a linear network through delay line (see appendix E RNNs encounter vanishing gradients for large delays d during training even with 'critical' initialization with 1. Our experiments show that gradient flossing can substantially improve the performance of RNNs on this task 3A, C). While Vanilla RNNs without gradient flossing fail to train reliably beyond 20, Vanilla RNNs with gradient flossing can be reliably trained for = 40 (Fig 3C). Note that we flossed here A = 40 Lyapunov exponents before training. We will later investigate the role of the number of flossed I Lyapunov exponents.
Second, we consider the temporal XOR task, which requires the RNN to perform a nonlinear input- output computation on sequential stream of scalar inputs, i.e., Vt- [z2-4/2 æt-a),, where denotes time delay of d time steps (For details see appendix [H). 3 demonstrates that gradient flossing helps to train networks on a substantially longer delay d. We found similar improvements through gradient flossing for RNNs initialized with orthogonal weights (see appendix G
7

---

7
Gradient Flossing During Training
delayed temporal XOR
delayed spatial XOR
A
100
100
flossing during training
90
preflossing
test accuracy (%)
90
test accuracy
no flossing
80
60
60
50
50
2000
4000
6000
8000
10000
0
2000
4000
6000
8000
10000
Epochs
Epochs
100
100
test accuracy (%)
test accuracy
80
flossing during training preflossing
70
70
no flossing
60
50
10
30
10
60
complexity (delay d)
complexity (delay
Figure 4: Gradient flossing during training further improves trainability
A) Test accuracy for Vanilla RNNs trained on delayed temporal binary XOR task yo = 2t-d/2 d with gradient flossing during training (green), preflossing (gradient flossing before training) (orange), and with no gradient flossing (blue) for d 70. Solid lines are mean across 20 network realizations, individual network realizations shown in transparent fine lines. B) Same as A for delayed spatial XOR task with ÿ =/__ x2- d 2 xr x²- Parameters (g 1, batch size b- 16). C) Test accuracy as a function of task difficulty (delay d) for delayed temporal XOR task. D) Test accuracy as a function of task difficulty (delay d) for delayed spatial XOR task. Parameters: g 1, batch size b = 16, N = 80, epochs = 10a, T = 300, gradient flossing for E 500 epochs on k = 75 before training and during training for green lines, and only before training for orange lines. Same plotting conventions as previous figure. Task loss: cross- cross-entropy between and ÿ.
We next investigate the effects of gradient flossing during the training and find that gradient flossing during training can further improve trainability. We trained RNNs on two more challenging tasks with variable temporal complexity and performed gradient flossing either both during and before training, only before training, not at all.
Fig AA shows the test accuracy for Vanilla RNNs training on the delayed temporal XOR task Vt Zt-d/2 @ Iz-d with random Berouulii process E 1). The accuracy of Vanilla RNNs falls to chance level for d 2 40 (Fig 4C). With gradient flossing before training, the trainability can be improved, but still goes to chance level for d 70. In contrast, for networks with gradient flossing during training, the accuracy is improved to > 80% at d = 70. In this case, we preflossed for 500 epochs before task training and again after 500 epochs of training D the networks have to perform the nonlinear XOR operation direediimessionll binary input signal x,, xx, and x² and generate the correct output with delay of d steps. While the solution of the task itself not difficult and could even be implemented by hand (see appendix). the task is challenging for backpropagation through time because nonlinear temporal associations bridging long time horizons have be formed. Again, we observe that gradient flossing before training improves the performance compared to baseline, but starts failing for long delays >> 60. In contrast, networks that are also flossed during training can solve even more difficult tasks (Fig (). We find that after gradient flossing, the norm of the error gradient with respect to initial conditions he amplified (appendix C). Interestingly, gradient flossing can also be detrimental to task performance if it is continued throughout all training epochs (appendix
the task. In Fig 4B. = =i_d @ xi_a @ x on a
We note that merely regularizing the spectral radius of the recurrent weight matrix W or the individual one-step Jacobians D. numerically or based on mean-field theory does not yield such training improvement. This suggests that taking the temporal correlations between Jacobians D, into account is important for improving trainability.
8

---

7.1 Gradient Flossing for Different Numbers of Flossed Lyapunov Exponents
We investigated how many Lyapunov exponents have to be flossed to achieve an improvement in training success (Fig 5.. We studied this in the binary temporal delayed XOR task with gradient flossing during training (same as Fig 3) and varied the task difficulty by changing the delay d.
We found that as the task becomes more difficult, networks where not enough Lyapunov exponents k are flossed begin to fall below 100% test accuracy (Fig Correspondingly, when measuring final test accuracy as a function of the number of flossed Lyapunov exponents, we observed that more Lyapunov exponent k have to be flossed to achieve 100% accuracy the tasks become more difficult BB) We also show the entire parameter plane of median test accuracy as a function of both number of flossed Lyapunov exponents k and task difficulty (delay d), and found the same trend (Fig Overall, we found that tasks with larger delay d require more Lyapunov exponents close to zero. We note that this might also partially be caused by the 'streaming" nature of the task: in our tasks, longer delays automatically imply that more values have to be stored as at any moment all the values in the 'delay line' have to be remembered to successfully solve the tasks. This is different from tasks where a single variable has to be stored and recalled after a long delay. It would be interesting to study tasks where the number of delay steps and the number of items in memory can be varied independently.
Finally, we did the same analysis on networks with only preflossing (gradient flossing before training) and found the same trend (supplement Fig (D), however, in that case even if all N Lyapunov exponents were flossed, thus N. they were not able to solve the most difficult tasks. This seems to indicate that gradient flossing during training cannot be replaced by just gradient flossing more Lyapunov exponents before training.
B
100
100
test accuracy (%)
test test accuracy (%)
k
delay
80
80
20
30
40
60
70
80
20
40
20
60
complexity (delay d)
number of flossed X k
Figure 5: Gradient flossing for different numbers of flossed Lyapunov exponents A) Test accuracy for delayed temporal XOR task as function of delay d with different numbers flossed Lyapunov exponents k. B) Same data as A but here test accuracy as function of number
of flossed Lyapunov exponents k. Parameters: g 1,, batch size b = 16, N = 80, epochs 10 for delayed temporal XOR, epochs- = 5000 for delayed spatial XOR, T = 300, gradient flossing for E, = 500 epochs before training and during training for A, B. Shaded areas are 25% and 75% percentile, solid lines are means, transparent dots are individual simulations, task loss: cross-entropy between y and ý.
8 Limitations
The mathematical connection between Lyapunov exponents and backpropagation through time exploited in gradient flossing is rigorously established only in the infinite-time limit. It would be interesting to extend our analysis to finite-time Lyapunov exponents.
Furthermore, the backpropagation through time gradient involves a sum over products of Jacobians of different time periods f т, but the Lyapunov exponent only considers the asymptotic longest product. Additionally, Lyapunov exponents characterize the asymptotic dynamics on the attractor of the dynamics, whereas RNNs often exploit transient dynamics from some initial conditions outside or towards the attractor.
Although our proposed method focuses on exploiting Lyapunov exponents, it neglects the geometry of covariant Lyapunov vectors (59], which could be used to improve training performance, speed, and reliability. Additionally, it is important to investigate how sensitive the method is to the choice of orthonormal basis employed because it is only guaranteed to become unique asymptotically [60].

---

Finally, the computational cost of our method scales with O(Nk²), where N is the network size and k is the number of Lyapunov exponents calculated. To reduce the computational cost, we suggest doing QR decomposition only sufficiently often to ensure that the orthonormal system is not ill-conditioned and using gradient flossing only intermittently or as pretraining. One could also calculate the Lyapunov spectrum for shorter time interval or use a cheaper proxy for the Lyapunov spectrum and investigate more efficient gradient flossing schedules.
9 Discussion
We tackle the problem of gradient signal propagation in recurrent neural networks through a dy- namical systems lens. We introduce novel method called gradient flossing that addresses the problem of gradient instability during training. Our approach enhances gradient signal stability both before and during training by regularizing Lyapunov exponents. By keeping the long-term Jacobian well-conditioned, gradient flossing optimizes both training accuracy and speed. To achieve this, we combine established dynamical systems methods for calculating Lyapunov exponents with an analytical pullback of the QR factorization. This allows us to establish and maintain gradient stability in in manner that is memory-efficient Our method is applicable to arbitrary RNN architectures, nonlinearities, and also neural ODEs [61]. Empirically, pre-training with gradient flossing enhances both training speed and accuracy. For difficult temporal credit assignment problems, gradient flossing throughout training further enhances signal propagation. We also demonstrate the versatility of our method on a set of synthetic tasks with controllable time-complexity and show that it can be combined with other approaches to tackle exploding and vanishing gradients, such as dynamic mean-field theory for initialization, orthogonal initialization and specialized single units, such as LSTMs.
numerically stable, and exact across long time horizons.
Prior research on exploding and vanishing gradients mainly focused on selecting network architectures that are less prone to exploding/vanishing gradients or finding parameter initializations that provide well-conditioned gradients least at the beginning of training. Our introduced gradient flossing can be seen as a complementary approach that can further enhance gradient stability throughout training. Compared to the work on picking good parameter initializations based on random matrix theory [46] and mean-fiele field heuristics (45), gradient flossing provides several improvements: First, mean-field theory only considers the gradient flow at initialization, while gradient flossing can maintain gradient flow and well-conditioned Jacobians throughout the training process. Second, random matrix theory and mean-field heuristics are usually confined the limit of large networks (58), while gradient flossing can be used for networks of any size. The link between Lyapunov exponents and the gradients of backpropagation through time has been described previously 9 12] and has been spelled out analytically and studied numerically [10]1] I1 62, 63 64]. In contrast, we use Lyapunov exponents here not only as a diagnostic tool for gradient stability but also to show that they can directly be part of the cure for exploding and vanishing gradients.
Future investigations could delve further into the roles of the second to Nthh Lyapunov exponents in trainability, and how it is related to the task at hand, the rank of the parameter update, the dimensional- ity of the solution space, as well as the network dynamics (sce also 37 65 66)). Our results suggest a trade-off between trainability across long time horizons and the nonlinear task demands that is worth exploring in more detail (appendix C. Applying gradient flossing to real-time recurrent learning and its biologically plausible variants is another avenue (67]. Extending gradient flossing to feedforward networks, space models and transformers is a promising avenue for future research (see also 70 71 72).. While Lyapunov exponents are only strictly defined for dynamical systems, such as maps or flows that are endomorphisms, the long-term Jacobian of deep feedforward networks could be treated similarly. This could also provide link between the stability of the network against adversarial examples and its dynamic stability, as measured by Lyapunov exponents. Given that time-varying input can suppress chaos in recurrent networks 9 12
31
47 68, 69
73
74 75 76), we
anticipate they may exacerbate vanishing gradients. Gradient flossing could also be applied in neural architecture search, to identify and optimize trainable networks. Finally, gradient flossing is applicable to other model parameters, as well. For instance, gradients of Lyapunov exponents with respect to single-unit parameters could optimize the activation function and single-neuron biophysics in biologically plausible neuron models.
10

---

Acknowledgments and Disclosure of Funding
I thank E. Izhikevich, F. Wolf, S. Goedeke, J. Lindner, L.F. Abbott, L. Logiaco, M. Schottdorf, G. Wayne and P. Sokol for fruitful discussions and J. Stone, L.F. Abbott, M. Ding, O. Marshall, S. Goedeke, S. Lippl, M.P. nouzel, J. Lindner and the reviewers for feedback on the manuscript. I thank Jinguo Liu for the Julia package BackwardsLinalge
and also the developers of Julia and Flux 78). Research supported by NSF NeuroNex Award (DBI-1707398) the Gatsby Charitable Foundation (GAT3708), the Simons Collaboration for the Global Brain (542939SPI), and the Swartz Foundation (2021-6).
References
[1] P. WERBOS. Beyond Regression: Ph. D. dissertation, Harvard University, 1974.
[2] DB Parker. Learning-logic (TR-47). Center for Computational Research in Economics and Management Science. MIT-Press, Cambridge, Mass, 1985.
[3] Y. LECUN. Une procedure d' d'apprentissage ponr rescau scuil asymetrique. Proceedings of Cognitiva 85, pages 599-604, 1985.
[4] David E. Rumelhart, Geoffrey F. Hinton, and Ronald J. Williams. Learning representations by back- propagating errors. Nature, (636688):533, October 1986.
[5] Sepp Hochreiter. Untersuchungen dynamischen neuronalen Netzen. PhD thesis, 1991.
[6] Sepp Hochreiter and Jürgen Schmidhuber. Long Short-Term Memory, Neural Computation, 978):1735- 1780, 1997.
[7] Y. Bengio, P. Simard, and P. Frasconi. Learning long-term dependencies with gradient descent difficult. IEEE Transactions on Neural Networks, 5(2):157-166, March 1994.
[8] Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. the difficulty of training Recurrent Neural Networks. urXiv:1211.506 [cx], November 2012. arXiv: 1211.5063.
Rainer Engelken, Fred Wolf, and I., F. Abbott. Lyapunow spectra of chaotic recurrent neural networks. arXiv:2006.0244 02427 [nlin, q-bio], June 2020. arXiv: 2006.02427.
[10] Jonas Mikhaeil, Zahra Monfared, and Daniel Durstewitz. On the difficulty of learning chaotic dynamics with RNNs. Advances Neural Information Processing Systems, 35:1 15:11129--1312, December 2022.
[11] Il Memming Park. Ábel Ságodi, and Piotr Aleksander Sokól. Persistent learning signals and working memory without continuous attractors. ArXiv, page arXiv:2308. 12585v1, August 2023.
[12] Rainer Engelken, Fred Wolf, and L., R. Abbott. Lyapumov spectra of chaotic recurrent neural networks. Physical Review Research, 5(4):043044, Octoher 2023.
[13] J. P. Eckmann and D. Ruelle. Ergodic theory of chaos and strange attractors. Reviews of Modern Physics, 57(3):617-656, July 1985.
[14] Arkady Pikovsky and Antonio Politi. Lyapunov Exponents: A Tool to Explore Complex Dynamics. Cambridge University Press, Cambridge, February 2016.
[15] Kyunghyun Cho, Bart van Merrienhoer, Dzmitry Bahdanau, and Yoshua Bengio. On the Propertics of Neural Machine Translation: Enender-Deendder Approaches. Technical report, September 2014. ADS Bibcode: 2014arXiv1400 12599 Type: article.
[16] Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. On the difficulty of training recurrent neural networks. In Proceedings of the 30th International Conference International Conference on Machine Learning Volume 28, ICML'1 13, pages 11I-1310-I1 111-1318, Atlanta, GA, USA. June 2013. JMLR.org.
[17] Tomáš Mikolov. Statistical language models based neural networks. Ph.D. thesis, Bmo University of Technology, Faculty of Information Technology, Brno, CZ, 2012.
[18] Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. MIT Press, November 2016. Google-Books-ID: omivDQAAQBAJ.
[19] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E. Hinton. Layer Normalization, July 2016. arXiv:1607006550 (cs, stat].

---

[20] Sergey loffe and Christian Szegedy. Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift. In Proceedings of the 32nd International Conference Machine Learning, pages 446. PMLR, June 2015.
[21] Yuxin Wu and Kaiming He. Group Normalization. pages 3 19, 2018.
[22] Bo Chang. Minmin Chen, Eldad Haber, and Bd H. Chi. AntisymmetricRNN: A Dynamical System View Recurrent Neural Networks. International Conference Learning Representations, December 2018.
[23] Andrew M. Saxe. James L. McClelland, and Surya Ganguli. Exact solutions to the nonlinear dynamics of learning in deep linear neural networks. arXiv:1312.6120 [cond-mat, q-bio, stat]. December 2013, arXiv: 1312.6120.
[24] Martin Arjovsky, Amar Shah, and Yoshua Bengio. Unitary Evolution Recurrent Neural Networks. In Proceedings of The 33rd International Conference on Machine Learning, pages 2016.
PMLR, June
[25] Kyle Helfrich, Devin Willmott, and Qiang Ye. Orthogonal Recurrent Neural Networks with Scaled Cayley Transform. In Proceedings of the 35th International Conference on Muchine Learning, pages 1969-1978. PMLR, July 2018.
[26] T. Konstantin Rusch and Siddhartha Mishra. Coupled Oscillatory Recurrent Neural Network (coRNN); An accurate and (gradient) stable architecture for learning long time dependencies. arXiv e-prints, page ar arXiv:2010.00951, October 2020.
[27] N. Benjamin Frichson, Omri Azencot, Alejandro Queiruga, Liam Hodgkinson, and Michael W. Mahoney. Lipschitz Recurrent Neural Networks. International Conference on Learning Representations, January 2021.
[28] Albert Gu, Karan Goel, and Christopher Ré. Efficiently Modeling Long Sequences with Structured State Spaces, August 2022. arXiv:2111.005: [es].
[29] Jimmy T. H. Smith, Andrew Warrington, and Scott W. Linderman. Simplified State Space Layers for Sequence Modeling, March 2023. 2eXiv:2208.04933 [es].
[30] Antonio Orvieto, Samuel L Smith, Albert Gu. Anushan Fernando, Caglar Gulcehre, Razvan Pascanu, and Soham De. Resurrecting Recurrntt Neural Networks for Long Sequences. Technical report, March 2023. ADS Bihcode: 0023arX1v2300063490
Type: article.
[31] Albert Gu and Tri Dao. Mamba: Lincar-Time Sequence Modeling with Selective State Spaces, December 2023, arXiv:2312.00752 [es].
[32] Herbert Jacger. The "ccho state" approach to analysing and training recurrent neural networks-with an erratum note. Bonn, Germany: German National Research Center for Information Technology GMD Technical Report, 148:34, 2001.
[33] Mustafa C. Ozturk, Dongming Xu, and José C. Príncipe. Analysis and Design of Echo State Networks. Neural Computation, 19(1):111111) January 2007.
[34] Rupesh Srivastava, Klaus Greff, and Jürgen Schmidhuber. Training Very Deep Networks. In Advances in Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015.
[35] Julian Georg Zilly, Rupesh Kumar Srivastava, Jan Koutnik, and Jürgen Schmidhuber. Recurrent Highway Networks. In Proceedings of the 34th International Conference on Machine Learning, pages 4189 4198. PMLR, July 2017.
[36] Piotr A. Sokół, Ian Jordan, Ehen Kadile, and II Memming Park, Adjoint Dynamics of Stahle Limit Cycle Neural Networks. In 2019 53rd Asilomar Conference on Signals, Systems, and Computers, pages 884-887, November 2019. ISSN: 2576-2303.
[37] Piotr A. Sokól. Geometry of Learning and Representations in Neural Networks. PhD thesis, Stony Brook University, May 2023.
[38] Xavier Glorot, Antoine Bordes, and Yoshua Bengio. Deep Sparse Rectifier Neural Networks. volume 15, pages 315-323, April 2011.
[39] Ben Poole, Subhaneil Lahiri, Maithra Raghu, Jascha bick-Dickstein, and Surya Ganguli. Exponential expressivity in deep neural networks through transient chaos. 1006::6666.033440 [cond-mat, statj, June 2016. arXiv: 1606.05340.
12

---

[40] Minmin Chen, Jeffrey Pennington, and Samuel S. Schoenholz. Dynamical Isometry and Mean Field Theory of RNNs: Gating Enables Signal Propagation in Recurrent Neural Networks. N.KYv:8806.05333 [cs, stat], August 2018, arXiv: 1806. 853944
[41] Boris Hanin and Mihai Nica. Products of Many Large Random Matrices and Gradients in Deep Neural Networks. December 2018.
[42] Samuel S. Schoenholz, Justin Gilmer, Surya Ganguli, and Jascha Sohl-Dickstein. Deep Information Propagation. November 2016.
[43] Boris Hanin and David Rolnick. How to Start Training: The Effect of Initialization and Architecture. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018
[44] Piotr A. Sokol and II Memming Park. Information Geometry of Orthogonal Initializations and Training. Technical report, October 2018. ADS Bibcode: 2018arXiv1 1110378511 Type: article.
[45] Dar Gilboa, Bo Chang, Minmin Chen, Greg Yang. Samucl S. Schoenholz, Ed H. Chi, and Jeffrey Pennington. Dynamical Isometry and a Mean Field Theory LSTMs and GRUs. January 2019.
[46] Tankut Can, Kamesh Krishnamurthy, and David J. Schwah. Gating creates slow modes and controls phase-space complexity GRUs and I.S.M. arXiv:2002.00 00025 [cond-mat, stat], January 2020. arXiv: 2002.00025.
[47] Darshil Doshi, Tianyu He, and Andrey Gromov. Critical Initialization Wide and Deep Neural Networks using Partial Jacobians: General Theory and Applications. November 2023.
[48] Valery lustinovich Oseledets. A multiplicative ergodic theorem. Characteristic jjapunov, exponents of dynamical systems. Trudy Moskovskogo Matematicheskogo Obshchestva, 19:179-210, 1968.
[49] Karlheinz Geist, Ulrich Parlitz, and Werner Lauterborn. Comparison of Different Methods for Computing Lyapunov Exponents. Progress of Theoretical Physics, 83(5)8755-83, May 1990.
[50] Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli, and Jean-Marie Strelcyn. Lyapunov Characteristic Exponents for smooth dynamical systems and for hamiltonian systems; A method for computing all of them. Part 2: Numerical application. Meccanica, 15(1):21-30, March 1980.
[51] Hai-Jun Liao, Jin-Guo Liu, Lei Wang, and Tao Xiang. Differentiable Programming Tensor Networks. Physical Review X, 933:033144,, September 2019. arXiv: 1903.09600 [cond-mat, physics.quntt-hh.
[52] S. F. Walter and L. Lehmann. Algorithmic Differentiation of Linear Algebra Functions with Application in Optimum Experimental Design (Extended Version). Technical report, January 2010. ADS Bibcode: 1010arXiv1001.
Type: article.
[53] Robert Schreiber and Charles Van Loan. A Storage-Efficient SWYS Representation for Products of Householder Transformations. SIAM Journal on Scientific and Statistical Computing, 10(1):53 57, January 1989.
[54] Matthias Sceger, Asmus Hetzel, Zhenwen Dai, Eric Meissner, and Neil D. Lawrence. Auto-Differentiiaiig Lincar Algebra. Technical report, October 2017. ADS Bibcode: 201 7ar77ai111000 Type: article.
[55] Jeffrey Pennington. Samuel Schoenholz, and Surya Ganguli. Resurrecting the sigmoid deep learning through dynamical isometry: theory and practice. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
[56] Jeffrey Pennington, Samuel S. Schocnholz, and Surya Ganguli. The Emergence of Spectral Universality in Deep Networks. 1eXiv:8002.09000 (cs, stat), February 2018. arXiv: 1802.09979.
[57] E. Cheney and David Kincaid. Numerical Mathematics and Computing. Cengage Learning, August 2007. Google-Books-ID: JUVZELIMMEC.
[58] Yasaman Bahri, Jonathan Kadmon, Jeffrey Pennington, Sam S. Schoenholz, Jascha Sohl- Dickstein, and Surya Ganguli. Statistical Mechanics of Deep Learning. Annual Review of Condensed Matter Physics, 11()):011-52,, 2020.
[59] F. Ginelli, P. Poggi, A. Turchi, H. Chaté, R. Livi, and A. Politi. Characterizing Dynamics with Covariant Lyapunov Vectors. Physical Review Letters, 130601, September 2007.
[60] Sergey V. Ershov and Alexei B. Potapov. On the concept of stationary Lyapunov basis. Physica D: Nonlinear Phenomena, 118(3):167-198, July 1998.
13

---

[61] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
[62] Javed Lindner. Investigating the exploding and vanishing gradients problem with Lyapunov exponents. Master's thesis, RWTH Aaachen, Aaachen/Juelich, 2021.
[63] Ryan Vogt. Maximilian Puelma Touzel, Eli Shlizerman, and Guillaume Lajoie. On Lyapunov Exponents for RNNs: Understanding Information Propagation Using Dynamical Systems Tools. Frontiers in Applied Mathematics and Statistics, 8, 2022.
[64] Kamesh Krishnamurthy, Tankut Can, and J. Schwab. Theory of Gating in Recurrent Neural Networks. Physical Review X. 12(1):01101,, J. 2022.
[65] L. Pontryagin. Mathematical Theory of Optimal Processes: The Mathematical Theory of Optimal Processes. Routledge, New York, 1st edition edition, March 1987.
[66] Daniel Liberzon. Calculus of Variations and Optimal Control Theory;
A Concise Introduction. Calculas of Variations and Optimal Control Theory. Princeton University Press, December 2011.
[67] Owen Marschall, Kyunghyun Cho, and Cristina Savin, unified framework of online learning algorithms for training recurrent neural networks. The Journal of Machine Learning Research, 21(1):135:5320- 135:5353, January 2020.
[68] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is All you Need. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
[69] Judy Hoffman, Daniel A. Roberts, and Sho Yaida. Robust Learning with Jacobian Regularization. Technical report, August 2019. ADS Bibcode: 2019arXiv1900027000 190802729H Type: article.
[70] Hongyi Zhang, Yann N. Dauphin, and Tengyu Ma. Fixup Initialization: Residual Learning Without Normalization. Technical report, January 2019. ADS Bibcode: 2019arXiv 1901093212 Type: article.
[71] Thomas Bachlechner, Bodhisattwa Prasad Majumder, Henry Mao, Gary Cottrell, and Julian McAuley. ReZero is all you need: fast convergence at large depth. In Proceedings the Thirty-Seventh Conference on Uncertainty in Artificial Intelligence, pages 1352- 1361. PMLR, December 2021. ISSN: 2640-3498.
[72] Tianyu He, Darshil Doshi, and Andrey Gromov. AutoInit: Automatic Initialization via Jacobian Tuning, June 2022. arXiv:2206. 13568 [cond-mat, stat].
[73] L. Molgedey, J. Schuchhardt, and H. G. Schuster. Suppressing chaos in neural networks by noise. Physical Review Letters, 69(26)3317-37799 December 1992.
[74] Kanaka Rajan, L. F. Abbott, and Haim Sompolinsky. Stimulus-dependent suppression of chaos in recurrent neural networks. Physical Review E, 82(1):01 1903, July 2010.
[75] Jannis Schuecker, Sven Goedeke, David Dahmen, and Moritz Helias. Functional methods for disordered neural networks. arXiv:1605.000 06758 [cond-mat, q-bio]. May 2016. arXiv: 1605.06758.
[76] Rainer Engelken, Alessandro Ingrosso, Ramin Khajeh, Sven Goedeke, and L.. F. Abbott. Input correlations impede suppression chaos and learning in balanced firing-rate networks. PLOS Computational Biology, 18(12)=10110500 December 2022.
[77] Bezanson, A. Edelman, S. Karpinski, and V. Shah. Julia: A Fresh Approach to Numerical Computing. SIAM Review, (-6-998999, January 2017.
[78] Innes. Flux: Flegant machine learning with Julia, Journal of Open Source Software, 3(25):602, May 2018.
[79] Edward OuL Chaos in Dynamical Systems. Cambridge University Press, August 2002. Google-Books-ID> PIXoAWAAQBAI.
[80] Angelo Vulpiani, Fabio Cecconi, and Massimo Cencini. Chaos: From Simple Models Complex Systems. World Scientific Pub Co Inc, Hackensack, NJ, September 2009.
14

---

A Backpropagation Through QR Decomposition
The backward pass of the QR decomposition is given by 51 52 53 54
Q Q + Q copyltu(MM) R
(11)
where M = RR Q' Q and the copyltu function generates symmetric matrix by copying the lower triangle of the input matrix to its upper triangle, with the element Copyltu(MM)s; M (min()))..mn((..))
52 53
52
54 Adjoint variable are written here as T = 0L/0T.
Using an analytical pullback is more memory-efficient and less computationally costly than directly doing automatic differentiation through the QR-dccomposition. Moreover, from practical perspective, for QR decomposition. often BLAS/LAPACK routines are utilized which programming frameworks like TensorFlow, PyTorch, JAX and Zygote. our implementation gradient flossing, we used the Julia package Backwards by Jinguo Liu available here
amenable to common differentiable
B Further Details and Analysis of Gradient Flossing
An example implementation of gradient flossing in Flux 78 machine learning library in Julia (77) available
hhttps//github..o httpgggthuu..ooRaaiiia wowm/ainerEngeleeeeeee wrreeeeeeeeeeeeennnnniiieeee GradientFlossing We actively developing implementations for other widely used diffcrentiable programming frameworks.
B.1
Gradient Flossing for recurrent LSTM and ReLU networks
LSTM
ReLU
A
c
target A
target À
actual A
actual À
4
0.0
2
-0.2
0.4
-2
20
40
60
80
100
20
60
80
100
Epochs
Epochs
8
D 10a
10 -1
10-2
10-1
10
10-
10-3
10-5
10-6
10-5
100
Epochs
Epochs
Figure 6: Gradient flossing for recurrent LSTM networks and recurrent ReLU networks A) First Lyapunov exponent of LSTM network as a function of training epochs. Minimizing the mean squared error between estimated first Lyapumov exponent and target Lyapunov exponent A1 0 by gradient descent. First Lyapunov exponent LSTM network (solid lines) converges to target value (thick dashed lines) within less than 100 epochs. 10 random LSTM RNNs
initialized with Gaussian recurrent weights, where standard deviations of weight scaling were drawn g ~ Unif(0, 1). B) Gradient flossing minimizes the square of the first Lyapunov exponent of random recurrent LSTM networks over epochs. C) Same as A for recurrent ReLU network. Here networks were initialized with Gaussian recurrent weights Wij N-(0.1, g²
/NN where values of g were drawn ~ Uni 1) D) B) for recurrent ReLU network. Parameters: network size N = 32 with 10 network realizations. Shaded regions in B. are 25% and 75% percentiles, solid line shows median.
15

---

We demonstrate that gradient flossing can also be applied to recurrent LSTM and ReLU networks in Fig 6 To this end, we generated random LSTM networks where the weights of all the different gates and biases were independently and identically distributed (i.i.d)) and sampled from Gaussian distributions of different variance. Our results show that gradient flossing can also constrain the Lyapunov exponent to be close to zero. The dynamics of each of the N LSTM units follows the map [
Se ao(Uphe-1 Wa by)
(12)
Oe dgUUchi-1 Wo + bo)
(13)
b ay(Udu--1 + W.r, + &c)
(14)
č = on(Uchi-1 + Watt + b.)
(15)
Ce = 0 is @ à
(16)
he or @ ((ee))
(17)
where denotes the Hadamard product, ox(x) 1-mxp(-)) the sigmoid function, (x((2) tanh(z) and
entries of the matrices U2 are drawn from U2 ~ NN / For simplicity, the bias terms b, are scalars. Subscripts f. o and denote respectively the forget gate, the output gate, the input gate, and i the cell state. In each LSTM unit, there are two dynamic variables and h, and three gates f. o, and that control the flow if signals into and out of the cell c. We set the values Jch, Gir. If2, by, John Jen, gor, Joz to be uniformly distributed between and and initialize b, br.bb as zero.
During gradient flossing. the actual Lyapumov exponents different random network realizations converge close to the target Lyapunov exponent X hinget = in fewer than 100 epochs as shown in Fig Fig 6 shows that the squared Lyapunov exponents converge towards zero. We note that for LSTM networks, a target Lyapunov exponent of Xurger -1 is achieved after 100 gradient flossing steps only for subset of random network realizations (not shown). We speculate that behavior influenced by the gating structure of LSTM units, which seems to naturally place the first Lyapunov exponent close to zero for certain initializations (See also
p 12 46 64
For the recurrent ReL.U networks, we considered the same Vanilla RNN dynamics as in the main manuscript in Eq 9
h,+1 [[hh,x-)) We(h.) + Vx.+1,
The initial entries of W are drawn independently from Gaussian distribution with negative mcan of -0.1 and variance g² /N, where g is gain parameter that controls the heterogeneity of weights. We use the transfer function (xx) max(z, 0). x, a sequence of inputs and V the input weight. x, is a stream of i.i.d. Gaussian imput xr ~ (0, 1) and the input weights V are (O, 1). Both W and V are trained during gradient flossing. We found that some ReLU network had initially unstable dynamics with positive Lyapunov exponents bC. However, during gradient flossing, these unstable networks were quickly stabilized. Fig shows that the squared Lyapunov exponents of ReLU networks converge towards zero
B.2 Additional Results for Different Numbers of Flossed Lyapunov Exponents
Additionally to the main Fig 5 we did the same analysis networks with only preflossing (gradient flossing before training) and found that more Lyapunov exponent & have to be flossed to achieve 100% accuracy as the tasks become more difficult P)), however, in that case even if all N Lyapunov exponents were flossed, thus k = N, they were not able to solve the most difficult tasks. This seems to indicate that gradient flossing during training cannot be replaced by just gradient flossing more Lyapunov exponents before training.
C Additional Results on Gradient Flossing Throughout Training
We now discuss some additional results gradient flossing throughout training. First, analyze how gradient flossing affects the gradients and find that during gradient flossing, norm of gradients that bridge many time steps are boosted. Morcover, subordinate singular values of the error norm of the recurrent weights are also boosted, indicating that gradient flossing can increase the effective rank of the parameter update. Additionally, we show that if gradient flossing is continued throughout training it can be detrimental to the accuracy. Finally, we show that Lyapunov exponents of successfully trained networks after training for the spatial delayed XOR task have simple relationship to the delay d.
D Gradient Flossing boosts the Gradient Norm for Long Time Horizons
In this section, we investigate the impact of gradient flossing on the norm and structure of the gradient. It is important to note that the complete error gradient of backpropagation through time is composed of summation of products of step Jacobians, reflecting the number of "loops" the error signal traverses through the recurrent
16

---

A
B
100
100
test accuracy (%)
test accuracy (%)
k
delay
80
80
20
30
40
50
60
70
80
20
40
60
0
20
40
60
80
complexity (delay d)
number of flossed A
C
100
D
100
70
test accuracy (%)
test accuracy (%)
(t (dets lelay
60
60
comp2 hopleextyyy (de
50
50
40
40
30
20
20
10
10
50
1 10 20
30 40 50 60 70 80
1 10 20 30 40 50 60 70 80
number of flossed A; k
number of flossed A; k
Figure 7: Gradient flossing for different numbers of flossed Lyapunov exponents A) Test accuracy for delayed temporal XOR task as function of delay with different numbers flossed I Lyapunov exponents k. B) Same data A but here test accuracy as function of number of flossed Lyapunov exponents k. C) Median test accuracy for delayed temporal XOR task as function of delay d and & for networks with gradient flossing during training (500 steps of gradient flossing at epochs e E {0, 100, 200, 300, 400)). D)Same as for preflossing only. Parameters: g = 1, batch size b 16, N 80, epochs = 10 for delayed temporal XOR, epochs 5000 for delayed spatial XOR, T = 300, gradient flossing for Ey = 500 epochs before training and during training for A, B, C, and only before training for C. Shaded areas are 25% and 75% percentiles, solid lines are means, transparent dots are individual simulations, task loss: cross-entropy btw. y, g.
dynamies before reaching its target. Consequently, when the singular values of the long-term Jacobian smaller than the influence of the shorter loops typically dominates the long-term Jacobian.
In tasks, have full control over the correlation structure of the task and thus know exactly which loop length of backpropagation through time necessary finding the correct association. We were moreover careful in our task design not to have any additional signals in our task that might help to bridge the long time scale. In the case of vanishing gradients, the gradient norm is predominantly influenced by the shorter loops. even though the actual signal in the gradient originates solely from the loop of length d in our task. To mitigate the contamination of spurious signals from shotter loops and effectively extract the gradient that spans long time horizons, we focus on the gradient with respect to the initial conditions ha-
aC. aC. Ph. Jh.
Σ
II
Sh oh
oh, ac. Σ Th.
I
Th Oh
8,
Th. Oh.
Tho))
(18)
Th.
We note that the sum conveniently drops as only the longest "loop", in other words, the only summand that contributes is the product of Jacobians going from 0 to f. By considering this gradient, can therefore ensure that no undesired signals stemming from shorter loops interfere with the analysis. Moreover, use the binary cross entropy loss which makes the derivative 2ct trivial. Th,
that
In Fig g we show that gradient flossing boosts the gradient with respect to the initial conditions. Specifically, we compare two identical networks trained on the binary delayed temporal XOR task with loop length of 70.
17

---

One network is trained with gradient flossing at epochs E {0, 200,1000,0000000000 200, 400)), while the other is trained without gradient flossing.
For the network without gradient flossing. the gradient norm of diminishes to extremely small values (< 10-9) and remains small throughout training. In contrast, for the network trained with gradient flossing, each episode of gradient flossing causes the norm 36 aho to spike, surpassing values larger than 10-2. These findings direct evidence that gradient flossing boosts the gradient norm, facilitating bridge long time horizons in challenging temporal credit assignment tasks. We observe that after several episodes of gradient flossing. the gradient latg of the networks stays around 10 and eventually rise up to values around 10-2 Subsequent in training, the test accuracy surpasses chance level (B) We observed this temporal relationshin between gradient norm de and training success consistently across numerous network realizations (Fig and D).
de
de can be a good predictor of learning success, sometimes These findings suggest that the gradient norm ang
hundreds of epochs before the accuracy exceeds the chance level of 50%. Indeed, when depicting the gradient norm aligned to the last epoch where accuracy was < 50%, we see for many network realizations a gradual growth of gradient norm oven epochs before accuracy surpasses chance level (.). Analogously, when plotting the accuracy as a function of epoch aligned with the last epoch with abo < 0.001, we observe for this ålg
task that the increase of gradient norm dé reliably precedes the epoch which the accuracy surpasses the chance level (See Fig(B)) We note that when measuring the overlap of the orientation of the gradient vector de with the first covariant Lyapunov vector of the forward dynamics, we found a significant increase in overlap Tho
around the training epoch where the accuracy surpasses the chance level both in networks with and without gradient flossing. This does not come a surprise as the covariant Lyapunov vector measure the most unstable (or least stable) direction in the tangent space of a trajectory and perturbations of tn that have to travel over many epochs align
D.1 Gradient Flossing Boosts Effective Dimension of Error Gradient
To further investigate the effect of flossing training. investigated the structure of the error gradient and how it is changed by gradient flossing. To this end, we decompose the recurrent weight gradient o; dC into in weighted sum of outer products using singular value decomposition (Fig TT).
As the Lyapunov exponents are the time-averaged logarithms the singular values of the asymptotic long-term Jacobian (h, ). this allows us to directly link the effect of pushing Lyapunov exponents toward zero during gradient flossing to the structure of the error gradient of the recurrent weights, as they are intimately linked:
ac, OW
aC,
T=(-1
Σ
II
0hr4+1 Th
Oh. JW
aL. Σ T (h+)
Oh, OW
(19)
Th.
Th.
r---1
We again note that different 'loops' contribute to the total gradient expression and the Lyapunov exponents only characterize the longest loop. Further, we note that in our controlled tasks, depending on delay d. only few of the summands are relevant for solving the task. We thus expect the relevant gradient summand that carries important signals about the task to be contaminated by summands of both shorter and longer chains, which contribute irrelevant fluctuations.
The singular values of the recurrent weight gradient 0, ( dC (aw)
function of training epoch reveal that the subordinate singular values subordinate singular value azo and au exhibit peaks the times of gradient flossing, while the first singular value di only shows slight peak (Fig [TP). This indicates that gradiens flossing increases the effective rank of the recurrent weight gradient dc In other words, gradient flossing facilitates high-dimensional parameter updates. Our interpretation is, gradient flossing pushes Lyapunov exponents zero, the different summands in the total gradient contribute more equitable as long loops have neither dominant contribution (which would happen for exploding gradients) nor a vanishing contribution (which would happen for vanishing gradients). This way. the sum of gradient terms has a higher effective rank.
zw
In contrast, without gradiens flossing, the subordinate singular values (in Fig 10 d'e0 and a rapidly diminish to extremely small values over training epochs and remain very small throughout training. Note however that the leading singular values d are of comparable size irrespective whether gradient flossing was performed or not.
We note that similar to the gradient norm of the loss with respect the initial condition de J, the subordinate singular values seem to predict when the test accuracy of networks with gradient flossing grows heyond chance level (Fig IOB). We confirmed this in multiple other network realizations and give here another example we the accuracy grows beyond chance only later during training (Fig IT).
D.2
Gradient Flossing Throughout Training Can Be Detrimental
We find that gradient flossing continued throughout all training epochs can be detrimental for performance We demonstrate this again in the binary delayed temporal XOR task. We compare three different
(Fig 12
18

---

10
0.05
10-4
0.04
10-9
0.03
86
10-8
0.02
10
with flossing during training without flossing
0.01
0.00
500
1000
1500 Epochs
2000
2500
3000
500
1000
1500
2000
2500
3000
Epochs
100
100
accuracy (%)
Ceuunnne (6)
with flossing during training without flossing
500
1000
2000 2500
500 1000
1500
2000
3000
Epochs
Epochs
Figure Gradient flossing boosts norm of long-term Jacobian A) Gradient norm of dC
dho as function of training epochs for networks without flossing (blue) and networks with flossing during training (orange). Error gradient norm is boosted after gradient flossing at epochs e E {0, 100, 200, 300, 400, 500)). In networks without gradient flossing, the gradient norms dC much smaller overall. One out of ten random network realizations with solid line, the other 9 with transparent line. B) Accuracy
dio are
function of epoch, same depiction and network realizations in A. Note that accuracy of networks with gradient flossing grows beyond chance level approximately when the gradient norm dhg becomes macroscopically large. C) Same as A in linear scale. Mean final test loss as a function of task difficulty (delay d) for delayed copy task. Different colors are different network realizations with gradient flossing during training. Black lines are without any gradient flossing. D) Accuracy function epochs, same colors as C. Note that for all network realizations the moments where gradient norm dho becomes macroscopically large coincides with 1, batch size b 16, N 80, epochs = 10 T = 300, flossing for E = 500 epochs on k = 75 Lyapunov exponents before training. Task: binary delayed XOR, delay d = 70, loss: cross entropy(y, 9)
the moment the accuracy beyond chance level. Parameters:
conditions: Either, we floss throughout the training every 100 training epochs for 500 flossing epochs (red), or we floss only early during training at training epochs e E (0, 100, 200, 300, 400) 400ee))reen) or we do not floss at all (bluc).
We observe that after every episode of gradient flossing. the accuracy drops down close to chance level of 50%
(Fig A2A red line). Between flossing, the accuracy quickly recovers but never reaches 100%. Simultaneously. when the accuracy drops the test error jumps up (Fig 12B). We also observed that the gradient
Who
stays close to indistinguishable once the gradient norm Tho macroscopically large (Fig [2C). This suggests that once gradient flossing facilitates signal propagation across long time horizons and the network picks up the relevant gradient signal, further gradient flossing be harmful to the actual task execution. We hypothesize that there might (at least for the Vanilla networks considered here), trade-off between the ability to bridge long time scales which seems to require one or several Lyapunov exponents of the forward dynamics close zero and nonlinear tasks requirements, which require at least fraction of the units be in the nonlinear regime of the nonlinearity é. where ø' (xx) < It would be an interesting future research avenue further investigate this potential trade-off also other network architectures.
initially boosted by gradient flossing.
becomes
D.3 Lyapunov Exponents after Training With and Without Gradient Flossing
B) and the tenth (Fig 13C. D) Lyapunov exponent after training on the spatial delayed XOR task both with and without gradient flossing. We find for successful networks with
In Fig we show the first (Fig 13
19

---

A
B
80
75
10
70
accuracy (%6
10
65
60
55
5
50
45
-3000
-2000
-1000
0
1000
0
100
200
300
400
Epochs
Epochs
Figure 9: Increase of gradient norm precedes epoch when accuracy exceeds chance level
de A) Gradient norm of dfg as function of training epochs for 20 network realizations with flossing during training. Epochs are aligned to last epoch where accuracy is < 50%. B) Same task and simulations as in A, here accuracy as a function of epoch, for 20 network realizations with flossing during training. Epochs are aligned to the last epoch with dC
Thg 0.001. Different colors are different network realizations. Parameters: = 1, batch size b = 16, N = 80, epochs = 10 T = 300, flossing for F. 500 epochs k = 75 Lyapunov exponents early during training training epochs e E {0, 100, 200, 300, 400). Task: binary delayed XOR, delay 70, loss: cross entropy(y, 9)
gradient flossing a systematic relationship between first Lyapunov exponent and delay, can fitted by approximately by A (d) = -..2eap. (-.00ddela). Unsuccessful networks with accuracy at chance level have a much smaller largest Lyapunov exponent. The same seems to hold true for the tenth Lyapunov exponent. In previous study [3 similar trend was observed, albeit in the context of task that did not possess an analytically tractable temporal correlation structure, which might partially explain the less conclusive results. is important to note that numerical evaluation of Lyapunov exponents recurrent LSTM networks 63 was based solely on the N x N Jacobian of the memory state. From a dynamical systems standpoint, a 2N x 2NV Jacobian matrix encompassing interactions between both memory and cell states into account is
required 9 12 64)
Gradient Flossing for Linear Network
We provide code for gradient flossing in linear networks We find that gradient flossing also helps train linear networks on tasks many time steps that can be solved by linear networks, for example the copy task, not tasks require a nonlinear input-output operation like temporally delayed XOR task. Full analytical description of gradient flossing for linear networks would be promising avenue for future research as networks with linear dynamics can still have nonlinear learning dynamics X However this is heyond the cope of the presented work.
Computational Complexity of Gradient Flossing
We present here main contributors to the computational cost [ First the RNN step, which has computational complexity of (N' b) per time step. where N is the dimension of the recurrent network state (which in case of Vanilla networks equals the numher of units) and b is the batch-size both in the forward and backward pass. Second, the Jacobian step which scales with O (N² per time step, where & is the number of flossed Lyapunov exponents. Third, the QR decomposition, which scales with O (N 12). where & is the number of Lyapunov exponents considered.
in in-depth scaling analysis of the computational cost of gradient flossing. There three
20

---

A
2
ide BLOR
0
500
1000
1500
2000
2500
3000
B
early training flossing
accuracy (%)
without flossing
50
0
500
1000
1500
2000
2500
3000
Epochs
Figure 10: Gradient flossing decreases condition number of recurrent weight error gradient A) Singular values of recurrent weight gradient ( function of training epochs for singular values E 1, 20, 40 for networks without gradient flossinn (blue) and early training gradient flossing (green). At epochs of gradient flossing, the subordinate singular value 2n and 4 are peaked, while the first singular value σj has only a slight peak. This indicates that gradient flossing increases the effective rank of the recurrent weight gradient ZW B) Accuracy as a function of training epochs. Note that accuracy of networks with gradient flossing grows beyond chance level approximately when the subordinate singular values singular value Ø and σ are peaked increase, which enables high-dimensional parameters updates. Parameters: g = 1, batch size 16, N 80, epochs 103. = 300, gradient flossing for E5 = 500 epochs on & = 75 Lyapunov exponents before training. Task: binary delayed XOR, delay d = 70, loss: cross entropy(y. 9)
Together, this results in a total amortized cost of (N² 61) per training epoch, where is the number of training time steps and total amortized costs per flossing epoch (N*T,(1 k/tons + k)) where T'y the number of flossing time steps.
In case of preflossing, thus, the total computation cost scale with (N2 + E, (1 (1 k/toNs + k))).
where E is the number of training epochs and Ep is the number of prefiossing epochs.
For gradient flossing during training (assuming that there is also prefiossing done), the amortized cost scale with O [N3 F.p T F, Ty k/toNs k)]), where F, is the total number of flossing epochs during training.
Empirically, find that both the number of preflossing epachs F.p and flossing episodes Fy necessary for training success is much smaller than the total number training epochs E. For example, the preflossing for 500 epochs in the numerical experiment of Fig book ~ 37 seconds, while the overall training on 10000 training epochs with batch size 6 16 took 1680 seconds. Thus only approximately 2.2% of the total training time spent gradient flossing. Moreover, Тр be smaller than T. it just has be long enough such that the temporal correlations in the task can bridged. In of the tasks discussed the manuscript, this would be the delay d. It remains an important challenge to infer the suitable number of flossing time steps Ty for tasks with unknown temporal correlation structure.
would also be interesting to investigate how the CPU hours/wall-clock hitee//loo//Julle/2-eemmmss spent gradient flossing vs training networks with larger N are trading off against each other. For this, we would suggest to first find the smallest network that median successfully trains binary temporal XOR task for a fixed given delay and measure
computational resources involved training it, e.g. in terms CPU hours. Then compare it to a network with g gradient flossing. This would be a promising analysis but is beyond our current computational budget. We will start such experiments might be able to provide results during the reviewer period.
21

---

A
2
ide BLOR
0
500
1000
1500
2000
2500
3000
B
early training flossing
accuracy (%)
without flossing
50
0
500
1000
1500 Epochs
2000
2500
3000
Figure 11: Gradient flossing decreases condition number of recurrent weight error gradient
Same as Fig 10 for different network realization.
G Additional controls
We also investigate the effects of gradient
flossiin flossing during the training with orthogonal weight initializations and confirm our finding that gradient flossing improves trainability on tasks that have long time horizons to bridge. Moreover, find that gradient flossing during training can further improve trainability. We replicated the challenging tasks from the main paper (Fig 4 for orthogonal initialization with variable temporal complexity and performed gradient flossing either both during and before training. only before training. or not
all.
Fig I4A shows the test accuracy for Vanilla RNNs with orthogonal initialization trained on the delayed temporal
XOR task = 2t-d/2 ( z with random Bernoulli process E {0. 1). The accuracy of orthogonal Vanilla RNNs falls to chance level for d > 40 (Fig [4 With gradient flossing before training, the trainability
be improved, but still falls close chance level for 70. In contrast, for initially orthogonal networks with gradient flossing during training, the accuracy improved to > 80% at 70. In this case, we preflossed for 500 epochs before task training and again after 500 epochs of training on the task. In Fig 14B, D the networks have to perform the nonlinear XOR operation signal wt, x², and n and generate the correct output with a delay of steps identical to Fig in the main text. Again, we observe similar to networks with Gaussian initialization that flossing before training improves the performance compared baseline, but starts failing for long delays d 60. In contrast, orthogonal networks that are also flossed during training solve more difficult tasks (Fig [4p). We note that for Fig [4] and D. we trained network only on 5000 epochs, compared to 10000 epochs networks with random Gaussian initialization because for 10000 epochs, both networks with gradient flossing only before training and with gradient flossing before and during training were able to bridge d = 70. These results suggest that orthogonal initialization does to slightly improve performance for tasks with long time horizons bridge and gradient flossing and additionally boost performance. Thus orthogonal initialization well together. It would be interesting study orthogonal initialization also reduces the number gradient flossing steps necessary to improve performance.
= ① _.i ① s²-a on three-dimensional binary input
gradient flossing seems to go
H Additional Details on Training Tasks
In this section, provide rigorous definition of the tasks used for training, discussed in Section 3:
22

---

A
accuracy (%)
100
A
with flossing throughout training early training flossing without flossing
50
0
500
1000
1500
2000
2500
3000
test error
100
10-1
500
1000
1500
2000
2500
3000
WM
10-3
with flossing throughout training early training flossing without flossing
de
10-7
0
500
1000
1500
2000
2500
3000
Epochs
Figure 12: Gradient flossing throughout training be detrimental to learning A) Accuracy a function of training epochs for binary temporal delayed XOR task for gradient flossing throughout training every 100 training epochs (red). Accuracy drops down close to chance level every time after gradient flossing but recovers quickly between. Same for only 5 episodes of gradient flossing at epochs e {0, 100, 200, 300, 400)) (green) and no flossing at all (blue). B) Test error as a function function of training epochs for networks without gradient flossing (blue) and networks with flossing throughout training (red) and early training gradient flossing (green). Error gradient norm is boosted after each gradient flossing. In networks without gradient flossing, the gradient norms dC are much smaller overall. Parameters: = 1, batch size b = 16, N = 80, epochs = 10 T 300, gradient flossing for F, 500 epochs on 75 Lyapunov exponents before training. Task: binary delayed XOR, delay d = 70, loss: cross entropy(y, 9).
of training epochs. C) Gradient norm of dC
dho
dho
H.1
Copy task
For the copy task, the target network readout at time i is yı -i-d, where denotes the delay. We chose the input to be sampled i...d. from a uniform distribution between 0 and 1.
H.2
Temporal XOR task
The temporal XOR task requires the target network readout yt at time to be computed as follows:
|rr- 22-al
(20)
where again d denotes a time delay of d time steps. In the case of æ E {0, 1}2 and E {0, 1}, the output y follows the truth table of the XOR digital logic gate (Table Thus, the function f(e., 2a.) = Iz. nal can be seen an analytical representation of the XOR gate. It important note that f(z, 0) æ only for implementation can easily be constructed analytically, for example.
this task requires nonlinearity. using two rectified linear units ((x) = max(x, 0) the outbut can constructed by
f( vb) Irr sal = Φ(za mτm) + d(0, .a).
(21)
Together with delay line transmit the signal æ:-á time, this solve the task.
23

---

no gradient flossing
gradient flossing during training
A
B
-0.1
-0.1
-0.2
-0.2
20
40
60
20
40
60
complexity (delay d)
complexity (delay d)
C
D
-0.1
0.1
Aι0
-0.2
Aιo
-0.3
-0.2
20
40
60
20
40 complexity (delay d)
complexity (delay d)
Figure 13: Lyapunov exponents of trained networks with and without gradient flossing A) First Lyapunov exponents Y for Vanilla networks trained on spatial delayed XOR task function of the delay with no gradient flossing. Colored-coded is accuracy at the end of training where red corresponds to 100% accuracy and blue to chance level (50%). B) Same as A for networks with gradient flossing during training. Black dashed line shows that Lyapunov exponents of successfully trained networks can be approximated by the empirical fit Y₃ (d) = - 0.2ezz. (-0.0delayy). (Proto- col for gradient moudientflossiinge flossing during training same as main text Fig B). ). C) Same as for tenth Lyapunov exponents λισ. D) Same as B for tenth Lyapunov exponents λιo. Same fit as in B also describes Å1o- Parameters: g = 1, batch size = 16, E 500 epochs on k = 75 Lyapunov exponents before training. Task: binary spatial delayed XOR, loss: cross entropy(y. 3)
80, epochs 10 T 300, gradient flossing for
I Additional Background on Lyapunov Exponents of RNNs
An autonomous dynamical system is usually defined by of ordinary differential equations dh/dt F(h), h E R the case of contimuous-time dynamics, or as a map h,++1 = f(h,) in the case of discrete-time dynamics. In the following. the theory is presented for discrete-time dynamical systems for ease of notation, but everything directly extends map forms a trajectory. As natural extension of linear stability analysis, one can ask how an infinitesimal perturbation há = ho + cuto evolves in time. Chaotic systems are sensitive to initial conditions; almost all infinitesimal perturbations eug of the initial condition grow exponentially with time |eu, AS exp(Art) ( Finite-size perturbations, therefore, may lead to drastically different subsequent behavior. The largest Lyapunov exponent Aı measures the average rate of exponential divergence convergence of nearby initial conditions:
antinuuus--ine systems 49 Together with initial condition ha, the
|eue]] cuu
A1 (ho) lim lim log
(22)
C-000
In dynamical systems that ergodic the attractor, the Lyapunov exponents do not depend the initial conditions as long as initial conditions are in basins attraction the attractor. Note to first take the limit e - 0 and then É -> co, as Aı (ho) would be trivially zero for a bounded attractor if the limits are exchanged, as lim log ||euell bounded finite perturbations even the system is chaotic. To Teuall
crucial
measure Lyapunov exponents, one to study the evolution of independent infinitesimal perturbations u. spanning tangent space:
Us+1 Daus
(23)
24

---

forward pass
backward pass
RNN dynamics
O ( b)
Jacobian step
0 (N2 k)
QR step
0 (N k²)
total amortized costs per training epoch
O (T)
total amortized costs per grudient floxsing epoch
o (N² k/LoNs 1:)
total amortized costs of preflossing
o (N2 + F.p (y + k/toNs + kjl)
total amortized costs flossing during training
O (N2(EbT + E, T + By
k))
Table 1: Computational cost for gradient flossing and training of RNNS
N denotes number of neurons, the batch size, T is the number of time steps in forward pass of training, T is the number of time steps in forward pass of flossing, loNs is the reorthonormalization interval, & the number of flossed Lyapunov exponents, E is the number of training epochs, Ep is the number of preflossing epochs, Ey is the number of flossing epochs during training. Empirically, we find that the necessary number of preflossing epochs Ep and flossing episodes E is much smaller than both the total number of training epochs E. Moreover, Tр can be smaller than T.
Table 2: XOR
input 2z-d
input Tt-26
target output y/+
0
0
1
I
1
1
where the N x N Jacohian D.(h.) df (h.)/hh characterizes the evolution of generic infinitesimal perturba- tions during step. Note that this Jacobian along the trajectory equivalent stability matrix only fixed point, i.e. when h f(h.) h..
We are interested in the asymptotic behavior, and therefore we study the long-term Jacobian
Te (bo) = Di-1 (hr- D. (hi Dooh))
(24)
Note that T, (bo) is a product of generally noncommuting matrices. The Lyapunov exponents À A2... are defined as the logarithms of the eigenvalues the Oseledets matrix
Ax
A(ho) = lim (To(ho) Tho)) t
(25)
where T denotes the transpose operation. The expression inside brackets is the Gram matrix of the long-term Jacobian T, (ho). Geometrically, the determinant of the Gram matrix the squared volume of the parallelotope spanned by the columns of T. (ho). Thus, the exponential volume growth rate is given by the sum of the logarithms of its first & (sorted) eigenvalues. Oscledets' multiplicative ergodic theorem guarantees the existence of the Oscledets matrix A A(ho) for almost all initial conditions ha [48). In ergodic systems, the Lyapunov exponents X, do not depend on the initial condition hy. However, for a numerical calculation of the Lyapunov spectrum, Eq cannot be used directly because the long-ten ten-termmaa Jacobian Tr (ho) quickly becomes ill cl-condiiioned, i.e., the ratio between its largest and smallest singular value diverges exponentially with time.
Algorithm for Calculating Lyapunov Spectrum of Rate Networks
For calculating the first & Lyapunov exponents, we exploit the fact that the growth rate of a k-dimensional infinitesimal volume element is given by (m) Σm A. Therefore, A -)), A2 ( A1, As
(3) d1 A2. 50 The volume growth rates can be obtained via decoomompnnntin.
25

---

orthogonal weight initialization
delayed temporal XOR
delayed spatial XOR
B
100
100
flossing during training
test accuracy (%)
test accuracy (%)
prefiossing
80
flassing
80
60
60
2000
4000
6000
10000
5000
Epochs
Epochs
100
100
test accuracy (%)
test accuracy (%)
80
80
flossing during training preflossing no flossing
10
30
50
complexity (delay d)
complexity (delay d)
Figure 14: Gradient flossing before and during training improves trainability for orthogonal nets
A) Test accuracy for orthogonally initialized vanilla RNNs trained delayed temporal binary XOR task y1 -t-4//2 O It-d with gradient flossing during training (green), preflossing (orange), and with no gradient flossing (blue) for 70. Solid lines are mean, transparent thin lines are individual network realizations B) Same as A for delayed spatial XOR task with yı = x @ x + x² t_d function of task difficulty (delay for delayed temporal XOR task. D) Test
C) Test aceuracy accuracy as a function of task difficulty (delay d) for delayed spatial XOR task. Parameters: 80, epochs 10 for delayed temporal XOR, epochs 5000 for delayed
1,
batch size b = 16, spatial XOR, T = flossing for E, 500 epochs on h = 75 Lyapunov exponents before training and during training for green lines, and only before training for orange lines. Shaded areas are 25% and 75% percentiles, solid lines are means, transparent dots are individual simulations, task loss is cross-entropy between y. ý-
First, we evolve an initially orthonormal system Q. = (q/, qe....q"" the tangent space along the trajectory using the Jacobian D.:
q."
Q++1= = D.Q.
(26)
continuous system can be transformed into a discrete system by considering a stroboscopic representation, where the trajectory is only considered at certain discrete time points. We use here the notation of discrete
dynamical systems where this corresponds to performing the product of Jacobians along the trajectory Q
D. Q.Q We study the discrete network dynamics in the limit of small time step At - 0 and for discrete time Al 1. The notation be readily extended continuous systems (49
Second, extract the exponential growth rates using the QR-decomposition,
Q.+1
Q++IR*1!
which uniquely decomposes Q.-1 into an orthonormal matrix Q.+1 size N x Q.+ 1mxm and to an upper triangular matrix R. of size & x & with positive diagonal elements. Geometrically, Q.+1 describes the rotation of Q. caused by D, and the diagonal entries of R describe the stretching and shrinking of columns of Q.. while the off-diagonal elements represent the shearing. Fig 15 visualizes D, and the QR deooppoopon for & = 2.
The Lyapumov exponents are given by time-averaged logarithms of the diagonal elements of R
A lim to
log
R s=1
lim <-000
log R
(27)
f
26

---

D2
QR
q
9
R
6+1
s++1
s+1
R2
Figure 15: Geometric illustration of Lyapunov spectrum calculation. An orthonormal matrix whose columns the axes of k-dimensional cube, is rotated and
Q. q ......."" ",",, distorted by the Jacobian D, into an k-dimensional parallelotope Q.+1 D,Q, embedded in RN, The figure illustrates this for = 2, in which case the columns Q.+1 span parallelogram, which can be divided into right triangle and trapezoid and rearranged into rectangle. Thus, the area of the gray parallelogram is the same as that of the orange rectangle. The QR-decomposition
reorthonormalizes Qs+1 1 by decomposing it into the product of an orthonormal matrix Qs+1 = [q 4 9 il and the upper-triangular matrix R*+1, Q++1 describes the rotation of Q. caused by D... The diagonal entries of R*+1 gives the stretching/shrinking along the columns of Qst1. thus the volume of the parallelotope formed by the first k columns of Qs+1 is given by Vm -m R The time-averaged logarithms of the diagonal elements of R give the Lyapunov spectrum: A lim
log IT. R lim 400 ==1 log R. Sin
Note that the QR QR-cccompoiition does not need to be performed at every simulation step, just sufficiently often, .. D. Q2 remains cell-conditioned
ie., once every sons steps such that a = DD++sons DD. ...
(50) An appropriate reorthonormalization interval soNs = tons / At thus depends on the condition number, the ratio of the smallest and largest singular value:
01 (R +
R FONS
K2g ac+noNs ) = K2 FONNS
(28)
0m ++toss
R
An initial transient should disregarded in the calculation the Lyapunov spectrum because first has to converge towards the attractor and Q has to converge to the unique eigenvectors of the Oseledets matrix (Eq25)
25
660 A simple example of this algorithm in pseudocode is:
Algorithm 2 Jacobian- based algorithm for Lyapunov spectrum
initialize h, Q
evolve h until it is attractor (avoid initial transient)
evolve Q until it converges to the eigenvectors of the backward Oseledets matrix
% for / = do
0
< f(h)
D df
Q < D Q
(mod SONS) then
Q R. + qr Q)
24 log(R.)
end if
end for
X; = T
is guaranteed that under general conditions initially random orthonormal systems will exponentially converge towards a unique basis that is given by the eigenvectors of the Oseledets matrix 25 [60] A minimal example of this algorithm in pseudocode is shown in appendix(3] A feasible strategy to determine the reorthonermalization time interval tons is to get first rough estimate of the Lyapunov spectrum using short simulation time fimm and longer simulation
a small tons and repeat estimate of the Lyapunov spectrum. Another strategy is, to first iteratively adapt tons on a short simulation run to get an acceptable condition number, It should be noted that there exists a diversity of other methods to
and tans based on the Lyapunov spectrum the rough
estimate the Lyapunov spectrum 44 49
79 80)
27

---

K
Convergence of Lyapunov Exponents of RNNs
we demonstrate the convergence of the Lyapunov exponents. We show the estimate of the Lyapunov exponents X for i 20, 60, 80 for different initial conditions but identical network realization.
In Fig. 16
A/1/steps)
-10
-20
30
-40
10a
10
10a
steps
Figure 16: Convergence of Lyapunov exponents Convergence of selected Lyapunov exponents
X; for ten identical network realizations with different initial conditions with simulation time (i 1, 20, 60, 80) for a = and 1. (Other parameters: N = 80, tsim = 100 steps, lons 1).
28


Reviews:
[{"review_id": "WLsgLrbnl7", "reviewer": "Reviewer_tH3Y", "summary": "The presented paper proposes a new method for tackling numerical instabilities during training for recurrent neural networks.\nThe proposed method exploits a theoretical link between Lyapunov exponents and the singular values of the long-term Jacobian.\nThe set of experiments is well-chosen to showcase the extent to which gradient flossing helps stabilize training for tasks involving long-range dependencies.\nThere is an in-depth coverage of the practical aspects of the methods in the appendix, which is very useful for both machine learning practitioners and researchers.\n\nBy steering the Lyapunov exponents to 0, either before and/or during training, the authors are able to show consistent improvements on a variety of synthetic learning tasks involving long-range dependencies in temporal sequences.", "questions": "- (line 197 - 201) How does fig 2.C shows that the estimated condition number based on Lyapunov exponents can predict differences in condition number originating from finite network size $N$ ?\nThe authors refer to appendix F for this matter, but it doesn't seem related to this claim.\n- Can the author explain in more details what are the \"the finite network size fluctuations of the Lyapunov exponents\" in figure 1 ?\n- Figure 11 in the appendix shows that too much gradient flossing during training can be detrimental, which means that the local minima of the gradient flossing procedure hardly coincides with the one corresponding to the actual learning task. It is also said in the introduction that appropriate initialization schemes can ensure well-behaved gradient at initialization, a property that can be lost during training. I suspect that the gradient flossing procedure prevents the weights from falling into an \"ill-conditioned\" region. Thus, how often should gradient flossing be used should depend on, given a fixed learning rate, how many iterations are needed to go from a \"well-behave\" region to an ill-condition one. Can the author elaborate on this specific intuition ?\nI am asking this question because I suspect that you already found regimes where the regularization objective conflicts with the actual learning task.", "limitations": "As emphasized by the author, the Lyapunov exponent characterizes the singularity of the long-term Jacobian, and thus cannot characterize local-in-time behavior.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- There are experiments validating basic hypotheses, such as the fact that the proposed algorithm is able to control the Lyapunov exponents, and that estimating the Lyapunov exponents via the QR decomposition matches the theory.\n- Practical aspects such as when to perform gradient flossing, how often should it be used during training, and also how often should the QR decomposition should be computed when unrolling the neural network dynamics.", "weaknesses": "- There are no evaluation on real world complex benchmark to estimate the impact of the proposed contribution. However, the paper is already quite dense.\n\n- Minor remark: I don't understand the \"the finite network size fluctuations of the Lyapunov exponents\" in figure 1.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We deeply appreciate the insightful feedback provided by the reviewer. Your comments and concerns have not only broadened our perspective but also greatly aided in refining our manuscript to a higher standard. Specifically, we acknowledge the importance of evaluating gradient flossing on real-world complex benchmarks and are currently in the process of addressing this. Further, we have included additional details and clarifications based on your feedback, particularly concerning the \"finite network size fluctuations of the Lyapunov exponents\".\n\nHere are our more detailed responses to the specific questions:\n\n> (line 197 - 201) How does Fig 2.C shows that the estimated condition number based on Lyapunov exponents can predict differences in condition number originating from finite network size $N$? The authors refer to Appendix F for this matter, but it doesn't seem related to this claim.\n\nFigure 2C shows that different network realizations of finite size (N=80) lead to different condition numbers $\\kappa_2$ (for instance, different $\\kappa_2$ among green dots). We note that these dots lie close to the diagonal, indicating that the differences in $\\kappa_2$ coming from different network realizations are well-captured by the theoretical prediction on the y-axis.\n\n> Can the author explain in more detail what are \"the finite network size fluctuations of the Lyapunov exponents\" in Figure 1 ?\n\nAbsolutely. For finite network size $N$, different realizations of the random recurrent weight matrix $\\mathbf{W}$ will lead to slightly different Lyapunov exponents. Such differences are expected to vanish for large $N$. We added a new supplementary figure to demonstrate that. We note that parts of the variations among Lyapunov spectra might also come from the interplay of the different network realizations and the dynamics of gradient flossing. We added this note to the updated manuscript.\n\n> Figure 11 in the appendix shows that too much gradient flossing during training can be detrimental, which means that the local minima of the gradient flossing procedure hardly coincides with the one corresponding to the actual learning task. It is also said in the introduction that appropriate initialization schemes can ensure a well-behaved gradient at initialization, a property that can be lost during training. I suspect that the gradient flossing procedure prevents the weights from falling into an \"ill-conditioned\" region. Thus, how often should gradient flossing be used should depend on, given a fixed learning rate, how many iterations are needed to go from a \"well-behaved\" region to an ill-condition one. Can the author elaborate on this specific intuition? I am asking this question because I suspect that you already found regimes where the regularization objective conflicts with the actual learning task.\n\nA thoughtful observation! Continued gradient flossing can indeed interfere with the actual learning task, as evident in Supplementary figure 11. We agree that the optimal gradient flossing schedule should depend on how many iterations are necessary to go from ill-conditioned to a well-behaved region in parameter space. After submission, we found that instead of having a fixed number of gradient flossing steps, \"adaptive gradient flossing\", that automatically stops once the Lyapunov exponents are close to zero is even better. It is conceivable that refining this technique might depend on the specific task, and further improvements are being explored.\n\nTo conclude, we sincerely thank you for your in-depth review, and assure you of our commitment to enhancing the quality of our work based on the valuable feedback received. We believe that the integration of gradient flossing will mark a significant advancement in the domain, and are keen to contribute effectively to the community's understanding.\n"}, {"review_id": "08MXE32orC", "reviewer": "Reviewer_ZKKk", "summary": "Authors propose gradient flossing for RNNs, which adds an additional regularization term that keeps the sum of Lyapunov exponents close to zero. This essentially encourages the singular values of the long-term Jacobian to be close to 1, hence addressing the vanishing/exploding gradient problem in RNNs. The authors point out that minimizing this regularization term again has the same gradient instability problem, and they purpose to apply QR decomposition once in a while to avoid the ill-conditioned long-term Jacobian pass. Finally, the authors show gradient flossing achieves faster convergence and a higher success rate on simple long-term tasks up to 500 steps. ", "questions": "\n* I’m still unclear how QR decomposition is able to avoid the long-term Jacobian being ill-conditioned. \n\n* Why can we not use QR decomposition directly on the original RNN objective in equation 1 to avoid the long-term gradient instability problem?\n", "limitations": "See above.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "* The idea of regularizing the Lyapunov exponent and using QR decomposition is novel and interesting.\n* The problem is well-motivated.\n", "weaknesses": "* Empirical evaluations are toy. Did not include results on standard long-time horizon benchmarks like Longe Range Arena [1].\n\n* The paper could have studied gradient flossing with architectures that are more suited for long horizon tasks such as state space models [2] and LRUs [3], instead of Vanilla RNNs and LSTMs.\n\n* Authors should compare to more recent methods as baselines such as HiPPO Matrices[2] and initialization methods used in LRU [3]. Those methods tend to be more competitive than a simple orthogonal initialization.\n\n* Clarity of the paper can be improved. For example, $\\tau$ in equation 1 is never defined, and equation 4 could be better explained. In addition, I did not quite understand the technical parts of the paper explaining how equation 4 is calculated in practice (see questions below). \n\n[1] Tay, Yi, et al. \"Long range arena: A benchmark for efficient transformers.\" arXiv preprint arXiv:2011.04006 (2020).\n\n[2] Gu, Albert, Karan Goel, and Christopher Ré. \"Efficiently modelling long sequences with structured state spaces.\" arXiv preprint arXiv:2111.00396 (2021).\n\n[3] Orvieto, Antonio, et al. \"Resurrecting recurrent neural networks for long sequences.\" arXiv preprint arXiv:2303.06349 (2023).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We deeply appreciate reviewer ZKKk for their thorough and insightful review of our manuscript. Your constructive feedback is invaluable, and we are committed to addressing each concern raised to enhance the quality and clarity of our work.\n\n> I’m still unclear how QR decomposition is able to avoid the long-term Jacobian being ill-conditioned.\n\nTo provide a more intuitive explanation before diving into technicalities: Imagine the Jacobians as a series of transformations of small perturbations that can stretch or compress the space along the trajectory of the RNN. Over time, when many such transformations are multiplied, they can cause vectors to become too aligned, leading to numerical instability. The QR decomposition offers a mathematically exact way to correct for this by ensuring that the space remains 'evenly spread' after each transformation.\n\nNow the technicalities: The core idea is that each individual Jacobian $D_s$ is not ill-conditioned, but their product $T_t(h_\\tau)$ is. Using QR decomposition, one can iteratively decompose the product of Jacobians sufficiently often such that the product is not ill-conditioned [1,2]. \n\nTo illustrate this geometrically, consider the evolution of initially random vectors $v^i$​ and $v^j$​ under the influence of the product of Jacobians $T_t(h_\\tau)$. As $t-\\tau$ is increased, the angle between $T_t(h_\\tau) v^i$ and $T_t(h_\\tau) v^j$​ become too small for numerical computation. This problem can be circumvented through the iterative application of QR decomposition, which repositions the vectors into orthogonal subspaces, thereby maintaining numerical stability [3].\n\nTo further illustrate this point dynamically for the forward pass of the calculation of Lyapunov exponents: Almost all perturbations of the initial recurrent network state $h$ will align over time with the fastest diverging (or slowest converging) tangent space direction, which is the first covariant Lyapunov vector [2]. The QR decomposition of a volume element in the tangent space will express the volume in terms of an orthonormal matrix $Q$ and an upper triangular matrix $R$, where the $Q$ matrix will converge to a unique basis when constraining the diagonal elements of $R$ to be positive [1,4]. \n\nIn summary, QR decomposition helps by iteratively repositioning vectors into orthogonal subspaces, thereby ensuring that numerical calculations remain accurate. The mechanics of this process align with the behavior of the Lyapunov exponents, as detailed in our earlier explanation.\n\n\n> Why can we not use QR decomposition directly on the original RNN objective in equation 1 to avoid the long-term gradient instability problem?\n\nGreat suggestion! In principle, one could implement the backward pass of backpropagation through time by hand and iteratively apply QR decomposition on the long-term Jacobians, directly regularizing the diagonal of the $R$ matrix. However, for gradient flossing, we chose here the much simpler solution of calculating the logarithm of the singular values of the long-term Jacobian in the forward pass numerically robustly and regularizing them using automatic differentiation. This has the advantage, that automatic differentiation takes care of the backward pass, which makes it easy to use gradient flossing, as only the forward pass has to be implemented. We will mention the elegant (and possibly more efficient) idea by the referee in the outlook of the manuscript.\n\nWe are genuinely grateful for the feedback, and we believe that addressing these points will significantly elevate the quality and impact of our manuscript. We hope that our revisions will meet your expectations and further advance the field's understanding of gradient flossing and its applications in RNNs. We are looking forward to answer follow-up questions and answer further comments.\n\n\n\n[1] K. Geist, U. Parlitz, and W. Lauterborn, Comparison of Different Methods for Computing Lyapunov Exponents, Prog. Theor. Phys. 83, 875 (1990).\n[2] A. Pikovsky and A. Politi, Lyapunov Exponents: A Tool to Explore Complex Dynamics (Cambridge University Press, Cambridge, 2016).\n[3] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Lyapunov Characteristic Exponents for Smooth Dynamical Systems and for Hamiltonian Systems - A Method for Computing All of Them. I - Theory. II - Numerical Application, Meccanica 15, 9 (1980).\n[4] S. V. Ershov and A. B. Potapov, On the Concept of Stationary Lyapunov Basis, Physica D: Nonlinear Phenomena 118, 167 (1998).\n"}, {"review_id": "1PZ4ycjYdI", "reviewer": "Reviewer_DGA6", "summary": "The authors propose a new technique for handling gradient instability (vanishing/exploding) in neural network models, especially sequence models such as RNNs. Specifically, by leveraging results from recent works that establish a connection between Lyapunov growth exponents (from dynamical systems) and the singular values of the model-Jacobian itself, the authors propose a (regularization) approach that they coin ‘gradient flossing’ to control the evolution of these singular values desirably. They demonstrate that, as a result, the gradient signal remains stable, thereby allowing for successful training of the model itself.\n\nAdditionally, the authors provide insights on the steps (schedule, pre/during training) to perform such a regularization leading to optimal gradient stabilization & computational benefits.", "questions": "- Eq. (2) does not necessarily need to refer to the ‘long-term Jacobian’ unless $\\Delta{t} \\equiv t - \\tau \\to \\infty$. That should be clarified.\n- Is the loss function from Eq. (4) used as an additional or the only loss to train?\n- Eq. (4) does not seem to contain a reference to the target $\\lambda$ values that the authors eventually aim to constrain them to\n- While the analytical $t\\to\\infty$ limit is understandable, can you provide an estimate of a reasonable sequence length for *practical data* that can serve as a finite size equivalent of this limit?\n- Is Eq. (5) subject to the assumption that parameter updates are gradient-based?\n- I believe the ’s’-index in Eq. (6) is supposed to run from $\\tau$ to be exact (?)\n- Clearly, in arriving at Eq. (6), the $Q$’s from *subsequent* time-steps $s$ and $s+1$ act on each other to reduce to unity ($Q_{s+1}^T Q_{s}$), leaving with the $R$’s only. Can the authors explain why that might be the case in general? Is there an assumption about a small learning rate here?\n- Since the authors claim the QR decomposition need not be applied at every time step, will the equivalent Eq. (6) thus arrived at still remain to be factually correct? The expression of $T_t$ contains an iterative product of $ D_s$, of which if one chooses to QR decompose a handful only, then that does not seem to yield $\\prod_s R^s$ for some period $t_{ons}$.\n- In Lines 140-141, you refer to using vanilla gradient descent for the experiments of Fig. 1 when you specify in Lines 132-133 that experiments were performed using ADAM.\n- Isn’t a target $\\lambda$ value of -1 indicative of vanishing gradients as per Eq. (3)? How is the model still reliably attaining that value regardless, then?\n- In Fig. 2 (A), the ‘theory’ and ‘simulation’ legends seem to be interchanged.\n- In Fig. 2 (A) & (C), the numerical and theoretical results should overlap in the long time-step regime, yet as time increases, the two appear to diverge in the former and maintain a uniform scaling (~10^2) in the latter. Can you explain why?\n- Can you explain what you mean by ‘flossing before training’? That does not seem to be well clarified in the text.", "limitations": "To the best of my knowledge, the authors act commendably in attributing to the various limitations of their approach.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The authors visit a longstanding problem in training NN models, esp. sequence models with long temporal dependencies. By leveraging recently established connections between Lyapunov exponents and the singular values of the (long-term) Jacobian, the authors propose a new technique in the form of:\n- Defining a loss function out of these Lyapunov exponents\n- Determine a clever strategy to evaluate this loss function through a QR decomposition\nto constrain these exponents directly to stable target values. This, in turn, allows control over the evolution of the Jacobian spectra itself.\n\nThe authors provide experimental evidence supporting this technique on two specific tasks that, by design, require long-time gradient signals to propagate efficiently.\n\nFurthermore, the authors provide practical insights into the (expensive) computational needs of the exact implementation of such a scheme and the place (pre/during training) of the same. They then propose a modified, intermittent regularization routine with preserved benefits yet having lower compute necessities.\n\nThe authors do an overall good job of satisfactorily articulating the problem statement, the proposed solution, observed benefits, and implementation limitations.", "weaknesses": "While I believe this is a good contribution, I feel the work can undergo a fair amount of extra polishing in terms of clarity. Please find below (Questions) my comments/questions on areas where this work can be further improved.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for their constructive feedback which has greatly contributed to improving the quality of our manuscript.\n\n> Eq. (2) does not necessarily need to refer to the ‘long-term Jacobian’ unless $\\Delta t \\equiv t-\\tau\\rightarrow \\infty$. That should be clarified.\n\nWe have updated the manuscript to only refer to $T$ as the long-term Jacobian when $t-\\tau$ is long, thereby clarifying this point.\n\n> Is the loss function from Eq. (4) used as an additional or the only loss to train?\n\nWe've clarified that while our current implementation utilizes the gradient flossing loss from Eq. (4) as the sole loss during its epochs, it can also be combined with a task-specific loss for computational efficiency.\n\n> Eq. (4) does not seem to contain a reference to the target $\\lambda$ values that the authors eventually aim to constrain them to.\n\nCorrect, we here set the target Lyapunov exponents to 0 and thus directly minimize the square of the Lyapunov exponents. It would be an exciting future direction to explore if other target values for the Lyapunov exponents can further improve trainability. For instance, in some tasks, one might want to have a certain number of contracting or growing directions.\n\n> While the analytical limit $t\\to\\infty$ is understandable, can you provide an estimate of a reasonable sequence length for practical data that can serve as a finite size equivalent of this limit?\n\nThank you for pointing this out. We added an extra figure to demonstrate how quickly the Lyapunov spectrum convergences with the number of time steps. We found that a relative error of 1% is already achieved after less than 100 steps. We note that even shorter sequence lengths improve training. We suspect that the sequence length should be on the order of the longest time horizon of that task that has to be bridged by backpropagation through time. \n\n\n> Is Eq. (5) subject to the assumption that parameter updates are gradient-based?\n\nEq. 5 is a general expression, and there's no assumption regarding the parameter update mechanism therein.\n\n> I believe the ’s’-index in Eq. (6) is supposed to run from $\\tau$ to be exact (?)\n\nWe are not exactly sure what this comment refers to. In the expression for the Lyapunov exponent,\n\n$\\lambda_{i}=\\lim_{t\\to\\infty}\\frac{1}{t}\\log\\prod_{s=1}^{t}R_{ii}^{s}=\\lim_{t\\to\\infty}\\frac{1}{t}\\sum_{s=1}^{t}\\log R_{ii}^{s}.$\nthe index $s$ correctly runs from s=1 to t, but one could equivalently write\n\n$\\lambda_{i}=\\lim_{(t-\\tau)\\to\\infty}\\frac{1}{t}\\log\\prod_{s=\\tau}^{t}R_{ii}^{s}=\\lim_{(t-\\tau)\\to\\infty}\\frac{1}{t-\\tau}\\sum_{s=\\tau}^{t}\\log R_{ii}^{s}.$\n\n> Clearly, in arriving at Eq. (6), the $Q$’s from subsequent time-steps $s$ and $s+1$ act on each other to reduce to unity ($Q_{s+t}^T Q_{s}$), leaving with the $R$’s only. Can the authors explain why that might be the case in general? Is there an assumption about a small learning rate here?\n\nNo, there is no assumption of small learning rates here. We apologize in case this was confusing, but $s$ is the time index, and learning happens across epochs. \n\n\n> Since the authors claim the QR decomposition need not be applied at every time step, will the equivalent Eq. (6) thus arrived at still remain to be factually correct? The expression of $T_t$ contains an iterative product of $D_s$, of which if one chooses to QR decompose a handful only, then that does not seem to yield $\\prod_s R^s$ for some period .\n\nYes, Eq. (6) remains mathematically correct even for infrequent QR decomposition. If one chooses to QR decompose less often, the resulting (fewer) R will big bigger:\nMore explicitly the result of \n$D_2 D_1 = D_2 Q_1 R_1 = Q_2 R_2 R_1$\nand \n$D_2 D_1 = \\tilde Q_2 \\tilde R_2$ is identical (as QR-decomposition is unique when constraining the signs of the diagonal of $R$ to be positive), thus $\\tilde R_2  = R_2 R_1$.\n\n> In Lines 140-141, you refer to using vanilla gradient descent for the experiments of Fig. 1 when you specify in Lines 132-133 that experiments were performed using ADAM.\n\nOur apologies for the inconsistency. We've revised the manuscript to clarify that ADAM was the optimizer used in all our experiments.\n\n> Isn’t a target value of -1 indicative of vanishing gradients as per Eq. (3)? How is the model still reliably attaining that value regardless, then?\n\nYou're right. A Lyapunov exponent of -1 does suggest vanishing gradients. However, we utilized the QR-algorithm to achieve this target, which inherently avoids the issue of vanishing gradients.\n\n> In Fig. 2 (A), the ‘theory’ and ‘simulation’ legends seem to be interchanged.\n\nThis oversight has been addressed in the revised manuscript.\n\n> In Fig. 2 (A) & (C), the numerical and theoretical results should overlap in the long time-step regime, yet as time increases, the two appear to diverge in the former and maintain a uniform scaling (~10^2) in the latter. Can you explain why?\n\nThanks for pointing this out. We've thoroughly investigated the differences and found the cause for the mismatches. For Figure 2A, discarding initial transients and increasing numerical precision resolves the discrepancy. In Figure 2C, there's no actual mismatch; the x-ticks and y-ticks were differently set, and we have updated the figure to rectify this.\n\n> Can you explain what you mean by ‘flossing before training’? That does not seem to be well clarified in the text.\n\nWe meant that our approach first involves several gradient flossing steps, after which the main task training begins without additional gradient flossing. This point has been elaborated upon for clarity in the revised manuscript.\n\nThank you once again for your valuable feedback. We believe these revisions have significantly strengthened our submission."}, {"review_id": "6ezueY32NR", "reviewer": "Reviewer_5Z2M", "summary": "The authors propose a novel approach called \"gradient flossing\" to tackle the instability of gradients in the training of recurrent neural networks (RNNs). The process of pushing the Lyapunov exponents toward zero is referred to as \"flossing\" the gradients. This stabilizes the gradients and improves network training.", "questions": "Can the authors comment that such a technique of pushing the Lyapunov exponents to zero could limit the capacity to learn certain tasks? Maybe enabling the network to amplify or forget things is important for some task contexts, which has been a motivation for specialized architectures like the LSTM. ", "limitations": "Experimental limitation mentioned above.", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "+ The authors show improvements in several toy sequence learning tasks: Delayed Copy Task, Temporal XOR Task, etc.\n\n+ Idea of using Lyapunov exponents towards zero to mitigate vanishing and exploding gradients is well justified and works as expected on toy tasks.\n\n", "weaknesses": "The lack of any experiments on non-simulated data is problematic, and the explanation given by the authors is not convincing. The paper's authors explicitly state that they deliberately did not use sequential MNIST or similar non-toy tasks commonly used to probe exploding/vanishing gradients. \n\n> We deliberately do not use sequential MNIST or similar toy tasks commonly used to probe exploding/vanishing gradients, because we want a task where the structure of long-range dependencies in the data is transparent and can be varied ad libitum\n\nWhile toy tasks are more transparent, real-world data is important to understand the proposed methods' potential failure modes and expected characteristics in *uncontrolled* settings. Avoiding running experiments on the basis that the results are more difficult to interpret than simulated data is not an acceptable excuse. Transparency in the work itself is hindered when experiments are consciously avoided.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\n\nThank you for your thoughtful review and feedback on our manuscript. We appreciate the acknowledgment and recognition of our novel approach involving the use of Lyapunov exponents towards zero to address the challenges of vanishing and exploding gradients in RNNs.\n\n> While toy tasks are more transparent, real-world data is important to understand the proposed methods' potential failure modes and expected characteristics in uncontrolled settings. Avoiding running experiments on the basis that the results are more difficult to interpret than simulated data is not an acceptable excuse. Transparency in the work itself is hindered when experiments are consciously avoided.\n\nWe recognize the importance of assessing our method on real-world datasets and thank the reviewer for emphasizing this aspect. Our rationale behind beginning with synthetic toy tasks was to \"open the black box\" and comprehensively understand the tangent space structure of gradient descent from a dynamical systems perspective and how it is improved by gradient flossing. This rigorous examination has allowed us to delve into not just vanishing and exploding gradients, but also:\n\n- Observe how the complexity of tasks impacts the number of required Lyapunov exponents to be flossed, as showcased in Supplementary figure 2.\n  \n- Study how gradient flossing enhances the norm of the long-term Jacobian (Supplementary figure 7).\n\n- Establish a relationship between the training epoch, where accuracy exceeds the chance level, and the increase in the gradient norm (Supplementary figure 8).\n\n- Analyze the positive effects of gradient flossing on the condition number of the recurrent weight gradient (Figure 2 and supplementary figure 8).\n\n- Formulate a mathematical correlation between the post-gradient-flossed Lyapunov exponents and task complexity, quantified by the delay \\(d\\) (Supplementary figure 9).\n\nTo add, our main contribution hinges on the connection between the condition number of the long-term Jacobian and the Lyapunov spectrum. Our work, especially in the section 5 on \"the condition number of the Long-Term Jacobian\" (Figure 2), provides novel insights into the dimensionality of the gradient update, offering an analytical perspective on the nature of signal propagation backward in time during RNN training.\n\n> Can the authors comment that such a technique of pushing the Lyapunov exponents to zero could limit the capacity to learn certain tasks? Maybe enabling the network to amplify or forget things is important for some task contexts, which has been a motivation for specialized architectures like the LSTM.\n\nYou've raised an excellent point. Our analysis in Supplementary figure 11 demonstrates potential conflicts of gradient flossing with certain task goals, where persistent gradient flossing over multiple Lyapunov exponents may inhibit the RNN's nonlinearity. We've addressed this challenge by limiting gradient flossing steps and subsequently focusing on task loss minimization without gradient flossing. For architectures like LSTM, this scenario could differ, as their latent variables can effectively carry gradient information over extended time frames, even amidst active nonlinearities on \\(h\\).\n\nIn conclusion, we deeply value the feedback and are keen to explore the applicability and nuances of our gradient flossing technique on real-world datasets beyond toy tasks, to make our contribution more comprehensive and robust in the field of RNN training.\n"}]

Author rebuttal:
We thank the reviewers for their time and valuable feedback on our manuscript, "Gradient Flossing: Improving Gradient Descent through Dynamic Control of Jacobians". We have carefully addressed each of the reviewers' comments in the subsequent sections. Here, we provide a concise summary to address the core critique raised by most reviewers.

We have taken note of the primary concerns raised, notably the application of our method to synthetic tasks instead of real-world tasks and the lack of comparison with other advanced architectures. We would like to highlight that the main objective of our manuscript was to bridge the understanding between the condition number of the long-term Jacobian and the Lyapunov spectrum, and "opening the black box" of gradient flossing in controlled synthetic tasks allowed us to rigorously investigate the tangent space structure. We utilized simple toy tasks because they allow a granular understanding by controlling every aspect of the problem, a strategy that has previously been employed effectively in foundational neural network research.

To address specific concerns raised and to provide a clearer picture of our results and methodology, we have included an additional page of figures:
*  **Updated Figure 2** provides a clearer representation of how gradient flossing reduces the condition number of the long-term Jacobian and improves numerical accuracy. We have also corrected an error in the legend.
* A new plot showcases the **convergence of the Lyapunov exponents** over time steps. With 20 distinct network realizations, it is evident that convergence is reached swiftly, within 100 time steps.
* We demonstrate how gradient flossing enhances the **dimensionality of the error gradient** of the recurrent weights, as calculated based on the SVD of these error gradients.

That being said, we acknowledge the significance of applying our findings to real-world tasks. Building on the foundational understanding we've achieved through synthetic tasks, we are poised to extend gradient flossing to more realistic tasks and diverse architectures in our subsequent research.

Now, let us address the specific points raised by each reviewer, including additional figures that further support our methodology and findings.

Meta-review: A practical method to deal with instabilities of gradients when training recurrent network is developed. The Lyapunov spectrum is estimated while training and the system dynamics is regularized through. It also provides theoretical insights into how to tune unstable networks at rest and increase the time horizon of “useful” recurrent dynamics. Despite mostly demonstrating on toy problems, this paper is of broad interest to the NeurIPS community.
Decision: Accept (poster)