Title: TEMPORALLY EFFICIENT DEEP LEARNING WITH SPIKES

Abstract: The vast majority of natural sensory data is temporally redundant. For instance, video frames or audio samples which are sampled at nearby points in time tend to have similar values. Typically, deep learning algorithms take no advantage of this redundancy to reduce computations. This can be an obscene waste of energy. We present a variant on backpropagation for neural networks in which computation scales with the rate of change of the data -not the rate at which we process the data. We do this by implementing a form of Predictive Coding wherein neurons communicate a combination of their state, and their temporal change in state, and quantize this signal using Sigma-Delta modulation. Intriguingly, this simple communication rule give rise to units that resemble biologically-inspired leaky integrate-and-fire neurons, and to a spike-timing-dependent weight-update similar to Spike-Timing Dependent Plasticity (STDP), a synaptic learning rule observed in the brain. We demonstrate that on MNIST, on a temporal variant of MNIST, and on Youtube-BB, a dataset with videos in the wild, our algorithm performs about as well as a standard deep network trained with backpropagation, despite only communicating discrete values between layers.

1. INTRODUCTION:

Currently, most algorithms used in Machine Learning work under the assumption that data points are independent and identically distributed, as this assumption provides good statistical guarantees for convergence. This is very different from the way data enters our brains. Our eyes receive a single, never-ending stream of temporally correlated data. We get to use this data once, and then it's gone. Moreover, most sensors produce sequential, temporally redundant streams of data. This can be both a blessing and a curse. From a statistical learning point of view this redundancy may lead to biased estimators when used to train models which assume independent and identically distributed input data. However, the temporal redundancy also implies that intuitively not all computations are necessary.

Online Learning is the study of how to learn in this domain -where data becomes available in sequential order and is given to the model only once. Given the enormous amount of sequential data, mainly videos, that are being produced nowadays, it seems desirable to develop learning systems that simply consume data on-the-fly as it is being generated, rather than collect it into datasets for offline-training. There is, however a problem of efficiency, which we hope to illustrate with two examples:

1. CCTV feeds. CCTV Cameras collect an enormous amount of data from mostly-static scenes.

The amount of new information in a frame, given the previous frame, tends to be low, i.e. the data tends to be temporally redundant. If we want to train a model from of this data (for example a pedestrian detector), we need to process a large amount of mostly-static frames.

If the frame rate doubles, so does the amount of computation. Intuitively, it feels that this should not be necessary. It would be nice to still be able to use all this data, but have the amount of computation scale with the amount of new information in each frame, not just the number of frames and dimensions of the data.

2. Robot perception. Robots have no choice but to learn online -their future input data (e.g. camera frames) are dependent on their previous predictions (i.e. motor actions). Not only does their data come in nonstationary temporal streams, but it typically comes from several sensors running at different rates. The camera may produce 1MB images at 30 frames/s, while the gyroscope might produce 1-byte readings at 1000 frames/s. It is not obvious, using current methods in deep learning, how we can integrate asynchronous sensory signals into a unified, trainable, latent representation, without undergoing the inefficient process of recomputing the function of the network every time a new signal arrives.

These examples point to the need for a training method where the amount of computation required to update the model scales with the amount of new information in the data, and not just the dimensionality of the data.

There has been a lot of work on increasing the computational efficiency of neural networks by quantizing neural weights or activations (see Section 4), but comparatively little work on exploiting redundancies in the data to reduce the amount of computation. O' Connor and Welling (2016b) , set out to exploit the temporal redundancy in video by having neurons only send their quantized changes in activation to downstream neurons, and having the downstream neurons integrate these changes over time. This approach (take the temporal difference, multiply by weights, temporally integrate) works for efficiently approximating the function of the network, but fails for training. The reason for this failure is that when the weights are functions of time, we no longer reconstruct the correct activation for the next layer. In other words, given a sequence of inputs x 0 ...x t with x 0 = 0 and weights w 1 ...w t : t τ =1 (x τ -x τ -1 ) • w τ = x t • w t unless w t = w 0 ∀t. Figure 2 describes the problem visually.

In this paper, we correct for this problem by encoding a mixture of two components of the layers activation x t : the proportional component k p x t , and the derivative component k d (x t -x t-1 ). When we invert this encoding scheme, we get a decoding scheme which corresponds to taking an exponentially decaying temporal average of past inputs. Interestingly, the resulting neurons begin to resemble models of biological spiking neurons, whose membrane potentials can approximately be modeled as an exponentially decaying temporal average of past inputs.

In this work, we present a scheme wherein the temporal redundancy of input data is used to reduce the computation required to train a neural network. We demonstrate this on the MNIST and Youtube-BB datasets. To our knowledge we are the first to create a neural network training algorithm which uses less computation as data becomes more temporally redundant.

2. METHODS:

We propose a coding scheme where neurons can represent their activations as a temporally sparse series of impulses. The impulses from a given neuron encode a combination of the value and the rate of change of the neuron's activation.

While our algorithm is designed to work efficiently with temporal data, we do not aim to learn temporal sequences in this work. We aim to efficiently approximate a function y t = f (x t ), where the current target y t is solely a function of the current input x t , and not previous inputs x 0 ...x t-1 . The temporal redundancy between neighbouring inputs x t-1 , x t will however be used to make our approximate computation of this function more efficient.

2.1. PRELIMINARY:

Throughout this paper we will use the notation

(f 3 • f 2 • f 1 )(x) = f 3 (f 2 (f 1 (x)

)) to denote function composition. We slightly abuse the notion of functions by allowing them to have an internal state which persists between calls. For example, we define the ∆ function in Equation 1as being the difference between the inputs in two consecutive calls (where persistent variable x last is initialized to 0). The Σ function, defined in Equation 2, returns a running sum of the inputs over calls. So we can write, for example, that when our composition of functions (Σ • ∆) is called with a sequence of input variables x 1 ...x t , then

