2.1. Linear f: Metropolis–Hastings with logarithmic Hamiltonian and Catoni’s energy transformation algorithm

In this subsection, we specialize to $f(u) = u$ . It turns out we can understand the landscape modification as if the Hamiltonian is on a logarithmic scale whenever $\mathcal{H}(x) > c$ .

For $x,y \in \{\mathcal{H}(y) > \mathcal{H}(x) \geqslant c\}$ , since

\begin{align*}\int_{\mathcal{H}(x)}^{\mathcal{H}(y)} \frac{1}{u-c + \epsilon}\,du = \ln \bigg(\dfrac{\mathcal{H}(y)-c+\epsilon}{\mathcal{H}(x)-c+\epsilon}\bigg),\end{align*}

substituting the expression into (2.4) gives

\begin{align*} \exp\!\big(\!-\!\big(\mathcal{H}_{\epsilon,c}^f(y)-\mathcal{H}_{\epsilon,c}^f(x)\big)_+\big) &= \begin{cases} 1 &\mbox{if } \mathcal{H}(y) \leqslant \mathcal{H}(x), \\[5pt] \exp\!\bigg(\!-\!\dfrac{1}{\epsilon}(\mathcal{H}(y) - \mathcal{H}(x))\bigg) & \mbox{if } c\geqslant \mathcal{H}(y) > \mathcal{H}(x), \\[9pt] \exp\!\bigg(\!-\!\dfrac{1}{\epsilon}(c - \mathcal{H}(x))\bigg) \dfrac{\epsilon}{\mathcal{H}(y)-c+\epsilon} & \mbox{if } \mathcal{H}(y) > c \geqslant \mathcal{H}(x), \\[9pt] \dfrac{\mathcal{H}(x)-c+\epsilon}{\mathcal{H}(y)-c+\epsilon} & \mbox{if } \mathcal{H}(y) > \mathcal{H}(x) > c. \end{cases}\end{align*}

This resulting dynamics $M^f$ coincides with the energy transformation method introduced by Catoni [Reference Catoni1, Reference Catoni2] on $\{\mathcal{H}(y) > \mathcal{H}(x) > c\}$ , which is based on a logarithmic Hamiltonian. We refer readers to Section 2.6 for a more detailed account of the connection between landscape modification and Catoni’s energy transformation.

2.2. Quadratic f: Metropolis–Hastings with $\arctan$ Hamiltonian

In this subsection, we take $f(u) = u^2$ . In this case, the effect of the landscape modification gives an inverse-tangent-transformed Hamiltonian whenever $\mathcal{H}(x) > c$ .

For $x,y \in \{\mathcal{H}(y) > \mathcal{H}(x) \geqslant c\}$ , using the inverse-tangent difference formula we obtain

\begin{align*} \int_{\mathcal{H}(x)}^{\mathcal{H}(y)} \frac{1}{(u-c)^2 + \epsilon}\,du &= \sqrt{\frac{1}{\epsilon}}\left(\arctan\!\left(\sqrt{\frac{1}{\epsilon}}(\mathcal{H}(y)-c)\right) - \arctan\!\left(\sqrt{\frac{1}{\epsilon}}(\mathcal{H}(x)-c)\right)\right) \\[5pt] &= \sqrt{\frac{1}{\epsilon}} \arctan\!\left(\dfrac{\sqrt{\frac{1}{\epsilon}}(\mathcal{H}(y)-\mathcal{H}(x))}{1+\frac{1}{\epsilon}(\mathcal{H}(y)-c)(\mathcal{H}(x)-c)}\right),\end{align*}

and substituting the above expression into (2.4) gives

\begin{align*} & \exp\!\big(\!-\!\big(\mathcal{H}_{\epsilon,c}^f(y)-\mathcal{H}_{\epsilon,c}^f(x)\big)_+\big) \\[5pt] &= \begin{cases} 1 &\mbox{if } \mathcal{H}(y) \leqslant \mathcal{H}(x), \\[5pt] \exp\!\left(\!-\!\dfrac{1}{\epsilon}(\mathcal{H}(y) - \mathcal{H}(x))\right) &\mbox{if } c\geqslant \mathcal{H}(y) > \mathcal{H}(x), \\[9pt] \exp\!\left(\!-\!\dfrac{1}{\epsilon}(c - \mathcal{H}(x)) - \sqrt{\frac{1}{\epsilon}}\arctan\!\left(\sqrt{\frac{1}{\epsilon}}(\mathcal{H}(y)-c)\right)\right) &\mbox{if } \mathcal{H}(y) > c \geqslant \mathcal{H}(x), \\[9pt] \exp\!\left(\sqrt{\dfrac{1}{\epsilon}}\left(\arctan\!\left(\sqrt{\frac{1}{\epsilon}}(\mathcal{H}(x)-c)\right) - \arctan\!\left(\sqrt{\frac{1}{\epsilon}}(\mathcal{H}(y)-c)\right)\right)\right) &\mbox{if } \mathcal{H}(y) > \mathcal{H}(x) > c. \end{cases}\end{align*}

2.3. Square root f: Metropolis–Hastings with sum of square root and logarithmic Hamiltonian

In the final example, we let $f(z) = \sqrt{z}$ for $z \geqslant 0$ . For $x,y \in \{\mathcal{H}(y) > \mathcal{H}(x) \geqslant c\}$ , consider the integral

\begin{align*}\int_{\mathcal{H}(x)}^{\mathcal{H}(y)} \frac{1}{\sqrt{u-c} + \epsilon}\,du = 2\bigg(\Big(\sqrt{\mathcal{H}(y)-c} - \sqrt{\mathcal{H}(x)-c}\Big) - \epsilon \ln \!\bigg(\dfrac{\sqrt{\mathcal{H}(y)-c} + \epsilon}{\sqrt{\mathcal{H}(x)-c} + \epsilon}\bigg)\bigg).\end{align*}

Substituting the expression into (2.4) gives

\begin{align*} & \exp\!\big(\!-\!\big(\mathcal{H}_{\epsilon,c}^f(y)-\mathcal{H}_{\epsilon,c}^f(x)\big)_+\big)\\[5pt] &= \begin{cases} 1 &\mbox{if } \mathcal{H}(y) \leqslant \mathcal{H}(x), \\[5pt] \exp\!\bigg(\!-\!\dfrac{1}{\epsilon}(\mathcal{H}(y) - \mathcal{H}(x))\bigg) &\mbox{if } c\geqslant \mathcal{H}(y) > \mathcal{H}(x), \\[9pt] \exp\!\bigg(\!-\!\dfrac{1}{\epsilon}(c - \mathcal{H}(x)) - 2\sqrt{\mathcal{H}(y)-c}\bigg) \bigg(\dfrac{\sqrt{\mathcal{H}(y)-c} + \epsilon}{ \epsilon}\bigg)^{2\epsilon} &\mbox{if } \mathcal{H}(y) > c \geqslant \mathcal{H}(x), \\[9pt] \exp\!\left(2\sqrt{\mathcal{H}(x)-c} - 2\sqrt{\mathcal{H}(y)-c}\right) \bigg(\dfrac{\sqrt{\mathcal{H}(y)-c} + \epsilon}{\sqrt{\mathcal{H}(x)-c} + \epsilon}\bigg)^{2\epsilon} &\mbox{if } \mathcal{H}(y) > \mathcal{H}(x) > c. \end{cases}\end{align*}

Remark 1. (Use of landscape modification for sampling.) This paper focuses on investigating the acceleration effect of landscape modification in stochastic optimization and simulated annealing. However, the technique of landscape modification can also be applied to sampling from multimodal distributions. The paper Zhang and Choi [Reference Zhang and Choi7] analyses the use of landscape modification for sampling. We now briefly describe the setting therein. Suppose that we are interested in sampling from a multimodal distribution that we denote by $\nu(x) \propto e^{-\mathcal{H}(x)}$ . Applying the idea of landscape modification, we then construct an MH chain with the transformed Hamiltonian function

(2.5) \begin{align} \mathcal{H}^f_{1,c,\alpha}(x) = \int_{\mathcal{H}_{\textrm{min}}}^{\mathcal{H}(x)} \dfrac{1}{\alpha f((u-c)_+) + 1}\,du, \end{align}

so that its stationary distribution is given by $\nu^f_{1,c,\alpha}(x) \propto e^{-\mathcal{H}^f_{1,c,\alpha}(x)}$ . We note that the parameter $\alpha \geqslant 0$ , introduced in (2.5), controls the bias between the distribution $\nu$ and its landscape-modified counterpart $\nu^f_{1,c,\alpha}$ . If we anneal this parameter by sending $\alpha_t \to 0$ as $t \to \infty$ , then $\nu^f_{1,c,\alpha_t}$ converges weakly to the target distribution $\nu$ . As a result, we are able to construct a non-homogeneous MH chain that converges to $\nu$ in the long run while enjoying the benefits of landscape modification.

2.4. A Metropolized Ehrenfest urn with landscape modification

