Simulated Annealing

Simulated Annealing in Health Economic Model Calibration

Simulated annealing is a random search method that health economists use to calibrate disease models, finding values for unobservable natural history parameters whose outputs best match observed data. At each step it proposes a nearby candidate, always moves when the fit improves, and sometimes moves to a worse fit with a probability that shrinks as a control parameter called the temperature falls. Those uphill moves let the search climb out of a poor valley in the fit score, which a method that only ever improves cannot do. It is a heuristic global method within numerical optimisation and needs only the model's score for each candidate, with no gradient. The price is many model runs and an answer that is the best candidate found in the time allowed, with no proof of global optimality and no measure of uncertainty.

Physical annealing and the Metropolis acceptance rule

Kirkpatrick, Gelatt and Vecchi took the idea from statistical mechanics. A crystal is grown by careful annealing, first melting the substance and then lowering the temperature slowly. If it is cooled too fast, the result has many defects or forms a glass with only metastable, locally optimal structures. They argued that ordinary iterative improvement, which accepts only rearrangements that lower the cost function, is like extremely rapid quenching to a temperature of zero, so its solutions are usually metastable.

The bridge between physics and optimisation is the Metropolis procedure, in which the cost function plays the role of energy. A proposed change that lowers the cost or leaves it unchanged is accepted. A change that raises it is accepted with a probability that falls as the rise grows and as the temperature drops.

