Year: 2021

Abstract:
Although spiking neural networks (SNNs) show strong advantages in terms of information encoding, power consumption, and computational capability, the underdevelopment of supervised learning algorithms is still a hindrance for training SNNs. We consider weight initialization to be a pivotal issue for efficient SNN training. It greatly influences gradient generation when using the method of backpropagation through time in the initial training stage. Focusing on the properties of spiking neurons, we first derive an asymptotic formula for their response curve, approximating the actual neuron response distribution. Then, we propose an initialization method based on the slant asymptote to overcome gradient vanishing. Finally, experiments with different coding schemes in classification tasks on the MNIST and CIFAR10 dataset show that our method can effectively improve the training speed and the model accuracy compared with traditional deep learning initialization methods and existing SNN initialization methods. Further validation on different neuron types and training hyperparameters shows comparably good versatility and superiority over other methods. Some suggestions are given for SNN training based on the analyses.

heading: 1 INTRODUCTION
text: To date, deep learning has contributed to outstanding performance on various tasks (Hinton et al., 2012; Hosu & Rebedea, 2016; He et al., 2016; LeCun et al., 2015). Most state-of-the-art models are based on analog neurons. Each analog value conveys activation; thus, the inference processes of analog neural networks (ANNs) consist of massive matrix manipulations. To overcome this, spiking neural networks (SNNs) are increasing attracting from researchers (Maass, 1997; Maass & Markram, 2004; Yang et al., 2019). The spiking neurons in an SNN emit sparse 0,1 spikes when activated. Compared with existing ANNs, SNNs are considered to encode spatiotemporal patterns and have better power efficiency. Hence, implementations of SNNs are highly appealing in powerlimited scenarios. In addition, SNNs have exhibited computational capabilities as powerful as those of conventional artificial neural networks (Lee et al., 2016).
Despite SNNs’ powerful capabilities, these networks lack efficient supervised learning algorithms. The discrete spikes prevent the direct use of backpropagation (BP) training. Recently proposed supervised algorithms have mainly focused on two methodologies: ANN-SNN conversion and BP through time (BPTT) with surrogate functions. Compared to ANN-SNN conversion (Diehl et al., 2015; Sengupta et al., 2019; Deng & Gu, 2021; Ho & Chang, 2020; Ding et al., 2021), BPTT is more friendly to the event-based datasets (Lee et al., 2016; Zheng et al., 2021; Neftci et al., 2019). However, Wu et al. (2018) has reported that BPTT is only feasible to shallow architectures. Therefore, BPTT still cannot be regarded as a mature SNN training algorithm.
We ascribe the immaturity of BPTT training to the influence of the gradients. One common reason is gradient vanishing and explosion. This is natural since SNNs can be regarded as a type of recurrent neural networks (RNNs) (He et al., 2020; Fang et al., 2020). In the wisdom of RNN methodology, researchers have managed to provide a simple solution by rethinking weight initialization (Le et al., 2015). The success of RNN training inspires consideration of whether there may also be a simple solution for SNNs. Another reason is early gradient generation, which has rarely been mentioned in the previous literature. Most surrogate functions produce gradients when the membrane potential is in the neighborhood of the threshold. Some produce gradients only when spikes are fired (Lee
et al., 2020a). Thus, the flow of the gradients strongly depends on the number of firing neurons and will influence the optimization of the network. A small experiment can help illustrate how early gradient generation affects training. As shown in Fig. 1, with an inadequate weight initialization method, the neurons receive inputs that are too small and cannot fire initially, leading to stall of the training process because gradients fail to be produced. Until some neurons fire, the loss does not begin to decrease, and the neurons must suffer a latter uncoordinated training schedule. Therefore, in general, good weight initialization is also important for efficient SNN training.
At present, initialization for SNN training is mainly achieved by following and adapting methods originally developed for ANNs. However, unlike the rectified linear unit (ReLU) activation function conventionally used in ANNs, spiking neurons respond very differently, and recent improvements in models have rarely focused on the unique properties of spiking neurons. In this paper, we investigate the theoretical response of spiking neurons to develop a better SNN initialization strategy. The main contributions of this paper are summarized as follows:
• We use iterative systems to model first-order integrate-and-fire neurons and theoretically investigate the response curve. A slant asymptote is derived to approximate the actual neuron response.
• We examine the difficulty of training deep SNNs and develop an initialization method based on the slant asymptote. The proposed initialization method relies on adequate scaled weights and nonzero biases to enable early gradient training.
• Experiments with rate coding and event coding in classification tasks show that our method can effectively improve the model training speed and accuracy. In addition, empirical analyses suggest that our asymptote initialization method offers versatility and superiority over a vast range of training hyperparameters, such as different neuron types, surrogate functions, and optimizers.

heading: 2 THEORETICAL ANALYSIS OF SNNS

heading: 2.1 SPIKING NEURONS AND SURROGATE FUNCTIONS
text: With the development of neuroscience, many neuron models have been proposed (Hodgkin & Huxley, 1952; Gerstner et al., 2014), commonly in the form of differential systems. However, in the field of computing, analyses are often based on first-order differential models because of computational efficiency. Typically, leaky and nonleaky integrate-and-fire (LIF and IF) models can be formulated in terms of the dynamics of the membrane potential v as shown in Eq. 1(a) and (b).
τ dv(t)dt = −(v(t)− vreset) +Ri(t) (a) τ is time constant. di(t) dt = 1 C i(t) (b)
(1)
By utilizing the Euler method, we can convert these dynamic systems into an iterative expression that computers prefer (Eq. 2). Note that vreset is usually 0 in the literature; hence, it is neglected in latter analyses. When an element of vt exceeds the firing threshold θ, the corresponding neuron will emit a spike of 1, and membrane potential will then be reset to 0. Without loss of generality, we use the following equations to describe the iterative systems for first-order integrate-and-fire neurons.
ut = kvt + λit
ot = ut > θ
vt+1 = ut(1− ot) (2)
In Eq. 2, k and λ are hyperparameters that control the neuron type and firing pattern (k ∈ (0, 1]). When k is set to 0, the iteration of the membrane potential over time is canceled, and the network is not a standard SNN. Many neuron models used in the existing literature on SNN training can be regarded as the model given in Eq. 2 with different settings of the two parameters. Table 1 illustrates some of these settings.
In Eq. 2, ot = ut > θ is the non-differentiable part. It can also be rewritten as ot = H(ut − θ) where H(·) is the Heaviside step function. H(x) = 0 when x ≤ 0, and H(x) = 1 when x > 0. The heaviside function has a gradient of 0 almost everywhere. Surrogate function methods first emerged in network quantization to cope with non-differentiable problems (Mairal, 2013; Hubara et al., 2016). Accordingly, four surrogate functions (sigmoid, tanh, hardtanh, arctan) are presented in this paper, parameterized by the parameter a, as listed in Table A3. These functions are the shifted and scaled version of the original functions, as shown in Fig. A7 in the appendix. The parameter a can be adjusted to rescale the output gradient. The gradient functions are typically symmetric w.r.t. an input of zero. An exception has been proposed by Lee et al. (2020a), such that only the production of spikes can trigger nonzero gradients.

