Analyzing and Improving Surrogate Gradient Training in Binary Neural
Networks Using Dynamical Systems Theory

Rainer Engelken 1 LF Abbott 1

2. Surrogate Lyapunov Exponents

Abstract
Training Binary Recurrent Networks on tasks thatspan long time horizons is challenging, as thediscrete activation function renders the error landscape non-differentiable. Surrogate gradient training replaces the discrete activation function witha differentiable one in the backward pass but still
suffers from exploding and vanishing gradients.

g = 0.5
g=1
g=2
g = 10

sgn

A

1
0
0

1

B

1
3

S
i (1/s)

Output

0

0

4
1
1

C

We leverage the connection between gradient stability and Lyapunov exponents to address thisissue from a dynamical systems perspective, extending our previous work on Lyapunov exponent
regularization to non-differentiable systems. We
use differentiable linear algebra to regularize surrogate Lyapunov exponents, a method we callsurrogate gradient flossing.

Tt (hτ ) =

0

x

Qt−1

1

∂hτ ′ +1
τ ′ =τ ∂hτ ′

2

i

1

D

70

log (|Dij|

g=10

2

1
neuron i

2

10

20
1

10

2

0

2

10

2

20

20

neuron j

Figure 1. Surrogate Lyapunov exponents quantify vanishing andexploding gradients of surrogate gradient training of binaryneural networks A) In binary networks hs+1 = W sgn(hs ) +Vxs+1 , the step-like activation function leads to discontinuous,
non-differentiable dynamics, creating a rugged error landscapein gradient-based learning. Surrogate gradient training replaces
the step-like activation function with a continuous function φ in
the backward pass, implemented here as a normalized derivativefor φ(x) = tanh(x) with φ′g (x) = sech(gx)2 . The subscript
denotes the sharpness g of the surrogate gradients. B) Surrogate
Lyapunov exponents λS
i for different values of the surrogate sharpness g. C) Vanishing and exploding gradients in discrete recurrent
networks arise from the attenuation or amplification of the prodQ∂hτ ′ +1
uct of surrogate Jacobians t−1
τ ′ =τ ∂h ′ D) Surrogate Jacobian

We show that surrogate gradient flossing enhances performance on temporally demanding
tasks.

1. Introduction
Quantizing neural networks reduces computational cost andmemory usage but complicates gradient-based training due
to non-differentiable discrete activation functions. Traditional approaches, such as surrogate gradients [1], mitigate
this issue by replacing the discrete nonlinearity with a differentiable expression in the backward pass. This technique
is widely used for training spiking neural networks [2].

τ

Dsurrogate = W ⊙ φ′ (hs ) for g = 10, which results in a columnsparse Jacobian matrix. Note the logarithmic scale.

This paper extends our previous work on Lyapunov exponent regularization [3] to non-differentiable models. We
introduce surrogate Lyapunov exponents to quantify surrogate gradient stability, provide analytical expressions forthese exponents, and analyze the impact of the sharpness ofthe surrogate gradient.

We introduce surrogate Lyapunov exponents as a tool for
analyzing and mitigating vanishing and exploding gradientsin discrete recurrent networks, described by
hs+1 = W sgn(hs ) + Vxs+1 .

1(1)

The core idea is to calculate surrogate Lyapunov exponents,
which quantify exploding and vanishing gradients, basedon the derivative of the surrogate activation function usedin the backward pass. This differs from the conventional
Lyapunov exponents that measure how small perturbationsin the forward pass grow or shrink on average [4].

Center for Theoretical Neuroscience, Zuckerman Institute,
Columbia University, New York, NY, USA. Correspondence to:Rainer Engelken <re2365@columbia.edu>.
Published at the 2 nd Differentiable Almost Everything Workshop
at the 41 st International Conference on Machine Learning, Vienna, Austria. July 2024. Copyright 2024 by the author(s).

1Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

In surrogate gradient training, network parameters θ aregradually changed to reduce a loss Lt by smoothing the nondifferentiable state transitions with a differentiable function
for the backward pass (see Fig 1A). For the discretized dynamics of a binary or spiking recurrent neural network, thesurrogate gradient of the loss Lt with respect to parameters θis computed by unrolling the network dynamics backwards
in time and propagating error signals across time steps (SeeFig 1C):τ =t−1
∂L X
∂Lt
=
∂θ
∂ht τ =0

t−1
Y

Dsurrogateτ′τ ′ =τ

!

∂hτ
,
∂θ

average, neither strongly attenuated nor amplified, facilitating reliable information propagation across many time
steps. We note that the surrogate Lyapunov exponents of
discretized spiking network dynamics are different from thetrue Lyapunov exponents of the forward dynamics [12]–[18].
Figure 1B shows the surrogate Lyapunov exponents for different values of surrogate sharpness g in a binary RNN. We
find that with increasing g, most surrogate Lyapunov exponents become very negative, leading to vanishing gradientsand an ill-conditioned product of surrogate Jacobians during
training (See Appendix G and F for details). Conversely,
for small values of g, we find, depending on W, and V,
many surrogate Lyapunov exponents close to zero, whichmore easily give rise to oscillations or exploding gradientsduring task training once the spectral radius of the recurrent
weights grows slightly. See supplementary Figures 7 and 8
in Appendix L for an analysis of the convergence of surrogate Lyapunov exponents with simulation time for different
values of g.

(2)

s+1
where Dsurrogate= ∂hs∂hs is the surrogate Jacobian that describes the propagation of infinitesimal state changes across
one time step using a differentiable approximation of theactivation function sgn. See Fig 1D for an example surrogate Jacobian of a binary recurrent network with surrogate
sharpness g = 10 and Appendix E for the derivation of
the analytical surrogate Jacobian
for a LIF network. The
Qt−1
product of surrogate Jacobians τ ′ =τ Dsurrogate
in Eq. 2 reτ′
flects the recursive dependence of network states on theirpast values and tends to give rise to vanishing and exploding gradients, making training across long time horizons
challenging in both discrete and continuous recurrent neuralnetworks.
Recent studies have shown that the singular values of the
product of Jacobians, which measure the rate of gradientexplosion, are mathematically closely related to Lyapunov
exponents [4]–[11]. Lyapunov exponents are a core concept in dynamical systems theory, quantifying stability andchaos. The largest Lyapunov exponent λ1 measures theaverage rate of divergence or convergence of nearby initialconditions. An N -dimensional system has N Lyapunov
exponents that quantify the average growth rates of volume
elements in the tangent space. To define surrogate Lyapunov
exponents, we adapt this concept for surrogate gradients byinstead considering the evolution of perturbations near thereference trajectory when evaluating the surrogate activation
function (see Fig. 1C). We thus linearize along a reference
trajectory of the forward dynamics of the discrete neuralnetwork but evaluate the surrogate Jacobians:1λsurrogate= lim log
i
t→∞
t

"

σi

t−1
Y

τ ′ =0Dsurrogate
τ′

!#

,

3. Surrogate Gradient Flossing
Our recent work built upon the mathematical relationshipbetween Lyapunov exponents and the challenge of exploding
and vanishing gradients in training of recurrent networks,
and introduced a novel regularization technique termed gradient flossing [3]. This method stabilizes the training process by constraining Lyapunov exponents to values close tozero, ensuring more trainable network dynamics througha better-conditioned product of Jacobians. Here, we apply
this idea for binary neural networks trained with surrogate
gradients. In our approach, we regularize the surrogate Lyapunov exponents, extending the concept of gradient flossing
from conventional Lyapunov exponents to the surrogate domain [3]. To mitigate exploding and vanishing gradients,
we regularize the mean of the squares of the first k largest
surrogate Lyapunov exponents λ1 , λ2 , . . . , λk , resulting inPkan additional loss term Lsurrogate flossing = k1 i=1 λ2i . To
calculate the surrogate Lyapunov exponents, we employ anestablished iterative orthonormalization method using QR
decomposition [4], [8], [19], rather than directly evaluating
the product of surrogate Jacobians (Eq.2). This method isnumerically robust, avoiding issues with the ill-conditioned
long-term surrogate Jacobian (See AppendixA for details
and pseudocode). Surrogate Lyapunov exponents are obtained by time-averaging the logarithms of the diagonal
during QR decompoentries of the Rs matrices
Pobtainedtsitions: λi = limt→∞ 1t s=1 log Rsii . To backpropagate
the gradient of the squared surrogate Lyapunov exponents,
we leverage recent progress in differentiable linear algebra,
which provided an analytical expression for the pullback of
the QR decomposition (See Appendix B). For clarity, we
described surrogate gradient flossing in terms of stochastic
gradient descent, but we actually implemented it with the(3)where σi denotes the ith singular value of the surrogate Jacobian product. These exponents reflect the average asymptotic growth rates of infinitesimal perturbations in the surrogate tangent space of the forward dynamics, thus constraining signal propagation during BPTT over long time horizons.
Specifically, modes with surrogate Lyapunov exponents nearzero indicate directions along which error signals are, on2Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

ADAM optimizer [20] using standard hyperparameters η,
β1 , and β2 . An example implementation in Julia [21] using
Flux [22] is available here. This algorithm is versatile andcan be applied to various neuron models and network architectures. The surrogate Jacobian matrix Dsurrogate can becalculated analytically, as done here, or through automaticdifferentiation.