In this section, we discuss a Metropolized Ehrenfest urn model; our exposition closely follows that of Deuschel and Mazza ([Reference Deuschel and Mazza8], Section 5.1.2). Let us first briefly fix the setting. We consider the state space $\mathcal{X} = \{0,1,\ldots,d\}$ with $d \in \mathbb{N}$ and take a linear Hamiltonian $\mathcal{H}(x) = x$ , where $\mathcal{H}_{\textrm{min}} = 0$ . The proposal birth–death chain has generator Q given by $Q(x,x+1) = 1 - x/d$ , $Q(x,x-1) = x/d$ , $Q(x,x) = -1$ and zero otherwise, and we note that the stationary measure of Q is $\mu(x) \propto 2^{-d} (\substack{d \\[5pt] x})$ . With these choices, the classical Metropolized dynamics is $M^0_{\epsilon}(x,x+1) = Q(x,x+1) e^{-\frac{1}{\epsilon}}$ , $M^0_{\epsilon}(x,x-1) = Q(x,x-1)$ , and $\pi^0(x) \propto e^{-\frac{1}{\epsilon} x} \mu(x)$ . It is shown in Deuschel and Mazza ([Reference Deuschel and Mazza8], Equation (5.1.1)) that

(2.6) \begin{align} \lambda_2\big(\!-\!M^0_{\epsilon}\big) \leqslant \lambda_2(\!-\!Q) \dfrac{1}{\sum_{x=0}^d e^{-\frac{1}{\epsilon} x} \mu(x)} = \dfrac{d}{2} \dfrac{2^d}{\big(1+e^{-\frac{1}{\epsilon}}\big)^d}.\end{align}

At a fixed temperature $\epsilon$ , we note that the upper bound in (2.6) is exponential in d.

Now, we consider the landscape-modified MH with $f(z) = z$ and $c = 1$ . With these parameters, we compute

\begin{align*}\mathcal{H}(x) = \int_0^x \dfrac{1}{(u-1)_+ + \epsilon} \, du = \frac{1}{\epsilon} + \ln\!\bigg(\dfrac{(x-1)_+ + \epsilon}{\epsilon}\bigg).\end{align*}

Using Deuschel and Mazza ([Reference Deuschel and Mazza8], Equation (5.1.1)) leads to

(2.7) \begin{align} \lambda_2\big(\!-\!M^f_{\epsilon,c=1}\big) \leqslant \lambda_2(\!-\!Q) \dfrac{1}{\sum_{x=0}^d e^{-\mathcal{H}(x)} \mu(x)} = \dfrac{d}{2} \dfrac{2^d e^{\frac{1}{\epsilon}}}{\frac{\epsilon}{\lfloor d/2 \rfloor - 1 + \epsilon} \left(\substack{d \\[5pt] {\lfloor d/2 \rfloor}}\right)} \sim d^3 \frac{1}{\epsilon} e^{\frac{1}{\epsilon}},\end{align}

where we use Stirling’s formula, which gives, for large enough d,

\begin{align*} \left(\substack{d \\[5pt] {\lfloor d/2 \rfloor}}\right) &\sim \dfrac{2^d}{\sqrt{\pi \lfloor d/2 \rfloor}}.\end{align*}

As a result, the upper bound in (2.7) gives a polynomial dependency on d, at the tradeoff of an exponential dependency on $\frac{1}{\epsilon}$ . In retrospect this result is not surprising, since we are working with a logarithmic Hamiltonian and hence the asymptotics is in the polynomial of d instead of $2^d$ . This section serves as a warm-up example to prepare us for discussing the more complicated CW model in Section 3.

2.5. Tuning strategies for f and c

In this section, we discuss tuning strategies for f and the threshold parameter c in the context of using the MH chain with landscape modification for stochastic optimization.

First, the function f controls how the landscape is transformed above the threshold parameter c. For example, taking a linear $f(z) = z$ gives a logarithmic transformation, while taking a quadratic $f(z) = z^2$ yields an $\arctan$ transformation, as shown in Sections 2.1 and 2.2. In general, we recommend using $f(z) = z^2$ , since the $\arctan$ transformation is uniformly bounded above by $\pi/2$ , which facilitates exploration on the part of the landscape above c by giving a higher transition rate compared with the choice of $f(z) = z$ . However, we suspect that in numerical investigations, and depending on the target Hamiltonian $\mathcal{H}$ , there are possibilities of using a linear f to yield improved convergence towards the global minimum.

Second, as we shall see in Section 4, the threshold parameter c controls the clipped critical height $c^*$ , and ideally it should be set as close to $\mathcal{H}_{\textrm{min}}$ as possible. This is possible if we have information about the value $\mathcal{H}_{\textrm{min}}$ , which is the case for some statistical physics and theoretical computer science models [Reference Nardi and Zocca9, Reference Zocca10]. In general, however, we may not have access to the value $\mathcal{H}_{\textrm{min}}$ , and one general method is to tune the threshold parameter adaptively by setting the value at time t to be the running minimum up to time t generated by the chain. However, the resulting process becomes non-Markovian because of the adaptive tuning. Another adaptive tuning strategy is to set the value of c at time t to be $c_t = \mathcal{H}(y_t) - d$ , where $y_t$ is the proposed state at time t generated by the proposal chain and $d \geqslant 0$ is a fixed number. In this way, the transition rate of the landscape-modified MH chain is always greater than or equal to that of the MH chain without landscape modification. We shall numerically investigate this tuning strategy and report positive results in Section 4.1.

As far as this paper is concerned, we assume a fixed c in all of the main theoretical results. We shall postpone to future work a systematic numerical study of investigating various choices of f and tuning strategies for c on benchmark functions, as well as a theoretical analysis of adaptively tuning c by the running minimum generated by the algorithm in an MH chain with landscape modification.

2.6. Connections between landscape modification and other acceleration techniques

In this section, we outline conceptual similarities and differences between MH with landscape modification and other common acceleration techniques in the literature for MH and simulated annealing.

2.6.1. Catoni’s energy transformation algorithm

Let $\alpha_1 \geqslant 0$ , $\alpha_2 \geqslant 0$ , $\alpha_3 > - \mathcal{H}_{\textrm{min}}$ be three parameters. In Catoni [Reference Catoni1, Reference Catoni2], the author introduces the energy transformation algorithm by transforming the Hamiltonian $\mathcal{H}$ to

\begin{align*}F_{\alpha_1,\alpha_2,\alpha_3}(x) \;:\!=\; \alpha_1 \mathcal{H}(x) + \alpha_2 \log \left(\mathcal{H}(x) + \alpha_3 \right).\end{align*}

Recall from Section 2.1 that MH with landscape modification can be considered as a state-dependent version of energy transformation if we take $f(z) = z$ , $\alpha_3 = -c + \epsilon$ , with $\alpha_1$ , $\alpha_2$ chosen in a state-dependent manner:

\begin{align*}\alpha_1(x) = \textbf{1}_{\{\mathcal{H}(x) \leqslant c\}}, \quad \alpha_2(x) = \textbf{1}_{\{\mathcal{H}(x) > c\}}.\end{align*}

Note that MH with landscape modification can give rise to other kinds of energy transformation through different choices of f; see for example the quadratic case or the square root case in Section 2.2 or 2.3, respectively. The idea of mapping or transforming the function from $\mathcal{H}$ to $F(\mathcal{H})$ with F being strictly increasing and concave can be dated back to R. Azencott.

2.6.2. Preconditioning of the Hamiltonian

Landscape modification can be understood as a state-dependent preconditioning of the Hamiltonian $\mathcal{H}$ . Recall that in (2.4) we compute the acceptance–rejection probability in MH by

\begin{align*} & \exp\!\big(\!-\!\big(\mathcal{H}_{\epsilon,c}^f(y)-\mathcal{H}_{\epsilon,c}^f(x)\big)_+\big)\\[5pt] &= \begin{cases} 1 &\mbox{if } \mathcal{H}(y) \leqslant \mathcal{H}(x), \\[5pt] \exp\!\left(\!-\!\dfrac{1}{\epsilon}(\mathcal{H}(y) - \mathcal{H}(x))\right) &\mbox{if } c\geqslant \mathcal{H}(y) > \mathcal{H}(x), \\[9pt] \exp\!\left(\!-\!\dfrac{1}{\epsilon}(c - \mathcal{H}(x)) - \int_{c}^{\mathcal{H}(y)} \frac{1}{f(u-c) + \epsilon}\,du\right) &\mbox{if } \mathcal{H}(y) > c \geqslant \mathcal{H}(x), \\[9pt] \exp\!\left(\!-\!\int_{\mathcal{H}(x)}^{\mathcal{H}(y)} \dfrac{1}{f(u-c) + \epsilon}\,du\right) &\mbox{if } \mathcal{H}(y) > \mathcal{H}(x) > c. \end{cases}\end{align*}

On the set $\{c\geqslant \mathcal{H}(y) > \mathcal{H}(x)\}$ , the acceptance–rejection probability is the same as the original MH to allow for exploitation, while on the set $\{\mathcal{H}(y) > c \geqslant \mathcal{H}(x)\}$ and $\{\mathcal{H}(y) > \mathcal{H}(x) > c\}$ , the acceptance–rejection probability is higher than that of the original MH to encourage exploration of the landscape. Therefore, there is a higher transition rate of moving to other states when the algorithm is above the threshold c.

2.6.4 Quantum annealing

In quantum annealing [Reference Wang, Wu and Zou3], given a target Hamiltonian $\mathcal{H}$ and an initial Hamiltonian $\mathcal{H}_{init}$ that is usually easy to optimize, we optimize a time-dependent function $Q_t$ defined, for $t \in [0,T]$ , by

\begin{align*}Q_t(x) \;:\!=\; A(t) \mathcal{H}_{\textrm{init}}(x) + B(t) \mathcal{H}(x),\end{align*}

