Title: LOCALIZED RANDOM PROJECTIONS CHALLENGE BENCHMARKS FOR BIO-PLAUSIBLE DEEP LEARNING

Abstract: Similar to models of brain-like computation, artificial deep neural networks rely on distributed coding, parallel processing and plastic synaptic weights. Training deep neural networks with the error-backpropagation algorithm, however, is considered bio-implausible. An appealing alternative to training deep neural networks is to use one or a few hidden layers with fixed random weights or trained with an unsupervised, local learning rule and train a single readout layer with a supervised, local learning rule. We find that a network of leaky-integrate-andfire neurons with fixed random, localized receptive fields in the hidden layer and spike timing dependent plasticity to train the readout layer achieves 98.1% test accuracy on MNIST, which is close to the optimal result achievable with errorbackpropagation in non-convolutional networks of rate neurons with one hidden layer. To support the design choices of the spiking network, we systematically compare the classification performance of rate networks with a single hidden layer, where the weights of this layer are either random and fixed, trained with unsupervised Principal Component Analysis or Sparse Coding, or trained with the backpropagation algorithm. This comparison revealed, first, that unsupervised learning does not lead to better performance than fixed random projections for large hidden layers on digit classification (MNIST) and object recognition (CI-FAR10); second, networks with random projections and localized receptive fields perform significantly better than networks with all-to-all connectivity and almost reach the performance of networks trained with the backpropagation algorithm. The performance of these simple random projection networks is comparable to most current models of bio-plausible deep learning and thus provides an interesting benchmark for future approaches.

1. INTRODUCTION:

While learning a new task, synapses deep in the brain undergo task-relevant changes (Hayashi-Takagi et al., 2015) . These synapses are often many neurons downstream of sensors and many neurons upstream of actuators. Since the rules that govern such changes deep in the brain are poorly understood, it is appealing to draw inspiration from deep artificial neural networks (DNNs) (LeCun et al., 2015) . DNNs and the cerebral cortex process information in multiple layers of many neurons (Yamins & DiCarlo, 2016; Kriegeskorte, 2015) and in both, the artificial and the biological neural networks, learning depends on changes of synaptic strengths (Hebbian theory, Hebb (1949) ). However, learning rules in the brain are most likely different from the backpropagation algorithm (Crick, 1989; Marblestone et al., 2016; Rumelhart et al., 1986) . Furthermore, biological neurons communicate by sending discrete spikes as opposed to real-valued numbers used in DNNs. Differences like these suggest that there exist other, possibly equally powerful, algorithms that are capable to solve the same tasks by using different, more biologically plausible mechanisms. Thus, an important question in computational neuroscience is how to explain the fascinating learning capabilities of the brain with bio-plausible network architectures and learning rules. On the other hand, from a pure machine learning perspective there is increasing interest in neuron-like architectures with local learning rules, mainly motivated by the current advances in neuromorphic hardware (Nawrocki et al., 2016) .

Image recognition is a popular task to test the proposed models. Because of its relative simplicity and popularity, the MNIST dataset (28×28-pixel grey level images of handwritten digits, LeCun (1998)) is often used for benchmarking. Typical performances of existing models are around 97-99% classification accuracy on the MNIST test set (see section 2 and Table 8 ). This value lies in the region of the benchmarks for a large class of classical DNNs trained with backpropagation but without data-augmentation or convolutional layers (see table in LeCun (1998) ). Thus, accuracies around this value are assumed to be an empirical signature of backpropagation-like deep learning (Lillicrap et al., 2016; Sacramento et al., 2017) . It is noteworthy, however, that several of the most promising approaches that perform well on MNIST have been found to fail on harder tasks (Bartunov et al., 2018) .

An alternative to supervised training of all layers with backpropagation are fixed random weights, as proposed by general approximation theory (Barron, 1993) and the extreme learning field (Huang et al., 2006) , or unsupervised training in the first layers, combined with supervised training of a readout layer. Unsupervised methods are appealing since they can be implemented with local learning rules, see e.g. "Oja's rule" (Oja, 1982; Sanger, 1989) for principal component analysis or algorithms in Olshausen & Field (1997) ; Rozell et al. (2008) ; Liu & Jia (2012) ; Brito & Gerstner (2016) for sparse coding. A single readout layer can also be implemented with a local delta-rule (also called "perceptron rule"), which may be implemented by pyramidal spiking neurons with dendritic prediction of somatic spiking (Urbanczik & Senn, 2014) . Since it is pointless to simply stack multiple fully connected layers trained with principal component analysis or sparse coding (Olshausen & Field, 1997) we investigate here networks with a single hidden layer.

The main objective of this study was to see how far we can go with a single hidden layer and local learning rules in networks of spiking neurons. To support the design choices of the spiking model, we compared the classification performance of different rate networks: networks trained with backpropagation, networks where the hidden layer is trained with unsupervised methods, and networks with fixed random projections in the hidden layer. Since sparse connectivity is sometimes superior to dense connectivity (Litwin-Kumar et al., 2017; Bartunov et al., 2018) and successful convolutional networks leverage local receptive fields, we investigated also sparse connectivity between input and hidden layer, where each hidden neuron receives input only from a few neighboring pixels of the input image.

2. RELATED WORK:

In recent years, many bio-plausible approaches to deep learning have been proposed (see e.g. Marblestone et al. (2016) for a review). For achieving performances similar to deep learning methods, existing approaches usually use either involved architectures or elaborate mechanisms to approximate the backpropagation algorithm. Examples include the use of convolutional layers (Tavanaei & Maida (2016) ; Lee et al. (2018); Kheradpisheh et al. (2018) and table therein), dendritic computations (Hussain et al., 2014; Guergiuev et al., 2016; Sacramento et al., 2017) or approximations of the backpropagation algorithm such as feedback alignment (Lillicrap et al., 2016; Baldi et al., 2016; Nøkland, 2016; Samadi et al., 2017; Kohan et al., 2018; Bartunov et al., 2018) equilibrium propagation (Scellier & Bengio, 2017) , membrane potential based backpropagation (Lee et al., 2016) , restricted Boltzmann machines and deep belief networks (O'Connor et al., 2013; Neftci et al., 2014) , (localized) difference target propagation (Lee et al., 2015; Bartunov et al., 2018) , reinforcementsignal models like AuGMEnT (Rombouts et al., 2015) or approaches using predictive coding (Whittington & Bogacz, 2017) . Many models implement spiking neurons to stress bio-plausibility (Liu et al. (2016) ; Neftci et al. (2017) ; Kulkarni & Rajendran (2018); Wu et al. (2018) ; Liu & Yue (2018) and table therein) or for coding efficiency (O'Connor et al., 2017) . The conversion of DNNs to spiking neural networks (SNN) after training with backpropagation (Diehl et al., 2015) is a common technique to evade the difficulties of training with spikes. Furthermore, there are models including recurrent activity (Spoerer et al., 2017; Bellec et al., 2018) or even starting directly from realistic circuits (Delahunt & Kutz, 2018) . We refer to Table 8 for a list of current bio-plausible MNIST benchmark models.

3. RESULTS:

We study networks that consist of an input (l 0 ), one hidden (l 1 ) and an output-layer (l 2 ) connected by weight matrices W 1 and W 2 (Figure 1 ). Training the hidden layer weights W 1 with standard supervised training involves (non-local) error backpropagation using the transposed weight matrix W T 2 (Figure 1A ). In the bio-plausible network considered in this paper (Figure 1B ), the input-tohidden weights W 1 are either learned with an unsupervised method (Principal Component Analysis or Sparse Coding) or are fixed random projections. The unsupervised methods assume recurrent inhibitory weights V 1 between hidden units to implement competition.

3.1. SPIKING LOCALIZED RANDOM PROJECTIONS:

We first present the results with networks of leaky integrate-and-fire (LIF) neurons. The network architecture is as in Figure 1B , but without the recurrent connections V 1 . For implementing localized Random Projections (l-RP) in the hidden layer weights W 1 , we first chose the centers of the localized receptive fields at random positions in the input space and then randomly chose the weights therein, see Figure 1C . The receptive field patches span p × p pixels around their center position (we used p = 10 for the 28×28-pixel MNIST data). The output layer weights W 2 are trained with a supervised spike timing dependent plasticity (STDP) rule.

3.1.1. LIF AND STDP DYNAMICS:

The spiking dynamics follow the usual LIF equations (see methods A.4) and the readout weights W 2 evolve according to a supervised STDP delta rule using post-synaptic spike-traces tr i (t) and a post-synaptic target trace tgt i (t)

EQUATION

Thus, for a specific readout weight w 2,ij , the post-synaptic trace is updated at every post-synaptic spike time t f i and the weight is updated at every pre-synaptic spike time t f j . The target trace is used for feeding in the one-hot coded, supervisory signal for the MNIST classification into the output layer (l 2 ).

For a proof-of-principle and efficient parameter search we first investigate an LIF rate model. This rate model mimics the LIF dynamics by using the LIF activation function ϕ LIF as nonlinearity,

EQUATION

)

where u is the membrane potential, ∆ abs the refractory period, τ m the membrane time constant and ϑ the firing threshold of the LIF model. Furthermore, it employs the rate-version of the STDP delta rule Equation 1 (see methods section A.4 for details)

EQUATION

where tgtrate post i is the post-synaptic target rate, corresponding to the post-synaptic target trace tgt i (t) in Equation 1. We obtained similar spiking and weight dynamics when the readout weights W 2 were either directly trained with STDP or trained with the LIF rate model and then plugged into the spiking LIF network (as done in e.g. Diehl et al. (2015) ).

To illustrate the LIF and STDP dynamics, a toy example consisting of one pre-connected to one post-synaptic neuron was integrated for 650 ms. The pre-and post-synaptic membrane potentials show periodic spiking (Figure 2A ) which induces post-synaptic spike traces and corresponding weight changes (Figure 2B ), according to Equation 1. For the MNIST task, Figure 2C shows a raster plot for an exemplary training and testing protocol. During activity transients after pattern switches, learning is disabled until regular spiking is recovered. This is done, first, to ensure stability during activity transients (see Naud et al. (2008) and references therein) and second, to achieve decorrelation between the activities of subsequent patterns, as needed for stochastic gradient descent (SGD). During the testing period, learning is shut off permanently (see methods section A.4 for more details). 

3.1.2. CLASSIFICATION RESULTS FOR LIF l-RP:

When directly trained with the STDP rule in Equation 1 the spiking LIF l-RP model (n h = 5000 hidden units and patch size p = 10) reaches 98.1% test accuracy on MNIST. The corresponding LIF rate model reaches 98.5% test accuracy. Transferring weights learned with the LIF rate model into the spiking LIF model resulted in similar accuracies as the LIF rate model. Table 1 compares the performances of the rate and spiking LIF l-RP models with the reference algorithm l-BP, which is a rate model trained with backpropagation, see subsection 3.2 and subsection 3.3 (for same hidden layer size n h and patch size p). We can see that the spiking LIF model almost reaches the performance of the corresponding rate model. The remaining gap (0.4%) between rate and spiking LIF model presumably stems from transients and the shorter training time of the spiking model (only 10 6 compared to 10 7 iterations due to long simulation times). Both, the rate and spiking LIF model of l-RP achieve accuracies close to the backpropagation reference algorithm l-BP and certainly lie in the range of current bio-plausible MNIST benchmarks, i.e. 97-99% test accuracy (see section 2 and Table 8 ). Based on these numbers we conclude that the spiking LIF model of localized random projections using STDP is capable of learning the MNIST task to a level that is competitive with known benchmarks for spiking networks.