manipulate one, several, or all surrogate Lyapunov exponents before or during network training, thereby shaping
both the norm and the condition number of the surrogate
gradients over long time horizons.

5. Surrogate Gradient Flossing EnhancesTrainability

0.2

g:

S
i (1/s)1

S (1/s)

0.00
0.0

0.25
0.50

k =

0.5
0.4

50

1

k = 160.75

1.0

C 10 1

k = 32

1.00
50

100
Epochs

150

0

10

20

i

30

40

50

Epochs3001|ddh0 |

0

0

Figure 2. Surrogate gradient flossing regularizes surrogate Lyapunov exponents and facilitates gradient signal propagation inbinary neural networks A) The first surrogate Lyapunov exponent of a recurrent binary network plotted as a function of training
epochs for different surrogate sharpness g. The square of the firstsurrogate Lyapunov exponent is minimized using gradient descent.
B) Surrogate Lyapunov spectrum of a recurrent binary network
after different numbers of Lyapunov exponents k have been driven
towards zero via surrogate gradient flossing for k ∈ {1, 16, 32}.
The gray lines show the surrogate Lyapunov spectra before surrogate gradient flossing. Parameters: network size N = 80, g = 1
for B. Input as in Fig. 3. The thin semitransparent lines in A and B
indicate nine network realizations; the full lines are their average.

10 13

100

50
D

10

15
20
complexity (delay)

25

10

15
20
complexity (delay)

25

100

|ddh0 |

B

0.2

B

adaptive flossing
no flossing

accuracy(%)

A

100

accuracy(%)

A

4. Control of Surrogate Lyapunov Exponents

0Epochs

3001

10 13

Figure 3. Gradient flossing improves binary RNN trainingA) Test accuracy for binary RNNs trained on the delayed temporalbinary XOR task yt = xt−d/2 ⊕ xt−d with adaptive gradient
flossing during training (orange) and without gradient flossing(blue) for d = 18. Solid lines are the median across 9 network
realizations, and individual network realizations are shown in transparent fine lines. B) Mean final test accuracy as a function of taskdifficulty (delay d) for delayed XOR task. C) Gradient norm withrespect to initial network state h0 . D) Gradient norm with respect
to initial network state as a function of temporal task complexity
T averaged over training epochs.

Figure 2 demonstrates that surrogate gradient flossing can
set one or several surrogate Lyapunov exponents to a target
value via gradient descent with the ADAM optimizer in
networks of binary neurons for different values of surrogate
gradient sharpness g. The initial entries of the recurrent
weight matrix W are independently drawn from a Gaussian
distribution with zero mean and variance 1/N . The recurrent binary network receives a random binary input stream,
and the input weights V are drawn from N (0, 1/Nin ). W,
g, and V are all modified by surrogate gradient flossing.
Figure 2A shows that for randomly initialized binary RNNs,
the surrogate Lyapunov exponent can be modified by surrogate gradient flossing to approach a desired target value.
During surrogate gradient flossing, the exponents quickly
approached the target value of zero for various surrogate
gradient sharpness values g. Notably, networks with large
g converge more slowly towards the target value of zero or
may not converge at all.
Figure 2B shows surrogate gradient flossing for different
numbers of flossed surrogate Lyapunov exponents k. During gradient descent, the mean of the squares of 1, 16, or
32 surrogate Lyapunov exponents is used as the loss. The
resulting surrogate Lyapunov spectrum after flossing shows
the corresponding number of exponents driven close to zero.
We conclude that surrogate gradient flossing can selectively

We present numerical results on a task with variable temporal complexity, demonstrating that adaptive surrogate
gradient flossing enhances the trainability of binary RNNs.
In adaptive surrogate gradient flossing, we initialize theinput weights V, recurrent weights W, and output weights
randomly, set g = 1, and first employ surrogate gradientPflossing for 500 epochs by minimizing Lflossing =1
i = 1k λ2i . We then alternated between surrogate grak
dient flossing and task training every 100 epochs until task
training epoch 500, and then continued with pure task training. Afterwards, we performed one episode of gradient
flossing adaptively whenever either λS1 > 0.2 or the gradiTent norm | ∂L
∂W | exceeded 1.
We consider a temporal XOR task that requires the binary
RNN to perform a nonlinear input-output computation on asequential stream of random binary input with a delay of dsteps. Specifically, the target output y is yt = xt−d/2 ⊕xt−d ,where x is a random Bernoulli process x ∈ {0, 1}. Thistask requires retaining d items in memory and performing anonlinear operation on them. While solving the task is trivial and the correct solution could be implemented manually,
RNNs trained on this task quickly encounter exploding or
vanishing gradients. We train the network by minimizing
3Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

the cross-entropy loss between the one-hot encoded target
output and a linear readout. To successfully solve this task,
the recurrent binary network must bridge the time horizonof d steps. We change the temporal complexity of the task
by scaling the delay d. Figures 3A and B show that adaptive
surrogate gradient flossing helps train networks for substantially longer delays d, while networks without flossing
fail on this temporally challenging task. Figures 3C and Dshow that binary networks suffer from vanishing gradients,
as revealed by measuring the gradient with respect to theinitial network state h0 . This issue is mitigated by adaptive
surrogate gradient flossing.

complex tasks. We adopted a dynamical systems perspective to analyze and improve the surrogate gradient trainingof binary and spiking neural networks.
We introduced surrogate Lyapunov exponents, a novel set ofindicators based on dynamical systems, to quantify the surrogate gradient flow in the backward pass, which is crucial
for addressing exploding and vanishing gradients. Surrogate
Lyapunov exponents account for the entire trajectory, revealing complex neuronal interactions previously neglected in
the analysis of surrogate gradients (see Appendix H for more
literature review). We called them surrogate Lyapunov exponents because, unlike conventional Lyapunov exponents thatstem from a linearization of the forward dynamics, thesearise from the surrogate dynamics of the backward pass.
Using differentiable linear algebra, we implement surrogate
gradient flossing to push surrogate Lyapunov exponents towards zero. This method slows collective timescales in thetangent space, facilitating error propagation during training. Unlike previous related work that introduced gradient
flossing to shape the Lyapunov exponents of the forward
dynamics and thus improve the trainability of differentiable
dynamical systems [3], we specifically address the surrogate gradients in the backward pass of models that havenon-differentiable error landscapes. We demonstrate thatsurrogate gradient flossing prevents the gradient norm from
exploding or vanishing during training. Empirically, surrogate gradient flossing enhances training on tasks that require
bridging long time horizons.
Prior research that addressed vanishing and exploding gradients for surrogate training focused mostly on optimizing
the surrogate functions of individual neurons [24], [25],or optimized naturalistic network initializations [26]. In
contrast, we studied the collective network state and thetangent dynamics in the backward pass that involves manyneurons and many time steps. We find that the sharpnessof the surrogate gradient g, which shapes how strongly theactivation function is smoothed during the backward pass,plays a crucial role in determining the norm and condition
number of the surrogate gradient. This is because g shapes
the sparsity of the surrogate Jacobian and thus the shape ofthe surrogate Lyapunov spectrum. Compared to a previous
method to control Lyapunov exponents of vanilla RNNs
with continuous activation function [3], here we addresseddiscrete networks and optimized the surrogate gradients relevant to the backward pass of gradient training.
Future studies should evaluate surrogate gradient flossing
in real-world applications and extend testing to complex architectures, including deep spiking networks, spiking transformers, and quantized neural networks.

6. Limitations
When adapting the input weights, recurrent weights, andg during surrogate gradient flossing, we observed that gbecame very small, causing the surrogate function ϕ to become close to linear. In that situation, the binary RNNs werestill able to solve linear tasks, such as a delayed copy task,
but not XOR. Additionally, binary networks with very small
g were sometimes susceptible to oscillations or diverging
gradients. We addressed this by using adaptive gradient
flossing and successively increasing g whenever the maximum surrogate Lyapunov exponent estimate λS1 > 0.2 orT
the gradient norm | ∂L
∂W | exceeded 1. In general, binary
RNNs with surrogate gradients seem to have a trade-off
between the ability to bridge long time horizons and theability to solve nonlinear tasks such as XOR.
We focused our analysis on recurrent binary and spiking
neural networks because their recurrent interactions makethem an interesting and powerful dynamical system. However, the insights from our study could also help address
vanishing and exploding gradients in deeply layered feedforward quantized neural networks, where these issues canbe linked to transient chaos or transient stability [23].
A limitation of the proposed surrogate gradient flossing implementation is that QR decomposition scales with O(N k 2 ),where N is the network size and k is the number of surrogate Lyapunov exponents being flossed. The number ofsurrogate Lyapunov exponents that need to be pushed towards zero√ depends on the task. For tasks requiring morethan O( N ) surrogate Lyapunov exponents close to zero,
the QR decomposition would become the computationalbottleneck. We previously suggested strategies to address
this computational bottleneck [3].

7. Discussion
Training discrete networks is challenging because of thenon-differentiable error landscapes created by discrete interactions. Surrogate gradient training is commonly used
for binary and spiking neural networks, but it often leads to
exploding or vanishing gradients, especially in temporally
4Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

