A First Runtime Analysis of Parallel Tempering in Quadratic Optimization

1 Introduction

Parallel tempering, also known as the replica exchange method, is a Markov Chain Monte Carlo (MCMC) technique for finding low-energy states in complex landscapes by executing multiple processes at different temperatures and allowing for states to migrate between parallel processes based on the Metropolis criterion. The intuition behind this approach is that configurations at low temperatures can efficiently exploit landscape features by stochastic hill-climbing, whereas configurations at high temperatures are more likely to escape local minima to allow for adequate exploration. Allowing processes at different temperatures to periodically exchange states results in a robust ensemble that balances intensification and diversification during the search process.

Originally introduced as a technique in statistical mechanics for simulating physical systems in order to efficiently compute thermodynamical properties such as specific heat [Swendsen and Wang, 1986; Geyer, 1991], parallel tempering has recently gained interest as a randomized search heuristic for optimization. For example, Zhu et al. 2020 pointed out the relationship between the statistical physics of Ising spin glasses and hard Boolean optimization problems and observed that parallel tempering can be used to tackle this kind of problem. Their borealis solver, based on parallel tempering, was ranked first for randomMaxSAT instances in the Incomplete Solver track at the Eleventh Max-SAT Evaluation at SAT 2016.

Implementations of parallel tempering have recently emerged as digital annealing hardware for solving quadratic unconstrained optimization (QUBO) problems [Aramon et al., 2019], as well as an API for a multi-core implementation [Almeida et al., 2025]. The approach has been successfully applied for decades on a wide range of problems such as overcoming local minima in numerical simulations of molecules with rough energy landscapes [Hansmann, 1997], Bayesian data analysis [Habeck et al., 2005], and makespan minimization on instances of various NP-hard scheduling problems [Almeida et al., 2025].

Despite a long and successful history in various applications and many solid experimental results, so far a rigorous theoretical foundation for parallel tempering is missing. In particular, it is not clear when and why parallel tempering is more powerful than the Metropolis algorithm maintaining a single state.

In this paper, we make a first step towards establishing a rigorous theoretical foundation for parallel tempering by conducting a first runtime analysis on a simple QUBO problem with a tunably rugged landscape. We prove that the classical Metropolis algorithm fails to efficiently optimize the problem when the temperature is set too low, in which case it cannot easily escape local minima, or when the temperature is set too high, in which case its behavior is too chaotic to effectively leverage the global structure of the energy landscape. We furthermore show that, for certain settings of energy barrier parameters in the landscape, these temperature regimes will overlap and consequently there is no ideal temperature setting that guarantees success. In contrast, we prove that a two-state parallel tempering approach using a low and high temperature and asymmetric migration rather than symmetric exchanges typically optimizes the problem in O(n2 log n) iterations.

Our analysis gives insights into the working principles of parallel tempering and reveals why the interplay of two states running at different temperatures is crucial to optimize the benchmark efficiently. Our results demonstrate the power of parallel tempering and lend credence to its design.

1.1 RelatedWork

Several runtime analyses have been conducted for the Metropolis algorithm, which maintains a single state and a single (fixed) temperature. Sasaki and Hajek [1988] showed that the Metropolis algorithm can approximate a maximum matching on graphs with n vertices within a (1 - ϵ) factor in expected time at most 6 ϵ2 22/ϵn3+2/ϵ . Jerrum and Sorkin [1998] proved that the Metropolis algorithm computes an optimal graph bisection in quadratic time.

Allowing the system temperature to vary with time results in the simulated annealing algorithm [Kirkpatrick et al., 1983; ˇCern´y, 1985], and a variety of runtime analyses investigate the utility of allowing temperature fluctuations in different search landscapes. Wegener [2005] gave the first example where simulated annealing outperforms the Metropolis algorithm with any fixed temperature. Our negative results for the Metropolis algorithm are inspired by his work. Jansen and Wegener [2007] proved that the Metropolis algorithm solves the ONEMAX benchmark function in O(n log n) time with a suitably low temperature setting. They also compared Metropolis, simulated annealing and a simple evolutionary algorithm called the (1+1) EA on a number of benchmark functions and proved the existence of functions that Metropolis solves in polynomial time, but the (1+1) EA needs at least exponential time. In contrast, Lissovoi, Oliveto and Warwicker [2023] proved that the Metropolis algorithm fails to efficiently solve the simple multi-modal CLIFF function, while Doerr et al. [2023] proved that the (1+1) EA is comparatively faster than Metropolis on CLIFF.

