Nelder-Mead Method

Nelder-Mead Method for Calibrating Health Economic Models

The Nelder-Mead method is a direct search routine that finds parameter values giving the lowest score on an objective, such as the distance between model outputs and calibration targets or a negative log-likelihood. It compares the score at a small set of trial points and keeps replacing the worst one, so it needs only the ability to run the model and read off a number. That property explains its long use in health economics, where natural history parameters of a Markov model or a microsimulation often cannot be observed directly and must be inferred by calibration. It is one of several forms of numerical optimisation and is unrelated to the simplex method of linear programming.

How the simplex moves towards a lower score

For a problem with $n$ unknown parameters, the method keeps a simplex of $n + 1$ points, a triangle when there are two. Nelder and Mead described it in 1965 as a method for minimising a function of $n$ variables in which the vertex with the highest value is replaced by another point, so that the simplex adapts to the local shape of the function. Lagarias, Reeds, Wright and Wright later stated the algorithm precisely, adding tie-breaking rules the original paper left open, and the rules below follow them.

Each iteration orders the vertices so that $x_1$ has the lowest score and $x_{n+1}$ the highest, then computes the centroid of all vertices except the worst.

$$\bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i$$

where $\bar{x}$ is the centroid, $x_i$ is the $i$th best vertex and $n$ is the number of parameters. Every candidate point in a normal iteration lies on the line through the worst vertex and this centroid.

$$\begin{aligned} x_r &= \bar{x} + \rho\,(\bar{x} - x_{n+1}) \\ x_e &= \bar{x} + \chi\,(x_r - \bar{x}) \\ x_c &= \bar{x} + \gamma\,(x_r - \bar{x}) \\ x_{cc} &= \bar{x} - \gamma\,(\bar{x} - x_{n+1}) \\ v_i &= x_1 + \sigma\,(x_i - x_1) \end{aligned}$$

where $x_r$ is the reflected point, $x_e$ the expanded point, $x_c$ the outside contraction point, $x_{cc}$ the inside contraction point and $v_i$ the shrunk position of vertex $i$. The coefficients are $\rho$ for reflection, $\chi$ for expansion, $\gamma$ for contraction and $\sigma$ for shrinkage, and Lagarias and colleagues describe $\rho = 1$, $\chi = 2$, $\gamma = 1/2$ and $\sigma = 1/2$ as the nearly universal standard choice.

The choice among these points follows the ordering of scores. Writing $f_1$, $f_n$ and $f_{n+1}$ for the best, next-worst and worst scores and $f_r$ for the score at $x_r$, the iteration runs as follows.

1. Reflect. If $f_1 \le f_r < f_n$, the reflected point replaces the worst vertex.
2. Expand. If $f_r < f_1$, the expanded point is kept if it beats $x_r$, and $x_r$ is kept if it does not.
3. Contract outside. If $f_n \le f_r < f_{n+1}$, the point $x_c$ is tried and kept if its score is no worse than $f_r$.
4. Contract inside. If $f_r \ge f_{n+1}$, the point $x_{cc}$ is tried and kept if it beats the worst vertex.
5. Shrink. If a contraction fails, every vertex except the best moves halfway towards $x_1$.

A normal iteration costs one or two model runs, and a shrink costs $n$ more. Each run can be a full simulation of a cohort or of many individuals, so evaluations matter more than iterations.

Uses in calibration and likelihood fitting

Stout and colleagues named the Nelder-Mead algorithm among the directed searches used in cancer simulation models. Taylor and colleagues compared three ways of calibrating a Markov model of human papillomavirus infection and cervical disease to US data. Nelder-Mead gave the best fit, a mean deviation of 7% from the targets against 10% for manual calibration, using substantially less analyst time but more computing, while random search reached 39%.

Nelder-Mead is the default method of the optim function in R, described in its documentation as using only function values and as relatively slow but reasonably effective for non-differentiable functions. The calibrate_model function in the heemod package for Markov models lists Nelder-Mead first among its methods and accepts a table of starting values, running once per row so that results from different starts can be compared.

The method also fits statistical models whose likelihood is awkward to differentiate, for example a parametric survival model estimated by maximum likelihood estimation with a custom likelihood. The 1965 paper already gave a procedure for estimating the Hessian matrix near the minimum, which statistical estimation needs for standard errors. In R, optim can return a numerically differentiated Hessian at the solution found, which is needed because the simplex itself carries no curvature information.

