Heath bad method for Ising model

[Gulli]

I'm trying to simulate a 100x100 Ising model in matlab. So far I've tried the metropolis method (which seems to work), but I'm having problems with the heat bath method. I'm not sure how it's supposed to work (I can't get the spins to all take the same +1 or -1 value even at really low temperatures, while I can easily get them to do that with the metropolis method): am I supposed to change spins based only on the state of its neighbors in the original configuration and am I supposed to choose the spins at random or can I go through them systematically, if it's the latter, am I supposed to repeat the whole process (which would kinda make it just as much work for a computer as the metropolis method)?

Anyway, can somebody tell me how it's supposed to work, roughly?

 

[Timo]

Isn't the "metropolis method" (that I personally refer to as "Metropolis choice" for the fact that it merely means fixing a free parameter to its optimum value) just a special (optimized) case of a "heat bath method"? In that case, your "heat bath method" would just be the same as before; only the probabilities to accept a proposed spin flip would change. I would not go through the spins systematically unless you can formally proof that it gives the correct results. Adding correlations between your Monte Carlo steps is usually not a desirable property of a simulation scheme.

 

[Gulli]

Isn't the "metropolis method" (that I personally refer to as "Metropolis choice" for the fact that it merely means fixing a free parameter to its optimum value) just a special (optimized) case of a "heat bath method"? In that case, your "heat bath method" would just be the same as before; only the probabilities to accept a proposed spin flip would change. I would not go through the spins systematically unless you can formally proof that it gives the correct results. Adding correlations between your Monte Carlo steps is usually not a desirable property of a simulation scheme.

Thanks for your reply.

