# STANLEY: Stochastic Gradient Anisotropic
Langevin Dynamics for Learning Energy-Based Models

Belhal Karimi, Jianwen Xie, Ping Li

Cognitive Computing Lab
Baidu Research
10900 NE 8th St. Bellevue, WA 98004, USA

<EMAIL_ADDRESS>
###### Abstract

We propose in this paper, STANLEY, a STochastic gradient ANisotropic LangEvin
dYnamics, for sampling high dimensional data. With the growing efficacy and
potential of Energy-Based modeling, also known as non-normalized probabilistic
modeling, for modeling a generative process of different natures of high
dimensional data observations, we present an end-to-end learning algorithm for
Energy-Based models (EBM) with the purpose of improving the quality of the
resulting sampled data points. While the unknown normalizing constant of EBMs
makes the training procedure intractable, resorting to Markov Chain Monte
Carlo (MCMC) is in general a viable option. Realizing what MCMC entails for
the EBM training, we propose in this paper, a novel high dimensional sampling
method, based on an anisotropic stepsize and a gradient-informed covariance
matrix, embedded into a discretized Langevin diffusion. We motivate the
necessity for an anisotropic update of the negative samples in the Markov
Chain by the nonlinearity of the backbone of the EBM, here a Convolutional
Neural Network. Our resulting method, namely STANLEY, is an optimization
algorithm for training Energy-Based models via our newly introduced MCMC
method. We provide a theoretical understanding of our sampling scheme by
proving that the sampler leads to a geometrically uniformly ergodic Markov
Chain. Several image generation experiments are provided in our paper to show
the effectiveness of our method.

##  1 Introduction

The modeling of a data generating process is critical for many modern learning
tasks. A growing interest in generative models within the realm of computer
vision has led to multiple interesting solutions. In particular, Energy-based
models (EBM) (Zhu et al., 1998; LeCun et al., 2006), are a class of generative
models that learns high dimensional and complex (in terms of landscape)
representation/distribution of the input data. EBMs have been used in several
applications including computer vision (Ngiam et al., 2011; Xie et al., 2016;
Du and Mordatch, 2019), natural language processing (Mikolov et al., 2013;
Deng et al., 2020), density estimation (Li et al., 2019; Song et al., 2019)
and reinforcement learning (Haarnoja et al., 2017). Formally, EBMs are built
upon an unnormalized log probability, called the energy function, that is not
required to sum to one as standard log probability functions. This noticeable
feature allows for more freedom in the way one parameterizes the EBM. For
instance, Convolutional Neural Network can be employed to parameterize this
function, see Xie et al. (2016).

The training procedure of such models consists of finding an energy function
that assigns to lower energies to observations than unobserved points. This
phase can be cast as an optimization task and several ways are possible to
solve it. In this paper, we will focus on training the EBM via Maximum
Likelihood Estimation (MLE). Particularly, while using MLE to fit the EBM on
observed data, the high non-convexity of the loss function leads to a non
closed form maximization step. In general, gradient based optimization methods
are thus used during that phase. Besides, given the intractability of the
normalizing constant of our model, the aforementioned gradient, which is an
intractable integral, needs to be approximated. A popular and efficient way to
conduct such approximation is to use Monte Carlo approximation where the
samples are obtained via Markov Chain Monte Carlo (MCMC) (Meyn and Tweedie,
2012). The goal of this embedded MCMC procedure while training the EBM is to
synthesize new examples of the input data and use them to approximate
quantities of interest.

Hence, the sampling phase is crucial for both the EBM training speed and its
final accuracy in generating new synthetic samples. The computational burden
of those MCMC transitions, at each iteration of the EBM training procedure, is
alleviated via different techniques in the literature. For instance, in
Nijkamp et al. (2019), the authors develop a short-run MCMC as a flow-based
generator mechanism despite its non convergence property. Other principled
approach, as in Hinton (2002), keeps in memory the final chain state under the
previous global model parameter and uses it as the initialization of the
current chain. The heuristic of such approach is that along the EBM
iterations, the conditional distributions, depending on the model parameter,
are more and more similar and thus using a good sample from the previous chain
is in general a good sample of the current one. Though, this method can be
limited during the first iterations of the EBM training since when the model
parameter changes drastically, the conditional distributions do change too,
and samples from two different chains can be quite inconsistent. Several
extensions modifying the way the chain is initialized can be found in Welling
and Hinton (2002); Gao et al. (2018); Du and Mordatch (2019).

An interesting line of work in the realm of MCMC-based EBM tackles the biases
induced by stopping the MCMC runs too early. Indeed, it is known, see Meyn and
Tweedie (2012), that before convergence, MCMC samples are biased and thus
correcting this bias while keeping a short and less expensive run is an
appealing option. Several contributions aiming at removing this bias for
improved MCMC training include coupling MCMC chains, see Qiu et al. (2020);
Jacob et al. (2020) or simply estimating this bias and correct the chain
afterwards, see Du et al. (2021).

Here, our work is in line with the context of – high-dimensional data, – EBM
parameterized by deep neural networks and – MLE-based optimization via MCMC,
which make our method particularly attractive to all of the above combined. We
also consider the case of a short-run MCMC for the training of an EBM. Rather
than focusing on debiasing the chain, we develop a new sampling scheme where
the goal is to obtain better samples from the target distribution in fewer
MCMC transitions. We consider that the shape of the target distribution, which
inspires our proposed method, is of utmost importance to obtain such negative
samples. Our contributions are summarized below:

  1. 1.

We develop STANLEY, an Energy-Based model training method that embeds a newly
proposed _convergent_ and _efficient_ MCMC sampling scheme, focusing on
curvature informed metrics of the target distribution one wants to sample
from.

  2. 2.

Based on an anisotropic stepsize, our method, which is an improvement of the
Langevin Dynamics, achieves to obtain negative samples from the Energy-Based
model data distribution and improves the overall optimization algorithm.

  3. 3.

We prove the geometric ergodicity uniformly on any compact set of our MCMC
method assuming some regularity conditions on the target distribution and on
the backbone model of the EBM.

  4. 4.

We empirically assess the effectiveness of our method on several image
generation tasks, both on synthetic and real datasets including the Oxford
Flowers 102 dataset, CIFAR-10 and CelebA. We conclude the work with an Image
inpainting experiment on a benchmark dataset.

Roadmap. Section 2 introduces important notations and related work. Section 3
develops the main algorithmic contribution of this paper, namely STANLEY.
Section 4 presents our main theoretical results focusing on the ergodicity of
the proposed MCMC sampling method. Section 5 presents several image generation
experiments on both synthetic and real datasets. The complete proofs of our
theoretical results can be found in the supplementary material.

##  2 On MCMC-based Energy Based Models

Given a stream of input data noted $x\in\mathcal{X}\subset\mathbb{R}^{p}$, the
EBM is a Gibbs distribution defined as follows:

$p_{\theta}(x)=\frac{1}{Z(\theta)}\mathrm{exp}(f_{\theta}(x))\;,$ (1)

where $\theta\in\Theta\subset\mathbb{R}^{d}$ denotes the global parameters
vector of our model and
$Z(\theta)\mathrel{\mathop{:}}=\int_{x}\mathrm{exp}(f_{\theta}(x))\textrm{d}x$
is the normalizing constant (with respect to $x$). In particular,
$f_{\theta}(x):\mathcal{X}\to\mathbb{R}$, is the energy function (up to a
sign) that can be parameterized by a Convolutional Neural Network for
instance. The natural way of fitting model (1) is to employ Maximum Likelihood
Estimation (MLE) maximizing the marginal likelihood $p(\theta)$, i.e., finding
the vector $\theta^{*}$ such that for any $x\in\mathcal{X}$,

$\theta^{*}=\arg\max\limits_{\theta\in\Theta}\mathsf{L}(\theta)=\arg\max\limits_{\theta\in\Theta}\mathbb{E}_{q(x)}[\log
p_{\theta}(x)]\;,$ (2)

where $q(x)$ denotes the true distribution of the input data $x$. The
optimization task (2) is not tractable in closed form and requires an
iterative procedure in order to be solved. The standard algorithm used to
train EBMs is Stochastic Gradient Descent (SGD), see Robbins and Monro (1951);
Bottou et al. (2007). SGD requires having access to the gradient of the
objective function $\log p(\theta)$ which requires computing an intractable
integral, due to the high nonlinearity of the generally utilized parameterized
model $f_{\theta}(x)$. Given the general form defined in (1), we have that:

$\begin{split}\nabla\mathsf{L}(\theta)=\mathbb{E}_{q(x)}[\nabla_{\theta}f_{\theta}(x)]-\mathbb{E}_{p_{\theta}(x)}[\nabla_{\theta}f_{\theta}(x)]\;,\end{split}$

and a simple Monte Carlo approximation of $\nabla\log p(\theta)$ yields the
following expression of the gradient

$\nabla\mathsf{L}(\theta)\approx\frac{1}{n}\sum_{i=1}^{n}\nabla_{\theta}f_{\theta}(x^{q}_{i})-\frac{1}{M}\sum_{m=1}^{M}\nabla_{\theta}f_{\theta}(z_{m})\;,$
(3)

where $\\{z_{m}\\}_{m=1}^{M}$ are samples obtained from the EBM
$p_{\theta}(x)$ and $\\{x^{q}_{i}\\}_{i=1}^{n}$ are drawn uniformly from the
true data distribution $q(x)$. While drawing samples from the data
distribution is trivial, the challenge during the EBM training phase is to
obtain samples from the EBM distribution $p_{\theta}(x)$ for any model
parameter $\theta\in\Theta$. This task is generally performed using MCMC
methods. State-of-the-art MCMC used in the EBM literature include Langevin
Dynamics, see Grenander and Miller (1994); Roberts and Rosenthal (1998);
Roberts and Tweedie (1996) and Hamiltonian Monte Carlo (HMC), see Neal (2011).
Those methods are detailed in the sequel and are important concepts of our
contribution.

Energy Based Models. Energy based models are a class of generative models that
leverage the power of Gibbs potential and high dimensional sampling techniques
to produce high quality synthetic image samples. EBMs are powerful tools for
generative modeling tasks, as a building block for a wide variety of tasks.
The main purpose of EBMs is to learn an energy function (1) that assigns low
energy to a stream of observation and high energy values to other inputs.
Learning of such models is done via MLE (Xie et al., 2016; Du and Mordatch,
2019) or Score Matching (Hyvärinen, 2005) or Noise Constrastive Estimation
(Gao et al., 2020). In several general applications, authors leverage the
power of EBMs for develop an energy-based optimal policy where the parameters
of that energy function are provided by the reward of the overall system.
Learning EBMs with alternative strategies include contrastive divergence (CD)
(Hinton, 2002; Tieleman, 2008), noise contrastive estimation (NCE) (Gutmann
and Hyvärinen, 2010; Gao et al., 2020), introspective neural networks (INN)
(Lazarow et al., 2017; Jin et al., 2017; Lee et al., 2018), cooperative
networks (CoopNets) (Xie et al., 2018a, 2020, 2021c, 2022a, 2022c; Zhao et
al., 2023), f-divergence (Yu et al., 2020), and triangle divergence (Han et
al., 2019, 2020). Recently, EBMs parameterized by modern neural networks have
drawn much attention from the computer vision and machine learning
communities. Successful applications with EBMs include image generation (Xie
et al., 2016; Gao et al., 2018; Du and Mordatch, 2019; Zhao et al., 2021;
Zheng et al., 2021), videos (Xie et al., 2017, 2021d), 3D volumetric shapes
(Xie et al., 2018b, 2022b), unordered point clouds (Xie et al., 2021a), texts
(Deng et al., 2020), molecules (Ingraham et al., 2019; Du et al., 2020), as
well as image-to-image translation (Xie et al., 2022a, 2021b), unpaired cross-
domain image translation (Xie et al., 2021b; Song et al., 2023), out-of-
distribution detection (Liu et al., 2020), inverse optimal control (Xu et al.,
2022), deep regression (Gustafsson et al., 2020), salient object detection
(Zhang et al., 2022) and latent space modeling (Pang et al., 2020; Zhang et
al., 2021, 2023). Yet, unlike VAE (Kingma and Welling, 2014) or GAN
(Goodfellow et al., 2014) EBMs enjoy from a single structure requiring
training (versus several networks) resulting in more stability. The use of
implicit sampling techniques, such as MCMC, as detailed in the sequel, allows
more flexibility by trading off sample quality for computation time. Overall,
the _implicit_ property of the EBM, seen as an energy function, makes it a
tool of choice as opposed to _explicit_ generators that are limited to some
design choices, such as the choice of the prior distribution for VAEs or both
neural networks design in GANs.