References
[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

Y. Bengio, N. Léonard, and A. Courville, Estimatingor Propagating Gradients Through Stochastic Neurons
for Conditional Computation, arXiv:1308.3432 [cs], Aug.
2013. DOI: 10.48550/arXiv.1308.3432. [Online].
Available: http://arxiv.org/abs/1308.3432(visited on 04/13/2024).
E. O. Neftci, H. Mostafa, and F. Zenke, “Surrogate gradientlearning in spiking neural networks: Bringing the powerof gradient-based optimization to spiking neural networks,”
IEEE Signal Processing Magazine, vol. 36, no. 6, pp. 51–
63, 2019.
R. Engelken, “Gradient Flossing: Improving Gradient Descent through Dynamic Control of Jacobians,” en, May
2023. [Online]. Available: https : / / openreview .
net / forum ? id = jEQRoJzDx8 & referrer =%5BAuthor % 20Console % 5D(%2Fgroup % 3Fid %3DNeurIPS . cc % 2F2023 % 2FConference %2FAuthors % 23your - submissions) (visited on
08/28/2023).
R. Engelken, F. Wolf, and L. F. Abbott, “Lyapunov spectra of chaotic recurrent neural networks,” Physical Review Research, vol. 5, no. 4, p. 043 044, Oct. 2023. DOI:
10.1103/PhysRevResearch.5.043044. [Online].
Available: https : / / link . aps . org / doi / 10 .
1103 / PhysRevResearch . 5 . 043044 (visited on
10/19/2023).
R. Engelken, F. Wolf, and L. F. Abbott, “Lyapunov spectraof chaotic recurrent neural networks,” arXiv:2006.02427[nlin, q-bio], Jun. 2020, arXiv: 2006.02427. [Online].
Available: http://arxiv.org/abs/2006.02427(visited on 06/14/2020).
R. Vogt, M. Puelma Touzel, E. Shlizerman, and G. Lajoie, “On Lyapunov Exponents for RNNs: Understanding
Information Propagation Using Dynamical Systems Tools,”
Frontiers in Applied Mathematics and Statistics, vol. 8,2022, ISSN: 2297-4687. [Online]. Available: https://www.frontiersin.org/articles/10.3389/
fams.2022.818799 (visited on 04/08/2023).
J. Mikhaeil, Z. Monfared, and D. Durstewitz, “On the
difficulty of learning chaotic dynamics with RNNs,”
en, Advances in Neural Information Processing Systems, vol. 35, pp. 11 297–11 312, Dec. 2022. [Online]. Available: https : / / proceedings .
neurips . cc / paper _ files / paper / 2022 /hash / 495e55f361708bedbab5d81f92048dcd Abstract - Conference . html (visited on
04/08/2023).
R. Engelken, “Gradient Flossing: Improving Gradient Descent through Dynamic Control of Jacobians,”
en, Advances in Neural Information Processing Systems, vol. 36, pp. 10 412–10 439, Dec. 2023. [Online]. Available: https : / / proceedings .
neurips . cc / paper _ files / paper / 2023 /hash / 214ce905bf2072535e34b3cf873cbbc8 Abstract - Conference . html (visited on
03/13/2024).
P. A. Sokół, I. Jordan, E. Kadile, and I. M. Park, “Adjoint
Dynamics of Stable Limit Cycle Neural Networks,” in2019 53rd Asilomar Conference on Signals, Systems, andComputers, ISSN: 2576-2303, Nov. 2019, pp. 884–887.
DOI : 10.1109/IEEECONF44664.2019.9049080.

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

5

P. A. Sokół, “Geometry of Learning and Representations
in Neural Networks,” eng, Ph.D. dissertation, Stony Brook
University, May 2023.
I. M. Park, Á. Ságodi, and P. A. Sokół, “Persistent learning
signals and working memory without continuous attractors,”
eng, ArXiv, arXiv:2308.12585v1, Aug. 2023, ISSN: 23318422.
M. Monteforte and F. Wolf, “Dynamical Entropy Production in Spiking Neuron Networks in the Balanced State,”
Physical Review Letters, vol. 105, no. 26, p. 268 104, Dec.
2010. DOI: 10.1103/PhysRevLett.105.268104.
[Online]. Available: http://link.aps.org/doi/
10.1103/PhysRevLett.105.268104 (visited on
06/02/2012).
M. Monteforte and F. Wolf, “Dynamic Flux Tubes Form
Reservoirs of Stability in Neuronal Circuits,” PhysicalReview X, vol. 2, no. 4, p. 041 007, Nov. 2012. DOI:
10 . 1103 / PhysRevX . 2 . 041007. [Online]. Available: http : / / link . aps . org / doi / 10 . 1103 /PhysRevX.2.041007 (visited on 11/17/2012).
A. Palmigiano, R. Engelken, and F. Wolf, “Boosting ofneural circuit chaos at the onset of collective oscillations,” en, eLife, vol. 12, Nov. 2023. DOI: 10 . 7554 /eLife . 90378. [Online]. Available: https : / /elifesciences . org / reviewed - preprints /90378 (visited on 12/28/2023).
R. Engelken, “Chaotic Neural Circuit Dynamics,” deu, Feb.
2018. [Online]. Available: https : / / ediss . uni goettingen . de / handle / 11858 / 00 - 1735 0000-002E-E349-9 (visited on 03/27/2018).
P. Manz, S. Goedeke, and R.-M. Memmesheimer, “Dynamics and computation in mixed networks containing neurons
that accelerate towards spiking,” eng, Physical Review. E,
vol. 100, no. 4-1, p. 042 404, Oct. 2019, ISSN: 2470-0053.
DOI : 10.1103/PhysRevE.100.042404.
R. Engelken, “SparseProp: Efficient Event-Based Simulation and Training of Sparse Recurrent Spiking
Neural Networks,” en, May 2023. [Online]. Available: https : / / openreview . net / forum ?
id = yzZbwQPkmP & referrer = %5BAuthor %20Console % 5D(%2Fgroup % 3Fid % 3DNeurIPS .
cc % 2F2023 % 2FConference % 2FAuthors %23your-submissions) (visited on 08/28/2023).
M. Puelma Touzel and F. Wolf, “Statistical mechanicsof spike events underlying phase space partitioning andsequence codes in large-scale models of neural circuits,”
Physical Review E, vol. 99, no. 5, p. 052 402, May 2019.
DOI : 10 . 1103 / PhysRevE . 99 . 052402. [Online]. Available: https : / / link . aps . org / doi /10 . 1103 / PhysRevE . 99 . 052402 (visited on
02/15/2022).
G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn,“Lyapunov Characteristic Exponents for smooth dynamical systems and for hamiltonian systems; A method forcomputing all of them. Part 2: Numerical application,” en,Meccanica, vol. 15, no. 1, pp. 21–30, Mar. 1980, ISSN:
0025-6455, 1572-9648. DOI: 10.1007/BF02128237.
[Online]. Available: http : / / link . springer .
com/article/10.1007/BF02128237 (visited on
04/14/2014).

Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory
[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

D. P. Kingma and J. Ba, “Adam: A Method for Stochastic
Optimization,” arXiv e-prints, vol. 1412, arXiv:1412.6980,
Dec. 2014. [Online]. Available: http : / / adsabs .
harvard.edu/abs/2014arXiv1412.6980K (visited on 05/10/2019).
J. Bezanson, A. Edelman, S. Karpinski, and V. Shah, “Julia:
A Fresh Approach to Numerical Computing,” SIAM Review,
vol. 59, no. 1, pp. 65–98, Jan. 2017, ISSN: 0036-1445. DOI:
10.1137/141000671. [Online]. Available: http://epubs . siam . org / doi / 10 . 1137 / 141000671(visited on 02/07/2018).
Innes, “Flux: Elegant machine learning with Julia.,” Journal of Open Source Software, vol. 3, no. 25, p. 602, May
2018. [Online]. Available: https : / / doi . org / 10 .
21105/joss.00602.
B. Poole, S. Lahiri, M. Raghu, J. Sohl-Dickstein, and S.
Ganguli, “Exponential expressivity in deep neural networks
through transient chaos,” arXiv:1606.05340 [cond-mat,stat], Jun. 2016, arXiv: 1606.05340. [Online]. Available:
http://arxiv.org/abs/1606.05340 (visited on
10/20/2016).
Y. Li, Y. Guo, S. Zhang, S. Deng, Y. Hai, and S. Gu, “Differentiable Spike: Rethinking Gradient-Descent for Training
Spiking Neural Networks,” in Advances in Neural Information Processing Systems, vol. 34, Curran Associates, Inc.,
2021, pp. 23 426–23 439. [Online]. Available: https :/ / proceedings . neurips . cc / paper / 2021 /hash / c4ca4238a0b923820dcc509a6f75849b Abstract.html (visited on 04/09/2024).
F. Zenke and T. P. Vogels, “The remarkable robustness ofsurrogate gradient learning for instilling complex functionin spiking neural networks,” Neural computation, vol. 33,
no. 4, pp. 899–925, 2021, Publisher: MIT Press One Rogers
Street, Cambridge, MA 02142-1209, USA journals-info. . . [Online]. Available: https://direct.mit.edu/
neco / article - abstract / 33 / 4 / 899 / 97482(visited on 04/09/2024).
J. Rossbroich, J. Gygax, and F. Zenke, “Fluctuation-driven
initialization for spiking neural network training,” Neuromorphic Computing and Engineering, vol. 2, no. 4,
p. 044 016, 2022, Publisher: IOP Publishing. [Online].
Available: https : / / iopscience . iop . org /article/10.1088/2634- 4386/ac97bb/meta(visited on 04/09/2024).
S. F. Walter and L. Lehmann, “Algorithmic Differentiationof Linear Algebra Functions with Application in Optimum
Experimental Design (Extended Version),” Tech. Rep.,
Jan. 2010, ADS Bibcode: 2010arXiv1001.1654W Type:article. [Online]. Available: https : / / ui . adsabs .
harvard.edu/abs/2010arXiv1001.1654W (visited on 04/08/2023).
H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, “Differentiable Programming Tensor Networks,” Physical Review
X, vol. 9, no. 3, p. 031 041, Sep. 2019, arXiv:1903.09650[cond-mat, physics:quant-ph], ISSN: 2160-3308. DOI: 10.
1103 / PhysRevX . 9 . 031041. [Online]. Available:
http://arxiv.org/abs/1903.09650 (visitedon 03/23/2023).
R. Schreiber and C. Van Loan, “A Storage-Efficient $WY$
Representation for Products of Householder Transformations,” SIAM Journal on Scientific and Statistical Computing, vol. 10, no. 1, pp. 53–57, Jan. 1989, ISSN: 0196-5204.
DOI : 10.1137/0910005. [Online]. Available: https:

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

[38]

6

/ / epubs . siam . org / doi / 10 . 1137 / 0910005(visited on 03/23/2023).
M. Seeger, A. Hetzel, Z. Dai, E. Meissner, andN. D. Lawrence, “Auto-Differentiating Linear
Algebra,” Tech. Rep., Oct. 2017, ADS Bibcode:
2017arXiv171008717S Type: article. [Online]. Available:
https : / / ui . adsabs . harvard . edu / abs /2017arXiv171008717S (visited on 04/08/2023).
F. Farkhooi and W. Stannat, “Complete Mean-Field Theory
for Dynamics of Binary Recurrent Networks,” PhysicalReview Letters, vol. 119, no. 20, p. 208 301, Nov. 2017.
DOI : 10.1103/PhysRevLett.119.208301. [Online]. Available: https : / / link . aps . org / doi /10.1103/PhysRevLett.119.208301 (visited on
01/24/2018).
H. Sompolinsky, A. Crisanti, and H. J. Sommers, “Chaos
in Random Neural Networks,” Physical Review Letters,
vol. 61, no. 3, pp. 259–262, Jul. 1988. DOI: 10.1103/
PhysRevLett . 61 . 259. [Online]. Available: http ://link.aps.org/doi/10.1103/PhysRevLett.
61.259 (visited on 01/08/2015).
J. Kadmon and H. Sompolinsky, “Transition to Chaos in
Random Neuronal Networks,” Physical Review X, vol. 5,
no. 4, p. 041 030, Nov. 2015. DOI: 10.1103/PhysRevX.
5.041030. [Online]. Available: http://link.aps.
org/doi/10.1103/PhysRevX.5.041030 (visitedon 11/08/2016).
L. Molgedey, J. Schuchhardt, and H. G. Schuster, “Suppressing chaos in neural networks by noise,” Physical Review Letters, vol. 69, no. 26, pp. 3717–3719, Dec. 1992.
DOI : 10 . 1103 / PhysRevLett . 69 . 3717. [Online]. Available: http : / / link . aps . org / doi /10 . 1103 / PhysRevLett . 69 . 3717 (visited on
12/23/2014).
J. Schuecker, S. Goedeke, and M. Helias, “Optimal Sequence Memory in Driven Random Networks,” Physical Review X, vol. 8, no. 4, p. 041 029, Nov. 2018. DOI:
10 . 1103 / PhysRevX . 8 . 041029. [Online]. Available: https://link.aps.org/doi/10.1103/
PhysRevX.8.041029 (visited on 11/20/2018).
A. Crisanti and H. Sompolinsky, “Path integral approach
to random neural networks,” Physical Review E, vol. 98,
no. 6, p. 062 120, Dec. 2018. DOI: 10.1103/PhysRevE.
98 . 062120. [Online]. Available: https : / / link .
aps.org/doi/10.1103/PhysRevE.98.062120(visited on 07/31/2019).
K. Rajan, L. F. Abbott, and H. Sompolinsky, “Stimulusdependent suppression of chaos in recurrent neural networks,” Physical Review E, vol. 82, no. 1, p. 011 903,
Jul. 2010. DOI: 10 . 1103 / PhysRevE . 82 . 011903.
[Online]. Available: http : / / link . aps . org /doi/10.1103/PhysRevE.82.011903 (visited on
11/09/2012).
R. Engelken, A. Ingrosso, R. Khajeh, S. Goedeke,
and L. F. Abbott, “Input correlations impede suppression of chaos and learning in balanced firing-rate networks,” en, PLOS Computational Biology, vol. 18,
no. 12, e1010590, Dec. 2022, ISSN: 1553-7358. DOI:
10 . 1371 / journal . pcbi . 1010590. [Online].
Available: https : / / journals . plos . org /ploscompbiol / article ? id = 10 . 1371 /journal.pcbi.1010590 (visited on 12/23/2022).

Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory
[39]

[40]

[41]

[42]

[43]

[44]

[45]

[46]

[47]

[48]

C. M. Newman, “The distribution of Lyapanov exponents:
Exact results for random matrices,” eng, Communications
in mathematical physics, vol. 103, no. 1, pp. 121–126,
1986, ISSN: 0010-3616. [Online]. Available: http : / /cat.inist.fr/?aModele=afficheN&cpsidt=
7875417 (visited on 10/22/2016).
M. Isopi and C. M. Newman, “The triangle law for Lyapunov exponents of large random matrices,” eng, Communications in mathematical physics, vol. 143, no. 3,
pp. 591–598, 1992, ISSN: 0010-3616. [Online]. Available:
http://cat.inist.fr/?aModele=afficheN&
cpsidt=5111053 (visited on 10/22/2016).
M. Bauer and W. Martienssen, “Lyapunov exponents anddimensions of chaotic neural networks,” en, Journal of
Physics A: Mathematical and General, vol. 24, no. 19,
p. 4557, 1991, ISSN: 0305-4470. DOI: 10.1088/03054470 / 24 / 19 / 019. [Online]. Available: http : / /stacks.iop.org/0305-4470/24/i=19/a=019(visited on 12/17/2016).
G. Curato and A. Politi, “Onset of chaotic dynamics in neural networks,” Physical Review E, vol. 88, no. 4, p. 042 908,
Oct. 2013. DOI: 10.1103/PhysRevE.88.042908.
[Online]. Available: http : / / link . aps . org /doi/10.1103/PhysRevE.88.042908 (visited on
11/28/2013).
A. Lerchner, C. Ursta, J. Hertz, M. Ahmadi, P. Ruffiot, and
S. Enemark, “Response Variability in Balanced Cortical
Networks,” Neural Computation, vol. 18, no. 3, pp. 634–
659, Mar. 2006, ISSN: 0899-7667. DOI: 10.1162/neco.
2006.18.3.634. [Online]. Available: https://doi.
org/10.1162/neco.2006.18.3.634 (visited on
06/20/2023).
R. Engelken, F. Farkhooi, D. Hansel, C. van Vreeswijk,
and F. Wolf, “A reanalysis of “Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons”,” en, F1000Research, vol. 5,
p. 2043, Aug. 2016, ISSN: 2046-1402. DOI: 10.12688/
f1000research.9144.1. [Online]. Available: http://f1000research.com/articles/5-2043/v1(visited on 09/28/2016).
C. van Vreeswijk and H. Sompolinsky, “Chaos in Neuronal
Networks with Balanced Excitatory and Inhibitory Activity,” Science, vol. 274, no. 5293, pp. 1724–1726, Dec. 1996.
DOI : 10 . 1126 / science . 274 . 5293 . 1724. [Online]. Available: http://www.sciencemag.org/
content/274/5293/1724.abstract (visited on
12/11/2011).
C. van Vreeswijk and H. Sompolinsky, “Chaotic Balanced
State in a Model of Cortical Circuits,” Neural Computation, vol. 10, no. 6, pp. 1321–1371, 1998, ISSN: 08997667. DOI: 10.1162/089976698300017214. [Online]. Available: http://dx.doi.org/10.1162/
089976698300017214.
F. Zenke and S. Ganguli, “SuperSpike: Supervised Learning in Multilayer Spiking Neural Networks,” Neural Computation, vol. 30, no. 6, pp. 1514–1541, Jun. 2018, ISSN:
0899-7667. DOI: 10.1162/neco_a_01086. [Online].
Available: https://doi.org/10.1162/neco_a_01086 (visited on 05/13/2023).
A. M. Saxe, J. L. McClelland, and S. Ganguli, “Exact
solutions to the nonlinear dynamics of learning in deep
linear neural networks,” arXiv:1312.6120 [cond-mat, q-bio,stat], Dec. 2013, arXiv: 1312.6120. [Online]. Available:

[49]

[50]

[51]

[52]

[53]

[54]

[55]

[56]

[57]

[58]

[59]

7

http://arxiv.org/abs/1312.6120 (visited on
11/12/2016).
J. Pennington, S. Schoenholz, and S. Ganguli, “Resurrecting the sigmoid in deep learning through dynamical isometry: Theory and practice,” in Advances in Neural Information Processing Systems, vol. 30, Curran Associates, Inc., 2017. [Online]. Available: https : / /proceedings . neurips . cc / paper / 2017 /hash / d9fc0cdb67638d50f411432d0d41d0ba Abstract.html (visited on 05/17/2023).
J. Pennington, S. S. Schoenholz, and S. Ganguli, “TheEmergence of Spectral Universality in Deep Networks,”
arXiv:1802.09979 [cs, stat], Feb. 2018, arXiv: 1802.09979.
[Online]. Available: http://arxiv.org/abs/1802.
09979 (visited on 08/02/2019).
E. Cheney and D. Kincaid, Numerical Mathematics andComputing, en. Cengage Learning, Aug. 2007, GoogleBooks-ID: ZUfVZELlrMEC, ISBN: 978-0-495-11475-8.
Y. Bahri, J. Kadmon, J. Pennington, S. S. Schoenholz,
J. Sohl-Dickstein, and S. Ganguli, “Statistical Mechanics of Deep Learning,” Annual Review of Condensed Matter Physics, vol. 11, no. 1, pp. 501–528, 2020. DOI: 10.
1146/annurev- conmatphys- 031119- 050745.
[Online]. Available: https://doi.org/10.1146/
annurev- conmatphys- 031119- 050745 (visitedon 05/17/2023).
X. Glorot, A. Bordes, and Y. Bengio, “Deep Sparse
Rectifier Neural Networks,” en, vol. 15, Apr. 2011,
pp. 315–323. [Online]. Available: https : / / hal .
archives-ouvertes.fr/hal-00752497 (visitedon 02/23/2019).
M. Chen, J. Pennington, and S. S. Schoenholz, “Dynamical Isometry and a Mean Field Theory of RNNs: Gating
Enables Signal Propagation in Recurrent Neural Networks,”
arXiv:1806.05394 [cs, stat], Aug. 2018, arXiv: 1806.05394.
[Online]. Available: http://arxiv.org/abs/1806.
05394 (visited on 06/01/2020).
B. Hanin and M. Nica, “Products of Many Large Random
Matrices and Gradients in Deep Neural Networks,” en,Dec. 2018. [Online]. Available: https://arxiv.org/
abs/1812.05994v1 (visited on 03/26/2019).
S. S. Schoenholz, J. Gilmer, S. Ganguli, and J. SohlDickstein, “Deep Information Propagation,” Nov. 2016.
[Online]. Available: https : / / openreview . net /forum?id=H1W1UN9gg (visited on 08/02/2019).
B. Hanin and D. Rolnick, “How to Start Training: The
Effect of Initialization and Architecture,” in Advancesin Neural Information Processing Systems, vol. 31, Curran Associates, Inc., 2018. [Online]. Available: https:/ / proceedings . neurips . cc / paper / 2018 /hash / d81f9c1be2e08964bf9f24b15f0e4900 Abstract.html (visited on 10/28/2023).
P. A. Sokol and I. Memming Park, “Information Geometry of Orthogonal Initializations and Training,” Tech. Rep.,
Oct. 2018, ADS Bibcode: 2018arXiv181003785S Type:article. [Online]. Available: https : / / ui . adsabs .
harvard.edu/abs/2018arXiv181003785S (visited on 05/17/2023).
D. Gilboa, B. Chang, M. Chen, et al., “Dynamical Isometry
and a Mean Field Theory of LSTMs and GRUs,” en, Jan.
2019. [Online]. Available: https : / / arxiv . org /abs/1901.08987v1 (visited on 03/26/2019).

Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory
[60]

