# Almost sure convergence of the accelerated weight histogram algorithm

Henrik Hult and Guo-Jhen Wu

###### Abstract

The accelerated weight histogram (AWH) algorithm is an iterative extended
ensemble algorithm, developed for statistical physics and computational
biology applications. It is used to estimate free energy differences and
expectations with respect to Gibbs measures. The AWH algorithm is based on
iterative updates of a design parameter, which is closely related to the free
energy, obtained by matching a weight histogram with a specified target
distribution. The weight histogram is constructed from samples of a Markov
chain on the product of the state space and parameter space. In this paper
almost sure convergence of the AWH algorithm is proved, for estimating free
energy differences as well as estimating expectations with adaptive ergodic
averages. The proof is based on identifying the AWH algorithm as a stochastic
approximation and studying the properties of the associated limit ordinary
differential equation.

##  1 Introduction

Consider a parameterized family of Gibbs measures which have densities
$\\{p_{X|\Lambda}(\cdot|\lambda),\lambda\in\mathcal{L}\\}$ with respect to a
common reference measure $\nu$ on a state space $\mathcal{X}$ given by

$\displaystyle p_{X|\Lambda}(x|\lambda)=\exp\\{-E(x,\lambda)+F(\lambda)\\},$
(1.1)

where $E(x,\lambda)$ is the energy, $F(\lambda)$
$=-\log\int\exp\\{-E(x,\lambda)\\}\nu(dx)$ is the free energy (i.e.,
$e^{-F(\lambda)}$ is the normalizing constant) and $\mathcal{L}$ is the
parameter space. Such models originate from statistical physics, but are
frequently encountered in various fields including Bayesian statistics,
artificial intelligence and machine learning, computational biology, etc. The
interests lie in computing free energy differences, $F(\lambda)-F(\sigma)$
with $\lambda,\sigma\in\mathcal{L}$, and expectations with respect to
$p_{X|\Lambda}(\cdot|\lambda)$ at a fixed parameter $\lambda\in\mathcal{L}$.
In classical models of statistical physics, $\mathcal{X}$ represents the state
space of a system of particles and the energy is commonly of the form
$E(x,\lambda)=\lambda E(x)$, with $\lambda$ as the inverse temperature. More
generally, the parameter $\lambda$ may represent a number of reaction
coordinates. In statistical machine learning models, the parameter space
$\mathcal{L}$ may be high-dimensional, consisting of millions of parameters.

Sampling configurations according to a Gibbs measure is challenging because
the free energy is unknown; the state space is often too large, which makes
summation/integration over the state space infeasible. A widely used approach
is Markov chain Monte Carlo (MCMC), where an ergodic Markov chain is
constructed with $p_{X|\Lambda}(\cdot|\lambda)$ as its invariant distribution,
see [14, 2, 27, 1]. The vanilla MCMC algorithms often suffer from slow mixing
(slow relaxation) that reduces the effective number of samples and even may
lead to biased results. In physical problems there are several potential
reasons that lead to slow mixing, including critical slow down near points of
second order phase transitions, nucleation associated with first order phase
transitions and trapping in metastable states around local minima in complex
energy landscapes. The latter problem is often encountered in models of spin
glasses [29, 6], interacting spins in random fields [24], and heteropolymers
[30, 15] as well as in complex statistical inference problems where models are
highly non-Gaussian, see e.g. [14, 27, 28].

Extended ensemble methods refer to computational methods where multiple copies
of the state and parameter spaces are considered and a Markov chain on the
extended space is constructed. Such methods have been developed since the mid
80’s to overcome the problems of slow mixing due to nucleation and
metastability. They include well known algorithms such as parallel
tempering/replica exchange, see [29], and related exchange Monte Carlo
methods, simulated tempering, see [24], multicanonical Monte Carlo (adaptive
umbrella sampling), see [25, 16, 7, 8, 5], as well as the more recent infinite
swapping algorithm, see [9], and the accelerated weight histogram (AWH)
algorithm introduced by Lidmar in [21]. More recent advances also include
generative models from machine learning that can be trained to approximate
samples from Gibbs distributions, see [26].

The main idea of extended ensemble algorithms is to sample from artificial
ensembles that are constructed as extensions or compositions of the original
ensemble. Fast mixing of the Markov chain in higher temperature (energy, etc.)
components of the artificial ensembles greatly facilitate the mixing in other
components. Averages over the original ensemble are calculated by
marginalization (e.g. parallel tempering), conditional sampling (e.g.
simulated tempering), or, a reweighting procedure (e.g. multicanonical Monte
Carlo, infinite swapping). Improved mixing is not the only benefit of extended
ensemble methods. They are particularly useful for estimating free energy
differences $F(\lambda)-F(\sigma)$ for $\lambda,\sigma\in\mathcal{L}$ and
expectations $\int\psi(x)p_{X|\Lambda}(x|\lambda)\nu(dx)$ for a given function
$\psi$ and parameter $\lambda\in\mathcal{L}$. They are also useful for
efficient sampling of rare events, and may be used to estimate the probability
of configurations with small probability or transition probabilities between
different metastable states.

The design of extended ensemble methods is intricate. There is a lot of
flexibility in selecting the extended variables, the target density on the
extended space, the transition kernel of the Markov chain as well as
particular design parameters such as temperature selection and spacings in the
parameter space. The current practice is to a large extent based on heuristics
and experimental success. Methods based on the spectral gap to study
convergence of ergodic Markov chains, as well as analysis of the asymptotic
variance of particular observables, are generally insufficient to fully
capture the effect of design parameters. Recent advances in the theory of
large deviations for empirical measures of Markov chains shows promising
results for studying the rate of convergence of weighted empirical averages
towards population averages in the context of parallel tempering and infinite
swapping, see [11, 12].

Although large deviations theory may provide useful information for algorithm
design, it is very challenging to make detailed calculations so that they
provide optimal designs in complex models. For example, the problem of optimal
temperature selection in an infinite swapping algorithm for the two-well
potential has recently been worked out by Dupuis and Wu, see [12], and even
such relatively simple model requires extensive and ingenious analysis. For a
general multi-well potential, in order to make the calculations and analysis
feasible, they have to comprise with a sub-optimal (but nearly optimal)
design. It is our belief that including a learning mechanism is necessary in
complex models, where design parameters may be updated adaptively during the
simulation.

An example of this learning mechanism is provided in the accelerated weight
histogram method, see [21], where an estimate of the free energy is updated
iteratively during the simulation and the invariant density of the extended
Markov chain is updated accordingly. Ultimately, this iterative procedure
results in a fixed point method that can be analyzed using the theory of
stochastic approximation (Robbins-Monro).

The AWH algorithm has connections to the Wang-Landau algorithm, see [31, 4],
which can be viewed as an adaptive importance sampling algorithm belonging to
the class of free energy biasing techniques, see e.g. [19]. The Wang-Landau
method was originally developed for simulations in the closely related
multicanonical ensemble, see [32], but is straightforward to adapt to extended
ensemble simulations, see [8]. In the Wang-Landau algorithm the biasing factor
is updated during the simulation to flatten the target probability, but the
design is delicate since the adaption mechanism must decrease slowly with the
number of simulations, as first observed in [4]. By representing the Wang-
Landau algorithm as a stochastic approximation, the almost sure convergence of
the Wang-Landau algorithm is proved in [13], using the general conditions
provided in [3].

In contrast to the Wang-Landau algorithm, the AWH algorithm allows for large
parameter steps by the use of a Gibbs sampler combined with a reweighting
procedure that takes advantage of the information collected during the
simulation. This allows for a rather densely spaced set of parameters, without
being limited by slow diffusion. Moreover, the free energy parameters are
updated based on a histogram of weights, rather than a histogram of visited
parameter values, which is combined with the information collected during
previous iterations.

In this paper the almost sure convergence of the AWH algorithm is proved by
identifying it as a stochastic approximation algorithm. Free energy
differences are considered jointly with ergodic averages for fixed parameters.
The stability issue is circumvented by assuming that the iterates take values
in a compact set and almost sure convergence is proved by verifying general
conditions for stochastic approximations with state-dependent noise in [18,
Chapter 6]. The main technical difficulties lie in the identification of the
limit ordinary differential equation (ODE) for the stochastic approximation
algorithm, and the construction of appropriate Lyapunov functions for the
limit ODE to characterize its limit set. Extensions to non-projected
algorithms may also be treated, using the conditions in [3].

The paper is organized as follows. In Section 2 the AWH algorithm is presented
in detail, explaining its use for estimating free energy differences and
ergodic averages. Section 3 contains the results on almost sure convergence of
estimates of free energy differences and Section 4 extends the results to
include the joint estimation of free energy differences and ergodic averages
at fixed parameters. The Appendix contains existing results on almost sure
convergence of stochastic approximations that are used to prove the main
results.

##  2 The Accelerated Weight Histogram Algorithm

In this section the AWH algorithm is presented in detail. It was first
introduced in [20] to estimate the free energy differences and further
developed for applications in computational biology in [21]. The design of an
appropriate target distribution for the AWH algorithm is analyzed in [22].

Let $\mathcal{X}$ be the configuration space and $\mathcal{L}$ be the
parameter space. For simplicity, it is assumed that $\mathcal{{X}}$ and
$\mathcal{L}$ are both finite sets. We are concerned with the setting where
$|\mathcal{X}|$ is much larger than $|\mathcal{L}|$ and summation over
$\mathcal{X}$ is infeasible, whereas summation over $\mathcal{L}$ is feasible.
Consider a family of Gibbs measures
$\\{p_{X|\Lambda}(\cdot|\lambda),\lambda\in\mathcal{L}\\}$ on $\mathcal{X}$,
where

$\displaystyle
p_{X|\Lambda}(x|\lambda)\doteq\exp\\{-E(x,\lambda)+F(\lambda)\\},$

$E(x,\lambda)$ is the energy and
$F(\lambda)=-\log\sum_{x\in\mathcal{X}}\exp\\{-E(x,\lambda)\\}$ is the unknown
free energy. In the accelerated weight histogram algorithm a joint
distribution on the extended state space $\mathcal{X}\times\mathcal{L}$ is
constructed as

$\displaystyle
p_{X,\Lambda|\Theta}(x,\lambda|\theta)\propto\exp\\{-E(x,\lambda)+\theta(\lambda)\\},$
(2.1)

where $\theta:\mathcal{L}\to\mathbb{R}$ represents a design parameter, which
ideally is equal to $F$ up to an additive constant. The joint distribution is
obtained by specifying the marginal distribution on $\mathcal{L}$ as

$\displaystyle
p_{\Lambda|\Theta}(\lambda|\theta)=\frac{e^{\theta(\lambda)-F(\lambda)}}{\sum_{\sigma\in\mathcal{L}}e^{\theta(\sigma)-F(\sigma)}}.$
(2.2)

Consequently, by Bayes’ theorem, the conditional distribution of $\Lambda$
given $X$ and $\Theta$ is formulated by

$\displaystyle
p_{\Lambda|X,\Theta}(\lambda|x,\theta)=\frac{e^{-E(x,\lambda)+\theta(\lambda)}}{\sum_{\sigma\in\mathcal{L}}e^{-E(x,\sigma)+\theta(\sigma)}}.$

Note that, by construction
$p_{X|\Lambda,\Theta}(x|\lambda,\theta)=p_{X|\Lambda}(x|\lambda)$ since

$p_{X|\Lambda,\Theta}(x|\lambda,\theta)=\frac{p_{X,\Lambda|\Theta}(x|\lambda,\theta)}{p_{\Lambda|\Theta}(\lambda|\theta)}=\exp\\{-E(x,\lambda)+F(\lambda)\\}=p_{X|\Lambda}(x|\lambda).$

Additionally, throughout the paper, $\theta$ is called the design parameter,
which can be thought of as an approximation of the free energy, not to be
confused with $\lambda$, which is a parameter that of the energy function and
denotes a point in the parameter space $\mathcal{L}$.

Let $\rho$ be a probability measure on $\mathcal{L}$, called the target
distribution, which is specified by the user. Throughout the paper, we only
consider $\rho$ which satisfies $\rho(\lambda)>0$ for all
$\lambda\in\mathcal{L}$. The aim of the AWH algorithm is to find the design
parameter $\theta_{*}$ that makes

$\displaystyle
p_{\Lambda|\Theta}(\lambda|\theta_{*})=\rho(\lambda),\quad\lambda\in\mathcal{L}.$

We call such $\theta_{*}$ the optimal design parameter. With such
$\theta_{*}$, free energy differences may be computed by

$\displaystyle
F(\lambda)-F(\sigma)=\theta_{*}(\lambda)-\theta_{*}(\sigma)-\log\rho(\lambda)+\log\rho(\sigma)$
(2.3)

since

$\frac{\rho(\lambda)}{\rho(\sigma)}=\frac{p_{\Lambda|\Theta}(\lambda|\theta_{*})}{p_{\Lambda|\Theta}(\sigma|\theta_{*})}=\frac{e^{\theta_{*}(\lambda)-F(\lambda)}}{e^{\theta_{*}(\sigma)-F(\sigma)}}$

for any $\lambda,\sigma\in\mathcal{L}$. However, computing $\theta_{*}$ is
equally difficult as computing the free energy, $F$.

In the AWH algorithm estimates of free energy differences are obtained from
estimates of $\theta_{*}$. An iterative scheme is introduced where an estimate
$\theta_{n}$ of $\theta_{*}$ is updated in each iteration using weight
histograms, as described in the following. Let
$W_{n,k}:\mathcal{L}\to[0,\infty)$, $n=0,\dots,N-1$, $k=1,\dots,N_{I}$ be
weight histograms with $W_{0,0}(\lambda)=0$, $\lambda\in\mathcal{L}$. The
algorithm is initiated with an initial estimate $\theta_{0}$ of the design
parameter. A Markov chain $\\{(X_{k},\Lambda_{k}),k\geq 0\\}$ on the extended
state space $\mathcal{X}\times\mathcal{L}$ is constructed to have
$p_{X,\Lambda|\theta_{0}}$ as its unique invariant distribution. The Markov
chain is simulated for $N_{I}$ time steps to produce updated weight histograms
on the parameter space $\Lambda$ as follows:

$\displaystyle
W_{0,k}(\lambda)=W_{0,k-1}(\lambda)+p_{\Lambda|X,\Theta}(\lambda|X_{k},\theta_{0}),\quad\lambda\in\mathcal{L},\;k=1,\dots,N_{I}.$

An updated estimate $\theta_{1}$ of the optimal design parameter $\theta_{*}$
is obtained by matching the normalized weight histogram $W_{0,N_{I}}/N_{I}$
with the target distribution $\rho$,

$\displaystyle\theta_{1}(\lambda)=\theta_{0}(\lambda)-\log\left(\frac{W_{0,N_{I}}(\lambda)}{N_{I}\rho(\lambda)}\right),\quad\lambda\in\mathcal{L}.$

The procedure is iterated until a termination criteria is fulfilled; in this
paper the algorithm terminates after $N$ updates of the estimate of the
optimal design parameter. A pseudo-code for the AWH algorithm is provided in
Algorithm 1 and additional details provided below.