MCMC procedures. Whether for sampling from a posterior distribution (Robert et
al., 2010; Han et al., 2017; Xie et al., 2019; Zhang et al., 2020; An et al.,
2021; Zhu et al., 2023; Xie et al., 2023), or in general intractable
likelihoods scenario (Doucet et al., 2000), various inference methods are
available. Approximate inference is a partial solution to the inference
problem and include techniques such as Variational Inference (VI) (Wainwright
and Jordan, 2008; de Freitas et al., 2001) or Laplace Approximation
(Wolfinger, 1993; Rue et al., 2009). Those methods allow the simplification of
the intractable quantities and result in the collection of good, yet
approximate, samples. As seen in (3), training an EBM requires obtaining
samples from the model itself. Given the nonconvexity of the structural model
$f_{\theta}(\cdot)$ with respect to the model parameter $\theta$, direct
sampling is not an option. Besides, in order to update the model parameter
$\theta$, usually through gradient descent type of methods (Bottou et al.,
2007), exact samples from the EBM are needed in order to compute a good
approximation of its (intractable) gradient, see (3). To do so, we generally
have recourse to MCMC methods. MCMC are a class of inference algorithms that
provide a principled iterative approach to obtain samples from any intractable
distribution. While being exact, the samples generally represent a larger
computation burden than methods such as VI. Increasing the efficiency of MCMC
methods, by obtaining exact samples, in other words constructing a chain that
converges faster, in fewer transitions is thus of utmost importance in the
context of optimizing EBMs. Several attempts have been proposed for the
standalone task of posterior sampling through the use of Langevin diffusion,
see the Unadjusted Langevin in Brosse et al. (2019), the MALA algorithm in
Roberts and Rosenthal (1998); Roberts and Tweedie (1996); Durmus et al. (2017)
or leveraging Hamiltonian Dynamics as in Girolami and Calderhead (2011). We
propose in the next section, an improvement of the Langevin diffusion with the
ultimate goal of speeding the EBM training procedure. Our method includes this
latter improvement in an end-to-end learning algorithms for Energy-Based
models.

##  3 Gradient Informed Langevin Diffusion

We now introduce the main algorithmic contribution of our paper, namely
STANLEY. STANLEY is a learning algorithm for EBMs, comprised of a novel MCMC
method for sampling samples from the intractable model (1). We provide
theoretical guarantees of our scheme in Section 4.

###  3.1 Preliminaries on Langevin MCMC based EBM

State-of-the-art MCMC sampling algorithm, particularly used during the
training procedure of EBMs, is the discretized Langevin diffusion, cast as
Stochastic Gradient Langevin Dynamics (SGLD), see Welling and Teh (2011). In
particular, several applications using EBM and SGLD have thrived in image
generation, natural language processing or even biology (Du et al., 2020).
Yet, the choice of the proposal, generally Gaussian, is critical for improving
the performances of both the sampling step (inner loop of the whole procedure)
and the EBM training. We recall the vanilla discretized Langevin diffusion
used in the related literature as follows:

$z_{k}=z_{k-1}+\frac{\gamma}{2}\nabla\log\pi_{\theta}(z_{k})+\sqrt{\gamma}B_{k}\;,$

where $\pi_{\theta}(\cdot):=p(\cdot,\theta)$ is the target potential one needs
samples from and defined in (1), $z_{k}$ represents the states of the chains
at iteration $k$, i.e., the generated samples in the context of EBM, $k$ is
the MCMC iteration index and $B_{k}$ is the Brownian motion, usually set as a
Gaussian noise and which can be written as
$B_{k}\mathrel{\mathop{:}}=\epsilon\,\xi_{k}$ where $\xi_{k}$ is a standard
Gaussian random variable and $\epsilon$ is a scaling factor for implementation
purposes. This method directs the proposed moves towards areas of high
probability of the stationary distribution $\pi_{\theta}$, for any
$\theta\in\Theta$ , using the gradient of $\log\pi_{\theta}$ and has been the
object of several studies (Girolami and Calderhead, 2011; Cotter et al.,
2013). In high dimensional and highly nonlinear settings, the burden of
computing this gradient for a certain number of MCMC transitions leads to a
natural focus: improving of the sampling scheme by assimilating information
about the landscape of the target distribution while keeping its ease of
implementation.

###  3.2 STANLEY, an Anisotropic Energy Based Modeling Approach

Given the drawbacks of current MCMC methods used for training EBMs, we
introduce a new sampler based on the Langevin updates presented above in Step
4 of Algorithm 1.

Algorithm 1 STANLEY for Energy-Based model

1: Input: Number of iterations $T$, number of MCMC steps $K$ and of samples
$M$, learning rates $\\{\eta_{t}\\}_{t>0}$, stepsize threshold th, initial
value $\theta_{0}$, MCMC initialization $\\{z_{0}^{m}\\}_{m=1}^{M}$ and data
$\\{x_{i}\\}_{i=1}^{n}$.

2: for $t=1$ to $T$ do

3: Compute the anisotropic stepsize as follows:
$\mathsf{\gamma}_{t}=\textsf{th}/\max(\textsf{th},|\nabla
f_{\theta_{t}}(z_{t-1}^{m})|)\;.$ (4)

4: Draw $M$ samples $\\{z_{t}^{m}\\}_{m=1}^{M}$ from the objective potential
(1) via Langevin diffusion:

5: for $k=1$ to $K$ do

6: Construct the Markov Chain with the Brownian motion $\mathsf{B}_{k}$ as
follows:  $z_{k}^{m}=z_{k-1}^{m}+\mathsf{\gamma}_{t}\nabla
f_{\theta_{t}}(z_{k-1}^{m})/2+\sqrt{\mathsf{\gamma}_{t}}\mathsf{B}_{k}\;,$ (5)

7: end for

8: Assign $\\{z_{t}^{m}\\}_{m=1}^{M}\leftarrow\\{z_{K}^{m}\\}_{m=1}^{M}$.

9: Sample $m$ positive observations $\\{x_{i}\\}_{i=1}^{m}$ from the empirical
data distribution.

10: Compute the gradient of the empirical log-EBM (1):
$\begin{split}\nabla\mathsf{L}(\theta_{t})=&\mathbb{E}_{p_{\text{data
}}}\left[\nabla_{\theta}f_{\theta_{t}}(x)\right]-\mathbb{E}_{p_{\theta}(x)}\left[\nabla_{\theta_{t}}f_{\theta}(z_{t})\right]\approx\frac{1}{n}\sum_{i=1}^{n}\nabla_{\theta}f_{\theta_{t}}\left(x_{i}\right)-\frac{1}{M}\sum_{m=1}^{M}\nabla_{\theta}f_{\theta_{t}}\left(z_{K}^{m}\right)\;.\end{split}$

11: Update the vector of global parameters of the EBM:
$\theta_{t+1}=\theta_{t}+\eta_{t}\nabla\mathsf{L}(\theta_{t})$

12: end for

13: Output: Vector of fitted parameters $\theta_{T+1}$.

Intuitions behind the efficacy of STANLEY: Some past modifications have been
proposed in particular to optimize the covariance matrix of the proposal of
the general MCMC procedure in order to better stride the support of the target
distribution. Langevin Dynamics is one example of those improvements where the
proposal is a Gaussian distribution where the mean depends on the gradient of
the log target distribution and the covariance depends on some Brownian
motion. For instance, in Atchadé (2006); Marshall and Roberts (2012), the
authors propose adaptive and geometrically ergodic Langevin chains. Yet, one
important characteristic of our EBM problem, is that for each model parameter
updated through the training iterations, the target distribution moves and the
proposal should take that adjustment into account. The techniques in Atchadé
(2006); Marshall and Roberts (2012) does not take the whole advantage of
changing the proposal using the target distribution. In particular, the
covariance matrix of the proposal is given by a stochastic approximation of
the empirical covariance matrix. This choice seems completely relevant as soon
as the convergence towards the stationary distribution is reached, in other
words it would make sense towards the end of the EBM training, as the target
distributions from a model parameter to the next one are similar. However, it
does not provide a good guess of the variability during the first iterations
since it is still very dependent on the initialization.

Moreover, in Girolami and Calderhead (2011), the authors consider the
approximation of a constant. Even though this simplification leads to ease of
implementation, the curvature metric chosen by the authors need to be
inverted, step that can be a computational burden if not intractable.
Especially in the case we are considering in our paper, i.e., ConvNet-based
EBM, where the high nonlinearity would lead to intractable expectations.
Therefore, in (4) and (5) of Algorithm 1, we propose a variant of Langevin
Dynamics, in order to sample from a target distribution, using a full
anisotropic covariance matrix based on the anisotropy and correlations of the
target distribution, see the $\sqrt{\mathsf{\gamma}_{t}}\mathsf{B}_{k}$ erm.

##  4 Geometric Ergodicity of STANLEY

We will present our theoretical analysis for the Markov Chain constructed
using Line 3-4 of Algorithm 1. Let $\Theta$ be a subset of $\mathbb{R}^{d}$
for some integer $d>0$. We denote by $\mathcal{Z}$ the measurable space of
$\mathbb{R}^{\ell}$ for some integer $\ell>0$. We define a family of
stationary distribution $\left(\pi_{\theta}(z)\right)_{\theta\in\Theta}$,
probability density functions with respect to the Lebesgue measure on the
measurable space $\mathcal{Z}$. This family of p.d.f. defines the stationary
distributions of our newly introduced sampler.

###  4.1 Notations and Assumptions

For any chain state $z\in\mathcal{Z}$, we denote by $\Pi_{\theta}(z,\cdot)$
the transition kernel as defined in the STANLEY update in Line 4 of Algorithm
1. The objective of this section is to rigorously show that each transition
kernel $\pi_{\theta}$, for any parameter $\theta\in\Theta$ is geometrically
ergodic and that this result holds for any compact subset
$\mathcal{C}\in\mathcal{Z}$. As a background note, a Markov chain, as
constructed Line 4, is said to be geometrically ergodic when $k$ iterations of
the same transition kernel is converging to the stationary distribution of the
chain with a geometric dependence on $k$.

As in Allassonniere and Kuhn (2015), we state the assumptions required for our
analysis. The first one is related to the continuity of the gradient of the
log posterior and the unit vectors pointing in the direction of the sample $z$
and in the direction of the gradient of the log posterior distribution at $z$:

######  H 1.

For all $\theta\in\Theta$, the structural model $f_{\theta}(\cdot)$ satisfies:

$\begin{split}\lim\limits_{|z|\to\infty}\frac{z}{|z|}\nabla
f_{\theta}(z)=-\infty\;,\;\;\lim\sup\limits_{|z|\to\infty}\frac{z}{|z|}\frac{\nabla
f_{\theta}(z)}{|\nabla f_{\theta}(z)|}<0\;.\end{split}$

Besides, we assume some regularity conditions of the stationary distributions
with respect to $\theta$:

######  H 2.

$\theta\to\pi_{\theta}$ and $\theta\to\nabla\log\pi_{\theta}$ are continuous
on $\Theta$.

For a positive and finite function noted $V:\mathcal{Z}\mapsto\mathbb{R}$, we
define the V-norm distance between two arbitrary transition kernels $\Pi_{1}$
and $\Pi_{2}$ as follows:

$\|\Pi_{1}-\Pi_{2}\|_{V}\mathrel{\mathop{:}}=\sup\limits_{z\in\mathcal{Z}}\frac{\|\Pi_{1}(z,\cdot)-\Pi_{2}(z,\cdot)\|_{V}}{V(z)}\;.$

