Expressive probabilistic sampling in recurrent neural networks
Authors: Shirui Chen
Linxing Preston Jiang
Rajesh P. N. Rao
Eric Todd SheaBrown
Venue: NeurIPS 2023 poster (2023)
Keywords: neural coding
probabilistic sampling
neural dynamics
recurrent neural network
TL;DR: Expressivity analysis of recurrent neural circuits that can sample from any probability distribution

Abstract: In sampling-based Bayesian models of brain function, neural activities are assumed to be samples from probability distributions that the brain uses for probabilistic computation. However, a comprehensive understanding of how mechanistic models of neural dynamics can sample from arbitrary distributions is still lacking. We use tools from functional analysis and stochastic differential equations to explore the minimum architectural requirements for $\textit{recurrent}$ neural circuits to sample from complex distributions. We first consider the traditional sampling model consisting of a network of neurons whose outputs directly represent the samples ($\textit{sampler-only}$ network). We argue that synaptic current and firing-rate dynamics in the traditional model have limited capacity to sample from a complex probability distribution. We show that the firing rate dynamics of a recurrent neural circuit with a separate set of output units can sample from an arbitrary probability distribution. We call such circuits $\textit{reservoir-sampler networks}$ (RSNs). We propose an efficient training procedure based on denoising score matching that finds recurrent and output weights such that the RSN implements Langevin sampling. We empirically demonstrate our model's ability to sample from several complex data distributions using the proposed neural dynamics and discuss its applicability to developing the next generation of sampling-based Bayesian brain models.

Expressive probabilistic sampling in recurrent neural networks
Shirui Chen"
Linxing Preston Jiang
Department of Applied Mathematics University of Washington, Seattle sc256@uw. edu
Paul G. Allen School of Computer Science & Engineering
University of Washington, Seattle prestonjecs. washington, edu
Rajesh N. Rao
Eric Shea-Brown Department of Applied Mathematics Computational Neuroscience Center University of Washington, Seattle etsbθuw. edu
Paul G. Allen School of Computer Science & Engineering and Center for Neurotechnology University of Washington, Scattle raoêcs,, washington edu
Abstract
In sampling-basee based Bayesian models of brain function, neural activities are assumed be samples from probability distributions that the brain uses for probabilistic computation. However, a comprehensive understanding of how mechanistic models of neural dynamics can sample from arbitrary distributions still lacking. We use tools from functional analysis and stochastic differential equations to explore the minimum architectural requirements for recurrent neural circuits to sample from complex distributions. We first consider the traditional sampling model consisting network of neurons whose outputs directly represent the samples (sampler-only network). We argue that synaptic current and ftring-rate dynamics in the traditional model have limited capacity to sample from a complex probability distribution. We show that the firing rate dynamics of recurrent neural circuit with a separate set of output units can sample from an arbitrary probability distribution. We call such circuits reservoir-sampler networks (RSNs). We propose efficient training procedure based on denoising score matching that finds recurrent and output weights such that the RSN implements Langevin sampling. We empirically demonstrate our model's ability to sample from several complex data distributions using the proposed neural dynamics and discuss its applicability to developing the next generation of sampling-bast based Bayesian brain models.
Introduction
There growing evidence that humans and other animals make decisions by representing uncertainty internally and carrying out probabilistic computations that approximate Bayesian inference [42, 30, 18, 23, 34]. How networks of neurons in the brain represent probability distributions for Bayesian inference has remained major open question. There exist two major theories: one assumes that the neural activities encode the parameters of the underlying posterior distributions over sensory stimuli [5, 35, 52]. The other is the sampling-bass based hypothesis, which assumes that the neural responses can be interpreted as samples from posterior distribution [28]. Under this hypothesis, recurrent neural circuits make use of their inherent stochasticity to produce samples from posterior distributions. The
"Corresponding author
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

sampling-basee theory has explained various experimental observations regarding neural variability [19, 20, 40]. perceptual decision maaking [24] and spontaneous cortical activity [7].
Many studies have proposed biologically plausible spiking rules and membrane dynamics models to implement campling-based probabilistic inference. However, most of these studies mainly consider sampling from discrete Boltzmann distributions [12] and multivariate Gaussian distributions [17, 38. 25], only match the first two moments of the distribution [19], or employ Monte-Carlo approximation [29]. Although these studies use algorithmic substrates that can sample from any distribution with a density function theory, it is not clear whether the underlying neural dynamics are capable of implementing sufficiently expressive version of these sampling methods. Furthermore, natural image statistics are strongly non Gaussian [39], and experimental evidence shows that humans use non- non- Gan-Gassiinn prior representations to support cognitive judgments [27, 22]. It is known that deep artificial neural networks can used to generate samples from complex data distributions [48, 49] using a "U-net" [43] backbone. However, the neural circuits in the cortex are highly coupled with an abundance of recurrent synapses. Therefore, an outstanding question for probabílistic computation in the brain is: what kind of recurrent neural network is capable of efficiently learning to produce samples from an arbitrarily complex posterior distribution?
In this paper, we study this question under the basic assumption that the dynamics of recurrent neural eircuits can be described by stochastic differential equations (SDEs). Note that this assumption is common to broad range of past research that uses rate-based neural dynamics to implement sampling- based coding[1. 19, 17]. Moreover, spike-based models [46, 45, 38] that implement balanced spiking networks (BSNs) [11, 13] essentially train spiking networks to follow underlying continuous-time SDEs, so our work applies this line of research as well (details in Appendix G).
The contributions of this paper are as follows (Figure l):
1. We establish the relationship between the sampling power of the neural dynamics and the ability of the dynamics to approximate the score function, which is the gradient of the log probability density function. (Section 3.1)
2. We show that the synaptic current dynamics of a network of neurons whose outputs directly represent the samples (traditional sampler-only network) is only able approximate score functions that are in finte-dimennional function space. (Section 3.2, Proposition 2)
3. We prove that the firing rate dynamics of our proposed reservoir-sample network can sample from a distribution whose score function can approximate that of arbitrary target distributions (with mild restrictions) arbitrary precision (Section 3.3, Theorem 3).
4. We derive a computationally efficient and biologically-plausibbe learning rule for our pro- posed model (Section 3.4) that sidesteps the demands of backpropagation through time, and we empirically demonstrate how our model can sample from several complex data distributions (Section 4).
And interpretation of our contributions in biological terms is as follows: 1. Flexible synaptic connectivity within circuit itself is not enough to allow that circuit to flexibly produce arbitrary patterns of variability involving all of its neurons. 2. In order for circuit to achieve full flexibility in its output patterns, there need to be hidden variables involved, e.g. states of neurons in an upstream brain area possibly non-synaptic signalling. If this condition is met, concrete and fairly efficient plasticity rules may be capable of shaping the output patterns as desired.
2 Background
2.1 Fokker-Planck equation and stationary distribution
We consider a general time-homogeneous SDE with drift vector b (X)) and diffusion matrix Σ =
goo":
dxX, = b(X))dt odB,
where E Raxm is the diffusion coefficient, and B, is an m-dimensional standard Wiener process. The Fokker-Planck equation of this SDE describes the time evolution of the probability density p(æ; I)
for given SDE, assuming the initial density p(x, 0) is known:

---

Reservoir- Sampler Network
Neural Dynamics in SDE
Neural Samples
Reservoir
Samplar
7
9
7
8
W
Firing dynamics:
dr
- Wou!"
di (W D) di odB,
1
7
7
5
4
C
9
9
1
2
1
Wee
Approximation power
5
4
7
4
9
# Funo. basis
Fbxo?
Closed?
1
7
9
Noise
neurons
711 neurons
90-FR
(Lm)
No
No
50-50
O(m)
Yea
Sampler-Ont Network
RS-FR
O(s)
R6-8C
Inractazle
Intractable
Intrentable
Sampler
SC: Semp-e--rr Only Archiactuum
SC: Syrapt
Noise
W
FSI Husaazee Sumpperr Arct Raciure
PV ringg Hate
Synaptic currant dynomics:
dv D [- d.v)
odB,
Firing rate dynamics: dr = D [-r 4(W n] di edB.
r: firing rate
neurons
Figure 1: Reservoir-Sampler Networks versus Traditional Sampler-Only Networks. The sampler neurons which produce samples from the target distribution are present for both the sampler-only (SO) network and the esseroor-sample (RS) network. The reservoir neurons (shown within the large black circle) are only present in the RS network. We explore the firing rate (FR) dynamics (Section 3.3) for both networks and synaptic current (SC) dynamics (Section 3.2) only for the sampler-only network because the stationary distribution of the output neurons is intractable in the RS-SC case. We evaluate the approximation power the function set represented by the drift term the neural dynamics, and list the number of basis functions that span the set, whether these basis functions are fixed, and whether the function set is closed under addition (shown the table).
0.p = V (EVp bp)
where 0, a/0t, V the divergence operator, and is the gradient operator. A stationary probability density function is one for which the right-hand-side of equation (2) is 0. Mild regularity conditions guaranteeing the existence and uniqueness of stationary density function are discussed in, for example, Cerrai [14]. We assume these are satisfied by the SDEs considered in this paper. Moreover, to ensure ergodicity and well-defined score function, we assume that p supported on RR i.e. p(x) ) > 0 for all x E R..
As special case, consider the following Langevin dynamics, in which the drift term is given by the gradient of the log stationary probability density Vy(r),
dx, V Vlogp(X,)ee + V2dB
It can be verified through the Fokker-Plance equation that p(x: in the dynamics above is indeed a stationary probability density function of the dynamics. Therefore the Langevin dynamics can sample from the distribution p(x) as t co. While the drift term in the Langevin dynamics is gradient vector field, this is not true general and is often not the case for recurrent neural dynamics (Proposition 5 in Appendix B).
2.2 Score-based generative modeling
would like a particular dynamics to implement Langevin dynamics, we need to the drift term of an SDE to the score sx()) = V log p(x) of the probability distribution that we are trying to sample from. This procedure of fitting
score function ag xx(x) is called score matching. section, we give brief summary of one of the major methods of score matching that we will use in this paper, denoising score matching (DSM) [51]. The general idea of DSM is to match the score of noise-perturbed version of the target distribution. More specifically, the explicit score matching loss
this

---