Result: $\theta_{N}$

Initialize: $n=0$, $X_{0}$, $\Lambda_{0}$, $\theta_{0}$, $\rho$, $W_{0,0}=0$ ;

while  _$n <N$_ do

for  _$k=1,\dots,N_{I}$_ do

Sample $X_{nN_{I}+k}$ from $q(X_{nN_{I}+k-1},x|\Lambda_{nN_{I}+k-1})$ ;

Sample $\Lambda_{nN_{I}+k}$ from
$p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})$;

Update $W_{n,k}(\lambda)\leftarrow
W_{n,k-1}(\lambda)+p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})$,
$\lambda\in\mathcal{L}$;

end for

Update
$\theta_{n+1}(\lambda)\leftarrow\theta_{n}(\lambda)-\log\left(\frac{W_{n,N_{I}}(\lambda)}{(n+1)N_{I}\rho(\lambda)}\right)$,
$\lambda\in\mathcal{L}$ ;

Set $W_{n+1,0}(\lambda)\leftarrow(n+1)N_{I}\rho(\lambda)$,
$\lambda\in\mathcal{L}$ ;

Update $n\leftarrow n+1$;

end while

Algorithm 1 Accelerated weight histogram

In the $n$th iteration the sampling steps aim to sample the Markov chain
$\\{(X_{nN_{I}+k},\Lambda_{nN_{I}+k}),k\geq 1\\}$ with invariant distribution
$p_{X,\Lambda|\Theta}(x,\lambda|\theta_{n})$. The two sampling steps in
Algorithm 1 form a Gibbs sampler with the desired invariant distribution. Note
that, in the first sampling step, one would ideally sample $X_{nN_{I}+k}$ from
$p_{X|\Lambda}(x|\Lambda_{nN_{I}+k-1})$, but this is difficult since the free
energy $F$ is unknown. Therefore the ideal step is replaced by sampling from a
transition probability $q(X_{nN_{I}+k-1},x|\Lambda_{nN_{I}+k-1})$ whose
invariant distribution is $p_{X|\Lambda}(x|\Lambda_{nN_{I}+k-1})$. That is,
$q(x,y|\lambda)$ satisfies

$\displaystyle
p_{X|\Lambda}(y|\lambda)=\sum_{x\in\mathcal{X}}q(x,y|\lambda)p_{X|\Lambda}(x|\lambda),\quad
y\in\mathcal{X},\lambda\in\mathcal{L}.$ (2.4)

Moreover, we assume that for any $\lambda\in\mathcal{L}$,

$\displaystyle\\{q(x,y|\lambda)\\}_{x,y\in\mathcal{X}},\;\text{is irreducible
and aperiodic.}$ (2.5)

This assumption will come in handy later on in the proof of almost sure
convergence, and the assumption can be fulfilled easily in an implementation.

In contrast, since we only consider the case when $|\mathcal{L}|$ is not too
large, the normalizing constant
$\sum_{\sigma\in\mathcal{L}}e^{-E(x,\sigma)+\theta_{n}(\sigma)}$ can be
explicitly computed in this case, which means that sampling
$\Lambda_{nN_{I}+k}$ from
$p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})$ is straightforward.

###  2.1 Estimation of free energy differences

For $n\geq 1$, the weight histogram can be written as

$\displaystyle
W_{n,k}(\lambda)=W_{n,0}(\lambda)+\sum_{j=1}^{k}p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+j},\theta_{n}),\quad\lambda\in\mathcal{L},k=1\dots,N_{I}.$

For large $N_{I}$ the normalized weight histogram approximates a convex
combination of the target density $\rho(\lambda)$ and the marginal density
$p_{\Lambda|\Theta}(\lambda|\theta_{n})$. Indeed,

$\displaystyle\frac{1}{(n+1)N_{I}}W_{n,N_{I}}(\lambda)$
$\displaystyle=\frac{n}{n+1}\rho(\lambda)+\frac{1}{(n+1)N_{I}}\sum_{k=1}^{N_{I}}p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})$
$\displaystyle\rightarrow\frac{n}{n+1}\rho(\lambda)+\frac{1}{n+1}\sum_{x\in\mathcal{X}}p_{\Lambda|X,\Theta}(\lambda|x,\theta_{n})p_{X|\Theta}(x|\theta_{n})$
$\displaystyle=\frac{n}{n+1}\rho(\lambda)+\frac{1}{n+1}p_{\Lambda|\Theta}(\lambda|\theta_{n}),\quad\lambda\in\mathcal{L},$

where the convergence is almost surely, as $N_{I}\rightarrow\infty$ , by the
ergodic theorem. The updating scheme of the design parameter is then
approximately

$\displaystyle\theta_{n+1}(\lambda)$
$\displaystyle=\theta_{n}(\lambda)-\log\left(\frac{W_{n,N_{I}}(\lambda)}{(n+1)N_{I}\rho(\lambda)}\right)$
$\displaystyle\approx\theta_{n}(\lambda)-\log\left(1-\frac{1}{n+1}\left(1-\frac{p_{\Lambda|\Theta}(\lambda|\theta_{n})}{\rho(\lambda)}\right)\right),\quad\lambda\in\mathcal{L},$

which has a fixed point at $\theta_{*}$ such that
$p_{\Lambda|\Theta}(\lambda|\theta_{*})=\rho(\lambda)$. Thus, $\theta_{n}$
computed by the AWH algorithm can be identified as a stochastic approximation
of the optimal design parameter $\theta_{*}$.

After termination of the AWH algorithm, e.g., after $N$ iterations, estimates
of free energy differences $F(\lambda)-F(\sigma)$ may be computed as

$\displaystyle\theta_{N}(\lambda)-\theta_{N}(\sigma)-\log\rho(\lambda)+\log\rho(\sigma),\quad\lambda,\sigma\in\mathcal{L},$

since (2.3) and $\theta_{N}$ is an estimate of $\theta_{*}$. This is one of
the main use of the AWH algorithm in [20].

###  2.2 Adaptive estimation of expectations with ergodic averages

In addition to estimating free energy differences, the AWH algorithm may also
be used to estimate expectations at a fixed parameter value
$\lambda\in\mathcal{L}$,

$\displaystyle\sum_{x\in\mathcal{X}}\psi(x)p_{X|\Lambda}(x|\lambda),$

for some function $\psi:\mathcal{X}\rightarrow\mathbb{R}$. At the $n$th
iteration, given design parameter $\theta_{n}$, the Markov chain
$\\{(X_{nN_{I}+k},\Lambda_{nN_{I}+k}),k=1,\dots,N_{I}\\}$ has invariant
distribution $p_{X,\Lambda|\Theta}(x,\lambda|\theta_{n})$ given by (2.1). By
the ergodic theorem, as $k\to\infty$,

$\displaystyle\frac{1}{k}\sum_{j=1}^{k}\psi(X_{nN_{I}+j})p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+j},\theta_{n})$
$\displaystyle\rightarrow\sum_{x\in\mathcal{X}}\psi(x)p_{\Lambda|X,\Theta}(\lambda|x,\theta_{n})p_{X|\Theta}(x|\theta_{n})$
$\displaystyle=\sum_{x\in\mathcal{X}}\psi(x)p_{X,\Lambda|\Theta}(x,\lambda|\theta_{n}),$

and

$\displaystyle\frac{1}{k}\sum_{j=1}^{k}p_{X|\Lambda,\Theta}(\lambda|X_{nN_{I}+j},\theta_{n})$
$\displaystyle\rightarrow\sum_{x\in\mathcal{X}}p_{\Lambda|X,\Theta}(\lambda|x,\theta_{n})p_{X|\Theta}(x|\theta_{n})=p_{\Lambda|\Theta}(\lambda|\theta_{n}).$

These two limits imply that, as $k\to\infty$,

$\displaystyle\frac{\sum_{j=1}^{k}\psi(X_{nN_{I}+j})p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+j},\theta_{n})}{\sum_{j=1}^{k}p_{X|\Lambda,\Theta}(\lambda|X_{nN_{I}+j},\theta_{n})}$
$\displaystyle\rightarrow\frac{\sum_{x\in\mathcal{X}}\psi(x)p_{X,\Lambda|\Theta}(x,\lambda|\theta_{n})}{p_{\Lambda|\Theta}(\lambda|\theta_{n})}$
$\displaystyle=\sum_{x\in\mathcal{X}}\psi(x)p_{X|\Lambda}(x|\lambda).$

Thus, it is in principle possible to estimate the expectation on the right
side in the last display by the left side in the last display for some large
$k\in\mathbb{N}$. Such estimate is inefficient since it would require fixing
$\theta_{n}$ and let the Markov chain run for a long time without updating
$\theta_{n}$, not taking advantage of the possible improved estimates of the
optimal design parameter along the simulation. A more efficient estimate is to
keep updating $\theta_{n}$ according to the AWH algorithm and use the estimate

$\displaystyle\frac{\sum_{n=0}^{N-1}\sum_{j=1}^{N_{I}}\psi(X_{nN_{I}+j})p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+j},\theta_{n})}{\sum_{n=0}^{N-1}\sum_{j=1}^{N_{I}}p_{X|\Lambda,\Theta}(\lambda|X_{nN_{I}+j},\theta_{n})}.$
(2.6)

As $\theta_{n}$ converges to the optimal $\theta_{*}$, it seems plausible that
$\theta_{n}$ will essentially remain constant for large values of $n$ and the
estimate in (2.6) will converge to
$\sum_{x\in\mathcal{X}}\psi(x)p_{X|\Lambda}(x|\lambda)$ almost surely. This is
indeed the case, as will be demonstrated in Section 4. The proof is based on
studying the convergence an extended stochastic approximation of $\theta_{*}$
and integrals of the form
$\sum_{x\in\mathcal{X}}\psi(x)p_{X,\Lambda|\Theta}(x,\lambda|\theta_{*})$.

##  3 Almost Sure Convergence of Free Energy Differences

In this section almost sure convergence of the AWH algorithm for estimating
free energy differences will be proved. The proof is based on the ODE approach
for proving almost sure convergence of projected stochastic approximation
algorithms, as presented in [18] and explained in more detail in the Appendix.
In this approach, originally introduced by Ljung [23], it is demonstrated that
the noise averages out so that, asymptotically, the stochastic approximation
algorithm may be interpreted as a discretization of a limit ODE. Consequently,
the set of limit points of the stochastic approximation algorithm is a subset
of the limit points of the limit ODE, which may be identified as a subset of
the level set of $0$ of a properly constructed Lyapunov function.

Our proof relies on identifying the AWH algorithm as a stochastic
approximation algorithm, identification of the associated limit ODE, the
construction of an associated Lyapunov function, and verifying the appropriate
set of conditions for almost sure convergence.

Let us begin by identifying the AWH algorithm as a stochastic approximation.
Recall that the recursion in the AWH algorithm is given by

$\displaystyle\theta_{n+1}(\lambda)$
$\displaystyle=\theta_{n}(\lambda)-\log\left(\frac{W_{nN_{I}}(\lambda)}{(n+1)N_{I}\rho(\lambda)}\right)$
$\displaystyle=\theta_{n}(\lambda)-\log\left(1-\frac{1}{n+1}+\frac{1}{n+1}\frac{\sum_{k=1}^{N_{I}}p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})}{N_{I}\rho(\lambda)}\right).$

Since $\log(1+x)=x+o(|x|)$ the recursion can be written as

$\displaystyle\theta_{n+1}(\lambda)=\theta_{n}(\lambda)+\frac{1}{n+1}\left(1-\frac{\sum_{k=1}^{N_{I}}p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})}{N_{I}\rho(\lambda)}+\beta_{n}(\lambda)\right),$

$\lambda\in\mathcal{L}$, where $\beta_{n}(\lambda)=o(1)$ almost surely. To
guarantee realistic estimates it is practical to project the updates on a
hyper-rectangle $H$ containing the origin. The projected algorithm circumvents
issues of stability and guarantees that estimates do not explode. The
recursion then becomes,

$\displaystyle\theta_{n+1}(\lambda)\\!=\\!\theta_{n}(\lambda)\\!+\\!\frac{1}{n+1}\left(\\!1-\frac{\sum_{k=1}^{N_{I}}p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})}{N_{I}\rho(\lambda)}\\!+\\!\beta_{n}(\lambda)\\!+\\!Z_{n}(\lambda)\\!\right),$
(3.1)

where $Z_{n}(\lambda)$ denotes the projection term due to the constraint to
$H$. The exact definition of the projection term can be found in [18, Section
4.3]. To sample $X_{nN_{I}+k}$, $k=1,\dots,N_{I}$, a Gibbs sampler is used by
alternating between sampling $X_{nN_{I}+k}$ from
$q(X_{nN_{I}+k-1},x|\Lambda_{nN_{I}+k-1})$ and $\Lambda_{nN_{I}+k}$ from
$p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})$, where
$q(x,y|\lambda)$ has the desired invariant distribution, that is,

$\displaystyle
p_{X|\Lambda}(y|\lambda)=\sum_{x\in\mathcal{X}}q(x,y|\lambda)p_{X|\Lambda}(x|\lambda),\quad
y\in\mathcal{X},\lambda\in\mathcal{L}.$

For $n\geq 1$, let $\theta_{n}=(\theta_{n}(\lambda):\lambda\in\mathcal{L})$,
$\varepsilon_{n}=1/(n+1)$

$\xi_{n}=(X_{(n-1)N_{I}+1},\dots,X_{nN_{I}},\Lambda_{(n-1)N_{I}+1},\dots,\Lambda_{nN_{I}}),$

and
$h:\mathbb{R}^{|\mathcal{L}|}\times(\mathcal{X}^{N_{I}}\times\mathcal{L}^{N_{I}})\to\mathbb{R}^{|\mathcal{L}|}$
be given by

$\displaystyle
h(\theta,\xi)(\lambda)=1-\frac{\sum_{k=1}^{N_{I}}p_{\Lambda|X,\Theta}(\lambda|\xi({k}),\theta)}{N_{I}\rho(\lambda)},$

where $\xi(k)$ denotes the $k$th component of the vector $\xi$. Notice that
$\xi_{n}(k)=X_{(n-1)N_{I}+k}$ and $\xi_{n}(N_{I}+k)=\Lambda_{(n-1)N_{I}+k}$
for $k=1,\dots,N_{I}$. The recursion (3.1) can then be stated as

$\displaystyle\theta_{n+1}$
$\displaystyle=\theta_{n}+\varepsilon_{n}\left(h(\theta_{n},\xi_{n+1})+\beta_{n}+Z_{n}\right)$
$\displaystyle=\theta_{n}+\varepsilon_{n}Y_{n}+\varepsilon_{n}Z_{n},$

where $Y_{n}=h(\theta_{n},\xi_{n+1})+\beta_{n}$. Let
$\mathcal{F}_{n}=\sigma(\theta_{0},Y_{0},\dots,Y_{n-1},\xi_{1},\dots,\xi_{n})$.
Then, $P(\xi_{n+1}\in\cdot|\xi_{i},\theta_{i},i\leq
n)=p(\xi_{n},\cdot|\theta_{n})$, where $p(\xi,\cdot|\theta)$ denotes the one-
step transition probability with starting point $\xi$, parameterized by
$\theta$, given by