The definition of this norm allows us to establish a convergence rate for our
sampling method by deriving an upper bound of
$\|\Pi_{\theta}^{k}-\pi_{\theta}\|_{V}$ where $k>0$ denotes the number of MCMC
transitions. We recall that $\Pi_{\theta}$ is the transition kernel defined by
Line 4 of Algorithm 1 and $\pi_{\theta}$ is the stationary distribution of our
Markov chain at a given EBM model $\theta$. This quantity characterizes how
close to the target distribution our chain is getting, after a finite time of
iterations and will eventually formalize the _V-uniform ergodicity_ of our
method. We specify that strictly speaking, $\pi_{\theta}$ is a probability
measure, and not a transition kernel. However
$\|\Pi_{\theta}^{k}-\pi_{\theta}\|_{V}$ is well-defined if we consider
$\pi_{\theta}$ as a kernel:

$\pi(z,\mathcal{C})\mathrel{\mathop{:}}=\pi(\mathcal{C})\quad\textrm{for}\quad\mathcal{C}\subset\mathcal{Z},\quad
z\in\mathcal{Z}\;.$

Here, for some $\beta\in]0,1[$ we define the $V_{\theta}$ function, also know
as the _drift_ , for all $z\in\mathcal{Z}$ as follows:

$V_{\theta}(z)\mathrel{\mathop{:}}=c_{\theta}\pi_{\theta}(z)^{-\beta}\;,$ (6)

where $c_{\theta}$ is a constant, with respect to the chain state $z$, such
that for all $z\in\mathcal{Z}$, $V_{\theta}(z)\geq 1$. Note that the V norm
depends on the chain state noted $z$ _and_ of the global model parameter
$\theta$ varying through the optimization procedure. Yet, in both main
results, the ergodicity and the convergence rate, including the underlying
drift condition, are established uniformly on the parameter space $\Theta$. We
also define the auxiliary functions, independent of the parameter $\theta$ as:

$V_{1}(z)\mathrel{\mathop{:}}=\inf\limits_{\theta\in\Theta}V_{\theta}(z)\quad\textrm{and}\quad
V_{2}(z)\mathrel{\mathop{:}}=\sup\limits_{\theta\in\Theta}V_{\theta}(z)\;,$
(7)

and assume the following:

######  H 3.

There exists a constant $a_{0}>0$ such that for all $\theta\in\Theta$ and
$z\in\mathcal{Z}$, the function $V_{2}^{a_{0}}(z)$, defined in (7), is
integrable against the kernel $\Pi_{\theta}(z,\cdot)$ and we have

$\lim\sup\limits_{a\to
0}\sup\limits_{\theta\in\Theta,z\in\mathcal{Z}}\Pi_{\theta}V_{2}^{a}(z)=1\;.$

###  4.2 Convergence Results

The result consists in showing V-uniform ergodicity of the chain, the
irreducibility of the transition kernels and their aperiodicity, following
Meyn and Tweedie (2012); Allassonniere and Kuhn (2015). We also prove a drift
condition which states that the transition kernels tend to bring back elements
into a small set. Then, V-uniform ergodicity of the transition kernels
$(\Pi_{\theta})_{\theta\in\Theta}$ boils down from the latter proven drift
condition.

Important Note: The stationary distributions depend on $\theta\in\Theta$ as
they vary at each model update during the EBM optimization phase. Thus uniform
convergence of the chain is important in order to characterize the sampling
phase _throughout the entire training phase_. Particularly at the beginning,
the shape of the distributions one needs to sample from varies a lot from a
parameter to another.

Theorem 4.2 shows two important convergence results for our sampling method.
First, it establishes the existence of a small set $\mathcal{O}$ leading to
the crucially needed aperiodicity of the chain and ensuring that each
transition moves towards a better state. Then, it provides a uniform
ergodicity result of our sampling method in STANLEY, via the so-called _drift
condition_ providing the guarantee that our transition kernels
$(\Pi_{\theta})_{\theta\in\Theta}$ attract the states into the small set
$\mathcal{O}$. Moreover, the independence on the EBM model parameter $\theta$
of $V$ in (9) leads to _uniform_ ergodicity as shown in the Corollary 4.2.

Assume H1-H3. For any $\theta\in\Theta$, there exists a drift function
$V_{\theta}$, a set $\mathcal{O}\subset\mathcal{Z}$, a constant
$0<\epsilon\leq 1$ such that

$\Pi_{\theta}(z,\mathcal{B})\geq\epsilon\int_{\mathcal{B}}\mathsf{1}_{\mathcal{X}}(z)\textrm{d}y\;.$
(8)

Moreover there exists $0<\mu<1$, $\delta>0$ and a drift function $V$,
independent of $\theta$ such that for all $z\in\mathcal{Z}$:

$\Pi_{\theta}V(z)\leq\mu V(z)+\delta\mathsf{1}_{\mathcal{O}}(z)\;.$ (9)

Assume H1-H3. A direct consequence of Theorem 4.2 is that the family of
transition kernels $(\Pi_{\theta})_{\theta\in\Theta}$ are uniformly
ergodic,i.e., for any compact $\mathcal{C}\subset\mathcal{Z}$, there exist
constants $\rho\in]0,1[$ and $e>0$ such for any MCMC iteration $k>0$,we have:

$\sup\limits_{z\in\mathcal{C}}\|\Pi_{\theta}^{k}u(\cdot)-\pi_{\theta}u(\cdot)\|_{V}\leq
e\rho^{k}\|u\|_{V}\;,$ (10)

where $V$ is the drift function in Theorem 4.2 and $u(\cdot)$ is any bounded
function we apply a transition to. While Theorem 4.2 is critical for proving
the aperiodicity and irreducibility of the chain, we establish the geometric
convergence speed of the chain. We do not only show the importance of the
_uniform_ ergodicity of the chain, which makes it appealing for the EBM
training since the model parameter $\theta$ is often updated, but we also
derive a geometrical rate in Corollary 4.2.

We encourage the readers to read through the sketch of the main Theorem of our
paper provided on the first page of the supplemental as we give the important
details leading to the desired ergodicity results. Those various techniques
are common in the MCMC literature and we refer the readers to several MCMC
handbooks such as Neal (2011); Meyn and Tweedie (2012) for more understanding.

##  5 Numerical Experiments

We conduct a collection of experiments to show the effectiveness of our
method, both on synthetic and real datasets. After verifying the advantage of
STANLEY on a Gaussian Mixture Model (GMM) retrieving the synthetic data
observations, we then investigate its performance when learning a distribution
over high-dimensional natural images such as pictures of flowers, see the
Flowers dataset in Nilsback and Zisserman (2008), or general concepts featured
in CIFAR-10 (Krizhevsky and Hinton, 2009). For both methods, we use the
Frechet Inception Distance (FID), as a reliable performance metrics as
detailed in Heusel et al. (2017). In the sequel, we tune the learning rates
over a fine grid and report the best result for all methods. For our method
STANLEY, the threshold parameter th, crucial for the implementation of the
stepsize (4) is tuned over a grid search as well. As mentioned above, we also
define a Brownian motion as $B_{k}\mathrel{\mathop{:}}=\epsilon\,\xi_{k}$, and
tune the scaling factor $\epsilon$ for better performances.

###  5.1 Toy Example: Gaussian Mixture Model

Datasets. We first demonstrate the outcomes of both methods including our
newly proposed STANLEY for low-dimensional toy distributions. We generate
synthetic 2D rings data and use an EBM to learn the true data distribution and
put it to the test of generating new synthetic samples.

Methods and Settings. We consider two methods. Methods are ran with
_nonconvergent_ MCMC, i.e.,, we do not necessitate the convergence to the
stationary distribution of the Markov chains. The number of transitions of the
MCMC is set to $K=100$ per EBM iteration. We use a standard deviation of
$0.15$ as in Nijkamp et al. (2020). Both methods have a constant learning rate
of $0.14$. The value of the threshold th for our STANLEY method is set to
$\textsf{th}=0.01$. The total number of EBM iterations is set to $T=10\,000$.
The global learning rate $\eta$ is set to a constant equal to $0.0001$.

Network architectures. For the backbone of the EBM model, noted
$f_{\theta}(\cdot)$ in (1), we chose a CNN of $5$ 2D convolutional layers and
Leaky ReLU activation functions, with the leakage parameter set to $0.05$. The
number of hidden neurons varies between $32$ and $64$.

Figure 1: (Rings Toy Dataset) Top: our method, namely STANLEY Bottom: vanilla
Langevin Dynamics. Methods are used with the same backbone architecture.
Generated samples are plotted through the iterations ever $2\,000$ steps.

Results. We observe Figure 1 the outputs using both methods on the toy
dataset. While they achieve a great representation of the truth after a large
number of iterations, we notice that STANLEY learns an energy that closely
approximates the true density during the first thousands of iterations if the
training process. The sharpness of the data generated by STANLEY in the first
iterations shows an empirically better ability to sample from the 2D dataset.

###  5.2 Image Generation

Datasets. We run our method and several baselines detailed below on the
CIFAR-10 dataset (Krizhevsky and Hinton, 2009) and the _Oxford Flowers 102_
dataset (Nilsback and Zisserman, 2008). CIFAR-10 is a popular computer-vision
dataset of $50\,000$ training images and $10\,000$ test images, of size
$32\times 32$. It is composed of tiny natural images representing a wide
variety of objects and scenes, making the task of self supervision supposedly
harder. The _Oxford Flowers 102_ dataset is composed of 102 flower categories.
Per request of the authors, the images have large scale, pose and light
variations making the task of generating new samples particularly challenging.

Figure 2: Left: Langevin Method. Right: STANLEY . After 100k iterations.

Methods and Settings for the Flowers dataset. Nonconvergent MCMC are also used
in this experiment and the number of MCMC transitions is set to $K=50$. Global
learning parameters of the gradient descent update is set to $0.001$ for both
methods. We run each method during $T=100\,000$ iterations and plot the
results using the final vector of fitted parameters.

Methods and Settings for CIFAR-10. We employ the same nonconvergent MCMC
strategies for this experiment. The value of the threshold th for our STANLEY
method is set to $\textsf{th}=0.0002$. The total number of EBM iterations is
set to $T=100\,000$. The global learning rate $\eta$ is set to a constant
equal to $0.0001$. In this experiment, we slightly change the last step of our
method in Algorithm 1. Indeed, Line 11 in Algorithm 1 is not a plain
Stochastic Gradient Descent here but we rather use the Adam optimizer (Kingma
and Ba, 2015). The scaling factor of the Brownian motion is $0.01$.

Network architectures for both. The backbone of the energy function for this
experiment is a vanilla ConvNet composed of $3\times 3$ convolution layers
with stride $1$. $5$ Convolutional Layers using ReLU activation functions are
stacked.

Results. (Flowers) Visual results are provided in Figure 2 where we have used
both methods to generate synthetic images of flowers. For each threshold
iterations number ($5\,000$ iterations) we sample $10\,000$ synthetic images
from the EBM model under the current vector of parameters and use the same
number of data observations to compute the FID similarity score as advocated
in Heusel et al. (2017). The evolution of the FID values are reported in
Figure 3 (Left) through the iterations. We note that our method outperforms
the other baselines for all iterations threshold, including the Vanilla
Langevin (in blue) which is an ablated form our STANLEY (no adaptive
stepsize).

Figure 3: (FID values per method against 100k iterations elapsed). Left:
Oxford Flowers dataset. Right: CIFAR-10.

(CIFAR-10) Visual results are provided in Figure 4 where we have used both
methods to generate synthetic images of flowers. The FID values are reported
in Figure 3 (Right) and have been computed using $10\,000$ synthetic images
from each model. The similarity score is then evaluated every $5\,000$
iterations. While the FID curves for the Flowers dataset exhibits a superior
performance of our method throughout the training procedure, we notice that in
the case of CIFAR-10, vanilla method seems to be slightly better than
STANLEYḋuring the first iterations, i.e., when the model is still learning the
representation of the images. Yet, after a certain number of iterations, we
observe that STANLEY leads to more accurate synthetic images. _This behavior
can be explained by the importance of incorporating curvature informed metrics
into the training process when the parameter reaches a neighborhood of the
solution_.

Figure 4: 1: Langevin 2: STANLEY 3: MH 4: HMC 5: GD without noise. After 100k
iterations.

###  5.3 Image Inpainting

The image inpainting experiment aims to fill missing regions of a damaged
image with synthesized content.

Figure 5: FID values per method against 50k iterations elapsed.

Datasets. We use the CelebA dataset (Liu et al., 2015) to evaluate our
learning algorithm, which contains more than 200k RGB color facial image. We
use 100k images for training and 100 images for testing.

Methods and Settings. Nonconvergent MCMC are also used in this experiment and
the number of MCMC transitions is set to $K=50$. Global learning parameters of
the gradient descent update is set to $0.01$. We run each method during
$T=50\,000$ iterations and plot the results using the final vector of fitted
parameters.

Results. Figure 5 displays the FID curves for all methods. We note that along
the iterations, the STANLEY outperforms the other baseline and is similar to
HMC, while only requiring first order information for the computation of the
stepsize whereas HMC computes second order quantity. Even with second order
information, the HMC samples does not lead to a better FID.

(a)

(b)

Figure 6: Image inpainting. Top: STANLEY Bottom: Vanilla Langevin.

Figure 6 shows the visual check on different samples between our method and
its ablated form, i.e., the vanilla Langevin sampler based EBM. In addition to
a measurable metrics comparison, visual check provide empirical insights on
the effectiveness of adding a curvature-informed stepsize in the sampler of
our generative model as in STANLEY.

##  6 Conclusion

We propose in this paper, an improvement of the so-called MCMC based Energy-
Based models. In the particular case of a highly nonlinear structural model of
the EBM, more precisely a Convolutional Neural Network in our paper, we tackle
the complex task of sampling negative samples from the energy function. The
multi-modal and highly curved landscape one must sample from inspire our
technique called STANLEY, and based on a Stochastic Gradient Anisotropic
Langevin Dynamics, that updates the Markov Chain using an anisotropic stepsize
in the vanilla Langevin update. We provide strong theoretical guarantees for
our novel method, including uniform ergodicity and geometric convergence rate
of the transition kernels to the stationary distribution of the chain. Our
method is tested on several benchmarks data and image generation tasks
including toy and real datasets.

## References

  * Allassonniere and Kuhn (2015) Stéphanie Allassonniere and Estelle Kuhn.  Convergent stochastic expectation maximization algorithm with efficient sampling in high dimension. application to deformable template model estimation.  _Computational Statistics & Data Analysis_, 91:4–19, 2015\.
  * An et al. (2021) Dongsheng An, Jianwen Xie, and Ping Li.  Learning deep latent variable models by short-run MCMC inference with optimal transport correction.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 15415–15424, virtual, 2021.
  * Atchadé (2006) Yves F Atchadé.  An adaptive version for the metropolis adjusted langevin algorithm with a truncated drift.  _Methodology and Computing in applied Probability_ , 8(2):235–254, 2006.
  * Bottou et al. (2007) Léon Bottou, Olivier Chapelle, Dennis DeCoste, and Jason Weston, editors.  _Large-Scale Kernel Machines_.  The MIT Press, Cambridge, MA, 2007.
  * Brosse et al. (2019) Nicolas Brosse, Alain Durmus, Éric Moulines, and Sotirios Sabanis.  The tamed unadjusted langevin algorithm.  _Stochastic Processes and their Applications_ , 129(10):3638–3663, 2019.
  * Cotter et al. (2013) Simon L Cotter, Gareth O Roberts, Andrew M Stuart, and David White.  MCMC methods for functions: modifying old algorithms to make them faster.  _Statistical Science_ , pages 424–446, 2013.
  * de Freitas et al. (2001) Nando de Freitas, Pedro A. d. F. R. Højen-Sørensen, and Stuart Russell.  Variational MCMC.  In _Proceedings of the 17th Conference in Uncertainty in Artificial Intelligence (UAI)_ , pages 120–127, Seattle, WA, 2001.
  * Deng et al. (2020) Yuntian Deng, Anton Bakhtin, Myle Ott, Arthur Szlam, and Marc’Aurelio Ranzato.  Residual energy-based models for text generation.  In _Proceedings of the 8th International Conference on Learning Representations (ICLR)_ , Addis Ababa, Ethiopia, 2020.
  * Doucet et al. (2000) Arnaud Doucet, Simon Godsill, and Christophe Andrieu.  On sequential monte carlo sampling methods for bayesian filtering.  _Statistics and Computing_ , 10(3):197–208, 2000\.
  * Du and Mordatch (2019) Yilun Du and Igor Mordatch.  Implicit generation and modeling with energy based models.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 3603–3613, Vancouver, Canada, 2019.
  * Du et al. (2020) Yilun Du, Joshua Meier, Jerry Ma, Rob Fergus, and Alexander Rives.  Energy-based models for atomic-resolution protein conformations.  In _Proceedings of the 8th International Conference on Learning Representations (ICLR)_ , Addis Ababa, Ethiopia, 2020.
  * Du et al. (2021) Yilun Du, Shuang Li, Joshua B. Tenenbaum, and Igor Mordatch.  Improved contrastive divergence training of energy-based models.  In _Proceedings of the 38th International Conference on Machine Learning (ICML)_ , pages 2837–2848, Virtual Event, 2021.
  * Durmus et al. (2017) Alain Durmus, Gareth O Roberts, Gilles Vilmart, and Konstantinos C Zygalakis.  Fast Langevin based algorithm for MCMC in high dimensions.  _The Annals of Applied Probability_ , 27(4):2195–2237, 2017.
  * Gao et al. (2018) Ruiqi Gao, Yang Lu, Junpei Zhou, Song-Chun Zhu, and Ying Nian Wu.  Learning generative convnets via multi-grid modeling and sampling.  In _Proceedings of the 2018 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 9155–9164, Salt Lake City, UT, 2018.
  * Gao et al. (2020) Ruiqi Gao, Erik Nijkamp, Diederik P. Kingma, Zhen Xu, Andrew M. Dai, and Ying Nian Wu.  Flow contrastive estimation of energy-based models.  In _Proceedings of the 2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 7515–7525, Seattle, WA, 2020.
  * Girolami and Calderhead (2011) Mark Girolami and Ben Calderhead.  Riemann manifold Langevin and Hamiltonian Monte Carlo methods.  _Journal of the Royal Statistical Society: Series B_ , 73(2):123–214, 2011.
  * Goodfellow et al. (2014) Ian J. Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron C. Courville, and Yoshua Bengio.  Generative adversarial nets.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 2672–2680, Montreal, Canada, 2014.
  * Grenander and 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.
  * Gustafsson et al. (2020) Fredrik K. Gustafsson, Martin Danelljan, Goutam Bhat, and Thomas B. Schön.  Energy-based models for deep probabilistic regression.  In _Proceedings of the 16th European Conference on Computer Viosion (ECCV), Part XX_ , pages 325–343, Glasgow, UK, 2020.
  * Gutmann and Hyvärinen (2010) Michael Gutmann and Aapo Hyvärinen.  Noise-contrastive estimation: A new estimation principle for unnormalized statistical models.  In _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 297–304, Chia Laguna Resort, Sardinia, Italy, 2010.
  * Haarnoja et al. (2017) Tuomas Haarnoja, Haoran Tang, Pieter Abbeel, and Sergey Levine.  Reinforcement learning with deep energy-based policies.  In _Proceedings of the 34th International Conference on Machine Learning (ICML)_ , pages 1352–1361, Sydney, Australia, 2017.
  * Han et al. (2017) Tian Han, Yang Lu, Song-Chun Zhu, and Ying Nian Wu.  Alternating back-propagation for generator network.  In _Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence (AAAI)_ , pages 1976–1984, San Francisco, CA, 2017.
  * Han et al. (2019) Tian Han, Erik Nijkamp, Xiaolin Fang, Mitch Hill, Song-Chun Zhu, and Ying Nian Wu.  Divergence triangle for joint training of generator model, energy-based model, and inferential model.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 8670–8679, Long Beach, CA, 2019.
  * Han et al. (2020) Tian Han, Erik Nijkamp, Linqi Zhou, Bo Pang, Song-Chun Zhu, and Ying Nian Wu.  Joint training of variational auto-encoder and latent energy-based model.  In _Proceedings of the 2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 7975–7984, Seattle, WA, 2020.
  * 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.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 6626–6637, Long Beach, CA, 2017.
  * Hinton (2002) Geoffrey E. Hinton.  Training products of experts by minimizing contrastive divergence.  _Neural Comput._ , 14(8):1771–1800, 2002.
  * Hyvärinen (2005) Aapo Hyvärinen.  Estimation of non-normalized statistical models by score matching.  _J. Mach. Learn. Res._ , 6:695–709, 2005.
  * Ingraham et al. (2019) John Ingraham, Adam J. Riesselman, Chris Sander, and Debora S. Marks.  Learning protein structure with a differentiable simulator.  In _Proceedings of the 7th International Conference on Learning Representationg (ICLR)_ , New Orleans, LA, 2019.
  * Jacob et al. (2020) Pierre E Jacob, John O Leary, and Yves F Atchadé.  Unbiased markov chain monte carlo methods with couplings.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 82(3):543–600, 2020.
  * Jarner and Hansen (2000) Søren Fiig Jarner and Ernst Hansen.  Geometric ergodicity of metropolis algorithms.  _Stochastic Processes and Their Applications_ , 85(2):341–361, 2000.
  * Jin et al. (2017) Long Jin, Justin Lazarow, and Zhuowen Tu.  Introspective classification with convolutional nets.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 823–833, Long Beach, CA, 2017.
  * Kingma and Ba (2015) Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In _Proceedings of the 3rd International Conference on Learning Representations (ICLR)_ , San Diego, CA, 2015.
  * Kingma and Welling (2014) Diederik P. Kingma and Max Welling.  Auto-encoding variational bayes.  In _Proceedings of the 2nd International Conference on Learning Representations (ICLR)_ , Banff, Canada, 2014.
  * Krizhevsky and Hinton (2009) Alex Krizhevsky and Geoffrey Hinton.  Learning multiple layers of features from tiny images.  _Technical Report, University of Toronto_ , 2009.
  * Lazarow et al. (2017) Justin Lazarow, Long Jin, and Zhuowen Tu.  Introspective neural networks for generative modeling.  In _Proceedings of the IEEE International Conference on Computer Vision (ICCV)_ , pages 2793–2802, Venice, Italy, 2017.
  * LeCun et al. (2006) Yann LeCun, Sumit Chopra, Raia Hadsell, M Ranzato, and F Huang.  A tutorial on energy-based learning.  _Predicting Structured Data_ , 1(0), 2006.
  * Lee et al. (2018) Kwonjoon Lee, Weijian Xu, Fan Fan, and Zhuowen Tu.  Wasserstein introspective neural networks.  In _Proceedings of the 2018 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 3702–3711, Salt Lake City, UT, 2018.
  * Li et al. (2019) Wenliang Li, Danica J. Sutherland, Heiko Strathmann, and Arthur Gretton.  Learning deep kernels for exponential family densities.  In _Proceedings of the 36th International Conference on Machine Learning (ICML)_ , pages 6737–6746, Long Beach, CA, 2019.
  * Liu et al. (2020) Weitang Liu, Xiaoyun Wang, John D. Owens, and Yixuan Li.  Energy-based out-of-distribution detection.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , virtual, 2020.
  * Liu et al. (2015) Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang.  Deep learning face attributes in the wild.  In _Proceedings of the 2015 IEEE International Conference on Computer Vision (ICCV)_ , pages 3730–3738, Santiago, Chile, 2015.
  * Marshall and Roberts (2012) Tristan Marshall and Gareth O Roberts.  An adaptive approach to Langevin MCMC.  _Statistics and Computing_ , 22(5), 2012.
  * Meyn and Tweedie (2012) Sean P Meyn and Richard L Tweedie.  _Markov Chains and Stochastic Stability_.  Springer Science & Business Media, 2012.
  * Mikolov et al. (2013) Tomás Mikolov, Ilya Sutskever, Kai Chen, Gregory S. Corrado, and Jeffrey Dean.  Distributed representations of words and phrases and their compositionality.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 3111–3119, Lake Tahoe, NV, 2013.
  * Neal (2011) Radford M Neal.  MCMC using hamiltonian dynamics.  _Handbook of Markov Chain Monte Carlo_ , 2(11):2, 2011.
  * Ngiam et al. (2011) Jiquan Ngiam, Zhenghao Chen, Pang Wei Koh, and Andrew Y. Ng.  Learning deep energy models.  In _Proceedings of the 28th International Conference on Machine Learning (ICML)_ , pages 1105–1112, Bellevue, WA, 2011.
  * Nijkamp et al. (2019) Erik Nijkamp, Mitch Hill, Song-Chun Zhu, and Ying Nian Wu.  Learning non-convergent non-persistent short-run MCMC toward energy-based model.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 5233–5243, Vancouver, Canada, 2019.
  * Nijkamp et al. (2020) Erik Nijkamp, Mitch Hill, Tian Han, Song-Chun Zhu, and Ying Nian Wu.  On the anatomy of MCMC-based maximum likelihood learning of energy-based models.  In _Proceedings of the Thirty-Fourth AAAI Conference on Artificial Intelligence (AAAI)_ , pages 5272–5280, New York, NY, 2020.
  * Nilsback and Zisserman (2008) Maria-Elena Nilsback and Andrew Zisserman.  Automated flower classification over a large number of classes.  In _Proceedings of the Sixth Indian Conference on Computer Vision, Graphics & Image Processing (ICVGIP)_, pages 722–729, Bhubaneswar, India, 2008.
  * Pang et al. (2020) Bo Pang, Tian Han, Erik Nijkamp, Song-Chun Zhu, and Ying Nian Wu.  Learning latent space energy-based prior model.  In _Advances in Neural Information Processing Systems (NeurIPS) 2020_ , virtual, 2020.
  * Qiu et al. (2020) Yixuan Qiu, Lingsong Zhang, and Xiao Wang.  Unbiased contrastive divergence algorithm for training energy-based latent variable models.  In _Proceedings of the 8th International Conference on Learning Representations (ICLR)_ , Addis Ababa, Ethiopia, 2020.
  * Robbins and Monro (1951) Herbert Robbins and Sutton Monro.  A stochastic approximation method.  _The Annals of Mathematical Statistics_ , pages 400–407, 1951.
  * Robert et al. (2010) Christian P Robert, George Casella, and George Casella.  _Introducing Monte Carlo methods with R_ , volume 18.  Springer, 2010.
  * Roberts and Rosenthal (1998) Gareth O Roberts and Jeffrey S Rosenthal.  Optimal scaling of discrete approximations to Langevin diffusions.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 60(1):255–268, 1998.
  * Roberts and Rosenthal (2004) Gareth O Roberts and Jeffrey S Rosenthal.  General state space Markov chains and MCMC algorithms.  _Probability Surveys_ , 1:20–71, 2004.
  * Roberts and Tweedie (1996) Gareth O Roberts and Richard L Tweedie.  Exponential convergence of Langevin distributions and their discrete approximations.  _Bernoulli_ , pages 341–363, 1996.
  * Rue et al. (2009) Håvard Rue, Sara Martino, and Nicolas Chopin.  Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 2009.
  * Song et al. (2023) Weinan Song, Yaxuan Zhu, Lei He, Ying Nian Wu, and Jianwen Xie.  Progressive energy-based cooperative learning for multi-domain image-to-image translation.  _arXiv preprint arXiv:2306.14448_ , 2023.
  * Song et al. (2019) Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon.  Sliced score matching: A scalable approach to density and score estimation.  In _Proceedings of the Thirty-Fifth Conference on Uncertainty in Artificial Intelligence (UAI)_ , pages 574–584, Tel Aviv, Israel, 2019.
  * Tieleman (2008) Tijmen Tieleman.  Training restricted boltzmann machines using approximations to the likelihood gradient.  In _Proceedings of the Twenty-Fifth International Conference on Machine Learning (ICML)_ , pages 1064–1071, Helsinki, Finland, 2008.
  * Wainwright and Jordan (2008) Martin J. Wainwright and Michael I. Jordan.  Graphical models, exponential families, and variational inference.  _Found. Trends Mach. Learn._ , 1(1-2):1–305, 2008\.
  * Welling and Hinton (2002) Max Welling and Geoffrey E. Hinton.  A new learning algorithm for mean field boltzmann machines.  In _Proceedings of the International Conference on Artificial Neural Networks (ICANN)_ , pages 351–357, Madrid, Spain, 2002.
  * Welling and Teh (2011) Max Welling and Yee Whye Teh.  Bayesian learning via stochastic gradient langevin dynamics.  In _Proceedings of the 28th International Conference on Machine Learning (ICML)_ , pages 681–688, Bellevue, WA, 2011.
  * Wolfinger (1993) Russ Wolfinger.  Laplace’s approximation for nonlinear mixed models.  _Biometrika_ , 80(4):791–795, 1993.
  * Xie et al. (2016) Jianwen Xie, Yang Lu, Song-Chun Zhu, and Ying Nian Wu.  A theory of generative ConvNet.  In _Proceedings of the 33nd International Conference on Machine Learning (ICML)_ , pages 2635–2644, New York City, NY, 2016.
  * Xie et al. (2017) Jianwen Xie, Song-Chun Zhu, and Ying Nian Wu.  Synthesizing dynamic patterns by spatial-temporal generative convnet.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 7093–7101, 2017.
  * Xie et al. (2018a) Jianwen Xie, Yang Lu, Ruiqi Gao, and Ying Nian Wu.  Cooperative learning of energy-based model and latent variable model via MCMC teaching.  In _Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, (AAAI)_ , pages 4292–4301, New Orleans, LA, 2018a. AAAI Press.
  * Xie et al. (2018b) Jianwen Xie, Zilong Zheng, Ruiqi Gao, Wenguan Wang, Song-Chun Zhu, and Ying Nian Wu.  Learning descriptor networks for 3D shape synthesis and analysis.  In _Proceedings of the 2018 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 8629–8638, Salt Lake City, UT, 2018b.
  * Xie et al. (2019) Jianwen Xie, Ruiqi Gao, Zilong Zheng, Song-Chun Zhu, and Ying Nian Wu.  Learning dynamic generator model by alternating back-propagation through time.  In _Proceedings of the Thirty-Third AAAI Conference on Artificial Intelligence (AAAI)_ , pages 5498–5507, Honolulu, HI, 2019.
  * Xie et al. (2020) Jianwen Xie, Yang Lu, Ruiqi Gao, Song-Chun Zhu, and Ying Nian Wu.  Cooperative training of descriptor and generator networks.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 42(1):27–45, 2020.
  * Xie et al. (2021a) Jianwen Xie, Yifei Xu, Zilong Zheng, Song-Chun Zhu, and Ying Nian Wu.  Generative PointNet: Deep energy-based learning on unordered point sets for 3D generation, reconstruction and classification.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 14976–14985, virtual, 2021a.
  * Xie et al. (2021b) Jianwen Xie, Zilong Zheng, Xiaolin Fang, Song-Chun Zhu, and Ying Nian Wu.  Learning cycle-consistent cooperative networks via alternating MCMC teaching for unsupervised cross-domain translation.  In _Proceedings of the Thirty-Fifth AAAI Conference on Artificial Intelligence (AAAI)_ , pages 10430–10440, Virtual Event, 2021b.
  * Xie et al. (2021c) Jianwen Xie, Zilong Zheng, and Ping Li.  Learning energy-based model with variational auto-encoder as amortized sampler.  In _Proceedings of the Thirty-Fifth AAAI Conference on Artificial Intelligence (AAAI)_ , pages 10441–10451, Virtual Event, 2021c.
  * Xie et al. (2021d) Jianwen Xie, Song-Chun Zhu, and Ying Nian Wu.  Learning energy-based spatial-temporal generative ConvNets for dynamic patterns.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 43(2):516–531, 2021d.
  * Xie et al. (2022a) Jianwen Xie, Zilong Zheng, Xiaolin Fang, Song-Chun Zhu, and Ying Nian Wu.  Cooperative training of fast thinking initializer and slow thinking solver for conditional learning.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 44(8):3957–3973, 2022a.
  * Xie et al. (2022b) Jianwen Xie, Zilong Zheng, Ruiqi Gao, Wenguan Wang, Song-Chun Zhu, and Ying Nian Wu.  Generative VoxelNet: Learning energy-based models for 3D shape synthesis and analysis.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 44(5):2468–2484, 2022b.
  * Xie et al. (2022c) Jianwen Xie, Yaxuan Zhu, Jun Li, and Ping Li.  A tale of two flows: Cooperative learning of langevin flow and normalizing flow toward energy-based model.  In _Proceedings of the Tenth International Conference on Learning Representations (ICLR)_ , Virtual Event, 2022c.
  * Xie et al. (2023) Jianwen Xie, Yaxuan Zhu, Yifei Xu, Dingcheng Li, and Ping Li.  A tale of two latent flows: Learning latent space normalizing flow with short-run langevin flow for approximate inference.  In _Proceedings of the Thirty-Seventh AAAI Conference on Artificial Intelligence (AAAI)_ , pages 10499–10509, Washington, DC, 2023.
  * Xu et al. (2022) Yifei Xu, Jianwen Xie, Tianyang Zhao, Chris Baker, Yibiao Zhao, and Ying Nian Wu.  Energy-based continuous inverse optimal control.  _IEEE Trans. Neural Networks Learn. Syst. (Early Access)_ , 2022.
  * Yu et al. (2020) Lantao Yu, Yang Song, Jiaming Song, and Stefano Ermon.  Training deep energy-based models with f-divergence minimization.  In _Proceedings of the 37th International Conference on Machine Learning (ICML)_ , pages 10957–10967, Virtual Event, 2020.
  * Zhang et al. (2020) Jing Zhang, Jianwen Xie, and Nick Barnes.  Learning noise-aware encoder-decoder from noisy labels by alternating back-propagation for saliency detection.  In _Proceedings of the 16th European Conference on Computer Viosion (ECCV), Part XVII_ , pages 349–366, Glasgow, UK, 2020.
  * Zhang et al. (2021) Jing Zhang, Jianwen Xie, Nick Barnes, and Ping Li.  Learning generative vision transformer with energy-based latent space for saliency prediction.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 15448–15463, virtual, 2021.
  * Zhang et al. (2022) Jing Zhang, Jianwen Xie, Zilong Zheng, and Nick Barnes.  Energy-based generative cooperative saliency prediction.  In _Proceedings of the Thirty-Sixth AAAI Conference on Artificial Intelligence (AAAI)_ , pages 3280–3290, Virtual Event, 2022.
  * Zhang et al. (2023) Jing Zhang, Jianwen Xie, Nick Barnes, and Ping Li.  An energy-based prior for generative saliency.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 45(11):13100–13116, 2023.
  * Zhao et al. (2021) Yang Zhao, Jianwen Xie, and Ping Li.  Learning energy-based generative models via coarse-to-fine expanding and sampling.  In _Proceedings of the 9th International Conference on Learning Representations (ICLR)_ , Virtual Event, 2021.
  * Zhao et al. (2023) Yang Zhao, Jianwen Xie, and Ping Li.  CoopInit: Initializing generative adversarial networks via cooperative learning.  In _Proceedings of the Thirty-Seventh AAAI Conference on Artificial Intelligence (AAAI)_ , pages 11345–11353, Washington, DC, 2023.
  * Zheng et al. (2021) Zilong Zheng, Jianwen Xie, and Ping Li.  Patchwise generative ConvNet: Training energy-based models from a single natural image for internal learning.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 2961–2970, virtual, 2021.
  * Zhu et al. (1998) Song Chun Zhu, Ying Nian Wu, and David Mumford.  Filters, random fields and maximum entropy (FRAME): towards a unified theory for texture modeling.  _Int. J. Comput. Vis._ , 27(2):107–126, 1998\.
  * Zhu et al. (2023) Yaxuan Zhu, Jianwen Xie, and Ping Li.  Likelihood-based generative radiance field with latent space energy-based model for 3D-aware disentangled image representation.  In _Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 4164–4180, Palau de Congressos, Valencia, Spain, 2023.

##  Appendix A Appendix: Proofs of the Theoretical Results

###  A.1 Sketch of the Proof of Theorem 4.2

Notations for the proof: We denote by $z\to\textsc{t}_{\theta}(z^{\prime},z)$,
the pdf of the Gaussian proposal of Line 3 for any current state of the chain
$z^{\prime}\in\mathcal{Z}$ and dependent on the EBM model parameter $\theta$.
The transition kernel from $z$ to $z^{\prime}$ is denoted by
$\Pi_{\theta}(z,z^{\prime})$. $\mathcal{Z}$ is a subset of $\mathbb{R}^{\ell}$
and $\mathcal{B}$ is a Borel set of $\mathbb{R}^{\ell}$.

The proof of our results are divided into two main parts. We first prove the
existence of a small set for our transition kernel $\Pi_{\theta}$, noted
$\mathcal{O}$ showing that for any state, the Markov Chain moves away from it.
It constitutes the first step toward proving its irreducibility and
aperiodicity. Then, we will establish the so-called _drift condition_ , also
known as the Foster-Lyapunov condition, crucial to proving the convergence of
the chain. The drift condition ensures the recurrence of the chain as the
property that a chain returns to its initial state within finite time, see
Roberts and Rosenthal (1998, 2004) for more details. Uniform ergodicity is
then established as a consequence of those drift conditions and thus proving
(9).

(i) Existence of a small set:  Let $\mathcal{O}$ be a compact subset of the
state space $\mathcal{Z}$. We recall the definition of the transition kernel
in the case of a Metropolis adjustment and for any model parameter
$\theta\in\Theta$ and state $z\in\mathcal{Z}$:

$\begin{split}\Pi_{\theta}(z,\mathcal{B})=&\int_{\mathcal{B}}\alpha_{\theta}(z,y)\textsc{t}_{\theta}(z,y)\textrm{d}y+\mathsf{1}_{\mathcal{B}(z)}\int_{\mathcal{Z}}(1-\alpha_{\theta}(z,y))\textsc{t}_{\theta}(z,y)\textrm{d}y\;,\end{split}$

where we have defined the Metropolis ratio between two states
$(z,y)\in\mathcal{Z}\times\mathcal{B}$ as
$\alpha_{\theta}(z,y)=\textrm{min}(1,\frac{\pi_{\theta}(z)\textsc{t}_{\theta}(z,y)}{\textsc{t}_{\theta}(y,z)\pi_{\theta}(y)})$.
Under H1 and due to the fact that the threshold th leads to a symmetric
positive definite covariance matrix with bounded non zero eigenvalues, then
the following holds:

$an_{\sigma_{1}}(z-y)\leq\textsc{t}_{\theta}(z,y)\leq bn_{\sigma_{2}}(z-y)\;,$
(11)

for all $\theta\in\Theta$ and where $\sigma_{1}$ and $\sigma_{2}$ are the
corresponding standard deviations of the two Gaussian distributions
$n_{\sigma_{1}}$ and $n_{\sigma_{2}}$. We denote by $\rho_{\theta}$ the ratio
$\frac{\pi_{\theta}(z)\textsc{t}_{\theta}(z,y)}{\textsc{t}_{\theta}(y,z)\pi_{\theta}(y)}$
and define the quantity

$\delta=\textrm{inf}(\rho_{\theta}(z,y),\theta\in\Theta,z\in\mathcal{O})>0\;,$
(12)

where we have used assumptions H1 and H2. Then,

$\begin{split}\Pi_{\theta}(z,\mathcal{B})&\geq\int_{\mathcal{B}\cap\mathcal{X}}\alpha_{\theta}(z,y)\textsc{t}_{\theta}(z,y)\textrm{d}y\geq\textrm{min}(1,\delta)m\int_{\mathcal{B}}\mathsf{1}_{\mathcal{X}}(z)\textrm{d}y\;.\end{split}$

According to (12), we can find a compact set $\mathcal{O}$ such that
$\Pi_{\theta}(z,\mathcal{B})\geq\epsilon$ where
$\epsilon=\textrm{min}(1,\delta)m\textbf{Z}$ where Z is the normalizing
constant of the pdf
$\frac{1}{\textbf{Z}}\mathbf{1}_{\mathcal{X}}(z)\textrm{d}y$ and the proposal
distribution is bounded from below by some quantity noted $m$. The
calculations above prove (8), i.e., the existence of a small set for our
family of transition kernels $(\Pi_{\theta})_{\theta\in\Theta}$.

(ii) Drift condition and ergodicity:  We begin by proving that
$(\Pi_{\theta})_{\theta\in\Theta}$ satisfies a drift property. For a given EBM
parameter $\theta\in\Theta$, we can see in Jarner and Hansen (2000) that the
drift condition boils down to proving that

$\sup\limits_{z\in\mathcal{Z}}\frac{\Pi_{\theta}V_{\theta}(z)}{V_{\theta}(z)}<\infty\quad\textrm{and}\quad\lim\sup\limits_{|z|\to\infty}\frac{\Pi_{\theta}V_{\theta}(z)}{V_{\theta}(z)}<1\;,$

where $V_{\theta}$ is the _drift function_ defined in (6) Let denote the
acceptation set, i.e., $\rho_{\theta}\geq 1$ by

$\mathcal{A}_{\theta}(z)\mathrel{\mathop{:}}=\\{y\in\mathcal{Z},\rho_{\theta}(z,y)\geq
1\\}$ (13)

for any state $y\in\mathcal{B}$ and its complementary set
$\mathcal{A}^{*}_{\theta}(z)$. The remaining of the proof is composed of three
main steps. Step (1) shows that for any $\theta\in\Theta$,

$\displaystyle\lim\sup\limits_{|z|\to\infty}\frac{\Pi_{\theta}V_{\theta}(z)}{V_{\theta}(z)}\leq
1-\lim\inf\limits_{|z|\to\infty}\int_{\mathcal{A}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\textrm{d}y\;.$

where the smoothness of a Gaussian pdf, assumption H1 and a collection of
inequalities based on (13), its complentary set and the interval in (11).
Then, using an important intermediary result, stated in Lemma A.1, that
initiates a relation between the set of accepted states noted
$\mathcal{A}_{\theta}(z)$ and the cone $\mathcal{P}(z)$ designed so that it
does not depend on the model parameter $\theta$.

Define
$\mathcal{P}(z)\mathrel{\mathop{:}}=\\{z-\ell\frac{z}{|z|}-\kappa\nu,\,\textrm{with}\quad\kappa<a-\ell,\nu\in\\{\nu\in\mathbb{R}^{d},\|\nu\|<1\\},|\nu-\frac{z-\ell\frac{z}{|z|}}{|z-\ell\frac{z}{|z|}|}\leq\frac{\epsilon}{2}\\}$
and
$\mathcal{A}_{\theta}(z)\mathrel{\mathop{:}}=\\{y\in\mathcal{Z},\rho_{\theta}(z,y)\geq
1\\}$. Then for $z\in\mathcal{Z}$,
$\mathcal{P}(z)\subset\mathcal{A}_{\theta}(z)$.

Noting the limit inferior as $\varliminf$, Step (2) establishes that
$1-\varliminf\limits_{|z|\to\infty}\int_{\mathcal{A}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\textrm{d}y\leq
1-c$ where $c$ is a constant, _independent of all the other quantities_
towards showing uniformity of the final result. Finally, Step (3) uses the
inequality
$\Pi_{\theta}V_{\theta}(z)\leq\bar{\mu}V_{\theta}(z)+\bar{\delta}\mathsf{1}_{\mathcal{O}}(z)$
dependent of $\theta$ and defines the V function, independent of $\theta$, as
$V(z)\mathrel{\mathop{:}}=V_{1}(z)^{\alpha}V_{2}(z)^{2\alpha}$ in order to
establish the main result of Theorem 4.2, i.e.,

$\displaystyle\Pi_{\theta}V(z)\leq\left(\frac{\bar{\mu}}{2\epsilon^{2}}+\frac{\epsilon^{2}}{1+\bar{\mu}}\right)V(z)+\frac{\bar{\delta}}{2\epsilon^{2}}\mathsf{1}_{\mathcal{O}}(z)\;.$

Setting
$\epsilon\mathrel{\mathop{:}}=\sqrt{\frac{\bar{\mu}(1+\bar{\mu})}{2}}$,
$\mu\mathrel{\mathop{:}}=\sqrt{\frac{2\bar{\mu}}{1+\bar{\mu}}}$ and
$\delta\mathrel{\mathop{:}}=\frac{\bar{\delta}}{2\epsilon^{2}}$ proves the
uniformity of the inequality (9).

The complete proof is deferred in the appendix and is also developed in
Allassonniere and Kuhn (2015) in the context of Bayesian Mixed Effect models
trained with the EM algorithm.

###  A.2 Proof of Theorem 4.2

######  Theorem.

Assume H1-H3. For any $\theta\in\Theta$, there exists a drift function
$V_{\theta}$, a set $\mathcal{O}\subset\mathcal{Z}$, a constant
$0<\epsilon\leq 1$ such that

$\Pi_{\theta}(z,\mathcal{B})\geq\epsilon\int_{\mathcal{B}}\mathsf{1}_{\mathcal{X}}(z)\textrm{d}y\;.$
(14)

Moreover there exists $0<\mu<1$, $\delta>0$ and a drift function $V$, now
independent of $\theta$ such that for all $z\mathcal{Z}$:

$\Pi_{\theta}V(z)\leq\mu V(z)+\delta\mathsf{1}_{\mathcal{O}}(z)\;.$ (15)

###### Proof.

Notations used throughout the proof are listed in Table 1

Table 1: Notations $\Pi_{\theta}$ | $\triangleq$ |  Transition kernel of the MCMC defined by (5)
---|---|---
$\mathcal{O}$ | $\triangleq$ |  Subset of $\mathbb{R}^{p}$ and small set for kernel $\Pi_{\theta}$
$B(z,a)$ | $\triangleq$ |  Ball around $z\in\mathcal{Z}$ of radius $a>0$
$\mathcal{A}_{\theta}(z)$ | $\triangleq$ |  Acceptance set at state $\in\mathcal{Z}$ such that $\rho_{\theta}\geq 1$
$\mathcal{A}^{*}_{\theta}(z)$ | $\triangleq$ |  Complementary set of $\mathcal{A}_{\theta}(z)$
$\textsc{t}_{\theta}(z^{\prime},z)$ | $\triangleq$ |  Probability density function of the Gaussian proposal
$\pi_{\theta}(\cdot)$ | $\triangleq$ |  Stationary/Target distribution under model $\theta\in\Theta$
$\Pi_{\theta}(z,z^{\prime})$ | $\triangleq$ |  Transition kernel from state $z$ to state $z^{\prime}$
$n_{\sigma}(z)$ | $\triangleq$ |  Pdf of a centered Normal distribution of standard deviation $\sigma>0$

The proof of our results are divided into two parts. We first prove the
existence of a set noted $\mathcal{O}$ as a small set for our family of
transition kernels $(\Pi_{\theta})_{\theta\in\Theta}$. Proving a small set is
crucial in order to show that for any state, the Markov Chain does not stay in
the same state, and thus help in proving its irreducibility and aperiodicity.

Then, we will prove the drift condition towards a small set. This condition is
crucial to prove the convergence of the chain since it states that the kernels
tend to attract elements into that set. finally, uniform ergodicity is
established as a consequence of those drift conditions.

(i) Existence of small set:  Let $\mathcal{O}$ be a compact subset of the
state space $\mathcal{Z}$. We also denote the probability density function
(pdf) of the Gaussian proposal of Line 3 as
$z\to\textsc{t}_{\theta}(z^{\prime},z)$ for any current state of the chain
$z^{\prime}\in\mathcal{Z}$ and dependent on the EBM model parameter $\theta$.
Given STANLEY’s MCMC update, at iteration $t$, the proposal is a Gaussian
distribution of mean $z_{t-1}^{m}+\mathsf{\gamma}_{t}/2\nabla
f_{\theta_{t}}(z_{t-1}^{m})$ and covariance
$\sqrt{\mathsf{\gamma}_{t}}\mathsf{B}_{t}$.

We recall the definition of the transition kernel in the case of a Metropolis
adjustment and for any model parameter $\theta\in\Theta$ and state
$z\in\mathcal{Z}$:

$\Pi_{\theta}(z,\mathcal{B})=\int_{\mathcal{B}}\alpha_{\theta}(z,y)\textsc{t}_{\theta}(z,y)\textrm{d}y+\mathsf{1}_{\mathcal{B}(z)}\int_{\mathcal{Z}}(1-\alpha_{\theta}(z,y))\textsc{t}_{\theta}(z,y)\textrm{d}y\;,$
(16)

where we have defined the Metropolis ratio between two states
$z\in\mathcal{Z}$ and $y\in\mathcal{B}$ as
$\alpha_{\theta}(z,y)=\textrm{min}(1,\frac{\pi_{\theta}(z)\textsc{t}_{\theta}(z,y)}{\textsc{t}_{\theta}(y,z)\pi_{\theta}(y)})$.
Thanks to Assumption H1 and due to the fact that the threshold th leads to a
symmetric positive definite covariance matrix with bounded non zero
eigenvalues implies that the proposal distribution can be bounded by two zero-
mean Gaussian distributions as follows:

$an_{\sigma_{1}}(z-y)\leq\textsc{t}_{\theta}(z,y)\leq
bn_{\sigma_{2}}(z-y)\quad\textrm{for all}\quad\theta\in\Theta\;,$ (17)

where $\sigma_{1}$ and $\sigma_{2}$ are the corresponding standard deviation
of the distributions and $a$ and $b$ are some scaling factors.

We denote by $\rho_{\theta}$ the ratio
$\frac{\pi_{\theta}(z)\textsc{t}_{\theta}(z,y)}{\textsc{t}_{\theta}(y,z)\pi_{\theta}(y)}$
and given the assumptions H1 and H2, define the quantity

$\delta=\textrm{inf}(\rho_{\theta}(z,y),\theta\in\Theta,\quad
z\in\mathcal{O})>0\;.$ (18)

Likewise, the proposal distribution is bounded from below by some quantity
noted $m$. Then,

$\Pi_{\theta}(z,\mathcal{B})\geq\int_{\mathcal{B}\cap\mathcal{X}}\alpha_{\theta}(z,y)\textsc{t}_{\theta}(z,y)\textrm{d}y\geq\textrm{min}(1,\delta)m\int_{\mathcal{B}}\mathsf{1}_{\mathcal{X}}(z)\textrm{d}y\;.$
(19)

Then, given the definition of (18), we can find a compact set $\mathcal{O}$
such that $\Pi_{\theta}(z,\mathcal{B})\geq\epsilon$ where
$\epsilon=\textrm{min}(1,\delta)m\textbf{Z}$ where Z is the normalizing
constant of the pdf
$\frac{1}{\textbf{Z}}\mathbf{1}_{\mathcal{X}}(z)\textrm{d}y$. The calculations
above prove (8), i.e., the existence of a small set for our family of
transition kernels $(\Pi_{\theta})_{\theta\in\Theta}$.

(ii) Drift condition and ergodicity:  We first need to prove the fact that our
family of transition kernels $(\Pi_{\theta})_{\theta\in\Theta}$ satisfies a
drift property.

For a given EBM model parameter $\theta\in\Theta$, we can see in Jarner and
Hansen (2000) that the drift condition boils down to proving that for the
drift function noted $V_{\theta}$ and defined in (6), we have

$\sup\limits_{z\in\mathcal{Z}}\frac{\Pi_{\theta}V_{\theta}(z)}{V_{\theta}(z)}<\infty\quad\textrm{and}\quad\lim\sup\limits_{|z|\to\infty}\frac{\Pi_{\theta}V_{\theta}(z)}{V_{\theta}(z)}<1\;.$
(20)

Throughout the proof, the model parameter is set to an arbitrary
$\theta\in\Theta$. Let denote the acceptation set, i.e., $\rho_{\theta}\geq 1$
by
$\mathcal{A}_{\theta}(z)\mathrel{\mathop{:}}=\\{y\in\mathcal{Z},\rho_{\theta}(z,y)\geq
1\\}$ for any state $y\in\mathcal{B}$ and its complementary set
$\mathcal{A}^{*}_{\theta}(z)$.

Step (1):  Following our definition of the drift function in (6) we obtain:

$\displaystyle\frac{\Pi_{\theta}V_{\theta}(z)}{V_{\theta}(z)}$
$\displaystyle=\int_{\mathcal{A}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\frac{V_{\theta}(y)}{V_{\theta}(z)}\textrm{d}y+\int_{\mathcal{A}^{*}_{\theta}(z)}\frac{\pi_{\theta}(y)\textsc{t}_{\theta}(y,z)}{\pi_{\theta}(z)\textsc{t}_{\theta}(z,y)}\textsc{t}_{\theta}(z,y)\frac{V_{\theta}(y)}{V_{\theta}(z)}\textrm{d}y$
(21)
$\displaystyle+\int_{\mathcal{A}^{*}_{\theta}(z)}(1-\frac{\pi_{\theta}(y)\textsc{t}_{\theta}(y,z)}{\pi_{\theta}(z)\textsc{t}_{\theta}(z,y)})\textsc{t}_{\theta}(z,y)\textrm{d}y$
$\displaystyle\overset{(a)}{\leq}\int_{\mathcal{A}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\frac{\pi_{\theta}(y)^{-\beta}}{\pi_{\theta}(z)^{-\beta}}\textrm{d}y+\int_{\mathcal{A}^{*}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\frac{\pi_{\theta}(y)^{1-\beta}}{\pi_{\theta}(z)^{1-\beta}}\textrm{d}y+\int_{\mathcal{A}^{*}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\textrm{d}y\;,$
(22)

where (a) is due to (6).

Furthermore, according to (17), we thus have that, for any state $z$ in the
acceptance set $\mathcal{A}_{\theta}(z)$:

$\int_{\mathcal{A}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\frac{\pi_{\theta}(y)^{-\beta}}{\pi_{\theta}(z)^{-\beta}}\textrm{d}y\leq
b\int_{\mathcal{A}_{\theta}(z)}n_{\sigma_{2}}(y-z)\textrm{d}y\;.$ (23)

For any state $z$ in the complementary set of the acceptance set, noted
$\mathcal{A}^{*}_{\theta}(z)$, we also have the following:

$\int_{\mathcal{A}^{*}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\frac{\pi_{\theta}(y)^{1-\beta}}{\pi_{\theta}(z)^{1-\beta}}\textrm{d}y\leq\int_{\mathcal{A}^{*}_{\theta}(z)}\textsc{t}_{\theta}(z,y)^{1-\beta}\textsc{t}_{\theta}(y,z)^{\beta}\textrm{d}y\leq
b\int_{\mathcal{A}^{*}_{\theta}(z)}n_{\sigma_{2}}(z-y)\textrm{d}y\;.$ (24)

While we can define the level set of the stationary distribution
$\pi_{\theta}$ as
$\mathcal{L}_{\pi_{\theta}(y)}=\\{z\in\mathcal{Z},\pi_{\theta}(z)=\pi_{\theta}(y)\\}$
for some state $y\in\mathcal{B}$, a neighborhood of that level set is defined
as
$\mathcal{L}_{\pi_{\theta}(y)}(p)=\\{z\in\mathcal{L}_{\pi_{\theta}(y)},z+t\frac{z}{|z|},|t|\leq
p\\}$. H1 ensures the existence of a radial $r$ such that for all
$z\in\mathcal{Z},|z|\geq r$, then $0\in\mathcal{L}_{\pi_{\theta}(y)}$ with
$\pi_{\theta}(z)>\pi_{\theta}(y)$. Since the function $y\to
n_{\sigma_{2}}(y-z)$ is smooth, it is known that there exists a constant $a>0$
such that for $\epsilon>0$, we have that

$\int_{B(z,a)}n_{\sigma_{2}}(y-z)\textrm{d}y\geq
1-\epsilon\quad\textrm{and}\quad\int_{B(z,a)\cap\mathcal{L}_{\pi_{\theta}(y)}(p)}n_{\sigma_{2}}(y-z)\textrm{d}y\leq\epsilon\;,$
(25)

for some $p$ small enough and where $B(z,a)$ denotes the ball around
$z\in\mathcal{Z}$ of radius $a$. Then combining (23) and (25) we have that:

$\int_{\mathcal{A}_{\theta}(z)\cap
B(z,a)\cap\mathcal{L}_{\pi_{\theta}(y)}(p)}\textsc{t}_{\theta}(z,y)\frac{\pi_{\theta}(y)^{-\beta}}{\pi_{\theta}(z)^{-\beta}}\textrm{d}y\leq
b\epsilon\;.$ (26)

Conversely, we can define the following set
$\mathcal{A}=\mathcal{A}_{\theta}(z)\cap B(z,a)\cap\mathcal{L}^{+}$ where
$u\in\mathcal{L}^{+}$ if $u\in\mathcal{L}_{\pi_{\theta}(y)}(p)$ and
$\phi_{\theta}(u)>\pi_{\theta}(p)$. Then using the second part of H1, there
exists a radius $r^{\prime}>r+a$, such that for $z\in\mathcal{Z}$ with
$|z|\geq r^{\prime}$ we have

$\int_{\mathcal{A}}(\frac{\pi_{\theta}(y)}{\pi_{\theta}(z)})^{1-\beta}\textsc{t}_{\theta}(y,z)\textrm{d}y\leq\mathsf{d}(p,r^{\prime})^{1-\beta}b\int_{\mathcal{A}_{\theta}(z)}n_{\sigma_{2}}(y-z)\textrm{d}y\leq
b\mathsf{d}(p,r^{\prime})^{1-\beta}\;,$ (27)

where
$\mathsf{d}(p,r^{\prime})=\sup\limits_{|z|>r^{\prime}}\frac{\pi_{\theta}(z+p\frac{z}{|z|})}{\pi_{\theta}(z)}$.
Note that H1 implies that $\mathsf{d}(p,r^{\prime})\to 0$ when
$r^{\prime}\to\infty$. Likewise with $\mathcal{A}=\mathcal{A}_{\theta}(z)\cap
B(z,a)\cap\mathcal{L}^{-}$ we have

$\int_{\mathcal{A}}(\frac{\pi_{\theta}(y)}{\pi_{\theta}(z)})^{-\beta}\textsc{t}_{\theta}(z,y)\textrm{d}y\leq
b\mathsf{d}(p,r^{\prime})^{\beta}\;.$ (28)

Same arguments can be obtained for the second term of (21), i.e.,
$\textsc{t}_{\theta}(z,y)\frac{\pi_{\theta}(y)^{1-\beta}}{\pi_{\theta}(z)^{1-\beta}}$
and we obtain, plugging the above in (21) that:

$\displaystyle\lim\sup\limits_{|z|\to\infty}\frac{\Pi_{\theta}V_{\theta}(z)}{V_{\theta}(z)}\leq\lim\sup\limits_{|z|\to\infty}\int_{\mathcal{A}^{*}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\textrm{d}y\;.$
(29)

Since $\mathcal{A}^{*}_{\theta}(z)$ is the complementary set of
$\mathcal{A}_{\theta}(z)$, the above inequality yields

$\displaystyle\lim\sup\limits_{|z|\to\infty}\frac{\Pi_{\theta}V_{\theta}(z)}{V_{\theta}(z)}\leq
1-\lim\inf\limits_{|z|\to\infty}\int_{\mathcal{A}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\textrm{d}y\;.$
(30)

Step (2):  The final step of our proof consists in proving that
$1-\lim\inf\limits_{|z|\to\infty}\int_{\mathcal{A}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\textrm{d}y\leq
1-c$ where $c$ is a constant, independent of all the other quantities.

Given that the proposal distribution is a Gaussian and using assumption H1 we
have the existence of a constant $c_{a}$ depending on $a$ as defined above
(the radius of the ball $B(z,a)$) such that

$\frac{\pi_{\theta}(z)}{\pi_{\theta}(z-\ell\frac{z}{|z|})}\leq
c_{a}\leq\inf\limits_{y\in
B(z,a)}\frac{\textsc{t}_{\theta}(y,z)}{\textsc{t}_{\theta}(z,y)}\quad\textrm{for
any}\,z\in\mathcal{Z},|z|\geq r^{*}\;.$ (31)

Then for any $|z|\geq r^{*}$, we obtain that
$z-\ell\frac{z}{|z|}\in\mathcal{A}_{\theta}(z)$. A particular subset of
$\mathcal{A}_{\theta}(z)$ used throughout the rest of the proof is the cone
defined as

$\mathcal{P}(z)\mathrel{\mathop{:}}=\\{z-\ell\frac{z}{|z|}-\kappa\nu,\,\textrm{with}\quad
i<a-\ell,\nu\in\\{\nu\in\mathbb{R}^{d},\|\nu\|<1\\},|\nu-\frac{z-\ell\frac{z}{|z|}}{|z-\ell\frac{z}{|z|}|}\leq\frac{\epsilon}{2}\\}\;.$
(32)

Using Lemma A.1, we have that $\mathcal{P}(z)\subset\mathcal{A}_{\theta}(z)$.
Then, we observe that

$\int_{\mathcal{A}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\textrm{d}y\overset{(a)}{\geq}\int_{\mathcal{A}_{\theta}(z)}an_{\sigma_{1}}(y-z)\textrm{d}y\overset{(b)}{\geq}a\int_{\mathcal{P}(z)}n_{\sigma_{1}}(y-z)\textrm{d}y\;,$
(33)

where we have used (17) in (a) and applied Lemma A.1 in (b).

If we define the translation of vector $z\in\mathcal{Z}$ by the operator
$\mathcal{I}\subset\mathbb{R}^{d}\to T_{z}(\mathcal{I})$, then

$\int_{\mathcal{A}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\textrm{d}y\geq
a\int_{\mathcal{P}(z)}n_{\sigma_{1}}(y-z)\textrm{d}y=\int_{T_{z}(\mathcal{P}(z))}n_{\sigma_{1}}(y-z)\textrm{d}y\;.$
(34)

Recalling the objective of Step (2) that is to find a constant $c$ such that
$1-\lim\inf\limits_{|z|\to\infty}\int_{\mathcal{A}_{\theta}(z)}\textsc{t}_{\theta}(z,y)\textrm{d}y\leq
1-c$, we deduce from (34) that since the set $\mathcal{P}(z)$ does not depend
on the EBM model parameter $\theta$ and that once translated by $z$ the
resulting set $T_{z}(\mathcal{P}(z))$ is independent of $z$ (but depends on
$\ell$, see definition (32), then the integral
$\int_{T_{z}(\mathcal{P}(z))}n_{\sigma_{1}}(y-z)\textrm{d}y$ in (34) is
independent of $z$ thus concluding on the existence of the constant $c$ such
that

$\lim\sup\limits_{|z|\to\infty}\frac{\Pi_{\theta}V_{\theta}(z)}{V_{\theta}(z)}\leq
1-c\;.$

Thus proving the second part of (20) which is the main drift condition we
ought to demonstrate. The first part of (20) can be proved by observing that
$\frac{\Pi_{\theta}V_{\theta}(z)}{V_{\theta}(z)}$ is smooth on $\mathcal{Z}$
according to H2 and by construction of the transition kernel. Smoothness
implies boundedness on the compact $\mathcal{Z}$.

Step (3):  We now use the main proven equations in (20) to derive the second
result (9) of Theorem 4.2.

We will begin by showing a similar inequality for the drift function
$V_{\theta}$, thus not having uniformity, as an intermediary step. The Drift
property is a consequence of Step (2) and (34) shown above. Thus, there exists
$0<\bar{\mu}<1$, $\bar{\delta}>0$ such that for all $z\mathcal{Z}$:

$\Pi_{\theta}V_{\theta}(z)\leq\bar{\mu}V_{\theta}(z)+\bar{\delta}\mathsf{1}_{\mathcal{O}}(z)\;,$
(35)

where $V_{\theta}$ is defined by (6). Using the two functions defined in (7),
we define for $z\in\mathcal{Z}$, the $V$ function independent of $\theta$ as
follows:

$V(z)=V_{1}(z)^{\alpha}V_{2}(z)^{2\alpha}\;,$ (36)

where $0<\alpha<\textrm{min}(\frac{1}{2\beta},\frac{a_{0}}{3})$, $a_{0}$ is
defined in H3 and $\beta$ is defined in (6). Thus for $\theta\in\Theta$,
$z\in\mathcal{Z}$ and $\epsilon>0$:

$\displaystyle\Pi_{\theta}V(z)$
$\displaystyle=\int_{\mathcal{Z}}\Pi_{\theta}(z,y)V_{1}(y)^{\alpha}V_{2}(y)^{2\alpha}\textrm{d}y$
$\displaystyle\;,$
$\displaystyle\overset{(a)}{\leq}\frac{1}{2}\int_{\mathcal{Z}}\Pi_{\theta}(z,y)(\frac{1}{\epsilon^{2}}V_{1}(y)^{2\alpha}+\epsilon^{2}V_{2}(y)^{4\alpha})\textrm{d}y\;,$
$\displaystyle\overset{(b)}{\leq}\frac{1}{2\epsilon^{2}}\int_{\mathcal{Z}}\Pi_{\theta}(z,y)V_{\theta}(y)^{2\alpha}+\frac{\epsilon^{2}}{2}\int_{\mathcal{Z}}\Pi_{\theta}(z,y)V_{2}(y)^{4\alpha}\textrm{d}y\;,$
(37)

where we have used the Young’s inequality in (a) and the definition of
$V_{1}$, see (7), in (b). Then plugging (35) in (37), we have

$\displaystyle\Pi_{\theta}V(z)$
$\displaystyle\leq\frac{1}{2\epsilon^{2}}(\bar{\mu}V_{\theta}(z)^{2\alpha}+\bar{\delta}\mathsf{1}_{\mathcal{O}}(z))+\frac{\epsilon^{2}}{2}\int_{\mathcal{Z}}\Pi_{\theta}(z,y)V_{2}(y)^{4\alpha}\textrm{d}y\;,$
(38)
$\displaystyle\leq\frac{\bar{\mu}}{2\epsilon^{2}}V(z)+\frac{\bar{\delta}}{2\epsilon^{2}}\mathsf{1}_{\mathcal{O}}(z)+\frac{\epsilon^{2}}{2}\int_{\mathcal{Z}}\Pi_{\theta}(z,y)V_{2}(y)^{4\alpha}\textrm{d}y\;,$
(39)
$\displaystyle\leq\frac{\bar{\mu}}{2\epsilon^{2}}V(z)+\frac{\bar{\delta}}{2\epsilon^{2}}\mathsf{1}_{\mathcal{O}}(z)+\frac{\epsilon^{2}}{2}\sup\limits_{\theta\in\Theta,z\in\mathcal{Z}}\int_{\mathcal{Z}}\Pi_{\theta}(z,y)V_{2}(y)^{4\alpha}\textrm{d}y\;,$
(40)
$\displaystyle\leq\frac{\bar{\mu}}{2\epsilon^{2}}V(z)+\frac{\bar{\delta}}{2\epsilon^{2}}\mathsf{1}_{\mathcal{O}}(z)+\frac{\epsilon^{2}}{1+\bar{\mu}}V(z)\;,$
(41)
$\displaystyle\leq\left(\frac{\bar{\mu}}{2\epsilon^{2}}+\frac{\epsilon^{2}}{1+\bar{\mu}}\right)V(z)+\frac{\bar{\delta}}{2\epsilon^{2}}\mathsf{1}_{\mathcal{O}}(z)\;,$
(42)

where we have used (36) and the assumption H3 in the last inequality, ensuring
the existence of such exponent $\alpha$.

Setting
$\epsilon\mathrel{\mathop{:}}=\sqrt{\frac{\bar{\mu}(1+\bar{\mu})}{2}}$,
$\mu\mathrel{\mathop{:}}=\sqrt{\frac{2\bar{\mu}}{1+\bar{\mu}}}$ and
$\delta\mathrel{\mathop{:}}=\frac{\bar{\delta}}{2\epsilon^{2}}$ proves the
uniform ergodicity in (9) and concludes the proof of Theorem 4.2. ∎

###  A.3 Proof of Lemma A.1

######  Lemma.

Define
$\mathcal{P}(z)\mathrel{\mathop{:}}=\\{z-\ell\frac{z}{|z|}-\kappa\nu,\,\textrm{with}\quad\kappa<a-\ell,\nu\in\\{\nu\in\mathbb{R}^{d},\|\nu\|<1\\},|\nu-\frac{z-\ell\frac{z}{|z|}}{|z-\ell\frac{z}{|z|}|}\leq\frac{\epsilon}{2}\\}$
and
$\mathcal{A}_{\theta}(z)\mathrel{\mathop{:}}=\\{y\in\mathcal{Z},\rho_{\theta}(z,y)\geq
1\\}$. Then for $z\in\mathcal{Z}$,
$\mathcal{P}(z)\subset\mathcal{A}_{\theta}(z)$.

###### Proof.

In order to show the inclusion of the set $\mathcal{P}(z)$ in
$\mathcal{A}_{\theta}(z)$ we start by selecting the quantity
$y=z-\ell\frac{z}{|z|}-\kappa\nu$ for $z\in\mathcal{Z}$ and $\kappa<a-\ell$
where $a$ is the radius of the ball used in (25) such that
$y\in\mathcal{P}(z)$. We will now show that $y\in\mathcal{A}_{\theta}(z)$.

By the generalization of Rolle’s theorem applied on the stationary
distribution $\pi_{\theta}$, we guarantee the existence of some $\kappa^{*}$
such that:

$\displaystyle\nabla\pi_{\theta}(z-\ell\frac{z}{|z|}-\kappa^{*}\nu)=\frac{\pi_{\theta}(y)-\pi_{\theta}(z-\ell\frac{z}{|z|})}{y-(z-\ell\frac{z}{|z|})}=-\frac{\pi_{\theta}(y)-\pi_{\theta}(z-\ell\frac{z}{|z|})}{\kappa\nu}\;.$
(43)

Expanding $\nabla\pi_{\theta}(z-\ell\frac{z}{|z|}-\kappa^{*}\nu)$ yields:

$\pi_{\theta}(y)-\pi_{\theta}(z-\ell\frac{z}{|z|})=-\kappa\nu\frac{z-\ell\frac{z}{|z|}-\kappa^{*}\nu}{|z-\ell\frac{z}{|z|}-\kappa^{*}\nu|}|\nabla\pi_{\theta}(z-\ell\frac{z}{|z|}-\kappa^{*}\nu)|\;.$
(44)

Yet, under assumption H1, there exists $\epsilon$ such that

$\frac{\nabla f_{\theta}(z)}{|\nabla
f_{\theta}(z)|}\frac{z}{|z|}\leq-\epsilon\;,$

and for any $y\in\mathcal{P}(z)$ we note that
$\frac{y}{|y|}-\frac{z}{|z|}|\leq\frac{\epsilon}{2}$, by construction of the
set. Thus,

$\frac{\nabla f_{\theta}(y)}{|\nabla f_{\theta}(y)|}\nu=\frac{\nabla
f_{\theta}(y)}{|\nabla
f_{\theta}(y)|}(\nu-\frac{z-\ell\frac{z}{|z|}}{|z-\ell\frac{z}{|z|}|})+\frac{\nabla
f_{\theta}(y)}{|\nabla
f_{\theta}(y)|}(\nu-\frac{z-\ell\frac{z}{|z|}}{|z-\ell\frac{z}{|z|}|}-\frac{y}{|y|})+\frac{\nabla
f_{\theta}(y)}{|\nabla f_{\theta}(y)|}\frac{y}{|y|}\leq 0\;,$ (45)

where $\nu$ is used in the definition of $\mathcal{P}(z)$. Additionally we let
$\frac{\nabla f_{\theta}(y)}{|\nabla f_{\theta}(y)|}\nu$ denote the vector
multiplication between the normalized gradient and $\nu$. Then plugging (45)
into (44) leads to $\pi_{\theta}(y)-\pi_{\theta}(z-\ell\frac{z}{|z|})\geq 0$
and $y\in\mathcal{P}(z)$ implies, using (31), that
$\pi_{\theta}(y)\geq\pi_{\theta}(z-\ell\frac{z}{|z|})\geq\frac{1}{c_{a}}\pi_{\theta}(z)$.
Finally $y\in\mathcal{P}(z)$ implies that $y\in\mathcal{A}_{\theta}(z)$,
concluding the proof of Lemma A.1.

∎