DOI HERE

Fei Ding et al.

[$\dagger$]Co-first author.

# Channelling Multimodality Through a Unimodalizing Transport: Warp-U Sampler
and Stochastic Bridge Sampling

Fei Ding  David E. Jones  Shiyuan He  Xiao-Li Meng  Department of Statistics,
Texas A&M University, College Station, 77840, TX, USA
Institute of Statistics and Big Data, Renmin University of China, Beijing,
100872, China
Department of Statistics, Harvard University, Cambridge, 02138, MA, USA.

###### Abstract

Monte Carlo integration is fundamental in scientific and statistical
computation, but requires reliable samples from the target distribution, which
poses a substantial challenge in the case of multi-modal distributions.
Existing methods often involve time-consuming tuning, and typically lack
tailored estimators for efficient use of the samples. This paper adapts the
Warp-U transformation (Wang et al., 2022) to form multi-modal sampling
strategy called Warp-U sampling. It constructs a stochastic map to transport a
multi-modal density into a uni-modal one, and subsequently inverts the
transport but with new stochasticity injected. For efficient use of the
samples for normalising constant estimation, we propose (i) an unbiased
estimation scheme based coupled chains, where the Warp-U sampling is used to
reduce the coupling time; and (ii) a stochastic Warp-U bridge sampling
estimator, which improves its deterministic counterpart given in Wang et al.
(2022). Our overall approach requires less tuning and is easier to apply than
common alternatives. Theoretically, we establish the ergodicity of our
sampling algorithm and that our stochastic Warp-U bridge sampling estimator
has greater (asymptotic) precision per CPU second compared to the Warp-U
bridge estimator of Wang et al. (2022) under practical conditions. The
advantages and current limitations of our approach are demonstrated through
simulation studies and an application to exoplanet detection.

###### keywords:

adaptive MCMC, astrostatistics, bridge sampling, Bayesian evidence, multi-
modal density, normalising constant estimation.

##  1 Monte Carlo Integration

###  1.1 Computing normalising constants — How hard can it be?

Computing or estimating normalising constants may sound like a routine
homework exercise, but it is a surprisingly wide-ranging and hard problem in
statistics and more broadly in many scientific applications, from computing
the free energy in physics (see Bennett, 1976) to cognitive studies in
psychology (see Gronau et al., 2019, 2020). In particular, Bayesian evidence
ratios (a.k.a., Bayes factors) are widely used for hypothesis testing and
model selection in many scientific fields. For example, Nelson et al. (2020)
and Pullen and Morris (2014) discussed the computation of Bayesian evidence in
the context of exoplanet detection and systems biology, respectively. Given
the many scientific uses of normalising constants, computationally and
statistically efficient methods for estimating them are of high practical
value, and many powerful algorithms have been (re)invented, as we shall review
briefly at the end of this article in Section 6.1. On the other hand, some
ubiquitous scenarios remain challenging, such as estimating normalising
constants for multi-modal target densities. For example, Nelson et al. (2020)
applied numerous strategies for estimating the Bayesian evidence for the
presence of an exoplanet orbiting a star and obtained somewhat divergent
estimates, even after substantial calibration efforts.

In general, the problem of computing normalising constants is exceedingly
challenging and requires methods tailored to the particular context because if
there were an effective generic method for it, then we would have a general
efficient method for computing integration, which is impossible with
contemporary software, hardware, and mindware111A term coined by psychologist
David Perkins, referring to the knowledge and problem solving techniques
available to human minds.. Specifically, for any (absolutely integrable)
function $f$ we have

$\int_{\Theta}f(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})=\int_{\Theta}f_{+}(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})-\int_{\Theta}f_{-}(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta}),$
(1)

where $\mu$ is a suitable baseline measure222See Kong et al. (2003) for the
essential role the baseline measure $\mu$ plays in forming an appropriate
likelihood theory for Monte Carlo integration., and both
$f_{+}=\max\\{f(\boldsymbol{\theta}),0\\}$ and
$f_{-}=-\min\\{f(\boldsymbol{\theta}),0\\}$ are non-negative functions and
hence can be viewed as un-normalised densities on (possibly a subset of)
$\Theta$. Consequently, although we will frame our approaches in the context
of Monte Carlo estimation of normalising constants, they are generally
applicable because of identity (1), and because of other ways of making the
integrand positive, such as using a lower bound on $f$ as a control variate.

We consider the following general problem: for an unnormalised probability
density $q$ with support $\Theta$, we seek to estimate the normalising
constant
$c=\int_{\Theta}q(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})$. Any
Monte Carlo approach for this problem involves two key tasks: (i) obtaining
samples
$\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2},\cdots,\boldsymbol{\theta}_{n}$
from $\pi=q/c$ or some related density (or densities), and (ii) constructing
an estimate of $c$ based on the Monte Carlo samples. These tasks may be
performed sequentially or in combination, depending on the specific strategy.
It is worth emphasising that when $\pi$ is multi-modal, substantial
inefficiencies can result from not addressing the difficulties created by
multi-modality in both (i) and (ii) (on the other hand addressing multi-
modality does not automatically guarantee an efficient estimator of $c$).

###  1.2 An Integrated Approach for Combating Multi-modality

We address the multi-modality challenge by an effective Warp-U sampler
designed for multi-modal target densities, followed by a complementary
strategy for estimating the corresponding target normalising constants. The
estimation step is complementary in the sense that much of the necessary
computation is completed during the sampling stage. Both our sampling and
estimation steps make use of the stochastic Warp-U transformation proposed by
Wang et al. (2022), which transforms multi-modal densities into approximately
uni-modal ones. The key conceptual step of constructing a Warp-U
transformation is to induce a mixture representation of the target
$\pi(\boldsymbol{\theta})=\sum_{k=1}^{K}w_{k}\pi_{k}(\boldsymbol{\theta})$,
where intuitively $\pi_{k}$ can be viewed as representing a “sub-population
component” that is responsible for the $k$th mode in the target density333It
is worthwhile, however, to remind ourselves that the relationship between
multi-modality and mixture distributions is a lose one. For example, a mixture
distribution can easily exhibit a single mode when the mixture components are
sufficiently close to each other.. We will review the construction details in
Section 2.1.

The remainder of Section 2 details the sampling step and estimation step of
our integrated approach. Section 2.3 develops a Markov chain Monte Carlo
(MCMC) sampling algorithm, the Warp-U sampler, which alternately applies the
stochastic Warp-U transformation and its inverse. The key to this sampler is
that the components (i.e. the $\pi_{k}$’s) of the target are stochastically
transported to an intermediate uni-modal density, then the latter gets
stochastically transported back to the original target. Since all the target
components are mapped to similar intermediate uni-modal densities, the step of
mapping back to the target density tends to have a high probability of mapping
any given draw to a component different from that which was originally
sampled. This property makes the proposed procedure effective for multi-mode
exploration, as initially demonstrated in the first version444 See Section 6
of the initial version at https://arxiv.org/pdf/1609.07690v1.pdf, which was
removed in the published version, following an editorial request to streamline
Wang et al. (2022). of Wang et al. (2022). Our full algorithm also
incorporates a classical Metropolis-Hastings update to ensure irreducibility,
and our simulation results given in Section 4 demonstrate that it can be
substantially more efficient than competing approaches for sampling from
multi-modal targets, such as parallel tempering (Geyer, 1991).

The quality of the Warp-U transformation obviously affects the convergence
behaviour of the Warp-U sampler. In the estimation strategy proposed by Wang
et al. (2022), a Gaussian mixture approximation is fit to the target density
$\pi$, which is then used to build a Warp-U map. But in the sampling context,
initially there are no Monte Carlo samples with which to fit any kind of
mixture. To address this issue, we adopt an adaptive Warp-U sampler which
starts by sampling from a uniform density, or any other over-dispersed crude
approximation to the target density, and then proceeds in stages, updating the
approximating density at each stage. Importantly, as with the Warp-U bridge
sampling estimation strategy of Wang et al. (2022), it is unnecessary for the
approximating density to be a precise representation of the target density in
order to substantially improve MCMC convergence speed; any reasonable
approximation is sufficient to generate practically meaningful acceleration.
This means that the number of stages needed in our algorithm is typically
small, e.g., fewer than ten. Later in Section 3.1, we prove that the adaptive
Warp-U sampler is ergodic, following the approach of Roberts and Rosenthal
(2007).

Given MCMC samples from $\pi$, we are often interested in estimating
$\mathbb{E}_{\pi}h(\boldsymbol{\theta})$ for a chosen function $h$, where the
expectation is taken with respect to $\boldsymbol{\theta}\sim\pi$. Traditional
methodology requires running a long Markov chain and performing convergence
diagnostics to guarantee accurate estimation. In Section 2.4, we devise an
unbiased estimation strategy for the Warp-U sampler, follow the strategy
proposed and tested by Jacob et al. (2020). The advantage of this unbiased
approach is its direct assessment of chain convergence, both theoretically and
practically, through the $L$-coupling method, as elaborated in Biswas et al.
(2019), and further investigated in Craiu and Meng (2022). Under this scheme,
two Warp-U sampling chains run concurrently and are prompted to meet. In
particular, we use optimal transport coupling to keep the two chains close
during their transitions between density modes. This implementation ensures
that we obtain genuinely unbiased moment estimates under multi-modal
distributions, regardless of the isolation of the density modes.

As another key ingredient to our estimation step, Section 2.5 proposes a
divide-and-conquer strategy to break a computational bottleneck that appears
in the Warp-U bridge estimators of Wang et al. (2022). Warp-U bridge
estimators apply standard bridge sampling estimators (Bennett, 1976; Meng and
Wong, 1996) to the (approximately) uni-modal $\tilde{q}$ resulting from the
Warp-U transformation, which preserves the normalising constant, i.e.,
$\int_{\Theta}\tilde{q}(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})=c=\int_{\Theta}q(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})$.
The bottleneck is created by the large number of “cross-component” function
evaluations needed for implementing bridging sampling estimators using
$\tilde{q}$ and all the Warp-U transformed draws. In Section 3.2 we prove
that, in comparison with the Warp-U bridge estimator, a substantial reduction
in asymptotic variance can be obtained by instead applying standard bridge
sampling to estimate the normalising constant of each mixture component
$q_{k}=c\pi_{k}$, for $k=1,\dots,K$, and then combining the results to
estimate $c$. The simulation studies reported in Section 4.3 show that, given
fixed computational resources, this divide-and-conquer approach incurs
substantially lower root mean squared error (RMSE) than the corresponding ones
given in Wang et al. (2022).

In Section 5, we apply our integrated approach to estimating the Bayesian
evidence for exoplanets based on radial velocity (RV) data, which demonstrates
its potential for outperforming existing sampling and estimation strategies
applied together in various combinations, including combining parallel
tempering with the Warp-U bridge estimator of Wang et al. (2022). In order to
report our main proposals as soon as possible, we defer the literature review
to Section 6, where we also discuss the limitations of our approach and
possible improvements and extensions, such as how the divide-and-conquer
method might be combined with the unbiased estimation strategy to reduce both
variance and bias. Proofs are given in the online Supplementary Material.

##  2 Warp-U Sampling and Estimation Methods

###  2.1 Bridge Sampling Estimation and Warp-U Transformations

Let $q_{1}$ and $q_{2}$ denote unnormalised densities with unknown normalising
constants $c_{1}$ and $c_{2}$, respectively, and for simplicity we assume they
share the same support $\Theta=\mathbb{R}^{d}$. We are interested in
estimating their ratio $r=c_{1}/c_{2}$ (e.g., a Bayes factor). For this
scenario, bridge sampling (Bennett, 1976; Meng and Wong, 1996) relies on the
following identity:

$r=\frac{c_{1}}{c_{2}}=\frac{\mathbb{E}_{p_{2}}[q_{1}(\boldsymbol{\theta})\alpha(\boldsymbol{\theta})]}{\mathbb{E}_{p_{1}}[q_{2}(\boldsymbol{\theta})\alpha(\boldsymbol{\theta})]},$
(2)

where $\alpha$ is the bridge function (discussed below), and
$\mathbb{E}_{p_{i}}$ denotes an expectation with respect to
$\boldsymbol{\theta}\sim p_{i}=q_{i}/c_{i}$, for $i=1,2$. This identity leads
naturally to the bridge sampling estimator:

$\hat{r}=\frac{n_{2}^{-1}\sum_{j=1}^{n_{2}}q_{1}(\boldsymbol{\theta}_{2,j})\alpha(\boldsymbol{\theta}_{2,j})}{n_{1}^{-1}\sum_{j=1}^{n_{1}}q_{2}(\boldsymbol{\theta}_{1,j})\alpha(\boldsymbol{\theta}_{1,j})},$
(3)

where
$\left\\{\boldsymbol{\theta}_{i,1},\dots,\boldsymbol{\theta}_{i,n_{i}}\right\\}$
are (possibly dependent) samples from $p_{i}$, and $n_{i}$ is the number of
samples from $q_{i}$, for $i=1,2$. When the samples are independent, Meng and
Wong (1996) showed that the optimal $\alpha$ is

$\alpha_{r}(\boldsymbol{\theta})\propto[s_{1}q_{1}(\boldsymbol{\theta})+rs_{2}q_{2}(\boldsymbol{\theta})]^{-1},$
(4)

where $s_{i}=n_{i}/(n_{1}+n_{2})\ (i=1,2)$, in the sense that this choice
yields the smallest asymptotic variance for the estimator $\hat{r}$. As
expected, the unknown $r$ appears in the optimal choice of $\alpha$, but this
issue can be easily addressed by introducing an iterative scheme
$\\{r^{(t)},t=1,\ldots,\\}$, where $r^{(t)}$ is given by (3) with
$\alpha=\alpha_{r^{(t-1)}}$. This approach was proposed by Meng and Wong
(1996) who showed that $\lim_{t\rightarrow\infty}r^{(t)}$ maintains the
asymptotic variance of the optimal bridge sampling estimator for $r$ that uses
the true value of $r$ in $\alpha_{r}$.555This turns out to be a consequence of
the fact that the fix-point equation implied by (3)-(4) is equivalent to the
score equation for the maximum likelihood estimator of $r$ from the likelihood
theory (for Monte Carlo integration) as formulated in Kong et al. (2003), and
hence the adaption does not lead to loss of information, at least
asymptotically.

In some cases, there is only one unnormalised density $q=c\pi$, but bridge
sampling is still applicable, and indeed preferred. We can set $q_{1}=q$ in
(3), and then have the freedom to choose the “pairing” density $q_{2}=p_{2}$
to be any normalised density (i.e., $c_{2}=1$), such as a Gaussian
distribution. When we are interested in estimating (ratios of) many
normalising constants, we can choose to apply bridge sampling to them
simultaneously, or separately by judiciously choosing a pairing density for
each (or some) of them (Meng and Schilling, 2002). The choice of $p_{2}$ is
important because the asymptotic variance of $\hat{c}$ decreases as the
separation between $p_{1}=\pi$ and $p_{2}$ decreases. Their separation is
measured by the sample-size adjusted harmonic divergence

$\displaystyle
H_{A}(p_{1},p_{2})=1-\int\big{[}\eta_{1}p_{1}^{-1}(\boldsymbol{\theta})+\eta_{2}p_{2}^{-1}(\boldsymbol{\theta})\big{]}^{-1}\mu(\mathrm{d}\boldsymbol{\theta}),$
(5)

where $\eta_{i}\propto s_{i}^{-1}$. Obviously, $p_{2}(\boldsymbol{\theta})$
should also be chosen to minimise computation in terms of both function
evaluation and sampling.

Given the above considerations, if $\pi=q/c$ is multi-modal, it is natural to
choose $p_{2}$ to be a Gaussian mixture distribution approximating $\pi$,
i.e.,

$\phi_{\text{mix}}(\boldsymbol{\theta})=\sum_{k=1}^{K}\phi^{(k)}(\boldsymbol{\theta})=\sum_{k=1}^{K}w_{k}|\mathcal{S}_{k}^{-1}|\phi(\mathcal{S}_{k}^{-1}(\boldsymbol{\theta}-\boldsymbol{\mu}_{k})),$
(6)

where $w_{k}$, $\boldsymbol{\mu}_{k}$, and $\mathcal{S}_{k}$ denote a weight,
mean vector, and scale matrix, respectively, and $\phi$ is the density of
$\mathcal{N}_{d}(\mathbf{0},\mathbf{I})$. In (6) and throughout, we use
$\phi^{(k)}$ to denote the $k$-th component of $\phi_{\text{mix}}$ including
its mixture weight $w_{k}$. Standard bridge sampling estimation would proceed
by applying (3) with densities $q_{1}=q$ (unnormalised) and
$q_{2}=\phi_{\text{mix}}$. The top left panel of Figure 1 provides an
illustrative example in which $\phi_{\text{mix}}$ (dashed curve) is a three
component Gaussian mixture and roughly approximates $\pi$ (solid curve). In
general, the quality of the approximation affects the accuracy of the
estimator $\hat{r}$ in (3).

Figure 1: Top left: the Gaussian mixture $\phi_{\text{mix}}$ (dashed line) is
a rough approximation to the target density $\pi$ (solid line). Top right: the
components $\pi^{(1)}$, $\pi^{(2)}$, and $\pi^{(3)}$ of $\pi$ sliced by
$\phi_{\text{mix}}$ are shown as red, green and blue regions, respectively.
Bottom left: the transformed density $\tilde{\pi}$ obtained by applying the
stochastic Warp-U transformation $\mathcal{F}_{\psi}$. Bottom right: applying
$\mathcal{G}_{\psi^{\prime},\psi}=\mathcal{H}_{\psi^{\prime}}\circ\mathcal{F}_{\psi}$
(i.e., the Warp-U transformation and its inverse) exactly recovers the
original $\pi$ (black solid line); the shaded colours indicate the proportion
of density transported from the original components of $\pi$. Algorithm 1
Warp-U estimation of $c$ (Wang et al., 2022)