$\displaystyle p(\xi,\eta|\theta)$
$\displaystyle=\prod_{k=1}^{N_{I}}q\left(\eta(k-1),\eta(k)|\hat{\eta}(N_{I}+k-1)\right)p_{\Lambda|X,\Theta}\left(\eta(N_{I}+k)|\eta(k),\theta\right),$
(3.2)

where $\eta(0)=\xi(N_{I})$, $\hat{\eta}(N_{I})=\xi(2N_{I})$ and
$\hat{\eta}(k)=\eta(k)$ for all $k\neq N_{I}$. Furthermore,

$\displaystyle
E_{n}[h(\theta_{n},\xi_{n+1})]=E[h(\theta_{n},\xi_{n+1})|\mathcal{F}_{n}]=\sum_{\eta}h(\theta_{n},\eta)p(\xi_{n},\eta|\theta_{n})=g(\theta_{n},\xi_{n}),$

where $E_{n}$ denotes the expectation conditioned on $\mathcal{F}_{n}$ and

$g(\theta,\xi)\doteq\sum_{\eta}h(\theta,\eta)p(\xi,\eta|\theta).$

Since $\rho(\lambda)>0$ for all $\lambda\in\mathcal{L}$, we find that $h$ is
bounded which implies $E|Y_{n}|<\infty$ and

$E_{n}Y_{n}=g(\theta_{n},\xi_{n})+\beta_{n}.$

For convenience and because $E_{n}\beta_{n}=o(1)$, here we use $\beta_{n}$ to
denote $E_{n}\beta_{n}$.

By introducing

$\displaystyle\delta M_{n}$
$\displaystyle=Y_{n}-E_{n}Y_{n}=h(\theta_{n},\xi_{n+1})-g(\theta_{n},\xi_{n})$
$\displaystyle=h(\theta_{n},\xi_{n+1})-\sum_{\eta}h(\theta_{n},\eta)p(\xi_{n},\eta|\theta_{n})$
$\displaystyle=\sum_{\eta}[h(\theta_{n},\xi_{n+1})-h(\theta_{n},\eta)]p(\xi_{n},\eta|\theta_{n}),$

the recursion (3.1) may be restated as,

$\theta_{n+1}=\theta_{n}+\varepsilon_{n}[E_{n}Y_{n}+\delta M_{n}+Z_{n}],$

or equivalently,

$\displaystyle\theta_{n+1}=\theta_{n}+\varepsilon_{n}[g(\theta_{n},\xi_{n})+\delta
M_{n}+\beta_{n}+Z_{n}],$ (3.3)

which is a standard form of a stochastic approximation algorithm with state-
dependent noise, see [18, Section 6.6].

###  3.1 Identification of the limit ODE

The limit ODE is a projected ODE of the form $\dot{\theta}=\bar{g}(\theta)+z$,
where $z$ is a projection term, that approximates the asymptotic behaviour of
the stochastic recursion (3.3). To identify $\bar{g}$ a first step is to
assume that the noise $\delta M_{n}$ and bias $\beta_{n}$ vanishes
asymptotically. Moreover, since
$g(\theta,\xi)=\sum_{\eta}h(\theta,\eta)p(\xi,\eta|\theta)$ and $\theta_{n}$
varies slowly with $n$, the intuition is that, for large $n$ and
$\theta_{n}\approx\theta$, the Markov chain $\xi_{n}$ is approximately
distributed according to the invariant distribution of the transition
probability $p(\xi,\cdot\,|\theta)$. Therefore, a candidate for $\bar{g}$ is
$E[h(\theta,\xi_{1})(\lambda)|\theta_{0}=\theta]$, where
$\xi_{1}=(X_{1},\dots,X_{N_{I}},\Lambda_{1},\dots,\Lambda_{N_{I}})$ depends on
$\theta_{0}$ and the expectation is taken under
$X_{0}\overset{d}{=}p_{X|\Theta}(\cdot|\theta)$.

Let $\bar{g}:\mathbb{R}^{|\mathcal{L}|}\to\mathbb{R}^{|\mathcal{L}|}$ be given
by

$\displaystyle\bar{g}(\theta)(\lambda)$
$\displaystyle=1-\frac{p_{\Lambda|\Theta}(\lambda|\theta)}{\rho(\lambda)}.$
(3.4)

The expression for $\bar{g}$ is motivated by the following lemma that
expresses $\bar{g}(\theta)$ as the ergodic average of $h(\theta,\xi_{n})$, for
fixed $\theta$.

###### Lemma 3.1

Let $\bar{g}$ be given by (3.4). For $\theta\in\mathbb{R}^{|\mathcal{L}]}$ and
$\lambda\in\mathcal{L}$,
$\bar{g}(\theta)(\lambda)=E[h(\theta,\xi_{1})(\lambda)|\theta_{0}=\theta]$.

Proof. For each $\lambda\in\mathcal{L}$, the expectation on the right hand
side can be written as

$\displaystyle E[h(\theta,\xi_{1})(\lambda)|\theta_{0}=\theta]$
$\displaystyle=\sum_{x_{0}\in\mathcal{X}}E[h(\theta,(X_{1},\dots,X_{N_{I}},\Lambda_{1},\dots,\Lambda_{N_{I}}))(\lambda)|\theta_{0}=\theta,X_{0}=x_{0}]p_{X|\Theta}(x_{0}|\theta)$
$\displaystyle=\sum_{x_{0}\in\mathcal{X}}\sum_{\lambda_{0}\in\mathcal{L}}E[h(\theta,(X_{1},\dots,X_{N_{I}},\Lambda_{1},\dots,\Lambda_{N_{I}}))(\lambda)|\theta_{0}=\theta,X_{0}=x_{0}]$
$\displaystyle\qquad\times
p_{X|\Lambda}(x_{0}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta)$
$\displaystyle=\sum_{x_{0}\in\mathcal{X}}\cdots\sum_{x_{N_{I}}\in\mathcal{X}}\sum_{\lambda_{0}\in\mathcal{L}}\cdots\sum_{\lambda_{N_{I}}\in\mathcal{L}}h(\theta,(x_{1},\dots,x_{N_{I}},\lambda_{1},\dots,\lambda_{N_{I}}))(\lambda)$
$\displaystyle\qquad\times\left(\prod_{k=1}^{N_{I}}q(x_{k-1},x_{k}|\lambda_{k-1})p_{\Lambda|X,\Theta}(\lambda_{k}|x_{k},\theta)\right)p_{X|\Lambda}(x_{0}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta)$
$\displaystyle=\sum_{x_{0}\in\mathcal{X}}\cdots\sum_{x_{N_{I}}\in\mathcal{X}}\sum_{\lambda_{0}\in\mathcal{L}}\cdots\sum_{\lambda_{N_{I}}\in\mathcal{L}}\left(1-\frac{\sum_{k=1}^{N_{I}}p_{\Lambda|X,\Theta}(\lambda|x_{k},\theta)}{N_{I}\rho(\lambda)}\right)$
$\displaystyle\qquad\times\left(\prod_{k=1}^{N_{I}}q(x_{k-1},x_{k}|\lambda_{k-1})p_{\Lambda|X,\Theta}(\lambda_{k}|x_{k},\theta)\right)p_{X|\Lambda}(x_{0}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta)$

We denote the last term in the previous display by $S_{N_{I}}$. It remains to
show that for any $n\in\mathbb{N}$,

$S_{n}=1-\frac{p_{\Lambda|\Theta}(\lambda|\theta)}{\rho(\lambda)}.$

For $n=1$,

$\displaystyle S_{1}$
$\displaystyle=\sum_{x_{0}\in\mathcal{X}}\sum_{x_{1}\in\mathcal{X}}\sum_{\lambda_{0}\in\mathcal{L}}\sum_{\lambda_{1}\in\mathcal{L}}\left(1-\frac{p_{\Lambda|X,\Theta}(\lambda|x_{1},\theta)}{\rho(\lambda)}\right)q(x_{0},x_{1}|\lambda_{0})$
$\displaystyle\qquad\times
p_{\Lambda|X,\Theta}(\lambda_{1}|x_{1},\theta)p_{X|\Lambda}(x_{0}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta)$
$\displaystyle=\sum_{x_{1}\in\mathcal{X}}\sum_{\lambda_{0}\in\mathcal{L}}\left(1-\frac{p_{\Lambda|X,\Theta}(\lambda|x_{1},\theta)}{\rho(\lambda)}\right)p_{X|\Lambda}(x_{1}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta)$
$\displaystyle=\sum_{x_{1}\in\mathcal{X}}\left(1-\frac{p_{\Lambda|X,\Theta}(\lambda|x_{1},\theta)}{\rho(\lambda)}\right)p_{X|\Theta}(x_{1}|\theta)$
$\displaystyle=1-\frac{p_{\Lambda|\Theta}(\lambda|\theta)}{\rho(\lambda)},$

where the second equation holds since we execute the summation over $x_{0}$
first with (2.4) and then we take summation over $\lambda_{1}$.

Next, if we assume that

$S_{n}=1-\frac{p_{\Lambda|\Theta}(\lambda|\theta)}{\rho(\lambda)}$

is true, then

$\displaystyle S_{n+1}$
$\displaystyle=\sum_{x_{0}\in\mathcal{X}}\cdots\sum_{x_{n+1}\in\mathcal{X}}\sum_{\lambda_{0}\in\mathcal{L}}\cdots\sum_{\lambda_{n+1}\in\mathcal{L}}\left(\frac{\sum_{k=1}^{n+1}[\rho(\lambda)-p_{\Lambda|X,\Theta}(\lambda|x_{k},\theta)]}{(n+1)\rho(\lambda)}\right)$
$\displaystyle\qquad\times\left(\prod_{k=1}^{n+1}q(x_{k-1},x_{k}|\lambda_{k-1})p_{\Lambda|X,\Theta}(\lambda_{k}|x_{k},\theta)\right)p_{X|\Lambda}(x_{0}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta)$
$\displaystyle=\bar{S}_{n+1}+\hat{S}_{n+1},$

where

$\displaystyle\bar{S}_{n+1}$
$\displaystyle\doteq\sum_{x_{0}\in\mathcal{X}}\cdots\sum_{x_{n+1}\in\mathcal{X}}\sum_{\lambda_{0}\in\mathcal{L}}\cdots\sum_{\lambda_{n+1}\in\mathcal{L}}\left(\frac{\rho(\lambda)-p_{\Lambda|X,\Theta}(\lambda|x_{1},\theta)}{(n+1)\rho(\lambda)}\right)$
$\displaystyle\qquad\times\left(\prod_{k=1}^{n+1}q(x_{k-1},x_{k}|\lambda_{k-1})p_{\Lambda|X,\Theta}(\lambda_{k}|x_{k},\theta)\right)p_{X|\Lambda}(x_{0}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta)$

and

$\displaystyle\hat{S}_{n+1}$
$\displaystyle\doteq\sum_{x_{0}\in\mathcal{X}}\cdots\sum_{x_{n+1}\in\mathcal{X}}\sum_{\lambda_{0}\in\mathcal{L}}\cdots\sum_{\lambda_{n+1}\in\mathcal{L}}\left(\frac{\sum_{k=2}^{n+1}[\rho(\lambda)-p_{\Lambda|X,\Theta}(\lambda|x_{k},\theta)]}{(n+1)\rho(\lambda)}\right)$
$\displaystyle\qquad\times\left(\prod_{k=1}^{n+1}q(x_{k-1},x_{k}|\lambda_{k-1})p_{\Lambda|X,\Theta}(\lambda_{k}|x_{k},\theta)\right)p_{X|\Lambda}(x_{0}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta).$

For $\bar{S}_{n+1}$, we first observe that

$\displaystyle\sum_{x_{2}\in\mathcal{X}}\cdots\sum_{x_{n+1}\in\mathcal{X}}\sum_{\lambda_{2}\in\mathcal{L}}\cdots\sum_{\lambda_{n+1}\in\mathcal{L}}\prod_{k=2}^{n+1}q(x_{k-1},x_{k}|\lambda_{k-1})p_{\Lambda|X,\Theta}(\lambda_{k}|x_{k},\theta)=1.$

It is not hard to check the equality holds by taking summation in the
following order: $\lambda_{n+1}\to x_{n+1}\to\lambda_{n}\to
x_{n}\to\cdots\to\lambda_{2}\to x_{2}$. This implies that

$\displaystyle\bar{S}_{n+1}$
$\displaystyle=\sum_{x_{0}\in\mathcal{X}}\sum_{x_{1}\in\mathcal{X}}\sum_{\lambda_{0}\in\mathcal{L}}\sum_{\lambda_{1}\in\mathcal{L}}\left(\frac{\rho(\lambda)-p_{\Lambda|X,\Theta}(\lambda|x_{n+1},\theta)}{(n+1)\rho(\lambda)}\right)q(x_{0},x_{1}|\lambda_{0})$
$\displaystyle\qquad\times
p_{\Lambda|X,\Theta}(\lambda_{1}|x_{1},\theta)p_{X|\Lambda}(x_{0}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta)$
$\displaystyle=\frac{1}{n+1}S_{1}.$

As for $\hat{S}_{n+1}$, since

$\displaystyle\sum_{x_{0}\in\mathcal{X}}\sum_{\lambda_{0}\in\mathcal{L}}q(x_{0},x_{1}|\lambda_{0})p_{\Lambda|X,\Theta}(\lambda_{1}|x_{1},\theta)p_{X|\Lambda}(x_{0}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta)$
$\displaystyle=\sum_{\lambda_{0}\in\mathcal{L}}p_{\Lambda|X,\Theta}(\lambda_{1}|x_{1},\theta)p_{X|\Lambda}(x_{1}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta)$
$\displaystyle=p_{\Lambda|X,\Theta}(\lambda_{1}|x_{1},\theta)p_{X|\Theta}(x_{1}|\theta)$
$\displaystyle=p_{X|\Lambda}(x_{1}|\lambda_{1})p_{\Lambda|\Theta}(\lambda_{1}|\theta),$

we have

$\displaystyle\hat{S}_{n+1}$
$\displaystyle=\sum_{x_{1}\in\mathcal{X}}\cdots\sum_{x_{n+1}\in\mathcal{X}}\sum_{\lambda_{1}\in\mathcal{L}}\cdots\sum_{\lambda_{n+1}\in\mathcal{L}}\left(\frac{\sum_{k=2}^{n+1}[\rho(\lambda)-p_{\Lambda|X,\Theta}(\lambda|x_{k},\theta)]}{(n+1)\rho(\lambda)}\right)$
$\displaystyle\qquad\times\left(\prod_{k=2}^{n+1}q(x_{k-1},x_{k}|\lambda_{k-1})p_{\Lambda|X,\Theta}(\lambda_{k}|x_{k},\theta)\right)p_{X|\Lambda}(x_{1}|\lambda_{1})p_{\Lambda|\Theta}(\lambda_{1}|\theta)$
$\displaystyle=\frac{n}{n+1}S_{n}.$