heading: 2.2 RESPONSE OF SPIKING NEURONS AND ASYMPTOTES
text: It is not easy to establish the relationship between the current input and the output based on an iterative equation. Considering that the difference between the iterative equation and the differential equation is small when stepsize = 1, and that the differential equation can be integrated to obtain the response curve, Theorem 1 gives the response of a first-order integrate-and-fire neurons. Theorem 1. For a first-order integrate-and-fire neuron as expressed in Eq. 2 (ut = kvt + λit, k ∈ (0, 1]), the model takes a constant current it = i as input and produces spikes as output. Let the discrete output firing rate w.r.t. i be expressed as f = TN , where T is the total inference time and N is the spike count. Then, the input-output response curve (f -i curve) is estimated as follows:
If k = 1 :
f =  0, i ≤ 0, λ θ i, 0 < i < θ λ ,
1, i ≥ θλ . If k ∈ (0, 1) :
f =  0, i ≤ 0, (k − 1)/ ln(1 + (k−1)θλi ), 0 < i < θ λ · k−1 ek−1−1 ,
1, i ≥ θλ · k−1 ek−1−1 .
(3)
The solid lines in Fig. 2 represent the response curves when k = 1 and k ∈ (0, 1) in Theorem 1. Fig. 2 also presents scatter plots of the actual response of a spiking neuron given the corresponding parameters (marked with circular dots). The actual neuron responses match the theoretical f − i curve for most inputs. Consider the region in which f ∈ (0, 1). When k ∈ (0, 1), it can be intuitively found that the f -i curve, which has a non-differentiable point, seems to exhibit a slant asymptote. Theorem 2 manifests the existence of this slant asymptote and provides a uniform expression for it. In Fig. 2, this asymptote approximately crosses the actual neuron response distribution. In this case, the slant asymptote of f -i can better describe the relationship between the expected input current and the output firing rate. Theorem 2. For a first-order integrate-and-fire neuron as expressed in Eq. 2 (ut = kvt + λit, k ∈ (0, 1]), the slant asymptote of the f -i curve satisfying f ∈ (0, 1) exists and can be uniformly expressed as follows:
fa = λ
θ i+ k − 1 2 . (4)

heading: 3 WEIGHT INITIALIZATION FOR SNN TRAINING
text: The gradient of the weights is related to the spike output. To better illustrate this, given an L-layered SNN, we denote the input current, membrane potential, intermediate variable, and spike output of the lth layer by il,t,vl,t,ul,t, and ol,t, respectively. The variables for the neurons in one layer satisfy Eq. 2. il,t is the weighted input consisting of the previous layer’s spike output (il,t = W lol−1,t+bl), and the first layer takes the encoded spike o0,t as input. When we backpropagate the loss ε through time T , the gradients w.r.t. the weights of the lth layer can be expressed as shown in Eq. 5. An inappropriate ol−1,t, such as an output of zero, will lead to a poor gradient.
∂ε
∂W l = T∑ t=1 ∂ε ∂ul,t ∂ul,t ∂il,t ∂il,t ∂W l = λ T∑ t=1 ∂ε ∂ul,t (ol−1,t)T (5)
Therefore, the main concern is how to produce spikes to trigger the learning process. Thankfully, the clipped version of the slant asymptote can replace the analysis of complex iterative functions, allowing us to focus on the scale of the parameters:
îl = W lôl−1 + bl
ôl = clip( λ
θ îl + k − 1 2 , 0, 1) (6)
where îi and ôi are the expected current and firing rate vector, respectively, of the ith layer with ni neurons and W i and bi are the corresponding parameters. Note that W i is an ni × ni−1 matrix. λ, k and θ represent the neuron parameters.
First, all weights in the SNN are initialized to be independent and identically distributed (i.i.d.). Let us assume that the input current variances are the same. Then, we have:
Var(̂i l ) = nl−1Var(W lôl−1) (7)
where î l , Wl and ôl−1 denote random variables of elements in îl, W l and ôl−1, respectively. Considering Var(Wlôl−1), Var(Wlôl−1) = Var(Wl)E[(ôl−1)2] + Var(ôl−1)E[Wl]2, let Wl be of zero mean; then:
Var(̂i l ) = nl−1Var(W l)E[(ôl−1)2]. (8)
In the function ô = clip(λθ î+ k−1 2 , 0, 1), i = θ(1−k) 2λ is the transition point between the nonspiking region and the linear response region. When bl is initialized as θ(1−k)2λ (l = 1, 2, · · · , L), a zerocentered W l × ôl−1 produces a symmetric distribution around zero. With bl initialized as θ(1−k)2λ ,
E[(ôl−1)2] = 1
2 Var(
λ θ î l−1 ) = 1 2
λ2 θ2 Var(̂i l−1 ). (9)
By plugging E[(ôl−1)2] into Eq. 8, we can obtain:
Var(̂i l ) =
1
2
λ2 θ2 nl−1Var(W l)Var(̂i l−1 ) (10)
To produce spikes is to maintain E[(ôl)2], which is related to Var(̂i l ) according to Eq. 9. Therefore, we determine the distribution of Wl as follows (for l = 2, 3, · · · , L):
Var(̂i l )
Var(̂i l−1 ) =
1
2
λ2 θ2 nl−1Var(W l) = 1 (11)
Accordingly, W l is initialized from a Gaussian distribution N (0, 2θ 2
λ2nl−1 ), where nl−1 is the number
of input neurons in layer l. As done in other initializers, we apply the initialization to the first layer for simplicity.
In previous work on ANNs, there is no need for the neurons to accumulate input. Thus, the biases are often set to zero in an ANN. For SNNs, because the response asymptote indicates a nonzero nondifferentiable point, it is natural for the biases in an SNN to be initialized to that non-differentiable value. We validate the impact to the mean firing rate of our initialization method in Fig. 3. Our method can effectively trigger spikes when compared with ANN-based methods.