[61]

[62]

[63]

[64]

[65]

[66]

[67]

[68]

[69]

[70]

T. Can, K. Krishnamurthy, and D. J. Schwab, “Gating
creates slow modes and controls phase-space complexity in GRUs and LSTMs,” arXiv:2002.00025 [cond-mat,stat], Jan. 2020, arXiv: 2002.00025. [Online]. Available:
http://arxiv.org/abs/2002.00025 (visited on
05/27/2020).
D. Doshi, T. He, and A. Gromov, “Critical Initialization ofWide and Deep Neural Networks using Partial Jacobians:
General Theory and Applications,” en, Nov. 2023. [Online].
Available: https : / / openreview . net / forum ?
id=wRJqZRxDEX (visited on 01/16/2024).
S.
Hochreiter,“Untersuchungen zu dynamischen neuronalen Netzen,” deu, Ph.D. dissertation, 1991. [Online]. Available: https : / /www . semanticscholar . org / paper /Untersuchungen - zu - dynamischen neuronalen - Netzen - Hochreiter /3f3d13e95c25a8f6a753e38dfce88885097cbd43(visited on 09/23/2021).
S. Hochreiter and Jürgen Schmidhuber, “Long Short-Term
Memory,” Neural Computation, vol. 9, no. 8, pp. 1735–
1780, 1997.
K. Cho, B. van Merrienboer, D. Bahdanau, and Y. Bengio,“On the Properties of Neural Machine Translation: EncoderDecoder Approaches,” Tech. Rep., Sep. 2014, ADS Bibcode: 2014arXiv1409.1259C Type: article. [Online]. Available: https://ui.adsabs.harvard.edu/abs/
2014arXiv1409.1259C (visited on 04/08/2023).
R. Pascanu, T. Mikolov, and Y. Bengio, “On the difficultyof training recurrent neural networks,” in Proceedings ofthe 30th International Conference on International Conference on Machine Learning - Volume 28, ser. ICML’13,
Atlanta, GA, USA: JMLR.org, Jun. 2013, pp. III–1310–III–
1318. (visited on 06/01/2020).
T. Mikolov, “Statistical language models based on neural networks,” Ph.D. thesis, Brno University of Technology, Faculty of Information Technology, Brno, CZ, 2012.
[Online]. Available: https://www.fit.vut.cz/
study/phd-thesis/283/.
I. Goodfellow, Y. Bengio, and A. Courville, Deep
Learning, en. MIT Press, Nov. 2016, Google-Books-ID:
omivDQAAQBAJ, ISBN: 978-0-262-33737-3.
B. Chang, M. Chen, E. Haber, and E. H. Chi, “AntisymmetricRNN: A Dynamical System View on Recurrent Neural Networks,” en, International Conference on Learning
Representations, Dec. 2018. [Online]. Available: https:/ / openreview . net / forum ? id = ryxepo0cFX(visited on 05/12/2023).
M. Arjovsky, A. Shah, and Y. Bengio, “Unitary Evolution Recurrent Neural Networks,” en, in Proceedings of
The 33rd International Conference on Machine Learning, PMLR, Jun. 2016, pp. 1120–1128. [Online]. Available: https://proceedings.mlr.press/v48/
arjovsky16.html (visited on 04/08/2023).
K. Helfrich, D. Willmott, and Q. Ye, “Orthogonal Recurrent Neural Networks with Scaled Cayley Transform,” en,in Proceedings of the 35th International Conference on Machine Learning, PMLR, Jul. 2018, pp. 1969–1978. [Online].
Available: https://proceedings.mlr.press/
v80/helfrich18a.html (visited on 04/08/2023).

