Sampling from Rough Energy Landscapes

TL;DR

This paper investigates the impact of multiscale roughness on sampling algorithms for Boltzmann distributions, showing that certain methods are sensitive to roughness while others are robust, and proposing new strategies to improve sampling performance.

Contribution

It demonstrates the limitations of the Metropolis Adjusted Langevin Algorithm with rough landscapes and introduces alternative methods that resist fine-scale roughness effects.

Findings

Metropolis Adjusted Langevin Algorithm performance deteriorates with increased roughness

Random Walk Metropolis remains insensitive to roughness

Proposed new strategies outperform Random Walk Metropolis in rough landscapes

Abstract

We examine challenges to sampling from Boltzmann distributions associated with multiscale energy landscapes. The multiscale features, or "roughness," corresponds to highly oscillatory, but bounded, perturbations of a smooth landscape. Through a combination of numerical experiments and analysis we demonstrate that the performance of Metropolis Adjusted Langevin Algorithm can be severely attenuated as the roughness increases. In contrast, we prove that Random Walk Metropolis is insensitive to such roughness. We also formulate two alternative sampling strategies that incorporate large scale features of the energy landscape, while resisting the impact of fine scale roughness; these also outperform Random Walk Metropolis. Numerical experiments on these landscapes are presented that confirm our predictions. Open questions and numerical challenges are also highlighted.

Figures

Tables

Table 1. Table 1: Ratios of the MSD at optimal σ 𝜎 \sigma for Modified MALA and Independence sampler to the RWM for the rough harmonic landscape, ( 4.1 ). As ϵ → 0 → italic-ϵ 0 \epsilon\to 0 , the amplification in performance saturates, but increases with dimension. The Independence sampler method appears to give better performance in moderate to high dimensions.

$ϵ$	$n=1$	$n=10$	$n=20$	$n=50$	$n=100$
$2^{-5}$	2.43	6.73	9.59	15.71	4.38
$2^{-6}$	2.43	6.74	9.61	15.88	15.16
$2^{-7}$	2.43	6.74	9.62	15.57	20.93
$2^{-8}$	2.43	6.73	9.60	15.40	25.07
$2^{-9}$	2.43	6.73	9.62	15.72	20.36
$2^{-10}$	2.43	6.74	9.61	15.48	22.51

Table 2. (a) Modified MALA Sampler

$ϵ$	$n=1$	$n=10$	$n=20$	$n=50$	$n=100$
$2^{-5}$	2.43	6.73	9.59	15.71	4.38
$2^{-6}$	2.43	6.74	9.61	15.88	15.16
$2^{-7}$	2.43	6.74	9.62	15.57	20.93
$2^{-8}$	2.43	6.73	9.60	15.40	25.07
$2^{-9}$	2.43	6.73	9.62	15.72	20.36
$2^{-10}$	2.43	6.74	9.61	15.48	22.51

Table 3. (b) Independence Sampler

$ϵ$	$n=1$	$n=10$	$n=20$	$n=50$	$n=100$
$2^{-5}$	2.20	10.15	17.91	40.36	12.59
$2^{-6}$	2.20	10.17	17.94	39.56	57.39
$2^{-7}$	2.20	10.17	17.88	40.07	70.85
$2^{-8}$	2.20	10.17	17.91	39.26	67.93
$2^{-9}$	2.20	10.17	17.94	39.98	67.61
$2^{-10}$	2.20	10.17	17.92	39.12	73.55

Table 4. Table 2: Ratios of the MSD at optimal σ 𝜎 \sigma for Modified MALA and Independence sampler to the RWM for the rough double well landscape, ( 4.3 ). As ϵ → 0 → italic-ϵ 0 \epsilon\to 0 , the amplification in performance saturates, but increases with dimension. The Independence sampler method appears to give better performance in moderate to high dimensions. In contrast to the harmonic potential, both the Modified MALA sampler and the Independence sampler show signs of decreased performance in high dimensions.

$ϵ$	$n=1$	$n=10$	$n=20$	$n=50$	$n=100$
$2^{-5}$	0.31	4.90	6.93	3.84	1.22
$2^{-6}$	0.30	4.89	6.98	11.45	1.80
$2^{-7}$	0.29	4.89	6.96	5.22	1.01
$2^{-8}$	0.29	4.89	6.97	11.69	15.40
$2^{-9}$	0.29	4.89	7.03	1.00	1.00
$2^{-10}$	0.29	4.90	6.94	11.64	17.30

Table 5. (a) Modified MALA Sampler

$ϵ$	$n=1$	$n=10$	$n=20$	$n=50$	$n=100$
$2^{-5}$	0.31	4.90	6.93	3.84	1.22
$2^{-6}$	0.30	4.89	6.98	11.45	1.80
$2^{-7}$	0.29	4.89	6.96	5.22	1.01
$2^{-8}$	0.29	4.89	6.97	11.69	15.40
$2^{-9}$	0.29	4.89	7.03	1.00	1.00
$2^{-10}$	0.29	4.90	6.94	11.64	17.30

Table 6. (b) Independence Sampler

$ϵ$	$n=1$	$n=10$	$n=20$	$n=50$	$n=100$
$2^{-5}$	6.30	330.89	589.73	441.69	62.95
$2^{-6}$	6.31	330.06	594.21	1317.64	262.12
$2^{-7}$	6.31	330.42	592.44	434.52	0.00
$2^{-8}$	6.32	330.32	594.40	1343.84	2150.46
$2^{-9}$	6.32	330.32	598.89	0.00	0.00
$2^{-10}$	6.31	330.81	592.33	1335.96	2287.34

Table 7. Table 3: Ratios of the MSD at optimal σ 𝜎 \sigma for Modified MALA and Independence Sampler to the RWM for the rough double well landscape, ( 4.3 ) with x 0 = ( 0 , 0 , … , 0 ) T ∈ ℝ n subscript 𝑥 0 superscript 0 0 … 0 𝑇 superscript ℝ 𝑛 x_{0}=(0,0,\ldots,0)^{T}\in\mathbb{R}^{n} . Compare with Table  2 .

$ϵ$	$n=1$	$n=10$	$n=20$	$n=50$	$n=100$
$2^{-5}$	0.31	4.91	6.95	3.88	1.22
$2^{-6}$	0.30	4.89	6.99	11.49	1.73
$2^{-7}$	0.29	4.89	6.98	11.54	7.48
$2^{-8}$	0.29	4.89	6.98	11.43	16.84
$2^{-9}$	0.29	4.88	6.99	11.67	17.11
$2^{-10}$	0.29	4.88	6.96	11.60	15.67

Table 8. (a) Modified MALA Sampler

$ϵ$	$n=1$	$n=10$	$n=20$	$n=50$	$n=100$
$2^{-5}$	0.31	4.91	6.95	3.88	1.22
$2^{-6}$	0.30	4.89	6.99	11.49	1.73
$2^{-7}$	0.29	4.89	6.98	11.54	7.48
$2^{-8}$	0.29	4.89	6.98	11.43	16.84
$2^{-9}$	0.29	4.88	6.99	11.67	17.11
$2^{-10}$	0.29	4.88	6.96	11.60	15.67

Table 9. (b) Independence Sampler

$ϵ$	$n=1$	$n=10$	$n=20$	$n=50$	$n=100$
$2^{-5}$	6.31	331.09	591.19	439.82	72.42
$2^{-6}$	6.31	330.56	594.66	1339.65	191.19
$2^{-7}$	6.32	330.24	594.54	1337.39	976.85
$2^{-8}$	6.32	330.13	594.19	1311.59	2562.82
$2^{-9}$	6.32	329.43	595.04	1355.55	2566.64
$2^{-10}$	6.31	329.40	593.38	1325.89	2199.09

Equations

$$\mu (dx)=Z^{-1}{e}^{-\beta V(x)}dx,Z=\int e^{-\beta V(x)}dx,V:{\mathbb{R}}^n\to \mathbb{R},$$

$$V_{ϵ}(x)=V_0(x)+V_1(x,x/ϵ),ϵ>0.$$

$$\text{RWM: }{X}_{k+1}^{\mathrm{p}}$$

$$\text{MALA: }{X}_{k+1}^{\mathrm{p}}$$

