Title: AN ADAPTIVE HOMEOSTATIC ALGORITHM FOR THE UNSUPERVISED LEARNING OF VISUAL FEATURES

Abstract: The formation of structure in the brain, that is, of the connections between cells within neural populations, is by large an unsupervised learning process: the emergence of this architecture is mostly self-organized. In the primary visual cortex of mammals, for example, one may observe during development the formation of cells selective to localized, oriented features. This leads to the development of a rough representation of contours of the retinal image in area V1. We modeled these mechanisms using sparse Hebbian learning algorithms. These algorithms alternate a coding step to encode the information with a learning step to find the proper encoder. A major difficulty faced by these algorithms is to deduce a good representation while knowing immature encoders, and to learn good encoders with a non-optimal representation. To address this problem, we propose to introduce a new regulation process between learning and coding, called homeostasis. Our homeostasis is compatible with a neuro-mimetic architecture and allows for the fast emergence of localized filters sensitive to orientation. The key to this algorithm lies in a simple adaptation mechanism based on non-linear functions that reconciles the antagonistic processes that occur at the coding and learning time scales. We tested this unsupervised algorithm with this homeostasis rule for a range of existing unsupervised learning algorithms coupled with different neural coding algorithms. In addition, we propose a simplification of this optimal homeostasis rule by implementing a simple heuristic on the probability of activation of neurons. Compared to the optimal homeostasis rule, we show that this heuristic allows to implement a more rapid unsupervised learning algorithm while keeping a large part of its effectiveness. These results demonstrate the potential application of such a strategy in machine learning and we illustrate this with one result in a convolutional neural network.
This supplementary information presents : first, the code to generate the figures from the paper, second, some control experiments that were mentionned in the paper, finally, some perspectives for future work inspired by the algorithms presented in the paper.

1. INTRODUCTION:

The neural architecture is a complex dynamic system that operates at different time scales. In particular, one of its properties is to succeed in representing quickly information (the coding step) while optimizing in the long term its encoding (the learning step). In the case of the mammalian primary visual cortex (V1) for instance, this rapid coding operation, of the order of 50 milliseconds in humans, is the key to the results of Hubel & Wiesel (1968) , who showed that some cells of V1 have relatively localized receptive fields which are predominantly selective at different orientations. As such, one can consider the rapid coding of the retinal image as a process of transforming the raw visual information into a rough "sketch" that represents the outlines of objects in the image by using elementary edge-like features. This internal representation and the visual information share the same property of being sparse: for most natural images, only a relatively small number of features are necessary to describe the input. Thus, the coding step consists in choosing the right encoder that selects as few features (called atoms) as possible among a collection of them (called the dictionary). Amazingly, Olshausen & Field (1996) have shown that when enforcing a sparse prior on the encoding step, such edge-like filters are emerging using a simple Hebbian unsupervised learning strategy.

Additionally, recent advances in machine learning, and especially on unsupervised learning have shed new light on the functioning of the underlying biological neural processes. By definition, unsupervised learning aims at learning the best dictionary to represent the input image autonomously, that is, without using other external knowledge such as in supervised or reinforcement learning.

Under review as a conference paper at ICLR 2019 Algorithms that combines such learning as the input to classical, supervised deep-learning show great success in tasks like image denoising (Vincent et al., 2008) or classification (Sulam et al., 2017) . A variant consists in forcing the generated representation to be sparsely encoded (Makhzani & Frey, 2013) , whether by adding a penalty term to the optimized cost function or by encoding each intermediate representation by a pursuit algorithm (Papyan et al., 2016) . Interestingly, (Papyan et al., 2016) proposes a model of Convolutional Sparse Coding (CSC) tightly connected with Convolutional Neural Network (CNN), so much that the forward pass of the CNN is equivalent to a CSC with a thresholding pursuit algorithm. These unsupervised algorithms are equivalent to a gradient descent optimization over an informational-type coding cost (Kingma & Welling, 2013) . This cost makes it then possible to quantitatively evaluate the joint exploration of new learning or coding strategies. As such, this remark shows us that unsupervised learning consists of two antagonistic mechanisms, a long time scale that corresponds to the learning and exploration of new components and a faster scale that corresponds to coding.

In particular, an aspect often ignored in this type of learning is the set of homeostasis mechanisms that control the average activity of neurons within a population. Indeed, there is an intrinsic complexity in unsupervised dictionary learning algorithms: how to adapt the regularization parameter of each atom to make sure no atoms are wasted because of improper regularization settings? In the original algorithms of sparse unsupervised learning (Olshausen & Field, 1997) , homeostasis is implemented as a heuristic that prevents the algorithm from diverging. In most unsupervised learning algorithms it takes the form of a normalization, that is, an equalization of the energy of each atom in the dictionary (Mairal et al., 2014) . However, the neural mechanisms of homeostasis are at work in many components of the neural code and are essential to the overall transduction of neural information. For example, the sub-networks of glutamate and GABA-type neurons may regulate the overall activity of neural populations (Marder & Goaillard, 2006) . In particular, such mechanisms could be tuned to balance the contribution of the excitatory populations with respect to that in inhibitory populations. As a consequence, this creates a so-called balanced network which may explain many facets of the properties of the primary visual cortex (Hansel & van Vreeswijk, 2012) . At the modelling level, these mechanisms are often implemented in the form of normalization rules (Schwartz & Simoncelli, 2001) which are considered as the basis of a normative theory to explain the function of the primary visual cortex (Carandini & Heeger, 2012) . However, when extending such model using unsupervised learning, most modelling effort is rather intended to show that the cells' selectivity that emerges have the same characteristics than those observed in neuro-physiology (Ringach, 2002; Rehn & Sommer, 2007; Loxley, 2017) . Other algorithms use non-linearities that implicitly implement homeostatic rules in neuro-mimetic algorithms (Brito & Gerstner, 2016) . These non-linearities are mainly used in the output of successive layers of deep learning networks that are nowadays widely used for image classification or artificial intelligence. However most of these non-linear normalization rules are based on heuristics mimicking neural mechanisms but are not justified as part of the global problem underlying unsupervised learning. Framing this problem in a probabilistic framework allows to consider in addition to coding and learning the intermediate time scale of homeostasis and allows us also to associate it to an adaptation mechanisms (Rao & Ballard, 1999) . Our main argument is that by optimizing unsupervised learning at different time scales, we allow for the implementation of fast algorithms compatible with the performance of biological networks and in comparison with classical (Olshausen & Field, 1997) or Deep Learning approaches.