My book "Introduction to Statistical Physics" (Kerson Huang) says 'It can be shown that one step in the heat bath method is equivalent to an infinite number of Metropolis steps.' It and other sources say the heat bath method is supposed to be much faster and only works for systems with a small number of possible spin values (the Ising model only has two: -1 and +1), unlike the Metropolis method. Also, what would be the physical justification that allows for the Metropolis method to be carried out with different sets of probabilities (as with a dice, you'd expect only one set of probabilities to be correct)? So I don't think the two methods are the same except for their probability sets.

My book advices me to go through the spins systematically for the Metropolis method (I think that is because the probability of selecting spins 1, 2 and 3 in sequence is the same as that of selecting spins 23, 48 and 71, so why bother with random selection that eats up computing resources?), it doesn't say how I should select spins for the heath bath method, some of my other sources say I should select at random.

 

[Timo]

Hi Gulli. I am not familiar with the "heat bath method". Here's a page I found that mentions and briefly describes is: http://cours-physique.lps.ens.fr/index.php/TD08_Heat-bath_algorithm_for_the_Ising_model_2011. I will assume that is what you books also talks about in the following.

According to this, you pick a random (random! see below) spin, and chose its new value according to a probability that depends on its neighbouring spins. Apart from some possible free parameters that is the same as flipping it with a probability that depends on its neighbouring spins and its current value.

My book "Introduction to Statistical Physics" (Kerson Huang) says 'It can be shown that one step in the heat bath method is equivalent to an infinite number of Metropolis steps.' It and other sources say the heat bath method is supposed to be much faster [...]

That is, at best, true for an infinite repeated Metropolis flipping of the same spin. And that is not the type of problems you want to solve with Monte Carlo. If you are interested in the problem in which only a single spin can flip, you can just solve that problem analytically. I don't see how the "heat bath" is supposed to be faster than the Metropolis move. But then, I don't know how your sources define "faster".

What would be the physical justification that allows for the Metropolis method to be carried out with different sets of probabilities (as with a dice, you'd expect only one set of probabilities to be correct)?

Generally, Monte-Carlo sampling is not physical itself (*), so neither is the reason. It is, however, quite easy to see: Assume that in the Metropolis algorithm you add a step where after having chosen a random lattice site, with 50% probability you either reject the flip immediately (and pick the next spin at random) or continue with the algorithm as normal. All that changes is that your algorithm slows down, because half of your proposed moves are rejected "for nothing". This is, of course, equivalent with accepting spin flips with only half of the original probability (note that in this case "always flip if the energy is reduced" is to be interpreted as a 100% flip probability).

So I don't think the two methods are the same except for their probability sets.

I still do, but admittedly didn't actually check the numbers.

My book advices me to go through the spins systematically for the Metropolis method (I think that is because the probability of selecting spins 1, 2 and 3 in sequence is the same as that of selecting spins 23, 48 and 71, so why bother with random selection that eats up computing resources?), it doesn't say how I should select spins for the heath bath method, some of my other sources say I should select at random.

Selecting the spins at random is correct for both algorithms (*). Selecting them sequentially is probably correct, too. But at least the mainstream proof for correctness I know doesn't work in this case. I am not convinced that it actually is faster, since you gain some CPU time for a single step, but introduce correlations between your steps that (in addition to putting the correctness of the whole algorithm in question) may well reduce the effectiveness of the algorithm globally. I would advice against sequential picking, but I understand that trading formal justification of your algorithm for a potential performance increase is tempting.

From what you wrote I have some doubts about the book you are using. The Amazon reviews and the fact that the author is over 80 years of age also support my feeling that the book is not a proper read for an introduction to Monte Carlo computer simulations. If you want to seriously go into Monte Carlo methods then a specialized book (i.e. one assuming you already know thermodynamics and statistical mechanics and only want to read about Monte Carlo methods as they are actually used in research in these fields) I can recommend Newman, Barkema: "Monte Carlo Methods in Statistical Physics" (disclaimer: as I am a 2nd generation "academic sibling" of one of the authors, this opinion is strongly biased by the "the book describes MC exactly the way I use it"-effect).

(*) Explaining these points in detail would be beyond the current scope. But I can gladly go into more detail if your are interested (for which it would help to know a bit about your background, though).

 

[Gulli]

Hi Gulli. I am not familiar with the "heat bath method". Here's a page I found that mentions and briefly describes is: http://cours-physique.lps.ens.fr/index.php/TD08_Heat-bath_algorithm_for_the_Ising_model_2011. I will assume that is what you books also talks about in the following.

According to this, you pick a random (random! see below) spin, and chose its new value according to a probability that depends on its neighbouring spins. Apart from some possible free parameters that is the same as flipping it with a probability that depends on its neighbouring spins and its current value.

Yes, that's it.

Generally, Monte-Carlo sampling is not physical itself (*), so neither is the reason. It is, however, quite easy to see: Assume that in the Metropolis algorithm you add a step where after having chosen a random lattice site, with 50% probability you either reject the flip immediately (and pick the next spin at random) or continue with the algorithm as normal. All that changes is that your algorithm slows down, because half of your proposed moves are rejected "for nothing". This is, of course, equivalent with accepting spin flips with only half of the original probability (note that in this case "always flip if the energy is reduced" is to be interpreted as a 100% flip probability).

I still do, but admittedly didn't actually check the numbers.

I see, if the numbers work out then you could be right.

Selecting the spins at random is correct for both algorithms (*). Selecting them sequentially is probably correct, too. But at least the mainstream proof for correctness I know doesn't work in this case. I am not convinced that it actually is faster, since you gain some CPU time for a single step, but introduce correlations between your steps that (in addition to putting the correctness of the whole algorithm in question) may well reduce the effectiveness of the algorithm globally. I would advice against sequential picking, but I understand that trading formal justification of your algorithm for a potential performance increase is tempting.

My assignment is to play around with the Ising model using Monte Carlo and simulate ferromagnetism. Testing whether the extra iterations required by correlation between steps negate time saved by not having to choose random spin sites, is exactly the sort of thing that fits the playing around part of my assignment (the formal proof isn't terribly important), thanks for the tip!

From what you wrote I have some doubts about the book you are using. The Amazon reviews and the fact that the author is over 80 years of age also support my feeling that the book is not a proper read for an introduction to Monte Carlo computer simulations. If you want to seriously go into Monte Carlo methods then a specialized book (i.e. one assuming you already know thermodynamics and statistical mechanics and only want to read about Monte Carlo methods as they are actually used in research in these fields) I can recommend Newman, Barkema: "Monte Carlo Methods in Statistical Physics" (disclaimer: as I am a 2nd generation "academic sibling" of one of the authors, this opinion is strongly biased by the "the book describes MC exactly the way I use it"-effect).

(*) Explaining these points in detail would be beyond the current scope. But I can gladly go into more detail if your are interested (for which it would help to know a bit about your background, though).

You might be right about my book: it doesn't explain much and I actually got stuck for an hour or so when programming the Metropolis method because the author had switched < and > signs somewhere.

I just started graduate school (I have a bsc in physics and know thermodynamics, etc...), taking introductory Monte Carlo courses, but no one ever explained the heat bath method to me and I guess my book really threw me of course because the language it uses to (very briefly) describe the heat bath method is very different from online sources. Anyway, I'll try implementing it now and report back what I find.

 

[Gulli]

I think I got the heat bath method working now, your replies has been very helpful and I found another mistake in the book (the author switched mixed things up again) that kept me from implementing the heat bath method succesfully. I'll never use that book again: it's the second edition but it still has two deadly mistakes on the four pages that deal with the Ising model.

 

[lww]

How about sweeping the lattice for more times?

 

[Gulli]

Well, I've now successfully implemented bith methods (Metropolis seems to be faster, I think because its steps don't always require computing an exponential).

Rather than randomly choosing spins both methods systematically go though all spins in my implementation because drawing random numbers doesn't work with the Metropolis method (because you can't eliminate double counting unless you keep a record of all the spins you've visited previously, which is tedious and time consuming) and for the heat bath method it's just very slow (computer takes too much time drawing random numbers).

Both methods are dramatically faster when I use an algorithm (based on the number of spins that are neighboring at least one spin with opposite spin, which is a measure for how big islands of +1 or -1 have become) that occassionally goes through the lattice systematically but only looks at spins that have at least one opposite neighbor (these spins tend to be located on the edge between two islands). When the islands become big enough "occassionally" changes into "always" and the program will only look at the edges between islands because that's where all the action would happen anyway (flipping a spin in the middle of an island would just result in noise since it's very likely to be forced to flip back by its neighbors the next round, so looking at these spins is a waste of time).

This is an algorithm by my own design and I think it doesn't change the outcome in any measurable way (since the odds against, a spin in the middle of a big island suddenly flipping and converting the entire island and then changing the outcome for the entire model, are astronomical), but of course there's always a chance I'm wrong.

 

[torquil]

Rather than randomly choosing spins both methods systematically go though all spins in my implementation because drawing random numbers doesn't work with the Metropolis method (because you can't eliminate double counting unless you keep a record of all the spins you've visited previously, which is tedious and time consuming) and for the heat bath method it's just very slow (computer takes too much time drawing random numbers).

It is not important if you are flipping at a given site more than one time during "one pass". With random site selection on e.g. an N\times N grid, it is customary to just define "one pass" as N^2 spin flips at random sites.

Both methods are dramatically faster when I use an algorithm (based on the number of spins that are neighboring at least one spin with opposite spin, which is a measure for how big islands of +1 or -1 have become) that occassionally goes through the lattice systematically but only looks at spins that have at least one opposite neighbor (these spins tend to be located on the edge between two islands). When the islands become big enough "occassionally" changes into "always" and the program will only look at the edges between islands because that's where all the action would happen anyway (flipping a spin in the middle of an island would just result in noise since it's very likely to be forced to flip back by its neighbors the next round, so looking at these spins is a waste of time).