$$\mathrm{Gap}(\mathcal{T})\equiv 1-\underset{\lambda \in \sigma (\mathcal{T})\setminus \{1\}}{\sup }\lambda \equiv 1-\Lambda .$$

$$V_n(x)=\sum _{i=1}^{n}v(x_i),x=(x_1,\dots ,x_n)\in {\mathbb{R}}^n.$$

$$\mu (dx)\propto \prod _{i=1}^n{e}^{-v(x_i)}d{x}_i.$$

$$V_{ϵ,n}(x)=\sum _{i=1}^n\underset{\equiv v_{ϵ}(x_i)}{\underbrace{v_0(x_i)+v_1(x_i,x_i/ϵ)}}.$$

$${\mathrm{MSD}}_n={\mathbb{E}}_{\mu }[|X_{k+1}-X_k{|}^2].$$

$${\mathrm{MSD}}_n/n={\ell }^2{n}^{-I}a(\ell ;v)+o(n^{-I}),I>0,$$

$${\sigma }^2={\ell }^2{n}^{-I}$$

$${\ell }_{\star }=\underset{\ell }{\mathrm{argmax}}{\ell }^{2}a(\ell ;v).$$

$$a(\ell ;v)=2\Phi \left(-\frac{{\ell }^{1/I}}{2}\sqrt{\mathcal{K}[v]}\right),$$

$$V_{ϵ}(x)=V_0(x_1)+V_1(x_1,x_2/ϵ),x=(x_1,x_2)\in {\mathbb{R}}^{n=n_1+n_2}.$$

$$V_{ϵ}(x)=V(x_1/ϵ,x_2,\dots ,x_n),x\in {\mathbb{R}}^n.$$

$$\underset{k\to \infty }{\lim }\frac{1}{k}\log \Vert {\rho }_0{T}^k-\mu {\Vert }_{\mathrm{TV}}$$

$$\underset{k\to \infty }{\lim }\frac{1}{k}{\mathrm{Var}}_{\mu }\left(\sum _{j=0}^{k-1}f(X_j)\right)$$

$$\underset{f\in L_{0,1}^2(\mu )}{\inf }\frac{1}{2}{\mathbb{E}}_{\mu }[|f(X_{k+1})-f(X_k){|}^2]=1-\underset{\lambda \in \sigma (\mathcal{T})\setminus \{1\}}{\sup }\lambda =\mathrm{Gap}(\mathcal{T}),$$

$$\mathrm{Gap}(\mathcal{T})=\underset{f\in L^2(\mu )}{\inf }\frac{{\mathbb{E}}_{\mu }[|f(X_{k+1})-f(X_k){|}^2]}{2{\mathrm{Var}}_{\mu }(f(X))}$$

$$\underset{ϵ\to 0}{\mathrm{liminf}}\mathrm{Gap}({\mathcal{T}}_{ϵ})>0,$$

$$\underset{ϵ\to 0}{\mathrm{liminf}}\frac{{\mathrm{MSD}}_n}{2{\mathrm{Var}}_{\mu }(X)}>0.$$

$$\sigma$$

$$X_{k+1}^{\mathrm{p}}=X_k-\frac{{\sigma }^2}{2}\nabla V_0(X_k)+\sigma {\xi }_{k+1}.$$

$$X_{k+1}^{\mathrm{p}}=Y_{k+1}\sim e^{-V_0(x)}dx.$$

$$R(x,y)=V(x)-V(y)+\log \frac{q(y\to x)}{q(x\to y)}.$$

$$F(r)=1\wedge e^r,$$

$$F(r)=(1+e^{-r}{)}^{-1}.$$

$$\mathrm{osc}{V}_1=\underset{x}{\sup }{V}_1(x,x/ϵ)-\underset{x}{\inf }{V}_1(x,x/ϵ)<\infty .$$

$$\begin{array}{rl}{e}^{-2\mathrm{osc}{V}_1}{\mathbb{E}}_{{\mu }_0}[|f(X_{k+1})-f(X_k){|}^2]& \le {\mathbb{E}}_{\mu }[|f(X_{k+1})-f(X_k){|}^2]\\ & \le e^{2\mathrm{osc}{V}_1}{\mathbb{E}}_{{\mu }_0}[|f(X_{k+1})-f(X_k){|}^2].\end{array}$$

$$e^{-3\mathrm{osc}{V}_1}\mathrm{Gap}({\mathcal{T}}_0)\le \mathrm{Gap}({\mathcal{T}}_{ϵ})\le e^{3\mathrm{osc}{V}_1}\mathrm{Gap}({\mathcal{T}}_0)$$

$$e^{-\mathrm{osc}{V}_1}{\mu }_0(dx)\le \mu (dx)\le e^{\mathrm{osc}{V}_1}{\mu }_0(dx)$$

$$Z^{-1}{e}^{-V(x)}\ge Z^{-1}{e}^{-V_0(x)}{e}^{-\sup V_1(x,x/ϵ)}$$

Full text

{AMS}

65C05, 65C40, 60J22

1 Introduction

In this work, we consider the task of sampling from a Boltzmann distribution,

[TABLE]

when $V$ is, in some sense “rough,” or “rugged.” We are particularly interested in multiscale landscapes of the form

[TABLE]

Here, $V_{0}$ is a smooth, long range, trapping potential ($V_{0}(x)\to\infty$ as $|x|\to\infty$) that is bounded from below, and $V_{1}$ is bounded with local short wavelength features. $V_{1}(x,y)$ will also be assumed to be smooth. An example of such a rough landscape, and its impact on the associated distribution, is shown in Figure 1. Model potentials like (1.2) serve as prototypes for rough and multiscale landscapes found in disordered media and soft matter, [1, 2, 3, 4, 5]. The goal of the present work is to assess how such roughness impacts the performance of well known Markov Chain Monte Carlo (MCMC) sampling strategies like Random Walk Metropolis (RWM) and Metropolis Adjusted Langevin (MALA).

We recall that RWM and MALA generate samples for $e^{-\beta V(x)}$ with proposals111Proposals will be denoted with a superscript $\mathrm{p}$.

[TABLE]

These proposals are then accepted or rejected with the appropriate rule to ensure detailed balance with respect to $\mu(dx)\propto e^{-\beta V(x)}dx$.

As an example, sample MALA paths with $\sigma=1$ for the landscapes in Figure 1 are shown in Figure 2. A path for the smooth landscape exhibits better mixing than the one for the rough landscape. On the rough landscape the trajectory stagnates.222Stagnation corresponds to persistent rejection of proposals. This begs the question of whether or not $\sigma=1$ was merely a poorly chosen value – perhaps with a different value the rough landscape would also be efficiently sampled. Large values of $\sigma$ result in proposals with greater magnitude, but few will be accepted, and the trajectory will move slowly. Conversely, small values of $\sigma$ produce more readily accepted proposals, but their size limits exploration of the state space. Consequently, an optimal choice of $\sigma$ is anticipated for each distribution.

Assuming we tune our sampler to the optimal $\sigma$ for each $V_{\epsilon}$ we seek to assess how $\epsilon$ impacts sampling performance. In this work, optimality, at a fixed value of $\epsilon$ and at a fixed dimension, will refer to maximization of some measurement of mixing, discussed below over the set of numerical parameters, i.e. $\sigma$ in the case of MALA. Ultimately, our work indicates that even at the optimal value of proposal variance, MALA will cease to be effective as $\epsilon\to 0$. In contrast, even a poorly tuned RWM sampler remains robust in the $\epsilon\to 0$ limit.

In a “global” sense, robustness refers to the stability of the mixing and asymptotic variance properties of the chain as $\epsilon$ vanishes. This can be quantified through the spectral gap of the transition operator $\mathcal{T}$ of a given MCMC method333Some authors refer to (1.5) as the interval, $\operatorname{Int}(\mathcal{T})$. In [6], $\operatorname{Gap}(\mathcal{T})$ is instead defined as $1-\sup|\lambda|$.

[TABLE]

The gap controls both the mixing of the process and the time averaged variance constant, [6]. In particular, if the gap vanishes ($\Lambda\to 1$), then the mixing collapses and the variance bound explodes. The relationship between the gap and these quantities is given below.