In this paper, we will first define a simple algorithm for controlling the selection of coefficients in sparse coding algorithms based on a set of non-linear functions similar to a generic neural gain normalization mechanisms. Such functions will be used to implement an homeostasis mechanism based on histogram equalization by progressively adapting these non-linear functions. In particular, this algorithm will extend an already existing algorithm of unsupervised sparse learning (Perrinet, 2010) to a more general setting. In particular, we will show quantitative results of this optimal algorithm by applying it to different pairs of coding and learning algorithms. Second, we will propose a simplification of this homeostasis algorithm based on the activation probability of each neuron and show that it yields similar quantitative results as the full homeostasis algorithm and that it converges more rapidly than classical methods (Olshausen & Field, 1997; Sandin & Martin-del Campo, 2017) . In particular, we focused in our architecture to be able to quantitatively cross-validate for every single hyper-parameters and all these scripts are available at https://github.com/XXX/ZZZ. Finally, we will conclude by showing an application of such an adaptive algorithm to CNNs and discuss on its development in real-world architectures. The compared algorithms are : None (using a simple simple normalization of the atoms), OLS (the method of (Olshausen & Field, 1997) ), HEH (using the optimal homeostasis described in this paper). (A) For each algorithm, we show 18 atoms from the 441 filters of the same size as the image patches (M = 22 × 22 = 484, circularly masked) and presented in a matrix (separated by a white border). The upper and lower row respectively show the least and most probably selected atoms. This highlights the fact that without proper homeostasis, dictionary learning leads to inhomogeneous representations. (B) Evolution of cost F (in bits, see Eq. 4) as a function of the number of iterations and cross-validated over 10 runs. While OLS provides a similar convergence than None, the HEH method provides a better final convergence.

2. UNSUPERVISED LEARNING AND THE OPTIMAL REPRESENTATION OF:



IMAGES:

Visual items composing natural images are often sparse, such that knowing a model for the generation of images, the brain may use this property to reconstruct images using only a few of these items. In the context of the representation of natural images1 y = (y k ) K k=1 ∈ R K×M represented in a matrix as a set of K vector samples (herein, we will use a batch size of K = 256) as images raveled along M pixels (each y k,j ∈ R are the corresponding luminance values), let us assume the generic Generative Linear Model, such that for any sample k the image was generated as y k = Φ T a k + , where by definition, the coefficients are denoted by a k = (a k,i ) N i=1 ∈ R N and the dictionary by Φ ∈ R N ×M . Finally, ∈ R M is a Gaussian iid noise which is Normal without loss of generality by scaling the norm of the dictionary's rows. Knowing this model, unsupervised learning aims at finding the least surprising causes (the parameters âk and Φ) for the data y k . In particular, the cost may be formalized in a probabilistic terms as (Olshausen & Field, 1997) :

EQUATION

Such hypothesis allows to retrieve the cost that is optimized in most of existing models of unsupervised learning. Explicitly, the representation is optimized by minimizing a cost defined on prior assumptions on representation's sparseness, that is on log p Φ (a k ). For instance, learning is accomplished in SPARSENET (Olshausen & Field, 1997 ) by defining a sparse prior probability distribution function for each coefficients in the factorial form log p Φ (a k ) ∼ -β i log(1 +

a 2 i σ 2 )

where β corresponds to the steepness of the prior and σ to its scaling (see Figure 13 .2 from (Olshausen, 2002) ). Then, knowing this sparse solution, learning is defined as slowly changing the dictionary using Hebbian learning.

Under review as a conference paper at ICLR 2019 Indeed, to compute the partial derivate of F with respect to Φ, we have ∀i:

EQUATION

This allows to define unsupervised learning as the gradient descent using this equation. Similarly to Eq. 17 in (Olshausen & Field, 1997) or to Eq. 2 in (Smith & Lewicki, 2006) , the relation is a linear "Hebbian" rule (Hebb, 1949) since it enhances the weight of neurons proportionally to the activity (coefficients) between pre-and post-synaptic neurons. Note that there is no learning for non-activated coefficients and also that we used a (classical) scheduling of the learning rate and a proper initialization of the weights (see Annex 2.5 & 2.6). The novelty of this formulation compared to other linear Hebbian learning rule such as (Oja, 1982) is to take advantage of the sparse (non-linear) representation, hence the name Sparse Hebbian Learning (SHL). In general, the parameterization of the prior in Eq. 1 has major impacts on results of the sparse coding and thus on the emergence of edge-like receptive fields and requires proper tuning. For instance, a L2-norm penalty term (that is, a Gaussian prior on the coefficients) corresponds to Tikhonov regularization (Tikhonov, 1977) and a L1-norm term (that is, an exponential prior for the coefficients) corresponds to the convex cost which is optimized by least-angle regression (LARS) (Efron et al., 2004) or FISTA (Beck & Teboulle, 2009) .

2.1. ALGORITHM: SPARSE CODING WITH A CONTROL MECHANISM FOR THE SELECTION OF:



ATOMS:

Concerning the choice of a proper prior distribution, the spiking nature of neural information demonstrates that the transition from an inactive to an active state is far more significant at the coding time scale than smooth changes of the firing rate. This is for instance perfectly illustrated by the binary nature of the neural code in the auditory cortex of rats (DeWeese et al., 2003) . Binary codes also emerge as optimal neural codes for rapid signal transmission (Bethge et al., 2003) . This is also relevant for neuromorphic systems which transmit discrete, asynchronous events (such as a network packet). With a binary event-based code, the cost is only incremented when a new neuron gets active, regardless to its (analog) value. Stating that an active neuron carries a bounded amount of information of λ bits, an upper bound for the representation cost of neural activity on the receiver end is proportional to the count of active neurons, that is, to the 0 pseudo-norm ||a k || 0 = |{i, a k,i = 0}|:

