SPIDE: A PURELY SPIKE-BASED METHOD FOR TRAINING FEEDBACK SPIKING NEURAL NETWORKS
Year: 2021

Abstract:
Spiking neural networks (SNNs) with event-based computation are promising brain-inspired models for energy-efficient applications on neuromorphic hardware. However, most high-performance supervised SNN training methods in machine learning research, such as conversion from artificial neural networks or direct training with surrogate gradients, require complex computation not supported by spiking neurons, which hinders them from spike-based energy-efficient training. Among them, the recently proposed method, implicit differentiation on the equilibrium state (IDE), for training feedback SNNs is a promising way possible for generalization to spike-based learning with flexible network structures. In this paper, we study spike-based implicit differentiation on the equilibrium state (SPIDE) that extends IDE for supervised local learning with spikes, which could be possible for energy-efficient training on neuromorphic hardware. Specifically, we first introduce ternary spiking neuron couples to realize ternary outputs with the common neuron model, and prove that implicit differentiation can be solved by spikes based on this design. With this approach, the whole training procedure can be made as event-driven spike computation and weights are updated locally with two-stage average firing rates. Then to reduce the approximation error of spikes due to the finite simulation time steps, we propose to modify the resting membrane potential. Based on it, the average firing rate, when viewed as a stochastic estimator, achieves an unbiased estimation of iterative solution for implicit differentiation and the variance of this estimator is reduced. With these key components, we can train SNNs with either feedback or feedforward structures in a small number of time steps. Further, the firing sparsity during training demonstrates the great potential for energy efficiency. Meanwhile, even with these constraints, our trained models can still achieve competitive results on MNIST, CIFAR-10, CIFAR-100, and CIFAR10-DVS. Our proposed method demonstrates the great potential for energy-efficient training of SNNs on neuromorphic hardware.

heading: 1 INTRODUCTION
text: Spiking neural networks (SNNs) are brain-inspired models that transmit spikes between neurons for event-driven energy-efficient computation. SNNs can be implemented with less energy on neuromorphic hardware (Akopyan et al., 2015; Davies et al., 2018; Pei et al., 2019; Roy et al., 2019), which can remedy the defects of large energy consumption of artificial neural networks (ANNs).
Different from ANNs, however, directly supervised training of SNNs is a hard problem due to the complex spiking neuron model which is discontinuous. To solve this problem, converting ANNs to SNNs (Hunsberger & Eliasmith, 2015; Rueckauer et al., 2017; Sengupta et al., 2019; Rathi et al., 2019; Deng & Gu, 2021; Yan et al., 2021), or many other direct SNN training methods (Wu et al., 2018; Bellec et al., 2018; Jin et al., 2018; Shrestha & Orchard, 2018; Wu et al., 2019; Neftci et al., 2019; Zhang & Li, 2019; Kim et al., 2020; Zheng et al., 2021; Bohte et al., 2002; Zhang & Li, 2020; Kim et al., 2020; Xiao et al., 2021) have been proposed. While these methods could partly tackle the problems of unsatisfactory performance or high latency, they require complex computation for gradient calculation or approximation, which cannot be implemented by common spiking neurons on neuromorphic hardware. They aim at training SNNs on commonly used computational units, e.g. GPU, and deploying trained models for energy-efficient inference. However they do not consider if the training procedure could leverage the same spike-based computation for gradient calculation and training to reduce the large energy consumption during training as well.
A few previous works try to train SNNs with spikes (Guerguiev et al., 2017; Neftci et al., 2017; Samadi et al., 2017; O’Connor & Welling, 2016; Thiele et al., 2019b;a). They either are based on direct feedback alignment (DFA) (Nøkland, 2016) and performs poorly, or require impractical special neuron models (Thiele et al., 2019b;a). Besides, they only focus on feedforward network structures imitated from ANNs, which ignores feedback connections that are ubiquitous in the human brain and enable neural networks to be shallower and more efficient (Kubilius et al., 2019; Xiao et al., 2021). Actually, feedback structures suit SNNs more since SNNs will naturally compute with multiple time steps, which could reuse representations and avoid uneconomical costs to unfold along time that ANNs suffer from (Xiao et al., 2021). So training algorithms for feedback SNNs, which may also degrade to feedforward structures by taking feedback as zero, is worth more exploration.
An ideal SNN training method should tackle the common problems, be suitable for flexible structures (feedforward or feedback) and be spike-based with high neuromorphic plausibility. The implicit differentiation on the equilibrium state (IDE) method (Xiao et al., 2021), which is recently proposed to train feedback spiking neural networks (FSNNs), is a promising method that may generalize to spike-based learning for requirement. They derive that the forward computation of FSNNs converges to an equilibrium state, which follows a fixed-point equation. Based on it, they propose to train FSNNs by implicit differentiation on this equation, which tackles the common difficulties for SNN training including non-differentiability and large memory costs, and has interesting local update properties. In their method, however, they leverage general root-finding methods to solve implicit differentiation, which requires complex computation on standard computation systems.
In this work, we extend the IDE method to spike-based IDE (SPIDE), which fulfills our requirements and has great potential for energy-efficient training of SNNs on neuromorphic hardware, by introducing ternary spiking neuron couples and proposing to solve implicit differentiation by spikes based on them. Our method is also applicable to feedforward structures by degrading the feedback connection as zero. In practice, however, it may require long time steps to stabilize the training with spikes due to approximation error for gradients. So we further dive into the approximation error from the statistical perspective, and propose to simply adjust the resting potential of SNNs to achieve an unbiased estimation of gradients and reduce the estimation variance of SNN computation. With these methods, we can train our models in a small number of time steps, which could further improve the energy efficiency as well as the latency. Our contributions include:
1. We propose the SPIDE method that is the first to train high-performance SNNs by spikes with common neuron models. Specifically, we propose ternary spiking neuron couples and prove that implicit differentiation for gradient calculation can be solved by spikes based on this design. Our method is applicable to both feedback and feedforward structures.
2. We theoretically analyze the approximation error of solving implicit differentiation by spikes, and propose to modify the resting potential to remove the approximation bias and reduce the estimation variance, which enables training in a small number of time steps.
3. Experiments show the low latency and firing sparsity during training, which demonstrates the great potential for energy-efficient training of SNNs on neuromorphic hardware. The performance on MNIST, CIFAR-10, CIFAR-100 and CIFAR10-DVS are also competitive.

heading: 2 RELATED WORK
text: Early works seek biologically inspired methods to train SNNs, e.g. spike-time dependent plasticity (STDP) (Diehl & Cook, 2015) or reward-modulated STDP (Legenstein et al., 2008). Since the rise of successful ANNs, several works try to convert trained ANNs to SNNs to obtain high performance (Hunsberger & Eliasmith, 2015; Rueckauer et al., 2017; Sengupta et al., 2019; Rathi et al., 2019; Deng & Gu, 2021; Yan et al., 2021). However, they suffer from extremely large time steps
and their structures are limited in the scope of ANNs. Others try to directly train SNNs by imitating backpropagation throught time (BPTT) and use surrogate derivative for discontinuous spiking functions (Lee et al., 2016; Wu et al., 2018; Bellec et al., 2018; Jin et al., 2018; Shrestha & Orchard, 2018; Wu et al., 2019; Zhang & Li, 2019; Neftci et al., 2019; Zheng et al., 2021) or compute gradient with respect to spiking times (Bohte et al., 2002; Zhang & Li, 2020; Kim et al., 2020). However, they suffer from approximation error and large memory costs. Xiao et al. (2021) propose the IDE method to train feedback spiking neural networks, which decouples the forward and backward procedures and avoids the common SNN training problems. However, all these methods require complex computation during training rather than spike-based. A few works focusing on training SNNs with spikes either are based on feedback alignment and limited in simple datasets (Guerguiev et al., 2017; Neftci et al., 2017; Samadi et al., 2017; O’Connor & Welling, 2016), or require impractical special neuron models that require consideration of accumulated spikes for spike generation (Thiele et al., 2019b;a), which is impractical on neuromorphic hardware. And they are only applicable to feedforward architectures. Instead, we are the first to leverage spikes with common neuron models to train SNNs with feedback or feedforward structures. The comparison of different methods is illustrated in Table 1.

heading: 3 PRELIMINARIES
text: We first introduce preliminaries about spiking neurons and the IDE training method. The basic thought of IDE (Xiao et al., 2021) is to identify the underlying equilibrium states of FSNN computation so that gradients can be calculated based on implicit differentiation on the equilibrium state. We will briefly introduce the conclusion of equilibrium states in Section 3.2 and the IDE method in Section 3.3. For more descriptions about the background please refer to Appendix A.