Experimentally accessible measurements of mixing can be found by looking at observables. In particular, we consider the mean square displacement of the chain in stationarity. If this remains positive in the $\epsilon\to 0$ limit, it provides a “local” (in the sense of a single observable) notion of robustness. In addition to being straightforward to estimate through simulations, the mean squared displacement provides an upper bound on the spectral gap.

To obtain better performance than RWM, we also formulate two related sampling strategies that incorporate information about the large scale (i.e. long wavelength) features of the energy landscape through $V_{0}$ in (1.2). Indeed, our results, particularly Theorem 1 and Corollary 2 show that for potentials that can be decomposed as in (1.2) with $V_{0}$ smooth and trapping and $V_{1}$ rough but bounded, if the proposal of the sampling strategy is $\epsilon$-independent, then the performance of the method will also be $\epsilon$-independent.

1.1 Review of prior work

The question of optimizing $\sigma$ to maximize performance was initially examined in [7, 8] and has been subsequently studied in other works, including [9, 10, 11, 12, 13, 14, 15, 16, 17]. Many of these works consider an energy landscape of the type

[TABLE]

This choice of the potentials induces Boltzmann distributions that are products (for brevity, we take $\beta=1$):

[TABLE]

Thus, coordinates only interact through the accept/reject step of the method. Some results for potentials other than (1.6) have also been obtained. In [9], the authors treat distributions which have a density $e^{-\Phi(x)}$ with respect to a product measure (1.7). In earlier works, [7, 8, 9], it was often assumed that whichever sampler is studied, the process is in stationarity. This has been relaxed in the more recent results, [11, 12, 13, 16, 17, 10].

Some of our results will also specialize to the rough analog of (1.6),

[TABLE]

As in the case of (1.2), we assume $v_{0}$ is trapping while $v_{1}$ is uniformly bounded. Both $v_{0}(x_{i})$ and $v_{1}(x_{i},y_{i})$ are assumed to be smooth functions.

In the case of (1.6), in stationarity, the performance can be measured by considering the means square displacement (MSD)

[TABLE]

The strategy of [9] is to use the MSD as a proxy for mixing, and $\sigma$ is selected to maximize it in the limit of $n\to\infty$. Indeed, this leads to the results that, as $n\to\infty$,

[TABLE]

where $I$ and $a$ both depend on the method, but not on the dimension $n$. The choice of $\sigma$ is then related to $\ell$ by

[TABLE]

The function $a$ is the mean acceptance rate in the $n\to\infty$ limit, [18].

Consequently, we can maximize this measure of performance as $n\to\infty$ by solving

[TABLE]

The optimal $\sigma_{\star}$ is inferred from (1.11), and there is an associated optimal acceptance rate, $a(\ell_{\star})$. This optimization provides a strategy for tuning the value of $\sigma$ to achieve the optimal acceptance rate approximately 23% for RWM. Analogously, one tunes MALA to have a 57% acceptance rate, [7, 8], and Hamiltonian Monte Carlo (HMC) to have a 65% acceptance rate, [15].

For RWM $I=1$ while it is $I=1/3$ for MALA. The function $a$ has the explicit form

[TABLE]

with $\mathcal{K}$ a functional involving derivatives of $v$; $\Phi$ is the standard normal ($N(0,1)$) cumulative distribution function. A similar result holds for HMC, [15].

There are caveats to applying these results in practical computations, [19]. In particular, the results are obtained as $n\to\infty$ for distributions of form (1.7), and often assume the process to be in stationarity. Many distributions of interest will not be of this form, so target acceptance rates, like 23% for RWM, may be inappropriate. However, in the recent work [20], it was demonstrated that for a more general class of distributions than (1.7), 23% remains the optimal acceptance rate for RWM.

We mention these results because they motivate certain aspects of this work, such as the examination of product measures and the examination of $\operatorname{MSD}_{n}$ as a proxy for mixing. However, we emphasize that this work is focused on problems at fixed $n$, letting $\epsilon\to 0$. The preceding results can provide guidance in the $\epsilon\to 0$ limit, but they would necessitate first taking $n\to\infty$.

We also highlight the recent work in [10] which studies “ridged” densities associated with potentials of the form

[TABLE]

Here, roughness is only present in a subset of the degrees of freedom ($x_{2}$). Examining RWM for such a problem, the authors are able to derive a limiting diffusion from which they can find an optimal step size. This limiting diffusion has a state dependent diffusion coefficient. Both the drift and diffusion coefficients are nontrivial requiring averaging against the rough degrees of freedom.

Related to the work on ridged densities, and the present work, is [21]. In this work, the authors consider the case that one of the coordinates is scaled differently than the others. This would correspond to a potential like

[TABLE]

In [21] the authors also look for algorithms that are less sensitive to length scale variations in the gradients. They obtain results on MALA and HMC showing poor behavior in the $\epsilon\to 0$ limit. An important tool that they use in their analysis is the Dirichlet form and its relationship to the spectral gap; we, too, make use of that approach.

Another relevant work is [22]. There, the authors sought to perform gradient based sampling on non-differentiable energy landscapes and proposed using a Moreau-Yosida regularization. This approach is related to one of the mechanisms that we propose in order to overcome roughness, though our potentials are smooth, but highly oscillatory.

1.2 Measures of performance and notions of robustness.

As mentioned, the key metrics that we use to assess performance are the spectral gap, (1.5), along with the MSD, (1.9). We recall the relationships amongst these quantities, [6, 21]. First, the relaxation to the equilibrium in the total variation (TV) and the time averaged variance constant (TAVC) are controlled by $\Lambda$ and $\operatorname{Gap}(\mathcal{T})$:

[TABLE]

for an initial distribution with $d\rho_{0}/d\mu\in L^{2}(\mu)$ and any function $f\in L^{2}(\mu)$. Consequently, if $\operatorname{Gap}(\mathcal{T})\to 0$, mixing ceases and there is no upper bound on the TAVC. On the other hand, if the gap remains positive then there is a priori bound on the TAVC. Due to the inequality in (1.16) a positive spectral gap does not imply a lower bound on mixing.

A key relationship is between spectral gap and the Dirichlet form

[TABLE]

where $L_{0,1}^{2}(\mu)$ is the subset of mean zero, unit variance functions in $L^{2}(\mu)$. An elementary computation reveals that this is equivalent to

[TABLE]

Consequently, for any square integrable function, if we can estimate the term on the right-hand side of (1.19) , we have obtained an upper bound on the spectral gap. This allows us to use (1.9) as a proxy for the gap; if we find, numerically, that $\operatorname{MSD}_{n}\to 0$ as $\epsilon\to 0$, that is strong empirical evidence that $\operatorname{Gap}(\mathcal{T}_{\epsilon})\to 0$ too. Again, due to the infimum in the Dirichlet form, positivity of the one step jumping distance of any observable, including $\operatorname{MSD}_{n}$, does not imply positivity of the gap.

We thus formalize two notions of robustness. Our robustness criteria for a method with a transition operator $\mathcal{T}_{\epsilon}$ for sampling the Boltzmann distribution of potential $V_{\epsilon}$ are

[TABLE]

Analogous robustness conditions can be constructed for other observables. Failure to be locally robust immediately implies the method cannot be globally robust. Likewise, a globally robust method automatically implies local robustness. However, there may be distributions and methods for which, on a particular observable, the method is locally robust, but for which it fails to be globally robust.

In (1.20) and (1.21), no mention is made of the choice of the numerical parameter $\sigma$. Some methods may fail to be robust if $\sigma$ is not chosen carefully. Thus, we introduce two conditional notions of robustness that depend on $\sigma$

[TABLE]

These notions generalize to methods with additional numerical parameters. Obviously, if a method fails to satisfy (1.22), it will fail to satisfy (1.20). Conversely, if the method satisfies (1.20), then it will also satisfy (1.22). Analogous relationships can be formulated for local robustness.

1.3 Results on performance in the presence of roughness.

Our key observations and results are:

1. (i)

At a fixed dimension $n$, the performance of RWM for potentials of the form (1.2), is globally robust in the sense of (1.20). This is a consequence of Corollary 2. Indeed, any method that uses $\epsilon$-independent proposals will similarly satisfy (1.20). Proposals with a sufficiently mild $\epsilon$-dependence will also be robust; see Corollary 9. The methods need not be optimally tuned for (1.20) to hold. 2. (ii)