EQUATION

This cost is similar with information criteria such as the Akaike Information Criteria (Akaike, 1974) or distortion rate (Mallat, 1998, p. 571) . For λ = log 2 N , it gives the total information (in bits) to code for the residual (using entropic coding) and the list of spikes' addresses. In general, the high inter-connectivity of neurons (on average approximately 10000 synapses per neurons) justifies such an informational perspective with respect to the analog quantization of information in the point-topoint transfer of information between neurons. However, Eq. 3 defines a harder cost to optimize (in comparison to convex formulations in Equation 1 for instance) since the hard 0 pseudo-norm sparseness leads to a non-convex optimization problem which is NP-complete with respect to the dimension M of the dictionary (Mallat, 1998, p. 418) .

Still, there are many solutions to this optimization problem and here, we will use a generalized version of the Matching Pursuit (MP) algorithm (Mallat, 1998, p. 422) . A crucial aspect of this algorithm is the arg max function as it produces at each step a competition among N neurons (that is, λ bits).

For this reason, we will introduce a mechanism to tune this competition. For any signal y k drawn from the database, we get the coefficients a k = S(y k ; Ψ = {Φ, z, N 0 }) thanks to Algorithm 1. The parameter N 0 > 0 controls the amount of sparsity that we impose to the coding. The novelty of this generalization of MP lies in the scalar functions z = {z i } i=1...N which control the competition for the best match across atoms. While an identical symmetric function is chosen in the original MP algorithm (that is, ∀i, z i (a k ) = |a k |), we will define these at a first attempt as the rescaled non-linear rectified linear unit (ReLU) with gain

γ i : ∀i, z i (a k,i ) = γ i * a k,i * δ(a k,i > 0)

where δ is Kronecker's indicator function.

We found as in (Rehn & Sommer, 2007) that by using an algorithm like Matching Pursuit (that is using the symmetric function or setting ∀i, γ i = 1 as in (Mairal et al., 2014) for instance), the Sparse Hebbian Learning algorithm could provide results similar to SPARSENET. An advantage is

Under review as a conference paper at ICLR 2019 Algorithm 1 Generalized Matching Pursuit: a k = S(y k ; Ψ = {Φ, z, N 0 }) 1: set the sparse vector a k to zero, 2: initialize āki = y k , Φ i for all i 3: while ||a k || 0 < N 0 do:

4: select the best match: i * = arg max i [z i ( āki )] 5:

update the sparse coefficient: Histogram Equalization Homeostasis and its role in unsupervised learning: (A) Nonlinear homeostatic functions z i , ∀i learned using Hebbian learning. These functions were computed for different homeostatic strategies (None, OLS or HEH) but only used in HEH. Note that for our choice of N 0 = 13, all cumulative functions start around 1 -N 0 /N ≈ .970. At convergence of HEH, the probability of choosing any filter is equiprobable, while the distribution of coefficients is more variable for None and OLS. As a consequence, the distortion between the distributions of sparse coefficients is minimal for HEH, a property which is essential for the optimal representation of signals in distributed networks such as the brain. (B) Effect of learning rate η (eta) and homeostatic learning rate η h (eta homeo) on the final cost as computed for the same learning algorithms but with different homeostatic strategies (None, OLS or HEH). Parameters were explored around a default value, on a logarithmic scale and over 4 octaves. This shows that HEH is robust across a wide range of parameters.

a k,i * = a k,i * + āk,i * , 6: update residual coefficients: ∀i, ā k,i ← ā k,i -a k,i * Φ i * , Φ i .

the non-parametric assumption on the prior based on this more generic 0 pseudo-norm sparseness. However, we observed that this class of algorithms could lead to solutions corresponding to a local minimum of the full objective function: Some solutions seem as efficient as others for representing the signal but do not represent edge-like features homogeneously (Figure 1 -A, None). Moreover, using other sparse coding algorithms which are implemented in the sklearn library, we compared the convergence of the learning with different sparse coding algorithms. In particular, we compared the learning as implemented with matching pursuit to that with orthogonal matching pursuit (OMP) (Pati et al., 1993) , LARS or FISTA (see Annex 2.4). For all these sparse coding algorithms, during the early learning step, some cells may learn "faster" than others. In particular, these cells have more peaked distributions of their activity and tend to be selected more often. There is the need for a homeostasis mechanism that will ensure convergence of learning. The goal of this work is to study the specific role of homeostasis in learning sparse representations and to propose a homeostasis mechanism based on the functions z i which optimizes the learning of an efficient representation.

2.2. ALGORITHM: HISTOGRAM EQUALIZATION HOMEOSTASIS:

Knowing a dictionary and a sparse coding algorithm, we may transform any data sample y k into a set of sparse coefficients using the above algorithm: a k = S(y k ; Ψ = {Φ, z, N 0 }) (see Algorithm 1). In particular, at any step during learning, dictionaries may not have been homogeneously learned and

Under review as a conference paper at ICLR 2019 may exhibit different distributions. However, this would not be taken into account in the original cost (see Eq. 3) as we assumed as in Olshausen & Field (1997) that p Φ ( âk ) is factorized, that is, that the components of the sparse vector are independent. As a consequence, we may use a deviation to this hypothesis as an additional component to the cost:

EQUATION

Where we used the mutual information MI as a proxy to measure the dependence between the components of the sparse vector. Indeed, as information is coded in the address of neurons, information transfer as computed through Shannon entropy, is optimized when the activity within the neural population is uniformly balanced, that is when each neuron is a priori selected with the same probability. In particular, a necessary (yet not sufficient) condition for minimizing this cost is that the prior probability of selecting coefficients are identical ∀(i, j), q Ψ (a k,i ) = q Ψ (a k,j ) to ensure the optimality of the choice of the 0 pseudo-norm and compare it to the representation in the primary visual cortex. As we have seen, we may use different transformation functions z to influence the choice of coefficients such that we may use these functions to optimize the objective cost defined by Eq. 4.