Hence,

$S_{n+1}=\bar{S}_{n+1}+\hat{S}_{n+1}=\frac{1}{n+1}S_{1}+\frac{n}{n+1}S_{n}=1-\frac{p_{\Lambda|\Theta}(\lambda|\theta)}{\rho(\lambda)}.$

We complete the proof by mathematical induction.

The limit ODE can now be identified by

$\displaystyle\dot{\theta}=\bar{g}(\theta)+z,$ (3.5)

where $z\in-C(\theta)$ is the reflection due to the projection onto the
hyperrectangle $H$. For $\theta$ in the interior of $H$, $C(\theta)=\\{0\\}$;
for $\theta$ on the boundary of $H$, $C(\theta)$ is the infinite convex cone
generated by the outer normals at $\theta$ of the faces on which $\theta$
lies, see [18, Section 4.3] for additional details.

To identify a set containing the limit points of the stochastic approximation
algorithm, it is sufficient to construct a Lyapunov function, $V$, for the
limit ODE and find the level set of $0$. If $\left\langle\nabla
V(\theta),\bar{g}(\theta)\right\rangle\leq 0$, with equality if and only if
$V(\theta)=0$, then the Lyapunov function decreases under the dynamics of the
limit ODE until it reaches $V(\theta(t))=0$, after which it remains constant
at $0$. Indeed,

$V(\theta(t_{1}))-V(\theta(t_{0}))=\int_{t_{0}}^{t_{1}}\langle\nabla
V(\theta(t)),\dot{\theta}(t)\rangle dt=\int_{t_{0}}^{t_{1}}\langle\nabla
V(\theta(t)),\bar{g}(\theta(t))\rangle dt<0,$

where the reflection term has been ignored for convenience. Consequently, the
limit set of the ODE is contained in the level set where $V(\theta)=0$.

To construct a Lyapunov function for the limit ODE of the AWH algorithm, let
$V:\mathbb{R}^{|\mathcal{L}|}\to[0,\infty)$ be given by

$\displaystyle
V(\theta)=\sum_{\sigma\in\mathcal{L}}\rho(\sigma)[\bar{g}(\theta)(\sigma)]^{2}=\sum_{\sigma\in\mathcal{L}}\rho(\sigma)\left(1-\frac{p_{\Lambda|\Theta}(\sigma|\theta)}{\rho(\sigma)}\right)^{2}.$
(3.6)

###### Lemma 3.2

Let $V$ be given by (3.6). Then $V$ is continuously differentiable and
$\left\langle\nabla V(\theta),\bar{g}(\theta)\right\rangle\leq 0$ with
equality if and only if $V(\theta)=0$.

Proof. It follows immediately from the definition that $V$ is continuously
differentiable. It remains to prove $\left\langle\nabla
V(\theta),\bar{g}(\theta)\right\rangle\leq 0$ with equality if and only if
$V(\theta)=0$. To determine the gradient, $\nabla V$, note first that since
$p_{\Lambda|\Theta}$ is given by (2.2),

$\displaystyle\frac{\partial
p_{\Lambda|\Theta}(\sigma|\theta)}{\partial\theta(\lambda)}=\left\\{\begin{array}[]{ll}-p_{\Lambda|\Theta}(\sigma|\theta)p_{\Lambda|\Theta}(\lambda|\theta),&\sigma\neq\lambda,\\\
p_{\Lambda|\Theta}(\lambda|\theta)(1-p_{\Lambda|\Theta}(\lambda|\theta)),&\sigma=\lambda.\end{array}\right.$

Consequently, the gradient of $V$ is given by

$\displaystyle\frac{\partial V(\theta)}{\partial\theta(\lambda)}$
$\displaystyle=\frac{\partial}{\partial\theta(\lambda)}\left(\sum_{\sigma\in\mathcal{L}}\rho(\sigma)\left(1-\frac{p_{\Lambda|\Theta}(\sigma|\theta)}{\rho(\sigma)}\right)^{2}\right)$
$\displaystyle=\sum_{\sigma\in\mathcal{L}}2\left(1-\frac{p_{\Lambda|\Theta}(\sigma|\theta)}{\rho(\sigma)}\right)\left(-\frac{\partial
p_{\Lambda|\Theta}(\sigma|\theta)}{\partial\theta(\lambda)}\right)$
$\displaystyle=2\sum_{\sigma\in\mathcal{L}}\left(1-\frac{p_{\Lambda|\Theta}(\sigma|\theta)}{\rho(\sigma)}\right)p_{\Lambda|\Theta}(\sigma|\theta)p_{\Lambda|\Theta}(\lambda|\theta)$
$\displaystyle\quad-2\left(1-\frac{p_{\Lambda|\Theta}(\lambda|\theta)}{\rho(\lambda)}\right)p_{\Lambda|\Theta}(\lambda|\theta)$
$\displaystyle=2\sum_{\sigma\in\mathcal{L}}\bar{g}(\theta)(\sigma)p_{\Lambda|\Theta}(\sigma|\theta)p_{\Lambda|\Theta}(\lambda|\theta)-2\bar{g}(\theta)(\lambda)p_{\Lambda|\Theta}(\lambda|\theta).$

Taking the inner product with $\bar{g}(\theta)$ yields,

$\displaystyle\frac{1}{2}\left\langle\nabla
V(\theta),\bar{g}(\theta)\right\rangle$
$\displaystyle=\sum_{\lambda\in\mathcal{L}}\left[\sum_{\sigma\in\mathcal{L}}\bar{g}(\theta)(\sigma)p_{\Lambda|\Theta}(\sigma|\theta)p_{\Lambda|\Theta}(\lambda|\theta)\right]\bar{g}(\theta)(\lambda)$
$\displaystyle\qquad-\sum_{\lambda\in\mathcal{L}}\bar{g}(\theta)(\lambda)p_{\Lambda|\Theta}(\lambda|\theta)\bar{g}(\theta)(\lambda)$
$\displaystyle=\left(\sum_{\lambda\in\mathcal{L}}\bar{g}(\theta)(\lambda)p_{\Lambda|\Theta}(\lambda|\theta)\right)^{2}-\sum_{\lambda\in\mathcal{L}}[\bar{g}(\theta)(\lambda)]^{2}p_{\Lambda|\Theta}(\lambda|\theta)$
$\displaystyle=-\text{Var}_{p_{\Lambda|\Theta}(\cdot|\theta)}(\bar{g}(\theta))$
$\displaystyle\leq 0.$

In the last equality it is observed that the desired expression can be
represented as the variance of $\bar{g}(\theta)$ under the distribution
$p_{\Lambda|\Theta}(\cdot|\theta)$ on $\mathcal{L}$. Consequently,
$\left\langle\nabla V(\theta),\bar{g}(\theta)\right\rangle=0$ if and only if
$\bar{g}(\theta)$ is a constant, i.e., there exists a constant
$c\in\mathbb{R}$ such that

$\bar{g}(\theta)(\lambda)=1-\frac{p_{\Lambda|\Theta}(\lambda|\theta)}{\rho(\lambda)}=c,\quad\text{
for all }\lambda\in\mathcal{L}.$

Multiplying both sides by $\rho(\lambda)$ and summing over $\lambda$, leads to

$\displaystyle c$
$\displaystyle=\sum_{\lambda\in\mathcal{L}}c\rho(\lambda)=\sum_{\lambda\in\mathcal{L}}\left(\rho(\lambda)-p_{\Lambda|\Theta}(\lambda|\theta)\right)=0.$

This implies that $c$ must be equal to $0$, and thus $\left\langle\nabla
V(\theta),\bar{g}(\theta)\right\rangle=0$ if and only if

$V(\theta)=\sum_{\lambda\in\mathcal{L}}\rho(\lambda)[\bar{g}(\theta)(\lambda)]^{2}=c^{2}=0.$

This completes the proof.

###  3.2 Convergence of estimates of free energy differences

In this subsection, the main theorem on the almost sure convergence of the
estimates of the design parameters $\\{\theta_{n}\\}$ of the AWH algorithm is
provided in Theorem 3.3. The subsequent Corollary 3.4 provides the result on
the almost sure convergence of estimates of free energy differences.

###### Theorem 3.3

Let $\\{\theta_{n}\\}$ be given by the AWH algorithm (3.3), with target
distribution $\rho$, $\rho(\lambda)>0,\lambda\in\mathcal{L}$, let $H$ a hyper-
rectangle in $\mathbb{R}^{|\mathcal{L}|}$ and $V$ be given by (3.6). Assume
that $H$ is large enough so that $\\{\theta\in H:V(\theta)=0\\}$ is non-empty.
Then $\\{\theta_{n}\\}$ converges almost surely to $\\{\theta\in
H:V(\theta)=0\\}$.

Proof. As noted previous, the recursion of the AWH can be written as a
stochastic approximation algorithm with state-dependent noise as in (3.3). It
follows from Theorem A.1 in the Appendix, once the conditions of that theorem
are verified, that $\\{\theta_{n}\\}$ converges almost surely to some limit
set of the limit ODE, given by (3.5). By Lemma 3.2 it follows that the limit
set is contained in $\\{\theta\in H:V(\theta)=0\\}$, see e.g. [18, Theorem
4.2.3].

To complete the proof, the following conditions of Theorem A.1 must be
verified: (5.1.1), (A1.1), (A1.2), (A1.4), (6.2), (A6.1), (A6.2), and
(A4.3.1).

Condition (5.1.1) is satisfied since $\varepsilon_{n}=1/(n+1)$, for $n\geq 0$,
and consequently, $\sum_{n=0}\varepsilon_{n}=\infty$, $\varepsilon_{n}\geq 0$,
and $\lim_{n\to\infty}\varepsilon_{n}=0$. To verify (A1.1), the upper bound

$\displaystyle\left|h(\theta,\xi)(\lambda)\right|=\left|1-\frac{\sum_{k=1}^{N_{I}}p_{\Lambda|X,\Theta}(\lambda|\xi_{k},\theta)}{N_{I}\rho(\lambda)}\right|\leq
1+\frac{1}{\rho(\lambda)},$ (3.7)

implies that $|h(\theta,\xi)(\lambda)|\leq C$ with $C\doteq
1+\max_{\lambda\in\mathcal{L}}(1/\rho(\lambda))<\infty$ for all $\theta,\xi$
and $\lambda$. Consequently,
$\sup_{n}E|Y_{n}|=\sup_{n}E|h(\theta_{n},\xi_{n+1})|\leq C<\infty.$ The
condition (A1.2) states that $g_{n}(\theta,\xi)=g(\theta,\xi)$ is continuous
in $\theta$ for each $\xi$. It follows immediately from the definitions, since
both $p_{\Lambda|X,\Theta}(\lambda|\xi(k),\theta)$ and $p(\xi,\eta|\theta)$
are continuous in $\theta$ for each $\xi,\eta$, $\lambda$ and
$k\in\\{1,2,\dots,N_{I}\\}$. To verify condition (A1.4) it is sufficient to
prove

$\displaystyle\text{(A2.1)}\quad\sup_{n}E|Y_{n}|^{2}<\infty\quad\text{and}\qquad\text{(A2.4)}\quad\sum_{n}\varepsilon_{n}^{2}<\infty.$

Indeed, according to the last paragraph on p. 137 in [18], (A2.1) and (A2.4)
implies the stronger statement (3.2) in [18], which says for some $T>0$,

$\lim_{n\rightarrow\infty}\left(\sup_{j\geq n}\max_{0\leq t\leq
T}\left|\sum_{i=m(jT)}^{m(jT+t)-1}\varepsilon_{i}\delta
M_{i}\right|\right)=0,\quad\text{ a.s}.$

Thus, for each $\mu>0$ and some $T>0$,

$\lim_{n\rightarrow\infty}P\left\\{\sup_{j\geq n}\max_{0\leq t\leq
T}\left|\sum_{i=m(jT)}^{m(jT+t)-1}\varepsilon_{i}\delta
M_{i}\right|\geq\mu\right\\}=0,$

which is condition (A1.4). The conditions (A2.1) and (A2.4) are satisfied by
the definition of $\varepsilon_{n}$ and since the function $h$ is bounded by
(3.7). Condition (6.2) states that
$P(\xi_{n+1}\in\cdot|\xi_{i},\theta_{i},i\leq n)=p(\xi_{n},\cdot|\theta_{n})$,
where $p(\xi,\cdot|\theta)$ denotes the one-step transition probability with
starting point $\xi$ and parameterized by $\theta$, which has already been
observed in connection to (3.2).

To verify conditions (A6.1) and (A6.2), recall that sufficient conditions for
them (provided in [18, Section 5.3]) are
(A3.1) For each $\mu>0$,

$\sum_{n}e^{-\mu/\varepsilon_{n}}<\infty.$

(A3.2) For some $T<\infty$, there is a $c_{1}(T)<\infty$ such that for all
$n$,

$\sup_{n\leq i\leq m(t_{n}+T)}\frac{\varepsilon_{i}}{\varepsilon_{n}}\leq
c_{1}(T).$

(A3.3) There is a real $K<\infty$ such that for small real $\gamma$, all $n$,
and each component $\alpha_{n,j}$ of $\alpha_{n}$, $\delta N_{n,j}$ of $\delta
N_{n}$,

$E_{n}e^{\gamma\alpha_{n,j}}\leq e^{\gamma^{2}K/2}\quad\mbox{ and }\quad
E_{n}e^{\gamma(\delta N_{n,j})}\leq e^{\gamma^{2}K/2}.$

One can easily verify that conditions (A3.1) and (A3.2) hold when
$\varepsilon_{n}=1/(n+1)$. It remains to check if (A3.3) holds for both
$\alpha_{n}$ and $\delta N_{n}$. To prove this, let $\bar{g}$ be given by
(3.4) and note that, by Lemma 3.1

$\bar{g}(\theta)=E[g(\theta,\xi_{1})|\theta_{0}=\theta],$

where the expectation is taken under the invariant distribution of $\xi_{1}$
given $\theta_{0}=\theta$. Since $p_{\Lambda|\Theta}(\lambda|\theta)$ is
continuous in $\theta$ for each $\lambda\in\mathcal{L}$ it follows that
$\bar{g}$ is continuous. Moreover, since the state space
$\mathcal{X}^{N_{I}}\times\mathcal{L}^{N_{I}}$ of $\\{\xi_{n}\\}$ is finite,
and $\\{\xi_{n}\\}$ is irreducible and aperiodic, $\\{\xi_{n}\\}$ is a
geometrically ergodic Markov chain. Indeed, the irreducibility and
aperiodicity of $\\{\xi_{n}\\}$ follows from the definition of its transition
probability $p(\xi,\eta|\theta)$ in (3.2), the assumption of irreducibility
and aperiodicity on $\\{q(x,y|\lambda)\\}_{x,y}$ and the definition of
$p_{\Lambda|X,\Theta}(\lambda|x,\theta)$. Consequently, the summands in the
definition of $v_{n}(\theta,\xi)$ converges to zero at a geometric rate. This
implies that $\alpha_{n}$ and $\delta N_{n}$ are bounded. As a result, they
both satisfy (A3.3) and the asymptotic rates of change of the processes
$A^{0}(t)$ and $N^{0}(t)$ are zero with probability one. Lastly, (A4.3.1)
holds due to our assumption on $H$. This completes the proof.

###### Corollary 3.4

Let $\\{\theta_{n}\\}$ be given by the AWH algorithm (3.3), with target
distribution $\rho$, $\rho(\lambda)>0,\lambda\in\mathcal{L}$ and $H$ a hyper-
rectangle in $\mathbb{R}^{|\mathcal{L}|}$ such that $\\{\theta\in
H:V(\theta)=0\\}$ is non-empty, where $V$ is given by (3.6). Then, for each
$\lambda,\sigma\in\mathcal{L}$, the sequence
$\\{\theta_{n}(\lambda)-\theta_{n}(\sigma)\\}$, converges almost surely to
$F(\lambda)-F(\sigma)+\log\rho(\lambda)-\log\rho(\sigma)$.

Proof. By Theorem 3.3 the sequence $\\{\theta_{n}\\}$ converges almost surely
to the set $\\{\theta\in H:V(\theta)=0\\}$. For any
$\theta\in\mathbb{R}^{|\mathcal{L}|}$ with $V(\theta)=0$ it follows that
$p_{\Lambda|\Theta}(\lambda|\theta)=\rho(\lambda)$, for each
$\lambda\in\mathcal{L}$, and consequently,

$\displaystyle\theta(\lambda)-F(\lambda)=\log\rho(\lambda)+\log\sum_{\tau\in\mathcal{L}}e^{\theta(\tau)-F(\tau)}.$

Hence, for $\lambda,\sigma\in\mathcal{L}$,

$\displaystyle\theta(\lambda)-\theta(\sigma)=F(\lambda)-F(\sigma)+\log\rho(\lambda)-\log\rho(\sigma).$

This completes the proof.

From Corollary 3.4 it follows that estimates of free energy differences can be
obtained by running the AWH algorithm for $N$ iterations, to obtain
$\theta_{N}$, and then estimate the free energy difference
$F(\lambda)-F(\sigma)$ by

$\displaystyle\theta_{N}({\lambda)}-\theta_{N}(\sigma)-\log\rho(\lambda)+\log\rho(\sigma).$

##  4 Almost Sure Convergence of Ergodic Averages

In addition to estimating free energy differences, the AWH algorithm may be
used to estimate expectations at a fixed parameter value. Namely, for a fixed
$\lambda\in\mathcal{L}$ and a function
$\psi:\mathcal{X}\rightarrow\mathbb{R}$, consider estimating

$\displaystyle\sum_{x\in\mathcal{X}}\psi(x)p_{X|\Lambda}(x|\lambda),$ (4.1)

which assumed to be finite.

At the $n$th iteration of the AWH algorithm, given the free energy estimate
$\theta_{n}$, the Markov chain
$\\{(X_{nN_{I}+k},\Lambda_{nN_{I}+k}),k=1,\dots,N_{I}\\}$ has invariant
distribution $p_{X,\Lambda|\Theta}(x,\lambda|\theta_{n})$ given by (2.1). The
estimate of (4.1) that we proposed in Subsection 2.2 is given by

$\displaystyle\frac{\sum_{n=0}^{N-1}\sum_{k=1}^{N_{I}}\psi(X_{nN_{I}+k})p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})}{\sum_{n=0}^{N-1}\sum_{k=1}^{N_{I}}p_{X|\Lambda,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})}.$
(4.2)