(Σ•∆)(x t ) = x t , because (x 1 -x 0 )+(x 2 -x 1 )+...+(x t -x t-1 )| x0=0 = x t .

In general, when we write y t = f (x t ), where f is a function with persistent state, it will be implied that we have previously called f (x τ ) for τ ∈ [1, .., t -1] in sequence. Variable definitions that are used later will be highlighted in blue. While all terms are defined in the paper, we encourage the reader to refer to Appendix A for a complete collection of definitions and notations.

2.2. POSITION-DERIVATIVE (PD) ENCODING:

∆ :x → y; Persistent:

EQUATION

EQUATION

Suppose a neuron has time-varying activation x 1 ..x t .

Taking inspiration from Proportional-Integral-Derivative (PID) controllers, we can "encode" this activation at each time step as a combination of its current activation and change in activation as a t enc(x t ) = k p x t + k d (x t -x t-1 ), (see Equation 4). The parameters k p and k d determine what portion of our encoding represents the value of the activation and the rate of change of that value, respectively. In Section 4, we discuss how this is a form of Predictive Coding and in Appendix E, we discuss the effect our choices for these parameters have on the network.

To derive our decoding formula, we can simply solve for x t as x t = at+k d xt-1 kp+k d (Equation 5), such that (dec • enc)(x t ) = x t . Notice that Equation 5 corresponds to decaying the previous decoder state by some constant k d /(k p + k d ) followed by adding the input a t /(k p + k d ). We can expand this recursively to see that this corresponds to a temporal convolution a * κ where κ is a causal exponential kernel κ τ =

1 kp+k d k d k d +kp τ if τ ≥ 0 otherwise 0 .

2.3. QUANTIZATION:

Our motivation for the aforementioned encoding scheme is that we want a sparse signal which can be quantized into a low bitrate discrete signal. This will later be used to reduce computation. We can quantize our signal a t into a sparse, integer signal s t Q(a t ), where the quantizer Q is defined in Equation 3. Equation 3 implements a form of Sigma-Delta modulation, a method widely used in signal processing to approximately communicate signals at low bit rates (Candy and Temes, 1962) . We can show that Q( 1 .

x t ) = (∆ • R • Σ)(x t ) (see Supplementary Material Section C), where ∆ • R • Σ indicates

2.3.1. SPECIAL CASES:

We can compactly write the entire reconstruction function as x = (dec . This is visible in the bottom row of Figure 1 . This was the encoding scheme used in O' Connor and Welling (2016b) .

• ∆ • R • Σ • enc)(x t ). k p = 0: When k p = 0, we get dec(x t ) = (k -1 d • Σ)(x t ) and enc(x t ) = (k d • ∆)(x t ), so our reconstruction reduces to x = (k -1 d • Σ • ∆ • R • Σ • k d • ∆)(x t ). Because Σ • k d • ∆

k d = 0: In this case, dec(x t ) = k -1 p x t and enc(x t ) = k p x t so our encoding-decoding process becomes x = (k -1 p • ∆ • R • Σ • k p )(x t ).

Neither our encoder nor our decoder have any memory, and we take no advantage of temporal redundancy.

2.4. SPARSE COMMUNICATION BETWEEN LAYERS:

The purpose of our encoding scheme is to reduce computation by sparsifying communication between layers of a neural network. Our approach is to approximate the matrix-product as a series of additions, where the number of additions is inversely proportional to the sparsity of the input data. Suppose we are trying to compute the pre-nonlinearity activation of the first hidden layer, z t ∈ R dout , given the input activation, x t ∈ R din . We approximate z t as:

z t x t • w t ≈ xt • w t dec(Q(enc(x t ))) • w t dec(s t ) • w t ≈ dec(s t • w t ) ẑt where: x t , xt ∈ R din ; s t ∈ Z din ; w ∈ R din×dout ; z t , ẑt ∈ R dout (7)

The first approximation comes from the quantization (Q) of the encoded signal, and the second from the fact that the weights change over time, as explained in Figure 2 . The effects of these approximations are further explored in Appendix E. 

s t = N n=1 sign(s t,in ) • γ in : i n ∈ [1..d in ]

where γ in is a length-d in one-hot vector with element (γ in ) in = 1. The matrix product s t • w can then be decomposed into a series of row additions:

s t • w = N n=1 sign(s t,in ) • γ in • w = N n=1 sign(s t,in )γ in • w = N n=1 sign(s t,in ) • w in,• (8)

If we include the encoding, quantization, and decoding operations, our matrix product takes a total of 2d in + 2d out multiplications, and n |s t,n | • d out + 3d in + d out additions. Assuming the n |s t,n | • d out term dominates, we can say that the relative cost of computing ẑt vs z t is:

EQUATION

2.5. A NEURAL NETWORK:

We can implement this encoding scheme on every layer of a neural network. Given a standard neural net f nn consisting of alternating linear (•w l ) and nonlinear (h l ) operations, our network function f pdnn can then be written as:

EQUATION

We can use the same approach to approximately calculate our gradients to use in training. If we define our layer activations as ẑl (dec

• •w l • Q • enc)(x) if l = 1 otherwise (dec • •w l • Q • enc)(ẑ l-1

), and L (f pdnn (x), y), where is some loss function and y is a target, we can backpropagate the approximate gradients as:

∂L ∂ ẑl = ∂L ∂z L if l = L h l (ẑ l ) • dec back l • •w T l+1 • Q back l+1 • enc back l+1 ( ∂L ∂ ẑl+1 ) otherwise (12)

This can be implemented by either executing a (sparse) forward and backward pass at each time-step, or in an "event-based" manner, where the quantizers fire "events" whenever incoming events push their activations past a threshold, and these events are in turn sent to downstream neurons. For ease of implementation, we opt for the former in our code. Note that unlike in regular backprop, computing these forward and backward passes results in changes to the internal state of the enc, dec, and Q components.

