Biologically plausible deep learning -- but how far can we go with shallow networks?
Bernd Illing Wulfram Gerstner Johanni Brea

1
.001
Biologically plausible deep learning – but how far can we go with shallow networks?
Illing et al.


mode = title]Biologically plausible deep learning – but how far can we go with shallow networks?                      

1]Bernd Illing[]
[1]

bernd.illing@epfl.ch




[1]School of Computer and Communication Science & School of Life Science, EPFL, 1015 Lausanne, Switzerland

1]Wulfram Gerstner[]
1]Johanni Brea[]

[cor1]Corresponding author

[S U M M A R Y]
Training deep neural networks with the error backpropagation algorithm 
 is considered implausible from a biological perspective. 
 Numerous recent publications suggest elaborate models for biologically plausible variants of deep learning, typically defining success as reaching around 98% test accuracy on the MNIST data set.
Here, we investigate how far we can go on digit (MNIST) and object (CIFAR10) classification with biologically plausible, local learning rules in a network with one hidden layer and a single readout layer. 
The hidden layer weights are either fixed (random or random Gabor filters) or trained with unsupervised methods (Principal/Independent Component Analysis or Sparse Coding) that can be implemented by local learning rules. 
The readout layer is trained with a supervised, local learning rule.
We first implement these models with rate neurons. 
This comparison reveals, first, that unsupervised learning does not lead to better performance than fixed random projections or Gabor filters for large hidden layers. 
Second, networks with localized receptive fields perform significantly better than networks with all-to-all connectivity and can reach backpropagation performance on MNIST.
We then implement two of the networks - fixed, localized, random & random Gabor filters in the hidden layer - with spiking leaky integrate-and-fire neurons and spike timing dependent plasticity to train the readout layer.  
These spiking models achieve > 98.2% test accuracy on MNIST,  
which is close to the performance of 
rate networks with one hidden layer trained with backpropagation.
The performance of our shallow network models is comparable to most current biologically plausible models of deep learning. 
Furthermore, our results with a shallow spiking network provide an important reference and suggest the use of datasets other than MNIST for testing the performance of future models of biologically plausible deep learning.



Deep learning
 Local learning rules
 Random Projections
 Unsupervised Feature Learning
 Spiking Networks
 MNIST
 CIFAR10