Inputs: samples
$\boldsymbol{\theta}_{1,1},\dots,\boldsymbol{\theta}_{1,n_{1}}\sim\pi$ and
$\boldsymbol{\theta}_{2,1},\dots,\boldsymbol{\theta}_{2,n_{2}}\sim\phi$, and
parameters $w_{k},\mathcal{S}_{k},\boldsymbol{\mu}_{k}$ for $k=1,\cdots,K$,
see (6).

1: Set
$\boldsymbol{\theta}_{1,j}^{*}=\mathcal{F}_{\psi_{j}}({\boldsymbol{\theta}_{1,j}})=\mathcal{S}_{\psi_{j}}^{-1}({\boldsymbol{\theta}_{1,j}}-\boldsymbol{\mu}_{\psi_{j}})$
with $\psi_{j}$’s independently drawn from (7) for $j=1,\dots,n_{1}$.

2: Apply (3) using the densities $q_{1}=\tilde{\pi}$ in (8) and $q_{2}=\phi$,
and the samples
$\boldsymbol{\theta}_{1,j}^{*},\dots,\boldsymbol{\theta}_{1,n_{1}}^{*}\sim\tilde{\pi}$
(from Step 1) and
$\boldsymbol{\theta}_{2,1},\dots,\boldsymbol{\theta}_{2,n_{2}}\sim\phi$ (from
the inputs).

Wang et al. (2022) proposed an improved approach called Warp-U bridge
estimation based on the idea of warp bridge sampling estimation (Meng and
Schilling, 2002). By the properties of $f$-divergences (Ali and Silvey, 1966),
for which the Harmonic diverge of (5) is a special case, any transformation
$\mathcal{F}$ will not increase the divergence between $\pi$ and
$\phi_{\text{mix}}$, i.e., $H_{A}(\pi,\phi_{\text{mix}})\geq
H_{A}(\mathcal{F}(\pi),\mathcal{F}(\phi_{\text{mix}}))$. Since the asymptotic
variance of $\hat{r}$ decreases with the divergence (5), a transformation of
the densities will generally decrease the asymptotic variance, as long as the
transformation does not alter the underlying normalising constant. In many
cases, a good choice of $\mathcal{F}$ can provide substantial efficiency
gains.

Given all the desirable properties discussed above, a natural choice of
transformation would be the distributional normalisation:
$\boldsymbol{\theta}^{*}=\mathcal{F}_{\psi}(\boldsymbol{\theta}):=\mathcal{S}_{\psi}^{-1}(\boldsymbol{\theta}-\boldsymbol{\mu}_{\psi})$,
where the index $\psi\in\\{1,\cdots,K\\}$ is drawn from

$\varpi(\psi|\boldsymbol{\theta})=\frac{\phi^{(\psi)}(\boldsymbol{\theta})}{\phi_{\text{mix}}(\boldsymbol{\theta})}.$
(7)

This conditional distribution is the weight of the $\psi$-th component of
$\phi_{\text{mix}}$ at $\boldsymbol{\theta}$. When
$\boldsymbol{\theta}\sim\phi_{\text{mix}}$ and $\psi$ is sampled from (7),
$\mathcal{F}_{\psi}(\boldsymbol{\theta})\sim\mathcal{N}_{d}(\mathbf{0},\mathbf{I})$.
That is, in addition to the usual normalisation that renders mean zero and
variance one, the stochastic map $\mathcal{S}_{\psi}$ normalises
$\phi_{\text{mix}}$ in distribution, since it becomes
$\mathcal{N}_{d}(\mathbf{0},\mathbf{I})$. This is the Warp-U transformation
given in Wang et al. (2022).

Since $\phi_{\text{mix}}$ approximates $\pi$, it is intuitive that applying
the stochastic Warp-U transformation to $\pi$ (with
$\boldsymbol{\theta}\sim\pi=cq$ and
$\psi\sim\varpi(\psi|\boldsymbol{\theta})$) also gives an approximately uni-
modal density. The transformed version of $q$ is given by

$\tilde{q}(\boldsymbol{\theta}^{*})=\sum_{k=1}^{K}\phi(\boldsymbol{\theta}^{*})\frac{q(\mathcal{S}_{k}\boldsymbol{\theta}^{*}+\boldsymbol{\mu}_{k})}{\phi_{\text{mix}}(\mathcal{S}_{k}\boldsymbol{\theta}^{*}+\boldsymbol{\mu}_{k})}w_{k}.$
(8)

Wang et al. (2022) show that $\tilde{q}$ maintains the same normalising
constant as $q$, and hence we can apply bridge sampling to $\tilde{q}$ for
estimating $c$, as detailed in Algorithm 1.

To visualise (8), the coloured regions in the top right panel of Figure 1 show
the components
$\pi^{(k)}(\boldsymbol{\theta})=\varpi(k|\boldsymbol{\theta})\pi(\boldsymbol{\theta})$
(as functions of $\boldsymbol{\theta}$) of $\pi$ (black solid curve) induced
by $\phi_{\text{mix}}$ (the dashed curve in the top left panel). The density
of $\pi^{(k)}$ gets transported towards the origin via the mapping
$\mathcal{F}_{\psi}$ with $\psi=k$. The bottom left panel of Figure 1 is the
resulting Warp-U transformed density $\tilde{\pi}$. The shaded regions
represent the proportion of density contributed from the original components
of $\pi$ with the corresponding colours.

###  2.2 A Distribution-Preserving Mass-swapping Transport

Recall that, when $\boldsymbol{\theta}\sim\pi$, the Warp-U transformed random
variable ${\boldsymbol{\theta}}^{*}=\mathcal{F}_{\psi}(\boldsymbol{\theta})$
follows $\tilde{\pi}=\tilde{q}/c$, where $\tilde{q}$ is given by (8). The
essence of our Warp-U sampler is to additionally apply the inverse Warp-U
transformation to obtain a draw $\theta^{\prime}$ that follows the original
target density $\pi$. That is, we set
$\boldsymbol{\theta}^{\prime}=\mathcal{H}_{\psi^{\prime}}(\boldsymbol{\theta}^{*}):=\mathcal{F}_{\psi^{\prime}}^{-1}(\boldsymbol{\theta}^{*})=\mathcal{S}_{\psi^{\prime}}\boldsymbol{\theta}^{*}+\boldsymbol{\mu}_{\psi^{\prime}}$,
where the random index $\psi^{\prime}$ is generated to be generally different
from $\psi$ but preserves $\boldsymbol{\theta}^{\prime}\sim\pi$. We can
achieve this easily and in general once we notice that as long as
$\mathcal{F}_{\psi}$ is bijective given $\psi$, we have

$\boldsymbol{\theta}=\mathcal{F}^{-1}_{\psi}(\boldsymbol{\theta}^{*}),\quad{\rm
and}\quad\boldsymbol{\theta}^{\prime}=\mathcal{F}^{-1}_{\psi^{\prime}}(\boldsymbol{\theta}^{*})$
(9)

by construction. This parallel expression immediately establishes the
following general result, which applies to any bijective $\mathcal{F}_{\psi}$
(e.g., it does not need to be a Warp-U transformation) and any stochastic
index $\psi$ (e.g., it does not need to be discrete).

######  Theorem 1.

Let $\boldsymbol{\theta}^{*}=\mathcal{F}_{\psi}(\boldsymbol{\theta})$ be a
stochastic map from $\Theta$, the support of a random variable
$\boldsymbol{\theta}$, to itself, where the index $\psi$ is a random variable
which has a well-defined joint distribution with $\boldsymbol{\theta}$.
Suppose $\mathcal{F}_{\psi}(\cdot)$ is bijective for any given value of
$\psi$. Let $\psi^{\prime}$ be a random draw from
$v(\psi|\boldsymbol{\theta}^{*})$, the conditional distribution of $\psi$
given $\boldsymbol{\theta}^{*}$, as implied by the joint distribution of
$(\psi,\boldsymbol{\theta})$. Then
$\boldsymbol{\theta}^{\prime}=\mathcal{F}^{-1}_{\psi^{\prime}}\left(\mathcal{F}_{\psi}(\boldsymbol{\theta})\right)$
and $\boldsymbol{\theta}$ are identically distributed. Furthermore, if
$\psi^{\prime}$ and $\psi$ are independent conditional on
$\boldsymbol{\theta}^{*}$, then $\boldsymbol{\theta}^{\prime}$ and
$\boldsymbol{\theta}$ are i.i.d. given $\boldsymbol{\theta}^{*}$.

Consequently, by drawing $\psi^{\prime}$ independently from

$\nu(\psi|\boldsymbol{\theta}^{*})\propto\varpi(\psi|\boldsymbol{\theta}=\mathcal{H}_{\psi}(\boldsymbol{\theta}^{*}))q(\mathcal{H}_{\psi}(\boldsymbol{\theta}^{*}))\times|\mathcal{H^{{}^{\prime}}}_{\psi}(\boldsymbol{\theta}^{*})|,$
(10)

where $|\mathcal{H^{{}^{\prime}}}_{\psi^{\prime}}(\boldsymbol{\theta}^{*})|$
is the Jacobian of $\mathcal{F}_{\psi}^{-1}$, we not only ensure
$\mathcal{G}_{\psi^{\prime},\psi}=\mathcal{H}_{\psi^{\prime}}\circ\mathcal{F}_{\psi}$
is distribution preserving, but also that its output
$\boldsymbol{\theta}^{\prime}$ is conditionally independent of
$\boldsymbol{\theta}$ given $\boldsymbol{\theta}^{*}$. This conditional
independence, which holds regardless of the distribution of
$\boldsymbol{\theta}^{*}$, is critical for preventing our algorithm from being
trapped by any particular mode, statistically speaking. Note that neither the
distribution preserving property nor the conditional independence requires
knowledge of how well $\phi_{\rm mix}$ approximates $\pi$.

The bottom right panel of Figure 1 illustrates the result of the two-step
stochastic transformation
$\mathcal{G}_{\psi^{\prime},\psi}=\mathcal{H}_{\psi^{\prime}}\circ\mathcal{F}_{\psi}$
applied to $\pi$, using the $\phi_{\text{mix}}$ in the top left panel. The
transformation exactly recovers the original density $\pi$ (the solid black
curve). More importantly, the transformation swaps probability masses among
the original components of $\pi$, which are coloured in the top right panel.
In the bottom right panel of Figure 1, the shaded regions indicate the
proportion of density contributed from each original component of $\pi$ with
the corresponding colour. At most $\boldsymbol{\theta}$ locations, there is a
considerable amount of density transported from each of the original
components of $\pi$. A significant consequence is that by repeatedly applying
the random transformation $\mathcal{G}$ to $\pi$ the mass between its
components gets quickly mixed.

Here the transported masses are calculated in the following way. Let

$\pi_{\psi^{\prime},\psi}(\boldsymbol{\theta})=\pi\big{(}\mathcal{G}_{\psi^{\prime},\psi}^{-1}(\boldsymbol{\theta})\big{)}\times|(\mathcal{G}_{\psi^{\prime},\psi}^{-1})^{\prime}(\boldsymbol{\theta})|\times
p(\psi^{\prime},\psi|\mathcal{G}_{\psi^{\prime},\psi}^{-1}(\boldsymbol{\theta})),$
(11)

where
$|(\mathcal{G}_{\psi^{\prime},\psi}^{-1})^{\prime}(\boldsymbol{\theta})|$ is
the Jacobian of $\mathcal{G}_{\psi^{\prime},\psi}^{-1}$, and
$p(\psi^{\prime},\psi|\xi)=\varpi\big{(}\psi|\xi)\times\nu(\psi^{\prime}|\mathcal{F}_{\psi}(\xi))$
is the probability of choosing the transformation
$\mathcal{G}_{\psi^{\prime},\psi}$ at $\boldsymbol{\theta}=\xi$ (with
$\boldsymbol{\theta}$ being generic notation), that is, the probability of
selecting $\psi$ and transitioning to $\psi^{\prime}$. Then by applying
Theorem 1 with discrete index $\psi$, we have

$\pi(\boldsymbol{\theta})=\sum_{\psi=1}^{K}\sum_{\psi^{\prime}=1}^{K}\pi_{\psi^{\prime},\psi}(\boldsymbol{\theta}),$
(12)

because Theorem 1 tells us that averaging over all the possible transitions
must recover $\pi$, since $\boldsymbol{\theta}$ and
$\boldsymbol{\theta}^{\prime}$ have the same distribution. Therefore, the
function

$f^{(\psi)}(\boldsymbol{\theta})=\sum_{\psi^{\prime}=1}^{K}\pi_{\psi^{\prime},\psi}(\boldsymbol{\theta})$
(13)

can be interpreted as the amount of density transported or
redistributed666Note that this redistribution includes the term
$\pi_{\psi,\psi}(\theta)$. from the $\psi$-th original component of $\pi$. The
densities $f^{(\psi)}(\boldsymbol{\theta}),\psi=1,2,3$ correspond to the three
shaded areas with different colours in the bottom right panel of Figure 1.

###  2.3 A Warp-U Sampler

Algorithm 2 Basic Warp-U MCMC sampler

Input: parameters $K$ and $w_{k},S_{k},\boldsymbol{\mu}_{k}$, for
$k=1,\dots,K$, see (6), proposal variance $\sigma^{2}$, initial value
$\boldsymbol{\theta}_{0}$, and number of samples to be collected $T$.

1: for $t=1,2,\dots,T$ do

2: (i) Generate ${\boldsymbol{\theta}}^{\text{\tiny{MH}}}$ using a Metropolis-
Hasting step with proposal
$\mathcal{N}(\boldsymbol{\theta}_{t-1},\sigma^{2}\mathbf{I})$.

3: (ii) Sample $\psi$ from
$\varpi(\psi|{\boldsymbol{\theta}}^{\text{\tiny{MH}}})=\phi^{(\psi)}({\boldsymbol{\theta}}^{\text{\tiny{MH}}})/\phi_{\text{mix}}({\boldsymbol{\theta}}^{\text{\tiny{MH}}}).$

4: (iii) Set
$\boldsymbol{\theta}^{*}=\mathcal{F}_{\psi}({\boldsymbol{\theta}}^{\text{\tiny{MH}}})$,
where
$\mathcal{F}_{\psi}({\boldsymbol{\theta}}^{\text{\tiny{MH}}})=\mathcal{S}_{\psi}^{-1}({\boldsymbol{\theta}}^{\text{\tiny{MH}}}-\boldsymbol{\mu}_{\psi})$.

5: (iv) Sample $\psi^{\prime}$ from
$\nu(\psi^{\prime}|\boldsymbol{\theta}^{*})\propto\varpi(\psi^{\prime}|\mathcal{H}_{\psi^{\prime}}(\boldsymbol{\theta}^{*}))q(\mathcal{H}_{\psi^{\prime}}(\boldsymbol{\theta}^{*}))|\mathcal{H^{{}^{\prime}}}_{\psi^{\prime}}(\boldsymbol{\theta}^{*})|$.

6: (v) Set
$\boldsymbol{\theta}_{t}=\mathcal{H}_{\psi^{\prime}}(\boldsymbol{\theta}^{*})=\mathcal{S}_{\psi^{\prime}}\boldsymbol{\theta}^{*}+\boldsymbol{\mu}_{\psi^{\prime}}$.

7: end for

The redistribution of masses via Warp-U transformations provides a candidate
MCMC sampler. Given one initial sample $\boldsymbol{\theta}_{0}$, we can
repeatedly apply the random transformation
$\boldsymbol{\theta}_{t}=\mathcal{G}_{\psi^{\prime},\psi}(\boldsymbol{\theta}_{t-1})$
to generate a sequence of $\boldsymbol{\theta}$’s which switches among the
target components. However, this switching is not quite sufficient to
construct a valid MCMC sampler, because the resulting Markov chain is not
guaranteed to be irreducible in general. For illustration, consider an example
where the auxiliary distribution (6) is set as

$\phi_{\text{mix}}(\boldsymbol{\theta})=\sum_{k=1}^{K}\phi(\boldsymbol{\theta}-\boldsymbol{\mu}_{k})/K.$
(14)

In this case, each
$\mathcal{G}_{\psi^{\prime},\psi}(\boldsymbol{\theta})=\boldsymbol{\theta}-\boldsymbol{\mu}_{\psi}+\boldsymbol{\mu}_{\psi^{\prime}}$
is a shift transformation. Given the initial $\boldsymbol{\theta}_{0}$, the
subsequent sequence $\boldsymbol{\theta}_{t}$ can only visit the set

$\Theta=\big{\\{}\boldsymbol{\theta}:\;\boldsymbol{\theta}=\boldsymbol{\theta}_{0}+j_{1}\boldsymbol{\mu}_{1}+\cdots+j_{K}\boldsymbol{\mu}_{K},\;j_{1},\cdots,j_{K}\in\mathbb{Z}\big{\\}}.$

where $\mathbb{Z}$ is the set of all integers. The set $\Theta$ is a countable
grid in $\mathbb{R}^{d}$, and the minimal spacing between the grid points is
$\min_{k\neq
k^{\prime}}\|\boldsymbol{\mu}_{k}-\boldsymbol{\mu}_{k^{\prime}}\|_{2}$. When
the target $\pi$ has a (continuous) compact support set $\Omega$, the chain
can only move among a finite number of candidates inside $\Omega\cap\Theta$,
and hence it cannot converge to the target density for any starting value
$\boldsymbol{\theta}_{0}$.

This issue is fully resolved by introducing a Metropolis-Hasting (MH) step
between consecutive $\mathcal{G}_{\psi^{\prime},\psi}$ transformations.
Algorithm 2 specifies the basic version of our Warp-U MCMC sampler. The
algorithm proceeds by executing a random walk MH step at the beginning of each
iteration followed by the stochastic transformations $\mathcal{F}_{\psi}$ and
$\mathcal{F}_{\psi^{\prime}}^{-1}$ with random indices $\psi,\psi^{\prime}$.
The ergodicity of this MCMC sampler will be verified in Section 3.1. It is
worth noting that the proposed sampler is still valid if the MH step in Line 2
of Algorithm 2 is replaced by other samplers, e.g., the Metropolis adjusted
Langevin or Hamiltonian Monte Carlo samplers.