heading: 4 RELATED WORK
text: From the perspective of the loss landscape, a good weight initialization is often located in a region with low loss values (Li et al., 2018). However, training deep neural networks is never an easy thing. After trails of layerwise training of deep nets (Hinton, 2007; Bengio et al., 2006), Glorot & Bengio (2010) provided a theoretical view on the prevention of exploding and vanishing gradients. They estimated the standard deviation of the weights from the number of layer neurons and initialized the biases to zero. He et al. (2015) extended their hypothesis to a piecewise linear case to cater to the growing prevalence of ReLU activation. Both methods are widely used in deep ANNs. There are also works for RNN. Saxe et al. (2014) examined the dynamics of learning and derived random
orthogonal initial conditions for deep nonlinear networks. Le et al. (2015) used an identity matrix to initialize the RNN weight matrix and thus maintain the gradients.
Recent SNN studies have also included some initialization methods. For example, one can regard ANN-SNN conversion as an initialization method (Rathi et al., 2020). In addition, Lee et al. (2016) set the weights and thresholds using a manually tuned hyperparameter. They also proposed a complex BP error normalization procedure. Wu et al. (2018) suggested normalizing weights sampled from a standard uniform distribution, with no further explanation. Lee et al. (2020b) followed the recommendation of He et al. (2015) with a manually selected constant κ.
With the similar motivation of this work, many recent studies have focused on regularizing the spiking response. Zheng et al. (2021) and Kim & Panda (2020) proposed batch normalization along channel axes and time axes, respectively, for BPTT; both of which are layerwise normalization methods. In addition, Kim & Panda (2021) utilized a global unsupervised firing rate normalization. The above methods lack direct attention to neuronal dynamics. Sengupta et al. (2019), however, used dynamics of simple IF neuron proposing Spike-Norm for ANN-SNN conversion. In contrast, our work aims to improve the training for general first-order integrate-and-fire neurons by analyzing neuronal dynamics.

heading: 5 EXPERIMENTAL VALIDATION

heading: 5.1 EXPERIMENTAL SETTING AND DATASETS
text: To verify the training acceleration performance and the model performance against various other initialization methods, we conducted experiments on classification tasks involving different codings for the MNIST dataset (Lecun et al., 1998) and the CIFAR10 dataset (Krizhevsky et al., 2009). Before being fed into an SNN, images should be converted into spike trains. We implemented constant coding, Poisson coding, retina-like coding, and event coding for both datasets. Poisson coding and constant coding are popular rate coding schemes for ANN-SNN conversion (Rueckauer et al., 2017; Sengupta et al., 2019). Retina-like coding mimics the integration process of a bipolar cell (Zhu et al., 2019). The N-MNIST dataset was produced by converting MNIST into spike events with dimensions of size 34× 34 using saccades (Orchard et al., 2015). The CIFAR10-DVS dataset comprises 10000 128× 128 event streams converted from CIFAR10 using a dynamic vision sensor (DVS) (Li et al., 2017).
Previous work has revealed that the accuracy of an SNN is sensitive to the choice of hyperparameters such as the neuron parameters (Fang et al., 2020). Therefore, we fixed the neuron parameters for a given dataset, as shown in Table 2, and compared our method with four other prevalent choices
(Glorot & Bengio (2010), He et al. (2015), Lee et al. (2016) and Wu et al. (2018)). An analysis of the hyperparameter sensitivity is presented in Sec. 6. For MNIST, we trained a 5-layer SNN. For the CIFAR10 dataset, we trained a VGG-9 model. For the CIFAR10-DVS dataset, we added one pooling layer and a fully-connected layer into the VGG-9 model. Time step in forwarding process is set to a constant value of 20. All the models were optimized using the Adam optimizer. The detailed training hyperparameters are shown in the appendix. Additionally, we conducted experiments for larger datasets, i.e., CIFAR100 and TinyImageNet, to compare to Batch Normalization (BN). See Appendix for the detailed descriptions and comments on experiments on larger datasets.

heading: 5.2 RESULTS
text: In general, Fig. 4 and Table 2 show that good parameter initialization can not only accelerate training but also lead to some performance improvement. On MNIST, the different initialization methods achieve similar performance in validation accuracy with all four encoding schemes listed in Table 2, and our method performs slightly better for three of them. On CIFAR10, our method achieves performance improvement, showing accuracy increases of 4.67%, 1.20%, 6.62%, and 0.5% compared with He’s method. As seen from the first row of Fig. 4, our method enables the model to generate early gradients in fewer iterations and cause the loss to drop faster. In more detail, our method allows the model to obtain a proper parameter distribution, which helps every layer coordinate with the input and output to better adapt to SNN training. This advantage enables network training to begin quickly, within a few iterations. In comparison, the other methods are slower than ours.
For MNIST, our method reaches 98% validation accuracy after only 1 epoch, while other methods reach this accuracy after about 2 or 3 epochs of training. For N-MNIST, our method reaches 98% accuracy after 2 epochs, while the others reach after 3 epochs. For CIFAR10, our method first achieves 60% accuracy at approximately the 4th epoch, whereas Wu’s, He’s, Lee’s, and Glorot’s methods reach this level later, respectively. For CIFAR10-DVS, our method and those of Wu and He reach 60% accuracy in order, while training fails under Lee’s and Glorot’s methods.
Due to the simplicity of MNIST, networks taking input encoded using each of the four schemes perform similarly. For CIFAR10, a model trained with constant coding achieves the highest accuracy because this data format is almost the same as the original data. Poisson coding performs worse than constant coding due to its more complicated coding rules. Retina-like coding leads to considerable degradation, probably due to the loss of some texture information. The results show that this coding scheme influences the other four initialization methods more than it does ours. The event data (CIFAR10-DVS) are sparse and complicated, and the amount of data is much smaller than that in CIFAR10, which increases the training difficulty and leads to performance degradation.

heading: 6 ANALYSIS OF EMPIRICAL RESULTS
text: To study whether the initialization strategy is sensitive to the hyperparameters of the model and the training process, we conducted many experiments. In these experiments, we used the MNIST dataset and the same 5-layer network as in Sec. 5. Training was performed for 3 epochs with a minibatch size of 64.