where A(t) and B(t) are smooth annealing schedules that satisfy $A(T) = B(0) = 0$ and T is the total annealing time. We also choose A(t) to be decreasing and B(t) to be increasing on the interval [0,T].

In simulated annealing with landscape modification, we also optimize a time-dependent function $\mathcal{H}^f_{\epsilon_t,c}$ which shares the same set of stationary points as the target $\mathcal{H}$ . In quantum annealing, $\mathcal{H}_{\textrm{init}}$ and $\mathcal{H}$ do not necessarily share the same set of stationary points. We mention the works Del Moral and Miclo [Reference Del Moral and Miclo11], Löwe [Reference Löwe12], and Frigerio and Grillo [Reference Frigerio and Grillo13] for simulated annealing with time-dependent energy function.

3.1. Main results

Without loss of generality, assume the external magnetic field is $h < 0$ . In the subcritical regime where $\frac{1}{\epsilon} > 1$ , it is known that there are two local minima of $g_{\epsilon}$ . We denote the global minimum of $g_{\epsilon}$ by $m_{-}^* < 0$ and the other local minimum by $m_+^* > 0$ , where $|m_-^*| > m_+^*$ , and let $z^*$ be the saddle point between $m_-^*$ and $m_+^*$ . We also write $m_-^*(N)$ (resp. $m_+^*(N)$ ) for the closest point in Euclidean distance to $m_-^*$ (resp. $m_+^*$ ) on $\Gamma_N$ .

Before we present the proof, we interpret the results in Theorem 1 intuitively: in Item 1, on the one hand we would like to choose c small enough so that the mapping

\begin{align*}m \mapsto \tanh\!\bigg(\dfrac{m+h}{f((E(m)-c)_+) + \epsilon}\bigg)\end{align*}

is flattened and only intersects the straight line $m \mapsto m$ at the global minimum $m_-^*$ . In this way the landscape of $g^f_{\epsilon,c}$ is transformed from a double-well to a single-well landscape, while the location of the global minimum at $m_-^*$ is preserved as that in the original landscape $g_{\epsilon}$ . This is illustrated in Figures 2 and 3. On the other hand, we cannot choose c to be too small if we are interested in seeking the ground state of $g_{\epsilon}$ , since otherwise, if $c < E(m_-^*)$ , then $m_-^*$ may no longer be the global minimum in the transformed free energy $g^f_{\epsilon,c}$ . This consequently yields a mean crossover time on the modified landscape that is subexponential in N, while the original mean crossover time is exponential in N, $1/\epsilon$ , and the original critical height $g_{\epsilon}(z^*) - g_{\epsilon}(m_+^*)$ . In Item 2 of Theorem 1, we choose a larger value of c compared with that in Item 1. Although the transformed free energy $g^f_{\epsilon,c}$ is not a convex function, it has a smaller critical height than the original free energy $g_{\epsilon}$ . This gives a reduced exponential dependence on the modified mean crossover time compared with the original mean crossover time.

Figure 2. Plots of the free energy $g_{\epsilon}$ and the modified free energy $g^{f}_{\epsilon,c} + E_{min}$ with $h = -0.05$ and $f(z) = z$ at temperature $\epsilon = 1/1.5$ , where $E_{min} = \min_{m \in [-1,1]} E(m)$ . We shift the modified free energy by $E_{min}$ so that it is on the same scale as the original free energy $g_{\epsilon}$ .

Figure 3. Plots of the mean-field equation (3.1) and the modified mean-field equation (3.4) with $h = -0.05$ and $f(z) = z$ at temperature $\epsilon = 1/1.5$ .

The power of landscape modification or energy transformation lies in tuning the parameter c appropriately. One way to tune c is to use the running minimum generated by the algorithm on the original free energy $g_{\epsilon}$ . Suppose we start in the well containing the local minimum $m_+^*$ ; setting c to be the running minimum eventually gives $c = E(m_+^*)$ , and hence Item 2 of Theorem 1 can be applied and the critical height on the modified landscape is reduced.

We illustrate Theorem 1 with a concrete numerical example in Figures 2 and 3, where we take $h = -0.05$ and $f(z) = z$ at temperature $\epsilon = 1/1.5$ . We numerically compute that $m_-^* = -0.8863$ , $m_+^* = 0.8188$ , and $z^* = 0.1524$ . As a result we have $E(m_-^*) = -0.4371$ , $E(m_+^*) = -0.2943$ , and $E(z^*) = -0.004$ . In the leftmost plots of Figures 2 and 3, we choose $c = -0.4 \in [E(m_-^*), E(m_+^*))$ . We numerically check that the conditions in Item 1 of Theorem 1 are satisfied, and we see that the location of the global minimum is the same for the blue curve as for the orange curve. In the rightmost plots of Figures 2 and 3, we choose $c = -0.2 \in [E(m_+^*), E(z^*)]$ . We see that for the blue curve and the red curve, the locations of the two local minima are the same, while the critical height is smaller than that in the original landscape $g_{\epsilon}$ .

3.2. Proof of Theorem 1

Before we give the proof, let us first recall the concept of critical height and the notation L as introduced in Section 2. In this section, we are interested in the CW model with and without landscape modification, with free energy $g^f_{\epsilon,c,N}$ and $g_{\epsilon,N}$ respectively. We therefore define the analogous concepts of critical height by inserting a subscript of N. This leads us to

\begin{align*} H^f_{\epsilon,c,N} &\;:\!=\; L\big(g^f_{\epsilon,c,N}\big) = \max_{x,y \in \mathcal{X}}\Big\{G\big(g^f_{\epsilon,c,N},x,y\big) - g^f_{\epsilon,c,N}(x) - g^f_{\epsilon,c,N}(y) \Big\} + \min_{x \in \Gamma_N} g^f_{\epsilon,c,N}(x), \\[5pt] H^0_{\epsilon,N} &\;:\!=\; L(g_{\epsilon,N}) = \max_{x,y \in \mathcal{X}}\Big\{G(g_{\epsilon,N},x,y) - g^f_{\epsilon,N}(x) - g_{\epsilon,N}(y) \Big\} + \min_{x \in \Gamma_N} g_{\epsilon,N}(x),\end{align*}

where we recall that G is introduced in (2.3).

We now proceed with the proof.

Proof of Theorem 1. First we prove Item 1. We observe that $\{E(m) \geqslant c\} = \{m \in [-h - \sqrt{h^2 - 2c},-h + \sqrt{h^2 - 2c}]\}$ . On this interval,

\begin{align*}\dfrac{d}{dm} g^f_{\epsilon,c}(m) = \textrm{arctanh}\,(m) - \dfrac{m+h}{f((E(m)-c)_+) + \epsilon} > 0,\end{align*}

and hence the modified free energy is strictly increasing on this interval. On the interval $\{m > -h + \sqrt{h^2 - 2c}\}$ ,

\begin{align*}\dfrac{d}{dm} g^f_{\epsilon,c}(m) = \dfrac{d}{dm} g_{\epsilon}(m) > 0,\end{align*}

as the original free energy is strictly increasing. On the interval $\{m < -h - \sqrt{h^2 - 2c}\}$ , we also have

\begin{align*}\dfrac{d}{dm} g^f_{\epsilon,c}(m) = \dfrac{d}{dm} g_{\epsilon}(m).\end{align*}

Thus, with these parameter choices, the only stationary point of $g^f_{\epsilon,c}$ is $m^*_-$ , which is the global minimum.

Next we proceed to proving (3.5). According to Löwe ([Reference Löwe12], Theorem 2.1), for $\xi_1(N)$ a polynomial function in N, we have

\begin{align*}\dfrac{1}{\lambda_2\Big(\!-\!M^f_{\epsilon,c,N}\Big)} \leqslant \xi_1(N) e^{N H^f_{\epsilon,c,N}}.\end{align*}

Using the random target lemma (Aldous and Fill, [Reference Aldous and Fill15], Section 4.2) and the above inequality leads to

(3.6) \begin{align} \pi^f_{\epsilon,c,N}(m^*_-(N)) \mathbb{E}_{m_+^*(N)} \Big(\tau^f_{m_-^*(N)}\Big) \leqslant \sum_{y \in \Gamma_N} \pi^f_{\epsilon,c,N}(y) \mathbb{E}_{m_+^*(N)} \left(\tau^f_{y}\right) &\leqslant (|\Gamma_N|-1) \dfrac{1}{\lambda_2\Big(\!-\!M^f_{\epsilon,c,N}\Big)} \nonumber \\[5pt] &\leqslant N \xi_1(N) e^{N H^f_{\epsilon,c,N}}. \end{align}

Now, let $\overline{m}(N) \;:\!=\; \arg \min g^f_{\epsilon,c,N}(m)$ and compute

\begin{align*} g^f_{\epsilon,c,N}(m^*_-(N)) - g^f_{\epsilon,c,N}(\overline{m}(N)) &= E^f_{\epsilon,c}(m^*_-(N)) - E^f_{\epsilon,c}(\overline{m}(N)) + I_N(m^*_-(N)) - I_N(\overline{m}(N)) \\[5pt] &= E^f_{\epsilon,c}(m^*_-(N)) - E^f_{\epsilon,c}(\overline{m}(N)) \\[5pt] &\quad + I_N(m^*_-(N)) - I(m^*_-(N)) + I(m^*_-(N)) - I_N(\overline{m}(N)) \\[5pt] &= E^f_{\epsilon,c}(m^*_-(N)) - E^f_{\epsilon,c}(\overline{m}(N)) \\[5pt] &\quad + [1+o(1)] \frac{1}{2 N} \ln \left(\frac{\pi N\left(1-m^*_-(N)^{2}\right)}{2}\right) \\[5pt] &\quad + I(m_-^*(N)) - I(\overline{m}(N)) \\[5pt] &\quad - [1+o(1)] \frac{1}{2 N} \ln \left(\frac{\pi N\left(1-\overline{m}(N)^{2}\right)}{2}\right) \\[5pt] &\rightarrow 0 \quad \text{as} \, N \to \infty, \end{align*}