Figure 2: One marginal of the true target density (solid red line) and the
estimated density (centre dash-dot blue line) obtained by running Algorithm 2.
The left and right panels show the case where $\phi_{\text{mix}}$ is set to be
the target density and a two Gaussian mixture approximation, respectively. The
shaded blue regions indicate pointwise $95\%$ confidence intervals. Algorithm
3 Adaptive Warp-U MCMC sampler

Input: proposal variance $\sigma^{2}$, initial value
$\boldsymbol{\theta}_{0}$, number of within stage samples $T$.

1: Obtain $T$ samples $\Theta_{0}$ from the uniform density with support
$\Theta$ (or another initial density with support $\Theta$).

2: Compute $\phi^{(0)}_{\text{mix}}$ by fitting $\phi_{\text{mix}}$ using the
samples $\Theta_{0}$ from Step 1.

3: for $s=1,2,\dots,M$ do

4: (i) Apply Algorithm 2 with the following inputs: the parameters of
$\phi^{(s-1)}_{\text{mix}}$, $\sigma^{2}$, $\boldsymbol{\theta}_{0}$, and $T$.
Thus obtain new samples $\Theta_{\text{new}}$. Then set
$\Theta_{s}=\Theta_{s-1}\cup\Theta_{\text{new}}$.

5: (ii) With probability $p_{s}$, obtain $\phi^{(s)}_{\text{mix}}$ by re-
fitting $\phi_{\text{mix}}$ using the samples $\Theta_{s}$ from Step 4,
otherwise set $\phi^{(s)}_{\text{mix}}=\phi^{(s-1)}_{\text{mix}}$.

6: end for

Figure 3: Log Wasserstein distance between the true target density (see
(LABEL:target_goodphimix) in Section LABEL:appendix:setting of the online
Supplementary Material) and a kernel density estimate computed from samples
generated via four versions of Algorithm 3, against the stage number of
Algorithm 3. The four algorithm versions are described in the main text and
Table 1.

The performance of Algorithm 2 depends on the quality of the approximation
$\phi_{\text{mix}}$. Figure 2 shows an example where samples obtained via
Algorithm 2 are employed to recover one marginal of a four-dimensional target
density $\pi$. The true density $\pi$ is a mixture of Gaussian distributions
specified by (LABEL:target_goodphimix) in Section LABEL:appendix:setting of
the online Supplementary Material. The left panel of Figure 2 shows the true
density (solid red line) and estimated density (centre dash-dot blue line) in
the case where $\phi_{\text{mix}}$ is the true target $\pi$. The right panel
shows the case where $\phi_{\text{mix}}$ is a two-mode Gaussian mixture
distribution given by (LABEL:target_twomode) in Section
LABEL:appendix:setting. The blue regions show pointwise 95% confidence regions
obtained by running Algorithm 2 fifty times with the same initial value and
$T=3,000$ and then applying the kernel density estimation function density in
R to each set of the $3,000$ samples obtained. In both panels, the confidence
regions approximately cover the true density, but as expected, the uncertainty
region is much narrower in the left panel where the underlying
$\phi_{\text{mix}}$ is a closer (exact) approximation to $\pi$. Note that we
deliberately did not impose any burn-in period in order to measure the quality
of the algorithms in their entirety.

In order to make our Warp-U MCMC sampler less dependent on the quality of
$\phi_{\text{mix}}$, we introduce an adaptive version in Algorithm 3, which
performs a periodic update of $\phi_{\text{mix}}$. In particular, Algorithm 3
adopts the strategy of using increasingly rare updates (Chimisov et al., 2018;
Pompe et al., 2020). After each round of $T$ iterations of Algorithm 2,
Algorithm 3 updates $\phi_{\text{mix}}$ from the accumulated samples
$\Theta_{s}$ with probability $p_{s}$, where $\lim_{s\to\infty}p_{s}=0$. (We
adopt $p_{s}=e^{1-s^{1/8}}$ in the numerical sections, but it is important to
note that any $p_{s}$ satisfying the condition $\lim_{s\to\infty}p_{s}=0$ can
be employed.) This strategy guarantees the ergodicity of the algorithm while
at the same time allowing the component number $K$ in $\phi_{\text{mix}}$ to
gradually increase when needed in order to approximate highly multi-modal
$\pi$.

In our numerical studies, we implement the EM algorithm to update
$\phi_{\text{mix}}$ based on the accumulated samples $\Theta_{s}$. The
complexity of updating the Gaussian mixture approximation $\phi_{\text{mix}}$
is $\mathcal{O}(t_{\max}KTs)$, where $t_{\max}$ is the maximum number of
iterations of the EM algorithm, and $s$ is the current stage number of
Algorithm 3. In practice, to save computational resources, instead of using
all $sT$ samples of $\Theta_{s}$ to update $\phi_{\text{mix}}$ in Step 5 of
Algorithm 3, we can use $T$ samples randomly selected from $\Theta_{s}$.
Alternatively, the update of $\phi_{\text{mix}}$ can be eventually stopped,
e.g., once $s>10$. Figure 3 shows the log Wasserstein distance between the
true target density and a kernel density estimate computed from samples
generated via the following four versions of Algorithm 3: (i) Algorithm 3 as
written above (dash-dot blue line), (ii) Step 5 is skipped for $s>10$ (dotted
orange line), (iii) the same as (i) except Step 5 uses $T$ samples randomly
selected from $\Theta_{s}$ to update $\phi_{\text{mix}}$, rather than all the
available samples (dashed purple line), and (iv) the same as (iii) except that
Step 5 is skipped for $s>10$ (solid green line). See Table 1 for a summary of
these four versions of Algorithm 3. The log Wasserstein distances were
computed using the R package ‘transport’. Figure 3 shows that all four
procedures perform similarly as the stage number $s$ increases.

Table 1: Four versions of Algorithm 3.

Step 5 not skipped Step 5 skipped for $s>10$ Step 5 using all samples (i) (ii)
Step 5 using randomly selected samples (iii) (iv)

###  2.4 An Unbiased Warp-U Estimator

Algorithm 4 Warp-U Sampler with Unbiased Estimator

1: Sample $\boldsymbol{\theta}_{1,0},\boldsymbol{\theta}_{2,0}$ from some
initial distribution $\pi_{0}$.

2: Generate ${\boldsymbol{\theta}}^{\text{\tiny{MH}}}_{1}$ using a Metropolis-
Hasting step with proposal
$\mathcal{N}(\boldsymbol{\theta}_{1,0},\sigma^{2}\mathbf{I})$.

3: Set
$\boldsymbol{\theta}^{*}_{1}=\mathcal{F}_{\psi_{1}}({\boldsymbol{\theta}}^{\text{\tiny{MH}}}_{1})$,
where
$\psi_{1}\sim\varpi(\psi_{1}|{\boldsymbol{\theta}}^{\text{\tiny{MH}}}_{1})$.

4: Set
$\boldsymbol{\theta}_{1,1}=\mathcal{H}_{\psi^{\prime}}(\boldsymbol{\theta}^{*}_{1})$
where
$\psi_{1}^{\prime}\sim\nu(\psi^{\prime}_{1}|\boldsymbol{\theta}^{*}_{1})$.

5: for $t=2,3,\dots,T$ do

6: Generate
${\boldsymbol{\theta}}^{\text{\tiny{MH}}}_{1},{\boldsymbol{\theta}}^{\text{\tiny{MH}}}_{2}$
jointly using a Metropolis-Hasting step with maximal coupling (or reflection
coupling). The marginals of the joint proposal are
$\mathcal{N}(\boldsymbol{\theta}_{1,t-1},\sigma^{2}\mathbf{I})$ and
$\mathcal{N}(\boldsymbol{\theta}_{2,t-2},\sigma^{2}\mathbf{I})$.

7: Set
$\boldsymbol{\theta}^{*}_{1}=\mathcal{F}_{\psi_{1}}({\boldsymbol{\theta}}^{\text{\tiny{MH}}}_{1})$
and
$\boldsymbol{\theta}^{*}_{2}=\mathcal{F}_{\psi_{2}}({\boldsymbol{\theta}}^{\text{\tiny{MH}}}_{2})$,
where $(\psi_{1},\psi_{2})$ are sampled with optimal transport coupling with
marginals $\varpi(\psi_{1}|{\boldsymbol{\theta}}^{\text{\tiny{MH}}}_{1})$ and
$\varpi(\psi_{2}|{\boldsymbol{\theta}}^{\text{\tiny{MH}}}_{2})$.

8: Set
$\boldsymbol{\theta}_{1,t}=\mathcal{H}_{\psi^{\prime}_{1}}(\boldsymbol{\theta}^{*}_{1})$
and
$\boldsymbol{\theta}_{2,t-1}=\mathcal{H}_{\psi^{\prime}_{2}}(\boldsymbol{\theta}^{*}_{2})$
where $(\psi^{\prime}_{1},\psi^{\prime}_{2})$ are sampled with optimal
transport coupling with marginals
$\psi_{1}^{\prime}\sim\nu(\psi^{\prime}_{1}|\boldsymbol{\theta}^{*}_{1})$ and
$\psi_{2}^{\prime}\sim\nu(\psi^{\prime}_{2}|\boldsymbol{\theta}^{*}_{2})$.

9: end for

10: Compute
$H_{l:m}(\boldsymbol{\theta}_{1,\cdot},\boldsymbol{\theta}_{2,\cdot})$ in
(17).

Given MCMC samples from $\pi$, we are often interested in obtaining an
estimate of $\mathbb{E}_{\pi}h(\boldsymbol{\theta})$ for some function
$h(\cdot)$, where the expectation is taken with respect to
$\boldsymbol{\theta}\in\pi(\boldsymbol{\theta})$. Jacob et al. (2020)
described a general strategy for unbiased estimation of
$\mathbb{E}_{\pi}h(\boldsymbol{\theta})$ that involves running two coupled
Markov chains $(\boldsymbol{\theta}_{1,t},\boldsymbol{\theta}_{2,t-1})$.
Marginally, each chain ($\\{\boldsymbol{\theta}_{1,t}\\}_{t=0}^{\infty}$ or
$\\{\boldsymbol{\theta}_{2,t}\\}_{t=0}^{\infty}$) is updated by the same
transition kernel which admits $\pi$ as its stationary distribution. In
addition, the coupled chains should meet the following three requirements:

>   1. R1.
>
> As
> $t\to\infty,\mathbb{E}\left[h\left(\boldsymbol{\theta}_{1,t}\right)\right]\rightarrow\mathbb{E}_{\pi}[h(\boldsymbol{\theta})]$.
> There exists an $\eta>0$ and $D<\infty$ such that
> $\mathbb{E}\left[\left|h\left(\boldsymbol{\theta}_{1,t}\right)\right|^{2+\eta}\right]\leq
> D$ for all $t\geq 0$.
>
>   2. R2.
>
> The meeting time of the two chains $\tau:=\inf\left\\{t\geq
> 1:\boldsymbol{\theta}_{1,t}=\boldsymbol{\theta}_{2,t-1}\right\\}$ satisfies
> $\mathbb{P}(\tau>t)\leq C\delta^{t}$ for all $t\geq 0$, for some constants
> $C<\infty$ and $\delta\in(0,1)$.
>
>   3. R3.
>
> The two chains stay together after meeting, i.e.
> $\boldsymbol{\theta}_{1,t}=\boldsymbol{\theta}_{2,t-1}$ for all $t\geq\tau$.
>
>

Given two coupled chains
$(\boldsymbol{\theta}_{1,t},\boldsymbol{\theta}_{2,t-1})$ satisfying R1–R3,
the unbiased estimation strategy of Jacob et al. (2020) is based on the
following relation for a fixed integer $l$:

$\displaystyle\mathbb{E}_{\pi}[h(\boldsymbol{\theta})]$
$\displaystyle=\mathbb{E}\left[h\left(\boldsymbol{\theta}_{1,l}\right)\right]+\sum_{t=l+1}^{\infty}\big{\\{}\mathbb{E}\left[h\left(\boldsymbol{\theta}_{1,t}\right)\right]-\mathbb{E}\left[h\left(\boldsymbol{\theta}_{1,t-1}\right)\right]\big{\\}}$
$\displaystyle=\mathbb{E}\left[h\left(\boldsymbol{\theta}_{1,l}\right)\right]+\sum_{t=l+1}^{\infty}\big{\\{}\mathbb{E}\left[h\left(\boldsymbol{\theta}_{1,t}\right)\right]-\mathbb{E}\left[h\left(\boldsymbol{\theta}_{2,t-1}\right)\right]\big{\\}}$
$\displaystyle=\mathbb{E}\Big{[}h\left(\boldsymbol{\theta}_{1,l}\right)+\sum_{t=l+1}^{\tau-1}\left(h\left(\boldsymbol{\theta}_{1,t}\right)-h\left(\boldsymbol{\theta}_{2,t-1}\right)\right)\Big{]}.$
(15)

An unbiased estimator of $\mathbb{E}_{\pi}[h(\boldsymbol{\theta})]$ based on
the coupled chains $(\boldsymbol{\theta}_{1,t},\boldsymbol{\theta}_{2,t-1})$
is

$H_{j}(\boldsymbol{\theta}_{1,\cdot}\,,\,\boldsymbol{\theta}_{2,\cdot})=h\left(\boldsymbol{\theta}_{1,j}\right)+\sum_{t=j+1}^{\tau-1}\left[h\left(\boldsymbol{\theta}_{1,t}\right)-h\left(\boldsymbol{\theta}_{2,t-1}\right)\right].$
(16)

In the above, the last summation term serves as a bias correct and it is zero
when $\tau\leq j+1$; and $j$ is an arbitrary fixed value. With a single run of
the coupled chains, we can obtain an average estimator (averaged across $j$):

$H_{l:m}(\boldsymbol{\theta}_{1,\cdot}\,,\,\boldsymbol{\theta}_{2,\cdot})=\frac{1}{m-l+1}\sum_{j=l}^{m}H_{j}(\boldsymbol{\theta}_{1,\cdot}\,,\,\boldsymbol{\theta}_{2,\cdot})$
(17)

as the final estimator for $\mathbb{E}_{\pi}h(\boldsymbol{\theta})$, for some
fixed integers $l$ and $m$ with $(l\leq m)$.

To facilitate the estimation of $\mathbb{E}_{\pi}h(\boldsymbol{\theta})$ when
$\pi$ is multi-modal, we adapt the strategy of Jacob et al. (2020) to the
basic Warp-U sampler in Algorithm 2. This unbiased Warp-U estimation method is
detailed in Algorithm 4. We initialise two chains via an initial distribution
$\boldsymbol{\theta}_{1,0},\boldsymbol{\theta}_{2,0}\sim\pi_{0}$ (Line 1,
Algorithm 4). A one-step update is then conducted for the first chain
$\boldsymbol{\theta}_{1,0}$ (Lines 2–4), which is essentially a single
iteration of Algorithm 2. Subsequently, we update the two chains
$(\boldsymbol{\theta}_{1,t},\boldsymbol{\theta}_{2,t-1})$ jointly (Lines 5–9).
Specifically, a Metropolis-Hasting step with maximal coupling or reflection
coupling (Jacob et al., 2020) is executed (Line 6), followed by the coupled
Warp-U transformation (Line 7) and its inverse (Line 8). The coupled Warp-U
transformation strategy is based on the optimal transport coupling (see
Supplementary Material, Section LABEL:appendix:unbiased:transport). This
strategy preserves the correct marginal distributions while encouraging
proximity of the two chains during their transitions between modes. Note this
strategy minimizes the expected squared distance between the two chains after
transition, as in Papp and Sherlock (2022).

Full algorithmic details are available in the online Supplementary Material
(Section LABEL:appendix:unbiased). The requirements R1–R3 are also verified
for the proposed algorithm in Section LABEL:appendix:unbiased:transport.
Meanwhile, a Rao-Blackwellization technique is developed in Section
LABEL:appendix:unbiased:raoblack to reduce the variance of (16) for estimating
$\mathbb{E}_{\pi}h(\boldsymbol{\theta})$. Specifically, the summands in (16)
are replaced by conditional expectations, where the conditioning is on one or
more of the intermediate samples
${\boldsymbol{\theta}}^{\text{\tiny{MH}}}_{1}$,
${\boldsymbol{\theta}}^{\text{\tiny{MH}}}_{2}$, $\boldsymbol{\theta}^{*}_{1}$
and $\boldsymbol{\theta}^{*}_{2}$ in Algorithm 4. Lastly, note that the
optimal transport couplings for the Warp-U transformation (Line 7) and its
inverse (Line 8) are calculated in two separate steps in Algorithm 4. By
integrating Lines 7–8 into a one-step coupling, we can reduce chain meeting
time $\tau$ but with higher computational cost. This version of the algorithm
is detailed in Algorithm LABEL:alg:unbiasedV2 in Section
LABEL:appendix:unbiased:transportV2.

###  2.5 Stochastic Warp-U Bridge Estimation

Algorithm 5 Stochastic Warp-U Bridge Sampling

Input: parameters $K$ and $w_{k},\mathcal{S}_{k},\boldsymbol{\mu}_{k}$, for
$k=1,\dots,K$, see (6).

1: Sample
$\boldsymbol{\theta}_{1,1},\dots,\boldsymbol{\theta}_{1,n_{1}}\sim\pi$.

2: for $i=1,\dots,n_{1}$ do

3: (i) sample $\psi_{i}$ from the probability
$\varpi(\psi|\boldsymbol{\theta})$ in (7).

4: (ii) Set
$\boldsymbol{\theta}^{*}_{1,i}=\mathcal{F}_{\psi_{i}}({\boldsymbol{\theta}}_{1,i})=\mathcal{S}_{\psi_{i}}^{-1}({\boldsymbol{\theta}_{1,i}}-\boldsymbol{\mu}_{\psi_{i}})$.

5: end for

6: for each component $k$ from $1$ to $K$ do

7: (A) Sample
$\boldsymbol{\theta}_{2,k,1},\dots,\boldsymbol{\theta}_{2,k,n_{2}}\sim
N(0,1)$;

8: (B) Apply all $\boldsymbol{\theta}^{*}_{1,i}$ with $\psi_{i}=k$ and the
$n_{2}$ samples from (A) to compute (20) (for whatever chosen $\alpha$)