Oliveto et al. [2018] investigated the Metropolis algorithm and a related model from population genetics called SSWM and proved that they are able to optimize the rugged VALLEY function more efficiently than simple evolutionary algorithms by exploiting their ability to accept worsening moves. Wang, Zheng and Doerr [2024] proved that the Metropolis algorithm solves the DECEPTIVELEADINGBLOCKS problem in quadratic time, while comparable algorithms that do not accept inferior moves require at least cubic time. Zheng et al. [2024] analyzed the Metropolis algorithm in the context of multi-objective optimization, considering a bi-objective variant of DECEPTIVELEADINGBLOCKS.

Migration is a key concept in parallel evolutionary algorithms such as island models, in which several populations evolve separately, but coordinate their search processes by migrating solutions between them [Luque and Alba, 2011; Alba et al., 2013; Sudholt, 2015]. For these algorithms, a theoretical foundation has been established, cf. [Sudholt, 2015] and the references therein. Most relevant to our work is [L¨assig and Sudholt, 2013], which showed that migration can yield exponential speedups for a constructed benchmark problem. In addition, heterogeneous island models, in which islands solve different aspects of a problem, were proven to be effective for the SETCOVER problem [Mambrini et al., 2012].

2 Preliminaries

A pseudo-Boolean function is a real-valued mapping f : {0, 1}n → R. A pseudo-Boolean function f is called quadratic when it can be written in the form of a quadratic multi-linear polynomial:

In quadratic unconstrained binary optimization (QUBO), given a quadratic pseudo-Boolean function f, the objective is to solve

We introduce a simple QUBO benchmark that allows for tuning the height of energy barriers, but is elementary enough so that the behavior of algorithms is better clarified.

Definition 1. Let n be a positive even integer and δL < δG be the energy barrier parameters with 1 ≤ δL = O(1). We define the quadratic pseudo-Boolean function RUGGEDLANDSCAPEn,δL,δG as follows. The quadratic terms following (1) are

The linear terms are ci = δL for i > 0 and the constant term is c0 = n2 (δG - δL).

The class RUGGEDLANDSCAPEn,δL,δG contains functions over length-n binary vectors where quadratic interactions are "partitioned" into n/2 disjoint bit-pairs (x1, x2), (x3, x4), . . . , (xn-1, xn), and no interactions occur between pairs. Moreover, bits within a pair exhibit socalled reciprocal sign epistasis [Poelwijk et al., 2011] in the sense that the contribution to the energy from a (0, 0)-pair is c0/(n/2) = δG - δL, the contribution from a (1, 0) or (0, 1) pair is c0/(n/2) + δL = δG, while the contribution from a (1, 1) pair is c0/(n/2) + 2δL - (δG + δL) = 0.

Let x ∈ {0, 1}n be a binary vector and π = (xi, xi+1) a bit pair such that n > i ≡ 1 (mod 2). When xi = xi+1 = 0, we say π is a metastable pair. When xi + xi+1 = 1, we say π is an unstable pair. Finally, when xi = xi+1 = 1, we say π is a stable pair.

To propose new states, the algorithms we study sample from the Hamming neighborhood, i.e., one bit is flipped at a time. This induces a topology on the energy function that makes it possible to define structural properties such as local minima and energy barriers. A state x is a local minimum with respect to Hamming neighbors if and only if it has no unstable pairs. If x is a local minimum with metastable pairs, it is a metastable state (non-global local minimum). If a local minimum has no metastable pairs, it is a global minimum with zero energy. The landscape of RUGGEDLANDSCAPE thus contains 2n/2 - 1 metastable states separated by states with unstable pairs. In order to escape a metastable state, a process must pass through an unstable state with higher energy. Let x be a metastable state and y a Hamming neighbor of y. If y has one fewer metastable pair than x, then f(y) = f(x)+δL > f(x). Similarly, if y has one fewer stable pair, then f(y) = f(x) + δG > f(x). In any case, to escape a metastable state, the processes must transition through an "energy barrier" of height at least δL.

Our proofs makes use of the multiplicative drift theorem, which connects the expected change in a random process to the first hitting time of that process to a certain state of interest. The multiplicative drift theorem can be stated as follows.

