Linear-Nonlinear-Poisson Neuron Networks Perform Bayesian Inference On Boltzmann Machines
Authors: Yuanlong Shao
ICLR 2013

Abstract: One conjecture in both deep learning and classical connectionist viewpoint is that the biological brain implements certain kinds of deep networks as its back-end. However, to our knowledge, a detailed correspondence has not yet been set up, which is important if we want to bridge between neuroscience and machine learning. Recent researches emphasized the biological plausibility of Linear-Nonlinear-Poisson (LNP) neuron model. We show that with neurally plausible settings, the whole network is capable of representing any Boltzmann machine and performing a semi-stochastic Bayesian inference algorithm lying between Gibbs sampling and variational inference.

arXiv:1210.8442v3  [cs.AI]  27 Jan 2013
Linear-Nonlinear-Poisson Neuron Networks Perform
Bayesian Inference On Boltzmann Machines
Louis Yuanlong Shao
Department of Computer Science & Engineering
The Ohio State University
shaoyu@cse.ohio-state.edu
Abstract
One conjecture in both deep learning and classical connecti onist viewpoint is
that the biological brain implements certain kinds of deep n etworks as its back-
end. However, to our knowledge, a detailed correspondence h as not yet been set
up, which is important if we want to bridge between neuroscie nce and machine
learning. Recent researches emphasized the biological pla usibility of Linear-
Nonlinear-Poisson (LNP) neuron model. We show that with neu rally plausible
settings, the whole network is capable of representing any Boltzmann machine and
performing a semi-stochastic Bayesian inference algorith m lying between Gibbs
sampling and variational inference.
1 Introduction
Classical connectionist viewpoint has long been inspired f rom how the brain works, such as the
invention of “perceptron” [1], and the “parallel distributed processing” approach [2]. Modern deep
learning and unsupervised feature learning also assume con nections between biological brain and
certain kinds of deep networks, either probabilistic or not . For example, properties of visual area
V2 are found to be comparable to those on the sparse autoencod er networks [3]; the sparse coding
learning algorithm [4] is originated directly from neurosc ience observations; also psychological
phenomenon such as end-stopping is observed in sparse coding experiments [5].
In addition, there are architectural similarity between neuroscience and machine learning. For exam-
ple, neuronal spike propagation may correspond to prediction or inference tasks in a learned model;
synaptic plasticity may correspond to parameter estimatio n given the network structure; neuroplas-
ticity may correspond to learning the network structure together with inventing new hidden units in
the network; axonal path ﬁnding in the help of guidance signals may be regarded as structural priors
making learning network structure easier.
For this reason, it may be helpful again to refer to how the bra in works when dealing with prob-
lems that are currently puzzling machine learning research ers. For example, training deep net-
works was made practical mostly after the breakthrough of de ep learning starting from 2006
(e.g. [6][7][8][9][10], etc.). The success may lie in a safe way of adapting the network structure,
such that the parameter estimation part can be effective eno ugh to reach non-trivial local optima.
However, since the brain is not exactly layer-wise, we may wonder is there any other way of choos-
ing the network structure along with parameter estimation? If we want to transfer neuroscience
knowledge for answering questions as such and build bridge b etween these two areas, we want to
ﬁrst ﬁgure out what exactly is the deep network that the brain is representing.
A good starting point of answering this may be to look at singl e neuron computation. Through
how they transfer information among them, we can reverse engineer both the representation and the
prediction / inference algorithm, which is the main purpose of this paper.
1
In theoretical neuroscience, single neuron models are divided into different levels of detail, (see [11]
for a ﬁve level survey). Detailed models can be quite realist ic and few simpliﬁcation assumptions
hold. Yet effective simpliﬁcation holds in different exten t. For example, integrate-and-ﬁre models
capture essential properties of Hodgkin-Huxley model [12] , and its variant exponential integrate-
and-ﬁre neuron models [13] are found to be capable of reprodu cing spike timing of several types
of cortical neurons. The community has also been investigat ing on the Linear-Nonlinear-Poisson
(LNP) models for long [14][15][16]. A recent study in [17] an alyzed how LNP models effectively
reproduce ﬁring rates in more realistic neurons. Based on th ese results, we start from formally
presenting the LNP model and then turn to presenting a semi-s tochastic inference algorithm on
Boltzmann machines, which is derived from combining Gibbs s ampling and variational inference.
We then make a detailed matching between LNP model and the inf erence algorithm. Since the
semi-stochastic inference algorithm has not been explored in the learning community, we also show
some experiments illustrating its computational property , including its stochastic convergence and
its similarity to variational inference.
2 Brief Review on Neural Plausibility of LNP Model
In this section, we brieﬂy review the Linear-Nonlinear-Poisson neuron model and its neural plausi-
bility, focusing on what modeling options we have to ﬁt it with useful inference algorithms.
The LNP neuron model [14][15][16] formalizes the spike trai n generated by each neuron as a non-
homogeneous Poisson point process, whose rate function over time depends on the input spike trains
from its presynaptic neurons, with certain form of short-term memory in the dependence.
Synaptic Transmission: Iij (τ) = Wij ·(Xj ∗α) (τ) (a)
Dendritic Integration: Ii (τ) = Multi - Linear
(
{Iij (τ)}j∈M(i)
)
(b)
Spike Generation:
{ λi (τ) = σ ((Ii ∗D) (τ)) (c)
Xi (τ) ∼PoissonPointProcess (λi (τ)) (d)
(1)
In all these equations, τ is the continuous time index, whose unit is millisecond. Sub script i is
the index of the neuron in question, M (i) is the index set of presynaptic neurons of neuron i. In
Eq. (1)(d), Xi (τ) = ∑
f δ(τ −τ(f)
i ) represents the spike train, with spikes as Dirac functions l o-
cated at time steps τ(f)
i . λi (τ) is the rate function of the Poisson point process. Eq. (1)(a)represents
the postsynaptic current at a particular synapse with efﬁcacy Wij . α is a certain non-negative func-
tion with certain time course. Evidence for “Poisson-like” statistics in cortical neurons can be found
in [18][19], which implies that spike counts over a certain interval may be Poisson distributed.
Eq. (1)(b) is the dendritic summation. There are different v iewpoints on whether linearity holds
true. Due to phenomenon like mutual inhibition on different dendritic shafts, it is commonly known
that nonlinear summation could happen (see [20][21] for exa mple), yet studies from [22] show
that when dendritic spines are presented, linear summation is true, and this happens frequently in
excitatory mammalian cortex. In addition, a survey on neuron models in [11] summarized different
computations that can be performed in dendritic processing, two notable examples are multiplication
and summation. We include these two and generalizes them int o the “Multi-Linear” function in
Eq. (1)(b), i.e., it is linear with respect to each of its arguments individua lly. Here we use the set
notation {·}j∈M(i) to denote the set of arguments. What is special in a “Multi-Li near” function is
that, when taking its expectation with respect to a certain argument, the expectation goes in and the
resulted formula simply replaces the argument by its mean. This is useful for mean-ﬁeld algorithms.
After dendritic summation, relationship between the summe d current and the ﬁring rate can be
well researched by in-vitro neuroscience experiments where input current can be extern ally con-
trolled. Various frequency-current function can also be derived from popular spike generation mod-
els, such as Leaky Integrate-and-Fire model (see [12] for an overview), Exponential Integrate-And-
Fire Model [13], Spike Response Model [23], or more realisti c models such as Hodgkin-Huxley
Model [24]. In particular, [14][15][16] formulated the LNP simpliﬁcation through a nonlinear map-
ping followed by a Poisson point process random sampling. Th e analytical study with simulation
in [17] show that the simpliﬁcation is effective. We put this in Eq. (1)(c).
2
−10 −5 0 5 10 15 200
0.2
0.4
0.6
0.8
1
Input Current I0
Output Frequency
Figure 1: Frequency-current curve
of leaky integrate-and-ﬁre neurons
for a particular set of parameters.
More from different neuron models
can be found in [12].
Since neural spikes has refractory period with at least 1ms,ﬁr-
ing rate will not exceed 1000Hz and there won’t be more than
one spike happening in 1ms, so λi (τ) ∈[0, 1]. In practice,
ﬁring rate is often even lower. In addition, discrete time steps
in the resolution of 1ms per each step seems quite ﬁne-grained
for most perceptual tasks (e.g., visual recognition emerges in
75-80ms [25]). For these reasons, discretizing time into 1m s
per each step may be reasonable. By LeCam’s theorem [26],
the discrete counterpart of nonhomogeneous Poisson point
process is the nonhomogeneous Bernoulli process, and the
Bernoulli probability at each discrete time step is approxi -
mated by the integration of rate function at that time period .
We use t ∈Z to denote the discrete time steps, whose unit
is millisecond. We use φ(t)
i ∈[0, 1] to denote the Bernoulli
probability and x(t)
i ∈{0, 1}to denote the spike indicator at
time t,
x(t)
i ∼Bernoulli
(
φ(t)
i
)
, φ (t)
i ≜
∫ t+1
t
λi (τ) dτ (2)
The corresponding discrete form of Eq. (1)(a,b,c) is by substituting Xi (τ) with x(t)
i and substituting
α, D functions with their discrete form similar to the deﬁnition of φ(t)
i in Eq. (2). There are two
ways to proceed with the formulation. In case D function is the Dirac function located at 0, or is
close to that, the convolution with D can be ignored, the discrete counterpart of Eq. (1)(a,b,c) are
I(t)
ij = Wij ·
∑ Kα
k=1
α (k) ·x(t−k)
i
φ(t)
i = σ
(
Multi-Linear
( {
I(t)
ij
}
j∈M(i)
)) (3)
Or if D is non-trivial, but the Multi-Linear function in Eq. (1)(b) reduces to a linear summation, by
the associativity of convolution, we can let ǫ = α ∗D, and the discrete counterpart of Eq. (1)(a,b,c)
reduces to the following, assuming ﬁnite time course of the function ǫ.
φ(t)
i = σ
( ∑
j∈M(i)
Wij ·
∑ K
k=1
ǫ (k) ·x(t−k)
i
)
(4)
In [17], D functions resulted from spiking neuron models are not Diracfunctions, but are very close.
In the latter part of this paper, we will focus on the case whenlinear summation is true (hence Eq. (4)
holds) although multi-linearity in Eq. (3) may generalize o ur formulation to high-order Boltzmann
machines. Another observation is that in Eq. (4), if D is exponential function (true in [17]) and α is
exponential function too (true in [12]), and assuming both functions have the same scale parameter,
ǫ becomes the classic α-function [27] used in postsynaptic modeling.
Another issue is whether linear summation in Eq. (4) or multi -linear integration in Eq. (3) contains
a constant bias term bi. In [28], there is such a bias term in the spiking response mod el. This will
make our latter part easier, but in our derivation later, we f ocus on the case when there is no such
bias. It will be trivial then to turn to the case where bias exists.
To sum up, the LNP neuron model we are going to deal with in the l atter part is as the following,
assuming His the index set of neurons whose activity is not bound to observations (external stimuli).
For all i ∈H, t ∈Z+,



x(t)
i ∼Bernoulli
(
φ(t)
i
)
φ(t)
i = σ
( ∑
j∈M(i) Wij ·∑ K
k=1 ǫ (k) ·x(t−k)
i
) (5)
σ function can be ﬁt from data or from realistic neuron models, and has some relationship with
the frequency-current curve. According to [12][17], such n onlinearity in different spiking neuron
models may be non-decreasing, close to zero towards the left , increase almost linearly in a certain
interval. Also, since ﬁring rate has upper bounds, σ won’t increase unboundedly. For this reason, σ
may be considered as a modiﬁcation of sigmoid function to meet certain demands of inference. An
example of frequency-current curve is shown in Fig. 1. See also [17] for another type of nonlinearity.
3
In Fig. 1, parameters of the leaky integrate-and-ﬁre model are set as τm = 20 ms, ∆ abs = 1 ms, R =
1 mΩ , υ = 1 mV (with notations following [12]).
In the biological brain, all neurons are computing in parall el. Thus we want to investigate on the
case when Eq. (5) is executed in parallel, vectors x(t) = ( x(t)
i )n
i=1 and φ(t) = ( φ(t)
i )n
i=1 over time
forms a Markov chain of orderK. The latter part focuses on showing its stochastic inference nature.
3 Semi-Stochastic Inference on Boltzmann Machines
In this section, ﬁrst we compare Gibbs sampling and variational inference on Boltzmann machines,
emphasizing their architectural similarity, then we prese nt the semi-stochastic inference algorithm
by combining them. Some notations we use in this section will overlap with the last section. This is
intended since we want to link the quantities in both sections.
Let Y ∈{A, B}n be a collection of random variables. Each of them has two possible values A and
B (e.g., A = 0 and B = 1 in the binary case). Let V, Hbe a partition of the index setI= {1, ..., n}.
We denote YV as the visible variables, YH as the hidden variables. We use lower case ofY to denote
its values, e.g., yi, yH or y. In this paper we use the Boltzmann machines with softmax units [29],
p (y) = 1
Z exp

 1
2 ·
∑
i,j:i̸=j;u,v∈{A,B}
[ [yi = u] ]·[ [yi = v] ]·Viu,jv −
∑
i∈I,u∈{A,B}
[ [yi = u] ]·ciu

 (6)
Here we denote [ [·] ]as the binary indicator. [ [·] ] = 1 if the statement is true and [ [·] ] = 0 otherwise.
V = ( Viu,jv)i,u,j,vis a four-dimensional tensor of size 4n2, c = ( ciu)i,u is a matrix of size 2n. Z
is the normalization constant (partition function), which depends on V and c. Note that this family
has no more capacity than the original Boltzmann machine [30] and the Ising model [31].
In Gibbs sampling without parallelism [32][33], we are give n an observed value yV , and we want
to ﬁnd a sequence
{
y(t)
H
} ∞
t=0
which converges in distribution to the true posterior p (YH|YV = yV ).
The algorithm proceeds as follows. First we initialize y(0)
H , and let y(0)
V take on the observed values.
Then at each iteration t ∈Z+ we pick an index i ∈H and do the following,
(a) Calculate proposal φ(t)
iu = σ
( ∑
j∈M(i),v∈{A,B}
Wiu,jv ·
[ [
y(t−1)
j = v
] ]
−biu
)
, ∀u ∈{A, B}
(b) Sample y(t)
i ∼Bernoulli
(
φ(t)
iA , φ(t)
iB
)
, pass down samples by y(t)
∼i = y(t−1)
∼i (7)
Here we denoteWiu,jv = Viu,jv−Vi¯u,jv and biu = ciu −ci¯u with ¯u deﬁned as {¯u}= {A, B}\{u}.
σ (x) = 1 / (1 + exp ( −x)) is the sigmoid function. M (i) is the index set of Markov blanket of Yi,
i.e., M (i) = {j ∈I|j ̸= i, ∃u, v ∈{A, B}, s.t. Viu,jv ̸= 0 }.
Variational inference with factorized approximation fami ly often takes the form of a mean-ﬁeld
version of Gibbs sampling ( e.g., [34]). We adopt the same family of variational distributio ns as
in [35], such that q (YH|Θ H) = ∏
i∈H q (Yi|θiA, θiB ), where each q (Yi|θiA, θiB ) is a Bernoulli
distribution such that θiu = q (Yi = u|θiA, θiB ) for u ∈{A, B}. 1 The objective of variational
inference is to ﬁnd Θ H by minimizing the loss function L (Θ H) = KL (q (Y|Θ ) ||p (Y)), here
KL (q||p) = Eq log (q/p) is the KL-divergence.
The algorithm proceeds as: ﬁrst we randomly initialize Θ (0)
H , then at each iteration t ∈Z+ we pick
an index i ∈H and do the following,
(a) Calculate update φ(t)
iu = σ
( ∑
j∈M(i),v∈{A,B}
Wiu,jv ·θ(t−1)
jv −biu
)
, ∀u ∈{A, B}
(b) Assign θ(t)
iu = φ(t)
iu , pass down θ(t)
∼iu = θ(t−1)
∼iu , ∀u ∈{A, B}. (8)
Our motivation for combining them is as follows. First, expo nential family distributions have sufﬁ-
cient statistics [36]. By accumulating the empirical expectation of sufﬁcient statistics, we can online
1For simplicity, when i ∈ V we write θiu and q (Yi|θiA, θiB ) as well, except that θiu for all iterations are
ﬁxed as the observed value. We denote Θ = (θiu)i∈I,u∈{A,B} and write q (Y|Θ ) =∏
i∈I q (Yi|θiA, θiB ).
4
estimate the parameter if data instances are completely observed, e.g., to online estimate a Gaussian
distribution N(X|µ, σ), we can incrementally calculate ∑ t
i=1 xi and ∑ t
i=1 x2
i , and estimate ˆµ and
ˆσ at each time t by only using these two statistics.
Now if we take the examples in Gibbs sampling as the online data instances and estimate a distribu-
tion from within the variational approximation family, we will end up getting a marginal probability
distribution for each random variable (when mixing is good). However, we can also change the way
how the sequence is generated, such that the online estimated variational distribution biases towards
whatever we want. In particular, we can take the calculated distribution in Eq. (8)(a) as the proposal
distribution, used to sample a new data example, and increme ntally update the variational distribu-
tion by using that. The updated variational distribution serve as the next input to Eq. (8). If decaying
as in [37] is used when accumulating sufﬁcient statistics, the variational distribution estimated in our
Bernoulli case, is simply a weighted sum of the most recent examples:
θ(t)
iu =
∑ −∞
k=0
ǫ (k) ·
[ [
y(t−k)
i = u
] ]
, for u ∈{A, B} (9)
The −∞in Eq. (9) simply means the starting time of accumulation. ǫ function is non-negative and
should sum up to one for all k in the effective time span. In case constant decaying ratio i s used on
the last accumulated sufﬁcient statistics as in [37], the ǫ function above decays exponentially. Some
choices of ǫ function may not have any online updating scheme that produces it. In practice, we can
either do online updating or let ǫ have a ﬁnite time course, up to some K ∈Z+.
In semi-stochastic inference, we choose a non-negative weight function ǫ (k) which sum up to 1 for
k ∈{1, 2, ..., K }. We initialize both y(0) and Θ (0). In each iteration t ∈Z+ we pick an index
i ∈H and do the following,
(a) φ(t)
iu = σ
( ∑
j∈M(i),v∈{A,B}
Wiu,jv ·θ(t−1)
jv −biu
)
, for u ∈{A, B}
(b) y(t)
i ∼Bernoulli
(
φ(t)
iA , φ(t)
iB
)
(10)
(c) θ(t)
iu =
∑ K
k=1
ǫ (k) ·
[ [
y(t−k+1)
i = u
] ]
, for u ∈{A, B}
In sequential inference, for every other random variable wh ich does not carry out these updating
steps, we will pass down θ and most recent K copies of y’s to the next time step. In step (c),
additional normalization is needed when t < K . Details are ignored for clarity of presentation.
In another way, we can also look at algorithm in Eq. (10) as a “r andom slowing-down” / stochastic
approximation / momentum version of variational inference . When an updated variational distri-
bution is computed, we don’t immediately turn to the updated distribution. Instead, we use it to
sample a new example, and use it to incrementally update the o riginal variational distribution. The
expectation of the new example will be the same as the updated variational distribution.
4 Making Detailed Correspondence
To connect LNP neuron model in Eq. (5) with semi-stochastic i nference algorithm in Eq. (10),
several issues needs to be resolved:
1. LNP neurons don’t have the bias terms (such as biu), although their nonlinear transform
function may have constant shift identical to all neurons.
2. In Eq. (10), two Linear-Nonlinear procedures in (a)(c) co rresponds to only one sampling
step in (b). This is an architectural difference.
3. Biological neurons can only be excitatory or inhibitory.
4. Biological neurons carry out the updates all in parallel.
To resolve issue 1, we notice that when the weights V and biases c in Eq. (6) are turned into W
and b in Eq. (10), for ﬁxed i, j ∈I, ﬁxed u, v ∈{A, B}, and any constant C ∈R the following
5
operation is invariant to the inference algorithm in Eq. (10) since θ(t)
iA and θ(t)
iB sum up to one.





bnew
iu ←bold
iu + C
W new
iu,jv ←W old
iu,jv + C
W new
iu,j¯v ←W old
iu,j¯v + C
(11)
For this reason, we can always modify the algorithm in Eq. (10 ) to make any bias we want. In our
experiments, when we want to remove a bias biu, we subtract biu/|M (i) |from both Wiu,jv and
Wiu,j¯v for each j ∈M (i).
To resolve issue 2, we propose the “Event-Network”. We ﬁrst s plit every sampling step (b) in
Eq. (10) into two, namely[ [y(t)
i = A] ]∼Bernoulli(φ(t)
iA ) and [ [y(t)
i = B] ]∼Bernoulli(φ(t)
iB ). Now we
can allocate two neurons to take care of the two parts in all three steps (a)(b)(c). Units of the resulted
network no longer corresponds to random variables, they are now probabilistic events. When we do
this, there is no guarantee thatθ(t)
iA and θ(t)
iB will sum up to 1, but this will be approximately true since
the weights Wiu,jv’s are compatible with each other. Whether this results in me aningful inference
algorithms needs further investigation. We show veriﬁcation in our experiment part.
To resolve issue 3, once we have a network containing neurons with both excitatory and inhibitory
outgoing synapses, we can duplicate each neuron into two, bo th have the same incoming synapses
with the original neural efﬁcacy. Then one of them take all positive outgoing synapses before split-
ting, and the other take all negative ones. Both of this and th e last modiﬁcation of the network are
exactly invariant for variational inference but may only be approximately true for semi-stochastic
inference. We show experimental veriﬁcation for this modiﬁcation as well.
If after resolving all these issues we denote x(t)
iu = [ [y(t)
i = u] ]for u ∈{A, B}, and re-index every-
thing, the resulted algorithm will be exactly the same as Eq. (5) if executed in parallel. Parallelism
is not only an issue in semi-stochastic inference. It is noti ced in [33] that when adjacent random
variables of a Boltzmann machine are updated in parallel in G ibbs sampling, the sample sequence
does not converge to the right distribution. This issue is al so raised very often recently ( e.g. [38]).
Variational inference is a ﬁxed point iteration. So it may co nverge to the same answer when up-
dated all in parallel, but we do observed experimentally that the variational parameters converges to
the ﬁxed point after severe oscillation. Whether the slowin g down version alleviates this problem
requires further investigation2.
5 Experimental Veriﬁcation
Figure 2: Samples of reconstruction results (selected randomly). Rows from top to bottom: original
images; VarO; SemiO; SemiEN; SemiB; SemiU.
In this section, we focus on illustrating that the semi-stoc hastic inference, as well as all the modi-
ﬁcations in section 4, are valid Bayesian inference algorit hms, which behave similar to variational
inference. In all experiments, we set ǫ(k) as exp(−k/2) normalized to have summation 1 over its
domain k ∈{1, ..., 30}. All inferences are carried out in synchronized full parallelism.
To avoid the complication of learning deep Boltzmann machin es, we take the learned Boltzmann
machine from code in [35]. Since we focus on inference on Boltzmann machine, we take the learned
model before back-propagation. The Boltzmann machine we us e consists of three layers, the ﬁrst
2We observe that different choices of ǫ functions result in different smoothness of the trajectori es, but it’s
not clear how. We leave this as future work.
6
0 20 40 60 80 100
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
# iterations
activation values θi
(t)
0 20 40 60 80 100
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
# iterations
activation values θi
(t)
0 20 40 60 80 100
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
# iterations
activation values θi
(t)
0 20 40 60 80 100
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
# iterations
activation values θi
(t)
Figure 3: Sample trajectory of θ(t)
i in VarO and SemiO. In all above, Blue line is for VarO, Green
line is for SemiO, and Red line is the moving average of Green line over window size 30.
layer for input image consists of 784 units, the second consists of 500 units and the third consists of
1000 units. In the following, we refer to θ(t)
i in both algorithms as the activation of neuron i.
To see the effect of modiﬁcations in section 4, we show result s of semi-stochastic inference after
each of them. Note that variational inference is exactly inv ariant on all these modiﬁcations, hence
there no need to verify them. We use the following abbreviations in this part:
•VarO: variational inference applied on the original deep Boltzmann machine.
•SemiO: semi-stochastic inference applied on the original deep Boltzmann machine.
•SemiEN: Same as SemiO except each sampling is duplicated and put on two neurons.
•SemiB: Same as SemiEN except biases are removed by transforming the network.
•SemiU: Same as SemiB except each neuron is duplicated into two, such that no neuro n
contains out-going synapses with different signs.
The experiment we did is similar to the reconstruction exper iment in [7]. We do the following: (i)
plug in an image on the input layer, (ii) apply the inference a lgorithm, (iii) frozen the activations
on the topmost layer and round up to 0 or 1, (iv) set the topmost layer as observed and other layers
as hidden, (v) infer back, (vi) read off the image as activati ons on the input layer. Both (ii) and (v)
are randomly initialized and have 100 iterations. For semi- stochastic inference, in steps (iii) and
(vi) the activations are taken as averages over the last 50 it erations, while for variational inference
we simply use the last activation of each neuron. This is because semi-stochastic inference has ran-
dom convergence, activations approach some ﬁxed point withrandom perturbation (see justiﬁcation
later). Some reconstruction examples are shown in Fig. 2. Ea ch algorithm has good and bad cases,
also note that because of the randomness, even for identical input image and identical initialization,
semi-stochastic inference may get different results in different trials.
In ﬁgure 3, we showed sample trajectory of activations in step (ii) of VarO and SemiO for the same
variables. Each trajectory are those θ(t)
i ’s of one neuron over all iterations. We see that activations
in variational inference converge to a stable solution with oscillation before convergence, while
semi-stochastic inference has their moving average conver ges to some similar values with random
perturbation. Figure 5(a) is a scatter plot of mean of trajec tories vs. the corresponding converged
activation in variational inference, extracted from 2000 trajectories randomly chosen from 50 trials.
We also observed that when the mean is close to extreme values (such as 0 and 1), the standard
variance will be smaller, as can be found in ﬁgure 3 and more cl early in ﬁgure 5(b). This means
that semi-stochastic inference converges randomly to the v ariational inference solution, with more
conﬁdent variables having less random perturbation.
Finally, we investigated on how well the splitting in sectio n 4 preserves identity relationships that
are expected on SemiU. In particular, if two neurons are resulted from splitting o ne sampling into
two in resolving issue 2, their summed activation should be c lose to 1. If two neurons are resulted
from splitting one neuron into two in resolving issue 3, thei r activation should be equal. Figure 4
illustrated three samples and histogram on 2000 random trajectories are given in ﬁgure 5(c)(d).
7
0 20 40 60 80 100 120
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
# iterations
activation values
0 20 40 60 80 100 120
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
# iterations
activation values
0 20 40 60 80 100 120
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
# iterations
activation values
(a) (b) (c)
Figure 4: Sample trajectory from SemiU. Each plot shows four neurons resulted from two splits
on a single neuron. Thick dash lines: neurons with positive o utgoing weights. Thin solid lines:
corresponding neurons with negative outgoing weights. Blue and Red lines are for the two neurons
from splitting one sampling into two. Green line is the average of Blue and Red. Ideally, lines with
same color should be equal, and green lines should be at 0.5.
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
converged activation values in VarO
mean activation values in SemiO
0 0.2 0.4 0.6 0.8 1
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
mean activation values in SemiO
std of activation values in SemiO
0 0.2 0.4 0.6 0.80
200
400
600
800
1000
0 0.5 10
200
400
600
800
1000
(a) (b) (c) (d)
Figure 5: Some statistics of the algorithms. (a) means of tra jectories in SemiO vs. converged
activation values in VarO for the same hidden variables. (b) standard deviation vs. me an, from the
trajectories in SemiO. (c) histogram of standard deviation(red line + blue line - 1.0). (d) histogram
of standard deviation(dash line - solid line). (Here red, blue, dash, solid lines refer to ﬁgure 4.)
6 Conclusions and Discussions
We pointed out the stochastic inference nature of the LNP neu ron models, which may serve as an
interpretation of neural coding. The stochastic convergen ce seems reasonable, e.g., when we see
ambiguous images, our perception will switch back and forth between plausible explanations.
There are many behavioral experiments showing statistical optimality of perception and learning
(e.g. [39]) and many calls for neuronal modeling that achieves this optimality (e.g. [40]). Arguably,
the mode-seeking nature of variational inference [41] make it natural for interpreting perception.
Since Boltzmann machines with hidden variables are universal approximators (e.g. [42]), modeling
knowledge representation by that is safe if the world can be a pproximately binarized. Yet we only
show that with particular choices of weights, neurons are ca pable of representing a Boltzmann ma-
chine and carry out inference. Then the question is what does this do if the network is a recurrent
neural network in general. Real biological neural networks may be learned by taking the Boltz-
mann machine representation and undergoing discriminativ e training or reinforcement learning to
optimize performance of inference directly (especially plausible when time dimension is involved),
which yields a model not necessarily convey consistent prob abilistic semantics but may do so in
the limit. The particular form of nonlinearity in biologica l neurons may also result from such opti-
mization. Another issue is the low average activity in biolo gical neurons. It may be a certain form
of sparse coding [4] on top of inference so that optimizing su ch inference procedure will not yield
degenerate solution. In this case, likelihood-based learn ing such as contrastive divergence may still
be applicable as a further reﬁnement.
Acknowledgement
I want to thank Geoffrey Hinton, Mikhail Belkin, DeLiang Wang, Brian Kulis for helpful discussion
and The Ohio Supercomputer Center for providing computation resource.
8
References
[1] Frank Rosenblatt. The perceptron: a probabilistic mode l for information storage and organiza-
tion in the brain. pages 89–114, 1988.
[2] David E Rumelhart, James L Mcclelland, and Geoffrey E Hin ton. Parallel Distributed Pro-
cessing: Explorations in the Microstructure of Cognition. MIT Press, August 1986.
[3] Honglak Lee, Chaitanya Ekanadham, and Andrew Y . Ng. Spar se deep belief net model for
visual area V2. In Advances in Neural Information Processing Systems 20 , pages 873–880.
2008.
[4] Bruno. A. Olshausen and David J. Field. Emergence of simp le-cell receptive ﬁeld properties
by learning a sparse code for natural images. Nature, 381:607–609, 1996.
[5] Honglak Lee, Alexis Battle, Rajat Raina, and Andrew Y . Ng . Efﬁcient sparse coding algo-
rithms. In Advances in Neural Information Processing Systems 19, pages 801–808. 2007.
[6] Geoffrey E. Hinton, Simon Osindero, and Yee-Whye Teh. A f ast learning algorithm for deep
belief nets. Neural Comput., 18(7):1527–1554, July 2006.
[7] Geoffrey E. Hinton and Ruslan R. Salakhutdinov. Reducin g the dimensionality of data with
neural networks. Science, 313(5786):504–507, 2006.
[8] Yoshua Bengio, Pascal Lamblin, Dan Popovici, and Hugo Larochelle. Greedy layer-wise train-
ing of deep networks. In B. Sch¨ olkopf, J. Platt, and T. Hoffman, editors, Advances in Neural
Information Processing Systems 19, pages 153–160. MIT Press, Cambridge, MA, 2007.
[9] Quoc Le, Marc’Aurelio Ranzato, Rajat Monga, Matthieu Devin, Kai Chen, Greg Corrado, Jeff
Dean, and Andrew Ng. Building high-level features using lar ge scale unsupervised learning.
In Proceedings of the Twenty-Ninth International Conferenceon Machine Learning, 2012.
[10] Honglak Lee, Roger Grosse, Rajesh Ranganath, and Andre w Y . Ng. Convolutional deep belief
networks for scalable unsupervised learning of hierarchical representations. In Proceedings of
the 26th Annual International Conference on Machine Learni ng, ICML ’09, pages 609–616,
New York, NY , USA, 2009. ACM.
[11] Andreas V . M. Herz, Tim Gollisch, Christian K. Machens,and Dieter Jaeger. Modelling single-
neuron dynamics and computations: A balance of detail and ab straction. Science, 314:80–84,
October 2006.
[12] Wulfram Gerstner and Werner M. Kistler. Spiking Neuron Models: Single Neurons, Popula-
tions, Plasticity. Cambridge University Press, 1 edition, 2002.
[13] Nicolas Fourcaud-Trocm, David Hansel, Carl van Vreesw ijk, and Nicolas Brunel. How spike
generation mechanisms determine the neuronal response to ﬂuctuating inputs. The Journal of
Neuroscience, 23(37):11628–11640, 2003.
[14] E. J. Chichilnisky. A simple white noise analysis of neuronal light responses. Network (Bristol,
England), 12(2):199–213, May 2001.
[15] Eero P. Simoncelli, Liam Paninski, Jonathan Pillow, and Odelia Schwartz. Characterization of
neural responses with stochastic stimuli. The cognitive neurosciences, 2004.
[16] Liam Paninski. Maximum likelihood estimation of casca de point-process neural encoding
models. Network: Comput. Neural Syst., 15(04):243–262, November 2004.
[17] Srdjan Ostojic and Nicolas Brunel. From spiking neuron models to Linear-Nonlinear models.
PLoS Comput Biol, 7(1):e1001056+, January 2011.
[18] David J. Tolhurst, J. A. Movshon, and A. F. Dean. The statistical reliability of signals in single
neurons in cat and monkey visual cortex. Vision Research, 23(8):775–785, 1983.
[19] Wei J. Ma, Jeffrey M. Beck, Peter E. Latham, and Alexandr e Pouget. Bayesian inference with
probabilistic population codes. Nature Neuroscience, 9(11):1432–1438, October 2006.
[20] Idan Segev. What do dendrites and their synapses tell the neuron? Journal of Neurophysiology,
95(3):1295–1297, March 2006.
[21] Idan Segev and Michael London. Dendritic processing. I n Michael A. Arbib, editor, The
handbook of brain theory and neural networks , pages 324–332. MIT Press, Cambridge, MA,
USA, 2003.
9
[22] Roberto Araya, Kenneth B. Eisenthal, and Rafael Yuste. Dendritic spines linearize the summa-
tion of excitatory potentials. Proceedings of the National Academy of Sciences of the Unite d
States of America, 103(49):pp. 18799–18804, 2006.
[23] W. Gerstner. Spike-response model. Scholarpedia, 3(12):1343, 2008.
[24] Alan L. Hodgkin and Andrew F. Huxley. A quantitative des cription of membrane current and
its application to conduction and excitation in nerve. The Journal of Physiology, 117(4):500–
544, 1952.
[25] Ruﬁn Vanrullen and Simon J. Thorpe. The time course of vi sual processing: From early
perception to decision-making. J. Cognitive Neuroscience, 13(4):454–461, May 2001.
[26] Lucien le Cam. An approximation theorem for the poisson binomial distribution. Paciﬁc
Journal of Mathematics, 10:1181–1197, 1960.
[27] Wilfrid Rall. Distinguishing theoretical synaptic po tentials computed for different soma-
dendritic distributions of synaptic input. Journal of neurophysiology , 30(5):1138–1168,
September 1967.
[28] Danilo J. Rezende, Daan Wierstra, and Wulfram Gerstner . Variational learning for recurrent
spiking networks. In J. Shawe-Taylor, R.S. Zemel, P. Bartlett, F.C.N. Pereira, and K.Q. Wein-
berger, editors, Advances in Neural Information Processing Systems 24, pages 136–144. 2011.
[29] Ruslan Salakhutdinov and Geoffrey Hinton. Replicated softmax: an undirected topic model.
In Y . Bengio, D. Schuurmans, J. Lafferty, C. K. I. Williams, and A. Culotta, editors, Advances
in Neural Information Processing Systems 22, pages 1607–1614. 2009.
[30] David H. Ackley, Geoffrey E. Hinton, and Terrence J. Sej nowski. A Learning Algorithm for
Boltzmann Machines. Cognitive Science 9, pages 147–169, 1985.
[31] Ernst Ising. Beitrag zur theorie des ferromagnetismus . Zeitschrift fr Physik A Hadrons and
Nuclei, 31:253–258, 1925. 10.1007/BF02980577.
[32] Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian
restoration of images. IEEE Transactions on Pattern Analysis and Machine Intellig ence,
6:721–741, 1984.
[33] Emile Aarts and Jan Korst. Simulated annealing and Boltzmann machines: a stochastic a p-
proach to combinatorial optimization and neural computing . John Wiley & Sons, Inc., New
York, NY , USA, 1989.
[34] Yee Whye Teh, David Newman, and Max Welling. A collapsed variational bayesian inference
algorithm for latent dirichlet allocation. In B. Sch¨ olkop f, J. Platt, and T. Hoffman, editors,
Advances in Neural Information Processing Systems 19 , pages 1353–1360. MIT Press, Cam-
bridge, MA, 2007.
[35] Ruslan Salakhutdinov and Geoffrey Hinton. Deep boltzmann machines. Artiﬁcial Intelligence,
5(2):448C455, 2009.
[36] Erich Leo Lehmann. Theory of point estimation / E.L. Lehmann . John Wiley, New York,
Singapore, 1983.
[37] Radford M. Neal and Geoffrey E. Hinton. A view of the EM algorithm that justiﬁes incremen-
tal, sparse, and other variants. In Learning in graphical models, pages 355–368. 1999.
[38] Arthur Gretton Carlos Guestrin Joseph E. Gonzalez, Yuc heng Low. Parallel gibbs sampling:
From colored ﬁelds to thin junction trees. In Proceedings of the 14th International Conference
on Artiﬁcial Intelligence and Statistics (AISTATS) 2011, pages 324–332, 2011.
[39] Tianming Yang and Michael N. Shadlen. Probabilistic re asoning by neurons. Nature, 2007.
[40] J´ ozsef Fiser, Pietro Berkes, Gerg˝ o Orb´ an, and M´ at´e Lengyel. Statistically optimal percep-
tion and learning: from behavior to neural representations . Trends in Cognitive Sciences ,
14(3):119–130, March 2010.
[41] Thomas Minka. Divergence measures and message passing . Technical report, Microsoft Re-
search, 2005.
[42] Nicolas Le Roux and Yoshua Bengio. Representational po wer of restricted boltzmann ma-
chines and deep belief networks. Neural Comput., 20(6):1631–1649, June 2008.
10
arXiv:1210.8442v3  [cs.AI]  27 Jan 2013
Supplemental Material
Several Justiﬁcations and Open Questions
In this supplemental material, we deal with two issues that are not made clear in the main text: (1) we
give the details of discrete time approximation, which even tually leads to the inference algorithm;
(2) we present preliminary results on stochastic stability of the neuron model. After that, we also
point out several open questions related to the Linear-Nonl inear-Poisson (LNP) neuron model and
its interpretation of Bayesian inference, with possible directions for solving these problems.
For now, we focus on the neuron model with a Poisson Point Process (PPP) when there is no linear
convolution with function D in the main text and the dendritic integration is linear, as the following
λi (τ) = σ
(∑
j∈M(i) Wij · (Xj ∗ α) (τ) + ei
)
,
Xi (τ) ∼ PoissonPointProcess (λi (τ)) , ∀i ∈ I
(1)
Here τ ∈ R is the continuous time index, which has unit millisecond. An external input term ei
is introduced so that we don’t need to explicitly distinguis h visible and hidden units. Every ei is
constant over time. For example, ei = 0 means the neuron corresponds to a hidden unit, while a
very positive and a very negativeei correspond to a visible unit with observation 1 and 0 respectively.
We also assume that τ = 0 is the beginning time and each λi (0) is chosen by the user
1 The Stochastic Differential Equation Formulation
In this section, we reformulate Eq. (1) as a stochastic diffe rential equation, which is for the ease of
the later parts.
If we choose the parametrization of α function as the exponential function α (τ) = a · exp (−a · τ),
and deﬁne Yi (τ) ≜ (Xi ∗ α) (τ), then Yi is the solution of the differential equation as the following.
dYi
dτ = −a · Yi + a · Xi (2)
Thus we can rewrite Eq. (1) in terms of Y = ( Yi)i∈I , and we get a system of coupled stochastic
differential equations (SDE) as the following
dY
dτ = −a · (Y − λ) + a · (X − λ) , where
∀i ∈ I ,
{
λi (τ) = σ
(∑
j∈M(i) Wij · Yj (τ) + ei
)
Xi (τ) ∼ PoissonPointProcess (λi)
(3)
The right hand side of Eq. (3) contains two terms, the determi nistic term −a · (Y − λ) and the
stochastic term a · (X − λ). The deterministic term reaches zero when Y = λ, and thus is the
equilibrium of the (recurrent) neural network, if there is one. The stochastic part, roughly speaking,
has short time integral close to zero, which is referred to as random perturbation. The deterministic
counterpart of Eq. (3) is deﬁned as the following.
dY
dτ = −a · (Y − λ) (4)
In the literature of stochastic differential equations, th e stochastic stability of Eq. (3) is deﬁned
in many sense analogous to the Lyapunov stability of Eq. (4). For example, in [1] the stability
1
under small random perturbation is introduced, which is deﬁned as a capability of the deterministic
system, in the sense that as long as the random force can be mad e sufﬁciently small, the combined
stochastic system will exhibit certain stability, although the requirement imposed upon the strength
of the random force may not satisfy It is likely that Eq. (4) sa tisﬁes this stability property. We
haven’t obtained a rigorous proof yet, but in the following p art, we give an analysis of this kind of
stochastic stability in the discrete time case. We ﬁrst give the detail about how the discrete case is
achieved and then proceed by analyzing it.
2 Details of Discrete Time Approximation
In this section, we give the details about how the PPP formula tion is discretized into the Bernoulli
process formulation, and how it is related to the ﬁxed-point iteration procedure, which are not made
clear in the main text.
Let t ∈ Z be the discrete time index such that discrete timet corresponds to continuous time τ = tǫ,
where ǫ ∈ R+ is the length of a small time step. The discrete time approxim ation of PPP in Eq. (3)
is as the following.
x(t)
i ∼ Bernoulli
(
ǫ · λ(t)
i
)
, where λ(t)
i = 1
ǫ
∫
[tǫ,tǫ+ǫ)
λi (τ) dτ (5)
By LeCam’s theorem [2], it can be shown that when ǫ → 0, for any ﬁxed time interval [τ1, τ2),
the sum of x(t)
i fall in that interval converges in distribution to a Poisson distribution with rate∫ τ2
τ1
λi (τ) dτ. That means the non-homogeneous Bernoulli process approac hes non-homogeneous
PPP We deﬁne y(t)
i as the following (in which we ignored the notation issue rega rding the starting
point)
y(t)
i =
∑ ∞
k=0
α(k)x(t−1−k)
i , where α(k) = a · (1 − aǫ)k (6)
With this deﬁnition, we have the following approximation formula
α(τ /ǫ) = a · (1 − aǫ)τ /ǫ → a · exp (−aτ) = α (τ) , as ǫ → 0 (7)
For this reason, y(t)
i deﬁned in Eq. (6) is an approximation of Yi (tǫ) when the Bernoulli process is
an approximation of PPP. Eq. (6) can also be rewritten recursively as the following.
y(t+1)
i = (1 − aǫ) · y(t)
i + a · x(t)
i (8)
In addition, for τ > 0, each Yi (τ) is right continuous with discontinuity located only at the spiking
points indicated by Xi (τ), and so is λi (τ). Thus when ǫ → 0, λ(⌈τ /ǫ⌉)
i → λi (τ), where ⌈τ/ǫ⌉ is
the smallest integer that is larger thanτ/ǫ. For this reason, λ(t)
i is an approximation of λi (tǫ). Since
λi (tǫ) is a function of Yi (tǫ) and y(t)
i approximates Yi (tǫ), λ(t)
i is approximately a function of y(t)
i .
As a conclusion, the discrete time approximation of Eq. (3) i s the following, where the ﬁrst line
is derived from Eq. (8). The approximation gets better (clos er to the continuous time counterpart)
when ǫ → 0
y(t+1)−y(t)
ǫ = −a ·
(
y(t) − λ(t))
+ a ·
(
x(t)/ǫ − λ(t))
, where
λ(t)
i = σ
(∑
j∈M(i) Wij · y(t)
j + e(t)
i
)
, x (t)
i ∼ Bernoulli
(
ǫ · λ(t)
i
) (9)
We can also replace the ﬁrst formula in Eq. (9) by Eq. (6). The r esulted formula matches the
neuron model in our main text. Again we can deﬁne the determin istic counterpart of Eq. (9) as the
following.
y(t+1) − y(t)
ǫ = −a ·
(
y(t) − λ(t)
)
(10)
Note that when aǫ = 1 , Eq. (10) is exactly the forward evaluation step of (recurre nt) neural net-
works. When aǫ < 1, it is the momentum version of the forward step. In case the ne ural network
represents the variational inference algorithm over Boltz mann machines, i.e., weights are symmet-
ric, the forward step is guaranteed to converge to a local optimum, no matter where the starting point
is.
2
3 Stochastic Stability of the Neuron Model
In this section, we illustrate the stochastic stability of E q. (9). Eventually, the stability of the
stochastic system will rely on the stability of the deterministic system as in Eq. (10). In general it is
not clear when the dynamical system, either discrete time in Eq. (10) or continuous time in Eq. (4)
has stable ﬁxed-point. Several known results are listed, fo r example, in [3]. An important case, as
what our main text was depending on is that, when the neural we ights are symmetric, the network
converges no matter where it starts. In case convergence is guaranteed with only one ﬁxed-point, the
corresponding dynamical system can be regarded as dissipat ive with its ﬁxed-point being globally
asymptotically stable. In this section, we denote y∗ as the position in the unit hypercube [0, 1]n that
satisﬁes the following.
y∗ = σ (Wy∗ + e) (11)
And unless speciﬁed otherwise, it is assumed that there is on ly one such ﬁxed-point. We further
assume that there is a Lyapunov function V of the continuous system, such that for y ∈ [0, 1]n, we
have
V (y) = 0 , iff y = y∗
V (y) > 0, if y ∈ [0, 1]n \ {y∗}⣨
∂V
∂y , −y + λ
⟩
⩽ 0, with eq. only aty∗
(12)
Here λ is deﬁned as in Eq. (9). We focus on the case when ǫ = 1 in which case Eq. (9) can be
rewritten as the following
y(t+1) = (1 − a) · y(t) + a · x(t)
x(t)
i ∼ Bernoulli
(
λ(t)
i
) (13)
Given the current position y(t) ∈ [0, 1]n within a unit hypercube, λ(t) as deﬁned in Eq. (9) deﬁnes
the tentative target where y(t) will evolve towards, in case of a deterministic system. Whil e in Eq.
(13), the new position y(t+1) is expressed as a convex combination of the current position y(t) and
the target λ(t), with combination ratios1−a and a respectively. The stability of the stochastic system
depends on two factors: the value a and where the ﬁxed-point y∗ is located. When a is closer to
1, y(t+1) has most of its part determined by the random sample rather th an the last position y(t).
Thus y along with time t evolves with large variation no matter wherey∗ is located. In particular, if
a > 0.5, we easily have
y(t)
i ∈ [0, 1 − a] ∪ [a, 1] , ∀i, t (14)
Thus no matter where y∗ is located, there is a certain region (1 − a, a)n centered at the center of
the hypercube where y(t) cannot reach at any time t. In case y∗ is inside that region, any trajectory
will be apart from y∗ with at least a ﬁxed distance, which means that none of the sta bility results
in [1] holds true. In another extreme, when a is very small, since we have that
⏐
⏐
⏐y(t+1)
i − y(t)
i
⏐
⏐
⏐ ⩽ a,
the next position y(t+1) will be within the range of a small hypercube of length 2a centered at y(t),
and the Lyapunov function V in this small hypercube will be approximately linear. We den ote the
linearized function V passing through y(t) as ¯V . Since the expectation of the next step is better in
terms of the Lyapunov function value, as the following
E
[
y(t+1)|y(t)
]
= (1 − a) · y(t) + a · λ(t) = y(t) − a ·
(
y(t) − λ(t)
)
(15)
When a is small, this approximately matches the evolution rule in t he deterministic system, which
by assumption in Eq. (12), guarantees to decrease the Lyapunov function value. Thus we have
E
[
¯V
(
y(t+1)
)
|y(t)
]
< ¯V
(
y(t)
)
(16)
When the linear approximation is good (in casea → 0), the above relationship with respect to¯V also
holds for V . This means that the trajectory, although being stochastic, will almost deterministically
converge to the ﬁxed-point and stay around there jumping ran domly with small perturbation. This
phenomenon can be easily seen by numerical simulation.
In particular, we found an interesting example. Suppose the parameters of the network are as the
following.
W =
[
0 20
15 0
]
, e =
[
−15
−10
]
, a = 0 .5, ǫ = 1 (17)
3
Figure. 1 illustrates the deterministic dynamical system deﬁned by Eq. (17). Table 1 (a) is∥y − λ∥2
plotted as color for each y, with red having high value and blue having low value. Table 1(b) is the
vector ﬁeld − (y − λ) showing how the deterministic system will evolve. The dynamical system has
two ﬁxed points y∗
1 = (0 .9922, 0.9925), and y∗
2 = (0 .0031e − 4, 0.4540e − 4). Starting from most
of the positions in the unit hypercube, the deterministic sy stem will go to the stronger ﬁxed-point
y∗
2. However, in the stochastic system as in Eq. (13), there are t wo situations. Most of the time the
trajectory goes to y∗
2 and stays there great stability. For very rare cases the trajectory goes to y∗
1 and
stays there for a while, but eventually jumps to y∗
2. It could takes long before such trajectory jumps
back to y∗
2, yet in our more than 10000 test trials, no single trajectory stays in y∗
1 forever.
0 0.2 0.4 0.6 0.8 1
0
0.2 
0.4 
0.6 
0.8 
1
0 0.2 0.4 0.6 0.8 1
0
0.2 
0.4 
0.6 
0.8 
1
(a) (b)
Figure 1: (a) Square of norm of r.h.s of Eq. (9); (b) The vector ﬁeld of r.h.s of Eq. (9).
The above example has very stable converged point, which is typical when the deterministic system
has one or more stable ﬁxed-points that are very strong ( i.e., located near the vertices of the hyper-
cube). In case the ﬁxed-point is not very strong, the situation is different. For example, if when both
∥W∥ and ∥e∥ are small, the following is also very small.


Wy(t) + e


 ⩽ ∥W∥ ·


y(t)


 + ∥e∥ (18)
In this case, λ(t) will have its entries close to 0.5 and thus each entry of x(t) will be very random.
Making y(t) evolves in changing directions. The stability of y(t) around y∗ still depends on how
large a is. In case the ﬁxed-pointy∗ has different strength in each of its entries, we found experimen-
tally that those stronger entries have more stable converged value, which is simply a combination of
what we illustrated in this section.
Experiments in our main text, especially Figure 3 and Figure5 (a) (b), demonstrated similar stability
in this section.
4 Open Questions
Presented in the following are several questions related to Bayesian inference by neural networks.
4.1 Conditional Speciﬁcation of Probability Models
The connection between neuron models and Bayesian inference algorithms as indicated in our main
text is based on the observation that the neural network structure, instead of representing the proba-
bility structure itself, can be used to represent the “infor mation ﬂow” of a certain message passing
algorithm on a probability model. In this case the probabili ty model refers to the Boltzmann ma-
chine and the message-passing algorithm refers to variatio nal inference and Gibbs sampling. As
mentioned in [4], many times variational inference has the mean-ﬁeld form of Gibbs sampling, thus
these two algorithms may share the same information ﬂow stru cture. However, not necessarily all
message-passing algorithms can be directly encoded in neural networks. An important requirement
for doing so is that the message each neuron sends to its downs tream neurons must be the same.
4
For example, belief-propagation (see [5]) and expectation propagation (see [6]) does not satisfy this
requirement, so it’s not clear at this point whether these algorithms can be directly encoded in neural
networks.
Assuming an inference algorithm can indeed be encoded in neu ral networks, the next question is
then, given an arbitrary neural network, when does it corres pond to a valid joint probability model
with the given inference algorithm? This is important, for example, if we are capable of carrying out
both inference and learning on the neural networks without a ssuming the joint probability model,
we can simply carry out that algorithm (which should be cheap er than actually learn probability
semantics), and very probably the learned model will exhibi t correct probability semantics in the
limit, in the sense that the network structure and parameter s under learning, converge to those that
do represent a joint probability model. For example, the pro of in [7] can be used to show that if
one carry out the generalized backpropagation algorithm [8] on Boltzmann machines, and preserve
the symmetric weights constraint, then error-driven learn ing can yield models with correct proba-
bility semantics, but what if the algorithm does not preserv e symmetricity, such as the Generalized
Recirculation Algorithm [9].
In the context of this paper, the question is that, suppose th e network structure is to represent local
conditional probability of a random variable, when does the speciﬁcation of a collection of local
conditional probability non-conﬂicting, and imply a uniqu ely determined joint probability model.
This is termed as the “conditional speciﬁcation model”, see[10] for a reference. By the methodology
in that book, it can be easily seen that, assuming each neuron corresponds to an individual binary
random variable and the input set of neurons represent the Ma rkov blanket, the only way to let
the neural network imply a joint probability model which is c ompatible to all speciﬁed conditional
probabilities is that the weights must be symmetric. In anot her way, the joint probability model
that the conditionals of a neural network represent has to be a Boltzmann machine. However, if we
allow a many-to-one mapping, i.e., different neurons can map to the same binary random variabl e,
the situation is different. As our main text implies, a neura l network without symmetric weights
can still represent a Boltzmann machine if two neurons speci fy the same binary random variable.
Furthermore, if we allow neurons to represent functions of random variables, or allow some hidden
neurons not representing any random variable, such that the ir purpose is only for the supporting of
the information ﬂow, then much more possibilities can happen. How do we even deﬁne this question
rigorously?
4.2 Universal Approximation by Distributions over Binarie s
Even though there are different ways of utilizing continuous variables in neural networks, it may still
be preferable to have all variables to be encoded in binary. Since Boltzmann machines are universal
approximator of distributions over binary variables [11],the natural question that follows is whether
it is also a good approximation to distributions over non-binary variables. If so, it would be sufﬁcient
for representing knowledge when developing intelligent systems.
The question is twofold: (1) Whether a continuous variable can be efﬁciently encoded by a collection
of binary variables, with the approximation error being controllable. (2) Whether distributions over
the continuous variables can be efﬁciently approximated by distributions over the encoded binary
variables.
When the continuous numbers are located on the unit interval , its binary approximation is natural.
We can simply take the binary digit representation and truncate it to have the ﬁrst set of digits. This
approach may be extended to domains where Fourier series functions can be well deﬁned, in which
case the sign of the Fourier series function evaluations can be used as a good binary approximation.
In case the continuous domain does not have Fourier series deﬁned or such approximation is inefﬁ-
cient, we may be able to look at the eigenfunctions of the Lapl acian operator. As indicated in [12],
there are relationships between Fourier series and eigenfu nctions of Laplacian operator. This may
provide a uniﬁed theoretical framework for investigating on binary approximation.
References
[1] Rafail Khasminskii. Stochastic stability of differential equations , volume 66 of Stochastic
Modeling and Applied Probability. Springer, 2nd edition, 2011.
5
[2] Lucien le Cam. An approximation theorem for the poisson b inomial distribution. Paciﬁc
Journal of Mathematics, 10:1181–1197, 1960.
[3] Barak A. Pearlmutter. Dynamic recurrent neural network s. 1990.
[4] Yee Whye Teh, David Newman, and Max Welling. A collapsed variational bayesian inference
algorithm for latent dirichlet allocation. In B. Sch¨ olkop f, J. Platt, and T. Hoffman, editors,
Advances in Neural Information Processing Systems 19 , pages 1353–1360. MIT Press, Cam-
bridge, MA, 2007.
[5] Judea Pearl. Probabilistic reasoning in intelligent systems: networks of plausible inference .
Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1988.
[6] Thomas Minka. Divergence measures and message passing. Technical report, Microsoft Re-
search, 2005.
[7] Joshua V . Dillon and Guy Lebanon. Stochastic composite l ikelihood. J. Mach. Learn. Res. ,
11:2597–2633, December 2010.
[8] Fernando J. Pineda. Generalization of back-propagatio n to recurrent neural networks. Phys.
Rev. Lett., 59:2229–2232, Nov 1987.
[9] Randall C. O’Reilly. Biologically plausible error-dri ven learning using local activation differ-
ences: The generalized recirculation algorithm. Neural Comput., 8(5):895–938, July 1996.
[10] B.C. Arnold, E. Castillo, and J.M. Sarabia. Conditional speciﬁcation of statistical models .
Springer, 1999.
[11] Nicolas Le Roux and Yoshua Bengio. Representational po wer of restricted boltzmann ma-
chines and deep belief networks. Neural Comput., 20(6):1631–1649, June 2008.
[12] Nathan Linial, Yishay Mansour, and Noam Nisan. Constan t depth circuits, fourier transform,
and learnability. J. ACM, 40:607–620, July 1993.
6


Peer review:
[Review]:
Thank you very much for the valuable reviews and references! I learned quite a lot from reading the suggested papers.

--> For Reviewer caa8: 

- Regarding the question raised in the end of your review, I think a somewhat related question is why neurons use spikes and whether we shall follow that in our computational model. I was previously approaching this question by exploring whether the combined semi-stochastic algorithm does something similar to simulated annealing, resulting in better local optima estimated in variational inference. But I failed to show this kind of superiority. In all the experiments I did, pure variational inference performs better than semi-stochastic inference in the context of classification with DBM on MNIST (probably because the DBM is well learned and the posteriors are well shaped with only one significant mode, also see the next paragraph). 

Another possibility of answering this is that, according to neural coding literature, such as works from Aurel Lazar, spikes are efficient ways of encoding time-varying real-valued functions if they are band limited (smooth in some sense). So neurons may be constrained by the energy they can use and have to choose the spiking approach. If that is the case, the 'randomness' of LNP is not that important, what's important of spikes are their property of reconstructing the function. And indeed more realistic neuron models such as Leaky Integrate-and-Fire (LIF) and Hodgkin-Huxley (HH) are not random (Chapter 5 of <6> provide a review about randomness vs. chaotic property of neurons). Are the spikes generated only to meet the reconstruction requirement? I think this worth further investigation. What I do found after submitting this paper is that, in the classification experiments I just mentioned, LIF networks work much better than LNP networks because the converged spiking pattern in LIF is periodic and the activation we get from convolving the recent spiking history is more stable, but if I replace the pseudo-random number generator in LNP by a quasi-random number generator (an effort to get more evenly distributed random samples), LNP and LIF behaves similarly. Another finding is that, in my recent GPU-cluster implementation of stochastic networks, the only way I can balance the data transmission with the computation is by transmitting packed binary spikes among computation nodes. So maybe LIF and HH are merely certain kinds of quasi-random generators used to do variational inference in an economic way. 

--> For Reviewer 4490

- I think the Poisson-unit BM direction is interesting, and I will definitely explore this in the future. Thanks for the review. In the supplemental material, I have a detailed justification about in what sense the discrete time Bernoulli model can be considered as approximations to the continuous time Poisson model. So for this paper, I think the foundation is OK.

- I agree to your comments that this paper is specific to the deep Boltzmann machine architecture. I didn't rule out the possibility that neurons can represent other models. But as long as they can represent Boltzmann machines with hidden variables, then the representation power is promised. As a compact universal approximator, even if behavior experiments reveals different types of probability models, they are not necessarily conflicting and may still be implemented by DBM. In addition, what's in my paper is extendable to high-order Boltzmann machines, as long as convolution with D function in section 2 of my paper can be eliminated (which I already validated on LIF and will put in my future publications). 

- As to whether it is enough to have positive weights only, I'm not sure. If the weights are symmetric, then the question is whether a Boltzmann machine with positive weights only is still universal for knowledge representation. I believe the answer is yes, because if yes, then the brain could be doing learning by using positive weights only in most of the normal time but add additional negative weights on demand when serious mistakes are made which need to be revised effectively. This way positive weights will be dominating but not exclusive. So yes I also think this issue needs to be delayed when we deal with neurally plausible learning. If Boltzmann machine with positive only weights and a constant bias is still universal, the original Boltzmann machine without the softmax manipulation in my paper could suffice for LNP modeling. 

- For whether the network will diverge to other parts of the posterior, I think this is very probable, since variational inference is a local optimum algorithm, and by my justification in the supplemental material, the semi-stochastic inference algorithm also approaches the local optima of variational inference. The good news is that by <3>, if a model is well learned, variational inference is good in the sense that the variational lower bound is close to the true likelihood, we can interpreting this result as the true posterior is highly single-moded. This alleviates the problem of local optima. 

- For low mean firing rate, I already considered this. A broader question is that if the nonlinear activation is not sigmoid (such as the Figure 1 in my paper), or if the maximum output of the activation function is low, can we still consider the LNP as a inference procedure. My answer is yes. In deriving the variational inference, we obtain the sigmoid activation from differentiating the KL-divergence loss function. If the activation is not sigmoid, we can easily reverse this derivation to obtain the new loss function, and by taking a difference to the original loss function, we get a regularization term. With different ways of fitting LNP to a Boltzmann machine, the regularization term can be different. Most of the time such regularization term is such that they favors low activity, which is reasonable. In this way, a different activation function can be regarded as a regularized variational inference. I can put more about this issue in the next revision if the paper is accepted. 

--> In the following I put a brief discussion of the reference papers provided by the reviewers. Please let me know if I made any mistake on my interpretation of these papers, as I read them in a hurry. 

- József Fiser, Pietro Berkes, Gergő Orbán, Máté Lengyel. Statistically optimal perception and learning: from behavior to neural representations', Trends Cogn Sci. 2010. 14(3):119-130

This paper previously inspired me a lot. Although there is no specific computational models proposed in this paper, most of the schemes which they discussed are what I'm following in my paper. For example, (1) they mentioned that sampling-based approaches, compared to parameter-based approach, has direct option for learning, but they lack the basis for experimental test. My work connects abstract computational models to LNP spiking neurons, allowing direct test (my on-going work right now extends this to Leaky Integrate-and-Fire and Hodgkin-Huxley as well). (2) They also mentioned that parameter-based approaches may suffer from the exponential explosion in the required size of neurons, while connectionist models such as Boltzmann machine rely on distributed coding which don't suffer from this problem. Furthermore, Boltzmann machines may be compact universal approximators. In this sense, modeling knowledge representation in terms of Boltzmann machines would be safe as long as the learnability issue can be addressed. (3) They also mentioned the 'spontaneous activity'. One of the issues when interpreting spikes as samples is that Monte-Carlo methods such as Gibbs sampling does not converge to the right distribution when all neurons sample together in parallel, while my approach provides an alternative viewpoint as a stochastic approximation of variational inference. The variational distribution one can get can be considered as an approximation of a mode of the joint posterior <4>, thus the activities in different neurons will be correlated according where the mode is located, this leads to an explanation of 'spontaneous activity'. (4) They also mentioned that inference and learning should be considered together when talking about representation. I didn't have it in my current paper, but my on-going work, built on this model, relates STDP to backpropagation, and by <1>, error-driven learning is hopeful to yield consistent probabilistic models with a properly chosen learning scheme. I will make these works available once the learning rule is tested in actual learning tasks. Also, if one favors likelihood-based learning such as contrastive divergence, the description on how hippocampus works in <2> implies that the positive phase of contrastive divergence, if implemented by variational inference as in <3>, can be preserved as short-term memory in the brain. The negative phase of contrastive divergence may be implemented by dreaming <5> (these early works about unlearning/reverse learning are about Hopfield Networks. But HN is a thresholded version of variational inference in Boltzmann machines with hidden variables, thus in terms of representation and learning, they are highly related to each other). 

- Reichert, D. Deep Boltzmann Machines as Hierarchical Generative Models of Perceptual Inference in the Cortex. Ph.D. Thesis. 2012.

For now as what I understand, I think this paper is more of an 'analogical model', the clues they use to relate the stochastic property of spiking to machine learning is through Gibbs sampling, which has the issue I mentioned above, while my focus is more on interpreting neural spiking as variational inference. 

- Shi, L., & Griffiths, T. L. (2009). Neural implementation of hierarchical Bayesian inference by importance sampling. Advances in Neural Information Processing Systems 22

To my understanding, this paper is about the implementation of importance sampler via one-hidden layer feedforward neural networks and then use it as building block to construct a hierarchical model for both top-down generative and bottom-up inference procedure. Thus if this is neurally plausible, it stands for other things that biological neural networks can do, which does not conflict to my work that neurons can perform approximate inference on DBM. The two lines of research can be separately proceeded.

- Ma WJ, Beck JM, Pouget A (2008) Spiking networks for Bayesian inference and choice. Current Opinion in Neurobiology 18, 217-22.

This paper is about the Probabilistic Population Code, which is another alternative for how neurons represent probability, belonging to the parameter-based approach discussed in the Fiser etc. 2010 paper above. 

- Pecevski D, Buesing L, Maass W (2011) Probabilistic Inference in General Graphical Models through Sampling in Stochastic Networks of Spiking Neurons. PLoS Comput Biol 7(12): e1002294

This paper is most interesting to me for now. I cannot comment on it before I read it carefully. I will come back with another post a bit later. 

Reference List: 

<1> Joshua V. Dillon, Guy Lebanon. Stochastic Composite Likelihood. Journal of Machine Learning Research 11 (2010) 2597-2633.

<2> O'Reilly, R. C., Bhattacharyya, R., Howard, M. D., & Ketz, N. (2011). Complementary learning systems. Cognitive Science

<3> Ruslan Salakhutdinov and Geoffrey Hinton. Deep boltzmann machines. Artiﬁcial Intelligence, 5(2):448C455, 2009.

<4> Thomas Minka. Divergence measures and message passing. Technical report, Microsoft Research, 2005.

<5> Francis Crick and Graeme Mitchison, The function of dream sleep, Nature 304, 111 - 114 (14 July 1983); doi:10.1038/304111a0.

<6> Wulfram Gerstner and Werner M. Kistler. Spiking Neuron Models: Single Neurons, Populations, Plasticity. Cambridge University Press, 1 edition, 2002.