Choices that change the answer

The analyst supplies a starting point and, directly or through defaults, the size of the first simplex. A small simplex searches only near the start, and a start in the wrong basin can end in a local minimum that looks like a good fit.

Scaling matters because the moves treat all directions alike. A transition probability near 0.01 and a scale parameter measured in years differ by orders of magnitude, so a simplex sized for one is badly sized for the other. R provides a parscale control for this, and many calibrations instead search on log or logit scales, which also keeps probabilities inside zero and one.

Stopping rules differ between implementations. R stops when a step cannot reduce the score by a relative tolerance, reltol, of about one part in a hundred million by default, or after 500 iterations for Nelder-Mead. SciPy uses absolute tolerances on the parameters and the score and defaults to 200 iterations per parameter. A run that ends on the iteration limit has not converged, and R reports it with code 1, while code 10 reports a degenerate simplex.

Bounds are handled crudely or not at all, because the basic method is for unconstrained problems. In R, supplying non-trivial bounds to optim switches the method to L-BFGS-B with a warning, and in SciPy bounds on the Nelder-Mead method simply clip the vertices. Transforming the parameters or a constrained optimisation method is usually cleaner. For a single parameter, R warns that the default method works poorly and points to the Brent method.

Worked example fitting a Weibull curve to two landmark survival figures

All numbers in this example are illustrative. A model needs a survival curve for a patient group whose only published evidence is overall survival of 80% at 5 years and 60% at 10 years. The analyst assumes a Weibull distribution and calibrates its two parameters to these two figures.

$$S(t) = \exp\left[-\left(\frac{t}{\lambda}\right)^{k}\right]$$

where $S(t)$ is the proportion alive at year $t$, $k$ is the shape parameter and $\lambda$ is the scale parameter in years. The score is the sum of squared gaps in percentage points.

$$f(k, \lambda) = \left[100\,(S(5) - 0.80)\right]^{2} + \left[100\,(S(10) - 0.60)\right]^{2}$$

where $f$ is the goodness of fit score to be minimised. The standard coefficients are used throughout.

1. Starting simplex. The vertices are A = (1.0, 20), B = (1.5, 20) and C = (1.0, 25), written as ($k$, $\lambda$). At A, $S(5) = \exp(-0.25) = 0.7788$ and $S(10) = \exp(-0.5) = 0.6065$, so the gaps are 2.12 and 0.65 points and the score is 4.49 plus 0.43, or 4.92. The same working gives 172.48 at B and 52.96 at C, so A is best, C next-worst and B worst.

2. Centroid and reflection. The centroid of A and C is (1.0, 22.5). Reflecting B through it gives $x_r = 2 \times (1.0, 22.5) - (1.5, 20) = (0.5, 25.0)$, with survival of 0.6394 and 0.5313 and a score of 305.12.

3. Inside contraction. The reflected score is worse than the worst vertex, so the inside point is tried, $x_{cc} = (1.0, 22.5) - 0.5 \times [(1.0, 22.5) - (1.5, 20)] = (1.25, 21.25)$. Its score of 83.49 beats 172.48, so it replaces B.

4. The next five iterations. Each row starts from the simplex left by the row before.

At iteration 5 the reflected score of 4.83 beat the best score of 4.92, so the expansion to (1.2422, 14.1797) was tried, but it scored 74.89 and the reflected point was kept. After six iterations and 14 model evaluations the best vertex gives survival of 79.5% at 5 years and 61.7% at 10 years.

5. Convergence. After 30 iterations and 60 evaluations the best vertex is $k = 1.1948$, $\lambda = 17.546$, with a score below one millionth. The exact solution follows from the two targets. Since $(10/5)^k = \ln 0.60 / \ln 0.80 = 2.2892$, the shape is $k = \ln 2.2892 / \ln 2 = 1.1949$, and $\lambda = 5 / (-\ln 0.80)^{1/k} = 5 / 0.2231^{1/1.1949} = 17.545$.

6. Interpretation. The search found the right answer quickly, but the problem is easy. Two parameters meet two targets exactly and the solution is unique. A zero score shows only that the curve passes through both points. Another two-parameter family through the same points would extrapolate differently beyond 10 years, and that choice belongs to model structure rather than to the optimiser.

What convergence does and does not show

