Exploiting Neuron and Synapse Filter Dynamics in Spatial Temporal Learning of Deep Spiking Neural Network

1 Introduction

Spiking neural networks have demonstrated their capability in signal processing and pattern detection by mimicking the behavior of biological neural systems. In SNNs, information is represented by sparse and discrete spike events. The sparsity of spike activities can be exploited by event driven implementation for energy efficiency. In a more bio-realistic neuron and synapse model, each neuron is a dynamic system, which is capable of spatial temporal information processing. The network made of such neurons can memorize and detect spatial temporal patterns with an ability superior to conventional artificial neural network (ANN) [Wu et al., 2018b].

The potential of SNNs has not been fully explored. First of all, due to the lack of unified and robust training algorithms, the performance of SNNs is still not comparable with deep neural networks (DNN). Directly adapting backpropagation is not feasible because their output is a sequence of Dirac delta functions, hence is non-differentiable. Secondly, most SNN models and training algorithms use rate coding, representing a numerical value in DNN by spike counts in a time window, and consider only the statistics of spike activities. Temporal structure of spike train and spike timing also convey information [Mohemmed et al., 2012]. Spike trains with similar rates may have distinct temporal patterns representing different information. To detect the temporal pattern in the spike train, novel synapse and neuron models with temporal dynamics are needed. However, synapse dynamics are often ignored in the computational models of SNNs.

To address the problem with non-differentiable neuron output, one approach is to train an ANN such as a multi-layer perceptron (MLP) and convert the model to an SNN. This method is straightforward, but it requires additional finetuning of weights and thresholds [Diehl et al., 2015]. There are also works that directly apply backpropagation to SNN training by approximating the gradient of the spiking function [Lee et al., 2016; Esser et al., 2015; Shrestha et al., 2019], or utilizing gradient surrogates [Wu et al., 2018b; Shrestha and Orchard, 2018]. Other approaches include using derivatives of soft spike [Neftci et al., 2019] or membrane potential [Zenke and Ganguli, 2018].

The ability of capturing temporal patterns relies on neuron and synapse dynamics [G¨utig and Sompolinsky, 2006]. Synapse function can be modeled as filters, whose states preserve rich temporal information. The challenge is how to capture the dependencies between the current SNN states and previous input spikes. This challenge has been addressed by some existing works. [G¨utig and Sompolinsky, 2006] and [G¨utig, 2016] train individual neuron to classify different temporal spike patterns. [Mohemmed et al., 2012] is capable to train neurons to associate an input spatial temporal pattern with a specific output spike pattern. However, the aforementioned works cannot be extended to multiple layers and therefore are not scalable. Some recent works utilize backpropagation through time (BPTT) to address the temporal dependency problems. [Wu et al., 2018b] proposed simplified iterative leaky integrate and fire (LIF) neuron model. [Gu et al., 2019] derived an iterative model from a current based LIF neuron. Based on the iterative model, network can be unrolled hence BPTT is possible. However, these works only consider the temporal dynamics of membrane potential, the synapse dynamics and the filter effect of SNN are ignored. There are also works that introduced the concept of IIR and FIR into Multi Layer Perceptron (MLP) [Back and Tsoi, 1991; Campolucci et al., 1999], which enabled MLP to model time series.

In this work, our contributions are summarized as follows:

1. The dynamic behavior of LIF neuron is formulated by infinite impulse response (IIR) filters. We exploit the synapse and neuron filter effect, derive a general representation of SNN as a network of IIR filters with neuron non-linearity.
2. A general algorithm is proposed to train such SNN to learn both rate-based and spatial temporal patterns. The algorithm does not only learn the synaptic weight, but is also capable to optimize the impulse response kernel of synapse filters to improve convergence. The similar learning behavior has been discovered in biological systems [Hennig, 2013]. Our training algorithm can be applied to train simple LIF, and neurons with more complex synapses such as alpha synapse, dual-exponential synapse etc.
3. Our algorithm is tested on various datasets including MNIST, neuromorphic MNIST, DVS128 gesture, TIDIGITS and Australian Sign Language dataset, and outperform state of the art approaches.

2 NeuronModel

