# Gaussian-Bernoulli RBMs Without Tears

Renjie Liao 1, Simon Kornblith2, Mengye Ren3, David J. Fleet2,4,5, Geoffrey
Hinton2,4,5
University of British Columbia1, Google Research, Brain Team2,
New York University3, University of Toronto4, Vector Institute5
<EMAIL_ADDRESS><EMAIL_ADDRESS>
{skornblith, davidfleet<EMAIL_ADDRESS>Work done partially as a
visiting faculty researcher at Google Brain.

###### Abstract

We revisit the challenging problem of training Gaussian-Bernoulli restricted
Boltzmann machines (GRBMs), introducing two innovations. We propose a novel
Gibbs-Langevin sampling algorithm that outperforms existing methods like Gibbs
sampling. We propose a modified contrastive divergence (CD) algorithm so that
one can generate images with GRBMs starting from noise. This enables direct
comparison of GRBMs with deep generative models, improving evaluation
protocols in the RBM literature. Moreover, we show that modified CD and
gradient clipping are enough to robustly train GRBMs with large learning
rates, thus removing the necessity of various tricks in the literature.
Experiments on Gaussian Mixtures, MNIST, FashionMNIST, and CelebA show GRBMs
can generate good samples, despite their single-hidden-layer architecture. Our
code is released at: https://github.com/lrjconan/GRBM

##  1 Introduction

Restricted Boltzmann machines (RBMs) (Smolensky, 1986; Freund & Haussler,
1991; Hinton, 2002) are energy-based generative models with stochastic binary
units. A variant of Boltzmann machines (Ackley et al., 1985), they have a
bipartite graphical structure that enables efficient probabilistic inference,
and they can be stacked to form deep belief networks (DBNs) (Hinton &
Salakhutdinov, 2006; Bengio et al., 2006; Hinton et al., 2006). Gaussian-
Bernoulli RBMs (GRBMs) (Welling et al., 2004; Hinton & Salakhutdinov, 2006)
extend RBMs to model continuous data by replacing the binary visible units of
the RBM with Gaussian random variables.

GRBMs remain challenging to learn, however, despite many proposed
modifications to the model or training algorithm. For instance, Lee et al.
(2007) add a regularization term to encourage sparsely activated binary hidden
units. Krizhevsky et al. (2009) attribute the difficulties in learning to
high-frequency noise present in natural images. Factorized high-order terms
were introduced in (Ranzato & Hinton, 2010; Ranzato et al., 2010) to allow
GRBMs to explicitly learn the covariance structure among pixels. Nair & Hinton
(2010) suggest that binary hidden units are problematic, and proposed model
variants with real-valued hidden units. Cho et al. (2011a; 2013) advocate the
use of parallel tempering sampling (Earl & Deem, 2005), adaptive learning
rate, and enhanced gradient (Cho et al., 2011b) to improve GRBM learning.
Melchior et al. (2017) conclude that difficulties in GRBM training are due to
training algorithms rather than the model itself; they advocate the use of
gradient clipping, specialized weight initialization, and contrastive
divergence (CD) (Hinton, 2002) rather than persistent CD (Tieleman, 2008).
Upadhya & Sastry (2021) propose a stochastic difference of convex functions
programming (S-DCP) algorithm to replace CD in training GRBMs.

An important motivation for seeking to improve GRBM learning is so that a GRBM
can be used to convert real-valued data to stochastic binary data. This would
make it easy for researchers to explore novel ways of implementing stochastic
binary Boltzmann machines to model real-valued data. To that end, we propose
improved GRBM learning methods for image data. Specifically,

  * •

We propose a hybrid Gibbs-Langevin sampling algorithm that outperforms
predominant use of Gibbs sampling. To the best of our knowledge this is the
first use of Langevin sampling for GRBM training (with or without Metropolis
adjustment).

  * •

We propose a modified CD algorithm so that one can generate images with
learned GRBMs starting from Gaussian noise. This enables a fair and direct
comparison of GRBMs with deep generative models, something beyond the reach of
existing GRBM learning methods.

  * •

We show that the modified CD with gradient clipping is sufficient to train
GRBMs, thus removing the need for heuristics that have been crucial for
existing approaches.

  * •

We empirically show that GRBMs can generate good samples on Gaussian Mixtures,
MNIST, FashionMNIST, and CelebA, despite they have a single hidden layer.

##  2 Related Work

Learning the variances  Learning the variance of visible units in GRBMs is
necessary for generating sharp and realistic images. But small variances tend
to cause the energy function and its gradient to have large values, thus
making the stochastic gradient estimates returned by CD numerically unstable.
Most existing methods fix the variance (_e.g_., to one) to avoid this issue.
Krizhevsky et al. (2009); Cho et al. (2011a) consider learning the variance
using a smaller learning rate than for other parameters, obtaining much better
reconstruction, thus supporting the importance of learning variances. However,
many of the learned filters are still noisy and point-like. Melchior et al.
(2017) learn a shared variance across all visible units, yielding improved
performance, especially with large numbers of hidden units. In this work, we
learn one variance parameter per visible unit and achieve much lower learned
variances than existing methods, _e.g_., approximately $1e^{-5}$ on MNIST.

Stochastic gradient estimation and learning rate  Due to the intractable log
partition function of GRBMs, one often estimates the gradients of the log
likelihood w.r.t. parameters via Monte Carlo. Gibbs sampling is predominant in
CD learning due to its simplicity, but it mixes slowly in practice. This
yields noisy gradient estimates which often cause training instabilities and
prohibits using large learning rates. Cho et al. (2011a) explore parallel
tempering with adaptive learning rates to obtain better reconstruction. Cho et
al. (2013) propose enhanced gradients that are invariant to bit-flipping in
hidden units. Melchior et al. (2017) show that gradient clipping and special
weight initialization support robust CD learning with large learning rates. We
advocate Langevin MC to improve gradients, and validate that gradient clipping
does enable training with large learning rates.

Model capacity  Theis et al. (2011) empirically show that GRBMs are
outperformed even by simple mixture models in estimating likelihoods for image
data. Wang et al. (2012); Melchior et al. (2017) demonstrate that GRBMs can be
expressed as either a product of experts or a constrained Gaussian mixture in
the visible domain, hinting that GRBMs need more hidden units than the true
number of components to fit additive mixture densities well. Krause et al.
(2013); Gu et al. (2022) provide theoretical guarantees on GRBMs for universal
approximation of mixtures and smooth densities. Although this shows that GRBMs
are expressive, they do not lead directly to practical GRBM learning
algorithms.

Model Evaluation  Like many deep generative models, evaluating GRBMs is
difficult, as the log likelihood is intractable. To date, GRBMs have been
evaluated by visually inspecting reconstructed images, filters and hidden
activation (_i.e_., features), and sampled images during CD training.
Quantitative metrics include reconstruction errors, and error rates of post-
hoc trained classifiers on learned features. However, these metrics do not
necessarily indicate if GRBMs are good generative models Melchior et al.
(2017). Unlike existing work, we sample from learned GRBMs, starting from
Gaussian noise, enabling direct comparisons with other generative models,
qualitatively (visually inspecting samples) and quantitatively (_e.g_.,
Frechet Inception distance (FID) (Heusel et al., 2017).

##  3 Gaussian-Bernoulli Restricted Boltzmann Machines

A Gaussian-Bernoulli Restricted Boltzmann Machine (GRBM) (Welling et al.,
2004; Krizhevsky et al., 2009; Cho et al., 2011a; Melchior et al., 2017) is a
Markov Random Field (MRF) with continuous stochastic visible units and binary
stochastic hidden units. Denoting $N$ visible units as
$\mathbf{v}\in\mathbb{R}^{N}$ and $M$ hidden units as
$\mathbf{h}\in\\{0,1\\}^{M}$, the energy function associated with a GRBM is
defined to be

$\displaystyle
E_{\theta}(\mathbf{v},\mathbf{h})=\frac{1}{2}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)^{\top}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)-\left(\frac{\mathbf{v}}{\bm{\sigma}^{2}}\right)^{\top}W\mathbf{h}-\mathbf{b}^{\top}\mathbf{h}\,\,,$
(1)

with weight matrix $W\in\mathbb{R}^{N\times M}$, bias
$\bm{b}\in\mathbb{R}^{M}$, mean $\bm{\mu}\in\mathbb{R}^{N}$, and variance
$\bm{\sigma}^{2}\in\mathbb{R}_{+}^{N}$, where, unless stated otherwise,
$\frac{\mathbf{x}}{\mathbf{y}}$ denotes element-wise division between vectors
$\mathbf{x}$ and $\mathbf{y}$, as is convention in the GRBM literature. We
denote the set of learnable parameters as
$\theta=\\{W,\bm{b},\bm{\mu},\bm{\sigma}^{2}\\}$. To ensure the variance
remains non-negative during learning, we adopt a reparameterization, directly
learning $\log\bm{\sigma}^{2}$ rather than $\bm{\sigma}^{2}$ or $\bm{\sigma}$.
Finally, given the energy function, one can define the Boltzmann distribution,
over visible and hidden states, as

$\displaystyle p_{\theta}(\mathbf{v},\mathbf{h})$
$\displaystyle=\frac{1}{Z}\exp\left(-E_{\theta}(\mathbf{v},\mathbf{h})\right)\,,~{}~{}\mbox{where}~{}~{}~{}Z=\int_{-\infty}^{+\infty}\sum_{\mathbf{h}}\exp\left(-E_{\theta}(\mathbf{v},\mathbf{h})\right)\mathrm{d}\mathbf{v}$
(2)

is the normalization constant, which is intractable for even moderately large
$M$.

The underlying graphical model, like an RBM, is a bipartite graph with edges
only connecting visible units to hidden units. This entails conditional
independence of the form
$p(\mathbf{v}|\mathbf{h})=\prod_{i}p(\mathbf{v}_{i}|\mathbf{h})$ and
$p(\mathbf{h}|\mathbf{v})=\prod_{j}p(\mathbf{h}_{j}|\mathbf{v})$. One can also
derive the following conditional distributions for GRBMs,

$\displaystyle p(\mathbf{v}|\mathbf{h})$
$\displaystyle=\mathcal{N}\left(\mathbf{v}|W\mathbf{h}+\bm{\mu},\mathrm{diag}(\bm{\sigma}^{2})\right)$
(3) $\displaystyle p(\mathbf{h}_{j}=1|\mathbf{v})$
$\displaystyle=\left[\text{Sigmoid}\left(W^{\top}\frac{\mathbf{v}}{\bm{\sigma}^{2}}+\mathbf{b}\right)\right]_{j},$
(4)

where
$\mathcal{N}\left(\mathbf{v}|W\mathbf{h}+\bm{\mu},\mathrm{diag}(\bm{\sigma}^{2})\right)$
is the multivariate Gaussian distribution with mean $W\mathbf{h}+\bm{\mu}$,
and the diagonal covariance matrix $\mathrm{diag}(\bm{\sigma}^{2})$. Here,
$\text{Sigmoid}(\mathbf{x})=1/(1+\exp(-\mathbf{x}))$ is applied to the vector
$\mathbf{x}$ in an element-wise manner, and $[\cdot]_{j}$ denotes the $j$-th
element of the corresponding vector.

Given the Boltzmann distribution, one can derive the marginal distribution
over visible units, _i.e_.,

$\displaystyle p(\mathbf{v})$
$\displaystyle=\frac{1}{Z}\exp\left(-\tilde{E}_{\theta}(\mathbf{v})\right),$
(5) $\displaystyle\tilde{E}_{\theta}(\mathbf{v})$
$\displaystyle=\frac{1}{2}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)^{\top}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)-\text{Softplus}\left(W^{\top}\frac{\mathbf{v}}{\bm{\sigma}^{2}}+\mathbf{b}\right)^{\top}\textbf{1}~{},$

where $\text{Softplus}(\mathbf{x})=\log(1+\exp(\mathbf{x}))$ is applied in an
element-wise manner, and 1 is the all-one vector of size $M$. We call
$\tilde{E}_{\theta}(\mathbf{v})$ the marginal energy to distinguish it from
the GRBM energy in Eq. (1). We leave the derivation to Appendix A.1. As shown
in Melchior et al. (2017), one can also rewrite the marginal distribution
$p(\mathbf{v})$ as a constrained Gaussian mixture.

Algorithm 1 Langevin Sampling for GRBMs

1:Input: $\mathbf{v}^{(0)}$, step size $\alpha_{0}$, total step $T$, burn-in
step $\tilde{T}$, adjust step $\eta$

2:For $t=1,\dots,T$

3: $\alpha_{t}=\text{CosineScheduler}(t,T,\alpha_{0})$

4:
$\mathbf{v}=\mathbf{v}^{(t-1)}-\alpha_{t}\frac{\partial\tilde{E}(\mathbf{v}^{(t-1)})}{\partial\mathbf{v}}+\sqrt{2\alpha_{t}}\bm{\xi}_{t}$
, $\bm{\xi}_{t}\sim\mathcal{N}(0,I)$ $\triangleright$ Use marginal energy in
Eq. (5)

5: If $t<=\eta$ or $\left(t>\eta\text{ and
}u\sim\mathcal{U}(0,1)<A(\mathbf{v},\mathbf{v}^{(t-1)})\right)$

6: $\mathbf{v}^{(t)}=\mathbf{v}$

7: Else

8: $\mathbf{v}^{(t)}=\mathbf{v}^{(t-1)}$

9:Return: $\\{\mathbf{v}^{(\tilde{T}+1:T)}\\}$ $\triangleright$ $i:j$ indexes
consecutive samples from $i$-th to $j$-th

###  3.1 Inference

When performing probabilistic inference, _e.g_., computing the marginal
distribution or the maximum a posterior (MAP) estimation, one often chooses
between variational inference (Hinton & Van Camp, 1993; Jordan et al., 1999)
and Markov chain Monte Carlo (MCMC) methods (Neal, 1993; Andrieu et al.,
2003). We focus on MCMC as common variational methods have been less effective
with RBMs and GRBMs (Gabrié et al., 2015; Takahashi & Yasuda, 2016). From the
generative modelling perspective, we wish to draw samples of visible units
during inference. There are two natural approaches to this: 1) sample from the
joint distribution in Eq. (2) and discard the samples of hidden units, or 2)
directly sample from the marginal distribution.

Gibbs sampling (Geman & Geman, 1984) is perhaps the predominant approach, due
to its simplicity. In the context of GRBMs, one alternates between sampling
hidden units given visible units, and sampling visible units given hidden
units. This produces samples from the joint distribution in Eq. (2). The
detailed Gibbs sampling algorithm is given in Appendix A.2.

##### Langevin Sampling