(left-hand side) and the denoising score matching loss (right-hand hand side) are related to each other as follows [51]:
1 [|0o(X) - Vx go 9x)(1 = Eg (x,x)
1
E (8)
2 [100()) - Vx 9o (x|x))) + C.
(4)
where X is noise-perturbed version of on a but not on 0. When the noise is 0, i.c x = x, the left-hand side of the equation above is the explicit score matching loss. Although in theory noise, and directly optimize the right hand of equation (4), empirically is beneficial to start with large Gaussian noise and gradually decrease the noise magnitude until do (x(x) p(x) and se (x) R Vx log Go (x) 2 Vx log p(x). Specifically, this has been shown to stabilize the training and improve score estimation [48].
(x(xx x)) (xx o), and C' is constant depending
start from extremely small Gaussian
3 Methods
3.1 Do we really need to match the score?
The Langevin dynamics provides an elegant way construct an SDE given a specific stationary distribution. However, as noted above, the neural network dynamics seldom have drift term that is gradient field (Appendix B). natural question is therefore whether an SDE with a drift term that not a gradient field (also known as irvversibilitty) also gives rise to a specific stationary distribution. The answer requires us to look at the Fokker-Planck equation:
0.p(v. = V PEVV - pFo(v))
(5)
where p(v,() () is the probability density function of the variable of interest v at time I, Fo(v) is the drift term of the neural dynamies parametrized by and Σ _o00=" the diffusion matrix. Since the right-hand side of (5) needs to be O) for given stationary distribution p(v). we have
(EVp pF)(v)) O.
(6)
Therefore == -EVp + pVv(v) needs to be a divergence-free (DF) vector field. In other words, pVa(v) is unique up to DF field given fixed stationary distribution p. Ma et al. [36] shows that there exists a akew-symmetric matrix Q such that the DF field can be written as G QVp+p> EE; 0v, Qij- however, this does not shed more light on how expressive {Fota needs be without more knowledge of Q and its derivative. We show below that under certain conditions, the DF field G can be regarded as component that orthogonal to the score function. Therefore the function space { (Fo }o needs to have enough basis functions so that (when projected) it is able to approximate the score function of the target distribution.
We first note that it would be convenient if (v could approximate the gradient fields シー log for any p. in which case G O. To find out if this is necessary condition for the dynamics to sample from an arbitrary target distribution we let Σ I and invoke the Helmholtz-Hoe Hodge decomposition (HHD) [8. 37J. The decomposition theorem states that any sufficiently smooth vector field in (R"; R") can be uniquely decomposed into a DF vector field and a pure gradient field. In other words, the function space of all DF fields and the function space of all gradient fields Vp are the orthogonal complement of each other. Therefore the projection of pFo(v) onto the subspace of smooth gradient fields still needs to be able to approximate Vp despite the freedom to choose arbitrary DF field G.
For example, in the 1-D case, if we assume that both pV'(v) and Vp are square-integrable (so they vanish at infinity), then the divergence-free vector field G (which must be constant in 1-D) is o, therefore Fv(v)= = logp. As result, {Fe), indeed needs be able to approximate Vlog p for every p. In higher dimensions, the same conclusion holds under the assumption of a strict orthogonality constraint:
Proposition Let p be the stationary distribution of the neural dynamics, and the diffusion matrix be the identity matrix. If the DF field G is strictly orthogonal to the gradient field Vp. meaning that G(v). Vp(v)= = 0 for all v, then the drift term Fo Fo(v) can be written as the sum of a divergence-free field G and gradient field V log p.

---

The proof and further detail is in Appendix A. The form of the above decomposition coincides with that of the HHD [8, 9]. Therefore if we enforce the normal-paralia boundary condition [15] for the gradient component, the HHD theorem [8] says that the orthogonal projection of the function space ( onto the space of gradient fields is the function space of gradient fields (Vlog logg), satisfying the boundary condition given the strict orthogonality constraint on G. The upshot that ( needs to admit enough basis functions (Appendix A). Note that the boundary condition will not be restrictive if we take sufficiently large bounded region. Therefore it is essential for the neural dynamics to have an expressive functional form that is able to approximate complex score functions, even if their dynamies are not gradient fields.
Previous work has rigorously established that the strict orthogonality constraint holds, in particular, when the nonlinear dynamics linearized around fixed points of v (where the drift term zero [2, 31, 41J). As a consequence, the conditions of Proposition I are true locally around fixed points.
In what follows, assume that the conditions of the Proposition 1 hold. Under this assumption, without loss of generality, we set to be 0 in the following text and explore whether E- Fo. which is determined by the specific neural dynamics (Figure 1), is able approximate complex score function V lo
3.2 Synaptic current dynamics: sampler-only networks with limited capacity
We consider the following stochastic synaptic current dynamics [16] (cf. eq. 7.39) that describe a recurrent neural network in terms of the synaptic current that each neuron receives:
(7) -mm E R is the synaptic current of the IR mmmm is a diagonal matrix where diagonal elements E Rmxm is the connection matrix, d(.) is a nonlinear transfer function?, I I is the external input, a E R.xx is the diffusion coefficient and B, an I-dimensional standard Wiener process. The diffusion term can be interpreted as input from other brain areas (due to the large number of incoming connections, the assumption of Gaussianity can be justified by the central limit theorem [21]). We assume that 0 DD, W, I are tunable parameters through biological learning. To show the limited expressivity of FSC we have the following corollary from the Hilbert projection theorem showing that F is only able to approximate functions in a dineesedmenii function space.
dv = DD-v + W Wø(v)= () + d)t + adB, == F5C (v)dt (v)(t adB; where (V) D W WV(v)= I). v [11,
m neurons in the recurrent network, are the decay constants, d, = T;
E
Proposition 2. Lei H RTm + Rm" a function in the Hilbert space (2(",, p)"::)) Let II be the orthogonal projection operator onto the vector subspace
E (Av + Bø(() v + I|A, B E R _nxm I [R]
If - III > 0, then info H(()) F (v) 2 (111 IH
Proof of the proposition is given in Appendix B.I. The proposition says that no matter how the parameter of F5 is tuned, the difference between FSC and the target function cannot approach 0, and the lower bound of the error is given by the norm of the component in the target function that is orthogonal to the finite-dimensionl function space F. Therefore, synaptic current dynamics have a limited ability to match the score function and hence limited ability to sample from complex probability distributions under the strict orthogonality constraint in Section 3.1. The conclusion holds even if we let the diffusion coefficient be tunable. Since Σ-1 is linear, [FSC] je share the same set of basis functions as { FSC1 o As we will see in the next section, the firing rate dynamics of recurrent neural circuit with a separate output layer (a reservoir-sample network) and a learnable diffusion coefficient a can sample from arbitrary stationary distributions.
3.3
Firing-rate dynamics could be a universal sampler
3.3.1 Sampler-only networks
In this section, we consider the firing rate dynamics [16] (cf. eq. 7.11) that describe a recurrent neural circuit in terms of firing rates of the neurons. We first consider the sampler-only network:
dr D + o(Wrect wer + I))dt + odB, == F (v)dt + adB,
(8)
We later need it be non non-nolyloiial in Theorem 3
5

---

Here, (v) = D (WWeer + I)) and D is a diagonal matrix with decay constants as its diagonal elements. The stationary solution of the corresponding Fokker- Planck equation satisfies
V VEVVp pF!
where Σ == 200" is symmetric positive definite (SPD) and FØ D o( + 1)). Here 0 = {D, Weee } are tunable parameters. Σ is invertible, equation (9) becomes V V. (E(Vp pΣ p Therefore if Σ Fg V log p is gradient field, then the stochastic dynamies of the recurrent neural network described by equation (8) have a stationary distribution p" such that the score of this distribution V Vlogp" pd E FF.R
Compared to the synaptic current dynamics where we could only have functional basis V, and d(v).. we can now freely choose the functional basis spanning { [FFR), Jo depending on W rec and I, but since there is no linear term before the nonlinear transformation Φ, the function set {F5 }o is not closed under addition. If we view FFR as neural network with one hidden layer, the number of hidden neurons must be the same as the input dimension, and the diffusion matrix > (hence -ー is restricted to be an SPD matrix. Therefore, we do not get universal approximation power from Σ-1 FFFR and combined with results in Section 3.2, we see that an RNN by itself does not intrinsically produce samples from arbitrary distributions. As we will see below, if we let population of output neurons receive inputs from large reservoir of recurrently connected neurons (a reservoir-sampler network), we are able to obtain samples from complex distributions from the output neurons.
3.3.2 Reservoir-sampler networks
Now we consider the reservoir-sampler network where there is a linear readout layer W E mx of the reservoir whose dynamics is given by equation (8) (see also the upper row of Figure 1). As special case of Ito's 's lemma, we have
W Wondr ddWout Wout F. di + oont
(10)
= (-t Dr + Wot Do(Weec + 1)) + Wo
Now we assume that Wree is the product of W rec and Wout, i.e. Wrec = WTec out and D al is scaled identity matrix. If we denote the output of the recurrent neural network as x == Wour E IR". we derive the following stochastic dynamics for output neurons:
dx = (-ax + aW Wou ol x + I))dt + Won == F (x)d( + ōdB,
(11)
where F (x) = (-ox- + &W WWrex + I)) and ā = Won Therefore in order for the output neurons to sample from a stationary distribution p, we need sx(x) (- FFr
(x)
(PR (x) to match the score V logg(x).. Here 3 = (Wr.:: Wout, 1,0] are tunable parameters.
Additionally, we assume that F3 is 0 outside reasonable range of x. This assumption is used to prevent (x) from behaving wildly outside the bounded region on which sx(x) has the expressivity to match the score. The following theorem proves that with a large enough number of reservoir neurons, the moore-natching loss can be arbitrarily small. The proof is given in Appendix C.
Theorem 3. Suppose that we are given a probability distribution with continuously differentiable density function p(x): : Rm + RR. and score function V llog p(x) for which there exist constants Mj, M2, a, k >> 0 such that
p(x) < Mie a||x||
(12)
| lVggg(()) < M2 ||x||
(13)
when ||x]| > for large enough L. Then for any E > O, there exists a recurrent neural network whose firing-rate dynamics are given by (11), whose recurrent weights, output weights and the diffusion coefficient are given by Wres Raxn of rank m, W R and paxm respectively, such that, for large enough n, the score of the stationary distribution of the output units se (x) satisfies Exxp(xx [I] log so(x) [|²] <.
This theorem says that for any realistic data distribution with a smooth positive density function, there always exists a reservoir of recurrenty-cooneece neurons whose output units give samples from

---

distribution whose score function approximates that of the data distribution to arbitrary precision. Given the bound on the score matching error, Block al. [10] (cf. Theorem 13) gives bounds in Wasserstein 2-distance between the stationary distribution of the trained recurrent dynamics and the true data distribution. It is also worth noting that the recurrent weight matrix of the neural circuit in the theorem is of low-rank, so regardless of how many neurons there are in the reservoir, we can always find a low-rank recurrent weight matrix such that the output neurons sample from a correspondingly low-dimensional distribution.
3.4
An efficient RNN weight learning algorithm
The training procedure derived from the proof of Theorem 3. The main idea is to first train an auxiliary neural network with one hidden layer, then transfer the weights of the auxiliary neural network to the weights of the recurrent network that we are considering. More specifically, we first optimize an auxiliary feedforward neural network with one hidden layer Wout (Wrecx X + +I) using backpropagation with the denoising score matching loss (4) such that
( oW & (Wres) rec I) x) V Vlgg((x).
(14)
Then we can calculate the diffusion coefficient = W out (Wou W and the real recurrent weights Wree Wre.Woo out accordingly. The noise magnitude added the data samples is decreased exponentially over the entire training period. Figure 2 illustrates how the network can gradually learn the score function during training. We refer the readers Appendix D for more details.
Note that this method of using an auxiliary neural network much more computationally efficient than directly matching the score of the stationary distribution of the dynamics (11), for which the score function, which involves the matrix inverse (ōo²)- needs to be recomputed at each optimization step. Furthermore, since the entire training procedure is equivalent to training feedforward network with one hidden layer, it sidesteps the often challenging temporal computations associated with the Backpropagation Through Time (BPTT) algorithm used to train deterministic RNNs.
Although we assumed that the divergence free field G for the purpose theoretical analysis, in practice, fast sampling is major concern for implementing sampling-based inference models of the brain [25, 1, 19, 38), and reversible stochastic dynamics, i.e. if O, are known for their slow sampling speed. Fortunately, there is straightforward way to extend our framework and train RSNs to implement irreversible dynamics with a non-zero divergence-fre field G. This results in improvements in sampling speed (see Appendix E for details).
4 Experimental results
In this section, we present the results from two tasks. First, we try to let the recurrent neural network learn and generate samples from 1-D double-mode Gaussian mixture distribution and 2-D mixture of heavy-tailed Laplace distributions. For the second task, we explore whether a reservoir-sample network with firing rate dynamics is able sample from the distribution of internal representations computed from PCA filtering of image inputs. All dynamics are simulated with the Euler-Maruyama method. See Appendix H for hyperparameters used and other training details.
4.1 Learning mixture distributions
We consider a I-D Gaussian mixture distribution whose density function is the average of two Gaussian distributions centered at +1, i.c. Pdata (x) (NV 0.255 + (1, 0.25)). We artificially generate 10000 data points from this distribution and minimize the denoising score-matching loss:
[7]
[20 (Wre + +) x)
ā x
-- σ²
C(W out, W Wree; I) Exaaaa (c) Ea-
(15)
We use Φ(.) Re)U(-) and set a = 1/2. See Figure 2 for the numerical results. It is worth noting that we use tanh as the transfer function, the sampler-only networks are able to learn the score function perfectly, as the score function of the Gaussian mixture distribution we considered is exactly spanned by (x) = x and /2(x)) tanh(2..). See Appendix F. for an example where sampler-only networks fail to learn the score function even if hyperbolic tangent nonlinearity is used.
7

---

SO-SC
SO-FR
RS-FR
a)
b)
c)
100
5.0
true score
5.0
true SCOTO
5.0
true
2.5
2.5
Noise
score
0.0
0.0
0.0
10-1
-2.
2.5
2.5
-5.0
-5.0
-5.0
2
d)
8)
f)
sampled points
200
200
Scaled density Sampled points
150
150
150
100
100
100
50
0
2 0
0
V
X
Figure 2: Bimodal distribution sampling results. The 3 tractable cases shown are Sampler-only (SO) networks with both synaptic current (SC) dynamics and firing rate (FR) dynamics and Reservoir- sampler (RS) networks with FR dynamics, which are named SO-SC, SO-FR and RS-FR respectively. a-c) The score function learned compared to the true score function (orange curve) decrease the noise level (the darker the line, the lower the noise level). We see that RS-FR is capable of perfectly fitting the score function, while SO-SC and SO-FR are only able to fit the score function with piccewise linear functions when using the ReLU transfer function. d-f) Histogram sampled points, and the (scaled) density function of the target distribution. Again the reservoir-sampler network is able to generate samples whose distribution matches the target distribution, while the sampler-only network is not able to do due to the incorrectly matched score function.
gradually
Next, we show that the model can learn mixtures of heavy-tailed distributions that are ev- ident in natural image statistics and the neural representations in the primary visual cortex [47, 39J. We trained the Reservoir- Reservoir-Sampler network with FR dynamics (RS-FR) model on 20000 sampled data points from 2-D Laplace mixture distribution, whose density is given by
0.9
1
-0.9]
Paata (x) Lap
0,
Lap 0.
0.9 -09
where Lap denotes the multi-
0.9
1
variate Laplace distribution. The model successfully learned the probability density of the mixture distribution (Figure 3 left vs. middle), and captured the heavy tails of the distribution measured by the kurtosis (Figure: 3 right).
4.2 MNIST generation task
We also tested the sampling ability of our model using the MNIST dataset [32] which contains 60,000 handwritten digits from 0 to 9. We projected MNIST images to a 300-D latent space spanned by the first 300 principal components, and trained the weights of the recurrent neural network as described in Section 3.4 so that the RNN can sample from the latent distribution. To test the model, we generated images by applying inverse PCA projection to samples generated by the model. The schematics and generated images are shown in Figure 4. Note that since we are essentially using a shallow network to match the score, we should not expect comparable performance to generative models that use deep ANNs. Our main goal to illustrate that the reservoir-sampler network using firing rate dynamics qualitatively more expressive than other traditional neural sampling models (Appendix G). Finally we also note that it highly nontrivial for recurrent neural dynamics to complete such generative task, and to the best of our knowledge, previous work has achieved such results.
8

---

True
Generated
sample
Generated sample
Kurtosis
X2
1.5
1,0
0.5
0.0
X1
X/2
Figure 3: RS-FR model learning a mixture of 2-D Laplace distribution. Left to right: sample density from the true distribution (brighter color denotes higher density); sample density from the learned distribution; marginalized kurtosis of each dimension from the true and learned distribution.
a)
Training
Reservoir
Sempler
Sampl ng
RS-FR
Training Samples
Wow
9
8
3
8
2
9
4
1
7
5
4
2
0
9
2
9 | / 8
9
2
3
8 6
8 M
S
4
7
9
3
Inx PCA
1
9
5
3
3
0
Ncise
b)
d)
f)
SO-SC
SO-FR
SO-FR
3
RIS-RR
72
M
Lo89
24
4
8
4
2
47
4
27
4
3
27
4
200
400
600
800
1000
Epoch
Figure 4: Learning to sample the MNIST image distribution. An MNIST image is projected to 300-D latent space (orange circle) spanned by the first 300 principal components using PCA. The sampler learns to sample from the distribution of this latent space and generates images by applying inverse PCA projection to these samples. The diagram illustrates the RS-FR model. b) The loss curves for different models during training. Every 100 epochs the noise level added the training samples is reduced (Appendix H), and the noise increases to higher value because the score matching loss magnitude depends on the noise level. As shown in the inset, the loss of the RS-FR model decreases throughout the training process when using the lowest fixed noise level. Meanwhile, the losses of the other two models remain unchanged. c-f) The images generated for the models compared to the digit images generated from latent training samples.
5 Discussion
From the perspective of functional analysis and SDE theory, we prove that under the strict orthog- onality constraint, it is essential for neural circuits to have drift term that has the expressivity to approximate complex score functions, despite the fact that dynamics do not have to exactly implement the Langevin dynamics. We investigated whether population of neurons can sample from an arbitrary distribution directly and proved that the drift term of the synaptic current dynamics only sample from dinee--oimssinnall function space. Although the drift term of the firing