To achieve this uniformity, we may define an homeostatic gain control mechanism based on histogram equalization, that is, by transforming coefficients in terms of quantiles by setting ∀i, z i (•) = P (• > a i ). Such a transform is similar to the inverse transform sampling which is used to optimize representation in auto-encoders (Doersch, 2016) and can be considered as a nonparametric extension of the "re-normalization trick" used in variational auto-encoders (Kingma & Welling, 2013) . Moreover, it has been found that such an adaptation mechanism is observed in the response of the retina to various contrast distributions (Laughlin, 1981) . However, an important point to note is that this joint optimization problem between coding and homeostasis is circular as we can not access the true posterior p Φ (a): Indeed, the coefficients depend on non-linear coefficients through a k = S(y k ; Ψ = {Φ, z i , N 0 }), while the non-linear functions depend on the (cumulative) distribution of the coefficients. We will make the assumption that such a problem can be solved iteratively by slowly learning the non-linear functions. Starting with an initial set of non-linear functions as in None, we will derive an approximation for the sparse coefficients. Then, the function z i for each coefficient of the sparse vector is calculated using an iterative moving average scheme (parameterized by time constant 1/η h ) to smooth its evolution during learning. At the coding level, this non-linear function is incorporated in the matching step of the matching pursuit algorithm, to modulate the choice of the most probable as that corresponding to the maximal quantile: i * = arg max i z i (a i ) (see Algorithm 1). We will coin this variant as Histogram Equalization Homeostasis (HEH). The rest of this Sparse Hebbian Learning algorithm is left unchanged. As we adapt the dictionaries progressively during Sparse Hebbian Learning, we may incorporate this HEH homeostasis during learning by choosing an appropriate learning rate η h . To recapitulate the different choices we made from the learning to the coding and the homeostasis, the unsupervised learning can be summarized using the following steps (see Algorithm 2).

Algorithm 2 Homeostatic Unsupervised Learning of Kernels: Φ = H(y; η, η h , N 0 ) 1: Initialize the point non-linear gain functions z i to similar cumulative distribution functions, 2: initialize atoms Φ i to random points on the K-unit sphere, 3: for T epochs do:

4:

draw a new batch y from the database of natural images, 5:

for each data point y k do: modify dictionary: ∀i,

Φ i ← Φ i + η • a i • (y k -Φa), 8: normalize dictionary: ∀i, Φ i ← Φ i /||Φ i ||, 9: update homeostasis functions: ∀i, z i (•) ← (1 -η h ) • z i (•) + η h • δ(a i ≤ •).

We compared qualitatively the set Φ of receptive filters generated with different homeostasis algorithms (see Fig. 1-A ). A more quantitative study of the coding is shown by comparing the decrease of the cost as a function of the iteration step (see Fig. 1-B ). This demonstrate that forcing the learning activity to be uniformly spread among all receptive fields results in a faster convergence of the representation error as represented by the decrease of the cost F . 

2.3. RESULTS: FAST UNSUPERVISED LEARNING USING HOMEOSTASIS:

We have shown above that we can find an exact solution to the problem of homeostasis during Sparse Hebbian Learning. However, this solution has several drawbacks. First, it is computationally intensive on a conventional computer as it necessitates to store each z i function to store the cumulative distribution of each coefficient. More importantly, it seems that biological neurons seem to rather use a simple gain control mechanism. This can be implemented by modifying the gain γ i of the slope of the ReLu function to operate a gradient descent on the cost based on the distribution of each coefficients. Such strategy can be included in the SHL algorithm by replacing line 8 in Algorithm 2.

For instance, the strategy of (Olshausen & Field, 1997 ) assumes a cost on the difference between the observed variance of coefficients V i as computed over a set of samples compared to a desired value σ 2 g (and assuming a multiplicative noise parameterized by α) :

EQUATION

This is similar to the mechanisms of gain normalization proposed by Schwartz & Simoncelli (2001) and which were recently shown to provide efficient coding mechanisms by Simoncelli (2017). However, compared to these methods which manipulate the gain of dictionaries based on the energy of coefficients, we propose to rather use a methodology based on the probability of activation. Indeed, the main distortion that occurs during learning is on higher statistical moments rather than variance, for instance when an atom is winning more at the earlier iterations, its pdf will typically be more kurtotic than a filter that has learned less.

Recently, such an approach was proposed by Sandin & Martin-del Campo (2017) . Based on the same observations, the authors propose to optimize the coding during learning by modulating the gain of each dictionary element based on the recent activation history. They base their Equalitarian Matching Pursuit (EMP) algorithm on a heuristics which cancels the activation of any filter that was more often activated than a given threshold probability (parameterized by 1 + α h ). In our setting, we may compute a similar algorithm using an evaluation of probability of activation followed by binary gates:

p i ← (1 -η h ) • p i + η h • 1/K k=1•••K δ(a i,k > 0) and γ i = δ(p i < N 0 /N * (1 + α h )) (6)

Under review as a conference paper at ICLR 2019 Interestingly, they reported that such a simple heuristic could improve the learning, deriving a similar result as we have shown in Figure 1 and Figure 2 . Again, such strategy can be included in Algorithm 2.

Similarly, we may derive an approximate homeostasis algorithm based on the current activation probability but using a gradient descent approach on gain modulation. Ideally, this corresponds to finding γ i such that we minimize the entropy -i=1•••N p i log p i . However, the sparse coding function S(y k ; Ψ = {Φ, z, N 0 }) is not differentiable. One possible heuristic is then to differentiate the change of modulation gain that would be necessary to achieve an equiprobable probability, that is when ∀i, p i = p 0 def.

= N 0 /N :

EQUATION

We will coin this variant of the algorithm Homeostasis on Activation Probability (HAP). Following these derivations, we quantitatively compared OLS, EMP and HAP to HEH (see Figure 3 ). This shows that while EMP slightly outperforms OLS (which itself is more efficient than None, see Figure 2-B ), HAP proves to be closer to the optimal solution given by HEH. In particular, we replicated in HAP the result of Sandin & Martin-del Campo (2017) that while homeostasis was essential in improving unsupervised learning, the coding algorithm (MP versus OMP) mattered relatively little (see Annex 2.4). Also, we verified the dependence of this efficiency with respect to different hyperparameters (as we did in Figure 2-B ). These quantitative results show that the HEH algorithm could be replaced by a simpler and more rapid heuristic, HAP, which is based on activation probability. This would generate similar efficiency for the coding of patches from natural images.

