Year: 2018

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 (CIFAR10); 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.

heading: 1 INTRODUCTION
text: While learning a new task, synapses deep in the brain undergo task-relevant changes (HayashiTakagi 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.

heading: 2 RELATED WORK
text: 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.

heading: 3 RESULTS
text: We study networks that consist of an input (l0), one hidden (l1) and an output-layer (l2) connected by weight matrices W1 and W2 (Figure 1). Training the hidden layer weights W1 with standard
supervised training involves (non-local) error backpropagation using the transposed weight matrix WT2 (Figure 1A). In the bio-plausible network considered in this paper (Figure 1B), the input-tohidden weights W1 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 V1 between hidden units to implement competition.

heading: 3.1 SPIKING LOCALIZED RANDOM PROJECTIONS
text: 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 V1. For implementing localized Random Projections (l-RP) in the hidden layer weights W1, 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 W2 are trained with a supervised spike timing dependent plasticity (STDP) rule.

heading: 3.1.1 LIF AND STDP DYNAMICS
text: The spiking dynamics follow the usual LIF equations (see methods A.4) and the readout weights W2 evolve according to a supervised STDP delta rule using post-synaptic spike-traces tri(t) and a post-synaptic target trace tgti(t)
τtr dtri(t) dt = −tri(t) + ∑ f δ ( t− tfi ) ∆w2,ij = α · ( tgtposti (t)− tr post i (t) ) δ ( t− tfj ) . (1)
Thus, for a specific readout weight w2,ij , the post-synaptic trace is updated at every post-synaptic spike time tfi 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 (l2).
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,
rate(u) = ϕLIF (u) = [ ∆abs − τm ln ( 1− ϑ
u
)]−1 , (2)
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)
∆wij = α̃ · rateprej · ( tgtrateposti − rate post i ) , (3)
where tgtrateposti is the post-synaptic target rate, corresponding to the post-synaptic target trace tgti(t) in Equation 1. We obtained similar spiking and weight dynamics when the readout weights W2 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 (nh = 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 nh 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 106 compared to 107 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.

heading: 3.2 BENCHMARKING RATE MODELS TRAINED WITH UNSUPERVISED LEARNING AND BACKPROPAGATION
text: 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 W1 (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 W2 and use fixed Random Projections (RP, as in subsection 3.1) for the hidden layer weights W1 (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 W2 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 nh (Figure 3). The green PCA curve in Figure 3 ends at the vertical line nh = d = 784 because the number of principal components (PCs), i.e. the number of hidden units nh, 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 nh and then saturates for larger nh, crossing the (dotted) Simple Perceptron line at nh = 25 PC hidden neurons. This intersection and other measures of effective dimensionality (see methods subsection A.1) suggest that the MNIST dataset lies mostly in a low-dimensional linear subspace with deff ≈ 25 d. SC performance (red curve) starts at a higher test error but improves as quickly with nh as PCA. With overcomplete representations (nh > d), the network achieves a remarkable classification performance of around 96 % test accuracy. This suggests that the sparse representation and the features extracted by SC are indeed useful for classification, especially in the overcomplete case.
The performance of RP (blue curve) for small numbers of hidden units (nh < d) is worse than for feature extractors like PCA and SC. Also for large hidden layers, performance improves only slowly with nh, 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 nh for both BP (O(1/nh)) and RP (O(1/n2/dh ), 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 deff/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.

heading: 3.3 LOCALIZED RANDOM RECEPTIVE FIELDS
text: 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 nh (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 nh = 5000 hidden neurons (blue arrow in
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 (lBP, 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 nh > d.
Since classification performances of l-RP and l-BP are very close for layer sizes above nh = 5000, we investigated the misclassified MNIST digits for both algorithms. We find that 75% of the (≈ 125) misclassified digits of l-BP (nh = 5000) are contained in the misclassified ones of l-RP (nh = 5000). This means that in roughly 75% of the cases l-RP fails, also the reference algorithm lBP fails, suggesting that these digits are particularly hard to recognize for networks with one hidden layer. We trained networks with up to nh = 100000 hidden neurons to test if (l-)RP can finally reach (l-)BP performance, since the latter saturates for large nh (see Figure 4B). Indeed for simulations with nh = 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, nh = 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 nh 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.

heading: 4 DISCUSSION
text: 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
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 nh 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.

heading: A METHODS

heading: A.1 RATE NETWORK MODEL
text: We use a 3-layer (input l0, hidden l1 = lh and output l2) feed-forward rate-based architecture with layer sizes (n0 for input), n1 (hidden) and n2 (output, with n2 = # classes). The layers are connected via weight matrices W1 ∈ Rn1×n0 and W2 ∈ Rn2×n1 and each neuron receives bias from the bias vectors b1 ∈ Rn1 and b2 ∈ Rn2 respectively (see Figure 1). The neurons themselves are nonlinear units with an element-wise, possibly layer-specific, nonlinearity ai = ϕl(ui). The feed-forward pass of this model thus reads
ul+1 = Wl+1ul + bl+1 al+1 = ϕl+1(ul+1). (4)
The simple perceptron (SP) only consists of one layer (l2, W2 ∈ Rn2×n0 , b2 ∈ Rn2 ). The sparse coding (SC) model assumes recurrent inhibition within the hidden layer l1. This inhibition is not modeled by an explicit inhibitory population, as required by Dale’s principle (Dale, 1935), but direct, plastic, inhibitory synapses V1 ∈ Rn1×n1 are assumed between neurons in l1. 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 deff of the MNIST data set can be obtained, e.g. via eigen-spectrum analysis, keeping 90% of the variance. We obtain values around deff ≈ 20. The measure proposed in Litwin-Kumar et al. (2017) gives the same value deff ≈ 20. Another measure is the crossing of the PCA curve with the Simple Perceptron line in Figure 3 at nh = 25(= deff). We checked that training a perceptron (1 hidden layer, nh = 1000, 107 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 deff ≈ 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.

heading: A.2 UNSUPERVISED TECHNIQUES

heading: A.2.1 PRINCIPAL COMPONENT ANALYSIS (PCA)
text: 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 ∈ Rn×d. This matrix can directly be used as feedforward matrix W1 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 l1 extracts another principal component of the data.
Since PCA is a linear model, biases b1 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:
a1 = u1 = W1 · a0 with W1 = P (5)
Since the maximum number of PCs that can be extracted is the dimensionality of the data, nmax = d, the number of neurons in the hidden layer n1 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.

heading: A.2.2 SPARSE CODING (SC)
text: For d-dimensional data, SC aims at finding a dictionary W ∈ Rh×d of features that lead to an optimal representation a1 ∈ Rh which is sparse, i.e. has as few non-zero elements as possible. The corresponding optimization problem reads:
Wopt, aopt1 = argmin L(W, a1)
L(W, a1) = 1
2 ‖a0 −W>a1‖22 + λ‖a1‖1. (6)
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 W1 and V1 and ask for optimal hidden activations a1. Because of the recurrent inhibition V1 the resulting equation for the hidden activities a1 is nonlinear and implicit. To solve this equation iteratively, we simulate the dynamics of a neural model with time-dependent internal and external variables u1(t) and a1(t) respectively. The dynamics of the system is then given by Zylberberg et al. (2011); Brito & Gerstner (2016):
τu du1(t) dt = −u1(t) + (W1a0(t)− V1a1(t))
a1(t) = ϕ(u1(t)) (7)
In practice the dynamics is simulated for Niter = 50 iterations, which leads to satisfying convergence (change in hidden activations < 5%).
2. Optimizing the weights: Now the activities a1 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):
∆W1,ji = αw · a0,i · a1,j ∆V1,jk = αv · a1,k · (a1,j − 〈a1,j〉) (8)
〈·〉 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 W1 are normalized after each update, however this can also be achieved by adding a weight decay term. Additionally the values of V1 are clamped to positive values after each update to ensure that the recurrent input is inhibitory. Also the diagonal of V1 is kept at zero to avoid self-inhibition.
During SC learning, at every iteration, the variabes u1(t) and a1(t) are reset (to avoid transients) before an input is presented. Then for every of the N iterations, equation 7 is iterated for Niter 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 W1 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 V1 are initialized and learned as in normal SC (full competition is kept). For a detailed parameter list, see Table 3.

heading: A.2.3 RANDOM PROJECTIONS (RP)
text: For RP, the weight matrix W1 between input and hidden layer is initialized randomly W1 ∼ N (0, σ2) with variance-preserving scaling: σ2 ∝ 1/n0. The biases b1 are initialized by sampling from a uniform distribution U([0, 0.1]) between 0 and 0.1. In practice we used the specific initialization
W1 ∼ N (0, σ2) σ2 = 1
100 n0 b1 ∼ U([0, 0.1]) (9)
for RP (keeping weights fixed), SC, SP and also BP & RF (both layers with W2,b2 and n1 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:
W1,rf ∼ N (0, σ2rf) σ2rf = c
100 p (10)
b1 ∼ U([0, 0.1]) (11) were the optimization factor c = 3 was found empirically through a grid-search optimization of classification performance. For exact parameter values, see Table 4.

heading: A.3 CLASSIFIER & SUPERVISED REFERENCE ALGORITHMS
text: The connections W2 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:
ul+1 = Wl+1ul + bl+1 (12) al+1 = ϕl+1(ul+1) (13)
The error ẽL is calculated from the comparison of activations in the last layer aL with the (one-hot encoded) target activations tgt, with respect to the chosen loss function: mean squared error (MSE),
ẽL = tgt− aL (14)
LMSE = 1
2 ‖tgt− aL‖22 (15)
or softmax/cross-entropy loss (CE),
p = softmax (aL) (16) ẽL = tgt− p (17)
LCE = − nL∑ i=1 tgti · log (pi) (18)
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 ϕ(ul) 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. ∆Wl,ij ∝ ∂L∂Wl,ij . The full BP algorithm for deep networks reads (Rumelhart et al., 1986):
eL = ϕ′L(uL) ẽL el−1 = ϕ′l−1(ul) W>l el
∆Wl = α · el ⊗ al−1 ∆bl = α · el (19)
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 Rl 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.

heading: A.4 SPIKING IMPLEMENTATION

heading: A.4.1 LIF MODEL
text: 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 dui(t)
dt = −ui(t) +RIi(t)
with Ii(t) = I ff i (t) + I ext i (t) = ∑ j,f wij ( t− tfj ) + Iexti (t)
and the spiking condition: if ui(t) ≥ ϑi: ui → ureset (20)
where ui(t) is the membrane potential, τm the membrane time-constant,R the membrane resistance, wij 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 ureset the reset potential after a spike. The input is split into a feed-forward (Iff (t)) and an external (Iext(t)) contribution. Each neuron in the input layer l0 (n0 = d) receives only external input Iext 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 tri(t) and a post-synaptic target-trace tgti(t).
τtr dtri(t) dt = −tri(t) + ∑ f δ ( t− tfi ) ∆wij = g ( trposti (t), tgt post i (t) ) δ ( t− tfj ) (21)
with the plasticity function g ( trposti (t), tgti(t) ) = α · ( tgtposti (t)− tr post i (t) ) . (22)
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 ampinp, the targets tgt(t) of the output layer are the one-hot-coded classes, scaled by the target amplitude amptgt. Additionally, every neuron receives a static bias input Iextbias ≈ ϑ to avoid silent units in the hidden layer. Every pattern is presented as fixed input for a time Tpat 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 Ttrans = 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 ampinp = 40, amptgt = 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 ampinp = 300, amptgt = 300, ϑmean = 20, σϑ= 0, τm = 50, α = 1.2 · 10−5 Tpat = 50 ms, Ttrans = 100 ms.

heading: A.4.2 LIF RATE MODEL
text: The LIF dynamics can be mapped to a rate model described by the following equations:
ul = Wlul−1 +RIext
al = ϕLIF (ul) ∆wij = g̃ ( aprej , a post i , tgt post i ) (23)
with the (element-wise) LIF-activation function ϕLIF(·) and the modified plasticity function g̃(·):
ϕLIF (uk) = [ ∆abs − τm ln ( 1− ϑk
uk
)]−1 (24)
g̃ ( aprej , a post i , tgt post i ) = α̃ · aprej · ( tgtposti − a post i ) (25)
The latter can 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 (see Tabs. 6 & 7). The learningrate α must be adapted since the LIF weight change depends on the presentation time of a pattern Tpat. In the limit of long pattern presentation times (Tpat τm, τtr), the transition from the learning rate of the LIF rate model (α̃) to the one of the spiking LIF model (α) is
α = 1000 ms Tpat [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. (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 106 image presentations (compared to 107 for the rate model) due to long simulation times.

heading: B PARAMETER TABLES
text: In the following tables we use scientific E-notation XeY = X · 10Y for better readability. For all simulations, we scaled the learning rate proportional to 1/nh for nh > 5000 to ensure convergence.

heading: C BIO-PLAUSIBLE MNIST BENCHMARKS

Review summaries:
In this work authors benchmark a biologically plausible network architecture for image classification. The employed architecture consists of one hidden layer, where input to hidden layer weights W1 are either trained with PCA or sparse coding, or are kept fixed after random initialization. The output layer units are modeled as leaky integrate-and-fire (LIF) neurons and hidden to output connections W2 are tuned using a rate model that mimics STDP learning dynamics in the LIF neurons. The authors compare classification results on MNIST and CIFAR10 datasets, where they also include results of an equivalent feed-forward network that is trained with standard error backpropagation.

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.