heading: 6.1 PERFORMANCE ANALYSIS ON DIFFERENT NEURONS
text: We performed an analysis of different neuron types, that is, different neuron parameter settings. While θ was fixed to 1.0, k was sampled from 0.1 to 1.0, and λ was sampled from 0.1 to 4.0. Fig. A8 compares the training speeds in this grid of experiments. For most (k, λ) settings, our method shows obvious acceleration relative to the other methods. In particular, our method can be used to effectively train networks for which training cannot be successfully started using other initialization methods.
We choose six representative groups to better illustrate this as shown in Fig. 5. For leaky neurons with different τ values, with decreasing τ , the training speed and accuracy of all initialization methods improve. In this case, our method maintains high accuracy and achieves a training acceleration of more than 11× compared to the baseline Glorot initialization. For nonleaky neurons, as shown in panels (a) and (e), the difference among methods is not obvious. As λ decreases, our method enables faster training. For neurons that leak quickly (k = 0.2), Glorot’s method performs poorly in terms of the start-up training speed but achieves comparable accuracy to our method. Among the faster methods, our method achieves the best accuracy.
Besides, instead of fixing θ, we plot the heatmap in Fig. A9 after sampling θ from {0.5, 1, 2, 4} . This figure illustrates that our method is suitable for almost every chosen parameter and can help address the stalling of training caused by extremely large threshold due to the compensation of the initial biases.

heading: 6.2 PERFORMANCE ANALYSIS WITH DIFFERENT OPTIMIZERS AND SURROGATE FUNCTIONS
text: We selected τ = 5 and 2 to test the sensitivity to the training hyperparameters. Two main factors that influence the gradient update are the optimizer and the surrogate function. Accordingly, Fig. 6 shows the training curves with different initialization methods, and the accuracy is marked in the figure in the order of the legends. Generally, if training can be successfully started via the gradient update with a given hyperparameter, the difference in accuracy is not huge. However, we can still observe an accuracy advantage with our method. For example, for (0.8, 0.2, 1.0), the accuracy of Wu’s and Lee’s methods is 1.8% lower than ours. In addition, roughly, the initialization methods can be roughly arranged in the following order by the training speeds observed from the experiments: ours > Wu’s > He’s ≥ Lee’s > Glorot’s. Moreover, better training is observed in experiments with a smaller τ .
Regarding the optimizers, SGD(lr=0.5) and ASGD(lr=2.0) fail to activate fast early training for other initialization methods, yet our methods still works. Compared with Adam(lr=1e-3,β =(0.9,0.99)), AdaGrad(lr=1e-2) seems less sensitive to the choice of the initialization methods.
With different surrogate functions, the order of the methods is maintained. Hardtanh produces coarse surrogate gradients and suppresses training under the other four methods. Surrogate functions that produce more fine-grained gradients all successfully activate training. Sigmoid(a=2.0) has some advantages in terms of training speed, but it is slightly inferior to tanh(a=2.0) and arctanh(a=1.0) in training accuracy.

heading: 7 DISCUSSIONS AND LIMITATIONS
text: The large number of experiments we conducted allows us to provide suggestions for SNN training. A good combination of hyperparameters is essential. Our initialization method can speed up training in most cases and adapt to more difficult neural parameters. In addition, we recommend the use of adaptive optimizers, as it helps accumulate the gradients. Finally, the arctan and tanh surrogate functions are suggested because they offer high accuracy and relatively high training speeds.
Apart from the advantages, our method still has some limitations in terms of neuronal dynamics. Although analyzing the response of spiking neurons and asymptotes is feasible in most cases, the theoretical results will be hampered if the input current carries more noise. The noise will have at least two impacts (Gerstner et al., 2014). One is that the neurons will fire even when the mean input is below zero. The other is that noise in the input current will cause the membrane potential to drift away from the theoretical asymptote. Furthermore, the noise will increase when we extend to wider and deeper networks. Therefore, for large datasets such as ImageNet, the asymptote should take the current noise into account.

heading: 8 CONCLUSION
text: In this work, we aim to enhance the efficient BPTT training of SNNs by proposing an initialization method that is consistent with the response of spiking neurons in initial training. Our method bridges the spiking neuron response with the wisdom of traditional deep learning training, which may have an influence on future research on topics such as ANN-SNN conversion with LIF neurons or other SNN training methods. Our experimental results validate that our initialization method shows superiority, versatility, and robustness over a vast training hyperparameter space.

