Neural Empirical Bayes

Introduction

This setup is denoted symbolically as $X\rightharpoonup Y$.

with $\epsilon\ll 1$ and $d_{\sharp}\ll d$. We revisit textbook calculations on concentration of measure and show that in high dimensions, the smoothed manifold $\mathscr{N}$ is approximately a sphere of dimension $d-1$:

This analysis paints a picture of “manifold disintegration-expansion” as a general phenomenon.

The first conceptual contribution of this paper is to expose a notion of manifold disintegration-expansion, with a general takeaway that in (very) high dimensions, gaussian smoothing pushes away (all) the probability masses near $\mathscr{M}$. This is due to the fact that in high dimensions, the Gaussian random variable

is not necessarily concentrated on $\mathscr{M}$, and there are many directions (asymptotically infinite) where samples from $Y$ “escape” the manifold. The thesis is that, in convolving $f_{X}$ and $f_{N}$, $\mathscr{M}$ would be mapped to a (much) higher dimensional manifold.

A seemingly unrelated theory is empirical Bayes, as formulated by Robbins (1956), which we use for pulling the probability mass back towards $\mathscr{M}$. In 1956, Robbins considered a scenario of a random variable $Y$ that depends “in a known way” on an “unknown” random variable $X$. The starting point in (Robbins, 1956) was the existence of a noisy random variable that was denoted by $X$. In that setup, one can only observe values of $X$. Here, we started with the “clean” i.i.d. sequence $X_{1},\dots,X_{n}$ and artifically created $Y=X+N(0,\sigma^{2}I_{d})$. This departure in starting points is related to our new take on empirical Bayes which will become clear in Section 3. He found that, given an observation $Y=y$, the least squares estimator of $X$ is the Bayes estimator, and quite remarkably, can be expressed purely in terms of the distribution of $Y$ (he showed this for Poisson, geometric, and Laplacian kernels). Robbins’ results were later extended to gaussian kernels by (Miyasawa, 1961) who derived the estimator

The least-squares estimator of $X$ above is derived for any $\sigma$, not necessarily infinitesimal. To signify this important fact, we also refer to this estimation as “Robbins jump” or jump for short. The empirical Bayes machinery is denoted symbolically as

The smoothing of kernel density estimation, $X\rightharpoonup Y$, and the denoising mechanism of empirical Bayes, $X\leftharpoondown Y$, come together to define the learning objective

(Throughout the paper, $\nabla$ is the gradient taken with respect to the inputs $y$, not parameters.)

When analyzing a finite i.i.d. sequence, $X_{1},\dots,X_{n}$, i.i.d. samples from $Y=X+N(0,\sigma^{2}I_{d})$ are generated as

and the learning objective $\mathcal{L}(\theta)$ is then approximated as

where $\sum$ is a shorthand for $\frac{1}{nm}\sum_{i=1}^{n}\sum_{j=1}^{m}$. In high dimensions, the samples $Y_{ij}$ are (approximately) uniformly distributed on a “thin spherical shell” around the sphere $\sigma\sqrt{d}S_{i}^{d-1}$ with its center at $X_{i}$ which leads to the following definition.

The sphere $\sigma\sqrt{d}S_{i}^{d-1}$ centered at $X_{i}$ is a useful abstraction and it is referred to as “$i$-sphere”. The $i$ in $i$-sphere refers to the index $i$ in $X_{i}$, and the index $j$ is reserved for the gaussian noise $\varepsilon_{j}$. Therefore, $Y_{ij}$ is the $j$th sample on the $i$-sphere. Fixing the index $i$, the samples $Y_{ij}$ are approximately uniformly distributed on a thin spherical shell around the $i$-sphere (see Figure 1).

DEEN: minimize $\mathcal{L}(\theta)$ with stochastic gradient descent and return $\theta^{*}$.

In most of the paper, we are interested in studying $\phi$ which is already at optimality with parameteres $\theta^{*}=\operatorname*{\rm{argmin}}\mathcal{L}(\theta)$. In these contexts, $\theta^{*}$ is dropped, and its presence is understood, e.g. at optimality, $\phi\approx-\log f$ (modulo a constant).

