# Improved annealing for sampling from multimodal distributions via landscape
modification

Michael C.H. Choi and Jing Zhang  Yale-NUS College, Singapore and School of
Data Science, The Chinese University of Hong Kong, Shenzhen and Shenzhen
Institute of Artificial Intelligence and Robotics for Society, Guangdong,
518172, P.R. China<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Given a target distribution $\mu\propto e^{-\mathcal{H}}$ to sample from with
Hamiltonian $\mathcal{H}$ on a finite state space, in this paper we propose
and analyze new Metropolis-Hastings sampling algorithms that target an
alternative distribution $\mu^{f}_{1,\alpha,c}\propto
e^{-\mathcal{H}^{f}_{1,\alpha,c}}$, where $\mathcal{H}_{1,\alpha,c}^{f}$ is a
landscape-modified Hamiltonian which we introduce explicitly. The advantage of
the Metropolis dynamics which targets $\mu_{1,\alpha,c}^{f}$ is that it enjoys
reduced critical height described by the threshold parameter $c$, function
$f$, and a penalty parameter $\alpha\geqslant 0$ that controls the state-
dependent effect.

First, we identify a convergence-bias tradeoff phenomenon: a larger value of
$\alpha$ gives faster convergence yet a larger bias from $\mu$. To mitigate
this bias, in the case of fixed $\alpha$ we propose a self-normalized
estimator that corrects for the bias and prove asymptotic convergence results
and Chernoff-type bound of the proposed estimator.

Next, we consider the case of annealing the penalty parameter $\alpha$ so that
the state-dependent effect diminishes to zero in the long run. We prove strong
ergodicity and bounds on the total variation mixing time of the resulting non-
homogeneous chain subject to appropriate assumptions on the decay of $\alpha$.

Finally, we illustrate the proposed algorithms by comparing their mixing times
with the original Metropolis dynamics on a variety of statistical physics
models including the ferromagnetic Ising model on the hypercube or the
complete graph and the $q$-state Potts model on the two-dimensional torus. In
these cases, the mixing times of the classical Glauber dynamics are at least
exponential in the system size as the critical height grows at least linearly
with the size, while the proposed annealing algorithm, with appropriate choice
of $f$, $c$ and annealing schedule on $\alpha$, mixes rapidly with at most
polynomial dependence on the size. The crux of the proof harnesses on the
important observation that the reduced critical height can be bounded
independently of the size that consequently gives rise to rapid mixing.

We emphasize that the techniques developed in this paper are not only limited
to Metropolis-Hastings or Glauber dynamics, and are broadly applicable to
other sampling algorithms.

AMS 2010 subject classifications: 60J27, 60J28, 82B20, 82B80

Keywords: Metropolis-Hastings; simulated annealing; Glauber dynamics;
landscape modification; Potts model; spectral gap; self-normalization;
Chernoff bound

###### Contents

  1. 1 Introduction
    1. 1.1 Notations
  2. 2 Preliminaries
    1. 2.1 Metropolis-Hastings with landscape-modified Hamiltonian
    2. 2.2 Landscape-modified Hamiltonian $\mathcal{H}_{1,\alpha,c}^{f}$ with $f(x)=x$
    3. 2.3 Landscape-modified Hamiltonian $\mathcal{H}_{1,\alpha,c}^{f}$ with $f(x)=x^{2}$
    4. 2.4 Landscape-modified Hamiltonian $\mathcal{H}_{1,\alpha,c}^{f}$ with $f(x)=e^{x}-1$
  3. 3 Spectral gap bounds in terms of the reduced critical height
  4. 4 Constant penalty $\alpha_{t}=\alpha$ for all $t\geqslant 0$
    1. 4.1 Proof of Theorem 4.1
    2. 4.2 Proof of Theorem 4.2
    3. 4.3 Proof of Theorem 4.3
  5. 5 Time-dependent penalty $(\alpha_{t})_{t\geqslant 0}$
    1. 5.1 Proof of Theorem 5.1
    2. 5.2 Proof of Theorem 5.2
    3. 5.3 Proof of Theorem 5.3
  6. 6 Applications
    1. 6.1 Ferromagnetic Ising model on the $n$-dimensional hypercube
      1. 6.1.1 Proof of Theorem 6.1
    2. 6.2 Ferromagnetic Ising model on the complete graph
      1. 6.2.1 Proof of Theorem 6.2
    3. 6.3 $q$-state Potts model on the two-dimensional torus
      1. 6.3.1 Proof of Theorem 6.3

##  1\. Introduction

Sampling lies at the heart of modern computational physics and Bayesian
statistics. Often we would like to generate samples efficiently from a given
target distribution $\mu\propto e^{-\mathcal{H}}$ with Hamiltonian
$\mathcal{H}$ on a finite state space $\mathcal{X}$, without knowing the
normalization constant or the partition function of $\mu$. In this setting,
the classical Metropolis-Hastings (MH) algorithm offers a way to construct
Markov chains that converge to $\mu$ in the long run: a proposal chain is
utilized to generate a possible move, followed by either accepting or
rejecting the proposal move based on a ratio of the target $\mu$. The
resulting MH chain is reversible with respect to $\mu$, and various results
have been developed concerning its long-time convergence. We refer readers to
Roberts and Rosenthal (2004); Brooks et al. (2011); Diaconis (2009) for
further pointers in the literature on the MH algorithm.

While the MH chain is known to mix quickly on unimodal target Johndrow and
Smith (2018), in reality, the target distribution is often multimodal and
high-dimensional, and as much the MH chain can possibly get stuck at a local
mode that consequently leads to torpid mixing. For instance, in Section 6 we
investigate the ferromagnetic Ising model on the $n$-dimensional hypercube.
The total variation mixing time of the classical MH algorithm in this model is
at least exponential in the number of vertices. It is, therefore, crucial to
develop acceleration techniques to tackle the challenges arising from sampling
multimodal $\mu$ and to possibly improve the convergence of Markov chain Monte
Carlo algorithms. A wide variety of tools and algorithms have been suggested
to overcome this bottleneck, for example, the equi-energy sampler Kou et al.
(2006); Andrieu et al. (2008); Ebbers and Löwe (2020), simulated tempering or
replica exchange methods Bhatnagar and Randall (2016); Madras and Zheng
(2003); Woodard et al. (2009), and more recently non-reversible Markovian
Monte Carlo algorithms Bierkens (2016); Hwang et al. (1993, 2005), to name but
a few.

The core of this paper centers around a technique that we call landscape
modification. The high-level idea is to run algorithms on an alternative
landscape with reduced critical height that hopefully accelerate the
convergence of the Markov chain. The starting point of this work originates
from Fang et al. (1997), in which the authors therein analyze a variant of the
overdamped Langevin diffusion with state-dependent diffusion coefficient for
simulated annealing. In Choi (2020), it is realized that the state-dependent
effect can be effectively cast to kinetic simulated annealing by simply
altering the target Hamiltonian function to optimize, and in Choi (2021a) the
same idea is applied in the discrete setting.

While the focus of the previous papers Choi (2020, 2021a) are on global
optimization, in this paper we propose and analyze new MH sampling algorithms
that target an alternative distribution $\mu^{f}_{1,\alpha,c}\propto
e^{-\mathcal{H}^{f}_{1,\alpha,c}}$, where $\mathcal{H}_{1,\alpha,c}^{f}$ is a
landscape-modified Hamiltonian which we shall define explicitly in Section 2.
One significant benefit of the Metropolis dynamics which targets
$\mu_{1,\alpha,c}^{f}$ is that it enjoys reduced critical height which depends
on the threshold parameter $c$, function $f$ and a penalty parameter
$\alpha\geqslant 0$ that controls the state-dependent effect.

We summarize the major contributions of this paper as follows:

  * •

Propose and analyze a self-normalized estimator in the case of fixed penalty:
To mitigate the bias of sampling from $\mu^{f}_{1,\alpha,c}$ instead of $\mu$,
we propose a self-normalized estimator that correct for the bias and prove
both asymptotic and non-asymptotic convergence results of the proposed
estimator.

  * •

Propose and analyze an annealing algorithm on the penalty parameter: The
second algorithm that we propose involves annealing the penalty parameter
$\alpha$. We prove strong ergodicity and bounds on the total variation mixing
time of the resulting non-homogeneous Markov chain subject to appropriate
assumptions on the decay of $\alpha$.

  * •

Rapid mixing of the proposed annealing algorithm on the ferromagnetic Ising
model on complete graph or hypercube and the $q$-state Potts model on the two-
dimensional torus: We illustrate the benefits of utilizing landscape
modification in the ferromagnetic Ising model on the hypercube or the complete
graph and the $q$-state Potts model on the two-dimensional torus. In these
cases, the mixing times of the classical Glauber dynamics are at least
exponential in the system size as the critical height grows at least linearly
with the number of vertices, while the proposed annealing algorithm, with
appropriate choice of $f$, $c$ and annealing schedule on $\alpha$, mixes
rapidly with at most polynomial dependence on the system size. The crux of the
proof relies on the observation that the reduced critical height can be
bounded independently of the size that consequently leads to rapid mixing.

We stress that the techniques developed in this paper are not only limited to
MH, and are in fact broadly applicable to other sampling algorithms by simply
altering the target distribution from $\mu$ to $\mu^{f}_{1,\alpha,c}$,
followed by either annealing the penalty parameter or to correct the bias via
self-normalization.

The rest of this paper is organized as follows. In Section 1.1, we fix some
commonly-used notations throughout the paper. In Section 2, we first recall
the classical MH and the notion of critical or communication height, followed
by defining the landscape-modified Hamiltonian
$\mathcal{H}_{1,\alpha_{t},c}^{f}$ and the MH chain on the modified landscape.
We conclude the Section by calculating explicitly
$\mathcal{H}_{1,\alpha_{t},c}^{f}$ when $f$ is a linear, quadratic or
exponential function. We then proceed to bound the spectral gap with reduced
critical height of the proposed chain in Section 3. In Section 4, we consider
the case of constant penalty and introduce a self-normalized estimator. We
prove related limit theorems as well as non-asymptotic Chernoff bounds of the
estimator. In Section 5, we investigate the situation with time-dependent
penalty and offer bounds on the total variation mixing time. Finally, we
illustrate the model-independent results obtained in previous Sections in the
ferromagnetic Ising model on hypercube or complete graph and the $q$-state
Potts model on the two-dimensional torus in Section 6.

###  1.1. Notations

Throughout this paper, we adopt the following notations. For
$x,y\in\mathbb{R}$, we write $x_{+}=\max\\{x,0\\}$ to denote the non-negative
part of $x$, and $x\wedge y=\min\\{x,y\\}$. We also write $\mathds{1}_{A}$ to
be the indicator function of the set $A$.

In order to compare the efficiency of algorithms, we now introduce a few
common asymptotic notations. For two functions
$g_{1},g_{2}:\mathbb{R}\to\mathbb{R}$, we say that $g_{1}=\mathcal{O}(g_{2})$
if there exists a constant $C>0$ such that for sufficiently large $x$, we have
$|g_{1}(x)|\leqslant Cg_{2}(x)$. On the other hand, we use $\Omega$ to denote
the tight lower bound, that is, we say that $g_{1}=\Omega(g_{2})$ if there
exists a constant $C>0$ such that for sufficiently large $x$, we have
$|g_{1}(x)|\geqslant Cg_{2}(x)$. The $\ell^{\infty}$-norm of a bounded $g$ is
defined as $\|g\|_{\infty}=\sup_{x}\\{|g(x)|\\}$.

We say that a sequence of random variables $(X_{n})_{n\in\mathbb{N}}$
converges almost surely to the random variable $X$, denoted as
$X_{n}\rightarrow^{a.s.}X$, if
$\mathbb{P}\left(\lim_{n\rightarrow\infty}X_{n}=X\right)=1$. Meanwhile, we say
that $(X_{n})_{n\in\mathbb{N}}$ converges in distribution to $X$, shown as
$X_{n}\rightarrow^{d}X$, if $\lim_{n\to\infty}\mathbb{P}(X_{n}\leqslant
x)=\mathbb{P}(X\leqslant x)$ for any $x\in\mathbb{R}$. Throughout the paper,
we denote $\mathbb{P}_{x}$ to be the probability measure when the Markov chain
is initialized at $x$, and $\mathbb{E}_{\mu}(X)$ denotes the expectation of
$X$ with respect to $\mu$. In Section 5, we introduce the time-dependent
penalty parameter $(\alpha_{t})_{t\geqslant 0}$ and write
$\alpha^{\prime}_{t}=d\alpha_{t}/dt$ as the derivative of $\alpha_{t}$ with
respect to $t$.

##  2\. Preliminaries

Let $\mathcal{H}:\mathcal{X}\rightarrow\mathbb{R}$ be a given Hamiltonian
function on a finite state space $\mathcal{X}$. The primary objective in this
paper is to sample from the target distribution $\mu\propto e^{-\mathcal{H}}$.
Among various Markov chain Monte Carlo algorithms to sample from $\mu$, we
shall be interested in the classical MH chain and its variants in the
continuous-time setting. The MH algorithm offers a way to construct Markov
chain that is reversible with respect to the target $\mu$, and now we recall
its definition:

######  Definition 2.1 (Classical MH with target distribution $\mu$).

Given a target distribution $\mu$ on $\mathcal{X}$ and a proposal continuous-
time Markov chain $X^{0}=(X^{0}(t))_{t\geqslant 0}$ with symmetric generator
$Q$, the MH chain has generator given by
$M^{0}=M^{0}(Q,\mu)=(M^{0}(x,y))_{x,y\in\mathcal{X}}$, where

$M^{0}(x,y)=M^{0}(Q,\mu)(x,y):=\begin{cases}\min\left\\{Q(x,y),\dfrac{\mu(y)}{\mu(x)}Q(y,x)\right\\}=Q(x,y)e^{-(\mathcal{H}(y)-\mathcal{H}(x))_{+}},&\mbox{if
}x\neq y,\\\ -\sum_{y:x\neq y}M^{0}(x,y),&\mbox{if }x=y.\end{cases}$

We shall also denote the transition semigroup to be $((P^{0})^{t})_{t\geqslant
0}$.

We remark that the superscript of $0$ that appears in Definition 2.1 will be
explained in Section 2.1 below: it turns out that in the special case when
$f=0$, the modified MH dynamics that we propose reduces to the classical MH.

We equip the Hilbert space $\ell^{2}(\mu)$ with the usual inner product
weighted by the distribution $\mu$: for $g_{1},g_{2}\in\ell^{2}(\mu)$,

$\langle
g_{1},g_{2}\rangle_{\mu}:=\sum_{x\in\mathcal{X}}g_{1}(x)g_{2}(x)\mu(x),~{}~{}~{}~{}\|g_{1}\|_{\mu}=\langle
g_{1},g_{1}\rangle_{\mu},~{}~{}\text{ and }~{}~{}\mu(g_{1})=\langle
g_{1},1\rangle_{\mu}.$

For any infinitesimal generator $M$ of an ergodic continuous-time Markov chain
$X$ with stationary distribution $\mu$, for instance $M^{0}$, we let
$\\{0=\lambda_{1},\lambda_{2},...,\lambda_{|\mathcal{X}|}\\}$ be the
eigenvalues, counted with their algebraic multiplicities and arranged in non-
descending order, of $-M$. In particular, $\lambda_{2}(-M)$ is the spectral
gap of $X$, that is,

(2.1)
$\displaystyle\lambda_{2}(-M):=\inf_{g\in\ell^{2}(\mu):\mu(g)=0}\dfrac{\langle-
Mg,g\rangle_{\mu}}{\langle g,g\rangle_{\mu}}.$

The spectral gap controls the rate at which a Markov chain converges to its
stationary distribution. In this paper, we shall be interested in bounding the
total variation mixing time. To this end, let us recall that for any
probability measure $\nu_{1},\nu_{2}$ with support on $\mathcal{X}$, the total
variation distance between $\nu_{1}$ and $\nu_{2}$ is

$||\nu_{1}-\nu_{2}||_{TV}:=\sup_{A\subset\mathcal{X}}|\nu_{1}(A)-\nu_{2}(A)|=\dfrac{1}{2}\sum_{x\in\mathcal{X}}|\nu_{1}(x)-\nu_{2}(x)|.$

The total variation mixing time of $X^{0}$ is then defined as:

(2.2) $t_{mix}^{0}(\varepsilon):=\inf\bigg{\\{}t\geqslant
0;~{}\max_{x\in\mathcal{X}}\left\|\left(P^{0}\right)^{t}(x,\cdot)-\mu\right\|_{TV}\leqslant\varepsilon\bigg{\\}}.$

To bound the spectral gap, we recall the concept of critical height or
communication height as in the metastability literature Holley and Stroock
(1987); Bovier and den Hollander (2015). For a given irreducible proposal
chain with generator $Q$, a path from $x$ to $y$ is any sequence of points
starting from $x_{0}=x,x_{1},x_{2},\ldots,x_{n}=y$ such that
$Q(x_{i-1},x_{i})>0$ for $i=1,2,\ldots,n$. For any $x\neq y$, such path exists
as $Q$ is irreducible. We write $\chi^{x,y}$ to be the set of paths from $x$
to $y$, and elements of $\chi^{x,y}$ are denoted by
$\omega=(\omega_{i})_{i=0}^{n}$. The highest value of the original Hamiltonian
function $\mathcal{H}$ along a path $\omega\in\chi^{x,y}$, known as the
elevation, is defined to be

$\mathrm{Elev}(\omega):=\max_{\omega_{i}\in\omega}\\{\mathcal{H}(\omega_{i})\\},$

and the lowest possible highest elevation along path(s) from $x$ to $y$ is

$H(x,y):=\min_{\omega\in\chi^{x,y}}\\{\mathrm{Elev}(\omega)\\}.$

For $X^{0}$, the associated critical height is defined to be

(2.3) $\displaystyle m^{0}$
$\displaystyle:=\max_{x,y\in\mathcal{X}}\\{H(x,y)-\mathcal{H}(x)-\mathcal{H}(y)+\min_{x}\mathcal{H}(x)\\}.$

###  2.1. Metropolis-Hastings with landscape-modified Hamiltonian

Instead of directly targeting $\mathcal{H}$ in the Gibbs distribution
$\mu(x)\propto e^{-\mathcal{H}(x)}$ in MH, we shall consider sampling from an
alternative Gibbs distribution

$\displaystyle\mu_{1,\alpha,c}^{f}(x)=\dfrac{e^{-\mathcal{H}_{1,\alpha,c}^{f}(x)}}{Z_{1,\alpha,c}^{f}}\propto
e^{-\mathcal{H}_{1,\alpha,c}^{f}(x)}$

in the MH chain, where

(2.4) $\displaystyle
Z_{1,\alpha,c}^{f}:=\sum_{x\in\mathcal{X}}e^{-\mathcal{H}_{1,\alpha,c}^{f}(x)}$

is the normalization constant. The modified or transformed Hamiltonian
function $\mathcal{H}_{1,\alpha,c}^{f}$ depends on three parameters: the
function $f$, the threshold parameter $c$ and the penalty parameter $\alpha$.
These parameters are assumed to be chosen in the following manner:

######  Assumption 2.1.

  1. (1)

The function $f:\mathbb{R}^{+}\to\mathbb{R}^{+}$ is differentiable, non-
negative and non-decreasing. Furthermore, $f$ satisfies

$f(0)=0.$

  2. (2)

The threshold parameter $c$ is chosen so that $\mathcal{H}_{\max}\geqslant
c\geqslant\mathcal{H}_{\min}$, where
$\mathcal{H}_{\max}:=\max_{x\in\mathcal{X}}\mathcal{H}(x)$ and
$\mathcal{H}_{\min}:=\min_{x\in\mathcal{X}}\mathcal{H}(x)$ are respectively
the global maximum and minimum value of $\mathcal{H}$.

  3. (3)

The penalty parameter $\alpha$ is non-negative, i.e. $\alpha\geqslant 0$.

In the rest of this paper, we suppose that Assumption 2.1 holds. We now state
the definition of $\mathcal{H}_{1,\alpha,c}^{f}$:

######  Definition 2.2 (Landscape-modified Hamiltonian
$\mathcal{H}_{\epsilon,\alpha,c}^{f}$).

Suppose Assumption 2.1 holds. At temperature $\epsilon>0$, the landscape-
modified function $\mathcal{H}_{\epsilon,\alpha,c}^{f}$ is

$\displaystyle\mathcal{H}^{f}_{\epsilon,\alpha,c}(x)$
$\displaystyle:=\int_{\mathcal{H}_{\min}}^{\mathcal{H}(x)}\dfrac{1}{\alpha
f((u-c)_{+})+\epsilon}\,du.$

In particular, if we take $\epsilon=1$, then

$\displaystyle\mathcal{H}_{1,\alpha,c}^{f}(x)=\int_{\mathcal{H}_{\min}}^{\mathcal{H}(x)}\dfrac{1}{\alpha
f((u-c)_{+})+1}\,du=\begin{cases}\mathcal{H}(x)-\mathcal{H}_{\min},\quad&\text{if}\,x\in\\{\mathcal{H}(x)<c\\},\\\
c-\mathcal{H}_{\min}+\int_{c}^{\mathcal{H}(x)}\dfrac{1}{\alpha
f(u-c)+1}\,du,\quad&\text{if}\,x\in\\{\mathcal{H}(x)\geqslant c\\}.\\\
\end{cases}$

If we further specialize into $f=0$ or $\alpha=0$, then we retrieve the
original target distribution as $\mu=\mu^{0}_{1,\alpha,c}=\mu^{f}_{1,0,c}$.

We shall see in Section 3 that one important benefit of running MH on the
modified landscape $\mathcal{H}_{1,\alpha,c}^{f}$ is that its associated
critical height is less than or equal to the original critical height $m^{0}$.

To understand $\mathcal{H}_{1,\alpha,c}^{f}$ intuitively, we recall that this
transformation is first introduced in Choi (2021a) for simulated annealing in
the discrete setting for global optimization, and is inspired by the
overdamped Langevin diffusion with state-dependent diffusion coefficient in
Fang et al. (1997). The parameter $c$ can be interpreted as a threshold
parameter: on $\\{\mathcal{H}(x)<c\\}$, $\mathcal{H}_{1,\alpha,c}^{f}$ shares
the same scale as the original $\mathcal{H}$; while on
$\\{\mathcal{H}(x)\geqslant c\\}$, the state-dependent effect is activated
since the integral

$\int_{c}^{\mathcal{H}(x)}\dfrac{1}{\alpha
f(u-c)+1}\,du\leqslant\mathcal{H}(x)-c.$

Thus, the function $f$ can be understood as determining the magnitude of the
state-dependent effect: if we choose a larger $f$, the above integral on the
left hand side and hence $\mathcal{H}_{1,\alpha,c}^{f}$ are smaller. On the
other hand, this effect is also controlled and penalized by $\alpha$. A larger
$\alpha$ gives smaller $\mathcal{H}_{1,\alpha,c}^{f}$ yet a larger bias
between $\mu_{1,\alpha,c}^{f}$ and $\mu$. The name penalty parameter is
largely motivated by the work in Karagulyan and Dalalyan (2020), where the
authors therein consider a penalized Langevin diffusion in which the negative
of the gradient of the potential function is perturbed with a penalty that
vanishes as time increases.

We now define formally the MH chain that targets $\mu_{1,\alpha,c}^{f}$:

######  Definition 2.3 (MH with target distribution $\mu_{1,\alpha,c}^{f}$).

Suppose Assumption 2.1 holds. Given the target distribution
$\mu_{1,\alpha,c}^{f}$ on $\mathcal{X}$ and a proposal continuous-time Markov
chain $X_{1,\alpha,c}^{f}=(X_{1,\alpha,c}^{f}(t))_{t\geqslant 0}$ with
symmetric generator $Q$, the MH chain has generator given by
$M_{1,\alpha,c}^{f}=M_{1,\alpha,c}^{f}(Q,\mu)=(M_{1,\alpha,c}^{f}(x,y))_{x,y\in\mathcal{X}}$,
where

$M_{1,\alpha,c}^{f}(x,y):=\begin{cases}\min\left\\{Q(x,y),\dfrac{\mu_{1,\alpha,c}^{f}(y)}{\mu_{1,\alpha,c}^{f}(x)}Q(y,x)\right\\}=Q(x,y)e^{-(\mathcal{H}_{1,\alpha,c}^{f}(y)-\mathcal{H}_{1,\alpha,c}^{f}(x))_{+}},&\mbox{if
}x\neq y,\\\ -\sum_{y:x\neq y}M_{1,\alpha,c}^{f}(x,y),&\mbox{if
}x=y.\end{cases}$

We shall also denote the transition semigroup to be
$((P_{1,\alpha,c}^{f})^{t})_{t\geqslant 0}$, and its spectral gap to be

$\lambda_{2}(-M_{1,\alpha,c}^{f}):=\inf_{g\in\ell^{2}(\mu_{1,\alpha,c}^{f}):\mu_{1,\alpha,c}^{f}(g)=0}\dfrac{\langle-
M_{1,\alpha,c}^{f}g,g\rangle_{\mu_{1,\alpha,c}^{f}}}{\langle
g,g\rangle_{\mu_{1,\alpha,c}^{f}}}.$

Note that if $f=0$ or $\alpha=0$, the above dynamics reduce to the classical
MH with target distribution $\mu$ since

$M^{0}=M^{0}_{1,\alpha,c}=M^{f}_{1,0,c}.$

######  Remark 2.1.

In view of Definition 2.3, the knowledge of $\mathcal{H}_{\min}$ is not
required in simulating $X_{1,\alpha,c}^{f}$ since the difference, for
$x,y\in\mathcal{X}$,

$\mathcal{H}_{1,\alpha,c}^{f}(y)-\mathcal{H}_{1,\alpha,c}^{f}(x)=\int_{\mathcal{H}(x)}^{\mathcal{H}(y)}\dfrac{1}{\alpha
f((u-c)_{+})+1}\,du$

does not depend on $\mathcal{H}_{\min}$.