l-RP spiking LIF l-RP rate LIF l-BP 98.1 ± 0.04 98.5 ± 0.16 98.8 ± 0.1

Table 1 : Test accuracies (%) on MNIST for n h = 5000 hidden neurons and receptive field size p = 10. The reference algorithm l-BP is a rate model trained with backpropagation, see subsection 3.2 and subsection 3.3. The rate models were trained with 10 7 iterations (pattern presentations), the spiking LIF model with 10 6 iterations.

3.2. BENCHMARKING RATE MODELS TRAINED WITH UNSUPERVISED LEARNING AND:



BACKPROPAGATION:

To justify the design choices of the spiking model, we systematically investigated rate models with different methods to initialize or learn the hidden layer weights W 1 (see Figure 1 and methods subsection A.1 for details). To set these hidden layer weights, we use either one of the unsupervised methods Principal Component Analysis (PCA) or Sparse Coding (SC), or train only the readout layer W 2 and use fixed Random Projections (RP, as in subsection 3.1) for the hidden layer weights W 1 (see Figure 1B ). All these methods can be implemented with local, bio-plausible learning rules (Oja, 1982; Olshausen & Field, 1997) . As a reference and upper performance bound, we train networks with the same architecture with standard backpropagation (BP, see Figure 1A ). As a more bio-plausible approximation of BP, we include Feedback Alignment (FA, Lillicrap et al. (2016) ) which uses fixed random feedback weights for error-backpropagation (see methods subsection A.3 for further explanation). A Simple Perceptron (SP) without a hidden layer serves as a simplistic reference, since it corresponds to direct classification of the input.

The hidden-to-output weights W 2 are trained with standard stochastic gradient descent (SGD), using a one-hot representation of the class label as target. Since no error-backpropagation is needed for a single layer, the learning rule is local ("delta" or "perceptron"-rule, similar to Equation 3 of the LIF rate model). Therefore the system as a whole is bio-plausible in terms of online learning and synaptic updates using only local variables. For computational efficiency, we train first the hidden layer and then the output layer, however, both layers could be trained simultaneously.

We compared the test errors on the MNIST digit recognition data set for varying numbers of hidden neurons n h (Figure 3 ). The green PCA curve in Figure 3 The performance of RP (blue curve) for small numbers of hidden units (n h < d) is worse than for feature extractors like PCA and SC. Also for large hidden layers, performance improves only slowly with n h , which is in line with theory (Barron, 1993) and findings in the extreme learning field (Huang et al., 2006) . However, for large hidden layers sizes, RP outperforms SC. This suggests that the high dimensionality of the hidden layers is more important for reaching high performance than the features extracted by PCA or SC. Tests on the object recognition task CIFAR10 lead to the same conclusion, indicating that this observation is not entirely task specific (see subsection 3.3 for further analysis on CIFAR10). For all tested methods and hidden layer sizes, performance is significantly worse than the one reached with BP (black curve in Figure 3 ). In line with (Lillicrap et al., 2016) , we find that FA (cyan curve) performs as well as BP on MNIST.

Universal function approximation theory predicts lower bounds for the squared error that follow a power law with hidden layer size n h for both BP (O(1/n h )) and RP (O(1/n 2/d h ), where d is the input dimension Barron et al. (1994) ; Barron (1993) ). In the log-log-plot in Figure 3 this would correspond to a factor d/2 = 784/2 = 392 between the slopes of the curves of BP and RP, or at least a factor d eff /2 ≈ 10 using an effective dimensionality of MNIST (see methods A.1). We find a much faster decay of classification error in RP and a smaller difference between RP and BP slopes than suggested by the theoretical lower bounds.

3.3. LOCALIZED RANDOM RECEPTIVE FIELDS:

There are good reasons to reduce the connectivity from all-to-all to localized receptive fields (Figure 1C ): local connectivity patterns are observed in real neural circuits (Hubel & Wiesel, 1962) , proven useful theoretically (Litwin-Kumar et al., 2017) and empirically (Bartunov et al., 2018) , and successfully used in convolutional networks (CNNs). Even though this modification seems well justified from both biological and algorithmic sides, it reduces the generality of the algorithm to input data such as images where neighborhood relations between pixels (i.e. input dimensions) are important.

For random projections with localized receptive fields (l-RP), the centers of the patches were chosen at random positions in the input space and their weights where randomly fixed (as in subsection 3.1, see Figure 1C ). We tested different patch sizes of p × p pixels and found an optimum around p ≈ 10 which is more pronounced for large hidden layer sizes n h (see Figure 4A ). Note that p = 1 corresponds to resampling the data with random weights, and p = 28 recovers fully connected RP performance.