Theorem 1 (Doerr and Goldberg, 2013). Let (Xt)t∈N be a non-negative stochastic process over a state space S ⊆ {0}∪ [smin,∞) where smin > 0, and let T := inf{t ≥ 0 | Xt = 0}. Suppose that X0 = s0 and there exists a δ > 0 such that for all s ∈ S \ {0} and all t ≥ 0

Then the following tail bound holds for the first hitting time of the process to zero:

We also use the negative drift theorem with self-loops presented in [Rowe and Sudholt, 2014] (an extension of the negative drift theorem [Oliveto and Witt, 2011; Oliveto and Witt, 2012] to stochastic processes with large self-loop probabilities). It is stated in the following for the sake of completeness.

Theorem 2 (Rowe and Sudholt, 2014). Consider a Markov process X0,X1, . . . on {0, . . . ,m} with transition probabilities pi,j and suppose there exist integers a, b with 0 < a < b ≤ m and ε > 0 such that for all a ≤ k ≤ b the drift towards 0 is

where pk,k is the self-loop probability at state k. Further assume there exist constants r, δ > 0 (i. e. they are independent of m) such that for all k ≥ 1 and all d ≥ 1

Let T be the first hitting time of a state at most a, starting from X0 ≥ b. Let ℓ = b - a. Then there is a constant c > 0 such that

3 MCMCAlgorithms for Optimization

Markov chain Monte Carlo (MCMC) sampling is a technique for exploring a state space X by defining a Markov process on X whose stationary distribution, under certain balance conditions, converges in the limit to a target distribution. Many of these techniques are useful in the context of optimization when considered as randomized search heuristics in which the state space is the set of candidate solutions and the target distribution is proportional to the objective function.

One of the first such methods is the Metropolis algorithm [Metropolis et al., 1953; Hastings, 1970]. This algorithm maintains a single state, and iteratively performs a biased random walk. In each iteration, a random neighbor y of the current state x is produced. If f(y) ≤ f(x), then the state y is accepted as the new current state. Otherwise, y is accepted with probability e-∆f/T where ∆f := f(y) - f(x) is the difference in objective function value and T is a temperature parameter that governs how likely a move to a higher energy state is. The Metropolis algorithm is listed in Algorithm 1.

The temperature T is one of the most important parameters in Metropolis and similar MCMC sampling techniques. The process can permeate an energy barrier with height at most T, but it cannot easily escape from a valley with depth greater than T . Furthermore, the process cannot easily follow gradients with magnitude smaller than T , so the precision of the result depends heavily on how the temperature is set.

One way to address this is to also adapt the temperature. Simulated annealing [Kirkpatrick et al., 1983; ˇCern´y, 1985] takes a "cooling" schedule and decreases the temperature monotonically during the run of the algorithm. Simulated tempering treats the temperature as another state to sample [Marinari and Parisi, 1992]. Given a set of temperatures T1 < · · · < Tκ, a single state (x, i) is maintained. The x update is a Metropolis update with temperature Ti, and i → i±1 is also updated stochastically. The idea is to use a wide range of temperatures without converging too soon into a low temperature suboptimal configuration.

Rather than maintaining a single state that explores the temperature space, parallel tempering also takes a sequence of temperatures T1 < · · · < Tκ, but additionally maintains a pool of κ parallel states x1, . . . , xκ, each performing a Metropolis update with its respective temperature value [Swendsen and Wang, 1986; Geyer, 1991]. Periodically, neighboring states are exchanged using a Metropolis criterion based on the temperature ratios. In particular, an exchange between xi and xi+1 is accepted with probability

The pseudocode for parallel tempering is listed in Algorithm 2. Varying approaches exist for implementing the frequency of exchanges, e.g., by performing an exchange every M steps, or in each step but only with a fixed probability q.

In the classical approach, the exchange operation is symmetric so no configuration is ever lost during an exchange. In particular if xi exchanges with xj in a time step, then xi and xj remain in the ensemble, but in different processes with different temperatures. In this paper, we employ a slight modification to this approach by allowing for asymmetric exchanges. This means it is possible for a state to be lost in an exchange if it is overwritten by an incoming state. For that reason, we call the procedure a migration rather than an exchange. This gives rise to Algorithm 3 in which, during a migration step, process i + 1 migrates to process i with probability

and process i migrates to process i + 1 with probability