9: end for

10: Calculate the final estimator $\hat{c}=\sum_{k=1}^{K}w_{k}\hat{c}_{k}$.

Having obtained samples from $\pi$ using the Warp-U sampler in Section 2.3, we
may want to use them to estimate the normalising constant $c$ of $\pi=q/c$,
which is a special case of the problem of estimating
$\mathbb{E}_{\pi}h(\boldsymbol{\theta})$. We propose a new estimation strategy
which adapts and improves the Warp-U bridge estimator proposed by Wang et al.
(2022). The main limitation of the Warp-U bridge sampling estimator is that it
has high computational cost if $q$ is expensive to evaluate, which is an
important case because many estimation approaches can work well if $q$ is
inexpensive to evaluate (indeed, exact integration can be performed in the
limit where evaluating $q$ is free). The high computational cost of the Warp-U
bridge estimator is due to the fact that each evaluation of $\tilde{q}$ in (8)
requires $K$ evaluations of $q$. Overall Algorithm 1 requires $K(n_{1}+n_{2})$
evaluations of $q$ to obtain an estimate $\hat{c}$. Thus, although the Warp-U
bridge estimator is statistically more efficient than the standard bridge
sampling applied to $\pi$ and $\phi_{\text{mix}}$, it is approximately $K$
times more expensive. Consequently, Wang et al. (2022) found that, given fixed
computational resources, Warp-U bridge sampling is only comparable to standard
bridge sampling (and is sometimes slightly worse).

To obtain a more computationally efficient way of estimating $c$, we first
note that because $\tilde{q}$ in (8) preserves the normalising constant $c$,
it implies that $c=\sum_{k=1}^{K}w_{k}c_{k}$, where $c_{k}$ is the normalising
constant of

$\tilde{q}_{k}(\boldsymbol{\theta}^{*})=\phi(\boldsymbol{\theta}^{*})\frac{q(\mathcal{S}_{k}\boldsymbol{\theta}^{*}+\boldsymbol{\mu}_{k})}{\phi_{\text{mix}}(\mathcal{S}_{k}\boldsymbol{\theta}^{*}+\boldsymbol{\mu}_{k})}.$
(18)

By the construction of
$\boldsymbol{\theta}^{*}=\mathcal{F}_{\psi}(\boldsymbol{\theta})$ and by (8),
$\tilde{q}_{k}(\boldsymbol{\theta}^{*})/c_{k}=p(\boldsymbol{\theta}^{*}|\psi=k)$.
In other words, the normalised counterpart of (8) is a mixture density with
the form

$\tilde{\pi}(\boldsymbol{\theta}^{*})=\frac{\tilde{q}(\boldsymbol{\theta}^{*})}{c}=\sum_{k=1}^{K}\frac{w_{k}c_{k}}{c}\frac{\tilde{q}_{k}(\boldsymbol{\theta}^{*})}{c_{k}}=\sum_{k=1}^{K}\tilde{w}_{k}\tilde{p}_{k}(\boldsymbol{\theta}^{*}),$
(19)

where $\tilde{w}_{k}=w_{k}c_{k}/c$ is the mixture weight for
$\tilde{p}_{k}(\boldsymbol{\theta}^{*})\equiv
p(\boldsymbol{\theta}^{*}|\psi=k)$, for $k=1,\dots,K$. This immediately
suggests that we can apply bridge sampling to estimate each $c_{k}$
separately, using the draws $\mathcal{F}_{\psi}(\boldsymbol{\theta})$ where
$\psi=k$, and pairing them with draws from the common auxiliary density
$\phi$, for all $k=1,\ldots,K$.

Our proposed bridge estimation procedure is detailed in Algorithm 5. The
algorithm computes $c_{k}$ from the $\phi$ samples
$\boldsymbol{\theta}_{2,k,j}$ and the $\pi$ samples
$\boldsymbol{\theta}_{1,j}^{*}$ for which $\psi_{j}=k$, i.e.,

$\hat{c}_{k}=\frac{n_{2}^{-1}\sum_{j=1}^{n_{2}}\tilde{q}_{k}(\boldsymbol{\theta}_{2,k,j})\alpha(\boldsymbol{\theta}_{2,k,j})}{n_{1k}^{-1}\sum_{j\in\mathcal{J}_{k}}\phi(\boldsymbol{\theta}^{*}_{1,j})\alpha(\boldsymbol{\theta}^{*}_{1,j})}$
(20)

where $\mathcal{J}_{k}=\\{j:\;\psi_{j}=k\\}$ and $n_{1k}=|\mathcal{J}_{k}|$ is
the cardinality of the index set $\mathcal{J}_{k}$. To estimate $c$, we use
the weighted sum $\hat{c}_{\text{SWB}}:=\sum_{k=1}^{K}w_{k}\hat{c}_{k}$. We
call our estimator $\hat{c}_{\text{SWB}}$ the stochastic Warp-U bridge (SWB)
estimator because it can be viewed as an estimator of the mean normalising
constant of a random unnormalised density. Algorithm 5 requires a total of
$n_{1}+Kn_{2}$ evaluations of $q$, which is $(K-1)n_{1}$ fewer than that
needed for Warp-U bridge estimation. In Section 3.2, we prove that, under some
conditions, $\hat{c}_{\text{SWB}}$ is more efficient than Warp-U bridge
estimation in terms of asymptotic variance and precision per CPU second. The
stochastic Warp-U bridge estimator can in principle be combined with the
unbiased estimation method described in Section 2.4, which is a topic for
future work and is briefly discussed in Section 6.2.

##  3 Theoretical Justification

###  3.1 Ergodicity Property

In this subsection, we establish that Algorithm 3 produces an ergodic Markov
chain with stationary distribution $\pi$. The technical proofs of the
theoretical results in this subsection can be found in the online
Supplementary Material, Section LABEL:appendix:ergodicity_proof.

We make use of Condition 1 below, which has previously been employed by Jarner
and Hansen (2000) and Atchadé (2006) to ensure the ergodicity of the random
walk Metropolis-Hastings algorithm. Let
$n(\boldsymbol{\theta})=\boldsymbol{\theta}/\|\boldsymbol{\theta}\|_{2}$ be a
unit vector and
$m(\boldsymbol{\theta})=\nabla\pi(\boldsymbol{\theta})/\|\nabla\pi(\boldsymbol{\theta})\|_{2}$
be the normalised gradient of $\pi$.

######  Condition 1.

The density $\pi$ is super-exponential, i.e.

$\limsup_{\|\boldsymbol{\theta}\|_{2}\to\infty}\;\langle
n(\boldsymbol{\theta}),\nabla\log\pi(\boldsymbol{\theta})\rangle=-\infty,$

where $\langle\cdot,\cdot\rangle$ is vector inner product. In addition, $\pi$
has a regular manifold contour, i.e.,
$\limsup_{\|\boldsymbol{\theta}\|_{2}\to\infty}\;\langle
n(\boldsymbol{\theta}),m(\boldsymbol{\theta})\rangle<0$.

A Gaussian density or a mixture of Gaussian densities is super-exponential. In
fact, densities of the form
$\pi(\boldsymbol{\theta})=h(\boldsymbol{\theta})\exp(-g(\boldsymbol{\theta}))$
belong to the super-exponential family for many choices of the polynomials
$h(\boldsymbol{\theta})$ and $g(\boldsymbol{\theta})$ (see Jarner and Hansen,
2000).

Recall that Algorithm 3 contains Algorithm 2 within a for loop. Let
$P_{\gamma}$ be the transition kernel for Lines 2–6 of Algorithm 2, where
$\gamma$ represents the adaptive parameters of $\phi_{\text{mix}}$.
Specifically,
$\gamma=\\{K,(\boldsymbol{\mu}_{k},\mathcal{S}_{k},w_{k})_{k=1}^{K}\\}$
includes the number of mixture components $K$, as well as the weight $w_{k}$,
the mean vector $\boldsymbol{\mu}_{k}$ and the scale matrix $\mathcal{S}_{k}$
of each component. The adaptive nature of Algorithm 3 means that $\gamma$ has
a chance to be updated periodically, but with decreasing probability. We
require each $\gamma$ to belong to a constraint set $\Gamma$, which is
specified in Condition 2 below.

######  Condition 2.

The adaptive parameter
$\gamma=\\{K,(\boldsymbol{\mu}_{k},\mathcal{S}_{k},w_{k})_{k=1}^{K}\\}$ always
lies inside the set $\Gamma$ satisfying

$0<K<K_{\max},\ w_{k}\geq 0,\ \sum_{k=1}^{K}w_{k}=1,\
\Lambda_{\min}\leq|\mathcal{S}_{k}|\leq\Lambda_{\max},\text{ and }\
\|\boldsymbol{\mu}_{k}\|_{2}\leq C_{\mu}.$

In the above, $K_{\max}>1$, $\Lambda_{\max}>\Lambda_{\min}>0$, and $C_{\mu}>0$
are constants which may depend on the underlying data (but not the Monte Carlo
samples).

Condition 2 states that the determinants of the scale matrices
$\mathcal{S}_{k}$ should be bounded from above and below. In numerical
implementation, this requirement can be easily satisfied via monitoring the
update of $\phi_{\text{mix}}$; components can be split or merged when
$|\mathcal{S}_{k}|$ gets too large or too small. Condition 2 also requires
that the maximum number of mixture components is bounded from above. In
applications, it suffices to set $K_{\max},C_{\mu}$ to be some large numbers
such that $\phi_{\text{mix}}$ can reasonably approximate $\pi$.

Based on the above conditions, we now establish the ergodicity of Algorithm 3
through three lemmas. Lemma 1 formalises the argument in Section 2.2 that the
target density $\pi$ is the stationary distribution of Algorithm 2.

######  Lemma 1.

For each $\gamma\in\Gamma$, the transition kernel $P_{\gamma}$ for Lines 2–6
of Algorithm 2 satisfies

$\int\pi(\boldsymbol{\theta})P_{\gamma}(\boldsymbol{\theta}^{\prime}|\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})=\pi(\boldsymbol{\theta}^{\prime}).$
(21)

To further prove that $\pi$ is also the stationary distribution of Algorithm
3, we need to establish that Algorithm 3 is ergodic, despite its adaptive
nature. To achieve this we apply the framework proposed by Roberts and
Rosenthal (2007). We will show that the diminishing adaptation and
simultaneous strongly aperiodic geometric ergodicity (SSAGE) conditions hold.
The diminishing adaptation condition is stated in Lemma 2 below and is
straightforward to verify. Indeed, Step 5 of Algorithm 3 automatically
guarantees the kernel update happens with probability decreasing to zero.

######  Lemma 2.

The sequence of the transition kernels within Algorithm 3 satisfies the
diminishing adaptation condition, i.e.,

$\lim_{s\to\infty}\mathbb{P}\big{(}\sup_{\boldsymbol{\theta}\in\Theta}\|P_{\hat{\gamma}^{s}}(\
\cdot\ |\boldsymbol{\theta})-P_{\hat{\gamma}^{s+1}}(\ \cdot\
|\boldsymbol{\theta})\|_{TV}\geq\delta_{1}\big{)}=0$ (22)

for any $\delta_{1}>0$, where TV indicates the total variation distance, and
$\hat{\gamma}^{s},\hat{\gamma}^{s+1}$ are the estimated parameters of
$\phi_{\text{mix}}$ used at stage $s$ and $s+1$, respectively.

Following Roberts and Rosenthal (2007), the family of transition kernels
$\\{P_{\gamma}\\}_{\gamma\in\Gamma}$ satisfies the simultaneous strongly
aperiodic geometric ergodicity (SSAGE) condition, if the following two
properties are satisfied:

  1. P1.

(Minorization) There exists a probability measure $\omega$, a compact Borel
set $C\subset\Theta$, and a constant $\delta\in(0,1]$, such that $\omega(C)>0$
and

$\inf_{\gamma\in\Gamma}P_{\gamma}(A\,|\boldsymbol{\theta})\geq\delta\cdot\omega(A)\cdot
I_{C}(\boldsymbol{\theta})$ (23)

for any Borel set $A\subset\Theta$ and any $\boldsymbol{\theta}\in\Theta$,
where $I_{C}$ is the indicator function of $C$.

  2. P2.

(Simultaneous Drift) There exists $b<\infty$, $\lambda\in(0,1)$, and a
function $V$ (with $V(\boldsymbol{\theta})\geq 1$ and
$\sup_{\boldsymbol{\theta}\in C}V(\boldsymbol{\theta})<\infty$), such that for
any $\boldsymbol{\theta}\in\Theta$ it holds that

$\sup_{\gamma\in\Gamma}\int
V(\boldsymbol{\theta}^{\prime})P_{\gamma}(\,d\boldsymbol{\theta}^{\prime}\,|\boldsymbol{\theta})\leq\lambda
V(\boldsymbol{\theta})+bI_{C}(\boldsymbol{\theta}).$ (24)

In the literature, the drift function $V$ is usually taken as
$V(\boldsymbol{\theta})=c_{\alpha}\pi^{-\alpha}(\boldsymbol{\theta})$ for some
constant $c_{\alpha}$ and $\alpha\in(0,1)$. The SSAGE condition implies that
the family of transition kernels $\\{P_{\gamma}\\}_{\gamma\in\Gamma}$ is
uniformly geometrically ergodic. Hence the $m$ step transition kernel
$P_{\gamma}^{m}$ of Algorithm 2 has the property that
$\sup_{\gamma\in\Gamma}\|P_{\gamma}^{m}(\boldsymbol{\theta},\cdot)-\pi(\cdot)\|_{TV}\leq
R\rho^{m}V(\boldsymbol{\theta})$, for some $\rho<1$ and $R<\infty$ (e.g., see
equation (2.14) of Atchadé, 2006).

The geometric drift property (P2) requires an additional restriction on
$\Gamma$ which regulates the behaviour of the Warp-U transformation of
$\boldsymbol{\theta}$ when $\|\boldsymbol{\theta}\|_{2}\to\infty$.
Specifically, when $\|\boldsymbol{\theta}\|_{2}\to\infty$, we impose an upper
bound on

$\displaystyle
U(\boldsymbol{\theta}):=\mathbb{E}V\big{(}\mathcal{G}_{\psi^{\prime},\psi}(\boldsymbol{\theta})\big{)}/V(\boldsymbol{\theta})=\sum_{\psi,\psi^{\prime}=1}^{K}\big{[}V(\mathcal{G}_{\psi^{\prime},\psi}(\boldsymbol{\theta}))/V(\boldsymbol{\theta})\big{]}\times
p(\psi^{\prime},\psi|\boldsymbol{\theta}),$ (25)

where the above expectation is taken with respect to the random transformation
$\mathcal{G}_{\psi^{\prime},\psi}$, and
$p(\psi^{\prime},\psi|\boldsymbol{\theta})$ is defined in (11). Let
$Q(\boldsymbol{\theta},\cdot)$ denote the measure of
$\mathcal{N}(\boldsymbol{\theta},\sigma^{2}\mathbf{I})$ which corresponds to
the MH proposal in Step (i) of Algorithm 2, and $p(\cdot|\boldsymbol{\theta})$
denote the density function of this proposal. Meanwhile, let
$\mathcal{R}(\boldsymbol{\theta})$ be set of $\boldsymbol{\theta}^{\prime}$
with acceptance probability less than $1$ when starting from
$\boldsymbol{\theta}$, i.e.,

$\mathcal{R}(\boldsymbol{\theta})=\big{\\{}\boldsymbol{\theta}^{\prime}\in\Theta:\
\pi(\boldsymbol{\theta}^{\prime})p(\boldsymbol{\theta}|\boldsymbol{\theta}^{\prime})<\pi(\boldsymbol{\theta})p(\boldsymbol{\theta}^{\prime}|\boldsymbol{\theta})\big{\\}}.$

Lastly, let
$\overline{Q}:=\limsup_{\|\boldsymbol{\theta}\|_{2}\to\infty}Q(\boldsymbol{\theta},\mathcal{R}(\boldsymbol{\theta}))$.
We use the following technical restriction on $\Gamma$.

######  Condition 3.

The adaptive parameter set $\Gamma$ satisfies

$\limsup_{\|\boldsymbol{\theta}\|_{2}\to\infty\;}\sup_{\gamma\in\Gamma}U(\boldsymbol{\theta})<1/\,\overline{Q}.$
(26)

We emphasise that $U(\boldsymbol{\theta})$ depends on the adaptive parameter
$\gamma$ through the random transformation $\mathcal{G}_{\psi^{\prime},\psi}$,
as seen in (25). Intuitively, Condition 3 asserts that the ratio
$\pi(\boldsymbol{\theta})/\pi(\mathcal{G}_{\psi^{\prime},\psi}(\boldsymbol{\theta}))$
is not expected to be too large when $\boldsymbol{\theta}$ is very far away
from the origin. In other words, when $\boldsymbol{\theta}$ is a substantial
distance from the origin, the stochastic transformation
$\mathcal{G}_{\psi^{\prime},\psi}$ is expected to push $\boldsymbol{\theta}$
towards the origin, where the target density is higher. This condition can be
easily verified in some simple settings. For example, when the parameter space
$\Theta=\mathbb{R}$ is one-dimensional, suppose the mixture components in
$\phi_{\text{mix}}$ have balanced weights, the scaling $\mathcal{S}_{k}$’s are
equal (i.e., $\mathcal{S}_{k}=\mathcal{S}_{k^{\prime}}$ for $k\neq
k^{\prime}$), and the component means are not arbitrarily close to each other
(i.e., $|\mu_{k}-\mu_{k^{\prime}}|\geq c_{\mu}$ for all $k\neq k^{\prime}$ and
for an arbitrarily small constant $c_{\mu}>0$). In this case, we can easily
see when $\theta$ is far away from 0, it is almost surely probable that the
stochastic transformation either keeps $\theta$ at its current position or
pushes it towards $0$. However, in general situations, the technical Condition
3 may not be straightforward to check.

######  Remark.

