Hello, bonjour, gruss gott, laba diena, Welcome to the the 9th week of Statistical Mechanics: Algorithms and Computations from the Physics Department of Ecole Normale Supérieure. During our course, we have mostly concentrated on equilibrium statistical mechanics Our algorithms allowed us to compute partition functions, energies and condensed fractions but they were rarely concerned with time-dependent phenomena. Physical time played a minor role except in the molecular dynamics calculations for the hard disks. Monte Carlo time really was of secondary importance for example the Wolff cluster algorithm that we discussed last week had a dynamics that had nothing in common with the physical dynamics of the Ising model. In this lecture, we explore time-dependent phenomena and strangely enough we are right back with the local Monte Carlo algorithm, as it provides often an excellent model for a time evolution in a physical system. So, what we are looking for is an efficient implementation of markov_ising.py in dynamical Monte Carlo, we are stuck with a given dynamical rule and we want to implement this rule as efficiently as possible. In this lecture we will treat markov_ising.py and study the faster-than-the-clock implementation and also its main limitation, that's called the futility. In this week's tutorial we plunge into the great and very general method called simulated annealing. Simulated annealing is a tool for solving difficult optimization problems with or without a relation with physics. It consists in reformulating the optimization problem in terms of an artificial physical system whose ground-state at zero temperature or at infinite pressure contains the solution to the original task. This state is then approached through a slow decrease in temperature or a slow increase in pressure. We will apply simulated annealing to the famous 13-spheres problem, that already Isaac Newton worked on. Simulated annealing gives a great solution. In this week's homework session - the final homework session of this course - we will also be concerned with simulated annealing in the 13-spheres problem and in the equally famous traveling salesman problem. So, let's get started with week 9 of Statistical Mechanics: Algorithms and Computations. To begin with, let us consider a magnetic field as shown here and a single Ising spin whose dynamics we want to study. We set the temperature equal to 1 and the energy of this spin is equal to -h if it is up and is equal to +h if it is down. In short, the energy of sigma is -h*sigma, where sigma can be plus or minus one. To simulate the dynamics, we use the Metropolis algorithm min(1, pi_new / p_old) So the probability to flip the spin from minus to plus is equal to one, because the statistical weight of the new configuration is larger than the statistical weight of the old configuration. And the probability to flip from +1 to -1 is equal to min(1, pi_new / pi_old) = exp(-beta h) / exp(beta h) = exp(- 2 beta h) In Python, this gives the following program naive_spin_dynamics.py, where we have set beta=1 Now, notice that if at time t the spin is down it will certainly be up at time t+1 This is the essence of Metropolis dynamics and we don't discuss here whether this is realistic. In any case, this algorithm satisfies detailed balance and in the long run - at long times - the two configurations up and down will be visited with their corresponding statistical weights: pi_up and pi_down. So now, let's look at the movie version of naive_spin_dynamics.py. Here is output of this program at high temperature: there are a lot of flips. And here is output at lower temperature you see: a long time nothing is happening and the spin is in the up position, then there is a flip, immediately followed by a back-flip, and then followed by a long desert of the spin being up and here is the dynamics of our system at very low temperature. Nothing going on, and just an occasional flip and our computer heating up without any result. Now remember: the probability to flip from down to up is equal to one and the probability from up to down is equal to exp(-2h) so now, to make progress, we choose a particular value of h. In fact, we choose log(6)/2 so now, the probability to flip from up to down is exp(-2 log(6) / 2) =exp(-log(6)) =1/6 So look at this flip-die it has 5 neutral faces and one face that says "flip" so, simulating the Metropolis algorithm at magnetic field log(6)/2 is like rolling this die so, let's roll it.. First trial: no result, the spin does not flip Second trial Third time we are getting quite tired.. no result.. Fifth attempt Sixth.. no.. Seventh.. so now we have to flip the spin In Python this gives the following program: naive_throw.py and here is output of this program so now, let us program this problem a little more intelligently than we have done before. So, the probability of throwing a blank face of our beautiful flip-die is equal to 5/6 The probability of throwing twice a blank face is equal to (5/6)^2 and the probability of throwing k times a blank face followed by once the flip face is given by (5/6)^k times 1/6 or (1 - 5/6) So this is equal to (5/6)^k - (5/6)^(k+1) so this can be made into a tower sampling problem between 0 and 1 with horizontal lines at 5/6 (5/6)^2, (5/6)^3 ... we throw a pebble into this tower and if it falls for example in the upper box we have a new flip right away if it falls into the second box, we have to wait once if it falls into the third box, we have to wait twice to find out the number of times we have to wait we have to compute the box, which means (5/6)^(k+1) smaller than the random number smaller than (5/6)^k Taking logarithms, we find that k is the integer part of the logarithm of the random number divided by the logarithm of 5/6. And this means - in our naive_throw program - that the time after which we have a new flip is given by t_flip=.. So, in Python this gives the program fast_throw.py its output is indistinguishable from the naive program but it advances one iteration per accepted flip and not per time, it is for this reason that it is called a faster-than-the-clock algorithm. Of course, a very simple one. We may now write the fast_spin_dynamics program, as shown here it will remain exactly one step in the down-spin configuration with spin equal to -1 and an average of 6 steps in the up-spin configuration so that the magnetization on average, will be 5/7 We can obtain the same result also by a pedestrian approach So now, before continuing please take a moment to download, run and modify the programs we just discussed So we had the die-throwing program naive_throw.py and the fast_throw.py, that implemented the faster-than-the-clock idea The output of the two programs is completely indistinguishable and then we had the related programs naive_spin_dynamics.py and fast_spin_dynamics.py that implemented the dynamics of a single spin in the magnetic field. From a single spin in a field, we now pass on to a full-fledged simulation of the Ising model on N sites in the faster-than-the-clock flavor We consider a spin sigma, all the spins that depends on time and we denote by sigma_k the configuration that comes out of sigma by flipping the spin number k The energy change that is associated with this spin flip is delta_E = E[sigma_k] - E[sigma] In the Metropolis algorithm the probability to go from sigma to sigma_k is given by (1/N) min(1, exp(-beta delta_E)) the 1/N stands for the probability of sampling the spin number k In the Ising model and its variances (for example spin glasses) at low temperature, most of the spin flips are rejected. So, let's roll back our sleeves and implement the faster-than-the-clock approach for the Ising model. Let's avoid all rejections. But there is a problem. In the die-rolling algorithm fast_throw.py we notice that the quantity that was important was not 1/6 the probability to have a flip but it was rather 5/6, the probability of doing nothing Now, in our N-sites Ising model we also have to compute the probability of doing nothing but as we have N choices for doing something, computing this probability takes a bit of time. So we can compute now lambda, the probability of doing nothing that was 5/6 before, is given by 1 minus the sum over all the probabilities we just wrote down a few seconds ago and we can implement the faster-than-the-clock approach through the tower sampling from lambda, lambda^2, lambda^3, ... In Python, this is implemented in this program dynamic_ising.py and you see that there are two tower samplings the first tower sampling is in lambda, it is from 0 to 1 and it computes after which time we accept a new spin flip then we have to decide which spin to flip and this again is a discrete probability and it is done in a second tower sampling from lambda to 1 and it is there that we compute k, the value of the spin that is going to flip. Output of dynamic_ising.py is completely indistinguishable from markov_ising.py and the only problem that we have is that doing one iteration of this algorithm takes of the order of N operations besides that, there is no rejection in dynamic_ising.py The complete recalculation of lambda - the probability of doing nothing - can be avoided because flipping a spin somewhere in the system does not change the local environment elsewhere. So the number of environments of a spin with respect to its neighbors fall into a finite number of classes: a number of up spins and down spins that are neighbors of a given spin. And flipping a spin somewhere in the system just changes the numbers of members of each class. This was first implemented by Bortz, Kalos, and Lebowitz in 1975 They called their algorithm the n-fold way, where n stands for the number of classes. Let's take as an example the Ising model in two dimensions with periodic boundary conditions As you see here, the number of classes is 10 from a up spin with four neighboring up spins to a down spin with four neighboring down spin The tower of N probabilities can thus be reduced to a tower of 10 classes of spin flips if only you know the number of spins per class So, in this example the central spin that we want to flip has just one plus neighbor, so it lives in class 3 flipping it brings it into class 8 and makes change also the classes of the neighboring spins this very elementary book-keeping is written up in dynamic_ising_patch.py With respect to the original version (dynamic_ising.py) where each spin flip involved of the order of N operations now we can do one spin flip even for a large system in of the order of one operation At small temperature, this algorithm dynamic_ising_patch.py is much faster than the original markov_ising.py but remember that it implements exactly the same dynamics moving from accepted spin flip to the next accepted spin flip. So, before continuing, please take a moment to download, run and modify the two programs we discussed in this section For once there was dynamic_ising.py that worked with a tower of N probabilities and the nice n-fold way algorithm first implemented by Kalos Bortz and Lebowitz in 1975 that uses the n-fold way and goes from the tower of N probabilities to a tower of 10 probabilities. So the complexity of this algorithm is O(1) per accepted spin flip. So in conclusion, we have given in this lecture an introduction to faster-than-the-clock algorithms, that implement efficiently the dynamics of our system. This algorithm spend for per accepted move and no longer per attempted move. But, going from the word "attempted" to the word "accepted" move involves book-keeping. This book-keeping makes also that in many applications they are not as good as they seem to be because to go from one configuration a to a configuration b you may have to do some book-keeping as we showed in the Ising model and then on the next move you will simply go back from b to a. This is called the futility problem. To overcome this futility problem many people have implemented quite elaborate algorithms that even increase the expense per accepted move In this week's tutorial we will take up the study of dynamics algorithms from another point of view, the one of simulated annealing. Simulated annealing is a great optimization algorithm. It also shows that the sampling approach is not only related to integration as we have said for many weeks now but it can also be used for searching and for orientation in complicated spaces. So have fun with this week's tutorial and homework sessions and see you again next week in the final week of Statistical Mechanics: Algorithms and Computations for a discussion of the alpha and the omega of Monte Carlo followed by a resume and by a big party..
