# Metropolis, simulated annealing and proposal distributions

I'm trying to understand the physical meaning of the proposal step generating function in Metropolis algorithm.

In the original paper, and most derivations I found, it seems that it's not much important as long as it's symmetric (for detailed-balance) and allows a good exploration of the configuration space (e.g. $\sim 0.5$ acceptance ratio). Usually it's either a uniform or a normal distribution.

Now when studying Simulated Annealing, most papers assume a common use of a normal distribution with $\sigma = T$.

Given simulated annealing is deeply rooted in physical analogy, if we put the temperature there, there must be some physical sense we're trying to model. I mean, we could use a fixed uniform distribution with the proper width and still get good enough parameter space exploration and boltzmann distribution in the result, so why we use a gaussian and why the temperature dependence?

Also, when Cauchy distribution is introduced (e.g. by Szu and Hartley, 1987) they say the temperature used as a control parameter is kind of artificial, implying that the one in the Gaussian was not.

So what's the physical phenomenon with Gaussian, temperature dependent, steps we're modeling?

I know it all can be better explained with markov chains, priors, posteriors and likelihood, but that way you lose the physical analogy...

• I know nothing about physics, but the temperature $T$ in the target $\pi_T$ and the proposal are two different things: the decreasing temperature induces the target $\pi_T$ to become more and more peaked, hence the maxima to become more and more attractive and sticky. The proposal is used to walk around the current target $\pi_T$ with a reasonable chance of success (e.g., an acceptance probability of $0.234$). Feb 13, 2016 at 8:22
• @Xi'an that's the role of the temperature in the target distribution, which is usually a Boltzmann one, the lower the temperature, the higher the rejection rate for distant steps. What I cannot understand though, is the role of the temperature in the proposal distribution which is usually something like $g(x) \simeq \exp(-x^2/T^2)$. Guess it comes from some thermal diffusion model, but I'd like to see a more formal derivation... Feb 13, 2016 at 8:34
• Again, I cannot give you any physical reason. From a probabilistic perspective, the scale of the random walk, whether normal or not, must adapt to the scale of the target $\pi_T$, according to theoretical developments (Roberts et al., 1997) that show the optimal variance of the random walk is twice the variance of the target. Hence the presence of $T$ in the Gaussian density. (I am unsure it is $1/T^2$ though.) The fact that it does not seem to matter for the Cauchy random walk is that the Cauchy distribution has no variance. Feb 13, 2016 at 8:40
• Will look into the paper, thanks. $1/T^2$ is taken from e.g. Szu and Hartely, 1989 and Ingber, 1989 Feb 13, 2016 at 8:49
• It all depends on the variance of the target distribution, $\pi_T$, so $T^2$ or $2T^2$ is not always the optimal choice for the random walk. Feb 13, 2016 at 9:07

I suppose one could draw a physical analogy, noting that the proposal length scale should be proportional to the kinetic energy of the particle, but as @Xi'an I too would argue that this is simply a practical choice to maintain sensible acceptance rates as the temperature goes down.

To see this, assume we want to optimize a typical statistical problem, where the target function is a likelihood, which is asymptotically normal. Assume for simplicity that the optimum is at 0, and the width (sd) is $\sigma_t$. We optimize the log Likelihood, which is proportional to

$$LL \propto - x^2 / \sigma_t^2$$

x being the distance from the maximum at 0. In the standard simulated annealing, our acceptance probability would thus be proportional to

$$p_a \propto e^{- \frac{x^2} {\sigma_t^2 \cdot T } }$$

Now, consider that the expected value of the random variable x depends on the proposal function

$$<x^2> \sim proposal$$

As this is tagged self-study, I'll leave it to you to see how you should scale the proposal function to get constant acceptance rates (I would argue that this is a sensible goal).