In the following subsections, we compute explicitly the landscape-modified
Hamiltonian $\mathcal{H}_{1,\alpha,c}^{f}$ with linear, quadratic or
exponential $f$.

###  2.2. Landscape-modified Hamiltonian $\mathcal{H}_{1,\alpha,c}^{f}$ with
$f(x)=x$

In this subsection, we compute explicitly $\mathcal{H}_{1,\alpha,c}^{f}$ when
we specialize into $f(x)=x$, which gives

$\displaystyle\mathcal{H}_{1,\alpha,c}^{f}(x)=\begin{cases}\mathcal{H}(x)-\mathcal{H}_{\min},\quad&\text{if}\,x\in\\{\mathcal{H}(x)<c\\},\\\
c-\mathcal{H}_{\min}+\dfrac{1}{\alpha}\ln(1+\alpha(\mathcal{H}(x)-c)),\quad&\text{if}\,x\in\\{\mathcal{H}(x)\geqslant
c\\}.\\\ \end{cases}$

Using the fact that

$\lim_{\alpha\to
0}\dfrac{1}{\alpha}\ln(1+\alpha(\mathcal{H}(x)-c))=\mathcal{H}(x)-c,$

we thus obtain, for $x\in\mathcal{X}$,

$\lim_{\alpha\to
0}\mathcal{H}_{1,\alpha,c}^{f}(x)=\mathcal{H}(x)-\mathcal{H}_{\min}.$

In other words, we retrieve the original target distribution
$\mu\propto{e^{-\mathcal{H}}}\propto e^{-(\mathcal{H}-\mathcal{H}_{\min})}$
when we take $\alpha\to 0$, as expected.

###  2.3. Landscape-modified Hamiltonian $\mathcal{H}_{1,\alpha,c}^{f}$ with
$f(x)=x^{2}$

In this subsection, we write down explicitly $\mathcal{H}_{1,\alpha,c}^{f}$
when we specialize into $f(x)=x^{2}$, which yields

(2.5)
$\displaystyle\mathcal{H}_{1,\alpha,c}^{f}(x)=\begin{cases}\mathcal{H}(x)-\mathcal{H}_{\min},\quad&\text{if}\,x\in\\{\mathcal{H}(x)<c\\},\\\
c-\mathcal{H}_{\min}+\dfrac{1}{\sqrt{\alpha}}\arctan(\sqrt{\alpha}(\mathcal{H}(x)-c)),\quad&\text{if}\,x\in\\{\mathcal{H}(x)\geqslant
c\\}.\\\ \end{cases}$

As a result, we observe that

$\max_{x}\mathcal{H}_{1,\alpha,c}^{f}(x)\leqslant
c-\mathcal{H}_{\min}+\dfrac{1}{\sqrt{\alpha}}\dfrac{\pi}{2}.$

Consequently, if we take $c=\mathcal{H}_{\min}+\delta$, where $\delta>0$, the
maximum of $\mathcal{H}_{1,\alpha,c}^{f}$ can be bounded above independently
of $\mathcal{H}_{\max}$ or $\mathcal{H}_{\min}$. This important observation
will be utilized in Section 6 on a few statistical physics models to
illustrate rapid mixing of the proposed annealing algorithm in Section 5.

###  2.4. Landscape-modified Hamiltonian $\mathcal{H}_{1,\alpha,c}^{f}$ with
$f(x)=e^{x}-1$

In this subsection, we calculate explicitly $\mathcal{H}_{1,\alpha,c}^{f}$
when we specialize into $f(x)=e^{x}-1$, which yields

(2.6)
$\displaystyle\mathcal{H}_{1,\alpha,c}^{f}(x)=\begin{cases}\mathcal{H}(x)-\mathcal{H}_{\min},\quad&\text{if}\,x\in\\{\mathcal{H}(x)<c\\},\\\
c-\mathcal{H}_{\min}+\dfrac{\log\left(\alpha\left(e^{\mathcal{H}(x)}-1\right)+1\right)-\log\left(\alpha\left(e^{c}-1\right)+1\right)-(\mathcal{H}(x)-c)}{\alpha-1},\quad&\text{if}\,x\in\\{\mathcal{H}(x)\geqslant
c\\}.\\\ \end{cases}$

Also, we know that

$\lim_{\alpha\to 0}\dfrac{\log(\alpha(e^{\mathcal{H}(x)}-1)+1)}{\alpha-1}=0,$

thus

$\lim_{\alpha\to
0}\dfrac{\log\left(\alpha\left(e^{\mathcal{H}(x)}-1\right)+1\right)-\log\left(\alpha\left(e^{c}-1\right)+1\right)-(\mathcal{H}(x)-c)}{\alpha-1}=H(x)-c,$

which implies that, for $x\in\mathcal{X}$,

$\lim_{\alpha\to
0}\mathcal{H}_{1,\alpha,c}^{f}(x)=\mathcal{H}(x)-\mathcal{H}_{\min}.$

In other words, we recover the original target distribution
$\mu\propto{e^{-\mathcal{H}}}\propto e^{-(\mathcal{H}-\mathcal{H}_{\min})}$
when we take $\alpha\to 0$, as expected.

##  3\. Spectral gap bounds in terms of the reduced critical height

In this Section, we bound the spectral gap of $X^{f}_{1,\alpha,c}$, the MH
dynamics with target distribution $\pi^{f}_{1,\alpha,c}$ as introduced in
Definition 2.3, in terms of its associated critical height. It turns out that
this critical height is reduced and is less than or equal to the original
critical height $m^{0}$, which demonstrates one advantage of running MH
algorithms on the modified landscape described by
$\mathcal{H}^{f}_{1,\alpha,c}$.

Recall that we write $\chi^{x,y}$ to be the set of paths from $x$ to $y$ as
specified by the proposal generator. The highest value of the landscape-
modified Hamiltonian function $\mathcal{H}^{f}_{1,\alpha,c}$ along a path
$\omega\in\chi^{x,y}$ is defined to be

$\mathrm{Elev}^{f}_{1,\alpha,c}(\omega):=\max_{\omega_{i}\in\omega}\\{\mathcal{H}^{f}_{1,\alpha,c}(\omega_{i})\\},$

and the lowest possible highest elevation along path(s) from $x$ to $y$ is

$H^{f}_{1,\alpha,c}(x,y):=\min_{\omega\in\chi^{x,y}}\\{\mathrm{Elev}^{f}_{1,\alpha,c}(\omega)\\}.$

With the above notations in mind, we define the concept of reduced critical
height:

######  Definition 3.1 (Reduced critical height).

Given the same proposal chain with generator $Q$ of both $X^{0}$ and
$X^{f}_{1,\alpha,c}$, Then the critical height of $X^{f}_{1,\alpha,c}$ is
defined to be

$m_{1,\alpha,c}^{f}:=\max_{x,y\in\mathcal{X}}\\{H^{f}_{1,\alpha,c}(x,y)-\mathcal{H}^{f}_{1,\alpha,c}(x)-\mathcal{H}^{f}_{1,\alpha,c}(y)\\}.$

######  Remark 3.1.

Suppose $(x^{*},y^{*})$ attains the critical height $m^{f}_{1,\alpha,c}$ of
$X^{f}_{1,\alpha,c}$ in the sense that it satisfies
$m^{f}_{1,\alpha,c}=\mathcal{H}^{f}_{1,\alpha,c}(x^{*})-\mathcal{H}^{f}_{1,\alpha,c}(y^{*})$.
Since $\alpha\geqslant 0$ and $f\geqslant 0$ we thus have

$\displaystyle
m_{1,\alpha,c}^{f}=\int_{\mathcal{H}\left(y^{*}\right)}^{\mathcal{H}\left(x^{*}\right)}\frac{1}{\alpha
f\left((u-c)_{+}\right)+1}\,du\leqslant\mathcal{H}(x^{*})-\mathcal{H}(y^{*})\leqslant
m^{0}.$

Thus, the reduced critical height $m_{1,\alpha,c}^{f}$ is bounded above by the
original critical height $m^{0}$. This illustrates one advantage of using the
chain $X^{f}_{1,\alpha,c}$.

Using (Holley and Stroock, 1988, Lemma $2.7$) with the Hamiltonian function
therein be replaced by $\mathcal{H}^{f}_{1,\alpha,c}$ and at temperature $1$,
we can readily bound the spectral gap $\lambda_{2}(-M^{f}_{1,\alpha,c})$ in
terms of $m^{f}_{1,\alpha,c}$:

######  Lemma 3.1.

$\tilde{C}e^{-m_{1,\alpha,c}^{f}}\geqslant\lambda_{2}(-M_{1,\alpha_{t},c}^{f})\geqslant
Ce^{-m_{1,\alpha,c}^{f}}\geqslant Ce^{-m^{0}},$

where $C$ and $\tilde{C}$ only depend on $\mathcal{X}$ and
$m_{1,\alpha,c}^{f}$ is the reduced critical height as introduced in
Definition 3.1.

Consequently, thanks to Lemma 3.1 we obtain

######  Corollary 3.1.

Suppose $\alpha=(\alpha_{t})_{t\geqslant 0}$ is time-dependent. We then have

(3.1) $\lambda_{2}(-M_{1,\alpha_{t},c}^{f})\geqslant Ce^{-m^{0}}$

where $Ce^{-m^{0}}$ is independent of $t$, which implies that

$\int_{0}^{\infty}\lambda_{2}(-M_{1,\alpha_{t},c}^{f})dt=\infty.$

We remark that Corollary 3.1 is useful in establishing strong ergodicity in
Section 5 when we consider the time-dependent penalty.

##  4\. Constant penalty $\alpha_{t}=\alpha$ for all $t\geqslant 0$

In this Section, we focus in the constant penalty case, i.e. the penalty
parameter takes on a fixed value with $\alpha_{t}=\alpha>0$ for all
$t\geqslant 0$. In this setting, we first identify a phenomenon that we call
the (total-variation) convergence-bias tradeoff: using triangle inequality, we
note that, for any $x\in\mathcal{X}$,

(4.1)
$\displaystyle\left\|\left(P_{1,\alpha,c}^{f}\right)^{t}(x,\cdot)-\mu\right\|_{TV}\leqslant\underbrace{\left\|\left(P_{1,\alpha,c}^{f}\right)^{t}(x,\cdot)-\mu_{1,\alpha,c}^{f}\right\|_{TV}}_{\text{convergence}}+\underbrace{\|\mu_{1,\alpha,c}^{f}-\mu\|_{TV}}_{\text{bias}}.$

On one hand, we would like to choose $\alpha$ to be large so that the reduced
critical height $m^{f}_{1,\alpha,c}$ is small, leading to fast convergence of
$\left(P_{1,\alpha,c}^{f}\right)^{t}(x,\cdot)$ towards $\mu_{1,\alpha,c}^{f}$.
On the other hand, we would like to choose $\alpha$ to be small so that the
bias between $\mu_{1,\alpha,c}^{f}$ and $\mu$ is small when measured in terms
of the total variation distance. We quantify the above statements into the
following theorem:

######  Theorem 4.1.

Suppose that Assumption 2.1 holds.

  1. (1)

(Large $\alpha$ leads to fast convergence): For $x\in\mathcal{X}$,

$\displaystyle\dfrac{1}{2}e^{-\lambda_{2}(-M_{1,\alpha,c}^{f})t}\leqslant\left\|\left(P_{1,\alpha,c}^{f}\right)^{t}(x,\cdot)-\mu_{1,\alpha,c}^{f}\right\|_{TV}\leqslant\dfrac{e^{\mathcal{H}_{\max}-\mathcal{H}_{\min}}|\mathcal{X}|}{2}e^{-\lambda_{2}(-M_{1,\alpha,c}^{f})t}.$

  2. (2)

(Small $\alpha$ gives small total variation bias):

$\|\mu_{1,\alpha,c}^{f}-\mu\|_{TV}=\mathcal{O}(\alpha).$

The proof of Theorem 4.1 is given in Section 4.1.

To mitigate the bias induced by running the Markov chain $X^{f}_{1,\alpha,c}$
to approximately sample from $\mu_{1,\alpha,c}^{f}$ instead of $\mu$, in this
Section, we consider a self-normalized estimator to correct the bias. Define
the weight function to be

$w(x)=w^{f}_{1,\alpha,c}(x):=\dfrac{e^{-(\mathcal{H}(x)-\mathcal{H}_{\min})}}{e^{-\mathcal{H}_{1,\alpha,c}^{f}(x)}},$

and for any $g\in\ell^{2}(\mu)$, the self-normalized estimator by

$\dfrac{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds}{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds}.$

We remark that the weight function $w(x)$ can be computed readily since

$\displaystyle w(x)$
$\displaystyle=\dfrac{e^{-(\mathcal{H}(x)-\mathcal{H}_{\min})}}{e^{-\mathcal{H}_{1,\alpha,c}^{f}(x)}}$
$\displaystyle=\dfrac{e^{-(\mathcal{H}(x)-\mathcal{H}_{\min})}}{\exp\bigg{\\{}-(c\wedge\mathcal{H}(x)-\mathcal{H}_{\min})-\int_{c\wedge\mathcal{H}(x)}^{\mathcal{H}(x)}\dfrac{1}{\alpha
f((u-c)_{+})+1}\,du\bigg{\\}}}$
$\displaystyle=\dfrac{e^{-\mathcal{H}(x)}}{\exp\bigg{\\{}-c\wedge\mathcal{H}(x)-\int_{c\wedge\mathcal{H}(x)}^{\mathcal{H}(x)}\dfrac{1}{\alpha
f((u-c)_{+})+1}\,du\bigg{\\}}}$