[
    [
    December 30, 2023
=====================






§ INTRODUCTION




While learning a new task, synapses deep in the brain undergo task-relevant changes <cit.>.
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) <cit.>.





DNNs and the cerebral cortex share that information is processed in multiple layers of many neurons <cit.> and that learning depends on changes of synaptic strengths <cit.>. 
However, learning rules in the brain are most likely different from the backpropagation algorithm <cit.>.
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 nearly 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 biologically plausible network architectures and learning rules.
Moreover 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 <cit.>.


Image recognition is a popular task to test the performance of neural networks.
Because of its relative simplicity and popularity, the MNIST dataset (28×28-pixel grey level images of handwritten digits, <cit.>) is often used for benchmarking.
Typical performances of existing models are around 97-99% classification accuracy on the MNIST test set (see <ref> and <ref>).
Since the performances of many classical DNNs trained with backpropagation (but without data augmentation or convolutional layers, see table in <cit.>) also fall in this region, accuracies around these values are assumed to be an empirical signature of backpropagation-like deep learning <cit.>.
It is noteworthy, however, that several of the most promising approaches that perform well on MNIST have been found to fail on harder tasks <cit.> or at least need major modifications to scale to deeper networks <cit.>.

There are two obvious alternatives to supervised training of all layers with backpropagation. 
The first one is to fix weights in the first layer(s) at random values , as proposed by general approximation theory <cit.> and the extreme learning field <cit.>. The second alternative is unsupervised training in the first layer(s). 
In both cases, only the weights of a readout layer are learned with supervised training.
Unsupervised methods are appealing since they can be implemented with local learning rules, see e.g. “Oja's rule” <cit.> for principal component analysis, nonlinear extensions for independent component analysis <cit.> or algorithms in <cit.> for sparse coding.
A single readout layer can be implemented with a local rule as well. A candidate is the delta-rule (also called “perceptron rule”), which may be implemented by pyramidal spiking neurons with dendritic prediction of somatic spiking <cit.>. 
Since straightforward stacking of multiple fully connected layers of unsupervised learning does not reveal more complex features <cit.> we focus here on networks with a single hidden layer (see also <cit.>).

The main objective of this study is to see how far we can go with networks with a single hidden layer and biologically plausible, local learning rules, preferably using spiking neurons. 
To do so we first compare the classification performance of different rate networks: networks trained with backpropagation, networks with fixed random projections or random Gabor filters in the hidden layer and networks where the hidden layer is trained with unsupervised methods (<ref>).
Since sparse connectivity is sometimes superior to dense connectivity <cit.> and successful convolutional networks leverage local receptive fields, we investigate sparse connectivity between input and hidden layer, where each hidden neuron receives input only from a few neighboring pixels of the input image (<ref>).
Finally we implement the simplest, yet promising and biologically plausible models - localized random projections and random Gabor filters - with spiking leaky integrate-and-fire neurons and spike timing dependent plasticity (<ref>). We discuss the performance and implications of this simplistic model with respect to current models of biologically plausible deep learning.





§ RELATED WORK








In recent years, many biologically plausible approaches to deep learning have been proposed, see e.g. <cit.> for reviews.
Existing approaches usually use either involved architectures or elaborate mechanisms to approximate the backpropagation algorithm. Examples include the use of convolutional layers 

<cit.> (and tables therein),
 dendritic computations

 <cit.> 
 or backpropagation approximations such as feedback alignment

<cit.> 
equilibrium propagation <cit.>, 
membrane potential based backpropagation <cit.>, 
restricted Boltzmann machines and deep belief networks <cit.>, 
(localized) difference target propagation <cit.>, 
 using reinforcement-signals <cit.> or approaches using predictive coding <cit.>.
Many models implement spiking neurons to stress bio-plausibility 

 <cit.> (and tables therein) or coding efficiency <cit.>. 
The conversion of DNNs to spiking neural networks (SNN) after training with backpropagation <cit.> is a common technique to evade the difficulties of training with spikes. 
Furthermore, there are models including recurrent activity  

<cit.>, starting directly from realistic circuits 

<cit.>, or combining unsupervised and supervised training <cit.> as in this paper. 
We refer to <ref> for an extensive list of current biologically plausible models tested on MNIST (see <ref> for abbreviations).





§ RESULTS



We study networks that consist of an input (l_0), one hidden (l_1) and an output-layer (l_2) of (nonlinear) units, connected by weight matrices W_1 and W_2 (<ref>).

Training the hidden layer weights W_1 with standard supervised training involves (non-local) error backpropagation using summation over output units, the derivative of the units' nonlinearity (φ'(·)) and the transposed weight matrix W^T_2 (<ref>a).
In the biologically plausible network considered in this paper (<ref>b & c), the input-to-hidden weights W_1 are either fixed random, random Gabor filters or learned with an unsupervised method (Principal/ Independent Component Analysis or Sparse Coding). 
The unsupervised learning algorithms assume recurrent inhibitory weights V_1 between hidden units to implement competition, i.e. to make different hidden units learn different features. For more model details we refer to <ref> - <ref>. Code for all (rate & spiking) models discussed below is publicly available at https://github.com/EPFL-LCN/pub-illing2019-nnetworkshttps://github.com/EPFL-LCN/pub-illing2019-nnetworks.




 §.§ Benchmarking biologically plausible rate models and backpropagation




To see how far we can go with a single hidden layer, we systematically investigate rate models using different methods to initialize or learn the hidden layer weights W_1 (see <ref> and methods <ref>-<ref> for details). 
We use two different ways to set the weights W_1 of the hidden layer: either using fixed Random Projections (RP) or Random Gabor filters (RG), see <ref>b & blue curves in <ref>, or using one of the unsupervised methods Principal Component Analysis (PCA), Independent Component Analysis (ICA) or Sparse Coding (SC), see <ref>c & red curves in <ref>.  
All these methods can be implemented with local, biologically plausible learning rules <cit.>. We refer to the methods <ref> for further details. 
As a reference, we train networks with the same architecture with standard backpropagation (BP, see <ref>a).
As a step from BP towards increased biologically plausibility, we include Feedback Alignment (FA, <cit.>) with fixed random feedback weights for error backpropagation (see methods <ref> for further explanation). 