Condition 3 is only technical and regulates the stochastic transformations
when $\boldsymbol{\theta}$ is extremely far away from the origin. This
condition can be completely removed by slightly modifying Algorithm 2 in the
following way. Given a pre-specified and arbitrarily large value
$C_{\theta}>0$, we executed Steps (ii)–(v) only when
$\|\boldsymbol{\theta}^{\text{MH}}\|_{2}\leq C_{\theta}$; otherwise, we set
$\boldsymbol{\theta}_{t}=\boldsymbol{\theta}^{\text{MH}}$. Correspondingly,
the index probability mass function (10) also needs to be changed to

$\nu(\psi^{\prime}|\boldsymbol{\theta}^{*})\propto\varpi(\psi^{\prime}|\mathcal{H}_{\psi^{\prime}}(\boldsymbol{\theta}^{*}))\pi(\mathcal{H}_{\psi^{\prime}}(\boldsymbol{\theta}^{*}))|\mathcal{H^{{}^{\prime}}}_{\psi^{\prime}}(\boldsymbol{\theta}^{*})|\cdot
I(\|\mathcal{H}_{\psi^{\prime}}(\boldsymbol{\theta}^{*})\|_{2}\leq
C_{\theta}).$

See Section LABEL:appendix:WarpUSamplerV2 of the online Supplementary Material
for more discussions.

######  Lemma 3.

Under Condition 1–3, the family of transition kernels
$\\{P_{\gamma}\\}_{\gamma\in\Gamma}$ for Algorithm 2 satisfies the SSAGE
condition.

Finally, Theorem 2 establishes that Algorithm 3 is ergodic. Specifically, the
total variation distance between $\pi$ and the $s$-step transition kernel of
Algorithm 3 converges to zero as $s\to\infty$.

######  Theorem 2.

Under Conditions 1–3, the adaptive Warp-U MCMC sampling algorithm (Algorithm
3) produces a Markov chain which is ergodic with the target distribution $\pi$
as its stationary distribution.

###  3.2 Computational and Statistical Efficiency of Stochastic Bridge
Estimation

In this subsection, we demonstrate that, under some conditions, our stochastic
bridge estimator has lower asymptotic variance and greater Precision per
Second (PpS) than Warp-U bridge estimation. The technical proofs for the
theoretical results in this subsection can be found in the online
Supplementary Material, Section LABEL:appendix:asymVariance.

We define PpS as $1/(\text{RMSE}\times\text{CPU seconds})$, and assume that
other computational costs are negligible compared with evaluating $q=c\pi$.
Therefore CPU seconds are given by $C\cdot E\cdot g(q)$, where $C$ is a
constant, $g(q)$ is the time taken to evaluate $q$ once, and $E$ is the number
of evaluations of $q$ used by the method under consideration. The values of
$E$ for the different methods we compare are shown in the right three columns
of Table 2, where $n_{1}$ and $n_{2}$ denote the number of the samples from
the target and auxiliary distribution, respectively. The auxiliary
distribution is the standard Gaussian $\phi$ for Warp-U bridge estimation and
stochastic bridge estimation, and $\phi_{\text{mix}}$ for classical bridge
sampling. The number of Warp-U bridge estimation target evaluations is lower
in the case of Warp-U MCMC sampling because some of the necessary evaluations
have already been computed during the sampling stage.

Table 2:  Number of evaluations of the unnormalised target density $q$ for
different sampling and estimation methods needed to obtain $n_{1}$ target
draws and $n_{2}$ auxiliary draws which are then used in the estimation step.
No. Iters represents the number of iterations (accept and reject), $n_{1}$ and
$\tilde{n}_{1}$ denote the number of samples and average number of samples at
each stage, respectively, $M$ is the number of stages, $M_{l}$ is the number
of temperature levels in Parallel Tempering, and $n_{2}$ is the number of
samples from the auxiliary distribution.

Sampling No. Iters Sampling Evals. Estimation Evals. Bridge Warp-U S. Warp-U
Warp-U MCMC $n_{1}M$ $(K+1)n_{1}M$ $n_{2}$ $Kn_{2}$ $Kn_{2}$ Generalised Wang-
Landau $\tilde{n}_{1}M$ $\tilde{n}_{1}M$ $n_{2}$ $(K-1)n_{1}+Kn_{2}$ $Kn_{2}$
Parallel Tempering $2n_{1}M_{l}$ $n_{1}M_{l}$ $n_{2}$ $(K-1)n_{1}+Kn_{2}$
$Kn_{2}$ Unknown Sampling $-$ $-$ $n_{1}+n_{2}$ $Kn_{1}+Kn_{2}$ $n_{1}+Kn_{2}$

Under the i.i.d. sample assumption, Meng and Wong (1996) showed that the
asymptotic variance of the logarithm $\hat{\lambda}=\log(\hat{r})$ of the
classical bridge estimator (3) is

$\left(n_{1}^{-1}+n_{2}^{-1}\right)\left[(1-H_{A}(p_{1},p_{2}))^{-1}-1\right]+o\left((n_{1}+n_{2})^{-1}\right),$
(27)

where $H_{A}(p_{1},p_{2})$ is the sample-size adjusted harmonic divergence
(5). We consider the setting when $p_{2}=\phi$ (with $c_{2}=1$) and $p_{1}$
corresponds to the Warp-U transformed density (19) for comparing the Warp-U
bridge estimator ($\hat{\lambda}_{\text{WB}}=\log(\hat{c}_{\text{WB}})$) and
the stochastic Warp-U bridge estimator
($\hat{\lambda}_{\text{SWB}}=\log(\hat{c}_{\text{SWB}})$). The asymptotic
variance (27) can not be directly applied to compare
$\hat{\lambda}_{\text{WB}}$ and $\hat{\lambda}_{\text{SWB}}$, because $H_{A}$
involves sample sizes and the stochastic Warp-U bridge estimator exploits
random subsamples to reduce the computational cost. Theorem 3 below provides
an alternative comparison of the asymptotic variance based on the Pearson’s
$\chi^{2}$ divergence,
$\chi_{P}^{2}(p_{1},p_{2})=\int\big{(}p_{2}(\boldsymbol{\theta})/p_{1}(\boldsymbol{\theta})-1\big{)}^{2}\times
p_{1}(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})$, which
approximates the sample-size adjusted harmonic divergence (5) as
$n_{1},n_{2}\to\infty$, see Nielsen and Nock (2013).

######  Theorem 3.

Suppose that the available samples are i.i.d. from $\pi=q/c$ and the standard
Gaussian distribution $\phi$. As $n_{1},n_{2}\to\infty$, suppose that
$n_{2}/n_{1}\to\beta(>0)$ and $\chi_{P}^{2}(\phi_{\text{mix}},\pi)\to 0$. Then

$\lim_{n_{1},n_{2}\to\infty}\frac{(n_{1}+n_{2})\mathrm{Var}(\hat{\lambda}_{\text{WB}})}{\chi_{P}^{2}(\phi,\tilde{\pi})}=1,$
(28)

and

$\lim_{n_{1},n_{2}\to\infty}\frac{(n_{1}+n_{2})\mathrm{Var}(\hat{\lambda}_{\text{SWB}})}{\sum_{k=1}^{K}\tilde{w}_{k}^{2}(1+\beta)/(\tilde{w}_{k}+\beta)\cdot\chi_{P}^{2}(\phi,\tilde{p}_{k})}=1.$
(29)

The notation $\tilde{\pi}$, $\tilde{p}_{k}$ and $\tilde{w}_{k}$ is defined in
(19).

MCMC draws generally have a complex dependence which is specific to the
sampling algorithm used. For simplification, Theorem 3 assumes i.i.d. samples.
Theorem 3 also assumes that the sample size ratio $n_{2}/n_{1}$ converges to a
constant, and that $\phi_{\text{mix}}$ approximates $\pi$ with increasing
accuracy as the sample size increases. The result suggests
$\mathrm{Var}(\hat{\lambda}_{\text{SWB}})$ can be expected to be smaller than
$\mathrm{Var}(\hat{\lambda}_{\text{WB}})$ in certain cases. Consider as an
illustrative case the scenario where $\beta\geq 1$, $K$ is a relatively large,
and the mixture components in $\tilde{\pi}$ have approximately equal weights
(i.e. $\tilde{w}_{k}\approx 1/K\approx 0$ such that
$\tilde{w}_{k}+\beta\approx\beta$). Then, the denominator in (29)
approximately equals

$\sum_{k=1}^{K}\tilde{w}_{k}^{2}(1+\beta)/(\tilde{w}_{k}+\beta)\cdot\chi_{P}^{2}(\phi,\tilde{p}_{k})\approx(1+1/\beta)\times\Big{[}\frac{1}{K}\sum_{k=1}^{K}\chi_{P}^{2}(\phi,\tilde{p}_{k})\Big{]}/K.$
(30)

The quantity inside the square brackets is the average of the component
divergences $\chi_{P}^{2}(\phi,\tilde{p}_{k})$, and this is further divided by
$K$ (which is large). Thus, intuitively, we can expect (30) to be smaller than
the denominator of (28) in this setting. In the following, we further exploit
Lemma 3 to rigorously characterise when
$\mathrm{Var}(\hat{\lambda}_{\text{SWB}})$ will be smaller than
$\mathrm{Var}(\hat{\lambda}_{\text{WB}})$ asymptotically.

Based on (28) and (29), the Warp-U bridge and stochastic Warp-U bridge
estimators have small asymptotic variance when $\tilde{\pi}$ or
$\tilde{p}_{k}$, respectively, is close to $\phi$. For our stochastic Warp-U
bridge estimator, the ideal case is that each mixture component
$\tilde{p}_{k}$ perfectly matches the transformation target $\phi$. In that
ideal case, the estimator would have zero asymptotic variance. In practice, we
can not expect such perfect agreement with $\phi$. We measure the discrepancy
by
$d_{k}(\boldsymbol{\theta})=\tilde{p}_{k}(\boldsymbol{\theta})/\phi(\boldsymbol{\theta})-1$,
and collect the weighted discrepancy for all densities $\tilde{p}_{k}$ in the
mixture (19) as a vector

$\displaystyle\boldsymbol{d}(\boldsymbol{\theta})=\big{(}\tilde{w}_{1}d_{1}(\boldsymbol{\theta}),\tilde{w}_{2}d_{2}(\boldsymbol{\theta}),\cdots,\tilde{w}_{K}d_{K}(\boldsymbol{\theta})\big{)}^{\top}.$
(31)

Note that the sign of $d_{k}(\boldsymbol{\theta})$ reflects the transformation
bias. When $d_{k}(\boldsymbol{\theta})>0$, the transformed density
$\tilde{p}_{k}$ overestimates $\phi$ at $\boldsymbol{\theta}$, whereas when
$d_{k}(\boldsymbol{\theta})<0$, it underestimates $\phi$ at
$\boldsymbol{\theta}$. The overall bias for the mixture distribution
$\tilde{\pi}$ at $\boldsymbol{\theta}$ is measured by
$\mathbf{1}^{\top}\boldsymbol{d}(\boldsymbol{\theta})$. Based on this
intuition, we decompose the discrepancy vector
$\boldsymbol{d}(\boldsymbol{\theta})$ as

$\boldsymbol{d}(\boldsymbol{\theta})=\mathbf{P}_{\mathbf{1}}\boldsymbol{d}(\boldsymbol{\theta})+(\mathbf{I}-\mathbf{P}_{\mathbf{1}})\boldsymbol{d}(\boldsymbol{\theta}).$
(32)

where $\mathbf{P}_{\mathbf{1}}=(1/K)\mathbf{1}\mathbf{1}^{\top}$ is the
orthonormal projection matrix onto $\text{span}\\{\mathbf{1}\\}$. Equation
(32) implies that the weighted sum of the component discrepancies
$\chi_{P}^{2}(\phi,\tilde{p}_{k})$ can be written as

$\displaystyle\sum_{k=1}^{K}\tilde{w}_{k}^{2}\chi_{P}^{2}(\phi,\tilde{p}_{k})=$
$\displaystyle\int\|\boldsymbol{d}(\boldsymbol{\theta})\|_{2}^{2}\phi(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})$
$\displaystyle=$
$\displaystyle\underbrace{\int\|\mathbf{P}_{\mathbf{1}}\boldsymbol{d}(\boldsymbol{\theta})\|_{2}^{2}\phi(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})}_{(\rm
I)}+\underbrace{\int\|(\mathbf{I}-\mathbf{P}_{\mathbf{1}})\boldsymbol{d}(\boldsymbol{\theta})\|_{2}^{2}\phi(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})}_{(\rm
II)}.$ (33)

In the above, the average magnitude of $\boldsymbol{d}(\boldsymbol{\theta})$
across all $\boldsymbol{\theta}$ is decomposed into two parts: term (I) equals
$\chi_{P}^{2}(\phi,\tilde{\pi})/K$ and reflects the bias of $\tilde{\pi}$
averaged across $\boldsymbol{\theta}$; term (II) measures the magnitude of the
remaining discrepancy between the Warp-U transformed components and $\phi$.
Term (I) may be interpreted as the “systematic bias”. The following theorem
reveals that when (I) is a non-negligible source of discrepancy, our
stochastic Warp-U bridge estimator will have better performance than the
Warp-U bridge estimator (Wang et al., 2022) in terms of asymptotic variance
and precision per second.

######  Theorem 4.

Suppose that $\beta\geq 1,K\geq 3$, and the conditions in Lemma 3 hold. Then,
when

$\int\|\mathbf{P}_{\mathbf{1}}\boldsymbol{d}(\boldsymbol{\theta})\|_{2}^{2}\phi(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})\geq\beta_{1K}\int\|(\mathbf{I}-\mathbf{P}_{\mathbf{1}})\boldsymbol{d}(\boldsymbol{\theta})\|_{2}^{2}\phi(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta}),$
(34)

with $\beta_{1K}=(\beta+1)/(\beta(K-1)-1)$, we have

$\limsup_{n_{1},n_{2}\to\infty}\big{[}\mathrm{Var}(\hat{\lambda}_{\text{SWB}})/\mathrm{Var}(\hat{\lambda}_{\text{WB}})\big{]}\leq
1;$ (35)

and

$\liminf_{n_{1},n_{2}\to\infty}\big{[}\mathrm{PpS}(\hat{\lambda}_{\text{SWB}})/\mathrm{PpS}(\hat{\lambda}_{\text{WB}})\big{]}\geq\frac{(1+\beta)K}{1+\beta
K}.$ (36)

When $\beta\geq 1$, we have $\beta_{1K}\leq 3/(K-1)$ and hence
$\beta_{1K}\|(\mathbf{I}-\mathbf{P}_{\mathbf{1}})\boldsymbol{d}(\boldsymbol{\theta})\|_{2}^{2}\leq
3/(K-1)\|(\mathbf{I}-\mathbf{P}_{\mathbf{1}})\boldsymbol{d}(\boldsymbol{\theta})\|_{2}^{2}$.
Based on this, a sufficient condition for (34) is that

$\int\|\mathbf{P}_{\mathbf{1}}\boldsymbol{d}(\boldsymbol{\theta})\|_{2}^{2}\phi(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta})\geq[3/(K-1)]\times\int\|(\mathbf{I}-\mathbf{P}_{\mathbf{1}})\boldsymbol{d}(\boldsymbol{\theta})\|_{2}^{2}\phi(\boldsymbol{\theta})\mu(\mathrm{d}\boldsymbol{\theta}).$
(37)

The right hand side of (37) measures the dimension-adjusted size of
$\boldsymbol{d}(\boldsymbol{\theta})$ inside $\text{span}(\mathbf{1})^{\perp}$
(multiplied by $3$), because
$\text{dim}\big{(}\text{span}(\mathbf{1})^{\perp}\big{)}=K-1$. The sufficient
condition (37) means that, adjusted by subspace dimension, the average size of
$\boldsymbol{d}(\boldsymbol{\theta})$ inside $\text{span}(\mathbf{1})$ should
be three times as large as the average size of
$\boldsymbol{d}(\boldsymbol{\theta})$ inside
$\text{span}(\mathbf{1})^{\perp}$.

In practice, condition (34) will be satisfied if the $\tilde{p}_{k}$’s tend to
jointly over-estimate (or under-estimate) $\phi$ across non-negligible regions
of $\Theta$. For example, if the target density $\pi$ is a mixture of
$t$-distribution while $\phi_{\text{mix}}$ is a mixture of Gaussians, then
each component of $\phi_{\text{mix}}$ could over-estimate around the modes and
under-estimate the tails of each corresponding component of the target
distribution. Systematic bias could potentially also arise from bias of the
estimation method used for fitting $\phi_{\text{mix}}$ to $\pi$. In all our
numerical studies, the stochastic Warp-U bridge estimator performs better than
the original Warp-U bridge estimator.

The precision per second comparison (36) takes the number of evaluations of
$q$ into account. Recall that the Warp-U bridge estimator requires
$(n_{1}+n_{2})K$ evaluations of $q$, while the stochastic Warp-U bridge
estimator requires $n_{1}+n_{2}K$ evaluations. As $n_{1}\to\infty$ and
$n_{2}/n_{1}\to\beta>0$, the ratio $(n_{1}+n_{2})K/(n_{1}+n_{2}K)$ converges
to $(1+\beta)K/(1+\beta K)$, which is the right hand side of (36).

##  4 Simulation Studies

We present simulation studies to illustrate the effectiveness of our
methodology: Section 4.1 focuses on the Warp-U MCMC sampling method described
by Algorithm 3, and Section 4.3 focuses on the stochastic Warp-U bridge
estimation method described by Algorithm 5 when combined with various sampling
strategies. Our main simulation studies include 4-dimensional, 10-dimensional,
and 30-dimensional examples. The empirical performance of unbiased Warp-U
estimation (Algorithm 4) is examined in Section 4.4. In Section 4.5, we
present some challenges and solutions in the context of the 50-dimensional
setting discussed in Wang et al. (2022), which has complex density contours
and heavy tails.

###  4.1 Comparison of Sampling Methods