---

rate dynamics can approximate functions spanned by different basis functions, the number of basis functions is limited. To address this problem, we proposed the reservoir-sample network for firing rate dynamics. We found that with learnable diffusion coefficients and a sufficiently large reservoir of hidden neurons, the output neurons described using the firing rate dynamics are able to sample from arbitrary data distributions. Our results partly answer the question of what architecture recurrent neural circuits need so that they are able to sample from complex data distributions.
Our analysis and empirical experiments affirm the universality of stochastic RNNs. However, this universality comes with limitations. First, we have only analytically shown the existence of weights that enable sampling from complex data distributions; there is no guarantee that one will find such weights through backpropagation. Additionally, in order to obtain the tunable diffusion coefficient during training. a matrix inverse is needed (likewise in the FORCE algorithm [50]). Further, the question of how biological circuits compute the specific gradient and implement the denoising score- matching algorithm remains an open question. Moreover, in our current formulation, we are only able to approximate the score function with a shallow network with one hidden layer. Our preliminary experiments show that one-hidden-layer RSNs cannot readily approximate high-dimensional heavy- tailed distributions (e.g., those of overcomplete sparse coding representations [39]). It is unclear this is because of insufficient number of reservoir neurons. Due to the limitation of the GPU memory, we did not try higher number of reservoir neurons.
Our model differs from recent diffusion models [26, 49], which can be seen as time- imm-innooooeeeeuus SDEs, and has the advantage of being able to run indefinitely in time, making a suitable candidate for modeling spontaneous activity in the brain. Moreover, while Song and Ermon [48] optimizes the denoising score matching loss at different noise levels jointly, we adopt a sequential learning procedure by gradually decreasing the noise level of the training samples. This procedure is more aligned with the developmental processes involved in forming visual representations in the infant brain, where the distribution of visual representations are thought to be noisier (less linearly separable) initially [4]. Our study therefore serves as a starting point for building a mechanistic model for probabilistic computation in the brain that has similar generative power to current Al generative models.
Biologically, there are multiple ways to interpret the reservoir neurons and sampler neurons in an RSN. First, reservoir and sampler neurons could be seen as different types of neurons a single brain area, where the dynamics of sampler neurons converge quickly to the equilibrium point. Second, even more straightforwardly, the sampler neurons could be seen as a more separate set of neurons located downstream of the reservoir. We also wish to suggest an alternative interpretation. Biological neural networks are known to have non-synaptic signaling networks (e.g. pervasive neuropeptidergic signaling [3], extensive aminergic signaling [6] or potential extrasynaptic signaling [53]) in addition to the synaptic connectivity that is typically modeled (i.e., via connection weights). We suggest that it is possible that the computations of the reservoir may be implemented by non-synaptic networks, and then "read out" by certain neurons' spikes. This possibility is supported by the recent finding of a low correlation between functional activity and the synaptic ("structural") connectome in C.elegans [53]. Morever, if we only take the structural connectome into consideration, then the resulting model of C. elegans would correspond to the sampler-only network, which, as theory predicts, will have limited sampling capability.
6 Conclusion
In this paper, we explore how recurrent neural circuit can sample from complex probability distributions, an important functional motif in probabilistic brain models. We start from a basic assumption that the recurrent neural circuit could be described as an SDE. We show that recurrently- connected neural population by itself has limited capability to implement stochastic dynamics that can sample from complex data distributions. In contrast, we prove that firing rate dynamies of the output units of a recurrent neural circuit (a reservoir-samme of probability distributions. These theoretical results, together with our preliminary experimental results, provide sufficient condition for neural sampling-base models to exhibit universal sampling capability. Our results therefore provide a foundation for the next generation of sampling-bases probabilistic brain models that can explain wider range of cognitive behaviors.
network) can sample from a richer range
10

---