• Thanks, I'll have to think about it, but it seems a promising and simple enough approach.About constant acceptance rates, I agree it should be a sensible goal but... Tipical rates start with a plateau at high temperatures, go roughly linearly as the annealing goes on, and reach a constant plateau $\sim 0$ at low temperatures. This is often used as a stopping criterion (IIRC Kirkpatrick used this very same criterion in his paper), so it doesn't seem a goal much shared in the Annealing literature. Feb 27, 2016 at 15:15
• Yes, OK, but the first part of the practical experience may also be because you haven't reached the optimum yet - in the code above I assume you are proposing to go away from the optimum. It's also my experience that acceptance goes to zero, but this may also be because people tune the proposal to shrink slower than the acceptance (which is what your $\sigma$ = T suggests to me). I wonder about the sense of that though, because there would be other ways to check convergence, so why would you waste your computation time on all these rejections. Feb 27, 2016 at 17:38
• In my practical experience with $\sigma = T$ the acceptance ratio grows again at small temperatures (small steps, small $\Delta E$, small $T$, proposal $\rightarrow 1$). Using $\sigma \propto T$ with a proper, big enough, constant, this effect is mitigated, but I didn't see any real gain (e.g. faster or better convergence) over constant $\sigma$ steps. If I may go a little off-topic, could you give me some pointers about this better methods of checking convergence? right now I'm monitoring energy variance but it's not that different in behaviour from acceptance rate. Feb 28, 2016 at 8:04
• better methods - not a big secret - just check if you still get better with time (by some epsilon), if not stop. Feb 28, 2016 at 19:01
• I finally had a little time to look into this again. $\sigma = \sqrt{T}$ actually leads to roughly constant acceptance rates as expected. This is not what Szu, Geman and Geman, and a lot of other literature on Bolzmann and Cauchy machines, use though: they use $\sigma = T$ with steps $\propto \exp{-x^2 / T^2}$ or $\propto T/(x^2+T^2)$. Same for Scipy annealing implementation (not that it really proves anything...) Mar 3, 2016 at 9:29

Here's an analogy from physics - consider the distribution in $n$ dimensions to have many peaks or valleys throughout. In a MCMC sampling problem, the 'height' of the peaks would be how improbable the state was - we would want to go to the valleys. Now, let us imagine this same surface to be a the map of potential energy for a physical system. Then, since a particle or object would want to find a minimum for the potential, it will try to go 'downhill'.

Now, in a classical sense, the particle can never pass through or 'tunnel' through the potential barriers. These barriers are the peaks in our distribution. However, in quantum theory, these particles always have a nonzero probability of tunneling through these peaks, and the probability of this tunneling behavior is dependent on the size of the barrier and the amount of energy the particle has.

Essentially, if the particle has a lot of energy, it is moving fast. If it is moving fast, it can tunnel through the potential barrier. If a particle's aggregate phase of matter has a high temperature, each particle has a lot of kinetic energy, on average. Let's reconsider our analogy with simulated annealing. If we have particles with high temperature, we allow them to escape local minima by going through otherwise unfavorable portions of state space. They are, in this analogy, 'tunneling' through the regions of low probability.

Now, as to how this relates to the distributions, I'm not sure.

• The thing is that as long as the temperature is non zero and the proposal distribution is such that the system is ergodic (symmetric uniform is enough) the particle will tunnel given enough time. Metropolis, at least in its original form, is good to describe the physics at equilibrium, how it gets there seems nothing more than a numeric method with little resemblance to the system dynamics. Reading Roberts et al, 1997 though it seems that Langevin diffusion is the key to physically understand the non-equilibrium dynamics, and other sources seem to confirm it. Feb 13, 2016 at 16:55
• If a particle has a lot of energy, it is passing above the barrier, not tunneling through it. Quantum tunneling does serve as a basis for some more advanced algorithms (see, e. g., those listed in the Wikipedia article in simulated annealing). Apr 17, 2022 at 8:14

In the original paper, and most derivations I found, it seems that it's not much important as long as it's symmetric (for detailed-balance) and allows a good exploration of the configuration space (e.g. ∼0.5

acceptance ratio). Usually it's either a uniform or a normal distribution.