This is an algorithm by my own design and I think it doesn't change the outcome in any measurable way (since the odds against, a spin in the middle of a big island suddenly flipping and converting the entire island and then changing the outcome for the entire model, are astronomical), but of course there's always a chance I'm wrong.

Have a look at "block" or "cluster" algorithms:

http://www.netlib.org/utk/lsi/pcwLSI/text/node292.html

 

[Gulli]

It is not important if you are flipping at a given site more than one time during "one pass". With random site selection on e.g. an N\times N grid, it is customary to just define "one pass" as N^2 spin flips at random sites.

That's how I defined a lattice sweep when trying out random selection. I allow for up to 1000 sweeps before I quit the algorith. So you're saying it really doesn't matter whether I sweep systematically or select randomly?

Have a look at "block" or "cluster" algorithms:

http://www.netlib.org/utk/lsi/pcwLSI/text/node292.html

I suspected I wasn't the first to think of such things, though what that link describes is still quite different from my algorithm. Sadly my algorithm doesn't seem to display a phase transition, but then again neither does the heat bath method which I copied from books and websites (only the pure Metropolis Method seems to have a phase transition in my implementation).

 

[torquil]

That's how I defined a lattice sweep when trying out random selection. I allow for up to 1000 sweeps before I quit the algorith. So you're saying it really doesn't matter whether I sweep systematically or select randomly?