As we will see later, this adaptation allows the algorithm to carefully leverage the different temperatures without losing valuable state information to an unlucky exchange. Moreover, a symmetric exchange occurs with the joint probability of both events, the product of which recovers the classical exchange probability in Equation (4).

Algorithm 1: Metropolis Algorithm

Input: state x Input: temperature T > 0 1 repeat 2 y ←neighbor(x); 3 p ← exp ( -(f(y)-f(x)) T ); 4 if random() ≤ p then x ← y; 5 until stop;

Algorithm 2: Parallel Tempering

Input: κ states x1, . . . , xκ Input: κ temperatures 0 < T1 < · · · < Tκ 1 repeat /* Perform Metropolis updates */ 2 for i ∈ {1, . . . , κ} do 3 y ← neighbor(xi); 4 p ← exp ( -(f (y)-f (xi)) Ti ); 5 if random() ≤ p then xi ← y; /* Pair exchange */ 6 periodically 7 Choose i ∈ {1, . . . , κ - 1} u.a.r.; 8 p ← exp (-(f(xi+1) - f(xi))( 1 Ti - 1 Ti+1 )); 9 if random() ≤ p then 10 (xi, xi+1) ← (xi+1, xi); 11 until stop;

4 Runtime Analysis

We now analyze the runtime of the Metropolis algorithm and parallel tempering with asymmetric exchanges on RUGGEDLANDSCAPE. We assume the initial state is selected uniformly at random from {0, 1}n , and the periodic exchange step occurs with constant probability 0 < q < 1. Given a temperature T > 0, we use the notation β = β(T ) = 1/T to denote the inverse temperature.

Algorithm 3: Parallel Tempering (Asymmetric Exchange)

Input: κ states x1, . . . , xκ Input: κ temperatures 0 < T1 < · · · < Tκ Input: migration probability q 1 repeat /* Perform Metropolis updates */ 2 for i ∈ {1, . . . , κ} do 3 y ← neighbor(xi); 4 p ← exp ( -(f (y)-f (xi)) Ti ); 5 if random() ≤ p then xi ← y; /* Asymmetric pair exchange */ 6 if random() ≤ q then 7 Choose i ∈ {1, . . . , κ - 1} u.a.r.; 8 (ˆxi, ˆxi+1) ← (xi, xi+1); 9 pi ← exp ( -(f(xi+1)-f(xi)) Ti ); 10 pi+1 ← exp ( -(f (xi)-f (xi+1)) Ti+1 ); 11 if random() ≤ pi then ˆxi ← xi+1; 12 if random() ≤ pi+1 then ˆxi+1 ← xi; 13 (xi, xi+1) ← (ˆxi, ˆxi+1); 14 until stop;

4.1 Metropolis Algorithm

Theorem 3. Consider the Metropolis algorithm (Algorithm 1) optimizing RUGGEDLANDSCAPE starting from a random initial state with inverse temperature β ≥ c δL lnn for any constant c > 0. With probability 1 - e-Ω(n) , the algorithm does not find the optimum before nc+1 /(2e) steps. Furthermore, when β = ω(logn), the algorithm requires nω(1) steps with probability 1 - e-Ω(n) .

Proof. During the run of the algorithm, any metastable pair must be transformed into an unstable pair at least once before it can become a stable pair. The probability that a metastable pair becomes unstable in any iteration is 2e-βδL /n, since the pair is changed with probability 2/n and the corresponding energy increase is accepted with probability e-βδL .

Let E be the event that the random initial state has at least n/16 metastable pairs. By a Chernoff bound [Mitzenmacher and Upfal, 2017, Theorem 4.5], Pr(E) ≥ 1-e-n/64 . Conditioning on this event, the optimum is not found until at least n/16 metastable pairs are removed. Applying a union bound, for any t > 0, the probability that a metastable pair becomes unstable before t iterations is at most 2te-βδL /n. Thus with probability at most (2te-βδL )n/16 , all of the metastable states have transitioned at least once in t iterations, which is a necessary condition for hitting the optimum. For t < nc+1 /(2e) this is at most exp (- n 16 (βδL - c ln n + 1)) ≤ e-n/16 when β ≥ c δL lnn. When β = ω(log n), for any t < 12 (n βδL 2 ln n ) = nω(1) , this probability is at most e-βδLn/32 = e-Ω(n) .