Langevin Monte Carlo (Grenander & Miller, 1994; Roberts & Tweedie, 1996;
Welling & Teh, 2011) is a class of MCMC methods that generate samples from a
probability distribution of continuous random variables by simulating Langevin
dynamics. Since GRBMs are hybrid graphical models, _i.e_., comprising
continuous and discrete random variables, we have at least two ways to
leverage Langevin sampling. One is to directly apply Langevin sampling to the
marginal distribution of visible units in Eq. (5). Suppose at time step $t-1$,
we have sample $\mathbf{v}_{t-1}$ and want to draw a new sample
$\mathbf{v}_{t}$. The proposal distribution corresponding to one-step Langevin
dynamics is given by

$\displaystyle
q(\mathbf{v}|\mathbf{v}^{(t-1)})=\mathcal{N}\left(\mathbf{v}\middle|\mathbf{v}^{(t-1)}-\alpha_{t}\frac{\partial\tilde{E}(\mathbf{v}^{(t-1)})}{\partial\mathbf{v}},2\alpha_{t}I\right),$
(6)

where the gradient the of marginal energy $\tilde{E}$ w.r.t. the visible units
is given in Appendix A.3. If we use the Metropolis-Hastings algorithm to
accept or reject proposed samples, the acceptance probability of a proposal
$\mathbf{v}_{t}$, given the previous state, $\mathbf{v}_{t-1}$, is (see
Appendix A.3 for derivation):

$\displaystyle
A(\mathbf{v}^{(t)},\mathbf{v}^{(t-1)})=\min\left(\\!1,\,\frac{\exp\left(-\tilde{E}_{\theta}(\mathbf{v}^{(t)})-\frac{1}{4\alpha_{t}}\left\|\mathbf{v}^{(t-1)}-\mathbf{v}^{(t)}+\alpha_{t}\frac{\partial\tilde{E}(\mathbf{v}^{(t)})}{\partial\mathbf{v}}\right\|^{2}\right)}{\exp\left(-\tilde{E}_{\theta}(\mathbf{v}^{(t-1)})-\frac{1}{4\alpha_{t}}\left\|\mathbf{v}^{(t)}-\mathbf{v}^{(t-1)}+\alpha_{t}\frac{\partial\tilde{E}(\mathbf{v}^{(t-1)})}{\partial\mathbf{v}}\right\|^{2}\\!\right)}\right).$
(7)

Alg. 1 shows the Metropolis-adjusted Langevin Algorithm (MALA) for the
marginal GRBM. Compared to generic MALA, it also includes an extra
hyperparameter, namely, the adjust step $\eta$. If $\eta$ is set to $0$, then
we perform a Metropolis adjustment at every sampling step, as prescribed in
the generic MALA. If $\eta$ is set to $K>0$, then we skip the Metropolis
adjustment for the first $K$ steps. The adjust step effectively controls a
trade-off between sampling accuracy111Here sampling accuracy means the
closeness between the underlying distribution of samples and the target
distribution measured in, _e.g_., total variation or Wasserstein distances.
and computational efficiency. Since we do not hope to see Gaussian noise in
our final-sampled images (_i.e_., beyond the level of intrinsic noise in the
observations), it is beneficial to decay the noise level, as in score-based
models (Song & Ermon, 2019). For certain step-size-annealing schedules and
energy functions, there are theoretical guarantees on the convergence of
Langevin sampling (Durmus & Moulines, 2019). For simplicity, we use the cosine
scheduler and find it works well in practice. More details about the scheduler
are provided in Appendix A.3.

##### Gibbs-Langevin Sampling

Algorithm 2 Gibbs-Langevin Sampling for GRBMs

1:Input: $\mathbf{v}^{(0)}$, $\mathbf{h}^{(0)}$, step size $\alpha_{0}$, total
step $T$, burn-in step $\tilde{T}$, adjust step $\eta$, Langevin step $K$

2:Function Langevin($\tilde{\mathbf{v}}^{(0)}$, $\mathbf{h}$, $\alpha_{0}$,
K):

3: For $k=1,\dots,K$

4: $\alpha_{k}=\text{CosineScheduler}(k,K,\alpha_{0})$

5:
$\tilde{\mathbf{v}}^{(k)}=\tilde{\mathbf{v}}^{(k-1)}-\alpha_{k}\frac{\partial
E(\tilde{\mathbf{v}}^{(k-1)},\mathbf{h})}{\partial\mathbf{v}}+\sqrt{2\alpha_{k}}\bm{\xi}_{k}$
, $\bm{\xi}_{k}\sim\mathcal{N}(0,I)$

6:Return $\tilde{\mathbf{v}}^{(K)}$

7:

8:For $t=1,\dots,T$

9:
$\mathbf{v}=\text{Langevin}(\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)},\alpha_{0},K)$

10: $\mathbf{h}\sim p(\mathbf{h}|\mathbf{v})$

11: If $t<=\eta$ or $\left(t>\eta\text{ and
}u\sim\mathcal{U}(0,1)<\tilde{A}\left((\mathbf{v},\mathbf{h}),(\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})\right)\right)$

12: $\mathbf{v}^{(t)},\mathbf{h}^{(t)}=\mathbf{v},\mathbf{h}$

13: Else

14: $\mathbf{v}^{(t)},\mathbf{h}^{(t)}=\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)}$

15:Return: $\\{(\mathbf{v}^{(\tilde{T}+1:T)},\mathbf{h}^{(\tilde{T}+1:T)})\\}$

We also introduce a new hybrid sampler for GRBMs (see Alg. 2). Like the Gibbs
sampler, it alternates between sampling hidden units conditioned on visible
units, and sampling visible units given the hidden units. Unlike generic
Gibbs, which directly samples from the Gaussian
$p(\mathbf{v}|\mathbf{h}^{(t)})$, we instead use Langevin MC to sample the
continuous visible units given the hidden state. The use of Langevin MC may
seem unnecessary because the Gaussian conditional permits a one-step sampling
algorithm. The subtlety comes from the fact that the finite-step Langevin
sampler explicitly depends on the initial sample. Specifically, the proposal
distribution of one complete outer-loop step in Alg. 2, _e.g_., at iteration
$t-1$, can be expressed as

$\displaystyle q(\mathbf{v},\mathbf{h}|\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})$
$\displaystyle\,=\,q(\mathbf{h}|\mathbf{v})~{}q(\mathbf{v}|\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})~{},$
(8)

where $q(\mathbf{h}|\mathbf{v})$ is given by Eq. (4), and
$q(\mathbf{v}|\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})$ is the proposal
distribution of a K-step Langevin sampler (_i.e_. from the inner loop). This
proposal distribution explicitly depends on the initial visible sample,
$\mathbf{v}^{(t-1)}$ from iteration $t\\!-\\!1$. By contrast, the generic
Gibbs sampler does not have such dependence, _i.e_.,
$q(\mathbf{v}|\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})=q(\mathbf{v}|\mathbf{h}^{(t-1)})$.
This dependence allows us to construct a persistent Markov chain in the space
of visible units. Moreover, the Langevin sampler leverages the informative
gradient of log density whereas Gibbs sampler does not. We find that our new
sampler performs significantly better than the vanilla Gibbs sampler in
practice.

The Metropolis adjustment for these Gibbs-Langevin proposals is, however,
somewhat more involved. Following Alg. 2, with
$\tilde{\mathbf{v}}^{(0)}=\mathbf{v}^{(t-1)}$ and
$\tilde{\mathbf{v}}^{(K)}=\mathbf{v}$, by marginalizing out the intermediate
states on the Markov chain, we obtain the proposal

$\displaystyle
q(\tilde{\mathbf{v}}^{(K)}|\tilde{\mathbf{v}}^{(0)},\mathbf{h}^{(t-1)})=\int\cdots\int\left(\prod_{k=1}^{K}q(\tilde{\mathbf{v}}^{(k)}|\tilde{\mathbf{v}}^{(k-1)},\mathbf{h}^{(t-1)})\right)\mathrm{d}\tilde{\mathbf{v}}^{(1)}\cdots\mathrm{d}\tilde{\mathbf{v}}^{(K-1)}.$
(9)

The integrand in Eq. 9 comprises $K$ one-step Langevin updates, each of which
is given by

$\displaystyle q(\mathbf{v}|\tilde{\mathbf{v}}^{(k-1)},\mathbf{h}^{(t-1)})$
$\displaystyle=\mathcal{N}\left(\mathbf{v}\middle|\tilde{\mathbf{v}}^{(k-1)}-\alpha_{k}\frac{\partial
E(\tilde{\mathbf{v}}^{(k-1)},\mathbf{h}^{(t-1)})}{\partial\mathbf{v}},\,2\alpha_{k}I\right)~{},$
(10)

for which the energy gradient is given in Appendix A.4. Although the multiple
integral in Eq. (9) appears intractable, one can use reparameterization to
derive the following analytical form,

$\displaystyle\\!q(\tilde{\mathbf{v}}^{(K)}|\tilde{\mathbf{v}}^{(0)},\mathbf{h}^{(t-1)})=\mathcal{N}\left(\bm{\beta}_{0}\tilde{\mathbf{v}}^{(0)}+\left(\sum_{k=1}^{K}\bm{\beta}_{k}\alpha_{k}\right)\frac{\bm{\mu}+W\mathbf{h}^{(t-1)}}{\bm{\sigma}^{2}},\,\mathrm{diag}\left(\sum_{k=1}^{K}2\alpha_{k}\bm{\beta}_{k}^{2}\right)\right),$
(11)

where
$\bm{\beta}_{k}=\prod_{j=k+1}^{K}\left(\bm{1}-\frac{\alpha_{j}}{\bm{\sigma}^{2}}\right)$,
$\forall k\in\\{0,\dots,K-1\\}$ and $\bm{\beta}_{K}=\bm{1}$. Based on this
result, one can show that the acceptance probability for the Metropolis
adjustment is

$\displaystyle\tilde{A}((\mathbf{v}^{(t)},\mathbf{h}^{(t)}),(\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)}))=$
$\displaystyle\min\left(1,\frac{\exp\left(-E_{\theta}(\mathbf{v}^{(t)},\mathbf{h}^{(t)})-\left\|\frac{\mathbf{v}^{(t-1)}-\bm{\beta}_{0}\mathbf{v}^{(t)}-\bm{a}(\bm{\mu}+W\mathbf{h}^{(t)})}{\sqrt{2}\tilde{\bm{\sigma}}}\right\|^{2}\right)q(\mathbf{h}^{(t-1)}|\mathbf{v}^{(t-1)})}{\exp\left(-E_{\theta}(\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})-\left\|\frac{\mathbf{v}^{(t)}-\bm{\beta}_{0}\mathbf{v}^{(t-1)}-\bm{a}(\bm{\mu}+W\mathbf{h}^{(t-1)})}{\sqrt{2}\tilde{\bm{\sigma}}}\right\|^{2}\right)q(\mathbf{h}^{(t)}|\mathbf{v}^{(t)})}\right),$
(12)

where
$q(\mathbf{h}_{j}=1|\mathbf{v})=\left[\text{Sigmoid}\left(W^{\top}\frac{\mathbf{v}}{\bm{\sigma}^{2}}+\mathbf{b}\right)\right]_{j}$,
$\bm{a}=\frac{\sum_{k=1}^{K}\bm{\beta}_{k}\alpha_{k}}{\bm{\sigma}^{2}}$, and
$\tilde{\bm{\sigma}}^{2}=\sum_{k=1}^{K}2\alpha_{k}\bm{\beta}_{k}^{2}$. We
leave derivations to Appendix A.4.

###  3.2 Learning

To learn GRBMs, we maximize the log likelihood of the observed data using
stochastic gradient-based methods, _e.g_., contrastive divergence (CD).
Depending on whether we use the joint (Eq. (2)) or the marginal (Eq. (5))
distribution, we have two possible gradient estimators.

Algorithm 3 Modified CD Learning Algorithm for GRBMs with Joint Density

1:Input: CD-step $K$, burn-in step $M$, learning Rate $\eta$, Langevin step
size $\alpha_{0}$, SGD step $T$

2:For $t=1,\cdots,T$

3: $\mathbf{v}^{+}=\mathbf{v}_{\text{data}}$

4: $\mathbf{h}^{+}\sim p(\mathbf{h}|\mathbf{v}^{+})$

5: $\nabla\theta^{+}=\left\langle\frac{\partial
E(\mathbf{v}^{+},\mathbf{h}^{+})}{\partial\theta}\right\rangle_{d}$
$\triangleright$ Compute Positive Gradient

6: $\mathbf{v}^{-}_{0}\sim\mathcal{N}(\mathbf{0},I)$, $\mathbf{h}^{-}_{0}\sim
p(\mathbf{h}|\mathbf{v}^{-}_{0})$ $\triangleright$ Modified CD to start with
noise

7:
$\\{\mathbf{v}^{-}_{M:K},\mathbf{h}^{-}_{M:K}\\}\sim\text{Sampler}(\mathbf{v}^{-}_{0},\mathbf{h}^{-}_{0},\alpha_{0}\bar{\bm{\sigma}})$
$\triangleright$ Alg. 5 or Alg. 2, $\bar{\bm{\sigma}}$ is current mean
variance

8: $\nabla\theta^{-}=\left\langle\frac{\partial
E(\mathbf{v}^{-},\mathbf{h}^{-})}{\partial\theta}\right\rangle_{m}$
$\triangleright$ Compute Negative Gradient

9: $\theta=\theta-\eta(\nabla\theta^{+}-\nabla\theta^{-})$ $\triangleright$
Compute Update

10:Return $\theta$

##### Learning with the Joint Distribution

When optimizing the GRBM with the joint distribution, one can express the
general form of the gradient of the log likelihood w.r.t. parameters $\theta$
as

$\displaystyle\nabla\theta=\left\langle-\frac{\partial
E_{\theta}(\mathbf{v},\mathbf{h})}{\partial\theta}\right\rangle_{d}-\left\langle-\frac{\partial
E_{\theta}(\mathbf{v},\mathbf{h})}{\partial\theta}\right\rangle_{m}.$ (13)

Here, following the notation in the RBM literature, we denote expectation
under the _data distribution_ , _i.e_.,
$p_{\theta}(\mathbf{h}|\mathbf{v})\,p_{\text{data}}(\mathbf{v})$, as
$\left\langle\cdot\right\rangle_{d}=\mathbb{E}_{p_{\theta}(\mathbf{h}|\mathbf{v})p_{\text{data}}(\mathbf{v})}\left[\cdot\right]$.
Similarly, we denote the expectation under the _model distribution_ ,
$p_{\theta}(\mathbf{v},\mathbf{h})$ as
$\left\langle\cdot\right\rangle_{m}=\mathbb{E}_{p_{\theta}(\mathbf{v},\mathbf{h})}\left[\cdot\right]$.
The expected gradients under the data and model distributions are called
_positive_ and _negative_ gradients respectively. Based on Eq. (13), we can
formulate the gradients of specific parameters as follows,