It doesn't matter much, and both methods should give the correct answer in the end. However, my comment was mostly to alleviate your concern about flipping at the same site multiple times within one sweep.

See pages 48-51 here.

I suspected I wasn't the first to think of such things, though what that link describes is still quite different from my algorithm. Sadly my algorithm doesn't seem to display a phase transition, but then again neither does the heat bath method which I copied from books and websites (only the pure Metropolis Method seems to have a phase transition in my implementation).

Are you sure your heat bath algorithm is correct? You should get the same result in both cases, as long as your numerical method converges to a good enough accuracy to see the phse transition. Perhaps your temperature is defined differently in the two cases, so that you are searching for a phase transition at the wrong temperatures when you are using the heat bath method?

 

[Timo]

A few random comments:

1) If you are calculating exponentials in a step of either method then your implementation is inefficient. There are only five possibilities for the total magnetization of the neighbors. Calculate your flip probabilities in advance and store them in a lookup table.

2) The Metropolis algorithm should still be more efficient, since it is, by definition, the most dynamic single spin flip algorithm with randomly chosen spins. Again, that statement may be impacted by the fact that I don't know the Heat Bath method, but from your description it shouldn't beat what you call the Metropolis Method.

3) Ignoring spins inside a large cluster may be faster. But it is wrong. And unless you can control how wrong you go (and verify that the error you introduce does not matter for your purpose) you shouldn't intentionally program something like that (unless you are an experimentalist and merely want "simulation results" to "validate" your experiments - then you can program anything because your referee won't know proper MC simulations, either ). If you want the T=0 state you do not need to program anything - just set all spins up or down and call it a result.

4) Consider it theorist nit-picking, but I cannot resist: You won't see a phase transition in a finite system in any case.

 

[Gulli]

It doesn't matter much, and both methods should give the correct answer in the end. However, my comment was mostly to alleviate your concern about flipping at the same site multiple times within one sweep.

Thanks.

Are you sure your heat bath algorithm is correct? You should get the same result in both cases, as long as your numerical method converges to a good enough accuracy to see the phse transition. Perhaps your temperature is defined differently in the two cases, so that you are searching for a phase transition at the wrong temperatures when you are using the heat bath method?

No, I didn't define the temperature any different, the heat bath method does get worse at high temperatures, but those temperatures are well above the critical temperature of the metropolis method, but maybe it has something to do with my lattice only being 100x100.

A few random comments:1) If you are calculating exponentials in a step of either method then your implementation is inefficient. There are only five possibilities for the total magnetization of the neighbors. Calculate your flip probabilities in advance and store them in a lookup table.

I know, at first I thought I couldn't do that because I wanted the simulation to work with more than just one value for the temperature and the magnetic field, but of course I could still calculate the probabilities once per sweep. Oh well, my deadline passed and sleep deprivation kept me from having that idea before the deadline so I handed my code in with this flaw.