7 Acknowledgements
We are thankful Profs. Hong Qian, Bamdad Hosseini and Edgar Walker for their guidance and insight with this project. We gratefully acknowledge the support of the grant NIH BRAIN R01 IRFIDA055669.
References
[1] Laurence Aitchison and Máté Lengyel. The Hamiltonian Brain: Efficient Probabilistic Inference with Eacitatory-Inhibitoo Neural Cireuit Dynamies. PLOS Computational Bi- ology, 12(12): 186, December 2016. ISSN 1553-7358. 1005186. URL https: j/jurna//// -aaenallslllss.....llm.ooo journal. Cournal.pcbi...
doi: 10.1371/journal.pebb
http:///sccppiooaatiitleeeeeeeee 1371/ Publisher: Public Library of Science.
[2] P. Ao. Potential in stochastic differential equations: novel construction. Journal of Physics A: Mathematical and General, 37(3):L25, January 2004. ISSN 0305-4470. doi: 10.1088/
0305-4470037//100L URL https: / ht////xx.doo// doi. org/ 10. 0tt5//003055444/7773////
[3] Cornelia Bargmann and Eve Marder. From connectome to brain function. Nature methods, 10(6):483-490. 2013.
[4] Laurie Bayet, Benjamin D. Zinszer, Emily Reilly, Julia K. Cataldo, Zoe Pruitt, Radoslaw M. Cichy, Charles A. Nelson, and Richard N. Aslin. Temporal dynamics of visual representations in the infant brain, August 2020. URL https //tp:://wwwww/ h/www.biorxxv biorxiv. wwaordvvoo//oooteeeeeee wetgg//oneet///tee 1101/2020. Pages: 2020.02.26.947911 Section: New Results.
02. 26. 1S9947911vv
[5] Jeffrey M. Beck, Wei Ma, Roozbeh Kiani, Tim Hanks, Anne K. Churchland, Jamie Roitman, Michael N. Shadlen, Peter E. Latham, and Alexandre Pouget. Probabilistic population codes for Bayesian decision making. Neuron. 60(6):1142-115: December 2008. ISSN 0896-6273. doi: 10.1016/j.neurrg wwtuuron.2000..oo 20.8.00.00000000 09.021. URL https:/ // //tp://wwww./.. wwbi.cbi.... nlm. nih wath.govpmm////m httpp//atales//////////m.essslet PMC2742921/.
[6] Barry Bentley, Robyn Branicky, Christopher I. Barnes, Yee Lian Chew, Eviatar Yemini, Ed- ward T Bullmore, Petra E Vértes, and William R Schafer. The multilayer connectome of caenorhabditis elegans. PLoS computational biology, 12(12):e100528, 2016.
[7] Pietro Berkes, Gergo Orban, Máté Lengyel, and József Fiser. Spontaneous Cortical Activity Reveals Hallmarks of an Optimal Internal Model of the Environment. Science, 36013))) 83-87, January 2011. ISSN 0036-8075, 1095-9203. doi: 10. 126/science. 1195870. URL
https: htopp::/www...e www.wscience..ommmmmmmmmmm /tt..////...11//////////// scwence..cc.ccccco..oo 1195870.
[8] Harsh Bhatia, Gregory Norgard, Valerio Pascucci, and Peer-Timo Bremer. The Helmholtz- Hodge Decomposition-A Survey. JEEE Transactions on Visualization and Computer Graphics, (988):386-1400 August 2013. ISSN 1941-0506. doi: 10.1 10.1109TTCC.20122. Conference Name: IEEE Transactions on Visualization and Computer Graphics.
[9] Harsh Bhatia, Valerio Pascucci, and Peer-Timo Bremer. The Natural Helmholtz-Hodge Decomposition for Open-Boundary Flow Analysis. IEEE Transactions on Visualization and Computer Graphics, 20(11):1566-1 1578, November 2014. ISSN 1941-0506. doi:
10. 111110CC...004420 20112 Conference Name: IEEE Transactions on Visualization and Computer Graphics.
[10] Adam Block, Youssef Mroueh, and Alexander Rakhlin. Generative Modeling with Denoising Auto-Encoders and Langevin Sampling, October 2022. URL http: http:::/axxxvv.oopppe wowpg////b//////// 00107. arXiv:2002.00100 [cs, math, stat].
[11] Martin Boerlin, Christian K. Machens, and Sophie Denève. Dynamical Variables in Balanced Spiking Networks. PLOS Computatioonal Biol-
Predictive Coding of
ogy, 9(1 10032258 November 2013. ISSN 1553-7358.
doi: 10. 10.1371
1003258. URL https:: https://journalss h/jjorrnalssppppp plarnaalsplos.......m plos. orgppleccoppoo//rtcicceetooo.
journal, pcbi 1003258. Publisher: Public Library of Science.

---

[12] Lars Buesing, Johannes Bill, Bernhard Nessler, and Wolfgang Maass. Neural Dynamics as Sampling: A Model for Stochastic Computation Recurrent Networks of Spiking Neu- rons. PLOS Computational Biology. 7(1 7(11)2e100221, November 2011. ISSN 1553-7358.
doi: po. 1002211 URL https:/ B jorn/all.lll..a// plos. org/ ppspppoccompbbll// Publisher: Public Library of Science.
article?id=10. 10loo.117ooooooooa pcbi 1002211,
[13] Camille E. Rullán Buxó and Jonathan W. Pillow. Poisson balanced spiking networks. PLOS Computationnal Biology. 16(1 t)ee 1008261, November 2020. ISSN 1553-7358. doi: 10.1371/ journal.pcbi. 1008261. URL https: https:///ooraiii journalss peremaassplo......... wtt/ploggeeeoooooooooooooooooooo httpa//ppp.o//pppoossoooooooooo.
10. 10.1371/ 1371/journal pcbi..pi...... 1008261. Publisher: Public Library of Science
[14] S. Cerrai. Second Order PDE's in Finite and Infinite Dimension: Probabilistic Approach. Number no. 1762 in Lecture Notes in Mathematics. Springer, 2001. ISBN 978-3-540-42136-88 URL https:
google e. cww/bookee waabokkksddrpDDeupkkk--.
[15] Alexandre J. Chorin and Jerrold E. Marsden. A Mathematical Introduction Fluid Mechanics, volume 4 of Texts in Applied Mathematics. Springer, New York, NY, 1993. ISBN 978-1- 4612-6934-2 978 1--612-0883-9 doi: 10.
URL http://link:
springer com 10. /978
[16] Peter Dayan and Laurence Abbott. Theoretical Neuroscience: Computational and Math- ematical Modeling of Neural Systems. MIT Press, August 2005. ISBN 978-0-262-54185555
Google-Books- ID; ILT4DWAAQBAJ,
[17] Xingsi Dong, Zilong Ji. Tianhao Chu, Tiejun Huang, Wenhao Zhang, and Si Wu. Adaptation Accelerating Sampling- based Bayesian Inference in Attractor Neural Networks. October 2022. URL
[18] Kenji Doya, Shin Ishii, Alexandre Pouget, and Rajesh Rao. Bayesian brain. Probabilistic approaches to neural coding. January 2007.
[19] Rodrigo Echeveste, Laurence Aitchison, Guillaume Hennequin, and Máté Lengyel. Cortical- like dynamics in recurrent circuits optimized for compling-baset probabilistic inference Nature Neuroscience. Nature Neuroscience, (3(9)1118-11444 September 2020. ISSN 1546-1726. doi: 10. 10000//4119300000007--11111 URL https: //www//wwww///// nature. com/ com/articles/
841593-020-0671--1
[20] Dylan Festa, Amir Aschner, Aida Davila, Adam Kohn, and Ruben Coen-Cagli. Neuronal variability reflects probabilistic inference tuned to natural image statistics. Nature Commu- nications, 12(1):3635, June 2021. ISSN 2041-1723. doi: 10. 10.108//444467-02 1467-021-22838-1 URL https: //tp:://www www.wr............m.a. wttp///rticceeee//11100000----20
x. Number: Publisher:
Nature Publishing Group.
[21] Wulfram Gerstner, Werner M. Kistler, Richard Naud, and Liam Paninski. Neuronal dynamics: From single neurons to networks and models of cognition. Cambridge University Press, 2014.
[22] Thomas L. Griffiths and Joshua B. Tenenbaum. Optimal predictions in everyday cognition. Psychological Science, 1709)767-773, September 2006. ISSN 0956-7976. doi: 10.1111/j.. wtt.99880.00001788000
[23] Thomas L. Griffiths, Nick Chater, Charles Kemp, Amy Perfors, and Joshua B. Tenenbaum. Probabilistic models of cognition: exploring representations and inductive biases. Trends in Cog- nitive Sciences, 14(8):357-364, August 2010. ISSN 1364-6613. doi: 10. 10.11116..icss.2100005.. URL
/www sci http:///wwaccccccooooooccccoaaaa
[24] Ralf M. Hacfner, Pietro Berkes, and József Fiser. Perceptual Decision as Probabilistic Inference by Neural Sampling. Neuron, 90(3):649-661) May 2016. ISSN 0896-6273, doi: 10.
1016//neeuron22016.0.00222 URL https: hwwp://wwwwww////oo///// www.sssciiiccceccccccmmmme p11/$08966223366000113
[25] Guillaume Hennequin, Laurence Aitchison, and Mate Lengyel. Fast Sampling-Based Inference in Balanced Neuronal Networks. In Advances in Newrul Information Processing Systems, vol- ume 27. Curran Associates, Inc., 2014. URI.
wtwpeeeeddinnnnunnnnnneeeeeeeeee neurips. meurips.cc/paper. cc/paerr/ html.
2014/ 11144//hsh/ah httaaa/ess..00eeeeeeeeee0eaaaaaa
12

---

[26] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising Diffusion Probabilistic Models. In Advances in Neural Information Processing Systems, volume 33, pages 6840-6851. Cur- ran Associates, Inc., 2020. URI htRt.htttss::/pooceeeiiii iegpeeedddnnnnnnnniiisseee
[27] Neil M.T. Houlsby, Ferenc Huszár, Mohammad M. Ghassemi, Gergo Orbán, Daniel M. Wolpert, and Máté Lengyel. Cognitive Tomography Reveals Complex, Task-Independent Mental Representations. Current Biology. 23(21)2216-22777 November 2013. ISSN 0960-9822. doi: 10. 1u..006/j.cbb2211..9900000 URL https: //www//wwwwwww/ wwbi ncbi nlm. nih, wtw.gg//////icccciie./ PMC3898796/.
[28] Patrik Hoyer and Aapo Hyvärinen. Interpreting Neural Response Variability Monte Carlo Sampling of the Posterior. In Advances in Neural Information Processing Systems, vol- ume 15. MIT Press, 2002. URL. hattpss:/proceeddnngg wttpssserer.....pppre////// http:sccccppppe/22////aaaea 00666oo766ac33222777762444433343 Abs.eecs---sttttt.ttt.....
[29] Yanping Huang and Rajesh PN Rao. Neurons as Monte Carlo Samplers: Bayesian Inference and Learning in Spiking Networks. In Advances in Neural Information Processing Systems, vol- ume 27. Curran Associates, Inc., 2014. URL https: /ttp//ppoceeeg...//ee/ neg www.nie...eeepeeeeeeeeee cc/paper/ 2014/ 2014/hash//oe httaaa//.caaaa..........7o07....
Abstract html.
[30] David C. Knill and Alexandre Pouget. The Bayesian brain: the role of uncertainty in neural coding and computation. Trends in Neurosciences, 27(12):712-719, December 2004. ISSN 0166-2236. doi: 10.1016/j..io wtw.:///..iis.2........
URL https://wwes wwwencedirect.......m..... com/
[31] Chulan Kwon, Ping Ao, and David J. Thouless. Structure of stochastic dynamics near fixed points. Proceedings of the National Academy of Sciences, (2237)):3029-1300 September 2005. doi: 10.107 1//////as..aaammmmm httpas..0000347/111.. 0506347102. URL https: /t////www pnas org/ /oo//0..... 10.10000000 107 3/ 1073/pnas. pnas. 0506347100 102.
Publisher: Proceedings of the National Academy of Sciences.
[32] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner. Gradient-based learning applied document recognition. Proceedings of the IEEE, (((1)))277-23224 November 1998. ISSN 1558-2256. doi: 10.1 (0.1109/5.72691, Conference Name: Proceedings of the IEEE.
[33] Moshe Leshno, Vladimir Ya. Lin, Allan Pinkus, and Shimon Schocken. Multilayer feedforward networks with a nonpolynomial activation function can approximate any function. Newral Ner- works, 6(6):861-866) January 1993. ISSN 0893-6080. doi: 10. 10.10016/00833--000005550031111
https: //wp://///www.//// www.scciieeecete.cccemeeeeeeeeee com coom//ccian htt.:////.aaaae/ee..///////////a
[34] Wei Ji Ma and Mehrdad Jazayeri. Neural Coding of Uncertainty and Probability. Annual Review of Neuroscience, 37(1):205-22 2014. doi: 10.1 10.1144//amure--eu00000111333330
URL attpp://don http::://ooooooooo org/10
.eprint:
https///doi.org/10. 46
[35] Wei Ji Ma, Jeffrey M. Beck, Peter E. Latham, and Alexandre Pouget. Bayesian inference with probabilistic population codes. Nature Neuroscience, 911): 1438, November 2006. ISSN 1546-1 1266. doi:
URL https://www nww.raateee.... Number: 11 Publisher: Nature Publishing Group.
[36] Yi-An Ma, Tianqi Chen, and Emily Fox. A complete recipe for stochastic gradient meme. Advances in neural information processing systems, 28, 2015.
[37] Andrew J. Majda and Andrea L. Bertozzi. Vorticity and incompressible flow. Camb. Texts Appl. Math. Cambridge University Press, Cambridge, 2002. ISBN 978-0-521-63057-3 978-0-521- 63948-4. doi: 11L10177CC009780511612222
[38] Paul Masset, Jacob A. Vavatone-Veth, J. Patrick Connor, Venkatesh N. Murthy, and Cengiz Pehlevan. Natural gradient enables fast sampling in spiking neural networks. October 2022.
13

---

[39] Bruno A. Olshausen and David J. Field. Emergence of simple-cell receptive field properties by learning a sparse code for natural images. Nature, 381(6583)6077-609) June 1996. ISSN 1476-4687. doi: 10. (0.1008//81600
http:: :/t::///wwwwww http://twtwt..eee//////////t nature. com com.arrtiiclees/31167aaaa Number: 6583 Publisher: Nature Publishing Group.
[40] Gergo Orban, Pietro Berkes, József Fiser, and Mäté Lengyel. Neural Variability and Sampling- Based Probabilistic Representations in the Visual Cortex. Neuron, 92(2):530-544 October 2016. ISSN 0896-6273. doi: 10. 10.10166 6.0.088. URL https: cell com/ neuron/ aaaroo/abbsract//0000 9063 O. Publisher: Elsevier.
6273 30639-0
[41] Hong Qian. The zeroth law of thermodynamies and volume- preserving conservative system in equilibrium with stochastic damping. Physics Letters A. 378(7):609 616, January 2014. ISSN 0375-9601. doi: 1...016f..hysseta.211..... science/artionn ccitemee//atticeee///////8//1//
URL https: htap::://www /ww. sciencedirect. com/
[42] Rajesh P. Rao, Bruno A. Olshausen, and Michael Lewicki, editors. Probabilistic models of the brain: Perception and neural function. Probabilistic models of the brain: Perception and neural function. The MIT Press, Cambridge, MA, US, 2002. ISBN 978-0-262-18224-9 doi: 10.755 / www.550000000000.11.1.
Pages: x, 324.
[43] Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-Net: Convolutional Networks for Biomedical Image Segmentation, May 2015. URL http:://rra httt:::/araii..og/aba//a/a/ orgg/bs//1055 04597. TrXiv:1505.0555 [cs].
[44] W. Rudin. Functional Analysis. International series in pure and applied mathematics. McGraw- Hill, 1991. ISBN 978-0-07-05423665. URL htttps:/www. hwt.ggle..ombboooo/eeeee Functional wtw.aaaa....aaaaaaaaaa..AAAA.AaA
[45] Camille Rullán Buxó and Cristina Savin. A sampling-based circuit for optimal decision making. In Advances in Neural Information Processing Systems, volume 34, pages 14163- 14175. Curran Associates, Inc., 2021. URL https::/ neurips httpsscccppppr//ppppppppp httm..
0011///aa/0444333312222227711111
[46] Cristina Savin and Sophie Denève. Spatio-temporal Representations of Uncertainty in Spiking Neural Networks. In Advances in Neural Information Processing Systems, volume 27. Cur- Associates, Inc., 2014. ht htt.:/5000300ee0000077----beeeee 42244ff8.ee777005550087777666666
wewpiinnnnnnnnnnneeeeeeeeeeeeeee nww.nnuursccc.m wwtp:ips.ccppppp/1e/4/11 html.
[47] Ecro P Simoncelli and Bruno A Olshausen. Natural Image Statistics and Neural Repre- sentation. Annual Review of Neuroscience, 24(1):1193-1255 2001. doi: 10.1 10.1146/annurev neuro.24.1.1 1193. URL https: //doi///ooi//// org/ //44/re neuro. 24. 1193. .eprint: http:////arrrrr////./////aaaeaae 140/////mmurrr.emuoo..11.1111
[48] Yang Song and Stefano Ermon. Generative Modeling by Estimating Gradients of the Data Distribution. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. URL https: http///peeeeenn.eeeeee////e
/ neurips cc/paper/2019/hat) Awwtract..tttt........aa. html.
[49] Yang Song. Jascha Sohl-Dickstein, Dicderik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-Based Generative Modeling through Stochastic Diffferential Equations, February 2021. URL http:://arrtiv.org//b////.1..
13456. 1:X1 [cs, stat].
[50] David Sussillo and L. F. Abbott. Generating Coherent Patterns of Activity from Chaotic Neural Networks. Neuron, 63(4))544-557, August 2009. ISSN 0896-6273, doi: 10160 J.neuron.209.07..00 07.018...111..1 URL ottpr//etttttttteoooooeaaaaeaaee 50896627309005479.
[51] Pascal Vincent. A Connection Between Score Matching and Denoising Autoencoders. Neural Computation, 23(7):1661-167 July 201 ISSN 0899-7667. doi: 10.1162/NFCO...000
Conference Name: Neural Computation.
14

---

[52] Eszter Vértes and Maneesh Sahani. Flexible and accurate inference and learning for deep generative models. Advances Neural Information Processing Systems, volume 31. Cur- ran Associates, Inc., 2018. URL https::/
wwweinng.....uiieeeee nwwrinuurncccom wwt..ee.....ppppeeee.eeeeeeeeeee wwt.cc/pppppe///////aaaaeee 45ccc557555844444222222667776-63 htmml.
[53] Eviatar Yemini, Albert Lin, Amin Nejatbakhsh, Erdem Varol, Ruoxi Sun, Gonzalo Mena, Aravinthan DT Samuel, Liam Paninski, Vivek Venkatachalam, and Oliver Hobert. Neuropal: a multicolor atlas for whole-brain neuronal identification in c. elegans. Cell, 184(1):272-288) 2021.
15