Without loss of generality, we consider a LIF neuron with dual exponential synapse for its biological plausibility. The neuron can be described as a hybrid system, i.e. the membrane potential and synapse status evolve continuously over time, depicted by ordinary differential equations (ODE), while a spike event triggers the update of the state variables as the following [Brette et al., 2007]:

Where xi is the state variable of the ith synapse, wi is the associated weight, andM is the total number of synapses. τm and τs are time constants, and η = τm/τs. v and vrest are the neuron membrane potential and rest potential. For simplicity, we set vrest = 0. Every synapse has its own potential, which is called postsynaptic potential (PSP). Neuron accumulates PSP of all input synapses. The membrane potential resets when an output spike is generated.

The ODE system is linear time invariant (LTI). It can also be interpreted as the convolution of an impulse response of a filter with the input spike train, which leads to the spike response model [Gerstner et al., 2014]. The relation between the v(t), O(t) and the historical spike input can clearly be seen in the spike response model. We denote the input spike trains as a sequence of time-shifted Dirac delta functions, Si(t) = ∑j δ(t - tji ), where tji denotes the jth spike arrival time from the ith input synapse. Similarly, output spike train can be defined as O(t) = ∑ δ(t - tf ), tf ∈ {tf : v(tf ) = Vth}. To simplify the discussion, we consider only one synapse. The impulse response kernel k(t) of a neuron described by above ODE system is obtained by passing a single spike at time 0 at the input, such that the initial conditions are x(0) = 1 and v(0) = 0. By solving equation 1a and 1b, we have k(t) = η η η-1 (e -t τm - e -t τs ). Given the general input S(t), PSP is the convolution of k(t) and S(t). For a neuron with M synapses, without reset, the sub-threshold membrane potential is the summation of all PSPs, such that v(t) = ∑Mi wi ∫ ∞ 0 k(s)Si(t - s)ds.

In hybrid model, the reset is modeled by simply setting v to vrest, and regarding the reset as the start of the next evaluation and discarding the neuron's history information. A more biological way is to treat reset as a negative current impulse applied to the neuron itself [Gerstner et al., 2014]. The reset impulse response is h(t) = -Vthe -t τr , where τr controls the decay speed of reset impulse. Such that the membrane potential is the summation of all PSPs and reset voltage:

Treating reset as a negative impulse enables adaptive threshold, which is observed in biological neurons. Neuron's threshold depends on its prior spike activity. With adaptation, frequent spike activity increases the reset voltage, which inhibits the neuron activity, preventing SNNs from over-activation. Such that additional tuning methods such as weight-thresholds balancing [Diehl et al., 2015] is not necessary.

Above equations reveal the filter nature of the biologically realistic neuron model. Each synapse act like a low pass filter. Synapse filter is causal, and the kernel is defined to decay over time, hence the current state of the PSP is determined by all previous input spikes up to current time. The temporal dependency calls for temporal error propagation in the training.

3 Neuron and Synapse as IIR Filters

In practice, for computational efficiency, spiking neural network are usually simulated in discrete time domain and network states are evaluated for every unit time. The discrete time version of equation 2 can be written as:

where t ∈ Z≥0. It is clear that v[t] is a combination of a reset filter and multiple synapse filters. However, the above form is not practical for implementation because of infinite convolution coefficients. We express the above system using Linear Constant-Coefficient Difference Equations (LCCD):

Introducing synapse dynamics could cause significantly large computation overhead because the number of synapses is quadratic to the number of neurons. Maintaining such large number of synaptic states is infeasible. In a biological system, spikes are transmitted through axons, an axon connects to multiple destination neurons through synapses. Therefore, the synapses that connect to the same axon have identical spike history hence same states. Based on this observation, tracking the states of synapses that have the same fan-in neuron is unnecessary as these synapses can share the same state and computation.

Neuron itself can also be a filter and v[t] may also rely on its previous states. We can extend equation 4a - 4b to a more general form, such that the SNN can be interpreted as a network of IIR filters with non-linear neurons:

Where l and i denote the index of layer and neuron respectively, and j denotes input index and t is the time, Nl is number of neurons in lth layer. V l i [t] is neuron membrane potential. I l i [t] is weighted input. Rli [t] is reset voltage, F j l [t] is PSP. OL i [t] is spike function, and U(x) is a Heaviside step function. P and Q denote the feedback and feed forward orders. λ, θ, αlj,p and βlj,q are coefficients of neuron filter, reset filter and synapse filter respectively. 5d is a general form of IIR filters, it allows PSP to be arbitrary shapes. The above formulation is not specific to neuron models and it provides a flexible and universal representation, it is capable of describing more complex spiking neuron models than LIF neuron. For example, by setting α1 = 2e -1 τ , α2 = -e -2 τ , αp = 0, p ∈ {2, 3, ..., P}, β1 = 1 τ e- 1 τ and βq = 0, q ∈ {0, 2, 3, ...,Q}, it models neuron with alpha synapse. By setting αp = 0, p ∈ {1, 2, ..., P}, β0 = 1, βq = 0, q ∈ {1, 2, ...Q}, the synapse filter is removed, the model becomes simple LIF neuron as in [Diehl et al., 2015; Gu et al., 2019]. Based on 5a - 5f, a general model of spiking neuron can be represented as a network of IIR as shown in Figure 1. Axonal delay is explicitly modeled in equation 5d by delayed input βlj,q Ol-1 j [t - q], hence it enables more complex and biologically plausible temporal behavior. Neurons can also have heterogeneous synapses, i.e. the synapses' feed forward order and feedback order can vary across layers. To avoid notation clutter, we assume that all neurons in this paper have homogeneous synapse types.

Equation 5a to 5f provide an explicitly iterative way to model synapse and neuron dynamics, hence it is possible to unfold the network over time and apply BPTT. The spatial and temporal data flow and unfolded network with second order synapse filter are shown in Figure 2. Similar formulations can be found in [Wu et al., 2018b; Gu et al., 2019]. However they are aimed at specific neuron models.

4 Spatial Temporal Error Propagation

We discuss the spatial temporal backpropagation in the context of two learning tasks. In the first, the neuron that fires most represent the correct result. Since this is a classification task, we use cross-entropy loss and spike count of the output neuron represents the probability. Loss is defined as:

where yi is the label, L is number of layers, OL i [t] denotes output of last layer.

In the second learning task, the goal is to train SNN to generate spikes at specified times such that the output spike pattern O[t] is spatially and temporally similar as the target spike pattern Starget[t]. We refer to it as temporal learning. The loss function of the learning is the distance between the actual output spike trains and the target spike trains. Inspired by Van Rossum distance, we pass the actual and target spike train through a synapse filter k[t], to convert them to continuous traces. The loss is defined as:

where Si

target [t] is the ith spike train of target spike patterns. For both tasks, we define: δl i [t] = ∂E ∂Oli [t] , ϵli [t] =

∂U(V l i [t]-Vth) ∂V l i [t] , κli [t] = ∂V l i [t+1] ∂(V l i [t]) . Please note that the spike activation function U(x) is not differentiable. Its approximation will be discussed in section 4.1. By unfolding the model into spatial path and temporal path as shown in Figure 2, BPTT can be applied to train the network. κli [t] can be computed as:

δl i [t] an be computed recursively as follows:

where

Where δl i [t+q] = 0 for t+q > T . Unlike LSTM/RNN, or SNN such as [Wu et al., 2018b; Gu et al., 2019], there may be dependency from layer l + 1 to layer l at multiple time steps due to axonal delay. Based on above equations, error can propagate recursively. By applying chain rule, we can obtain the gradient with respect to weight:

In real biological system, synapses may respond to spike differently. The PSP kernels can be modulated by input as part of the synaptic plasticity [Hennig, 2013]. It is possible to employ gradient descent to optimize the filter kernels in equation 5d [Campolucci et al., 1999]. The gradients of L with respect to αlj,p and βlj,p are:

where ∂I l i [t] ∂αlj,p and ∂I l i [t] ∂βlj,q are:

Above learning rule assumes the SNN to be an LTI system. The loss calculation, error propagation, filter coefficients and synaptic weights update are performed at the end of each training iteration. Therefore, within one iteration, the SNN is still linear time-invariant.

4.1 Spike Function Gradient Approximation