A rigorous result is established that, subject to certain assumptions, MALA fails to be globally robust when $\sigma$ is inadequately scaled for potentials like (1.8). This is the content of Theorem 11. Specifically, when $\sigma$ is too large relative to $\epsilon$, the spectral gap will close. 3. (iii)

Numerical experiments and an explicit example indicate that for $n=1$, the optimal scaling of MALA is $\sigma\propto\sqrt{\epsilon}$, so that $\operatorname{MSD}_{1}\propto\epsilon$. In contrast, for $n$ sufficiently large, the empirical optimal scaling is $\sigma\propto\epsilon$, so that $\operatorname{MSD}_{n}/n\propto\epsilon^{2}$.

The experiments also indicate that MALA is not locally robust even at an empirically determined optimal $\sigma$, (1.23). Thus, the method suffers generically in the $\epsilon\to 0$. 4. (iv)

We formulate two alternative methods for potentials of type (1.2) that use large scale information contained in $V_{0}$. The first method, which we call Modified MALA, uses the proposals (at $\beta=1$)

[TABLE]

This method fits in the class of sampling methods studied in [15] where the authors rigorously established that for the path-wise accuracy of Metropolized integrators it is sufficient to accurately simulate the diffusion term, which (1.24) does.

The second method, we call an Independence Sampler, uses proposals (at $\beta=1$)

[TABLE]

The $Y_{k}$ samples are assumed to be independent, generated by an auxiliary process. Both of these methods are also insensitive to the roughness and outperform RWM.

As both of these methods use $\epsilon$-independent proposals Corollary 2 allows us to conclude they will also be globally robust in the $\epsilon\to 0$ limit. 5. (v)

For potentials that do not admit an obvious decomposition like (1.2) we propose using local entropy approximation, [23, 24] to extract the large scale information needed for either the Modified MALA method or the Independence Sampler.

In Section 2 we identify a bound on the performance of RWM and other algorithms, and we consider the asymptotic behavior of MALA as $\epsilon\to 0$. In Section 3 we present alternative methods that are also robust to $\epsilon\to 0$. Numerical experiments are presented in Section 4, and we conclude with a discussion in Section 5.

Acknowledgements

The work of PP was partially supported by DARPA project W911NF-15-2-0122 and GS was supported by US National Science Foundation Grant DMS-1818716. The authors thank UCLA IPAM for hosting them during the beginning of this project. The authors also thank M. Luskin, S. Osher, J. Mattingly, and N. Bou-Rabee for helpful discussions.

2 Bounds on performance.

In this section we present bounds on performance with respect to the robustness criteria. For any MCMC method, let $q(x\to y)$ denote the associated proposal kernel, and define

[TABLE]

Consequently, the proposal $X_{n}^{\mathrm{p}}$ is accepted with probability $F(R(X_{n},X_{n}^{\mathrm{p}}))$. The two forms of $F$ that we consider here are

[TABLE]

2.1 Roughness independent bounds.

For potentials of type (1.2) one can obtain $\epsilon$-independent upper and lower bounds on a variety of quantities. Indeed, by the boundedness assumption in the introduction, we are assured that

[TABLE]

Our main results in this subsection are the following:

Theorem 1**.**

Let $V_{\epsilon}$ be a potential of type (1.2) and assume $V_{1}$ has uniform in $\epsilon$ bounded oscillation in the sense of (2.3). If the sampling strategy uses proposals $q(x\to y)$ that are $\epsilon$-independent, then for any $f$, the one step jumping distance is bounded by $\epsilon$ independent constants:

[TABLE]

In (2.4), $\mu_{0}(dx)\propto e^{-V_{0}(x)}dx$.

A corollary to this result provides $\epsilon$ independent bounds on the spectral gap:

Corollary 2**.**

Under the same assumptions as Theorem 1

[TABLE]

where $\mathcal{T}_{0}$ is the transition operator of the method with proposals generated by $q$ sampling $\mu_{0}$.

Remark 3**.**

We emphasize that these results, and the results in the rest of Section 2.1, are all at a fixed dimension $n$. Constants, like $e^{\pm 2\operatorname{osc}V_{1}}$ in (2.4), may depend in an unfavorable way on $n$. Additionally, Theorem 1 and Corollary 2 hold independently of the choice of any numerical parameters, like $\sigma$ in RWM.

To prove Theorem 1 and its corollary, we first prove the following bounds on the distribution.

Lemma 4**.**

Let $V_{\epsilon}$ be a potential of type (1.2), and assume $V_{1}$ has uniform in $\epsilon$ bounded oscillation in the sense of (2.3). Then

[TABLE]

Proof 2.1**.**

The proof of this follows from direct estimates on the densities. First,

[TABLE]

while

[TABLE]

As an immediate consequence we have bounds for the mean.

Corollary 5**.**

Let $V_{\epsilon}$ satisfy the same assumptions as in Lemma 4. Then for any non-negative observable, $f$, which may depend on a small parameter $\epsilon$,

[TABLE]

Similarly we obtain bounds on the variance.

Lemma 6**.**

Let $V_{\epsilon}$ satisfy the same assumptions as in Lemma 4. Then $L^{2}(\mu)$ and $L^{2}(\mu_{0})$ are equivalent as sets, and

[TABLE]

Proof 2.2**.**

For $f\in L^{2}(\mu)$, Lemma 4 ensures

[TABLE]

Hence, $f\in L^{2}(\mu_{0})$ too and $L^{2}(\mu)\subset L^{2}(\mu_{0})$. A similar computation shows the reverse inclusion. Next, for $f\in L^{2}(\mu)$, let $f_{c}(x)=f(x)-\mathbb{E}_{\mu_{0}}[f(X)]$. Then, using Lemma 4 again,

[TABLE]

Reversing the roles of $\mu$ and $\mu_{0}$ establishes the analogous lower bound in (2.8).

Remark 7**.**

For potentials satisfying the assumptions of the preceding results, it will often be sufficient to examine $\mathbb{E}_{\mu}[|f(X_{k+1})-f(X_{k})|^{2}]$ since a prior bound on the variance is provided by Lemma 6.

Proof 2.3**.**

Proof of Theorem 1.

Given potential $V=V_{0}+V_{1}$ of the form (1.2) satisfying the assumptions:

[TABLE]

Since $\operatorname{osc}V_{1}$ is bounded,

[TABLE]

Consequently, for either choice of (2.2),

[TABLE]

Since the proposal, $q(x\to y)$, is also assumed $\epsilon$-independent, then, for any $X_{0}$,

[TABLE]

In the preceding upper and lower bounds on $\mathbb{E}[|X_{1}-X_{0}|^{2}]$, no $\epsilon$ is present. If we now apply Lemma 4 to the upper and lower bounds in (2.12), we obtain (2.4).

Corollary 2 follows from the theorem and Lemma 6.

Proof 2.4**.**

Proof of Corollary 2.

For any non constant $f$,

[TABLE]

Taking the infinum over $L^{2}(\mu)$,

[TABLE]

This last equality is due to $L^{2}(\mu)$ and $L^{2}(\mu_{0})$ being equivalent as sets. An analogous computation establishes the lower bound.

Theorem 1 and Corollary 2 immediately apply to RWM, as it has a proposal independent of $\epsilon$. Consequently, for RWM,

[TABLE]

In contrast, as MALA proposals include gradients of the potential, for potentials of the form (1.2), these results will not apply.

The results can be modified to allow for proposals that have some $\epsilon$ dependence:

Theorem 8**.**

Let $V_{\epsilon}$ be a potential of the type (1.2), and assume $V_{1}$ has bounded oscillation, uniformly in $\epsilon$, in the sense of (2.3). Assume the sampling strategy’s proposal kernel, $q(x\to y)$, satisfies the inequality

[TABLE]

where $q_{0}(x\to y)$ is an $\epsilon$-independent proposal kernel. Then the performance, as measured by $\operatorname{MSD}_{n}$, is $\epsilon$-independent:

[TABLE]

Proof 2.5**.**

As in the proof of Theorem 1, we begin by writing

[TABLE]

Analogously,

[TABLE]