---

A Helmholtz-Hodge decomposition and proof of Proposition 1
First we give a proof for Proposition 1:
Proposition. Let be the stationary distribution of the neural dynamics, and the diffussin matrix be the identity matrix. If the DF field G is strictly orthogonal the gradient field Vp, meaning that G(v) V p(v) = 0 for all v. then the drift term Fo(v) can be written as the sum of divergence-free field p G and gradient field logp.
Proof. Since is the stationary distribution of the neural dynamics, from the Fokker-Planck equation we have
Fo(v) = p G V logp
(16)
where G is divergence free. Now
V G) = Vp-
G - G-
Vp-G
G is divergence free G and Vp is pointwise orthogonal
p2
0.
Therefore, P is also divergence free.
Now we give statement of the Helmholtz-Hove Hodge decomposition, which is an application of the Hodge decomposition theorem to vector calculus on R" [8, 37].
Theorem 4 (Helmoholtz-Hodgg Decomposition). A smooth vector field V on bounded domain of IR" can be uniquely decomposed into the sum of a gradient field Vu that is normal to the boundary, divergence-free field v that is parallel to the boundary and a harmonic field h, i.e. V == Vu+ h. Moreover, the three components are orthogonal to each other in the L2 sense.
We note that the form of the decomposition in Proposition coincides with the Helmoholtz-Hor Hodgge decomposition with zero harmonic field. Therefore, we can define the orthogonal projection operator that projects Fg onto its gradient component log Under the strict orthogonality constraint that requires G to be orthogonal to Vp pointwise, to he able to sample from broad range of distributions with V log orthogonal to the boundary, we need that for all p, there exists Gp. and such that (v = p-ppp- Gp + log p. Applying the projection operator on both sides, we have
(F&(v) P(p GG + V log logp.
(17)
Moreover, the function space { admits a set of basis functions {eih ie. F6 = E, 6ee, then
v log PFa c((.P..)
(18)
This means that the space of V log needs to be spanned by {Pei), so { should also be spanned by an infinite number of basis functions {e,), ; if we want sample from complex class of V log p.
B Neural dynamics not following gradient fields
For the meaning of notations, see Section 3.2. Proposition 5. When > 2, 1SCC (v) D(-v wo(v) + is gradient field if and only if W a zero matrix.
Proof.
If (v) VG(v) is gradient field. then
a FSC (v);
a
PC (v), =
G(v)
(19)
Fr.
Fu,du,
16

---

because (v) is smooth, and we can switch the order of partial differentiation. Written in the form of (7), we have for all vs and bi (where # j).
a
( -dju + d,
d + d, Σ Wad(vn) + +d,1, Bwy -d,e,
Σ W W)L((os) + tdzzs) km
dv,
d W (w))) Joj
do: d, (i(((())))
d &,,, = (
since u, and Vj are arbitrary and is nonlinear (so the derivative d,, is not the same for all v), this implies di aj = dy W = 0. But the decay constants cannot be O. Therefore, Wiz = W =
(=)
If W a zero matrix, then F (v) = VD (- v.v + +. is gradient field.
B.1
Hilbert projection theorem and proof of Proposition 2
Here we directly state the Hilbert Projection Theorem. For detailed proof, we refer readers to standard functional analysis textbooks [44].
Theorem 6. (Hilbert Projection Theorem) Given a Hilbert space H, and a nonempty closed convex
set E H. for any E H, there exists unique II E E, such that |[ z|l = inf || - cll-
Here II is called the orthogonal projection operator. Next, we give the proof of Proposition 2.
Proof. Here we consider special case of a Hilbert space L2(R",, R""; :R."":) where the inner product defined by (f.g)
Jem If 9) p(x)xx and we let E {Av + Bé(v) [|A, B Rm E [_xi]. The set functions E is a closed convex vector sub- space L2 : because closed under addition and dinee--imeniional. Hence by Hilbert projection theorem, there exists unique ITH E E such that [[ - IT) | = inl H Fl > 0.
Because {FSC), 10 C if the decay matrix D is diagonal, we have
inf inf
H
SC
2 inf | F = [1 H >
(20)
C Proof of Theorem 3
First, we give statement of the universal approximation theorem as in [33].
Theorem C.1. (Universal Approximation Theorem) The continuous transfer function is not polynomial if and only for every N, m N, compact set U RRt,, C(U, Rm), and there exists N. W2 Baxk b R. W1 E Rmxn such that
sup Il 9(x))l(l1]110)0) <
(21)
=EU
where g(x) = W1d(W2r + b). Or equivalently, Σn = span (Ø(w: x b) w R dense in C(R)) if and only if & is not polynomial.
b R.
is
The theorem says that for any transfer function that is not polynomial (tanh for example), there is a multi-layer neural network with the transfer function that able to approximate any continuous function to arbitrary precision on closed and bounded subset of R The proof of Theorem proceeds in three steps:
1. Choose a large enough compact support U for so(x) such that
p(x) logp(() s0(x) (x)ll dx = -ammp p(x) V dx
(22)
Bm
2. Choose large enough n and appropriate 0 such that the universal approximation theorem can be applied to show that
[ p(x) V Vlogp(x) - so(x)]1 dx < 8
(23)
17

---

3. Combine equation (22) and (23) so that we have
(xx log p(x) - s0(x)]1 {
p(x) V log p(x) - so(x))1 <x
<E
(24)
We begin with two lemmas that will be helpful in the proof of the theorem. Both of them standard results but we include here for completeness. Let the upper incomplete Gamma function be defined by [(s, C) -c y eddr.
Lemma C.2. If's is a positive integer, then I'ss, C) = 1)le
Li=o
C
Proof.
This can be proved by induction. For = 1,
-
(25)
Therefore the conclusion holds. For >> 1, we assume that for = k. the conclusion holds. Then for k + 1, Through integration by parts, have
T(h 1,,+ C) C
(26)
Therefore
r(k +1,0 = k(k 1)le
C
+
C
i!
1=0
C
C
= kee
(27)
k!
i!
0=
O
kle
i!
==0
Therefore the conclusion holds 8 = k:
as well. By induction, holds for all positive integer
Lemma C.3. C)s, 0 as C
Proof. When C + xc, ['(s, C) < le -l Therefore W.L.O.G., we can assume that 1 is a positive integer. By lemma C.2,
C
lim C(s,C)= (8 1)! C-
lim
(28)
C
i!
1=0
since the exponential term e C beyond all orders.
Now give formal proof of Theorem 3.
Theorem. Suppose that are given probability distribution with continuously differentiable density function p(x) : Rm + RR+ and score function log p(x) for which there exist constants M1, M2, a, k >> such that
p(x) < Mie a||]||
(29)
V log p(x) < M2 ||x||
(30)
when ||x]| > for large enough L. Then for any e >> O, there exists a recurrent neural network whose firing-rate dynamics given by (11), whose recurrent weights, output weights and the diffusion coefficient are given by Wiee E Raxx of rank m, Wont E Rmix and a E Raxm respectively. such that, for large enough n, the score of the stationary distribution of the output units se (x) satisfies Exx<p(xill logp(x) - sp (x) <
see http:://dlm..nstt.ov......
18

---

Proof.
Step 1 First of all, we would like to find a compact set Uc = [x||| I C} such that fgmque ||V log)(x) sx)]l] p(x)dx < We can upper bound the integral by the following calculation using spherical coordinates when C > L.
J Jan Vlogp(xx> - sx(x]] p(x)dxx=
/ JRmAUe
VVlogp())) p(x)dx
(31)
| JRTAUU
M2 a||x||
|xx|| dx
(32)
m E e-r
(33)
((
2M, m/
m e =dddu
(34)
am+
JaC
2M, M2x
I'm + k. aC)
(35)
amthr(
Note that we assume sx(x) = O outside of UC 4. By lemma C.3, there exists large enough constant Ö such that I'(m + k, a,a <
m TY [())s
1M, M2 Therefore
2M, 7/2 m/2
am*kT(m)) m
|V so(x)]l p(x) <
(36)
RmAUa
a=*)
42
2
Step 2 Next we write down the specific form of sx(x), the score of the stationary distribution of x for the reservoir-sampler dynamics (11)
s0(x) (x) = (WWonooo "'''' ( -2 (-ax + (Wrea
rec + I)).
(37)
Since a is learnable, we can define = WT. (W WoT.) such that W oo WT. identity matrix 5:
out is
W Wonoo "WT out Wout WT (W out WT. ) - (Won WT - Wout WT
out
(38)
As p is continuously differentiable, by Theorem C.I, there exists W Wout, W such that Vlogp(x)
Wout (onneeWWeece x -(x+
sup x
(39)
2a
Sa
Therefore with those choices of Wout. Wree and 0,, we have
|ogg((x sx(x)]l p(x)dx < sup xEU
VVoog(( -so(x)s))))) (x))
V logp(x) 2(-axx + aW oWWon (exx 1))
sup xEU
V logg((x)
4a²
sup
+x Wout W rec + 1)
xEUa
2a
4a²
8a.
2
(40)
'Equality (31) is the only place where this assumption is used. The equality also holds if sx) grows at similar rate logg.. p.
"A technical subtlety is that Wout W is not necessarily invertible. However, since every transformation considered here is continuous, we could add small perturbation to Wost such that Wont W out is invertible, and result still holds.
19

---

Step 3
+ -(mm
(((e) p(x) VV logp(x) - se(xx]] dx
- s0(
Rmmwa
2 + E
=E
2
D Pseudocode for training the reservoir-sampler network
Algorithm 1:
Training RSN
Input: Training samples {x:), pseudo recurrent weights W output weights Wout, start noise level øi, noise decay factor C, number of noise level N Output: trained recurrent weights WTe. trained output weights Wout, trained diffussion coefficient
for i 1....,, .N do
while not done do
Sample from the training set
Perturb the samples with Gaussian noise x x, +E, N(0,0
Compute the gradient of the score-mase matching loss
V Wines
|| 2ou re 1) -)x) x)
Update Wrec, Wout with gradient descent
0++1 =Co;
W Wrec Wout a f- W (W oot WoL)
out
return W rec, Wont
E Accelerated sampling with reservoir-sampler network
To construct an irreversible dynamics that accelerates sampling, we let the divergence-free (DF) field be of the form JVp, where = -J" can be any skew-symmetric matrix that is either learned through training or prescribed. Then from equation (6) we see that the drift term Fa satisfies
Fo(v) = VV log p(v) + JV log = + JJV log p(v).
(41)
Therefore the corresponding score matching problem is axE + JJ (W @W rec + x) R log p(x). With a slight modification, the pseudo code for an accelerated RSN training algorithm is given in algorithm 2. We repeat the experiment in Section 4.1 with the accelerated RSN. The results are shown in Figure 5, We see that the accelerated RSN is able to sample from the distribution faster than the RSN.
F
Additional experiments
F1
Failure cases for sampler-only networks with hyperbolic tangent nonlinearity
Figure shows that with tanh nonlinearity, sampler-only network able to sample from the distri- bution we consider in Section 4.1. However, this is because the score function of the distribution considered can be spanned by just two basis functions equipped by the drift term of the sampler-only network. Therefore, our theory predicts that the sampler-only network will not be able to approximate score functions that are not exactly spanned by these two basis functions. Indeed, if we consider slightly different Gaussian mixture distribution, Pdata (z) = (W1,,, 0.0625) WN(1, 0.25)), then the sampler-only network will not able recover the score function as shown in Figure 7, while RSN still succeeds in recovering the score function.
20

---

Algorithm 2: Training accelerated RSN
Input: Training samples {x pseudo recurrent weights Wee. output weights W Wott pseudo diffusion coefficient o, skew-symmetric matrix J, start noise level 01, noise decay factor C, number of noise level N
Output: trained recurrent weights Wree trained output weights Won., trained diffussion coefficient
for = =....., NN do while not done do
Sample x from the training set Perturb the samples with Gaussian noise x x; + ~~ ~N(10,07 Compute the gradient of the score-matching loss
1 W... a + (Won + 1) x) x-
Update Wree, Wont; with gradient descent Update J with gradient descent (Optional)
0++1= =C,;
W WrecWoor out a g W (W out Wol! - Ō
return W rec; Wont
X1
X2
1.8
1.75
Accelerated
Accelerated
Wasserstiin distance
1.6
Langevin
Wasserstein distance
1.50
Lengevin
1.4
1.25
1.2
1.00
1.0
0.75
0.8
0.6
0.50
0
50
100
150
200
50
100
150
200
time (ms)
time (ms)
Figure 5: Sampling speed of accelerated RSN We repeat the experiment described in Section 4. where we try to sample from 2-D Laplacian mixture distribution. The figure shows the Wasserstein distance between the marginalized distributions of generated samples and true samples for both dimensions. We see that compared to vanilla RSN, accelerated RSN is able to sample from the stationary distribution with faster relaxation time.
21

---

SO-SC
SO-FR
RS-FR
b)
C)
100
true
5.0
5.0
5
2.5
2.5
Noise
score
0.0
0.0
10-1
2.5
22.5
-5
-5.0
-5.0
d)
2
e)
-2
f)
-2
# sampled points
200
200
200
Scaled density Sampled points
150
150
100
100
100
50
-2
2
-2
0
2
V
r
X
Figure 6: Sampler-only networks with hyperbolic tangent nonlinearity are able to sample from a I-D Gaussian mixture distribution with a symmetric score function.. See Figure for an explanation of the abbreviations. Here we use tanh nonlinearity for sampler-only networks (SO-SC and SO-FR). All types of networks are able to match the score function and sample from the distribution.
SO-SC
SO-FR
RS-FR
a)
b)
C)
100
20
true score
true score
score
score
10
10
10
Noise
0
10-1
-10
-10
2
2
-2
0
2
d)
f)
400
400
400
Scaled density
Sampled points
# sampled points
300
200
200
200
100
100
100
-2
0
2
-2
-2
0
2
Figure 7: Sampler-only network with hyperbolic tangent nonlinearity fails to sample from 1-D Gaussian mixture distribution with a non-symmetric score function. See Figure 2 for an explanation of the abbreviations. Here we use tanh nonlinearity for sampler-only networks (SO-SC and SO-FR). Only RSN is able to match the score function and sample from the distribution.
22