heading: A APPENDIX
text: A.1 APPENDIX FOR THEORETICAL ANALYSIS
Here, we provide theoretical proofs of the theorems in Sec. 2. These are the foundation of the analyses in this paper. Additionally, Fig. A7 is presented here to illustrate the surrogate functions and symmetric gradients.
−2 0 2 Input
(a)
0.0
0.5
1.0
1.5
2.0
2.5
Ou tp
ut
Heaviside Sigmoid(a=2) Tanh(a=2) Hardtanh(a=1) Arctan(a=1)
−2 0 2 Input (b)
0.0
0.5
1.0
1.5
2.0
2.5
Gr ad
ien t
Dirac Sigmoid(a=2) Tanh(a=2) Hardtanh(a=1) Arctan(a=1)
Figure A7: Surrogate functions and corresponding gradients.
Table A3: Surrogate functions for the non-differentiable model (a is the parameter).
Surrogate function h h(x, a) h′(x, a)
Sigmoid-like 1/(1 + e−ax) a · h(x, a)(1− h(x, a)) Tanh-like 12 tanh(ax) + 1 2 1 2a · (1− tanh(ax)
2) Hardtanh-like clip(ax+ 12 , 0, 1) a ·Heaviside( 1 2a − |x|) Arctan-like 1πarctan(πax) + 1 2 a/(1 + (πax) 2)
Theorem 1. For a first-order integrate-and-fire neuron as expressed in Eq. 2 (ut = kvt + λit, k ∈ (0, 1]), the model takes a constant current it = i as input and produces spikes as output. Let the discrete output firing rate w.r.t. i be expressed as f = TN , where T is the total inference time and N is the spike count. Then, the input-output response curve (f -i curve) is estimated as follows:
If k = 1 :
f =  0, i ≤ 0, λ θ i, 0 < i < θ λ ,
1, i ≥ θλ . If k ∈ (0, 1) :
f =  0, i ≤ 0, (k − 1)/ ln(1 + (k−1)θλi ), 0 < i < θ λ · k−1 ek−1−1 ,
1, i ≥ θλ · k−1 ek−1−1 .
(A1)
Proof. From Eq. 2, the iterative functions are first converted into differential equations. Here, we use the scalar variable form because we are interested in the time interval between the resting state and the firing state for a constant input. Our goal is to determine the time t∗ of firing (u(t∗) = θ)
with u(0) = 0. The discrete output firing rate f = TN can then be approximated as the inverse of the spiking interval t∗.
du dt = (k − 1)u+ λi (A2)
One can use integration to obtain the expression for u given a constant input. When k = 1, the integral is merely an accumulation of the scaled input. Considering u(0) = 0, we have the following expression:
u(t) = λi · t. (A3)
Then, with u(t∗) = θ, it is calculated that t∗ = θλi . Thus, when k = 1,
f = T N ≈ 1 t∗ = λ θ i. (A4)
For i ≥ θλ , the neuron constantly fires, and f reaches its upper bound of 1. When i ≤ 0, the neuron accumulates negative input and produce zero spikes, f = 0. When k ∈ (0, 1), Eq. A2 is an ordinary differential equation. By separating the integration variables u and t, we obtain:
du
(k − 1)u+ λi = dt.
Integrating both sides of the equation yields
ln((1− k)u− λi) = (k − 1)t+ C (1− k)u− λi = C ′e(k−1)t,
where C and C ′ are the constants of integration. Considering u(0) = 0, we have C ′ = −λi, and thus,
u(t) = λi
1− k (1− e−(1−k)t). (A5)
Then, by solving u(t∗) = θ, we obtain t∗ = 1k−1 · ln(1 + (k−1)θ λi ). Hence,
f = T N ≈ 1 t∗ = (k − 1)/ ln(1 + (k − 1)θ λi ). (A6)
When i = θλ k−1 ek−1−1 , f = 1. Thus, when i > θ λ k−1 ek−1−1 , f is clipped to 1.
Theorem 2. For a first-order integrate-and-fire neuron as expressed in Eq. 2 (ut = kvt + λit, k ∈ (0, 1]), the slant asymptote of the f -i curve satisfying f ∈ (0, 1) exists and can be uniformly expressed as follows:
fa = λ
θ i+ k − 1 2 . (A7)
Proof. From Theorem 2, the relationship between i and f satisfying f ∈ (0, 1) is given as follows:
f =
{ λ θ i, k = 1,
(k − 1)/ ln(1 + (k−1)θλi ), k ∈ (0, 1). (A8)
Suppose that the asymptote is formulated as fa = K · i+B. Then, for k = 1, K = λθ , and B = 0. For k ∈ (0, 1), if a slant asymptote exists, then limi→∞ fi → K and limi→∞[f −Ki] → B.
lim i→∞
f i = lim i→∞ k − 1 i ln(1 + θ(k−1)λi )
= lim i→∞
k−1 i
ln(λi+ θ(k − 1))− ln(λi)
= lim i→∞
−k−1i2 λ
λi+θ(k−1) − λ λi
(L’Hospital’s Rule)
= lim i→∞ λi+ θ(k − 1) θi
= λ
θ (L’Hospital’s Rule)
Thus, a slant asymptote exists for k ∈ (0, 1), and K = λθ for all k ∈ (0, 1]. For k ∈ (0, 1),
B = lim i→∞
[ k − 1
ln(1 + θ(k−1)λi ) − λ θ i
]
= lim i→∞ k − 1− λθ i[ln(λi+ θ(k − 1))− ln(λi)] ln(λi+ θ(k − 1))− ln(λi)
= lim i→∞
−λθ ln(λi+ θ(k − 1))− λ θ λi λi+θ(k−1) + λ θ ln(λi) + λ θ
λ λi+θ(k−1) − λ λi
(L’Hospital’s Rule)
= lim i→∞
λ
θ
ln(λi)− ln(λi+ θ(k − 1)) + 1− λiλi+θ(k−1) λ
λi+θ(k−1) − λ λi
= lim i→∞
λ
θ
1 i − λ(λi+θ(k−1))+λθ(k−1) (λi+θ(k−1))2
− λ2(λi+θ(k−1))2 + 1 i2
(L’Hospital’s Rule)
= lim i→∞
λ
θ θ2(k − 1)2i 2λθ(k − 1)i+ θ2(k − 1)2
= λ
θ θ(k − 1) 2λ (L’Hospital’s Rule)
= k − 1 2
When k = 1, B = 0. Thus, fa = λθ i+ k−1 2 is the uniform slant asymptote for all k ∈ (0, 1].
A.2 APPENDIX FOR EXPERIMENTAL RESULTS
The detailed hyperparameters are presented in Table A4. Note that the data in CIFAR10-DVS and N-MNIST take the form of address event representation (AER). To split the events into T slices that match the input streams of our network, we utilized the technique referenced in Appendix F of Fang et al. (2020) to preprocess the data.
Performance Analysis under Different Neuron Configurations
Eq. 2 shows that the neuron types and firing patterns are controlled by k and λ, which may vary in the biological brain. To cover as many of the neuron models used in the existing neuroscience literature as possible, we discretely sampled k from 0.1 to 1.0 and λ from 0.1 to 4.0. In each experiment, we compared four other methods with our initialization method in two aspects: training speed and accuracy. The four methods considered for comparison were proposed by Glorot & Bengio (2010), He et al. (2015), Lee et al. (2016) and Wu et al. (2018).
Table A4: Detailed hyperparameter settings in the experimental validation. The networks were designed to fit the different shapes of input data.
MNIST N-MNIST CIFAR10 CIFAR10-DVS
Input size 28× 28 T×32× 32 32× 32 T×128× 128 Timestep 20 20 20 20
Network
conv 3× 3× 32 LIF neuron pool 2× 2 conv 3× 3× 32 LIF neuron pool 2× 2
fc 500 LIF neuron
fc 100 LIF neuron
fc 10
conv 3× 3× 32 LIF neuron pool 2× 2 conv 3× 3× 32 LIF neuron pool 2× 2
fc 500 LIF neuron
fc 100 LIF neuron
fc 10
conv 3× 3× 128 LIF neuron conv 3× 3× 128 LIF neuron pool 2× 2 conv 3× 3× 256 LIF neuron conv 3× 3× 256 LIF neuron pool 2× 2 conv 3× 3× 256 LIF neuron conv 3× 3× 256 LIF neuron conv 3× 3× 256 LIF neuron pool 2× 2
fc 1024 LIF neuron
dropout fc 10
conv 3× 3× 64 LIF neuron
conv 3× 3× 128 LIF neuron pool 2× 2 conv 3× 3× 128 LIF neuron conv 3× 3× 128 LIF neuron pool 2× 2 conv 3× 3× 128 LIF neuron pool 2× 2 conv 3× 3× 128 LIF neuron pool 2× 2 conv 3× 3× 128 LIF neuron pool 2× 2
fc 1024 LIF neuron
dropout fc 1024
LIF neuron dropout
fc 10
Epoch 50 50 150 150
Batch size 32 32 32 16
LR 0.001 0.001 0.001 0.001
Optimizer Adam(0.9,0.99) Adam(0.9,0.99) Adam(0.9,0.99) Adam(0.9,0.99)
Scheduler Cosine(Tmax=50) Cosine(Tmax=50) Multistep (T=[50,90,130]) Multistep (T=[50,90,130])
Surrogate Arctan(a=2.0) Arctan(a=2.0) Arctan(a=2.0) ArcTan(a=2.0)
The exact iteration step at which the training accuracy first reached a certain expected level (60%) was recorded to represent the training speed. In Fig. A8, green bars represent our method, while bars in other colors represent other methods. It is intuitive that our method reaches the expected accuracy faster or at least at the same performance level for most neuron types compared to other methods; this is more obvious when λ falls in the range of (−∞, 0.6) or (2.0,+∞). In particular, there are no valid records at all using other methods in some cases in which λ is set to 0.1, 0.2, 8.0, and so on. However, even in these cases, the training accuracy of our method still reaches the expected level in a very short time.
In addition to our method’s superiority in terms of training speed, we measured its performance in terms of the validation accuracy to check its robustness on classification tasks under various configurations. In Fig. A9, the classification accuracies of five initialization methods under different neuron configurations are shown as heatmaps. Our method can be used to effectively and efficiently
Table A5: Optimizers and their parameters used to test the versatility.
Optimizer Learning rate Other Parameters
Adam 0.001 betas=(0.9, 0.99) SGD 0.5 momentum=0.9 ASGD 2.0 lambda=0.000001, alpha=0.75 AdaGrad 0.01 -
train SNNs with almost all neurons types, while other methods fail to train networks with some neuron types.
To further analyze whether the proposed method exhibits versatility for various choices of optimizers and surrogated funtions, eight representative groups of (k, λ, θ) were selected: (0.3, 0.1, 1.0), (0.8, 0.2, 1.0), (0.8, 0.5, 1.0), (0.2, 0.8, 1.0), (1.0, 1.0, 1.0), (0.5, 4.0, 1.0), (0.6, 0.3, 4.0) and (0.4, 0.4, 4.0).
k = 0.1
37
28 11 28 11 28 11 28 11
23
16 35
28 11
45 0
28 11
20
36 0
71 4
16 0
71 3
17
16 0
38 5
82 19
2
16 10
3
27 3
60 11
9
15 69
21 1
38 86
15 50
18 3
31 58
Ours He's Glorot's Wu's Lee's
17 36
15 2
24 43 17 32
13 7
22 3 0 20 2 6
11 6
22 2 7
31 30 41 38
22 36
86 88 11
7
73
k = 0.3
36
28 11 28 11 28 11 28 11
22
71 3
13 17
25 6
16 75
17
20 2
42 7
10 5
24 8
16
11 6
28 1
60 13
1
14 74
21 0
43 86
16 51
17 2
28 57 14 37
14 3
26 42 14 30
12 0
23 3 6 15 2 4
10 1
21 2 6 15 2 0
90
21 2 5
20
37 29 37 30
29
12 7 13
4
85 83
k = 0.5
30
28 11 28 11
13 93
28 11
20
33 5
59 1
16 1
43 6
16
15 4
29 8
60
15 9
14
88
20 9
54 90
14 56
16 0
31 61
14 38
12 9
28 4 2 14 30
10 5
23 36 16 2 3
89
22 2 6
16 18
79
19 2 4 15 15
67
17 2 2 2
2
47
11
34
22
28
17 5
16 5
12 0
84
k = 0.7
30
92 5
15 26
48 6
19 80
20
19 5
32 4
11 8
20 2
15
10 4
19 9
58 11
2
14 56
14 8
43 69 15 43
10 8
31 54
15 30
85
29 3 6
15 2 3
73
27 27 14 2 0
57
22 23 14 15
48
17 2 3
14 12
43
14 1 7 22
48
14
37 29
26
33 8
18 0
12 4 13
6
k = 0.9
31
33 8
49 7
28 2
41 1
15
12 0
18 3
11 1
15 2
15
69
11 5
60 91
15 43
81
38 56
15 30
56
31 3 7
16 22
43
29 28 15 15 30 23 23 15 13 23 19 18 15 14 21 17 17 17 12 1 5 21 16
22
79
28 55 47
30
16 4 1
97 18 0
14 7
= 0.1
k = 1.0
29
27 1 3
25
28 1 30
2
= 0.2 21
10 2
14 3
11 1 11 2
= 0.3
15
56 83 57 69
= 0.4
14
36 54 38 4 3
= 0.5
15
29 30 30 30
= 0.6
17 15
29 27 28
= 0.7
19 14
21 25 21
= 0.8
17 13 1
5 22 17
= 0.9
15 15 14 22 17
= 1.0
23 23
11
30
17
= 2.0
28
86
28 50 48
= 4.0 33
21 1 22
1
15 9
77
Figure A8: Training iterations at which the accuracy first reaches the expected level (60%) for firstorder integrate-and-fire neurons with different (k, λ) settings and initialization methods. Note that θ is set to 1.0. The total number of iterations of 3 training epochs is 2811. A value of 2811 in this figure means that in the corresponding setting, the network cannot be trained to reach the expected training accuracy level.
Versatility Across Optimizers
During model training, different optimizers may be used to minimize the loss function. To test the versatility of our method across various optimizers, we chose four classic and prevalent optimizers for model training, i.e., Adam, SGD, ASGD, and AdaGrad. The detailed optimizer parameters are listed in Table A5.
We present the accuracy of training in Table A7. In Fig. A10, we show the training curves during the first 500 iterations to illustrate the results in more detail. As shown in this figure, when SGD or ASGD is chosen to optimize the network, our method demonstrates an advantage in accelerating the training in the start-up phase with the accuracy reaching a good level within 100 iterations in most cases. For example, when k = 0.8, λ = 0.5, and θ = 1.0, Wu’s and He’s methods allow training to begin early but with a slower speed, achieving an accuracy of approximately 90% after 500 iterations, while Glorot’s and Lee’s methods show much worse performance. When Adam or AdaGrad is chosen to optimize the network, despite the similar performance of all methods in some cases (for example the 3rd-5th rows in the figure), our method outperforms the others in accelerating
Figure A9: Training accuracy heatmaps for first-order integrate-and-fire neurons with different (k, λ, θ) settings and initialization methods. The color represents the validation accuracy in the given setting.
the training process. In addition, our weight initialization for SNNs enables valid training in some difficult cases in which training totally fails with other initialization methods, such as (k = 0.4, λ = 0.4, θ = 0.4).
Versatility Across Surrogate Functions
Various surrogate functions may be chosen for SNN training w.r.t. different needs and tasks. As listed in Table A3, four surrogate functions (sigmoid, tanh, hardtanh, and arctan) are popularly used in BPTT. As described in Sec. 2.1, the input is subjected to the Heaviside function to perform the forward pass, while the backward gradient is computed using surrogate functions. We trained the 5-layer SNN using surrogate functions consistent with those in Table A3 by adjusting the parameter a.
The validation accuracy after 3 epochs of training is presented in Table A7, and the training curves from the first 220 iterations with the four surrogate functions are shown in Fig. A11. The curves in green show our method’s effectiveness. Our method reaches acceptable accuracy within relativelt few iterations (fewer than 50 in most cases), while in comparison, the decline in the network loss starts more slowly under other initialization methods, taking parameter sets of (0.8, 0.2, 1.0), (0.8, 0.5, 1.0), (0.2, 0.8, 1.0) and (0.5, 4.0, 1.0) as examples. Especially when k and λ have small values or θ has a large value, other initialization methods cannot generate effective gradients for BP, leading to training failure. For almost all sets of surrogate functions and neuronal parameters, our method enables neuron firing and early generation of network gradients. Fig. A11 reveals that our method enables models to start training faster and to end training with considerably greater classification accuracy.
A.3 COMPARISON TO BATCH NORMALIZATION (BN)
As Eq. 11 implies, the variances of input current from the adjacent layers have been scale to be identical when using our initialization method, which can be viewed as a kind of normalization of current before the start of training. The difference lies in the term of validity. Normalization methods such as Kim & Panda (2020), Kim & Panda (2021) are normalizing current throughout training. To clarify this, we examined the validation accuracy with or without normalization methods on CIFAR100 and TinyImageNet datasets. We chose vanilla BN (Ioffe & Szegedy, 2015) and BNTT
0.0 0.3 0.6 0.9 A cc
Adam SGD ASGD
k = 0.3 = 0.1 = 1.0
AdaGrad
0.0 0.3 0.6 0.9 A cc k = 0.8 = 0.2 = 1.0
0.0 0.3 0.6 0.9 A cc k = 0.8 = 0.5 = 1.0
0.0 0.3 0.6 0.9 A cc k = 0.2 = 0.8 = 1.0
0.0 0.3 0.6 0.9 A cc k = 1.0 = 1.0 = 1.0
0.0 0.3 0.6 0.9 A cc k = 0.5 = 4.0 = 1.0
0.0 0.3 0.6 0.9 A cc k = 0.6 = 0.3 = 4.0
0 50 100 150 200 250 300 350 400 450 Iterations
0.0 0.3 0.6 0.9 A cc
0 50 100 150 200 250 300 350 400 450 Iterations
Glorot's He's Wu's Lee's Ours
0 50 100 150 200 250 300 350 400 450 Iterations 0 50 100 150 200 250 300 350 400 450 Iterations
k = 0.4 = 0.4 = 4.0
Figure A10: Training accuracy curves obtained using four different optimizers, different (k, λ, θ) settings and different initialization methods.
(batch normalization through time) (Kim & Panda, 2020) for comparison. When the network is without BN, a dropout is implemented in place of BN. (k = 0.99, λ = 1.0) is the default setting inKim & Panda (2020). The validation accuracy curve is plotted in Fig. A12. The use of BNTT make the effect of initialization insignificant as it works throughout training. Table A8 shows that with BN
Table A6: Validation accuracy achieved when training on MNIST for 3 epochs using four different optimizers, different (k, λ, θ) settings and different initialization methods. A hyphen in this table indicates invalid training in a given setting.
k = 0.3, λ = 0.1, θ = 1.0 k = 0.8, λ = 0.2, θ = 1.0
Ours Glorot’s He’s Wu’s Lee’s Ours Glorot’s He’s Wu’s Lee’s
Adam 0.979 - - - - 0.982 0.972 0.972 0.973 0.971
SGD 0.776 - - - - 0.929 - - - -
ASGD 0.636 - - - - 0.887 - - - -
AdaGrad 0.966 0.878 0.898 0.876 0.740 0.975 0.959 0.967 0.956 0.957
k = 0.8, λ = 0.5, θ = 1.0 k = 0.2, λ = 0.8, θ = 1.0
Ours Glorot’s He’s Wu’s Lee’s Ours Glorot’s He’s Wu’s Lee’s
Adam 0.984 0.980 0.986 0.983 0.984 0.984 0.978 0.988 0.982 0.989
SGD 0.975 0.921 0.959 0.964 0.945 0.968 - 0.950 0.956 0.884
ASGD 0.958 0.165 0.910 0.918 0.866 0.951 - 0.905 0.926 0.612
AdaGrad 0.984 0.978 0.984 0.984 0.982 0.977 0.978 0.986 0.985 0.986
k = 1.0, λ = 1.0, θ = 1.0 k = 0.5, λ = 4.0, θ = 1.0
Ours Glorot’s He’s Wu’s Lee’s Ours Glorot’s He’s Wu’s Lee’s
Adam 0.987 0.986 0.987 0.986 0.986 0.986 0.987 0.984 0.985 0.986
SGD 0.985 0.982 0.985 0.983 0.985 - - - - -
ASGD 0.986 0.981 0.986 0.981 0.984 0.981 - - - -
AdaGrad 0.986 0.984 0.986 0.985 0.986 0.983 - - 0.977 0.984
k = 0.6, λ = 0.3, θ = 4.0 k = 0.4, λ = 0.4, θ = 4.0
Ours Glorot’s He’s Wu’s Lee’s Ours Glorot’s He’s Wu’s Lee’s
Adam 0.980 - - - - 0.984 - - - -
SGD 0.895 - - - - 0.923 - - - -
ASGD 0.841 - - - - 0.876 - - - -
AdaGrad 0.968 0.882 0.889 0.910 0.743 0.972 0.906 0.924 0.926 0.752
or BNTT, the gaps of the accuracy of different initialization methods are not obvious. But when the network is without BN, our method can also work under (0.5, 0.5) configuration. Therefore, we are not stating that our initialization method can outperform some well-designed learning algorithms. Instead, we believe our method is a useful supplement to make SNNs trainable on larger datasets.
Table A7: Validation accuracy achieved when training on MNIST for 3 epochs using four different surrogate functions, different (k, λ, θ) settings and different initialization methods. A hyphen in this table indicates invalid training in a given setting.
k = 0.3, λ = 0.1, θ = 1.0 k = 0.8, λ = 0.2, θ = 1.0
Ours Glorot’s He’s Wu’s Lee’s Ours Glorot’s He’s Wu’s Lee’s
Hardtanh(a=1.0) 0.982 - - 0.310 - 0.985 0.974 0.980 0.980 0.977
Arctan(a=1.0) 0.983 - - - - 0.984 0.972 0.980 0.980 0.976
Tanh(a=2.0) 0.984 - - - - 0.986 0.976 0.979 0.981 0.978
Sigmoid(a=2.0) 0.984 - - - - 0.985 - - - -
k = 0.8, λ = 0.5, θ = 1.0 k = 0.2, λ = 0.8, θ = 1.0
Ours Glorot’s He’s Wu’s Lee’s Ours Glorot’s He’s Wu’s Lee’s
Hardtanh(a=1.0) 0.987 0.985 0.987 0.984 0.986 0.987 0.983 0.990 0.986 0.988
Arctan(a=1.0) 0.988 0.983 0.986 0.987 0.989 0.988 0.986 0.991 0.987 0.991
Tanh(a=2.0) 0.989 0.987 0.987 0.988 0.990 0.988 0.984 0.990 0.989 0.991
Sigmoid(a=2.0) 0.988 - - - - 0.986 - - - -
k = 1.0, λ = 1.0, θ = 1.0 k = 0.5, λ = 4.0, θ = 1.0
Ours Glorot’s He’s Wu’s Lee’s Ours Glorot’s He’s Wu’s Lee’s
Hardtanh(a=1.0) 0.988 0.988 0.988 0.984 0.989 0.985 0.987 0.980 0.984 0.986
Arctan(a=1.0) 0.986 0.987 0.986 0.981 0.984 0.987 0.987 0.973 0.973 0.978
Tanh(a=2.0) 0.988 0.988 0.988 0.990 0.987 0.988 0.983 0.985 0.984 0.980
Sigmoid(a=2.0) 0.987 0.986 0.987 0.887 0.989 0.199 0.197 - - 0.176
k = 0.6, λ = 0.3, θ = 4.0 k = 0.4, λ = 0.4, θ = 4.0
Ours Glorot’s He’s Wu’s Lee’s Ours Glorot’s He’s Wu’s Lee’s
Hardtanh(a=1.0) 0.979 - - - - 0.983 - - - -
Arctan(a=1.0) 0.917 - - - - 0.953 - - - -
Tanh(a=2.0) 0.968 - - - - 0.976 - - - -
Sigmoid(a=2.0) 0.887 - - - - 0.951 - - - -
Table A8: Accuracy comparison with different normalization technique training with VGG11 on large datasets. A hyphen in this table indicates invalid training in a given setting.
Dataset (k, λ) Technique Initialization method Glorot’s He’s Wu’s Lee’s Ours
CIFAR100 (0.5, 0.5) w/o BN, Dropout - - - - 54.94 CIFAR100 (0.99, 1.0) w/o BN, Dropout - 54.03 54.46 - 53.79 CIFAR100 (0.99, 1.0) w/ BN 57.18 57.19 56.45 56.64 57.02 CIFAR100 (0.99, 1.0) w/ BNTT 65.89 65.55 64.91 65.65 65.44 TinyImageNet (0.5, 0.5) w/o BN - - - - 44.16 TinyImageNet (0.99, 1.0) w/ BN 45.97 47.18 47.60 45.03 47.51 TinyImageNet (0.99, 1.0) w/ BNTT 57.03 57.09 56.90 57.19 56.80

Review summaries:
The paper investigates the effects of weight initialization for deep spiking neural networks (SNNs) on training speed and final accuracy. The main finding is that popular initialization methods for conventional ANNs and RNNs do not match the specific dynamics of SNNs, in particular they ignore the need to have sufficient firing early on to generate good gradients. The paper derives a theoretical first-order approximation on the response curve, and from it a weight initialization scheme that initializes the weights in a region where neurons fire from the first epoch on, but also avoid explosion or reduction of activities throughout the network. In various experiments on standard and spiking variants of MNIST and CIFAR, the initialization scheme shows good performance and fast convergence for different neuron types (e.g. non-leaking vs. fast leaking), and different encoding schemes. For CIFAR10 the results are significantly better than for other initialization schemes.

 This paper focuses on the parameter initialization problem of training SNNs. The authors derive the theoretical response of spiking neurons, and propose an initialization method based on slant asymptote, which can overcome the gradient vanishing problem. The results show that the proposed method can effectively improve training speed and accuracy.

The authors tackle the problem of improving spiking neural net training with better initialization. They observe that error backpropagation efficiency in SNNs depends on availability of responding neurons and show that traditional surrogate backprop methods are inefficient as they take many iterations to move weights to the regime supporting neuronal responses.  The authors approximate the SNN neuron model (LIF, IF etc) I/O function using a piecewise-linear  iterative expression and show that clipping the output to (0,1) and normalizing the variance of the random weight initialization to insure that neurons produce substantial responses even with initial random weights results improves SNN training substantially. The improvement is to a large extent due to early onset of convergence. The authors demonstrate the superiority of their approach on MNIST, N-MNIST, DVS-MNIST and CIFAT datasets and present extensive experiments using different optimization algorithms and response functions.  

This paper enhances the efficient BPTT training of SNNs by proposing an initialization method to match the response of spiking neurons in initial training. Their method bridges the spiking neuron response to the wisdom of traditional deep learning training, which may have an influence on future research like ANN-to-SNN conversions with LIF neurons or other SNN training methods. The authors conduct experiments on CIFAR10, MNIST, and neuromorphic datasets to show the efficacy of their technqiue. One of the main contribution is their theorteical analysis of iterative systems to model first-order integrate-and-fire neurons and investigate the response curve.