Differentiable Spike: Rethinking Gradient-Descent for Training Spiking Neural Networks
Year: 2021

Abstract:
Spiking neural networks (SNNs) have emerged as a biology-inspired method mimicking the spiking nature of brain neurons. This biomimicry derives SNNs’ energy efficiency of inference on neuromorphic hardware. However, it also causes an intrinsic disadvantage in training high-performing SNNs from scratch since the discrete spike prohibits the gradient calculation. To overcome this issue, the surrogate gradient (SG) approach has been proposed as a continuous relaxation. Yet the heuristic choice of SG leaves it vacant how the SG benefits the SNN training. In this work, we first theoretically study the gradient descent problem in SNN training and introduce finite difference gradient to quantitatively analyze the training behavior of SNN. Based on the introduced finite difference gradient, we propose a new family of Differentiable Spike (Dspike) functions that can adaptively evolve during training to find the optimal shape and smoothness for gradient estimation. Extensive experiments over several popular network structures show that training SNN with Dspike consistently outperforms the state-of-the-art training methods. For example, on the CIFAR10-DVS classification task, we can train a spiking ResNet-18 and achieve 75.4% top-1 accuracy with 10 time steps.

heading: 1 Introduction
text: Recently, spiking neural networks (SNNs) have received increasing attention due to their biologyinspired neural behavior and efficient computation. SNNs deal with binary spike information and therefore enjoy the advantage of multiplication-free inference. Neuromorphic hardware such as TrueNorth [1] and Loihi [2] demonstrates that SNNs can save energy by orders of magnitude. Hybrid architecture like Tianjic [3] suggests its potential power for general intelligence when integrated with traditional artificial infrastructure. However, the bio-mimicry also causes an intrinsic disadvantage in training high-performing SNNs from scratch due to the discrete spike. As a result, obtaining a high-performing SNN has long been a critical problem limiting SNN’s deployment in practice.
To be more specific, the appealing advantage of the multiplication-free computation in SNN is not recognized in its training [4]. Modern ML frameworks such as Pytorch [5], JAX [6], and Tensorflow [7] do not provide efficient and general instructions or high-level functions to accelerate the convolution with 0/1 spike activation. Moreover, the binary activation in SNN produces allor-nothing gradients [8] that are incompatible with typical gradient-based optimization methods for efficient neural network training. To circumvent this problem, many works adopt surrogate gradients [8] (typically the derivative of a soft-relaxed function) to replace the spike non-linear gradient. Although the heuristic choice of SG is beneficial to enable the gradient descent in SNN training, the theoretical soundness of SG has not been justified so far. What’s more, the gradient descent in SNN training may follow different updating rules from ANN. In fact, since the true gradient
∗Equal Contributions. Corresponding author. Primarily supported by NSF-61876032.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
value of spike activation does not exist, the adoption of SG raises a new question of how to replace the∞-gradient adaptively through the training process for a good updating strategy of parameters. In this work, we study this specialized problem of adaptive gradient estimation in SNN training. We first revisit the gradient descent algorithm with the finite difference method to accurately depict the loss landscape of adopting a surrogate gradient for the non-differentiable spike firing function. Then we further propose the Differentiable Spike (Dspike) function, which can adaptively change its shape and capture the direction of finite difference gradients. With such adaptive surrogate gradients, we establish an efficient algorithm of training SNN with large-scale models. The overall workflow of our algorithm is visualized in Fig. 1.
Our contributions can be summarized as follows:
• We identify the problem of adaptively estimating gradient in SNN training and propose a novel analysis framework for SNN’s surrogate gradient, through which we are able to not only explain why it works well but also indicate how we can improve its adaptivity to data through training.
• We also introduce the novel Dspike function built from the hyperbolic tangent function. The Dspike function can learn its shape accordingly by matching it with the finite difference gradient.
• Extensive experiments on both static and dynamic datasets show that our method is highly effective and efficient. For example, our model can achieve 75.4% top-1 accuracy on the CIFAR-DVS. Our SNN also only consumes 7.8% energy of the corresponding ANN on CIFAR-10.