---

F.2 Sampling using autoencoder
We show that we can replace the fixed linear PCA projection used in Section 4.2 with pretrained nonlinear autoencoder. We test the algorithm on both the MNIST dataset (Figure 8) and the CIFAR-10 dataset (Figure 9).
3
6
8
5
8
8
8
6
8
7
3
9
8
3
6
6
6
7
S
6
8
8
6
9
6
8
5
8
8
8
8
8
6
8
6
8
9
8 0
6
B
9
6
8
6
9
4
4
8
6
6
9
3
7
9
8
9
8
8
4
8
6
3
8
8
8
6
8
3
5
8
0
4
8
8
3
4
8
2
9
3
6
8
8
6
9
9
8
8
0
5
8
8
6
8
Figure 8: Samples produced by our reservoir sampler network when trained to fit the latent distribution produced by non-linear convolution autoencoder from MNIST. The distribution 32-dimensional. The reservoir sampler network had 500 reservoir neurons.
G Related works
In this section, we examine the neural dynamics considered [17, 1, 38, 19] and discuss how our work is related.
Dong al. [17] uses Hamiltonian Langevin dynamics accelerate sampling in continuous attractor networks. The Hamiltonian dynamics below are considered in the paper (cf. eq (19) and (20)):
Ta dt
ay
(42)
dy = ---a-11- + A(s" -) VTa
de
The overdamped Langevin dynamics that samples from the same marginal stationary distribution over as (42):
di = A(s - 8) + 22c.
(43)
where s is an observation generated by a latent feature s according to a Gaussian distribution, and the linear term Ass s) comes from the external input. So the network dynamics in Dong et al. [17] themselves do not generate samples without the external input, and the study does not clarify how to supply the network with nonlincar score functions.
23

---

CIFAR 10 Images
Reconstructed Images
Samples
Figure 9: Samples produced by our reservoir sampler network when trained to fit the latent distribution produced by non-linear convolution autoencoder from CIFAR-10. We show both the original images and the reconstructed images to evaluate the sampling quality. The distribution was 128 28-dimensional. The reservoir sampler network had 20000 reservoir neurons.
Aitchison and Lengyel [1] have shown that Hamiltonian Monte Carlo (HMC) naturally maps onto the dynamics of excitatory-inhiiitrry neural networks. However, the neural dynamics considered entirely linear in terms of the network variables (cf. eq (23) and (24)). Therefore, the dynamics are again equivalent to the synaptic current dynamics with an identity transfer function Section 3.2. The same applies the analysis in Masset et al. [38], which also only considers linear neural dynamics. By Proposition 2 and the decomposition discussed in Kwon et al. [31], these dynamics only sample from probability distributions whose score functions are linear, i.e. Gaussians.
The work by Echeveste et al. [19] provides a unifying model for several dynamical phenomena in sensory cortices by training nonlinear recurrent excitatory-inbibioory neural circuit model to perform sampling-bass probabilistic inference. The synaptic neural dynamics below (cf. equation (8)) is used, where the nonlinear transfer function is supralinear, i.e. o(u) = =[u]""
N
du; dl
U;
+
Wi ( ())))
(44)
As discussed in Section 3.2, while they are richer than linear dynamics, dynamics of this form still have limited ability to sample from complex probability distributions. Moreover, we note that the training criterion considered Echeveste et al. [19] only matches the mean and the variance of the stationary distribution. These points are in contrast to our work here, where we show that the reservoir-sampler network, which extends the dynamics (44), can be trained via alternate methods to sample from complex probability distributions by matching their score functions.

---

H Details of numerical experiments
All experiments were run on one NVIDIA Quadro RTX 6000 GPU. We always use the model saved at the last training iteration.
H.1 Training details for Section 4.1
For the 1-D bimodal Gaussian distribution experiment (Figure 2), used 1000 reservoir neurons
(only in the reservoir sampler network) and sampler neuron. We trained all 3 networks (SO-FR, SO-SC, RS-FR) for 79000 iterations with the Adam optimizer. The learning rate was 0.0001, and the batch size was 128. We started by adding Gaussian noise with standard deviation (std) 01 to the training samples. We decreased the std every 7900 iterations by constant factor, i.c. 0++1/0 C, such that 0'0 .01. When simulating the stochastic neural dynamics, we used step size of 10 10-4, and ran 10000 steps.
For the mixture 2-D Laplace distributions experiment (Figure 3), we used 1000 reservoir neurons (only in the reservoir sampler network) and sampler neurons. We trained the RS-FR network for 62,500 iterations with the Adam optimizer. The learning rate was 0.0001, and the batch size was 128. We started by adding Gaussian noise with std σ1 to the training samples. We decreased the std every 6250 iterations by a constant factor, i.c. 0.+1/o; = C, such that 010 = ,01. When simulating the stochastic neural dynamics, used step size of 10 4 and ran 20000 steps.
H.2 Training details for Section 4.2
For the MNIST experiment (Figure 2), we used 20000 reservoir neurons (only in the reservoir sampler network) and 300 sampler neurons. We trained all 3 networks (SO-FR, SO-SC, RS-FR) for 1000 epochs with the Adam optimizer. The learning rate was 0.001, and the batch size was 64, We started by adding Gaussian noise with std σ1 1 to the training samples. We decreased the std every 100 epochs by a constant factor, i.e. 0.41/a; C, such that = .01. When simulating the stochastic neural dynamics, we used step size of 10 and ran 30000 steps.
I Statement of code availability
All code is available on Github
25