We now consider a 4-dimensional multivariate Gaussian mixture target
distribution $\pi$ and compare the performance of our Warp-U MCMC sampling
method to Parallel Tempering (PT) and the Metropolis-Hastings (MH) algorithm
with a fixed proposal distribution. PT (Geyer, 1991) is a widely used method
designed to sample from multimodal target distributions, but it is necessary
to tune the temperature levels, which can be time-consuming. For tuning the
temperature levels, we apply the method of Atchadé et al. (2011) which
iteratively selects the inverse temperatures by the stochastic approximation
algorithm given in Robbins and Monro (1951) and Andrieu and Robert (2001). For
MH, we set the proposal density to be $\phi_{\text{mix}}$, i.e., the same
density used by our Warp-U MCMC sampler. The unnormalised target distribution
is given by

$\displaystyle q(\boldsymbol{\theta})$
$\displaystyle=\frac{1}{15}\exp(-\frac{1}{2}(\boldsymbol{\theta}+\textbf{11})^{\top}(\boldsymbol{\theta}+\textbf{11}))+\frac{2}{15}\exp(-\frac{1}{2}(\boldsymbol{\theta}-\textbf{12})^{\top}(\boldsymbol{\theta}-\textbf{12}))$
$\displaystyle+\frac{3}{15}\exp(-\frac{1}{2}(\boldsymbol{\theta}+\textbf{8})^{\top}(\boldsymbol{\theta}+\textbf{8}))+\frac{4}{15}\exp(-\frac{1}{2}(\boldsymbol{\theta}-\textbf{7})^{\top}(\boldsymbol{\theta}-\textbf{7}))$
$\displaystyle+\frac{5}{15}\exp(-\frac{1}{2}(\boldsymbol{\theta}+\textbf{2})^{\top}(\boldsymbol{\theta}+\textbf{2})),$

where bold numbers indicate the vectors in $\mathbb{R}^{4}$, e.g.,
$\textbf{2}=(2,2,2,2)^{\top}$. Note that some of the modes (especially
$-\textbf{2}$ and 7) are well separated, and therefore a simple Metropolis-
Hastings algorithm with a Gaussian proposal is not expected to work well.

Figure 4: The top left, top right, and bottom left panels show the true target
density (red line) and the density estimated from samples obtained using our
adaptive Warp-U MCMC sampler (middle dash-dot blue line) before the first,
second, and eighth stage of Algorithm 3, respectively. The blue shaded regions
show pointwise $95\%$ confidence intervals (marked by the upper and lower
dash-dot blue lines). The bottom right panel shows the square of the 1/log
Wasserstein distance (solid blue line) between the true marginal density (the
first dimension of the target) and the estimated marginal density at different
stages of the algorithm, as well as $95\%$ confidence intervals (solid orange
lines) and a linear fit (dashed red line).

We apply Algorithm 3 with $T=4000$ (number of samples at each stage), $M=11$
(number of stages), and $K=10$ (number of components). We repeat the
simulation $50$ times, and use kernel density estimation to infer the target
density for each simulation, and thereby obtain approximate pointwise $95\%$
confidence intervals. Figure 4 shows the inferred density and pointwise $95\%$
confidence interval (middle dash-dot blue line and blue shaded regions) based
on samples obtained at different stages of our Warp-U MCMC sampling algorithm:
(i) based on the initial samples (top left panel), (ii) based on samples
obtained in stage $1$ (top right panel), and (iii) based on stage $7$ samples
(bottom left panel), respectively. The solid red lines show the true target
density. We can see that for the initial samples, the confidence interval is
highly inconsistent with the target density (since the initial sampling is
from a uniform distribution), but that the confidence interval at stage $7$
covers the target density well ($88.3$% coverage). We found that the square of
the inverse of the log Wasserstein distance between the true marginal density
and the estimated marginal is well approximated by a linear function of the
stage number, see the bottom right of Figure 4.

Figure 5 shows the log Wasserstein distance between the true target density
and the density estimated from samples, as a function of target evaluations,
when the sampling method is our adaptive Warp-U MCMC method (solid blue lines)
or PT (dashed orange lines). Figure 5 shows that the samples obtained via our
adaptive Warp-U MCMC sampler better approximate the true target density than
those obtained via PT, for a given number of target evaluations. In other
words, in this example, our approach has lower computational cost than PT for
a given level of accuracy.

Figure 5: Log Wasserstein distance between the true target density and
estimated density constructed using our adaptive Warp-U MCMC method (centre
solid blue line) and PT (centre dashed orange line), as well as pointwise
$95\%$ confidence intervals (shaded regions).

Figure 6 shows the trace plot and auto-correlation function (ACF) for our non-
adaptive Warp-U MCMC sampler specified in Algorithm 2 (left column) and the
Metropolis-Hastings (MH) algorithm with proposal density $\phi_{\text{mix}}$
(right column). We use the non-adaptive version of our sampler to avoid having
to change the proposal density of MH algorithm that we are comparing to at
each stage and simplify the comparison. The Gaussian mixture distribution
$\phi_{\text{mix}}$ used for this non-adaptive version is given by
(LABEL:target_twomode) in Section LABEL:appendix:setting of the online
Supplementary Material. The trace and ACF plots show that the Warp-U
transformation used in Algorithm 2 greatly reduces the dependence between
samples compared with the Metropolis-Hastings approach, and facilitates jumps
between modes leading to a larger effective sample size for the Warp-U MCMC
sampler.

Figure 6: Trace and autocorrelation plots for our Warp-U MCMC method (left
column) and MH with proposal density $\phi_{\text{mix}}$ (right column). The
y-axis of the panels in the first row is the first dimension of the samples
from the target density, i.e., $\boldsymbol{\theta}_{[1]}$.

###  4.2 Comparison of Sampling Methods in High Dimension

A notable limitation of PT is that it converges slowly for high dimensional
mixture distributions whose components have unequal variances, as Woodard et
al. (2009) demonstrated theoretically. Woodard et al. (2009) provided a
theoretical demonstration of PT’s slow convergence in certain high-dimensional
contexts, especially when the target distribution is a mixture of Gaussians
with distinct variances. We examine the performance of PT ans Warp-U sampler
in two related settings.

In the first setting, we adopt the target distribution in Woodard et al.
(2009):

$\displaystyle q(\boldsymbol{\theta})$
$\displaystyle=\frac{1}{2\bar{\sigma}_{1}^{d}}\exp(-\frac{1}{2\bar{\sigma}_{1}^{2}}(\boldsymbol{\theta}+\textbf{1})^{\top}(\boldsymbol{\theta}+\textbf{1}))+\frac{1}{2\bar{\sigma}_{2}^{d}}\exp(-\frac{1}{2\bar{\sigma}_{2}^{2}}(\boldsymbol{\theta}-\textbf{1})^{\top}(\boldsymbol{\theta}-\textbf{1})),$

and we compare the performance of PT and Warp-U sampler as the dimension $d$
increases. Specifically, in this experiment, we set $\bar{\sigma}_{1}^{2}=0.8$
and $\bar{\sigma}_{2}^{2}=0.2$, and vary the dimension
$d\in\\{5,10,30,50,80,100,300,1000\\}$.

For PT, we run 20 chains with equally spaced inverse temperature. On the other
hand, recall the Warp-U sampler in Algorithm 2 requires the specification of
the auxiliary density $\phi_{\rm mix}$. The adaptive version in Algorithm 3 is
less sensitive to the particular form of $\phi_{\rm mix}$, which is only used
for initialisation. However, for simplicity, in this experiment, we test the
Warp-U sampler in Algorithm 2 with two choices of $\phi_{\rm mix}$ as follows.

Given $q$, a simplistic choice of the auxiliary distribution would nbe

$\phi_{\text{mix},1}(\boldsymbol{\theta})=\frac{1}{2}\phi(\boldsymbol{\theta}-\boldsymbol{\mu}_{1})+\frac{1}{2}\phi(\boldsymbol{\theta}-\boldsymbol{\mu}_{2}),$
(38)

where $\boldsymbol{\mu}_{1}=-\textbf{1}$ and
$\boldsymbol{\mu}_{2}=\textbf{1}$. Here the modes
($\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2}$) can be estimated from random
initial points. It is a simplistic choice because it ignores the variances. A
more principled choice would attempt to model the unknown variances, such as

$\phi_{\text{mix},2}(\boldsymbol{\theta})=\int\frac{1}{2\sigma_{1}^{d}}\phi((\boldsymbol{\theta}-\boldsymbol{\mu}_{1})/\sigma_{1})p(\sigma_{1}^{2})\,\mathrm{d}\sigma_{1}^{2}+\int\frac{1}{2\sigma_{2}^{d}}\phi((\boldsymbol{\theta}-\boldsymbol{\mu}_{2})/\sigma_{2})p(\sigma_{2}^{2})\,\mathrm{d}\sigma_{2}^{2},$
(39)

which enables us to account for some uncertainty regarding the variances. In
this experiment, we set $p(\sigma_{1}^{2})$ and $p(\sigma_{2}^{2})$ to be the
inverse Gamma distribution, $\text{Inv-Gamma}(a,b)$, with shape parameter
$a=2.25$ and scale parameter $b=1.25$. Note, when $\sigma_{1}^{2}\sim
p(\sigma_{1}^{2})$, we have $\mathbb{E}(\sigma_{1}^{2})=1$ and
$\mathrm{Var}(\sigma_{1}^{2})=4$. This reflects our weak knowledge of the
variances in the target $q$. When the Warp-U sampler runs with (39), in
addition to the primiary mixture component index $\psi\in\\{1,2\\}$, we also
need to sample the variance $\sigma_{\psi}^{2}$. The implementation details
for the Warp-U sampler with (39) are provided in Section
LABEL:appendix:numerical:pt of the Online Supplements.

Figure 7: Comparison of sampling methods in high dimension. The Wasserstain
distance between MCMC samples and true target i.i.d samples for Setting 1
(left panel) and Setting 2 (right panel). The distance is computed for the
joint marginal $(\theta_{1},\theta_{d})$ and is averaged across 20 simulation
replicates.

Figure 8: Comparison of PT and Warp-U sampler in Setting 1. Each panel shows
the density plot for $(\theta_{1},\theta_{d})$, which is drawn by one run of a
MCMC sampler. Each panel also shows the marginal histogram for the Markov
chain samples, and the black curves are the true marginal density. The three
columns are for PT, Warp-U sampler with $\phi_{\text{mix},1}$ (see (38)), and
Warp-U sampler with $\phi_{\text{mix},2}$ (see (39)), respectively. The four
rows correspond to dimension $d=10$, 30, 100, and 1000, respectively.

Figure 9: Comparison of PT and Warp-U sampler in Setting 2. Each panel shows
the density plot for $(\theta_{1},\theta_{d})$, which is drawn by one run of a
MCMC sampler. Each panel also shows the marginal histogram for the Markov
chain samples, and the black curves are the true marginal density. The three
columns are for PT, Warp-U sampler with $\phi_{\text{mix},1}$ (see (41)), and
Warp-U sampler with $\phi_{\text{mix},2}$ (see (42)), respectively. The four
rows correspond to dimension $d=10$, 30, 100, and 1000, respectively.

For each specified value of $d$, we run each algorithm 20 times to understand
their convergence properties. In each run, the initial sample is drawn from
the standard Gaussian distribution, and the first 1000 samples are discarded,
and the subsequent 4000 samples are retained. The simulation outcomes from one
run are depicted in Figure 8. In the figure, each panel shows the density plot
for $(\theta_{1},\theta_{d})$, i.e., the joint density plot for the first and
the last dimension. Note each 2-dimensional density plot is for MCMC samples.
Each panel also shows the marginal histogram for MCMC samples, and the black
curves represent the target marginal densities. The three columns are for PT,
Warp-U sampler with $\phi_{\text{mix},1}$ in (38), and Warp-U sampler with
$\phi_{\text{mix},2}$ in (39), respectively. The four rows correspond to
dimension $d=10$, 30, 100, and 1000, respectively. When the dimension equals
$d=10$, all three algorithms accurately detects two local modes of the target
density. When $d=30$, the Warp-U samples depicts the two modes, while PT fails
to do so. When $d\geq 100$, both PT and the Warp-U sampler with
$\phi_{\text{mix},1}$ fail to discover two local modes, while the samples from
Warp-U sampler with $\phi_{\text{mix},2}$ maintain a good representation of
the target.

The left panel of Figure 7 presents the Wasserstein distance between i.i.d
samples from the target joint density and the Markov chain samples, for the
two dimensions $(\theta_{1},\theta_{d})$. The distance is averaged across 20
replicates. As the dimension $d$ increases, the Warp-U sampler using (39)
consistently achieves better accuracy than PT. While the performance of the
Warp-U sampler with (38) does decline as $d$ increases, it still outperforms
PT in moderate dimensions.

We note that the average computing time for each run of PT, Warp-U sampler
with $\phi_{\text{mix},1}$ in (38), and Warp-U sampler with
$\phi_{\text{mix},2}$ are 7.44 min, 0.34 min, 47.68 min, respectively. Warp-U
sampler with $\phi_{\text{mix},2}$ took substantially longer because we need
to use metropolis algorithm to draw the inverse Warp-U transformation. But
this significantly longer running time is worthwhile in the sense that running
the first two algorithms for this long would not get them out of their trapped
regions.

In the second simulation setting, we consider the target distribution:

$\displaystyle q(\boldsymbol{\theta})$
$\displaystyle=\sum_{k=1}^{2}\frac{1}{2|\overline{\mathbf{D}}_{k}|^{1/2}}\exp\Big{(}-\frac{1}{2}(\boldsymbol{\theta}-\boldsymbol{\mu}_{k})^{\top}\overline{\mathbf{D}}_{k}^{-1}(\boldsymbol{\theta}-\boldsymbol{\mu}_{k})\Big{)}.$
(40)

The above is also a Gaussian mixture distribution with two components. The
elements of the first mean vector $\boldsymbol{\mu}_{1}$ are drawn from
Uniform$(-2.5,-1.5)$; the elements of the second mean vector
$\boldsymbol{\mu}_{2}$ are drawn from Uniform$(1.5,2.5)$. Besides, for
$k=1,2$, the covariance matrix
$\overline{\mathbf{D}}_{k}=\mathrm{diag}(\overline{\mathbf{D}}_{k1},\overline{\mathbf{D}}_{k2},\ldots,\overline{\mathbf{D}}_{k5})\in\mathbb{R}^{d\times
d}$ is diagonal of size $d$. For $j=1,\ldots,5$, the block
$\overline{\mathbf{D}}_{kj}=\mathrm{diag}(\bar{\sigma}_{kj}^{2},\bar{\sigma}_{kj}^{2},\ldots,\bar{\sigma}_{kj}^{2})\in\mathbb{R}^{(d/5)\times(d/5)}$
has exactly the same value for its diagonal elements. Specifically, we use

$\displaystyle(\bar{\sigma}_{11}^{2},\bar{\sigma}_{12}^{2},\bar{\sigma}_{13}^{2},\bar{\sigma}_{14}^{2},\bar{\sigma}_{15}^{2})$
$\displaystyle=(0.25,0.3,0.35,0.4,0.45),\text{ and }$
$\displaystyle(\bar{\sigma}_{21}^{2},\bar{\sigma}_{22}^{2},\bar{\sigma}_{23}^{2},\bar{\sigma}_{24}^{2},\bar{\sigma}_{25}^{2})$
$\displaystyle=(1,0.95,0.9,0.85,0.8),$

for the two mixture components, respectively.

In this setting, we run Warp-U sampler with the auxiliary distribution:

$\phi_{\text{mix},1}(\boldsymbol{\theta})=\frac{1}{2}\phi(\boldsymbol{\theta}-\boldsymbol{\mu}_{1})+\frac{1}{2}\phi(\boldsymbol{\theta}-\boldsymbol{\mu}_{2}),$
(41)

where the mean vectors $\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2}$ are the
same as in (40). Meanwhile, we will also run Warp-U sampler with another
auxiliary distribution:

$\phi_{\text{mix},2}(\boldsymbol{\theta})=\sum_{k=1}^{2}\int\frac{1}{2|\mathbf{D}_{k}|^{1/2}}\phi(\mathbf{D}_{k}^{-1/2}(\boldsymbol{\theta}-\boldsymbol{\mu}_{k}))\prod_{j=1}^{5}p(\sigma_{kj}^{2})\,\mathrm{d}\sigma_{k1}^{2}\ldots\,\mathrm{d}\sigma_{k5}^{2}.$
(42)

In the above,
$\mathbf{D}_{k}=\mathrm{diag}(\mathbf{D}_{k1},\mathbf{D}_{k2},\ldots,\mathbf{D}_{k5})\in\mathbb{R}^{d\times
d}$ has exactly the same diagonal structure as $\overline{\mathbf{D}}_{k}$ in
(40), but with unknown values of $\sigma_{kj}^{2}$ on the diagonal of
$\mathbf{D}_{kj}$. Each $p(\sigma_{kj}^{2})$ is an inverse Gamma distribution
with with shape parameter $a=2.25$ and scale parameter $b=1.25$.

The result for one simulation replicate is shown in Figure 9. The
interpretation is similar as that for Figure 8. For PT in the first column,
when $d=10$, though it discovers two local modes of the target density, the
Markov chain samples under-represent the mode with smaller variances. When
$d\geq 30$, PT samples only discover the mode with larger variance. For Warp-U
sampler with $\phi_{\text{mix},1}$ in the second column, its samples are no
longer an accurate portrayal of the two local modes when $d\geq 100$. For
Warp-U sampler with $\phi_{\text{mix},2}$ in the third column, its samples
consistently maintain a good representation of the target. In the right panel
of Figure 7, we show the Wasserstain distance between the MCMC sampled
$(\theta_{1},\theta_{d})$ and the i.i.d samples from the target density. The
distance is averaged across 20 simulation replicates. Warp-U sampler with
$\phi_{\text{mix},2}$ consistently generates high-quality samples across
distinct specifications of the dimension.

###  4.3 Comparison of Combined Sampling and Estimation Methods