8

[71]

[72]

[73]

[74]

[75]

[76]

[77]

[78]

[79]

[80]

T. Konstantin Rusch and S. Mishra, “Coupled Oscillatory
Recurrent Neural Network (coRNN): An accurate and (gradient) stable architecture for learning long time dependencies,” en, arXiv e-prints, arXiv:2010.00951, Oct. 2020.
DOI : 10 . 48550 / arXiv . 2010 . 00951. [Online].
Available: https://ui.adsabs.harvard.edu/
abs / 2020arXiv201000951K / abstract (visitedon 04/08/2023).
N. B. Erichson, O. Azencot, A. Queiruga, L. Hodgkinson, and M. W. Mahoney, “Lipschitz Recurrent Neural
Networks,” en, International Conference on Learning Representations, Jan. 2021. [Online]. Available: https://openreview.net/forum?id=-N7PBXqOUJZ (visited on 05/12/2023).
A. Gu, K. Goel, and C. Ré, Efficiently Modeling Long Sequences with Structured State Spaces, arXiv:2111.00396[cs], Aug. 2022. DOI: 10 . 48550 / arXiv . 2111 .
00396. [Online]. Available: http : / / arxiv . org /abs/2111.00396 (visited on 01/16/2024).
J. T. H. Smith, A. Warrington, and S. W. Linderman,
Simplified State Space Layers for Sequence Modeling,
arXiv:2208.04933 [cs], Mar. 2023. DOI: 10 . 48550 /arXiv . 2208 . 04933. [Online]. Available: http :/ / arxiv . org / abs / 2208 . 04933 (visited on
01/16/2024).
A. Orvieto, S. L. Smith, A. Gu, et al., “Resurrecting Recurrent Neural Networks for Long Sequences,” Tech. Rep.,
Mar. 2023, ADS Bibcode: 2023arXiv230306349O Type:article. [Online]. Available: https : / / ui . adsabs .
harvard.edu/abs/2023arXiv230306349O (visited on 04/08/2023).
A. Gu and T. Dao, Mamba: Linear-Time Sequence Modeling with Selective State Spaces, arXiv:2312.00752 [cs],Dec. 2023. DOI: 10 . 48550 / arXiv . 2312 . 00752.
[Online]. Available: http://arxiv.org/abs/2312.
00752 (visited on 01/16/2024).
H. Jaeger, “The “echo state” approach to analysing andtraining recurrent neural networks-with an erratum note,”Bonn, Germany: German National Research Center for Information Technology GMD Technical Report, vol. 148,
p. 34, 2001. [Online]. Available: http : / / minds .
jacobs - university . de / sites / default /files/uploads/papers/EchoStatesTechRep.
pdf (visited on 10/30/2016).
M. C. Ozturk, D. Xu, and J. C. Prı́ncipe, “Analysis andDesign of Echo State Networks,” Neural Computation,
vol. 19, no. 1, pp. 111–138, Jan. 2007, ISSN: 0899-7667.
DOI : 10.1162/neco.2007.19.1.111.
R. K. Srivastava, K. Greff, and J. Schmidhuber, “Training
Very Deep Networks,” in Advances in Neural Information Processing Systems, vol. 28, Curran Associates, Inc.,
2015. [Online]. Available: https://proceedings.
neurips . cc / paper _ files / paper / 2015 /hash / 215a71a12769b056c3c32e7299f1c5ed Abstract.html (visited on 05/17/2023).
J. G. Zilly, R. K. Srivastava, J. Koutnik, and J. Schmidhuber, “Recurrent Highway Networks,” en, in Proceedingsof the 34th International Conference on Machine Learning, PMLR, Jul. 2017, pp. 4189–4198. [Online]. Available: https://proceedings.mlr.press/v70/
zilly17a.html (visited on 05/17/2023).

Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory
[81]