The main finding here is the significant improvement in performance using l-RP: the optimum around p ≈ 10 almost reaches BP performance for n h = 5000 hidden neurons (blue arrow in

B A:

Figure 4 : Effect of localized connectivity on MNIST. A Test error for localized Random Projections (l-RP), dependent on receptive field size p for different hidden layer sizes n h . The optimum at p = 10 is more pronounced for large hidden layer sizes. Full connectivity is equivalent to p = 28. Note the log-lin scale. B Localized receptive fields decrease test errors of Sparse Coding (SC), Random Projections (RP) and Backpropagation (BP). The effect is most significant for l-RP, which outperforms l-SC and almost approaches (l-)BP performance for large n h and p = 10 (blue arrow). All other reference lines as in Figure 3 . Note the log-log scale.

Figure 4B ). As expected l-RP and the LIF rate model of l-RP in subsection 3.1 perform equally well. To achieve a fair comparison BP and SC were also tested with localized receptive fields (l-BP, l-SC, see Figure 4B ). Also these algorithms seem to benefit from localized connectivity (also with an optimum for patch size p = 10), however, not as much as RP. This makes l-RP a strong competitor of SC (and also FA, see Figure 3 ) as a bio-plausible algorithm in the regime of large, overcomplete hidden layers n h > d.

Since classification performances of l-RP and l-BP are very close for layer sizes above n h = 5000, we investigated the misclassified MNIST digits for both algorithms. We find that 75% of the (≈ 125) misclassified digits of l-BP (n h = 5000) are contained in the misclassified ones of l-RP (n h = 5000). This means that in roughly 75% of the cases l-RP fails, also the reference algorithm l-BP fails, suggesting that these digits are particularly hard to recognize for networks with one hidden layer. We trained networks with up to n h = 100000 hidden neurons to test if (l-)RP can finally reach (l-)BP performance, since the latter saturates for large n h (see Figure 4B ). Indeed for simulations with n h = 100000 and p = 10, l-BP and l-RP performance was not significantly different any more, both being at 1.2% test error.

To test whether l-RP only works for the relatively simple MNIST data set (centered digits, noninformative margin pixels, no clutter, uniform features and perspective etc.) or generalizes to more difficult tasks, we applied it to the CIFAR10 data set (Krizhevsky, 2013) . We first reproduced a typical benchmark performance of a fully connected network with one hidden layer trained with standard BP (≈ 56% test accuracy, n h = 5000, see also Lin & Memisevic (2016) ). Again, l-RP outperforms the unsupervised methods PCA and l-SC in the case of large, overcomplete hidden layers (see Table 2 ). Furthermore, as on MNIST, classification performance increases for increasing hidden layer size n h and localized receptive fields perform better than full connectivity for all methods. Also on CIFAR10, l-RP comes close to the performance of the reference algorithm l-BP, however, the difference between l-RP and l-BP is larger than on MNIST. Given that state-of-the-art performance on the CIFAR10 dataset with deep convolutional neural networks is close to 98% (e.g. Real et al. (2018) ), the limitations of l-RP and the difference in difficulty between MNIST and CIFAR10 become apparent.

4. DISCUSSION:

The rules that govern plasticity of synapses deep in the brain remain elusive. In contrast to bioplausible deep learning based on approximations of the backpropagation algorithm, we focused 2500) PCs for MNIST (CIFAR10). All other methods use n h = 5000 hidden neurons and receptive field size p = 10. Note that CIFAR10 has d = 32×32×3 = 3072 input channels (the third factor is due to the color channels), MNIST only d = 28×28 = 784. The models were trained for 10 7 iterations (≈ 167 epochs).

here on training a readout layer with a supervised, local learning rule combined with a single hidden layer with either fixed random weights or trained with unsupervised, local learning rules.

To our surprise, randomly initialized fixed weights (RP) of large hidden layers lead to better classification performance than training them with unsupervised methods like PCA or sparse coding (SC). This implies that the inductive bias of PCA and sparse coding is not well suited for the task of digit classification and object recognition. It may be interesting to search for alternative unsupervised, local learning rules with a stronger inductive bias.

Replacing all-to-all connectivity with localized input filters is such an inductive bias that was already seen to be useful in other models (Bartunov et al., 2018) and proved to be particularly useful in conjunction with randomly initialized static weights. Already for a hidden layer size of 5000 neurons the performance of l-RP almost reaches the performance of backpropagation on MNIST. Furthermore, performance scaling with the number of hidden units n h was found to be orders of magnitudes better than the lower bound suggested by universal function approximation theory (Barron, 1993) .

Since we wanted to keep our models as simple as possible, we used online (no mini-batches) stochastic gradient descent (SGD) with a constant learning rate in all our experiments. There are many known ways to further tweak the final performance, e.g. with adaptive learning rate schedules or data augmentation, but our goal here was to demonstrate that even the simple model with localized random projections and spike timing dependent plasticity with a constant learning rate achieves results that are comparable with more elaborate approaches that use e.g. convolutional layers with weight sharing (Panda & Roy, 2016) , backpropagation approximations (Lee et al., 2016) , multiple hidden layers (Lillicrap et al., 2016) , dendritic neurons (Sacramento et al., 2017) , recurrence (Diehl & Cook, 2015) or conversion from rate to spikes (Diehl et al., 2015) .

Above 98% accuracy we have to take into account a saturating effect of the network training: better models will only lead to subtle improvements in accuracy. It is not obvious whether improvements are really a proof of having achieved deep learning or just the result of tweaking the models towards the peculiarities of the MNIST dataset (centered digits, non-informative margin pixels, no clutter, uniform features and perspective etc.). We observed that more challenging data sets such as CIFAR10 clearly highlight the limitations of l-RP and thus are better suited to test deep learning capabilities. We are aware that state-of-the-art deep learning has moved from MNIST to harder datasets, such as ImageNet (Deng et al., 2009) , long ago. Yet MNIST seems to be the current reference task for most bio-plausible deep learning models (see section 2 and Table 8 ).