2.6. PARAMETER UPDATES:

There is no use in having an efficient backward pass if the parameter updates are not also efficient. In a normal neural network trained with backpropagation and simple stochastic gradient descent, the 

x t = h l-1 (z l -1), its quantized version, xt = (Q • enc)(x t )

, and its reconstruction xt = dec(x t ). Middle: Another signal, representing the postsynaptic gradient of the error e = ∂L ∂z l , along with its quantized (ē) and reconstructed (ê) variants. Bottom: The true weight gradient ∂L ∂wt and the reconstruction gradient ∂L ∂wt . At the time of the spike in ēt , we have two schemes for efficiently computing the weight gradient that will be used to increment weight (see Section 2.6). The past scheme computes the area under x • ê since the last spike, and the future scheme computes the total future additional area due to the current spike.

parameter update for weight matrix w has the form w ← w -η ∂L ∂w where η is the learning rate. If w connects layer l -1 to layer l, we can write ∂L ∂w = x t ⊗ e t where x t h l-1 (z l-1,t ) ∈ R din is the presynaptic (layer l -1) activation, e t ∂L ∂z l,t ∈ R dout is the postsynaptic (layer l) backpropagating gradient and ⊗ is the outer product. So we require d in • d out multiplications to update the parameters for each sample.

We want a more efficient way to compute this product, which takes advantage of the sparsity of our encoded signals to reduce computation. We can start by applying our encoding-quantizing-decoding scheme to our input and error signals as

xt (Q • enc)(x t ) ∈ Z din and ēt (Q • enc)(e t ) ∈ Z dout ,

and approximate our true update as ∂L ∂w recon,t xt ⊗ êt where xt dec(x t ) and êt dec(ē t ). This does not do any good by itself, because the vectors involved in the outer product, xt and êt , are still not sparse. However, we can exactly compute the sum of this value over time using one of two sparse update schemes -past updates and future updates -which are depicted in Figure 3 . We give the formula for the Past and Future update rules in Appendix D, but summarize them here:

Past Updates: For a given synapse w i,j , if either the presynaptic neuron spikes (x ti = 0) or the postsynaptic neuron spikes (ē ti = 0), we increment the w i,j by the total area under xτ,i êτ,j since the last spike. We can do this efficiently because between the current time and the time of the previous spike, xτ,i êτ,j is a geometric sequence. Given a known initial value u, final value v, and decay rate r, a geometric sequence sums to u-v 1-r . The area calculated is shown in pink on the bottom row of Figure 3 , and one algorithm to calculate it is in Equation 19 in Appendix D.

Future Updates: Another approach is to calculate the Present Value of the future area under the integral from the current spike. This is depicted in the blue-gray area in Figure 3 , and the formula is in Equation 20 in Appendix D.

Finally, because the magnitude of our gradient varies greatly over training, we create a scheme for adaptively tuning our k p , k d parameters to match the running average of the magnitude of the data. This is described in detail in Appendix E.

2.7. RELATION TO STDP:

An extremely attentive reader might have noted that Equation 20 has the form of an online implementation of Spike-Timing Dependent Plasticity (STDP). STDP (Markram et al., 2012) emerged from neuroscience, where it was observed that synaptic weight changes appeared to be functions of the Using the quantized input signal x and error signal ē, and their reconstructions xt and êt as defined in the last section, we define a causal convolutional kernel

κ t = k β (k α ) t if t ≥ 0 otherwise 0 ,

where

k α = k d kp+k d , k β 1 kp+k d .

We can then define a "cross-correlation kernel"

g t = {κ t if t ≥ 0 otherwise κ -t } = k β (k α ) |t| : t ∈ Z

which defines the magnitude of a parameter update as a function of the difference in timing between pre-synaptic spikes from the forward pass and post-synaptic spikes from the backward pass. The middle plot of Figure 4 is a plot of g. We define our STDP update rule as:

EQUATION

We note that while our version of STDP has the same double-exponential form as the classic STDP rule observed in neuroscience (Markram et al., 2012) , our "presynaptic" spikes come from the forward pass while our "postsynaptic" spikes come from the backwards pass. STDP is not normally used to as a learning rule networks trained by backpropagation, so the notion of forward and backward pass with a spike-timing-based learning rule are new. Moreover, unlike in classic STDP, we do not have the property that sign of the weight change depends on whether the presynaptic spike preceded the postsynaptic spike.

In Section D in the supplementary material we show experimentally that while Equations ∂L ∂w recon , ∂L ∂w past , ∂L ∂w f uture , ∂L ∂w stdp may all result in different updates at different times, the rules are equivalent in that for a given set of pre/post-synaptic spikes x, ē, the cumulative sum of their updates over time converges exactly.

3. EXPERIMENTS:



3.1. TEMPORAL MNIST:

To evaluate our network's ability to learn, we train it on the standard MNIST dataset, as well as a variant we created called "Temporal MNIST". Temporal MNIST is simply a reshuffling of the MNIST dataset so that so that similar inputs (in terms of L2-pixel distance), are put together. Figure 6 shows several snippets of consecutive frames in the temporal MNIST dataset. We compare our Proportional-Derivative Net against a conventional Multi-Layer Perceptron with the same architecture (one hidden layer of 200 ReLU hidden units and a softmax output). The results are shown in Figure 5 . Somewhat surprisingly, our predictor slightly outperformed the MLP, getting 98.36% on the test set vs 98.25% for the MLP. We assume this improvement is due to the regularizing effect of the quantization. On Temporal MNIST, our network was able to converge with less computation than it required for MNIST (it used 32 • 10 12 operations for MNIST vs 15 • 10 12 for Temporal MNIST), but We now put the number of computational operations on the x-axis. We see that as a result our PDNN shifts to the left. Right: Because our network computes primarily with additions rather than multiplications. When we multiply our operation counts with the estimates of Horowitz (2014) for the computational costs of arithmethic operations (0.1pJ for 32-bit fixed-point addition vs 3.2pJ for multiplication), we can see that our algorithm would be at an advantage on any hardware where arithmetic operations were the computational bottleneck. ended up with a slightly worse test score when compared with the MLP (the PDNN achieved 97.99% vs 98.28% for the MLP). The slightly higher performance of the MLP on Temporal MNIST may be explained by the fact that the gradients on Temporal MNIST tend to be correlated across time-steps, so weights will tend to move in a single direction for a number of steps, which will interfere with the PDNN's ability to accurately track layer activations (see Figure 2 ). Appendix F contains a table of results with varying hyperparameters.