where we use Bovier and den Hollander ([Reference Bovier and den Hollander14], Equation (13.2.5)) in the third equality, and $m_-^*(N), \overline{m}(N) \to m_-^*$ as $N \to \infty$ . The above computation combined with (3.6) yields

(3.7) \begin{align} \limsup_{N \to \infty} \mathbb{E}_{m_+^*(N)} \Big(\tau^f_{m_-^*(N)}\Big) \leqslant \lim_{N \to \infty} \dfrac{1}{N} \log H^f_{\epsilon,c,N} = 0. \end{align}

On the other hand, as the magnetization $(m^f(t))_{t \geqslant 0}$ is a birth–death process, using Equation (13.2.2) of Bovier and den Hollander [Reference Bovier and den Hollander14], the mean hitting time can be calculated explicitly as

\begin{align*} & \mathbb{E}_{m_{+}^{*}(N)}\Big(\tau^f_{m_{-}^{*}(N)}\Big) \\[5pt] &\quad= \sum_{\substack{m, m^{\prime} \in \Gamma_{N}, m \leqslant m^{\prime} \\ m_{-}^{*}(N)<m \leqslant m_{+}^{*}(N)}} \frac{\pi^f_{\epsilon,c, N}\left(m^{\prime}\right)}{\pi^f_{\epsilon,c, N}(m)} \frac{1}{M^f_{\epsilon,c, N}\left(m, m-2 N^{-1}\right)} \\ &\quad\geqslant \dfrac{1}{M^f_{\epsilon,c, N}\left(m_{+}^{*}(N), m_{+}^{*}(N)-2 N^{-1}\right)} = \dfrac{2}{1+m_{+}^{*}(N)} e^{ N(E^f_{\epsilon,c}(m_{+}^{*}(N)-2 N^{-1})-E^f_{\epsilon,c}(m))_+} \\ &\quad\geqslant \dfrac{2}{1+m_{+}^{*}(N)}. \end{align*}

As a result, as $m_+^*(N) \to m_+^*$ we have

(3.8) \begin{align} \liminf_{N \to \infty} \dfrac{1}{N} \log \mathbb{E}_{m_+^*(N)} \Big(\tau^f_{m_-^*(N)}\Big) \geqslant 0. \end{align}

Using both (3.8) and (3.7) gives (3.5).

Next, we prove Item 2, the argument for which closely follows that of Bovier and den Hollander ([Reference Bovier and den Hollander14], Theorem 13.1). Since we choose $c \in [E(m_+^*), E(z^*)]$ , we have $E(m_-^*) - c < E(m_+^*) - c \leqslant 0$ , and hence the landscape-modified mean-field equation (3.4) has at least two solutions $m_+^*$ and $m_-^*$ , which are exactly the same as those of the original mean-field equation (3.1). As the landscape-modified mean-field equation is continuous in m, there exists $\textbf{z}^* = \arg \max_{m_-^* \leqslant m \leqslant m_+^*} g^f_{\epsilon,c}(m)$ which also satisfies (3.4). Now, for $m \in \Gamma_N$ we consider

\begin{align*} N \big(E^f_{\epsilon,c}(m - 2N^{-1}) - E^f_{\epsilon,c}(m)\big) &= N \int_{E(m)}^{E(m) + 2N^{-1}(m+h - N^{-1})} \dfrac{1}{f((u-c)_+) + \epsilon} \, du \\[5pt] &\rightarrow \dfrac{2(m+h)}{f((E(m)-c)_+) + \epsilon} \quad \text{as } N \to \infty. \end{align*}

If we take $N \to \infty$ and $m \rightarrow \textbf{z}^*$ , we obtain

\begin{align*} \frac{1}{M^f_{\epsilon,c, N}\big(m, m-2 N^{-1}\big)} \to \dfrac{2}{1+\textbf{z}^*} \exp \bigg(\dfrac{2(\textbf{z}^*+h)_+}{f((E(\textbf{z}^*)-c)_+) + \epsilon}\bigg) = \dfrac{2}{1 - |\textbf{z}^*|}, \end{align*}

since, if $\textbf{z}^* > 0$ , then $\textbf{z}^* + h > 0$ and satisfies the mean-field equation (3.4). Using the mean hitting time formula again leads to the following, for any $\delta > 0$ :

(3.9) \begin{align} \mathbb{E}_{m_{+}^{*}(N)}\Big(\tau^f_{m_{-}^{*}(N)}\Big) &= \sum_{\substack{m, m^{\prime} \in \Gamma_{N}, m \leqslant m^{\prime} \\ m_{-}^{*}(N)<m \leqslant m_{+}^{*}(N)}} \frac{\pi^f_{\epsilon,c, N}\left(m^{\prime}\right)}{\pi^f_{\epsilon,c, N}(m)} \frac{1}{M^f_{\epsilon,c, N}\left(m, m-2 N^{-1}\right)} \nonumber \\ &= e^{N\left[g^f_{\epsilon,c,N}\left(\textbf{z}^{*}\right) - g^f_{\epsilon,c,N}\left(m_{+}^{*}\right)\right]} \frac{2}{1 - |\textbf{z}^*|}[1+o(1)] \nonumber \\ &\quad \times \sum_{\substack{m, m^{\prime} \in \Gamma_{N} \\ \left|m-\textbf{z}^{*}\right|<\delta,\left|m^{\prime}-m_{+}^{*}\right|<\delta}} \mathrm{e}^{N\left[g^f_{\epsilon,c, N}(m)-g^f_{\epsilon,c, N}\left(\textbf{z}^{*}\right)\right]- N\left[g^f_{\epsilon,c, N}\left(m^{\prime}\right)-g^f_{\epsilon,c, N}\left(m_{+}^{*}\right)\right]}. \end{align}

Using Stirling’s formula and the same argument as in Bovier and den Hollander ([Reference Bovier and den Hollander14], Equation (13.2.7)) yields

\begin{align*} e^{N\left[g^f_{\epsilon,c,N}\left(\textbf{z}^{*}\right) - g^f_{\epsilon,c,N}\left(m_{+}^{*}\right)\right]} &= (1+o(1)) e^{N\left[g^f_{\epsilon,c}\left(\textbf{z}^{*}\right) - g^f_{\epsilon,c}\left(m_{+}^{*}\right)\right]} \sqrt{\dfrac{1-\textbf{z}^{*2}}{1-m_+^{*2}}}. \end{align*}

Substitution into (3.9) gives

(3.10) \begin{align} & \mathbb{E}_{m_{+}^{*}(N)}\left(\tau^f_{m_{-}^{*}(N)}\right) \nonumber \\[5pt] &= e^{N\left[g^f_{\epsilon,c}\left(\textbf{z}^{*}\right) - g^f_{\epsilon,c}\left(m_{+}^{*}\right)\right]} \frac{2}{1 - |\textbf{z}^*|} \sqrt{\dfrac{1-\textbf{z}^{*2}}{1-m_+^{*2}}}[1+o(1)] \nonumber \\[5pt] &\quad \times \sum_{\substack{m, m^{\prime} \in \Gamma_{N}\\ \left|m-\textbf{z}^{*}\right|<\delta,\left|m^{\prime}-m_{+}^{*}\right|<\delta}}\sqrt{\dfrac{1-m^{2}}{1-\textbf{z}^{*2}}} \sqrt{\dfrac{1-m^{*2}_+}{1-m^{\prime 2}}} \mathrm{e}^{N\left[g^f_{\epsilon,c}(m)-g^f_{\epsilon,c}\left(\textbf{z}^{*}\right)\right]- N\left[g^f_{\epsilon,c}\left(m^{\prime}\right)-g^f_{\epsilon,c}\left(m_{+}^{*}\right)\right]}. \end{align}

We use a Laplace method argument to handle the sum in (3.10). Note that as f is assumed to be twice differentiable with $f(0) = f^{\prime}(0) = f^{\prime \prime}(0) = 0$ , this implies that $g^f_{\epsilon,c}$ is three times differentiable, and applying the third-order Taylor expansion gives

(3.11) \begin{align} g^f_{\epsilon,c}(m) - g^f_{\epsilon,c}\big(\textbf{z}^{*}\big) &= \dfrac{(m - \textbf{z}^{*})^2}{2} \big(g_{\epsilon,c}^{f}\big)^{\prime \prime}( \textbf{z}^*) + \mathcal{O}\left((m - \textbf{z}^{*})^3\right),\end{align}

(3.12) \begin{align} g^f_{\epsilon,c}\left(m^{\prime}\right)-g^f_{\epsilon,c}\left(m_{+}^{*}\right) &= \dfrac{(m^{\prime} - m_{+}^{*})^2}{2} \big(g_{\epsilon,c}^{f}\big)^{\prime \prime}\big( m_{+}^{*}\big) + \mathcal{O}\left((m^{\prime} - m_{+}^{*})^3\right), \end{align}

