The Remarkable Robustness of Surrogate Gradient Learning for Instilling Complex Function in Spiking Neural Networks

The computational power of deep neural networks (LeCun, Bengio, & Hinton, 2015; Schmidhuber, 2015) has reinvigorated interest in using in-silico systems to study information processing in the brain (Barrett, Morcos, & Macke, 2019; Richards et al., 2019). For instance, performance-optimized artificial neural networks bear striking representational similarity with the visual system (Maheswaranathan et al., 2018; McClure & Kriegeskorte, 2016; McIntosh, Maheswaranathan, Nayebi, Ganguli, & Baccus, 2016; Pospisil, Pasupathy, & Bair, 2018; Tanaka et al., 2019; Yamins & DiCarlo, 2016; Yamins et al., 2014) and serve to formulate hypotheses about their mechanistic underpinnings. Similarly, the activity of artificial recurrent neural networks optimized to solve cognitive tasks resembles cortical activity in prefrontal (Cueva et al., 2019; Mante, Sussillo, Shenoy, & Newsome, 2013), medial frontal (Wang, Narain, Hosseini, & Jazayeri, 2018), and motor areas (Michaels, Schaffelhofer, Agudelo-Toro, & Scherberger, 2019; Stroud, Porter, Hennequin, & Vogels, 2018), thus providing us with new vistas for understanding the dynamic properties of computation in recurrent neural networks (Barrett et al., 2019; Sussillo & Barak, 2012; Williamson, Doiron, Smith, & Yu, 2019).

All of these studies rely on conventional artificial neural networks with graded activation functions as commonly used in machine learning. The recipe for building a deep neural network is straightforward. The value of a scalar loss function defined at the output of the network is decreased through gradient descent. Deep neural networks differ from biological neural networks in important respects. For instance, they lack cell type diversity and do not obey Dale's law while ignoring the fact that the brain uses spiking neurons. We generally accept these flaws because we do not know how to construct more complicated networks. For instance, gradient descent works only when the involved system is differentiable. This is not the case for spiking neural networks (SNNs).

Surrogate gradients have emerged as a solution to build functional SNNs capable of solving complex information processing problems (Bellec, Salaj, Subramoney, Legenstein, & Maass, 2018; Cramer, Stradmann, et al., 2020; Esser et al., 2016; Hunsberger & Eliasmith, 2015; Lee, Delbruck, & Pfeiffer, 2016; Neftci, Mostafa, & Zenke, 2019; Pfeiffer & Pfeil, 2018; Shrestha & Orchard, 2018). To that end, the actual derivative of a spike, which appears in the analytic expressions of the gradients, is replaced by any well-behaved function. There are many possible choices of such surrogate derivatives, and consequently, the resulting surrogate gradient is, unlike the true gradient of a system, not unique. A number of studies have successfully applied different instances of surrogate derivatives to various problem sets (Bellec et al., 2018; Esser et al., 2016; Huh & Sejnowski, 2018; Shrestha & Orchard, 2018; Woźniak, Pantazi, Bohnstingl, & Eleftheriou, 2020; Zenke & Ganguli, 2018). While this suggests that the method does not crucially depend on the specific choice of surrogate derivative, we know relatively little about how the choice of surrogate gradient affects the effectiveness and whether some choices are better than others. Previous studies did not address this question because they solved different computational problems, thus precluding a direct comparison. In this letter, we address this issue by providing benchmarks for comparing the trainability of SNNs on a range of supervised learning tasks and systematically vary the shape and scale of the surrogate derivative used for training networks on the same task.

To systematically evaluate the performance of surrogate gradients, we sought to repeatedly train the same network on the same problem while changing the surrogate gradient. Toward this end, we required a demanding spike-based classification problem with a small computational footprint to serve as benchmark. There are only few established benchmarks for SNNs. One approach is to use analog-valued machine learning data sets as input currents directly (Hunsberger & Eliasmith, 2015) or to first convert them to Poisson input spike trains (Lee et al., 2016; Pfeiffer & Pfeil, 2018). These input paradigms, however, do not fully capitalize on the ability to encode information in spike timing, an important aspect of spiking processing. Gütig (2016) addressed this point with the Tempotron by classifying randomly generated spike timing patterns in which each input neuron fires a single spike. Yet completely random timing precludes the possibility of assessing generalization performance, that is, the ability to generalize to previously unseen inputs.