3.2. YOUTUBE VIDEO DATASET:

Next, we want to simulate the setting of CCTV cameras, discussed in Section 1, where we have a lot of data with only a small amount of new information per frame. In the absence of large enough public CCTV video datasets, we investigate the surrogate task of frame-based object classification on wild YouTube videos from the large, recently released Youtube-BB dataset Real et al. (2017) . Our subset consists of 358 Training Videos and 89 Test videos with 758,033 frames in total. Each video is labeled with an object in one of 24 categories.

We start from a VGG19 network (Simonyan and Zisserman, 2014): a 19-layer convolutional network pre-trained on imagenet. We replace the top three layer with three of our own randomly initialized layers, and train the network both as a spiking network, and as a regular network with backpropagation.

While training the entire spiking network end-to-end works, we choose to only train the top layers, in order to speed up our training time.

We compare our training scores and computation between a spiking and non-spiking implementation.

The learning curves in Figure 7 show that our spiking network performs comparably to a non-spiking network, and Figure 8 shows how the computation per frame of our spiking network decreases as we increase the frame rate (i.e. as the input data becomes more temporally redundant). Because our spiking network uses only additions, while a regular deep network does multiply-adds, we use the estimated energy-costs per op of Horowitz (2014) to compare computations to a single scale, which estimates the amount of energy required to do multiplies and adds in fixed-point arithmetic. Right: Learning curves with respect to computational energy using the conversion of Horowitz (2014) . The spiking network slightly outperforms the non-spiking baseline -we suspect that this is because the added noise of spiking acts as a regularizer. Figure 8 : We simulate different frame-rates by selecting every n'th frame. This plot shows our network's mean computation over several snippets of video, at varying frame rates. As our frame rate increases, the computation per-frame of our spiking network goes down, while with a normal network, it remains fixed.

4. RELATED WORK:

Noise-Shaping is a quantization technique that aims to increase the fidelity of signal reconstructions, per unit of bandwidth of the encoded signal, by quantizing the signal in such a way that the quantization noise is pushed into a higher frequency band which is later filtered out upon decoding. Sigma-Delta (also known as Delta-Sigma) quantization is a form of noise-shaping. Shin (2001) first suggested that biological neurons may be performing a form of noise shaping, and Yoon (2017) found standard spiking neuron models actually implement a form of Sigma-Delta modulation.

The encoding/decoding scheme we use in this paper can be seen as a form of Predictive Coding. Predictive coding is a lossless compression technique wherein the predictable parts of a signal are subtracted away so that just the unpredictable parts are transmitted. The idea that biological neurons may be doing some form of predictive coding was first proposed by Srinivasan et al. (1982) . In a predictive-coding neuron (unlike neurons commonly used in Deep Learning), there is a distinction between the signal that a neuron represents and the signal that it transmits. The neurons we use in this paper can be seen as implementing a simple form of predictive coding where the "prediction" is that the neuron maintains a decayed form of its previous signal -i.e. that pred(x t ) Chklovskii and Soudry (2012) suggest that the biological spiking mechanism may be thought of as consisting of a sigma-delta modulator within a predictive-coding circuit.

k d kp+k d x t-1 (See Appendix B for detail).

To our knowledge, none of the aforementioned work has yet been used in the context of deep learning.

There has been sparse but interesting work on merging the notions of spiking neural networks and deep learning. Diehl et al. (2015) found a way to efficiently map a trained neural network onto a spiking network. Lee et al. ( 2016) devised a method for training integrate-and-fire spiking neurons with backpropagation -though their neurons did not send a temporal difference of their activations. O'Connor and Welling (2016a) created a method for training event-based neural networks -but their method took no advantage of temporal redundancy in the data. Binas et al. (2016) and (O'Connor and Welling, 2016b) both took the approach of sending quantized temporal changes to reduce computation on temporally redundant data, but their schemes could not be used to train a neural network. Bohte et al. (2000) showed how one could apply backpropagation for training spiking neural networks, but it was not obvious how to apply the method to non-spiking data. Zambrano and Bohte (2016) developed a spiking network with an adaptive scale of quantization (which bears some resemblance to our tuning scheme described in Appendix E), and show that the spiking mechanism is a form of Sigma-Delta modulation, which we also use here. Courbariaux et al. (2015) showed that neural networks could be trained with binary weights and activations (we just quantize activations). Bengio et al. (2015) found a connection between the classic STDP rule (Figure 4 , right) and optimizing a dynamical neural network, although the way they arrived at an STDP-like rule was quite different from ours (they frame STDP as a way to minimze an objective based on the rate of change of the real-valued state of the network, whereas we use it approximately compute gradients based on spike-encodings of layer activations).

5. DISCUSSION:

We set out with the objective of reducing the computation in deep networks by taking advantage of temporal redundancy in data. We described a simple rule (Equation 4) for sparsifying the communication between layers of a neural network by having our neurons communicate a combination of their temporal change in activation, and the current value of their activation. We show that it follows from this scheme that neurons should behave as leaky integrators (Equation 5). When we quantize our neural activations with Sigma-Delta modulation, a common quantization scheme in signal processing, we get something resembling a leaky integrate-and-fire neuron. We derive efficient update rules for the weights of our network, and show these to be related to STDP -a learning rule first observed in neuroscience. Finally, we train our network, verify that it does indeed compute more efficiently on temporal data, and show that it performs about as well as a traditional deep network of the same architecture, but with significantly reduced computation. Finally, we show that our network can train on real video data.