Now when studying Simulated Annealing, most papers assume a common use of a normal distribution with σ=T .From a physicist point of view simulated annealing seeks the most probable state of a system, described by the Boltzmann distribution: $$w(p,q|T) \propto \exp\left( -\frac{H(p,q)}{T}\right)$$ where $$H(p,q)$$ is the Hamiltonian, also referred to as the energy of the system, whereas $$T$$ is temperature. Generalized momenta and coordinates, $$p$$ and $$q$$ are related via Hamiltonian evolution (heavily exploited, e.g., in the Hamiltonian Monte Carlo); these need not he continuous - e.g., they could by discrete up/down states of magnetic moments in ferromagnet or placement of atoms A and B in a binary alloy. MCMC is used to sample different variable configurations not connected via Hamiltonian evolution - physically it mostly means sampling states of different energies. This does not imply a gaussian proposal density and largely depends on the nature of the variables and the Hamiltonian.

The essential part is that at the temperature of interest, $$T_0$$ the system us likely to be found in its lowest energy state. Finding this state corresponds to maximizing the probability/likelihood $$w(p,q|T_0)$$.

In case of a multimodal likelihood we are not guaranteed to find the global maximum - physically it means that there may exist multiple global maxima (e g., corresponding to different directions of magnetization in a ferromagnet) or there may be some metastable/long-living "glass" states, the transitions from which to the global maximum are unlikely, since they require passing through highly improbable intermediate configurations.

As the probability of a state is proportional to $$\exp\left( -\frac{H(p,q)}{T}\right)$$, intermediate states are easier reached, if the temperature $$T$$ is higher. The solution is then to explore the energy states at a high temperature, and then use the result to redo the exploration at lower and lower temperature, till we reach the desired one.

Probabilistically, instead of exploring distribution $$p(x)\propto \exp\left(-l(x)\right),$$ we work with a rescaled negative log-probability $$p_\alpha(x)\propto \exp\left(-\alpha l(x)\right),\alpha \leq 1,$$ which is more easily amenable to exploration, e.g., via the Metropolis-Hastings algorithm.

One particular class of physical problems is particles undergoing diffusive/Langevin dynamic. These can be modeled as a random walk with a gaussian proposal density. They are referred to sometimes as drift-diffusion problems, but they are rather ubiquitous in science and may appear under many different names. The easy "physical" way of thinking of such systems is water on a rugged surface, which is being shaked. When the holes in the surface are deep, while the shaking is weak, the water remains where it is - the shaking explores only the surface near the minima. Shaking harder would allow the water to leak to other places. Annealing achieves it by rescaling the rugged landscape, allowing the water flow, and then gradually restoring the original size of the roughness.

Update

In the original paper, and most derivations I found, it seems that it's not much important as long as it's symmetric (for detailed-balance) and allows a good exploration of the configuration space (e.g. ∼0.5 acceptance ratio). Usually it's either a uniform or a normal distribution.

Now when studying Simulated Annealing, most papers assume a common use of a normal distribution with σ=T

Near the probability maximum (energy minimum) log-likelihood can be considered as parabolic, i.e., we have $$p(x)\propto \exp\left(-\frac{kx^2}{T}\right),$$ which means that the characteristic extent of the distribution is $$\sigma=\sqrt{T/k}$$. This gives the idea for the suitable step of the Metropolis chain to get the acceptance rate around 0.5: if the proposal distribution is much wider, most steps will be rejected, while if it is much narrower, most of them will be accepted.

As it is seen from the above argument, we are not talking here literally about $$\sigma=\sqrt{T}$$, but rather about the proposal distribution width scaling with temperature as $$\sqrt{T}$$ - there is still a factor that depends on the curvature of the probability distribution near its peak.

References

1. Simulated annealing: Theory and applications (1987) by Laarhoven and Aarts is an early book taking simulated annealing from physics domain to general statistical applications.
2. Reaction rate: fifty years after Kramers (1990) by Hänggi et al. is a review article describing many physical, chemical and other systems, where diffusion-like dynamics leads to transitions between local energy minima (probability maxima).
3. Probabilistic inference using Markov chain Monte Carlo methods by Neal specifically discusses the size of the annealing step.