# A Stochastic Tunneling Approach for Global Minimization of Complex Potential Energy Landscapes

\[

## Abstract

We investigate a novel stochastic technique for the global optimization of complex potential energy surfaces (PES) that avoids the freezing problem of simulated annealing by allowing the dynamical process to tunnel energetically inaccessible regions of the PES by way of a dynamically adjusted nonlinear transformation of the original PES. We demonstrate the success of this approach, which is characterized by a single adjustable parameter, for three generic hard minimization problems.   

PACS: 02.70.Pn, 02.70.Lq, 02.50.Ey, 02.70.Ln

\]

The development of methods that efficiently determine the global minima of complex and rugged energy landscapes remains a challenging problem with applications in many scientific and technological areas. In particular for NP-hard problems stochastic methods offer an acceptable compromise between the reliability of the method and its computational cost. Branch-and-bound techniques offer stringent error estimates but scale exponentially in their computational effort. In many stochastic approaches the computational cost to determine the global minimum with a given probability grows only as a power-law with the number of variables .

In such techniques the global minimization is performed through the simulation of a dynamical process for a “particle” on the multi-dimensional potential energy surface. Widely used is the simulated annealing (SA) technique  where the PES is explored in a series of Monte-Carlo simulations at successively decreasing temperatures. Its success depends often strongly on the choice of the cooling schedule, yet even the simplest geometric cooling schedule is characterized by three parameters (starting temperature, cooling rate and number of cooling steps) that must be optimized to obtain adequate results. For many difficult problems with rugged energy landscapes SA suffers from the notorious “freezing” problem, because the escape rate from local minima diverges with decreasing temperature. To ameliorate this problem many variants of the original algorithm have been proposed. Unfortunately these proposals often increase the number of parameters even further, which complicates their application for practical problems.

In this letter we investigate the stochastic tunneling method, a generic physically motivated generalization of SA. This approach circumvents the “freezing” problem, while reducing the number of problem dependent parameters to one. In this investigation we demonstrate the success of this approach for three hard minimization problems: the Coulomb spin-glass (CSG), the traveling salesman problem (TSP) and the determination of low autocorrelation binary sequences(LABS) in comparison with other techniques.

Method: The freezing problem in stochastic minimization methods arises when the energy difference between “adjacent” local minima on the PES is much smaller than the energy of intervening transition states separating them. As an example consider the dynamics on the model potential in Figure (1)(a). At high temperatures a particle can still cross the barriers, but not differentiate between the wells. As the temperature drops, the particle will get eventually trapped with almost equal probability in any of the wells, failing to resolve the energy difference between them. The physical idea behind the stochastic tunneling method is to allow the particle to “tunnel” forbidden regions of the PES, once it has been determined that they are irrelevant for the low-energy properties of the problem. This can be accomplished by applying a non-linear transformation to the PES:

$$f_{\mathrm{STUN}}(x)=1\mathrm{exp}\left[\gamma (f(x)f_0)\right]$$
(1)   
where $`f_0`$ is the lowest minimum encountered by the dynamical process so far (see Figure (1)(b) + (c)). The effective potential preserves the locations of all minima, but maps the entire energy space from $`f_0`$ to the maximum of the potential onto the interval $`[0,1]`$. At a given finite temperature of O(1), the dynamical process can therefore pass through energy barriers of arbitrary height, while the low energy-region is still resolved well. The degree of steepness of the cutoff of the high-energy regions is controlled by the tunneling parameter $`\gamma >0`$. Continously adjusting the reference energy $`f_0`$ to the best energy found so far, successively eliminates irrelevant features of the PES that would trap the dynamical process.

