Note: Moment approximations for Bernoulli neurons with sigmoidal nonlinearity

Consider a stochastic, binray, linear-nonlinear unit, with spiking output $s$, synaptic inputs $\mathbf{x}$, weights $\mathbf{w}$, and bias (threshold) $b$:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}s& \sim \mathrm{Bernoulli}[p=\mathrm{\Phi }(a)]\\ a& ={\mathbf{w}}^{\mathrm{\top }}\mathbf{x}+b,\end{array}\end{array}$$

where $\mathrm{\Phi }(\cdot )$ is the cumulative distribution function of a standard normal distribution. Note that $\mathrm{\Phi }(\cdot )$ can be rescaled to closely approximate the logistic sigmoid if desired. Assuming the mean $\mu$ and covariance $\mathrm{\Sigma }$ of $\mathbf{x}$ are known, can we obtain the mean and covariance of $s$?

The wrong way: A small-variance approximation

It is common to model the stochastic response in terms of the mean-field (deterministic) transfer function $\mathrm{\Phi }(a)$, plus a small correction assuming that the variance in $a$ is small. This works acceptable for weakly-nonlinear units, like Poisson, but fails for Bernoulli units.

One can approximate the variance of the spiking $s$ output as the sum of the variance from the Bernoulli sampling, and the variance in the Bernoulli rate itself. The variance of a Bernoulli variable is $p(1-p)$, and the variance in $p$ itself can be obtained from a locally-linear Gaussian approximation as the variance of the activation, multiplied by the slope of the effective transfer function $f^{\prime }({\mu }_a)$:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\sigma }_s^2& =p(1-p)+f^{\prime }({\mu }_a{)}^2\cdot {\sigma }_a^2\end{array}\end{array}$$

This approximation is convenient, and works for generic nonlinearities (not just $\mathrm{\Phi }(a)$). More generally, we can use the linear-Gaussian approximation to estimate the covariance (${\mathrm{\Sigma }}_s$) of a population of output neurons driven by shared (noisy) inputs:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\mu }_s& =f(W{\mu }_x+B)\\ {\mathrm{\Sigma }}_s& =J{\mathrm{\Sigma }}_x{J}^{\mathrm{\top }}+Q\\ J& ={\partial }_x\mu (⟨x⟩)\\ Q& =\mathrm{Diag}[\mu (x)(1-\mu (x))],\end{array}\end{array}$$

Where $W$ denotes a matrix of input weights, $B$ denotes a vector of per-unit biases, and $\overrightarrow{\mu }(x)$ reflects a (mean-field estimate) of the vector of output mean-rates.

This approximation arises from a locally-linear approximation of the nonlinear transfer function $f$, which transforms correlated inputs ${\mathrm{\Sigma }}_x$ according to the jacobian $J$ of the firing-rate nonlinearity. This also includes a Gaussian (diffusive) approximation of the noise arising from Bernoulli sampling, $Q$, which is equal to $p(1-p)$. A similar result holds for the Poisson (low firing-rate limit), for which $Q=\mathrm{Diag}[\overrightarrow{\mu }(x)]\mathrm{\Delta }t$.

Additional corrections can be added, for example accounting for the effect of variance on $\overrightarrow{\mu }(x)$ for estimating the noise source term $Q$, or providing additional corrections based on higher-order moments or moment-closure approximations thereof.

Similar locally-quadratic approximations for noise have been advanced in the context of chemical reaction networks (Ale et al. 2013). Rule and Sanguinetti (2018) and Rule et al (2019) use this approach in spiking neuron models, and Keeley et al. (2019) explored spiritually similar quadratic approximations for point-processes.

The small-variance correction is essentially the first term in a family of series expansions, which use the Taylor expansion of the firing-rate nonlinearity to capture how noise is transformed from inputs to outputs. In the case of linear-nonlinear-Bernoulli neurons, approximations based on series expansions like this have poor convergence. Polynomial approximations to $\mathrm{\Phi }$ diverge as the activation becomes very large or very small, whereas the sigmoidal nonlinearity is bounded. Global approximations are therefore desirable when the variance of the input is large.