In either case, by the law of total probability, all of the metastable pairs have transitioned within the first t iterations with probability at most e-Ω(n) + (1 - Pr(E)) = e-Ω(n) ,

which proves the claim.

The following theorem shows that high temperatures lead to chaotic behavior and exponential optimization times.

Theorem 4. Consider the Metropolis algorithm (Algorithm 1) optimizing RUGGEDLANDSCAPE starting from a random initial state with inverse temperature β ≤ c/(δG - δL) for an arbitrary constant c > 0. The number of steps until the algorithm finds the optimum is at least 2dn with probability 1-2-Ω(n) for some constant d > 0 that may depend on c.

Proof. Theorem 2 will be applied to the count of bit pairs that are not in a stable state when this count is in the interval [0, αn] for a small constant 0 < α < 3/4 to be determined later. Note that the optimum is only found when the number of bit pairs that are not stable has reached 0. Again applying a Chernoff bound, the Metropolis algorithm starts with at least αn non-optimal bit pairs with probability 1 - 2-Ω(n) . We assume in the following that this happens; the claimed probability bound then follows from a union bound of the failure probability 2-Ω(n) and a failure probability 2-Ω(n) which will result from an application of the drift theorem.

Note that every unstable pair (i.e., of the form (0, 1) or (1, 0)) will eventually transition into a metastability (0, 0) or a stability (1, 1), and this transition will always happen with probability one since unstable pairs contribute the highest energy. Conditioning on this transition, the probability of reaching a metastability is 1/2 and the probability of reaching a stability is 1/2.

We consider a modified process to simplify the analysis: we assume that every unstable pair immediately transitions to a metastability or a stability, with the aforementioned conditional probabilities. One way to envision this is that for a particular pair, the random decision about its next transition is made in advance and stored. When the time comes that an unstable pair is mutated, we use the stored transition. Since the algorithm flips only one bit at a time, all pairs can be considered independently, which justifies this early-decision mechanism. Intuitively, this speeds up a transition that would be made in the future. As a consequence, in the modified process, all bit pairs are always metastabilities or stabilities. Since we eliminate the time a bit pair spends in the unstable state, we claim that the optimization time of the modified process is stochastically dominated by that of the original process (though we do not give a formal proof).

Fix one bit pair. If that bit pair is in a metastable state, the probability that it will mutate into a stable state in one step is (1/n)e-βδL (since the original process moves to an unstable state with probability (2/n)e-βδL and subsequently moves to a stable state with probability 1/2). Similarly, a pair in a stable state will transition to a metastable state with probability (1/n)e-βδG .

Consider the stochastic process (Xt)t∈N that counts the number of pairs that are not in a stable state at time t and denote as Et the event {Xt ≤ αn}. The self-loop probability pk,k in the statement of Theorem 2 is the probability that the state count remains unchanged in a single iteration. Thus we have Pr (Xt+1 ̸= Xt | Xt, Et) = 1-pk,k, and we bound this as follows.

since δG > δL. The drift of the Xt process below αn can be bounded as

Choose α small enough so that 2 (α(1 + e-c ) - e-c /2) < 0. By (7), we have -e-βδL < -2(1 - pk,k), so

satisfying the the first condition (2) of Theorem 2 for ε := -2 (α(1 + e-c ) - e-c /2). The second condition (3) follows trivially since the number of non-optimal components only changes by at most 1 in any step. Invoking the negative drift theorem implies the claim.

We now show that the low temperature and high temperature regime can overlap, and the performance of the Metropolis algorithm can degrade significantly as a function of energy barrier parameters in RUGGEDLANDSCAPE.

Corollary 1. Let δG = δL+ln-1 n. For any temperature setting, the Metropolis algorithm requires Ω(nδL+1 ) steps to optimize RUGGEDLANDSCAPE with probability exponentially close to one.

Proof. For all β ≤ lnn = 1/(δG - δL), Theorem 4 with c = 1 yields an exponential optimization time for Algorithm 1 with probability 1 - 2-Ω(n) . For all β ≥ ln n, Theorem 3 with c = δL yields an optimization time of Ω(nδL+1 ) with probability 1 - e-Ω(n) .

4.2 Parallel Tempering