To illustrate the physical content of the transformation we consider a Monte-Carlo (MC) process at some fixed inverse temperature $`\beta `$ on the STUN-PES. A MC-step from $`x_1`$ to $`x_2`$ with $`\mathrm{\Delta }=f(x_2)f(x_1)`$ is accepted with probability $`\stackrel{~}{w}_{12}=\mathrm{exp}\left[\beta (f_{\mathrm{STUN}}(x_2)f_{\mathrm{STUN}}(x_1))\right]`$. In the limit $`\gamma \mathrm{\Delta }1`$ this reduces to $`\stackrel{~}{w}_{12}\mathrm{exp}(\stackrel{~}{\beta }\mathrm{\Delta })`$ with an effective, energy dependent temperature $`\stackrel{~}{\beta }=\beta \gamma e^{\gamma (f_0f(x_1))}\beta \gamma .`$ The dynamical process on the STUN potential energy surface with fixed temperature can thus be interpreted as an MC process with an energy dependent temperature on the original PES. In the latter process the temperature rises rapidly when the local energy is larger than $`f_0`$ and the particle diffuses (or tunnels) freely through potential barriers of arbitrary height. As better and better minima are found, ever larger portions of the high-energy part of the PES are flattened out. In analogy to the SA approach this behavior can be interpreted as a self-adjusting cooling schedule that is optimized as the simulation proceeds.

Since the transformation in equation (1) is bounded, it is possible to further simplify the method: On the fixed energy-scale of the effective potential one can distinguish between phases corresponding to a local search and “tunneling” phases by comparing $`f_{\mathrm{STUN}}`$ with some fixed, problem independent predefined threshold $`f_\mathrm{t}`$ (see Fig. 1(c)). For the success of the method it is essential that the minimization process spends some time tunne-

ling and some time searching at any stage of the minimization process. We therefore adjust the parameter $`\beta `$ accordingly during the simulation: If a short-time moving average of $`f_{\mathrm{STUN}}`$ exceeds the threshold $`f_{\mathrm{thresh}}0.03`$, $`\beta `$ is reduced by some fixed factor, otherwise it is increased. Following this prescription the method is characterized by the single problem-dependent parameter ($`\gamma `$).

Applications: In order to test the performance of this algorithm we have investigated three families of complicated NP-hard minimization problems. For each problem we have determined either the exact ground-state energy or a good estimate thereof. We computed the average error of the various optimization methods as a function of the computational effort to determine the computational effort required to reach a prescribed accuracy. For the applications presented here we have fixed the functional form of the transformation and the “cooling schedule” for $`\beta `$ in order to demonstrate that these choices are sufficient to obtain adequate results. Obviously this does not guarantee that these choices are optimal.

(CSG) The determination of low-energy configurations of glassy PES is a notoriously difficult problem. We have verified by direct comparison that the method converges quickly to the exact ground states  for two-dimensional short-range Ising spin-glasses of linear dimension $`10`$ to $`30`$ with either discrete or Gaussian distributions of the coupling parameters. Next we turned to the more demanding problem of the Coulomb spin-glass, where classical charges $`\{s_i\}`$ with $`s_i=\pm 1`$ are placed on fixed randomly chosen locations within the unit cube. The energy of the system

$$E(\{s_i\})=\underset{ij}{\overset{N}{}}\frac{s_is_j}{|\stackrel{}{r}_i\stackrel{}{r}_j|},$$
(2)   
is minimized as a function of the distribution of the $`\{s_i\}`$.

The results of grand-canonical simulations for ten replicas of $`N=100`$ and $`N=500`$ charges are shown in Figure 2. We first conducted twenty very long STUN runs for each replica to determine upper bounds for the true ground-state energy. For the same charge distributions we then averaged the error of STUN, SA, parallel tempering (PT)  and simulated tempering(ST)  for twenty runs per replica as function of the numerical effort. We found that the average STUN energy converged in $`10^6`$ MC-steps to within 1% of the estimated true ground-state energy. Over two decades of the numerical effort we found a consistent significant advantage of the proposed method over the SA approach. Fitting the curves in the figure with a power-law dependence we estimate that STUN is two orders of magnitude more efficient than SA.

We found no consistent ranking of ST and PT relative to SA for the two system sizes considered. Both methods offer alternative routes to overcome the freezing problem in SA. In PT the configurations of concurrent simulations at a variety of temperatures are occasionally exchanged. In ST only a single simulation is undertaken, but its temperature is considered to be a dynamical variable. Temperature and configuration are distributed according to: $`p(s,T)=e^{E(s)/Tg(T)}`$ and the weights $`g(T)`$ are optimized for a discretized temperature distribution, such that all temperatures are visited with equal probability. In both methods, a configuration can escape a local minimum when the instantaneous temperature is increased. The choice of the temperature set (along with values for $`g(T)`$) is system dependent and must be optimized much like the annealing schedule in SA. In accordance with other studies our results indicate that ST performs significantly better than SA for long simulation times. PT was successful only for the larger system (N=500), where it reached the same accuracy as STUN for $`10^6`$ steps. STUN converged faster than any of the competing methods, but showed a tendency to level off at high accuracy. In the limit of large computational its accuracy was matched by ST for N=100 and PT for N=500.

