A detailed comparison of MIGRAD and Optim.BFGS() turned out to be one prompt away. I found the result interesting and instructive, so I am posting it as a follow-up of our discussion.
How this note was made
Model: Codex Sol 5.6, reasoning level High
Project context: NativeMinuit.jl + Optim.jl
Prompt:
I have been fascinated by the numerical algorithm of Minuit2 and stability it has. Optim.BFGS() should use a similar algorithm, but it seems not to be well tuned for the problems I solve. I would like a compact explanation of how the algorithms work, how similar they are, and where they diverge. Study the implementations and check the existing setting mappings in FitterHEP.jl.
The analysis used the ROOT Minuit2, Optim.jl, NativeMinuit.jl, and FitterHEP.jl sources. The matrix-update identities were also checked numerically. The principal local versions were Optim 2.2.1 and NativeMinuit 0.6.2.
Common quasi-Newton iteration
Both algorithms maintain an approximation M_k to the inverse Hessian,
M_k \approx [\nabla^2 f(x_k)]^{-1},
and form the search direction
p_k=-M_k g_k,\qquad g_k=\nabla f(x_k).
A line search selects \alpha_k, after which
x_{k+1}=x_k+\alpha_kp_k.
The curvature update uses
s=x_{k+1}-x_k,\qquad y=g_{k+1}-g_k,
and attempts to satisfy the secant equation
M_{k+1}y=s.
This is the common algorithmic structure of MIGRAD and Optim’s BFGS implementation.
BFGS and MIGRAD matrix updates
Optim’s BFGS implementation uses the inverse-BFGS update
M_{k+1}
=
M_k+
\frac{\delta+\gamma}{\delta^2}ss^\mathsf T
-
\frac{M_ky\,s^\mathsf T+s\,y^\mathsf TM_k}{\delta},
where
\delta=s^\mathsf Ty,\qquad
\gamma=y^\mathsf TM_ky.
Optim applies this update only when \delta>0; otherwise it retains the previous matrix.
Default MIGRAD first constructs the DFP update
M_{\mathrm{DFP}}
=
M_k+\frac{ss^\mathsf T}{\delta}
-\frac{(M_ky)(M_ky)^\mathsf T}{\gamma}.
When \delta>\gamma, it additionally applies
\gamma
\left(
\frac{s}{\delta}-\frac{M_ky}{\gamma}
\right)
\left(
\frac{s}{\delta}-\frac{M_ky}{\gamma}
\right)^\mathsf T.
Expanding this term gives exactly the inverse-BFGS formula. Therefore, default MIGRAD uses:
- the BFGS update when s^\mathsf Ty>y^\mathsf TM_ky;
- the DFP update otherwise.
This behaviour is explicit in ROOT’s DavidonErrorUpdator.cxx. ROOT also provides a separate pure-BFGS updater, but the standard MnMigrad constructor uses the Davidon updater.
I checked the identity numerically using non-collinear random secant pairs. For \delta/\gamma=3, the MIGRAD and BFGS matrices agreed to 1.9\times10^{-15}. For \delta/\gamma=0.5, the MIGRAD update agreed with DFP and differed from BFGS.
Thus, some MIGRAD iterations use exactly the BFGS matrix update. The standard MIGRAD algorithm as a whole is a DFP/BFGS hybrid.
Initial scaling
Optim’s default initial inverse-Hessian approximation is the identity matrix. It can be replaced through initial_invH or scaled through initial_stepnorm, as shown in the BFGS constructor and initialization code.
Minuit derives its initial metric from the user-supplied parameter steps. For an unbounded parameter with initial step e_i, Minuit estimates the diagonal curvature as
g_{2,i}=\frac{2\,\mathrm{Up}}{e_i^2},
where Up is errordef. The corresponding diagonal inverse-Hessian estimate is
(M_0)_{ii}\approx\frac{e_i^2}{2\,\mathrm{Up}}.
This construction appears in ROOT’s InitialGradientCalculator.cxx.
Minuit’s initial parameter errors therefore have two roles: they influence numerical differentiation and define the initial optimization metric. This matters when parameters have different units or characteristic scales.
An equivalent Optim initialization is
initial_invH =
_ -> Diagonal(step_sizes.^2 ./ (2errordef))
For a negative-log-likelihood objective with errordef = 0.5, this reduces to
initial_invH = _ -> Diagonal(step_sizes.^2)
For a \chi^2 objective with errordef = 1, the diagonal is instead step_sizes.^2 / 2.
Line search and numerical derivatives
Optim’s default BFGS method uses a Hager–Zhang line search. It tests sufficient-decrease and curvature conditions and may evaluate both the objective and directional derivative at several trial points.
MIGRAD uses a parabolic line search. It first evaluates the nominal quasi-Newton step and then uses two- and three-point quadratic interpolation. ROOT’s MnLineSearch.cxx limits this procedure to 12 line-search iterations.
This difference affects the cost of numerical gradients. A directional-derivative evaluation in Optim may require a complete numerical gradient. MIGRAD primarily evaluates the objective during its line search and recalculates the gradient after choosing a new point.
Minuit’s numerical gradient calculation also refines the finite-difference step separately for every parameter. The number of refinement cycles and their tolerances depend on the strategy level; the corresponding settings are defined in MnStrategy.cxx.
With accurate analytic or automatic gradients, the Hager–Zhang search provides stronger general-purpose line-search conditions. With an expensive objective and numerical derivatives, the Minuit procedure may require fewer objective evaluations per line search.
Stopping criteria
Optim normally tests the infinity norm of the gradient,
\lVert g\rVert_\infty \leq g_{\mathrm{abstol}},
with a default threshold of 10^{-8}. Changes in the parameters or objective can also terminate minimization. These conditions are implemented in Optim’s convergence assessment.
MIGRAD instead uses the Expected Distance to Minimum,
\mathrm{EDM}=\frac12g^\mathsf TMg.
Under the local quadratic model, EDM estimates the remaining decrease in the objective.
The effective MIGRAD target is
\mathrm{EDM}_{\mathrm{goal}}
=
0.002\,
\texttt{tolerance}\,
\mathrm{errordef},
subject to a machine-precision floor. During the variable-metric iteration, Minuit corrects EDM using its covariance-change estimate Dcovar,
\mathrm{EDM}_{\mathrm{corrected}}
=
\mathrm{EDM}\,(1+3\,\mathrm{Dcovar}).
Depending on the strategy and Dcovar, MIGRAD may recompute the Hessian and continue minimization from the refined state. This control flow is implemented in ROOT’s VariableMetricBuilder.cxx.
The EDM callback proposed earlier in this thread,
edm = dot(g, invH * g) / 2
therefore reproduces the principal Minuit convergence quantity when invH is Optim’s current BFGS approximation. It does not reproduce the Dcovar correction, Hessian recomputation, or subsequent MIGRAD restart.
There is no parameter-independent conversion from Minuit’s tolerance to Optim’s g_tol. A raw gradient threshold depends on parameter units and scaling, whereas EDM includes the local inverse-Hessian metric.
Matrix validation and recovery
Optim preserves positive definiteness by applying its BFGS update only when s^\mathsf Ty>0. If the search direction is not a descent direction, the inverse-Hessian approximation can be reset.
MIGRAD additionally:
- checks whether -Mg is a descent direction;
- repairs a non-positive-definite metric by shifting its eigenvalues;
- tracks the estimated covariance change through
Dcovar;
- can run HESSE when the variable-metric covariance remains uncertain;
- can resume the variable-metric iteration after HESSE.
These operations belong to the VariableMetricBuilder control flow rather than to the Davidon update itself.
Bounds
Minuit transforms bounded external parameters into unconstrained internal coordinates and performs MIGRAD in the internal space.
Optim normally handles box constraints through Fminbox, which adds a logarithmic barrier around an unconstrained inner optimizer. FitterHEP can alternatively use logistic and softplus transformations.
These methods enforce the same physical limits but do not define the same optimization problem in the working coordinates. A comparison involving bounds therefore includes differences in the transformation or barrier method, in addition to differences between MIGRAD and BFGS.
FitterHEP setting mappings
The FitterHEP Optim backend initializes Optim’s BFGS metric using
Diagonal(step_sizes .^ 2)
and supplies the corresponding diagonal Hessian preconditioner to LBFGS.
This agrees with Minuit for errordef = 0.5. The general mapping should include errordef:
initial_invH =
Diagonal(step_sizes.^2 ./ (2errordef))
and, for the LBFGS Hessian preconditioner,
P = Diagonal((2errordef) ./ step_sizes.^2)
For example, with steps [2, 3]:
errordef = 0.5 gives [4, 9];errordef = 1 gives [2, 4.5];errordef = 2 gives [1, 2.25].
Fixed parameters are treated consistently in intent: FitterHEP removes them from Optim’s active vector, while Minuit marks them as fixed.
The stopping options are not equivalent mappings. g_tol, x_tol, and f_tol retain their Optim meanings and do not correspond directly to Minuit’s EDM tolerance.
The current benchmarks also do not initialize the backends identically. In the mass-fit benchmark, Minuit receives the explicit steps
[0.01, 0.05, 0.05, 0.1]
while Optim receives generic default steps.
In addition, MinuitBackend enables its SIMPLEX fallback by default, whereas the Optim BFGS backend has no equivalent fallback. These are reasonable backend defaults, but they should be equalized when comparing the underlying minimizers.
Summary
The comparison gives the following conclusions:
- MIGRAD and
Optim.BFGS() use the same quasi-Newton iteration structure.
- Default MIGRAD switches between DFP and BFGS matrix updates.
- EDM can be calculated directly from Optim’s BFGS inverse-Hessian approximation.
- Minuit differs from plain BFGS through its initialization, numerical derivatives, line search, EDM correction, covariance monitoring, matrix repair, Hessian recomputation, and bound transformations.
- Minuit’s initial steps map to
(M_0)_{ii}=\frac{e_i^2}{2\,\mathrm{errordef}},
rather than to e_i^2 for every objective convention.
- A controlled comparison should use the same initial metric, derivative information, working coordinates, stopping rule, fallback policy, and objective-call budget.
I saw your messages on the