In this section we show that parallel tempering can be highly effective, and already two parallel processes suffice to demonstrate an advantage over a single-temperature process. Our setting employs one process at a low temperature 1/β1 and one at a high temperature 1/β2. The low-temperature process ensures stability, as it is unlikely to accept worsenings, whereas the high-temperature process facilitates exploration by traversing energy barriers. We show that improvements are discovered through exploration in the hightemperature process and transferred to the low-temperature process, where they are permanently retained. This guarantees global monotonicity.

Theorem 5. Consider a run of a κ = 2 state parallel tempering process with asymmetric exchanges (Algorithm 3) and a constant migration probability q ∈ (0, 1) on the quadratic RUGGEDLANDSCAPE function with any inverse temperatures β1 ≥ (2 + ε) lnn and β2 = O(1) for any constant ε > 0. The algorithm optimizes the function in O(n2 logn) steps with probability 1 - o(1).

Proof. We consider the parallel stochastic processes (X(1) t ,X(2) t )t∈N on R that correspond to the value of f(x1) and f(x2) in the t-th step, respectively. We begin by arguing that for all times t < n2+ε/2 , X(1) t ≤ X(1) t-1 with probability 1 - e-Ω(β1) . At any time step t, X(1) t > X(1) t-1 only if a higher energy state was accepted in line 4 of Algorithm 3, or if a higher energy state was swapped in from x2 in line 10 (or both). This happens with probability at most 2e-β1δL . Taking a union bound, the probability that X(1) t does not increase during any of the first n2+ε/2 steps is at least 1 - 2n2+ε/2 e-β1δL = 1 - e-Ω(β1) since we have assumed β1 ≥ (2 + ε) lnn. For the remainder of the proof, we condition on the event that the energy does not increase at state x1 before n2+ε/2 = ω(n2 log n) steps.

We now argue that Pr(X(1) t+2 < X(1) t ) ≥ 2q2 (1-q)e-β2δL n2 . We make a case distinction on the relative values of X(1) t and X(2) t . If X(2) t < X(1) t , a sufficient condition for X(1) t+1 < X(1) t is that the lower energy state from x2 was swapped in to x1 at line 10. This occurs with probability q, and since we assume that the energy does not increase in state x1, the new lower energy is preserved, and the claim also holds for X(1) t+2 .

The case X(2) t ≥ X(1) t is slightly more involved. A sufficient condition to decrease the energy at x1 in three steps is for exactly one metastable or unstable pair in xt to transform into a stable pair. This requires relying on the higher temperature process x2 to traverse an energy barrier. Fix an index pair π that corresponds to a pair that is not stable at time t in x1. If π is unstable, then it becomes stable at time t + 1 with probability (1 - q)/n, which results in X(1) t+1 < X(1) t and the claim holds as before.

Otherwise, assume that π is a metastable pair. In order for π to become stable, it suffices that three consecutive events occur. Let E1 denote the event that x1 is swapped into x2 at line 11 at time t. Let E2 denote the event that π transitions frommetastable to unstable at line 4 in x2 at time t+1, and is not overwritten by an exchange during that iteration. Finally, let E3 denote the event that π transitions from unstable to stable during iteration t + 2 and is subsequently copied back to x1 by an exchange (or that the energy at x1 has already decreased, in which case we are done).

It remains to bound the probability of this event sequence. Since X(2) t ≥X(t) t , we have p2 ≥ 1 at line 11, and it suffices that the exchange occurs, so we have Pr(E1) = q. For E2, the Metropolis update at time t + 1 for x2 must flip one of the zero bits in π to a one with probability 2/n, after which the resulting state must be accepted, which happens with probability e-β2δL . This update is not overwritten by a subsequent exchange during step t + 1 with probability at least 1 - q. Thus we have Pr(E2 | E1) ≥ (1 - q)2e-β2δL /n. For E3, the probability that π is transformed into a stable pair in the update step is 1/n, and this proposal results in a strict energy decrease, which is accepted with probability 1. The resulting state is copied back to x1 in the exchange step of iteration t with probability q, yielding Pr(E3 | E2 ∩ E1) ≥ q/n. Event E1 ∩ E2 ∩ E3 is sufficient for {X(1) t+2 < X(1) t } so we have Pr(X(1) t+2 < X(1) t ) ≥ 2q2 (1-q)e-β2δL n2 as claimed.

To complete the proof, consider the stochastic process (Yt)t∈N where Yt counts the pairs in x1 that are not stable at time t. Since the preceding arguments hold independently for all such pairs, it follows that