where we use $(g_{\epsilon,c}^{f})^{\prime}(\textbf{z}^*) = (g_{\epsilon,c}^{f})^{\prime}(m_{+}^{*}) = 0$ . Now we observe that the sum in (3.10) is

\begin{align*} &\quad (1+o(1)) \dfrac{N}{4} \int_{\mathbb{R}} \int_{\mathbb{R}} \exp \left[\frac{1}{2} \big(g_{\epsilon,c}^{f}\big)^{\prime \prime}\left(\textbf{z}^{*}\right) u^{2}-\frac{1}{2} \big(g_{\epsilon,c}^{f}\big)^{\prime \prime}\left(m_{+}^{*}\right) u^{\prime 2}\right]\, du du^{\prime} \\[5pt] &= (1+o(1)) \dfrac{N}{4} \frac{2 \pi}{\sqrt{\left[-\big(g_{\epsilon,c}^{f}\big)^{\prime \prime}\left(\textbf{z}^{*}\right)\right] \big(g_{\epsilon,c}^{f}\big)^{\prime \prime}\left(m_{+}^{*}\right)}}, \end{align*}

since $\big(g_{\epsilon,c}^{f}\big)^{\prime \prime}\left(\textbf{z}^{*}\right) < 0$ and $\big(g_{\epsilon,c}^{f}\big)^{\prime \prime}\left(m_{+}^{*}\right) > 0$ .

3.3. Extension to the random-field Curie–Weiss model

In Section 3.1, we discussed the classical CW model with landscape modification under a fixed magnetic field. In this section, we consider the random-field CW model with landscape modification. We discuss related metastability results and the ground-state free energy in such a setting, with the aim of illustrating that landscape modification can also be applied in the setting of a random energy landscape. Let us begin by recalling the random-field CW model. We shall adapt the setting of Mathieu and Picco [Reference Mathieu and Picco16]. Let $(h_i)_{i \in \mathbb{N}}$ be a sequence of independent and identically distributed random variables with $\mathbb{P}(h_i = 1) = \mathbb{P}(h_i = -1) = 1/2$ . We consider the random Hamiltonian function given, for a fixed $\theta > 0$ and $\sigma \in \{-1,1\}^N$ , by

\begin{align*}\textbf{H}_N(\sigma) = \textbf{H}_N(\sigma, \omega) &\;:\!=\; - \dfrac{N}{2} \textbf{m}_N(\sigma)^2 - \theta \sum_{i=1}^N h_i(\omega) \sigma_i \\[5pt] &= - \dfrac{N}{2} (\textbf{m}_N^+(\sigma) + \textbf{m}_N^-(\sigma))^2 - \theta N (\textbf{m}_N^+(\sigma) - \textbf{m}_N^-(\sigma)),\end{align*}

where

\begin{align*}\textbf{m}_N(\sigma) \;:\!=\; (1/N) \sum_{i=1}^N \sigma_i, \qquad\textbf{m}_N^+(\sigma) \;:\!=\; (1/N) \sum_{i=1;\; h_i = 1}^N \sigma_i, \qquad\textbf{m}_N^-(\sigma) \;:\!=\; (1/N) \sum_{i=1;\; h_i = -1}^N \sigma_i.\end{align*}

In the sequel, we shall suppress the dependency on $\omega$ . Denote the Gibbs distribution at temperature $\epsilon$ on $\{-1,1\}^N$ by

\begin{align*}\boldsymbol{\nu}_N(\sigma) \propto \exp\bigg\{-\frac{1}{\epsilon} \textbf{H}_N(\sigma)\bigg\}.\end{align*}

Let $N^+ \;:\!=\; |\{i;\;h_i = +1\}|$ , $N^- \;:\!=\; |\{i;\;h_i = -1\}|$ , and define the random set

\begin{align*}\boldsymbol{\Gamma}_N \;:\!=\; \bigg(\!-\!\frac{N^{+}}{N},-\frac{N^{+}}{N} + \frac{2}{N}, \ldots, \frac{N^{+}}{N}\bigg) \times\bigg(\!-\!\frac{N^{-}}{N},-\frac{N^{-}}{N} + \frac{2}{N}, \ldots, \frac{N^{-}}{N}\bigg).\end{align*}

For $\textbf{m} = (\textbf{m}^+, \textbf{m}^-) \in \boldsymbol{\Gamma}_N$ , with slight abuse of notation we write

\begin{align*}\textbf{H}_N(\textbf{m}) = - \dfrac{N}{2} (\textbf{m}^+ + \textbf{m}^-)^2- \theta N (\textbf{m}^+ - \textbf{m}^-) \;=\!:\; N \cdot \textbf{E}(\textbf{m}),\end{align*}

where $\textbf{m}^+$ and $\textbf{m}^-$ are the magnetization among the sites i where respectively $h_i = 1$ and $h_i = -1$ . Let $\boldsymbol{\pi}_{\epsilon,N}^0$ denote the image Gibbs distribution of $\boldsymbol{\nu}_N$ by $\boldsymbol{\Gamma}_{N}$ , where

\begin{align*} \boldsymbol{\pi}_{\epsilon,N}^0(\textbf{m}) &\propto \exp\bigg\{- \frac{1}{\epsilon} N \textbf{g}_{\epsilon,N}(\textbf{m})\bigg\}, \\[5pt] \textbf{g}_{\epsilon,N}(\textbf{m}) &\;:\!=\; - \dfrac{1}{2} (\textbf{m}^+ + \textbf{m}^-)^2- \theta (\textbf{m}^+ - \textbf{m}^-) - \frac{1}{\frac{1}{\epsilon} N} \log \left(\begin{array}{c} N^{+} \\[5pt] \frac{N^{+}}{2}+\textbf{m}^{+} \frac{N}{2} \end{array}\right)\left(\begin{array}{c} N^{-} \\[5pt] \frac{N^{-}}{2}+\textbf{m}^{-} \frac{N}{2} \end{array}\right).\end{align*}

As $N \to \infty$ , by the strong law of large numbers, $\textbf{g}_{\epsilon,N}$ converges almost surely to the free energy given by

\begin{align*}\textbf{g}_{\epsilon}(\textbf{m}) \;:\!=\; - \dfrac{1}{2} (\textbf{m}^+ + \textbf{m}^-)^2- \theta (\textbf{m}^+ - \textbf{m}^-) + \frac{1}{2\frac{1}{\epsilon}}\left(I(2\textbf{m}^+) + I(2\textbf{m}^-)\right),\end{align*}

where I(m) is the Cramér rate function as introduced in Section 3. The critical points of $\textbf{g}_{\epsilon}$ satisfy

(3.13) \begin{align} \textbf{m}^+ &= \dfrac{1}{2} \tanh\!\bigg(\frac{1}{\epsilon}(\textbf{m}^+ + \textbf{m}^- + \theta)\bigg),\end{align}

(3.14) \begin{align} \textbf{m}^- &= \dfrac{1}{2} \tanh\!\bigg(\frac{1}{\epsilon}(\textbf{m}^+ + \textbf{m}^- - \theta)\bigg).\end{align}

In this section, we shall only consider the subcritical regime where $\frac{1}{\epsilon} > \cosh^2(\frac{1}{\epsilon} \theta)$ . It can be shown (see e.g. [Reference Mathieu and Picco16]) that there are exactly three critical points. Let $\textbf{m}_{*} > 0$ be the unique positive solution to the mean-field equation

\begin{align*}\textbf{m}_{*} = \dfrac{1}{2}\bigg(\tanh\!\bigg(\frac{1}{\epsilon}(\textbf{m}_{*} + \theta)\bigg) + \tanh\!\bigg(\frac{1}{\epsilon}(\textbf{m}_{*} - \theta)\bigg)\bigg).\end{align*}

The three critical points of $\textbf{g}_{\epsilon}$ are given by

\begin{align*} \textbf{m}_{0} &= \bigg(\frac{1}{2} \tanh \bigg(\frac{1}{\epsilon} \theta\bigg),-\frac{1}{2} \tanh \bigg(\frac{1}{\epsilon} \theta\bigg)\bigg), \\[5pt] \textbf{m}_{1} &= \bigg(\frac{1}{2} \tanh \bigg(\frac{1}{\epsilon} \textbf{m}_{*}+\frac{1}{\epsilon} \theta\bigg), \frac{1}{2} \tanh \bigg(\frac{1}{\epsilon} \textbf{m}_{*}-\frac{1}{\epsilon} \theta\bigg)\bigg), \\[5pt] \textbf{m}_{2} &= \bigg(\frac{1}{2} \tanh \bigg(\!-\!\frac{1}{\epsilon} \textbf{m}_{*}+\frac{1}{\epsilon} \theta\bigg),-\frac{1}{2} \tanh \bigg(\frac{1}{\epsilon} \textbf{m}_{*}+\frac{1}{\epsilon} \theta\bigg)\bigg),\end{align*}

where $\textbf{m}_0$ is the saddle point and $\textbf{m}_1, \textbf{m}_2$ are the two global minima. Consider the continuized Glauber dynamics $(\sigma_N(t))_{t \geqslant 0}$ by picking a node uniformly at random and changing the sign of the selected spin, while targeting the Gibbs distribution $\boldsymbol{\nu}_N$ at temperature $\epsilon$ . Denote by $\textbf{m}_N(t) \;:\!=\; \textbf{m}_N(\sigma_N(t))$ the induced dynamics on the magnetization, and denote its infinitesimal generator by $\textbf{M}^0_{\epsilon,N}$ . This is proven to be a Markov chain in Mathieu and Picco [Reference Mathieu and Picco16], with stationary measure $\boldsymbol{\pi}^0_{\epsilon,N}$ .