Reviews:
[{"review_id": "suOiwr1zzB", "reviewer": "Reviewer_sqZn", "summary": "The paper investigates architectural requirements for recurrent neural circuits to sample from complex distributions using diffusion models. It presents a model where traditional sampler-only networks are enhanced with additional firing-rate dynamics and a set of separate output units, called reservoir-sampler networks (RSNs). An efficient training procedure based on denoising score matching is proposed. Empirical experiments are presented to demonstrate the model's ability to sample from complex data distributions.", "questions": "See weaknesses.", "limitations": "I am satisfied with the discussion of limitations.", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "\n- The research addresses a significant issue in computational neuroscience and Bayesian learning – how neural dynamics can sample a complex distribution. \n- Biological plausibility and mathematical computation efficiency are particularly commendable aspects of this work.\n- The theoretical analysis is well constructed, providing strong support for the model design, including the choice of firing-rate dynamics and the reservoir-sampler network.\n- The empirical experiments presented provide evidence of the effectiveness of the proposed method.", "weaknesses": "- The experiments conducted need improvements. Despite using score matching for model training, the quality of the generated images leaves room for improvement. The incorporation of Unet-like or transformer backbones might have enhanced the results.\n- The paper lacks a clear, diagrammatic representation to explain the utilization of diffusion models in the design of recurrent neural circuits architecture.\n- The results and conclusions seems to have limited impact to the general AI/ML areas, giving limited experimental results and comparison (e.g. precision and efficiency ) with existing diffusion probabilistic models is somewhat lacking.  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the comments and suggestions. Below we address the reviewer’s concerns:\n\nWeakness:\n- Image quality is poor: It is an accurate point that the fidelity of the generated images is substantially lower than that would be realized by, for example, a U-net.  However, our goal here is to study probabilistic sampling in the class of RNNs that is broadly used as mechanistic and functional models in neurobiology.  Specifically, we ask: 1) whether the RNN that are considered in the current literature can sample from complex distributions, and 2) If not, what form of RNN is able to achieve this? While we believe our paper gives significant answers to both questions, ur, we fully acknowledge the subpar performance of using the RSN to generate complex images in contrast to U-nets.  We will acknowledge this explicitly in our revision and point the reader to the discovery of biological components that enable better score matching as an important area of further research. \n\n- We thank the reviewer for this important comment.   We wish to underscore that our approach, while partially inspired by the diffusion model, also holds differences. While the diffusion model can be conceptualized as a time-inhomogeneous SDE with a finite time horizon, our work delves into a time-homogeneous SDE (the neural dynamics) with an infinite time horizon. The entire training and sampling pipeline is diagrammatically described in Figure 1 and Figure 4a). Recognizing the need to make this more clear, in our revision, we will (1) more clearly refer to this diagram in the main text, and (2) add to the diagram a clear set of pseudocode.  We are confident that these two changes will both more clearly explain our methodology and illustrate to the readers the similarities and differences with the classical diffusion approach. \n\n- We appreciate that our emphasis is on implications for neuroscience, and as such we have limited implications for the cutting-edge AI/ML area in the short term. However, we do note that reviewers are starting to look at how sampling/diffusion can play a role in RL, robotics, and modeling human behavior [1-4]. Our paper explores the neural underpinnings of those applications and has the potential to inspire more AI/ML work in this regard. \n\n[1] Ajay, Anurag, Yilun Du, Abhi Gupta, Joshua Tenenbaum, Tommi Jaakkola, and Pulkit Agrawal. “Is Conditional Generative Modeling All You Need for Decision-Making?” arXiv, July 10, 2023. http://arxiv.org/abs/2211.15657.\n\n[2] Chi, Cheng, Siyuan Feng, Yilun Du, Zhenjia Xu, Eric Cousineau, Benjamin Burchfiel, and Shuran Song. “Diffusion Policy: Visuomotor Policy Learning via Action Diffusion.” arXiv, June 1, 2023. http://arxiv.org/abs/2303.04137.\n\n[3] Janner, Michael, Yilun Du, Joshua B. Tenenbaum, and Sergey Levine. “Planning with Diffusion for Flexible Behavior Synthesis.” arXiv, December 20, 2022. https://doi.org/10.48550/arXiv.2205.09991.\n\n[4] Pearce, Tim, Tabish Rashid, Anssi Kanervisto, Dave Bignell, Mingfei Sun, Raluca Georgescu, Sergio Valcarcel Macua, et al. “Imitating Human Behaviour with Diffusion Models.” arXiv, March 3, 2023. http://arxiv.org/abs/2301.10677.\n"}, {"review_id": "EVFK999q4Q", "reviewer": "Reviewer_Y31x", "summary": "The authors proposed a revervoir-sampler network whose firing rate dynamics can sample from an arbitrary probability distribution. They first established the relationship between the sampling power of the neural dynamics and the ability of the dynamics to approximate the score function. Then they showed that the synaptic current dynamics of the traditional sampler-only networks is only able to approximate score functions that are in a finite-dimensional function space. Importantly they proved that the  revervoir-sampler network can sample from an arbitrary probability distribution to arbitrary precision. Finally, they developped a biologically-plausible learning rule for the proposed model.", "questions": "Q1: Why did the authors choose to use the reservoir network as the building block? \n\nQ2: Following the first question, does the RSN model work because of the linear readout layer, which rescales the distributions? Can the authors explaine more about this？\n\nQ3: Can the current model be extended to the multiplicative noise scenario which could be the feature of the poisson spike neurons, or (at the network level) be the feature of the E-I balanced spiking net.\n\nQ4: How does the RSN model compare to the coupled attractor model which also implement Langevin sampling [1]?\n\n[1]: Zhang, W. H., Lee, T. S., Doiron, B., & Wu, S. (2020). Distributed sampling-based Bayesian inference in coupled neural circuits. bioRxiv, 2020-07.\n", "limitations": "The authors addressed some limitations.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "This paper is very clearly presented and the scientific question is quite interesting. ", "weaknesses": "Comparison with previous models should be addressed in the experiment part.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer for the thorough summary of our paper, as well as the comments and suggestions. \n\nWeakness:  We appreciate the importance of stronger comparisons with previous models, and in our revision, we will significantly expand our treatment of these in the main text. In the current version of the paper, the sampler-only architecture serves the role of previous models, and we discussed how previous works can be seen as special cases of either SO-SC or SO-FR model in Appendix D. In the revised version, we will include and expand the contents of Appendix D in the main text.\n\nQuestions:\n\nQ1: Although the brain exhibits a hierarchical structure (e.g., the visual hierarchy), within each hierarchy neurons are recurrently connected, i.e. they may form a reservoir. This is our motivation for using the reservoir network as a basic building block of our model.\n\nQ2: Intuitively, the linear readout layer linearly combines different basis functions, and the reservoir is responsible for supplying a sufficient number of basis functions. Therefore the RSN model works because of 1) a large enough reservoir and 2) the linear readout layer. \n\nQ3: Thank you for this great question.  Yes, with a slight modification of the Fokker-Planck equation in the paper. If we take the diffusion coefficient to be $\\sigma(x) = \\sigma \\sqrt{\\operatorname{diag}(x)}$ and the diffusion matrix $\\Sigma = \\frac{1}{2} \\sigma(x) \\sigma(x)^T = \\frac{1}{2} \\sigma^2 \\operatorname{diag}(x) $, it can be shown that the stationary distribution $p$ satisfies $\\nabla \\cdot(\\frac{1}{2} \\sigma^2 p + \\Sigma \\nabla p - p F) = 0$ where $F$ is the drift term of the neural dynamics. Again if we take the divergence field to be 0 and divide both sides by $p$, then the corresponding score matching problem would be $\\nabla \\log p \\approx \\Sigma^{-1} (F - \\frac{1}{2} \\sigma^2)$. \n\nQ4: Zhang et al. propose that coupled CANNs can implement Langevin sampling on the space of latent stimulus features. To see the connection between their work and ours, it is helpful to consider the functional form of the Langevin dynamics used by these authors (cf. eq. 4 in their paper). The score function is set directly by means of a formula to obtain the desired target distribution, and since only Gaussian distributions are considered, the stimulus feature needs only be linearly transformed (recall that the score function of Gaussian is linear). Our work, by contrast, (1) treats general target distributions, beyond the Gaussian case, (2) proposes a way to iteratively learn the score function by just looking at the samples from the target distribution, and (3) specifies conditions in which the score functions corresponding to general target distributions can and cannot be successfully learned. \n"}, {"review_id": "G4n2IyRG5k", "reviewer": "Reviewer_wwgi", "summary": "- This study solves an outstanding problem involving arbitrary density sampling using stochastic neural networks.\n- The authors propose a a Reservoir Sampling architecture, whereby an auxiliary recurrently-connected population facilitates the sampling from a non-trivial distribution.\n- The work seems solid and elegant (although SDEs are not my domain of expertise), and the numerical experiments are convincing.\n\nUPDATE: Sep 1, 2023. I have read the rebuttal, it addressed my Qs and I maintain my score.", "questions": "- This auxiliary RNN population approach reminds me of a recent poster/paper at cosyne by Duong et al. (\"Adaptive whitening in neural populations with gain-modulating interneurons\"; ICML 2023 https://openreview.net/forum?id=cEWB5hABV5), where they show that if the auxiliary population is greater than a specific size, it can represent and whiten a multivariate gaussian densities exactly. In your framework, is there a relationship between the dimensionality of reservoir and the expressivity of the network in terms of what kinds of densities it can represent? Could this provide testable predictions for an experimentalist to search for? \n- Suggestion: I did enjoy reading this manuscript, but there was more emphasis on technical parts than necessary to get your point across. If the text were reframed with less emphasis on the math, I think it would flow much nicer and would be much more approachable.", "limitations": "There was adequate discussion of limitations.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The math looks solid, although I did not have time to go through the supplement.\n- There are numerical experiments on complex distributions (MNIST PCs), demonstrating that their procedure can operate on non-trivial densities.", "weaknesses": "- This work seems to address an _existence_ problem with arbitrary density sampling with RNNs, but how would one validate this kind of model in physiological experiments?\n- In fact, I would appreciate more discussion on the state of affairs in the neurophysiology of density representation.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer’s critiques and suggestions. Our point-by-point replies follow, together with changes we will make to the revision in light of the reviewer’s points. \n\nWeaknesses:\n\n- Validation in physiology experiments: This is a great question and one which we now see we should have addressed fully in our submission; our revised manuscript will add to the discussion to improve the paper in this way.  In particular, we will describe the following general methodology: A first step in analyzing physiological data in our framework is to identify windows where the distribution of recorded neural activity is roughly stationary over time.  The next step is to train different network models using our framework to check if these models are capable of generating samples from the same stationary distribution as measured. Validation would involve not only asking whether this procedure succeeds, but evaluating the size of the “unobserved” reservoir network (i.e., number of hidden neurons) that is necessary, as well as whether different features of single-neuron physiology (modifications to the single-neuron and coupling terms that define the network) and priors on network connectivity lead to the better fitting of the experimental distribution.  We note that the outcome here is an answer to the question of what the conditions are on network size, together with the set of weights/other physiological parameters, that can reproduce the distribution of physiological data.  We do caution, however, that it will not be a complete answer, as failure to match these data could also arise from limitations in the underlying learning algorithm and/or optimization process.\n\n- Regarding the state of affairs in the neurophysiology of density representation:  We agree that providing more here would improve our paper, and in our revision, we will add text to the introduction and discussion to relay at least the following:  Strong evidence shows that both humans and animals use some form of uncertainty information to guide their behavior; moreover, multiple brain areas have been identified as neural correlates of this uncertainty [2]. Neurophysiology data and analyses that directly verify how neural activity represents probability densities remain in their early stages. However, neural sampling theories, to which our work contributes, have made empirical predictions of certain physiological properties such as noise correlations and Fano factors, and these have been verified in, e.g., [1]. We see our work as opening the door to an allied study of higher-order statistics that more fully describe complete probability distributions. \n\nQuestions:\n\n- The higher the dimensionality is, the closer the stationary distribution of the sampler is to the target distribution because in this case, one has more basis functions at one’s disposal to match the score function. For testable predictions, a hallmark of whether enough neurons are considered is whether the distribution can indeed be captured; if not, our work suggests that might require a reservoir of neurons that supply input signals. Please also see further experimental implications under the response to your first “weakness” point above.\n\n- We are glad that the reviewer enjoyed reading our work, and appreciate the technical nature of much of our presentation as well.  For the revision, we will go through the paper top to bottom, and identify (1) at least several places where technical details can be moved to the supplement, and (2) clear topic sentences that can be added in less technical language that describe the import of a mathematical argument that follows, with the latter set off by a phrase such as “The technical details of this are as follows,” to accomplish the same objective.  \n\n[1] Orbán G, Berkes P, Fiser J, Lengyel M. Neural Variability and Sampling-Based Probabilistic Representations in the Visual Cortex. Neuron. 2016 Oct 19;92(2):530-543. doi: 10.1016/j.neuron.2016.09.038. PMID: 27764674; PMCID: PMC5077700.\n\n[2] Rullán Buxó, Camille, and Cristina Savin. \"A sampling-based circuit for optimal decision making.\" Advances in Neural Information Processing Systems 34 (2021): 14163-14175.\n\n"}, {"review_id": "GgiZY5K9gj", "reviewer": "Reviewer_kc8J", "summary": "In this paper, the authors address the question of how the dynamics of recurrent neural networks (RNNs) can generate samples from probability distributions of interest. This is of interest to the neuroscience community, since it has been proposed that biological networks can use such sampling-based inference to represent and compute with probability distributions in order to account for their noisy surroundings. The authors prove that standard RNNs have limited expressivity when used directly for sampling-based inference.\nThey proceed to show that RNNs with separate (linear) readouts in contrast have the necessary expressivity to approximate any probability distribution through Langevin sampling. Finally, it is demonstrated empirically that such a 'reservoir' formulation allows RNNs to generate samples from a variety of complex distributions in practice, including multimodal, heavy-tailed, and high-dimensional distributions.", "questions": "L60: The authors state that they derive a 'biologically plausible' learning rule, but the training procedure described in Eq 14 does not seem particularly biologically plausible, requiring both backprop and matrix inverses (as the authors also highlight themselves in L288-291). Perhaps the introduction could be reworded to make this initial statement more in tune with the subsequent methods/results?\n\nL73: it would be useful to briefly state what these regularity conditions entail.\n\nL114: Helmholtz typo\n\nFigure 2 legend: \"SO-SC and SO-FR is (_sic_) only able to fit the score function with piecewise linear function (_sic_) when using ReLU transfer function (_sic_)\". Why the contrast between 'piecewise linear' in the legend and what looks like globally linear approximations to the score function in the figure?\n\nRef 5 & 37: formatting errors in author name (\"GergHo Orban\").\n\nI am slightly confused by the authors' network dynamic equations (Eq 7 & Eq 8). When modeling neural firing rates, a positive (e.g. ReLU) nonlinearity is commonly used since firing rates are strictly positive. However, would the Brownian noise term added to the rates in Eq 8 not allow for negative neural firing rates?\n\nI was also surprised by the synaptic current equations, since RNN dynamics are usually formulated in terms of either firing rates or membrane potentials. Potentials and currents are of course equivalent for ohmic units, but I was wondering if there was a reason that the authors decided to describe Eq 7 in terms of currents rather than membrane potentials?\n\nThe authors have discussed the case of a non-linear RNN with a linear readout and shown that it can approximate arbitrary probability distributions. Is the same true for non-linear RNNs with non-linear readouts? and if so, could the reservoir sampler network be considered a special case of a non-linear RNN where only a subset of units are used to represent samples from the target distribution, and the remainder are used for additional computational power/expressivity?", "limitations": "The authors have adequately addressed the limitations of their study.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "The paper contains several interesting theoretical results related to the topic of sampling-based inference using recurrent neural networks. This is an important open problem in computational neuroscience, where it remains unknown how the brain performs (approximate) probabilistic computations and represents the necessary probability distributions for such computations. By providing a firmer theoretical grounding for our understanding of sampling-based inference, this work has the potential to inspire both new computational and experiment studies of probabilistic computations in neural circuits. The paper is also generally well written and demonstrates impressive empirical results, although only on relatively 'toy' problems.", "weaknesses": "Much of the motivation and discussion of the paper is written with reference to the neuroscience literature on sampling-based inference, especially in sensory cortices. However, the empirical results all use artificial toy examples with no obvious relation to questions of interest to the neuroscience community. It would have been interesting to e.g. look at natural images or similar and delve a bit deeper into the dynamics learned by the network to make it more interesting to a neuroscience audience.\n\nMany of the figure axis labels and legends are too small to be easily readable, and they would benefit from an increased font size.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer’s recognition of the significance of our work and the comments and suggestions. Below we address the reviewer’s concerns and questions.\n\nWeaknesses:\n\n- Regarding the application of our model to natural images:\nWe thank the reviewer for bringing up this point. We agree that training our reservoir sampler on natural image patches would be more relevant for neuroscience applications. We plan to train the model on the CIFAR 10 dataset and evaluate the Fréchet inception distance\n(FID) of the model in the revised version.\n\n- We thank the reviewer for the accurate critique regarding our figure formatting and will fix the axis labels and legends in the revised version. \n\nQuestions:\n\n- L60: The reviewer is absolutely correct that there is still considerable distance between our proposed algorithm and a biological implementation, and our revision will clearly note this.  Our meaning here was indeed more narrow than we now see we expressed, in that the learning rule “sidesteps the demands of backpropagation through time (BPTT)”.  We will carefully rephrase our usage of “biologically plausible” to make sure this qualification is clear throughout our revision.  \n\n- L73: please see Cerrai, 2001 Hypothesis 2.1 & 2.2. Roughly this indicates that 1) the drift term b(x) does not increase too fast (slower than exponential rate) as x goes to infinity, 2) the drift term is locally Lipschitz, and 3) the diffusion coefficient is invertible as we assumed in our work. \n\n- L114: the typo will be fixed in the camera-ready version.\n\n- Figure 2 legend: We appreciate this astute point and the opportunity to clarify here; in addition, our revision will include a revised version of this figure that does not have the confusing linearity of the score function (please see Fig S1).  In reply, we first recall that the solutions are found by minimizing the denoising score-matching loss. Because we are using a distribution whose score function is symmetric with respect to the origin, it makes sense for the solution to have the same property.  With a piecewise linear activation function, SO-FR and SO-SC have the ability to “bend” the score function at a single point, but due to this symmetry they will not benefit from this in terms of the score matching loss, hence the seemingly “global linear approximation”. We believe that the new  Figure S1, which essentially conducts the same experiment using tanh nonlinearity instead of the ReLU function, presents the underlying situation much more clearly, and thank the reviewer for inspiring it.\n\n- Ref 5 & 37: will be fixed in the camera-ready version\n\n- Noise term can produce negative values within firing rate dynamics: We thank the reviewer for another on-target and careful point.  Yes indeed, this is true, and an imprecision (or flaw) that comes along with our basic SDE formulation of the rate-RNN on an unrestricted domain. In our revision, we will carefully note this.  We will also state that one quick remedy would be to add a high energy barrier at 0, i.e. add $C \\cdot \\operatorname{ReLU}(-x)$ to the score function where $C$ is large compared to the noise magnitude so that the dynamics would almost never output negative values. Since this only adds one more basis function, it would not affect our main argument.\n\n- Nomenclature of synaptic current dynamics: We absolutely agree that this terminology can be confusing!  After surveying a number of sources in the literature without identifying a clear consensus, we simply settled with the present choice as it follows the convention in the seminal text of Abbott & Dayan 2001 “Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems” (cf. eq. 7.39). \n\n- What about non-linear readouts: If the nonlinear function following the linear layer is invertible (e.g. tanh), then the result still holds, and the density function would be scaled by the determinant of the jacobian. Otherwise, the result may not hold. For example, if we just apply ReLU to the sampler neurons, then we can only obtain positive samples. And yes, if we ignore the transients of the sampler neurons, then we can say that RSN is a special case of a non-linear RNN where only a subset of units are used to represent samples from the target distribution. \n"}, {"review_id": "AYrNONrlBo", "reviewer": "Reviewer_pX3M", "summary": "This paper studies circuit algorithms for sampling-based Bayesian inference in continuous-time rate-based recurrent neural networks. Specifically, it argues that using a linear readout from a noisy reservoir provides a substantial expressivity benefit relative to representing the sampled distribution directly in the recurrent population. ", "questions": "- The figure axis labels are too small, making them hard to read. Figure 4b is also too small. Please increase the font and panel sizes to improve legibility. \n\n- Lines 32-33: Along with the other references on Gaussian sampling, Hennequin, Aitchison, and Lengyel, \"Fast Sampling-Based Inference in Balanced Neuronal Networks\" (NeurIPS 2014), should be cited. \n\n- Line 225: I do not see how the proof of Theorem 3 is \"constructive,\" since it is based on an existence-type universal approximation theorem. \n\n- Lines 303-304: \"This procedure is more aligned with the developmental processes involved in forming visual representations in the infant brain, where visual representations are thought to be noisier (less linearly separable) initially [3]\" Please elaborate; I don't quite follow. \n\n- Line 317-318: This is quite vague. Having a clearer proposal for what this model might help explain in neuroscience would make for a much stronger paper. \n\n- It would be helpful to repeat the theorem statements in the Supplement. \n\n- Lemmas C.2 and C.3 can be replaced by using standard results on the asymptotic expansion of the upper incomplete Gamma function for large argument, see https://dlmf.nist.gov/8.11. \n\n- In Line 562, in the proof of Theorem 3, should $m$ be replaced by $n$?", "limitations": "The authors provide a generally adequate discussion of the limitations of their work (modulo the weaknesses discussed above), but fail to show the data for some points (see point 3 under Weaknesses). ", "rating": 5, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. The question of the minimal recurrent circuit architectures that allow sampling from complex distributions is clearly relevant to the neuroscience community at NeurIPS, and I think the main result is interesting (albeit unsurprising).\n\n2. The experiments show to a reasonable degree that the proposed method can sample from non-trivial non-Gaussian distributions (MNIST). This is an improvement on most comp neuro papers on sampling, which largely focus on Gaussian distributions. My enthusiasm for this point is, however, dampened a bit by the fact that the networks here are rate-based and the authors aren't trying to optimize functionals of the dynamics (e.g. as in work from the Lengyel group), so it's not very surprising that this can be made to work. ", "weaknesses": "1. How does Proposition 1 differ from the results of Ma, Chen, and Fox, \"A Complete Recipe for Stochastic Gradient MCMC\" (NeurIPS 2015)? Perhaps I've missed something, but I do not see a novel contribution in Section 3.1. \n\n2. The present manuscript does not address the speed of relaxation to the stationary distribution. The ability to obtain accurate samples quickly is a necessary feature of any model of sampling-based inference in the brain, and indeed the question of how to achieve fast sampling in biologically-plausible implementations has been a major focus of past work (e.g., Hennequin et al. 2014, Aitchison and Lengyel 2016, Echeveste et al. 2020, Masset et al. 2022). In my view, this is a key shortcoming of the learning algorithm proposed in Section 3.4: no attention is paid to whether the resulting network enjoys favorable convergence properties. Of course, one could modify the resulting drift and diffusion terms (as in Hennequin et al. 2014 or Masset et al. 2022), but this would require an additional optimization step. The authors should at least discuss this point; adding more experiments to probe this issue would make for a stronger paper. \n\n3. In Lines 293-295 of the Discussion, the authors state that \"Our preliminary experiments show that one-hidden-layer RSNs cannot readily approximate high-dimensional heavy-tailed distributions (e.g., those of overcomplete sparse coding representations [36]).\" If the authors have such results, they should show them. In particular, related to point (2) above, they should show how the convergence speed and number of reservoir neurons required to obtain accurate samples scales with the dimension of the target distribution. As the authors themselves discuss, their approach is currently limited to approximating the score with a single-hidden-layer MLP, and such Cybenko-Hornik-type approximation results may require immensely wide networks. This is also a limitation of the MNIST experiment: MNIST is intrinsically not that high-dimensional, and already here the authors used a network with 20000 reservoir neurons (Line 638)! \n\n4. I think it's important to mention that the striking result of Figure 2 depends critically on the choice of ReLU activation function, particularly as much of the paper is concerned with existence results. If instead one used $\\tanh$, then the performance gap between the trained sampler-only firing rate and reservoir firing rate models should be small. In particular, a sampler-only firing rate network can implement exact Langevin sampling in this case. The score of $p(x) = \\frac{1}{2} \\mathcal{N}(\\mu , \\sigma^{2}) + \\frac{1}{2} \\mathcal{N}(-\\mu,\\sigma^2)$ is $\\frac{d}{dx} \\log p(x) = \\frac{1}{\\sigma^{2}} \\left[ - x + \\mu \\tanh\\left(\\frac{\\mu}{\\sigma^{2}} x\\right) \\right]$. This gives the Langevin dynamics $dx(t) = \\frac{d \\log p}{dx} dt + \\sqrt{2} dB(t) = \\frac{1}{\\sigma^{2}} \\left[ - x + \\mu \\tanh\\left(\\frac{\\mu}{\\sigma^{2}} x\\right) \\right] dt + \\sqrt{2} dB(t)$. With $\\mu = 1$ as is used in Figure 2, this is exactly of the form of a single-neuron sampler-only firing rate network in eq (8) of the submitted manuscript, with $D = 1/\\sigma^{2}$, $W_{rec} = 1/\\sigma^{2}$, $I=0$, and $\\phi(x) = \\tanh(x)$. One can easily verify numerically that this network produces accurate samples from the desired bimodal distribution. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer’s critiques and suggestions. Below we address the reviewer’s concerns point by point:\n\nStrengths:\n\nWith respect to the two enthusiasm-limiting aspects, we believe we understand, and in our revision will add text specifically pointing these limitations out, with citations to the relevant literature (e.g. Echeveste et al. 2020 and allied papers) as well as the explicit suggestions that future work should extend our findings to these settings.  In more detail, we will explain: First, it is true that we study the ability to produce samples from given “target” distributions, as opposed to studying (input-dependent) networks that additionally solve an inference problem (as in Echeveste et al. 2020).  Second, we study this problem in the context of continuous-valued “rate” networks as opposed to jump-discontinuous spiking systems.   \n\nThis said, we do believe that our results, which are novel in their scope of general distributions and fairly complex technically, are best presented in terms of sampling from fixed, general (beyond-gaussian)  target distributions and in terms of the broadly used, and mathematically tractable rate networks.  (With respect to a tractable partial extension towards spiking networks, please also see our response to Q3 of Reviewer Y31x.)    \n\nWeakness:\n\n1. We appreciate the reference to Ma, Chen, and Fox 2015, and in our revised Section 3.1 we will both cite this paper and explain the distinction of our results from their Proposition.  Specifically, these authors assume a particular form of the key divergence-free field (G in our paper), and to do so they introduce a second “gamma’’ term to keep G divergence-free (see equation 3 therein). However, with this approach to the gamma term, it is not easy to relate the expressivity of a neural network’s “drift” (i.e., its dynamics, given by decay plus weights times activations) term with the score function (specifically, without knowledge of their Q matrix and its derivative).  In section 3.1 and appendix A, we take a novel and distinct approach, based on the Helmholtz-Hodge decomposition, that allows us to directly relate the score and required aspects of the network dynamics. Specifically, Section 3.1 shows that even if the neural dynamics are allowed to be irreversible (G is not 0), the function class of all dynamics still needs to have enough basis functions to approximate the score function part of the drift term that determines the stationary distribution. \n\n2. We thank the reviewer for this critique and their important idea here. This has led us to an explicit new treatment of sampling speed (relaxation) to be included in our revision -- please see the new figure in the .pdf attached to our global reply.  We will also add text covering this important topic and the literature the reviewer notes in our revision. In more detail, following the reviewer's critique we have identified a simple way to increase sampling speed without additional optimization steps.  This is by letting the divergence-free (DF) field G in our paper be of the form $J\\nabla p$ (treat Q as constant in Ma, Chen, and Fox 2015), where $J = -J^T$ can be any skew-symmetric matrix that is either learned through training or prescribed. Then the equation 14 in our paper would be $2\\alpha (I + J)^{-1}(W_{out} \\phi(\\tilde{W}_{rec} x + I) - x) \\approx \\nabla\\log p(x)$. The rest of the training procedure stays the same. Note that we do not need to alter the form of the dynamics in Equation 11. This improvement to the score matching loss gives a significant increase in sampling speed as we show in Figure S2. \n\n3. As requested, we show the preliminary results of our reservoir sampler trained on the high-dimensional latent representations produced by a sparse coding model on MNIST (Fig. S4). Here, the latent space dimension is 3136 (four times overcomplete) and we used 30000 hidden neurons in the sampler.  We do note that this result is highly preliminary, and based on our work with the PCA projection we believe that it is possible that future work (on more powerful GPU systems) with larger reservoirs (though, of course, not yet approaching millions of units available in local circuits of the mammalian cortex) and/or other learning rules and hyperparameters are quite likely to succeed.\nIn addition, we argue that MNIST is low-dimensional only if it undergoes a nonlinear transformation.  The first 300 PC components from MNIST just explain ~90% of the variance. Using a nonlinear autoencoder, we can indeed learn the distribution with a much smaller number (500) of reservoir neurons; see also Figure S3. \n\n4. We thank the reviewer for this insightful comment: they are absolutely right that a tanh nonlinearity can lead to sampling of bimodal distributions.  However, another limitation does arise in this case: since tanh(x) and (x) are both odd functions, any distribution whose score function does not have a related symmetry property will not be able to be approximated (See Figure S1). \n\nQuestions:\n\n- We thank the reviewer for pointing out the missing literature, it will be added in the camera-ready version. \n\n- The proof of theorem 3 is constructive in the sense that it constructs the weight of the RSN from a single-hidden-layer MLP.  We will add text to the revision to clarify this meaning. \n\n- Line 303-304 is referring to the fact that denoising score matching learns from gradually less noisy samples from a target distribution. Here we notice the striking similarity of this procedure to how the infant brain learns less noisy distribution of vision representations over time. \n\n- The main application of our work is to mechanistically model neural circuits that produce varied, non-Gaussian prior distributions that support Bayes-optimal behavior (reference 21, 26 in our paper)  \n  \n- With thanks, we will make the implied suggestions in the other careful reviewer questions as well.\n"}]

Author rebuttal:
We thank all reviewers for their expert review and insightful critiques, suggestions, and comments.  We have written a detailed, point by point reply to each reviewer's issue or concern in the replies to individual reviewers.  

In the pdf attached to this general response, we illustrate the results requested by the reviewers.  Specifically, we show more results on the double peak experiment using tanh nonlinearity (Figure S1), accelerated sampling with irreversible dynamics (Figure S2), improved MNIST generation using autoencoder (Figure S3), and the preliminary result on the sparse coding setting (Figure S4). 


Meta-review: Recurrent neural circuits and algorithms are studied in the broad context of Bayesian sampling hypothesis of the brain. Authors demonstrated non-trivial distribution can be sampled. Reviewer-Author discussions have generated potential improvements and ideas and I suggest the authors put an extended amount of efforts in making the final version addressing the discussion on new comparisons and clarifications.
Decision: Accept (poster)