Next, we define a notion of interactions between $i$-spheres that we find useful in thinking about the goal of approximating the score function through optimizing $\mathcal{L}(\theta)$. However, this can be skipped until our discussions of the associative memory in Sections 6 and 7.

The interaction between $i$-spheres is defined as the competition that emerge due to the presence of $\nabla\phi$ in the problem of minimizing $\mathcal{L}$. In other words, it refers to “the set of constraints” that $\nabla\phi$, the gradient of the globally defined energy function, has to satisfy locally on $i$-spheres to arrive at the optimality such that $\nabla\phi\approx-\nabla\log f.$ It is indeed an abstract notion in this work, but it is a useful abstraction to have which we come back to in thinking about the gradient flow $y^{\prime}=-\nabla\phi(y).$ It should be stated that there is also a notion of interaction that already exists for the kernel density in the sense of “mass aggregation” (see Section 7 for references), but in this framework, $i$-spheres interact even when they do not overlap in that they “communicate” via $\nabla\phi$.

The third conceptual contribution of this paper is introducing the physical picture of interactions between i.i.d. samples. This is a useful abstraction in this work in thinking about the problem of approximating the score function with the gradient of an energy function. The interaction is a code for “the set of constraints” that $\nabla\phi$ (evaluated at i.i.d. samples $Y_{ij}$) has to satisfy to arrive at optimality, $\nabla\phi\approx-\nabla\log f$. These interactions could also give rise to some collective phenomena (Anderson, 1972).

Next, we outline the technical contributions of this paper. Approximating $\nabla\log f$ in the empirical Bayes setup leads to a novel sampling algorithm, a new notion of associative memory and the emergence of “creative memories”. First, we define the matrix $\chi$ that was introduced to quantify $i$-sphere overlaps and was used for designing experiments.

Concentration of measure leads to a geometric picture for the kernel density as a “mixture of $i$-spheres”, and we define the matrix $\chi$ to quantify the extent to which the spheres in the mixture overlap; the entries in $\chi$ are essentially pairwise distances, scaled by $2\sqrt{d}$ (where the scaling is related to the concentration of isotropic Gaussians in high dimensions):

The $i$-sphere and the $i^{\prime}$-sphere do not overlap if $\sigma<\chi_{ii^{\prime}}$, and they overlap if $\sigma>\chi_{ii^{\prime}}$ (see Figure 2). Of special interest is $\sigma_{c}$, defined as

We define “extreme noise” as the regime $\sigma>\sigma_{c}$; it is the regime where all $i$-spheres have some degree of overlap.

“Walk-jump sampling” is an approximative sampling algorithm that first draws exact samples from the density $\exp(-\phi)/Z$ ($Z$ is the unknown normalizing constant) by Langevin MCMC:

where $\delta\ll 1$ is the step size and here $t$ is discrete time. At an arbitrary time $\tau$, approximative samples “close” to $\mathscr{M}$ are generated with the least-squares estimator of $X$—the jumps:

A signature of walk-jump sampling is that the walks and the jumps are decoupled in that the jumps can be made at arbitrary times.

“Neural empirical Bayes associative memory”—named NEBULA, where the “LA” refers to à la Hopfield (1982)—is defined as the flow to strict local minima of the energy function $\phi$. In continuous time, the memory dynamics is governed by the gradient flow:

In retrieving a “memory”, one flows deterministically to an attractor. NEBULA is an intriguing construct as the deterministic dynamics is governed by a probability density function, since $\nabla\phi\approx-\nabla\log f$—this is very far from true for Hopfield networks.

Memories in NEBULA are believed to be formed due to $i$-sphere interactions (Definition 3). We provide evidence for the presence of this abstract notion of interaction in this problem by reporting the emergence of very structured memories in a regime of highly overlapping $i$-spheres. They are named “creative memories” because they have intuitively appealing structures, while clearly being new instances, in fact quite different from the training data.

This paper is organized as follows. Manifold disintegration-expansion is discussed in Section 2. In Section 3, we review empirical Bayes and give a derivation of the least squares estimator of $X$ for gaussian kernels. We then talk about a Gedankenexperiment—in the school of Robbins (1956) and Robbins and Monro (1951)—to bring home the unification scheme $X\rightleftharpoons Y$ and shed some light on the inner-working of the learning objective. In Section 4, we define the notion of extreme noise and demonstrate two extreme denoising experiments. In Section 5, we present the walk-jump sampling algorithm. NEBULA is defined in Section 6, where the abstract notion of $i$-sphere interactions is grounded to some extent; we present two sets of experiments with qualitatively different behaviors. In Section 7, we report the emergence of creative memories for two values of $\sigma$ close to $\sigma_{c}$. We finish with a summary.