$$P(\text{accept}) = \min\left\{1, \exp\left(-\frac{\Delta E}{T}\right)\right\}, \quad \Delta E = f(x') - f(x)$$

where $x$ is the current candidate, $x'$ is the proposed neighbour, $f$ is the objective to be minimised, $\Delta E$ is the change in the objective, and $T$ is the temperature, a control parameter in the same units as the cost function. In practice a random number is drawn uniformly between 0 and 1 and the worse candidate is kept if the number is less than the acceptance probability.

Kirkpatrick and colleagues observed that the temperature sorts moves by size, with large changes in the objective at high temperatures and small ones deferred until low temperatures.

How the search runs in a calibration

In model calibration, the candidate is a vector of parameter values and the objective is a goodness-of-fit score comparing model outputs with calibration targets. A run then repeats four actions until a stopping rule is met.

1. Propose a neighbour by a small random change within the allowed ranges.
2. Run the model and score it, which costs one full model run.
3. Accept or reject with the Metropolis rule.
4. Lower the temperature by the cooling schedule, recording the best set seen so far.

A geometric schedule multiplies the temperature by a fixed ratio after each step or each stage.

$$T_{k+1} = \alpha T_k, \quad 0 < \alpha < 1$$

where $T_k$ is the temperature at step or stage $k$ and $\alpha$ is the cooling ratio. Kirkpatrick and colleagues, splitting the circuits of a microprocessor between two chips, started at a temperature of 10, at which nearly all proposed moves were accepted, and cooled with a ratio of 0.9 between stages.

What convergence theory guarantees

Hajek gave a necessary and sufficient condition on the cooling schedule for the state of the algorithm to converge in probability to the set of global minima. For a logarithmic schedule the condition takes a simple form.

$$T(t) = \frac{c}{\log(1 + t)}$$

where $T(t)$ is the temperature at step $t$ and $c$ is a constant. Hajek showed that convergence to the global minima requires $c$ to be at least the depth of the deepest local minimum that is not a global minimum, depth being roughly how far the score must rise before a better state can be reached.

The guarantee is asymptotic and the schedule is very slow. With an illustrative $c = 2.5$, the temperature is still 0.18 after a million steps, because $c / \log(1 + 10^6) = 2.5 / 13.8 = 0.18$. Calibrations limited to a few thousand model runs use much faster geometric cooling, outside the theorem, so repeated runs from different starting points are the practical check.

Calibrating a lung cancer microsimulation

Kong and colleagues applied simulated annealing to the Lung Cancer Policy Model, a microsimulation built to evaluate screening and other lung cancer control interventions. They searched 28 natural history parameters at once against 11 registry and cohort targets combined into one weighted sum. Each iteration simulated 500,000 people, and 1,000 iterations took about 400 hours of processor time.

The search was started from 10 random points, each run limited to 1,000 iterations or to reaching a preset stopping score. Over 10 runs, simulated annealing reached a mean best score of 86.5 (range 80.7 to 90.6), against 91.3 (85.2 to 94.1) for a genetic algorithm, a significant difference. The authors judged simulated annealing the better choice for automated calibration.

Calibrating transition probabilities for an alcohol model

Barbosa and colleagues used simulated annealing to estimate long-term drinking patterns for treatment-seeking people with alcohol use disorder, so that cost-effectiveness analyses of alcohol treatment could extrapolate beyond short follow-up. Short-term transition probabilities between abstinence and four World Health Organization drinking risk levels were estimated with multinomial logit models fitted to a national survey. Simulated annealing then calibrated the transition probabilities to a matched sample from a separate 16-year study of alcohol consumption.

Each proposal picked one transition at random and raised or lowered it by one percentage point in all 24 demographic transition matrices, adjusting the rest of the row proportionally so every row still summed to 1 with no negative entries. The temperature started at 10 and was multiplied by 0.995 at each iteration, and the search stopped when it fell below 0.01. The distance from the target, a sum of squared differences between predicted and observed state distributions, fell from 84.69 to 0.01 for the second set of transition probabilities and from 39.60 to 0.01 for the third. The design shows how the neighbourhood can be built to keep every candidate a valid Markov model input.

Worked example with a local minimum in a calibration score

All numbers in this example are illustrative. A cohort model has one unknown annual probability of progression, $p$, searched over seven values. The goodness-of-fit score for each value is shown below, lower being better. There is a local minimum at 0.03, a barrier at 0.04 and the global minimum at 0.06.

A neighbour is the next value up or down, chosen with equal probability.

1. Greedy search stalls. A search that only accepts improvements, started at 0.03, rejects both neighbours (15.0 and 13.5 are worse than 11.0) and stops at a score of 11.0.

2. Uphill move at high temperature. At $T = 5$, the move to 0.04 has $\Delta E = 13.5 - 11.0 = 2.5$, so the acceptance probability is $e^{-2.5/5} = e^{-0.5} = 0.607$. The uniform random draw is 0.42, which is below 0.607, so the worse value is accepted.

3. Downhill moves are always accepted. From 0.04, the move to 0.05 has $\Delta E = 9.5 - 13.5 = -4.0$, and the next move to 0.06 has $\Delta E = 6.0 - 9.5 = -3.5$. Both are accepted and the score reaches 6.0, the global minimum.

4. Cooling locks the search in. At $T = 0.5$, an uphill move from 0.06 to 0.07 has probability $e^{-2.0/0.5} = e^{-4} = 0.018$, and a move to 0.05 has $e^{-3.5/0.5} = e^{-7} = 0.0009$. The search now rarely leaves the valley it is in.

5. Timing matters. The same barrier crossing from 0.03 to 0.04 at $T = 0.5$ would have probability $e^{-2.5/0.5} = e^{-5} = 0.0067$. A search still in the local valley when the temperature is this low will almost certainly stay there.

6. Count the iterations. With a cooling ratio of 0.995 and a starting temperature of 5, the temperature first drops below 0.5 after 460 iterations, since $\log(0.1) / \log(0.995) = 459.4$. With the settings Barbosa and colleagues used, a start at 10 and a stop below 0.01, the run lasts 1,379 iterations, since $\log(0.001) / \log(0.995) = 1378.1$. In a microsimulation each of those iterations is a full model run.

7. Interpret. The search reports $p = 0.06$, a point estimate. Other values with scores close to 6.0, here 0.07 at 8.0, may be almost as consistent with the targets, and the run says nothing about how uncertain $p$ is.

Choices that change the result

Settings fixed before the search shape the result. Kirkpatrick and colleagues noted that a schedule may be developed by trial and error or set by warming the system until it is clearly melted and then cooling in slow stages. A starting temperature that is too low turns the method into a greedy search, while one that is too high wastes model runs on a random walk. A cooling ratio close to 1 explores more thoroughly at higher cost.

The neighbourhood and the bounds shape what can be found. Kong and colleagues noted as a limitation that their best scores were found within estimated allowable ranges and that lower scores might lie outside them. The weights used to combine several targets into one score also change which parameter set counts as best, so the weighting is part of the model specification.

Stochastic models add noise. Menzies and colleagues point out that optimisation is particularly difficult for stochastic models, whose outputs for a given parameter set carry Monte Carlo error, so an accepted move may reflect noise more than a real change in fit. A larger simulated sample per run reduces that error, and common random numbers across candidates reduce the noise in score differences.

Reading a calibrated parameter set in an economic evaluation

Simulated annealing returns one best-fitting parameter set. Menzies and colleagues note that results calculated from parameter point estimates provide no uncertainty estimates for modelled outcomes, unlike probabilistic sensitivity analysis. A cost-effectiveness model that carries forward only the annealing result treats calibrated parameters as known and leaves their parameter uncertainty out. Carrying forward several good-fitting sets, or using a Bayesian calibration, keeps that uncertainty in.

Menzies and colleagues add that if optimisation finds different best-fitting sets from different starting points, this suggests the posterior is multimodal or otherwise difficult to sample. Repeat runs from random starts, as in Kong's 10 searches, are therefore a diagnostic as well as a check on the optimiser. The ISPOR Optimization Methods Emerging Good Practices Task Force classes simulated annealing among heuristic approaches, which give approximate solutions when an exact approach is unavailable or computationally expensive. A calibration report should give the schedule, neighbourhood, bounds, stopping rule, number of runs and spread of their results.

Simulated annealing beside other search methods

The acceptance rule is the Metropolis rule, the symmetric-proposal case of the Metropolis-Hastings algorithm, but the purpose differs. Markov chain Monte Carlo holds the target distribution fixed and keeps sampling to describe a posterior, whereas simulated annealing lowers the temperature so the chain concentrates on the minimum and the sample path is discarded.

Among optimisers, the Nelder-Mead method moves a simplex of points downhill and can stall in a local minimum. A genetic algorithm and particle swarm optimisation move a population of candidates, while simulated annealing moves a single one. Grid search scores a fixed lattice that grows quickly with the number of parameters, and Bayesian optimisation fits an emulator to choose each next run. Problems with linear objectives and constraints have exact methods in linear programming, so a heuristic earns its place when no exact method fits, as when the model must be run to score each candidate.

Limitations and common misreadings

The commonest misreading treats the result as the global optimum. Hajek's guarantee needs a very slow schedule, such as logarithmic cooling with a large enough constant, run over an unbounded number of steps, which no practical calibration uses. A low score from one run shows a good fit was found, not that no better fit exists.

A second misreading treats a good fit as identification. Several parameter sets may match the targets almost equally well, and the single best set hides that. The calibrated values are only as credible as the targets, the weights and the model structure behind them.

Sources

• Barbosa C, Dowd WN, Aldridge AP, Timko C, Zarkin GA. Estimating long-term drinking patterns for people with lifetime alcohol use disorder. Medical Decision Making. 2019;39(7):765-780. https://doi.org/10.1177/0272989X19873627
• Crown W, Buyukkaramikli N, Thokala P, et al. Constrained optimization methods in health services research: an introduction. Report 1 of the ISPOR Optimization Methods Emerging Good Practices Task Force. Value in Health. 2017;20(3):310-319. https://doi.org/10.1016/j.jval.2017.01.013
• Hajek B. Cooling schedules for optimal annealing. Mathematics of Operations Research. 1988;13(2):311-329. https://doi.org/10.1287/moor.13.2.311
• Kirkpatrick S, Gelatt CD, Vecchi MP. Optimization by simulated annealing. Science. 1983;220(4598):671-680. https://doi.org/10.1126/science.220.4598.671
• Kong CY, McMahon PM, Gazelle GS. Calibration of disease simulation model using an engineering approach. Value in Health. 2009;12(4):521-529. https://doi.org/10.1111/j.1524-4733.2008.00484.x
• Menzies NA, Soeteman DI, Pandya A, Kim JJ. Bayesian methods for calibrating health policy models: a tutorial. PharmacoEconomics. 2017;35(6):613-624. https://doi.org/10.1007/s40273-017-0494-4