In this paper we presented a new MNIST benchmark for bio-plausible spiking networks. Using localized random projections (l-RP) and STDP learning, our spiking LIF model reached 98.1% test accuracy on MNIST which lies within the range of current benchmarks for bio-plausible models for deep learning (see section 2 and Table 8 ). Our network model is particularly simple, i.e. it has only one trainable layer and does not depend on sophisticated architectural or algorithmic features (e.g. to approximate backpropagation). Instead it relies on the properties of high-dimensional localized random projections. We suggest that novel, progressive approaches to bio-plausible deep learning should significantly outperform the benchmark presented here.

A METHODS:



A.1 RATE NETWORK MODEL:

We use a 3-layer (input l 0 , hidden l 1 = l h and output l 2 ) feed-forward rate-based architecture with layer sizes (n 0 for input), n 1 (hidden) and n 2 (output, with n 2 = # classes). The layers are connected via weight matrices W 1 ∈ R n1×n0 and W 2 ∈ R n2×n1 and each neuron receives bias from the bias vectors b 1 ∈ R n1 and b 2 ∈ R n2 respectively (see Figure 1 ). The neurons themselves are nonlinear units with an element-wise, possibly layer-specific, nonlinearity a i = ϕ l (u i ). The feed-forward pass of this model thus reads

u l+1 = W l+1 u l + b l+1 a l+1 = ϕ l+1 (u l+1 ).

(4)

The simple perceptron (SP) only consists of one layer (l

2 , W 2 ∈ R n2×n0 , b 2 ∈ R n2

). The sparse coding (SC) model assumes recurrent inhibition within the hidden layer l 1 . This inhibition is not modeled by an explicit inhibitory population, as required by Dale's principle (Dale, 1935) , but direct, plastic, inhibitory synapses V 1 ∈ R n1×n1 are assumed between neurons in l 1 . Classification error variances in Figure 3 & Figure 4 are displayed as shaded, semi-transparent areas with the same colors as the corresponding curves. Their lower and upper bounds correspond to the 25% and 75% percentiles of at least 10 independent runs.

An effective dimensionality d eff of the MNIST data set can be obtained, e.g. via eigen-spectrum analysis, keeping 90% of the variance. We obtain values around d eff ≈ 20. The measure proposed in Litwin-Kumar et al. (2017) gives the same value d eff ≈ 20. Another measure is the crossing of the PCA curve with the Simple Perceptron line in Figure 3 at n h = 25(= d eff ). We checked that training a perceptron (1 hidden layer, n h = 1000, 10 7 iterations, ReLU, standard BP) on the first 25 PCs of MNIST leads to 1.7% test error (vs 1.5% test error on the full MNIST data). Together, these findings suggest that the MNIST dataset lies mostly in a low-dimensional linear subspace with d eff ≈ 25 d. The MNIST (& CIFAR10) data was rescaled to values in [0,1] and mean centered, which means that the pixel-wise average over the data was subtracted from the pixel values of every image. The code for the implementation of our rate network model will be available online upon acceptance.

A.2 UNSUPERVISED TECHNIQUES:



A.2.1 PRINCIPAL COMPONENT ANALYSIS (PCA):

In this paper we do not implement PCA learning explicitly as a neural learning algorithm but by a standard PCA algorithm (https://github.com/JuliaStats/MultivariateStats. jl). For d-dimensional data such algorithms output the values of the n ≤ d first principal components as well as the principal subspace projection matrix P ∈ R n×d . This matrix can directly be used as feedforward matrix W 1 in our network since the lines of P correspond to the projections of the data onto the single principal components. In other words each neuron in the hidden layer l 1 extracts another principal component of the data.

Since PCA is a linear model, biases b 1 were set to 0 and the nonlinearity was chosen linear, i.e. ϕ 1 (u) = u. With this, we can write the (trained) feed-forward pass of the first layer of our PCA model as follows:

a 1 = u 1 = W 1 • a 0 with W 1 = P (5)

Since the maximum number of PCs that can be extracted is the dimensionality of the data, n max = d, the number of neurons in the hidden layer n 1 is limited by d. This makes PCA unusable for overcomplete hidden representations as investigated for SC and RP.

Consistency between the used standard algorithm and neural implementations of PCA ("Sanger's" rule Sanger (1989) ) was checked by comparing the extracted PCs and visualizing the learned projections (lines of P) for the case of 30 extracted PCs, i.e. n = 30.

A.2.2 SPARSE CODING (SC):

For d-dimensional data, SC aims at finding a dictionary W ∈ R h×d of features that lead to an optimal representation a 1 ∈ R h which is sparse, i.e. has as few non-zero elements as possible. The corresponding optimization problem reads:

EQUATION

)

Since this is a nonlinear optimization problem with latent variables (hidden layer) it cannot be solved directly. Usually an iterative two step procedure is applied (akin to the expectation-maximization algorithm) until convergence: First optimize with respect to the activities a with fixed weights W.

Second, assuming fixed activities, perform a gradient step w.r.t to weights.

We implement a biologically plausible SC model using a 2-layer network with recurrent inhibition and local plasticity rules similar to the one in Brito & Gerstner (2016). For a rigorous motivation (and derivation) that such a network architecture can indeed implement sparse coding we refer to Olshausen & Field (1997) ; Zylberberg et al. (2011) ; Pehlevan & Chklovskii (2015) ; Brito & Gerstner (2016). We apply the above mentioned two step optimization procedure to solve the SC problem given our network model. The following two steps are repeated in alternation until convergence of the weights:

1. Optimizing the hidden activations:

We assume given and fixed weights W 1 and V 1 and ask for optimal hidden activations a 1 . Because of the recurrent inhibition V 1 the resulting equation for the hidden activities a 1 is nonlinear and implicit. To solve this equation iteratively, we simulate the dynamics of a neural model with time-dependent internal and external variables u 1 (t) and a 1 (t) respectively. The dynamics of the system is then given by Zylberberg et al. (2011); Brito & Gerstner (2016) :

EQUATION

In practice the dynamics is simulated for N iter = 50 iterations, which leads to satisfying convergence (change in hidden activations < 5%).

2.. Optimizing the weights::

Now the activities a 1 are kept fixed and we want to update the weights following the gradient of the loss function. The weight update rules are Hebbian-type local learning rules (Brito & Gerstner, 2016) :

EQUATION

• is a moving average (low-pass filter) with some time constant τ mav . At the beginning of the simulation (or after a new pattern presentation) τ mav is increased starting from 0 to τ mav during the first τ mav . The values of the lines of W 1 are normalized after each update, however this can also be achieved by adding a weight decay term. Additionally the values of V 1 are clamped to positive values after each update to ensure that the recurrent input is inhibitory. Also the diagonal of V 1 is kept at zero to avoid self-inhibition.

During SC learning, at every iteration, the variabes u 1 (t) and a 1 (t) are reset (to avoid transients) before an input is presented. Then for every of the N iterations, equation 7 is iterated for N iter steps and the weights are updated according to equation 8.

For comparison with localized RP (l-RP, see subsubsection A.2.3), a localized version of SC was implemented with the same initialization of W 1 as in l-RP. The usual SC learning rule equation 8 is applied and the localized connectivity is kept by clamping weights outside the receptive fields to zero. Lateral inhibition weights V 1 are initialized and learned as in normal SC (full competition is kept). For a detailed parameter list, see Table 3 .

A.2.3 RANDOM PROJECTIONS (RP)

For RP, the weight matrix W 1 between input and hidden layer is initialized randomly W 1 ∼ N (0, σ 2 ) with variance-preserving scaling: σ 2 ∝ 1/n 0 . The biases b 1 are initialized by sampling from a uniform distribution U([0, 0.1]) between 0 and 0.1. In practice we used the specific initialization

EQUATION

for RP (keeping weights fixed), SC, SP and also BP & RF (both layers with W 2 , b 2 and n 1 respectively). The initialization of the biases b was found to be uncritical in the range of [0,0.1]

For localized RP (l-RP), neurons in the hidden layer receive input only from a fraction of the input units called a receptive field. Receptive fields are chosen to form a compact patch over neighbouring pixels in the image space. For each hidden neuron a receptive field of size p × p (p ∈ N) input neurons is created at a random position in the input space. The weight values for each receptive field (rf) and the biases are initialized as:

EQUATION

were the optimization factor c = 3 was found empirically through a grid-search optimization of classification performance. For exact parameter values, see Table 4 .

A.3 CLASSIFIER & SUPERVISED REFERENCE ALGORITHMS:

The connections W 2 from hidden to output layer are updated by a simple delta-rule which is equivalent to BP in a single-layer network and hence is bio-plausible. For having a reference for our bio-plausible models (Figure 1B ), we compare it to networks with the same architecture (number of layers, neurons, connectivity) but trained in a fully supervised way with standard backpropagation (Figure 1A ). The forward pass of the model reads:

EQUATION

The error ẽL is calculated from the comparison of activations in the last layer a L with the (one-hot encoded) target activations tgt, with respect to the chosen loss function: mean squared error (MSE), ẽL = tgta L (14)

EQUATION

or softmax/cross-entropy loss (CE),

EQUATION

EQUATION

Classification results (on the test set) for MSE-and CE-loss were found to be not significantly different. Rectified linear units (ReLU) were used as nonlinearity ϕ(u l ) for all layers (MSE-loss) or for the first layer only (CE-loss).

In BP the weight and bias update is obtained by stochastic gradient descent, i.e. ∆W l,ij ∝ ∂L ∂W l,ij . The full BP algorithm for deep networks reads (Rumelhart et al., 1986) :

EQUATION

)