does not depend on $\mathcal{H}_{\min}$. In the special case where $f=0$ or
$\alpha=0$, we note that $w=1$. Recall that $Z_{1,\alpha,c}^{f}$ is the
normalization constant of $\mu_{1,\alpha,c}^{f}$ (2.4), and we adopt the
notation

$Z^{0}:=Z^{0}_{1,\alpha,c}=\sum_{x\in\mathcal{X}}e^{-(\mathcal{H}(x)-\mathcal{H}_{\min})}$

to be the normalization constant when $f=0$. We then have the following
asymptotic convergence results for the proposed self-normalized estimator:

######  Theorem 4.2.

For $g\in\ell^{2}(\mu)$, as $t\to\infty$ we have

  1. (1)

(Strong consistency of the self-normalized estimator)

$\dfrac{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds}{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds}\rightarrow^{a.s.}\mu(g).$

  2. (2)

(Central limit theorem of the self-normalized estimator)

$\sqrt{t}\left(\dfrac{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds}{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds}-\mu(g)\right)\rightarrow^{d}\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}N(0,\mathrm{Var}_{\mu_{1,\alpha,c}^{f}}(wg)).$

where the limiting variance is defined to be

$\mathrm{Var}_{\mu_{1,\alpha,c}^{f}}(wg):=\lim_{t\to\infty}\mathbb{E}_{\mu}\left(\dfrac{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds}{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds}-\mu(g)\right)^{2}.$

In particular, in the special case when $f=0$ or $\alpha=0$ so that $w=1$, the
above results reduce to the classical Markov chain law of large number and
central limit theorem for the MH chain respectively.

The proof of Theorem 4.2 is delayed to Section 4.2.

Our next result concerns the finite-time non-asymptotic deviation of the self-
normalized estimator from the target mean $\mu(g)$. To this end, for any
$\eta,\varepsilon>0$, $g\in\ell^{2}(\mu)$, we define the first time such that
the deviation probability is less than or equal to $\varepsilon$ to be:

(4.2) $\displaystyle t_{conc}(\eta,\varepsilon):=\inf\bigg{\\{}t\geqslant
0;~{}\sup_{g\in\ell^{2}(\mu);~{}\left\lVert g\right\rVert_{\infty}\leqslant
a}\sup_{x\in\mathcal{X}}\mathbb{P}_{x}\left(\dfrac{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds}{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds}-\mu(g)\geqslant\eta\right)\leqslant\varepsilon\bigg{\\}}.$

######  Theorem 4.3 (Chernoff’s bound of the self-normalized estimator).

Given $\eta,\varepsilon>0$, $g\in\ell^{2}(\mu)$ such that $\left\lVert
g\right\rVert_{\infty}\leqslant a$, if we have

$t\geqslant\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}\dfrac{32a^{2}(1+\eta)^{2}}{\lambda_{2}(-M_{1,\alpha,c}^{f})\eta^{2}}\ln\left(\dfrac{1}{\varepsilon\min_{x\in\mathcal{X}}\mu_{1,\alpha,c}^{f}(x)}\right),$

then for any $x\in\mathcal{X}$,

$\mathbb{P}_{x}\left(\dfrac{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds}{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds}-\mu(g)\geqslant\eta\right)\leqslant\varepsilon.$

In particular,

$\displaystyle
t_{conc}(\eta,\varepsilon)=\mathcal{O}\left(|\mathcal{X}|\dfrac{a^{2}(1+\eta)^{2}}{\lambda_{2}(-M_{1,\alpha,c}^{f})\eta^{2}}\ln\left(\dfrac{|\mathcal{X}|\exp\\{\max_{x\in\mathcal{X}}\mathcal{H}_{1,\alpha,c}^{f}(x)\\}}{\varepsilon}\right)\right).$

The proof of Theorem 4.3 can be found in Section 4.3, which utilizes the
results in Lezaud (2001).

###  4.1. Proof of Theorem 4.1

First, we prove item (1). For the lower bound, let $\phi$ be an eigenfunction
of $-M_{1,\alpha,c}^{f}$ with eigenvalue $\lambda\neq 0$. For
$x\in\mathcal{X}$, we have

$\displaystyle e^{-\lambda
t}|\phi(x)|=\left|\left(P_{1,\alpha,c}^{f}\right)^{t}\phi(x)\right|$
$\displaystyle=\left|\sum_{y\in\mathcal{X}}\left[\left(P_{1,\alpha,c}^{f}\right)^{t}(x,y)\phi(y)-\mu_{1,\alpha,c}^{f}(y)\phi(y)\right]\right|$
$\displaystyle\leqslant\|\phi\|_{\infty}2\left\|\left(P_{1,\alpha,c}^{f}\right)^{t}(x,\cdot)-\mu_{1,\alpha,c}^{f}\right\|_{TV}.$

Taking $x$ that attains $|\phi(x)|=\|\phi\|_{\infty}$ and
$\lambda=\lambda_{2}(-M_{1,\alpha,c}^{f})$ yield the desired result since

$e^{-\lambda_{2}(-M_{1,\alpha,c}^{f})t}\leqslant
2\left\|\left(P_{1,\alpha,c}^{f}\right)^{t}(x,\cdot)-\mu_{1,\alpha,c}^{f}\right\|_{TV}.$

As for the upper bound, it follows from (Levin et al., 2009, equation
$(20.21)$) that

$\displaystyle
2\left\|\left(P_{1,\alpha,c}^{f}\right)^{t}(x,\cdot)-\mu_{1,\alpha,c}^{f}\right\|_{TV}$
$\displaystyle\leqslant
e^{-\lambda_{2}(-M_{1,\alpha,c}^{f})t}\sum_{y\in\mathcal{X}}\frac{\mu_{1,\alpha,c}^{f}(y)}{\sqrt{\mu_{1,\alpha,c}^{f}(y)\mu_{1,\alpha,c}^{f}(x)}}$
$\displaystyle\leqslant\frac{e^{-\lambda_{2}(-M_{1,\alpha,c}^{f})t}}{\min_{x\in\mathcal{X}}\mu_{1,\alpha,c}^{f}(x)}$
$\displaystyle=\frac{Z_{1,\alpha,c}^{f}}{e^{-\max_{x\in\mathcal{X}}\mathcal{H}_{1,\alpha,c}^{f}(x)}}e^{-\lambda_{2}(-M_{1,\alpha,c}^{f})t}$
$\displaystyle\leqslant
e^{\mathcal{H}_{\max}-\mathcal{H}_{\min}}Z_{1,\alpha,c}^{f}e^{-\lambda_{2}(-M_{1,\alpha,c}^{f})t}$
$\displaystyle\leqslant
e^{\mathcal{H}_{\max}-\mathcal{H}_{\min}}|\mathcal{X}|e^{-\lambda_{2}(-M_{1,\alpha,c}^{f})t}.$

Next, we prove item (2). First, we compute the difference of
$Z_{1,\alpha,c}^{f}$ and $Z^{0}$:

$\displaystyle Z_{1,\alpha,c}^{f}-Z^{0}$
$\displaystyle=\sum_{x\in\mathcal{X}}e^{-\mathcal{H}_{1,\alpha,c}^{f}(x)}-e^{-(\mathcal{H}(x)-\mathcal{H}_{\min})}$
$\displaystyle=\sum_{x:\mathcal{H}(x)>c}e^{-\mathcal{H}_{1,\alpha,c}^{f}(x)}-e^{-(\mathcal{H}(x)-\mathcal{H}_{\min})}$
$\displaystyle\leqslant\sum_{x:\mathcal{H}(x)>c}\left|e^{-\mathcal{H}_{1,\alpha,c}^{f}(x)}-e^{-(\mathcal{H}(x)-\mathcal{H}_{\min})}\right|$
$\displaystyle\leqslant\sum_{x:\mathcal{H}(x)>c}\left|{\mathcal{H}_{1,\alpha,c}^{f}(x)}-{(\mathcal{H}(x)-\mathcal{H}_{\min})}\right|$
$\displaystyle=\sum_{x:\mathcal{H}(x)>c}\left|\int_{\mathcal{H}_{\min}}^{\mathcal{H}(x)}\left(\frac{1}{\alpha
f\left((u-c)_{+}\right)+1}-1\right)du\right|$
$\displaystyle\leqslant\sum_{x:\mathcal{H}(x)>c}\int_{\mathcal{H}_{\min}}^{\mathcal{H}(x)}{\alpha
f\left((u-c)_{+}\right)}du\leqslant\alpha{\xi}$

where
$\xi:=|{x:\mathcal{H}(x)>c}|({\mathcal{H}_{\max}-\mathcal{H}_{\min}}){f((\mathcal{H}_{\max}-c)_{+})}$.
Using the bound above, we calculate the total variation distance between
$\mu_{1,\alpha,c}^{f}$ and $\mu$:

$\displaystyle\|\mu_{1,\alpha,c}^{f}-\mu\|_{TV}$
$\displaystyle=\frac{1}{2}\sum_{x\in\mathcal{X}}|\mu_{1,\alpha,c}^{f}(x)-\mu(x)|$
$\displaystyle=\frac{1}{2}\sum_{x\in\mathcal{X}}\left|\frac{e^{-\mathcal{H}_{1,\alpha,c}^{f}(x)}}{Z_{1,\alpha,c}^{f}}-\frac{e^{-(\mathcal{H}(x)-\mathcal{H}_{\min})}}{Z^{0}}\right|$
$\displaystyle=\frac{1}{2}\sum_{x\in\mathcal{X}}\left|\frac{e^{-\mathcal{H}_{1,\alpha,c}^{f}(x)}}{Z_{1,\alpha,c}^{f}}-\frac{e^{-(\mathcal{H}(x)-\mathcal{H}_{\min})}}{Z_{1,\alpha,c}^{f}}+\frac{e^{-(\mathcal{H}(x)-\mathcal{H}_{\min})}}{Z_{1,\alpha,c}^{f}}-\frac{e^{-(\mathcal{H}(x)-\mathcal{H}_{\min})}}{Z^{0}}\right|$
$\displaystyle\leqslant\frac{1}{2}\frac{Z_{1,\alpha,c}^{f}-Z^{0}}{Z_{1,\alpha,c}^{f}}+\frac{1}{2}\frac{Z_{1,\alpha,c}^{f}-Z^{0}}{Z_{1,\alpha,c}^{f}}\leqslant\dfrac{\alpha\xi}{|\mathcal{X}|}e^{\max_{x}\mathcal{H}_{1,\alpha,c}^{f}(x)}\leqslant\dfrac{\alpha\xi}{|\mathcal{X}|}e^{\mathcal{H}_{\max}-\mathcal{H}_{\min}}=\mathcal{O}(\alpha).$

###  4.2. Proof of Theorem 4.2

First, we prove item (1). Using the Markov chain law of large number, we have

$\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}\dfrac{1}{t}\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds\rightarrow^{a.s.}\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}\mu_{1,\alpha,c}^{f}(wg)=\mu(g).$

Applying the above result to $g=1$ yields

$\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds\rightarrow^{a.s.}\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}\mu_{1,\alpha,c}^{f}(w)=\mu(1)=1.$

Collecting the two results above together with the continuous mapping theorem
gives the desired result.

Next, we prove item (2). Using the Markov chain central limit theorem, we have

$\sqrt{t}\left(\frac{1}{t}\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds-\mu_{1,\alpha,c}^{f}(wg)\right)\rightarrow^{d}N(0,\mathrm{Var}_{\mu_{1,\alpha,c}^{f}}(wg)).$

By Slutsky’s theorem, as
$\frac{1}{t}\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds\rightarrow^{a.s.}\dfrac{Z^{0}}{Z_{1,\alpha,c}^{f}}$,
the left-hand-side can be expressed as

$\sqrt{t}\left(\dfrac{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds}{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds}-\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}\mu_{1,\alpha,c}^{f}(wg)\right)\rightarrow^{d}\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}N(0,\mathrm{Var}_{\mu_{1,\alpha,c}^{f}}(wg)).$

Finally, using
$\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}\mu_{1,\alpha,c}^{f}(wg)=\mu(g)$ gives the
desired result since

$\sqrt{t}\left(\dfrac{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds}{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds}-\mu(g)\right)\rightarrow^{d}\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}N(0,\mathrm{Var}_{\mu_{1,\alpha,c}^{f}}(wg)).$

###  4.3. Proof of Theorem 4.3

Using equation $(1.2)$ in Lezaud (2001), for $g$ such that $\left\lVert
g\right\rVert_{\infty}\leqslant a$, if we have

$t\geqslant\dfrac{32a^{2}}{\lambda_{2}(-M_{1,\alpha,c}^{f})\gamma^{2}}\ln\left(\dfrac{1}{\varepsilon\min_{x\in\mathcal{X}}\mu_{1,\alpha,c}^{f}(x)}\right),$