(TSP) The traveling salesman problem is another ubiquitous NP-hard minimization problem. We have investigated the problem in its simplest incarnation: i.e. as a minimization of the euclidian distance along a closed path of N cities. Using long-range updates, i.e. the reversal and exchange of paths of arbitrarily length, we found that both SA and STUN perform about equally well and reach the global optimum for $`N=20,50`$ and $`100`$ very quickly (see right side of Table (I)).

Nevertheless it is instructive to analyze this model somewhat further as it provides insight into the interplay of move-construction and complexity of the minimization problem. The unconstrained TSP is a rare instance among NP-hard minimization problems, where it is possible to construct efficient “long-range” hops on the PES. In most practical applications of minimization problems related to the TSP, the construction of global moves is severely complicated by the existence of “hard constraints” on the routes taken. For such problems, as well as the other examples reported here, the alteration of just a few variables of the configurations leads to unacceptably high energies in almost all cases. As a result, the construction of global moves is not an efficient way to facilitate the escape from local minima. When only local moves, i.e. transpositions of two adjacent cities, are considered high barriers that were circumvented in the presence of global moves hamper the progress of SA. The results on the left side of Table (I) demonstrate that in this scenario SA performs significantly worse than STUN.

(LABS) Finally we turn to the construction of low-autocorrelation binary sequences. The model can be cast as a ground-state problem for a one-dimensional classical spin-1/2 chain with long-range four-point interactions

$$E=\frac{1}{N}\underset{k=1}{\overset{N1}{}}\left[\underset{j=1}{\overset{Nk}{}}s_js_{j+k}\right]^2$$
(3)   
and is one of the hardest discrete minimization problems known. Even highly sophisticated and specialized optimization algorithms have failed to find configurations anywhere near (within 20%) the ground-state energy that can be extrapolated from exact enumeration studies for small systems ($`N<50`$). The reason for this difficulty has been attributed to the “golf-course” character of the energy landscape and there is convincing evidence that SA will fail to converge to the ground-state energy even in the limit of adiabatic cooling. The situation is significantly improved if the original potential energy surface is replaced by a piecewise constant energy surface that is obtained by a local minimization of the original PES at each point. Obviously the latter surface preserves all ground-state configurations and energies of the original PES, but eliminates many “plateaus” of the “golf-course” landscape. Using the modified energy surface we are able to compare SA to STUN, since SA can now determine the ground state energy of medium

size systems (N=49) with a large, but finite computational effort. Table II summarizes the results for the average error of 20 SA and STUN runs for system sizes N=49 and N=101 as a function of the computational effort. In direct comparison we find that STUN is two orders of magnitude more efficient than SA. Both methods are at least a dozen orders of magnitude more efficient than SA on the original PES.

Discussion: Using three NP-hard minimization problems with high-barriers separating local minima we have demonstrated that the stochastic tunneling approach offers a reliable, generic and efficient route for the determination of low-energy configurations. One chief advantage of the method lies in the fact that only a single parameter must be adjusted to adapt it for a specific problem. One of the drawbacks of STUN is that in contrast to e.g. PT, no thermodynamic expectation values for the system can be obtained from the simulation. Secondly, because the non-linear transformation will map any unbounded PES onto an interval bounded from above, the dynamical process in STUN will experience “tunneling” phases at any finite temperature. For PES that do not contain high barriers, or in the presence of efficient global moves that circumvent such barriers, STUN may therefore be less efficient than competing methods. In many realistic optimization problems where the construction of global moves is exceedingly difficult or very expensive the tunneling approach can ameliorate the difficulties associated with the existence of high energy barriers that separate local minima of the PES. 

We gratefully acknowledge stimulating discussions with C. Gros and U. Hansmann.