heading: 2 Related Work
text: ANN-SNN Conversion. Generally, the algorithms to obtain a high-performance can be divided into two routes: (1) ANN-SNN conversion and (2) direct training SNN from scratch. Conversion-based methods [9, 10] utilize the knowledge from ANN and convert the ReLU activation to spike activation mechanism. This type of method can produce an SNN in a short time. For example, Rueckauer et al. [11] propose to traverse the training dataset to find the maximum activation. Then, the SNN can be converted by directly replacing several modules. In [12], maximum activation is abandoned since it may out-lie the central range of activation. The authors choose to find the percentile numbers for threshold voltage to improve the robustness. TSC [13] adopts the soft-reset (also known as reset by subtraction) mechanism and no leakage for conversion, which proves to be effective in conversion. Deng & Gu [14] decompose the conversion error to each layer, then propose to reduce the error by shifting the bias term. [15] first propose the calibration to reduce the conversion accuracy. ANN-SNN conversion is fast and accurate. However, achieving near lossless conversion requires a considerable amount of time steps (>200) to accumulate the spikes, which may significantly add up the latency.
Direct Training SNN. In contrast to the conversion method, direct training of SNNs can lower the time step and may further reduce energy consumption. Direct training leverages gradient-descent algorithms. And many works have been developed to tackle the non-differentiability of spike activation. This non-differentiable activation function is similar to the quantization or binarization
neural network [16, 17, 18] or adder neural network [19]. [20] achieve high performance on CIFAR10 and other neuromorphic datasets. Some works use probabilistic models to approximate the gradient [21], which is also developed in binary stochastic model [22]. Another common approach is to use rate-encoding networks for gradients calculation [23, 24]. They quantize the spike rate as the magnitude of information transmitted. (Wu et al. 2019 [25] improves the leaky integrate-and-fire (LIF) model to an iterative LIF model and develops STBP learning algorithm. Zheng et al. [26] further propose the tdBN algorithm to balance the gradient norm and smooth the loss function in SNN training, enabling large-scale model optimization.

heading: 3 Preliminary

heading: 3.1 Leaky Integrate-and-Fire Model
text: We use bold italic letters for vector representations. For example, x and y denote the input and target output variables. Bold capital letters like W denote the matrices (or tensors as clear from the text). Constants are denoted by small upright letters, e.g., a. With bracketed superscript and subscript, we can denote the time dimension and the element indices, respectively. For example, u(t)i means the i-th membrane potential at time step t. We adopt the well-known Leaky Integrate-and-Fire (LIF) neuron model for a spiking neural network (SNN). Given a membrane potential u(t), the pre-synaptic potential is updated as
u(t+1),pre = τu(t) + c(t+1),where c(t+1) = Wx(t+1). (1)
Here, τ is a constant leaky factor within (0, 1) and c(t+1) is the pre-synaptic input at time step t+ 1. The input can be calculated from the dot-product between the weights and the spike from the previous layer. Let Vth be the firing threshold, the membrane potential will fire a spike when it exceeds the threshold, and then hard-reset to 0, given by
y(t+1) = { 1 if u(t+1),pre > Vth 0 otherwise , u(t+1) = u(t+1),pre · (1− y(t+1)) (2)
After firing, the spike output y(t+1) will propagate to the next layer and become the input x(t+1) of the next layer, (note that we omit the layer index for simplicity). Vth is set to 0.5 in this work.
Softmax in SNN. In the general classification model, the final network output will be leveraged to compute the softmax and then to find the winner class. However, the network output is highly variant per input sample, and unlike ReLU activation, the output could be positive or negative. For this reason, we choose to only integrate the network output, and do not fire them across time, as did in recent practices [26, 27].
ynet = 1
T T∑ i=1 c (t) net = 1 T T∑ i=1 Wx(t), (3)
Then, we can compute the cross-entropy loss based on the ture label and Softmax(ynet).

heading: 3.2 Time-dependent Batch Normalization
text: Batch Normalization (BN) [28] is used to reduce internal covariate shift and smooth the loss landscape [29] during the optimization of an ANN. It enables the stable training of ANN even with more than 100 layers on the large-scale dataset, e.g. ResNet-152 [30]. However, BN layers are not designed to normalize spatial-temporal data at first. And simply transplanting BN in SNN training will lead to undesired results. To solve this problem, Zheng et al. [26] propose time-dependent Batch Normalization (tdBN), which normalizes the data in both spatial and temporal paradigm. tdBN will gather the pre-synaptic input along the time dimension, denoted by C = [c(1), c(2), . . . , c(T )], and then normalize it by
C̄ = αVth(C− E[C])√
VAR(C) + # normalize (4)
Ĉ = γC̄ + β, # affine transformation (5)
where E and VAR compute the mean and variance in channel dimension. γ and β are learnable parameters. Note that in tdBN the normalized tensor will multiply αVth, which can prevent over-fire and under-fire by choosing α wisely. After tdBN, the pre-synaptic input in each time step will be fed into Eq. (1). tdBN stabilizes the training of SNN and achieves good results. We also adopt this training technique in our framework. See [26] for more details of tdBN.

heading: 4 Methodology

heading: 4.1 Non-Differentiable Spike
text: One of the most notorious problems in SNN training is the non-differentiability of the firing function Eq. (2). The firing function could be viewed as a variant of sign function, whose derivative is identified as the Dirac delta function. This derivative is 0 almost everywhere except for the threshold that grows to infinity. Seemingly, directly using the Dirac delta function is not an ideal choice for gradient descent-based learning. To demonstrate this problem, let L denote the task loss objective, we can compute the gradient of weights parameters using the spatial-temporal backpropagation (STBP):
∂L ∂W = ∑ t ∂L ∂y(t) ∂y(t) ∂u(t),pre ∂u(t),pre ∂c(t) ∂c(t) ∂W . (6)
In above equation, ∂y (t)
∂u(t),pre is the gradient of firing function and is 0 almost everywhere, as aforemen-
tioned. As a consequence, the gradient descent (W←W− η ∂L∂W ) either freezes weights or updates the weights to infinity.
To mitigate this problem, prior work seeks an alternative gradient (also known as surrogate gradient) to circumvent this problem. A popular surrogate gradient may refer to the rectangular function proposed in [4], given by
∂y(t)
∂u(t),pre =
1 a sign (∣∣∣u(t),pre − Vth∣∣∣ < a 2 ) , (7)
where a is a hyper-parameter and usually set to 1. When Vth = 0.5, a = 1, the gradient is simplified to 1 if 0 ≤ u(t),pre ≤ 1 otherwise 0. This is known as the Straight-Through Estimator in Bengio et al. [22] and is broadly applied in binary networks [31]. Another family of surrogate gradient may refer to triangular gradients [32], given by
∂y(t)
∂u(t),pre = γmax
( 0, 1− ∣∣∣∣u(t),preVth − 1 ∣∣∣∣) . (8)
This surrogate gradient is computed based on the distance between the potential and threshold. Despite various surrogate gradients achieve acceptable accuracy, there lacks a theoretical justification for its reason and effectiveness.

heading: 4.2 Rethinking Gradient-Descent in Spiking Neural Networks
text: Gradient-descent (GD) algorithm is a first-order iterative optimization method. Theoretically, for each iteration, GD aims to find the steepest direction for the descent that minimizes the task loss:
min ||∆W||≤ε L(W + ∆W) ≈ min ||∆W||≤ε
L(W) + ∆W ∂L
∂W (9)
Note that the right-hand side of the minimization is approximated by the first-order Taylor expansion. Given that L(W) is constant at the current iteration, therefore we only have to find the right ∆W that minimizes the second term. Since the 2nd term is a linear in ∆W, it can be made as negative as we like by taking ∆W large. Therefore the minimum is ∆W = −ε ∇WL||∇WL|| , (where∇WL = ∂L ∂W ).
Note that the above condition is only sufficient for ε → 0 to correctly approximate the loss with Taylor expansion. However, ε → 0 means ∆W = −ε ∇WL||∇WL|| → 0, but in practice GD updates the weights by a small step (determined by learning rate), which, indicates a mismatch between the practical implementation and the theoretical assumption. In other words, gradient only inspects
the change rate of loss around an infinitely small neighborhood, however, we are interested in the change rate around a tiny (but not infinitely small) area. In the case of SNN, the gradient will always be 02, yet for a tiny neighborhood the loss function will change, leading to non-zero gradients. To penetrate the non-differentiability of the firing function, we adopt a technique called finite difference. It estimates the gradient by evaluating the loss objective near the coordinate:
∇̂ε,wL = ( L(w1 + εe1)− L(w1 − εe1)
2ε , . . . , L(wk + εek)− L(wk − εek) 2ε
) , (10)
where w is the flattened vector of W and ei refers to standard basis vector with 1 at i-th index. ε is some small constant, e.g. 0.001. We call ∇̂ε,wL as finite difference gradient (FDG). By the definition of gradient, we know FDG becomes the true gradient when ε → 0. FDG can be a useful tool in GD-based learning, however, evaluating finite difference for each element in the weight vector is tedious, which greatly slows down the training process. Nevertheless, FDG still offers us a choice to evaluate the gradient mismatch problem quantitatively.
In order to illustrate this, we construct a toy problem: a simple 2-layer CNN on the MNIST dataset. The first layer is a 3x3 convolutional layer (output channel number is 6) with ReLU activation (or LIF spike activation) and average pooling for down-sampling. Then, we directly add a fully connected layer to predict the classification results. We train the toy model on MNIST for 5 epochs, with batch size 128 and learning rate 0.01 decayed to 0. In an effort to compute the FDG of the first layer, we have to perform 2× 1× 6× 3× 3 = 108 times forward pass 3. In contrast, one has to perform 22 million times forward passes on a ResNet-18 [30] to compute the FDG of all parameters.
Experiment I: FDG in ANN. We first train an ANN on MNIST and compare the disparity between real gradients and FDG. We use cosine similarity (i.e. cos_sim(x1,x2) = x1x2||x1||2||x2||2 ) to determine the angle between these two gradient vectors. For each epoch, we compute the mean cosine similarity in the first 100 training iterations. The results are demonstrated in Fig. 2(a). The FDG and real gradients consistently have an extremely similar direction (> 0.99 cos_sim) with ε varying from 0.1 to 0.001. In conclusion, we show that (1) FDG is a valid signal and (2) the real gradients can effectively find the change rate within a small neighborhood in ANN optimization.
Experiment II: FDG in SNN. Then, we compare the cosine similarity between the surrogate gradients in SNN and the FDG. In Fig. 2(b), we can tell that surrogate gradients share similar gradient direction with FDG when ε = 0.1 (i.e. more than 0.8 cosine similarity). However, decreasing ε will lead to lower cosine similarity. In the case of ε = 0.001, the angle between FDG and the surrogate gradients is nearly 90 degrees, (for comparison, in high-dimensional space, two arbitrary random Gaussian vectors are very likely to be orthogonal). This observation concludes that surrogate gradient can predict the change rate of loss in a relatively larger neighborhood, which justifies the soundness of surrogate gradients. However, SG fails to predict the FDG in small step size ε. To this end, we need to figure out whether a large step size is better or a small step size is better.
Experiment III: Training SNN with FDG. In SNN optimization, a natural question arises: What ε is optimal for gradient descent? First, the ε has to be small enough to compute the non-linearity
2Assuming the potential is not equivalent to Vth. 3According to Eq. (10), for each element we evaluate the loss twice, the first layer has 54 weight elements.
of the loss surface, meanwhile, it has to be large enough to ignore the noise from the rough step function. For this reason, we explore the optimal ε by applying FDG in SNN optimization. The results are demonstrated in the Table 1. Note that all methods are trained on MNIST for only 10 epochs here. We show that the rectangular surrogate gradient with a = 1 reaches 94.9% accuracy. For FDG optimization, we evaluate ε ranging from 0.001 to 0.5, and we find 0.1 is the best choice for optimizing SNN. The results confirm our assumption that ε has to be neither too big nor too small. It also tells us using the instant change rate (very small ε) will lead to undesired performance degradation.

heading: 4.3 Differentiable Spike
text: In this section, we introduce the Dspike function, a novel surrogate gradient estimator, which can cover a large range of choices for SG. We choose the hyperbolic tangent function tanh(x) = e
x−e−x ex+e−x
as our basic building function for gradient estimation:
Dspike(x) =  1 if x > 1 a · tanh(b(x− c)) + d if 0 ≤ x ≤ 1 0 if x < 0
(11)
When x < 0 or x > 1, we choose to clip x as did in other practice. The vanilla tanh function maps (−∞,+∞) to (−1, 1), thus we have four parameters to control its shape to alter the function mapping: from [−1, 1] to [1, 1]. We first use c = Vth = 0.5 to set the symmetry center to (0.5, 0.5). Parameter b(b > 0) is called the temperature factor, it controls the shape of the Dspike function. Parameter a and d can then be determined by setting Dspike(0) = 0,Dspike(1) = 1. As a result, we can rewrite the function as:
Dspike(x, b) = tanh(b(x− 0.5)) + tanh(b/2)
2(tanh(b/2)) , if 0 ≤ x ≤ 1. (12)
Dspike function has the capability to adjust the shape of the function, leading to different surrogate gradient estimations. To demonstrate that, we have the following lemma: Lemma 4.1. By changing the temperature b, the shape of Dspike function can be altered. Specifically, we have limb→0+ Dspike(x)→ x and limb→+∞Dspike(x)→ sign(x− 0.5).
Proof. We first prove the case b→ 0+. In this case, applying Taylor expansion to the tanh we have:
tanh(bx) ≈ bx− b 3x3 3 + · · ·+ 2 (2n)(2(2n) − 1)B2n(bx)2n−1 (2n)! , (13)
where Bn is the n-th Bernoulli number. When b→ 0, we can safely ignore higher-order terms and only keep bx. Thus, we can rewrite Eq. (12) to Dspike(x, b) = bx/b = x. Then, we prove the case b→ +∞. Since limx→+∞ tanh(x)→ 1, we have:
lim b→+∞
Dspike(x, b)←  2∗tanh(+∞) 2∗tanh(+∞) = 1 if x > 0.5 tanh(+∞) 2∗tanh(+∞) = 0.5 if x = 0.5 0
2∗tanh(+∞) = 0 if x < 0.5 (14)
With this temperature factor, our algorithm can simulate a wide range of surrogate gradients with different smoothness. We visualize the transition in Fig. 3 from both forward perspective and backward perspective. As aforementioned, surrogate gradient can mimic the FDG in a relatively large neighborhood. Our intuition is that the optimal surrogate gradient can be found by choosing temperature wisely. Thus our task becomes to optimize temperature and network parameters simultaneously during training. As did in former toy experiments (Sec. 4.2), we verify the cosine similarity between our Dspike gradients and the FDG at the beginning of each training epoch and update the temperature to better match the FDG. This algorithm is detailed in Algo. 1.
As abovementioned, evaluating FDG could be time-consuming, due to the two times evaluation of the final loss for each weight element. Therefore, we choose to only compute the FDG in the first layer for two reasons. First, the first layer has much fewer weights. In ResNet-18 (for CIFAR10 dataset), the first layer only has 3× 3× 3× 64 = 576 elements. Second, the gradient propagation may accumulate the error in the former layers, which effectively measures the gradient mismatch.
Algorithm 1: Training SNN with Dspike gradient Input: SNN to be trained; Training dataset, total training epoch E, total training iteration in one epoch I , FDG step size ε, temperature update step size ∆b, initialized temperature b0. for all e = 1, 2, . . . , E-th epoch do
Collect input data and labels, computing the FDG (∇̂ε,wL) of the first layer using ε; for all b in {be−1, be−1 + ∆b, be−1 −∆b} do
Compute Dsipke surrogate gradient∇b,wL and cos_sim(∇̂ε,wL,∇b,wL) Find the optimal temperature with the highest cosine similarity and update it to be; for all i = 1, 2, . . . , I-iteration do
Get training data and labels, compute Dspike surrogate gradient using optimal b; Descend loss function and update weights;
return trained SNN.

heading: 4.4 Time Inheritance Training
text: We also introduce a novel training pipeline that can effectively reduce the time for training multiple SNNs. Unlike ANN-SNN conversion, direct training SNN is performed with a fixed time step T . As a result, we must learn a new SNN from scratch if we decide to use another different T . We propose Time Inheritance Training (TIT) pipeline to ameliorate this problem. First, we train an SNN to converge with high T , e.g. 20. Then, we use this well-trained model to initialize an SNN with lower T , say T = 10. Consequently, the initialized model has exactly the same output as the former model in the first 10 time steps, enabling better initialization and faster convergence. In other words, TIT progressively reduces the time step finetuning from a higher time step SNN first. This idea resembles the curriculum learning proposed in Bengio et al. [33], which suggests learning the easy task first. Armed with TIT, we can reduce the training time of subsequent models by 75%, while sharing the same task performances.

heading: 5 Experiments
text: In this section, we report experimental results for our Dspike algorithm on CIFAR [34] and ImageNet dataset [35]. For static image dataset, we encode the images to spike using the first layer in the SNN, as adopted in recent works [27]. For the neuromorphic image dataset, the images are already converted to 0/1 spike format. We also replace the max-pooling layer with the average-pooling layer in the network architecture. The detailed implementation is described in each subsection.

heading: 5.1 Network Architecture Choice
text: For CIFAR series experiments. We mainly adopt the ResNet-family architectures [30]. Unlike traditional ResNet structure [30], we use a modified version with 3 stem layers and avg-downsample
block, as introduced in [36]. Our ResNet-18 architectures contains 4 stages with 2 basic residual block in each stage. Our initial channel is set to 64 and doubled in each stage. In comparison with [26], they use ResNet-19 with 3 stages, with 3 blocks in the first 2 stages and 2 blocks in the last stage. However, they choose 128 channel and they has larger receptive field than our model. In general, their ResNet-19 has ∼10 times more operations than our ResNet-18. For ImageNet, we use the same VGG-16 as used in [27] and the same ResNet-34 as used in [26].

heading: 5.2 CIFAR
text: The CIFAR10 (100) dataset contains 50K training and 10K test images with 32×32 pixels. We use AutoAugment [37] and Cutout [38] for data augmentation. We use ResNet-18 architecture for running experiments. We train the corresponding ANN first and use it to initialize the first SNN with (T = 6), then we use Time Inheritance Training to gradually reduce the time step to 2. For FDG computation, we use the best practice ε = 0.1 in the toy experiments. The temperature is initialized to 1, and we set the update step size to ∆b = 0.2. We adopt SGD optimizer with 0.9 momentum. In the first round of TIT, we train the model for 300 epochs with a learning rate of 0.01 cosine decayed [39] to 0. In the rest rounds of TIT, we only train the model for 50 epochs with a learning rate of 0.004. Weight decay (L2 regularization) is set to 0.0001. We run the model with 4 GTX 1080Tis. For each run, we report the mean accuracy as well as the standard deviation with 3 trials.
We summarize the overall results in Table 2. On CIFAR10, our work can outperform state-of-the-art by 1.09% accuracy when T = 6. It is also worth noting that the ResNet-19 used in [26] has 2x channel numbers and a larger feature-map, translating to 10x computation cost. Nevertheless, our ResNet-18 with Dspike still achieves higher accuracy. Notably, our ResNet-18 (T = 2, 4) is trained with TIT for only 50 epochs, yet still reaches a new state-of-the-art. We also conduct an ablation study: training a 2 time step spiking ResNet-18 without TIT for 300 epochs. The accuracy is only 93.13%, which is subpar with the TIT method. We also verify Dspike on CIFAR100, with the same training pipeline here and we find our model can achieve much better performances than Diet-SNN [27].
To further demonstrate the effectiveness of our method, we apply the Dspike function on the eventstream dataset CIFAR10-DVS [45]. CIFAR10-DVS contains 10k 128×128 images converted from CIFAR10. We split the dataset into 9k training images and 1k test images following the prior practice [26]. For data pre-processing and augmentation, we first resize the training image frames to 48×48, then apply random horizontal flip. Additionally, we apply random roll within 5 pixels. The test images are directly resized to 48×48. Training hyper-parameters are kept the same with
the CIFAR10 experiments. Our proposed method can achieve 75.4% averaged top-1 accuracy on CIFAR10-DVS, with 7.6% absolute accuracy improvement compared with the existing state of the art.

heading: 5.3 ImageNet
text: In this section we verify our algorithm on the ImageNet dataset. ImageNet dataset has more than 1250k training images and 50k test images. We use the standard pre-processing and augmentation for training data [30]. The test data is directly centered cropped to 224×224. We train ResNet-34 and VGG-16 on ImageNet using the same architecture in the CIFAR experiments. Each network is trained for 100 epochs with SGD optimizer with 0.9 momentum. The learning rate is set to 0.04 and weight decay is set to 0.00004. We train each model on 32 GTX 1080Tis with batch size of 10 per GPU. Following [26], to avoid the biased statistics, we synchronize the batch mean and variance across each device. Results are demonstrated in Table 4, Dspike training obtain 4.5% improvement compared to STBP-tdBN [26] and 2.2% improvement compared to Diet-SNN [27].

heading: 5.4 Sparsity, Efficiency, and Energy Cost
text: In this section, we measure the sparsity and efficient computation as well as the energy consumption of SNN. For ANN, the dot-product performs multiply-accumulate (MAC) operation, which is composed of an addition and a multiplication. In the case of SNN, its event-based operation and binary activation enable multiplication-free (except for the first rate-encoding layer) inference. Only one addition is needed if a spike is activated. In Fig. 4 we visualize the mean sparsity of activation of our spiking ResNet-18 on the whole CIFAR10 test set. We can see that each LIF layer has no more than 32% sparsity rate, which greatly saves energy. We also measure the operation number in Table 3. The addition count is calculated by s ∗ T ∗A, where s is the mean sparsity, T is the time step and A is the addition number in ANN. For multiplication in SNN, we set it to the FLOPs of the first layer and scale it by T . Our spiking ResNet-18 with 2 time-step has the lowest operation need for the inference. Then, we measure the energy consumption. Following [27], the energy is measured in 45nm CMOS technology. The MAC operation in ANN costs 4.6pJ energy and the accumulation in SNN costs 0.9pJ energy. Based on this, we compute the energy cost and put it in the table. Our network only costs 80.3mJ for a single forward, and consumes 12.8× lower energy compared with ANN.

heading: 5.5 Ablation Study
text: Our experiments use the same ε for dataset and network architectures. In this section, we use different
value for ε and compare their results. The results is included in Table 5. Generally, we find a low temperature will decrease the performance, as we explained in Sec. 4.2, the desired optimization algorithm should ensure each update will decrease the loss function. Thus, the instantaneous derivative in SNN cannot measure such curvature information since the instantaneous
derivative is always 0. Therefore, we find a large epsilon can help improve the performance. Indeed, the optimal ε requires tuning, but our empirical results show that 0.1 can have good results.

heading: 5.6 Visualization
text: In this section, we provide further experimental visualization. In Fig. 5, we plot the temperature curve for CIFAR10-DVS experiments. The spiking ResNet-18 was trained for 300 epochs with cosine learning rate decay. We initialize b to 1.0 and use ∆b=0.2 to optimize both temperature for gradient calculation and the weights parameters for network learning. For the initial 25 epochs, the temperature remains less-than-1. We think this is primarily due to the complex loss landscape at the starting stage of the training. During the 20 to 80 epoch, the temperature continues to be increasing. In this stage of learning the FDG desires to learn a sharp, narrow gradient with low bias but high variance. We also find the accuracy curve becomes more unstable in this stage of the training, as a consequence of the sharp gradient. After 80 epochs, the temperature becomes decreasing. Although the decreasing process is not aggressive as the former increasing stage, the trend is still obvious. From 80-250, the temperature decreases from 4.0 to 1.6 with some fluctuations. Then, at the end of the training stage, the temperature again becomes increasing.
The intuition behind this tuning mechanism may remain unexplainable right now, however, we think it is the intrinsic nature of the gradient estimation in SNN. Through FDG, we are able to optimize the gradient calculation along with the weights learning.

heading: 6 Discussion
text: In this work, we inspect the training process of SNN when using the surrogate gradient. The proposed finite difference gradient and Dspike function can be assembled with other optimization methods as a fundamental module of gradient estimation for training SNN. This suggests that the performance of SNN can be further increased when more sophisticated training approaches are integrated with the currently proposed methods in the future. Although the improvement on model performance is very significant, be that as it may, our work is still an extension to the conventional SG. How to efficiently find the best FDG gradient to directly deal with the non-linearity of spiking function still remains a problem. In addition, we only consider the LIF neurons for the analyses here, the efficiency of the proposed approaches on other types of spiking neurons still needs examination in practice.

Review summaries:
This paper studies the differentiable training problem of spiking neural network (SNN), which is a promising alternative to the widely used ANN. The authors first identify the problem of the current surrogate gradient of SNN and propose a new surrogate gradient estimator. I really like the author's research method of thinking about problems from experiments, but the result of the experiment is just a tuning experiment of some parameters, and I did not see the author's insight into the existing SNN training problems. Furthermore, the Spike method proposed by the author lacks clear motivation. The proposed Dspike method is just a differentiable piecewise function, which does not have much methodological contribution. The author exaggerated the experimental settings and did not present the experimental results of neuromorphic images. In summary, I think this article does not meet the standards of the NeurIPS conference.  The author’s reply addressed most of my concerns. I agree that AdderNet is another line of work with different purpose. I can see that the authors are very thoughtful in the model design, however it is not well-described in the main paper. I have raised my rating by 1.

The authors outline a principled approach--finite difference gradients (FDG)--for choosing a surrogate gradient for spiking neural networks. FDG is demonstrated to approximate gradient descent in classic rate-based ANNs very well, approximate classic rectangular surrogate gradient function descent in SNNs for medium neighborhoods and thus theoretically justify the method, and be useful in finding an optimal neighborhood in which to approximate a gradient for SNNs.  This last step is then used to optimize a continuously varying surrogate spike gradient estimator (Dspike), modeled as a tanh approximation with the gain as a free parameter. Using the chosen Dspike function (which depends on the time step), the authors then demonstrate the appreciable accuracy and energy efficiency improvements of the resultant SNN.

Summary.  This paper attempts to facilitate the surrogate gradient descent in training SNNs. The proposed FDG quantitatively analyzes the gradients in the training procedure. And then, the authors propose a Dspike function, which is claimed as a proper searching function for the optimal surrogate gradients.  Numerical experiments demonstrate the superiority of the proposed models to several contenders in CIFAR and ImageNet data sets.  Contributions.  I have to admit that it's hard to understand and follow the objective and soundness of this paper. At first glance at the Tables in Section 5, the proposed methods achieve a commendable result compared to other SNN models. However, I can't entirely agree with the author's claimed contributions, especially the first and second ones.  As claimed that "through which (FDG) we are able to not only explain why it (surrogate gradient) works well but also indicate how we can improve its (gradient) adaptivity to data through training", so what are the answers of why and how? In detail,  [1] The derivation of FGD can be traced back to the directional derivative and difference in mathematics. It is reasonable to use Eq. (10) to approximate the instantaneous derivative at a certain point in ANNs, which have continuous and differentiable activations. The high-similarity curves in experiment $\mathrm{I}$ verify my conjectures. It's a good point.  [2] However, the finite difference technique is not suitable for estimating the derivative of holder-continuous and non-differentiable functions. There may be a large gap between the instantaneous and average derivatives. Thus, the similarity between FDG and surrogate gradients in SNNs cannot support the "soundness of surrogate gradients" in Line 175 due to the imprecision of FDG.  [3] Further, how to conclude that "This observation concludes that surrogate gradient can predict the change rate of loss in a relatively larger neighborhood" between Lines 175-176? Doesn't dissimilarity mean invalidity? I hope the authors can clarify this point.  [4] My conjecture in [2] can be verified by Figure 2. At the initial state that epoch = 0 or 1, parameters $\mathbf{W}$ may be far away from the optimal solution. In this case, FDG is not similar to the surrogate gradient if employing small $\epsilon$. As the setting in the experiment that "weight decay (L2 regularization) is set to 0.0001", the parameters may converge to the optimal local minima as weights decay to 0.0001. Thus, FDG still works, although $\epsilon$ is small. So it may be a practical trick but not holds as the authors theoretically claimed, right?  About the claim that "optimal surrogate gradient can be found by choosing temperature wisely" in Lines 210-211. In other words, the Dspike function is effective since optimal surrogate gradient can be found by choosing temperature wisely. How to evaluate this point? I would have appreciated more clarity on the experimental results not as omitted by "As did in former toy experiments (Sec. 4.2), we verify the cosine similarity between our Dspike gradients and the FDG at the beginning of each training epoch and update the  temperature to better match the FDG."  In summary, I find the interests of this work and the proposed model may be effective in practice. However, this conclusion should be testified either theoretically or experimentally. Of course, this paper still needs improvement in these two aspects. So the current version is improper for acceptance.