[82]

K. Geist, U. Parlitz, and W. Lauterborn, “Comparison ofDifferent Methods for Computing Lyapunov Exponents,”
en, Progress of Theoretical Physics, vol. 83, no. 5, pp. 875–
893, May 1990, ISSN: 0033-068X, 1347-4081. DOI: 10.
1143 / PTP . 83 . 875. [Online]. Available: http : / /ptp . oxfordjournals . org / content / 83 / 5 /
875 (visited on 08/05/2016).
V. I. Oseledets, “A multiplicative ergodic theorem. Characteristic Ljapunov, exponents of dynamical systems,”
ru, Trudy Moskovskogo Matematicheskogo Obshchestva,
vol. 19, pp. 179–210, 1968. [Online]. Available: http://mi.mathnet.ru/eng/mmo/v19/p179.

9

Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

A. Algorithm and Pseudocode for Surrogate Gradient Flossing
We implement Surrogate Gradient Flossing by making the established method for calculating Lyapunov exponents
differentiable and combining it with analytical surrogate Jacobians. In the established reorthonormalization method,an orthonormal system Qs = [q1s , q2s , . . . , qks ] is evolved in the tangent space along the network trajectory using thesurrogates+1e
Qs at each time step. .
surrogate Jacobian Dsurrogate= ∂hs∂hs . This process involves the calculation of: Qs+1 = Ds
s+1e
Subsequently, the exponential growth rates are extracted via QR decomposition: Qs+1 = Qs+1 R , which decomposes
e s+1 into an orthonormal matrix Qs+1 and an upper triangular matrix Rs+1 with positive diagonal elements.
Q
Surrogate Lyapunov exponents are then determined by the time-averaged logarithms of the diagonal entries of Rs :The following pseudocode implements the surrogate gradient flossing algorithm, with surrogate gradient function steps in
red and conventional forward dynamics in blue.

B. Backpropagation Through QR Decomposition
The adjoint of the QR decomposition is given by [8], [27]–[30]
Q = Q + Q copyltu(M) R−T ,
T

(4)

T

where M = RR − Q Q and the copyltu function generates a symmetric matrix by copying the lower triangle of the
input matrix to its upper triangle, with the element [copyltu(M )]ij = Mmax(i,j),min(i,j) [27]–[30]. Adjoint variable are
written here as T = ∂L/∂T .
Using an analytical pullback is more memory-efficient and less computationally costly than directly performing automatic differentiation through the QR decomposition. We provide implementations of surrogate gradient flossing
both in Julia (using the Julia package BackwardsLinalg.jl by Jinguo Liu available here) and in PyTorch available athttps://github.com/RainerEngelken/SurrogateGradientFlossing.

C. Details on Surrogate Gradient Flossing for Binary Networks
We consider binary recurrent networks [31], with parallel update of the N units that follows:
hs+1 = W sgn(hs ) + Vxs+1 .

(5)where sgn denotes the sign function. The initial entries of W are drawn independently from a Gaussian distribution withzero mean and variance 1/N . In the backward pass of binary networks, we replace the pullback of sgn with the surrogate
10Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theorys+1′
function ϕ′g (x) := sech(gx)2 . The surrogate Jacobians are given by Dsurrogate= ∂hs∂hs = W ⊙ ϕ (hs ), where ⊙ denotes
the Hadamard product.

For large random networks, the largest Lyapunov exponent λmax can be obtained from dynamic mean-field theory, as
described earlier in both the discrete and continuous-time cases with constant input and additive Gaussian white noise drive[32]–[36] or other time-varying inputs [35], [37], [38]. These insights can be leveraged for surrogate gradient training. For
large binary networks, not only the first but also the full surrogate Lyapunov spectrum can be approximated analytically at
initialization for arbitrary g [4] based on results on product of uncorrelated Gaussian matrices, whose eigenvalue distribution
approximately follows a triangle law [39], [40]. This results in:(i
λmax + 21 log 1 − pN
, i ≤ pN
λi =
(6)
−∞,
i > pN
1where p = erf √2g∆
is the average fraction of neurons that are not in the saturated regime of the surrogate gra0
dient. This fraction depends on both the sharpness g and the variance of h, which can be obtained using dynamic
mean-field theory [4], [34]. This extends earlier work [4], [41], [42], which used the triangular law for the complete
Lyapunov spectrum of the forward dynamics of discrete-time recurrent neural networks, to an analysis of the gradientnorm and condition number in the backward pass. We provide code for surrogate gradient flossing in binary networks athttps://github.com/RainerEngelken/SurrogateGradientFlossing. We find that surrogate gradient flossing also helps train
binary networks on tasks with many time steps, such as the copy task and the temporally delayed XOR task. A full analyticaldescription of surrogate gradient flossing for binary networks would be a promising avenue for future research, as networks
with binary dynamics can still have complex nonlinear learning dynamics. Moreover, quantized networks are increasinglyof interest because of their memory efficiency. However, this is beyond the scope of the work presented here.

D. Additional Details on Spiking Network Model and Parameter Initialization
In the case of spiking networks, we considered a time-discretized implementation of a network of N leaky integrate-and-fireneurons with exponentially decaying synaptic currents, following the notation of [2], [26]:Ui [n + 1] = (βUi [n] + (1 − β)Ii [n])(1 − Sirec [n])

(7)



∆t
where β ≡ exp − τmem
reflects the decay due to leak, Ui is the membrane potential of neuron i, Ii is the synaptic current.
Similarly, the synaptic currents are updated as follows:
X
X
Wij Sjin [n](8)Vij Sjrec [n] +Ii [n + 1] = αIi [n] +jjUnless noted otherwise, we used the default parameters ∆t = 2 ms, τmem = 10 ms, and τmem = 5 ms. We chosethese time constants for both the recurrent neurons and the readout neurons. Our implementation in Julia [21] using Flux[22] is available at https://github.com/RainerEngelken/SurrogateGradientFlossing. We initialized the recurrent networkin the fluctuation-driven regime following [26] by adjusting the input weights Wij and recurrent weights Vij to induce afluctuation-driven regime at initialization.
We note that during task training, LIF networks can develop unphysiological negative voltage excursions if the input variancebecomes large, inducing bursty network activity states with a coefficient of variation of the interspike interval distribution
larger than one, indicating super-Poissonian spiking statistics [43], [44]. It will be important to further study surrogate
gradient flossing and spiking network training when constraining networks to a biologically plausible activity regime andwhen imposing further biological constraints like Dale’s principle and dynamical balance of excitatory and inhibitory input
currents [45], [46].

E. Analytical Surrogate Jacobian of the Spiking Neural Networkn+1
To calculate surrogate Lyapunov exponents, we first need to obtain the surrogate Jacobian, Dsurrogate= ∂hn∂hn , either throughautomatic differentiation or analytically. The surrogate Jacobians evaluated along the forward trajectory describe how an11Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

infinitesimal perturbation of the network state (e.g., here all membrane potentials and synaptic currents) evolves in one time
step.
Here, we give the analytical expression of the 2N × 2N surrogate Jacobian for the discretized dynamics of a recurrent
network of leaky integrate-and-fire neurons with exponential synapses. We only consider the dynamics of the recurrent
neurons and not of the readout neurons, as these are nonspiking and only passively driven and therefore have trivial surrogate
1.
Lyapunov exponents of the relaxation time constants − τ1syn and − τmem
δhn+1 = Dsurrogate
δhn
n

(9)

The surrogate Jacobian matrix, Dsurrogate, is block-structured, according to the partial derivatives of the membrane potentialsnand synaptic currents between subsequent time steps:!
∂U[n+1]∂U[n+1]
Dsurrogate
=
n

∂U[n]
∂I[n+1]
∂U[n]

∂I[n]
∂I[n+1]
∂I[n]

(10)

Often in surrogate gradient training, the voltage reset term (1 − Si ) is ignored in the backward pass [2]. Here, we provide thesurrogate Jacobian entries both including and excluding the voltage reset term. If the voltage reset is ignored, the Jacobian
entries are given by:∂Ui [n + 1]
= βδij + Vij ϕ′g (Uj [n] − ϑ)∂Uj [n]∂Ui [n + 1]
= αδij∂Ij [n]∂Ii [n + 1]
= Vij ϕ′g (Uj [n] − ϑ)∂Uj [n]∂Ii [n + 1]
= αδij∂Ij [n]

(11)
(12)
(13)
(14)

1
ϕ′g (x) is the surrogate gradient function. We use here the derivative of the fast sigmoid ϕ′g (x) = (g|x|+1)
2 [47]. Note that
we absorbed the term (1 − β) from Eq 7 into the input weights and recurrent weights.

If the voltage reset term is not ignored in the derivative, the Jacobian entries are:
∂Ui [n + 1]
= δij β + Vij ϕ′g (Uj [n] − 1) (1 − Si [n])∂Uj [n]− (βUi [n] + αIi [n] + xi [n] + Iirec [n]) ϕ′g (Ui [n] − 1),∂Ui [n + 1]
= δij α (1 − Si [n]) ,∂Ij [n]∂Ii [n + 1]
= Vij ϕ′g (Uj [n] − 1),∂Uj [n]∂Ii [n + 1]
= δij α.
∂Ij [n]P
P
where Iirec [n] = j Vij ϕg (Uj [n] − ϑ) and xi [n] = j Wij Sjin [n].

(15)
(16)
(17)
(18)

To corroborate these analytical expressions, we compare them in the following to expressions of the Jacobian obtained
through automatic differentiation:

F. Linking Surrogate Lyapunov Exponents & Gradient Norm
We here show in more detail that surrogate Lyapunov exponents give access to the norm and structure of the surrogate
gradients for long time horizons, by utilizing a recently discovered link between Lyapunov exponents of RNNs and the
exploding and vanishing gradient problem for the case of surrogate gradient training of spiking networks. [4], [6]–[11].
12

10

20

15
1

10

DijAD|)