Consequently, for a function $f$,

[TABLE]

Using Lemma 4 yields the result.

An analog of Corollary 2 can then be established, which we present without proof:

Corollary 9**.**

Under the same assumptions as Theorem 8,

[TABLE]

A method where Theorem 8 would apply is the tamed (or truncated) MALA, [25, 26, 27]. Given a parameter $\delta>0$, one form of the tamed MALA (at $\beta=1$) is

[TABLE]

This has the effect of ensuring that the gradient term is never larger than $1/\delta$, mitigating stiff regimes of the state space.

Corollary 10**.**

The performance, as measured by the spectral gap, of the tamed MALA on potentials of the type (1.2), is insensitive to $\epsilon$.

Proof 2.6**.**

For the tamed MALA as implemented in (2.16), let

[TABLE]

Then its proposal density is

[TABLE]

Since $|f_{\delta}(x)|\leq\delta^{-1}$,

[TABLE]

Setting

[TABLE]

we have that the tamed MALA proposal kernel satisfies

[TABLE]

Thus, Corollary 9 applies.

2.2 Optimal Proposal Variance for MALA.

Next, we perform a formal calculation on the performance of MALA with respect to proposal variance $\sigma$ and roughness $\epsilon$. Though this does not lead to any immediate conclusions, it may guide future analysis.

We define $M(\sigma)$ to be the $\operatorname{MSD}_{n}$ obtained with $\sigma$, assuming stationarity, $x\sim\mu$, with the Barker rule, (2.2b). Thus, the optimal value of $\sigma$ can be obtained by solving $M^{\prime}(\sigma)=0$. Proceeding with this strategy we have

[TABLE]

where $g$ is the Gaussian density of the normal distribution $N(x-\frac{\sigma^{2}}{2}\nabla V(x),\sigma^{2}I)$ and $\mu(dx)$ is the Boltzmann distribution.

Proposition 2.7**.**

At a critical point of (2.21), where $M^{\prime}(\sigma)=0$,

[TABLE]

or

[TABLE]

While this is a formal calculation, it does not require any particular assumptions on the potential beyond the associated distribution being normalizable and having the necessary moments.

Proof 2.8**.**

See Appendix A for the derivation of (2.22). The expression (2.23) follows by solving for $\sigma^{2}$.

2.3 Insights from the high dimensional limit.

We recall from the introduction that for product distributions such as (1.7), as $n\to\infty$, the $\operatorname{MSD}_{n}$ per degree of freedom is given by (1.10). Here, we consider the $n\to\infty$ limit at fixed $\epsilon$, and then examine the impact of the roughness through $\epsilon$. This corresponds to a consideration of the order $\lim_{\epsilon\to 0}\lim_{n\to\infty}$. Elsewhere in this manuscript, we focus on fixed $n$ while letting $\epsilon\to 0$.

The optimal $\ell$ for (1.10) can be found using (1.12) and (1.13). From (1.11) we can then find the optimal choice of $\sigma$ and obtain the scaling of the $\operatorname{MSD}_{n}/n$ in (1.10). After changing coordinates to $\lambda=-\ell^{1/I}\sqrt{\mathcal{K}[v]}/2$, and optimizing in this transformed coordinate, one can deduce that the optimal $\ell^{2}\sim\mathcal{K}^{-I}$. Consequently, for potentials like (1.8), the optimal scaling, as $n\to\infty$, is $\sigma^{2}\sim\left({n\mathcal{K}}\right)^{-I}$ and

• (i)

For RWM, $I=1$ and $\mathcal{K}[v_{\epsilon}]=\mathbb{E}_{\mu}[(v_{\epsilon}^{\prime})^{2}]$, so that $\sigma^{2}\sim\epsilon^{2}/n$;

• (ii)

For MALA, $I=1/3$ and $\mathcal{K}[v_{\epsilon}]=\mathbb{E}_{\mu}[5(v_{\epsilon}^{\prime\prime\prime})^{2}+3(v_{\epsilon}^{\prime\prime})^{3}]/48\sim\epsilon^{-6}$ so that $\sigma^{2}\sim\epsilon^{2}/n^{1/3}$.

In this limit, both RWM and MALA will fail to be locally robust at optimal $\sigma$ (in the sense of (1.23)) as $\epsilon\to 0$. Thus, they will also suffer when a suboptimal value is selected for any observable. This does not contradict the analysis of Section 2.1. There, in the product case, our lower bound would have the prefactor $e^{-n\operatorname{osc}v_{1}}$, which vanishes more rapidly than the $\epsilon^{2}/n$ we have here.

2.4 Scaling in dimension one.

We consider the question of sampling from the distribution associated with $V(x)=\tfrac{1}{2}\epsilon^{-1}x^{2}$ in $n=1$, using MALA. While this may seem odd, as this corresponds to a Gaussian distribution, it reveals a distinct scaling in low dimensions. It can be interpreted as the necessary scaling to efficiently sample individual modes of the highly multi-modal landscape in Figure 1(b), or alternatively, as the analog of the stiff differential equation problem $\dot{x}=-\epsilon^{-1}x$, $0<\epsilon\ll 1$.

In this case,

[TABLE]

and the term in the accept/reject rule is

[TABLE]

With the aid of a computer algebra system (see, also, Appendix B) we conclude

[TABLE]

A consequence of this computation is that it provides an explicit example for which MALA will not satisfy the local robustness criterion, even with optimal $\sigma$, (1.23).

The function $m(\delta)$ has a single maximum (see Figure 12), and the maximum value is achieved at $\delta_{\star}=1.27797$, and $m(\delta_{\star})=1.8494$. Thus, the optimal scaling for the time step is

[TABLE]

This scaling is a different from the scaling found in Section 2.2. We note that this optimal choice is smaller than the Euler-Maruyama stability threshold which requires $\sigma\leq 2\sqrt{\epsilon}$. Additionally, the acceptance rate at this optimal value is $0.70$, which is quite different than the $n\to\infty$ acceptance rate of $0.57$.

2.5 Poorly scaled proposals.

We show rigorously that, under particular assumptions, if $\sigma$ is improperly scaled, MALA will fail to be robust in the sense of (1.20). Our approach is based on the method of [21] (see, also, [6, 28]) using the relationship between the Dirichlet form and the spectral gap, (1.5).

As in [21], the interval $\Lambda$ is bounded in terms of the conductance of a set, $K$, with $\mu(K)\in(0,1/2)$. Using the observable

[TABLE]

we immediately have for a method with the transition operator $\mathcal{T}$

[TABLE]

Proposition 2.9**.**

Assume that

The potential is $V(x)=V_{0}(x)+V_{1}(x)$, where $V_{0}$ is trapping and $V_{1}(x)=V_{1}(x,x/\epsilon)$ is bounded. 2. 2.

A concentration inequality holds for $\mu$

[TABLE]

where $C_{0}$, $\kappa$, and $p$ are independent of $\epsilon$. 3. 3.

There exists $m>0$ such that for all $\epsilon$ sufficiently small, the set $K=\{x\,|\,|x|>m\}$ satisfies

[TABLE]

where $\underline{\mu}_{K}$ and $\overline{\mu}_{K}$ are independent of $\epsilon$. 4. 4.

The gradient $\nabla V_{0}$ satisfies a growth bound with an exponent $\delta>0$, such that for all $x$

[TABLE] 5. 5.

Let $\alpha\geq 0$, $\theta\in\mathbb{R}$, and $\gamma>0$ satisfy

[TABLE]

Let $\sigma=\epsilon^{\alpha}$, and define the sets

[TABLE]

Then it holds

[TABLE]

In (2.35), the relation $\lesssim$ is introduced. Generically, if $a\lesssim b$, then there is a constant, $C$, independent of $b$, such that $a\leq Cb$.

Before proving the result, a few remarks are in order.

1. (i)

The concentration inequality for $\mu$ will hold provided $\mu_{0}$, the distribution associated with $V_{0}$, has a concentration inequality. This is a consequence of Corollary 5. 2. (ii)

The bound (2.30) will similarly hold provided such a condition holds for $V_{0}$. Again, this is a consequence of Corollary 5 by selecting a set $K$ for which $0<\mu_{0}(K)<\tfrac{1}{2}e^{-\operatorname{osc}V_{1}}$. 3. (iii)