The efficiency of our approach hinges on the temporal redundancy of our input data and neural activations. There is an interesting synergy here with the concept of slow-features (Wiskott, 1999) . Slow-Feature learning aims to discover latent objects that persist over time. If the hidden units were to specifically learn to respond to slowly-varying features of the input, the layers in a spiking implementation of such a network would have to communicate less often. In such a network, the tasks of feature-learning and reducing inter-layer communication may be one and the same.

Code is available at github.com/petered/pdnn.

A NOTATION:

Here we present a legend of notation used throughout this paper. While the paper is intended to be self-contained, the reader may want to consult this list if ever there is any doubt about the meaning of a variable used in the paper. Here we indicate the section in which each symbol is first used.

Section 2.1 ∆: A "temporal difference" operator. See Equation 1Σ: A "temporal integration" operator. See Equation 2Q: Sigma-Delta quantization. See Equation 3φ: The internal state variable of the quantizer Q. enc: An "encoding" operation, which takes a signal and encodes it into a combination of the signal's current value and its change in value since the last time step. See Equation 4. dec: A "decoding" operation, which takes an encoding signal and attempts to reconstruct the original signal that it was encoded from. If there was no quantization done on the encoded signal, the reconstruction will be exact, otherwise it will be an approximation. See Equation 5.

R::

The "rounding" operation, which simply rounds an input to the nearest integer value.

x t : Used throughout the paper to represent the value of a generic analog input signal at time t. In Sections 2.1, 2.2, and 2.3 it represents a scalar, and thereafter it represents a vector of inputs.

Section 2.2 k p , k d ∈ R + : Positive scalar coefficients used in the encoder and decoder, controlling how the extent to which the encoding is proportional to the input (k p ) vs proportional to the temporal difference of the input (k d ). xt dec(s t ): Used to represent the reconstruction of input x t " obtained by encoding, quantizing, and decoding x t .

Section 2.4 w t ∈ R din×dout is the value of a weight matrix at time t. z t

x t • w t ∈ R dout is the value of a pre-nonlinearity hidden layer activation in a non-spiking network at time t.

ẑt dec(Q(enc(x t )) • w t ) ∈ R dout

is the value of a pre-nonlinearity hidden layer activation in the spiking network at time t. It is an approximation of z t .

Section 2.5 (•w l ) indicates applying a function which takes the dot-product of the input with the l'th weight matrix: (•w l )(x) x • w l h l indicates an elementwise nonlinearity (e.g. a ReLU).

Q l indicates the quantization step applied at the l'th layer (because quantization has internal state, φ, and an associated layer dimension, we use the subscript to distinguish quantizers at different layers.) dec l , enc l are likewise the (stateful) encoding/decoding functions applied before/after layer l. ∂L ∂z l is the derivative of the loss with respect to the (pre-nonlinearity) activation of layer l. (•w T l ) indicates the dot product with the transpose of w l . This is simply backpropagation across a weight matrix: If u x • w l , then ∂L ∂x = ∂L ∂u • w T l ẑl is the approximation to the (pre-nonlinearity) activation to layer l (ie the output of dec l ), computed by the spiking network.

( h l (ẑ l )) is a function that performs an elementwise multiplication of the input with the derivative of nonlinearity h l evaluated at ẑl . This is simply backpropagation across a nonlinearity: If u h l (x),

then ∂L ∂x = ∂L ∂u h l (x) dec back l , enc back l , Q back l

serve the same functions as dec l , enc l , Q l , but for the backward pass.

∂L ∂ ẑl ∈ R d l Is our approximation to the derivative of the loss of our network with respect to ẑl , which is itself an approximation of the activation z l in a non-spiking implementation of the network.

Section 2.6

In the updates section we describe how we calculate the weight gradients in layer l. Because this description holds for any arbitrary layer, we get rid of the layer subscript and use the following notation: 9 ). Section 2.7 ∂L ∂w stdp,t is the gradient approximation taken using the STDP-type update. It also converges to the same value as ∂L ∂w recon,t when averaged out over time.

x t h l-1 (z l-1,t ) ∈ R

k α k d kp+k d ∈ [0, 1], k β 1 kp+k d ∈ R + : A reparametrization of k p , k d in

terms of the memory in our decoder k α and the scaling of our encoded signal (k β ). This reparametrization is also used when discussing the automatic tuning of k p , k d to match the dynamic range of our data in Appendix E B RELATION TO PREDICTIVE CODING Our encoding/decoding scheme is an instance of predictive coding -an idea imported from the signal processing literature into neuroscience by Srinivasan et al. (1982) wherein the power of a transmitted signal is reduced by subtracting away the predictable component of this signal before transmission, then reconstructing it after (This requires that the encoder and decoder share the same prediction model). Bharioke and Chklovskii (2015) formulate feedforward predictive coding as follows (with variables names changed to match the conventions of this paper):

a t x t -C f eedf orward (x t-1 , x t-2 , ...) (14) = x t - ∞ τ =1 ω τ x t-τ

In the case of Linear Predictive Coding (15)

Where the reconstruction is done by:

EQUATION

They go on to define "optimal" liner filter parameters [w 1 , w 2 , ...] that minimize the average magnitude of a t in terms of the autocorrelation and signal-to-noise ratio of x.

Our scheme defines:

EQUATION

)

So it is identical to feedforward predictive coding with

ω τ = k d kp+k d if τ = 1 0 otherwise up to a scaling constant of (k p + k d ).