$\displaystyle\nabla W_{ij}$
$\displaystyle=\left\langle\frac{\mathbf{v}_{i}}{\bm{\sigma}_{i}^{2}}\mathbf{h}_{j}\right\rangle_{d}-\left\langle\frac{\mathbf{v}_{i}}{\bm{\sigma}_{i}^{2}}\mathbf{h}_{j}\right\rangle_{m}$
(14) $\displaystyle\nabla\mu_{i}$
$\displaystyle=\left\langle\frac{\mathbf{v}_{i}-\mu_{i}}{\bm{\sigma}_{i}^{2}}\right\rangle_{d}-\left\langle\frac{\mathbf{v}_{i}-\mu_{i}}{\bm{\sigma}_{i}^{2}}\right\rangle_{m}$
(15) $\displaystyle\nabla\log\bm{\sigma}_{i}^{2}$
$\displaystyle=\left\langle\frac{(\mathbf{v}_{i}-\mu_{i})^{2}}{2\bm{\sigma}_{i}^{2}}-\frac{\sum_{j}\mathbf{v}_{i}W_{ij}\mathbf{h}_{j}}{\bm{\sigma}_{i}^{2}}\right\rangle_{d}-\left\langle\frac{(\mathbf{v}_{i}-\mu_{i})^{2}}{2\bm{\sigma}_{i}^{2}}-\frac{\sum_{j}\mathbf{v}_{i}W_{ij}\mathbf{h}_{j}}{\bm{\sigma}_{i}^{2}}\right\rangle_{m}$
(16) $\displaystyle\nabla\mathbf{b}_{i}$
$\displaystyle=\langle\mathbf{h}_{i}\rangle_{d}-\langle\mathbf{h}_{i}\rangle_{m}~{}.$
(17)

Since the expectations in these gradients are generally intractable, we use
Monte Carlo methods to approximate them. To sample from the joint density, we
can use Gibbs or Gibbs-Langevin samplers as described in Sec. 3.1. The overall
learning algorithm is outlined in Alg. 3. An important detail is that we
multiply the initial Langevin step size by the average variance at each
gradient update step and then feed it to the sampler. Since the variance is
decreasing (the energy function and its gradient are increasing) as learning
goes on, keeping the step size roughly invariant to such scaling would make
the sampling more effective.

##### Learning with the Marginal Distribution

Now we turn to learning the model under the marginal distribution in Eq. (5).
Since we have the marginal distribution of visible units, we can directly get
the gradients of log likelihood w.r.t. model parameters $\theta$ as,

$\displaystyle\nabla\theta=\left\langle-\frac{\partial\tilde{E}_{\theta}(\mathbf{v})}{\partial\theta}\right\rangle_{d}-\left\langle-\frac{\partial\tilde{E}_{\theta}(\mathbf{v})}{\partial\theta}\right\rangle_{m}.$
(18)

Since the gradient
$\frac{\partial\tilde{E}_{\theta}(\mathbf{v})}{\partial\theta}$ does not
depend on $\mathbf{h}$ anymore, we have

$\displaystyle\left\langle-\frac{\partial\tilde{E}_{\theta}(\mathbf{v})}{\partial\theta}\right\rangle_{d}=\mathbb{E}_{p_{\text{data}}(\mathbf{v})}\left[-\frac{\partial\tilde{E}_{\theta}(\mathbf{v})}{\partial\theta}\right],\qquad\left\langle-\frac{\partial\tilde{E}_{\theta}(\mathbf{v})}{\partial\theta}\right\rangle_{m}=\mathbb{E}_{p_{\theta}(\mathbf{v})}\left[-\frac{\partial\tilde{E}_{\theta}(\mathbf{v})}{\partial\theta}\right].$
(19)

Based on above results, we can work out the detailed gradients which are the
same as those in Eq. (14) to Eq. (17) but with $\mathbf{h}$ replaced with
$\text{Sigmoid}\left(W^{\top}\frac{\mathbf{v}}{\bm{\sigma}^{2}}+\mathbf{b}\right)$.
More details are left to Appendix A.5. We use the Langevin sampler in Sec. 3.1
to sample from the marginal density to approximate the intractable
expectation. The overall learning algorithm is outlined in Alg. 4.

Algorithm 4 Modified CD Learning Algorithm for GRBMs with Marginal Density

1:Input: CD-step $K$, burn-in step $M$, learning rate $\eta$, Langevin step
size $\alpha_{0}$, SGD step $T$

2:For $t=1,\cdots,T$

3: $\mathbf{v}^{+}=\mathbf{v}_{\text{data}}$

4:
$\nabla\theta^{+}=\left\langle\frac{\partial\tilde{E}(\mathbf{v}^{+})}{\partial\theta}\right\rangle_{d}$
$\triangleright$ Compute Positive Gradient

5: $\mathbf{v}^{-}_{0}\sim\mathcal{N}(\mathbf{0},I)$ $\triangleright$ Modified
CD to start with noise

6:
$\\{\mathbf{v}^{-}_{i}|i=M,\cdots,K\\}\sim\text{Sampler}(\mathbf{v}^{-}_{0},\alpha_{0}\bar{\bm{\sigma}})$
$\triangleright$ Alg. 1, $\bar{\bm{\sigma}}$ is current mean variance

7:
$\nabla\theta^{-}=\left\langle\frac{\partial\tilde{E}(\mathbf{v}^{-})}{\partial\theta}\right\rangle_{m}$
$\triangleright$ Compute Negative Gradient

8: $\theta=\theta-\eta(\nabla\theta^{+}-\nabla\theta^{-})$ $\triangleright$
Compute Update

9:Return $\theta$

##### Modified Contrastive Divergence

The above two learning algorithms resemble CD if one ignores the specific
sampler used. There exists a subtle yet important difference however. For most
deep generative models one generates samples starting from noise. But this
does not work well for models trained with CD, where sampling starts from
observed data. This discrepancy of the starting sample between training and
testing would be a significant issue if the Markov chain does not mix
sufficiently quickly. We therefore modify CD by running two Markov chains to
collect samples for positive and negative gradients respectively. The positive
Markov chain is the same as in CD, _i.e_., starting from observed data. The
negative Markov chain now starts from a sample of standard Normal noise rather
than the reconstructed data222The reconstructed data is typically obtained by
running one complete step of Gibbs sampler from the observed data, thus being
highly likely close to observed data.. Since the positive chain starting from
data will usually stay close to the data distribution, this modification
pushes the negative Markov chain, starting from noise, toward the data
distribution. Moreover, the discrepancy between training and testing ceases to
be important as we can start from standard Normal noise while sampling from
the learned model.

##  4 Experiments

We examine the empirical behavior of our new GRBM algorithms on benchmark
image datasets, namely, MNIST, Fashion-MNIST (Xiao et al., 2017), and CelebA
(Liu et al., 2015).

|  |  |  |  |
---|---|---|---|---|---
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
(a) | (b) | (c) | (d) | (e) | (f)
Figure 1:  Density modelling using GRBMs on data from a Gaussian mixtures with
isotropic (rows 1 and 2) and anisotropic variances (rows 3 and 4). Rows 1 and
3 show normalized GMM densities and (unnormalized) negative energy values for
GRBMs. Rows 2 and 4 show samples drawn under different models and methods;
_i.e_., (a) Ground Truth; (b) Gibbs; (c) Langevin wo. Adjust; (d) Langevin w.
Adjust; (e) Gibbs-Langevin wo. Adjust; (f) Gibbs-Langevin w. Adjust.

##### Implementation Details

We found that training with modified CD alone occasionally diverges,
necessitating careful tuning of the learning rate. However, adding gradient
clipping (_e.g_., clip gradient norm to 10) enables stable training with all
aforementioned sampling methods. We therefore set learning rate to $0.01$ for
all experiments. Such a large learning rate almost never works in the
literature. Melchior et al. (2017) used gradient clipping and similarly large
learning rates, but they had to set the learning rate for the variances 100
times smaller than that for the weights and biases during CD training. But
thanks to the modified CD and gradient clipping, we found this special
treatment of variances is unnecessary. We do not use momentum, weight decay,
PCD, or other tricks.

###  4.1 Modeling Gaussian Mixture Densities

We first evaluate density modelling by GRBMs when the data density is known,
_i.e_., Gaussian mixture models (GMMs) in our case. This is challenging for
GRBMs as the marginal distribution of visible units of GRBMs is essentially a
constrained Gaussian mixture, _i.e_., the weights of mixture components depend
on one another (Melchior et al., 2017). As such, the mixture components in
GRBMs can not be freely placed in the visible domain so one actually needs
more hidden units than the log of the number of mixture components to fit GMMs
well. We consider the 2D case for simplicity and better visibility. We
generate 1,000 samples from two types (isotropic and anisotropic variances) of
GMMs with $3$ components as shown in Fig. 1, and learn GRBMs using our
modified CD with different sampling algorithms, from which we can draw
samples. Here all samplers run for 100 steps during both CD training and
testing (see Appendix B.1 for more detail). Density plots and samples are
shown in Fig. 1. Notice that Gibbs manages to recover the three modes in the
isotropic case but fails in the anisotropic case. Both Langevin and Gibbs-
Langevin sampling collapse when the adjustment is absent. We believe the
cosine step size schedule contributes to the collapse as it removes more
stochasticity of Langevin dynamics with small step sizes, thus making sampling
more similar to gradient descent. But as we will see later, in image
modelling, this may not be so severe; there are more modes so that the
sampling may collapse to different modes, and the diversity of images remains
acceptable. Finally, both Langevin and Gibbs-Langevin do recover all three
modes with the adjustment, which shows the adjustment helps the mixing in this
synthetic case.

###  4.2 Image Generation

We learn GRBMs to fit image datasets including MNIST, FashionMNIST, and
CelebA. To the best of our knowledge, this is the first time that GRBMs have
been shown to (unconditionally) generate good images. We provide the ablation
study in Appendix B.2 and more results in Appendix B.3.

|  |  |  |  |  |  |
---|---|---|---|---|---|---|---
|  |  |  |  |  |  |
|  |  |  |  |  |  |
|  |  |  |  |  |  |
|  |  |  |  |  |  |
|  |  |  |  |  |  |

Figure 2: Intermediate samples from Gibbs-Langevin sampling.

Methods | FID
---|---
VAE | 16.13
2sVAE (Dai & Wipf, 2019) | 12.60
PixelCNN++ (Salimans et al., ) | 11.38
WGAN (Arjovsky et al., 2017) | 10.28
NVAE (Vahdat & Kautz, 2020) | 7.93
GRBMs |
Gibbs | 47.53
Langevin wo. Adjust | 43.80
Langevin w. Adjust | 41.24
Gibbs-Langevin wo. Adjust | 17.49
Gibbs-Langevin w. Adjust | 19.27

Table 1: Results on MNIST dataset.

##### MNIST

We train GRBMs with hidden size 4096 and 100 sampling steps on MNIST. We
compare FID scores of GRBMs with other deep generative models in Table 1. From
the table, we can see that Gibbs-Langevin family works significantly better
than the Langevin family. The Metropolis adjustment improves Langevin slightly
but degrades Gibbs-Langevin slightly, which is different from what we observed
on synthetic data. This is likely because the image distribution is so
complicated (_e.g_., having significantly more modes) that the adjustment
rejects proposed moves more frequently than before. Some sophisticated
strategy may be needed to increase the acceptance probability. Nevertheless,
GRBMs trained with Gibbs-Langevin without adjustment achieve FID scores
comparable to other deep generative models, which is impressive given the
single-hidden-layer architecture. The learning curve of (natural) log variance
is shown in Fig. 3(a). The learned variance converges to around $1e^{-5}$
which is significantly smaller than those reported in the literature. The
learned filters are shown in Fig. 3(b). Although some point-like filters still
exist, stroke-like filters are common, thus indicating GRBMs indeed learn
meaningful features. We show samples drawn from the best GRBM in Fig. 3(c).
The intermediate samples from Gibbs-Langevin are shown in Fig. 2. Since Gibbs-
Langevin without adjustment works the best, we use it for remaining
experiments.

##### FashionMNIST

We then train GRBMs on FahsionMNIST which is more challenging than MNIST. We
set hidden size to 10,000 and the sampling step to 100. Samples drawn from
learned GRBMs are shown in Fig. 4(a). GRBMs successfully learn the shapes of
clothes, shoes, bags, and so on. However, they fail to capture fine textures.
Since many images in this dataset look similar in shape but differ in texture,
the resulting samples look similar to each other.

(a)

(b)

(c)

Figure 3: (a) Learning curve of (natural) log variances, (b) learned filters,
and (c) samples on MNIST.

##### CelebA

Last, we consider the even more challenging CelebA dataset. In particular, we
explore two versions of this dataset: 1) CelebA-32 where we center-crop (140
$\\!\times\\!$ 140) and downsample images to 32 $\\!\times\\!$ 32; 2)
CelebA-2K-64 where randomly select 2,000 images from the original CelebA and
apply the same center crop and downsampling to 64 $\\!\times\\!$ 64\. We set
hidden size to 10,000 and explore the number of 100 and 200 sampling steps.
Generated samples are shown in Fig. 4(b) and 4(c). From the figure, we can see
that GRBMs can learn to generate reasonably good face images.

(a)

(b)

(c)

Figure 4:  Samples from GRBMs on (a) FashionMNIST, (b) CelebA-32, and (c)
CelebA-2K-64.

##  5 Conclusion

In this paper, we revisit learning Gaussian-Bernoulli restricted Boltzmann
machines. We investigate Langevin Monte Carlo and propose a novel Gibbs-
Langevin sampling method. Furthermore, we modify the contrastive divergence
(CD) algorithm so that one can sample data from learned GRBMs starting from
noise. Modified CD along with gradient clipping enables robust training of
GRBMs with large learning rates. Finally, we show that GRBMs can
unconditionally generate images with good qualities, despite its single-
hidden-layer architecture. In the future, it would be beneficial to extend the
current GRBMs to convolutional GRBMs which should be able to learn better
localized filters. Meanwhile, it would be interesting to explore Gaussian deep
belief networks (GDBNs), which are deeper than GRBMs and should be superior.
At last, investigating our Gibbs-Langevin sampling for hybrid deep energy
based models could be a fruitful direction.

#### Acknowledgments

This work was funded, in part, by the NSERC Discovery Grant. Resources used in
preparing this research were provided, in part, by Google, Vector Institute,
and UBC.