Theory offers little reassurance. Lagarias and colleagues proved convergence to a minimiser for strictly convex functions of one variable and only limited results for two. They cite a counterexample by McKinnon, a family of strictly convex functions of two variables and starting simplices for which the method converges to a point that is not a minimiser. For general nonconvex functions, they conclude, it does not converge to a minimiser, and the open question is why it tends to work well in practice by producing a rapid initial fall in the score.

A run that stops should be restarted from its best point with a fresh simplex, since an answer that moves on restart was not a minimum, and runs from several dispersed starting points show whether the search keeps finding the same basin. The convergence of the algorithm should also be kept apart from identification, since a model with more free parameters than informative targets can fit equally well at many parameter sets.

Stochastic models add a further problem. In a microsimulation each score carries Monte Carlo error, and because the method acts only on the ranking of vertex scores, noise can reverse a ranking and send the simplex the wrong way. Common random numbers for every evaluation, or enough simulated individuals that the noise is small beside real differences in fit, reduce this.

Nelder-Mead beside other search methods

Methods that use derivatives usually reach a smooth minimum in fewer evaluations. Gradient descent steps against the slope, the Newton-Raphson method adds curvature and the BFGS method in R uses function values and gradients. Nelder-Mead suits problems with a handful of parameters where gradients are unavailable or unreliable, and SciPy presents its adaptive option, which sets the coefficients by dimension, as useful for high-dimensional minimisation.

Grid search evaluates a fixed lattice and grows rapidly with the number of parameters. Simulated annealing and a genetic algorithm use random moves to escape local minima, and Bayesian optimisation fits a surrogate to earlier runs when each run is slow. The Solver add-in for Excel does not offer Nelder-Mead. Microsoft lists GRG Nonlinear for smooth nonlinear problems, Simplex LP for linear problems and Evolutionary, based on genetic algorithms, for non-smooth problems, and Simplex LP is the linear programming method that Lagarias and colleagues distinguish from Nelder-Mead.

Common misreadings

Most errors come from treating the point where the search stopped as more than it is. The analyst still has to show that the stopping point is a credible answer.

1. Reading a success code as proof of a global minimum. The code means only that a stopping rule was met, so restarts and multiple starting points are still needed.
2. Treating a near-zero score as a validated model. A good fit to calibration targets leaves extrapolation and structural choices untested.
3. Ignoring scale. Parameters of very different magnitude distort the simplex unless they are rescaled or transformed.
4. Expecting standard errors from the simplex. Uncertainty needs a separate Hessian estimate or a probabilistic or Bayesian calibration.

Sources

• Lagarias JC, Reeds JA, Wright MH, Wright PE. Convergence properties of the Nelder-Mead simplex method in low dimensions. SIAM Journal on Optimization. 1998;9(1):112-147. https://doi.org/10.1137/S1052623496303470
• Microsoft. Define and solve a problem by using Solver. Microsoft Support, accessed 26 October 2026. https://support.microsoft.com/en-us/office/define-and-solve-a-problem-by-using-solver-5d1a388f-079d-43ac-a7eb-f63e45925040
• Nelder JA, Mead R. A simplex method for function minimization. The Computer Journal. 1965;7(4):308-313. https://doi.org/10.1093/comjnl/7.4.308
• R Core Team. optim: General-purpose Optimization. R documentation, package stats, accessed 26 October 2026. https://stat.ethz.ch/R-manual/R-devel/library/stats/html/optim.html
• SciPy community. minimize(method='Nelder-Mead'). SciPy v1.18.0 manual, accessed 26 October 2026. https://docs.scipy.org/doc/scipy/reference/optimize.minimize-neldermead.html
• Stout NK, Knudsen AB, Kong CY, McMahon PM, Gazelle GS. Calibration methods used in cancer simulation models and suggested reporting guidelines. PharmacoEconomics. 2009;27(7):533-545. https://doi.org/10.2165/11314830-000000000-00000
• Taylor DC, Pawar V, Kruzikas D, Gilmore KE, Pandya A, Iskandar R, Weinstein MC. Methods of model calibration: observations from a mathematical model of cervical cancer. PharmacoEconomics. 2010;28(11):995-1000. https://doi.org/10.2165/11538660-000000000-00000
• Zarca K, Filipovic-Pierucci A, et al. heemod: Markov Models for Health Economic Evaluations. R package version 1.1.0, CRAN, 2025. Reference manual, calibrate_model. https://cran.r-project.org/web/packages/heemod/refman/heemod.html