then

$\displaystyle\mathbb{P}_{x}\left(\dfrac{1}{t}\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds-\mu_{1,\alpha,c}^{f}(wg)\geqslant\gamma\right)$
$\displaystyle\leqslant\varepsilon,$
$\displaystyle\mathbb{P}_{x}\left(\dfrac{1}{t}\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds-\dfrac{Z^{0}}{Z_{1,\alpha,c}^{f}}\mu(g)\geqslant\gamma\right)$
$\displaystyle\leqslant\varepsilon.$

Specializing into $g=-1$ gives, for

$t\geqslant\dfrac{32}{\lambda_{2}(-M_{1,\alpha,c}^{f})\gamma^{2}}\ln\left(\dfrac{1}{\varepsilon\min_{x\in\mathcal{X}}\mu_{1,\alpha,c}^{f}(x)}\right),$

then

$\displaystyle\mathbb{P}_{x}\left(\dfrac{Z^{0}}{Z_{1,\alpha,c}^{f}}-\gamma\geqslant\dfrac{1}{t}\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds\right)$
$\displaystyle\leqslant\varepsilon.$

From now on, without loss of generality, we assume that $g\geqslant 0$, and so

$\displaystyle\mathbb{P}_{x}$
$\displaystyle\left(\dfrac{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds}{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds}-\mu(g)\geqslant\dfrac{\gamma}{\dfrac{Z^{0}}{Z_{1,\alpha,c}^{f}}-\gamma}\right)$
$\displaystyle\leqslant\mathbb{P}_{x}\left(\dfrac{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds}{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds}-\mu(g)\geqslant\dfrac{\gamma+\gamma\mu(g)}{\dfrac{Z^{0}}{Z_{1,\alpha,c}^{f}}-\gamma}\right)\leqslant\varepsilon.$

The desired result follows by taking

$\eta:=\dfrac{\gamma}{\dfrac{Z^{0}}{Z_{1,\alpha,c}^{f}}-\gamma},$

and rearranging gives

$\gamma=\dfrac{Z^{0}}{Z_{1,\alpha,c}^{f}}\dfrac{\eta}{1+\eta}.$

##  5\. Time-dependent penalty $(\alpha_{t})_{t\geqslant 0}$

In Section 4, we focus on the constant penalty parameter case: we run a
homogeneous Markov chain
$X_{1,\alpha,c}^{f}=(X_{1,\alpha,c}^{f}(t))_{t\geqslant 0}$ at a fixed
$\alpha$, and correct the bias by proposing a self-normalized estimator. In
this Section, instead of fixing $\alpha$, we anneal the penalty parameter by
sending $\alpha_{t}\to 0$ as the time $t\to\infty$. This results in running a
non-homogeneous Markov chain
$X_{1,\bm{\alpha},c}^{f}=(X_{1,\bm{\alpha},c}^{f}(t))_{t\geqslant 0}$ to
approximately sample from $\mu$, where $\bm{\alpha}=(\alpha_{t})_{t\geqslant
0}$. Unlike simulated annealing where the temperature parameter is annealed,
we instead anneal the penalty parameter so that the state-dependent effect
diminishes to zero in the long run, which allows for approximate sampling from
the target distribution $\mu$.

In the following, we shall denote, for $x\in\mathcal{X}$ and
$A\subseteq\mathcal{X}$,

$(P_{1,\bm{\alpha},c}^{f})^{t}(x,A):=e^{\int_{0}^{t}M^{f}_{1,\alpha_{s},c}\,ds}(x,A)=\mathbb{P}_{x}(X_{1,\bm{\alpha},c}^{f}(t)\in
A).$

Our first main result states that the non-homogeneous chain
$X_{1,\bm{\alpha},c}^{f}$ is weakly ergodic as long as $\alpha^{\prime}_{t}\to
0$. In other words, if the change of $\alpha_{t}$ approaches $0$ as
$t\to\infty$, the non-homogeneous chain tracks closely its instantaneous
stationary distribution $\mu_{1,\alpha_{t},c}^{f}$:

######  Theorem 5.1 (Weak ergodicity of $X_{1,\bm{\alpha},c}^{f}$).

If $\alpha^{\prime}_{t}\to 0$ as $t\to\infty$, then for any $x\in\mathcal{X}$,

(5.1)
$\left\|\left(P_{1,\bm{\alpha},c}^{f}\right)^{t}(x,\cdot)-\mu_{1,\alpha_{t},c}^{f}\right\|_{TV}\to
0\text{ as }t\to\infty.$

The proof of Theorem 5.1 is given in Section 5.1.

If we impose the additional assumption that $\alpha_{t}\downarrow 0$, then the
non-homogeneous chain $X_{1,\bm{\alpha},c}^{f}$ is strongly ergodic and
converges to the desired target distribution $\mu$ in total variation
distance:

######  Theorem 5.2 (Strong ergodicity of $X_{1,\bm{\alpha},c}^{f}$ ).

If $\alpha_{t}\downarrow 0,~{}\alpha^{\prime}_{t}\to 0$ as $t\to\infty$, then
for any $x\in\mathcal{X}$,

$\left\|\left(P_{1,\bm{\alpha},c}^{f}\right)^{t}(x,\cdot)-\mu\right\|_{TV}\to
0\text{ as }t\to\infty.$

Unlike simulated annealing where logarithmic cooling schedule of the
temperature is commonly employed to ensure strong ergodicity Choi (2021b);
Miclo (1992); Monmarché (2018), in this context the penalty parameter can
decay to zero at a faster rate than logarithmic cooling while maintaining
strong ergodicity. For instance, we can adopt the exponential annealing
schedule by taking $\alpha_{t}=e^{-t}$. With this choice, we note that the
conditions in Theorem 5.2 are satisfied.

The proof of Theorem 5.2 is deferred to Section 5.2.

In our next result, we quantify the convergence rate of the total variation
mixing by specializing into the case where $f(x)=x^{2}$. The primary reason
for considering quadratic $f$ is that the reduced critical height can be
bounded independently of the system size, leading to rapid mixing in the
examples to be presented in Section 6.

We would be interested in estimating the total variation mixing time
$t_{mix}^{f}(\varepsilon)$ of the non-homogeneous chain
$X_{1,\bm{\alpha},c}^{f}$, defined to be for $\varepsilon>0$,

(5.2) $\displaystyle t_{mix}^{f}(\varepsilon):=\inf\bigg{\\{}t\geqslant
0;~{}\max_{x\in\mathcal{X}}\left\|\left(P_{1,\bm{\alpha},c}^{f}\right)^{t}(x,\cdot)-\mu\right\|_{TV}\leqslant\varepsilon\bigg{\\}}.$

Note that the dependence on $\bm{\alpha}$ and $c$ are suppressed on the left
hand side.

######  Theorem 5.3 (Total variation mixing time with $f(x)=x^{2}$ and
logarithmic annealing schedule).

Suppose that $f(x)=x^{2}$, where we recall $\mathcal{H}_{1,\alpha,c}^{f}$ is
computed in (2.5). For a large enough universal constant $C>0$ that is
independent of $\mathcal{X}$ or $\mathcal{H}$, denote

$\beta_{t}:=\frac{1}{\alpha_{t}}=\frac{2}{3\pi}\ln(e^{\frac{3\pi}{2}}+pt),\quad
p:=\frac{\pi
K}{\mathcal{M}},\quad\mathcal{M}:=C\left((\ln|\mathcal{X}|)\vee\left(\mathcal{H}_{\max}-\mathcal{H}_{\min}\right)\right).$

For $t\geqslant 0$, we employ the following logarithmic annealing schedule on
the penalty parameter with

$\alpha_{t}:=\frac{3\pi}{2}\dfrac{1}{\ln(e^{\frac{3\pi}{2}}+pt)}.$

Assume further that there exists a positive constant $K=K(c)$ such that

(5.3) $\lambda_{2}(-M_{1,\alpha_{t},c}^{f})\geqslant
Ke^{-\frac{3\pi}{2}\sqrt{\beta_{t}}}\geqslant Ke^{-\frac{3\pi}{2}\beta_{t}},$

and we choose $c$ such that $c-\mathcal{H}_{\min}=\delta=\mathcal{O}(1)$. For
any $x\in\mathcal{X}$, we have

$\left\|\left(P_{1,\bm{\alpha},c}^{f}\right)^{t}(x,\cdot)-\mu\right\|_{TV}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{O}(|\mathcal{X}|e^{\mathcal{M}}e^{-\mathcal{M}\beta_{t}}).}$

In particular,

$\displaystyle t_{mix}^{f}(\varepsilon)$
$\displaystyle=\mathcal{O}\left(\dfrac{1}{p}\left(\exp\bigg{\\{}\dfrac{3\pi(\ln|\mathcal{X}|+\mathcal{M}+\log(1/\varepsilon))}{2\mathcal{M}}\bigg{\\}}-\exp\\{3\pi/2\\}\right)\right)$
$\displaystyle=\mathcal{O}\left(\dfrac{(\ln|\mathcal{X}|)\vee\left(\mathcal{H}_{\max}-\mathcal{H}_{\min}\right)}{K}\left(\dfrac{1}{\varepsilon}\right)^{\frac{3\pi}{2\left((\ln|\mathcal{X}|)\vee\left(\mathcal{H}_{\max}-\mathcal{H}_{\min}\right)\right)}}\right).$

The proof of Theorem 5.3 is given in Section 5.3.

###  5.1. Proof of Theorem 5.1

In view of (Gidas, 1985, Theorem $A$), it suffices to find
$(\gamma_{t})_{t\geqslant 0}$ that satisfies

(5.4)
$\left|\frac{d}{dt}\mu_{1,\alpha_{t},c}^{f}(x)\right|\leqslant\gamma_{t}\cdot\mu_{1,\alpha_{t},c}^{f}(x),$

where $\gamma_{t}$ is independent of $x$ and fulfills

(5.5)
$\lim_{t\rightarrow\infty}\frac{\gamma_{t}}{\lambda_{2}(-M_{1,\alpha_{t},c}^{f})}=0$

with

(5.6)
$\displaystyle\int_{0}^{\infty}\lambda_{2}(-M_{1,\alpha_{t},c}^{f})dt=+\infty.$

First, (5.6) is satisfied due to Corollary (3.1). We proceed to show (5.4).
Denote the negative time derivative of $\mathcal{H}_{1,\alpha_{t},c}^{f}$ as

$G_{1,\alpha_{t},c}^{f}(x):=-\frac{d}{dt}\mathcal{H}_{1,\alpha_{t},c}^{f}(x)\leqslant
0.$

We observe that

$G_{1,\alpha_{t},c}^{f}(x)=\int_{\mathcal{H}_{\min}}^{\mathcal{H}(x)}\frac{\alpha^{\prime}_{t}f((u-c)_{+})}{(\alpha_{t}f\left((u-c)_{+}\right)+1)^{2}}du\geqslant\alpha_{t}^{\prime}f(\mathcal{H}_{\max}-c)(\mathcal{H}_{\max}-\mathcal{H}_{\min})$

and

$\displaystyle\frac{d}{dt}\mu_{1,\alpha_{t},c}^{f}(x)$
$\displaystyle=\frac{1}{Z_{1,\alpha_{t},c}^{f}}\frac{d}{dt}\left(e^{-\mathcal{H}_{1,\alpha_{t},c}^{f}(x)}\right)+e^{-\mathcal{H}_{1,\alpha_{t},c}^{f}(x)}\frac{d}{dt}\left(\frac{1}{Z_{1,\alpha_{t},c}^{f}}\right)$
$\displaystyle=G_{1,\alpha_{t},c}^{f}(x)\cdot\mu_{1,\alpha_{t},c}^{f}(x)-e^{-\mathcal{H}_{1,\alpha_{t},c}^{f}(x)}\left(\frac{1}{Z_{1,\alpha_{t},c}^{f}}\right)^{2}\sum_{x\in\mathcal{X}}e^{-\mathcal{H}_{1,\alpha_{t},c}^{f}(x)}\cdot
G_{1,\alpha_{t},c}^{f}(x)$
$\displaystyle=G_{1,\alpha_{t},c}^{f}(x)\cdot\mu_{1,\alpha_{t},c}^{f}(x)-\mu_{1,\alpha_{t},c}^{f}(x)\sum_{x\in\mathcal{X}}\mu_{1,\alpha_{t},c}^{f}(x)\cdot
G_{1,\alpha_{t},c}^{f}(x),$

which implies

$\displaystyle\left|\frac{d}{dt}\mu_{1,\alpha_{t},c}^{f}(x)\right|$
$\displaystyle\leqslant\mu_{1,\alpha_{t},c}^{f}(x)\left(\max_{x}G_{1,\alpha_{t},c}^{f}(x)-\min_{x}G_{1,\alpha_{t},c}^{f}(x)\right)$
$\displaystyle\leqslant-\alpha_{t}^{\prime}f(\mathcal{H}_{\max}-c)(\mathcal{H}_{\max}-\mathcal{H}_{\min})\cdot\mu_{1,\alpha_{t},c}^{f}(x).$