We are interested in the first hitting time of (Yt)t∈N to zero, T := inf{t ≥ 0 | Yt = 0}, which corresponds to the first time that f(x1) = 0, i.e., the optimal solution. Clearly Y0 ≤ n/2 and Theorem 1 can be adapted to supply a tail bound on T by setting s0 and smin accordingly. Note that the timestep lag in the drift condition in (8) only slows the drift by a factor of two. We therefore conclude that for all r > 0 the probability that the process requires more than n2 q2 (1-q) eβ2δL (ln(n/2) + r) steps until the process hits the optimal solution is at most e-r . Since this outcome is conditioned on the event that the energy of x1 never increases during a polynomial number of steps, the final tail bound is 1 - e-Ω(β1) - e-r . As β1 ≥ (2 + ε) lnn, setting r = (2 + ε) ln n yields the claimed result.

5 Experiments

To complement our theoretical results, we report on the results of a selection of experimental studies. In particular, we are interested in the behavior of the algorithms on RUGGEDLANDSCAPE n,δL,δG for concrete parameter settings. We fix the problem size at n = 100. Our proofs show that when δG is only marginally larger than δL, the Metropolis algorithm fails to distinguish metastable and stable states; thus we set δL = 100 and δG = 105. We consider a sweep of 5000 temperature values T ∈ {0.01, 0.02, . . . , 50.00} and run the Metropolis algorithm 100 times with the corresponding temperature setting for 106 iterations each. The mean final energy discovered at each temperature value is plotted in Figure 1. The minimum energy found over all 500000 runs was 35 (7 metastable states) discovered in two individual runs at temperatures 12.87 and 18.32. In general, runs with an energy less than 50 occurred only for T ∈ [10, 30].

To study a typical Metropolis run at a promising temperature, we track the count of metastable, unstable and stable pairs at temperature T = 12 in Figure 2. At the beginning of the run, each pair-type occupies roughly a 1/4 fraction of the total count of n/2 = 50 pairs, as expected from a uniform random initial state. In the first 100 iterations, the two unstable pair-types (0, 1) and (1, 0) quickly vanish. Subsequently, the count of metastable and stable pairs transition into an equilibrium around 2/5 and 3/5, respectively.

We repeat this for a two-state parallel tempering process with temperatures T1 = 1/1000 and T2 = 100 and migration probability q = 0.1 and plot the results in Figure 3. We observe similar behavior early in the process, whereas the number of metastable pairs vanish steadily later in the process. Finally, we compare the energy as a function of time step for both runs in Figure 4. We observe that the Metropolis algorithm stagnates at energy values around 100, whereas parallel tempering continues to make progress and eventually finds the global optimum.

6 Conclusions

We gave the first proof that parallel tempering can be more effective than a single instance of the Metropolis algorithm in the context of optimization. Our proposed QUBO class RUGGEDLANDSCAPE is highly rugged as it contains an exponential number of local optima. We showed that a singletemperature Metropolis algorithm fails if the temperature is too small or too large. In the former case, the algorithm is unable to escape from local optima, and in the latter case, the search dynamics are too chaotic to approach the global optimum. We furthermore showed that the temperature regimes of these two failure modes can overlap, precluding any ideal temperature setting that admits efficient behavior.

In contrast, we proved that a parallel tempering algorithm with one low-temperature process and one high-temperature Time process is effective on RUGGEDLANDSCAPEn,δL,δG . The low-temperature process typically retains the best solution found so far. Improvements can be discovered by migrating this solution to the high-temperature process, which can traverse an energy barrier and find improvements that are migrated back to the low-temperature system. Hence, the interplay between these parallel processes of different temperatures is crucial for optimization. Our work lends credence to the algorithmic design of parallel tempering and takes a first step towards establishing a rigorous theoretical foundation that explains when and why parallel tempering is efficient, and how it can be used most effectively.

Figure 4: Energy over time for both algorithms on RUGGEDLANDSCAPE.

We also introduced a minor modification to the algorithm in which exchanges are allowed to be asymmetric. This adjustment may be particularly useful in cases where the algorithm is applied to solving optimization problems, as it allows the process a reasonable chance to maintain promising candidate solutions without the risk of being overwritten during an exchange. Two natural future directions are to extend the analysis to many parallel states, and to investigate other energy landscapes, such as those induced by combinatorial optimization or multi-objective problems.

Acknowledgments

This research was funded by NSF grant 2144080.