A Simple Perceptron (SP) without a hidden layer serves as a further reference, 
since it corresponds to direct classification of the input. 
We expect any biologically plausible learning algorithm to achieve results somewhere between SP (“lower”) and BP (“upper performance bound”)



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). 
Therefore the two-layer network as a whole is biologically plausible in terms of online learning and synaptic updates using only local variables. 
For computational efficiency, we first train the hidden layer and then the output layer, however, both layers could be trained simultaneously.


We compare the test errors on the MNIST digit recognition data set for varying numbers of hidden neurons n_h (<ref>).
The PCA (red dashed) and ICA (red dotted) curves in <ref> end at the vertical line n_h = d = 784 because the number of principal/independent components (PCs/ICs), 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. 
ICA does not seem to discover significantly more useful features than PCA, leading to similar classification performance.

SC (red solid line) extracts sparse representations that can be overcomplete (n_h > d), leading to 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.

As expected, the performance of RP (blue solid) for small numbers of hidden units (n_h < d) is worse than for feature extractors like PCA, ICA or SC.
Also for large hidden layers, performance improves only slowly with n_h, which is in line with theory <cit.> and findings in the extreme learning field <cit.>.
However, for large hidden layers sizes, RP outperforms SC.

As a reference, we also studied fixing the hidden layer weights to Gabor filters of random orientation, phase and size, located at the image center (RG, blue dashed, see <ref>). For hidden layers with more than 1000 neurons, SC is only marginally better than the network with fixed random Gabor filters.




For all tested methods and hidden layer sizes, performance is significantly worse than the one reached with BP (black solid in <ref>).
In line with <cit.>, we find that FA (black dashed) 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 (𝒪(1/n_h)) and RP (𝒪(1/n_h^2/d), where d is the input dimension <cit.>). In the log-log-plot in <ref> 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_/2 ≈ 10 using an effective dimensionality of MNIST (see methods <ref>).
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.

Taken together, these results show that the high dimensionality of the hidden layers is more important for reaching high performance than the global features extracted by PCA, ICA or SC. Tests on the object recognition task CIFAR10 lead to the same conclusion, indicating that this observation is not entirely task specific (see <ref> for further analysis on CIFAR10).





 §.§ Localized receptive fields boost performance




There are good reasons to reduce the connectivity from all-to-all to localized receptive fields (<ref>e & f): local connectivity patterns are observed in real neural circuits <cit.>, useful theoretically <cit.> and empirically <cit.>, and successfully used in convolutional neural 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.

To obtain localized receptive fields (called “l-” methods in the following) patches spanning p× p pixels in the input space are assigned to the hidden neurons. The centers of the patches are chosen at random positions in the input space, see <ref>e & f. For localized Random Projections (l-RP) and localized random Gabor filters (l-RG) the weights within the patches are randomly drawn from the respective distribution and then fixed. 
For the localized unsupervised learning methods (l-PCA, l-ICA & l-SC) the hidden layer is split into 500 independent populations. Neurons within each population compete with each other while different populations are independent, see <ref>f. This split implies a minimum number of n_h = 500 hidden neurons for these methods. For l-PCA and l-ICA a thresholding nonlinearity was added to the hidden layer to leverage the local structure (otherwise PCA/ICA act globally due to their linear nature, see methods <ref>).

We test l-RP for different patch sizes p and find an optimum around p≈ 10 (see <ref>a).
Note that p = 1 corresponds to resampling the data with random weights, and p = 28 recovers fully connected RP performance.
The other methods show similar optimal values around p = 10 (not shown).
The main finding here is the significant improvement in performance using localized receptive fields. 
All tested methods improve by a large margin when switching from full image to localized patches and some methods (l-RG and l-ICA) even reach BP performance for n_h = 5000 hidden neurons (see <ref>b). 
To achieve a fair comparison BP is also implemented with localized receptive fields (l-BP) which leads to a minor improvement compared to global BP. 
This makes local random projections or local unsupervised learning strong competitors to BP as biologically plausible algorithms in the regime of large, overcomplete hidden layers n_h > d - at least for MNIST classification.