Manifold disintegration-expansion in high dimensions

Our low-dimensional intuitions break down spectacularly in high dimensions. This is in large part due to the exponential growth of space in high dimensions (this abundance of space also underlies the curse of dimensionality, well known in statistics and machine learning). Our focus in this section is to develop some intuitions on the effect of gaussian smoothing, $X\rightharpoonup Y$, on the data manifold $\mathscr{M}$. We start by a summary of results on concentration of measure. The textbooks (Ledoux, 2001; Tao, 2012; Vershynin, 2018) should be consulted to fill in details. We then extend the textbook calculations for a “gaussian manifold” and make a case for manifold disintegration-expansion.

Start with the isotropic gaussian. In 2D, samples from the gaussian form a “cloud of points”, centered around its mode as illustrated in Figure 1a. Next, take the isotropic Gaussian in $d$ dimensions, $(X_{1},\dots,X_{d})=N(0,I_{d})$, and ask where the random variable is likely to be located. Before stating a well-known result on the concentration of measure, we can build intuition by observing the following identities for the expectation and the variance of the squared norm:

Since the components of $X=(X_{1},\dots,X_{d})$ are jointly independent, and since $d\gg 1$, the following holds with high probability,

The concentration of norm is visualized in Figure 1b; also see Figure 3.2 in (Vershynin, 2018). In contradiction to our low-dimensional intuition, in high dimensions, the Gaussian random variable $N(0,I_{d})$ is not concentrated close to the mode of its density $f_{N}$ (at the origin here), but in a “thin spherical shell” of width $O(1)$ around the sphere of radius $\sqrt{d}$. The concentration of norm suggests an even deeper phenomenon, the concentration of measure, captured by a non-asymptotic result for the deviation of the norm $\|X\|$ from $\sqrt{d}$, where the deviation has a sub-gaussian tail (Tao, 2012; Vershynin, 2018):

In above expressions, $C$ and $c$ are absolute constants. In high dimensions, the Gaussian random variable is thus approximated with the uniform distribution on the sphere $\sigma\sqrt{d}S^{d-1}$,

The approximation becomes equality in distribution, as $d\rightarrow\infty$.

where $\epsilon\ll 1$. The manifold $\mathscr{M}$ is viewed as a “gaussian manifold” where, in an abuse of notation, $\text{dim}(\mathscr{M})\approx d_{\sharp}$. Now consider $Y=X+N(0,\sigma^{2}I_{d})$ which means

We repeat the concentration of norm calculations:

the “manifold terms”, $d_{\sharp}(1+\sigma^{2})$ and $d_{\sharp}(1+\sigma^{2})^{2}$, become negligible, and the following holds with high probability:

The calculations in the previous page are reproduced by setting $\sigma=1$, $\epsilon=0$, $d_{\sharp}=0$.

In summary, under the convolution $f_{Y}=f_{X}*f_{N}$, the “gaussian manifold” $\mathscr{M}$ is transformed/mapped to

therefore $\text{dim}(\mathscr{N})\approx d-1$. In our terminology, $\mathscr{M}$ has been disintegrated-expanded to $\mathscr{N}$. We expect this phenomenon to hold for any $\mathscr{M}$ asymptotically, i.e. as $d\rightarrow\infty$, $\mathscr{N}=\sigma\sqrt{d}S^{d-1}$, but (unfortunately) the limit $d\rightarrow\infty$ is not practical. However, the picture should stay: in high dimensions, the convolution of $f_{X}$ and $f_{N}$ has severe side effects on the manifold $\mathscr{M}$ itself. Smoothing is indeed desired, but also a mechanism to restore $\mathscr{M}$. In the next section, we discuss such a mechanism using empirical Bayes, where the “restoration” is in least squares.

Neural empirical Bayes