Now let us consider the landscape-modified Hamiltonian on $\boldsymbol{\Gamma}_N$ ,

(3.15) \begin{align} \textbf{E}^f_{\epsilon,c}(\textbf{m}) &\;:\!=\; \int_{d}^{\textbf{E}(\textbf{m})} \dfrac{1}{f((u-c)_+) + \epsilon} \, du,\end{align}

(3.16) \begin{align} \textbf{H}^f_{\epsilon,c,N}(\textbf{m}) &\;:\!=\; N \cdot \textbf{E}^f_{\epsilon,c}(\textbf{m}),\end{align}

where $d \in \mathbb{R}$ can be chosen arbitrarily since we are only interested in the difference of $\textbf{E}^f_{\epsilon,c}$ . The transformed image Gibbs distribution is therefore

\begin{align*} \boldsymbol{\pi}_{\epsilon,c,N}^f(\textbf{m}) &\propto \exp\{- N \textbf{g}^f_{\epsilon,c,N}(\textbf{m})\}, \\[5pt] \textbf{g}^f_{\epsilon,c,N}(\textbf{m}) &\;:\!=\; \textbf{E}^f_{\epsilon,c}(\textbf{m}) - \frac{1}{N} \log \left(\begin{array}{c} N^{+} \\[5pt] \frac{N^{+}}{2}+\textbf{m}^{+} \frac{N}{2} \end{array}\right)\left(\begin{array}{c} N^{-} \\[5pt] \frac{N^{-}}{2}+\textbf{m}^{-} \frac{N}{2} \end{array}\right).\end{align*}

The strong law of large numbers yields that as $N \to \infty$ , $\textbf{g}^f_{\epsilon,c,N}$ converges almost surely to the transformed free energy

\begin{align*}\textbf{g}^f_{\epsilon,c}(\textbf{m}) \;:\!=\; \textbf{E}^f_{\epsilon,c}(\textbf{m}) + \frac{1}{2}\left(I(2\textbf{m}^+) + I(2\textbf{m}^-)\right).\end{align*}

The critical points of $\textbf{g}^f_{\epsilon,c}$ satisfy the following landscape-modified mean-field equations:

(3.17) \begin{align} \textbf{m}^+ &= \dfrac{1}{2} \tanh\!\bigg(\dfrac{\textbf{m}^+ + \textbf{m}^- + \theta}{f((\textbf{E}(\textbf{m})-c)_+) + \epsilon}\bigg),\end{align}

(3.18) \begin{align} \textbf{m}^- &= \dfrac{1}{2} \tanh\!\bigg(\dfrac{\textbf{m}^+ + \textbf{m}^- - \theta}{f((\textbf{E}(\textbf{m})-c)_+) + \epsilon}\bigg).\end{align}

Note that (3.17) and (3.18) reduce to the classical case (3.13) and (3.14) if we take $f = 0$ . Consider the continuized Glauber dynamics $(\sigma_N^f(t))_{t \geqslant 0}$ by picking a node uniformly at random and changing the sign of the selected spin, while targeting the Gibbs distribution with Hamiltonian $\epsilon \textbf{H}^f_{\epsilon,c,N}$ at temperature $\epsilon$ . Denote by $\textbf{m}_N^f(t) \;:\!=\; \textbf{m}^f_N(\sigma_N^f(t))$ the induced dynamics on the magnetization, and denote by $\textbf{M}^f_{\epsilon,c,N}$ its infinitesimal generator, which is a Markov chain with stationary measure $\boldsymbol{\pi}^f_{\epsilon,c,N}$ .

In the following, we shall consider the case where $c \in [\textbf{E}(\textbf{m}_1), \textbf{E}(\textbf{m}_0)]$ . It can be seen that the two global minima of $\textbf{g}^f_{\epsilon,c}$ remain at $\textbf{m}_1, \textbf{m}_2$ with this choice of c. For any path $\gamma^{\textbf{m}_1,\textbf{m}_0}$ connecting $\textbf{m}_1$ and $\textbf{m}_0$ , we define

\begin{align*} \textbf{m}_3(\gamma^{\textbf{m}_1,\textbf{m}_0}) &\;:\!=\; \arg \max \{\textbf{g}^f_{\epsilon,c}(\gamma_i);\;\gamma_i \in \gamma^{\textbf{m}_1,\textbf{m}_0}\}, \\[5pt] \Delta \textbf{g}^f_{\epsilon,c} &\;:\!=\; \min_{\gamma^{\textbf{m}_1,\textbf{m}_0}} \textbf{m}_3(\gamma^{\textbf{m}_1,\textbf{m}_0}) - \textbf{g}^f_{\epsilon,c}(\textbf{m}_1) = \min_{\gamma^{\textbf{m}_1,\textbf{m}_0}} \textbf{m}_3(\gamma^{\textbf{m}_1,\textbf{m}_0}) - \textbf{g}^f_{\epsilon,c}(\textbf{m}_0),\end{align*}

where $\Delta \textbf{g}^f_{\epsilon,c}$ is the critical height on the modified landscape. We also write $\Delta \textbf{g}_{\epsilon}$ to denote the critical height on the original landscape. Suppose that $\Delta \textbf{g}_{\epsilon}$ is attained at $\textbf{m}_4$ , so that $\Delta \textbf{g}_{\epsilon} = \textbf{g}_{\epsilon}(\textbf{m}_4) - \textbf{g}_{\epsilon}(\textbf{m}_0)$ ; we deduce

\begin{align*} \Delta \textbf{g}^f_{\epsilon,c} \leqslant \textbf{g}^f_{\epsilon,c}(\textbf{m}_4) - \textbf{g}^f_{\epsilon,c}(\textbf{m}_0) &= \int_{\textbf{E}(\textbf{m}_0)}^{\textbf{E}(\textbf{m}_4)} \dfrac{1}{f((u-c)_+) + \epsilon} \, du \\[5pt] &\quad + \frac{1}{2}\big(I\big(2\textbf{m}_4^+\big) + I\big(2\textbf{m}_4^-\big)\big) - \frac{1}{2}\big(I\big(2\textbf{m}_0^+\big) + I\big(2\textbf{m}_0^-\big)\big) \\[5pt] &\leqslant \frac{1}{\epsilon} \left(\textbf{g}_{\epsilon}(\textbf{m}_4) - \textbf{g}_{\epsilon}(\textbf{m}_0)\right) = \frac{1}{\epsilon} \Delta \textbf{g}_{\epsilon}.\end{align*}

In other words, the critical height of the free energy in the modified landscape is bounded above by $\frac{1}{\epsilon}$ times the critical height of the free energy in the original landscape.

A direct application of Mathieu and Picco ([Reference Mathieu and Picco16], Theorem 2.7) yields the following result on the asymptotics of the spectral gap.

Theorem 2. (Asymptotics of the spectral gap.) Suppose $\theta > 0$ , $\frac{1}{\epsilon} > \cosh^2\left(\frac{1}{\epsilon} \theta\right)$ are fixed, and $\textbf{m}_0$ is the saddle point, while $\textbf{m}_1, \textbf{m}_2$ are the two global minima on the original free energy landscape $\textbf{g}_{\epsilon}$ . For $c \in [\textbf{E}(\textbf{m}_1), \textbf{E}(\textbf{m}_0)]$ , we have $\mathbb{P}$ -almost surely that

\begin{align*} \lim_{N \to \infty} \dfrac{1}{N} \log \lambda_2(\!-\!\textbf{M}^f_{\epsilon,c,N}) &= - \Delta \textbf{g}^f_{\epsilon,c} \geqslant - \frac{1}{\epsilon} \Delta \textbf{g}_{\epsilon} = \lim_{N \to \infty} \dfrac{1}{N} \log \lambda_2(\!-\!\textbf{M}^0_{\epsilon,N}). \end{align*}

In essence, the relaxation time in the mean-field limit of the transformed generator $\textbf{M}^f_{\epsilon,c,N}$ is asymptotically less than or equal to that of the original generator $\textbf{M}^0_{\epsilon,N}$ .

This yields a reduced exponential dependence of the relaxation time on the modified landscape, compared with the relaxation time on the original landscape.

4.1. Numerical illustrations

Before discussing the proofs of the main results above, we illustrate the convergence performance of simulated annealing with landscape modification, which we call improved simulated annealing (ISA), and compare it to the performance of the classical simulated annealing algorithm (SA) for the travelling salesman problem (TSP). We first state the parameters that we used to generate the numerical results:

• • TSP and its objective function. We select 50 nodes uniformly at random on the grid $[0,100] \times [0,100]$ . The objective is to find a configuration that minimizes the total Euclidean distance with the same starting and ending point. Each node can only be visited once.

• • Initial configuration. Both ISA and SA have the same initialization. They are initialized using the output of the nearest-neighbour algorithm: a node is randomly chosen as the starting point, which is then connected to the closest unvisited node. This is repeated until every node has been visited, and subsequently the last node is connected back to the starting node.

• • Proposal chain. Both ISA and SA share the same proposal chain: at each step, a proposal move is generated using the 2-OPT algorithm [Reference Croes18].

• • Acceptance–rejection mechanism. In SA, the proposed move is accepted with probability

  \begin{align*}\min \left\{ 1,e^{\beta(\mathcal{H}(x)-\mathcal{H}(y))} \right\},\end{align*}
  while in ISA, the acceptance probability is computed as in Section 2.1. In other words, we use a linear f in ISA. Both SA and ISA have the same source of randomness.