To test whether localized receptive fields only work for the relatively simple MNIST data set (centered digits, uninformative margin pixels, no clutter, uniform features and perspective etc.) or generalizes to more difficult tasks, we apply it to the CIFAR10 data set <cit.>.
We first reproduce 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 <cit.>).
Again, classification performance increases for increasing hidden layer size n_h and localized receptive fields perform better than full connectivity for all methods.
Furthermore, as on MNIST, we can see similar performances for local feature learning methods (l-PCA, l-ICA & l-SC) and local random features (l-RP, l-RG) in the case of large, overcomplete hidden layers (see <ref>).

Also on CIFAR10, localized random filters and local feature learning reach the performance of biologically plausible models of deep learning <cit.> and come close to the performance of the reference algorithm l-BP. However, the difference remains statistically significant here.

Given that the state-of-the-art performance on CIFAR10 with deep convolutional neural networks is close to 98% (e.g. <cit.>), the limitations of our shallow local network and the well-known differences in difficulty between MNIST and CIFAR10 become apparent.

In summary, the main message of this section is that unsupervised methods, as well as random features, perform significantly better when applied locally. Equipped with local receptive fields our shallow network can outperform many current models of biologically plausible deep learning (see <ref>). On MNIST some models (l-RG & l-ICA) even reach backpropagation performance, while on CIFAR10 large differences to state-of-the-art deep convolutional networks remain.





 §.§ Spiking localized random projections


Real neural circuits communicate with short electrical pulses, called spikes, instead of real numbers such as rates. We thus extend our shallow network model to networks of leaky integrate-and-fire (LIF) neurons.  
The network architecture is the same as in <ref>b. 
To keep it simple we implement the two models with fixed random weights with LIF neurons: fixed localized Random Projections (l-RP) and fixed localized random Gabor filters (l-RG) with patches of size p× p - as in <ref>. 
The output layer weights W_2 are trained with a supervised spike timing dependent plasticity (STDP) rule.




The spiking dynamics follow the usual LIF equations (see methods <ref>) and the readout weights W_2 evolve according to a supervised delta rule via spike timing dependent plasticity (STDP) using post-synaptic spike-traces 
_i(t) and a post-synaptic target trace _i(t)

    τ_d _i(t)/dt   =    -_i(t) + ∑_f δ(t - t_i^f)
    Δ w_2,ij   =   α·(^_i(t) - ^_i(t))δ(t - t_j^f),

where α is the learning rate. 
Thus, for a specific readout weight w_2,ij, the post-synaptic trace is updated at every post-synaptic spike time t_i^f and the weight is updated at every pre-synaptic spike time t_j^f.
The target trace is constant while a pattern is presented and uses a standard one-hot coding for the supervisor signal in the output layer (l_2).

To illustrate the LIF and STDP dynamics, a toy example consisting of one pre- connected to one post-synaptic neuron is integrated for 650 ms. The pre- and post-synaptic membrane potentials show periodic spiking (<ref>a) which induces post-synaptic spike traces and corresponding weight changes (<ref>b), according to <ref>. 
For the MNIST task, <ref>c shows a raster plot for an exemplary training and testing protocol. 
During activity transients after a switch from one pattern to the next, learning is disabled until regular spiking is recovered. 
We experienced that without disabling learning during these transient phases the networks never reached a low test error. 
This is not surprising, since in this phase the network activities carry information both about the previously presented pattern and the current one, but the learning rule is designed for network activities in response to a single input
pattern. 
It is also known that LIF neurons differ from biological neurons in response to step currents (see <cit.> and references therein). 
During the testing period, learning is shut off permanently (see methods section <ref> for more details).
The LIF and STDP dynamics can be mapped to a rate model (see e.g. <cit.> and <ref> for details). However all following results are obtained with the fully spiking LIF/STDP model.



When directly trained with the STDP rule of <ref>, the spiking LIF models closely approach the performance of their rate counterparts.