In a seminal work from 1956, titled an empirical Bayes approach to statistics, Robbins considered an intriguing scenario of a random variable $Y$ having a probability distribution that depends “in a known way” on an “unknown” random variable $X$ (Robbins, 1956). Observing $Y=y$, the least-squares estimator of $X$ is the Bayes estimator (see Remark 1):

If $f_{X}$ is known to “the experimenter”, $\widehat{x}$ is a computable function, but what if the prior $f_{X}$ is unknown? It is quite remarkable that for a large class of kernels, the least-squares estimator can be derived in closed form purely in terms of the distribution of $Y$. Informally speaking, there is an “abstraction barrier” (Abelson and Sussman, 1985) between $X$ and $Y$ where the knowledge of $f_{X}$ is not needed to estimate $X$. This functional dependence of $\widehat{x}$ on $f_{Y}$ was achieved for the Poisson, geometric, and Laplacian kernels in (Robbins, 1956), and it was later extended to the gaussian kernel by (Miyasawa, 1961).

Multiply the expression above by $f(x)$ and integrate,

which then leads to the expression $\widehat{x}(y)=y+\Sigma~{}\nabla\log f(y).$ For the isotropic case $\Sigma=\sigma^{2}I_{d}$, the estimator takes the form

To sum up, the estimator above is obtained in a setup where only the corrupted data (the random variable $Y$) are observed, with the knowledge of the measurement (gaussian) kernel alone and without any knowledge of the prior $f_{X}$. The remarkable result is the fact the least-squares estimator of $X$ is written in closed form purely as a functional of $\nabla\log f$, also known as the score function (Hyvärinen, 2005). The expression above for the least squares estimator of $X$ is the basis for an algorithm to approximate $\nabla\log f$ which is discussed next.

In empirical Bayes, the random variable $X$ is estimated in least squares from the noisy observations $Y$, without any knowledge of the prior $f_{X}$. But how did we get here? We started the paper with the i.i.d. sequence $X_{1},\dots,X_{n},$ where $Y$ did not even exist!

The idea is a Gedankenexperiment of sort, where we first construct $Y=X+N(0,\sigma^{2}I_{d})$ by taking samples $Y_{ij}\sim f_{X}*f_{N}$; in turn, the experimenter (in the school of empirical Bayes) “observes” $Y_{ij}$ and estimates $X$. But in this artificially supervised setup (the “target” is $X_{i}$), the error signal $\|X_{i}-\widehat{x}(Y_{ij})\|^{2}$ can also be measured—this is not the case in (Robbins, 1956)—and it drives the learning. The learning is achieved with stochastic gradient descent as the experimenter has parameterized the energy function with a neural network and can compute the stochastic gradients (Robbins and Monro, 1951),

Approximating $\nabla\log f$ in this manner is the essence of neural empirical Bayes.

For finite samples, the best we can do is to approximate $\nabla\log f$. Given $X_{1},\dots,X_{n}$, the i.i.d. samples $Y_{ij}\sim f_{X}*f_{N}$ are given by $Y_{ij}=X_{i}+\varepsilon_{j},\ \varepsilon_{j}\sim N(0,\sigma^{2}I_{d}).$ These are indeed i.i.d. samples from density estimated by the kernel density estimator

At optimality (achieved with stochastic gradient descent), $\theta^{*}=\operatorname*{\rm{argmin}}\mathcal{L}(\theta)$,

Extreme noise/denoising

In empirical Bayes, the least-squares estimator of $X$,

is derived for any $\sigma$, which inspires us to call this a “jump” to capture the fact $\sigma$ may in fact be “large”. Denoising is not our actual goal, but the expression above can be viewed as “denoising $y$ to $\widehat{x}(y)$”. The denoising performance is important since in the machinery of neural empirical Bayes, the goal of approximating the score function is formulated by the “denoising objective”

As the least-squares estimator of $X$ is derived for any $\sigma$, the denoising performance of DEEN may also be tested to its limits. Here, we report such experiments for MNIST, where $d=784$ and $X$ is in the hypercube $^{d}$ (see Remark 8). Next, we elaborate on the geometric meaning of the noise levels, where in addition we define a notion of “extreme noise”.