The non-differentiable spike activation is a major road-block for applying backpropagation. One solution is to use a gradient surrogate [Neftci et al., 2019]. In the forward path, a spike is still generated by a hard threshold function, while in the backward path, the gradient of the hard threshold function is replaced by a smooth function. One of such surrogates can be spike probability [Esser et al., 2015; Neftci et al., 2019]. Although the LIF neuron is deterministic, stochasticity can be obtained from noise [Stevens and Zador, 1996]. Under Gaussian noise of mean 0 and variance σ, in a short interval, LIF neuron can behave like a Poisson neuron such that the spike probability is a function of the membrane potential v as follows:

where erfc(x) represents a complementary error function. With this replacement, the gradient of U(x) can be approximated as:

5 Experiments

Proposed model and algorithm are implemented in PyTorch 1 . We demonstrate the effectiveness using three experiments; the first experiment is a non-trivial generative task using associative memory; the second is vision classification, and the third is to classify temporal patterns. In following experiments, we use Adam optimizer, learning rate is set to 0.0001, batch size is 64. We employ synapse model depicted by equation 4c, in which τm = 4, τs = 1, λ = 0, θ = e -1 τm , Vth = 1.

5.1 Associative Memory

An associative memory network retrieves stored patterns that most closely resembles the one presented to it. To demonstrate the capability of our approach to learn complex spatial temporal spike patterns, we train a network of structure 300x500x200x500x300. We generate 10 spatial temporal spike train patterns, each contains 300 spike trains of length 300, samples of these patterns are shown in Figure 3a. Each dot corresponds to a spike event, the x-axis represents the time, and the y-axis represents the spike train index. The SNN is trained to reconstruct the pattern. First column of 3b shows two noisy sample inputs. Noisy samples are formed by adding random noise, which includes obfuscation and deletion of some part of the patterns, jitter in input spikes' timing following a Gaussian distribution and random background spikes. After 50 epochs of training, the network is able to reconstruct the original patterns and remove background noise. Corresponding outputs at epoch 5 and 50 are shown in 3b. Such a task is difficult for rate-based training methods as they are not capable of capturing temporal dependencies. It is noteworthy that the intermediate layer has 200 neurons, which is smaller than the input layer. And the intermediate layer is learning the spatial and temporal representation of the input patterns. Thus, this network also acts like a spatial temporal auto-encoder.

We drove the input of the network with 64 different testing samples and record output of 200 neurons in the intermediate layer. Figure 4a color codes the spiking rate of those neurons. The x-axis gives the index of the neurons, and y-axis gives the index of different testing samples. Those samples belong to 10 different classes, and are sorted so that data of the same class are placed close to each other vertically. The 10 different colors on the left side bar indicate each of the 10 classes. The pixel (x, y) represents the spiking rate of neuron x given testing sample y. Spiking rate of any neuron is almost a constant regardless of which class the testing sample belongs to. Figure4b shows the Van Rossum distances between the 200 neurons' output spike train. x-axis and y-axis give the input sample index. The color intensity of pixel (x, y) is proportional to the Van Rossum distance between the 200 neurons' output when given input sample x and y respectively. Similar as 4a, the color bar on left side indicates the class of each sample. It can clearly be seen from figure 4b that the temporal structure of these 200 neurons' outputs are significantly different. The fact that our model is able to take those 64 sets of spike trains with almost the same firing rate and generate

1 Code is available at: https://github.com/Snow-Crash/snn-iir

10 different classes indicates that it is capable of utilizing features in the temporal distribution of the spikes in addition to the spike rates.

5.2 Vision Tasks

We evaluated our method on three vision datasets. Results and comparisons with state-of-the-art works in the SNN domain are shown in Table 1. For MNIST, we utilize ratebased encoding to convert the input image into 784 spike trains where number of spikes in the spike train is proportional to the pixel value. With a convolutional SNN with the structure 32C3-32C3-64C3-P2-64C3-P2-512-10, our model achieves state-of-the-art accuracy in the SNN domain. The work next in terms of accuracy (99.42 %) [Jin et al., 2018] employs ensemble learning of 64 spiking CNNs. Compared to conversion-based approaches that require hyper-parameter search and fine tuning [Diehl et al., 2015], our approach does not require post-training processing. It directly trains SNN using BPTT and obtains models with comparable quality as DNN.