The set $L_{\epsilon,\theta}$ captures the points in space where $\nabla V_{1}$ will be large. To obtain a rate, it is necessary to know the measure of $L_{\epsilon,\theta}^{c}$, requiring additional details on the structure of the potential. In particular, it is essential to estimate

[TABLE]

which is clearly related to the conductance of set $K$ with the additional constraints that $X_{k}$ and $X_{k+1}$ both reside in $L_{\epsilon,\theta}^{c}$, the set where the gradient is “small”. 4. (iv)

To see how the set (2.36) can become small, consider the case that $V_{0}=x^{2}/2$ and $V_{1}=\cos(x/\epsilon)$. Then $\nabla V_{1}(x)=-\sin(x/\epsilon)/\epsilon$. For any $\theta>0$ and all $\epsilon>0$ sufficiently small

[TABLE]

Then, recalling that $g(y;x,\sigma)$ is the Gaussian density of the proposal

[TABLE]

For any $x\in K\cap L_{\epsilon,\theta}^{c}$

[TABLE]

Additionally,

[TABLE]

Combining (2.37) and (2.38), (2.36) is thus bounded as

[TABLE]

For $\theta>1/2$ and $\gamma$ and $\alpha$ sufficiently small the above expression is $\ll\epsilon$. 5. (v)

As the result requires $\alpha<1/2$, $\sigma\gg\sqrt{\epsilon}$ in this regime.

Proof 2.10**.**

Proof of Proposition 2.9.

Since $\mu(K)$ is, by assumptions, uniformly bounded away from [math] and $1/2$ as $\epsilon\to 0$, it is sufficient to study $\mathbb{P}_{\mu}(X_{k+1}\in K^{c},X_{k}\in K)$. For all sufficiently small $\epsilon$ we can write

[TABLE]

where we have used the assumed concentration inequality at the end.

Next, we divide the set depending on whether $X_{k}$ and $X_{k+1}$ lie in the high or low gradient sets

[TABLE]

We consider the first of the three terms on the right hand side of (2.39):

[TABLE]

Conditioned on $X_{k}\in K_{\epsilon,\gamma}\cap L_{\epsilon,\theta}$,

[TABLE]

therefore

[TABLE]

By our assumptions on the exponents, for all $\epsilon$ small enough,

[TABLE]

where we have a tail inequality on $N(0,I_{n})$.

For the second term in (2.39) we use reversibility to get

[TABLE]

Using the same approach as above we have

[TABLE]

Again, by our assumptions on the exponents and a tail inequality on the Gaussian

[TABLE]

Combining these estimates we obtain (2.35).

The next proposition provides an estimate of the last term in (2.35) for a particular class of potentials.

Theorem 11**.**

Assume that

There exists $p\geq 1$ and constants $V_{\min}$ and $C_{p}>0$

[TABLE] 2. 2.

The growth bound of Proposition 2.9 holds for $\nabla V_{0}$. 3. 3.

The potential $V_{1}$ takes the form

[TABLE] 4. 4.

$\alpha$, $\gamma$, $\delta$, and $\theta$ satisfy (2.32).

Then there exists a constant $\kappa$ such that

[TABLE]

Before proving this result we make a few remarks:

1. (i)

If $V_{0}$ is harmonic it will satisfy these assumptions with $p=2$ and $\delta=1$. 2. (ii)

As long as $\alpha<2\theta$, the performance will deteriorate as $\epsilon\to 0$. By (2.32), this necessitates $\alpha<2/5$. 3. (iii)

In the numerical experiments discussed below, $\operatorname{MSD}_{1}$ was found to be maximized when $\sigma\sim\sqrt{\epsilon}$, while in higher dimensions, $\operatorname{MSD}_{n}$ was maximized when $\sigma\sim\epsilon$. Thus, there is a “gap” between values of $\alpha$ for which Theorem 11 predicts poor performance, and values for which we numerically observe peak performance. We interpret this as a gap in the analysis because, even at the optimal choice of $\sigma$, our results show $\operatorname{MSD}_{n}\to 0$ as $\epsilon\to 0$, implying the method is not locally (or globally) robust in the sense of (1.23). 4. (iv)

As $n$ increases, the gap vanishes more and more rapidly as $\epsilon\to 0$.

Proof 2.11**.**

Proof of Theorem 11.

A concentration inequality will obviously hold for $\mu$, as such an inequality holds for $\mu_{0}$, with $\mathbb{P}_{\mu_{0}}(|X|>t)\lesssim e^{-\frac{1}{2}\kappa t^{p}}$; this is a consequence of (2.40). The constant $\kappa$ is related to $C_{p}$, $p$, and the dimension.

Since the concentration inequality holds, and $V_{0}$ is trapping while $V_{1}$ is bounded, Corollary 5 ensures that a set $K=\{|x|>m\}$ exists and satisfies (2.30). Thus, the assumptions of Proposition 2.9 hold.

Since

[TABLE]

it will be sufficient to estimate $\mu(K\cap L_{\epsilon,\theta}^{c})$ and the integral $\int_{K^{c}\cap L_{\epsilon,\theta}^{c}}g(y;x,\sigma)dy$ for $x\in K\cap L_{\epsilon,\theta}^{c}$.

We first bound the set $L_{\epsilon,\theta}^{c}$. For all $\epsilon$ small enough

[TABLE]

Next, we bound the set $K$ by

[TABLE]

Consequently, for all $\epsilon$ small enough,

[TABLE]

Similarly

[TABLE]

and

[TABLE]

Now, using (2.40) and (2.43),

[TABLE]

The first term in the above sum is as in the example (2.38)

[TABLE]

The other terms are treated the same way

[TABLE]

Consequently

[TABLE]

For $x\in L_{\epsilon,\theta}^{c}\cap K$, using (2.44),

[TABLE]

Combining this with (2.45) yields the result.

3 Alternative algorithms for rough landscapes.

Having concluded that there are impediments to MALA for sampling on rough landscapes, but desiring a method that may perform better than RWM, we propose two methods here. Both these methods require a smoothed version of the potential which may not be accessible. Indeed, as the potential of interest is unlikely to admit a simple decomposition like (1.2), we also propose a method for approximating a smoothed landscape.

3.1 Modified MALA.

The first method, which we call Modified MALA, involves generating proposals from an auxiliary potential, $U$. The algorithm generates the proposals

[TABLE]

and $R(x,y)$ is modified to be

[TABLE]

If the scale separation in (1.2) holds, we might take $U=V_{0}$ to get (1.24). Since the proposal is $\epsilon$-independent, Corollary 2 applies to this method.

3.2 Independence sampler.

The second method we propose is to generate auxiliary samples, $Y_{k}$, from another distribution $\propto e^{-U(x)}dx$, and perform independence sampling against $\mu(dx)\propto e^{-V(x)}dx$. More explicitly,

[TABLE]

This approach is agnostic as to how samples from the auxiliary distributions are generated – any strategy which produces (approximately) independent samples will be satisfactory. When the assumptions about (1.2) hold, then we would take $U=V_{0}$ giving (1.25). Corollary 2 again applies to this method. Given the choice of $U$ and the mechanism for sampling from its distribution, this algorithm has no free parameters – there is no $\sigma$ that can be tuned.

In our numerical examples below, where the potentials take the form (1.8), we numerically approximate an inverse function sampler on each degree of freedom, generating i.i.d. samples $\propto e^{-v_{0}(x_{i})}$ for $i=1,\ldots,n$. This approach, of course, will not be practical for general landscapes, and we return to this issue in the discussion.

3.3 Finding smoothed landscapes.

A question that remains is what to do when (1.2) fails to hold. One option is to use physical intuition about the problem to identify a potential $U(x)$ that has suitable properties. More systematically, we can use the local entropy approach formulated in [24, 23], or, equivalently the Moreau-Yosida approximation to estimate a smoothed version of $V(x)$.

Given $\gamma>0$, we define $V_{\gamma}$ as

[TABLE]

which corresponds to a Gaussian filter of the Boltzmann distribution. The associated gradient is then

[TABLE]

The density, $\rho_{\infty}(y\mid x)$, can be defined by first introducing

[TABLE]