In our case, the function of this additional constant is to determine the coarseness of the quantization.

From this relationship it is clear that this work could be extended to come up with more efficient predictive coding schemes which could further reduce computation by learning the temporal characteristics of the input signal.

C SIGMA-DELTA UNWRAPPING

Here we show that Q = ∆ • R • Σ, where Q, ∆, R, Σ are defined in Equations 3, 2, 6, 1, respectively.

From Equation 3 (Q) we can see that

y t ← round(x t + φ t-1 ) ∈ Z φ t ← φ t-1 + x t -y t ∈ R

Now we can unroll for y t and use the fact that if s ∈ Z then round(a + s) = round(a) + s:

y t = round(x t + φ t-1 ) = round(x t + φ t-2 + x t-1 -y t-1 ) = round(x t + x t-1 + φ t-2 ) -y t-1 = round(x t + x t-1 + φ t-2 ) -round(x t-1 + φ t-2 ) = round( t τ =1 x τ + 0 φ 0 ) - t-2 τ =0 y τ -round( t-1 τ =1 x τ + 0 φ 0 ) - t-2 τ =0 y τ = round( t τ =1 x τ ) -round( t-1 τ =1 x τ ) (18)

At which point it is clear that Q is identical to a successive application of a temporal summation, a rounding, and a temporal difference. That is why we say

Q = ∆ • R • Σ.

D UPDATE ALGORITHMS:

In Section 2.6, we visually describe what we call the "Past" and "Future" parameters updates. Here we present the algorithms for implementing these schemes.

To simplify our expressions in the update algorithms, we re-parametrize our k p , k d coefficients as Figure 9 : In Section 2.6 and 2.7, we described 4 different update rules ("Reconstruction", "Past", "Future", and "STDP"), and stated that while they do not necessarily produce the same updates at the same times, they produce the same result in the end. Here we demonstrate this empirically. We generate two random spike-trains representing the presynaptic input and the postsynaptic error signal to a single synapse and observe how the weight changes according to the different update rules. Top: A randomly generated presynaptic quantized signal x, along with its reconstruction x. Middle: A randomly generated postsynaptic quantized error signal ē, along with its reconstruction ê. Bottom:

EQUATION

The cumulative weight update arising from our four updates methods. "recon" is just t τ =1 xτ êτ , "past" and "future" are described in Section 2.6 and "STDP" is described in Section 2.7. Note that by the end all methods arrive at the same final-weight value. E TUNING k p , k d

E.1 CAUSES OF APPROXIMATION ERROR:

Equation 7 shows how we make two approximations when approximating z

t = x t • w t with ẑt = (dec • w • Q • enc)(x t ).

The first is the "nonstationary weight" approximation, arising from the fact that w changes in time. The second is the "quantization" approximation, arising from the quantization of x. Here we do a small experiment in which we multiply a time-varying scalar signal x t with a time-varying weight w t for many different values of k p , k d to understand the effects of k p , k d on our approximation error. The bottom-middle plot in Figure 10 shows that we enter a high-reconstructionerror regime (blue on plot) when k d is small (high quantization error), or when k d >> k p (high nonstationary-weight error). The bottom-right plot shows that blindly increasing k p and k d leads to representing the signal with many more spikes. Thus we need to tune hyperparameters to find the "sweet spot" where reconstruction error is fairly low but our encoded signal remains fairly sparse, keeping computational costs low. The cosine distance between the "true" signal and the approximation due to the quantization of x. Bottom Middle: The Cosine Distance between the "true" signal and the full approximation described in Equation 7. This shows why we need both k p and k d to be nonzero. Bottom Right: The Number of spikes in the encoded signal. In a neural network this would correspond to the number of weight-lookups required to compute the next layer's activation: dec(Q(enc(x)) w).

E.2 AN AUTO-TUNING SCHEME FOR k p , k d:

The smaller the magnitude of a signal, the more severely distorted it is by our quantizationreconstruction scheme. We can see that scaling a signal by K has the same effect on the quantized version of the signal, s t , as scaling k p and k d by K:

s t = (Q • enc kp,k d )(Kx t ) = Q(k p Kx t + k d (Kx t -Kx t-1 )) = Q(Kk p x t + Kk d (x t -x t-1 )) = (Q • enc Kkp,Kk d )(x t ).

The fact that the reconstruction quality depends on the signal magnitude presents a problem when training our network, because the error gradients tend to change in magnitude throughout training (they start large, and become smaller as the network learns). To keep our signal within the useful dynamic range of the quantizer, we apply a simple scheme to heuristically adjust k p and k d for the forward and backward passes separately, for each layer of the network. Instead of directly setting k p , k d as hyperparameters, we fix the ratio k α k d kp+k d , and adapt the scale k β 1 kp+k d to the magnitude of the signal. Our update rule for k β is: 

EQUATION

H INSTABILITY IN NEURAL NETWORK REPRESENTATIONS:

Figure 8 seems to show that computation doesn't quite approach zero as our frame-rate increases, but flat-lines at a certain point. We think this may have to do with the fact that hidden layer activations are not necessarily smoother in time than the input. We demonstrate this by taking 5 video snippets from the Youtube-BB dataset and running them through a (non-spiking) 19 layer VGGNet architectures (Simonyan and Zisserman, 2014) , which was pre-trained on ImageNet.

Given these 5 snippets, we measures how much the average relative change in layer activation |at-at-1| 2(|at|+|at-1|) varies as we increase our frame-rate, at various layer-depths. We simulate lower frame rates by skipping every N'th frame of video. (so for example to get a 10FPS frame rate we simply select every 3rd frame of the 30FPS video). For each selected frame rate, and for the given layers, we measure the average inter-frame change at various layers: Where: S = 5 is the number of video snippets we average over T is the number of frames in each snippet a t is the activation of a layer at time t n is the number of frames we are skipping over.