2) The Metropolis algorithm should still be more efficient, since it is, by definition, the most dynamic single spin flip algorithm with randomly chosen spins. Again, that statement may be impacted by the fact that I don't know the Heat Bath method, but from your description it shouldn't beat what you call the Metropolis Method.

It was just my book saying that "one step in the heat bath method was equivalent to an infinite number of steps in the metropolis method", in my implementation the metropolis method was much faster.

3) Ignoring spins inside a large cluster may be faster. But it is wrong. And unless you can control how wrong you go (and verify that the error you introduce does not matter for your purpose) you shouldn't intentionally program something like that (unless you are an experimentalist and merely want "simulation results" to "validate" your experiments - then you can program anything because your referee won't know proper MC simulations, either ). If you want the T=0 state you do not need to program anything - just set all spins up or down and call it a result.

I'm not totally ignoring them: the algorithm ignores them with a probability based on the number of spins that are not on the egde of a cluster at the beginning of the sweep, this probability is very low at the beginning and stays below 75% almost until the end. I figured that ignoring them doesn't introduce a large error, but I did see that the metropolis lost its phase transition and did include in my report that this indicates an error. It wasn't a thesis, btw, it was the final graded assignment of one subject and the whole idea was to screw around with Monte Carlo techniques and try out your own ideas and algorithms.

4) Consider it theorist nit-picking, but I cannot resist: You won't see a phase transition in a finite system in any case.

What do you mean by that? Are you saying the phase transition is more blurry in finite systems or that I shouldn't see one at all.

 

[torquil]

Thanks.
maybe it has something to do with my lattice only being 100x100.

No, 100x100 should be plenty to see obvious signs of a phase transition.

What do you mean by that? Are you saying the phase transition is more blurry in finite systems or that I shouldn't see one at all.

Timo means that strictly speaking, a phase transition involves discontinuities and singularities of the thermodynamic quantities as functions of temperature, and this can only happen for infinite systems.

However, for finite systems you will instead get a "blurry" transition as you call them. E.g. instead of a singular peak of the heat capacity at the critical temperature, you would get only a nonsingular peak for a finite system, but still a clearly visible peak for a 100x100 system.

 

[Gulli]

If you are calculating exponentials in a step of either method then your implementation is inefficient. There are only five possibilities for the total magnetization of the neighbors. Calculate your flip probabilities in advance and store them in a lookup table.

I did a quick test in Matlab to be sure, but as I suspected the exponentials don't take much time to calculate, drawing random numbers takes up >10x more time, so this inefficiency with exponentials didn't matter too much.

P.S. I passed the course.

 

[Timo]

That's an interesting finding. Perhaps even more interestingly, I found exactly the opposite: a single exp() taking ten times as long as drawing a random number in a simple test (although I am sceptical of my results, too - a factor of ten seems a bit too much).

Here's the code I used:

#include <iostream>
#include <cmath>
#include "MTRand.h"

using namespace std;


const int N=100000000;
mccore::MTRand RND(31);

double randoms() {
  double result(0.); 
  for (int i=0; i<N; ++i)
    result += RND();
  return result;
}

double expos() {
  double result(0.); 
  double d(-1.e-5);

  for (int i=0; i<N; ++i, d*=1.1) {
    result += exp(d);
    if (d < -10.) d=-1.e-5; // keep d in sensible range, CPU time may explode, otherwise.
  }
  return result;
}


int main() {

  cout<<randoms()<<endl;
  // cout<<expos()<<endl;
}

Compiled with g++ -O3, random number generator is a Mersenne Twister generator. Code ran 3 seconds for random numbers and 25 seconds for exponentials. Granted, a lot of the performance of the exponentials depends on the number range you put in, but I don't see how you'd come from 25 seconds to 0.3 seconds. Perhaps you were putting integers into the argument of exp() and Matlab internally used a lookup table?

Anyways, nice to hear you passed your course. Most of my comments here were actually targetted at a more advanced audience, say a PhD student starting to work in computer simulations. So don't worry if you found them too picky for your case.