3. DISCUSSION AND CONCLUSION:

One core advantage of sparse representations is the efficient coding of complex signals using compact codes. Inputs are thus represented as combination of few elements drawn from a large dictionary of atoms. As a consequence, a common design for unsupervised learning rules relies on a gradient descent over a cost measuring representation quality with respect to sparseness. This constraint introduces a competition between atoms. In the context of the efficient processing of natural images, we proposed here that such strategies can be optimized by including a proper homeostatic regulation enforcing a fair competition between the elements of the dictionary. We implemented this rule by introducing a non-linear gain normalization similar to what is observed in biological neural networks.

We validated this theoretical insight by challenging this adaptive unsupervised learning algorithm with alternate heuristics for homeostasis. Simulations show that at convergence, while the coding accuracy did not vary much, including homeostasis changed qualitatively the learned features. In particular, homeostasis results in a more homogeneous set of orientation selective filters, which is closer to what is found in the visual cortex of mammals (Ringach, 2002; Rehn & Sommer, 2007; Loxley, 2017) . To further validate these results, we quantitatively compared the efficiency of the different variants of the algorithms, both at the level of homeostasis (homeostatic learning rate, parameters of the heuristics), but also to the coding (by changing M , N or N 0 ) and to the learning (by changing the learning rate, the scheduling or M ). This demonstrated that overall, this neuro-inspired homeostatic algorithm provided with the best compromise between efficiency and computational cost.

In summary, this biologically-inspired learning rule demonstrates that principles observed in neural computations can help improve real-life machine learning algorithms. Indeed, by developing this fast learning algorithm, we hope for its rapid application in artificial intelligence algorithms. This type of architecture is economical, efficient and fast. It makes it possible to be transferred to most deep learning algorithms. Along with this, we hope that this new type of rapid unsupervised learning algorithm can provide a normative theory for the coding of information in low-level sensory processing, whether it is visual or auditory, for example. Moreover, by its nature, this algorithm can easily be extended to convolutional networks such as those used in deep learning neural networks. This extension is possible by extending the filter dictionary by the hypothesis of invariances to the translation of representations. Our results on different databases show the stable and rapid emergence of characteristic filters on these different bases (see Figure 4 and Annex 3.1). This result shows a probable prospect of extending this representation and for which we hope to obtain classification results superior to the algorithms existing in the state-of-the-art. As such, empirical evaluations of the proposed algorithms should be extended. For instance, it would be very useful to test for image classification results on standard benchmark datasets. We extend the HAP algorithm to a single layered CNN with 20 kernels and using the ATT face database. We show here the kernels learned without (None, top row) and with (HAP, bottom row) homeostasis (note that we used the same initial conditions). As for the simpler case, we observe a heterogeneity of activation counts without homeostasis, that is, in the case which simply normalizes the energy of kernels (see (A)).

With homeostasis, we observe the convergence of the activation probability for the different kernels (see (B)). This demonstrates that this heuristic extends well to a CNN architecture.

Figures for "An adaptive algorithm for unsupervised learning":