heading: 3.1 SPIKING NEURAL NETWORK MODELS
text: Spiking neurons draw inspirations from the human brain to communicate with each other by spikes. Each neuron integrates information from input spike trains by maintaining a membrane potential through a differential equation, and generates an output spike once the membrane potential exceeds a threshold, following which the membrane potential is reset to the resting potential. We consider the commonly used integrate and fire (IF) model and simple current model, whose discretized computational form is:  ui [t+ 0.5] = ui[t] + ∑ j wijsj [t] + b,
si[t+ 1] = H(ui [t+ 0.5]− Vth), ui[t+ 1] = ui [t+ 0.5]− (Vth − urest)si[t+ 1],
(1)
where ui[t] is the membrane potential of neuron i at time step t, si[t] is the binary output spike train of neuron i, wij is the connection weight from neuron j to neuron i, b is bias,H is the Heaviside step function, Vth is the firing threshold, and urest is the resting potential. We use subtraction as the reset operation. urest is usually taken as 0 in previous work, while we will reconsider it in Section 4.3.

heading: 3.2 EQUILIBRIUM STATES OF FEEDBACK SPIKING NEURAL NETWORKS
text: Xiao et al. (2021) derive that the (weighted) average rate of spikes during FSNN computation with common neuron models would converge to an equilibrium state following a fixed-point equation given convergent inputs. We focus on the conclusions with the discrete IF model under both singlelayer and multi-layer feedback structures. The single-layer structure has one hidden layer of neurons with feedback connections on this layer. The update equation of membrane potentials is:
u[t+ 1] = u[t] +Ws[t] + Fx[t] + b− (Vth − urest)s[t+ 1], (2) where u[t] and s[t] are the vectors of membrane potentials and spikes of these neurons, x[t] is the input at time step t, W is the feedback weight matrix, and F is the weight matrix from inputs to these neurons. The average input and average firing rate are defined as x[t] = 1t+1 ∑t τ=0 x[τ ] and
α[t] = 1t ∑t τ=1 s[τ ], respectively. Define σ(x) = min(1,max(0, x)).
The equilibrium state of the single-layer FSNN is described as (Xiao et al., 2021): If the average inputs converge to an equilibrium point x[t]→ x∗, and there exists γ < 1 such that ‖W‖2 ≤ γVth, then the average firing rates of FSNN with discrete IF model will converge to an equilibrium point α[t] → α∗, which satisfies the fixed-point equation α∗ = σ ( 1 Vth (Wα∗ + Fx∗ + b) )
. Note that they take urest = 0 in this conclusion, if we consider nonzero urest, the constraint and the fixedpoint equation should be ‖W‖2 ≤ γ(Vth − urest) and α∗ = σ ( 1
Vth−urest (Wα ∗ + Fx∗ + b)
) .
The multi-layer structure incorporates more non-linearity into the equilibrium fixed-point equation, which has multiple layers with feedback connections from the last layer to the first layer. The update equations of membrane potentials are expressed as:{
u1[t+ 1] = u1[t] +W1sN [t] + F1x[t] + b1 − (Vth − urest)s1[t+ 1], ul[t+ 1] = ul[t] + Flsl−1[t+ 1] + bl − (Vth − urest)sl[t+ 1], l = 2, · · · , N.
(3)
The equilibrium state of the multi-layer FSNN with urest is described as (Xiao et al., 2021): If the average inputs converge to an equilibrium point x[t] → x∗, and there exists γ < 1 such that ‖W1‖2‖FN‖2 · · · ‖F2‖2 ≤ γ(Vth − urest)N , then the average firing rates of multi-layer FSNN with discrete IF model will converge to equilibrium points αl[t] → αl∗, which satisfy the fixedpoint equations α1∗ = f1 ( fN ◦ · · · ◦ f2(α1 ∗ ),x∗ ) and αl+1∗ = fl+1(αl ∗ ), where f1(α,x) =
σ (
1 Vth−urest (W
1α+ F1x+ b1) ) and fl(α) = σ (
1 Vth−urest (F
lα+ bl) ) .

heading: 3.3 IDE TRAINING METHOD
text: Based on the equilibrium states in Section 3.2, we can train FSNNs by calculating gradients with implicit differentiation (Xiao et al., 2021). Let α = fθ(α) denote the fixed-point equation of the equilibrium state which is parameterized by θ, gθ(α) = fθ(α) − α, and let L(α∗) denote the objective function with respect to the equilibrium state α∗. The implicit differentiation satisfies( I − ∂fθ(α
∗) ∂α∗
) dα∗
dθ = ∂fθ(α ∗) ∂θ (Bai et al., 2019) (we follow the numerator layout convention for
derivatives). Therefore, the differentiation of L(α∗) for parameters can be calculated as: ∂L(α∗) ∂θ = −∂L(α ∗) ∂α∗ ( J−1gθ |α∗ ) ∂fθ(α∗) ∂θ , (4) where J−1gθ |α∗ is the inverse Jacobian of gθ evaluated at α ∗. The calculation of inverse Jacobian can be avoided by solving an alternative linear system (Bai et al., 2019; 2020; Xiao et al., 2021):( J>gθ |α∗ ) β + ( ∂L(α∗) ∂α∗ )> = 0. (5)
Note that a readout layer after the last layer of neurons will be constructed for output (Xiao et al., 2021), which is equivalent to a linear transformation on the approximate equilibrium state, i.e. the output would be o = Woα∗, and the loss will be calculated between o and labels y with a common criterion such as cross-entropy. Then the gradient on the equilibrium state could be calculated. For the solution of implicit differentiation, Xiao et al. (2021) follow Bai et al. (2019; 2020) to leverage root-finding methods, while we will solve it by spike dynamics, as will be derived in Section 4. We treat the forward computation of average firing rates α[T ] of FSNNs at time step T roughly reach the equilibrium state. Then by substituting α∗ by α[T ] in the above equations, gradients for the parameters can be calculated only with α[T ] and the equation, and we calculate them based on spikes. With the gradients, first-order optimization methods such as SGD (Rumelhart et al., 1986) and its variants can be applied to update parameters.

heading: 4 SPIKE-BASED IMPLICIT DIFFERENTIATION ON THE EQUILIBRIUM STATE
text: In this section, we present our SPIDE method that calculates the whole training procedure based on spikes. We first introduce ternary spiking neuron couples in Section 4.1 and how to solve implicit differentiation in Section 4.2. Then we theoretically analyze the approximation error and propose the improvement in Section 4.3. Finally, a summary of the training pipeline is presented in Section 4.4.

heading: 4.1 TERNARY SPIKING NEURON COUPLES
text: The common spiking neuron model only generates spikes when the membrane potential exceeds a positive threshold, which limits the firing rate from representing negative information. To enable approximation of possible negative values for implicit differentiation calculation in Section 4.2, we require negative spikes, whose expression could be:
si[t+ 1] = T (ui[t+ 0.5], Vth) =  1, ui[t+ 0.5] > Vth
0, |ui[t+ 0.5]| ≤ Vth − 1, ui[t+ 0.5] < −Vth , (6)
and the reset is the same as usual: ui[t+1] = ui[t+0.5]−(Vth−urest)si[t+1]. Direct realization of such ternary output, however, may be not supported by common neuromorphic hardware for SNNs.
We propose to leverage two coupled common neurons to realize this computation. As illustrated in Figure 1, the two coupled neurons with the common IF model (Eq. (1)) receive opposite inputs and output opposite spikes, which aim to deal with positive information and negative information with spikes, respectively. They should share a reset operation in order to accord with Eq. (6), which can be realized by the connection between them: as we use subtraction as the reset operation, the connection whose weight equals Vth − urest enables one neuron to reset another equivalently. To see how this works, consider the condition that the accumulated membrane potential of neuron 1 reaches Vth, then neuron 1 would generate a spike and reset, and the output is this positive spike. At the same time, the membrane potential of neuron 2 is −Vth and the neuron will not fire and reset, but the spike from neuron 1 will reset it to −urest, which accords with our desired reset for ternary output. Similarly, if the inputs are negative,
neuron 2 will generate a spike which will be treated negative as output, and both neurons are reset. For the operation of taking negative, one solution is to enable the reverse operation on hardware, another is to reconnect neuron 2 with other neurons while taking the weight negative to that of neuron 1. Therefore, such kind of coupled neurons with the common IF model could realize ternary output.
We note that the SpikeGrad algorithm (Thiele et al., 2019a) also requires neurons for ternary output. However, they do not consider how such kind of operation can be implemented with the common neuron model on neuromorphic hardware, and moreover, they propose another modified ternary model in practice that requires consideration of accumulated spikes for spike generation, which is further impractical on neuromorphic hardware. Differently, our method can be realized with the common neuron model suitable for neuromorphic hardware.

heading: 4.2 SOLVING IMPLICIT DIFFERENTIATION WITH SPIKES
text: Based on the coupled neurons in Section 4.1, we can solve implicit differentiation with spikes. For notation simplicity, we directly use Eq. (6) as a ternary neuron without detailing coupled neurons below. Our main focus is on solving Eq. (5) with spikes. The brief outline for the derivation is: we first derive the update equation of membrane potentials in SNN computation, then we derive the equivalent equation of the rate of spikes with eliminating perturbation, finally, we could prove that the rate of spikes converges to the solution of Eq. (5).
We first consider the single-layer condition. Let α[TF ] denote the average firing rate of these neurons after the forward computation with time steps TF as an approximate equilibrium state (we treat
the forward procedure as the first stage), g = (
∂L ∂α[TF ]
)> denote the gradient of the loss function
on this approximate equilibrium state, and m = σ′(α[TF ]),M = Diag(m) denote a mask indicator based on the firing condition in the first stage, where σ′(x) = { 1, 0 < x < 1
0, else . We will have another
TB time steps in the second backward stage to calculate implicit differentiation. We set the input to these neurons as g at all time steps, which can be viewed as input currents (Zhang & Li, 2020; Xiao et al., 2021). Then along the inverse connections of neurons and with a mask on neurons or weights and an output rescale, the computation of FSNN with ternary neurons is calculated as:
u[t+ 1] = u[t] + 1
Vth − urest (MW)>s[t] + g − (V bth − ubrest)s[t+ 1], (7)
where Vth, urest and V bth, u b rest are the threshold and resting potential during the first and second stage, respectively. Define the ‘average firing rate’ at this second stage as β[t] = 1t ∑t τ=1 s[τ ], and u[0] = 0, s[0] = 0, then through summation, we have:
β[t+ 1] = 1
V bth − ubrest
( t
t+ 1
1 Vth − urest (MW)>β[t] + g − u[t+ 1] t+ 1
) . (8)
Since there could be at most t spikes during t time steps, β would be bounded in the range of [−1, 1]. The membrane potential ui[t] will maintain the exceeded terms, i.e. define vi[t] = ( t t+1 1 Vth−urest (MW) >β[t] + g ) i , we can divide ui[t] as uEi [t] + u B i [t], where
uEi [t] = max ( vi[t]− V bth, 0 ) +min ( vi[t] + V b th, 0 ) is the exceeded term while uBi [t] is a bounded term (Xiao et al., 2021) which is typically bounded in the range of [−V bth, V bth]. Then, Eq. (8) turns into:
β[t+ 1] = φ
( 1
V bth − ubrest
( t
t+ 1
1 Vth − urest (MW)>β[t] + g − u
B [t+ 1]
t+ 1
)) , (9)
where φ(x) = min(1,max(−1, x)). Note that if the input g and weight (MW)> are in an appropriate range, there would be no exceeded term and therefore φ will not take effect. Indeed we will rescale the loss to control g in an appropriate range, as will be indicated in Section 4.3. With this consideration, we could derive that the average firing rate β[t] converges to the solution of Eq. (5). Theorem 1. If there exists γ < 1 such that ‖(MW)>‖2 ≤ γ(Vth − urest)(V bth − ubrest), then the average firing rate β[t] will converge to an equilibrium point β[t] → β∗. When V bth − ubrest = 1, and there exists λ < 1 such that ‖(MW)>‖∞ ≤ λ(Vth− urest) and ‖g‖∞ ≤ 1−λ, then β∗ is the solution of Eq. (5).
The proof and discussion of assumptions are in Appendix B. With Theorem 1, we can solve Eq. (5) by simulating this second stage of SNN computation to obtain the ‘firing rate’ β[TB ] as the approximate solution. Plugging this solution to Eq. (4), the gradients can be calculated by: ∇WL =
1 Vth−urestMβ[TB ]α[TF ] >,∇FL = 1Vth−urestMβ[TB ]x[TF ] >,∇bL = 1Vth−urestMβ[TB ].
Note that in practice, even if the data distribution is not properly in the range of φ, we can still view φ as a kind of clipping for improperly large numbers, which could be similar to empirical techniques like “gradient clipping” to stabilize the training.
Then we consider the extension to the multi-layer condition. Let αl[TF ], l = 1, 2, · · · , N denote the average firing rate of neurons in layer l after the forward computation, g = (
∂L ∂αN [TF ]
)> denote
the gradient of the loss function on the approximate equilibrium state of the last layer, and ml = σ′(αl[TF ]),M
l = Diag(ml) denote the mask indicators. Similarly, we will have another TB time steps in the second stage to calculate implicit differentiation. We set the input to the last layer as g at all time steps. Then along the inverse connections of neurons and with a mask on neurons or weights and an output rescale, the computation of FSNN with ternary neurons is calculated as: uN [t+ 1] = uN [t] + 1 Vth − urest (M1W1)>s1[t] + g − (V bth − ubrest)sN [t+ 1], ul[t+ 1] = ul[t] + 1
Vth − urest (Ml+1Fl+1)>sl+1[t+ 1]− (V bth − ubrest)sl[t+ 1], l = N − 1, · · · , 1.
(10) The ‘average firing rates’ βl[t] are similarly defined for each layer, and the equivalent form can be similarly derived as: βN [t+ 1] = φ ( 1 V bth − ubrest ( t t+ 1 1 Vth − urest (M1W1)>β1[t] + g − u NB [t+ 1] t+ 1 )) , βl[t+ 1] = φ ( 1
V bth − ubrest
( 1
Vth − urest (Ml+1Fl+1)>βl+1[t+ 1]− u
lB [t+ 1]
t+ 1
)) .
(11) The convergence of the ‘firing rate’ at the last layer to the solution of Eq. (5) can be similarly derived as Theorem 1. However, we need to calculate gradients for each parameter as Eq. (4), which is more complex than the single layer condition. Actually, we can derive that the ‘firing rates’ at each layer converge to equilibrium points, based on which the gradients can be easily calculated with information from the adjacent layers. Theorem 2 gives a formal description. Theorem 2. If there exists γ < 1 such that ‖(M1W1)>‖2‖(MNFN )>‖2 · · · ‖(M2F2)>‖2 ≤ γ(Vth − urest)N (V bth − ubrest) N , then the average firing rates βl[t] will converge to equilibrium points βl[t] → βl∗. When V bth − ubrest = 1, and there exists λ < 1 such that ‖(M1W1)>‖∞ ≤ λ(Vth − urest), ‖(MlFl)>‖∞ ≤ λ(Vth − urest), l = 2, · · · , N and ‖g‖∞ ≤ 1 − λN , then βN ∗
is the solution of Eq. (5), and βl ∗ = ( ∂hN (α
N∗) ∂hl(αN ∗) )> βN ∗ , l = N − 1, · · · , 1, where hl(αN ∗ ) =
fl ◦ · · · ◦ f2 ( f1(α N ∗,x∗) ) , l = N, · · · , 1.
The functions fl are defined in Section 3.2. For the proof please refer to Appendix C. With Theorem 2, by plugging the solutions explicitly into Eq. (4), the gradients can be calculated by ∇FlL =
1 Vth−urestM lβl[TB ]α l−1[TF ] >, l = 2, · · · , N,∇F1L = 1Vth−urestM 1β1[TB ]x[TF ] >,∇W1L = 1 Vth−urestM lβl[TB ]α N [TF ] >,∇blL = 1Vth−urestM lβl[TB ].
Note that the gradient calculation shares an interesting local property, i.e. it is proportional to the firing rates of the two neurons connected by it: ∇Fli,jL = 1 Vth−urestm l iβ l iα l−1 j . During calculation, since we will have the firing rate of the first stage before the second stage, this calculation can also be carried out by event-based calculation triggered by the spikes in the second stage. So it would be plausible on neuromorphic hardware as well.
Also, note that the theorems still hold if we degrade our feedback models to feedforward ones by setting feedback connections as zero. In this setting, the dynamics and equilibriums degrade to direct functional mappings, and the implicit differentiation degrades to the explicit gradient. We can still approximate gradients with this computation.
In the following, we take V bth−ubrest = 1 by default to fulfill the assumption of theorems (it may take other values if we correspondingly rescale the outputs and we set 1 for simplicity). Other techniques like dropout can also be integrated into the calculation. Please refer to Appendix D for details.

heading: 4.3 REDUCING APPROXIMATION ERROR
text: Section 4.2 derives that we can solve implicit differentiation with spikes, as the average firing rate will gradually converge to the solution. In practice, however, we will simulate SNNs for finite time steps, and a smaller number of time steps is better for lower energy consumption. This will introduce approximation error which may hamper training. In this subsection, we theoretically study the approximation error and propose to adjust the resting potential to reduce it. Inspired by the theoretical analysis on quantized gradients (Chen et al., 2020), we will analyze the error from the statistical perspective.
For the ‘average firing rates’ βl[t] in Eq. (8) and the multilayer counterparts, the approximation error e to the equilibrium states consists of three independent parts ee, er and ei: the first is ul E [t + 1] that is the exceeded term due to the limitation of spike number, the second is ulB [t + 1] which can be viewed as a bounded random variable, and the third is the convergence error of the iterative update scheme without ul[t + 1], i.e. let bl[t] denote the iterative sequences for solving βl∗ as bl[t+1] = tt+1 1 Vth−urest (M
l+1Fl+1)>bl[t], the convergence error is ‖bl[t]−βl∗‖. The second part er can be again decomposed into two independent components er = eq + es: eq is the quantization effect due to the precision of firing rates ( 1T for T time steps) if we first assume the same average inputs at all time steps, and es is due to the random arrival of spikes rather than the average condition, as there might be unexpected output spikes, e.g. the average input is 0 and the expected output should be 0, but two large positive inputs followed by one larger negative input at the last time would generate two positive spikes while only one negative spike. So the error is divided into: e = ee + eq + es + ei. Since the iterative formulation is certain for ei, we focus on ee, eq and es.
Firstly, the error eq due to the quantization effect is influenced by the input scale and time steps TB . To enable proper input scale and smaller time steps, we will rescale the loss function by a factor sl, since the magnitude of gradients considering the direct cross-entropy loss function is relatively small. We scale the loss to an appropriate range so that information can be propagated by SNNs in smaller time steps, and most signals are in the range of φ as analyzed in Section 4.2. The base learning rate will be scaled by 1sl correspondingly. This is also adopted by Thiele et al. (2019a).
Then given the scale and number of time steps, eq , ee and es can be treated as random variables from statistical perspective, and we view βl[t] as stochastic estimators for the equilibrium states with ei. For the stochastic optimization algorithms, the expectation and variance of the gradients are important for convergence and convergence rate (Bottou, 2010), i.e. we hope an unbiased estimation of gradients and smaller estimation variance. As for ee and es, they depends on the input data and the expectations are E[ee] = 0,E[es] = 0 (the positive and negative parts have the same probability). While for eq , it will depend on our hyperparameters V bth and u b rest. Since the remaining terms in ul B [t + 1] caused by the quantization effect is in the range of [ubrest, V b th] for positive terms while [−V bth,−ubrest] for negative ones, given V bth − ubrest and considering the uniform distribution, only
when ubrest = −V bth, the expectation E[eq] = 0 for both positive and negative terms. Therefore, we should adjust the resting potential from commonly used 0 (Wu et al., 2018; Sengupta et al., 2019; Xiao et al., 2021) to −V bth for unbiased estimation, as described in Proposition 1. Proposition 1. For fixed V bth − ubrest and uniformly distributed inputs and eq , only when ubrest = −V bth, βl[t] are unbiased estimators for bl[t].
Also, taking ubrest = −V bth achieves the smallest estimation variance for the quantization effect eq , considering the uniform distribution on [ubrest, V b th]∪ [−V bth,−ubrest]. Since the effects of ee, es and ei are independent of eq and their variance is certain given inputs, it leads to Proposition 2.
Proposition 2. Taking ubrest = −V bth reduces the variance of estimators βl[t].
With this analysis, we will take V bth = 0.5, u b rest = −0.5 in the following to stabilize the training. For Vth and urest during the first forward stage, we will also take urest = −Vth.

heading: 4.4 DETAILS AND TRAINING PIPELINE
text: The original IDE method (Xiao et al., 2021) leverages other training techniques including modified batch normalization (BN) and restriction on weight spectral norm. Since the batch statistical information might be hard to obtain for calculation on neuromorphic hardware and we seek for algorithms that could be possible on it, we drop the BN component in our SPIDE method. The restriction on the weight norm, however, is necessary for the convergence of feedback models, as indicated in theorems. We will adjust it for a more friendly calculation, please refer to Appendix D for details.
We summarize our training pipeline as follows (we also provide detailed pseudocodes in Appendix E). There are two stages for forward and backward procedures respectively. In the first stage, SNNs will receive inputs and perform the calculation as Eq. (1,2,3) for TF time steps, after which we get the output from the readout layer, and save the average inputs as well as the average firing rates and masks of each layer for the second stage. In the second stage, the last layer of SNNs will receive gradients for outputs and perform calculation along the inverse connections as Eq. (6,7,10) for TB time steps, after which we get the ‘average firing rates’ of each layer. Based on the firing rates from two stages, the gradients for parameters can be calculated as in Section 4.2 and then the first-order optimization algorithm is applied to update the parameters.

heading: 5 EXPERIMENTS
text: In this section, we conduct experiments to demonstrate the effectiveness of our method and the great potential for energy-efficient training. We simulate the computation on common computational units. Please refer to Appendix D for implementation details and descriptions.
We first evaluate the effectiveness of our method in a small number of time steps. As shown in Table 2, we can train high-performance models with low latency (TF = 30) in a very small number time steps during training (e.g. TB = 50), which indicates the low latency and high energy efficiency. Note that the ANN-SNN methods usually require hundreds to thousands of time steps just for satisfactory inference performance, and direct training methods show that relatively small time steps are enough for inference, while we are the first to demonstrate that even training of SNNs can be carried out with spikes in a very small number of time steps. This is due to our analysis and improvement to reduce the approximation error, as illustrated in the ablation study in Appendix F.1.
Then we analyze the firing rate statistics to demonstrate the potential of energy efficiency. Since the energy consumption on event-driven neuromorphic hardware is proportional to the number of spikes, we present the average firing rates for forward and backward stages (for backward, both positive and negative spikes are considered as firing) in Figure 2. It shows the firing sparsity of our method, and spikes are sparser in the backward stage with around only 3%. Combined with the small number of time steps, this demonstrates the great potential for the energy-efficient training of SNNs based on our method on neuromorphic hardware.
Finally we evaluate the performance of our method on MNIST (LeCun et al., 1998), CIFAR-10 and CIFAR-100 (Krizhevsky & Hinton, 2009). We compare our method to several ANN-SNN methods (Hunsberger & Eliasmith, 2015; Sengupta et al., 2019; Deng & Gu, 2021), direct SNN training methods (Wu et al., 2018; Xiao et al., 2021), and SpikeGrad (Thiele et al., 2019a). As shown in Table 3, we can train models with a small number of time steps and our trained models achieve
competitive results on MNIST and CIFAR-10. Compared with the original IDE method (Xiao et al., 2021), since we discard the BN component, our generalization performance is poorer (a detailed discussion is in Appendix F.2). Compared with SpikeGrad (Thiele et al., 2019a), we can use fewer neurons and parameters due to flexible network structure choices, and a small number of time steps while they do not report this important feature. Besides, we use common neuron models while they require impractical models, as indicated in Section 4.1. The results and discussion on CIFAR-100 and CIFAR10-DVS are in Appendix F.3 and F.4 due to the space limit, and our model could achieve 64.07% and 60.7% accuracy respectively. The result on CIFAR-100 is competitive for networks without BN, though it is poorer than IDE with BN. And the result on CIFAR10-DVS is competitive among results of common SNN models. It shows the effectiveness of our method even with constraints of purely spike-based training. Future work could seek normalization techniques friendly for neuromorphic computation and our desired algorithm to further improve the performance.

heading: MNIST

heading: CIFAR-10

heading: 6 CONCLUSION
text: In this work, we propose the SPIDE method that generalize the IDE method to enable the whole training of SNNs with either feedback or degraded feedforward structures to be based on spikes with common neuron models. We prove that the implicit differentiation can be solved with spikes by our coupled neurons. We also analyze the approximation error due to finite time steps, and propose to adjust the resting potential of SNNs. Experiments show that we could achieve competitive performance with a small number of training time steps and sparse spikes, which demonstrates the great potential of our method for energy-efficient training of SNNs on neuromorphic hardware.

heading: A MORE BACKGROUND ABOUT THE IDE TRAINING METHOD
text: Due to the complex spiking neuron model which is discontinuous, directly supervised training of SNNs is a hard problem, since the explicit computation is non-differentiable and therefore backpropagation along the forward computational graph could be problematic. The IDE training method (Xiao et al., 2021) considers another approach to calculating gradients that does not rely on the exact reverse of the forward computation, which avoids the problem of non-differentiability as well as large memory costs by BPTT-like methods with surrogate gradients. Specifically, the IDE training method first derives that the (weighted) average firing rate of FSNN computation with common neuron models would gradually evolve to an equilibrium state along time, which follows a fixed-point equation. Then by viewing the forward computation of FSNN as a black-box solver for this equation, and applying implicit differentiation on the equation, gradients can be calculated only based on this equation and the (weighted) average firing rate during forward computation rather than the exact forward procedure. Therefore, the forward and backward procedures are decoupled and the non-differentiability is avoided.
As briefly introduced in Section 3.2, the IDE method defines the average firing rates of spikes during forward computation, i.e. α[t]. Then IDE could derive an equivalent update equation for α[t], based on the integrated update equations of membrane potentials Eq. (2,3). With the equivalent equation for α[t], IDE proves that under certain conditions, α[t] converges to an equilibrium state following a fixed-point equation. With the assumption that α[T ] after simulation of T time steps roughly follows the fixed-point equation of the equilibrium state, gradients of the loss function for parameters can be calculated by implicit differentiation, as introduced in Section 3.3. So the training pipeline for the IDE method can be summarized as: first simulate FSNN computation for T time steps to obtain the rate of spikes α[T ], then solve the implicit differentiation by root-finding methods for gradient calculation based on α[T ] and the derived fixed-point equation of equilibrium states, finally apply gradient-based optimizers to update parameters.
Contribution of this work compared with IDE. Compared with IDE, this work extends the thought of equilibrium of spikes to solving implicit differentiation, which enables the whole training procedure to be based on spike computation with common neuron models and provides the potential for energy-efficient training of SNNs on neuromorphic hardware. IDE only derives the equilibrium states of forward FSNN computation and requires general root-finding methods with complex computation to solve implicit differentiation for training. It remains unclear if it could be solved by common spiking neuron models with only positive firing. This work design ternary spiking neuron couples and prove that the equilibrium of spike computation can be leveraged to solve implicit differentiation based on the design, and we propose to modify the resting membrane potential to make it practical in a relatively small number of time steps. This enables the proposed SPIDE method to be the first to train high-performance SNNs with low latency and firing sparsity by spikes with common neuron models, demonstrating the great potential for energy-efficient training of high-performance SNNs on neuromorphic hardware.

heading: B PROOF OF THEOREM 1
text: Proof. We first prove the convergence of β[t]. Let Vu and V bu denote Vth − urest and V bth − ubrest respectively. Consider ‖β[t+ 1]− β[t]‖, it satisfies:
‖β[t+ 1]− β[t]‖
= ∥∥∥∥φ( 1V bu ( t t+ 1 1 Vu (MW)>β[t] + g − u B [t+ 1] t+ 1 )) −φ ( 1
V bu
( t− 1 t 1 Vu (MW)>β[t− 1] + g − u B [t] t ))∥∥∥∥ ≤
∥∥∥∥φ( 1V bu ( 1 Vu (MW)>β[t] + g )) − φ ( 1 V bu ( 1 Vu (MW)>β[t− 1] + g ))∥∥∥∥ +
∥∥∥∥φ( 1V bu ( t t+ 1 1 Vu (MW)>β[t] + g − u B [t+ 1] t+ 1 )) − φ ( 1 V bu ( 1 Vu (MW)>β[t] + g ))∥∥∥∥
+ ∥∥∥∥φ( 1V bu ( t− 1 t 1 Vu (MW)>β[t− 1] + g − u B [t] t )) − φ ( 1 V bu ( 1 Vu (MW)>β[t− 1] + g ))∥∥∥∥ ≤
∥∥∥∥φ( 1V bu ( 1 Vu (MW)>β[t] + g )) − φ ( 1 V bu ( 1 Vu (MW)>β[t− 1] + g ))∥∥∥∥ + 1
V bu (∥∥∥∥ 1t+ 1 1Vu (MW)>β[t] ∥∥∥∥+ ∥∥∥∥uB [t+ 1]t+ 1 ∥∥∥∥+ ∥∥∥∥1t 1Vu (MW)>β[t− 1] ∥∥∥∥+ ∥∥∥∥uB [t]t ∥∥∥∥) . (12)
As ‖(MW)>‖2≤ γVuV bu , γ < 1, and |uBi [t]| is bounded, we have ‖β[t+ 1]‖ ≤ γ ‖β[t]‖ + 1Vu ( ‖g‖+ ∥∥∥uB [t+1]t+1 ∥∥∥) ≤ γ ‖β[t]‖ + c, where c is a constant. Therefore ‖β[t]‖ is bounded. Then ∀ > 0,∃T1 such that when t > T1, we have:
1
V bu (∥∥∥∥ 1t+ 1 1Vu (MW)>β[t] ∥∥∥∥+ ∥∥∥∥uB [t+ 1]t+ 1 ∥∥∥∥+ ∥∥∥∥1t 1Vu (MW)>β[t− 1] ∥∥∥∥+ ∥∥∥∥uB [t]t ∥∥∥∥) ≤ (1− γ)2 . (13)
And since ‖(MW)>‖2≤ γVuV bu , we have:∥∥∥∥φ( 1V bu ( 1 Vu (MW)>β[t] + g )) − φ ( 1 V bu ( 1 Vu (MW)>β[t− 1] + g ))∥∥∥∥ ≤ γ ‖β[t]− β[t− 1]‖ . (14)
Therefore, when t > T1, it holds that:
‖β[t+ 1]− β[t]‖ ≤ γ ‖β[t]− β[t− 1]‖+ (1− γ) 2 . (15)
By iterating the above inequality, we have ‖β[t + 1] − β[t]‖ ≤ γt−T1‖β[T1 + 1] − β[T1]‖ + (1−γ)
2
( 1 + γ + · · ·+ γt−T1−1 ) < γt−T1‖β[T1 +1]−β[T1]‖+ 2 . There exists T2 such that when
t > T1 + T2, γt−T1‖β[T1 + 1] − β[T1]‖ ≤ 2 , and therefore ‖β[t + 1] − β[t]‖ < . According to Cauchy’s convergence test, the sequence {β[t]}∞i=0 converges to β∗. Considering the limit, it satisfies β∗ = φ ( 1 V bu ( 1 Vu (MW)>β∗ + g )) .
When V bth − ubrest = 1, and there exists λ < 1 such that ‖(MW)>‖∞ ≤ λ(Vth − urest) and ‖g‖∞ ≤ 1− λ, the equation turns into β∗ = φ ( 1 Vu (MW)>β∗ + g ) . We have:
‖β∗‖∞ = ∥∥∥∥φ( 1Vu (MW)>β∗ + g )∥∥∥∥ ∞ ≤ ∥∥∥∥ 1Vu (MW)>β∗ ∥∥∥∥ ∞ + ‖g‖∞ ≤ λ ‖β ∗‖∞ + ‖g‖∞ .
(16)
Therefore, ‖β∗‖∞ ≤ ‖g‖∞ 1−λ ≤ 1, and ∥∥∥ 1Vu (MW)>β∗ + g∥∥∥∞ ≤ ∥∥∥ 1Vu (MW)>β∗∥∥∥∞ + ‖g‖∞ ≤ λ+ (1− λ) = 1. It means β∗ = φ ( 1 Vu (MW)>β∗ + g ) = 1Vu (MW) >β∗ + g.
Taking fθ(α∗) = σ (
1 Vth−urest (Wα
∗ + Fx∗ + b) )
(i.e. the fixed-point equation at the equilibrium state as in Section 3.2) explicitly into Eq. (5), the linear equation turns into(
1 Vu
(MW)> − I ) β + g = 0, where Vu,M,g are previously defined. Therefore, β∗ satisfies
this equation. And since ‖(MW)>‖2≤ γVu, γ < 1, the equation has the unique solution β∗.
Remark 1. As for the assumptions in the theorem, firstly, when V bth − ubrest = 1 as we will take, the assumption for the convergence is weaker than that for the convergence in the forward stage (in Section 3.2), because ‖(MW)>‖2 ≤ ‖W‖2 as M is a diagonal mask matrix. We will restrict the
spectral norm of W following Xiao et al. (2021) to encourage the convergence of the forward stage (in Appendix D), then this backward stage would converge as well.
The assumptions for the consistency of the solution is a sufficient condition. In practice, the weight norm will be partially restricted by weight decay and our restriction on Frobenius norm (in Appendix D), as well as the diagonal mask matrix M which would be sparse if the forward firing events are sparse, and we will rescale the loss so that the input g is in an appropriate range, as indicated in Section 4.3. Even if these assumptions are not satisfied, we can view φ as a kind of empirical clipping techniques to stabilize the training, as indicated in Section 4.2. The discussion is similar for the multi-layer condition (Theorem 2) in the next section.

heading: C PROOF OF THEOREM 2
text: Proof. We first prove the convergence of βl[t]. Let Vu and V bu denote Vth − urest and V bth − ubrest respectively. Let gt+1N (β,g,u B) = φ (
1 V bu ( t t+1 1 Vu (M1W1)>β + g − u B t+1 )) ,
gtl (β,u B) = φ ( 1 V bu ( 1 Vu (Ml+1Fl+1)>β − u B t )) , l = 1, · · · , N − 1,
gN (β,g) = φ (
1 V bu ( 1Vu (M 1W1)>β + g)
) ,
gl(β) = φ (
1 V bu ( 1 Vu (Ml+1Fl+1)>β )) , l = 1, · · · , N − 1.
Then βN [t+ 1] = gt+1N ( gt1 ( · · · gtN−1 ( βN [t],uN−1 B [t] ) · · · ,u1B [t] ) ,g,uN B [t+ 1] ) .
We have:
∥∥βN [t+ 1]− βN [t]∥∥ =
∥∥∥gt+1N (gt1 (· · · gtN−1 (βN [t],uN−1B [t]) · · · ,u1B [t]) ,g,uNB [t+ 1]) −gtN ( gt−11 ( · · · gt−1N−1 ( βN [t− 1],uN−1B [t− 1] ) · · · ,u1B [t− 1] ) ,g,uN B [t] )∥∥∥
≤ ∥∥gN (g1 (· · · gN−1 (βN [t]) · · · ) ,g)− gN (g1 (· · · gN−1 (βN [t− 1]) · · · ) ,g)∥∥ + ∥∥∥gt+1N (gt1 (· · · gtN−1 (βN [t],uN−1B [t]) · · · ,u1B [t]) ,g,uNB [t+ 1])− gN (g1 (· · · gN−1 (βN [t]) · · · ) ,g)∥∥∥
+ ∥∥∥gtN (gt−11 (· · · gt−1N−1 (βN [t− 1],uN−1B [t− 1]) · · · ,u1B [t− 1]) ,g,uNB [t]) −gN ( g1 ( · · · gN−1 ( βN [t− 1] ) · · · ) ,g )∥∥
≤ ∥∥gN (g1 (· · · gN−1 (βN [t]) · · · ) ,g)− gN (g1 (· · · gN−1 (βN [t− 1]) · · · ) ,g)∥∥ + 1
V bu (∥∥∥∥ 1t+ 1 1Vu (M1W1)>gt1 ( · · · gtN−1 ( βN [t],uN−1 B [t] ) · · · ,u1B [t] )∥∥∥∥ +
∥∥∥∥ 1Vu (M1W1)> ( gt1 ( · · · gtN−1 ( βN [t],uN−1 B [t] ) · · · ,u1B [t] ) − g1 ( · · · gN−1 ( βN [t] ) · · · ))∥∥∥∥︸ ︷︷ ︸
A
+ ∥∥∥∥1t 1Vu (M1W1)>gt−11 ( · · · gt−1N−1 ( βN [t− 1],uN−1B [t− 1] ) · · · ,u1B [t− 1] )∥∥∥∥ +
∥∥∥∥ 1Vu (M1W1)> ( gt−11 ( · · · gt−1N−1 ( βN [t− 1],uN−1B [t− 1] ) · · · ,u1B [t− 1] ) − g1 ( · · · gN−1 ( βN [t− 1] ) · · · ))∥∥∥∥︸ ︷︷ ︸
B
+ ∥∥∥∥∥uN B [t+ 1] t+ 1 ∥∥∥∥∥+ ∥∥∥∥∥uN B [t] t ∥∥∥∥∥ ) . (17)
For the term A and B, they are bounded by:
A ≤ 1 V bu (∥∥∥∥ 1Vu (M1W1)> 1Vu (MNFN )> ( gt2 ( · · · gtN−1 ( βN [t],uN−1 B [t] ) · · · ,u2B ) − g2 ( · · · gN−1 ( βN [t] ) · · · ))∥∥∥∥
+ ∥∥∥∥∥ 1Vu (M1W1)>u 1B [t] t ∥∥∥∥∥ )
≤ · · · · · ·
≤ 1 V bu ∥∥∥∥∥ 1Vu (M1W1)>u 1B [t] t ∥∥∥∥∥+ · · ·+ 1V buN−1 ∥∥∥∥∥ 1VuN−1 (M1W1)>(MNFN )> · · · (M3F3)>u N−1B [t] t ∥∥∥∥∥ , (18)
and B has the same form as A by substituting t with t− 1.
Since ‖(M1W1)>‖2‖(MNFN )>‖2 · · · ‖(M2F2)>‖2 ≤ γV Nu V bu N , we have:∥∥gN (g1 (· · · gN−1 (βN [t]) · · · ) ,g)− gN (g1 (· · · gN−1 (βN [t− 1]) · · · ) ,g)∥∥
≤ ∥∥∥∥ 1V bu 1Vu (M1W1)> (g1 (· · · gN−1 (βN [t]) · · · )− g1 (· · · gN−1 (βN [t− 1]) · · · )) ∥∥∥∥ ≤ · · · · · ·
≤ ∥∥∥∥∥ 1V buN 1V Nu (M1W1)>(MNFN )> · · · (M2F2)> ( βN [t]− βN [t− 1] )∥∥∥∥∥ ≤ γ ∥∥βN [t]− βN [t− 1]∥∥ .
(19)
And since uli B [t] is bounded, then ∀ > 0,∃T1 such that when t > T1, we have:∥∥βN [t+ 1]− βN [t]∥∥ ≤ γ ∥∥βN [t]− βN [t− 1]∥∥+ (1− γ)
2 . (20)
Then ‖βN [t + 1] − βN [t]‖ < γt−T1‖βN [T1 + 1] − βN [T1]‖ + 2 , and there exists T2 such that when t > T1 + T2, ‖βN [t + 1] − βN [t]‖ < . According to Cauchy’s convergence test, βN [t] converges to βN ∗, which satisfies βN ∗ = gN ( g1 ◦ · · · ◦ gN−1(βN ∗ ),g ) . Considering the limit, βl[t] converges to βl∗, which satisfies βl∗ = gl(βl+1 ∗ ).
When V bth − ubrest = 1, and there exists λ < 1 such that ‖(M1W1)>‖∞ ≤ λ(Vth − urest), ‖(MlFl)>‖∞ ≤ λ(Vth − urest), l = 2, · · · , N and ‖g‖∞ ≤ 1− λN , we have:∥∥∥βN ∗∥∥∥ ∞ = ∥∥∥gN (g1 ◦ · · · ◦ gN−1(βN ∗),g)∥∥∥ ∞ ≤ ∥∥∥∥ 1Vu (M1W1)>g1 ◦ · · · ◦ gN−1(βN ∗) ∥∥∥∥ ∞ + ‖g‖∞
≤ λ ∥∥∥g1 ◦ · · · ◦ gN−1(βN ∗)∥∥∥
∞ + ‖g‖∞ ≤ · · · · · · ≤ λ
N ∥∥∥βN ∗∥∥∥
∞ + ‖g‖∞ (21)
Therefore, ∥∥∥βN ∗∥∥∥
∞ ≤ ‖g‖∞
1−λN ≤ 1, and ∥∥∥g̃N−1(βN ∗)∥∥∥
∞ ≤ λ ∥∥∥βN ∗∥∥∥ ∞ ≤ λ, · · · · · · ,∥∥∥g̃1 ◦ · · · ◦ g̃N−1(βN ∗)∥∥∥
∞ ≤ λN−1, ∥∥∥g̃N (g̃1 ◦ · · · ◦ g̃N−1(βN ∗),g)∥∥∥ ∞ ≤ λN + (1 − λN ) =
1, where g̃N (β,g) = 1Vu (M 1W1)>β + g, g̃l(β) = 1 Vu (Ml+1Fl+1)>β, l = 1, · · · , N − 1, (i.e. g̃l is gl without the function φ). It means βN ∗ = gN ( g1 ◦ · · · ◦ gN−1(βN ∗ ),g ) =
g̃N ( g̃1 ◦ · · · ◦ g̃N−1(βN ∗ ),g ) and βl∗ = gl(βl+1 ∗ ) = g̃l(β l+1∗).
Taking α1∗ = f1 ( fN ◦ · · · ◦ f2(α1 ∗ ),x∗ ) and αl+1∗ = fl+1(αl ∗ ) (i.e. the fixed-point equation at the equilibrium state as in Section 3.2) explicitly into Eq. (5), the linear equation turns into g̃1 ◦ · · · ◦ g̃N−1(β) − β + g = 0. Therefore, βN
∗ satisfies this equation. And since ‖(M1W1)>‖2‖(MNFN )>‖2 · · · ‖(M2F2)>‖2 ≤ γV Nu , γ < 1, the equation has the unique
solution βN ∗. Further, because g̃l(β) = ( ∂hl+1(α N∗) ∂hl(αN ∗) )> β, where hl(αN ∗ ) = fl ◦ · · · ◦
f2 ( f1(α N ∗,x∗) ) , l = N, · · · , 1, we have βl∗ = ( ∂hN (α N∗) ∂hl(αN ∗) )> βN ∗ , l = N − 1, · · · , 1.

heading: D TRAINING DETAILS

heading: D.1 DROPOUT
text: Dropout is a commonly used technique to prevent over-fitting, and we follow Bai et al. (2019; 2020); Xiao et al. (2021) to leverage variational dropout, i.e. the dropout of each layer is the same at different time steps. Since applying dropout on the output of neurons is a linear operation with a mask and scaling factor, it can be integrated into the weight matrix without affecting the conclusions of convergence. The detailed computation with dropout is illustrated in the pseudocode in Appendix E.

heading: D.2 RESTRICTION ON WEIGHT NORM
text: As indicated in the theorems, a sufficient condition for the convergence to equilibrium states in both forward and backward stages is the restriction on the weight spectral norm. Xiao et al. (2021) leverages re-parameterization to restrict the spectral norm, i.e. they re-parameterize W as W = α W‖W‖2 , where ‖W‖2 is computed as the implementation of Spectral Normalization and α is a learnable parameter to be clipped in the range of [−c, c] (c is a constant). However, the computation of spectral norm and re-parameterization may be hard to realize on neuromorphic hardware. We adjust it for a more friendly calculation as follows.
First, the spectral norm is upper-bounded by the Frobenius norm: ‖W‖2 ≤ ‖W‖F . We can alternatively restrict the Frobenius norm which is easier to compute. Further, considering that connection weights may not be easy for readout compared with neuron outputs, we can approximate
‖W‖F by ‖W‖F = √ tr(WW>) = √ E ∈N (0,Id) [‖ >W‖22], according to the Hutchinson estimator (Hutchinson, 1989). It can be viewed as source neurons outputting noises and target neurons accumulating signals to estimate the Frobenius norm. We will estimate the norm based on the Monte-Carlo estimation (we will take 64 samples), which is similarly adopted by Bai et al. (2021) to estimate the norm of their Jacobian matrix. Then based on the estimation, we will restrict W as W = α W‖W‖F where α = min(c, ‖W‖F ), c is a constant for norm range. This estimation and calculation may correspond to large amounts of noises in our brains, and a feedback inhibition on connection weights based on neuron outputs.
Following Xiao et al. (2021), we only restrict the norm of feedback connection weight W1 for the multi-layer structure, which works well in practice.

heading: D.3 OTHER DETAILS
text: For SNN models with feedback structure, we set Vth = 1, urest = −1 in the forward stage to form an equivalent equilibrium state as Xiao et al. (2021). The constant for restriction in Appendix D.2 is c = 2. Following Xiao et al. (2021), we train models by SGD with momentum for 100 epochs. The momentum is 0.9, the batch size is 128, and the initial learning rate is 0.05. For MNIST, the learning rate is decayed by 0.1 every 30 epochs, while for CIFAR-10 and CIFAR-100, it is decayed by 0.1 at the 50th and 75th epoch. We apply linear warmup for the learning rate in the first 400 iterations for CIFAR-10 and CIFAR-100. We apply the weight decay with 5× 10−4 and variational dropout with rate 0.2 for AlexNet-F and 0.25 for CIFARNet-F. The initialization of weights follows Wu et al. (2018), i.e. we sample weights from the standard uniform distribution and normalize them on each output dimension. The scale for the loss function (as in Section 4.3) is 100 for MNIST, 400 for CIFAR-10, and 500 for CIFAR-100.
For SNN models with degraded feedforward structure, our hyperparameters mostly follow Thiele et al. (2019a), i.e. we set Vth = 0.5, urest = −0.5, train models by SGD with momentum 0.9 for
60 epochs, set batch size as 128, and the initial learning rate as 0.1 which is decayed by 0.1 every 20 epochs, and apply the variational dropout only on the first fully-connected layer with rate 0.5.
The notations for our structures mean: ‘64C5’ represents a convolution operation with 64 output channels and kernel size 5, ‘s’ after ‘64C5’ means convolution with stride 2 (which downscales 2×) while ‘u’ after that means a transposed convolution to upscale 2×, ‘P2’ means average pooling with size 2, and ‘F’ means feedback layers. The network structures for CIFAR-10 are:
AlexNet (Wu et al., 2019): 96C3-256C3-P2-384C3-P2-384C3-256C3-1024-1024,
AlexNet-F (Xiao et al., 2021): 96C3s-256C3-384C3s-384C3-256C3 (F96C3u),
CIFARNet (Wu et al., 2019): 128C3-256C3-P2-512C3-P2-1024C3-512C3-1024-512,
CIFARNet-F (Xiao et al., 2021): 128C3s-256C3-512C3s-1024C3-512C3 (F128C3u).
We simulate the computation on commonly used computational units. The code implementation is based on the PyTorch framework (Paszke et al., 2019), and experiments are carried out on one NVIDIA GeForce RTX 3090 GPU.

heading: E PSEUDOCODE OF THE SPIDE ALGORITHM
text: Our algorithm consists of two-stage SNN computation, as explained in Section 4.4. The detailed computation for both stages are illustrated in Algorithm 1 and Algorithm 2 respectively.
Algorithm 1 Forward procedure of SPIDE training - Stage 1. Input: Network parameters F 1, b1, · · · , FN , bN ,W 1,W o, bo; Input data x; Time steps TF ; Forward threshold Vth; Dropout rate r; Output: Output of the readout layer o.
1: Initialize ui[0] = 0, i = 1, 2, · · · , N 2: If use dropout, randomly generate dropout masks Di(i = 1, 2, · · · , N), Df with rate r // Di, Df are saved for backward 3: for t = 1, 2, · · · , TF do 4: if t == 1 then 5: u1[t] = u1[t− 1] +D1 (F 1x[t] + b1) 6: else 7: u1[t] = u1[t− 1] +D1 (F 1x[t] + b1) +Df (W 1sN [t− 1]) 8: s1[t] = H(u1[t] ≥ Vth) 9: u1[t] = u1[t]− 2Vths1[t] // urest = −Vth, the same below 10: for l = 2, 3, · · · , N do 11: ul[t] = ul[t− 1] +Dl (F lsl−1[t] + bl) 12: sl[t] = H(ul[t] ≥ Vth) 13: ul[t] = ul[t]− 2Vthsl[t] 14: o = o+W osN [t]+bo // o can be accumulated here, or calculated later by o =W oαN+bo
15: o = oT 16: αi = ∑T t=1 s i[t]
T , i = 1, 2, · · · , N // Save for backward, firing rate in Stage 1 17: mi = H(αi > 0) ∧H(αi < 1) // Save for backward, mask 18: If x is not constant, save x = ∑T t=1 x[t]
T for backward 19: return o

heading: F MORE EXPERIMENTAL COMPARISONS

heading: F.1 ABLATION STUDY
text: In this section, we conduct ablation study on our improvement to reduce the approximation error by setting the resting potential as negative threshold. To formulate equivalent equilibrium states, we
Algorithm 2 Backward procedure of SPIDE training - Stage 2. Input: Network parameters F 1, b1, · · · , FN , bN ,W 1,W o, bo; Forward output o; Label y; Time steps TB ; Forward threshold Vth; Backward threshold V bth = 0.5; Other hyperparameters and saved variables; Output: Trained network parameters F 1, b1, · · · , FN , bN ,W 1,W o, bo.
1: Calculate g = ∂L(o,y)∂o // for CE loss, ∂L(o,y) ∂o = softmax(o)− y, in practice we will scale the
loss by a factor sl, then ∂L(o,y) ∂o = sl (softmax(o)− y)
2: Initialize ui[0] = 0, i = 1, 2, · · · , N 3: for t = 1, 2, · · · , TB do 4: if t == 1 then 5: uN [t] = uN [t− 1] +W o>g 6: else 7: uN [t] = uN [t − 1] +W o>g + 12VthW
1>(Df m1 s1[t − 1]) // mi is the saved mask in Stage 1
8: sN [t] = T (uN [t], 0.5) // realized by two coupled neurons 9: uN [t] = uN [t]− sN [t] // realized by two coupled neurons
10: for l = N − 1, N − 2, · · · , 1 do 11: ul[t] = ul[t− 1]+ 12VthF
l>(Dl ml+1 sl+1[t]) // mi is the saved mask in Stage 1 12: sl[t] = T (ul[t], 0.5) // realized by two coupled neurons 13: ul[t] = ul[t]− sl[t] // realized by two coupled neurons 14: βi = ∑T t=1 s i[t]
T , i = 1, 2, · · · , N // “firing rate” in Stage 2 15: Calculate gradients: 16: (1) ∇F 1L = 12Vth (m
1 β1)x> // mi, x are the saved mask and average input in Stage 1 17: (2)∇F iL = 12Vth (m
i βi)αi−1>, i = 2, 3, · · · , N // mi, αi are the saved mask and firing rate in Stage 1
18: (3) ∇biL = 12Vth (m i βi), i = 1, 2, · · · , N 19: (4) ∇W 1L = 12Vth (m 1 β1)αN> // mi, αi are the saved mask and firing rate in Stage 1
20: (5) ∇W oL = αN ( ∂L(o,y) ∂o )> // αi is the saved firing rate in Stage 1
21: (6) ∇boL = ( ∂L(o,y) ∂o )> 22: Update F 1, b1, · · · , FN , bN ,W 1,W o, bo based on the gradient-based optimizer // SGD
learning rate η + momentum α & weight decay µ, the base learning rate is scaled by the factor sl of the loss, i.e. η = ηsl
23: (1) Update the momentum Mθ = α ∗Mθ + (1− α) ∗ ∇θL, θ ∈ {F i, bi,W 1,W o, bo} 24: (2) Update parameters θ = (1− µ) ∗ θ + η ∗Mθ, θ ∈ {F i, bi,W 1,W o, bo} 25: (3) Restrict the norm of W 1 26: return F 1, b1, · · · , FN , bN ,W 1,W o, bo
take the same Vth − urest = Vu and the same V bth − ubrest = V bu , and we consider the following settings: (1) both forward and backward stages apply our improvement, i.e. urest = −Vth, ubrest = −V bth; (2) remove the improvement on the backward stage, i.e. V bth = V bu , ubrest = 0; (3) remove the improvement on both forward and backward stages, i.e. Vth = Vu, urest = 0 and V bth = V b u , u b rest = 0. The latter two setting are denoted by “w/o B” and “w/o F&B” respectively.
The models are trained on CIFAR-10 with AlexNet-F structure and 30 forward time steps. The training and testing curves under different settings and backward time steps are illustrated in Figure 3 and Figure 4 respectively. It demonstrates that without our improvement, the training can not perform well within a small number of backward time steps, probably due to the bias and large variance of the estimated gradients. When the backward time steps are large, the performance gap is reduced since the bias of estimation is reduced. It shows the superiority of our improvement to training SNNs within a small number of backward time steps.