Unlike MNIST, which consists of static images, Neuromorphic MNIST (N-MNIST) is a dynamic dataset which consists of spike events captured by DVS camera and is a popular dataset for SNN evaluation. An N-MNIST sample is obtained by mounting the DVS camera on a moving platform to record MNIST image on the screen. The pixel change triggers spike event. Thus, this dataset contains more temporal information. With a convolutional network of size 32C3-32C3-64C3-P264C3-P2-256-10, our model outperforms the current state-ofthe-art. The results are shown in Table 1. [Lee et al., 2016] introduced additional winner-take-all (WTA) circuit to improve performance. [Wu et al., 2019] gets 99.35% accuracy with a very large network, the structure is 128C3-256C3-AP2512C3-AP2-1024C3-512C3-1024FC-512FC-Voting. There is also additional voting circuit at output layer. We use a significantly smaller network to achieve the same accuracy, and no additional voting layer or WTA circuits are required.

DVS128 Gesture Dataset contains 10 hand gestures such as hand clapping, arm rolling etc. collected from 29 individuals under 3 illumination conditions using DVS camera. The network is trained to classify these hand gestures. This dataset contains rich temporal information. For this task, we utilize a network with 64C7-32C3-32C3-256-10 structure. The advantage of our work is clearly seen in the third column of Table 1. We achieved 96.09 % accuracy, which is state-of-theart in the spiking domain, while other works, such as [Amir et al., 2017], requires additional filters for data preprocessing and WTA circuit at the output layer. Our model and learning algorithm doesn't need specialized neuron circuits or any data preprocessing techniques as the spike streams are directly fed into the network.

We also studied the effect of training the synapse response kernels. The learned synapse kernels are shown in Figure 5. The solid red line represents the original kernel. The decay speed of synapse response of the learned kernel diverges from original kernel. Slower decay speed indicates the synapses are capable of remembering information for a longer time. Such behavior is similar to the gates in an LSTM. The accuracy with and without training synapse filter kernel are shown in Figure 6. No improvements are observed for MNIST dataset, the accuracy with training and without trained kernel are 99.46% and 99.43% respectively. This is because MNIST is a static dataset, hence no temporal information. There is slight improvement in NMNIST by training synapse filter kernel, the accuracy increases from 99.24% to 99.39%. In DVS 128 dataset, the advantage of training the synapse filter kernel is clearly seen, the model not only converges faster, the accuracy also increases from 94.14% to 96.09%.

5.3 Time Series Classification

Our work also shows advantages in time series classification. We evaluated our work in TIDIGITS and Australian Sign Language [Kadous and others, 2002] dataset. TIDIGITS is a speech dataset that consists of more than 25,000 digit sequences spoken by 326 individuals. For training and testing, we extracted MFCC from each sample, resulting 20 time series of length 90. The Australian sign language dataset [Kadous and others, 2002] is a multivariate time series dataset, collected from 22 data glove sensors that track acceleration and hand movements such as roll, pitch etc. Each recorded hand sign is a sequence of sensor readings. The average duration of a hand sign is 45 samples. The dataset has 95 classes of hand signs, To convert time series into spike trains, we use current-based LIF neuron as encoder. It accumulates input data as current and converts time varying continuous values to time varying spike patterns.

Networks to classify TIDIGITS and Australian Sign Language have a structure 300-300-11 and 300-300-95 respectively. We trained two-layer stacked LSTM of unit size 300 as baseline. Results are shown in Table 2. The best accuracy in TIDIGITS is achieved by [Abdollahi and Liu, 2011], however, it is a non-spiking approach. In Australian Sign Language dataset, we outperformed vanilla LSTM and DNN based approaches. [Shrestha et al., 2019] uses EMSTDP to train an SNN to classify 50 classes of the hand signs, the network size is 990-150-150-50. It buffers the entire sequence and flattened the time series into a vector. While our work is trained to classify all 95 classes, and it processes the time series in a more efficient and natural way, the input data is converted into spikes on the fly. Since flattening is no longer necessary, the input dimension is also reduced.

6 Conclusion

In this work, we proposed a general model to formulate SNN as network of IIR filters with neuron non-linearity. The model is independent of neuron types and capable to model complex neuron and synapse dynamics. Based on this model, we derived a learning rule to efficiently train synapse weights and synapse filter impulse response kernel. The proposed model and method are evaluated on various tasks, including associative memory, MNIST, NMNIST, DVS 128 gesture, TIDIGITS etc. and achieved state-of-the-art accuracy.