Table 3: Summary statistics of the log target density evaluated at $10^{6}$ samples, where the target density is a 10-dimensional mixed skew-$t$ distribution. $\log(q)$ | Min. | 1st Qu. | Median | Mean | 3rd Qu. | Max.
---|---|---|---|---|---|---
Target samples | -171.72 | -27.86 | -23.47 | -24.67 | -20.11 | -10.26
Uniform samples | -53.07 | -45.23 | -43.71 | -43.52 | -42.01 | -26.65

We now compare our stochastic bridge estimation approach given by Algorithm 5
with bridge sampling (Meng and Wong, 1996) and Warp-U bridge estimation (Wang
et al., 2022). The target density we consider is that used in Wang et al.
(2022), and in particular is a $10$ dimensional mixture of $25$ multivariate
skew-$t$ distributions whose mass mainly lies on the support $[-20,20]^{10}$.
A similar target density was first introduced in a simulation study designed
by Azzalini (2013). The target density presents a challenge in that large
regions of its support have low density. To illustrate this, we generated
$10^{6}$ samples first from the target density itself and then from a uniform
distribution on $[-20,20]^{10}$. In the case of the uniform samples, we found
that $99.99\%$ had target density value lower than one-hundredth of the median
target density value of the samples directly obtained from the target density,
see Table 3.

The three estimation methods we consider are all designed to estimate
normalising constants given samples from the target density, and we
investigate their performance when the samples are obtained via three
different approaches: (i) our Warp-U MCMC sampler given by Algorithm 3; (ii)
the Generalised Wang-Landau (GWL) algorithm (Liang, 2005); and (iii) PT. The
GWL algorithm does not directly sample from the target density and therefore
in the case of (ii) we apply the additional sub-sampling step suggested in
Section 5.2 of Wang et al. (2022); see Section 6.1 for further discussion. We
set the number of samples from the target density to be the same for all three
sampling methods; specifically, we collect $26,000$ samples under each method,
as in Wang et al. (2022).

Figure 10 shows the root mean square error (RMSE) of the estimator as a
function of the number of target evaluations for the classical bridge sampling
(BS) estimator (green circles, and dash-dot line), the Warp-U bridge (WB)
estimator (red triangles, and dotted line), and our stochastic Warp-U bridge
(SWB) estimator (blue squares, and solid line), when the samples are provided
by our adaptive Warp-U MCMC sampler (left), GWL (middle), and PT (right). Our
stochastic bridge estimator has the smallest RMSE regardless of the sampling
method used, except when the number of target evaluations is small. Whether
using GWL or Warp-U MCMC samples, the results are similar and slightly better
than with PT samples. On the other hand, GWL does not directly generate
samples from the target distribution, and additional sub-sampling steps (see
Section 5.2 of Wang et al. (2022)) are needed to get the target samples, so
the total computational cost is higher than the other sampling methods (only
the cost of estimation is illustrated in Figure 10).

Figure 10: Root mean square error (RMSE) of three normalising constant
estimators as a function of the number of target evaluations, where the
samples are obtained using our adaptive Warp-U MCMC sampler (left), GWL
(middle), and PT (right). The three estimators are the classical bridge
sampling (BS) estimator (green circles, and dash-dot line), the Warp-U bridge
(WB) estimator (red triangles, and dotted line), and our stochastic Warp-U
bridge (SWB) estimator (blue squares, and solid line). The error bars show $2$
times the standard errors. The $x$-axis gives the number of estimation stage
target evaluations as $\log_{10}$ of the proportion of target evaluations
needed for an 11th stage of the GWL algorithm.

We now return to Table 2 to compare the computational cost of the three
estimation methods. In Table 2, $M$ denotes the number of stages in our
adaptive Warp-U MCMC sampler and the GWL algorithm, and $M_{l}$ denotes the
number of temperature levels for PT. For the GWL algorithm, the number of
samples collected at each stage is different, so we use $\tilde{n}_{1}$ to
denote the average number of samples across all stages (in practice,
$\tilde{n}_{1}$ is relatively large compared with $n_{1}$). In addition, burn-
in is needed for all three sampling methods. Note that the number of
iterations given in Table 2 is not necessarily the number of unique samples,
because some proposed samples are rejected. Table 2 shows that classical
bridge sampling and stochastic Warp-U bridge estimation both require no more
computation than Warp-U bridge estimation, even when some estimation stage
computation is saved by the use our Warp-U MCMC sampling method. The last row
of Table 2 shows the computation cost of the three estimation methods given
the samples, i.e., without considering the computation saved in the sampling
stage. It shows that stochastic Warp-U bridge estimation required no more
computation than Warp-U bridge estimation when all samples are given.

###  4.4 Warp-U Sampler Combined with an Unbiased Estimator

In Section 2.4 we introduced genuinely unbiased estimation strategies for the
multi-modal distribution context which do not rely on checks for the
convergence of the MCMC sampler. We now examine the performance of these
strategies, specifically the Warp-U unbiased estimation methods presented in
Algorithm 4 (separate coupling for the Warp-U transformation and its inverse)
and Algorithm LABEL:alg:unbiasedV2 (combined coupling). We also consider the
three estimators (LABEL:eqn:unbiased:raoLevel0)–(LABEL:eqn:unbiased:raoLevel2)
in the online Supplementary Material, which correspond to distinct levels of
Rao-Blackwellization. The three estimators are denoted as L0, L1 and L2
respectively. L0 corresponds to the original estimator in (16)-(17) without
Rao-Blackwellization. L1 exploits the conditional expectations
$\mathbb{E}\big{[}h(\boldsymbol{\theta}_{1,t})|\boldsymbol{\theta}_{1}^{*}\big{]}$
and
$\mathbb{E}\big{[}h(\boldsymbol{\theta}_{2,t})|\boldsymbol{\theta}_{2}^{*}\big{]}$
for Rao-Blackwellization, while L2 uses the conditional expectations
$\mathbb{E}\big{[}h(\boldsymbol{\theta}_{1,t})|\boldsymbol{\theta}^{\text{\tiny{MH}}}_{1}\big{]}$
and
$\mathbb{E}\big{[}h(\boldsymbol{\theta}_{2,t})|\boldsymbol{\theta}^{\text{\tiny{MH}}}_{2}\big{]}$.
See Section LABEL:appendix:unbiased of the online Supplementary Material for
more details. In this subsection we compare the six candidate methods
corresponding to the two algorithms (Algorithm 4 and Algorithm
LABEL:alg:unbiasedV2) and three levels of Rao-Blackwellization (L0, L1, and
L2).

We consider two settings which are of higher dimension and more challenging
than those in the previous subsections. In these two simulation settings, the
target distribution is a mixture of multivariate Gaussian distributions and a
mixture of multivariate skewed $t$-distributions (Gupta, 2003), respectively.
In both settings there are 5 mixture components and the dimension of
$\boldsymbol{\theta}$ is 30. The multivariate skewed $t$-distributions have 12
degree of freedom, and the elements of the skewness vector take values between
$-100$ to $200$. For the two mixture settings, density contours for the first
4 dimensions are plotted in Figure LABEL:fig:unbiased:contour01 and Figure
LABEL:fig:unbiased:contour02 in the online Supplementary Material,
respectively.

The high dimension of $\boldsymbol{\theta}$ makes it unlikely that the coupled
chains will meet in the random walk Metropolis-Hastings step of Algorithm 4
and Algorithm LABEL:alg:unbiasedV2, which is consistent with the similar
observations made by Jacob et al. (2020). Therefore, in our simulation study,
all the random walk Metropolis-Hasting steps are replaced by the Hamiltonian
Monte Carlo step proposed by Heng and Jacob (2019). The simulation is repeated
200 times for each setting.

Table 4: Average value and 10%, 30%, 50%, 70%, 90% quantile for the meeting
time $\tau$ in the two mixture settings. Two Warp-U coupling schemes are
compared: separate (Algorithm 4) and combined (Algorithm
LABEL:alg:unbiasedV2).

Setting | Mixture Component | Warp-U Coupling | Average | 10% | 30% | 50% | 70% | 90%
---|---|---|---|---|---|---|---|---
1 | Gaussian | Separate | 21.9 | 7 | 11 | 17 | 24 | 46
|  | Combined | 20.8 | 7 | 10 | 16 | 24 | 40
2 |  Skewed $t$-distribution | Separate | 70.1 | 26 | 42 | 58 | 82 | 128
|  | Combined | 62.9 | 26 | 36 | 52 | 74 | 115

Table 4 shows that the median meeting time ($\tau:=\inf\left\\{t\geq
1:\boldsymbol{\theta}_{1,t}=\boldsymbol{\theta}_{2,t-1}\right\\}$) is less
than 20 and 60 for the two settings, respectively. The combined Warp-U
coupling scheme (Algorithm LABEL:alg:unbiasedV2) can further reduce the
meeting time slightly, see rows two and four in Table 4. If we directly adopt
the method of Jacob et al. (2020) without our Warp-U sampler, then the two
chains generally get stuck in two different local modes and in practice are
unlikely to meet at all. In conclusion, the proposed algorithms are effective
for encouraging the meeting of the two chains in a multi-modal distribution
context. Fast meeting of the two chains is important for the unbiased
estimation scheme of Jacob et al. (2020) both theoretically and practically.
Theoretically, the scheme requires the tail of the random meeting time to
decay geometrically (see requirement R2 listed in Section 2.4). The practical
side is that, after the two chains meet, we only need to run a single chain
which saves computational resources.

The results of the first simulation setting (target distribution $\pi$ is a
mixture of Gaussian distribution) are shown in Figure
LABEL:fig:unbiased:s1h1Iter. Specifically, Figures
LABEL:fig:unbiased:s1h1Iter:a–LABEL:fig:unbiased:s1h1Iter:b present the
results for unbiased estimation of $\mathbb{E}_{\pi}h(\boldsymbol{\theta})$ by
$H_{l:m}(\boldsymbol{\theta}_{1,\cdot}\,,\,\boldsymbol{\theta}_{2,\cdot})$ in
(17), where $h(\boldsymbol{\theta})=\sum_{j=1}^{50}\theta_{j}$. The true value
is $\mathbb{E}_{\pi}h(\boldsymbol{\theta})\approx-4.786$. Figures
LABEL:fig:unbiased:s1h2Iter:a–LABEL:fig:unbiased:s1h2Iter:b present the
results for $h(\boldsymbol{\theta})=\sum_{j=1}^{50}\theta_{j}^{2}$ with true
value $\mathbb{E}_{\pi}h(\boldsymbol{\theta})\approx 182.939$.

In Figure LABEL:fig:unbiased:s1h1Iter:a, the two rows correspond to $l=500$
and $l=1000$, where $l$ is the subscript of
$H_{l:m}(\boldsymbol{\theta}_{1,\cdot}\,,\,\boldsymbol{\theta}_{2,\cdot})$.
The first column contains the results of Algorithm 4, while the second column
shows the results of Algorithm LABEL:alg:unbiasedV2. Each panel plots root
mean squared error (RMSE) versus the number of iteration (i.e. $m$ in the
subscript of $H_{l:m}$). In each panel, distinct curves represent different
levels of Rao-Blackwellization (L0, L1 and L2). From the figure, we can find
the estimation achieves high accuracy. Moreover, Rao-Blackwellization (L1 and
L2) improves the accuracy over the original estimator (L0). The performance
between L1 and L2 is similar. Meanwhile, the combined Warp-U coupling in
Algorithm LABEL:alg:unbiasedV2 has accuracy close to the separate Warp-U
coupling in Algorithm 4.

Figure LABEL:fig:unbiased:s1h1Iter:b displays how RMSE decreases with
computational time in seconds. The RMSE is computed for the estimator
$H_{l:m}$ with a proper value of $m$, such that the computational time falls
within a certain time interval on the horizontal axis. Recall from Figure
LABEL:fig:unbiased:s1h1Iter:a that L1 and L2 levels of Rao-Blackwellization
have similar performance as the number of iteration increases. However, in
Figure LABEL:fig:unbiased:s1h1Iter:b, RMSE for L2 decreases slower as the
computation time increases. This is due to the heavier computational burden
associated with L2, which requires $K^{2}$ evaluations of the target $q$. On
the other hand, L1 does not significantly increase the computational cost
compared with the original estimator L0. This is because most of the necessary
computation, which involves $K$ evaluations of the target $q$, has already
been completed during the Warp-U sampling phase. Our findings suggest that in
most cases the L1 Rao-Blackwellization is as good or better than the L0 and L2
versions in terms of the computation time needed to reduce RMSE below a given
threshold.

Lastly, Figure LABEL:fig:unbiased:s1h2Iter:a–LABEL:fig:unbiased:s1h2Iter:b
carry a similar message for the estimation of
$\mathbb{E}_{\pi}h(\boldsymbol{\theta})$ with
$h(\boldsymbol{\theta})=\sum_{k=1}^{50}\theta_{k}^{2}$. The results of the
second simulation setting (where the target distribution $\pi$ is a mixture of
multivariate skewed $t$-distribution) are shown in Figures
LABEL:fig:unbiased:s2h1Iter. The numerical results have a similar
interpretation as those of the first setting.

###  4.5 Sampling and Estimation in a Challenging 50-Dimensional Setting

To test the limit of the proposed method, we consider the 50 dimensional
simulation setting in Section 4.3 of Wang et al. (2022). The target
distribution $\pi$ is a mixture of 30 distributions, including normal
distributions, $t$-distributions (including Cauchy distributions), and
multivariate distributions with gamma and/or exponential marginal
distributions and normal copulas. In the work of Wang et al. (2022),
independent and identically distributed samples are drawn from $\pi$ and fed
to the Warp-U bridge estimator. In this subsection, we endeavour to draw
samples via our proposed adaptive Warp-U MCMC sampler (Algorithm 3) and then
estimate the normalising constant.

The left panel of Figure 11 shows marginal two-dimensional histograms of i.i.d
samples from the target distribution. Because the target distribution contains
Cauchy distributions and $t$-distributions as its mixture components, some
i.i.d. samples take extreme values. The plot only shows the density around the
origin. Since the target distribution has heavy tails and complex density
contours in 50 dimensions, obtaining samples from the target distribution
using MCMC is very challenging. Direct application of the adaptive Algorithm 3
with the basic sampler in Algorithm 2 will not work effectively. One obstacle
is that in 50 dimensions the random walk Metropolis Hasting step in Algorithm
2 is ineffective, and therefore needs to be replaced by a Hamiltonian Monte
Carlo (HMC) step. Even with this replacement, the Warp-U sampler can still get
stuck in local modes.

For illustration, we set the auxiliary distribution $\phi_{\text{mix}}$ to
have $K=200$ components. We initialise
$\phi_{\text{mix}}(\boldsymbol{\theta})$ by drawing its mean vector
$\boldsymbol{\mu}_{k}$ from a uniform distribution
$\text{Uniform}\big{(}[-15,15]^{50}\big{)}$, and setting the covariance matrix
of each mixture component of $\phi_{\text{mix}}(\boldsymbol{\theta})$ to be
$3I$. Lastly, we set the mixture components of
$\phi_{\text{mix}}(\boldsymbol{\theta})$ to have equal weights, i.e.
$w_{k}=1/K$. This is a typical $\phi_{\text{mix}}(\boldsymbol{\theta})$ to
choose when we have little prior knowledge about the target distribution
before the adaptive algorithm begins. The resulting Warp-U MCMC samples are
shown in the right panel of Figure 11. The figure indicates the Warp-U sampler
with this initial $\phi_{\text{mix}}(\boldsymbol{\theta})$ fails to explore
the full sample space.

Figure 11: Left panel: two-dimensional histogram of i.i.d. samples from the
target distribution for some selected pairs of dimensions. Right panel: Warp-U
MCMC samples (blue points) with a $\phi_{\text{mix}}$ whose mixture component
mean vectors are drawn from a uniform distribution.

Figure 12:  The Warp-U transformation and its inverse tries to shift the
current sample to one of
$\mathcal{H}_{\psi^{\prime}}\circ\mathcal{F}_{\psi}\big{(}\boldsymbol{\theta}^{\text{HMC}}\big{)}$
($\psi,\psi^{\prime}=1,2,\cdots,200$). Left panel: the black point is the
current sample $\boldsymbol{\theta}^{\text{HMC}}$, and the blue points
represent 200 candidates with a fixed $\psi$ and $\psi^{\prime}=1,\cdots,200$.
The background is the two-dimensional histogram for i.i.d. samples. Right
panel: the heatmap of
$\log\big{(}\pi\big{(}\mathcal{H}_{\psi^{\prime}}\circ\mathcal{F}_{\psi}\big{(}\boldsymbol{\theta}^{\text{HMC}}\big{)}\big{)}\big{)}$
for all the 400,00 candidate points with $\psi,\psi^{\prime}=1,2,\cdots,200$.

The reason that Warp-U transformation fails to shift samples across modes can
be explained by Figure 12. In the left panel, the black point is a sample
$\boldsymbol{\theta}^{\text{HMC}}$ generated by the HMC step. The Warp-U
transform and its inverse essentially propose to shift
$\boldsymbol{\theta}^{\text{HMC}}$ to one of the following points

$\boldsymbol{\theta}_{\psi,\psi^{\prime}}=\mathcal{H}_{\psi^{\prime}}\circ\mathcal{F}_{\psi}\big{(}\boldsymbol{\theta}^{\text{HMC}}\big{)}=\boldsymbol{\theta}^{\text{HMC}}-\boldsymbol{\mu}_{\psi}+\boldsymbol{\mu}_{\psi^{\prime}},\qquad\psi,\psi^{\prime}=1,2,\cdots,K.$

Among these $K^{2}=40,000$ candidates, 200 candidates (with a fixed $\psi$ and
$\psi^{\prime}=1,\cdots,200$) are shown as blue points in the right panel. The
logarithm of the target density
$\log\pi(\boldsymbol{\theta}_{\psi,\psi^{\prime}})$ at all 40,000 candidates
(with $\psi,\psi^{\prime}=1,\cdots,200$) is shown in the right panel of Figure
12. We can see the density at $\psi=\psi^{\prime}$ is much larger than for
$\psi\not=\psi^{\prime}$. This indicates that for $\psi\not=\psi^{\prime}$ the
Warp-U transformation usually shifts the current
$\boldsymbol{\theta}^{\text{HMC}}$ to a low density region. Therefore, with
$\phi_{\text{mix}}(\boldsymbol{\theta})$ constructed as described above, the
Warp-U sampler only chooses the identity transformation
$\mathcal{H}_{\psi}\circ\mathcal{F}_{\psi}\big{(}\boldsymbol{\theta}^{\text{\tiny
HMC}}\big{)}=\boldsymbol{\theta}^{\text{\tiny HMC}}$ (with
$\psi=\psi^{\prime}$), and so the resulting Markov chain fails to jump
efficiently between modes.