heading: F.2 COMPARISON TO IDE
text: Table 3 in Section 5 shows that the SPIDE method performs poorer than the original IDE method (Xiao et al., 2021). We further investigate the training and testing curves during optimization to analyze this phenomenon. As shown in Figure 5, the SPIDE method could achieve the same training accuracy as the IDE method, while the generalization performance is poorer. Since the hyperparameters are the same for experiments, except that we drop the modified BN component (as explained in Section 4.4), the performance gap may be caused by the implicit regularization effect of BN. Therefore, the optimization ability of the SPIDE method should be similar to that of IDE, while future work would be to investigate how to realize operations similar to BN with possible computation that could be friendly on neuromorphic hardware.
Also, we note that models trained by the IDE method (Xiao et al., 2021) have sparser spikes, as their forward average firing rates are only around 1% (Xiao et al., 2021), while ours are around 7% as shown in Figure 2. It is again probably due to BN which would subtract the statistical mean value of neuron inputs, therefore regularizing the weights so that neurons will generate sparser spikes. An interesting future work is to further reduce the number of spikes considering this phenomenon.

heading: F.3 RESULTS ON CIFAR-100
text: In this section, we present the results on CIFAR-100. As shown in Table 4, our model could achieve 64.07% accuracy. Compared with IDE, the performance is poorer, and the main reason is probably again the absence of BN which could be important for alleviating overfitting on CIFAR-100 with relatively small number of images per class. The training accuracy of SPIDE is similar to IDE (around 93% v.s. around 94%) while the generalization performance is poorer. Despite this, the performance of our model is competitive for networks without BN and our model is with fewer neurons and parameters and a small number of time steps. Compared with SpikeGrad (Thiele et al., 2019a), we can use fewer neurons and parameters due to flexible network structure choices, and we leverage common neuron models while they require complex impractical models. Future work could investigate more suitable structures and if there are normalization techniques friendly for neuromorphic computation and our desired algorithm to further improve the performance.