To study the almost sure convergence of the estimator in the last display, it
is sufficient to study the almost sure convergence of

$\displaystyle\hat{\phi}_{N}=\frac{1}{N\,N_{I}}\sum_{n=0}^{N-1}\sum_{k=1}^{N_{I}}\phi(X_{nN_{I}+k})p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n}),$

which is an estimate of

$\displaystyle\sum_{x\in\mathcal{X}}\phi(x)p_{\Lambda|X,\Theta}(\lambda|x,\theta_{*})p_{X|\Theta}(x|\theta_{*})=\sum_{x\in\mathcal{X}}\phi(x)p_{X|\Lambda}(x|\lambda)p_{\Lambda|\Theta}(\lambda|\theta_{*}),$

where the optimal design parameter $\theta_{*}$ satisfying

$\displaystyle
p_{\Lambda|\Theta}(\lambda|\theta_{*})=\rho(\lambda),\quad\lambda\in\mathcal{L}.$

With $\phi=\psi$ and $\phi\equiv 1$ it follows that the estimator (4.2) is
given by $\hat{\psi}_{N}/\hat{1}_{N}$, which estimates the desired quantity

$\displaystyle\frac{\sum_{x\in\mathcal{X}}\psi(x)p_{X|\Lambda}(x|\lambda)p_{\Lambda|\Theta}(\lambda|\theta_{*})}{\sum_{x\in\mathcal{X}}p_{X|\Lambda}(x|\lambda)p_{\Lambda|\Theta}(\lambda|\theta_{*})}=\sum_{x\in\mathcal{X}}\psi(x)p_{X|\Lambda}(x|\lambda).$

In order to show the almost sure convergence of the sequence
$\\{\hat{\phi}_{n}\\}_{n\geq 1}$, observe that $\\{\hat{\phi}_{n}\\}_{n\geq
1}$ can be written as an iteration, with $\hat{\phi}_{0}=0$ and for $n\geq 0$,

$\displaystyle\hat{\phi}_{n+1}=\hat{\phi}_{n}+\frac{1}{n+1}\left(\Phi(\xi_{n},\theta_{n})-\hat{\phi}_{n}\right),$

where

$\Phi(\xi,\theta)=\frac{1}{N_{I}}\sum_{k=1}^{N_{I}}\phi\left(\xi(k)\right)p_{\Lambda|X,\Theta}\left(\lambda|\xi(k),\theta\right),\;\theta\in\mathbb{R}^{|\mathcal{L}|}$

and, as in the previous section,

$\xi_{n}=(X_{(n-1)N_{I}+1},\dots,X_{nN_{I}},\Lambda_{(n-1)N_{I}+1},\dots,\Lambda_{nN_{I}}).$

Consider the extended design parameter
$\vartheta=(\theta,\zeta)\in\mathbb{R}^{|\mathcal{L}|}\times[-C,C]$ with
$C\doteq\max_{x\in\mathcal{X}}|\phi(x)|<\infty$, and without loss of
generality, we assume $0\not\in\mathcal{L}$ and define
$\mathcal{L}_{0}=\mathcal{L}\cup\\{0\\}$. For $u\in\mathcal{L}$ we interpret
$\vartheta(u)=\theta(u)$, and for $u=0$, $\vartheta(u)=\zeta$. The optimal
value of extended design parameter is $\vartheta_{*}=(\theta_{*},\zeta_{*})$,
where
$\zeta_{*}=\sum_{x\in\mathcal{X}}\phi(x)p_{X|\Lambda}(x|\lambda)p_{\Lambda|\Theta}(\lambda|\theta_{*})\in[-C,C]$.
The estimator of the optimal extended design parameter is

$\vartheta_{n}=(\theta_{n},\zeta_{n})=((\theta_{n}(\sigma):\sigma\in\mathcal{L}),\hat{\phi}_{n}),$

$\varepsilon_{n}=1/(n+1)$ and the extended function
$h:(\mathbb{R}^{|\mathcal{L}|}\times[-C,C])\times(\mathcal{X}^{N_{I}}\times\mathcal{L}^{N_{I}})\to\mathbb{R}^{|\mathcal{L}_{0}|}$
given by