so that

[TABLE]

The potential $V_{\gamma}$ could be estimated with the simple Monte Carlo scheme

[TABLE]

Likewise, (3.5) could be estimated by integrating the auxiliary diffusion

[TABLE]

If $V$ is assumed to grow at infinity, then for sufficiently small $\gamma$, $U_{\gamma}(y\mid x)$ will be convex about $y=x$, and (3.8) will converge to equilibrium exponentially fast. The gradient $\nabla V_{\gamma}$ can then be estimated by Monte Carlo

[TABLE]

where each $Y^{(j)}\sim\rho_{\infty}(\cdot\mid x)$. Of course, $Y^{(j)}$ must now be appropriately sampled. Running short i.i.d. samples of (3.8) using, for instance, MALA, the algorithm introduces three additional numerical parameters: $N_{\mathrm{s}}$, the number of samples; $\delta t$, the fine scale time step; and $N_{\delta t}$, the number of fine scale time steps. We must also specify initial conditions.

With estimates of $V_{\gamma}$ and $\nabla V_{\gamma}$, we then use them as the smoothed landscapes in our coupled independence sampler scheme or Modified MALA scheme.

4 Numerical experiments.

In this section we demonstrate, numerically, how roughness can impede MALA, and how RWM, along with the methods discussed in Section 3, can resist the roughness. In these computational examples, we use the Metropolis rule, (2.2a). By Remark 7, it is sufficient to study $\operatorname{MSD}_{n}$ without including the variance factor in (1.21) or (1.23).

4.1 Rough harmonic potential.

As a first example, we consider a potential of the type (1.8), with

[TABLE]

The potential is depicted in Figure 4 at various values of $\epsilon$.

We explore this potential by varying both $\epsilon$ and the dimension $n$. A priori, we do not know the optimal value of $\sigma$ for each $\epsilon$ and $n$ pair, for each algorithm. Thus, we examine a range of $\sigma$ values, running long trajectories for each, and then empirically estimate the one with the maximum $\operatorname{MSD}_{n}$. See Figure 5 for examples of such a computations. In this way we are able to compare performance across methods. In these experiments inverse temperature is set to $\beta=5$ and $10^{8}$ iterations are performed in each run. The starting point for these runs is $x_{0}=(0,0,\ldots,0)^{T}\in\mathbb{R}^{n}$. For the independence sampler, we combine numerical quadrature and interpolation to approximate the inverse cumulative distribution function of the distribution $e^{-\beta v_{0}(x_{i})}$. This allows us to perform inverse function sampling. Consequently, there are no free parameters to tune for this strategy.

Having estimated the optimal $\sigma$, we can then compare performance between the algorithms, as a function of $\epsilon$ and $n$. This is shown in Figure 6. This indicates that $\sigma\propto\epsilon$ for MALA, provided $n$ is large enough. In this scaling, Theorem 11 does not apply. However, the numerical experiments reveal that mixing, as measured by MSD, is impeded because of the size of $\sigma$. We thus conclude from the numerical experiment that MALA fails to satisfy (1.23), implying that it will also fail to be globally robust, even with optimal $\sigma$. Also note that the mean acceptance rates across a range of $\epsilon$ and $n$ deviate from the idealized $n\to\infty$ RWM and MALA values.

Since RWM, Modified MALA, and the Independence sampler are all globally robust to $\epsilon\to 0$ in the sense of condition (1.20), we examined amplification in performance at different $\epsilon$ and $n$, as measured by

[TABLE]

This is shown in Table 1. The alternative methods always beat RWM, and, there is a greater improvement in higher dimension, though the performance improvement saturates as $\epsilon\to 0$. The independence sampler method typically outperforms the Modified MALA method.

4.2 Rough double well potential.

We repeat the experiment from Section 4.1 for a more challenging problem of sampling from the distribution given by the potential

[TABLE]

with parameters otherwise the same as in the harmonic case. The starting point for these runs is $x_{0}=(-1,-1,\ldots,-1)^{T}\in\mathbb{R}^{n}$. Results here are shown in Figure 7 and Table 2. These are similar to the experiments for the harmonic potential, but there are some serious differences and some cases where there is little or no performance gain.

First, in comparison to the harmonic problem, the performance amplification of the independence sampler is significantly larger in the double well problem in dimensions $n=10$ and $n=20$. We attribute this to the ability of the independence sampler to switch between the left and right super basins, a behavior that is entirely absent from the harmonic problem and one which RWM and the Modified MALA method are incapable of reproducing.

The second difference is that both the Modified MALA sampler and the independence sampler have significant performance issues in dimensions $n=50$ and $n=100$ values of $\epsilon$ including $2^{-7}$ and $2^{-9}$. Examining the actual trajectories, in the case of the Modified MALA, the optimal $\sigma$ is so small that the contribution of $\tfrac{\sigma^{2}}{2}\nabla V_{0}$ is negligible. Consequently, it is trending towards RWM resulting in no performance amplification. This is why the ratio is unity. For the independence sampler, the trajectory is stagnant. This is partially a function of the choice of initial conditions for the trajectories. If, we instead start the independence sampler at $x_{0}=(0,0,\ldots,0)^{T}\in\mathbb{R}^{n}$, instead of Table 2 (b), we obtain Table 3. This variation requires a more detailed discussion which we provide in Section 5. Despite these problems in the higher dimension, there is clearly improvement upon RWM in more modest dimension over a broad range of $\epsilon$.

In Figure 8, we repeat the experiment from Figure 2 with $n=1$, but for the Modified MALA and the Independence sampler schemes. Computing at $\sigma=1$, we see better mixing than in Figure 2, for the same rough energy landscape.

4.3 Results in dimension one.

We briefly consider the behavior in dimension one for the harmonic potential and the double well. For both energy landscapes, as shown in Figure 9, the optimal $\sigma\propto\sqrt{\epsilon}$, not $\sigma\propto\epsilon$. Thus, while the performance of MALA will also degrade in $n=1$, a different scaling appears here. This is consistent with the harmonic problem analyzed in Section 2.4.

4.4 Local entropy approximations.

We briefly consider the possibility of using local entropy, introduced in Section 3.3 with (3.4) and (3.5). This approach may be of use in problems where no straightforward scale separation, of the type found in (1.2), is present in the energy landscape. As motivation we consider the potential in dimension one

[TABLE]

where the $c_{j}$ and $k_{j}$ are random, from a particular distribution. Indeed, taking $M=10$ and $c_{j}\sim U(-.1,.1)$, $\log(k_{j})\sim U(10^{1},10^{3})$, we obtain the landscape shown in Figure 10, along with the leading contribution, $V_{0}$, and $V_{\gamma}$ obtained through a numerical quadrature at $\beta=5$ and $\gamma=0.05$. Clearly, the local entropy approximation eliminates the fine scale roughness found in the original potential.

In general, Monte Carlo approximations of $V_{\gamma}$ and $\nabla V_{\gamma}$ will be needed, as a quadrature will be impractical in high dimensions. An example of this is shown in Figure 11. In Figure 11(a), we compare $V_{0}$, $V_{\gamma}$ and a Monte Carlo estimate of $V_{\gamma}$ computed using (3.7) with $N_{\mathrm{s}}=10^{2}$.

In Figure 11(b), we compare $\nabla V_{0}$, $\nabla V_{\gamma}$, and a Monte Carlo estimate of $V_{\gamma}$, computed using (3.9). In this latter figure we take $N_{\mathrm{s}}=10^{4}$, and each sample is obtained by taking $N_{\rm\delta t}=4$ time steps with $\delta t=1$ in a variant of MALA that exactly linearly integrates the Ornstein-Uhlenbeck component of (3.8). Obviously, there are many options for how Monte Carlo estimates of $V_{\gamma}$ and $\nabla V_{\gamma}$ can be obtained. We will return to this in the discussion.

5 Discussion.

We have examined a class of rough energy landscapes where the performance of MALA can be driven to zero at a fixed dimension. When $\sigma$ is inadequately scaled with $\epsilon$, Proposition 2.9 and Theorem 11 reveal that MALA will fail to be globally robust, (1.20). Even if $\sigma$ is optimally scaled, according to the empirical estimates, the numerical simulations indicate that it also suffers.