A better way: Dichotomized Gaussian (probit) moment approximation

Global approximations have been presented elsewhere for other types of firing nonlinearity. Echeveste et al. (2019) used exact solutions for propagation of moments for rectified-polynomial nonlinearities (Hennequin and Lengyel 2016). Rule and Sanguinetti (2018) also illustrate an example with exponential nonlinearities.

These approaches fall under the umbrella of "moment-based methods", and entail solving for the propagation of means and correlations under some distributional anstaz (often Gaussian, although see Byrne et al. 2019 for an important application using circular distributions). In general, there are few guarantees of accuracy for these methods (Schnoerr et al. 2014, 2015, 2017), although they are often empirically useful.

Moment approximations fair poorly for the linear-nonlinear-Bernoulli neuron. However, when one takes the firing-rate nonlinearity to be the CDF of the standard normal distribution, global approximations are possible. This yields suitable approximations for other sigmoidal nonlinearities, provided that these nonlinearities can be approximated by the normal CDF under a suitable change of variables.

The variance and covariances in a population of dichotomized Gaussian neurons can be expressed in terms of the multivariate normal CDF. To derive the population covariance, consider a single entry which reflects the covariance between a pair units.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathrm{\Sigma }}_{12}& =⟨(s_1-p_1)(s_1-p_2)⟩\\ & =⟨s_1{s}_2⟩-p_1{p}_2\\ ⟨s_1{s}_2⟩& =Pr(s_1=s_2=1)\end{array}\end{array}$$

If $a=w^{\mathrm{\top }}x+b$ is the activation, and $u=a+\xi$ is the activity combined with zero-mean unit-variance threshold noise $\xi$, we can evaluate $⟨s_1{s}_2⟩$ by considering the joint distribution of $u_1$ and $u_2$ as Gaussian (for a numerical recipe see Drezner and Wesolowsky 1989; numerical implementations are provided in standard computing packages, e.g. Matlab or Scipy in Python):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}⟨s_1{s}_2⟩& =Pr(u_1>0\text{ and }{u}_2>0).\\ u& \sim \mathcal{N}({\mu }_u,{\mathrm{\Sigma }}_u)\\ {\mu }_u& ={\mu }_a\\ {\mathrm{\Sigma }}_u& ={\mathrm{\Sigma }}_a+\mathrm{I}.\end{array}\end{array}$$

A numerical solution in terms of the bivariate Gaussian CDF is useful for propagating activity, but challenging for building a differentiable model suitable for optimization. However, practical approximations exist.

Faster approximations to dichotomized Gaussian moment approximation

For a single neuron, the mean and variance of the spiking output are those of a $\mathrm{Bernoulli}(p)$ distribution, with probability $p=Pr(s=1)$. The mean rate ${\mu }_s$ is equal to the probability of firing ($p$), and the variance ${\sigma }_s^2$ is equal to $p(1-p)$ (Fig a).

Binary spiking units with a Gaussian CDF nonlinearity $\mathrm{\Phi }(\cdot )$ can be modeled as a thresholded Gaussian noise source. When this noise is above a certain threshold ($-a$), the unit emits a "1", otherwise, a "0". This makes it easy to model the effect of additional noise (variance) in the synatpic activation "a". This extra noise simply sums with the existing Gaussian noise in the model of the stochastic spiking.

The spiking probability $p$ of a dichotomized-Gaussian unit being driven by noisy, Gaussian inputs can be obtained by treating the effect of noise in activation (${\sigma }_a^2$) as a decrease in gain:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}a& \sim \mathcal{N}({\mu }_a,{\sigma }_a^2)\\ p& =⟨\mathrm{\Phi }(a)⟩=\mathrm{\Phi }\left(\gamma {\mu }_a\right)\\ \gamma & =\frac{1}{\sqrt{1+{\sigma }_a^2}}\\ {\mu }_s& =p\\ {\sigma }_s^2& =p(1-p)\end{array}\end{array}$$