10

58

1

60

10

62

64

20

66

log2 (|

20

5

postsynaptic neuron i

10

C

1

log2 (|

10

Dijanalytical|)5

postsynaptic neuron i

B1

log2 (|

postsynaptic neuron iA

DijAD Dijanalytical|)Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

15

20

1

presynaptic neuron j

10

20

1

presynaptic neuron j

1020

presynaptic neuron j

10

20

15
1

10

DijAD|)

10

1
60

10

62

64

20

log2 (|

20

5

postsynaptic neuron i

10

C

1

log2 (|

10

Dijanalytical|)5

postsynaptic neuron i

B1

log2 (|

postsynaptic neuron iA

DijAD Dijanalytical|)Figure 4. Comparison of analytical and numerical surrogate Jacobians A) Analytical surrogate Jacobian for g = 1. B) surrogate
Jacobian obtained from automatic differentiation for g = 10. C) log2 of absolute value of difference between AD and analytical surrogate
analyticalAD
Jacobian. The difference is close to machine precision. log2 (|Dij
− Dij
|)

15

20

1

presynaptic neuron j

10

20

1

presynaptic neuron j

1020

presynaptic neuron j

DijAD|)

10

10

15
20
20

25
1

10

20

presynaptic neuron j

1
62.5

65.0
10
67.5

70.0
20

72.5