To assess if SNNs could learn to categorize spike patterns and generalize to unseed patterns, we created a number of synthetic classification data sets with added temporal structure. Specifically, we created spike rasters for a given set of input afferents. Each afferent only fired one spike, and the spike times of all afferents were constrained to lie on a low-dimensional smooth, random manifold in the space of all possible spike timings. All data points from the same manifold were defined as part of the same input class, whereas different manifolds correspond to different classes.

We first confirmed that learning is poor when we used the actual gradient. To that end, we computed it using the derivative of the hard threshold nonlinearity of the spikes. As expected, the hard threshold nonlinearity prevented gradient flow into the hidden layer (Neftci et al., 2019) and consequently led to poor performance (see Figures 2b and 2c). In contrast, when we used surrogate gradients to train the same network, the problem disappeared. Learning took place in both the hidden and output layers and resulted in a substantial reduction of the loss function (see Figures 2b to 2e).

Next we tested different surrogate derivative shapes, namely a standard sigmoid (Sigmoid$'$⁠) and piece-wise linear function (Esser et al., 2016; see Figure 3a; Bellec et al., 2018). This manipulation led to a reduction of the size of the parameter regime in $β$ in which the network was able to perform the task, which is presumably due to vanishing gradients (Hochreiter, 1998). However, there was no substantial reduction in maximum performance (see Figures 3b and 3c). Using a piece-wise linear surrogate derivative (Esser et al.) led to a further reduction of viable parameters $β$ (see Figure 3b), but did not affect maximum performance, regardless of whether we used one or two hidden layers (see Figure 3c). To check whether a surrogate derivative was required at all for solving the random manifold problem, we assayed the learning performance for $β=0$⁠, which corresponds to setting the function to 1. This change resulted in a significant drop in performance comparable to a network without hidden units (see Figure 3c) suggesting that a nonlinear voltage dependence is crucial to learn useful hidden-layer representations. Finally, we confirmed that these findings were robust to different initial network parameters and data sets (see Figure 3d) apart from only a few outliers with low performance in the two-hidden-layer case (⁠$nh=2$⁠). These outliers point at the vital role of proper initialization (He, Zhang, Ren, & Sun, 2015; Mishkin & Matas, 2016).

To gain a better understanding of how derivative scales larger than one affect surrogate gradient learning, we trained networks on a fixed random manifold task (see Figures 3a to 3c), using an asymptotic version of the SuperSpike surrogate (aCtl; see Figure 3) and our well-tested standard SuperSpike function (sCtl; see Figure 3) as a control (sCtl). As we expected, the difference in scale to manifest itself primarily in the presence of recurrence, we compared networks in which we treated the spike reset as differentiable (DR) with networks in which its contribution was ignored by detaching it from the computational graph. Technically speaking, we prevented PyTorch's auto-differentiation routines (Paszke et al., 2019) from considering connections in the computational graph that correspond to the spike reset when computing the gradient with backpropagation through time (BPTT) (see equation 4.1). While both the normalized (sCtl) and the asymptotic surrogate (aCtl) performed equally well when the reset was detached, combining a differentiable reset with an asymptotic surrogate derivative (aDR) led to impaired performance (see Figures 5b and 5c). This adverse effect on learning was amplified in deeper networks (see Figure 5d). Thus, the scale of the surrogate derivative plays an important role in learning success if implicit recurrence, as contributed by the spike reset, is present in the network.

Since the spike reset constitutes a specific form of implicit recurrence (cf. see Figure 4), we were wondering whether we would observe a similar phenomenon for explicit recurrence through recurrent synaptic connections. To that end, we repeated the performance measurements in networks with recurrent connections but kept the spike reset term detached to prevent gradient flow. We observed a small but measurable reduction in accuracy for the best-performing networks (see Figures 5c and 5e). Importantly, however, there was a substantial decrease in classification performance when gradients were allowed to flow through recurrent connections and the asymptotic SuperSpike variant was used (see Figure 5e). Unlike in the differentiable reset (DR) case, the effect was severe enough to drop network performance to chance level even for a network with a single hidden layer (see Figure 5f). In summary, surrogate gradients are sensitive to the scale of the surrogate derivative. More specifically, when the scale of the surrogate derivative is too large, and either implicit or explicit recurrence is present in the network, the effect on learning can be detrimental.

Further, to study the effect of explicit recurrence, we ran separate experiments for recurrently connected networks (RC). Importantly, we did not observe any substantial performance differences between strictly feedforward and recurrently connected networks either (see Figures 6a and 6b).

We speculated that the absence of an effect may be due to the short duration of the input spike trains considered so far (approximately 50 ms). Recurrent connections are typically thought to provide neural networks with longer timescale dynamics, effectively giving the network a working memory. Thus, the beneficial effects of recurrent connections may only emerge when using stimuli of longer duration, with potentially multiple spikes from any given input neuron.

To test this hypothesis, we trained networks on the Spiking Heidelberg Digits (SHD) data set (Cramer, Stradmann et al., 2020), which consists of simulated input spikes from the auditory pathway of varied duration between 0.6 and 1.4 s (see Figure 7c). Indeed, we found that the best-performing models in this case were recurrent and achieving state-of-the-art classification accuracy of $0.82±0.02%$ (see Figure 6c). These data are consistent with the notion that working memory plays a vital role in classifying longer input patterns.

Until now, we have considered spike-based data sets. While spiking inputs are arguably the most natural input to SNNs, they come with an important caveat. All spiking data sets assume a specific encoding model, which is used to convert analog input data into a spiking representation. The chosen model, however, may not be optimal and thus adversely affect classification performance. To avoid this issue, we sought to learn the spike encoding by directly feeding current-based input to a set of spiking units (Zimmer, Pellegrini, Singh, & Masquelier, 2019). To test this idea, we converted the raw audio data of the Heidelberg Digits (Cramer, Stradmann et al., 2020) to Mel-spaced spectrograms (see Figure 7d and section 4). To reduce overfitting, we decreased the number of channels and time steps to values commonly used in artificial speech recognition systems. Specifically, we used 40 channels and 80 time frames, corresponding to compression in time by a factor of about five (see section 4). The networks trained on this RawHD data set showed reduced overfitting yet still benefited from recurrent connections compared to strictly feedforward networks (see Figure 6d). Concretely, recurrent networks reached $94±2%$ test accuracy, whereas the feedforward networks reached only $85±3%$⁠. In agreement with the results on Randman and MNIST, there were no significant differences between the Sum and Max readout configurations (see Figure 6).

We wanted to know whether the discrepancy between feedforward and recurrent networks would increase for more challenging data sets. This was made possible by the performance gain from reducing the input dimension to 40 channels, which allowed us to train SNNs on the larger Speech Command data set (Warden, 2018; see section 4). This data set contains over 100,000 utterances from 35 classes, including “yes,” “no,” and “left.” In contrast to the original intended use for keyword spotting, here we assayed its top-1 classification accuracy over all classes, a more challenging problem than accurate detection of only a subset of words. The best SNNs achieved $85.3±0.3%$ on this challenging benchmark (see Figure 6e). On the same task, the spiking feedforward network performed at $70±2%$⁠. There was a clear benefit from adding recurrent connections to the network, with the performance of the Max (⁠$85.3±0.3%$⁠) being slightly better than the Sum (⁠$80.7±0.4%$⁠) readout configuration.

These findings illustrate that surrogate gradient learning is robust to changes in the input paradigm, including spiking and nonspiking data sets. For more complex data sets, recurrently connected networks performed better than strictly feedforward networks. Finally, in the majority of cases, surrogate gradient learning was robust to the details of how the output loss was defined.

We investigated whether surrogate gradients could instantiate SNNs in this biologically plausible, sparse activity regime. To that end, we trained SNNs with added activity regularization that penalized high spiking activity (see Figure 8b and section 4) and recorded their hidden-layer activity. While neurons in the unregularized network frequently emityed the maximally possible number of 500 spikes per trial (see Figure 8c), regularization drastically reduced both the maximal spike count per event and the overall probability of high spike count events (see Figure 8d).

In the case of MNIST, fewer than 10 spikes were sufficient on average to achieve the point of diminishing returns beyond which additional spiking activity did not improve classification performance.

This trend could be replicated for all other data sets, with varying degrees of spike reduction. Adding a hidden layer generally required more spikes for the same performance. Recurrency generally did not have a great effect on the minimum number of spikes and did not improve performance except on RawSC. On RawSC, 80% classification accuracy was achieved only by a few feedforward networks with more than 2000 spikes on average. In the recurrent network, this level was already attained with about 150 spikes (see Figure 9).

In all cases, the transition from chance level to maximum accuracy occurred in less than one order of magnitude change in the mean number of spikes. On all data sets we tested, the addition of a second hidden layer led to an overall increase in mean spiking activity, which did not yield a notable performance change on Randman and MNIST but resulted in small improvements on RawHD and RawSC.

These results illustrate that activity-regularized SNNs can perform with high accuracy down to some critical activity threshold at which their performance degrades rapidly. Importantly, we found several network configurations that showed competitive performance with an average number of spikes substantially lower than the number of hidden units. For instance, to classify MNIST with high accuracy, an average of 10 to 20 action potentials was sufficient. Such low activity levels are more consistent with the sparse neuronal activity observed in biological neural circuits and illustrate that surrogate gradients are well suited to building SNNs that use such plausible sparse activity levels for information processing.

Surrogate gradients offer a promising way to instill complex functions in artificial models of spiking networks. This step is imperative for developing brain-inspired neuromorphic hardware and using SNNs as in silico models to study information processing in the brain. In this letter, we have focused on two aspects of surrogate gradient learning in SNNs. We showed, using a range of supervised classification problems, that surrogate gradient learning in SNNs is robust to different shapes of surrogate derivatives. In contrast, inappropriate choice of scale adversely affected learning performance. Our results imply that for practical applications, surrogate derivatives should be appropriately normalized. Second, by constraining their activity through regularization, we showed that surrogate gradients could produce SNNs capable of efficient information processing with sparse spiking activity.

Surrogate gradients have been used by a number of studies to train SNNs (Neftci et al., 2019), solve small-scale toy problems with fractionally predictive neurons (Bohte, 2011), train convolutional SNNs on challenging neuromorphic (Amir et al., 2017; Orchard, Jayawant, Cohen, & Thakor, 2015) and vision benchmarks (Esser et al., 2016), or train recurrent SNNs on temporal problems requiring working memory (Bellec et al., 2018; Shrestha & Orchard, 2018). These studies used different surrogate derivatives ranging from exponential (Shrestha & Orchard, 2018), piece-wise linear (Bellec et al., 2018), or tanh (Woźniak et al., 2020), sometimes with a nonstandard neuron model with a constant leak term (Esser et al., 2016), but due to the different function choices and data sets, they are not easily comparable. Here we provide such a comprehensive comparison.

Towards this end, we had to make some compromises. Our study, like previous ones, is limited to supervised classification problems, because supervised learning offers a well-defined and intuitive quantification of computational performance. To keep the number of model parameters tractable, we focused on current-based LIF neurons. Moreover, we entirely dispensed with Dale's law and relied solely on all-to-all connectivity. However, we believe that most of our findings will carry over to more realistic neuronal, synaptic, and connectivity models. Our study thus provides a set of blueprints and benchmarks to accelerate the design of future studies.

An alternative way of training SNNs, which does not rely on surrogate gradients, was introduced by Huh and Sejnowski (2018), who put forward a differentiable formulation of the neuronal spiking dynamics, thus allowing gradient descent with exact gradients. However, thus far, the approach has only been demonstrated for theta neurons, and nonleaky integrate-and-fire neurons and an extension to neuron models with nondifferentiable reset dynamics, like the LIF neurons we used here, is still pending.

Although considered biologically implausible (Crick, 1989), we have limited our study to training SNNs with backpropagation, the de facto standard for computing gradients in systems involving recurrence and hidden neurons. Although there exist more plausible forward-in-time algorithms like real-time recurrent learning (RTRL) (Williams & Zipser, 1989) they are prohibitively expensive or require additional approximations that affect learning performance (Bellec et al., 2019; Murray, 2019; Neftci et al., 2019; Zenke & Ganguli, 2018). Instead, using BPTT enabled the targeted manipulation of gradient flow through different elements of the network, which allowed us to dis-entwine some of the complexity of surrogate gradient learning. Finally, our study was purely numerical, and several important questions remain open. For instance, how does one optimally initialize hidden-layer weights? How well do these findings translate to convolutional neural networks? And why do the surrogate gradients work so well, despite ignoring the spike reset in their evaluation? Answering these questions requires additional numerical experiments and a rigorous theoretical understanding of surrogate gradient learning in SNNs, both of which remain for future work.

In summary, surrogate gradients allow translating the success of deep learning to biologically inspired SNNs by optimizing their connectivity toward functional complexity through end-to-end optimization. The in-depth study of both surrogate gradients and the resulting functional SNNs will occupy scientists for years to come and will likely prove transformative for neural circuit modeling.

In this study, we used a number of synthetic and real-world learning tasks with the overarching aim of balancing computational feasibility and practical relevance. To that end, we focused on synthetic data sets generated from random manifolds and real-world auditory datasets.

Concretely, we generated spike trains by randomly and uniformly sampling points from a $D$-dimensional hypercube with a side length of one. We then interpreted these points as intrinsic manifold coordinates and converted them to the corresponding embedding space coordinates. We next standardized these values along all embedding dimensions to lie between 0 and $τrandman$ and interpreted these $M$-dimensional coordinates as the firing times of $M$ distinct input neurons belonging to the same class of the classification problem (see Figure 1c). Example code to generate smooth random manifolds following this procedure is available at https://github.com/fzenke/randman.

We chose a default parameter set for most of our experiments that struck a good balance between minimizing both the embedding dimension and the number of samples needed to solve the problem while simultaneously not being solvable by a two-layer network without hidden units. Specifically, we chose a 10-way problem with $D=α=1$⁠, $M=20$⁠, and $τrandman=50ms$ and fixed the random seed in situations in which we reused the same data set. In practice, the spike time $tic$ of input neuron $i$ in a given trial of class $c$ is $tic=fic(X)$⁠, where the $fic$ are fixed random functions and $X$ is a uniformly distributed random number between zero and one. For all simulation experiments, we generated 1000 data points for each class, out of which we used 800 for training and two sets of 100 each for validation and testing purposes.

We used both spiking and nonspiking auditory data sets of digit and word utterances. Specifically, we used the SHDs without any further preprocessing (Cramer, Stradmann et al., 2020). For performance reasons and to dispense with the spike conversion process, we ran additional simulations with nonspiking auditory inputs. Specifically, we worked with the raw Heidelberg Digits (RawHD) and Pete Warden's Speech Commands data set (RawSC) (Warden, 2018) which were preprocessed as follows. We first applied a preemphasis filter to the raw audio signal $x(t)$ by computing $y(t)=x(t)-0.95x(t-1)$⁠. We then computed 25 ms frames with a 10 ms stride from the resulting signal and applied a Hamming window to each frame. For each frame, we computed the 512-point fast Fourier transform to obtain its power spectrum. From the power spectrum, we further computed the filter banks by applying 40 triangular filters on a Mel-scale (Huang, Acero, Hon, & Reddy, 2001). After cropping or padding to 80 (RawHD) or 100 (RawSC) steps by repeating the last frame, the analog-valued filter banks were fed directly to the SNNs.

To train SNN models with surrogate gradients, we implemented them in PyTorch (Paszke et al., 2019). To that end, all models were explicitly formulated in discrete time with time step $Δt$⁠.

We used all-to-all connectivity in all simulations without bias terms unless mentioned explicitly. The weights were initialized from a uniform distribution $U(-k,k)$ with $k=1ninputs$ where $ninputs$ is the number of afferent connections.

We minimized the loss $Ltot$ by adjusting the parameters $W$ and $V$ in the direction of the negative surrogate gradients using Adam with default parameters (Kingma & Ba, 2014). Surrogate gradients were computed with backpropagation through time using PyTorch's automatic differentiation capabilities. To deal with the nondifferential spiking nonlinearity of the hidden-layer neurons, we approximated their derivatives $S'Uiμ[n]=Θ'Uiμ[n]-1$ with suitable surrogates $σ'Uiμ[n]=hUiμ[n]-1$⁠. Throughout, we used the following functions $h(x)$⁠:

• SuperSpike:$h(x)=1βx+12$

• Sigmoid$'$⁠:$h(x)=s(x)1-s(x)$ with the sigmoid function $s(x)=11+exp(-βx)$

• Esser et al.:$h(x)=max(0,1.0-βx)$

F.Z. was supported by the Wellcome Trust (110124/Z/15/Z) and the Novartis Research Foundation. T.P.V. was supported by a Wellcome Trust Sir Henry Dale Research fellowship (WT100000), a Wellcome Trust Senior Research Fellowship (214316/Z/18/Z), and an ERC Consolidator Grant SYNAPSEEK.