Then, (5.4) is satisfied with
$\gamma_{t}=-\alpha_{t}^{\prime}f(\mathcal{H}_{\max}-c)(\mathcal{H}_{\max}-\mathcal{H}_{\min})$
and (5.5) is satisfied since $\alpha^{\prime}_{t}\to 0$ as $t\to\infty$ and
the time-independent lower bound of the spectral gap as in Corollary (3.1).
Thus, (5.1) is proved by applying (Gidas, 1985, Theorem A).

###  5.2. Proof of Theorem 5.2

Using the triangle inequality, we have

$\left\|\left(P_{1,\bm{\alpha},c}^{f}\right)^{t}(x,\cdot)-\mu\right\|_{TV}\leqslant\left\|\left(P_{1,\bm{\alpha},c}^{f}\right)^{t}(x,\cdot)-\mu_{1,\alpha_{t},c}^{f}\right\|_{TV}+\|\mu_{1,\alpha_{t},c}^{f}-\mu\|_{TV}\to
0\text{ as }t\to\infty$

where the first term converges to $0$ due to Theorem 5.1, and we recall from
Theorem 4.1 that the second term is

$\|\mu_{1,\alpha_{t},c}^{f}-\mu\|_{TV}=\mathcal{O}(\alpha_{t}).$

###  5.3. Proof of Theorem 5.3

We follow the same roadmap as the proof of (Deuschel and Mazza, 1994, Appendix
A.1). First, using the Cauchy-Schwartz inequality, we have

