Optimizer Selection¶
Bumps has a number of different optimizers available, each with its own control parameters:
Random Lines [experimental]
Particle Swarm [experimental]
Parallel Tempering [experimental]
In general there is a trade-off between convergence rate and robustness, with the fastest algorithms most likely to find a local minimum rather than a global minimum. The gradient descent algorithms (Levenberg-Marquardt, Quasi-Newton BFGS) tend to be fast but they will find local minima only, while the population algorithms (DREAM, Differential Evolution) are more robust and likely slower. Nelder-Mead Simplex is somewhere between, with a small population keeping the search local but more robust than the gradient descent algorithms.
Each algorithm has its own set of control parameters for adjusting the search process and the stopping conditions. The same option may mean slightly different things to different optimizers. The Bumps package provides a dialog box for selecting the optimizer and its options when running the fit wx application. This only includes the common options for the most useful optimizers. For full control, the fit will need to be run from the command line interface or through a python script.
For parameter uncertainty, most algorithms use the covariance matrix at the optimum to estimate an uncertainty ellipse. This is okay for a preliminary analysis, but only works reliably for weakly correlated parameters. For full uncertainty analysis, DREAM uses a random walk to explore the parameter space near the minimum, showing pair-wise correlations amongst the parameter values. In order for DREAM to return the correct uncertainty, the function to be optimized should be a conditional probability density, with nllf as the negative log likelihood function of seeing point \(x\) in the parameter space. Other functions can be fitted, but uncertainty estimates will be meaningless.
Most algorithms have been adapted to run in parallel at least to some degree. The implementation is not heavily tuned, either in terms of minimizing the overhead per function evaluation or for distributing the problem across multiple processors. If the theory function is implemented in parallel, then the optimizer should be run in serial. Mixed mode is also possible when running on a cluster with a multi-threaded theory function. In this case, only one theory function will be evaluated on each cluster node, but the optimizer will distribute the parameters values to the cluster nodes in parallel. Do not run serial algorithms (Levenberg-Marquardt, Quasi-Newton BFGS) on a cluster.
We have included a number of optimizers in Bumps that did not perform particularly well on our problem sets. However, they may be perfect for your problem, so we have left them in the package for you to explore. They are not available in the GUI selection.
Levenberg-Marquardt¶
The Levenberg-Marquardt [1] [2] algorithm has been the standard method for non-linear data fitting. As a gradient descent trust region method, it starts at the initial value of the function and steps in the direction of the derivative until it reaches the minimum. Set up as an explicit minimization of the sum of square differences between theory and model, it uses a numerical approximation of the Jacobian matrix to set the step direction and an adaptive algorithm to set the size of the trust region.
When to use¶
Use this method when you have a reasonable fit near the minimum, and you want to get the best possible value. This can then be used as the starting point for uncertainty analysis using DREAM. This method requires that the problem definition includes a residuals method, but this should always be true when fitting data.
When modeling the results of an experiment, the best fit value is an accident of the measurement. Redo the same measurement, and the slightly different values you measure will lead to a different best fit. The important quantity to report is the credible interval covering 68% (1-\(\sigma\)) or 95% (2-\(\sigma\)) of the range of parameter values that are somewhat consistent with the data.
This method uses lmfit from scipy, and does not run in parallel.
Options¶
Steps is the number of gradient steps to take. Each step requires
a calculation of the Jacobian matrix to determine the direction. This
needs \(2 m n\) function evaluations, where \(n\) is the number of parameters and
each function is evaluated and \(m\) data points (assuming center point
formula for finite difference estimate of the derivative). The resulting
linear equation is then solved, but for small \(n\) and expensive function
evaluation this overhead can be ignored. Use --steps=n from
the command line.
f(x) tolerance and x tolerance are used to determine when
the fit has reached the point where no significant improvement is expected.
If the function value does not improve significantly within the step, or
the step is too short, then the fit will terminate. Use --ftol=v and
--xtol=v from the command line.
From the command line, --starts=n will automatically restart the algorithm
after it has converged so that a slightly better value can be found. If
--keep_best is included then restart will use a value near the minimum,
otherwise it will restart the fit from a random point in the parameter space.
Use --fit=lm to select the Levenberg-Marquardt fitter from the command line.
Notes¶
v0.8.2 Changed from scipy.leastsq to mpfit for better bounds handling. Use
--fit=scipy.leastsq to restore the previous behaviour.
References¶
Nelder-Mead Simplex¶
The Nelder-Mead [3] downhill simplex algorithm is a robust optimizer which does not require the function to be continuous or differentiable.
It uses the relative values of the function at the corners of a simplex (an n-dimensional triangle) to decide which points of the simplex to update. It will take the worst value and try moving it inward or outward, or reflect it through the centroid of the remaining values stopping if it finds a better value. If none of these values are better, then it will shrink the simplex and start again. The name amoeba comes from the book Numerical Recipes [4] wherein they describe the search as acting like an amoeba, squeezing through narrow valleys as it makes its way down to the minimum.
When to use¶
Use this method as a first fit to your model. If your fitting function is well behaved with few local minima this will give a quick estimate of the model, and help you decide if the model needs to be refined. If your function is poorly behaved, you will need to select a good initial value before fitting, or use a more robust method such as Differential Evolution or DREAM.
The uncertainty reported comes from a numerical derivative estimate at the minimum.
This method requires a series of function updates, and does not benefit much from running in parallel.
Options¶
Steps is the simplex update iterations to perform. Most updates
require one or two function evaluations, but shrinking the simplex evaluates
every value in the simplex. Use --steps=n from the command line.
Starts tells the optimizer to restart a given number of times.
Each time it restarts it uses a random starting point. Use
--starts=n from the command line.
Simplex radius is the initial size of the simplex, as a portion of
the bounds defining the parameter space. If a parameter is unbounded, then
the radius will be treated as a portion of the parameter value. Use
--radius=n from the command line.
x tolerance and f(x) tolerance are used to determine when the
fit has reached the point where no significant improvement is expected.
If the simplex is tiny (that is, the corners are close to each other) and
flat (that is, the values at the corners are close to each other),
then the fit will terminate. Use --xtol=v and --ftol=v from
the command line.
From the command line, use --keep_best so that restarts are centered on a
value near the minimum rather than restarting from a random point within the
parameter bounds.
Use --fit=amoeba to select the Nelder-Mead simplex fitter from the
command line.
References¶
Quasi-Newton BFGS¶
Broyden-Fletcher-Goldfarb-Shanno [5] is a gradient descent method which uses the gradient to determine the step direction and an approximation of the Hessian matrix to estimate the curvature and guess a step size. The step is further refined with a one-dimensional search in the direction of the gradient.
When to use¶
Like Levenberg-Marquardt, this method converges quickly to the minimum. It does not assume that the problem is in the form of a sum of squares and does not require a residuals method.
The \(n\) partial derivatives are computed in parallel.
Options¶
Steps is the number of gradient steps to take. Each step requires
a calculation of the Jacobian matrix to determine the direction. This
needs \(2 m n\) function evaluations, where \(n\) is the number of parameters and
each function is evaluated and \(m\) data points (assuming center point
formula for finite difference estimate of the derivative). The resulting
linear equation is then solved, but for small \(n\) and expensive function
evaluation this overhead can be ignored.
Use --steps=n from the command line.
Starts tells the optimizer to restart a given number of times.
Each time it restarts it uses a random starting point.
Use --starts=n from the command line.
f(x) tolerance and x tolerance are used to determine when
the fit has reached the point where no significant improvement is expected.
If the function is small or the step is too short then the fit
will terminate. Use --ftol=v and --xtol=v from the command line.
From the command line, --keep_best uses a value near the previous minimum
when restarting instead of using a random value within the parameter bounds.
Use --fit=newton to select BFGS from the command line.
References¶
Differential Evolution¶
Differential evolution [6] is a population based algorithm which uses
differences between points as a guide to selecting new points. For each member
of the population a pair of points is chosen at random, and a difference vector
is computed. This vector is scaled, and a random subset of its components are
added to the current point based on crossover ratio. This new point is
evaluated, and if its value is lower than the current point, it replaces
it in the population. There are many variations available within DE that
have not been exposed in Bumps. Interested users can modify
bumps.fitters.DEFit and experiment with different crossover and
mutation algorithms, and perhaps add them as command line options.
Differential evolution is a robust directed search strategy. Early in the search, when the population is disperse, the difference vectors are large and the search remains broad. As the search progresses, more of the population goes into the valleys and eventually all the points end up in local minima. Now the differences between random pairs will often be small and the search will become more localized.
The population is initialized according to the prior probability distribution for each each parameter. That is, if the parameter is bounded, it will use a uniform random number generate within the bounds. If it is unbounded, it will use a uniform value in [0,1]. If the parameter corresponds to the result of a previous measurement with mean \(\mu\) and standard deviation \(\sigma\), then the initial values will be pulled from a Gaussian random number generator.
When to use¶
Convergence with differential evolution will be slower, but more robust.
Each update will evaluate \(k\) points in parallel, where \(k\) is the size of the population.
Options¶
Steps is the number of iterations. Each step updates each member
of the population. The population size scales with the number of fitted
parameters. Use --steps=n from the command line.
Population determines the size of the population. The number of
individuals, \(k\), is equal to the number of fitted parameters times the
population scale factor. Use --pop=k from the command line.
Crossover ratio determines what proportion of the dimensions to update
at each step. Smaller values will likely lead to slower convergence, but
more robust results. Values must be between 0 and 1. Use --CR=v from
the command line.
Scale determines how much to scale each difference vector before adding
it to the candidate point. The selected mutation algorithm chooses a scale
factor uniformly in \([0,F]\). Use --F=v from the command line.
f(x) tolerance and x tolerance are used to determine when the
fit has reached the point where no significant improvement is expected.
If the population is flat (that is, the minimum and maximum values are
within tolerance) and tiny (that is, all the points are close to each
other) then the fit will terminate. Use ftol=v and xtol=v from the
command line.
Use --fit=de to select differential evolution from the command line.
References¶
DREAM¶
DREAM [7] is a population based algorithm like differential evolution, but instead of only keeping individuals which improve each generation, it will sometimes keep individuals which get worse. Although it is not fast and does not give the very best value for the function, we have found it to be a robust fitting engine which will give a good value given enough time.
The progress of each individual in the population from generation to generation can considered a Markov chain, whose transition probability is equal to the probability of taking the step times the probability that it keeps the step based on the difference in value between the points. By including a purely random stepper with some probability, the detailed balance condition is preserved, and the Markov chain converges onto the underlying equilibrium distribution. If the theory function represents the conditional probability of selecting each point in the parameter space, then the resulting chain is a random draw from the posterior distribution.
This means that the DREAM algorithm can be used to determine the parameter uncertainties. Unlike the hessian estimate at the minimum that is used to report uncertainties from the other fitters, the resulting uncertainty need not Gaussian. Indeed, the resulting distribution can even be multi-modal. Fits to measured data using theory functions that have symmetric solutions have shown all equivalent solutions with approximately equal probability.
When to use¶
Use DREAM when you need a robust fitting algorithm. It takes longer but it does an excellent job of exploring different minima and getting close to the global optimum.
Use DREAM when you want a detailed analysis of the parameter uncertainty.
Like differential evolution, DREAM will evaluate \(k\) points in parallel, where \(k\) is the size of the population.
Options¶
Samples is the number of points to be drawn from the Markov chain.
To estimate the 68% interval to two digits of precision, at least
1e5 (or 100,000) samples are needed. For the 95% interval, 1e6
(or 1,000,000) samples are needed. The default 1e4 samples
gives a rough approximation of the uncertainty relatively quickly.
Use --samples=n from the command line.
Burn-in Steps is the number of iterations to required for the Markov
chain to converge to the equilibrium distribution. If the fit ends
early, the tail of the burn will be saved to the start of the steps.
Use --burn=n from the command line.
Population determines the size of the population. The number of
individuals, \(k\), is equal to the number of fitted parameters times the
population scale factor. Use --pop=k from the command line.
Initializer determines how the population will be initialized. The options are as follows:
eps (epsilon ball), in which the entire initial population is chosen at random from within a tiny hypersphere centered about the initial point
lhs (latin hypersquare), which chops the bounds within each dimension in \(k\) equal sized chunks where \(k\) is the size of the population and makes sure that each parameter has at least one value within each chunk across the population.
cov (covariance matrix), in which the uncertainty is estimated using the covariance matrix at the initial point, and points are selected at random from the corresponding Gaussian ellipsoid
random (uniform random), in which the points are selected at random within the bounds of the parameters
Use --init=type from the command line.
Thinning is the amount of thinning to use when collecting the
population. If the fit is somewhat stuck, with most steps not improving
the fit, then you will need to thin the population to get proper
statistics. Use --thin=k from the command line.
Convergence gives a cutoff value \(\alpha\) for determining when
the Markov chain has converged. The default is --alpha=0.00 for no
convergence tests. Various tests are used, such as comparing the distribution
of points in the first part of the chain to the last part and looking for
trends in the log-likelihood values. You may need to use smaller \(\alpha\) for
shorter sequences (samples over variables times population) since the test
statistics will have higher variance. Convergence is tested every \(n\) steps.
Outliers is the test to use to check for outlier chains. Default is
--outliers=none for no outlier test. Options are iqr, which uses
the inter-quartile range on the likelihoods, grubbs, which uses a t-test
on the likelihoods, and mahal which looks at the distance from the best
chain in parameter space. Outlier removal occurs every \(2n\) steps where
\(n\) is #samples/(#pars #pop), or when the convergence test indicates the
chains are stable. Outliers are replaced by non-outlier chains at random.
These new chains need at least \(n\) steps to mix before being used. If the
MCMC exploration stops due to time, some of the chains may not be properly
mixed.
Burn-in trim is used to clear spurious samples from the Markov chains.
If --trim=true then Bumps finds the “burn point” after which the
chains appear to have converged. Samples before this point are ignored
when computed statistics and making plots. The trimmed samples are still
written to the MCMC output files so they will be available when the fit
is resumed.
Calculate entropy, if true, computes the entropy for the fit. This is
an estimate of the amount of information in the data. Use --entropy=method
from the command line, where method is one of llf (default), gmm, mvn
or wnn. See below for details.
Steps, if not zero, determines the number of iterations to use for
drawing samples after burn in. Each iteration updates the full population,
which is (population x number of fitted parameters) points. This option
is available for compatibility; it is more useful to set the number of
samples directly. Use --steps=n from the command line.
Use --fit=dream to select DREAM from the command line. Consider using
--parallel and --checkpoint as well. When running in a batch queue,
add --batch and use --mpi rather than --parallel.
Output¶
DREAM produces a number of different outputs, and there are a number of
things to check before using its reported uncertainty values. The main
goal of selecting --burn=n is to wait long enough to reach the
equilibrium distribution.
Fig. 142 This DREAM fit is incomplete, as can be seen on all four plots. The Convergence plot is still decreasing, the Parameter Trace plot shows a reduction in the mixing of Markov chain values, the Correlation plots are fuzzy and mostly empty, or show obvious correlations, and the Uncertainty plot shows black histograms (indicating that there are a few stray values far away from the best) whilst the green maximum likelihood spikes do not match the histogram (indicating that the region around the best value has not been adequately explored).¶
Fig. 143 This DREAM fit completed successfully. The Convergence plot is flat, the Parameter Trace plot is flat and messy indicating good mixing of the Markov chain values, the Correlation plots show nice defined blobs (with a bit of correlation between the M1.radius parameter and the M1.radius. width parameter), and the uncertainty plots show a narrow range of -log(P) values in the mostly brown histograms well-matched to the green constrained maximum likelihood line.¶
For each parameter in the fit, DREAM finds the mean, median and best value, as well as the 68% and 95% credible intervals. The mean value is defined as \(\int x P(x) dx\), which is just the expected value of the probability distribution for the parameter. The median value is the 50% point in the probability distribution, and the best value is the maximum likelihood value seen in the random walk. The credible intervals are the central intervals which capture 68% and 95% of the parameter values respectively. You need approximately 100,000 samples to get two digits of precision on the 68% interval, and 1,000,000 samples for the 95% interval. [9]
# |
Parameter |
mean |
median |
best |
[ 68% interval] |
[ 95% interval] |
|---|---|---|---|---|---|---|
1 |
M1.background |
0.059925(41) |
0.059924 |
0.059922 |
[0.05988 0.05997] |
[0.05985 0.06000] |
2 |
M1.radius |
2345.3(15) |
2345.234 |
2345.174 |
[2343.83 2346.74] |
[2342.36 2348.29] |
3 |
M1.radius.width |
0.00775(41) |
0.00774 |
0.00777 |
[ 0.0074 0.0081] |
[ 0.0070 0.0086] |
4 |
M1.scale |
0.21722(20) |
0.217218 |
0.217244 |
[0.21702 0.21743] |
[0.21681 0.21761] |
The Convergence plot shows the range of \(\chi^2\) values in the population for each iteration. The band shows the 68% of values around the median, and the solid line shows the minimum value. If the distribution has reached equilibrium, then convergence graph should be roughly flat, with little change in the minimum value throughout the graph. If there is no convergence, then the remaining plots don’t mean much.
The Correlations plot shows cross correlation between each pair of parameters. If the parameters are completely uncorrelated then the boxes should contain circles. Diagonals indicate strong correlation. Square blocks indicate that the fit is not sensitive to one of the parameters. The range plotted on the correlation plot is determined by the 95% interval of the data. The individual correlation plots are too small to show the range of values for the parameters. These can instead be read from the Uncertainty plot for each parameter, which covers the same range of values and indicates 68% and 95% intervals. If there are some chains that are wandering around away from the minimum, then the plot will look fuzzy, and not have a nice blob in the center. If a correlation plot has multiple blobs, then there are multiple minima in your problem space, usually because there are symmetries in the problem definition. For example, a model fitting \(x + a^2\) will have identical solutions for \(\pm\,a\).
The Uncertainty plot shows histograms for each fitted parameter generated from the values for that parameter across all chains. Within each histogram bar the values are sorted and displayed as a gradient from black to copper, with black values having the lowest \(\chi^2\) and copper values having the highest. The resulting histogram should be dark brown, with a black hump in the center and light brown tips. If there are large lumps of light brown, or excessive black then its likely that the optimizer did not converge. The green line over the histogram shows the best value seen within each histogram bin (the maximum likelihood given \(p_k == x\)). With enough samples and proper convergence, it should roughly follow the outline of the histogram. The yellow band in the center of the plot represents the 68% interval for the data. The histogram cuts off at 95%. These values along with the median are shown as labels along the x axis. The green asterisk represents the best value, the green E the mean value and the vertical green line the median value. If the fit is not sensitive to a parameter, or if two parameters are strongly correlated, the parameter histogram will show a box rather than a hump. Spiky shapes (either in the histogram or the maximum likelihood line) indicate lack of convergence or maybe not enough steps. A chopped histograms indicates that the range for that parameter is too small.
The Parameter Trace plot is diagnostic for models which have poor mixing. In this cases no matter how the parameter values are changing, they are landing on much worse values for the \(\chi^2\). This can happen if the problem is highly constrained with many tight and twisty values.
The Data and Theory plot should show theory and data lining up pretty well, with the theory overlaying about 2/3 of the error bars on the data (1-\(\sigma\) = 68%). The Residuals plot shows the difference between theory and data divided by uncertainty. The residuals should be 2/3 within [-1, 1], They should not show any structure, such as humps where the theory misses the data for long stretches. This indicates some feature missing from the model, or a lack of convergence to the best model.
If entropy is requested, then Bumps will show the total number of bits of information in the fit, where entropy is defined as:
Since we already have a sample from the posterior distribution \(p(\Theta)\) the Monte Carlo integral should be \(S \approx \sum_{k} \log_2 p(\theta_k)\). However, we do not know \(p(\theta_k)\), especially when we are integrating over nuisance parameters and only computing entropy for the parameters of interest. There are numerous methods in the literature for performing this calculation, and we have implemented the following:
gmm fits the MCMC sample to a Gaussian mixture model (GMM) and then estimates the entropy of the GMM through Monte Carlo integration.
llf finds the average ratio between the unnormalized negative log likelihood (NLLF) and a kernel density estimate (sklearn KernelDensity with default options), then estimates the entropy from the normalized likelihood through Monte Carlo integration. [8] This technique will not work for marginal likelihood estimates.
mvn fits the MCMC sample to a multivariate Gaussian and returns the entropy of that Gaussian. This is fast and accurate when the sample is well behaved (i.e., the uncertainty distribution is approximately Gaussian).
wnn estimates entropy from nearest-neighbour distances in the sample. [10]
Using entropy and simulation we hope to be able to make experiment planning decisions in a way that maximizes information, by estimating whether it is better to measure more precisely or to measure different but related values and fit them with shared parameters.
References¶
Particle Swarm¶
Inspired by bird flocking behaviour, the particle swarm [11] algorithm is a population-based method which updates an individual according to its momentum and a force toward the current best fit parameter values. We did not explore variations of this algorithm in any detail.
When to use¶
Particle swarm performed well enough in our low dimensional test problems, but made little progress when more fit parameters were added.
The population updates can run in parallel, but the tiny population size limits the amount of parallelism.
Options¶
--steps=n is the number of iterations. Each step updates each member
of the population. The population size scales with the number of fitted
parameters.
--pop=k determines the size of the population. The number of
individuals, \(k\), is equal to the number of fitted parameters times the
population scale factor. The default scale factor is 1.
Use --fit=ps to select particle swarm from the command line.
Add a few more lines
References¶
Random Lines¶
Most of the population based algorithms ignore the value of the function when choosing the points in the next iteration. Random lines [12] is a new style of algorithm which fits a quadratic model to a selection from the population, and uses that model to propose a new point in the next generation of the population. The hope is that the method will inherit the robustness of the population based algorithms as well as the rapid convergence of the newton descent algorithms.
When to use¶
Random lines works very well for some of our test problems, showing rapid convergence to the optimum, but on other problems it makes very little progress.
The population updates can run in parallel.
Options¶
--steps=n is the number of iterations. Each step updates each member
of the population. The population size scales with the number of fitted
parameters.
--pop=k determines the size of the population. The number of
individuals, \(k\), is equal to the number of fitted parameters times the
population scale factor. The default scale factor is 0.5.
--CR=v is the crossover ratio, determining what proportion of the
dimensions to update at each step. Values must be between 0 and 1.
--starts=n tells the optimizer to restart a given number of times.
Each time it restarts it uses a random starting point.
--keep_best uses a value near the previous minimum when restarting
instead of using a random value within the parameter bounds. This option is
not available in the options dialog.
Use --fit=rl to select random lines from the command line.
References¶
Parallel Tempering¶
Parallel tempering [13] is an MCMC algorithm for uncertainty analysis. This version runs at multiple temperatures simultaneously, with chains at high temperature able to more easily jump between minima and chains at low temperature to fully explore the minima. Like DREAM it has a differential evolution stepper, but this version uses the chain history as the population rather than maintaining a population at each temperature.
This is an experimental algorithm which does not yet perform well.
When to use¶
When complete, parallel tempering should be used for problems with widely spaced local minima which dream cannot fit.
Options¶
--steps=n is the number of iterations to include in the Markov
chain. Each iteration updates the full population. The population size
scales with the number of fitted parameters.
--burn=n is the number of iterations to required for the Markov
chain to converge to the equilibrium distribution. If the fit ends
early, the tail of the burn will be saved to the start of the steps.
--CR=v is the differential evolution crossover ratio to use when
computing step size and direction. Use a small value to step through the
dimensions one at a time, or a large value to step through all at once.
-nT=k, -Tmin=v and --Tmax=v specify a log-spaced initial
distribution of temperatures. The default is 25 points between
0.1 and 10. DREAM runs at a fixed temperature of 1.0.
Use --fit=pt to select parallel tempering from the command line.