<ref> compares the performances of the rate and spiking LIF l-RP & l-RG models with the reference algorithm l-BP (for same hidden layer size n_h and patch size p, see <ref>).
The remaining gap (< 0.3%) between rate model and spiking LIF model presumably stems from noise introduced by the spiking approximation of rates and the activity transients mentioned above. 
Both, the rate and spiking LIF model of l-RP/l-RG achieve accuracies close to the backpropagation reference algorithm l-BP and fall in the range of performance of prominent, biologically plausible models, i.e. 98-99% test accuracy (see <ref> and <ref>).
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.



§ DISCUSSION



In contrast to biologically plausible deep learning algorithms that are derived from approximations of the backpropagation algorithm <cit.>, we focus here on shallow networks with only one hidden layer. 
The weights from the input to the hidden layer are either learned by unsupervised algorithms with local learning rules; or they are fixed. 
If fixed, they are drawn randomly or represent random Gabor filters.
The readout layer is trained with a supervised, local learning rule.

When applied globally, randomly initialized fixed weights/ Gabor filters (RP/RG) of large hidden layers lead to better classification performance than training them with unsupervised methods like Principal/Independent Component Analysis (PCA/ICA) or Sparse Coding (SC). Such observations also occur in different contexts, e.g. <cit.> showed that (sparse) random projections, combined with dimensionality expansion outperform known  algorithms for locality-sensitive hashing.

It may be interesting to search for alternative unsupervised, local learning rules with an inductive bias that is better adapted to image processing tasks than the one of SC.

Replacing all-to-all connectivity with localized input filters is such an inductive bias that already proved useful in supervised models <cit.> but turns out to be particularly powerful in conjunction with unsupervised learning (l-PCA, l-ICA & l-SC). Interestingly, non of the local unsupervised methods could significantly outperform localized random Gabor filters (l-RG). 

Furthermore, we find that the performance scaling with the number of hidden units n_h is orders of magnitudes better than the lower bound suggested by universal function approximation theory <cit.>.


To move closer to realistic neural circuits we implement our shallow, biologically plausible network with spiking neurons and spike timing dependent plasticity to train the readout layer. 
Spiking localized random projections (l-RP) and localized Gabor filters (l-RG) reach >98% test accuracy on MNIST which lies within the range of current benchmarks for biologically plausible models for deep learning (see <ref> and <ref>).
Our network model is particularly simple, i.e. it has only one trainable layer and does not depend on sophisticated architectural or algorithmic features typically necessary to approximate backpropagation <cit.>.
Instead it only relies on the properties of high-dimensional localized random projections.

Since we want to keep our models as simple as possible, we use online stochastic gradient descent (SGD, no mini-batches) with a constant learning rate. 
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 is to demonstrate that even a simple model with constant learning rate achieves results that are comparable with more elaborate approaches that use e.g. convolutional layers with weight sharing <cit.>, backpropagation approximations <cit.>, multiple hidden layers <cit.>, dendritic neurons <cit.>, recurrence <cit.> or conversion from rate to spikes <cit.>.

Above 98% accuracy we also 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.
Localized random filters or local unsupervised feature learning perform remarkably well compared to fully-connected backpropagation in shallow networks, even on more challenging data sets such as CIFAR10. This makes our model an important benchmark for future, biologically plausible models but also clearly highlights the limitations of our shallow two-layer model.  
A long time ago state-of-the-art deep learning has moved from MNIST to harder datasets, such as CIFAR10 or ImageNet <cit.>.  
Yet MNIST seems to be the current reference task for most biologically plausible deep learning models (see <ref> and <ref>).
We suggest that novel, progressive approaches to biologically plausible deep learning should significantly outperform the results presented here. 
Furthermore, they should be tested on tasks other than MNIST, where real deep learning capabilities become necessary.



§ ACKNOWLEDGMENTS


This research was supported by the Swiss National Science Foundation (no. 200020_165538 and 200020_184615) and by the European Union Horizon 2020 Framework Program under grant agreement no. 785907 (HumanBrain Project, SGA2).





§ GENERAL RATE MODEL DETAILS



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 = 10 = number of classes). The layers are connected via weight matrices W_1 ∈ℝ^n_1 × n_0 and W_2 ∈ℝ^n_2 × n_1 and each neuron receives bias from the bias vectors b_1 ∈ℝ^n_1 and b_2 ∈ℝ^n_2 respectively (see <ref>). 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+1a_l + b_l+1
    a_l+1   =   φ_l+1 (u_l+1).