To overcome the above issue, we consider a further modification of the basic
sampler (Algorithm 2) to encourage sample space exploration in the early
stages of adaptive Warp-U (Algorithm 3). The modified version of Algorithm 2
is presented in Algorithm LABEL:alg:highAdaptSampler in the online Supplement
Material. This version of algorithm draws the inverse Warp-U transformation
$\mathcal{H}_{\psi^{\prime}}$ with probability

$\widetilde{\nu}(\psi^{\prime}|\boldsymbol{\theta}^{*})\propto\exp\big{\\{}\log\big{(}\nu(\psi^{\prime}|\boldsymbol{\theta}^{*})\big{)}/C^{(s)}\big{\\}},$

instead of $\nu(\psi^{\prime}|\boldsymbol{\theta}^{*})$ in Line 5 of Algorithm
2. The constant $C^{(s)}(\geq 1)$ forces the probability towards a uniform
distribution, which is similar to the simulated annealing strategy. We
investigate combining Algorithm 3 with Algorithm LABEL:alg:highAdaptSampler,
instead of Algorithm 2. As the adaptive stage $s$ increases in Algorithm 3,
the constant $C^{(s)}(\geq 1)$ is gradually decreased to $1$. In the final
stages, we set $C^{(s)}=1$ and Algorithm LABEL:alg:highAdaptSampler reduces to
the original version of basic Warp-U sampler (i.e. Algorithm 2 except with an
HMC proposal). See Section LABEL:appendix:numerical of the online
Supplementary Material for more details of Algorithm
LABEL:alg:highAdaptSampler .

To further assist local mode exploration in the adaptive stage, we adopted the
following auxiliary distribution

$\phi_{\rm
mix}=0.5\phi_{\text{mix}}^{\text{(adapt.)}}+0.5\phi_{\text{mix}}^{\text{(opt.)}}.$
(43)

In the above, $\phi_{\text{mix}}^{\text{(opt.)}}$ is a mixture distribution
with 11 mixture components. Its centres are found via optimisation:
specifically, we ran an optimisation algorithm with 2,000 random initial
points and 11 unique local maxima were detected. As before, each mixture
component of $\phi_{\text{mix}}^{\text{(opt.)}}$ has equal weight and diagonal
covariance matrix $3I$. The distribution $\phi_{\text{mix}}^{\text{(adapt)}}$
is fitted adaptively to MCMC samples as Algorithm 3 proceeds, while
$\phi_{\text{mix}}^{\text{(opt.)}}$ is fixed. This choice of $\phi_{\rm mix}$
and its mixture weights is arbitrary and not optimised, but is sufficient to
illustrate the benefits of this approach in our empirical studies, which we
now describe.

Table 5: Settings for the bridge estimator. The auxiliary
$\phi_{\text{mix}}^{\text{(i.i.d.)}}$ is fitted by i.i.d. samples by EM
algorithm, $\phi_{\text{mix}}^{\text{(opt.)}}$ is a mixture of Gaussian whose
centres are local modes detected by optimisation. The auxiliary
$\phi_{\text{mix}}^{\text{(adapt)}}$ is fitted to the adaptive Warp-U samples.

Setting | $\phi_{\rm mix}$ | Samples for bridge estimators
---|---|---
i.i.d. | $\phi_{\text{mix}}^{\text{(i.i.d.)}}$ | i.i.d samples
HMC+WarpU(i.i.d.) | $0.5\phi_{\text{mix}}^{\text{(i.i.d.)}}+0.5\phi_{\text{mix}}^{\text{(opt.)}}$ | Warp-U MCMC samples with HMC proposal
HMC+WarpU(adaptive) | $0.5\phi_{\text{mix}}^{\text{(adapt)}}+0.5\phi_{\text{mix}}^{\text{(opt.)}}$ | Warp-U MCMC samples with HMC proposal

Figure 13: Warp-U bridge and stochastic Warp-U bridge estimation of the
normalising constant under 100 simulation replicates. The true value of the
normalising constant is $1$. From left to right, the boxplots are divided into
three groups. See Table 5 for the specific setting of each group.

For the 50 dimensional target distribution, we compare the Warp-U bridge and
stochastic Warp bridge estimation under three specific settings based on the
discussion above. The settings are summarized in Table 5. In the i.i.d.
setting, the auxiliary distribution $\phi_{\text{mix}}^{\text{(i.i.d.)}}$ is
obtained by applying the EM algorithm to 30,000 i.i.d. samples; meanwhile, an
independent set of 10,000 i.i.d. samples are fed to the bridge estimators. In
the setting HMC+WarpU(i.i.d.), the auxiliary distribution is chosen as
$0.5\phi_{\text{mix}}^{\text{(i.i.d.)}}+0.5\phi_{\text{mix}}^{\text{(opt.)}}$,
where $\phi_{\text{mix}}^{\text{(i.i.d.)}}$ is the Gaussian mixture fitted to
i.i.d. samples, and $\phi_{\text{mix}}^{\text{(opt.)}}$ is the Gaussian
mixture whose centres are found by optimisation. In the setting
HMC+WarpU(adaptive), we use
$0.5\phi_{\text{mix}}^{\text{(adapt.)}}+0.5\phi_{\text{mix}}^{\text{(opt.)}}$
in (43) as the auxiliary distribution. For both HMC+WarpU(i.i.d) and
HMC+WarpU(adaptive), we run Algorithm LABEL:alg:highAdaptSampler with
$C^{(s)}=1$ to generate 10,000 MCMC samples during the last adaptive stage.
These MCMC samples are then supplied to the bridge estimators.

For HMC+WarpU(i.i.d) and HMC+WarpU(adaptive), the samples generated by one run
of the modified algorithm are shown in Figure LABEL:fig:dim50:iidPhiMix and
Figure LABEL:fig:dim50:adaptivePhiMix in the online Supplement Material,
respectively. With a properly fitted $\phi_{\rm mix}$, the Warp-U sampler is
capable of exploring multiple modes in this complex high dimensional setting.

The estimation of the normalising constants under 100 replicates is shown in
Figure 13. The true value of the normalising constant is 1. The leftmost
boxplot indicates the Warp-U bridge estimator with i.i.d sample provides a
very accurate estimation, which is consistent with the work of Wang et al.
(2022). The next boxplot shows the stochastic Warp-U bridge estimator with
i.i.d. samples has a slightly larger estimation error. This is because the
stochastic Warp-U bridge estimator requires that each transformed density
matches the standard Gaussian distribution, and this requirement is difficult
to meet for the current target distribution with complex density contours and
heavy tails. Under the HMC+WarpU(i.i.d) and HMC+WarpU(adaptive) settings, the
estimation of the normalising constant is close to 1 in most simulation
replicates. A few replicates produce estimates close to 100.

The computational findings presented in this subsection demonstrate that our
proposed adaptive sampler might struggle to perform efficiently when faced
with high-dimensional distributions with complex density contours. During the
initial adaptive phase, $\phi_{\rm mix}$ is initialized with random guessing,
which can result in Warp-U transformations failing to shift the Markov chain
to a different high-density region. To mitigate this issue, we combined the
concept of simulated annealing with the Warp-U transformations, which
facilitates sample space exploration in the early stages. Despite the tactics
applied within this subsection, the derived $\phi_{\rm mix}$ does not
precisely emulate the target distribution. A potential avenue for future
research is to develop other variational inference techniques (Arenz et al.,
2020; Yan et al., 2023) to fit a high-quality $\phi_{\rm mix}$ from the onset.
Additionally, the efficiency of local exploration around individual local
modes is compromised due to the inefficiency of Hamiltonian Monte Carlo (HMC)
when sampling densities with highly irregular density contours and heavy tails
(Hoffman et al., 2019; Livingstone et al., 2019). These limitations could be
the reason why the Warp-U samples used in the bridge estimators have failed to
deliver highly accurate estimates of the true normalizing constant.

##  5 Exoplanet Detection Using Radial Velocity (RV) Data

In astronomy, one of the most successful approaches for detecting exoplanets
is the radial velocity (RV) method (exoplanets are planets outside our Solar
System). The radial velocity of a star is its velocity towards or away from
the Earth in meters per second (m/s). When an exoplanet orbits a star, the
gravitational force of the planet impacts the RV of the star, and RV data can
therefore be used to detect exoplanets. Consider a candidate model
$\mathcal{M}$ for capturing the physical system and noise, e.g., $\mathcal{M}$
might be a Keplerian model for a single exoplanet orbiting a star with
Gaussian measurement noise. From a Bayesian perspective, it is natural to
compute the Bayesian evidence of the model $\mathcal{M}$:

$\mathcal{Z}\equiv p(\boldsymbol{d}|\mathcal{M})=\int
p(\boldsymbol{d}|\boldsymbol{\theta},\mathcal{M})p(\boldsymbol{\theta}|\mathcal{M})\mu(\mathrm{d}\boldsymbol{\theta}),$
(44)

where $\boldsymbol{d}$ is the RV data and $\boldsymbol{\theta}$ denotes the
parameters of $\mathcal{M}$. Here
$p(\boldsymbol{d}|\boldsymbol{\theta},\mathcal{M})$ is the likelihood function
and $p(\boldsymbol{\theta}|\mathcal{M})$ as the prior. The Bayesian evidence
$\mathcal{Z}$ is the normalising constant of the posterior distribution of
$\boldsymbol{d}$ given up to proportionality by $\propto
p(\boldsymbol{d}|\boldsymbol{\theta},\mathcal{M})p(\boldsymbol{\theta}|\mathcal{M})$,
and represents the evidence in support of the model $\mathcal{M}$. Consider
the following model for the observed RV $v_{i}$ at time $t_{i}$:

$v_{i}=v_{\text{pred}}(t_{i}|\boldsymbol{\theta})+\epsilon_{i},$ (45)

where $v_{\text{pred}}$ is a Keplerian model for the planetary system (see
Danby, 1988; Loredo et al., 2012), $\boldsymbol{\theta}$ denotes the physical
parameters, and $\epsilon_{i}$ is a noise term. Further details of the
physical model are given in Section LABEL:appendix:physical_model of the
online Supplementary Material. One important statistical feature of RV data is
that the noise exhibits correlation across observations. In accordance with
Rajpaul et al. (2015) and Jones et al. (2022), we make the assumption that
this noise follows a Gaussian process, i.e.,
$\boldsymbol{\epsilon}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma})$. We
model the covariance matrix $\boldsymbol{\Sigma}$ by

$\boldsymbol{\Sigma}_{ij}=\kappa_{ij}+\delta_{ij}(\sigma_{i}^{2}+\sigma_{L}^{2}),$
(46)

where $\kappa_{ij}$ is a quasi-periodic kernel, $\delta_{ij}$ is the Kronecker
delta, $\sigma_{i}^{2}$ is the variance due to measurement error, and
$\sigma_{L}^{2}$ captures additional variation. The quasi-periodic kernel is
given by

$\kappa_{ij}=\alpha^{2}\exp\left[-\frac{1}{2}\left\\{\frac{\sin[\pi(t_{i}-t_{j})/\tau]}{\lambda_{p}^{2}}+\frac{t_{i}-t_{j}}{\lambda_{e}^{2}}\right\\}\right],$
(47)

and the kernel hyperparameters treated as known in our dataset (described
below) are fixed at $\alpha=\sqrt{3}\text{ meters}/\text{second}$,
$\lambda_{e}=50.0$ days, $\lambda_{p}=0.5$ (unit-less), and $\tau=20$ (days).

Our dataset consists of $n=200$ simulated observations from the Extremely
Precise Radial Velocities (EPRV3) Evidence Challenge which compared Bayesian
evidence estimates produced by different statistical methods in the context of
RV exoplanet detection; see Nelson et al. (2020) for details. Each observation
is accompanied by the time at which the RV is measured and a known measurement
error (i.e., standard deviation). The data are plotted in the left panel of
Figure 14.

Figure 14: The left panel shows the true radial velocity signal as a function
of time (solid red line), and observed measurements and their associated
measurement errors $\sigma_{i}^{2}$ (blue triangles and vertical bars). Note
that the measurement errors represent errors typically reported by an
observing telescope, and as can be seen in (46) and the left panel they do not
account for all uncertainties. The right panel shows the marginal posterior
distribution of the mean anomaly parameter. It compares the estimated
densities using the samples obtained by Warp-U MCMC (dash-dot blue line), PT
(dotted orange line) and HMC (dashed green line). The solid red line is the
estimated target density by the numerical integral.

###  5.1 Comparison with Parallel Tempering and Hamiltonian MC

We compare the performance of our adaptive Warp-U MCMC sampler to that of PT
and Hamiltonian Monte Carlo (HMC) (see Betancourt and Girolami, 2015; Neal et
al., 2011), in terms of its ability to recover the marginal posterior
distribution of the Mean Anomaly physical parameter used in the Keplerian
model introduced in (45) and Section LABEL:appendix:physical_model of the
online Supplementary Material. The Hamiltonian Monte Carlo sampling is
performed using the RStan package; see Stan Development Team (2019) and Stan
Development Team (2020). We focus on the Mean Anomaly because its posterior
distribution is expected to be multimodal. The temperature grid needed for PT
was tuned via the algorithm proposed in Atchadé et al. (2011). To obtain a
baseline for comparison, we also applied brute-force numerical integration (at
a high computational cost) to obtain the true posterior distribution of the
model parameters. The right panel of Figure 14 compares the marginal posterior
distribution of the Mean Anomaly obtained under the three different sampling
methods. The figure demonstrates that our adaptive Warp-U MCMC sampler more
accurately recovers the marginal posterior distribution than PT and HMC. It is
also worth noting that PT required substantial tuning, whereas our Warp-U MCMC
method did not (the Stan package automatically ran tuning required for HMC).

###  5.2 Estimation of the Bayesian Evidence

Next, we compare the performance of stochastic Warp-U bridge estimation with
standard bridge sampling estimation and Warp-U bridge estimation. We first
applied PT and our adaptive Warp-U MCMC sampler to get the target samples
needed by the estimation methods. To make the comparison fair, the total
number of evaluations of the target density across the sampling and estimation
steps is set to be the same for all the methods.

Table 6: RMSE (and associated SE) when estimating the $\log_{10}$ Bayesian evidence for a planet using bridge sampling, Warp-U bridge estimation, and stochastic Warp-U bridge estimation. | Bridge | Warp-U | S. Warp-U
---|---|---|---
Sampling | RMSE | SE | RMSE | SE | RMSE | SE
PT | $0.238$ | $0.015$ | $0.239$ | $0.017$ | $0.115$ | $0.008$
Warp-U | $0.136$ | $0.004$ | $0.137$ | $0.008$ | $0.059$ | $0.002$

Following Nelson et al. (2020), we use the median value of the Bayesian
evidence obtained across all the methods investigated in the Extremely Precise
Radial Velocities (EPRV3) Evidence Challenge as the quasi true value (i.e.,
$\log_{10}(\hat{c})=-193.71$). Table 6 compares the root mean square error
(RMSE) for bridge sampling, Warp-U bridge estimation, and stochastic Warp-U
bridge estimation. The RMSE is the smallest under our stochastic Warp-U bridge
estimation method. Table 6 also shows that the RMSE is lower when the target
samples are obtained by our adaptive Warp-U MCMC sampler than when they are
obtained via PT, which again illustrates the benefits of our sampling method.
Our closest estimate to the quasi-true value was $\log_{10}(\hat{c})=-193.74$,
and was obtained by stochastic Warp-U bridge estimation using the samples from
Warp-U MCMC sampler. Based on the bias shown in Nelson et al. (2020), which
ignores the variance and therefore represents the methods in Nelson et al.
(2020) favourably, our stochastic Warp-U bridge estimator utilising samples
from the Warp-U MCMC sampler has an RMSE that is comparable to the best
methods investigated in the EPRV3 Evidence Challenge. Lastly, note that the
standard errors of the RMSE in Table 6 are smaller under our stochastic Warp-U
bridge estimation strategy.

##  6 From Past to the Future

###  6.1 A brief Overview of Comparable Methods

Many algorithms have been developed to sample from multi-modal densities, a
number of which simultaneously perform sampling and estimation of normalising
constants. A leading example of the latter category of techniques is the
Generalised Wang-Landau (GWL) algorithm proposed by Liang (2005), which is an
energy based adaptive importance sampling method. The multi-stage approach
used in our adaptive method was inspired by the GWL algorithm, and earlier
adaptive importance sampling strategies such as Liang (2002), Berg and Neuhaus
(1991), and Wang and Landau (2001). There have been several extensions to the
GWL algorithm, including Liang et al. (2007) and Bornn et al. (2013), but also
some concerns about its convergence properties. Jacob et al. (2014) showed
that only some variations reach the so-called flat histogram convergence
criterion in finite time, whereas other variations do not. Furthermore, Wang
et al. (2022) illustrated that the GWL normalising constant estimator is
sometimes inefficient, and the alternative strategy of applying Warp-U bridge
estimation to the GWL draws (after weighted resampling) can substantially
reduce RMSE (for fixed computational resources).

Indeed, although it is conceptually appealing to combine sampling and
estimation in a single step, existing techniques for performing these tasks
separately are in some ways more developed. Some existing algorithms also
apply the idea of transporting the mass of the target density, e.g., Parno and
Marzouk (2018) constructs transport maps to match the target distribution and
a reference distribution for more efficient Metropolis-Hasting proposal. Some
theory regarding the optimal transformation has been developed by Villani
(2003). Pompe et al. (2020) proposed a sampling method that begins by finding
the modes of the target distribution via optimization; then, based on this
knowledge of the mode locations, they augment the parameter space with a mode