• • Cooling schedule. Both ISA and SA use the same logarithmic cooling schedule, of the form

  \begin{align*}\epsilon_t = \dfrac{\sqrt{50}}{\ln(t+1)}.\end{align*}

• • Choice of c in ISA. In this experiment, if we denote the proposal configuration at time t by $y_t$ , we set $c = c_t$ to be

  \begin{align*}c_t = \mathcal{H}(y_t) - 5.\end{align*}
  This tuning strategy has been discussed in Section 2.5.

• • Number of iterations. We run both ISA and SA for 100,000 iterations.

We generate 1000 random TSP instances. For each instance we compute what we call the improvement percentage ( $\textrm{IP}$ ) of ISA over SA, defined by

\begin{align*}\textrm{IP} \;:\!=\; 100 \dfrac{\min_{t \in [0,100,000]} \mathcal{H}(X^0(t)) - \min_{t \in [0,100,000]} \mathcal{H}(X^f(t))}{\min_{t \in [0,100,000]} \mathcal{H}(X^0(t))}.\end{align*}

The summary statistics of $\textrm{IP}$ are provided in Table 1, while its histogram over these 1000 instances can be found in Figure 5. The code for reproducing these results can be found at https://github.com/mchchoi/Improved-discrete-simulated-annealing.

The summary statistics in Table 1 and the histogram in Figure 5 offer empirical evidence for using ISA over SA: out of the 1000 TSP instances, there are 798 instances in which the $\textrm{IP}$ is non-negative. The sample mean of the $\textrm{IP}$ is approximately 1.87%, while its sample median is 1.47%.

Next, in Figure 6, we investigate the difference between SA and ISA in a particular instance. We see that SA (the blue curve) is stuck at a local minimum, while ISA (the orange curve) is able to escape the local minimum, owing to the increased acceptance probability compared with SA, and it reaches regions where the objective value is smaller than that of SA.

Table 1. Summary statistics of $\textrm{IP}$ on 1000 random TSP instances.

Figure 5. Histogram showing IP of ISA over SA on 1000 randomly generated TSP instances.

Figure 6. TSP objective value against iteration of ISA and SA.

The rest of this section is devoted to the proofs of Theorem 3, Corollary 1, Lemma 1, and Theorem 4.

4.2 Proof of Theorem 3

First, using the classical result of Holley and Stroock ([Reference Holley and Stroock19], Theorem 2.1), it is immediate that

\begin{align*} C_2 e^{-H^f_{\epsilon,c}} \leqslant \lambda_2(\!-\!M^f_{\epsilon,c}) \leqslant C_3 e^{-H^f_{\epsilon,c}}.\end{align*}

For any arbitrary $x_1, x_2 \in \{\mathcal{H}(x_1) \geqslant \mathcal{H}(x_2)\}$ , we deduce the following upper bound:

\begin{align*}& \mathcal{H}^f_{\epsilon,c}(x_1) - \mathcal{H}^f_{\epsilon,c}(x_2)\\[5pt] &= \int_{\mathcal{H}(x_2)}^{\mathcal{H}(x_1)} \dfrac{1}{f((u-c)_+) + \epsilon}\,du \\[5pt] &= \begin{cases} \dfrac{1}{\epsilon}(\mathcal{H}(x_1) - \mathcal{H}(x_2)) &\mbox{if } c\geqslant \mathcal{H}(x_1) > \mathcal{H}(x_2), \\[9pt] \dfrac{1}{\epsilon}(c - \mathcal{H}(x_2)) + \int_{c}^{\mathcal{H}(x_1)} \dfrac{1}{f(u-c) + \epsilon}\,du &\mbox{if } \mathcal{H}(x_1) > c \geqslant \mathcal{H}(x_2), \\[9pt] \int_{\mathcal{H}(x_2)}^{\mathcal{H}(x_1)} \dfrac{1}{f(u-c) + \epsilon}\,du &\mbox{if } \mathcal{H}(x_1) > \mathcal{H}(x_2) > c, \end{cases} \\[9pt] &\leqslant \begin{cases} \dfrac{1}{\epsilon}(\mathcal{H}(x_1) \wedge c - \mathcal{H}(x_2) \wedge c) &\mbox{if } c\geqslant \mathcal{H}(x_1) > \mathcal{H}(x_2), \\[9pt] \dfrac{1}{\epsilon}(\mathcal{H}(x_1) \wedge c - \mathcal{H}(x_2)\wedge c) + \ln\!\bigg(1+\dfrac{1}{\epsilon}(\max \mathcal{H}-c)\bigg) &\mbox{if } \mathcal{H}(x_1) > c \geqslant \mathcal{H}(x_2), \\[9pt] \dfrac{1}{f(\delta)} \left(\max \mathcal{H} - \min \mathcal{H}\right) &\mbox{if } \mathcal{H}(x_1) > \mathcal{H}(x_2) > c. \end{cases}\end{align*}

As a result, $C_2^{-1}e^{H^f_{\epsilon,c}} \leqslant e^{\frac{1}{\epsilon} c^*} C_1(\epsilon)$ . On the other hand, we have the following lower bound:

\begin{align*} \mathcal{H}^f_{\epsilon,c}(x_1) - \mathcal{H}^f_{\epsilon,c}(x_2) &\geqslant \begin{cases} \dfrac{1}{\epsilon}(\mathcal{H}(x_1) \wedge c - \mathcal{H}(x_2) \wedge c) &\mbox{if } c\geqslant \mathcal{H}(x_1) > \mathcal{H}(x_2), \\[9pt] \dfrac{1}{\epsilon}(\mathcal{H}(x_1) \wedge c - \mathcal{H}(x_2)\wedge c) &\mbox{if } \mathcal{H}(x_1) > c \geqslant \mathcal{H}(x_2), \\[9pt] \dfrac{1}{f(\mathcal{H}(x_1) - c)} \left(\mathcal{H}(x_1) - \mathcal{H}(x_2)\right) &\mbox{if } \mathcal{H}(x_1) > \mathcal{H}(x_2) > c, \end{cases}\end{align*}

and hence $e^{H^f_{\epsilon,c}} \geqslant e^{\frac{1}{\epsilon} c^*} C_4^{-1}$ .

4.3. Proof of Corollary 1

We first prove Item 1. For a continuous-time reversible Markov chain, by Levin and Peres ([Reference Levin and Peres20], Theorems 12.5 and 20.6), we can bound the total-variation mixing time by relaxation time via

\begin{align*}\dfrac{1}{\lambda_2(\!-\!M^f_{\epsilon,c}\Big)} \log 2 \leqslant t_{mix}\Big(M^f_{\epsilon,c},1/4\Big) \leqslant \dfrac{1}{\lambda_2(\!-\!M^f_{\epsilon,c})} \log\!\left(\dfrac{4}{\pi^f_{min}}\right),\end{align*}

where

\begin{align*}\pi_{min}^f \;:\!=\; \min_x \pi^f_{\epsilon,c}(x) = \dfrac{\mu(x^*)}{Z^f_{\epsilon,c}}\end{align*}

for some $x^* \in S_{\textrm{min}}$ , and

\begin{align*}Z^f_{\epsilon,c} \;:\!=\; \sum_{x \in \mathcal{X}} e^{-\mathcal{H}_{\epsilon,c}^f(x)}\mu(x)\end{align*}

is the normalization constant. Note that since $ \ln Z^f_{\epsilon,c} \to \ln \mu(S_{\textrm{min}})$ , we have

\begin{align*}\lim_{\epsilon \to 0} \epsilon \ln Z^f_{\epsilon,c} = 0.\end{align*}

Item 1 follows from collecting the above results together with Theorem 3.

Next, we prove Item 2. First, using the random target lemma ([Reference Aldous and Fill15], Section 4.2), we have

\begin{align*}\pi^f_{\epsilon,c}(\eta) \mathbb{E}_{\sigma}\big(\tau^f_\eta\big) \leqslant \sum_{x \in \mathcal{X}} \pi^f_{\epsilon,c}(x) \mathbb{E}_{\sigma}\big(\tau^f_x\big) \leqslant (|\mathcal{X}|-1) \dfrac{1}{\lambda_2(\!-\!M^f_{\epsilon,c})}.\end{align*}

Since $\mathcal{H}(\eta) = 0$ and $Z^f_{\epsilon,c} \leqslant 1$ , upon rearranging and using Theorem 2 we obtain

\begin{align*}\limsup_{\epsilon \to 0} \epsilon \log \mathbb{E}_{\sigma}\big(\tau^f_\eta\big) \leqslant \lim_{\epsilon \to 0} \epsilon \log \dfrac{1}{\lambda_2\Big(\!-\!M^f_{\epsilon,c}\Big)} = c^*.\end{align*}

Define the equilibrium potential and capacity of the pair $(\sigma,\eta)$ as in Bovier and den Hollander ([Reference Bovier and den Hollander14], Chapter 7.2) to be respectively

\begin{align*} h_{\sigma,\eta}^{f}(x) &\;:\!=\; \mathbb{P}_x\big(\tau^f_{\sigma} < \tau^f_{\eta}\big), \\[5pt] \mathrm{cap}^{M^f_{\epsilon,c}}(\sigma,\eta) &\;:\!=\; \inf_{f: f|_{A} = 1, f|_{B} = 0} \langle -M^f_{\epsilon,c}f,f \rangle_{\pi^f_{\epsilon,c}} = \langle -M^f_{\epsilon,c} h_{\sigma,\eta}^{f}, h_{\sigma,\eta}^{f} \rangle_{\pi^f_{\epsilon,c}}.\end{align*}