The Simple Perceptron (SP) only consists of one layer (l_2, W_2 ∈ℝ^n_2 × n_0, b_2 ∈ℝ^n_2). 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 <cit.>, but direct, plastic, inhibitory synapses V_1 ∈ℝ^n_1 × n_1 are assumed between neurons in l_1. Classification error variances in <ref> & <ref> 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_ of the MNIST data set can be obtained, e.g. via eigen-spectrum analysis, keeping 90% of the variance. We obtain values around d_≈ 20. The measure proposed by <cit.> gives the same value d_≈ 20. 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 instead of the full data set leads to a comparable MNIST performance (1.7% vs 1.5% test error respectively).
Together, these findings suggest that the MNIST dataset lies mostly in a low-dimensional linear subspace with d_≈ 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. Simulations were implemented and performed in the https://julialang.org/-language. The code for the implementation of our rate network models is publicly available at https://github.com/EPFL-LCN/pub-illing2019-nnetworkshttps://github.com/EPFL-LCN/pub-illing2019-nnetworks.






§ UNSUPERVISED METHODS (PCA, ICA & SC)

In this paper we do not implement PCA/ICA learning explicitly as a neural learning algorithm but by a standard PCA/ICA algorithm (https://github.com/JuliaStats/MultivariateStats.jlMultivariateStats.jl) since biologically plausible online algorithms for both methods are well known <cit.>. For d-dimensional data such algorithms output the values of the n ≤ d first principal/ independent components as well as the corresponding subspace projection matrix P∈ℝ^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/independent principal components. In other words each neuron in the hidden layer l_1 extracts another principal/independent component of the data. ICA was performed with the usual pre-whitening of the data.

Since PCA/ICA is a linear model, biases b_1 were set to 0 and φ_1(u) = u. With this, we can write the (trained) feed-forward pass of the first layer of our PCA/ICA model as follows:

    a_1 = u_1 = W_1·a_0W_1 = P


Since the maximum number of principal/independent components that can be extracted is the dimensionality of the data, n_ = d, the number of neurons in the hidden layer n_1 is limited by d. This makes PCA/ICA unusable for overcomplete hidden representations as investigated for SC and RP. 

In the localized version of PCA/ICA we assume the hidden layer to consist of independent populations, each extracting PCs/ICs of its respective localized receptive field (see <ref>). The hidden layer was divided into 500 of those populations, resulting in a minimum number of n_h = 500 hidden neurons (1 PC/IC per population) for these methods (and up to 10 PCs/ICs per population for n_h = 5000). The classifier was then trained on the combined activations of all populations of the hidden layer. Because PCA/ICA are linear methods the localized PCA/ICA version would not extract significantly different features unless we introduce a nonlinearity in the hidden units. This was done by simply thresholding the hidden activations (ReLU with threshold 0). No further optimization in terms of nonlinearity- and threshold-tuning was performed.



Sparse coding (SC) aims at finding a feature dictionary W∈ℝ^h × d (for d-dimensional data) that leads to an optimal representation a_1 ∈ℝ^h which is sparse, i.e. has as few non-zero elements as possible. The corresponding optimization problem reads:


    W^opt, a_1^opt   =   ℒ(W,a_1)
    ℒ(W,a_1)    =   1/2‖a_0 - W^⊤a_1 ‖^2_2 + λ‖a_1‖_1.

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 <cit.>. For a rigorous motivation (and derivation) that such a network architecture can indeed implement sparse coding we refer to <cit.>. 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:




  * 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 <cit.>:

    τ_u du_1(t)/dt   =    - u_1(t) + (W_1a_0(t) - V_1 a_1(t))
    a_1(t)    =   φ(u_1(t))

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



  * Optimizing the weights:

Now the activities a_1 are kept fixed and we update the weights following the gradient of the loss function. The weight update rules are Hebbian-type local learning rules <cit.>:

    Δ W_1,ji   =   α_w ·  a_0,i· a_1,j
    Δ V_1,jk   =   α_v ·  a_1,k·(a_1,j - ⟨ a_1,j⟩)

⟨·⟩ is a moving average (low-pass filter) over several past hidden representations (after convergence of the recurrent dynamics) with some time constant τ_, e.g. τ_ = 100 patterns. 
At the beginning of the simulation (or after a new pattern presentation) τ_ is increased starting from 0 to τ_ during the first τ_. 
The values of the rows 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 variables u_1(t) and a_1(t) are reset (to avoid transients) before an input is presented. 
Then for every of the N iterations, (<ref>) is iterated for N_ steps and the weights are updated according to (<ref>).


Similar to localized PCA/ICA, the localized version of SC uses independent populations in the hidden layer (see <ref>). The SC algorithm above was applied to each population and its respective receptive field independently. The classifier was then trained on the combined activations of all populations of the hidden layer.




§ FIXED RANDOM FILTERS (RP & RG)


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

    W_1    ∼   𝒩(0,σ^2)   σ^2 = 1/100 n_0
    b_1    ∼   𝒰([0,0.1])

for RP (keeping weights fixed), SC, SP and also BP & RF (both layers with W_2, b_2 and n_1 respectively).

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 ∈ℕ) 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:

    W_1,   ∼   𝒩(0,σ_^2)   σ_^2 = c/100 p
    b_1    ∼   𝒰([0,0.1])

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