## References

  * Ackley et al. (1985) David H Ackley, Geoffrey E Hinton, and Terrence J Sejnowski.  A learning algorithm for boltzmann machines.  _Cognitive science_ , 9(1):147–169, 1985.
  * Andrieu et al. (2003) Christophe Andrieu, Nando De Freitas, Arnaud Doucet, and Michael I Jordan.  An introduction to mcmc for machine learning.  _Machine Learning_ , 50(1):5–43, 2003.
  * Arjovsky et al. (2017) Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein generative adversarial networks.  In _International Conference on Machine Learning_ , pp. 214–223. PMLR, 2017.
  * Bengio et al. (2006) Yoshua Bengio, Pascal Lamblin, Dan Popovici, and Hugo Larochelle.  Greedy layer-wise training of deep networks.  _Advances in neural Information Processing Systems_ , 19, 2006.
  * Cho et al. (2013) Kyung Hyun Cho, Tapani Raiko, and Alexander Ilin.  Gaussian-bernoulli deep boltzmann machine.  In _The 2013 International Joint Conference on Neural Networks (IJCNN)_ , pp. 1–7. IEEE, 2013.
  * Cho et al. (2011a) KyungHyun Cho, Alexander Ilin, and Tapani Raiko.  Improved learning of gaussian-bernoulli restricted boltzmann machines.  In _International Conference on artificial Neural networks_ , pp. 10–17. Springer, 2011a.
  * Cho et al. (2011b) KyungHyun Cho, Tapani Raiko, and Alexander Ilin.  Enhanced gradient and adaptive learning rate for training restricted boltzmann machines.  In _ICML_ , 2011b.
  * Dai & Wipf (2019) Bin Dai and David Wipf.  Diagnosing and enhancing vae models.  _arXiv preprint arXiv:1903.05789_ , 2019.
  * Durmus & Moulines (2019) Alain Durmus and Eric Moulines.  High-dimensional bayesian inference via the unadjusted langevin algorithm.  _Bernoulli_ , 25(4A):2854–2882, 2019.
  * Earl & Deem (2005) David J Earl and Michael W Deem.  Parallel tempering: Theory, applications, and new perspectives.  _Physical Chemistry Chemical Physics_ , 7(23):3910–3916, 2005.
  * Freund & Haussler (1991) Yoav Freund and David Haussler.  Unsupervised learning of distributions on binary vectors using two layer networks.  _Advances in Neural Information Processing Systems_ , 4, 1991.
  * Gabrié et al. (2015) Marylou Gabrié, Eric W Tramel, and Florent Krzakala.  Training restricted boltzmann machine via the thouless-anderson-palmer free energy.  _Advances in Neural Information Processing Systems_ , 28, 2015.
  * Geman & Geman (1984) Stuart Geman and Donald Geman.  Stochastic relaxation, gibbs distributions, and the bayesian restoration of images.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , (6):721–741, 1984.
  * Grenander & Miller (1994) Ulf Grenander and Michael I Miller.  Representations of knowledge in complex systems.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 56(4):549–581, 1994.
  * Gu et al. (2022) Linyan Gu, Lihua Yang, and Feng Zhou.  Approximation properties of gaussian-binary restricted boltzmann machines and gaussian-binary deep belief networks.  _Neural Networks_ , 2022.
  * Heusel et al. (2017) Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Hinton (2002) Geoffrey E Hinton.  Training products of experts by minimizing contrastive divergence.  _Neural Computation_ , 14(8):1771–1800, 2002\.
  * Hinton & Salakhutdinov (2006) Geoffrey E Hinton and Ruslan R Salakhutdinov.  Reducing the dimensionality of data with neural networks.  _science_ , 313(5786):504–507, 2006.
  * Hinton & Van Camp (1993) Geoffrey E Hinton and Drew Van Camp.  Keeping the neural networks simple by minimizing the description length of the weights.  In _Proceedings of the sixth annual Conference on Computational Learning theory_ , pp. 5–13, 1993.
  * Hinton et al. (2006) Geoffrey E. Hinton, Simon Osindero, and Yee Whye Teh.  A fast learning algorithm for deep belief nets.  _Neural Computation_ , 18:1527–1554, 2006.
  * Jordan et al. (1999) Michael I Jordan, Zoubin Ghahramani, Tommi S Jaakkola, and Lawrence K Saul.  An introduction to variational methods for graphical models.  _Machine Learning_ , 37(2):183–233, 1999.
  * Krause et al. (2013) Oswin Krause, Asja Fischer, Tobias Glasmachers, and Christian Igel.  Approximation properties of dbns with binary hidden units and real-valued visible units.  In _International Conference on Machine Learning_ , pp. 419–426. PMLR, 2013.
  * Krizhevsky et al. (2009) Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * Lee et al. (2007) Honglak Lee, Chaitanya Ekanadham, and Andrew Ng.  Sparse deep belief net model for visual area v2.  _Advances in Neural Information Processing Systems_ , 20, 2007.
  * Liu et al. (2015) Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang.  Deep learning face attributes in the wild.  In _Proceedings of International Conference on Computer Vision (ICCV)_ , December 2015.
  * Melchior et al. (2017) Jan Melchior, Nan Wang, and Laurenz Wiskott.  Gaussian-binary restricted boltzmann machines for modeling natural image statistics.  _PloS one_ , 12(2):e0171015, 2017.
  * Nair & Hinton (2010) Vinod Nair and Geoffrey E Hinton.  Rectified linear units improve restricted boltzmann machines.  In _Icml_ , 2010.
  * Neal (1993) Radford M Neal.  _Probabilistic inference using Markov chain Monte Carlo methods_.  Department of Computer Science, University of Toronto Toronto, ON, Canada, 1993.
  * Ranzato & Hinton (2010) Marc’Aurelio Ranzato and Geoffrey E Hinton.  Modeling pixel means and covariances using factorized third-order boltzmann machines.  In _2010 IEEE Conference on Computer Vision and Pattern Recognition_ , pp. 2551–2558. IEEE, 2010.
  * Ranzato et al. (2010) Marc’Aurelio Ranzato, Alex Krizhevsky, and Geoffrey Hinton.  Factored 3-way restricted boltzmann machines for modeling natural images.  In _Proceedings of the thirteenth International Conference on Artificial Intelligence and Statistics_ , pp. 621–628. JMLR Workshop and Conference Proceedings, 2010.
  * Roberts & Tweedie (1996) Gareth O Roberts and Richard L Tweedie.  Exponential convergence of langevin distributions and their discrete approximations.  _Bernoulli_ , pp. 341–363, 1996.
  * (32) Tim Salimans, Andrej Karpathy, Xi Chen, and Diederik P Kingma.  Pixelcnn++: A pixelcnn implementation with discretized logistic mixture.  _ICLR_.
  * Smolensky (1986) P Smolensky.  Information processing in dynamical systems: foundations of harmony theory.  In _Parallel distributed processing: explorations in the microstructure of cognition, vol. 1: foundations_ , pp. 194–281. 1986.
  * Song & Ermon (2019) Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Takahashi & Yasuda (2016) Chako Takahashi and Muneki Yasuda.  Mean-field inference in gaussian restricted boltzmann machine.  _Journal of the Physical Society of Japan_ , 85(3):034001, 2016.
  * Theis et al. (2011) Lucas Theis, Sebastian Gerwinn, Fabian Sinz, and Matthias Bethge.  In all likelihood, deep belief is not enough.  _The Journal of Machine Learning Research_ , 12:3071–3096, 2011.
  * Tieleman (2008) Tijmen Tieleman.  Training restricted boltzmann machines using approximations to the likelihood gradient.  In _Proceedings of the 25th International Conference on Machine Learning_ , pp. 1064–1071, 2008.
  * Upadhya & Sastry (2021) Vidyadhar Upadhya and PS Sastry.  Learning gaussian-bernoulli rbms using difference of convex functions optimization.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2021\.
  * Vahdat & Kautz (2020) Arash Vahdat and Jan Kautz.  Nvae: A deep hierarchical variational autoencoder.  _Advances in Neural Information Processing Systems_ , 33:19667–19679, 2020.
  * Wang et al. (2012) Nan Wang, Jan Melchior, and Laurenz Wiskott.  An analysis of gaussian-binary restricted boltzmann machines for natural images.  In _ESANN_ , 2012.
  * Welling & Teh (2011) Max Welling and Yee W Teh.  Bayesian learning via stochastic gradient langevin dynamics.  In _Proceedings of the 28th International Conference on Machine Mearning (ICML-11)_ , pp. 681–688. Citeseer, 2011.
  * Welling et al. (2004) Max Welling, Michal Rosen-Zvi, and Geoffrey E Hinton.  Exponential family harmoniums with an application to information retrieval.  _Advances in Neural Information Processing Systems_ , 17, 2004.
  * Xiao et al. (2017) Han Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms, 2017.

##  Appendix A Derivations

###  A.1 Marginal Probability Distribution of Visible Units of GRBMs

We derive the marginal distribution of visible units as follows,

$\displaystyle p(\mathbf{v})$
$\displaystyle=\sum_{\mathbf{h}}p(\mathbf{v},\mathbf{h})$
$\displaystyle=\frac{1}{Z}\sum_{\mathbf{h}}\exp(-E_{\theta}(\mathbf{v},\mathbf{h}))$
$\displaystyle=\frac{1}{Z}\exp\left(-\frac{1}{2}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)^{\top}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)\right)\sum_{\mathbf{h}}\exp\left(\left(\frac{\mathbf{v}}{\bm{\sigma}^{2}}\right)^{\top}W\mathbf{h}+\mathbf{b}^{\top}\mathbf{h}\right)$
$\displaystyle=\frac{1}{Z}\exp\left(-\frac{1}{2}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)^{\top}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)\right)\prod_{i}\left(1+\exp\left(\left(\left(\frac{\mathbf{v}}{\bm{\sigma}^{2}}\right)^{\top}W\right)_{i}+\mathbf{b}_{i}\right)\right)$
$\displaystyle=\frac{1}{Z}\exp\left(-\frac{1}{2}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)^{\top}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)\right)\prod_{i}\exp\left(\text{Softplus}\left(\left(\left(\frac{\mathbf{v}}{\bm{\sigma}^{2}}\right)^{\top}W\right)_{i}+\mathbf{b}_{i}\right)\right)$
$\displaystyle=\frac{1}{Z}\exp\left(-\frac{1}{2}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)^{\top}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)\right)\exp\left(\text{Softplus}\left(W^{\top}\frac{\mathbf{v}}{\bm{\sigma}^{2}}+\mathbf{b}\right)^{\top}\textbf{1}\right)$
$\displaystyle=\frac{1}{Z}\exp\left(-\frac{1}{2}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)^{\top}\left(\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}}\right)+\text{Softplus}\left(W^{\top}\frac{\mathbf{v}}{\bm{\sigma}^{2}}+\mathbf{b}\right)^{\top}\textbf{1}\right).$
(20)

###  A.2 Gibbs Sampling

Algorithm 5 Gibbs Sampling for GRBMs

1:Input: number of steps $T$, burn-in step $\tilde{T}$

2:$\mathbf{v}^{(0)}\sim\mathcal{N}(\mathbf{0},I)$

3:For $t=1,\dots,T$

4: $\mathbf{h}^{(t)}\sim p(\mathbf{h}|\mathbf{v}^{(t-1)})$ $\triangleright$
following Eq. (4)

5: $\mathbf{v}^{(t)}\sim p(\mathbf{v}|\mathbf{h}^{(t)})$ $\triangleright$
following Eq. (3)

6:Return: $\\{(\mathbf{v}^{(t)},\mathbf{h}^{(t)})|t=\tilde{T}+1,\cdots,T\\}$

Gibbs sampling (Geman & Geman, 1984) is perhaps the most popular approach due
to its simplicity. In the context of GRBMs, we can alternate between sampling
hidden units given visible units and sampling visible units given hidden
units. Alg. 5 is a blocked Gibbs sampler; it samples all visible units (a
block of random variables) at once given all hidden units (the other block)
and vice versa. Given the conditional independence in the bipartite graphical
model, this block Gibbs sampler is equivalent to a univariate Gibbs sampler
that updates one variable at a time given the others following some schedule.
In fact, any schedule comprising a sequence of all hidden units followed by
all visible units or vice versa would make the equivalence hold. In other
words, it preserves the convergence of the original univariate Gibbs sampler
and runs as fast as a blocked Gibbs sampler. Relying on this Gibbs sampler, we
can get samples of visible and hidden units from the joint distribution in Eq.
(2). We can then discard the samples within the burn-in stage and treat
remaining ones as the final set of samples.

###  A.3 Langevin Sampling

The gradient of the marginal energy w.r.t. visible units is,

$\displaystyle\frac{\partial\tilde{E}(\mathbf{v})}{\partial\mathbf{v}}=\frac{\mathbf{v}-\bm{\mu}}{\bm{\sigma}^{2}}-\frac{W\,\text{Sigmoid}\left(W^{\top}\frac{\mathbf{v}}{\bm{\sigma}^{2}}+\mathbf{b}\right)}{\bm{\sigma}^{2}}.$
(21)

The cosine scheduler for annealing the step size is,

$\displaystyle\alpha_{k}$
$\displaystyle=\text{CosineScheduler}(k,K,\alpha_{0})=\frac{1}{2}\alpha_{0}\left(1+\cos\left(\frac{k}{K}\pi\right)\right)$
(22)

where $\alpha_{k}$ is the $k$-th step size, $\alpha_{0}$ is the initial step
size, and $K$ is the total number of steps.

The derivation of the Metropolis adjustment for Langevin sampling is as
follows,

$\displaystyle\tilde{A}(\mathbf{v}^{(t)},\mathbf{v}^{(t-1)})$
$\displaystyle=\min\left(1,\frac{p(\mathbf{v}^{(t)})q(\mathbf{v}^{(t-1)}|\mathbf{v}^{(t)})}{p(\mathbf{v}^{(t-1)})q(\mathbf{v}^{(t)}|\mathbf{v}^{(t-1)})}\right)$
$\displaystyle=\min\left(1,\frac{\exp\left(-\tilde{E}_{\theta}(\mathbf{v}^{(t)})\right)\exp\left(-\frac{1}{4\alpha_{t}}\left\|\mathbf{v}^{(t-1)}-\mathbf{v}^{(t)}+\alpha_{t}\frac{\partial\tilde{E}(\mathbf{v}^{(t)})}{\partial\mathbf{v}}\right\|^{2}\right)}{\exp\left(-\tilde{E}_{\theta}(\mathbf{v}^{(t-1)})\right)\exp\left(-\frac{1}{4\alpha_{t}}\left\|\mathbf{v}^{(t)}-\mathbf{v}^{(t-1)}+\alpha_{t}\frac{\partial\tilde{E}(\mathbf{v}^{(t-1)})}{\partial\mathbf{v}}\right\|^{2}\right)}\right)$
$\displaystyle=\min\left(1,\frac{\exp\left(-\tilde{E}_{\theta}(\mathbf{v}^{(t)})-\frac{1}{4\alpha_{t}}\left\|\mathbf{v}^{(t-1)}-\mathbf{v}^{(t)}+\alpha_{t}\frac{\partial\tilde{E}(\mathbf{v}^{(t)})}{\partial\mathbf{v}}\right\|^{2}\right)}{\exp\left(-\tilde{E}_{\theta}(\mathbf{v}^{(t-1)})-\frac{1}{4\alpha_{t}}\left\|\mathbf{v}^{(t)}-\mathbf{v}^{(t-1)}+\alpha_{t}\frac{\partial\tilde{E}(\mathbf{v}^{(t-1)})}{\partial\mathbf{v}}\right\|^{2}\right)}\right).$
(23)