$(\sigma>\sigma_{c})$. The value $\sigma_{c}=\max_{ii^{\prime}}\chi_{ii^{\prime}}$ is defined as the onset of extreme noise. Called “extreme”, because for $\sigma>\sigma_{c}$, all $i$-spheres in the dataset overlap (see Figure 2b). For MNIST, $\sigma_{c}\approx 0.2884$ was estimated from $10^{7}$ handwritten digit pairs (see Table 1). The results for $\sigma=0.3$ on the test set is presented in Figure 3. It should be stated that, this value of “extreme noise”, just above $\sigma_{c}$, is not visually perceived as extreme.

$(\sigma>2\sigma_{c})$. In this regime, due to geometry, the whole dataset is in inside of each $i$-sphere (see Figure 2c). This is very extreme! For MNIST, we experimented with $\sigma=0.7$ which is in this regime. The results on the test set are presented in Figure 4.

The following comments are in order regarding the network architecture we used in our experiments.

(activation function) The denoising results are a significant improvement over (Saremi et al., 2018). This was due to the use of ConvNets, instead of a fully connected network (more on that below) and the use of a “smooth ReLU” activation function

where the default was $\beta=1$. The activation function above converges uniformly to ReLU in the limit $\beta\rightarrow\infty$. It has been studied independently from different angles (Elfwing et al., 2018; Ramachandran et al., 2017). ReLU consistently gave a higher loss, which we believe is due to the term $\nabla_{y}\phi$ in the learning objective, which must be computed first before computing $\nabla_{\theta}\mathcal{L}$ to update the parameters $\theta\leftarrow\theta-\epsilon\nabla_{\theta}\mathcal{L}.$

(wide convnets) All experiments reported in the paper were performed in a fixed wide ConvNet architecture with the expanding $\rm channels=(256,512,1024)$, without pooling, and with a bottleneck layer of size $10$. All hidden layers were activated with the activation function $\sigma(\cdot,\beta=1)$, and the readout layer was linear.

(readout overparametrization) We observed a slightly faster convergence to lower loss (especially, very early in the training) by overparameterizing the linear output layer. These experiments were inspired by the “acceleration effects” studied for linear neural networks (Arora et al., 2018), but we did not study it thoroughly.

We used the Adam optimizer (Kingma and Ba, 2014). The optimization was stable over a wide range of mini-batch sizes and learning rates, and as expected, we did not observe the validation/test loss going up in the experiments. This stability is due to the fact that the inputs to $\phi$ are noisy samples. The automatic differentiation (Baydin et al., 2018) was implemented in PyTorch (Paszke et al., 2017). DEEN is “memory hungry” but the computational costs are not addressed here in detail.

Walk-jump sampling: Langevin walk, Robbins jump

than $f_{X}$ (putting aside the harder problem of learning $\nabla\log f_{X}$). At any time, the estimator $\widehat{x}$ can be used to jump near $\mathscr{M}$. “Near” is intuitive, and it should be stated that neither $f_{X|Y}$ nor $f_{X}$ is exactly sampled during jumps; the jump samples must be seen as heuristic.

In what follows, it is understood that $\phi$ is at optimality. We first give a description of the algorithm and then express the equations.

(walks) The Langevin MCMC is used to draw statistical samples $y_{t}\sim\exp(-\phi)/Z.$ Langevin MCMC is based on discretizing the Langevin diffusion (van Kampen, 1992) and the updates are based solely on $\nabla\phi$ (the equations are coming next), therefore not knowing the normalizing constant $Z$ is of no concern.

(jumps) Given $y_{\tau}$, ideally an exact sample (MacKay, 2003) from $\exp(-\phi)/Z$, the jumps are made with the Bayes estimator of $X$. The spirit of empirical Bayes is fully present here: the jump samples are generated without knowing $\nabla f_{X}$ (or its approximation), and without running Langevin MCMC on its “rougher terrain”.

What emerges is a sampling algorithm in two parts. The Langevin MCMC is the engine that should run continously. By contrast, the jumps can be made at any time. That is:

Sample from $\exp(-\phi)/Z$ with Langevin MCMC,

where $\delta\ll 1$ is the step size and here $t$ is discrete time.

At any time $\tau$, use the Bayes estimator of $X$ to jump from $y_{\tau}$ to $\widehat{x}(y_{\tau})$,