$\displaystyle
h(\vartheta,\xi)(u)=\left\\{\begin{array}[]{ll}1-\frac{\sum_{k=1}^{N_{I}}p_{\Lambda|X,\Theta}(u|\xi({k}),\theta)}{N_{I}\rho(u)},&\text{
if }u\in\mathcal{L},\\\ \Phi(\xi,\theta)-\zeta,&\text{ if
}u=0.\end{array}\right.$

The recursion of the AWH algorithm for the extended design parameter can be
stated as

$\displaystyle\vartheta_{n+1}$
$\displaystyle=\vartheta_{n}+\varepsilon_{n}\left(h(\vartheta_{n},\xi_{n+1})+\beta_{n}+Z_{n}\right)$
$\displaystyle=\vartheta_{n}+\varepsilon_{n}Y_{n}+\varepsilon_{n}Z_{n},$

where $Y_{n}=h(\vartheta_{n},\xi_{n+1})+\beta_{n}$. Let
$\mathcal{F}_{n}=\sigma(\theta_{0},Y_{0},\dots,Y_{n-1},\xi_{1},\dots,\xi_{n})$,
$p(\xi,\cdot|\theta)$ be as in (3.2), and

$g(\vartheta,\xi)\doteq\sum_{\eta}h(\vartheta,\eta)p(\xi,\eta|\theta).$

Since $h$ is bounded, it follows that $E|Y_{n}|<\infty$ and

$E_{n}Y_{n}=g(\vartheta_{n},\xi_{n})+\beta_{n}.$

By introducing

$\displaystyle\delta M_{n}$
$\displaystyle=Y_{n}-E_{n}Y_{n}=h(\vartheta_{n},\xi_{n+1})-g(\vartheta_{n},\xi_{n})$
$\displaystyle=h(\vartheta_{n},\xi_{n+1})-\sum_{\eta}h(\vartheta_{n},\eta)p(\xi_{n},\eta|\theta)$
$\displaystyle=\sum_{\eta}[h(\vartheta_{n},\xi_{n+1})-h(\vartheta_{n},\eta)]p(\xi_{n},\eta|\theta)$

the recursion may be restated as

$\vartheta_{n+1}=\vartheta_{n}+\varepsilon_{n}[E_{n}Y_{n}+\delta M_{n}+Z_{n}]$

or equivalently,

$\displaystyle\vartheta_{n+1}=\vartheta_{n}+\varepsilon_{n}[g(\vartheta_{n},\xi_{n})+\delta
M_{n}+\beta_{n}+Z_{n}],$ (4.3)

which, again, is a standard form of a stochastic approximation algorithm with
state-dependent noise, see [18, Section 6.6].

###  4.1 Identification of the limit ODE

To derive the limit ODE we proceed as in the previous section with the
extended design parameter $\vartheta=(\theta,\zeta)$. Let
$\bar{g}:\mathbb{R}^{|\mathcal{L}|}\times[-C,C]\to\mathbb{R}^{|\mathcal{L}_{0}|}$
be given by

$\displaystyle\bar{g}(\vartheta)(u)$
$\displaystyle=\left\\{\begin{array}[]{ll}1-\frac{p_{\Lambda|\Theta}(u|\theta)}{\rho(u)},&\text{if
}u\in\mathcal{L},\\\
\sum_{\mathcal{X}}\phi(x)p_{X,\Lambda|\Theta}(x,\lambda|\theta)-\zeta,&\text{if
}u=0.\end{array}\right.$ (4.6)

Similar to Section 3 the expression for $\bar{g}$ is the ergodic average of
$h(\vartheta,\xi_{1})$, for fixed $\vartheta_{0}=(\theta,\zeta)$. That is, the
expectation over
$\xi_{1}=(X_{1},\dots,X_{N_{I}},\Lambda_{1},\dots,\Lambda_{N_{I}})$ where
$X_{0}$ is sampled from $p_{X|\Theta}(\cdot|\theta)$.

###### Lemma 4.1

Let $\bar{g}$ be given by (4.6). For
$\vartheta\in\mathbb{R}^{|\mathcal{L}|}\times[-C,C]$ and
$u\in\mathcal{L}_{0}$,
$\bar{g}(\vartheta)(u)=E[h(\vartheta,\xi_{1})(u)|\vartheta_{0}=(\theta,\zeta)]$.

Proof. For $u\in\mathcal{L}$ the statement reduces to that of Lemma 3.1, so it
is sufficient to consider $u=0$. The expectation on the right hand side can be
written

$\displaystyle E[h(\vartheta,\xi_{1})(u)|\vartheta_{0}=(\theta,\zeta)]$
$\displaystyle=E[h(\vartheta,(X_{1},\dots,X_{N_{I}},\Lambda_{1},\dots,\Lambda_{N_{I}}))(u)|\vartheta_{0}=(\theta,\zeta)]$
$\displaystyle=\sum_{x_{0}\in\mathcal{X}}\cdots\sum_{x_{N_{I}}\in\mathcal{X}}\sum_{\lambda_{0}\in\mathcal{L}}\cdots\sum_{\lambda_{N_{I}}\in\mathcal{L}}\left(\frac{1}{N_{I}}\sum_{k=1}^{N_{I}}\phi(x_{k})p_{\Lambda|X,\Theta}(\lambda|x_{k},\theta)-\zeta\right)$
$\displaystyle\qquad\times\left(\prod_{k=1}^{N_{I}}q(x_{k-1},x_{k}|\lambda_{k-1})p_{\Lambda|X,\Theta}(\lambda_{k}|x_{k},\theta)\right)p_{X|\Lambda}(x_{0}|\lambda_{0})p_{\Lambda|\Theta}(\lambda_{0}|\theta)$

By applying the mathematical induction argument as in the proof of Lemma 3.1
with minor modifications, we can find that the last term in the previous
display equals

$\sum_{\mathcal{X}}\phi(x)p_{X,\Lambda|\Theta}(x,\lambda|\theta)-\zeta.$

This completes the proof.

The limit ODE is given by

$\displaystyle\dot{\vartheta}=\bar{g}(\vartheta)+z,$ (4.7)

where $z\in-C(\vartheta)$ is the reflection. To construct a Lyapunov function
for the limit ODE, recall that $\theta_{*}$ denotes the optimal design
parameter and let, for $\delta>0$,
$V_{\delta}:\mathbb{R}^{|\mathcal{L}|}\times[-C,C]\to[0,\infty)$ be given by

$\displaystyle
V_{\delta}(\vartheta)\\!=\\!\sum_{\sigma\in\mathcal{L}}\rho(\sigma)\left(1\\!-\\!\frac{p_{\Lambda|\Theta}(\sigma|\theta)}{\rho(\sigma)}\right)^{2}\\!\\!+\delta\left(\sum_{x\in\mathcal{X}}\phi(x)p_{X,\Lambda|\Theta}(x,\lambda|\theta_{*})\\!-\\!\zeta\right)^{2}\\!\\!.$
(4.8)

###### Lemma 4.2

Let $V_{\delta}$ be defined as in (4.8) with $\delta=m^{2}/(2C^{2})$, where
$C\doteq\max_{x\in\mathcal{X}}|\phi(x)|$ and

$m\doteq\inf_{\theta\in
H,\sigma\in\mathcal{L}}p_{\Lambda|\Theta}(\sigma|\theta)=\inf_{\theta\in
H,\sigma\in\mathcal{L}}\left(\frac{e^{\theta(\sigma)-F(\sigma)}}{\sum_{\tau\in\mathcal{L}}e^{\theta(\tau)-F(\tau)}}\right)>0.$

Then $V_{\delta}$ is continuously differentiable and, for $\delta>0$
sufficiently small, $\left\langle\nabla
V_{\delta}(\vartheta),\bar{g}(\vartheta)\right\rangle\leq 0$ with equality if
and only if $V_{\delta}(\vartheta)=0$.

Proof. It follows immediately from the definition that $V_{\delta}$ is
continuously differentiable. It remains to prove $\left\langle\nabla
V_{\delta}(\vartheta),\bar{g}(\vartheta)\right\rangle\leq 0$ with equality if
and only if $V(\vartheta)=0$.

From the proof of Lemma 3.2 it follows that for $u\in\mathcal{L}$

$\displaystyle\frac{\partial V_{\delta}(\vartheta)}{\partial\vartheta(u)}$
$\displaystyle=\sum_{\sigma\in\mathcal{L}}2\left(1-\frac{p_{\Lambda|\Theta}(\sigma|\theta)}{\rho(\sigma)}\right)\left(-\frac{\partial
p_{\Lambda|\Theta}(\sigma|\theta)}{\partial\theta(u)}\right)$
$\displaystyle=2\sum_{\sigma\in\mathcal{L}}\bar{g}(\vartheta)(\sigma)p_{\Lambda|\Theta}(\sigma|\theta)p_{\Lambda|\Theta}(u|\theta)-2\bar{g}(\vartheta)(u)p_{\Lambda|\Theta}(u|\theta).$

For $u=0$ we have

$\displaystyle\frac{\partial V_{\delta}(\vartheta)}{\partial\vartheta(u)}$
$\displaystyle=-2\delta\bar{g}((\theta_{*},\zeta))(u).$

Taking the inner product with $\bar{g}(\vartheta)$ yields,

$\displaystyle\frac{1}{2}\left\langle\nabla
V_{\delta}(\vartheta),\bar{g}(\vartheta)\right\rangle$
$\displaystyle=\sum_{\tau\in\mathcal{L}}\left[\sum_{\sigma\in\mathcal{L}}\bar{g}(\vartheta)(\sigma)p_{\Lambda|\Theta}(\sigma|\vartheta)p_{\Lambda|\Theta}(\tau|\vartheta)\right]\bar{g}(\vartheta)(\tau)$
$\displaystyle\qquad-\sum_{\tau\in\mathcal{L}}\bar{g}(\vartheta)(\tau)p_{\Lambda|\Theta}(\tau|\vartheta)\bar{g}(\vartheta)(\tau)-\delta\bar{g}((\theta_{*},\zeta))(0)\cdot\bar{g}((\theta,\zeta))(0)$
$\displaystyle=\left(\sum_{\tau\in\mathcal{L}}\bar{g}(\vartheta)(\tau)p_{\Lambda|\Theta}(\tau|\vartheta)\right)^{2}-\sum_{\tau\in\mathcal{L}}[\bar{g}(\vartheta)(\tau)]^{2}p_{\Lambda|\Theta}(\tau|\vartheta)$
$\displaystyle\qquad-\delta\bar{g}((\theta_{*},\zeta))(0)\cdot\bar{g}((\theta,\zeta))(0)$
$\displaystyle\leq\left(\sum_{\tau\in\mathcal{L}}\bar{g}(\vartheta)(\tau)p_{\Lambda|\Theta}(\tau|\vartheta)\right)^{2}-\sum_{\tau\in\mathcal{L}}[\bar{g}(\vartheta)(\tau)]^{2}p_{\Lambda|\Theta}(\tau|\vartheta)$
$\displaystyle\qquad+\frac{\delta}{2}\left[\bar{g}((\theta_{*},\zeta))(0)-\bar{g}((\theta,\zeta))(0)\right]\left[\bar{g}((\theta_{*},\zeta))(0)-\bar{g}((\theta,\zeta))(0)\right]$
$\displaystyle\qquad-\frac{\delta}{2}[\bar{g}((\theta,\zeta))(0)]^{2}.$

The inequality holds since for any $x,y\in\mathbb{R}$,
$-2xy\leq(x-y)^{2}-y^{2}$.

Note that the first and third terms on the right of the inequality are always
non-positive. The middle term can be written as

$\displaystyle\frac{\delta}{2}\left(\sum_{x\in\mathcal{X}}\phi(x)p_{X|\Lambda}(x|\lambda)\right)^{2}\left(p_{\Lambda|\Theta}(\lambda|\theta_{*})-p_{\Lambda|\Theta}(\lambda|\theta)\right)^{2}.$
(4.9)

The sum appearing in the square is bounded by
$C\doteq\max_{x\in\mathcal{X}}|\phi(x)|<\infty$, and since
$p_{\Lambda|\Theta}(\lambda|\theta_{*})=\rho(\lambda)$ and $\rho(\lambda)\leq
1$, the expression in (4.9) can be bounded above by

$\displaystyle\delta
C^{2}\left(1-\frac{p_{\Lambda|\Theta}(\lambda|\theta)}{\rho(\lambda)}\right)^{2}.$

Consequently, to show that $\frac{1}{2}\left\langle\nabla
V_{\delta}(\vartheta),\bar{g}(\vartheta)\right\rangle\leq 0$ for all
$\vartheta$, it suffices to show that for all $\vartheta$

$\displaystyle\left(\sum_{\tau\in\mathcal{L}}\bar{g}(\vartheta)(\tau)p_{\Lambda|\Theta}(\tau|\vartheta)\right)^{2}-\sum_{\tau\in\mathcal{L}}[\bar{g}(\vartheta)(\tau)]^{2}p_{\Lambda|\Theta}(\tau|\vartheta)+\delta
C^{2}\left(1-\frac{p_{\Lambda|\Theta}(\lambda|\theta)}{\rho(\lambda)}\right)^{2}$
(4.10)

is always non-positive.

To prove the last statement, recall that

$\bar{g}(\vartheta)(\tau)=1-\frac{p_{\Lambda|\Theta}(\tau|\theta)}{\rho(\tau)},\text{for
}\tau\in\mathcal{L},$

and in order to make the discussion easier, the following simplified notation
will be used. Suppose that $\mathcal{L}=\\{\lambda_{1},\dots,\lambda_{n}\\}$
for some $n\in\mathbb{N}$, denote $p_{\Lambda|\Theta}(\lambda_{i}|\theta)$ by
$p_{i}$, $\rho(\lambda_{i})$ by $\rho_{i}$, and let $c_{i}=p_{i}/\rho_{i}$ for
all $i$. To show that (4.10) is non-positive, it is equivalent to show that

$\left(\sum_{i=1}^{n}p_{i}(1-c_{i})\right)^{2}-\sum_{i=1}^{n}p_{i}(1-c_{i})^{2}+\delta
C^{2}(1-c_{k})^{2}\leq 0,$

for some $k\in\\{1,\dots,n\\}$. Without loss of generality, assume $k=1$.
Suppose $c_{1}=1$, i.e. $\rho_{1}=p_{1}$. Then, by Jensen’s inequality,

$\displaystyle\left(\sum_{i=1}^{n}p_{i}(1-c_{i})\right)^{2}-\sum_{i=1}^{n}p_{i}(1-c_{i})^{2}+\delta
C^{2}(1-c_{1})^{2}$
$\displaystyle\qquad\left(\sum_{i=1}^{n}p_{i}(1-c_{i})\right)^{2}-\sum_{i=1}^{n}p_{i}(1-c_{i})^{2}\leq
0.$

Suppose $c_{1}\neq 1$. If $c_{1}>1$, then there must be some $j\neq 1$ such
that $c_{j}<1$. Otherwise, $1=\sum_{i=1}^{n}p_{i}>\sum_{i=1}^{n}\rho_{i}=1$,
which is a contradiction. Similarly, if $c_{1}<1$, then there must be some
$j\neq 1$ such that $c_{j}>1$ by the same argument. In either case, we may
assume that $j=2$, and it follows that $(c_{1}-c_{2})^{2}>(1-c_{1})^{2}$. With
this inequality and the use of Lemma 4.3, it follows that

$\displaystyle\left(\sum_{i=1}^{n}p_{i}(1-c_{i})\right)^{2}-\sum_{i=1}^{n}p_{i}(1-c_{i})^{2}+\delta
C^{2}(1-c_{1})^{2}$
$\displaystyle\qquad=-\sum_{j=2}^{n}\frac{p_{j}\sum_{k=1}^{j-1}p_{k}}{\sum_{k=1}^{j}p_{k}}\left((1-c_{j})-\frac{\sum_{k=1}^{j-1}(1-c_{k})p_{k}}{\sum_{k=1}^{j-1}p_{k}}\right)^{2}+\delta
C^{2}(1-c_{1})^{2}$
$\displaystyle\qquad\leq-\frac{p_{2}p_{1}}{p_{1}+p_{2}}((1-c_{2})-(1-c_{1}))^{2}+\delta
C^{2}(1-c_{1})^{2}$
$\displaystyle\qquad=-\frac{p_{2}p_{1}}{p_{1}+p_{2}}(c_{1}-c_{2})^{2}+\delta
C^{2}(1-c_{1})^{2}$
$\displaystyle\qquad\leq\left(-\frac{p_{2}p_{1}}{p_{1}+p_{2}}+\delta
C^{2}\right)(1-c_{1})^{2}$ $\displaystyle\qquad\leq(-m^{2}+\delta
C^{2})(1-c_{1})^{2}\leq 0.$

Equality can only occur when $c_{i}$ are equal to $1$ for all $i$. This
implies that $p_{\Lambda|\Theta}(\tau|\theta)=\rho(\tau)$ for each
$\tau\in\mathcal{L}$, from which the conclusion $\theta=\theta_{*}$ follows.
Consequently, the term in (4.9) is also zero and therefore
$\frac{1}{2}\left\langle\nabla
V_{\delta}(\vartheta),\bar{g}(\vartheta)\right\rangle=0$ if and only if
$V_{\delta}(\vartheta)=0$. This completes the proof.

The following elementary result is used in the proof of Lemma 4.2.

###### Lemma 4.3

For any integer $n\geq 2$, real numbers
$\\{x_{i}\\}_{i=1}^{n}\subset\mathbb{R}$ and probability
$\\{p_{i}\\}_{i=1}^{n}$ on $\\{1,\dots,n\\}$,

$\displaystyle\left(\sum_{i=1}^{n}x_{i}p_{i}\right)^{2}-\sum_{i=1}^{n}x_{i}^{2}p_{i}=-\sum_{j=2}^{n}\frac{p_{j}\sum_{k=1}^{j-1}p_{k}}{\sum_{k=1}^{j}p_{k}}\left(x_{j}-\frac{\sum_{k=1}^{j-1}x_{k}p_{k}}{\sum_{k=1}^{j-1}p_{k}}\right)^{2}.$
(4.11)

Proof. When $n=2$, since $p_{1}+p_{2}=1$, it follows that the left-hand-side
of (4.11) is given by,

$\displaystyle(x_{1}p_{1}+x_{2}p_{2})^{2}-[x_{1}^{2}p_{1}+x_{2}^{2}p_{2}]$
$\displaystyle\quad=(x_{1}p_{1}+x_{2}(1-p_{1}))^{2}-[x_{1}^{2}p_{1}+x_{2}^{2}(1-p_{1})]$
$\displaystyle\quad=x_{1}^{2}p_{1}^{2}+2x_{1}x_{2}p_{1}(1-p_{1})+x_{2}^{2}(1-p_{1})^{2}-x_{1}^{2}p_{1}-x_{2}^{2}(1-p_{1})$
$\displaystyle\quad=x_{1}^{2}p_{1}(p_{1}-1)+2x_{1}x_{2}p_{1}(1-p_{1})-x_{2}^{2}p_{1}(1-p_{1})$
$\displaystyle\quad=p_{1}(1-p_{1})[-x_{1}^{2}+2x_{1}x_{2}-x_{2}^{2}]=-p_{1}(1-p_{1})(x_{1}-x_{2})^{2}.$

Additionally, the right-hand-side of (4.11) equals

$-\frac{p_{2}p_{1}}{\sum_{k=1}^{2}p_{k}}\left(x_{2}-\frac{x_{1}p_{1}}{p_{1}}\right)^{2}=-p_{1}p_{2}(x_{1}-x_{2})^{2}=-p_{1}(1-p_{1})(x_{1}-x_{2})^{2}.$

Thus, (4.11) holds when $n=2$.

Now suppose that (4.11) holds for some $n\geq 2$, then

$\displaystyle\left(\sum_{i=1}^{n+1}x_{i}p_{i}\right)^{2}$
$\displaystyle=\left(x_{n+1}p_{n+1}+(1-p_{n+1})\sum_{i=1}^{n}x_{i}\frac{p_{i}}{1-p_{n+1}}\right)^{2}$
$\displaystyle=x_{n+1}^{2}p_{n+1}+\left(\sum_{i=1}^{n}x_{i}\frac{p_{i}}{1-p_{n+1}}\right)^{2}(1-p_{n+1})$
$\displaystyle\qquad-
p_{n+1}(1-p_{n+1})\left(x_{n+1}-\sum_{i=1}^{n}x_{i}\frac{p_{i}}{1-p_{n+1}}\right)^{2},$

where the second equality comes from applying (4.11) with two terms.

On the other hand, by assumption, (4.11) holds for $n$, which implies that

$\displaystyle\left(\sum_{i=1}^{n}x_{i}\frac{p_{i}}{1-p_{n+1}}\right)^{2}$
$\displaystyle\qquad=\sum_{i=1}^{n}x_{i}^{2}\frac{p_{i}}{1-p_{n+1}}-\sum_{j=2}^{n}\frac{\frac{p_{j}}{1-p_{n+1}}\sum_{k=1}^{j-1}p_{k}}{\sum_{k=1}^{j}p_{k}}\left(x_{j}-\frac{\sum_{k=1}^{j-1}x_{k}p_{k}}{\sum_{k=1}^{j-1}p_{k}}\right)^{2}.$

Consequently,

$\displaystyle\left(\sum_{i=1}^{n}x_{i}\frac{p_{i}}{1-p_{n+1}}\right)^{2}(1-p_{n+1})$
$\displaystyle\qquad=\sum_{i=1}^{n}x_{i}^{2}p_{i}-\sum_{j=2}^{n}\frac{p_{j}\sum_{k=1}^{j-1}p_{k}}{\sum_{k=1}^{j}p_{k}}\left(x_{j}-\frac{\sum_{k=1}^{j-1}x_{k}p_{k}}{\sum_{k=1}^{j-1}p_{k}}\right)^{2}.$

Using the equality in the last display, we conclude that

$\displaystyle\left(\sum_{i=1}^{n+1}x_{i}p_{i}\right)^{2}$
$\displaystyle=x_{n+1}^{2}p_{n+1}+\sum_{i=1}^{n}x_{i}^{2}p_{i}-\sum_{j=2}^{n}\frac{p_{j}\sum_{k=1}^{j-1}p_{k}}{\sum_{k=1}^{j}p_{k}}\left(x_{j}-\frac{\sum_{k=1}^{j-1}x_{k}p_{k}}{\sum_{k=1}^{j-1}p_{k}}\right)^{2}$
$\displaystyle\qquad-
p_{n+1}(1-p_{n+1})\left(x_{n+1}-\sum_{i=1}^{n}x_{i}\frac{p_{i}}{1-p_{n+1}}\right)^{2}$
$\displaystyle=\sum_{i=1}^{n+1}x_{i}^{2}p_{i}-\sum_{j=2}^{n+1}\frac{p_{j}\sum_{k=1}^{j-1}p_{k}}{\sum_{k=1}^{j}p_{k}}\left(x_{j}-\frac{\sum_{k=1}^{j-1}x_{k}p_{k}}{\sum_{k=1}^{j-1}p_{k}}\right)^{2}.$

The proof is complete by mathematical induction.

###  4.2 Convergence of ergodic averages

In this subsection, the main theorem on the almost sure convergence of the
estimates of the extended design parameters $\\{\vartheta_{n}\\}$ of the AWH
algorithm is provided in Theorem 4.4.

###### Theorem 4.4

Let $\\{\vartheta_{n}\\}$ be given by (4.3), with target distribution $\rho$,
$\rho(\lambda)>0,\lambda\in\mathcal{L}$, $H$ a hyper-rectangle in
$\mathbb{R}^{|\mathcal{L}_{0}|}$, and $V_{\delta}$ be given by (4.8), where
$\delta>0$ is sufficiently small as in Lemma 4.2. Assume that $H$ is large
enough so that $\\{\vartheta\in H:V_{\delta}(\vartheta)=0\\}$ is non-empty.
Then $\\{\vartheta_{n}\\}$ converges almost surely to $\\{\vartheta\in
H:V_{\delta}(\vartheta)=0\\}$.

Proof. The proof of the statement in the theorem follows from almost identical
arguments to those in the proof of Theorem 3.3, with Lemma 3.1 and 3.2
replaced by Lemma 4.1 and 4.2, respectively. The main difference is that the
upper bound (3.7) of $h$ needs to be modified to the extended definition of
$h$. To this end (3.7) is replaced by

$\displaystyle|h(\vartheta,\xi)(u)|=\left\\{\begin{array}[]{ll}\left|1-\frac{\sum_{k=1}^{N_{I}}p_{\Lambda|X,\Theta}(u|\xi(k),\theta)}{N_{I}\rho(u)}\right|\leq
1+\frac{1}{\rho(u)},&\text{ if }u\in\mathcal{L},\\\
|\Phi(\xi,\theta)-\zeta|\leq 2\max\limits_{x\in\mathcal{X}}|\phi(x)|,&\text{
if }u=0.\end{array}\right.$

###### Corollary 4.5

Let $\\{\vartheta_{n}\\}$ be given by (4.3), with target distribution $\rho$,
$\rho(\lambda)>0,\lambda\in\mathcal{L}$, $H$ a hyper-rectangle in
$\mathbb{R}^{|\mathcal{L}_{0}|}$, and $V_{\delta}$ be given by (4.8), where
$\delta>0$ is sufficiently small as in Lemma 4.2. Assume that $H$ is large
enough so that $\\{\vartheta\in H:V_{\delta}(\vartheta)=0\\}$ is non-empty.
Then for a fixed $\lambda\in\mathcal{L}$ and a function
$\psi:\mathcal{X}\rightarrow\mathbb{R}$, the adaptive ergodic average

$\displaystyle\frac{\sum_{n=0}^{N-1}\sum_{k=1}^{N_{I}}\psi(X_{nN_{I}+k})p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})}{\sum_{n=0}^{N-1}\sum_{k=1}^{N_{I}}p_{X|\Lambda,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})}$

