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]:
This paper proposes a scheme for utilizing LNP model neurons to perform inference in Boltzmann Machines.  The contribution of the work is to map a Boltzmann Machine network onto a set of LNP model units and to demonstrate inference in this model.

The idea of using neural spiking models to represent probabilistic inference is not new (see refs. at end).  The primary contribution of this work is to take a learned deep Boltzmann machine from the literature, and to implement this network using LNP neurons, with the necessary modifications.  Therefore, the contribution is specific to the deep Boltzmann machine architecture.  The existing work in the literature often takes a different approach: taking realistic neural models and asking how these models can represent variables probabilistically.

Pros:
Developing mappings between machine learning algorithms and neural responses is an important direction.
To my knowledge, the implementation of a deep-BM with spiking neurons is novel.

Cons:
The clarity of the text and presentation of the mathematics needs improvement.
The resulting model suffers from some non-biological phenomenology.
The empirical results are not very compelling.
I would liked to have seen a comparison to the existing approaches for using spiking neurons to implement inference.  Particularly: [2-5]. Is there not a mapping from those models to the deep BM?  Why is the proposed mapping necessary, or what are the limitations of those existing proposals for a deep BM?


Other comments:
The paper provides a lengthy introduction to LNP and inference.  I would encourage the author to justify the various details that are introduced, and those that are not directly relevant for the proposed network, should be left out.  In general, the exposition needs clarification.

The proposed network seems like a logical series of steps, but the end result leads to a biologically implausible network, (at least when considering known properties of cortex).  I think a broad approach might be warranted for this problem. For example, starting from the LNP model and using this model as an element in a Boltzmann machine.

A related note: Isn't it just more plausible to estimate a deep-network with positive only weights? (to deal with Dale's law) There is likely some work to be done there, but it seems this direction wouldn't require the paired neurons you have here.  Or a network with realistic excitatory-inhibitory ratio?

Why not start with a Poisson-unit Boltzmann machine, and examine its properties? see (Welling et al. 2005)

I found the empirical evaluation to be weak.  I don't understand how running the network is a demonstration of correct inference.  Wouldn't we expect each of these networks to diverge and sample different parts of the posterior?

The statistics in Figure 5 need more justification.  I did not understand why these are relevant, or what degree of variability should be acceptable.

There are, of course, a variety of biological issues that seem to be incongruent with the proposal.
In cortex the distribution of excitatory to inhibitory neurons is 4:1.  The current proposal seems to require 1:1.
The pairing of neuron weights seems unlikely, but maybe this could be solved through learning?
What about mean firing rates?, are these consistent between model and cortical responses?

The title is a little misleading. I might suggest something more like:
Networks of LNP neurons are capable of performing Bayesian inference on Boltzmann machines.

The work of Shi and Griffiths 2009 seems highly relevant, and addresses some of the questions posed by the author.

Note that Dale's law is not generally applicable, but I am not sure about any refuttaion in cortex, which I assume is where you would imagine the deep network. (see co-transmission)

[1] Welling, M., Rosen-Zvi, M., and Hinton, G. E. (2005). Exponential family harmoniums with an application to information retrieval. Advances in Neural Information Processing Systems 17, pages 1481-1488. MIT Press, Cambridge, MA.

[2] Shi, L., & Griffiths, T. L. (2009). Neural implementation of hierarchical Bayesian inference by importance sampling. Advances in Neural Information Processing Systems 22

[3] Ma WJ, Beck JM, Pouget A (2008) Spiking networks for Bayesian inference and choice. Current Opinion in Neurobiology 18, 217-22.

[4] 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

[5] 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