heading: F.4 RESULTS ON CIFAR10-DVS
text: In this section, we supplement some results on the spiking dataset CIFAR10-DVS (Li et al., 2017). The CIFAR10-DVS dataset is the neuromorphic version of the CIFAR-10 dataset converted by a Dynamic Vision Sensor (DVS), which is composed of 10,000 samples, one-sixth of the original CIFAR-10. It consists of spike trains with two channels corresponding to ON- and OFF-event spikes. The pixel dimension is expanded to 128× 128. Following the common practice, we split the dataset into 9000 training samples and 1000 testing samples. As for the data pre-processing, we reduce the time resolution by accumulating the spike events (Fang et al., 2021) into 30 time steps, and we reduce the spatial resolution into 48× 48 by interpolation. We apply the random crop augmentation as CIFAR-10 to the input data. We leverage the network structure: 512C9s (F512C5), where the notations follow Appendix D.3. We train the model by SGD with momentum for 70 epochs. The momentum is 0.9, the batch size is 128, the weight-decay is 5× 10−4, and the initial learning rate is 0.05 which is decayed by 0.1 at the 50th epoch. No dropout is applied. The initialization of weights follows the widely used Kaiming initialization. The constant for restriction in Appendix D.2 is c = 10 due to the large channel size, and the scale for the loss function as well as the firing thresholds and resting potentials are the same as the CIFAR-10 experiment.
As shown in Table 5, our model could achieve 60.7% accuracy, which is competitive among results of common SNN models, demonstrating the effectiveness of our method. Fang et al. (2021) leverages many techniques such as learnable membrane time constant, batch normalization, and max pooling to achieve better performance. We do not aim at outperforming the state-of-the-art results, but demonstrate that a competitive performance could be achieved even with our constraints of purely spike-based training in a relatively small number of time steps, verifying the effectiveness.

Review summaries:
In this paper, the authors proposed a method to train Spiking Neural Networks (SNN) with spike-based implicit differentiation on the equilibrium state. Main idea is to use a spike-triggered event instead of average firing rate to approximate implicit differentiation of Feedback Spiking Neural Networks (FSNN). To enable such idea and to further reduce the approximation errors, the authors proposed several techniques such as adopting ternary spiking neuron couples and shifting resting potential. The experimental results showed that the proposed method can achieve high accuracy in several tasks such as MNIST, CIFAR-10, and CIFAR-100 with fewer time steps for training compared to existing methods.

The paper aims at porting the IDE method  into a spike-fbased and more bio-plausible version. The previous IDE used firing rates rather than spikes for computation, although reference Xiao et al in NeurIPS 2021 had already addressed implementations in spiking neural networks. The authors analyze the approximation error that results from solving implicit differentiation by spikes and report a solution based on ternary spiking neurons, that can be implemented with pairs of standard spiking neurons. They achieve in this way quite good performance for MNIST and CIFAR10.