We tested the algorithm on the handwritten digit database for $\sigma=0.3$ and $\sigma=0.15$. The results are shown in Figure 5 and 6 respectively. For $\sigma=0.3$ (in the regime of highly overlapping spheres), there was mixing between styles/classes. For $\sigma=0.15$ (in the regime of mostly non-overlapping spheres), samples stayed within a class while the styles changed.

Some side-by-side comparisons with denoising autoencoders are in order.

In walk-jump sampling, the key step is sampling from the density $\exp(-\phi)/Z$. The least-squares estimation of $X$—the jumps—are decoupled from Langevin MCMC. This is in contrast to the chain $X_{t}\rightarrow Y_{t}\rightarrow X_{t+1}$ constructed in generalized denoising autoencoders (Bengio et al., 2013b), which does not suffer from the problems we mentioned for the jumps.

Here, the learning objective is derived for any $\sigma$. By contrast, denoising autoencoders approximate the score function only in the limit $\sigma\rightarrow 0$ (Alain and Bengio, 2014). Generalized denoising autoencoders (Bengio et al., 2013b) and generative stochastic networks (Alain et al., 2016) were devised as a remedy for those $\sigma\rightarrow 0$ limitations. This is avoided altogether here.

In this work, there is a clear geometrical notion for large and small $\sigma$, which is based on the concentration of measure phenomenon. This picture is lacking in the references cited, and it is not clear which ranges of noise values should be considered there.

Most importantly, denoising autoencoders suffer from “limited parameterization” as expressed by Alain and Bengio (2014), and independently observed in (Saremi et al., 2018). To summarize, in denoising autoencoders, one must learn a curl-free encoder-decoder to be able to properly approximate the score function and this is problematic beyond one-hidden-layer architectures. In our work, this problem is avoided due to the energy parameterization (see Remark 7).

Neural empirical Bayes associative memory (NEBULA)

In this section, the neural empirical Bayes machinery is used to define a new notion of associative memory. Associative memory (also known as content-addressable memory) is a deep subject in neural networks with a rich history and with roots in psychology and neuroscience (long-term potentiation and Hebb’s rule). The depiction of this computation is even present in the arts and literature, championed by Marcel Proust, in the form of “stream of consciousness”, and the constant back and forth between perception and memory.

In 1982, Hopfield brought together earlier attempts to formulate associative memory, and showed that the collective computation of a system of neurons with symmetric weights, in the form of asynchronous updates of McCulloch-Pitts neurons (McCulloch and Pitts, 1943), minimize an energy function (Hopfield, 1982). The energy function was constructed in a Hebbian fashion. The associative memory was then formulated as the “flow” (the neurons were binary) to the local minima of the energy function.

However, Hopfield’s energy function is not the energy functionIn addition, Hopfield networks have severe limitations in memory capacity, addressed most recently in (Hillar and Tran, 2014; Krotov and Hopfield, 2016; Chaudhuri and Fiete, 2017).: it does not approximate (learn) the negative log probability density function of its stored memories. The Boltzmann machine (Hinton and Sejnowski, 1986) was developed in fact inspired by that observation, in which learning the energy function was achieved by introducing hidden units. Regarding the associative memory and the phase space flow, the problem is that in Boltzmann machines, hidden units need to be inferred first before having any notion of flow for the visible units. And inference is indeed computationally very expensive—“the curse of inference” (but not nearly as fundamental as the curse of dimensionality)—in probabilistic graphical models (Wainwright and Jordan, 2008; Koller and Friedman, 2009).

What we have achieved so far is to learn a function $\phi$ which approximates the energy function of $Y$, $\phi\approx-\log f$ (modulo a constant). The key here is that $\phi$ is a function (a “computer”) that computes $f(y)$ for any $y$; the “hidden units” in $\phi$ are not inferred. In other words, the hidden units (and the parameters) in $\phi$ are there solely for universal approximation, $-\nabla\phi\approx\nabla\log f.$ In Boltzmann machines, the hidden units are there to think! And this is the big difference. The definition below should be viewed against this backdrop.

We define the neural empirical Bayes associative memory (NEBULA), à la Hopfield (1982), as the flow to strict local minima of the energy function $\phi$. In continuous time, the memory dynamics is governed by the gradient flow:

NEBULA is identified by its attractors, the set of all the strict local minima of $\phi$:

What is so special about the local minima of $\phi$? Start with a single sample $X_{1}$,