There are several outstanding questions on MALA that merit investigation. First, it would be desirable to develop a rigorous understanding of why the optimal $\sigma$ scaling is $\sigma\propto\sqrt{\epsilon}$ in dimension one, while it is $\sigma\propto\epsilon$ for $n$ sufficiently large. Next, there is an analysis of why the spectral gap collapses at the optimal scaling as $\epsilon\to 0$; recall that our result in Theorem 11 does not apply to either scaling. Another, related, challenge is how to quantify the size of the small gradient set, (2.36), outside of the separable case, as $\epsilon\to 0$, or, alternatively, how to avoid analyzing the small set. Separately, there is the question of how $\sigma$ should be scaled for MALA before it has reached stationarity.

We have also demonstrated that RWM along with the Modified MALA and the Independence sampler are insensitive to the roughness of the landscape and are globally robust. While modified MALA and the Independence sampler require a smoothed energy landscape for proposals, RWM is a viable option without additional information. Indeed, Corollary 2 tells us that for rough energy landscapes like (1.2), with roughness bounded uniformly in $\epsilon$, if the proposal of the MCMC scheme is $\epsilon$-independent, then the performance will be insensitive to $\epsilon$. Corollary 9 indicates that other methods that have a weak $\epsilon$-dependence, such as the tamed MALA, will also exhibit robustness to the roughness.

As the results in Section 4 show, there still appear to be challenges in high dimensions. The degradation of the Independence sampler shown in Table 2 is a consequence of an unusually difficult starting point, $x_{0}=(-1,-1,\ldots,-1)^{T}\in\mathbb{R}^{n}$, as the following calculation shows.

For $\epsilon=2^{-7}$, $\cos(-1/\epsilon)=-0.69$ and at $\epsilon=2^{-9}$, $\cos(-1/\epsilon)=-1.00$. Consequently, every single component of the initial guess is close to a global minimum of the energy landscape and almost every proposal will be at a position in state space in which every coordinate has higher energy. Starting at any $X_{0}=x_{0}$,

[TABLE]

Under our assumptions on $U=V_{0}$, the fourth moment term is bounded by a constant independent of both $n$ and $\epsilon$. In the above expressions the expectation is over the proposal. The exponent in the acceptance probability is

[TABLE]

With $x_{1}=x_{2}=\ldots x_{n}=-1$, the $r_{i}$’s are thus i.i.d. Furthermore, since

[TABLE]

Suppose $\mu_{r}\equiv\mathbb{E}[r_{i}]<0$; $\mathbb{E}[R_{n}]=n\mu_{r}$. Then using Hoeffding’s inequality, [29],

[TABLE]

Consequently the mean acceptance and the sampling performance will both be driven to zero exponentially fast as $n\to\infty$ once $\mu_{r}<0$. Checking this numerically with $x_{1}=-1$, $\beta=5$, $\epsilon=2^{-9}$, and $v_{0}=\beta\cdot(x^{2}-1)^{2}$, $v_{1}=\beta\cdot\tfrac{1}{8}\cos(x)$,

[TABLE]

In fact, due to the oscillatory nature of the integral, $\mu_{r}\approx v_{1}(x_{1}/\epsilon)$. With such a value of $\mu_{r}$, it is straightforward to see, from (5.4), that the performance will rapidly be driven to zero for large enough $n$. In contrast, in Table 3, we see better results when $x_{1}=0$. Repeating the above computation, $\mu_{r}=0.625>0$, making it resistant to the previous pathology.

This phenomenon does not plague RWM because, in contrast to the independence sampler, RWM has a parameter, $\sigma$, that can be tuned to maintain an $\operatorname{O}(1)$ mean acceptance probability as $n\to\infty$. Indeed, the scaling $\sigma^{2}=\operatorname{O}(n^{-1})$ for RWM, discussed in Section 2.3, ensures this. We conjecture that a Metropolis within Gibbs sampler, by which only a subset of the $n$ coordinates are altered in each step of the sampler, will alleviate this problem in the Independence sampler.

Generating the samples from the smooth distribution for the Independence sampler is also a challenge in the general case. We believe this can be easily accomplished using MALA or HMC. These other samplers should be well behaved on the smooth energy landscape, allowing for the straightforward construction of approximately independent proposals for the rough landscape.

Smoothed energy landscapes might be available through a known decomposition like (1.2). We also conjecture that the optimal choice of approximate landscapes for potentials like (1.2) corresponds to the homogenized energy landscapes discussed in [3, 4, 5]. Unfortunately, computing such smoothed energies requires solving an elliptic PDE in a space of the same dimension as the considered state space; a Monte Carlo estimator of the solution may partially overcome this difficulty. Alternatively, physical knowledge of the system may motivate some choice for a surrogate smoothed landscape.

When these options are not available, the local entropy approximation is another possibility. The challenge to using local entropy, which we do not further develop here, is that unless the problem is in a very low dimension, auxiliary sampling algorithms must be formulated and tuned to first estimate $V_{\gamma}$ and $\nabla V_{\gamma}$. This task would involve determining a sample size, a sampling strategy, and some form of parallelization in order to outperform simpler alternatives like RWM.

Finally, the weakness of MALA in the presence of roughness can be seen as the MCMC manifestation of stiffness. We conjecture that it is a generic problem in gradient based MCMC methods, including HMC. Indeed, the magnitude of $\nabla V$ in HMC will constrain the time step of, for instance, the Verlet method used in the Hamiltonian flow subroutine. Thus, the number of force calls per HMC step will tend to increase with roughness degrading the overall performance. This has been partially addressed in [21], where the authors considered potentials of the form (1.15) and rigorously established that HMC suffers from scaling issues.

Appendix A Details of the mean square displacement computation

In this section, we give details of the derivation of (2.22). Differentiating (2.21) with respect to $\sigma$

[TABLE]

Consequently,

[TABLE]

Assuming that the optimal $\sigma$ occurs at a finite value, the first order condition $M^{\prime}=0$ will hold. The expression (A.5), at the optimal value, can then be expressed as

[TABLE]

This calculation makes use of the identity

[TABLE]

Since $M=\mathbb{E}[|x-y|^{2}F]$, we can re-write (A.6) to get (2.22).

Next, note that for the Barker proposal, $F(r)=(1+e^{-r})^{-1}=1-F(-r)$ so

[TABLE]

Additionally, recall that the Metropolis method satisfies detailed balance. Therefore,

[TABLE]

and we conclude

[TABLE]

Using this result in (A.6) gives us (2.22).

Appendix B Details of computations in dimension one

In this section, we provide a derivation of (2.26). We denote $g_{x}(x)$ the Gaussian density for $N(0,\epsilon)$ and $g_{y}(y|x)$ the Gaussian density for $N((1-\sigma^{2}\epsilon^{-1}/2)x,\sigma^{2})$. Then

[TABLE]

and we observe that

[TABLE]

Furthermore, we note that the set $R(x,y)>0$ corresponds to $|x|>|y|$ and $R(x,y)<0$ corresponds to $|x|<|y|$. We can thus use the symmetry $(x,y)\mapsto(y,x)$ to reduce (B.1) to

[TABLE]

Since the integrand, $2(y-x)^{2}g_{x}(x)g_{y}(y|x)$ is invariant to $(x,y)\mapsto(-x,-y)$, $I=IV$ and $II=III$. Thus we have

[TABLE]

Using Mathematica and making the change of variables, $\xi=x/\sqrt{\epsilon}$ and $\delta=\epsilon^{-1}\sigma^{2}/2$,

[TABLE]

Analogously,

[TABLE]

As in the case of the computation of $A_{1}$, we use the symmetry

[TABLE]

to reduce (B.6) to

[TABLE]

This is the split into four integrals, as in (B.3). Since the integrand, $2(y-x)^{2}g_{x}(x)g_{y}(y|x)$ is invariant to $(x,y)\mapsto(-x,-y)$, $I=IV$ and $II=III$. Thus

[TABLE]

Using Mathematica with $\xi=x/\sqrt{\epsilon}$ and $\delta=\epsilon^{-1}\sigma^{2}/2$, $\operatorname{MSD}_{n}=\epsilon m(\delta)$ we have

[TABLE]

This function is plotted in Figure 12.