given in (4.2) converges almost surely to

$\sum_{x\in\mathcal{X}}\psi(x)p_{X|\Lambda}(x|\lambda).$

Proof. Recall that $\vartheta_{n}=(\theta_{n},\zeta_{n})$. By Theorem 4.4 the
sequence $\\{\vartheta_{n}\\}$ converges almost surely to the set
$\\{\vartheta\in H:V_{\delta}(\vartheta)=0\\}$. For any
$\vartheta\in\mathbb{R}^{|\mathcal{L}|}\times[-C,C]$ with
$V_{\delta}(\vartheta)=0$ it follows that

$\zeta=\sum_{x\in\mathcal{X}}\phi(x)p_{X,\Lambda|\Theta}(x,\lambda|\theta_{*}).$

Consequently, for any $\phi:\mathcal{X}\rightarrow\mathbb{R}$,

$\zeta_{N}=\hat{\phi}_{N}\doteq\frac{1}{N\,N_{I}}\sum_{n=0}^{N-1}\sum_{k=1}^{N_{I}}\phi(X_{nN_{I}+k})p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})$

converges almost surely to

$\zeta=\sum_{x\in\mathcal{X}}\phi(x)p_{X,\Lambda|\Theta}(x,\lambda|\theta_{*}).$

Moreover, by taking $\phi=\psi$ and $\phi\equiv 1$, and since the ratio of two
almost surely convergent sequences is still an almost surely convergent
sequence, we know that

$\displaystyle\frac{\sum_{n=0}^{N-1}\sum_{k=1}^{N_{I}}\psi(X_{nN_{I}+k})p_{\Lambda|X,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})}{\sum_{n=0}^{N-1}\sum_{k=1}^{N_{I}}p_{X|\Lambda,\Theta}(\lambda|X_{nN_{I}+k},\theta_{n})}=\frac{\hat{\psi}_{N}}{\hat{1}_{N}}$

converges almost surely to

$\frac{\sum_{x\in\mathcal{X}}\psi(x)p_{X,\Lambda|\Theta}(x,\lambda|\theta_{*})}{\sum_{x\in\mathcal{X}}p_{X,\Lambda|\Theta}(x,\lambda|\theta_{*})}=\sum_{x\in\mathcal{X}}\psi(x)p_{X|\Lambda}(x|\lambda).$

This completes the proof.

##  Appendix A Convergence of Stochastic Approximations

In order to make this paper self-contained, some notations and statements used
in [18] are introduced in this Appendix. Additionally, Theorem 6.1 from [18,
Chapter 6] is stated, which is the main result that will be applied to show
the almost sure convergence of the AWH algorithm in the proof of Theorem 3.3
in this paper.

Consider a stochastic approximation with projection onto a constraint set $H$
defined as follows:

$\displaystyle\theta_{n+1}=\theta_{n}+\varepsilon_{n}Y_{n}+\varepsilon_{n}Z_{n},$
(A.1)

where $Y_{n}$ is the observation at time $n$ and $Z_{n}$ is the projection
term. Define $\mathcal{F}_{n}$ as the $\sigma$-algebra determined by the
initial condition $\theta_{0}$ and observations $Y_{i},i<n$, and assume that
there exists measurable functions $g_{n}(\cdot,\cdot)$ such that

$\displaystyle E_{n}Y_{n}=g_{n}(\theta_{n},\xi_{n})+\beta_{n},$ (A.2)

where $E_{n}$ denotes the expectation conditioned on $\mathcal{F}_{n}$,
$\xi_{n}$ denotes the correlated driving noise, and $\beta_{n}$ represents a
bias which is asymptotically unimportant.

Next, we define martingale differences $\delta M_{n}\doteq Y_{n}-E_{n}Y_{n}$
and use (A.2) to write (A.1) as

$\displaystyle\theta_{n+1}=\theta_{n}+\varepsilon_{n}[g_{n}(\theta_{n},\xi_{n})+\delta
M_{n}+\beta_{n}+\varepsilon_{n}Z_{n}].$ (A.3)

Moreover, define $t_{0}=0$ and $t_{n}=\sum_{i=0}^{n-1}\varepsilon_{i}$. For
$t\geq 0$, let $m(t)$ denote the unique value of $n$ such that $t_{n}\leq
t<t_{n+1}$.

With these notations, we can now state Theorem 6.1 from [18, Chapter 6]. We
only state the part that is relevant for the present paper.

###### Theorem A.1

Consider $\\{\theta_{n}\\}$ defined by (A.3). Assume

  * •

(5.1.1) $\sum_{n=0}\varepsilon_{n}=\infty,\varepsilon_{n}\geq
0,\varepsilon_{n}\rightarrow 0,\mbox{ for }n\geq 0.$

  * •

(A1.1) $\sup_{n}E|Y_{n}|<\infty.$

  * •

(A1.2) $g_{n}(\theta,\xi)$ is continuous in $\theta$ for each $\xi$ and $n$.

  * •

(A1.4) For each $\mu>0$ and some $T>0$,

$\lim_{n\rightarrow\infty}P\left\\{\sup_{j\geq n}\max_{0\leq t\leq
T}\left|\sum_{i=m(jT)}^{m(jT+t)-1}\varepsilon_{i}\delta
M_{i}\right|\geq\mu\right\\}=0.$

  * •

(6.2) $P(\xi_{n+1}\in\cdot|\xi_{i},\theta_{i},i\leq
n)=p(\xi_{n},\cdot|\theta_{n})$, where $p(\xi,\cdot|\theta)$ denotes the one-
step transition probability with starting point $\xi$ and parameterized by
$\theta$.

  * •

(A6.1) There is a continuous function $\bar{g}(\cdot)$ such that for
$\theta\in H$, the expression

$v_{n}(\theta,\xi_{n})=\sum_{i=n}^{\infty}\varepsilon_{i}E_{n}[g_{i}(\theta,\xi_{i}(\theta))-\bar{g}(\theta)]$

is well defined when the initial condition for $\\{\xi_{i},i\geq n\\}$ is
$\xi_{n}(\theta)=\xi_{n}$, and $v_{n}(\theta_{n},\xi_{n})\rightarrow 0$ w.p.1.

  * •

(A6.2) The asymptotic rates of change of the processes

$A^{0}(t)=\sum_{i=0}^{m(t)-1}\alpha_{i}\mbox{ and
}N^{0}(t)=\sum_{i=0}^{m(t)-1}\delta N_{i}$

are zero with probability one, where
$\alpha_{n}=v_{n+1}(\theta_{n+1},\xi_{n+1})-v_{n+1}(\theta_{n},\xi_{n+1})$ and
$\delta
N_{n}=v_{n+1}(\theta_{n},\xi_{n+1})-E_{n}v_{n+1}(\theta_{n},\xi_{n+1})$. In
other words, for $\omega$ not in some null set and for any positive $T$,

$\lim_{n}\sup_{j\geq n}\max_{0\leq t\leq
T}|A^{0}(\omega,jT+t)-A^{0}(\omega,jT)|=0.$

and

$\lim_{n}\sup_{j\geq n}\max_{0\leq t\leq
T}|N^{0}(\omega,jT+t)-N^{0}(\omega,jT)|=0.$

and the following constraint set condition

  * •

(A4.3.1) $H$ is a hyperrectangle. In other words, there are real numbers
$a_{i}<b_{i},i=1,\dots,r$, such that $H=\\{x:a_{i}\leq x_{i}\leq b_{i}\\}$.

Then for almost all $\omega$, $\\{\theta_{n}(\omega)\\}$ converges to some
limit set of the ODE

$\dot{\theta}=\bar{g}(\theta)+z,z\in-C(\theta),$

where $z$ is the minimum force needed to keep the solution in $H$.

###### Remark A.2

According to the last paragraph on p.137, for some positive number $T$,

$\displaystyle\lim_{n}\sup_{j\geq n}\max_{0\leq t\leq
T}|M^{0}(jT+t)-M^{0}(jT)|=0\text{ w.p.1.}$ (A.4)

is guaranteed by (A2.1) $\sup_{n}E|Y_{n}|^{2}<\infty$ and (A2.4)
$\sum_{i}\varepsilon_{i}^{2}<\infty$, where

$M^{0}(t)=\sum_{i=0}^{m(t)-1}\varepsilon_{i}\delta M_{i},\quad\delta
M_{n}=Y_{n}-E_{n}Y_{n}.$

Observe that (A.4) can be written as

$\lim_{n}\sup_{j\geq n}\max_{0\leq t\leq
T}\left|\sum_{i=m(jT)}^{m(jT+t)-1}\varepsilon_{i}\delta M_{i}\right|=0\text{
w.p.1},$

which implies (A1.4).

## Acknowledgements

Financial support for Guo-Jhen Wu from the Brummer & Partners MathDataLab and
the Swedish e-Science Research Center (SeRC) Data Science MCP is gratefully
acknowledged. The authors would like to thank Jack Lidmar for introducing us
to the accelerated weight histogram algorithm and numerous discussions on
related topics.