where stands for element-wise multiplication, ⊗ is the outer (dyadic) product, ϕ l (•) is the derivative of the nonlinearity and α is the learning rate. FA (Lillicrap et al., 2016) uses a fixed random matrix R l instead of the transpose of the weight matrix W l for the error backpropagation step in equation 19.

To allow for a fair comparison with l-RP, BP and FA were implemented with full connectivity and with localized receptive fields with the same initialization as in l-RP. During training with BP (or FA), the usual weight update equation 19 was applied to the weights in the receptive fields, keeping all other weights at zero. The exact parameter values can be found in Table 5 .

A.4 SPIKING IMPLEMENTATION:



A.4.1 LIF MODEL:

The spiking simulations were performed with a custom-made event-based leaky integrate-and-fire (LIF) integrator written in the Julia-language. For large network sizes, the exact, event-based integration can be inefficient due to a large frequency of events. To alleviate dramatic slow-down, an Euler-forward integration was added to the framework. For sufficiently small time discretization (e.g. ∆t ≤ 5 • 10 -2 ms for the parameters given in Table 6 ) the error of this approximate integration does not have negative consequences on the learning outcome. Consistent results were obtained using event-based and Euler-forward integration. The code of this framework will be available online upon acceptance.

The dynamics of the LIF network is given by:

τ m du i (t) dt = -u i (t) + RI i (t) with I i (t) = I f f i (t) + I ext i (t) = j,f w ij t -t f j + I ext i (t)

and the spiking condition:

EQUATION

)

