Every time you call
.fit() on a statsmodels GLM, fit a survival model, calibrate a mixture model, or write your own custom likelihood in Python, there's an optimizer quietly grinding away underneath — usually scipy.optimize.minimize (we minimize the negative log-likelihood, since there's no maximize). Most of the time nobody thinks about which algorithm it's actually running. It defaults to BFGS, the fit converges, the numbers look plausible, and everyone moves on.The problem is that "converged" and "found the maximum likelihood estimate" are not the same statement. An optimizer converges when it stops moving — not necessarily when it has found the best point in the whole parameter space. For convex, well-behaved likelihoods (most GLMs, most simple exponential-family models) this distinction never bites you. For anything with latent structure — mixture models, hidden Markov models, random-effects models, some survival and copula models — the likelihood surface can be genuinely non-convex, riddled with saddle points and inferior local optima, and the algorithm you pick (and where you start it) determines whether you get the real answer or a plausible-looking wrong one.
This is worth being deliberate about, in the same spirit as Wilson confidence intervals and Cohen's h — the tool you default to isn't automatically the right one for the shape of the problem in front of you.
The Cast of Characters
scipy.optimize.minimize(fun, x0, method=...) supports a long list of algorithms. Five come up constantly in MLE work, and they split into two families.Gradient-based (need the surface to be smooth):
- BFGS — a quasi-Newton method. It doesn't compute the true Hessian (the matrix of second derivatives); it builds up an approximation to the inverse Hessian iteratively from the gradients it observes, then takes Newton-like steps using that approximation. Superlinear convergence near the optimum, and very few function evaluations for smooth, well-scaled problems.
- L-BFGS(-B) — the limited-memory version of BFGS. Instead of storing a full approximate Hessian, it keeps only the last update vectors (typically 5–20) and reconstructs the search direction from those on the fly. For a handful of parameters this barely matters; for a model with hundreds or thousands of parameters (a large hierarchical model, an embedding layer, a big design matrix), storing the full Hessian is simply not possible, and L-BFGS is what makes quasi-Newton optimization feasible at all. The
-Bvariant also supports simple box constraints, which is exactly what you need to keep a variance or scale parameter positive during MLE without reparameterizing.
Derivative-free (don't need — or trust — a gradient):
- Nelder-Mead — the simplex method. It maintains points forming a simplex and moves it downhill by reflecting, expanding, or contracting through the worst point. No derivatives, no line search, nothing but function evaluations. Robust to a noisy or non-smooth likelihood (e.g. one involving a numerical integral or a simulation step where gradients aren't available or aren't reliable), but it scales poorly beyond a handful of parameters and has no real convergence guarantees on non-convex surfaces.
- Powell's method — also derivative-free, but structured differently: it does a sequence of one-dimensional line searches along a set of directions, then updates the direction set based on the net progress made each cycle (a simplified version of conjugate-direction search). Often more direct than Nelder-Mead on smoother derivative-free problems, at the cost of more function evaluations per iteration.
Global (don't trust that any one local search will find the right hill):
- Basin-hopping — not a local optimizer at all, but a wrapper around one. It perturbs the current best point at random, runs a local optimizer (Nelder-Mead or BFGS, your choice) from the perturbed point, and accepts or rejects the new local optimum with a Metropolis criterion — similar in spirit to simulated annealing, except the "state" being annealed is which basin of attraction you're standing in. It is much more expensive (many local optimizations instead of one), but it's the closest thing on this list to a genuine defence against local optima.
None of these algorithms know or care that your objective function is a negative log-likelihood. They are general-purpose numerical optimizers. What makes them behave differently on MLE problems specifically is that likelihood surfaces for latent-variable models (mixtures, HMMs, random effects) are frequently multimodal, with symmetric "label-switching" optima and degenerate saddle points — properties that trip up local optimizers in very characteristic ways.
A Concrete Trap: The Mixture Model Saddle Point
The textbook example is a two-component Gaussian mixture. Say your data actually come from two clusters, and you're maximizing the likelihood over the two component means (holding the mixing weight and variance fixed for simplicity). This surface has:
- Two genuine global optima that are mirror images of each other — and give identical likelihood, because nothing in the model distinguishes "component 1" from "component 2". This is the well-known label-switching symmetry.
- A degenerate saddle point at (both components collapsed onto the overall sample mean) — a stationary point of the likelihood that a purely gradient-based method can get stuck at or near, depending on the starting point, even though it's a dramatically worse fit than the real answer.
Start a local optimizer near the middle of the data and there's a real chance it walks straight into that saddle point and calls it converged, rather than separating into the two clusters that are actually there.
Try It Yourself — Watch Five Optimizers Race
Below is exactly this setup. The heatmap is a negative log-likelihood surface over ; the gold diamond is the true global optimum found by brute-force grid search. Pick a starting point and watch where BFGS, L-BFGS, Nelder-Mead, Powell, and basin-hopping each end up.
Switch the surface to "Easy — two independent means (convex)" — a genuinely convex NLL — and every method lands in the same spot almost immediately, because there's only one hill to descend. It's when you flip to the "Hard — 2-component Gaussian mixture" surface and start near the middle that the algorithms diverge in where they end up: some walk straight into the saddle point, some find one of the two symmetric global optima, and basin-hopping — at the cost of an order of magnitude more function evaluations — reliably escapes the trap by restarting from random perturbations.
What This Actually Means for Picking a Method
If your likelihood is smooth and (you believe) unimodal: use BFGS, or L-BFGS-B if you have many parameters or need box constraints (e.g. keeping a variance ). Fastest convergence, fewest function evaluations, and gradient information — even numerically approximated — genuinely helps.
If your likelihood involves something non-smooth (a hard threshold, a simulation step, a numerical integral with limited precision) so gradients are unreliable or unavailable: use Nelder-Mead or Powell. They're slower and don't scale to high dimensions, but they don't depend on a gradient that might be garbage.
If you have real reason to suspect multiple local optima — mixture models, HMMs, random-effects models, anything with latent structure or label-switching symmetry — a single call to any local optimizer, from a single starting point, is not enough. Either run basin-hopping, or do the cheaper version yourself: multi-start. Fit from a handful of different, sensibly spread-out starting points (e.g. from a k-means initialization for a mixture model) and keep the best result. This is standard practice in the mixture-model and EM literature for exactly this reason.
Notice that this is the same underlying lesson as Cohen's h: the default statistical tool answers a narrower question than the one you actually care about. Wilson intervals tell you "is this detectable," not "is this material." BFGS tells you "I found a stationary point," not "I found the global maximum likelihood estimate." In both cases the fix is the same shape — pair the fast default with a second, cheap check (an effect size; a second starting point) before you trust the output.
Practical Takeaway
scipy.optimize.minimize's default (BFGS, or L-BFGS-B if you set bounds) is the right choice for the large majority of everyday MLE fitting, and there's no need to reach for anything fancier when your likelihood is smooth and well-behaved. But the moment your model has latent structure — mixtures, hidden states, random effects, anything where components could plausibly swap identities — treat a single optimizer run as a starting hypothesis, not an answer. Multi-start it, or wrap it in basin-hopping, and check that the different runs actually agree before you report the fitted parameters to anyone.