## References

  * [1] Asmussen, S. and Glynn, P.W.: _Stochastic simulation_ , Springer, New York, 2007.
  * [2] Andrieu, C., de Freitas, N., Doucet, A., and Jordan, M.I.: An introduction to MCMC for machine learning. _Machine Learning_ , 50(1), 5–43, 2003.
  * [3] Andrieu, C., Moulines, E., and Priouret, P.: Stability of stochastic approximations under verifiable conditions. _SIAM J. Control Optim._ 44(1), 283-312, 2005.
  * [4] Belardinelli, R.E. and V. D. Pereyra, V.D.: Fast algorithm to calculate density of states. _Phys. Rev. E_ 75, 046701, 2007.
  * [5] Berg, B.A.: Introduction to multicanonical Monte Carlo simulations, _Monte Carlo Methods_ , Fields Inst. Commun. 26, N. Madras (ed.), (American Mathematical Society) p.1-24, 2000.
  * [6] Berg, B.A. and Celik, T.: New approach to spin-glass simulations, _Phys. Rev. Lett._ 69, 2292, 1992.
  * [7] Berg, B.A. and Neuhaus, T.: Multicanonical algorithms for first order phase transitions, _Phys. Lett._ B267, p. 249, 1991.
  * [8] Berg, B.A. and Neuhaus, T.: Multicanonical ensemble: a new approach to simulate first-order phase transitions, _Phys. Rev. Lett._ 68, 9-12, 1992.
  * [9] Dupuis, P., Liu, Y., Plattner, N., and Doll, J.D.: On the infinite swapping limit for parallel tempering, _Multiscale Model. Simul._ , 10(3), 986–1022, 2012.
  * [10] Doll, J.P., Plattner, N., Freeman, D.L., Liu, Y., and Dupuis, P.: Rare-event sampling: occupation-based performance measures for parallel tempering and infinite swapping, _J. Chem. Phys._ , 137, 2012.
  * [11] Dupuis, P., Wu, G.-J.: Large deviation properties of the empirical measure of a metastable small noise diffusion, _J. Theor. Probab._ , 2021.
  * [12] Dupuis, P., Wu, G.-J.: Analysis and optimization of certain parallel Monte Carlo methods in the low temperature limit, _Submiited_ , 2021.
  * [13] Fort, G., Jourdain, B., Kuhn, E., Leliévre, T., and Stoltz, G.: Convergence of the Wang-Landau algorithm, _Mathematics of Computation_ , 84(295), 2297-2327, 2015.
  * [14] Gilks, W.R., Richardson, S., and Spiegelhalter, D.J.: _Markov chain Monte Carlo in practice_ , Chapman and Hall, London, 1996.
  * [15] Hansmann, U.H.E., and Okamoto, Y.: New Monte Carlo algorithms for protein folding, _Current Opinion in Structural Biology_ 9, 177-183, 1999.
  * [16] Hooft, R.W.W., van Eijck, B.P., and Kroon, J.: An adaptive umbrella sampling procedure in conformational analysis using molecular dynamics and its application to glycol, _J. Chem. Phys._ 97, 6690-6694, 1992.
  * [17] Kushner, H. and Dupuis, P.: Constrained stochastic approximation via the theory of large deviations. _Adaptive statistical procedures and related topics_ , Institute of Mathematical Statistics, Hayward, CA, pp. 253–265, 1986.
  * [18] Kushner, H. and Yin, G.: _Stochastic approximation and recursive algorithms and applications_ , Springer, New York, 2003.
  * [19] Lelievre, T., Rousset, M., and Stoltz, G.: Computation of free energy profiles with adaptive parallel dynamics, _J. Chem. Phys._ 126, 134111, 2007.
  * [20] Lidmar, J.: Improving the efficiency of extended ensemble simulations : The accelerated weight histogram method, _Physical Review E. Statistical, Nonlinear, and Soft Matter Physics_ , 85(5), 056708, 2012.
  * [21] Lindahl, V., Lidmar, J., and Hess, B.: Accelerated weight histogram method for exploring free energy landscapes. _J. Chem. Phys._ , 141. 10.1063/1.4890371, 2014.
  * [22] Lindahl, V., Lidmar, J., and Hess, B.: Riemann metric approach to optimal sampling of multidimensional free-energy landscapes, _Physical review. E_ , 98(2), 2018.
  * [23] Ljung, L.: Analysis of recursive stochastic algorithms, _IEEE Trans. Automatic Control_ , 22, 551-575, 1977.
  * [24] Marinari, E. and Parisi, G.: Simulated tempering: a new Monte Carlo scheme, _Europhys. Lett._ 19, 451, 1992.
  * [25] Mezei, M.: Adaptive umbrella sampling: Self-consistent determination of the non-Boltzmann bias, _J. Comp. Phys._ 68, 237-248, 1987.
  * [26] Noé, F., Olsson, S., Köjler, J., and Wu, H.: Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning, _Science,_ Vol. 365, Issue 6457, eaaw1147, 06 Sep 2019.
  * [27] Robert, C.P. and Casella, G.: _Monte Carlo statistical methods_ , Springer, New York, 2004.
  * [28] Robert, C.P. and Rosenthal, J.: General state space Markov chains and MCMC algorithms, _Probab. Surv._ , 1, 20-71, 2004.
  * [29] Swendsen, R.H. and Wang, J.S.: Replica Monte Carlo simulation of spin-glasses, _Phys. Rev. Lett._ 57, 2607–2609, 1986.
  * [30] Urakami, N. and Takasu, M.: Multicanonical Monte Carlo simulation of a polymer with stickers, _J. Phys. Soc. Jpn_ 65, p. 2694, 1996.
  * [31] Wang, F. and Landau, D.P.: Efficient, multiple-range random walk algorithm to calculate the density of states. _Phys. Rev. Lett._ 86(2050), 2001.
  * [32] Zhang, C. and Ma, J.: Simulation via direct computation of partition functions, _Physical Review E_ , 76, 036708, 2007.


#  $J/\psi$ production at high $p_{T}$ in $p+p$ and $A+A$ collisions at STAR

Zebo Tang111The author was supported in part by the National Natural Science
Foundation of China under Grant No. 10610285, 10610286, 10575101 and the
Knowledge Innovation Project of Chinese Academy of Sciences under Grant No.
KJCX2-YW-A14. (for the STAR collaboration)  Dept. of Modern Physics,
University of Science and Technology of China, Hefei, Anhui, China, 230026;
Brookhaven National Laboratory, Upton, NY, 11973, USA<EMAIL_ADDRESS>
###### Abstract

The preliminary results of $J/\psi$ spectra at high transverse momentum
($5<p_{T}<14$ GeV/c) in $p+p$ and Cu+Cu collisions at $\sqrt{s_{\mathrm{NN}}}$
= 200 GeV are reported. The nuclear modification factor is measured to be
$0.9\pm 0.2$ at $p_{T}>5$ GeV/c. The correlations between $J/\psi$ and charged
hadrons are also studied in $p+p$ collisions to understand the $J/\psi$
production mechanism at high $p_{T}$.

##  1 Introduction

$J/\psi$ dissociation from color-screening of Quantum Chromodynamics (QCD) in
a Quark-Gluon Plasma (QGP) is one of the major signatures of QCD de-
confinement in relativistic heavy-ion collisions. Recent calculations which
assuming the AdS/CFT duality is valid for QCD expect that a heavy fermion pair
bound state (an analog of quarkonium in QCD) will have an effective
dissociation temperature decreasing with $p_{T}$ [1]. This requires a
measurement of $J/\psi$ extending to $p_{T}>5$ GeV/c where the effective
$J/\psi$ dissociation temperature is expected to decrease to the temperature
reached at RHIC ($\sim$ 1.5 $T_{c}$) . In this paper, we report the $J/\psi$
spectra at high transverse momentum ($5<p_{T}<14$ GeV/c) in $p+p$ and Cu+Cu
collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200 GeV. In addition, we performed an
analysis of $J/\psi$-hadron correlations in $p+p$ collisions to understand the
$J/\psi$ production mechanism at high $p_{T}$. The technique is similar to
that used by UA1 [2] and dihadron correlations analyzed by STAR [3].

UA1 simulated $J/\psi$-hadron correlation and found two cases: When a $J/\psi$
originated from $\chi_{c}$ there as no visible near-side correlation, whereas
$J/\psi$’s originating from $B$ meson decays showed a strong near-side
correlation. The large acceptance of STAR Time Projection Chamber (TPC) [4]
and the Barrel Electromagnetic Calorimeter (BEMC) [4] covering $|\eta|<1$ are
very well suited for such analyzes.

##  2 Data analysis and Results

At STAR, both the TPC and BEMC can provide electron identification [4]. At
high $p_{T}$, the BEMC is very powerful for electron identification and can
also be used to set up a fast trigger to enrich the electron sample. At
moderate $p_{T}$, the TPC can identify electrons efficiently. In this
analysis, the high $p_{T}$ $J/\psi$ was reconstructed through the dielectron
channel, one electron at high $p_{T}$ identified by combing the BEMC and TPC
and the other electron at lower $p_{T}$ identified by the TPC only. We used
the BEMC triggered data in $p+p$ and Cu+Cu collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 200 GeV. The integrated luminosity is $\sim$ 2.8
(11.3) $pb^{-1}$ for $p+p$ collisions collected in year 2005 (2006) with
transverse energy threshold $E_{T}>$ 3.5 (5.4) GeV, and $\sim$ 860 $\mu
b^{-1}$ for Cu+Cu collisions collected in year 2005 with $E_{T}>$ 3.75 GeV.

Figure 1: The dielectron invariant mass distributions in $p+p$ (left) and
Cu+Cu (right) collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200 GeV.

Figure 1 shows the high $p_{T}$ $J/\psi$ signal in $p+p$ (left) and Cu+Cu
(right) collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200 GeV. The background is
represented by the dashed lines from like-sign technique. We applied a cut of
$p_{T}$ $>2.5-4$ GeV/c to the EMC triggered electrons and the cut of $p_{T}$
$>1.2-1.5$ GeV/c for lower $p_{T}$ electrons. This ensured clean $J/\psi$
identification. The signal/background (S/B) ratio in the analysis is 22/2
(40/14) in $p+p$ collisions using year 2005 (2006) data and 17/23 in Cu+Cu
collisions. The $p_{T}$ coverage in $p+p$ and Cu+Cu collisions taken in year
2005 is $5<p_{T}<8$ GeV/c, while in $p+p$ collisions taken in year 2006, the
$J/\psi$ $p_{T}$ can reach 14 GeV/c due to higher recorded luminosity and full
BEMC coverage. The $J/\psi$ invariant cross section $B_{ee}\times
Ed^{3}\sigma/dp^{3}$, after efficiency correction, are shown as symbols in
Figure 2 (left).

Figure 2: Left: $J/\psi$ invariant cross section as a function of $p_{T}$ in
$p+p$ and Cu+Cu collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200 GeV. Errors shown
are statistical only. Middle: $x_{T}$ scaling of pions, protons and $J/\psi$s.
The data from other measurements can be found in references [8, 9, 10, 12, 5,
2, 11]. Right: $J/\psi$ $R_{AA}$ as a function of $p_{T}$. The dot-dashed line
represents the fit by constant to all the data points at $5<p_{T}<10$ GeV/c.
The boxes on the right show the normalization uncertainty.

The invariant cross section of inclusive pion and proton production in high
energy $p+p$ collisions have been found to follow the $x_{T}$ scaling law:
$E\frac{d^{3}\sigma}{dp^{3}}=\frac{g(x_{T})}{\sqrt{s}^{n}}$, where
$x_{T}=2p_{T}/\sqrt{s}$. The value of the power $n$ depends on the quantum
exchanged in the hard scattering and is related to the number of point-like
constituents taking an active role in parton model. It reaches 8 in the case
of a diquark scattering and reaches 4 in more basic scattering processes (as
in QED). Figure 2 (middle) shows the $x_{T}$ scaling of $J/\psi$, pion and
proton. The power $n$ was found to be $6.5\pm 0.8$ for pion and proton [5] and
$5.6\pm 0.2$ for $J/\psi$, which indicates that the high $p_{T}$ $J/\psi$
production mechanism is closer to parton-parton scattering.

Figure 2 (right) shows the $J/\psi$ nuclear modification factor $R_{AA}$ as a
function of $p_{T}$ in 0-20% and 0-60% Cu+Cu from PHENIX [6] and STAR
measurements. $R_{AA}$ tends to increase from low to high $p_{T}$, although
the error bars currently do not allow to draw strong conclusions. One can
nevertheless do a combined fit to all the high-pt data and find that
$R_{AA}=0.9\pm 0.2$. This result is in contrast to the expectation from
AdS/CFT-based models [1] and from the Two-Component model [7] which predict a
decreasing $R_{AA}$ with increasing $p_{T}$. This result could indicate that
other $J/\psi$ production mechanisms such as virtual photons or formation time
[13] play a role at high $p_{T}$.

With large S/B ratios, the $J/\psi$-hadron correlations were also measured in
$p+p$ collisions. Figure 3 (left) shows the azimuthal angle correlations
between high $p_{T}$ $J/\psi$ ($p_{T}>5$ GeV/c) and charged hadrons. No
significant near side correlations were observed, which is in contrast to the
dihadron correlation measurements [3]. Since the Monte Carlo simulation
results show a strong near side correlation if the $J/\psi$ is produced from
$B$-meson decay, these results can be used to constrain the $B$-meson
contribution to $J/\psi$ production. Figure 3 (right) shows the associated
charged hadron $p_{T}$ distribution on the near side and away side with
respect to $J/\psi$ triggers and charged hadron triggers. On the away side,
the yields of the associated charged hadrons with respect to both kinds of
triggers are consistent with each other, which indicates that the hadrons on
the away side of $J/\psi$ triggers are from light quark or gluon
fragmentation. On the near side, the associated charged hadron yields with
respect to $J/\psi$ triggers are significantly lower than those with respect
to charged hadron triggers. This indicates that the $B$-meson is not a
dominant contributor to the inclusive high $p_{T}$ $J/\psi$.

Figure 3: Left: $J/\psi$-hadron correlations after background subtraction in
$p+p$ collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200 GeV. Right: Associated
charged hadron $p_{T}$ distribution on the near and away side with respect to
$J/\psi$ triggers and charged hadron triggers.

##  3 Summary

We reported the STAR preliminary results of $J/\psi$ spectra from 200 GeV
$p+p$ and Cu+Cu collisions at high $p_{T}$ ($5<p_{T}<14$ GeV/c) at mid-
rapidity through the dielectron channel. The high $p_{T}$ $J/\psi$ production
was found to follow the $x_{T}$ scaling with a beam energy dependent factor
$\sim$ $\sqrt{s_{\mathrm{NN}}}$ 5.6±0.2. The $J/\psi$ nuclear modification
factor $R_{AA}$ in Cu+Cu increases from low to high $p_{T}$ which challenges
some models. The average of $R_{AA}$ at $p_{T}$ $>$ 5 GeV/c is $0.9\pm 0.2$,
consistent with no $J/\psi$ suppression. It implies that high $p_{T}$ $J/\psi$
may be produced from virtual photon or formed outside of the hot interaction
region [13]. The $J/\psi$-hadron correlations were also discussed. We observed
an absence of charged hadrons accompanying high $p_{T}$ $J/\psi$ on the near
side which indicates that the $B$-meson is not a dominant contributor to the
inclusive high $p_{T}$ $J/\psi$.

## References

## References

  * [1] H. Liu, K. Rajagopal and U.A.Wiedemann, Phys. Rev. Lett. 98 (2007) 182301.
  * [2] UA1 Collaboration, Phys. Lett. B 256 (1991) 112.
  * [3] J. Adams, et al., (STAR Collaboration), Phys. Rev. Lett. 95 (2005) 152301.
  * [4] M. Anderson, et al., Nucl. Instrum. Meth. A 499 (2003) 659; M. Beddo, et al., Nucl. Instrum. Meth. A 499 (2003) 725; B.I. Abelev, et al., (STAR Collaboration) Phys. Rev. Lett. 98 (2007) 192301.
  * [5] J. Adams, et al. (STAR Collaboration), Phys. Lett. B 616, (2005) 8; J. Adams, et al., (STAR Collaboration), Phys. Lett. B 637 (2006) 161.
  * [6] A.Adare, et al., (PHENIX Collaboration), arXiv:0801.0220.
  * [7] X. Zhao and R. Rapp, arXiv:0712.2407, private communication for Cu+Cu.
  * [8] F. Abe, et al., (CDF Collaboration), Phys. Rev. Lett. 69 (1992) 3704; D. Acosta et al., (CDF Collaboration), Phys. Rev. D 71 (2005) 032001.
  * [9] D. Abtreasyan, et al., Phys. Rev. D 19 (1979) 764.
  * [10] B.Alper, et al., (British-Scandinavian Collaboration) Nucl. Phys. B 100 (1975) 237; C. Kourkoumelis, et al., Phys. Lett. 91B (1980) 481.
  * [11] M. Banner, et al., (UA2 Collaboration), Phys. Lett. B 115 (1982) 59.
  * [12] A.Adare, et al., (PHENIX Collaboration), Phys. Rev. Lett. 98 (2007), 232002.
  * [13] F. Karsch and R. Petronzio, Phys. Lett. B 193 (1987), 105; J.P. Blaizot and J.Y. Ollitrault, Phys. Lett. B 199 (1987), 499