The (localized) random Gabor filters in RG have the same receptive field structure as in l-RP (see <ref>) but instead of choosing the weights within the receptive field as random values, they are choosen according to Gabor filters W_1 ∝ g(x,y). Here, x and y denote the pixel coordinates within the localized receptive field relative to the patch center. The Gabor filters have the following functional form:


    g(x,y; λ,Θ,ψ,σ,γ)    =   
    exp( -x'^2+γ^2 y'^2/2σ^2)    ·   cos(2πx'/λ+ψ)
    [ x'; y' ]   =   [  cosΘ  sinΘ; -sinΘ  cosΘ ]·[ x; y ]


To obtain diverse, random receptive fields we draw the parameters λ,Θ,ψ,σ,γ of the Gabor functions from uniform distributions over some intervals. The bounds of the sampling interval are optimized using Bayesian optimization (https://github.com/jbrea/BayesianOptimization.jlBayesianOptimization.jl) with respect to classification accuracy on the training set.



§ CLASSIFIER & SUPERVISED REFERENCE ALGORITHMS (BP, FA & SP)

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 biologically plausible. 
For having a reference for our biologically plausible models (<ref>b & c), we compare it to networks with the same architecture (number of layers, neurons, connectivity) but trained in a fully supervised way with standard backpropagation (<ref>a). The forward pass of the model reads:

    u_l+1   =   W_l+1a_l + b_l+1
    a_l+1   =   φ_l+1 (u_l+1)


Given the one-hot encoded target activations tgt, the error ẽ_L is 

    ẽ_L    =   tgt - a_L
 
 when minimizing mean squared error (MSE)

    ℒ_   =   1/2‖tgt - a_L ‖_2^2

 or 

    p   =   (a_L)
    ẽ_L    =   tgt - p
 
for the softmax/cross-entropy loss (CE),

    ℒ_   =    - ∑_i=1^n_L_i ·log(p_i).


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∝∂ℒ/∂ W_l,ij. The full BP algorithm for deep networks reads <cit.>:

    e_L    =   φ '_L(u_L) ⊙ẽ_L
    e_l-1   =   φ '_l-1(u_l) ⊙W_l^⊤e_l
    ΔW_l    =   α·e_l ⊗a_l-1
    Δb_l    =   α·e_l








where ⊙ stands for element-wise multiplication, ⊗ is the outer (dyadic) product, φ '_l(· ) is the derivative of the nonlinearity and α is the learning rate. FA <cit.> uses a fixed random matrix R_l instead of the transpose of the weight matrix W_l^⊤ for the error backpropagation step in (<ref>). 


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 (<ref>) was applied to the weights within the receptive fields. The exact parameter values can be found in <ref>.







§ SPIKING IMPLEMENTATION OF RP & RG 

The spiking simulations were performed with a custom-made event-based leaky integrate-and-fire (LIF) integrator written in the https://julialang.org/-language. 
Code is available at https://github.com/EPFL-LCN/pub-illing2019-nnetworkshttps://github.com/EPFL-LCN/pub-illing2019-nnetworks.
 For large network sizes, the exact, event-based integration can be inefficient due to a large frequency of events. 
We thus also added an Euler-forward integration mode to the framework. 
For sufficiently small time discretization (e.g. Δ t ≤ 5· 10^-2 ms for the parameters given in <ref>) the error of Euler-forward integration does not have negative consequences on the learning outcome. 
The dynamics of the LIF network is given by:

    τ_m du_i(t)/dt   =    -u_i(t) + RI_i(t)
     I_i(t)    =    I^ff_i(t) + I^ext_i(t) 
       =   ∑_j,f w_ijϵ(t - t_j^f)+  I^ext_i(t)

and the spiking condition:


u_i(t) ≥ϑ_i: u_i → u_


, 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 is the post-synaptic potential evoked by a pre-synaptic spike arrival, ϑ_i is the spiking threshold and u_ the reset potential after a spike.

The input is split into a feed-forward (I^ff(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 
_i(t) and a post-synaptic target-trace _i(t). 

    τ_d _i(t)/dt   =    -_i(t) + ∑_f δ(t - t_i^f)
    
    Δ w_ij   =    g(^_i(t),^_i(t))δ(t - t_j^f)

with the plasticity function 

    g(^_i(t),_i(t)) = α·(^_i(t) - ^_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 _, the targets tgt(t) of the output layer are the one-hot-coded classes, scaled by the target amplitude _. 
Additionally, every neuron receives a static bias input I^_≈ϑ to avoid silent units in the hidden layer. 
Every pattern is presented as fixed input for a time T_ and the LIF dynamics as well as the learning evolves according to (<ref>) and (<ref>) respectively. Learning is disabled after pattern switches for a duration of T_ = 4τ_m since the noise introduced by these transient phases was found to deteriorate learning progress.
With the parameters we used for the simulations (see <ref>), firing rates of single neurons in the whole network stayed below 1 kHz which was considered as a biologically plausible regime. 
For the toy example in <ref>a& b we used static input and target with the parameters amp_ = 40, amp_ = 5 (i.e. target trace = 0.005), ϑ_ = 20,  σ_ϑ= 0, τ_m = 50, α = 1.2· 10^-5. For the raster plot in <ref>c we used amp_ = 300, amp_ = 300, ϑ_ = 20,  σ_ϑ= 0, τ_m = 50, α = 1.2· 10^-5, T_ = 50 ms, T_ = 100 ms. The LIF dynamics can be mapped to a rate model described by the following equations:

    u_l    =   W_l u_l-1 + RI^ext
    a_l    =   φ_(u_l)
    Δ w_ij   =   g̃(a^_j,a^_i,^_i)

with the (element-wise) LIF-activation function φ_(·) and the modified plasticity function g̃(·):

    φ_(u_k)    =   [Δ_ - τ_m ln(1 - ϑ_k/u_k)]^-1
    g̃(a^_j,a^_i,^_i)    =   α̃· a^_j ·(^_i - a^_i)

The latter can be obtained by integrating the STDP rule of <ref> and taking the expectation over spike times. 
Most of the parameters of the spiking- and the LIF rate models can be mapped to each other directly (see <ref>). 
The learningrate α must be adapted since the LIF weight change depends on the presentation time of a pattern T_. In the limit of long pattern presentation times (T_≫τ_m , τ_), the theoretical transition from the learning rate of the LIF rate model (α̃) to the one of the spiking LIF model (α) is 


α = 1000/T_· 1000·α̃,



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. 
The reasons for the remaining difference in performance presumably lie in transients and single-spike effects that cannot be captured by the rate model. 
Furthermore the new target was presented immediately after a pattern switch even though the activity obviously needs at least a couple time constants (τ_ or τ_m) to propagate through the network. 
Removing this asynchrony between input and target should further shrink the discrepancy between rate and spiking models.













unsrtnat

  










gobble%

auxoutThis is a test line.