This shows something interesting. While our deeper layers do indeed show less relative change in activation over frames than our input/shallow layers, we note that as frame-rate increases, this seems to approach zero much more slowly than our input/shallow layers. This is a problem for our method, which relies on temporal smoothness in all layers (especially those hidden layers with large amounts of downstream units) to save computation. It suggests that methods for learning slow feature detectors -layers that are trained specifically to look for slowly varying features of the input, may be helpful to us. For increasing network depth (red=shallow, violet=deep) 

:



Figures and Tables:

t ) = sign(e t ) tpost tpre sign(x t ) sign(e t ) tpost tpre Classic STDP Rule

Figure4: Left: Our STDP rule, when both the input and error spikes have the same sign. Middle: Our STDP rule, when the input and error spikes have opposite signs. Right: The classic STDP ruleMarkram et al. (2012), where the weight update is positive when a presynaptic spike preceeds a postsynaptic spike, and negative otherwise.

Figure 5: Top Row: Results on MNIST. Bottom Row: Results on Temporal MNIST. Left Column: the training and test scores as a function of epoch. Middle:We now put the number of computational operations on the x-axis. We see that as a result our PDNN shifts to the left. Right: Because our network computes primarily with additions rather than multiplications. When we multiply our operation counts with the estimates ofHorowitz (2014) for the computational costs of arithmethic operations (0.1pJ for 32-bit fixed-point addition vs 3.2pJ for multiplication), we can see that our algorithm would be at an advantage on any hardware where arithmetic operations were the computational bottleneck.

Figure 6: Some samples from the Temporal-MNIST dataset. Each column shows a snippet of adjacent frames.

Figure7: Left: Learning Curves on the Youtube Dataset. Right: Learning curves with respect to computational energy using the conversion ofHorowitz (2014). The spiking network slightly outperforms the non-spiking baseline -we suspect that this is because the added noise of spiking acts as a regularizer.

t enc(x t ): Used to represent the encoded version of x t . Section 2.3 s t Q(a t ): Used to represent the quantized, encoded version of x t .

w, dec(enc(x) w)) w, dec(Q(enc(x))) w) w, dec(Q(enc(x)) w))

Figure 10: Top Left: A time varying signal x t , the quantized signal Q(enc(x t )), and the time-varying "weight" w t . Bottom Left: Compare the true product of these signals x t •w t with the dec(enc(x t )•w t ), which shows the effects of the non-stationary weight approximation, and dec(Q(enc(x t )) • w) which shows both approximations. Top Middle: The Cosine distance between the "true" signal x w and the approximation due to the nonstationary w, scanned over a grid of k p , k d values. Top Right:The cosine distance between the "true" signal and the approximation due to the quantization of x. Bottom Middle: The Cosine Distance between the "true" signal and the full approximation described in Equation7. This shows why we need both k p and k d to be nonzero. Bottom Right: The Number of spikes in the encoded signal. In a neural network this would correspond to the number of weight-lookups required to compute the next layer's activation: dec(Q(enc(x)) w).

Figure 11: 16 Frames from the Youtube-BB dataset. Each video annotated as having one of 24 objects in it. It also comes with annotated bounding-boxes, which we do not use in this study.

nt -a (n-1)t | 2(|a nt | + |a (n-1)t |)

Figure12: The average relative change in layer activation between frames, as frame-rate increases. For increasing network depth(red=shallow, violet=deep)

applying a temporal summation, a rounding, and a temporal difference, in series. If x t is temporally redundant and we set k p to be small, then |a t | 1∀t, and we can expect s t to consist of mostly zeros with a few 1's and -1's.We can now approximately reconstruct our original signal x t as xt dec(s t ) by applying our decoder, as defined in Equation 5. As our coefficients k p , k d become larger, our reconstructed signal xt should become closer to the original signal x t . We illustrate examples of encoded signals and their reconstructions for different k p , k d in Figure

Table content: <table/>

