EVALUATION OF SIMULATED TEMPERING FOR OPTIMIZATION PROBLEMS THOMAS ROSE* Research Apprentice Northeast Parallel Architectures Center Syracuse University Mathematics/Computer Science undergraduate University of Texas of the Permian Basin PAUL CODDINGTON Research Scientist Syracuse University ENZO MARINARI Professor Syracuse University ABSTRACT The application of massively parallel computers in simulating disordered physical systems and in the search for global minima of discrete optimization problems is addressed. Simulated tempering is a new algorithm which is presented as an alternative to the Monte Carlo procedure known as simulated annealing. A sequential version of this algorithm has been shown to greatly improve the simulation of a specific spin model when compared to conventional Monte Carlo techniques. In the present work, an efficient parallel program which performs simulated tempering on Ising spin models is produced and implemented on a MasPar system, using the MasPar Programming Language. The performance of the parallel version of the algorithm for finding the ground state of an Ising spin glass is then analyzed and compared to the performance of a parallel version of the simulated annealing algorithm. It is found that methods for tuning the variable parameters in the tempering algorithm must be discovered before simulated tempering can be as effective as simulated annealing for general optimization problems. *Research Apprentice, 1992 SCCS REU Program; currently, Computer Science and Mathematics major, University of Texas of the Permian Basin, Odessa, TX. 1. INTRODUCTION Frequently in physics, systems involving many degrees of freedom arise. To derive properties of such systems using the normal approaches of classical mechanics would prove impossible because of the complexity of the system, so in many cases statistical mechanics is used to give a probabilistic description of the physical system. Monte Carlo techniques employing pseudo-random number generators can be used to simulate such systems which are probabilistic or statistical in nature. These simulations are useful for describing the thermal, electronic, and magnetic properties of matter, metallic alloys, diffusion, percolation, quantum chemistry, field theory, etc. [1]. Simulated tempering is a new Monte Carlo technique which has been proposed for effectively simulating physical systems which have different co-existing phases at some non-zero temperature. This procedure has been shown to greatly improve over other Monte Carlo methods for simulating the Random Field Ising Model [2]. Because of its recent development, however, the potential of this algorithm has not yet been fully investigated. It is believed that simulated tempering might overcome the difficulties encountered by conventional Monte Carlo techniques in simulating spin glasses and other disordered systems. In addition to the potential benefits offered by simulated tempering for simulating physical systems, the algorithm might also be applicable to general optimization problems. This belief is founded on simulated tempering's relationship to the well known Monte Carlo technique called simulated annealing. Simulated annealing has demonstrated its ability to find approximate solutions to general optimization problems [3]. Because of the similarity between the two algorithms, it is believed that simulated tempering should do at least as well as simulated annealing in general optimization problems. 2. FORMULATION OF SIMULATED TEMPERING ALGORITHM Models in which the many degrees of freedom are located on a lattice and interact locally are frequently used to simulate typical physical systems. In statistical thermodynamics the state of the system is characterized by some temperature (T). The spin glass is one such system whose simulation begins by assigning certain values to the microscopic degrees of freedom (atom positions, electron spins). In this case the possible spin values are called "spin up" and "spin down" which are assigned values of +1 and -1, respectively. All possible configurations of the system contribute to the thermal average at non-zero temperature and are weighted according to the Boltzmann factor exp( -B * E ), where B is proportional to 1/T and E is the internal energy for the configuration. For spin glass the equation for E is: E = - sum_ J_{i,j} s_i s_j , where means sum over nearest neighbor sites i and j, s_i and s_j are the spin values, and J_{i,j} denotes the strength of interaction between spins. For Monte Carlo simulation of the spin glass using the standard Metropolis algorithm [4], one starts with an initial configuration of spins and repeats the following steps over and over again: 1) Select one spin to flip (i.e. change from "up" to "down", or vice versa). 2) Compute the change in energy resulting from this flip. 3) Calculate the "transition probability" W = exp( - B * (change in E) ) for this flip. 4) Generate a random number 'x' beween zero and one. 5) If x < W flip the spin, otherwise, leave it unchanged. 6) Analyze the resulting configuration as desired. >From these simulations one is able to determine the average magnetization, the average energy, and the specific heat of the system, and numerous other quantities which are experimentally observable. Simulated annealing is a useful Monte Carlo technique which is used for simulations similar to this. In annealing, one "heats" up the system to some high temperature and slowly "cools" it, following some predetermined schedule, to a low temperature. This is done for the same reason annealing is used in manufacturing steel and other metals, which is because "quenching" the system (i.e. cooling it very rapidly) results in a metastable configuration of the system, whereas slow cooling results in the approximate minimization of the energy, and consequently, a more stable configuration. A straightforward generalization of simulated annealing for a general problem is not possible. This is because we are often interested in simulating systems at non-zero temperature, where it is not the energy of the system which must be minimized but the free energy. Simulated tempering is similar to simulated annealing but offers an effective way to minimize free energy. The difference between tempering and annealing is that where every change in temperature in annealing drives the system out of equilibrium, tempering consists of changing the temperature while remaining at equilibrium [2]. The complete simulated tempering algorithm for simulating a spin glass includes the following steps which are done after a predetermined number of Metropolis iterations: 1) Generate a random number 'y' between zero and one. 2) If y < 0.5 attempt to decrease the temperature, otherwise, attempt to increase the temperature. 3) Accept the temperature change with a "transition probability" proportional to exp( - (change in B) * E ). Some footnotes are necessary to further clarify the above three steps: i) The different temperatures are stored in an array in decreasing order. When the temperature to be updated is an extreme value, there is only one possible temperature change. In this case, the attempted change in temperature mentioned in step 2) is predetermined. ii) The amount by which the temperatures are changed, and hence, the values in array of temperatures, is determined before the actual execution of the code. This choice is made such that the acceptance of a temperature change is neither too high nor too low, and is dependent on the system being simulated and the temperature under investigation. For a discussion on effectively choosing the changes in temperature see ref. [2]. 3. USING SIMULATED TEMPERING FOR OPTIMIZATION Simulated annealing was originally used to find the ground state of a spin glass, which is a hard (NP-complete) optimization problem. It was later realized that this method could be applied to a general optimization problem by simply replacing the equation for the energy of the system by the cost function to be minimized, and introducing a fictitious temperature parameter [3]. This has proven to be a very effective optimization technique for many difficult optimization problems. One problem with simulated annealing is that the optimal cooling schedule is not known, and must be determined by trial and error. With simulated tempering, the changes in temperature are at least partly determined by the system, with a Metropolis-like update. Also, in simulated annealing there is only one cooling step, which finds a single local minimum that is hoped to be close to the global minimum. With simulated tempering the system is constantly being heated and cooled, so many local minima are obtained, and the smallest of these can be chosen as an approximation to the global minimum. It is therefore likely that simulated tempering will do at least as well as simulated annealing in general optimization problems. To check this possibility, simulated tempering was used to find the ground state of a spin glass. This requires setting one of the temperatures in the temperature array to be T = 0. There are still some free parameters in this problem, such as the number of possible temperatures and their values, and the number of Metropolis sweeps performed at each temperature. Methods for choosing these parameters also need to be investigated. 4. IMPLEMENTATION OF THE ALGORITHM ON A PARALLEL COMPUTER The main computational cost of the simulated tempering algorithm occurs during system updates, specifically, from the multitude of near-neighbor calculations which must be made for each site. As mentioned previously, these updates are accomplished using the standard Metropolis Monte Carlo algorithm. This process is straightforward to parallelize using domain decomposition. Because the algorithm is local and regular, perfect load balancing and good efficiency was achieved. The simulated tempering algorithm was implemented on a SIMD-architecture MasPar computer using the MasPar Programming Language (MPL), a data parallel version of the C programming language. The MasPar used consisted of a processor array composed of 4K processing elements. Domain decompostion was accomplished by mapping each lattice point, or subset of lattice points, to each processing element. During Metropolis updates of the system, attempts to flip spins on each processing element were done simultaneously. MPL's reduceAdd library function was then used to calculate configuration sums. While perfect load-balancing can be easily achieved, this is only possible if the lattice size is at least twice as large as the number of processing elements in the processor array. This is because the energy at a site, and thus the "transition probability" of flipping the spin at this site, is dependent on the "spin values" of this site's nearest-neighbors. Because of this, if each processing element contains only one lattice point spin updates must be done in a checkerboard fashion. This, in turn, requires that two passes be made to update the whole lattice, with half of the processing elements remaining idle during each pass. If there are at least two lattice points on each processing element, then it is possible to attempt to update one spin on each processing element simultaneously without encountering the above complication. However, if the lattice size is not a multiple of the number of processing elements, perfect load-balancing is again impossible for at least one pass during a system update. This problem necessarily arises because some processing elements must store more lattice points than others, requiring a different number of spin updates on certain processing elements. The boundary conditions of the lattice were handled by imposing a toroidal topology on the system, which is common practice for problems of this type [5]. This corresponds nicely with the topology of the MasPar's processor array, which is also toroidal. In addition, the simulated tempering algorithm is a local algorithm, requiring only nearest-neighbor communication. For these reasons, the parallel version of the tempering algorithm was efficient. For general optimization problems, such as the Travelling Salesman Problem (TSP), the parallel update algorthm may be more complicated. One cause of this complexity is the long range communication often required to obtain good approximations to the optimal solution of such problems. Another cause involves the many different possible topologies for general optimization problems. For example, the topology of the TSP is a ring, which must be taken into consideration during domain decomposition. Because of the similarities between simulated annealing and simulated tempering, the parallel update algorithm would need to be incorporated into parallel versions of both algorithms, and neither algorithm would be more advantageous to use should this be a concern. Fortunately, many highly efficient load-balancing techniques for general optimization problems have been developed and continue to be greatly researched. Some of these techniques are mentioned in ref [6]. 5. RESULTS AND DISCUSSION The data from representative runs of the parallel tempering and annnealing programs is reproduced graphically in figures 1 and 2. Both figures plot time (measured in Metropolis sweeps) versus energy. These results are for a lattice consisting of 256 x 256 spin sites. The decreasing energy in both graphs corresponds to the decreasing of the temperature in the systems. When the temperature becomes zero, the systems simulated by both codes "freeze" at approximately the same energy. However, the annealing code produces slightly better results. While simulated annealing relies on this constant gradual decrease in the temperature of the system to obtain a good approximation to the minimum energy, this is not the case with the simulated tempering algorithm. Simulated tempering's credibility relies on its ability to both increase and decrease the temperature of the system under investigation. With this in mind, the goal of tempering is not to converge to an energy indefinitely at T=0, as in figure 2; rather, the goal of tempering is to converge to an energy at T=0 and then begin to increase the temperature of the system. It is then hoped that this increase in temperature will lead to the decorrelation of the system from its previous configuration, and eventually lead to the convergence of the system to a new minimum energy as the temperature decreases to zero once again. Figure 2 does show some evidence of tempering's ability to increase the temperature of the system in the area above 40-80 iterations. While this does demonstrate to a very small extent what was desired, it is threoretically useless for the system under investigation. The failure of the tempering algorithm was a result of the inability to develope an appropriate tempering schedule. Initially, methods described in ref. [2] were used to theoretically calculate the values of the previously mentioned constants, which are used to modify the acceptance rate of temperature changes. After this failed to produce the desired results, a method of trial and error was employed. The tempering schedule used to produce the graph in figure 2 gave the best results out of a multitude of simulations utilizing different tempering schedules. On the other hand, the annealing schedule used to produce figure 1 was produced solely by trial and error. This proved to be a simple task when compared to the effort required to produce even the least succesful of tempering schedules. One possible explanation for the difficulty encountered while developing the tempering schedule, concerns the use of an increased lattice size made possible by the parallelization. Simulated tempering is based on the fact that there is an overlap in the energies at contiguous temperatures in the tempering schedule [2]. However, this overlap in energies for extremely large lattices may be negligable enough to hamper tempering's performance. 6. CONCLUSION The simulated tempering algorithm is especially well suited for implementation on the data parallel MasPar computer system. This is because data parallelization is easily accomplished and perfect load-balancing is easily achieved. Also, the parallel code is efficient because simulated tempering is a local algorithm, requiring only near-neighbor communication. For general optimization problems such as the Travelling Salesman Problem, the parallel update algorithm may be more complex. An efficient parallel program running on the a MasPar computer system which performs simulated tempering on Ising spin models was produced. However, for this program to run effectively, especially for near zero temperature in optimization problems, procedures for tuning the variable parameters in the tempering algorithm must be developed. These parameters are the temperature changes, and the constants which can be added to the energy at each temperature. Estimates for these quantities can be obtained if the energies are known. However, to get reasonable acceptance rates for the Monte Carlo update on the temperature, these variables will need to be fine tuned, probably by some dynamical procedure. Good methods of tuning these variables to obtain good acceptances will need to be found before simulated tempering can be as effective as simulated annealing for general optimization problems. 7. ACKNOWLEDGEMENTS Work on this project was funded in part by Grant CDA-9200577 from the National Science Foundation for a Research Experiences for Undergraduates (REU) Site Program at the Syracuse Center for Computational Science (SCCS), Syracuse University, Syracuse, NY. Additional support for the SCCS REU Program was provided by the Office of the Vice President for Research and Computing, Syracuse University. Hardware and software support was provided by the Northeast Parallel Architectures Center at Syracuse University. Additional computing resources were provided by the University of Texas of the Permian Basin and Louisiana State University. The first author would like to thank the second author, Prof. Edward Bogucz, and Dr. Marcin Paprzycki for their continuous motivation and help throughout the course of this project. Without the help and support of those mentioned, this project would not have been possible. REFERENCES [1] Koonin, Steven E., Computational Physics, Addison-Wesley, 1986. [2] E. Marinari and G. Parisi, "Simulated Tempering: A New Montecarlo Technique", submitted to Europhys. Lett., February 1992. [3] van Laarhoven, P. J. M. and E. H. L. Aarts, Simulated Annealing: Theory and Applications, D. Riedel Publishing, 1987. [4] N. Metropolis et al., J. Chem. Phys. 21, 1087 (1953). [5] Binder, K., Application of the Monte Carlo Method in Statistical Physics, Springer-Verlag, 1984. [6] Fox, Geoffrey C. et al., Solving Problems on Concurrent Processors, Prentice-Hall, 1988. FIGURE CAPTIONS 1 Energy as a function of time measured in terms of Metropolis iterations (as produced by the parallel version of the simulated annealing algorithm) 2 Energy as a function of time measured in terms of Metropolis iterations (as produced by the parallel version of the simulated tempering algorithm)