To see this in more depth, observe that the variance in the spiking output ${\sigma }_s^2$ is a combination of the average spiking variance ${\sigma }_{\text{noise}}^2=⟨p(1-p)⟩$, plus whatever input noise in the firing rate (${\sigma }_p^2$) is passed through the nonlinearity. In the dichotomized Gaussian model of a linear-nonlinear-Bernoulli neuron, we find that ${\sigma }_s^2\approx {\mu }_p(1-{\mu }_p)$:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\sigma }_{\text{noise}}^2& =⟨p(1-p)⟩\\ & =⟨p⟩-⟨p^2⟩\\ & ={\mu }_p-({\mu }_p^2+{\sigma }_p^2)\\ & ={\mu }_p(1-{\mu }_p)-{\sigma }_p^2\\ \\ {\sigma }_p& \approx {\sigma }_a\cdot ⟨{\partial }_a\mathrm{\Phi }(a)⟩\\ & ={\sigma }_a\cdot {\partial }_a⟨\mathrm{\Phi }(a)⟩\\ & ={\sigma }_a\cdot {\partial }_a\mathrm{\Phi }({\mu }_a\gamma )\\ & ={\sigma }_a\cdot \varphi ({\mu }_a\gamma )\cdot \gamma \\ \\ {\sigma }_s^2& ={\sigma }_p^2+{\sigma }_{\text{noise}}^2\\ & ={\sigma }_p^2+[{\mu }_p(1-{\mu }_p)-{\sigma }_p^2]\\ & ={\mu }_p(1-{\mu }_p).\end{array}\end{array}$$

This generalizes to the multivariate case, and provides an approximation for how correlations in inputs propagate to correlations in the output:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\mathrm{\Sigma }}_s& ={\mathrm{\Sigma }}_p+{\mathrm{\Sigma }}_{\mathrm{noise}}\\ {\mathrm{\Sigma }}_p& \approx J{\mathrm{\Sigma }}_a{J}^{\mathrm{\top }}\\ {\mathrm{\Sigma }}_{\mathrm{noise}}& \approx \mathrm{Diag}[p(1-p)-{\sigma }_p^2]\\ J& =\mathrm{Diag}[\varphi (\gamma {\mu }_a)\cdot \gamma ]\\ p& =\mathrm{\Phi }(\gamma {\mu }_a),\\ \gamma & ={(1+\mathrm{Diag}[{\mathrm{\Sigma }}_a])}^{-\frac{1}{2}}\end{array}\end{array}$$

The accompanying figure shows a toy example of variance approximation, using a network of three neurons (Fig. b). Compared to the small-variance approximation, the approximation derived for the dichotomized Gaussian case provides a better approximation of the moments of the output, and accounts for how noise in the input propagates to the output (Fig. c).

Figure: variance propagation in the dichotomized Gaussian neuron (a) For a single neuron, the effect of input variability (${\sigma }_a^2$) can be viewed as a modulation of the gain of the nonlinear transfer function. The output variance is then similar to that of a Bernoulli distribution. (b) In a feed-forward network of nonlinear stochastic neurons, noise propagates to downstream neurons, affecting the computational properties of the circuit. (c) The output (blue) of this circuit is stochastic, and noise in the first layer (black, red) propagates to the output (left panel: Monte-Carlo samples, shaded = 5-95${}^{th}$ percentile), but can be modeling in a differentiable way using moment approximation. The small variance approximation (linear noise approximation or LNA, in this case: middle ) loses some accuracy for small circuits, since the is very little averaging to attenuate spiking noise. The moment approximation using a dichotomized Gaussian (DG) model is more accurate (right).