If we prove that

(4.4) \begin{align} \mathrm{cap}^{M^f_{\epsilon,c}}(\sigma,\eta) \leqslant \dfrac{1}{Z^f_{\epsilon,c}} \left(\sum_{x,y} \mu(x) Q(x,y)\right) e^{-\frac{1}{\epsilon} G^0(\sigma,\eta)},\end{align}

then the mean hitting time formula combined with the equilibrium potential and capacity leads to

\begin{align*} \mathbb{E}_{\sigma}(\tau^f_{\eta}) &= \dfrac{1}{\mathrm{cap}^{M^f_{\epsilon,c}}(\sigma,\eta)} \sum_{y \in \mathcal{X}} \pi_{\epsilon,c}^f(y) h_{\sigma,\eta}^{f}(y) \\[5pt] &\geqslant \dfrac{1}{\mathrm{cap}^{M^f_{\epsilon,c}}(\sigma,\eta)} \pi_{\epsilon,c}^f(\sigma) h_{\sigma, \eta}^{f}(\sigma) \geqslant \dfrac{\mu(\sigma)}{\sum_{x,y} \mu(x) Q(x,y)} e^{\frac{1}{\epsilon} c^*},\end{align*}

and the desired result follows, since

\begin{align*}\liminf_{\epsilon \to 0} \epsilon \log \mathbb{E}_{\sigma}\big(\tau^f_\eta\big) \geqslant c^*.\end{align*}

It therefore remains to prove (4.4). Define

\begin{align*}\Phi(\sigma,\eta) \;:\!=\; \{x\in\mathcal{X};\; G^0(x, \sigma) \leqslant G^0(x,\eta)\}.\end{align*}

Let $\textbf{1}_{A}$ denote the indicator function of the set A; the Dirichlet principle of capacity gives

\begin{align*} \mathrm{cap}^{M^f_{\epsilon,c}}(\sigma,\eta) \leqslant \langle -M^f_{\epsilon,c}\textbf{1}_{\Phi(\sigma,\eta)},\textbf{1}_{\Phi(\sigma,\eta)}\rangle_{\pi^f_{\epsilon,c}} &= \dfrac{1}{Z^f_{\epsilon,c}}\sum_{x \in \Phi(\sigma,\eta), y \notin \Phi(\sigma,\eta)} e^{-\frac{1}{\epsilon} \left( H(x) \vee H(y) \right)} \mu(x) Q(x,y) \\[5pt] &\leqslant \dfrac{1}{Z^f_{\epsilon,c}} \left(\sum_{x,y} \mu(x) Q(x,y)\right) e^{-\frac{1}{\epsilon} G^0(\sigma,\eta)},\end{align*}

where in the last inequality we use the fact that $G^0(\sigma,\eta)$ is the lowest possible highest elevation between $\sigma$ and $\eta$ .

4.4. Proof of Lemma 1

We first prove Item 1. The lower bound is immediate, while the upper bound can be deduced via

\begin{align*}\epsilon_t \mathcal{H}^f_{\epsilon_t,c} \leqslant \epsilon_t \int_{\mathcal{H}_{\textrm{min}}}^{\mathcal{H}(x)} \dfrac{1}{\epsilon_t} \,du \leqslant M.\end{align*}

Next, we prove Item 2. We consider

\begin{align*} \dfrac{\partial}{\partial t} \epsilon_t \mathcal{H}^f_{\epsilon_t,c}(x) &= \mathcal{H}^f_{\epsilon_t,c}(x) \left(\dfrac{\partial}{\partial t} \epsilon_t\right) + \epsilon_t \dfrac{\partial}{\partial t} \mathcal{H}^f_{\epsilon_t,c}(x) \\[5pt] &= \mathcal{H}^f_{\epsilon_t,c}(x) \dfrac{-c^* - \epsilon}{(\ln(t+1))^2} \dfrac{1}{t+1} + \epsilon_t \epsilon_t^{\prime} \int_{\mathcal{H}_{\textrm{min}}}^{\mathcal{H}(x)} - \dfrac{1}{(f((u-c)_+) + \epsilon_t)^2} \, du.\end{align*}

This leads to

\begin{align*} \left|\dfrac{\partial}{\partial t} \epsilon_t \mathcal{H}^f_{\epsilon_t,c}(x)\right| \leqslant \dfrac{2M}{(\ln(t+1))(t+1)}.\end{align*}

Thirdly, we prove Item 3. Using Item 2, we calculate that

\begin{align*} \beta_t^{\prime} M + \beta_t R_t &\leqslant \dfrac{M}{(c^*+\epsilon)(t+1)} + \dfrac{2M}{(c^*+\epsilon)(1+t)} = \dfrac{3M}{(c^*+\epsilon)(1+t)}.\end{align*}

Finally, we prove Item 4. Following exactly the same calculation as in the proof of Löwe ([Reference Löwe12], Lemma 3.5), we see that

\begin{align*} \left\lVert g - \pi_{\epsilon_t,c}^f(g)\right\rVert_{\ell^p(\pi^f_{\epsilon_t,c})}^2 &\leqslant A \left(1 + \beta_t(\max \mathcal{H} - c)\right) e^{\beta_t(c^* + M \frac{p-2}{p})}\langle -M^f_{\epsilon_t,c}g,g \rangle_{\pi^f_{\epsilon_t,c}} \\[5pt] &\leqslant A e^{\beta_t(c^* + M \frac{p-2}{p} + \max \mathcal{H} - c)}\langle -M^f_{\epsilon_t,c}g,g \rangle_{\pi^f_{\epsilon_t,c}}.\end{align*}

The desired result follows if we let

\begin{align*}\epsilon = \dfrac{M(p-2)}{p} + \max \mathcal{H} - c,\end{align*}

so that

\begin{align*}p = \dfrac{2M}{M + \max \mathcal{H} - \epsilon - c}.\end{align*}

4.5 Proof of Theorem 4

We would like to invoke the results in Löwe [Reference Löwe12] for a time-dependent target function in simulated annealing. In Lemma 1, Items 1, 2, 3, and 4, we verify that Equations (11) and (12) and Assumptions (A1) and (A2), respectively, from Löwe [Reference Löwe12] hold. Consequently, if we let $h_t(y) \;:\!=\; \mathbb{P}_x(X^f_{\epsilon_t,c}(t) = y)/\pi^f_{\epsilon_t,c}(y)$ , then according to Holley and Stroock ([Reference Holley and Stroock19], Lemma 1.7), its $\ell^2$ norm is bounded by

\begin{align*}\left\lVert h_t\right\rVert_{\ell^2(\pi^f_{\epsilon_t,c})} \leqslant 1+K^{\frac{1}{2\overline{\epsilon}}},\end{align*}

for $t \geqslant e^{1/\overline{\epsilon}}-1$ . The desired result follows from exactly the same argument as in Löwe ([Reference Löwe12], Theorem 3.8).

Now we calculate $\pi^f_{\epsilon_t,c}(\mathcal{X} \backslash S_{\textrm{min}})$ . For $c \geqslant \underline{d}$ , we compute that

\begin{align*} \pi^f_{\epsilon_t,c}(\mathcal{X} \backslash S_{\textrm{min}}) &= \dfrac{\sum_{x \in \mathcal{X} \backslash S_{\textrm{min}}} e^{-\mathcal{H}_{\epsilon_t,c}^f(x) }\mu(x)}{\sum_{x \in \mathcal{X}} e^{-\mathcal{H}_{\epsilon_t,c}^f(x)}\mu(x)} \\[5pt] &\leqslant \dfrac{\sum_{x;\; \mathcal{H}(x) \geqslant \underline{d}} \mu(x)}{\sum_{x; \mathcal{H}(x) \leqslant \underline{d}} \exp\{\int_{\mathcal{H}(x)}^{\underline{d}} \frac{1}{\epsilon_t}\,du\}\mu(x)} \\[5pt] &\leqslant e^{-\frac{1}{\epsilon_t}(\underline{d} - \mathcal{H}_{\textrm{min}})} \dfrac{1}{\mu(S_{\textrm{min}})}.\end{align*}

On the other hand, for $\mathcal{H}_{\textrm{min}} \leqslant c < \underline{d}$ ,

\begin{align*} \pi^f_{\epsilon_t,c}(\mathcal{X} \backslash S_{\textrm{min}}) &= \dfrac{\sum_{x \in \mathcal{X} \backslash S_{\textrm{min}}} e^{-\mathcal{H}_{\epsilon_t,c}^f(x) }\mu(x)}{\sum_{x \in \mathcal{X}} e^{-\mathcal{H}_{\epsilon_t,c}^f(x)}\mu(x)} \\[5pt] &\leqslant \dfrac{1}{\mu(S_{\textrm{min}})}\exp\bigg\{-\int_{\mathcal{H}_{\textrm{min}}}^{\underline{d}} \dfrac{1}{f((u-c)_+) + \epsilon_t}du\bigg\} \\[5pt] &\leqslant \dfrac{1}{\mu(S_{\textrm{min}})}\exp\bigg\{-\int_{c}^{\underline{d}} \dfrac{1}{f(u-c) + \epsilon_t}du\bigg\}.\end{align*}