where u i (t) is the membrane potential, τ m the membrane time-constant, R the membrane resistance, w ij are the synaptic weights, (t) = δ(t)/τ m (with τ m in seconds) is the post-synaptic potential evoked by a pre-synaptic spike arrival, ϑ i is the spiking threshold and u reset the reset potential after a spike. The input is split into a feed-forward (I f f (t)) and an external (I ext (t)) contribution. Each neuron in the input layer l 0 (n 0 = d) receives only external input I ext proportional to one pixel value in the data. To avoid synchrony between the spikes of different neurons, the starting potentials and parameters (e.g. thresholds) for the different neurons are drawn from a (small) range around the respective mean values.

We implement STDP using post-synaptic spike-traces tr i (t) and a post-synaptic target-trace tgt i (t). To train the network, we present patterns to the input layer and a target-trace to the output layer. The MNIST input is scaled by the input amplitude amp inp , the targets tgt(t) of the output layer are the one-hot-coded classes, scaled by the target amplitude amp tgt . Additionally, every neuron receives a static bias input I ext bias ≈ ϑ to avoid silent units in the hidden layer. Every pattern is presented as fixed input for a time T pat and the LIF dynamics as well as the learning evolves according to equation 20 and equation 21 respectively. To ensure stability during transients (see Naud et al. (2008) and references therein), learning is disabled after pattern switches for a duration of about T trans = 4τ m . With the parameters we used for the simulations (see Table 6 ), firing rates of single neurons in the whole network stayed below 1 kHz which was considered as a bio-plausible regime. For the toy example in Figure 2A & B we used static input and target with the parameters amp inp = 40, amp tgt = 5 (i.e. target trace = 0.005), ϑ mean = 20, σ ϑ = 0, τ m = 50, α = 1.2 • 10 -5 . For the raster plot in Figure 2C we used amp inp = 300, amp tgt = 300, ϑ mean = 20, σ ϑ = 0, τ m = 50, α = 1.2 • 10 -5 T pat = 50 ms, T trans = 100 ms.

A.4.2 LIF RATE MODEL:

The LIF dynamics can be mapped to a rate model described by the following equations:

EQUATION

)

with the (element-wise) LIF-activation function ϕ LIF (•) and the modified plasticity function g(•): 2015)). The reasons for the remaining difference in performance presumably lie in transients and single-spike effects that cannot be captured by the rate model. Also, the spiking network was only trained with 10 6 image presentations (compared to 10 7 for the rate model) due to long simulation times.

EQUATION

B PARAMETER TABLES:

In the following tables we use scientific E-notation XeY = X • 10 Y for better readability. For all simulations, we scaled the learning rate proportional to 1/n h for n h > 5000 to ensure convergence.

:



Figures and Tables:

Figure 1: The network model. A Training with Backpropagation (BP) through one hidden layer. B Architecture with unsupervised feature learning or fixed Random Projections (RP) in the first layer and a supervised classifier in the second layer. W stands for feed-forward, V for recurrent, inhibitory weights (Only used for unsupervised feature learning of the first layer weights W 1 ). C Illustration of fully connected and localized receptive fields of W 1 .

Figure 2: Spiking LIF and STDP dynamics. A Dynamics of the pre-and postsynaptic membrane potentials, spike-traces and the weight value (B) of a toy example with two neurons and one synapse. The weight decreases when the post-trace is above the post-target-trace (c.g. Equation 1 and subsection A.4). Both neurons receive static external input: I ext pre I ext post ≈ ϑ (spiking threshold). C Rasterplot of a network trained on MNIST, where every spike is marked with a dot. The background color indicates the corresponding layers: input (blue, n 0 = 100 neurons), hidden (green, n h = 100) and output (red, n 2 = 10). Bold vertical lines indicate pattern switches, thin lines indicate ends of transient phases (indicated by semi-transparency), during which learning is disabled. Left: Behaviour at the beginning of the training phase. Right: Testing period (learning off) after 10 4 iterations (presented patterns), which is 1/6 of an epoch. The output layer has started to learn useful, 1-hot encoded class predictions. A downsampled (12 × 12) version of MNIST was used for improved visibility. Other parameters, see appendix Appendix B.

τ

