Choosing an Optimizer for Econometric Estimation

2.1 Newton: use the observed curvature

Newton solves $Hs=-g$. For a positive-definite Hessian, this is the minimizer of the quadratic model. Near a well-behaved optimum it can converge rapidly; for a positive-definite quadratic, the exact unmodified step reaches the optimum in one iteration. An indefinite Hessian can instead produce an ascent direction, and a good local direction can still be too long for the nonlinear criterion.

The tested reference implementation is safeguarded Newton. It replaces Hessian eigenvalues below $10^{-5}$ by that floor and uses a line search to choose how far to move. The floor supplies a descent direction, but changes the step along flat directions. Its effects matter for the ill-conditioned OLS results; those results describe this implementation rather than every Newton solver.

2.2 BHHH: estimate curvature from individual scores

BHHH replaces the observed Hessian by an outer product of individual likelihood scores. For a mean negative log-likelihood, its curvature matrix is $C=n^{-1}\sum_i s_i s_i'$. The information identity connects this matrix to expected likelihood curvature at the correctly specified truth. They need not match in a finite sample or far from the optimum.

This avoids differentiating the score, at the cost of forming a score cross-product at each step. The reference BHHH implementation uses the same eigenvalue floor and line search as safeguarded Newton. Gaussian OLS and NLS require the variance rescaling implied by their squared-error normalization; ridge adds the penalty’s curvature separately. BHHH is used for these likelihood-based problems, logit and Poisson. Fixed-weight IV/GMM receives no BHHH label.

2.3 DFP and BFGS: learn curvature from successive gradients

Quasi-Newton methods build an inverse-curvature approximation $M$ and use $d=-Mg$. After moving by $s$, they observe the gradient change $y$ and update $M$ to satisfy the inverse secant equation $M^+y=s$.

DFP updates the inverse approximation directly:

$$M^+=M+\frac{ss'}{s'y}-\frac{Myy'M}{y'My}.$$

BFGS uses a different update with the same secant target:

$$M^+=(I-\rho sy')M(I-\rho ys')+\rho ss',\qquad \rho=(s'y)^{-1}.$$

Both retain positive definiteness when the curvature condition $s'y>0$ holds. The tested implementations start at $M=I$, skip numerically degenerate updates, and reset a non-descent direction. Their line search is consequential: changing its curvature requirement materially changes DFP’s success rate in this experiment.

2.4 L-BFGS: keep a short curvature history

L-BFGS uses a limited history of step/gradient-change pairs to apply an approximate inverse Hessian without storing a dense matrix. The tested implementations retain seven pairs. This is attractive when parameters are numerous or Hessians are expensive. Here it also provides a useful first-derivative-only comparator at modest dimension. ensmallen and CrabbyMetrics supply separate implementations with their own line-search and stopping rules.

2.5 TRU: minimize the quadratic inside a trust region

GALAHAD TRU solves a local step problem of the form

$$\min_{\|s\|\le\Delta}\; g's+\tfrac12s'Hs.$$

The radius $\Delta$ limits the region over which the model is trusted. The algorithm compares the actual reduction in $f$ with the model’s predicted reduction, accepting useful steps and adjusting the radius. The constraint also makes the step problem well posed when the Hessian is indefinite.

At a global solution of the Euclidean subproblem, the step obeys $(H+\lambda I)s=-g$ with $H+\lambda I\succeq0$. Thus the method can change the direction through a curvature shift, as well as limiting its length. The benchmark uses exact dense Hessians and direct subproblem solves.

2.6 ARC: penalize long steps with a cubic term

GALAHAD ARC solves a regularized model of the form

$$\min_s\;g's+\tfrac12s'Hs+\tfrac{\sigma}{3}\|s\|^3.$$

The cubic term discourages large steps and dominates negative quadratic curvature at large distances. Agreement between the model and the evaluated criterion determines how strongly subsequent steps are regularized. Its step also satisfies a shifted system, with $\lambda=\sigma\|s\|$ in Euclidean coordinates. TRU and ARC therefore share an adaptive curvature-shifting mechanism, while controlling the step through a radius and a penalty, respectively.

2.7 Step acceptance in the reference methods

Reference Newton, BHHH, DFP and BFGS share a strong-Wolfe line search. It asks for sufficient decrease in the objective and a reduction in the directional slope. The experiment compares the usual curvature parameter $c_2=.9$ with the stricter $c_2=.1$. The latter requires a slope closer to zero along the search direction; it can require more trial evaluations and does not simply mean taking a shorter step. Package solvers retain their own controls and are timed once per repetition, without artificial duplication across the two reference settings.

2.8 The constrained synthetic-control methods

Synthetic control minimizes a quadratic over nonnegative donor weights that sum to one. The comparison therefore uses three constrained methods:

• GALAHAD QPB applies a primal–dual interior-point method to the quadratic program.
• CrabbyMetrics’ native solver uses an active-set method, identifying weights on the boundary and solving for the remaining weights.
• ensmallen Frank–Wolfe moves toward a simplex vertex chosen from the current gradient, using convex combinations to preserve feasibility.

These methods exploit the actual constraint set. The unconstrained smooth solvers are outside this track. CrabbyMetrics’ other native estimators likewise use their existing fitting routines and initializers; their timing covers a complete estimator API rather than isolating a generic optimizer.

4.1 Reliability depends on the implementation and the problem

TRU, ARC and reference BFGS pass all 72 smooth instances. The remaining methods show why a single success-rate number requires context: safeguarded Newton’s misses are concentrated in ill-conditioned OLS, while BHHH has difficulty with distant NLS starts. Stronger Wolfe conditions help DFP substantially.

