# Iteratively refine bounds on exp for Metropolis criterion

In Monte Carlo simulations, using the Metropolis criterion, one often has to compare a random number $a$, $0 \leq a < 1$, to the Boltzmann distribution $exp(-\beta\Delta E)$, where $\Delta E$ is the energy cost of a transition between two states in the simulation.

A common optimization, is to check if $\Delta E > 0$, in which case, $a$ is certainly smaller than $exp(-\beta\Delta E)$.

For my simulation, it seems, the calculation of the exponential function, seems to take the most time and I'm looking into methods of reducing the cost.

One method I thought of, is based on the fact that for even $K$, the Taylor series $\sum_{k = 0}^{K}x^k/k!$, is a lower bound to $exp(x)$, while for odd $K$, it is an upper bound. Thus, one can write a function like the following,


static inline int fast_metropolis(double a, double x){
double est = 1.0 + x;
double f = x;
int i = 1;
while(true){
if(rand_num < est) return 1;
f *= x / (i + 1);
est += f;
if(rand_num > est) return 0;
f *= x / (i + 2);
est += f;
i += 2;
}
return 0;
}



Unfortunately, in hindsight, the above code doesn't work very negative arguments, since the absolute value of the bounds increases.

Are there any alternatives to the method described above i.e. could we iteratively refine some upper and lower bounds while making sure their absolute value is decreasing?

Or am I better of using some 'fast' estimate of the exponential function?

• Typically, the "fast" exponential functions trade off speed for accuracy; if you're willing to accept a less accurate result, you can compute it faster, as you've intuited above. Can you decrease the precision of your random numbers without compromising the physics of your simulation? If so, then you might also be able to limit the accuracy of your exponential to the precision of the random numbers you generate. – Geoff Oxberry Jun 30 '14 at 7:52

The simplest solution to my question was quite obvious, too bad I didn't see it before posting the question.

You can split the argument of $exp(x)$ into an integer part and a fractional part. For the integer part one can precompute a lookup table, while for the fractional part you can use the algorithm described in the question.

Below I show the final code I used,


static inline int fast_metropolis(double rand_num, double x){
int exp_i = -(int)x;
if(exp_i > 256) return rand_num < exp(x);
double prefactor = (exp_i > 0)? exp_table[exp_i - 1]: 1.0;
double k = x + exp_i;
double f = k;
double est = 1.0 + k;
int i = 1;
while(1){
if(rand_num < prefactor * est) return 1;
f *= k / (i + 1);
est += f;
if(rand_num > prefactor * est) return 0;
f *= k / (i + 2);
est += f;
i += 2;
}
return 0;
}



The algorithm described, was quite accurate in my simulations, with only a $8.4 \cdot 10^{-6}\%$ of the cases being different from using the standard library's exp function. In comparison, the fast exponential function, found here, was different in $0.3\%$ of the cases, although it had slightly better performance in simulation. Performance tests with a uniform distribution for the argument of $exp(x)$ seem to indicate a 3x performance gain.

I will not mark my answer correct for a couple of days, in case someone comes up with a faster/more accurate algorithm.