be obtained by integrating the STDP rule Equation 21 and taking the expectation. Most of the parameters of the spiking-and the LIF rate models can be mapped to each other directly (seeTabs. 6 & 7). The learningrate α must be adapted since the LIF weight change depends on the presentation time of a pattern T pat . In the limit of long pattern presentation times (T pat τ m , τ tr ), the transition from the learning rate of the LIF rate model (α) to the one of the spiking LIF model (α) isα = 1000 ms T pat [ms] • 1000 • α,(26)where the second factor comes from a unit change from Hz to kHz. It is also possible to train weight matrices computationally efficient in the LIF rate model and plug them into the spiking LIF model afterwards (as in e.g. Diehl et al. (

ends at the vertical line n h = d = 784 because the number of principal components (PCs), i.e. the number of hidden units n h , is limited by the input dimension d. Since the PCs span the subspace of highest variance, classification performance quickly improves when adding more PCs for small n h and then saturates for larger n h , crossing the (dotted) Simple Perceptron line at n Figure3: MNIST classification with rate networks, according to Figure1. The test error decreases for increasing hidden layer size n h for Principal Component Analysis (PCA), Sparse Coding (SC), fixed Random Projections (RP) and the fully supervised reference algorithms Backpropagation (BP) and Feedback Alignment (FA). The dashed dotted line at 90 % is chance level, the dotted line around 14 % is the performance of a Simple Perceptron without hidden layer. The vertical line marks the input dimension d = 784, i.e. the transition from under-to overcomplete hidden representations. Note the log-log scale.

Table content: <table/>

Test accuracies (%) on MNIST and CIFAR10. PCA uses full connectivity and 700 (

Table content: <table><tr><td/><td>PCA</td><td>l-SC</td><td>l-RP</td><td>l-BP</td></tr><tr><td>MNIST</td><td>90.3 ± 0.02</td><td>97.3 ± 0.03</td><td>98.4 ± 0.1</td><td>98.8 ± 0.1</td></tr><tr><td colspan="5">CIFAR10 22.6 ± 0.003 35.3 ± 0.004 52.0 ± 0.004 58.3 ± 0.002</td></tr></table>

(t), tgt i (t) = α • tgt post i (t)tr post i (t) .

Table content: <table><tr><td>with the plasticity function</td><td/></tr><tr><td colspan="2">g tr post i (22)</td></tr><tr><td>j</td><td>(21)</td></tr></table>

Peer review (iclr2019):
 This article compares different methods to train a two-layer spiking neural network (SNN) in a bio-plausible way on the MNIST dataset, showing that fixed localized random connections that form the hidden layer, in combination with a supervised local learning rule on the output layer can achieve close to state-of-the-art accuracy compared to other SNN architectures. The authors investigate three methods to train the first layer in an unsupervised way: principal component analysis (PCA) on the rates, sparse coding of activations, and fixed random local receptive fields. Each of the methods is evaluated on the one hand in a time-stepped simulator, using LIF neurons and on the other hand using a rate-approximated model which allows for faster simulations. Results are compared between each other and as reference with standard backpropagation and feedback alignment. The main finding is that localized random projections outperform other unsupervised ways of computing first layer features, and with many hidden neurons approaches backpropagation results. These results are summarized in Table 8, which compares results of the paper and other state-of-the-art and bio-plausible SNNs. PCA and sparse coding work worse on MNIST than local random projections, regardless if the network is rate-based, spike-based or a regular ANN trained with the delta rule. Feedback Alignment, although only meant for comparison, performs best of the algorithms investigated in this paper. In general the question how to train multi-layer spiking neural networks in a bio-plausible way is very relevant for computational neuroscience, and has attracted some attention from the machine learning community in recent years (e.g. Bengio et al. 2015, Scellier & Bengio 2016, Sacramento et al. 2018). It is therefore a suitable topic for ICLR. Of course the good performance of single-layer random projections is not surprising, because it is essentially the idea of the Extreme Learning Machine, and this concept has been well studied also for neuromorphic approaches (e.g. Yao & Basu, 2017), and versions with local receptive fields exist as well (Huang et al. 2015 "Local Receptive Fields Based Extreme Learning Machine"). While the comparison of different unsupervised methods on MNIST is somehow interesting, it fails to show any deeper insights because MNIST is a particularly simple task, and already the CIFAR 10 results are far away from the state-of-the-art (which is >96% using CNNs). Another interesting comparison that is missing is with clustering weights, which has shown good performance for CNNs e.g. in (Coates & Ng, 2012) or (Dundar et al. 2015), and is also unsupervised. The motivation is not 100% clear because the first experiment uses spikes, and shows a non-negligible difference to rate models (the authors claim it*s almost the same, but for MNIST differences of 0.5% are significant). All later results are purely about rate models. The authors apparently did not explore e.g. conversion techniques as in (Diehl et al. 2015) to make the spiking results match the rate versions better e.g. by weight normalization. I would rate the significance to the SNN community as average, and to the entire ICLR community as low. The significance would be higher if it was shown that this method scales to deeper networks or at least can be utilized in deeper architectures. Scrutinizing the possibilitites with random projections on the other hand could lead to more interesting results. But the best results here are obtained with 5000 neurons with 10x10 receptive fields on images of size 28x28, thus the representation is more than overcomplete, and of higher complexity than a convolution layer with 3x3 kernels and many input maps. Because the results provide only limited insights beyond MNIST I can therefore not support acceptance at ICLR. Pros: + interesting comparison of unsupervised feature learning techniques + interesting topic of bio-plausible deep learning Cons: - only MNIST, no indications if method will scale - results are not better than state-of-the-art Minor comments: The paper is generally well-written and structured, although some of the design choices could have been explained in more detail. Generally, it is not discussed if random connections have any advantage over other spiking models in terms of accuracy, efficiency or speed, besides the obvious fact that one does not have to train this layer. The title is a bit confusing. While it*s not wrong, I had to read it multiple times to understand what was meant. The first sentence in the caption for Fig. 2 is also confusing, mixing the descriptions of panel A and B. Also, in A membrane potentials are shown, but the post-membrane potential seems to integrate a constant current instead of individual spikes. Is this already the rate approximation of Eq. 2? Or is it because of the statement in the caption that they both receive very high external inputs. In general, the figures in panel A and B do not make the dynamics of the network or the supervised STDP much clearer. Principal Component Analysis and Sparse Coding are done algorithmically instead of using a sort of nonlinear Hebbian Learning as in Lillicrap 2016. It would have been interesting to see if this changes the comparatively bad results for PCA and SC. In Fig. 3, the curve in the random projections case is not saturated, maybe it would have been interesting to go above n_h = 5000. As there are 784 input neurons, a convolutional neural network with 7 filter banks already would have around 5000 neurons, but in this case each filter would be convolved over the whole image, while with random projections the filter only exists locally. In Eq. 1, the notation is a bit ambigous: The first delta-function seems to be the Dirac-delta for continuous t, while the second delta is a Kronecker-delta with discrete t. In A.1 and A.4.2 it is stated that the output of a layer is u_{t+1} = W u_t + b but I think in both cases it should be W a_t + b where a_t = phi(u_t). Otherwise, you just have a linear model and no activations. In Table 3, a typo: "eq. equation"