For many samples, the learning objective is optimized such that $-\nabla\phi$ evaluated on $i$-spheres point to $X_{i}$ (see Figure 1c), but there are conflicts of interests, and these “conflicts” are the essence of $i$-sphere interactions, stated in Definition 3; in its practical summary, larger $\sigma$ means more $i$-sphere overlaps, and therefore larger “interaction couplings”. For NEBULA, the presence of $i$-sphere interactions implies that $X_{a}^{*}$ are not statistical samples from $f_{X}$. In other words, given $X_{i}\sim f_{X}$, $\mathscr{A}(X_{i})\not\sim f_{X}$. But what is the “right metric” $d(X_{i},\mathscr{A}(X_{i})),$ to measure the distance between $X_{i}\sim f_{X}$ and $\mathscr{A}(X_{i})$? The problem is that the flow is the gradient flow of $\phi\approx-\log f$ but the attractors are roughly speaking “close to $\mathscr{M}$”. The total length of the path taken by the gradient flow is a natural choice but we leave that for furture studies. Our expectation (based on $i$-sphere interactions) is that $d(X_{i},\mathscr{A}(X_{i}))$ will be larger for larger $\sigma$, visualized in Figures 7 and 8, with pronounced qualitative differences.

Emergence of “creative memories”

We finish the paper with some surprising results on the emergence of highly-structured non-digit memories, named “creative memories”. The experiments designed in this section are a natural continuation of the experiments from the previous section. In its summary, NEBULA is “well behaved” when it is initialized at $X_{i}\sim f_{X}$. But what if it is not initialized as such? Here, we explore $i$-sphere interactions (see Definition 5) from a different angle by initializing NEBULA at a random point,

The results in Figures 9 and 10 are the strongest evidence for $i$-sphere interactions as being a good abstraction. In viewing Figures 9 and 10, consider Figure 2, but with 60000 highly overlapping $i$-spheres, and $-\nabla\phi$ on each $i$-sphere that should “point towards” $X_{i}$ in expectation, where in addition, there are other complexities (touched upon in Remark 6) regarding the $i$-spheres themselves.

The images shown in Figures 9 and 10 are not statistical samples from $f_{X}$ nor $f_{Y}$.

They are not the result of mass aggregation of the kernel density estimator $\widehat{f}_{Y}$. Mode seeking is a deep topic in the kernel density literature around the mean shift algorithm (fukunaga1975estimation; Cheng, 1995; Comaniciu and Meer, 2002), also around the intriguing topic of “ghost modes” (Carreira-Perpiñán and Williams, 2003). Of course, the attractors reported here has not been reported in the kernel density literature (see (Carreira-Perpiñán, 2015) and the examples therein for MNIST).

The results presented are from all the random initializations in the unit hypercube from which we ran the algorithm, without any hand-picking!

Summary

In neural empirical Bayes, the smoothing of a random variable $X$ to $Y$ ($f=f_{X}*f_{N}$), denoising, and the (unnormalized) density estimation of the smoothed density $f$ were unified in a single machinery. Its inner workings were captured symbolically by $X\rightleftharpoons Y$, as well as by a Gedankenexperiment with an “experimenter” in the school of Robbins (1956) and Robbins and Monro (1951).In our presentation, the smoothing part was viewed from the angle of kernel density estimation, with an important distinction that the kernel bandwidth was a hyperparameter. But on reflection, this is not necessary. Fundamentally, setting aside the “neural energy function”, one can go on with our program by just knowing the (deep) machineries of empirical Bayes (Robbins, 1956) and stochastic gradients (Robbins and Monro, 1951). In this machinery, the energy function is parametrized with a neural network and SGD becomes the engine for learning, whose end result is captured by $\nabla\phi(\cdot,\theta^{*})\approx-\nabla\log f(\cdot)$. We proposed an approximative walk-jump sampling scheme which produces samples which are particularly appealing visually due to the denoising jump used. The energy function can further be used as an associative memory (NEBULA), which even has some “creative” capacities.

Acknowledgement

This work was supported by the Canadian Institute for Advanced Research, the Gatsby Charitable Foundation, and the National Science Foundation through grant IIS-1718991. We especially thank Bruno Olshausen, Eero Simoncelli, and Francis Bach for discussions.

References