###  A.4 Gibbs-Langevin Sampling

We now derive the Metropolis Adjustment for Gibbs-Langevin sampling. At time
step t-1, the proposal distribution in Alg. 2 is

$\displaystyle q(\mathbf{v},\mathbf{h}|\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})$
$\displaystyle=q(\mathbf{h}|\mathbf{v})q(\mathbf{v}|\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)}),$
(24)

where

$\displaystyle q(\mathbf{h}|\mathbf{v})$
$\displaystyle=\text{Sigmoid}\left(W^{\top}\left(\frac{\mathbf{v}}{\bm{\sigma}^{2}}\right)+\mathbf{b}\right).$
(25)

Denoting $\mathbf{v}^{(t-1)}=\tilde{\mathbf{v}}^{(0)}$ and
$\mathbf{v}=\tilde{\mathbf{v}}^{(K)}$, we have,

$\displaystyle q(\mathbf{v}|\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})$
$\displaystyle=q(\tilde{\mathbf{v}}^{(K)}|\tilde{\mathbf{v}}^{(0)},\mathbf{h}^{(t-1)})$
$\displaystyle=\int\cdots\int\left(\prod_{k=1}^{K}q(\tilde{\mathbf{v}}^{(k)}|\tilde{\mathbf{v}}^{(k-1)},\mathbf{h}^{(t-1)})\right)\mathrm{d}\tilde{\mathbf{v}}^{(1)}\cdots\mathrm{d}\tilde{\mathbf{v}}^{(K-1)},$
(26)

where

$\displaystyle q(\mathbf{v}|\tilde{\mathbf{v}}^{(k-1)},\mathbf{h}^{(t-1)})$
$\displaystyle=\mathcal{N}\left(\mathbf{v}\middle|\tilde{\mathbf{v}}^{(k-1)}-\alpha_{k}\frac{\partial
E(\tilde{\mathbf{v}}^{(k-1)},\mathbf{h}^{(t-1)})}{\partial\mathbf{v}},2\alpha_{k}I\right)$
(27) $\displaystyle\frac{\partial
E(\mathbf{v},\mathbf{h})}{\partial\mathbf{v}}$
$\displaystyle=\frac{\mathbf{v}-\bm{\mu}-W\mathbf{h}}{\bm{\sigma}^{2}}.$ (28)

The key question here is how to derive the analytical form of
$q(\tilde{\mathbf{v}}^{(K)}|\tilde{\mathbf{v}}^{(0)},\mathbf{h}^{(t-1)})$. The
most straightforward way is to compute the multiple integral directly. By
fixing all variables except for $\tilde{\mathbf{v}}^{(k)}$ in Eq. (8), we can
integrate out $\tilde{\mathbf{v}}^{(k)}$ analytically via the Gaussian
integral trick, _i.e_.,
$\int_{-\infty}^{\infty}\exp(-ax^{2}+bx+c)\mathrm{d}x=\sqrt{\frac{\pi}{a}}\exp(\frac{b^{2}}{4a}+c)$.
Then by applying the same trick recursively, one can ideally integrate out all
$\tilde{\mathbf{v}}^{(1)},\dots,\tilde{\mathbf{v}}^{(K-1)}$ in an analytical
manner. However, this process is quite involved due to the fact that the
integral of $\tilde{\mathbf{v}}^{(k)}$ depends on both
$\tilde{\mathbf{v}}^{(k+1)}$ and $\tilde{\mathbf{v}}^{(k-1)}$.

We instead resort to the reparameterization trick. In particular, at the outer
loop step $t$, the $k$-th inner loop step of Langevin sampling is as follows,

$\displaystyle\tilde{\mathbf{v}}^{(k)}$
$\displaystyle=\tilde{\mathbf{v}}^{(k-1)}-\alpha_{k}\frac{\partial
E(\tilde{\mathbf{v}}^{(k-1)},\mathbf{h}^{(t-1)})}{\partial\mathbf{v}}+\sqrt{2\alpha_{k}}\bm{\xi}_{k}$
$\displaystyle=\tilde{\mathbf{v}}^{(k-1)}-\alpha_{k}\frac{\tilde{\mathbf{v}}^{(k-1)}-\bm{\mu}-W\mathbf{h}^{(t-1)}}{\bm{\sigma}^{2}}+\sqrt{2\alpha_{k}}\bm{\xi}_{k}$
$\displaystyle=\left(\bm{1}-\frac{\alpha_{k}}{\bm{\sigma}^{2}}\right)\tilde{\mathbf{v}}^{(k-1)}+\alpha_{k}\frac{\bm{\mu}+W\mathbf{h}^{(t-1)}}{\bm{\sigma}^{2}}+\sqrt{2\alpha_{k}}\bm{\xi}_{k},$
(29)

where $\forall k\in\\{1,\dots,K\\}$, $\bm{\xi}_{k}\sim\mathcal{N}(0,I)$. This
discretization of Langevin dynamics gives a sample path of the distribution
$q(\tilde{\mathbf{v}}^{(K)}|\tilde{\mathbf{v}}^{(0)},\mathbf{h}^{(t-1)})$. We
now show that this sample path could be reparameterized as a simpler one which
gives the desirable analytical form of
$q(\tilde{\mathbf{v}}^{(K)}|\tilde{\mathbf{v}}^{(0)},\mathbf{h}^{(t-1)})$. To
simplify the derivation, we introduce
$\bm{\beta}_{k}=\prod_{j=k+1}^{K}\left(\bm{1}-\frac{\alpha_{j}}{\bm{\sigma}^{2}}\right)$,
$\forall k\in\\{0,\dots,K-1\\}$ and $\bm{\beta}_{K}=\bm{1}$. Therefore, after
$K$ steps, we have,

$\displaystyle\tilde{\mathbf{v}}^{(K)}$
$\displaystyle=\left(\bm{1}-\frac{\alpha_{K}}{\bm{\sigma}^{2}}\right)\tilde{\mathbf{v}}^{(K-1)}+\alpha_{K}\frac{\bm{\mu}+W\mathbf{h}^{(t-1)}}{\bm{\sigma}^{2}}+\sqrt{2\alpha_{K}}\bm{\xi}_{K}$
$\displaystyle=\left(\prod_{k=1}^{K}\left(\bm{1}-\frac{\alpha_{k}}{\bm{\sigma}^{2}}\right)\right)\tilde{\mathbf{v}}^{(0)}+\sum_{k=1}^{K}\left(\prod_{j=k+1}^{K}\left(\bm{1}-\frac{\alpha_{j}}{\bm{\sigma}^{2}}\right)\right)\left(\alpha_{k}\frac{\bm{\mu}+W\mathbf{h}^{(t-1)}}{\bm{\sigma}^{2}}+\sqrt{2\alpha_{k}}\bm{\xi}_{k}\right)$
$\displaystyle=\bm{\beta}_{0}\tilde{\mathbf{v}}^{(0)}+\left(\sum_{k=1}^{K}\bm{\beta}_{k}\alpha_{k}\right)\frac{\bm{\mu}+W\mathbf{h}^{(t-1)}}{\bm{\sigma}^{2}}+\sum_{k=1}^{K}\bm{\beta}_{k}\sqrt{2\alpha_{k}}\bm{\xi}_{k}$
(30)

Here $\\{\bm{\xi}_{k}|k=1,\dots,K\\}$ are independent random variables from
the standard Normal distribution. Since we know that the linear combination of
several independent Gaussian random variables leads to another Gaussian random
variable, we have

$\displaystyle
q(\tilde{\mathbf{v}}^{(K)}|\tilde{\mathbf{v}}^{(0)},\mathbf{h}^{(t-1)})=\mathcal{N}\left(\bm{\beta}_{0}\tilde{\mathbf{v}}^{(0)}+\left(\sum_{k=1}^{K}\bm{\beta}_{k}\alpha_{k}\right)\frac{\bm{\mu}+W\mathbf{h}^{(t-1)}}{\bm{\sigma}^{2}},\sum_{k=1}^{K}2\alpha_{k}\bm{\beta}_{k}^{2}\right).$
(31)

We can compute the acceptance probability,

$\displaystyle
A((\mathbf{v}^{(t)},\mathbf{h}^{(t)}),(\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)}))$
$\displaystyle=\min\left(1,\frac{p(\mathbf{v}^{(t)},\mathbf{h}^{(t)})q(\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)}|\mathbf{v}^{(t)},\mathbf{h}^{(t)})}{p(\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})q(\mathbf{v}^{(t)},\mathbf{h}^{(t)}|\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})}\right)$
$\displaystyle=\min\left(1,\frac{p(\mathbf{v}^{(t)},\mathbf{h}^{(t)})q(\mathbf{h}^{(t-1)}|\mathbf{v}^{(t-1)})q(\mathbf{v}^{(t-1)}|\mathbf{v}^{(t)},\mathbf{h}^{(t)})}{p(\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})q(\mathbf{h}^{(t)}|\mathbf{v}^{(t)})q(\mathbf{v}^{(t)}|\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})}\right)$
$\displaystyle=\min\left(1,\frac{\exp\left(-E_{\theta}(\mathbf{v}^{(t)},\mathbf{h}^{(t)})-\left\|\frac{\mathbf{v}^{(t-1)}-\bm{\beta}_{0}\mathbf{v}^{(t)}-\bm{a}(\bm{\mu}+W\mathbf{h}^{(t)})}{\sqrt{2}\tilde{\bm{\sigma}}}\right\|^{2}\right)q(\mathbf{h}^{(t-1)}|\mathbf{v}^{(t-1)})}{\exp\left(-E_{\theta}(\mathbf{v}^{(t-1)},\mathbf{h}^{(t-1)})-\left\|\frac{\mathbf{v}^{(t)}-\bm{\beta}_{0}\mathbf{v}^{(t-1)}-\bm{a}(\bm{\mu}+W\mathbf{h}^{(t-1)})}{\sqrt{2}\tilde{\bm{\sigma}}}\right\|^{2}\right)q(\mathbf{h}^{(t)}|\mathbf{v}^{(t)})}\right),$
(32)

where
$q(\mathbf{h}_{j}=1|\mathbf{v})=\left[\text{Sigmoid}\left(W^{\top}\frac{\mathbf{v}}{\bm{\sigma}^{2}}+\mathbf{b}\right)\right]_{j}$,
$\bm{a}=\frac{\sum_{k=1}^{K}\bm{\beta}_{k}\alpha_{k}}{\bm{\sigma}^{2}}$, and
$\tilde{\bm{\sigma}}^{2}=\sum_{k=1}^{K}2\alpha_{k}\bm{\beta}_{k}^{2}$.

###  A.5 Learning

We derive the detailed gradients of the marginalized log likelihood of visible
units w.r.t. model parameters as below.

$\displaystyle\nabla W_{ij}$
$\displaystyle=\left\langle\frac{\mathbf{v}_{i}}{\bm{\sigma}_{i}^{2}}\left[\text{Sigmoid}\left(W^{\top}\frac{\mathbf{v}}{\sigma^{2}}+\mathbf{b}\right)\right]_{j}\right\rangle_{d}-\left\langle\frac{\mathbf{v}_{i}}{\bm{\sigma}_{i}^{2}}\left[\text{Sigmoid}\left(W^{\top}\frac{\mathbf{v}}{\sigma^{2}}+\mathbf{b}\right)\right]_{j}\right\rangle_{m}$
(33) $\displaystyle\nabla\mu_{i}$
$\displaystyle=\left\langle\frac{\mathbf{v}_{i}-\mu_{i}}{\bm{\sigma}_{i}^{2}}\right\rangle_{d}-\left\langle\frac{\mathbf{v}_{i}-\mu_{i}}{\bm{\sigma}_{i}^{2}}\right\rangle_{m}$
(34) $\displaystyle\nabla\log\bm{\sigma}_{i}^{2}$
$\displaystyle=\left\langle\frac{(\mathbf{v}_{i}-\mu_{i})^{2}}{2\bm{\sigma}_{i}^{2}}-\frac{\sum_{j}\mathbf{v}_{i}W_{ij}\left[\text{Sigmoid}\left(W^{\top}\frac{\mathbf{v}}{\sigma^{2}}+\mathbf{b}\right)\right]_{j}}{\bm{\sigma}_{i}^{2}}\right\rangle_{d}$
$\displaystyle\quad-\left\langle\frac{(\mathbf{v}_{i}-\mu_{i})^{2}}{2\bm{\sigma}_{i}^{2}}-\frac{\sum_{j}\mathbf{v}_{i}W_{ij}\left[\text{Sigmoid}\left(W^{\top}\frac{\mathbf{v}}{\sigma^{2}}+\mathbf{b}\right)\right]_{j}}{\bm{\sigma}_{i}^{2}}\right\rangle_{m}$
(35) $\displaystyle\nabla\mathbf{b}_{i}$
$\displaystyle=\langle\left[\text{Sigmoid}\left(W^{\top}\frac{\mathbf{v}}{\sigma^{2}}+\mathbf{b}\right)\right]_{i}\rangle_{d}-\langle\left[\text{Sigmoid}\left(W^{\top}\frac{\mathbf{v}}{\sigma^{2}}+\mathbf{b}\right)\right]_{i}\rangle_{m}.$
(36)

##  Appendix B More Experimental Results

For all experiments, we set the initial variances of GRBMs to be $1$, clip the
gradient norm to be no larger than $10$, and use SGD with neither momentum nor
weight decay. We divide the total energy of a mini-batch by the batch size so
that we are minimizing the average negative log likelihood. We also decay the
learning rate of SGD from the initial value $0.01$ to $0$ using the same
cosine scheduler as described in Eq. (22). The burn-in step in CD learning is
set to $0$, _i.e_., we do not discard any samples from any Markov chains. For
all experiments involving images, we standardize the input image by
subtracting the pixel-wise mean and dividing by the pixel-wise standard
deviation. For color images, the subtraction and division is performed
channel-wise.

###  B.1 Gaussian Mixture Densities

The batch size and the hidden size are set to $100$ and $256$ respectively. We
adjust at every step whenever Metropolis adjustmentment is used as the
experiments with Gaussian mixture densities are fast. Although smaller hidden
size could work, but we found this size makes learning converges stably for
all sampling algorithms. To ensure a fair comparison, we train all GRBMs for
50K epochs and use the last model to draw the density plots and samples,
despite the learning processes with most of inference algorithms converge
within 5K to 10K epochs.

###  B.2 Ablation Study on MNIST