$\displaystyle\left\|\left(P_{1,\bm{\alpha},c}^{f}\right)^{t}(x,\cdot)-\mu\right\|_{TV}$
$\displaystyle=\frac{1}{2}\sum_{y\in\mathcal{X}}\left|\left(P_{1,\bm{\alpha},c}^{f}\right)^{t}(x,y)-\mu(y)\right|$
$\displaystyle\leqslant\frac{1}{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sqrt{\sum_{y\in\mathcal{X}}\left(\frac{\left(P_{1,\bm{\alpha},c}^{f}\right)^{t}(x,y)-\mu(y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\right)^{2}\mu_{1,\alpha_{t},c}^{f}(y)}=:\frac{1}{2}\sqrt{u_{t}(x)}}$

where

$u_{t}(x):=\sum_{y\in\mathcal{X}}\left(g_{t}(x,y)-\frac{\mu(y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\right)^{2}\mu_{1,\alpha_{t},c}^{f}(y),\quad
g_{t}(x,y):=\frac{\left(P_{1,\bm{\alpha},c}^{f}\right)^{t}(x,y)}{\mu_{1,\alpha_{t},c}^{f}(y)}.$

Next, we calculate the time derivative of $u_{t}$:

$\displaystyle\frac{du_{t}(x)}{dt}$
$\displaystyle=\frac{d}{dt}\left(\sum_{y\in\mathcal{X}}\left(g_{t}(x,y)-\frac{\mu(y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\right)^{2}\mu_{1,\alpha_{t},c}^{f}(y)\right)$
$\displaystyle=\sum_{y\in\mathcal{X}}2\left(g_{t}(x,y)-\frac{\mu(y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\right)\left(\mu_{1,\alpha_{t},c}^{f}(y)\frac{dg_{t}(x,y)}{dt}+\frac{\mu(y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\frac{d\mu_{1,\alpha_{t},c}^{f}(y)}{dt}\right)$
$\displaystyle\qquad\qquad\qquad\qquad+\sum_{y\in\mathcal{X}}\left(g_{t}(x,y)-\frac{\mu(y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\right)^{2}\frac{d\mu_{1,\alpha_{t},c}^{f}(y)}{dt}$
$\displaystyle=\sum_{y\in\mathcal{X}}2\left(g_{t}(x,y)-\frac{\mu(y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\right)\left(M_{1,\alpha_{t},c}^{f}g_{t}(x,y)-\frac{g_{t}(x,y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\frac{d\mu_{1,\alpha_{t},c}^{f}(y)}{dt}\right)\mu_{1,\alpha_{t},c}^{f}(y)$
$\displaystyle\qquad+\sum_{y\in\mathcal{X}}2\left(g_{t}(x,y)-\frac{\mu(y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\right)\frac{\mu(y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\cdot\frac{d\mu_{1,\alpha_{t},c}^{f}(y)}{dt}+\sum_{y\in\mathcal{X}}\left(g_{t}(x,y)-\frac{\mu(y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\right)^{2}\frac{d\mu_{1,\alpha_{t},c}^{f}(y)}{dt}$
(5.7) $\displaystyle=-2\langle-
M_{1,\alpha_{t},c}^{f}g_{t},g_{t}\rangle_{\mu_{1,\alpha_{t},c}^{f}}-\sum_{y\in\mathcal{X}}\left(g_{t}(x,y)-\frac{\mu(y)}{\mu_{1,\alpha_{t},c}^{f}(y)}\right)^{2}\frac{d\mu_{1,\alpha_{t},c}^{f}(y)}{dt},$

where the third equality follows from (Deuschel and Mazza, 1994, Lemma A.1.1)
and the reversibility of $M^{f}_{1,\alpha_{t},c}$ with respect to
$\mu_{1,\alpha_{t},c}^{f}$.

In the special case where $f(x)=x^{2}$, $\mathcal{H}_{1,\alpha,c}^{f}$ is
computed as in (2.5), we can then calculate the derivative of
$\mu_{1,\alpha_{t},c}^{f}$ with respect to $\beta=1/\alpha$:

(5.8)
$\displaystyle\frac{d}{d\beta}\mu_{1,\alpha,c}^{f}(y)=\mu_{1,\alpha,c}^{f}(y)\left(-\frac{d\mathcal{H}_{1,\alpha,c}^{f}(y)}{d\beta}\right)-\mu_{1,\alpha,c}^{f}(y)\sum_{x\in\mathcal{X}}\mu_{1,\alpha,c}^{f}(x)\left(-\frac{d\mathcal{H}_{1,\alpha,c}^{f}(x)}{d\beta}\right)$

with

$-\frac{d\mathcal{H}_{1,\alpha,c}^{f}(y)}{d\beta}=\left\\{\begin{array}[]{ll}0,&H(y)\leqslant
c,\\\
\dfrac{\mathcal{H}(y)-c}{2(\mathcal{H}(y)-c)^{2}+2\beta}-\dfrac{\text{arctan}\left(\sqrt{\frac{1}{\beta}}\left(\mathcal{H}(y)-c\right)\right)}{2\sqrt{\beta}},&H(y)>c.\end{array}\right.$

We also have the bounds that

$\frac{-\pi}{4\sqrt{\beta}}\leqslant-\frac{d\mathcal{H}_{1,\alpha,c}^{f}(y)}{d\beta}\leqslant\frac{1}{2\beta}\left(\mathcal{H}_{\max}-\mathcal{H}_{\min}\right),$

which when combine with (5.8) implies

(5.9)
$\displaystyle\left|\frac{d}{d\beta}\ln\mu_{1,\alpha,c}^{f}(y)\right|\leqslant\max_{y\in\mathcal{X}}\left(-\frac{d\mathcal{H}_{1,\alpha,c}^{f}(y)}{d\beta}\right)-\min_{x\in\mathcal{X}}\left(-\frac{d\mathcal{H}_{1,\alpha,c}^{f}(x)}{d\beta}\right)=\mathcal{O}\left(\frac{1}{\sqrt{\beta}}\left(\mathcal{H}_{\max}-\mathcal{H}_{\min}\right)\right).$

With the choice of $K,\beta_{t}$ and $\mathcal{M}$ in the theorem, we have

$Ke^{-\frac{\pi}{2}\beta_{t}}-\frac{\mathcal{M}}{2}\frac{d\beta_{t}}{dt}=\mathcal{M}\dfrac{d\beta_{t}}{dt},$

and using (5.7) and (5.3) give

$\frac{du_{t}(x)}{dt}\leqslant-2\left(Ke^{-\frac{\pi}{2}\beta_{t}}-\frac{\mathcal{M}}{2}\frac{d\beta_{t}}{dt}\right)u_{t}(x)=-2\mathcal{M}\cdot\frac{d\beta_{t}}{dt}\cdot
u_{t}(x)$

which implies that

$u_{t}(x)\leqslant u_{0}(x)e^{-2\mathcal{M}(\beta_{t}-1)}$

with

$\displaystyle u_{0}(x)$
$\displaystyle=\left(\frac{1-\mu(x)}{\mu_{1,\alpha_{0},c}^{f}(x)}\right)^{2}\mu_{1,\alpha_{0},c}^{f}(x)+\sum_{y\neq
x}\left(\frac{\mu(y)}{\mu_{1,\alpha_{0},c}^{f}(y)}\right)^{2}\mu_{1,\alpha_{0},c}^{f}(y)$
$\displaystyle\leqslant\frac{1}{\min_{x\in\mathcal{X}}\mu_{1,\alpha_{0},c}^{f}(x)}+\max_{x\in\mathcal{X}}\left(\frac{\mu(x)}{{\mu}_{1,\alpha_{0},c}^{f}(x)}\right)^{2}$
$\displaystyle\leqslant
e^{c-\mathcal{H}_{\min}+\frac{1}{\sqrt{\alpha_{0}}\cdot}\frac{\pi}{2}}+\left(\frac{Z_{1,\alpha_{0},c}^{f}}{Z^{0}}\right)^{2}\leqslant\mathcal{O}(|\mathcal{X}|^{2}),$

where the last inequality follows from the assumption that
$c-\mathcal{H}_{\min}=\delta=\mathcal{O}(1)$. The desired result follows.

##  6\. Applications

In this Section, we shall apply the model-independent mixing time results that
we obtain in Section 4 and Section 5 to some common statistical physics models
where the Glauber dynamics is utilized for sampling from the target Gibbs
distribution $\mu$. In particular, we will be interested in comparing the
following mixing parameters: for any $\eta,\varepsilon>0$, $g\in\ell^{2}(\mu)$

  * •

(Large deviation from the target mean) Recall that in (4.2), we introduce

$t_{conc}(\eta,\varepsilon)=\inf\bigg{\\{}t\geqslant
0;~{}\sup_{g\in\ell^{2}(\mu);~{}\left\lVert g\right\rVert_{\infty}\leqslant
a}\sup_{x\in\mathcal{X}}\mathbb{P}_{x}\left(\dfrac{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))g(X_{1,\alpha,c}^{f}(s))\,ds}{\int_{0}^{t}w(X_{1,\alpha,c}^{f}(s))\,ds}-\mu(g)\geqslant\eta\right)\leqslant\varepsilon\bigg{\\}}.$

  * •

(Total variation mixing time of the modified dynamics with $f(x)=x^{2}$ and
logarithmic annealing schedule of the penalty parameter) Recall that in (5.2),
we introduce

$t_{mix}^{f}(\varepsilon)=\inf\bigg{\\{}t\geqslant
0;~{}\max_{x\in\mathcal{X}}\left\|\left(P_{1,\bm{\alpha},c}^{f}\right)^{t}(x,\cdot)-\mu\right\|_{TV}\leqslant\varepsilon\bigg{\\}}.$

  * •

(Total variation mixing time of the original dynamics) Recall that in (2.2),
we define the total variation mixing time of the MH dynamics that target $\mu$
to be

$t_{mix}^{0}(\varepsilon)=\inf\bigg{\\{}t\geqslant
0;~{}\max_{x\in\mathcal{X}}\left\|\left(P^{0}\right)^{t}(x,\cdot)-\mu\right\|_{TV}\leqslant\varepsilon\bigg{\\}}.$

We shall investigate the ferromagnetic Ising model on the $n$-dimensional
hypercube in Section 6.1, the ferromagnetic Ising model on the complete graph
in Section 6.2, and finally the $q$-state Potts model on the two-dimensional
torus in Section 6.3.

###  6.1. Ferromagnetic Ising model on the $n$-dimensional hypercube

In this Section, we specialize into the case of the Ising model on a finite
connected graph $G=(V,E)$. Let $\mathcal{X}=\\{-1,1\\}^{V}$ be the set of
possible configurations $\sigma=\\{\sigma(v);~{}v\in V\\}$ in which we
associate each vertex $v\in V$ with a spin value $\sigma(v)\in\\{-1,1\\}$. The
Hamiltonian is given by

(6.1) $\displaystyle\mathcal{H}(\sigma)=-\dfrac{J}{2}\sum_{(v,w)\in
E}\sigma(v)\sigma(w)-\dfrac{h}{2}\sum_{v\in V}\sigma(v),$

where $J=\mathcal{O}(1)>0$ is the ferromagnetic pair potential and
$h=\mathcal{O}(1)>0$ is the external magnetic field. We also denote by
$\mathbf{+1}$ (resp. $\mathbf{-1}$) to be the all-plus (resp. all-minus)
configuration, where we recall that $\mathbf{+1}$ is the ground-state or
stable configuration while $\mathbf{-1}$ is the metastable configuration of
$\eqref{eq:Hamhypercube}$.

For the Glauber dynamics, we consider a simple random walk proposal with
transition matrix
$P^{SRW}=(P^{SRW}(\sigma_{1},\sigma_{2}))_{\sigma_{1},\sigma_{2}\in\mathcal{X}}$
given by

$\displaystyle P^{SRW}(\sigma_{1},\sigma_{2})$
$\displaystyle:=\dfrac{1}{n}\mathds{1}_{\\{\textrm{there exists }i\textrm{
such that }\sigma_{1}(i)=-\sigma_{2}(i)\textrm{ and
}\sigma_{1}(j)=\sigma_{2}(j)\textrm{ for all }j\neq i\\}}.$

The stationary distribution of $P^{SRW}$ is the uniform distribution on
$\mathcal{X}$.

######  Theorem 6.1.

Suppose that $G=(V,E)$ is an $n$-dimensional hypercube. Consider the Glauber
dynamics on $\mathcal{X}=\\{-1,1\\}^{V}$ with proposal generator $P^{SRW}-I$,
target distribution $\mu_{1,\alpha,c}^{f}$ with Hamiltonian given in (6.1),
$f(x)=x^{2}$ and $c=\mathcal{H}(\mathbf{+1})+\delta$, where $\delta>0$. Given
$\eta,\varepsilon>0$, we have

  1. (1)

for $g\in\ell^{2}(\mu)$ such that $\left\lVert g\right\rVert_{\infty}\leqslant
a$,

$t_{conc}(\eta,\varepsilon)=\mathcal{O}\left(n^{2}2^{n}e^{3\delta+\frac{3}{\sqrt{\alpha}}\frac{\pi}{2}}\dfrac{a^{2}(1+\eta)^{2}}{\eta^{2}}\ln\left(\dfrac{1}{\varepsilon}\right)\right).$

  2. (2)

[Rapid mixing: mixing time is at most polynomial in the number of vertices in
the modified Glauber dynamics with time-dependent penalty] Assume the
logarithmic annealing schedule
$\alpha_{t}=\frac{3\pi}{2}\dfrac{1}{\ln(e^{\frac{3\pi}{2}}+pt)}$, where $p$ is
chosen as in Theorem 5.3, then

$t_{mix}^{f}(\varepsilon)=\mathcal{O}\left(n^{2}2^{n}e^{3\delta}\left(\dfrac{1}{\varepsilon}\right)^{\frac{3\pi}{2n\ln
2}}\right).$

  3. (3)

[Torpid mixing: mixing time is at least exponential in the number of vertices
in the original Glauber dynamics] Assume in addition that $0\leqslant
h/J\leqslant n$. For constant $\mathbb{M}=\mathbb{M}(J,h)$ that does not
depend on $n$,

$t_{mix}^{0}(\varepsilon)=\Omega\left(e^{\mathbb{M}2^{n}}\ln\left(\dfrac{1}{2\varepsilon}\right)\right).$

######  Remark 6.1 (On the choice of $c$).

We notice that in both item (1) and item (2), the upper bounds there depend on
$\delta$ through $\mathcal{O}(e^{3\delta})$. Ideally, the parameter $c$ should
be chosen to be as close as $\mathcal{H}(\mathbf{+1})$ so that $\delta$ is
small enough. We also utilize the fact that the ground-state configuration in
this setting is the all-ones configuration $\mathbf{+1}$.

We observe that the original Glauber dynamics has an exponentially large
mixing time in the system size $2^{n}$, while the proposed annealing algorithm
mixes rapidly in the sense that the mixing time is at most polynomial in
$2^{n}$.

####  6.1.1. Proof of Theorem 6.1

In this subsection, we give the proof of Theorem 6.1. We first state a
spectral gap lower bound result that will be used in subsequent sections:

######  Lemma 6.1.

Given a graph $G=(V,E)$, consider the Glauber dynamics on
$\mathcal{X}=\\{-1,1\\}^{V}$ with proposal generator $P^{SRW}-I$, target
distribution $\mu_{1,\alpha,c}^{f}$ and denote $n=|V|$, then

$\displaystyle\lambda_{2}(-M_{1,\alpha,c}^{f})\geqslant\dfrac{2}{n}e^{-3\max_{\sigma}\mathcal{H}_{1,\alpha,c}^{f}(\sigma)}.$

In particular, if we take $f(x)=x^{2}$ and
$c=\mathcal{H}(\mathbf{+1})+\delta$, where $\delta>0$, we have

$\displaystyle\lambda_{2}(-M_{1,\alpha,c}^{f})\geqslant\dfrac{2}{n}e^{-3\delta-\frac{3}{\sqrt{\alpha}}\frac{\pi}{2}}.$

[Proof of Lemma 6.1. ] To prove the desired result, we shall invoke the
following classical stability result of the spectral gap by Holley and Stroock
(1987):

######  Lemma 6.2.

Let $(\mathcal{X},d,\mu)$ be a finite metric measure space and suppose
$d\nu=\frac{e^{-U(x)}}{Z}d\mu$. Suppose that $L_{1}$ is a reversible generator
with respect to $\nu$, and that $L_{2}$ is another reversible generator with
respect to $\mu.$ Suppose furthermore that $AL_{2}(x,y)\leqslant L_{1}(x,y)$
for all $x,y\in\mathcal{X}$. Then if $L_{2}$ has spectral gap
$\lambda_{2}(-L_{2})$, the spectral gap of $L_{1}$ satisfies

$Ae^{-2(\max U-\min U)}\lambda_{2}(-L_{2})\leqslant\lambda_{2}(-L_{1}).$

The desired result follows if we take $L_{1}=M_{1,\alpha,c}^{f}$,
$L_{2}=P^{SRW}-I$ and
$A=\exp\bigg{\\{}-\max_{\sigma}\mathcal{H}_{1,\alpha,c}^{f}(\sigma)\bigg{\\}},$
and recall that in (Diaconis and Saloff-Coste, 1996, Example $3.2$) the
spectral gap of $P^{SRW}$ is given by $\lambda_{2}(-(P^{SRW}-I))=2/n$.

In the special case when $f(x)=x^{2}$ and $c=\mathcal{H}(\mathbf{+1})+\delta$,
we utilize Section 2.3 and the fact that $\sup_{x}\arctan(x)=\pi/2$.

We now give the proof of Theorem 6.1.

[Proof of Theorem 6.1. ] First, we prove item (1). The desired result follows
from Theorem 4.3 and the spectral gap lower bound Lemma 6.1 as well as the
inequality that

$\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}\leqslant Z_{1,\alpha,c}^{f}\leqslant
2^{n}.$

Next, we prove item (2), which simply follows from Theorem 5.3 and the
observation from Lemma 6.1 that

$\mathcal{H}_{\max}-\mathcal{H}_{\min}=\mathcal{O}(n2^{n}),\quad
K=\dfrac{2}{n}e^{-3\delta}.$

Finally, we prove item (3). Using (Levin and Peres, 2017, Lemma $20.11$), we
see that

$\displaystyle
t_{mix}^{0}(\varepsilon)\geqslant\dfrac{1}{\lambda_{2}(-M^{0})}\ln\left(\dfrac{1}{2\varepsilon}\right).$

According to Jovanovski (2017); Dommers et al. (2017), the critical height
$m^{0}$ is bounded below by

$m^{0}\geqslant\mathbb{M}2^{n},$

and the desired result follows from (Holley and Stroock, 1988, Lemma $2.3$)
since

$\lambda_{2}(-M^{0})=\mathcal{O}(e^{\mathbb{M}2^{n}}/4^{n})=\mathcal{O}(e^{\mathbb{M}2^{n}}).$

Note that in the above equations the value of $\mathbb{M}$ may possibly be
changed in each line.

###  6.2. Ferromagnetic Ising model on the complete graph

In the second example, we consider the ferromagnetic Ising model on the
complete graph $G=K_{n}$ with $n=|V|$ vertices. The Hamiltonian remains to be
(6.1), except now the underlying graph structure is $K_{n}$.

######  Theorem 6.2.

Suppose that $G=K_{n}$ is the complete graph on $n$ vertices. Consider the
Glauber dynamics on $\mathcal{X}=\\{-1,1\\}^{V}$ with proposal generator
$P^{SRW}-I$, target distribution $\mu_{1,\alpha,c}^{f}$ with Hamiltonian given
in (6.1), $f(x)=x^{2}$ and $c=\mathcal{H}(\mathbf{+1})+\delta$, where
$\delta>0$. Given $\eta,\varepsilon>0$, we have

  1. (1)

For $g\in\ell^{2}(\mu)$ such that $\left\lVert g\right\rVert_{\infty}\leqslant
a$,

$t_{conc}(\eta,\varepsilon)=\mathcal{O}\left(n^{2}2^{n}e^{3\delta+\frac{3}{\sqrt{\alpha}}\frac{\pi}{2}}\dfrac{a^{2}(1+\eta)^{2}}{\eta^{2}}\ln\left(\dfrac{1}{\varepsilon}\right)\right).$

  2. (2)

[Rapid mixing: mixing time is at most polynomial in the number of vertices in
the modified Glauber dynamics with time-dependent penalty] Assume the
logarithmic annealing schedule
$\alpha_{t}=\frac{3\pi}{2}\dfrac{1}{\ln(e^{\frac{3\pi}{2}}+pt)}$, where $p$ is
chosen as in Theorem 5.3, then

$t_{mix}^{f}(\varepsilon)=\mathcal{O}\left(n^{3}e^{3\delta}\left(\dfrac{1}{\varepsilon}\right)^{\frac{3\pi}{2n\ln
2}}\right).$

  3. (3)

[Torpid mixing: mixing time is at least exponential in the number of vertices
in the original Glauber dynamics] Assume in addition that $h/J$ is not an
integer. For constant $\mathbb{M}=\mathbb{M}(J,h)$ that does not depend on
$n$,

$t_{mix}^{0}(\varepsilon)=\Omega\left(e^{\mathbb{M}n^{2}}\ln\left(\dfrac{1}{2\varepsilon}\right)\right).$

######  Remark 6.2 (On the choice of $c$).

Similar to Remark 6.1, ideally the parameter $c$ should be chosen to be as
close as $\mathcal{H}(\mathbf{+1})$ so that $\delta$ is small enough, since in
both item (1) and item (2), the upper bounds there depend on $\delta$ through
$\mathcal{O}(e^{3\delta})$.

Similar to the previous example where we consider the ferromagnetic Ising
model on the hypercube, we observe that the original Glauber dynamics has an
exponentially large mixing time in the system size $n$, while the proposed
annealing algorithm mixes rapidly with its mixing time being at most
polynomial in $n$. As both the hypercube and the complete graph belong to the
class of expander graphs, the original critical height grows at least linearly
with $n$. The advantage of the proposed annealing algorithm is that the
reduced critical height is bounded independently of $n$, and consequently this
leads to rapid mixing.

####  6.2.1. Proof of Theorem 6.2

The proof is similar to the proof of Theorem 6.1. First, we prove item (1).
Making use of Theorem 4.3, the spectral gap lower bound Lemma 6.1 and the
inequality that

$\dfrac{Z_{1,\alpha,c}^{f}}{Z^{0}}\leqslant Z_{1,\alpha,c}^{f}\leqslant
2^{n},$

the desired result follows.

Next, we prove item (2), which readily follows from Theorem 5.3 and Lemma 6.1
that

$\mathcal{H}_{\max}-\mathcal{H}_{\min}=\mathcal{O}(n^{2}),\quad
K=\dfrac{2}{n}e^{-3\delta}.$

Finally, we prove item (3). Using (Levin and Peres, 2017, Lemma $20.11$), we
see that

$\displaystyle
t_{mix}^{0}(\varepsilon)\geqslant\dfrac{1}{\lambda_{2}(-M^{0})}\ln\left(\dfrac{1}{2\varepsilon}\right).$

According to Dommers et al. (2017), the critical height $m^{0}$ is bounded
below by

$m^{0}\geqslant\mathbb{M}n^{2},$

and the desired result follows from (Holley and Stroock, 1988, Lemma $2.3$)
since

$\lambda_{2}(-M^{0})=\mathcal{O}(e^{\mathbb{M}n^{2}}/4^{n})=\mathcal{O}(e^{\mathbb{M}n^{2}}).$

Note that in the above equations the value of $\mathbb{M}$ may possibly be
changed in each line.

###  6.3. $q$-state Potts model on the two-dimensional torus

The $q$-state Potts model is a generalization of the Ising model that we saw
in the previous two sections, which is a spin system characterized by a set
$S=\\{1,\ldots,q\\}$ of spin values on a finite graph $G=(V,E)$ with
$q\in\mathbb{N}$. In this Section, we shall specialize into the two-
dimensional torus $G=\left(\mathbb{Z}/n\mathbb{Z}\right)^{2}$. Each
configuration $\sigma\in\mathcal{X}=S^{V}$ assigns a spin value $\sigma(v)\in
S$ to each vertex $v\in V$, and the Hamiltonian is given by

(6.2) $\displaystyle\mathcal{H}(\sigma):=-J\sum_{(v,w)\in
E}\mathds{1}_{\\{\sigma(v)=\sigma(w)\\}},\quad\sigma\in\mathcal{X},$

where $J=\mathcal{O}(1)>0$ is the coupling constant. In subsequent analysis,
we shall be interested in monochromatic configurations. Denote $\mathbf{m}$ to
be the spin configuration in which $\sigma(v)=m$ for all $v\in V$, where $m\in
S$. Unlike the previous two examples of ferromagnetic Ising model where
$\mathbf{+1}$ is the only ground-state configuration, we note that the $q$
monochromatic configurations are the ground-state configurations in the
$q$-state Potts model.

For the Glauber dynamics, we consider a simple random walk proposal with
transition matrix
$P^{Potts}=(P^{Potts}(\sigma_{1},\sigma_{2}))_{\sigma_{1},\sigma_{2}\in\mathcal{X}}$
given by

$\displaystyle P^{Potts}(\sigma_{1},\sigma_{2})$
$\displaystyle:=\dfrac{1}{q|V|}\mathds{1}_{\\{|v\in
V;~{}\sigma_{1}(v)\neq\sigma_{2}(v)|=1\\}}.$

The stationary distribution of $P^{Potts}$ is the uniform distribution on
$\mathcal{X}$. We now state the main result of this Section:

######  Theorem 6.3.

Suppose $G=\left(\mathbb{Z}/n\mathbb{Z}\right)^{2}=(V,E)$, consider the
Glauber dynamics on $\mathcal{X}=\\{1,2,\ldots,q\\}^{V}$ with proposal
generator $P^{Potts}-I$, target distribution $\mu_{1,\alpha,c}^{f}$ with
Hamiltonian given in (6.2), $f(x)=x^{2}$ and
$c=\mathcal{H}(\mathbf{1})+\delta$. Given $\varepsilon>0$, we have

  1. (1)

[Rapid mixing: mixing time is at most polynomial in $n$ in the modified
Glauber dynamics with time-dependent penalty] Assume the logarithmic annealing
schedule $\alpha_{t}=\frac{3\pi}{2}\dfrac{1}{\ln(e^{\frac{3\pi}{2}}+pt)}$,
where $p$ is chosen as in Theorem 5.3, then

$t_{mix}^{f}(\varepsilon)=\mathcal{O}\left(n^{4}e^{3\delta}\left(\dfrac{1}{\varepsilon}\right)^{\frac{3\pi}{2n^{2}\ln
q}}\right).$

  2. (2)

[Torpid mixing: mixing time is at least exponential in $n$ in the original
Glauber dynamics] Assume in addition that $2\leqslant q\leqslant 4$. For
constant $\mathbb{M}=\mathbb{M}(q)$ that does not depend on $n$,

$t_{mix}^{0}(\varepsilon)=\Omega\left(e^{\mathbb{M}n}\ln\left(\dfrac{1}{2\varepsilon}\right)\right).$

######  Remark 6.3.

For $m\in\\{1,\ldots,q\\}$, we can alternatively take
$c=\mathcal{H}(\textbf{m})+\delta=\mathcal{H}(\textbf{1})+\delta$, where we
recall that there are $q$ ground-state configurations for the $q$-state Potts
model.

####  6.3.1. Proof of Theorem 6.3

We first state a spectral gap lower bound that generalizes Lemma 6.1:

######  Lemma 6.3.

Given $G=\left(\mathbb{Z}/n\mathbb{Z}\right)^{2}=(V,E)$, consider the Glauber
dynamics on $\mathcal{X}=\\{1,2,\ldots,q\\}^{V}$ with proposal generator
$P^{Potts}-I$ and target distribution $\mu_{1,\alpha,c}^{f}$, then

$\displaystyle\lambda_{2}(-M_{1,\alpha,c}^{f})\geqslant\dfrac{q}{(q-1)|V|}e^{-3\max_{\sigma}\mathcal{H}_{1,\alpha,c}^{f}(\sigma)}.$

In particular, if we take $f(x)=x^{2}$ and $c=\mathcal{H}(\mathbf{1})+\delta$,
where $\delta>0$, we have

$\displaystyle\lambda_{2}(-M_{1,\alpha,c}^{f})\geqslant\dfrac{q}{(q-1)|V|}e^{-3\delta-\frac{3}{\sqrt{\alpha}}\frac{\pi}{2}}.$

[Proof of Lemma 6.3. ] To prove the desired result, we again invoke Lemma 6.2:
the desired result follows if we take $L_{1}=M_{1,\alpha,c}^{f}$,
$L_{2}=P^{Potts}-I$ and
$A=\exp\bigg{\\{}-\max_{\sigma}\mathcal{H}_{1,\alpha,c}^{f}(\sigma)\bigg{\\}}$
therein, and recall that in (Levin and Peres, 2017, Corollary $12.13$) and
(Aldous and Fill, 2002, Chapter $5$ Example $9$) the spectral gap of
$P^{Potts}$ can be computed to be
$\lambda_{2}(-(P^{Potts}-I))=\frac{q}{(q-1)|V|}$.

In the special case when $f(x)=x^{2}$ and $c=\mathcal{H}(\mathbf{1})+\delta$,
we make use of Section 2.3 and $\sup_{x}\arctan(x)=\pi/2$.

We proceed to give the proof of Theorem 6.3.

[Proof of Theorem 6.3. ] First, we prove item (1), which simply follows from
Theorem 5.3 and the observation from Lemma 6.3 that

$\mathcal{H}_{\max}-\mathcal{H}_{\min}=\mathcal{O}(n^{2}),\quad
K=\dfrac{q}{(q-1)4n^{2}}e^{-3\delta}.$

Next, we prove item (2). Using (Levin and Peres, 2017, Lemma $20.11$), we see
that

$\displaystyle
t_{mix}^{0}(\varepsilon)\geqslant\dfrac{1}{\lambda_{2}(-M^{0})}\ln\left(\dfrac{1}{2\varepsilon}\right).$

The desired result follows from (Gheissari and Lubetzky, 2018, Theorem $4$)
since

$\lambda_{2}(-M^{0})=\mathcal{O}(e^{\mathbb{M}n}).$

## Acknowledgements

The first author Michael Choi acknowledges the financial support from the
startup grant and the special pocket research grant from Yale-NUS College,
National University of Singapore. The second author Jing Zhang acknowledges
the financial support from AIRS - Shenzhen Institute of Artificial
Intelligence and Robotics for Society Project 2019-INT002, NSFC Project number
12001460, and from SRIBD - Shenzhen Research Institute of Big Data.

## References

  * Aldous and Fill (2002) D. Aldous and J. A. Fill.  Reversible Markov Chains and Random Walks on Graphs, 2002.  Unfinished monograph, recompiled 2014, available at http://www.stat.berkeley.edu/~aldous/RWG/book.html.
  * Andrieu et al. (2008) C. Andrieu, A. Jasra, A. Doucet, and P. Del Moral.  A note on convergence of the equi-energy sampler.  _Stoch. Anal. Appl._ , 26(2):298–312, 2008.
  * Bhatnagar and Randall (2016) N. Bhatnagar and D. Randall.  Simulated tempering and swapping on mean-field models.  _J. Stat. Phys._ , 164(3):495–530, 2016.
  * Bierkens (2016) J. Bierkens.  Non-reversible Metropolis-Hastings.  _Stat. Comput._ , 26(6):1213–1228, 2016.
  * Bovier and den Hollander (2015) A. Bovier and F. den Hollander.  _Metastability_ , volume 351 of _Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]_.  Springer, Cham, 2015.  A potential-theoretic approach.
  * Brooks et al. (2011) S. Brooks, A. Gelman, G. L. Jones, and X.-L. Meng, editors.  _Handbook of Markov chain Monte Carlo_.  Chapman & Hall/CRC Handbooks of Modern Statistical Methods. CRC Press, Boca Raton, FL, 2011.
  * Choi (2020) M. C. Choi.  On the convergence of an improved and adaptive kinetic simulated annealing.  _arXiv preprint arXiv:2009.00195_ , 2020.
  * Choi (2021a) M. C. H. Choi.  On the convergence of an improved discrete simulated annealing via landscape modification, 2021a.
  * Choi (2021b) M. C. H. Choi.  An improved variant of simulated annealing that converges under fast cooling.  _Markov Process. Related Fields_ , 27(1):123–154, 2021b.
  * Deuschel and Mazza (1994) J.-D. Deuschel and C. Mazza.  $L^{2}$ convergence of time nonhomogeneous Markov processes. I. Spectral estimates.  _Ann. Appl. Probab._ , 4(4):1012–1056, 1994.
  * Diaconis (2009) P. Diaconis.  The Markov chain Monte Carlo revolution.  _Bull. Amer. Math. Soc. (N.S.)_ , 46(2):179–205, 2009.
  * Diaconis and Saloff-Coste (1996) P. Diaconis and L. Saloff-Coste.  Logarithmic Sobolev inequalities for finite Markov chains.  _Ann. Appl. Probab._ , 6(3):695–750, 1996.
  * Dommers et al. (2017) S. Dommers, F. den Hollander, O. Jovanovski, and F. R. Nardi.  Metastability for Glauber dynamics on random graphs.  _Ann. Appl. Probab._ , 27(4):2130–2158, 2017\.
  * Ebbers and Löwe (2020) M. Ebbers and M. Löwe.  Equi-energy sampling does not converge rapidly on the mean-field Potts model with three colors close to the critical temperature.  _J. Phys. A_ , 53(11):115001, 14, 2020.
  * Fang et al. (1997) H. Fang, M. Qian, and G. Gong.  An improved annealing method and its large-time behavior.  _Stochastic Process. Appl._ , 71(1):55–74, 1997\.
  * Gheissari and Lubetzky (2018) R. Gheissari and E. Lubetzky.  Mixing times of critical two-dimensional Potts models.  _Comm. Pure Appl. Math._ , 71(5):994–1046, 2018\.
  * Gidas (1985) B. Gidas.  Global optimization via the Langevin equation.  In _1985 24th IEEE Conference on Decision and Control_ , pages 774–778, Dec 1985.
  * Holley and Stroock (1987) R. Holley and D. Stroock.  Logarithmic Sobolev inequalities and stochastic Ising models.  _J. Statist. Phys._ , 46(5-6):1159–1194, 1987\.
  * Holley and Stroock (1988) R. Holley and D. Stroock.  Simulated annealing via Sobolev inequalities.  _Comm. Math. Phys._ , 115(4):553–569, 1988.
  * Hwang et al. (1993) C.-R. Hwang, S.-Y. Hwang-Ma, and S. J. Sheu.  Accelerating Gaussian diffusions.  _Ann. Appl. Probab._ , 3(3):897–913, 1993.
  * Hwang et al. (2005) C.-R. Hwang, S.-Y. Hwang-Ma, and S.-J. Sheu.  Accelerating diffusions.  _Ann. Appl. Probab._ , 15(2):1433–1444, 2005\.
  * Johndrow and Smith (2018) J. Johndrow and A. Smith.  Fast mixing of Metropolis-Hastings with unimodal targets.  _Electron. Commun. Probab._ , 23:Paper No. 71, 9, 2018.
  * Jovanovski (2017) O. Jovanovski.  Metastability for the Ising model on the hypercube.  _J. Stat. Phys._ , 167(1):135–159, 2017.
  * Karagulyan and Dalalyan (2020) A. G. Karagulyan and A. S. Dalalyan.  Penalized langevin dynamics with vanishing penalty for smooth and log-concave targets.  In _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020_ , 2020.
  * Kou et al. (2006) S. C. Kou, Q. Zhou, and W. H. Wong.  Equi-energy sampler with applications in statistical inference and statistical mechanics.  _Ann. Statist._ , 34(4):1581–1652, 2006.  With discussions and a rejoinder by the authors.
  * Levin and Peres (2017) D. A. Levin and Y. Peres.  _Markov chains and mixing times_.  American Mathematical Society, Providence, RI, 2017.
  * Levin et al. (2009) D. A. Levin, Y. Peres, and E. L. Wilmer.  _Markov chains and mixing times_.  American Mathematical Society, Providence, RI, 2009.
  * Lezaud (2001) P. Lezaud.  Chernoff and Berry-Esséen inequalities for Markov processes.  _ESAIM Probab. Stat._ , 5, 2001.
  * Madras and Zheng (2003) N. Madras and Z. Zheng.  On the swapping algorithm.  _Random Structures Algorithms_ , 22(1):66–97, 2003.
  * Miclo (1992) L. Miclo.  Recuit simulé sur ${\bf R}^{n}$. Étude de l’évolution de l’énergie libre.  _Ann. Inst. H. Poincaré Probab. Statist._ , 28(2):235–266, 1992.
  * Monmarché (2018) P. Monmarché.  Hypocoercivity in metastable settings and kinetic simulated annealing.  _Probab. Theory Related Fields_ , 172(3-4):1215–1248, 2018.
  * Roberts and Rosenthal (2004) G. O. Roberts and J. S. Rosenthal.  General state space Markov chains and MCMC algorithms.  _Probab. Surv._ , 1:20–71, 2004.
  * Woodard et al. (2009) D. B. Woodard, S. C. Schmidler, and M. Huber.  Conditions for rapid mixing of parallel and simulated tempering on multimodal distributions.  _Ann. Appl. Probab._ , 19(2):617–640, 2009.