Figure 1: An example signal x t (blue), encoded with k p varying across rows and k d varying across columns. s t (black) is the quantized signal produced by the successive application of encoding (Equation4) and quantization (Equation3, where N indicates the total number of spikes. xt (orange) is the reconstruction of x t produced by applying Equation 5 to s t . One might, after a careful look at this figure, ask why we bother with the proportional (k p ) term at all? Figure2anticipates this question and answers it visually. Top: we generate random signals for a single scalar activation x t and scalar weight w t . Row 2: We efficiently approximate z t by taking the temporal difference, multiplying by w t then temporally integrating, to produce ẑt , as described in Section 2.4. As the weight w t changes over time, our estimate ẑ diverges from the correct value. Rows 3, 4: Introducing k p allows us to bring our reconstruction back in line with the correct signal.

Table content: <table><tr><td/><td>2.5</td><td/><td/><td/><td/><td/><td/><td/><td>x t</td><td>s t</td><td>x t</td></tr><tr><td>k p = 2.0</td><td>2.5 0.0</td><td/><td colspan="4">|x t x t | t = 0.16, N = 1031</td><td/><td colspan="2">|x t x t | t = 0.09, N = 1033</td><td>|x t x t | t = 0.042, N = 1039</td></tr><tr><td>k p = 0.01 k p = 0.1</td><td>2.5 0.0 2.5 0.0</td><td/><td colspan="4">|x t x t | t = 1.8, N = 5 |x t x t | t = 1.6, N = 52</td><td/><td colspan="2">|x t x t | t = 0.24, N = 10 |x t x t | t = 0.2, N = 53</td><td>|x t x t | t = 0.065, N = 40 |x t x t | t = 0.056, N = 69</td></tr><tr><td>k p = 0.0</td><td>0.0</td><td>0</td><td colspan="4">200 k d = 0 |x t x t | t = nan, N = 0 400</td><td>0</td><td colspan="2">200 k d = 1.0 |x t x t | t = 0.26, N = 12 400</td><td>0</td><td>200 k d = 4.0 |x t x t | t = 0.057, N = 37 400</td></tr><tr><td/><td colspan="3">kp = 0</td><td>0 2.5 0.0 2.5</td><td colspan="2">In all plots, kd = 1.0</td><td colspan="2">xt</td><td>wt</td><td>zt = xt wt</td><td>zt = deck pkd (Q(enck pkd (xt)) wt)</td></tr><tr><td/><td colspan="3">kp = 0.01 kp = 0.1</td><td>10 5 0 5 2.5 0.0 2.5</td><td/><td/><td/><td/></tr><tr><td/><td/><td/><td/><td/><td>0</td><td>250</td><td>500</td><td>750</td><td>1000 t</td><td>1250</td><td>1500</td><td>1750</td><td>2000</td></tr><tr><td colspan="10">Figure 2: The problem with only sending changes in activation (i.e. k p = 0) is that during training,</td></tr><tr><td colspan="10">weights change over time. xt = round(x t • k d )/k d , with no dependence on x t-1</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td/><td>all commute</td></tr><tr><td colspan="10">with one another, we can simplify this to xt = (k -1 d • R • k d )(x t ). so our decoded signal is</td></tr></table>

1. Computing z t takes d in • d out multiplications and (d in -1) • d out additions. The cost of computing ẑt , on the other hand, depends on the contents of s t . If the data is temporally redundant, s t ∈ Z din should be sparse, with total magnitude N i |s t,i |. s t can be decomposed into a sum of one-hot vectors

Table content: <table/>

Figure 3: A visualization of our efficient update schemes from Section 2.6. Top: A scalar signal representing a presynaptic neuron activation

Table content: <table><tr><td>x x x</td><td/></tr><tr><td>e e e</td><td>spike</td></tr><tr><td>w t = xtet w t = xtet</td><td>w t, past w t, future</td></tr><tr><td/><td>t</td></tr></table>

din here is defined as a shorthand for "the input to layer l". ∈ R dout is simply a shorthand for "the approximate backpropagated gradient at layer l" (x t and ēt ) are the encoded and quantized versions of signals (x t and e t ) (x t and êt ) are the reconstructed versions of signals (x t and e t ), taken from the quantized (x t and ēt ) xt ⊗ êt ∈ R din×dout is the approximate gradient of weight matrix w, as calculated by taking the outer product of the (input, error) reconstructions, x, ê. ∂L ∂w f uture,t are the gradient approximations at time t taken using the (past, future) approximation methods, defined in Appendix D. They are more efficient to calculate than ∂L ∂w recon,t but converge to the same value when averaged over time (i.e. lim T →∞

Table content: <table><tr><td colspan="3">∂L ∂ ẑl,t ∂L e t</td></tr><tr><td colspan="3">∂w recon,t</td></tr><tr><td/><td colspan="2">∂L ∂w past,t , 1 T</td><td>T t</td><td>∂L ∂w f uture,t =</td></tr><tr><td>1 T</td><td>T t</td><td>∂L ∂w recon,t (see Figure</td></tr></table>

Peer review (iclr2018):
 This paper presents a novel method for spike based learning that aims at reducing the needed computation during learning and testing when classifying temporal redundant data. This approach extends the method presented on Arxiv on Sigma delta quantized networks (Peter O’Connor and Max Welling. Sigma delta quantized networks. arXiv preprint arXiv:1611.02024, 2016b.). Overall, the paper is interesting and promising; only a few works tackle the problem of learning with spikes showing the potential advantages of such form of computing. The paper, however, is not flawless. The authors demonstrate the method on just two datasets, and effectively they show results of training only for Feed-Forward Neural Nets (the authors claim that “the entire spiking network end-to-end works” referring to their pre-trained VGG19, but this paper presents only training for the three top layers). Furthermore, even if suitable datasets are not available, the authors could have chosen to train different architectures. The first dataset is the well-known benchmark MNIST also presented in a customized Temporal-MNIST. Although it is a common base-line, some choices are not clear: why using a FFNN instead that a CNN which performs better on this dataset; how data is presented in terms of temporal series – this applies to the Temporal MNIST too; why performances for Temporal MNIST – which should be a more suitable dataset — are worse than for the standard MNIST; what is the meaning of the right column of Figure 5 since it’s just a linear combination of the GOps results. For the second dataset, some points are not clear too: why the labels and the pictures seem not to match (in appendix E); why there are more training iterations with spikes w.r.t. the not-spiking case. Overall, the paper is mathematically sound, except for the “future updates” meaning which probably deserves a clearer explanation. Moreover, I don’t see why the learning rule equations (14-15) are described in the appendix, while they are referred constantly in the main text. The final impression is that the problem of the dynamical range of the hidden layer activations is not fully resolved by the empirical solution described in Appendix D: perhaps this problem affects CCNs more than FFN. Finally, there are some minor issues here and there (the authors show quite some lack of attention for just 7 pages): - Two times “get” in “we get get a decoding scheme” in the introduction; - Two times “update” in “our true update update as” in Sec. 2.6; - Pag3 correct the capital S in 2.3.1 - Pag4 Figure 1 increase font size (also for Figure2); close bracket after Equation 3; N (number of spikes) is not defined - Pag5 “one-hot” or “onehot”; - in the inline equation the sum goes from n=1 to S, while in eq.(8) it goes from n=1 to N; - Eq(10)(11)(12) and some lines have a typo (a .cdot) just before some of the ws; - Pag6 k_{beta} is not defined in the main text; - Pag7 there are two “so that” in 3.1; capital letter “It used 32x10^12..”; beside, here, why do not report the difference in computation w.r.t. not-spiking nets? - Pag7 in 3.2 “discussed in 1” is section 1? - Pag14 Appendix E, why the labels don’t match the pictures; - Pag14 Appendix F, explain better the architecture used for this experiment.