In plt.show() ----------Non ---------size of dictionary = (number of filters, size of imagelets) = ( 44 'eta', 'alpha_homeo', 'eta_homeo', ' panel A: plotting some dictionaries In [38] : pname = '/tmp/panel_A' #pname = fname + '_A' ----------OLS ---------size of dictionary = (number of filters, size of imagelets) = (44 1, 324) average of filters = -4.089243933727358e-06 +/-0.001241096006797 0878 average energy of filters = 1.0 +/-3.9562611248144994e-17 ----------HEH ---------size of dictionary = (number of filters, size of imagelets) = (44 1, 324) average of filters = -6.6112572753952305e-06 +/-0.00121065448870 92556 average energy of filters = 1.0 +/-3.700743415417188e-17 ----------EMP ---------size of dictionary = (number of filters, size of imagelets) = (44 1, 324) average of filters = 4.993730484951632e-05 +/-0.0012218228270885 788 average energy of filters = 1.0 +/-3.700743415417188e-17 ----------HAP ---------size of dictionary = (number of filters, size of imagelets) = (44 1, 324) average of filters = -2.429586935952582e-05 +/-0.001195729444507 5826 average energy of filters = 1.0 +/-3.775513461943296e-17

In [39]: subplotpars = dict( left=0.042, right=1., bottom=0., top=1., wspace =0.05, hspace=0. 

getting help:

In [ ]: help(shl)

In [ ]: help(dico)

loading a database:

Loading patches, with or without mask:

In 

Testing different number of neurons and sparsity:

As suggested by AnonReviewer3, we have tested how the convergence was modified by changing the number of neurons. By comparing different numbers of neurons we could re-draw the same figures for the convergence of the algorithm as in our original figures. In addition, we have also checked that this result will hold on a range of sparsity levels. In particular, we found that in general, increasing the l0_sparseness parameter, the convergence took progressively longer. Importantly, we could see that in both cases, this did not depend on the kind of homeostasis heuristic chosen, proving the generality of our results. This is shown in the supplementary material that we have added to our revision ("Testing different number of neurons and sparsity") . This useful extension proves the originality of our work as highlighted in point 4, and the generality of these results compared to the parameters of the network.

In 

Version used:

In [ ]: %load_ext version_information %version_information numpy, shl_scripts

version control:

In [ ]: !git status

In [ ]: !git pull

In [ ]: !git commit -am' {tag} : re-running notebooks'

In [ ]: !git push

exporting the notebook:

In [ ]: !jupyter nbconvert Annex.ipynb

In [ ]: #!jupyter-nbconvert --template report --to pdf Annex.ipynb

In [ ]: !pandoc Annex.html -o Annex.pdf

In [ ]: !zip Annex.zip Annex.html Done. Thanks for your attention!

:

We use image patches drawn from large images of outdoor scenes, as provided in the kodakdb database which is available in the code's repository.

Figures and Tables:

Figure 1: Role of homeostasis in learning sparse representations: We show the results of Sparse Hebbian Learning using different homeostasis algorithms at convergence (1024 learning epochs).The compared algorithms are : None (using a simple simple normalization of the atoms), OLS (the method of(Olshausen & Field, 1997)), HEH (using the optimal homeostasis described in this paper). (A) For each algorithm, we show 18 atoms from the 441 filters of the same size as the image patches (M = 22 × 22 = 484, circularly masked) and presented in a matrix (separated by a white border). The upper and lower row respectively show the least and most probably selected atoms. This highlights the fact that without proper homeostasis, dictionary learning leads to inhomogeneous representations. (B) Evolution of cost F (in bits, see Eq. 4) as a function of the number of iterations and cross-validated over 10 runs. While OLS provides a similar convergence than None, the HEH method provides a better final convergence.

Figure2: Histogram Equalization Homeostasis and its role in unsupervised learning: (A) Nonlinear homeostatic functions z i , ∀i learned using Hebbian learning. These functions were computed for different homeostatic strategies (None, OLS or HEH) but only used in HEH. Note that for our choice of N 0 = 13, all cumulative functions start around 1 -N 0 /N ≈ .970. At convergence of HEH, the probability of choosing any filter is equiprobable, while the distribution of coefficients is more variable for None and OLS. As a consequence, the distortion between the distributions of sparse coefficients is minimal for HEH, a property which is essential for the optimal representation of signals in distributed networks such as the brain. (B) Effect of learning rate η (eta) and homeostatic learning rate η h (eta homeo) on the final cost as computed for the same learning algorithms but with different homeostatic strategies (None, OLS or HEH). Parameters were explored around a default value, on a logarithmic scale and over 4 octaves. This shows that HEH is robust across a wide range of parameters.

representation vector a = S(y k ; Ψ = {Φ, z, N 0 }) (see Algorithm 1), 7:

Figure3: Homeostasis on Activation Probability (HAP) and a quantitative evaluation of homeostatic strategies: (A) 18 from the 441 dictionaries learned for the two heuristics EMP and HAP and compared to the optimal homeostasis (see Figure1-A, HEH). Again, the upper and lower row respectively show the least and most probably selected atoms. (B) Comparison of the cost F during learning and cross-validated over 10 runs: The convergence of OLS is similar to EMP. The simpler HAP heuristics gets closer to the more demanding HEH homeostatic rule, demonstrating that this heuristic is a good compromise for fast unsupervised learning.

Figure 4: Extension to Convolutional Neural Networks (CNNs):We extend the HAP algorithm to a single layered CNN with 20 kernels and using the ATT face database. We show here the kernels learned without (None, top row) and with (HAP, bottom row) homeostasis (note that we used the same initial conditions). As for the simpler case, we observe a heterogeneity of activation counts without homeostasis, that is, in the case which simply normalizes the energy of kernels (see (A)). With homeostasis, we observe the convergence of the activation probability for the different kernels (see (B)). This demonstrates that this heuristic extends well to a CNN architecture.

figure 2: Histogram Equalization Homeostasis

figure 4: Convolutional Neural Network



Table content: <table><tr><td colspan="2">spine.set_color('none') # don't draw spine</td></tr><tr><td>In [6]: shl?</td><td/></tr><tr><td colspan="2"># turn off ticks where there is no spine</td></tr><tr><td>if 'left' in spines:</td><td/></tr><tr><td colspan="2">ax.yaxis.set_ticks_position('left')</td></tr><tr><td>else:</td><td/></tr><tr><td># no yaxis ticks</td><td/></tr><tr><td>ax.yaxis.set_ticks([])</td><td/></tr><tr><td>if 'bottom' in spines:</td><td/></tr><tr><td colspan="2">ax.xaxis.set_ticks_position('bottom')</td></tr><tr><td>else:</td><td/></tr><tr><td># no xaxis ticks</td><td/></tr><tr><td>ax.xaxis.set_ticks([])</td><td/></tr><tr><td colspan="2">[1]: %load_ext autoreload import matplotlib In [7]: print('number of patches, size of patches = ', data.shape)</td></tr><tr><td colspan="2">%autoreload 2 pylab_defaults = { print('average of patches = ', data.mean(), ' +/-', data.mean(axis</td></tr><tr><td>'font.size': 10, =1).std())</td><td/></tr><tr><td colspan="2">'xtick.labelsize':'medium', SE = np.sqrt(np.mean(data**2, axis=1)) In [2]: import numpy as np # 'font.sans-serif' : ['Helvetica'], np.random.seed(seed) # 'font.family' : 'sans-serif', seed = 42 'text.usetex': False, np.set_printoptions(precision=2, suppress=True) 'ytick.labelsize':'medium', print('average energy of data = ', SE.mean(), '+/-', SE.std())</td></tr><tr><td>}</td><td/></tr><tr><td colspan="2">In [3]: # some overhead for the formatting of figures</td></tr><tr><td colspan="2">import matplotlib.pyplot as plt #matplotlib.rcParams.update({'font.size': 18, 'font.family': 'STIXG</td></tr><tr><td colspan="2">eneral', 'mathtext.fontset': 'stix'})</td></tr><tr><td colspan="2">matplotlib.rcParams.update(pylab_defaults) In [8]: #!ls -l {shl.cache_dir}/{tag}*</td></tr><tr><td colspan="2">#matplotlib.rcParams.update({'text.usetex': True}) !ls {shl.cache_dir}/{tag}*lock*</td></tr><tr><td colspan="2">!rm {shl.cache_dir}/{tag}*lock*</td></tr><tr><td>#!rm {shl.cache_dir}/{tag}*</td><td/></tr><tr><td colspan="2">fig_width_pt = 318.670 # Get this from LaTeX using \showthe\column #!ls -l {shl.cache_dir}/{tag}*</td></tr><tr><td>width</td><td/></tr><tr><td colspan="2">fig_width_pt = 450 # Get this from LaTeX using \showthe\columnwidt</td></tr><tr><td>h DEBUG = True</td><td/></tr><tr><td>#fig_width_pt = 1024 #221 DEBUG = False</td><td># Get this from LaTeX using \showthe\</td></tr><tr><td colspan="2">columnwidth / x264 asks for a multiple of 2 hl, hs = 10*'-', 10*' '</td></tr><tr><td colspan="2">ppi = 72.27 # (constant) definition of the ppi = points per inch</td></tr><tr><td colspan="2">inches_per_pt = 1.0/ppi # Convert pt to inches In [4]: tag = 'ICLR' #inches_per_cm = 1./2.54 datapath = '../../SparseHebbianLearning/database' fig_width = fig_width_pt*inches_per_pt # width in inches # different runs grid_fig_width = 2*fig_width #opts = dict(datapath=datapath, verbose=0) phi = (np.sqrt(5) + 1. ) /2 #opts = dict(cache_dir='cache_dir_cluster', datapath=datapath, verb #legend.fontsize = 8 ose=0) #fig_width = 9 #opts = dict(cache_dir='cache_dir_ICLR', datapath=datapath, verbose fig_height = fig_width/phi =0) figsize = (fig_width, fig_height) opts = dict(cache_dir='cache_dir_cluster25', eta=0.002, eta_homeo=0</td></tr><tr><td colspan="2">.005, datapath=datapath, verbose=0)</td></tr><tr><td>def adjust_spines(ax, spines):</td><td/></tr><tr><td colspan="2">In [5]: from shl_scripts.shl_experiments import SHL</td></tr><tr><td>shl = SHL(**opts)</td><td/></tr><tr><td colspan="2">data = shl.get_data(matname=tag)</td></tr></table>



Table content: <table><tr><td colspan="2">In [10]: homeo_methods = ['None', 'OLS', 'HEH']</td></tr><tr><td colspan="2">Type: list_figures = ['show_dico', 'time_plot_error', 'time_plot_logL', ' SHL</td></tr><tr><td colspan="2">String form: &lt;shl_scripts.shl_experiments.SHL object at 0x10919412 time_plot_MC', 'show_Pcum']</td></tr><tr><td colspan="2">8&gt; list_figures = []</td></tr><tr><td colspan="2">File: dico = {}</td></tr><tr><td colspan="2">~/science/SparseHebbianLearning/shl_scripts/shl_experiments.py for i_cv in range(N_cv):</td></tr><tr><td colspan="2">Docstring: dico[i_cv] = {}</td></tr><tr><td colspan="2">Base class to define SHL experiments: for homeo_method in homeo_methods:</td></tr><tr><td/><td>-initialization shl = SHL(homeo_method=homeo_method, seed=seed+i_cv, **opts</td></tr><tr><td>)</td><td>-coding and learning</td></tr><tr><td/><td>-visualization dico[i_cv][homeo_method] = shl.learn_dico(data=data, list_f</td></tr><tr><td colspan="2">-quantitative analysis igures=list_figures, matname=tag + '_' + homeo_method + '_seed=' +</td></tr><tr><td colspan="2">str(seed+i_cv))</td></tr><tr><td colspan="2">list_figures = ['show_dico']</td></tr><tr><td colspan="2">for i_cv in [one_cv]:</td></tr><tr><td/><td>for homeo_method in homeo_methods:</td></tr><tr><td/><td>print(hl + hs + homeo_method[:3] + hs + hl)</td></tr><tr><td/><td>shl = SHL(homeo_method=homeo_method, seed=seed+i_cv, **opts</td></tr><tr><td colspan="2">number of patches, size of patches = (65520, 324) )</td></tr><tr><td colspan="2">average of patches = 5.0641928164665185e-19 +/-0.0095770518654 shl.learn_dico(data=data, list_figures=list_figures, matnam</td></tr><tr><td colspan="2">37931 e=tag + '_' + homeo_method + '_seed=' + str(seed+i_cv))</td></tr><tr><td colspan="2">average energy of data = 0.29851622590347293 +/-0.08935954499531</td></tr><tr><td>101</td><td>print('size of dictionary = (number of filters, size of ima</td></tr><tr><td colspan="2">gelets) = ', dico[i_cv][homeo_method].dictionary.shape)</td></tr><tr><td/><td>print('average of filters = ', dico[i_cv][homeo_method].di</td></tr><tr><td colspan="2">ctionary.mean(axis=1).mean(),</td></tr><tr><td/><td>'+/-', dico[i_cv][homeo_method].dictionary.mean(axis</td></tr><tr><td colspan="2">=1).std())</td></tr><tr><td/><td>SE = np.sqrt(np.sum(dico[i_cv][homeo_method].dictionary**2,</td></tr><tr><td colspan="2">axis=1))</td></tr><tr><td colspan="2">ls: cache_dir_cluster25/ICLR*lock*: No such file or directory print('average energy of filters = ', SE.mean(), '+/-', SE.</td></tr><tr><td colspan="2">rm: cache_dir_cluster25/ICLR*lock*: No such file or directory std())</td></tr><tr><td colspan="2">In [9]: fname = 'figure_map'</td></tr><tr><td colspan="2">N_cv = 10</td></tr><tr><td colspan="2">one_cv = 0 # picking one to display intermediate results</td></tr></table>



Table content: <table><tr><td>In [ ]: shl = SHL(one_over_F=False, **opts)</td></tr><tr><td>dico_w = shl.learn_dico(data=data, matname=tag + '_WHITE', list_fig</td></tr><tr><td>ures=[])</td></tr><tr><td>shl = SHL(one_over_F=True, **opts)</td></tr><tr><td>dico_1oF = shl.learn_dico(data=data, matname=tag + '_OVF', list_fig</td></tr><tr><td>ures=[])</td></tr><tr><td>fig_error, ax_error = None, None</td></tr><tr><td>fig_error, ax_error = shl.time_plot(dico_w, variable='F', fig=fig_e</td></tr><tr><td>rror, ax=ax_error, color='blue', label='white noise')</td></tr><tr><td>fig_error, ax_error = shl.time_plot(dico_1oF, variable='F', fig=fig</td></tr><tr><td>_error, ax=ax_error, color='red', label='one over f')</td></tr><tr><td>#ax_error.set_ylim((0, .65))</td></tr><tr><td>ax_error.legend(loc='best')</td></tr><tr><td>[ ]: N_patches = 12 from shl_scripts.shl_tools import show_data Testing two different learning rates strategies</td></tr><tr><td>opts_ = opts.copy() We use by defaut the strategy of ADAM, see https://arxiv.org/pdf/1412.6980.pdf opts_.update(verbose=0) for i, (do_mask, label) in enumerate(zip([False, True], ['Without m (https://arxiv.org/pdf/1412.6980.pdf)</td></tr><tr><td>ask', 'With mask'])):</td></tr><tr><td>data_ = SHL(DEBUG_DOWNSCALE=1, N_patches=N_patches, n_image=1, In [ ]: shl = SHL(beta1=0., **opts) do_mask=do_mask, seed=seed, **opts_).get_data() dico_fixed = shl.learn_dico(data=data, matname=tag + '_fixed', list fig, axs = show_data(data_) _figures=[]) axs[0].set_ylabel(label); shl = SHL(**opts) plt.show() dico_default = shl.learn_dico(data=data, matname=tag + '_default',</td></tr><tr><td>list_figures=[])</td></tr><tr><td>fig_error, ax_error = None, None</td></tr><tr><td>fig_error, ax_error = shl.time_plot(dico_fixed, variable='F', fig=f Testing different algorithms ig_error, ax=ax_error, color='blue', label='fixed')</td></tr><tr><td>fig_error, ax_error = shl.time_plot(dico_default, variable='F', fig</td></tr><tr><td>=fig_error, ax=ax_error, color='red', label='ADAM')</td></tr><tr><td>#ax_error.set_ylim((0, .65))</td></tr><tr><td>ax_error.legend(loc='best')</td></tr><tr><td>mp']</td></tr><tr><td>: # 'threshold', 'lasso_cd',</td></tr><tr><td>opts_ = opts.copy()</td></tr><tr><td>opts_.update(homeo_method=homeo_method, learning_algorithm=</td></tr><tr><td>algorithm, verbose=0)</td></tr><tr><td>shl = SHL(**opts_)</td></tr><tr><td>dico= shl.learn_dico(data=data, list_figures=[],</td></tr><tr><td>matname=tag + ' -algorithm={}'.format(algor</td></tr><tr><td>ithm) + ' -homeo_method={}'.format(homeo_method))</td></tr><tr><td>fig, ax = shl.time_plot(dico, variable='F', fig=fig, ax=ax,</td></tr><tr><td>label=algorithm +'_' + homeo_method)</td></tr><tr><td>ax.legend()</td></tr><tr><td>White Noise Initialization + Learning</td></tr></table>

for i_ax, variable in enumerate(variables): for color, homeo_method in zip(colors, homeo_methods): opts_ = opts.copy() opts_.update(homeo_method=homeo_method, datapath=datapath) experiments = SHL_set(opts_, tag=tag + '_' + homeo_method) fig, axs[i_ax] = experiments.scan(variable=variable, list_f igures=[], display='final', fig=fig, ax=axs[i_ax], color=color, dis play_variable='F', verbose=0) #, label=homeo_metho axs[i_ax].set_xlabel('') #variable axs[i_ax].text(.1, .8, variable, transform=axs[i_ax].trans

Table content: <table><tr><td>[ ]: from shl_scripts.shl_experiments import SHL_set</td></tr><tr><td>homeo_methods = ['None', 'OLS', 'HEH']</td></tr><tr><td>homeo_methods = ['None', 'EMP', 'HAP', 'HEH', 'OLS']</td></tr><tr><td>variables = ['l0_sparseness', 'n_dictionary']</td></tr><tr><td>list_figures = []</td></tr><tr><td>#n_dictionary=21**2</td></tr><tr><td>for homeo_method in homeo_methods:</td></tr><tr><td>opts_ = opts.copy()</td></tr><tr><td>opts_.update(homeo_method=homeo_method, datapath=datapath)</td></tr><tr><td>experiments = SHL_set(opts_, tag=tag + '_' + homeo_method)</td></tr><tr><td>experiments.run(variables=variables, n_jobs=1, verbose=0)</td></tr><tr><td>fig, axs = plt.subplots(len(variables), 1, figsize=(fig_width/2, fi</td></tr><tr><td>g_width/(1+phi)), gridspec_kw=subplotpars, sharey=True)</td></tr><tr><td>Axes)</td></tr><tr><td>#axs[i_ax].get_xaxis().set_major_formatter(matplotlib.ticke</td></tr><tr><td>r.ScalarFormatter())</td></tr></table>

Peer review (iclr2019):
 This paper discusses the addition of a regularizer to a standard sparse coding/dictionary learning algorithm to encourage the atoms to be used with uniform frequency. I do not think this work should be accepted to the conference for the following reasons: 1: The authors show no benefit of this scheme except perhaps faster convergence. If faster training of dictionary learning models was a bottleneck in practical applications, this might be of interest, but it is not. SPAMS (http://spams-devel.gforge.inria.fr/) can train a model on image patches as the authors do here in a few tens of seconds on a modern computer. On the other hand, the authors give no evidence, empirical or otherwise, that their method is useful on any downstream tasks. In my view, they do not even show that the distribution of atom usage will be better with their algorithm after the learning has converged, as at least according to their learning curves, the baselines have not finished converging. It is not even clear that the final compression of the baselines would not be better. Even if they did show these convincingly, it is not obvious to me that it is valuable; the authors need to *show* that uniform usage is desirable. 2: The authors should compare against several costs/algorithms (e.g. l_0 with OMP, l_1 with LARS, etc.), and across various N_0/sparsity penalties, and across several datasets. The empirical evaluation is quite weak- one sparsity setting, two baselines, one dataset. Even without the "train to convergence" question above, I don*t think the authors have demonstrated that their claims on the properties of their algorithms/formulations are generally true.