Table 2. Passed instances / attempted instances. BHHH’s denominator is the 60 likelihood-based instances; the other rows cover 72 smooth instances. Native estimator APIs are reported separately in the supplement.

BHHH passes every logit and Poisson instance under both settings, but only four of twelve NLS instances. In the latter model, its score cross-product contains squared residuals as well as derivatives of the exponential mean. Far from a good fit, this curvature approximation can differ greatly from the observed Hessian. Its information-identity motivation provides no uniform guarantee over these starting points.

For DFP, stricter Wolfe search raises acceptance from 54 to 68 of 72 instances. NLS improves from four to twelve passes; IV/GMM improves from ten to twelve. The reference BFGS implementation passes every instance under either setting, while the existing CrabbyMetrics BFGS call passes 49. The label “BFGS” leaves important choices about step acceptance, initialization and termination open. These are comparisons of the stated implementations and tolerances.

4.2 Speed is meaningful only after the accuracy check

Figure 1 compares each solver with the fastest accepted solver on the same instance. Failed fits receive an infinite time ratio, so an early failure cannot appear as a speed advantage. The height at ratio one is the fraction of wins; the eventual height is the fraction of all instances passed.

With reference Wolfe .9, safeguarded Newton is fastest on 34 of the 72 smooth instances and ensmallen L-BFGS on 25. With Wolfe .1, those counts are 41 and 19. TRU and ARC provide broader reliability than either of these configurations, with extra subproblem work. The dimensions here are small enough that forming and factoring a Hessian can repay its cost through fewer iterations; the same balance can change at much larger parameter dimension.

4.3 Scaling and distant starts expose different weaknesses

In Poisson, changing regressor units leaves the statistical problem unchanged, yet ensmallen L-BFGS passes only eight of twelve instances and each generic CrabbyMetrics method passes four. The reference methods, TRU and ARC pass all twelve. This is a coordinate-scaling result rather than a deterioration in the data’s information about the mean.

NLS instead holds the observations fixed and changes initialization. TRU, ARC, safeguarded Newton and reference BFGS pass all twelve instances. ensmallen L-BFGS passes eleven; the two generic CrabbyMetrics methods pass four each. Figure 2 shows a difficult starting point alongside the badly scaled Poisson case. The full timing panels below retain DFP, both Wolfe settings and all package paths, including failures.

4.4 Timing and convergence by estimator

The following panels report every estimator family’s timing and convergence paths. Each uses the larger, stress-2, seed-0 instance: 20,000 observations and 24 coefficients for the smooth problems, or 500 periods and 50 donors for synthetic control. The 84-case atlas contains every size, stress and seed.

Within each figure, the left column measures joint objective/gradient evaluations and the right column measures elapsed time during the instrumented diagnostic warm-up. Smooth-problem rows use reference Wolfe .9 and .1; the package runs are shared across rows. Timing tables and Figure 1 use the separate, uninstrumented repetitions. Lines include trial points and trace the best evaluated objective. Circles rescore the returned vector; open markers fail the common accuracy test. QPB and native estimator APIs expose no callback path, so their time panels show returned endpoints only.

4.4.1 OLS

Collinearity creates flat directions; a small objective error can coexist with slow coefficient adjustment.

4.4.2 Ridge

The slope penalty improves curvature on the same samples used for OLS.

4.4.3 Logit

The strong-signal design weakens likelihood curvature while retaining a finite MLE.

4.4.4 Poisson

The highest-stress case changes regressor units while preserving the true mean function.

4.4.5 IV/GMM

The weak-first-stage case combines a quadratic sample criterion with limited statistical identification.

4.4.6 Nonlinear least squares

The distant initial coefficients test step control in the exponential-mean criterion.

4.4.7 Synthetic control

All three methods enforce the simplex; the Frank–Wolfe gap supplies the convergence certificate.

4.5 Solving the criterion does not remove identification problems

OLS and ridge clarify the distinction between numerical and statistical accuracy. All twelve ridge instances pass for the reference methods and both GALAHAD methods. Without the ridge penalty, collinearity makes a small gradient compatible with appreciable coefficient disagreement. Among accepted OLS fits, reference BFGS differs from the SVD coefficients by as much as 13.8% in relative Euclidean norm, and ARC by 8.8%. Their largest relative prediction disagreements are only about $1.2\times10^{-4}$ and $9.3\times10^{-5}$, respectively. The acceptance target controls objectives and gradients; applications requiring particular coefficients may need a stronger target.

Sampling error is a separate issue. Across the four IV instances at each strength level, the sample-reference estimate’s median raw coefficient RMSE rises from 0.010 to 0.213 as instruments weaken. An optimizer can closely reproduce this sample estimate while remaining far from the generating coefficient. OLS collinearity produces an even sharper coefficient/prediction contrast. Figure 10 separates recovery against truth from numerical acceptance. The full seven-family recovery plots are in the supplement.

4.6 Synthetic control rewards specialized constrained solvers

Both QPB and CrabbyMetrics’ native active-set solver satisfy the tight simplex certificate in all twelve instances. The native solver is faster in these small dense problems. Classic Frank–Wolfe misses the required gap within its fixed budget, despite preserving feasibility. A method’s simplicity and low per-step cost therefore need to be assessed alongside the accuracy it actually reaches.

Table 3. Median time across the twelve instances, after taking the median of three repetitions within each instance. Failed attempts are included in the time column; Frank–Wolfe’s time is therefore a cost to reach its stopping point, without achieving the required solution accuracy.

As donor series become more alike, the reference weight RMSE rises while its prediction RMSE stays small. Several weight vectors can describe nearly the same fitted path. This is the constrained analogue of the flat-direction problem in OLS, and is visible in the full recovery diagnostics.