{'fraction_non_alphanumeric': 0.035840938105820605, 'fraction_numerical': 0.0078490708780557, 'mean_word_length': 4.596718708653083, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Early Inference in Energy-Based Models Approximates Back-Propagation

We show that Langevin MCMC inference in an energy-based model with latent variables has the property that the early steps of inference, starting from a stationary point, correspond to propagating error gradients into internal layers, similarly to back-propagation. The error that is back-propagated is with respect to visible units that have received an outside driving force pushing them away from the stationary point. Back-propagated error gradients correspond to temporal derivatives of the activation of hidden units. This observation could be an element of a theory for explaining how brains perform credit assignment in deep hierarchies as efficiently as back-propagation does. In this theory, the continuous-valued latent variables correspond to averaged voltage potential (across time, spikes, and possibly neurons in the same minicolumn), and neural computation corresponds to approximate inference and error back-propagation at the same time.


Introduction
It has been hypothesized numerous times (Hinton and Sejnowski, 1986;Friston and Stephan, 2007;Berkes et al., 2011), that, given a state of sensory information (current and past inputs), neurons are collectively performing inference, i.e., moving towards configurations that better "explain" the observed sensory data. We can think of the configuration of internal neurons (hidden units or latent variables) as an "explanation" for the observed sensory data.
Under this hypothesis, when an unexpected signal arrives at visible neurons 1 , the other neurons will change their state (from one that is near a stochastic equilibrium) so as to better reflect that change in input. If this error signal is slowly driving visible neurons towards the externally observed value, this creates a perturbation for the whole network. In this paper, we consider what happens early in this process, when the perturbation due to the incorrect prediction of visible neurons is propagated towards inner areas of the brain (hidden units and hidden layers). We show how the propagation of this perturbation is mathematically equivalent to the propagation of activations gradients by the back-propagation algorithm in a deep neural network.
This result assumes that latent variables are continuousvalued (unlike in the Boltzmann machine), making the system an energy-based model similar to the continuous Hopfield net (Hopfield, 1984), defined by an energy function over both continuous visible and hidden units. It also assumes that the neurons are leaky integrators with injected noise. As a consequence neural computation corresponds to gradually going down the energy while adding noise, which corresponds to running a Langevin Monte-Carlo Markov chain (MCMC) as inference mechanism.
This work is only a stepping stone. Several points still need to be elucidated before a complete theory of learning, inference and credit assignment is elaborated that is both biologically plausible and makes sense from a machine learning point of view for whole networks (with global optimization of the whole network and not being limited to learning of visible neurons that receive a target). In particular, energy based models require symmetry of connections, but note that hidden units in the model need not correspond exactly to actual neurons in the brain (it could be groups of neurons in a cortical microcircuit, for example). It remains to be shown how a form of symmetry could arise from the learning procedure itself, for example because of an unsupervised learning objective similar to the one used to train autoencoders. As shown in this paper, synaptic changes approximately proportional to the stochastic gradient with respect to the prediction error could be achieved if the synaptic updates follow the learning rule studied by Bengio et al. (2015b) in order to mimic spike-timing dependent plasticity (STDP).

Neural computation does inference: going down the energy
We consider the hypothesis that a central interpretation of neural computation (what neurons do, on a short time scale at which weights can be considered fixed) is that it is performing iterative inference. Iterative inference means that the hidden units h of the network are gradually changed towards configurations that are more probable, with the given sensory input x and according to the current "model of the world" associated with the parameters of the model. In other words, they are approximately moving towards configurations more probable under P (h|x), and eventually sampling from P (h|x).
Before making the connection between Boltzmann machines or energy-based models and back-propagation, let us see in more mathematical detail how neural computation could be interpreted as inference.

Leaky integrator neuron as Langevin MCMC
For this purpose, consider the classical leaky integrator neural equation. Let s t represent the state of the system at time t, a vector with one element per unit, where s t,i is a real valued quantity associated with the i-th unit, corresponding to a time-integrated voltage potential.
Let us denote x t for the visible units, a subset of the elements of s t (i.e., the externally driven input) and h t for the hidden units, i.e., s t = (x t , h t ). Let f be the function that computes the new value of the complete state s t given its previous value, with f = (f x , f h ) to denote the parts of f that respectively outputs the predictions on the clamped (visible) units and on the unclamped (hidden) units. The time evolution of the unclamped units is assumed to follow a leaky integration equation, i.e., where R(s) = (R x (s), R h (s)) represents the networkgenerated pressure on neurons, i.e., R i (s) is what the rest of the network asks neuron i to move towards, and the corrupted states is results from synaptic noise and spiking effects. Here we roughly model this corruption by simple additive noise:s and we see that the above equation corresponds to the discretization of a differential equation which brings h exponentially fast towards the "target" value R h (s), along with the random walk movements brought by the noise η. We assume that R i (s) is a weighted sum of input signals coming from the neurons connected into neuron i, although the derivation below does not depend on the specific form of R, only on the assumption that it corresponds to an energy gradient, an idea developed below.

Machine Learning Interpretation
but injecting noise is important in order to find not just a single local mode of P (h|x) but explore the full distribution.
We now draw an interesting link with recent work on unsupervised learning using denoising auto-encoders and denoising score matching (Vincent, 2011;Alain and Bengio, 2013). If R(s) is the linear combination of input rates ρ(s), the above papers make a link between R(s) − s and the energy of a probabilistic model P (s) ∝ e −E(s) with energy function E(s), i.e., they find that With this interpretation, the leaky integration neural computation of Eq. 1 seems to follow a Langevin Monte-Carlo Markov chain (Andrieu et al., 2003): where for the last line we used Eqs. 2 and 3 and we see that we are going down the gradient froms t . Hence from the point of view of the noisy statess, we see that the update equation corresponds tõ which we recognize as going down the gradient of the energy with "learning rate" ǫ and adding "noise" η t+1 − (1 − ǫ)η t .

A Possible Energy Function
To fix ideas and illustrate the possibility of a driving function R corresponding to the gradient of an energy function, we propose the following energy function, closely related to the Boltzmann machine energy function, but with continuous non-linearities inserted.
where W i,j is the weight between unit j to unit i, and ρ is the neural non-linearity, some kind of monotonic bounded function which outputs a value between 0 and 1 corresponding to a firing rate, as for example with thresholds β 1 and β 2 such that β 2 − β 1 = 1 and β 1 < 0 < β 2 .
With this energy function, the driving function R would be To obtain this, we have assumed that Otherwise, we would get that , which automatically symmetrizes the weight matrix.
This formula for the driving input R is similar to the usual weighted sum of firing rates, except for the novel factor ρ ′ (s i ), which would suggest that when a neuron is saturated (either being shut off or firing at the maximal rate), the external inputs have no impact on its state. The only term that remains in the neural update equation (Eq. 1) is the one that drives the state towards 0, i.e., bringing it out of the saturation region and back into a regime where the neuron is sensitive to the outside feedback, so long as ρ (0) is not a saturated value. This idea is developed further below.

Fixed Point Behavior
In particular, it is interesting to note what happens around a fixed point of the state dynamics.
which means that Let us consider the hypothesis where the unit is saturated, i.e., ρ ′ (s i ) ≈ 0. In that case, R i (s) = 0 and the neural update becomes which converges towards s i = 0. If the origin corresponds to a region where the derivative is significantly nonzero, |ρ ′ (0)| > 0, we get that when ρ ′ (s i ) = 0, the network cannot be at a fixed point. Otherwise the state would move towards 0 and thus it could not have been at a fixed point.

Link to Back-propagation
We are now ready to present the main result of this paper, i.e., a link between neural computation as inference in an energy-based model and back-propagation of prediction error gradients.

Propagation of Pertubations
Consider what happens when a network such as described above sits near equilibrium, i.e., near a fixed point as per Eq. 9. At that point, as per that equation, the average gradient of the energy is 0 and weight updates are also 0 in average.
To make the link to supervised back-propagation simpler, let us consider two kinds of visible units: input units x and output units y, and s = (x, y, h). Suppose that we start by letting the network settle to a fixed point with x clamped to the observed input values. Then we obtain an output y at the fixed pointŝ, where R y (ŝ) =ŷ, R h (ŝ) =ĥ.
Equivalently, we have that i.e., the "free" units (hidden and output) have settled to a value that is in agreement with the clamped input units. Now suppose that a target value y is observed and gradually drives the output units from their fixed point valueŷ, towards y. This happens because the output units are also leaky integrator neurons, meaning that their state gradually changes based on the input they receive, in direction of the driving signal (now y rather than R y (ŝ). Let us denote that initial change ofŷ when going from time step 0 to time step 1 (following the neural update equation 1, Let us consider as training objective the squared prediction error but with R y (ŝ) replaced by y). That would push the global stateŝ away from the equilibrium where it was sitting, and into a region where ∂E(s) ∂s is non-zero.
which corresponding to the mismatch between the prediction R y (ŝ) and the target value y driving the output units, or equivalently because at equilibrium R y (ŝ) =ŷ. Note that where ǫ can be seen as a learning rate if we were trying to do SGD onŷ directly. Now, how would the rest of the network react to this external perturbation? Each hidden neuron would approximately move in the direction of the gradient of C, but only those (call them h 1 ) that are directly connected to the output would initially feel the pressure to minimize C. That perturbation (in the form of a volley of additional spikes, for real neurons) would then travel to the next circle of neurons, those directly connected to h 1 but not to the output, etc.
Let us look at this in more detail. Consider a typical multi-layer architecture with connections between the output layer (yielding valueŷ at equilibrium) and the top hidden layer (yielding valueĥ 1 ), between the top and the next hidden layer (yieldingĥ 1 andĥ 2 ), etc. The change ∆y would propagate toĥ 1 via the neural update, which, when we ignore the effect of the injected noise, would yield a change inĥ 1 ∆h 1 = ǫ(R h1 ((x,ŷ + ∆y,ĥ)) −ĥ 1 ) .
With ∆y small (arising out of our assumption that the visible units only gradually move towards their target), we can approximate the above by taking the Taylor expansion of R aroundŷ, R h1 ((x,ŷ + ∆y,ĥ)) =ĥ 1 + ∂R h1 (ŝ) ∂ŷ ∆y + o(ǫ) . (17) exploiting R h1 (ŝ) =ĥ 1 at the fixed point, yielding Hence we have that Note, that with the assumed layer-wise neural network structure (i.e. having no connections between neurons insight one layer) o(ǫ) (and o(ǫ 2 ) respectively) is zero, if the neural non-linearity ρ fulfills ρ ′′ (s i ) = 0, which is the case for the non-linearity given in Eqs. 7. This is a consequence of having all the higher derivatives of ρ being zero, since ρ is piecewise linear.
In order to obtain backprop, what we would like to get, though is not ∂R h 1 (ŝ) ∂ŷ ∆y but ∂Ry(ŝ) ∂ĥ1 ∆y since that would correspond to an application of the chain rule and we would have ∆h 1 ∝ ∂C ∂ĥ1 . The good news is that this equality is true because R is a first derivative of a function related to the energy function: where and R y (ŝ) = ∂L(ŝ) ∂ŷ , R h1 (ŝ) = ∂L(ŝ) ∂ĥ1 , so that ∂R h 1 (ŝ) ∂ŷ and ∂Ry(ŝ) ∂ĥ1 are cross-derivatives of L. As we know that cross derivatives are symmetric, Now note that since at the fixed pointŷ = R y (ŝ) and we are ready to exploit that to rewrite ∆h 1 (Eq. 19) in a form that equates it with a backpropagated gradient: Similarly, the perturbation ∆h 1 will be transmitted at the next time step to the unitsĥ 2 that are directly connected tô h 1 (but not toŷ), and yield

Stochastic Gradient Descent Weight Update
The above result is consistent with and inspired by the idea previously proposed by Hinton (2007) that temporal change can encode back-propagated gradients. What would it take for the ∆ĥ k at layer k to turn into a stochastic gradient descent (SGD) weight update with respect to the prediction error ||y −ŷ|| 2 ? Since the state changeṡ represents the gradient of the prediction error with respect to s, SGD on W i,j would require the weight change ∆W i,j being proportional to the rate of change of the state of the post-synaptic neuron,ṡ i and proportional to the gradient of the fixed point state s with respect to W i,j , i.e., to ∂si ∂Wi,j . To first approximation, this is the firing rate of the pre-synaptic neuron, ρ(s j ). This would yield It turns out that such a learning rule allows to simulate the relationship between spike-timing and synaptic change according to the STDP (spike-timing dependent plasticity), as shown via simulations by Bengio et al. (2015a). See Xie and Seung (2000) for a similar learning rule, also related to STDP via a different analysis.
Thus, with this STDP-compatible learning rule, the change in weights due to the initial perturbation would be approximately proportional to the back-propagation update, since it corresponds to ∆h ∂h ∂W . However, note the multiplicative factors ǫ k+2 for unitsĥ k at layer k, that make the initial changes much slower for the more remote layers. This is because the leaky integration neurons have not had time to integrate the information yet, so practically it will take on the order of the time constant times k for the change inĥ k to become significant, unless we adjust the per-layer learning rates accordingly.
Although we see that the proposed neural dynamics and weight updates will behave approximately like backpropagation, there are differences, especially when we consider what happens after more time steps. But maybe the most important take-home message from this link with back-propagation is the following. We know that backpropagation works extremely well to train both supervised and unsupervised networks. We see here that backpropagation essentially corresponds to a variational update when the inference is infinitesimal, i.e., we only allow a single step of inference corresponding to small moves in the direction of reducing the energy function.

Related work, contributions and future work
An important inspiration for this work is the idea proposed by Hinton (2007) that brains could implement backpropagation by using temporal derivatives to represent activation gradients, and the suggestion that combining this assumption with STDP would approximately yield SGD on the synaptic weights.
The idea of neural computation corresponding to a form of stochastic relaxation towards lower energy configurations is of course very old, for example with the Boltzmann machine (Hinton and Sejnowski, 1986) and its Gibbs sampling procedure. For more recent work in this direction, see also (Berkes et al., 2011). What differs here from the Boltzmann machine is that we consider the state space to be continuous (associated with the expected voltage potential, integrating out the random effects due to spikes), rather than discrete, and that we consider very small steps (going down the gradient of the energy), which is more like a Langevin MCMC, rather than allowing each neuron to stochastically jump with higher probabibility to its optimal state, given its neighbors configuration, which is what Gibbs sampling does.
There are of course many other papers on theoretical interpretations of STDP, and the reader can find many references in Markram et al. (2012), but more work is needed to explore the connection of STDP to an unsupervised learning objective that could be used to train not just a single layer network (like PCA and traditional Hebbian updates) but also a deep unsupervised model. Many approaches (Fiete and Seung, 2006;Rezende and Gerstner, 2014) rely on variants of the REINFORCE algorithm (Williams, 1992) to estimate the gradient of a global objective function (basically by correlating stochastic variations at each neuron with the changes in the global objective). Although this principle is simple, it is not clear that it will scale to very large networks due to the linear growth of the variance of the estimator with the number of neurons. It is therefore tempting to explore other avenues, and we hope that the building blocks introduced here and the links made with energy-based approaches with variational inference for unsupervised learning can form useful material for a more efficient unsupervised learning principle for deep networks that is also consistent with STDP.
If the energy function proposed here is closer to a biological truth than the energy defined for continuous Hopfield networks (Hopfield, 1984), we should see that (a) firing rate returns to its baseline when the neuron is saturated (completely turned off or maximally turned on) and (b) synaptic weight changes should also vanish under this condition, as seen by inspection of Eq. 26. It would clearly be interesting to test these predictions in actual biological experiments.
Much remains to be done to obtain a complete probabilistic theory of unsupervised learning that is consistent with STDP, but we believe that we have put interesting ingredients in place. One aspect that requires a lot more development is how the proposed STDP update helps to fit the sensory observations x. If, as hypothesized above, neural computation is approximately doing inference (e.g. Langevin MCMC), then each step of inference, in average, brings us towards an equally likely or even more likely configuration of h, given x, according to the model. Hence each step is approximately pointing down the energy of P (h|x). Now, in an EM or variational EM context such as discussed in Neal and Hinton (1999); Kingma and Welling (2014); Bengio et al. (2015b), with x fixed, the distribution we want to model and consider as a target towards which parameters should be updated is precisely the joint of h ∼ P (h|x) and x ∼ the training data, which we now call Q(h, x) (the inference distribution), following the above papers. We would like the parameters to move in the direction that makes the model more consistent with Q(h, x), which is what is required to maximize the variational EM bound on the data likelihood P (x). The idea is that we change the inference process so that it would reach its final state faster, which corresponds to a configuration of h that fits well the observed x.
Another open question is how to reconcile the need for symmetric weights when we introduce an energy function and the fact that W i,j and W j,i are stored at two physically different places in biological neurons. An encouraging observation is that earlier work on auto-encoders empirically showed that even when the forward and backward weights are not tied, they tend to converge to symmetric values, and in the linear case the minimization of reconstruction error automatically yields symmetric weights (Vincent et al., 2010). Another encouraging piece of evidence, also linked to autoencoders, is the theoretical result from Arora et al. (2015), showing that the symmetric solution minimizes the autoencoder reconstruction error between two successive layers of rectifying (ReLU) units.