log2 (|

20
20

5

postsynaptic neuron i

15

C

1

log2 (|

10

10

Dijanalytical|)5

postsynaptic neuron i

B1

log2 (|

postsynaptic neuron iA

DijAD Dijanalytical|)Figure 5. Same as 4 for g = 10

25
1

10

20

1

presynaptic neuron j

1020

presynaptic neuron j

Figure 6. Same as 4 for g = 100

Note that the complete error gradient of backpropagation through time in Eq 2 is composed of a summation of productsof surrogate Jacobians, reflecting the number of ”loops” the error signal traverses through the recurrent spiking dynamics
before reaching its target. Consequently, when the leading singular values of the surrogate Jacobians are smaller than 1, theinfluence of the shorter loops typically dominates the surrogate gradients.
In the case of vanishing surrogate gradients, the gradient norm is dominated by the shorter loops, even though the actual
signal in the surrogate gradient originates from the loop of duration T , where T is the maximum temporal delay relevant forsolving the temporal credit assignment of the task. To mitigate the contamination of spurious signals from shorter loops andeffectively extract the gradient that spans long time horizons, we focus on the gradient at the last time point LT with respect
to the initial conditions h0 .
!
τ =T −1 TY
−1
T −1
∂LT
∂hτ ′ +1 ∂hτ∂LT X∂LT Y ∂hτ ′ +1
=
=
(19)
∂h0
∂hT τ =0
∂hτ ′
∂h0
∂hT ′
∂hτ ′
′
τ =τ

τ =0

We note that the sum conveniently drops the only summand that contributes is the product of surrogate Jacobians goingtfrom 0 to T . We note that we considered the binary cross entropy loss which makes the derivative ∂L∂ht trivial. Using the13Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

definition of surrogate Lyapunov exponents in Eq 3 and the definition of the spectral norm, we obtainTY−1τ ′ =0∂hτ ′ +1∂hτ ′2

= sup ∥∥x∥2 =1

TY−1

τ ′ =0∂hτ ′ +1x∥2 = σ1
∂hτ ′

TY
−1

τ ′ =0Dsurrogateτ′!

≈ exp(λS1 T )ST
we can thus approximate the surrogate gradient norm | ∂L
∂h0 | ∝ exp(λ1 T ).

G. Linking Surrogate Lyapunov Exponents with Gradient Dimensionality
We pointed out in the main text that a surrogate Lyapunov exponents correspond to exponentially exploding gradient modes,
while negative surrogate Lyapunov exponents correspond to exponentially vanishing surrogate gradient modes.
A well-conditioned surrogate Jacobian is essential for efficient and fast learning [48]–[50]. Surrogate gradient flossing
improves the condition number of the long-term Jacobian which constrains the error signal propagation across long time
horizons in backpropagation [3].The condition number κ2 of a linear map A measures how close the map is to being(A)
singular and is given by the ratio of the largest singular value σmax and the smallest singular values σmin , so κ2 (A) = σσmax.
min (A)p
According to the rule of thumb given in [51], if κ2 (A) = 10 , one can anticipate losing at least p digits of precision whensolving the equation Ax = b. Note that the long-term Jacobian Tt is composed of a product of surrogate Jacobians, which
generically makes it ill-conditioned.
Our theoretical estimate of the condition number κ2 of an orthonormal system Q of size N × m that is temporally evolvedby the long-term Jacobian Tte t+τ ) = κ2 Tt (hτ )Qt = σ1 (Tt (hτ )) ≈ exp ((λ1 − λm )(t − τ )) .
κ2 ( Qσm (Tt (hτ ))(20)

where σ1 (Tt (hτ )) and σm (Tt (hτ )) are the first and mth singular value of the long-term Jacobian. We note that thistheoretical estimate of the condition number follows from the asymptotic definition of Lyapunov exponents and should be
exact in the limit of long times. Given that gradient flossing reduces the condition number by a factor whose magnitude
increases exponentially with T , we can expect that surrogate gradient flossing has a stronger effect on problems with a longtime horizon to bridge.
Moreover, surrogate Lyapunov exponents enable the estimation of the number of gradient dimensions available for the
backpropagation of error signals. Generally, the long-term Jacobian is ill-conditioned, however, the Lyapunov spectrum
provides for a 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 a fixed acceptable condition
number—determined, for example, by noise level or floating-point precision—we observe that gradient flossing increases
the number of usable tangent space dimensions for backpropagation.
We stress that these insights go beyond mean-field methods, which usually only consider uncorrelated i.i.d. weight matrices
and become exact only in the large-network limit N → ∞ [52].

H. Further Relation to Previous Literature
Generally, most previous analyses of the notorious vanishing and exploding gradient problem addressed differentiableneural networks and considered different strategies against vanishing and exploding gradients, for instance suitable choiceof weights that guarantee well-conditioned Jacobian at initialization [23], [48], [53]–[61], specialized neuron models whose
gating interactions shield the latent memory state and which can therefore transport information across long time horizons
steps [62]–[64] or pragmatic training tricks like gradient clipping, which re-scales the gradient norm [65] or their individual
elements [66] if they become too large [67]. Moreover, specialized network architectures were introduced to tackle gradientissues, for example, antisymmetric networks [68], orthogonal/unitary initializations [48], [69], [70], coupled oscillatory
RNNs [71], Lipschitz RNNs [72], state-space models [73]–[76], echo state networks [77], [78], (recurrent) highway networks[79], [80], and stable limit cycle neural networks [9]–[11]. The work presented here can be seen as complementary to thelatter.
14

Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

I. Further Details and Analysis of Surrogate Gradient Flossing
An example implementation of surrogate gradient flossing in PyTorch and Flux [22], a machine learning library in Julia [21]is available at https://github.com/RainerEngelken/SurrogateGradientFlossing. We also provide code to directly reproduce
Fig 1.

J. Computational Complexity of Surrogate Gradient Flossing
We present here a more in-depth scaling analysis of the computational cost of surrogate gradient flossing. There arethree main contributors to the computational
cost (table 1): First, the forward dynamics of the spiking RNN, which has
a computational complexity of O N 2 b per time step, where N is the dimension of the recurrent network state (whichin the case of leaky integrate-and-fire neurons with exponentially decaying synapses is twice the number of units) and bis the batch size both in the forward and backward pass. Second, the Jacobian step which scales with O N 2 k per time
step, where k is the number of flossed surrogate Lyapunov exponents. Third, the QR decomposition, which scales with
O N k 2 , where k is the number of surrogate Lyapunov exponents considered.

Together, this results in a total amortized cost of O N 2 b T per training epoch,where T is the number of training time
steps and a total amortized costs per flossing epoch of O N 2 Tf (1 + k/tONS + k) where Tf is the number of flossing time
steps.
In the case of surrogate preflossing, thus, the total computation cost scale with O N 2 [Eb T + Ep Tf (1 + k/tONS + k)] ,where E is the number of training epochs and Ep is the number of surrogate preflossing epochs.
For surrogate gradient flossing during training (assuming that there
is also surrogate preflossing done), the amortizedcost scale with O N 2 [Eb T + Ep Tp + Ef Tf (1 + k/tONS + k)] , where Ef is the total number of flossing epochs during
training.
√
We note that as long as k < N the QR decomposition is not the computational bottleneck.
Moreover, we find empirically that both the number of surrogate preflossing epochs Ep and surrogate flossing episodes Efnecessary for training success is much smaller than the total number of training epochs E. It remains an important challenge
to infer the suitable number of surrogate flossing time steps Tf for tasks with unknown temporal correlation structure.
It would also be interesting to investigate how the CPU hours/wall-clock time/flops/Joule/CO2-emission spent on surrogate
gradient flossing vs on training networks with larger N are trading off against each other. For this, we would suggest firstfinding the smallest network that on median successfully trains on a binary temporal XOR task for a fixed given delay T andmeasuring the computational resources involved in training it, e.g. in terms of CPU hours. Then compare it to a networkwith surrogate gradient flossing.

K. Additional Background on Lyapunov Exponents of RNNs
An autonomous dynamical system is usually defined by a set of ordinary differential equations dh/dt = F(h), h ∈ RN in
the case of continuous-time dynamics, or as a map hs+1 = f (hs ) in the case of discrete-time dynamics. In the following, thetheory is presented for discrete-time dynamical systems for ease of notation, but everything directly extends to continuoustime systems [81]. Together with an initial condition h0 , the map forms a trajectory. As a natural extension of linear
stability analysis, one can ask how an infinitesimal perturbation h′0 = h0 + ϵu0 evolves in time. Chaotic systems are
sensitive to initial conditions; almost all infinitesimal perturbations ϵu0 of the initial condition grow exponentially withtime |ϵut | ≈ exp(λ1 t)|ϵu0 |. Finite-size perturbations, therefore, may lead to a drastically different subsequent behavior.
The largest Lyapunov exponent λ1 measures the average rate of exponential divergence or convergence of nearby initialconditions:1||ϵut ||λ1 (h0 ) = lim lim log(21)t→∞ t ϵ→0||ϵu0 ||
In dynamical systems that are ergodic on the attractor, the Lyapunov exponents do not depend on the initial conditions aslong as the initial conditions are in the basins of attraction of the attractor. Note that it is crucial to take the limit ϵ → 0 first||ϵut ||and then t → ∞, as λ1 (h0 ) would be trivially zero for a bounded attractor if the limits are exchanged, as limt→∞ log ||ϵu0 ||is bounded for finite perturbations even if the system is chaotic. To measure k Lyapunov exponents, one has to study theevolution of k independent infinitesimal perturbations us spanning the tangent space:
15

Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

RNN dynamics
Jacobian step
QR step
total amortized costs
per training epoch

forward pass
O N2 b
O N2 k
O N k2
O N2 b T

backward pass””””

total amortized costs
per surrogate gradient flossing epochO N 2 Tf (1 + k/tONS + k)

total amortized costsof surrogate preflossing

O N 2 [Eb T + Ep Tf (1 + k/tONS + k)]

total amortized costs
flossing during training

O N 2 [Eb T + Ep Tp + Ef Tf (1 + k/tONS + k)]



”



”


”

Table 1. Computational cost for surrogate gradient flossing and training of RNNs
N denotes number of neurons, b is the batch size, T is the number of time steps in forward pass of training, Tf is the number of time steps
in forward pass of flossing, tONS is the reorthonormalization interval, k is the number of flossed Lyapunov exponents, E is the number oftraining epochs, Ep is the number of surrogate preflossing epochs, Ef is the number of flossing epochs during training. Empirically, we
find that the necessary number of surrogate preflossing epochs Ep and flossing episodes Ef is much smaller than both the total number oftraining epochs E. Moreover, Tp can be smaller than T .

us+1 = Ds us(22)where the N × N Jacobian Ds (hs ) = df (hs )/dh characterizes the evolution of generic infinitesimal perturbations during
one step. Note that this Jacobian along the trajectory is equivalent to a stability matrix only at a fixed point, i.e., whenhs+1 = f (hs ) = hs .
We are interested in the asymptotic behavior, and therefore we study the long-term JacobianTt (h0 ) = Dt−1 (ht−1 ) . . . D1 (h1 )D0 (h0 ).

(23)Note that Tt (h0 ) is a product of generally noncommuting matrices. The Lyapunov exponents λ1 ≥ λ2 · · · ≥ λN are definedas the logarithms of the eigenvalues of the Oseledets matrix1

Λ(h0 ) = lim [Tt (h0 )⊤ Tt (h0 )] 2t ,
t→∞

(24)

where ⊤ denotes the transpose operation. The expression inside the brackets is the Gram matrix of the long-term JacobianTt (h0 ). Geometrically, the determinant of the Gram matrix is the squared volume of the parallelotope spanned by thecolumns of Tt (h0 ). Thus, the exponential volume growth rate is given by the sum of the logarithms of its first k (sorted)eigenvalues. Oseledets’ multiplicative ergodic theorem guarantees the existence of the Oseledets matrix Λ(h0 ) for almost
all initial conditions h0 [82]. In ergodic systems, the Lyapunov exponents λi do not depend on the initial condition h0 .
However, for a numerical calculation of the Lyapunov spectrum, Eq 24 cannot be used directly because the long-term
Jacobian Tt (h0 ) quickly becomes ill-conditioned, i.e., the ratio between its largest and smallest singular value diverges
exponentially with time.

L. Convergence of Surrogate Lyapunov Exponents of discretized Spiking Neural Network
In Fig 7, we demonstrate the convergence of the surrogate Lyapunov exponents. We show ing Fig 7A the estimate of thesurrogate Lyapunov exponents λSi for i = 1, 20, 60, 80 for different random seeds for both the network weights, the input
weights, the random input, the orthonormal system and the initial conditions of the network state. B Shows the same, but for
one fixed network weight realization, where the seed of input, the initial conditions, and the orthonormal system are varied.
We note that with fixed network parameters, the different estimates of the surrogate Lyapunov exponent converge tighter.
16Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

A

B
0

0

100
(1/s)
i

i

(1/s)

100

200

200

300

10

300

3

10

2

10

1

10

0

10

t (s)

3

10

2

10

1

10

0

t (s)

Figure 7. Convergence of Lyapunov exponents for g = 1 A Convergence of selected Lyapunov exponents λi for ten identical network
realizations with different initial conditions with simulation time (i = 1, 26, 52, 77, 103, 128) for g = 1. B Same as A but for fixed
network parameters. (Other parameters: N = 64, tsim = 1.5 seconds, tONS = 1).

17

Analyzing and Improving Surrogate Gradient Training in Binary Neural Networks Using Dynamical Systems Theory

B
80

100

100

120

120

140

140

i

(1/s)

80

i

(1/s)

A

160

180

180

200

200

220

220

10

160

3

10

2

10

1

10

0

10

t (s)

3

10

2

10

1

10

0

t (s)

Figure 8. Convergence of Lyapunov exponents for g = 20 A Convergence of selected Lyapunov exponents λi for ten identical network
realizations with different initial conditions with simulation time (i = 1, 26, 52, 77, 103, 128) for g = 1. B Same as A but for fixed
network parameters. (Other parameters: N = 64, tsim = 1.5 seconds, tONS = 20).

18