In this part, we perform ablation study on MNIST to investigate the effect of
several important factors. In all experiments on MNIST, we set batch size to
$512$ and the number of epochs to $3000$. First, we vary the CD step and the
hidden size while fixing the other hyperparameters. The results are shown in
Table 2. We found that $4096$ hidden size and $100$ CD steps work the best on
MNIST. More CD steps would potentially be better but take longer time to
train. Then we turn to study the number of Langevin sampling steps, the adjust
step size, the initial Langevin step size, and its annealing. Here annealing
means we decay the initial Langevin step size to $0$ following the cosine
scheduler as training goes on. The results are shown in Table 3. We can see
that the larger the initial Langevin step size, the better the performance.
But values larger than $0.04$ would sometimes make the sampling numerically
fail. The more the Langevin steps, the better the performance would be. Again,
it comes with more computational cost with more Langevin steps. We also find
that it may not be necessary to adjust at every step and annealing the initial
step size slightly improves the performance of Gbbis-Langevin with adjustment.

###  B.3 More Visual Results

We train 3K epochs for experiments on both FashionMNIST and CelebA-32
datasets. For CeleA-2K-64, we train 4K epochs. The batch size on FashionMNIST
and CelebA-32 is $512$ whereas the batch size on CeleA-2K-64 is $100$.

We show the samples drawn from the best GRBMs learned with different sampling
methods in Fig. 5. It is clear that samples corresponding to Gibbs-Langevin
have better visual qualities than those from Langevin and Gibbs. We also show
more results of GRBMs learned with Gibbs-Langevin in Fig. 6, Fig. 7, Fig. 8,
and Fig. 9.

Methods | CD Step | Hidden Size | FID
---|---|---|---
Gibbs-Langevin wo. Adjust | 50 | 2048 | 32.33
Gibbs-Langevin wo. Adjust | 50 | 4096 | 21.02
Gibbs-Langevin wo. Adjust | 50 | 8192 | 22.05
Gibbs-Langevin wo. Adjust | 100 | 4096 | 17.49
Table 2: Ablation study of the hidden size and the number of CD steps on MNIST dataset. Methods |  | Langevin
---
Step $K$
| Langevin
---
Step Size $\alpha_{0}$
Anneal $\alpha_{0}$ |  | Adjust
---
Step $\eta$
FID
Gibbs-Langevin wo. Adjust | 1 | 20 | ✗ | - | 35.08
Gibbs-Langevin wo. Adjust | 5 | 20 | ✗ | - | 19.00
Gibbs-Langevin wo. Adjust | 10 | 20 | ✗ | - | 17.49
Gibbs-Langevin wo. Adjust | 10 | 10 | ✗ | - | 21.05
Gibbs-Langevin wo. Adjust | 10 | 5 | ✗ | - | 25.67
Gibbs-Langevin w. Adjust | 10 | 20 | ✗ | 0 | 21.31
Gibbs-Langevin w. Adjust | 10 | 20 | ✗ | 25 | 21.25
Gibbs-Langevin w. Adjust | 10 | 20 | ✗ | 50 | 20.64
Gibbs-Langevin w. Adjust | 10 | 20 | ✓ | 50 | 19.27
Table 3: Ablation study of the number of Langevin steps $K$, the initial
Langevin step size $\alpha_{0}$, annealing of the initial Langevin step size,
and the Metropolis adjust step $\eta$ on MNIST dataset. All runs use 100 CD
steps.

(a) Gibbs

(b) Langevin

(c) Langevin
w. adjust

(d) Gibbs-Langevin

(e) Gibbs-Langevin w. adjust

Figure 5: Samples from GRBMs learned with different sampling algorithms on
MNIST.

Figure 6: More samples from the learned GRBM (Gibbs-Langevin) on MNIST.

Figure 7: More samples from the learned GRBM (Gibbs-Langevin) on FashionMNIST.

Figure 8: More samples from the learned GRBM (Gibbs-Langevin) on CelebA-32.

Figure 9: More samples from the learned GRBM (Gibbs-Langevin) on CelebA-2K-64.


# Minimal Model for Reservoir Computing

Yuzuru Sato  RIES / Department of Mathematics, Hokkaido University, Kita 12
Nishi 7, Kita-ku, Sapporo, Hokkaido 060-0812, Japan  London Mathematical
Laboratory, 8 Margravine Gardens, Hammersmith, London W6 8RH, UK  Miki U.
Kobayashi  Faculty of Economics, Rissho University, 4-2-16 Osaki, Shinagawa,
Tokyo 141-8602, Japan

###### Abstract

A minimal model for reservoir computing is studied. We demonstrate that a
reservoir computing exists that emulates given coupled maps by constructing a
modularised network. We describe a possible mechanism for collapses of the
emulation in the reservoir computing. Such transitory behaviour is caused by
either (i) an escape from a finite-time stagnation near an unstable chaotic
set, or (ii) a critical transition driven by the effective parameter drift.
Our approach reveals the essential mechanism for reservoir computing and
provides insights into the design of reservoir computing for practical
applications.

##  1 Introduction

Machine learning, including the reservoir computing, has proven useful for
model-free and data-driven predictions in a variety of scientific fields in
recent years (see [1, 2] and references therein). In studies on nonlinear
phenomena, the reservoir computing is typically applied to large-scale
phenomena such as chemical reaction [3], fluid turbulence [4, 5], and climate
and environmental dynamics [6, 7]. The purpose of the reservoir computing is
to emulate the dynamics $u(t+1)=g(u(t))$ behind time-series data $\\{u(t)\\}$
using the dynamics $\\{r(t)\\}$ that is generated by a large recurrent neural
network $r(t+1)=f_{A}(r(t))$ with a network matrix $A$ [8, 9, 10]. Assuming
that $f_{A}$ approximates any continuous functions by tuning $A$ when a
sufficiently large network is used [11], the trivial solution to this problem
is $r(0)=u(0),~{}f_{A}\simeq g$.

As opposed to solving the non-convex optimisation problem to approximates $g$
with $f_{A}$ using a statistical optimisation method such as gradient descent,
we may solve a simple convex optimisation problem using linear regression to
obtain a projection of the network dynamics $P\\{r(t)\\}$, which is close to
$\\{u(t)\\}$. The computation cost of the projection $P\\{r(t)\\}$ is clearly
substantially smaller than that of the function approximation $f_{A}\simeq g$;
however, we need to search for an appropriate initial network matrix $A$ for
the successful emulation. The matrix $A$ is typically provided as a random
matrix with a spectral radius $\rho(A)<1$ [8] (See Appendix A). One can search
a number of the quenched random dynamics $\\{r(t)\\}$ by changing the random
matrix $A$ to minimise $\lVert P\\{r(t)\\}-\\{u(t)\\}\rVert$. However, the
design principle of such random matrices that satisfy appropriate quality,
size, topology, and the order of the search cost is largely unknown.

As a statistical interpretation, it is pointed out that the reservoir
computing is some kind of nonlinear regression [12]. Indeed, a reservoir
computing performed by a neural network without hidden layers $f_{O}$ is a
simple perceptron, that is equivalent to logistic regression. It is also known
that random neural networks also approximate $C^{1}$ functions arbitrary well
[13], which corresponds to the classical universal approximation theorem [11],
and that the non-autonomous framework of reservoir computing is analogous to
Takens’s embedding theorem, that works even under presence of noise [14, 15,
13]. See also Appendix B for the simplest example of a classical delay
coordinate embedding performed by a non-autonomous reservoir computing.

In this study, instead of a Monte Carlo search for the network matrix $A$ and
a validation of the quality of the reservoir computing with randomly selected
$A$, we construct the smallest possible network that can be included in $A$ to
emulate the given dynamical systems. We analyse the behaviour of the
constructed dynamical systems, including transitory behaviour that induces the
collapse of the emulation [16, 17]. Our approach reveals the essential
mechanism for reservoir computing and provides insights into the design of
reservoir computing for practical applications.

##  2 Reservoir Computing

We consider a reservoir computing at time $t$ with $M$-dimensional input
neurons $\bm{x}(t)=(x^{(1)}(t),\ldots,x^{(M)}(t))^{T}$, output neurons
$\bm{y}(t)=(y^{(1)}(t),\ldots,y^{(M)}(t))^{T}$, and a reservoir $\cal{A}$. The
reservoir $\cal{A}$ includes $K$ hidden layer neurons
$\bm{r}(t)=(r^{(1)}(t),\ldots,r^{(K)}(t))^{T}$, a $K\times K$ random matrix
$A=(a_{ij})$ with a spectral radius $\rho(A)<1$, and an activation function
$\sigma(\cdot)=\tanh(\cdot)$. The weights for the input and output are
expressed as $W_{\rm{in}}=(w_{km})$ with $|w_{km}|<1$ and
$W_{\rm{out}}=(v_{mk})$, where $m=1,\ldots,M,~{}k=1,\ldots,K$, respectively.
The bias is given as $\bm{b}=(b_{k})$ with $(k=1,\ldots,K)$. We assume that
the training data are generated using an $M$-dimensional continuous map
$\bm{u}(t)=\bm{g}(\bm{u}(t))$.

In the non-autonomous training phase $t=-N,\ldots,0$, the entire system
dynamics is determinied as follows:

$\left\\{\begin{array}[]{l}\bm{x}(t)=\bm{u}(t)\\\
\bm{y}(t)=W_{\rm{out}}\bm{r}(t)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(t=-N,\ldots,0)\\\
r^{(k)}(t+1)=(1-\gamma)r^{(k)}(t)+\gamma\tanh\left[(A\bm{r}(t))_{k}+(W_{\rm{in}}\bm{x}(t))_{k}+b_{k}\right]\\\
\end{array}\right.,$ (1)

where $\gamma$ denotes the decay rate. At t=0, we minimise
$\lVert\\{\bm{y}(t)\\}-\\{\bm{u}(t)\\}\rVert$ based on the generated
$\\{\bm{r}(t)\\}$ to determine the optimal $W_{\rm{out}}=W_{\rm{out}}^{*}$.
With a sufficiently large training data, linear regression is used to
determine the optimal $W_{\rm{out}}^{*}$. At $t=0$, we replace $\bm{u}(t)$
with $\bm{y}(t)=W_{\rm{out}}^{*}\bm{r}(t)$, and express the system state as
$\bm{x}(0)=\bm{u}(0)=\bm{y}(0)$.

In the autonomous prediction phase $t=1,2,\ldots$, we obtain

$\left\\{\begin{array}[]{l}\bm{x}(t)=\bm{y}(t)=W_{\rm{out}}^{*}\bm{r}(t)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(t=1,2,\ldots)\\\
r^{(k)}(t+1)=(1-\gamma)r^{(k)}(t)+\gamma\tanh\left[(A\bm{r}(t))_{k}+(W_{\rm{in}}\bm{y}(t))_{k}+b_{k}\right]\end{array}\right.,$
(2)

and expect the emulation $\\{\bm{u}(t)\\}\simeq\\{\bm{y}(t)\\}$.

##  3 Emulating coupled maps

To construct a minimal example, we assume that the training data are generated
by a class of one-dimensional maps given as

$u(t+1)=(1-\gamma)u(t)+\gamma\sum_{k=1,2}v_{k}^{*}\tanh\left[\beta_{k}u(t)+\alpha_{k}\right].$
(3)

Eq. (3) shows chaotic behaviour in a broad range of parameters. We adopt
$\gamma=0.9,~{}~{}\beta_{1}=0.75,~{}~{}\beta_{2}=0.25,~{}~{}\alpha_{1}=\alpha_{2}=0,~{}~{}v^{*}_{1}=-12$,
and $v_{2}^{*}=16$ as an example of a generator of chaotic training data. We
emulate the dynamics of Eq. (3) using a reservoir $\cal{B}$ with $K=2$ (Fig. 1
(a)) to set

$B=\begin{pmatrix}\mu_{1}v_{1}^{*}&\mu_{1}v_{2}^{*}\\\
\mu_{2}v_{1}^{*}&\mu_{2}v_{2}^{*}\\\ \end{pmatrix},~{}~{}0<\mu_{1},\mu_{2}\ll
1,~{}\rho(B)<1,$ (4) $W_{\rm{in}}=\begin{bmatrix}\beta_{1}-\mu_{1}\\\
\beta_{2}-\mu_{2}\\\
\end{bmatrix},~{}~{}W_{\rm{out}}=\begin{bmatrix}v_{1}&v_{2}\\\
\end{bmatrix},~{}~{}\bm{b}=\begin{bmatrix}\alpha_{1}\\\ \alpha_{2}\\\
\end{bmatrix}.$ (5)

As $\mu_{1},\mu_{2}$ are free parameters, we may set $B$ satisfying the
condition $\rho(B)<1$. Following linear regression with a sufficiently large
training data $\\{u(t)\\}_{t=-N,\ldots,0}$, $W_{\rm{out}}$ is optimised to

$\displaystyle W_{\rm{out}}^{*}=[v_{1}^{*}~{}v_{2}^{*}]$ (6)

and we obtain

$\left\\{\begin{array}[]{l}x(t)=y(t)=v_{1}^{*}r^{(1)}(t)+v_{2}^{*}r^{(2)}(t)\\\
r^{(1)}(t+1)=(1-\gamma)r^{(1)}(t)+\gamma\tanh\left[\beta_{1}y(t)+\alpha_{1}\right]\\\
r^{(2)}(t+1)=(1-\gamma)r^{(2)}(t)+\gamma\tanh\left[\beta_{2}y(t)+\alpha_{2}\right]\\\
\end{array}\right.,$ (7)

resulting in

$y(t+1)=(1-\gamma)y(t)+\gamma\sum_{k=1,2}v_{k}^{*}\tanh\left[\beta_{k}y(t)+\alpha_{k}\right],$
(8)

which is equivalent to Eq. (3). Recall that at $t=0$ we obtain
$x(0)=u(0)=y(0)$, and the dynamics is on
$x=y=W_{\rm{out}}^{*}\bm{r}=v_{1}^{*}r^{(1)}+v_{2}^{*}r^{(2)}$, which may
include the generalised synchronisation manifold in the reservoir $\cal{B}$
for the input $\\{u(t)\\}$ [18, 19]. Thus, we achieve a perfect emulation
$\\{y(t)\\}=\\{u(t)\\}$ for $t=0,1,2,\ldots$ by using the $K=2$ reservoir
$\cal{B}$. We confirm that the above emulation is achieved in numerical
experiments with linear regression (Fig. 2 (a)). Many other possible
constructions are available for given one-dimensional maps.

Next, we emulate a class of coupled maps, as follows:

$\left\\{\begin{array}[]{l}u^{(1)}(t+1)=(1-\epsilon)f(u^{(1)}(t))+\epsilon
f(u^{(2)}(t))\\\ u^{(2)}(t+1)=(1-\epsilon)f(u^{(2)}(t))+\epsilon
f(u^{(1)}(t))\\\
f(u)=(1-\gamma)u+\gamma\sum_{k=1,2}v_{k}^{*}\tanh\left[\beta_{k}u+\alpha_{k}\right]\\\
\end{array}\right..$ (9)

To emulate the dynamics of Eq. (9) using a reservoir $\cal{C}$ with $K=4$
(Fig. 1 (b)), we set

$\begin{array}[]{l}C=\begin{pmatrix}\mu_{1}(1-\epsilon)v_{1}^{*}&\mu_{1}(1-\epsilon)v_{2}^{*}&\mu_{1}\epsilon
v_{1}^{*}&\mu_{1}\epsilon v_{2}^{*}\\\
\mu_{2}(1-\epsilon)v_{1}^{*}&\mu_{2}(1-\epsilon)v_{2}^{*}&\mu_{2}\epsilon
v_{1}^{*}&\mu_{2}\epsilon v_{2}^{*}\\\ \mu_{3}\epsilon
v_{1}^{*}&\mu_{3}\epsilon
v_{2}^{*}&\mu_{3}(1-\epsilon)v_{1}^{*}&\mu_{3}(1-\epsilon)v_{2}^{*}\\\
\mu_{4}\epsilon v_{1}^{*}&\mu_{4}\epsilon
v_{2}^{*}&\mu_{4}(1-\epsilon)v_{1}^{*}&\mu_{4}(1-\epsilon)v_{2}^{*}\\\
\end{pmatrix},\\\ \\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}0<\mu_{1},\mu_{2},\mu_{3},\mu_{4}\ll
1,~{}~{}\rho(C)<1\\\ \end{array},$ (10)
$W_{\rm{in}}=\begin{bmatrix}\beta_{1}-\mu_{1}&0\\\ \beta_{2}-\mu_{2}&0\\\
0&\beta_{1}-\mu_{3}\\\ 0&\beta_{2}-\mu_{4}\\\
\end{bmatrix},~{}~{}W_{\rm{out}}=\begin{bmatrix}v_{11}&v_{12}&v_{13}&v_{14}\\\
v_{21}&v_{22}&v_{23}&v_{24}\\\
\end{bmatrix},~{}~{}\bm{b}=\begin{bmatrix}\alpha_{1}\\\ \alpha_{2}\\\
\alpha_{1}\\\ \alpha_{2}\\\ \end{bmatrix}.$ (11)

Following linear regression, $W_{\rm{our}}$ is optimised to

$W_{\rm{out}}^{*}=\begin{bmatrix}(1-\epsilon)v_{1}^{*}&(1-\epsilon)v_{2}^{*}&\epsilon
v_{1}^{*}&\epsilon v_{2}^{*}\\\ \epsilon v_{1}^{*}&\epsilon
v_{2}^{*}&(1-\epsilon)v_{1}^{*}&(1-\epsilon)v_{2}^{*}\\\ \end{bmatrix},$ (12)

and we obtain

$\left\\{\begin{array}[]{l}x^{(1)}(t)=y^{(1)}(t)=(W_{\rm{out}}^{*}\bm{r}(t))_{1}\\\
x^{(2)}(t)=y^{(2)}(t)=(W_{\rm{out}}^{*}\bm{r}(t))_{2}\\\
r^{(1)}(t+1)=(1-\gamma)r^{(1)}(t)+\gamma\tanh\left[\beta_{1}y^{(1)}(t)+\alpha_{1}\right]\\\
r^{(2)}(t+1)=(1-\gamma)r^{(2)}(t)+\gamma\tanh\left[\beta_{2}y^{(1)}(t)+\alpha_{2}\right]\\\
r^{(3)}(t+1)=(1-\gamma)r^{(3)}(t)+\gamma\tanh\left[\beta_{1}y^{(2)}(t)+\alpha_{1}\right]\\\
r^{(4)}(t+1)=(1-\gamma)r^{(4)}(t)+\gamma\tanh\left[\beta_{2}y^{(2)}(t)+\alpha_{2}\right]\\\
\end{array}\right.,$ (13)

resulting in

$\left\\{\begin{array}[]{l}y^{(1)}(t+1)=(1-\epsilon)f(y^{(1)}(t))+\epsilon
f(y^{(2)}(t))\\\ y^{(2)}(t+1)=(1-\epsilon)f(y^{(2)}(t))+\epsilon
f(y^{(1)}(t))\\\
f(y)=(1-\gamma)y+\gamma\sum_{k=1,2}v_{k}^{*}\tanh\left[\beta_{k}y+\alpha_{k}\right]\\\
\end{array}\right.,$ (14)

which is equivalent to Eq. (9). Noting that $\bm{x}(0)=\bm{u}(0)=\bm{y}(0)$,
and the dynamics is on $\bm{x}=\bm{y}=W_{\rm{out}}^{*}\bm{r}$, we achieve a
perfect emulation $\\{\bm{y}(t)\\}=\\{\bm{u}(t)\\}$.

We confirm that the above emulation is achieved in numerical experiments with
linear regression (Fig. 2 (b)). Many other possible constructions are
available for given coupled maps.

Figure 1: (a) Model of $K=2$ reservoir for one-dimensional map, (b) model of
$K=4$ reservoir for $2$-coupled maps, and (c) model of emulation collapse
induced by excess neurons. The dotted circle indicates the collective excess
neurons. Figure 2: (a) Left: map provided in Eq. (3) and return plots of
$\\{y(t)\\}$ generated by constructed reservoir with $K=2$, and right:
attractor in $(r^{1)},r^{(2)},y)$, and the plane $y=W_{\rm{out}^{*}}\bm{r}$.
(b) Map provided in Eq. (9) and return plots of $\\{\bm{y}(t)\\}$ generated by
constructed reservoir with $K=4$. The parameters are
$\gamma=0.9,~{}~{}\beta_{1}=0.75,~{}~{}\beta_{2}=0.25,~{}~{}\alpha_{1}=\alpha_{2}=0,~{}~{}v^{*}_{1}=-12,~{}~{}v_{2}^{*}=16,\epsilon=0.05$,
and $\mu_{1}=\mu_{2}=\mu_{3}=\mu_{4}=0.01$.

In general, according to the property of neural networks as universal function
approximators, we may emulate any one-dimensional maps

$u(t+1)=(1-\gamma)u(t)+\gamma h(u(t)),$ (15)

where $h$ denotes an arbitrary continuous functions. For a large $K$, the
projected dynamics $\\{y(t)\\}$ of a reservoir with $K$ neurons can be
expressed as follows:

$\begin{array}[]{l}y(t+1)=(1-\gamma)y(t)+\gamma\sum_{k=1}^{K}v_{k}^{*}\tanh\left[\beta_{k}y(t)+\alpha_{k}\right]\\\
~{}~{}\stackrel{{\scriptstyle
K\to\infty}}{{\longrightarrow}}~{}~{}y(t+1)=(1-\gamma)y(t)+\gamma h(y(t)).\\\
\end{array}$ (16)

For example, the projection $\\{y(t)\\}$ may emulate the logistic map
$y(t+1)=a-y(t)^{2}$ when $h(y)\simeq-(1-\gamma)y+a-y^{2}$. Similarly, we can
emulate arbitrary $M$-coupled maps with a linear coupling. To do that we use
an $M$-dimensional projected dynamics of a reservoir with $KM$ neurons and a
network matrix that consists of $M$ blocks of the size $K$. For example, with
the construction in Eq. (9) where $f(y)\simeq a-y^{2}$, and the adequately
optimised $W_{\rm{out}}^{*}$, the $M$-dimensional projected dynamics
$\\{\bm{y}(t)\\}$ can emulate the dynamics of $M$-globally coupled logistic
map
$y^{(m)}(t+1)=(1-\epsilon)f(y^{(m)}(t))+\frac{\epsilon}{M}\sum_{i=1}^{M}f(y^{(i)}(t))$.
This construction for coupled maps partially explains the successful emulation
of PDEs by reservoir computing [16, 3]. Using a reservoir with a sufficiently
large network, it is certainly possible to emulate the dynamics of a wide
class of nonlinear PDEs, which are implemented on computers in the form of
coupled maps. However, the essential problem is rather determining why a
finite-size reservoir computing with a randomly selected network matrix
efficiently emulates such dynamical systems with large degrees of freedom.

##  4 Emulation in large networks

Although the matrix $B$ or $C$ may be a part of a larger network matrix $A$,
in certain cases in large networks, the system may degenerate and the
effective dynamics may depend only on a smaller part of the network
represented by $B$ or $C$. We consider such degeneration to the reservoir
$\cal{B}$.

In general, some of neurons may not contribute to the emulation. We refer to a
neuron with index $k$, where $w_{km}=v_{mk}=0$ holds for all $m$, as an
”excess neuron.” Assuming that $r^{(k)}~{}(k=3,\ldots,N)$ are excess neurons,
and that their dynamics vanishes as $r^{(k)}\to 0$ when without the bias, the
degeneration to the perfect emulation is achieved. We set a reservoir
$\cal{A}$ with $K$ neurons, as follows:

$A=\begin{pmatrix}B&R_{1}\\\ O&R_{2}\\\ \end{pmatrix},~{}~{}\rho(A)<1,$ (17)
$W_{\rm{in}}=\begin{bmatrix}\beta_{1}-\mu_{1}\\\ \beta_{2}-\mu_{2}\\\ 0\\\
\vdots\\\ 0\\\ \end{bmatrix},~{}~{}W_{\rm{out}}=\begin{bmatrix}v_{1}\\\
v_{2}\\\ v_{3}\\\ \vdots\\\ v_{K}\\\
\end{bmatrix}^{T},~{}~{}\bm{b}=\begin{bmatrix}\alpha_{1}\\\ \alpha_{2}\\\ 0\\\
\vdots\\\ 0\\\ \end{bmatrix}.$ (18)

When $R_{1}=O$ and $R_{2}$ is a random matrix, $W_{\rm{out}}$ is optimised to
$W_{\rm{out}}^{*}=\begin{bmatrix}v_{1}^{*}&v_{2}^{*}&0&\cdots&0\\\
\end{bmatrix}$, and Eq. (3) is emulated. In other cases, the projected
dynamics $\\{\bm{y}(t)\\}$ may deviate from the perfect emulation by the
perturbations from the collective excess neurons, which we discuss later in
the paper.

The system may also degenerate by synchronisation of neurons. For example, we
set a reservoir $\cal{D}$ as follows;

$\begin{array}[]{l}D=\begin{pmatrix}\mu_{1}\frac{v_{1}^{*}}{2}&\mu_{1}\frac{v_{2}^{*}}{2}&\mu_{1}\frac{v_{1}^{*}}{2}&\mu_{1}\frac{v_{2}^{*}}{2}\\\
\mu_{2}\frac{v_{1}^{*}}{2}&\mu_{2}\frac{v_{2}^{*}}{2}&\mu_{2}\frac{v_{1}^{*}}{2}&\mu_{2}\frac{v_{2}^{*}}{2}\\\
\mu_{3}\frac{v_{1}^{*}}{2}&\mu_{3}\frac{v_{2}^{*}}{2}&\mu_{3}\frac{v_{1}^{*}}{2}&\mu_{3}\frac{v_{2}^{*}}{2}\\\
\mu_{4}\frac{v_{1}^{*}}{2}&\mu_{4}\frac{v_{2}^{*}}{2}&\mu_{4}\frac{v_{1}^{*}}{2}&\mu_{4}\frac{v_{2}^{*}}{2}\\\
\end{pmatrix},~{}~{}0<\mu_{1},\mu_{2},\mu_{3},\mu_{4}\ll 1,~{}\rho(D)<1\\\
\end{array},$ (19) $W_{\rm{in}}=\begin{bmatrix}\beta_{1}-\mu_{1}\\\
\beta_{2}-\mu_{2}\\\ \beta_{1}-\mu_{3}\\\ \beta_{2}-\mu_{4}\\\
\end{bmatrix},~{}~{}W_{\rm{out}}^{*}=\begin{bmatrix}v_{1}^{*}&v_{2}^{*}&v_{3}^{*}&v_{4}^{*}\\\
\end{bmatrix},~{}~{}\bm{b}=\begin{bmatrix}\alpha_{1}\\\ \alpha_{2}\\\
\alpha_{1}\\\ \alpha_{2}\\\ \end{bmatrix}.$ (20)

and obtain

$\left\\{\begin{array}[]{l}\bm{x}(t)=\bm{y}(t)=W_{\rm{out}}^{*}\bm{r}(t)\\\
r^{(1)}(t+1)=(1-\gamma)r^{(1)}(t)+\gamma\tanh\left[\beta_{1}y(t)+\alpha_{1}\right]\\\
r^{(2)}(t+1)=(1-\gamma)r^{(2)}(t)+\gamma\tanh\left[\beta_{2}y(t)+\alpha_{2}\right]\\\
r^{(3)}(t+1)=(1-\gamma)r^{(3)}(t)+\gamma\tanh\left[\beta_{1}y(t)+\alpha_{1}\right]\\\
r^{(4)}(t+1)=(1-\gamma)r^{(4)}(t)+\gamma\tanh\left[\beta_{2}y(t)+\alpha_{2}\right]\\\
\end{array}\right..$ (21)

The dynamics of $(r^{(3)},r^{(4)})$ in Eq. (21) is eventually synchronised to
those of $(r^{(1)},r^{(2)})$ because

$\left\\{\begin{array}[]{l}r^{(1)}(t)-r^{(3)}(t)=(1-\gamma)(r^{(1)}(t)-r^{(3)}(t))\\\
r^{(2)}(t)-r^{(4)}(t)=(1-\gamma)(r^{(2)}(t)-r^{(4)}(t))\end{array}\right..$
(22)

As a result, the dynamics of $y(t)$ prefectly emulates Eq. (3). In large
matrices, when the $K$-dimensional vector $\bm{r}$ is degenerated to a
$2$-dimensional vector $(r^{(1)},r^{(2)})$ based on the synchronization to $2$
clusters $r^{(1)}=r^{(i)},r^{(2)}=r^{(j)},~{}(i\neq j,~{}3\leq i,j\leq K)$,
again, the perfect emulation of Eq. (3) is achieved. Many other constructions
are available for the degeneration by synchronisation.

While in gradient descent learning in multi-layer perceptrons, the
degeneration of the network weights may cause vanishing gradient and
decelerate learning [20, 21], the degeneration of the hidden layer dynamics in
the reservoir computing may enhance precise emulations.

##  5 Collapse of emulation

When excess neurons exist, the dynamics $\\{\bm{y}(t)\\}$ in the subspace
$\bm{x}=\bm{y}=W_{\rm{out}}^{*}\bm{r}$ can be perturbed, and the emulation may
collapse in the long run. We create a model for the collapse of the emulation
in Eq. (3) using a reservoir with $K=4$ and

$\begin{array}[]{l}E=\begin{pmatrix}\mu_{1}v_{1}^{*}&\mu_{1}v_{2}^{*}&0&0\\\
\mu_{2}v_{1}^{*}&\mu_{2}v_{2}^{*}&q_{1}&q_{2}\\\
0&0&\nu_{1}q_{1}&\nu_{1}q_{2}\\\ 0&0&\nu_{2}q_{1}&\nu_{2}q_{2}\\\
\end{pmatrix},~{}~{}0<\mu_{1},\mu_{2}\ll 1,~{}\rho(E)<1\\\ \end{array},$ (23)
$W_{\rm{in}}=\begin{bmatrix}\beta_{1}-\mu_{1}\\\ \beta_{2}-\mu_{2}\\\ 0\\\
0\\\
\end{bmatrix},~{}~{}W_{\rm{out}}^{*}=\begin{bmatrix}v_{1}^{*}&v_{2}^{*}&0&0\\\
\end{bmatrix},~{}~{}\bm{b}=\begin{bmatrix}\alpha_{1}\\\ \alpha_{2}\\\
\alpha_{3}\\\ \alpha_{4}\\\ \end{bmatrix}.$ (24)

For the entire system dynamics, we obtain

$\left\\{\begin{array}[]{l}x(t)=y(t)=v_{1}^{*}r^{(1)}(t)+v_{2}^{*}r^{(2)}(t)\\\
\eta(t+1)=q_{1}r^{(3)}(t)+q_{2}r^{(4)}(t)\\\
r^{(1)}(t+1)=(1-\gamma)r^{(1)}(t)+\gamma\tanh\left[\beta_{1}y(t)+\alpha_{1}\right]\\\
r^{(2)}(t+1)=(1-\gamma)r^{(2)}(t)+\gamma\tanh\left[\beta_{2}y(t)+\alpha_{2}+\eta(t)\right]\\\
r^{(3)}(t+1)=(1-\gamma)r^{(3)}(t)+\gamma\tanh\left[\nu_{1}\eta(t)+\alpha_{3}\right]\\\
r^{(4)}(t+1)=(1-\gamma)r^{(4)}(t)+\gamma\tanh\left[\nu_{2}\eta(t)+\alpha_{4}\right]\\\
\end{array}\right.,$ (25)
$\begin{array}[]{l}y(t+1)=(1-\gamma)y(t)+\gamma(v_{1}^{*}\tanh\left[\beta_{1}y(t)+\alpha_{1}\right]+v_{2}^{*}\tanh\left[\beta_{2}y(t)+\alpha_{2}+\eta(t))\right]\\\
\end{array},$ (26)
$\eta(t+1)=(1-\gamma)\eta(t)+\gamma\sum_{k=1,2}q_{k}\tanh\left[\nu_{k}\eta(t)+\alpha_{k}\right].$
(27)

When $\eta=0$, linear regression yields $W=W_{\rm{out}}^{*}$ and Eq. (26) is
equivalent to Eq. (3). Thus, the collective excess neurons provides external
perturbations to the projected dynamics $\\{y(t)\\}$. We present two examples
of collapse of the emulation.

In the first example, we set $\gamma=0.9,q_{1}=-12,q_{2}=16,\nu_{1}=0.75$, and
$\nu_{2}=0.5$ satisfying $\rho(E)<1$. When $\alpha_{3}=\tanh^{-1}(2/3)$ and
$\alpha_{4}=\tanh^{-1}(1/2)$, Eq. (27) has a tangency at the origin and a
stable fixed point $\eta=\eta^{*}\simeq 3.773341$. In this case, the perfect
emulation can be performed as a neutral solution on the centre manifold
$\eta=0$. The chaotic saddle is approached when we start at an initial point
in the effective basin defined by $\eta(0)<0$ and
$\eta(0)\notin(-2.197527,-0.985279)$; otherwise, the dynamics converges to
$\eta=\eta^{*}$ and the linear regression fails. If we are not at the onset of
the bifurcation point but very close to it with a narrow channel, the dynamics
of $\eta(t)$ remains near the origin for a long time. If the waiting time near
the origin, which typically follows a power law, is longer than the training
time $N$, we start at $\eta(0)\simeq\eta^{*}$ at time $t=0$ and the projected
dynamics $\\{y(t)\\}$ emulates $\\{u(t)\\}$ for a finite time. However, in the
long run, the dynamics of $\\{\eta(t)\\}$ escapes from the narrow channel near
the origin and converges to the stable fixed point $\eta=\eta^{*}$, and the
emulation collapses (FIG 3 (a)).

Figure 3: (a) Return plots and dynamics of the emulation collapse caused by a
tangency. The parameters are
$q_{1}=-12,q_{2}=16,\nu_{1}=0.75,\nu_{2}=0.5,\alpha_{3}=\tanh^{-1}(2/3)$, and
$\alpha_{4}=\tanh^{-1}(1/2)$. The black line indicates Eqs. (26) and (27). The
origin is neutrally stable and the fixed point $\eta=\eta^{*}\simeq 3.773341$
is stable in Eq. (27). (b) Return plots and dynamics of the emulation collapse
caused by a critical transition. The parameters are
$q_{1}=1,\nu_{1}=1-10^{-8}$, and $\alpha_{3}=0.001$. The black line indicates
Eqs. (26) and (28). A narrow channel near the origin and a stable fixed point
$\eta=\eta^{*}\simeq 0.143626$ exist.

The second example is based on only one excess neuron. We set
$\gamma=0.9,q_{1}=1,\nu_{1}=1-10^{-8}$, and $\alpha_{3}=0.001$ satisfying
$\rho(E)<1$ and obtain a hyperbolic tangent map

$\eta(t+1)=(1-\gamma)\eta(t)+\gamma\tanh\left[\nu_{1}\eta(t)+\alpha_{3}\right].$
(28)

In this case, the perfect emulation can be performed on an unstable chaotic
set. The dynamics of Eq. (28) does not have tangencies but may have a narrow
channel near the origin. It has a stable fixed point $\eta=\eta^{*}\simeq
0.143626$. The origin can be approached when we start at $\eta(0)<0$;
otherwise, the dynamics converges to the stable fixed point $\eta=\eta^{*}$
and the linear regression fails. A very slow uniform motion
$\eta(t)\simeq\gamma\tanh(\alpha_{3})t\simeq 10^{-4}t$ occurs near the origin,
which functions as a parameter drift in Eq. (26). Thus, considering that other
excess neurons function as external noise in a large network, bifurcation that
is similar to the critical transition may occur [22, 17]. We observe orbits
that are close to the training data for a finite time. However, before $\eta$
reaches to the stable fixed point $\eta=\eta^{*}$, it may arrive at the
tipping point of Eq. (26) and the emulation collapses (FIG 3 (b)). In these
two scenarios, the initial matrices and the initial conditions of $\bm{r}(0)$
have a positive measures in the parameter and state space. Thus, the above
phenomena can be observed in a quenched random dynamical systems with a random
matrix $A$ and with an initial condition $\bm{r}(0)$.

##  6 Conclusion

We have demonstrated that a reservoir computing exists that emulates given
coupled maps by constructing a modularised network. We have proposed a
possible mechanism for the collapses of the emulation in reservoir computing.
Such transitory behaviour is caused by either (i) a finite-time stagnation
near an unstable chaotic set or (ii) a critical transition by the effective
parameter drift. The essential problem in reservoir computing is determining
why a finite-size reservoir computing with a randomly selected network can
efficiently emulate dynamical systems with large degrees of freedom. Our
approaches provides a minimal model for understanding reservoir computing,
thereby providing better insights into the design of reservoir computing for
practical applications. Problems of quenched random dynamical systems analysis
for bifurcations, generalised synchronisation, and collapse of the emulation
will be studied elsewhere.

Acknowledgements We acknowledge H. Suetani for useful discussion and comments.
The research leading to these results has been partially funded by the Grant
in Aid for Scientific Research (B) No. 21H01002, JSPS, Japan, and the London
Mathematical Laboratory External Fellowship, United Kingdom.

##  Appendix A Echo state property

The echo state property is a characteristic of a given reservoir dynamics with
a network matrix $A$, where the reservoir dynamics converge to the same
reservoir behavior as time tends to infinity, regardless of the choice of
initial reservoir states [23]. As an empirical measure of the echo state
property, the condition of the spectral radius of the network matrix
$\rho(A)<1$ has been discussed. However, in our construction, we can set
arbitrary spectral radii for perfect emulations of any coupled maps. In
general, if multiple attractors exist, the echo state property does not hold
[24].

##  Appendix B Delay coordinate embedding by reservoir computing

Assuming that the training data is provided as one-dimensional time series
$\\{u(t)\\}$ and the attractor of the dynamics generating the training data
$\bm{y}(t+1)=\bm{g}(\bm{y}(t))$ is a compact manifold of dimension $d$, the
non-autonomous reservoir computing given by Eq. (1) with $2d+1$ output
neurons, and with a smooth activation function is embedding by the embedding
theorem [14, 15]. Both embedding and the reservoir computing unfold the
training data into a high-dimensional space (or a reservoir), and take a
projection of the high-dimensional structure. To see that with a simple
example, we construct a $K=3$ reservoir computing $\cal{Q}$ with
$1$-dimensional input neurons, $3$-dimensional output neurons, that embeds
$1$-dimensional time series to $3$-dimensional state space. We assume that the
training data are provided as $\\{u(t)\\}$ that is generated by a dynamical
system whose effective dimension of stationary attractor is less than $3$.
When the training data $\\{u(t)\\}$ is scaled to the dynamics very close to
the origin, the activation function $\sigma(x)=\tanh(x)$ is an almost
identical function as $\tanh(x)\simeq x~{}(x\simeq 0)$. We set $\gamma=1$, and

$Q=\begin{pmatrix}0&0&0\\\ 1&0&0\\\ 0&1&0\\\
\end{pmatrix},~{}~{}W_{\rm{in}}=\begin{bmatrix}1\\\ 0\\\ 0\\\
\end{bmatrix},~{}~{}W_{\rm{out}}^{*}=\begin{bmatrix}1&0&0\\\ 0&1&0\\\ 0&0&1\\\
\end{bmatrix},~{}~{}\bm{b}=\begin{bmatrix}0\\\ 0\\\ 0\\\ \end{bmatrix},$ (29)

and obtain

$\left\\{\begin{array}[]{l}x(t)=u(t)\\\
\bm{y}(t)=W_{\rm{out}}^{*}\bm{r}(t)=\bm{r}(t)\\\
r^{(1)}(t+1)=\tanh(x(t))\simeq x(t)\\\ r^{(2)}(t+1)=\tanh(r^{(1)}(t))\simeq
r^{(1)}(t)\\\ r^{(3)}(t+1)=\tanh(r^{(2)}(t))\simeq r^{(2)}(t)\\\
\end{array}\right..$ (30)

Thus, $\bm{y}(t)$ is the classical delay coordinate embedding of the training
data $\\{u(t)\\}$

$\left\\{\begin{array}[]{l}y^{(1)}(t+3)=r^{(1)}(t+3)\simeq x(t+2)=u(t+2)\\\
y^{(2)}(t+3)=r^{(2)}(t+3)\simeq x(t+1)=u(t+1)\\\
y^{(3)}(t+3)=r^{(3)}(t+3)\simeq x(t)=u(t)\\\ \end{array}\right.$ (31)

via a reservoir $\cal{Q}$. See [13] for the general embedding theorem for non-
autonomous echo state networks.

## References

  * Lukoševičius and Jaeger [2009] M. Lukoševičius and H. Jaeger, Computer science review 3, 127 (2009).
  * Tanaka et al. [2019] G. Tanaka, T. Yamane, J. B. Héroux, R. Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose, Neural Networks 115, 100 (2019).
  * Pathak et al. [2018] J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, Physical review letters 120, 024102 (2018).
  * Nakai and Saiki [2018] K. Nakai and Y. Saiki, Phys. Rev. E 98, 023111 (2018).
  * Kobayashi et al. [2021] M. U. Kobayashi, K. Nakai, Y. Saiki, and N. Tsutsumi, Phys. Rev. E 104, 044215 (2021).
  * Paulin [2010] C. Paulin, ournal of Hydrology 381, 76 (2010).
  * Mitsui and Boers [2021] T. Mitsui and N. Boers, Environmental Research Letters 16, 074024 (2021).
  * Jaeger [2001] H. Jaeger, GMD Report 148, 13 (2001).
  * Jaeger and Haas [2004] H. Jaeger and H. Haas, Science 304, 78 (2004).
  * Maass et al. [2002] W. Maass, T. Natschläger, and H. Markram, Neural Computation 14, 2531 (2002).
  * Cybenko [1989] G. Cybenko, Mathematics of control, signals and systems 2, 303 (1989).
  * Gauthier et al. [2021] D. J. Gauthier, E. Bollt, A. Griffith, and W. A. Barbosa, Nature communications 12, 5564 (2021).
  * Hart et al. [2020] A. Hart, J. Hook, and J. Dawes, Neural Networks 128, 234 (2020).
  * Takens [2006] F. Takens, in _Dynamical Systems and Turbulence, Warwick 1980: proceedings of a symposium held at the University of Warwick 1979/80_ (Springer, 2006), pp. 366–381.
  * Sauer et al. [1991] T. Sauer, J. A. Yorke, and M. Casdagli, Journal of statistical Physics 65, 579 (1991).
  * Pathak et al. [2017] J. Pathak, Z. Lu, B. R. Hunt, M. Girvan, and E. Ott, Chaos: An Interdisciplinary Journal of Nonlinear Science 27 (2017).
  * Kong et al. [2021] L.-W. Kong, H.-W. Fan, C. Grebogi, and Y.-C. Lai, Physical Review Research 3, 013090 (2021).
  * Lu et al. [2018] Z. Lu, B. R. Hunt, and E. Ott, Chaos: An Interdisciplinary Journal of Nonlinear Science 28 (2018).
  * Lymburn et al. [2019] T. Lymburn, D. M. Walker, M. Small, and T. Jüngling, Chaos: An Interdisciplinary Journal of Nonlinear Science 29 (2019).
  * Fukumizu and Amari [2000] K. Fukumizu and S.-I. Amari, Neural networks 13, 317 (2000).
  * Sato et al. [2022] Y. Sato, D. Tsutsui, and A. Fujiwara, Physica D: Nonlinear Phenomena 430, 133095 (2022).
  * Scheffer et al. [2009] M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter, V. Dakos, H. Held, E. H. Van Nes, M. Rietkerk, and G. Sugihara, Nature 461, 53 (2009).
  * Yildiz et al. [2012] I. B. Yildiz, H. Jaeger, and S. J. Kiebel, Neural networks 35, 1 (2012).
  * Ceni et al. [2020] A. Ceni, P. Ashwin, L. Livi, and C. Postlethwaite, Physica D: Nonlinear Phenomena 412, 132609 (2020).