Four Optimizers for a Logistic Mixture

BHHH and modern likelihood computation

Robert E. “Bob” Hall shared the 2026 Nobel Memorial Prize in Economic Sciences with Michael Woodford “for microeconomic foundations of macroeconomic stabilisation.” The award provides an occasion to revisit an earlier computational contribution: Estimation and Inference in Nonlinear Structural Models (1974), by Ernst Berndt, Bronwyn Hall, Robert Hall, and Jerry Hausman. The initials give the BHHH algorithm its name. Their treatment of likelihood optimization fits naturally with today’s automatic differentiation and vectorized computation.

Economic models often characterize choices through first-order conditions. Estimating their parameters introduces a related computational problem: for an identified model with a smooth criterion $f$ and an interior optimum, find $\widehat\theta$ satisfying

$$g(\widehat\theta)=\nabla f(\widehat\theta)=0.$$

Identification comes from the model and the distribution of observables; curvature and the objective value help distinguish among stationary points. Computational convenience has long favored familiar criteria such as least squares and standard logit, whose gradients and Hessians are readily available. Newton–Raphson solves the first-order condition by linearizing the gradient: $g(\theta+d)\approx g(\theta)+H(\theta)d$, giving $d=-H(\theta)^{-1}g(\theta)$. Richer heterogeneity makes hand-derived curvature more involved. The mixture below, for example, adds terms from changing posterior probabilities of latent types.

BHHH exploits the statistical structure of maximum likelihood. Let $\ell_i(\theta)$ be observation $i$’s log likelihood and $s_i(\theta)=\nabla\ell_i(\theta)$ its score. Under correct specification and regularity, the information-matrix equality at the true parameter $\theta_0$ is

$$\mathbb E_{\theta_0}[s_i(\theta_0)s_i(\theta_0)^\top] =-\mathbb E_{\theta_0}[\nabla^2\ell_i(\theta_0)].$$

For $f=-n^{-1}\sum_i\ell_i$, the sample matrix $B=n^{-1}\sum_i s_i s_i^\top$ therefore estimates the curvature needed near a consistent likelihood solution. BHHH uses $d=-B^{-1}g$ with a line search. Its curvature calculation requires only individual first derivatives. The outer-product matrix is positive semidefinite by construction and supplies a descent direction when positive definite. Finite-sample differences from the Hessian and behavior away from the solution make the line search important.

Autodiff makes this recipe especially convenient: write an observation’s log likelihood, differentiate it, vectorize over observations, and accumulate the score outer products. This lowers the derivative-bookkeeping cost of fitting flexible economic models. The worked example uses a two-type logistic mixture with one million observations. BHHH works particularly well here; Newton, DFP, and BFGS provide comparisons on the same likelihood and starting values.

Model and simulated sample

Newton, BHHH, DFP, and BFGS differ in where they obtain curvature. Newton evaluates the Hessian, BHHH constructs an outer product of individual scores, and DFP and BFGS learn curvature from changes in the gradient. A two-type logistic mixture illustrates these differences through its smooth, nonconvex likelihood.

Let $Z_i\in\{1,2\}$ be an unobserved type, independently drawn for each observation, with $\Pr(Z_i=1)=\pi=\sigma(a)$ and $\sigma(t)=(1+e^{-t})^{-1}$. Conditional on the type,

$$Y_i\mid X_i,Z_i=j\sim\operatorname{Bernoulli}(p_{ij}), \qquad p_{ij}=\sigma(X_i^\top\beta_j).$$

Each observation contributes one binary outcome. With five regressors including the intercept, the parameter vector $\theta=(\beta_1^\top,\beta_2^\top,a)^\top$ has eleven elements. The following code creates a million observations; the latent labels are used only to simulate outcomes.

import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import torch
from torch.nn import functional as F
from IPython.display import display, Markdown

torch.set_default_dtype(torch.float64)
torch.set_num_threads(4)

rng = torch.Generator().manual_seed(20261012)
n = 1_000_000
X = torch.cat((torch.ones(n, 1), torch.randn(n, 4, generator=rng)), dim=1)
X[:, 2] = 0.5 * X[:, 1] + np.sqrt(0.75) * X[:, 2]
beta1 = torch.tensor([-1.2, 1.6, -1.0, 0.8, 0.4])
beta2 = torch.tensor([1.1, -0.6, 1.8, 0.5, -1.2])
pi = 0.4
truth = torch.cat((beta1, beta2, torch.tensor([np.log(pi / (1 - pi))])))
type1 = torch.rand(n, generator=rng) < pi
p = torch.where(type1, (X @ beta1).sigmoid(), (X @ beta2).sigmoid())
y = (torch.rand(n, generator=rng) < p).double()
del type1, p

The marginal success probability is $q_i=\pi p_{i1}+(1-\pi)p_{i2}$. We minimize the mean negative log likelihood,

$$f(\theta)=-\frac1n\sum_i\log\left[ \pi p_{i1}^{Y_i}(1-p_{i1})^{1-Y_i} +(1-\pi)p_{i2}^{Y_i}(1-p_{i2})^{1-Y_i}\right].$$

The implementation uses log-sum-exp for the mixture and autodiff for the gradient. losses also accepts one observation, which will be useful for BHHH. Every derivative and curvature matrix below uses the same $1/n$ normalization.

def losses(theta, X, y):
    eta = X @ theta[:10].reshape(2, 5).T
    log_weights = torch.stack((F.logsigmoid(theta[-1]),
                               F.logsigmoid(-theta[-1])))
    log_components = y.unsqueeze(-1) * eta - F.softplus(eta)
    return -torch.logsumexp(log_components + log_weights, dim=-1)

def fg(theta):
    theta = theta.detach().requires_grad_(True)
    f = losses(theta, X, y).mean()
    g, = torch.autograd.grad(f, theta)
    return f.item(), g.detach()

Swapping the two coefficient vectors and replacing $a$ by $-a$ leaves the likelihood unchanged. Estimates must therefore be label-aligned before comparing parameters. Coincident components and nearly zero mixing shares can produce degeneracy.

Newton: the observed Hessian

Write $c_{ij}=Y_i X_i^\top\beta_j-\log(1+e^{X_i^\top\beta_j})$. The posterior probability of type 1 is

$$r_i=\frac{\pi e^{c_{i1}}}{\pi e^{c_{i1}}+(1-\pi)e^{c_{i2}}}.$$

The complete-data scores are

$$u_{i1}=\begin{pmatrix}(Y_i-p_{i1})X_i\\0\\1-\pi\end{pmatrix}, \qquad u_{i2}=\begin{pmatrix}0\\(Y_i-p_{i2})X_i\\-\pi\end{pmatrix}.$$

Thus the observed score is $s_i=r_i u_{i1}+(1-r_i)u_{i2}$. With $v_i=u_{i1}-u_{i2}$, the Hessian of the mean negative log likelihood is

$$H=\frac1n\sum_i\left[ \operatorname{blockdiag}\left( r_i p_{i1}(1-p_{i1})X_iX_i^\top, (1-r_i)p_{i2}(1-p_{i2})X_iX_i^\top, \pi(1-\pi)\right) -r_i(1-r_i)v_iv_i^\top\right].$$

This is complete information minus missing information. The subtraction accounts for the fact that changing $\theta$ also changes the responsibilities. This term can make $H$ indefinite away from a minimum.

The code accumulates small matrix products in chunks to limit temporary memory. It computes the exact analytic Hessian, before any safeguard used by the solver.

@torch.no_grad()
def hessian(theta, X, y):
    H = torch.zeros(11, 11)
    pi = theta[-1].sigmoid()
    log_weights = torch.stack((F.logsigmoid(theta[-1]),
                               F.logsigmoid(-theta[-1])))
    for A, b in zip(X.split(100_000), y.split(100_000)):
        eta = A @ theta[:10].reshape(2, 5).T
        p = eta.sigmoid()
        r = (b[:, None] * eta - F.softplus(eta) + log_weights).softmax(1)[:, 0]
        w1 = r * p[:, 0] * (1 - p[:, 0])
        w2 = (1 - r) * p[:, 1] * (1 - p[:, 1])
        H[:5, :5] += A.T @ (w1[:, None] * A)
        H[5:10, 5:10] += A.T @ (w2[:, None] * A)
        H[-1, -1] += len(b) * pi * (1 - pi)
        v = torch.cat(((b - p[:, 0])[:, None] * A,
                       -(b - p[:, 1])[:, None] * A,
                       torch.ones(len(b), 1)), dim=1)
        H -= v.T @ ((r * (1 - r))[:, None] * v)
    return (H + H.T) / (2 * len(y))

Newton’s direction is $d=-H^{-1}g$. To ensure descent in nonconvex regions, the solver clips eigenvalues of $H$ below $10^{-5}$ before solving. This is safeguarded Newton, becoming ordinary Newton when the floor is inactive. Second-order autodiff provides an alternative way to compute the same Hessian.

BHHH: curvature from individual scores

BHHH replaces $H$ by $B=n^{-1}\sum_i s_is_i^\top$ and takes $d=-B^{-1}g$. Vectorized autodiff supplies the observation-level scores used in this sum; their sign is irrelevant in the outer product.

individual_grad = torch.func.vmap(torch.func.grad(losses), in_dims=(None, 0, 0))

@torch.no_grad()
def opg(theta, X, y):
    B = torch.zeros(11, 11)
    for A, b in zip(X.split(100_000), y.split(100_000)):
        G = individual_grad(theta, A, b)
        B += G.T @ G
    return B / len(y)

Under correct specification and regularity, the expected OPG and negative log-likelihood Hessian agree at the true parameter. Their realized sample values generally differ. BHHH uses first derivatives to estimate curvature; its local convergence rate depends on the accuracy of that approximation.

DFP and BFGS: curvature from successive gradients

Both methods take $d=-Mg$, starting from $M=I$. After an accepted step, define $s=\theta_+-\theta$ and $z=g_+-g$. The secant condition is $M_+z=s$. The inverse-Hessian updates are

$$M_+^{\mathrm{DFP}}=M+\frac{ss^\top}{s^\top z} -\frac{(Mz)(Mz)^\top}{z^\top Mz},$$

$$M_+^{\mathrm{BFGS}}=(I-\rho sz^\top)M(I-\rho zs^\top)+\rho ss^\top, \qquad \rho=(s^\top z)^{-1}.$$

Both preserve positive definiteness when $s^\top z>0$. These updates depend on the change in the gradient along the optimization path. BHHH uses variation in scores across observations.

def update(M, s, z, method):
    if method == "DFP":
        Mz = M @ z
        out = M + torch.outer(s, s) / (s @ z) - torch.outer(Mz, Mz) / (z @ Mz)
    else:  # BFGS
        rho = 1 / (s @ z)
        V = torch.eye(11) - rho * torch.outer(s, z)
        out = V @ M @ V.T + rho * torch.outer(s, s)
    return (out + out.T) / 2

A common line search and stopping rule

Every method uses the same strong-Wolfe conditions,

$$f(\theta+\alpha d)\le f(\theta)+10^{-4}\alpha g^\top d, \qquad |\nabla f(\theta+\alpha d)^\top d|\le c_2|g^\top d|.$$

The search first doubles the trial step to find a bracket, then bisects it. It returns only a step satisfying both conditions; forty unsuccessful trials cause an explicit failure. Smaller $c_2$ demands a slope closer to zero along the search direction.

def wolfe(theta, f, g, d, c2):
    slope0 = float(g @ d)
    assert slope0 < 0
    calls = 0

    def evaluate(alpha):
        nonlocal calls
        calls += 1
        ft, gt = fg(theta + alpha * d)
        return ft, gt, float(gt @ d)

    def zoom(lo, hi, flo):
        while calls < 40:
            alpha = (lo + hi) / 2
            ft, gt, slope = evaluate(alpha)
            if not np.isfinite(ft + slope) or ft > f + 1e-4 * alpha * slope0 or ft >= flo:
                hi = alpha
            else:
                if abs(slope) <= -c2 * slope0:
                    return alpha, ft, gt, calls
                if slope * (hi - lo) >= 0:
                    hi = lo
                lo, flo = alpha, ft
            if abs(hi - lo) < 1e-16:
                break
        raise RuntimeError("line search failed")

    previous, previous_f, alpha = 0., f, 1.
    while calls < 40:
        ft, gt, slope = evaluate(alpha)
        if (not np.isfinite(ft + slope) or ft > f + 1e-4 * alpha * slope0
                or (calls > 1 and ft >= previous_f)):
            return zoom(previous, alpha, previous_f)
        if abs(slope) <= -c2 * slope0:
            return alpha, ft, gt, calls
        if slope >= 0:
            return zoom(alpha, previous, ft)
        previous, previous_f, alpha = alpha, ft, 2 * alpha
    raise RuntimeError("line search failed")

The following loop is the complete optimizer. It stops when $\|g\|_\infty\le10^{-7}$ or after 2,000 steps. The same eigenvalue floor is available for Newton and BHHH. Numerically degenerate secant updates are skipped; a non-descent direction is reset to steepest descent. A trace retains elapsed time, objective, and gradient norm at every accepted iterate.

def fit(method, initial, c2):
    theta, M = initial.clone(), torch.eye(11)
    start = time.perf_counter()
    f, g = fg(theta)
    evaluations, clipped, skipped, resets = 1, 0, 0, 0
    trace = [(time.perf_counter() - start, f, float(g.abs().max()))]
    status = "iteration limit"
    for _ in range(2000):
        if g.abs().max() <= 1e-7:
            break
        if method in ("Newton", "BHHH"):
            C = hessian(theta, X, y) if method == "Newton" else opg(theta, X, y)
            eigenvalues, Q = torch.linalg.eigh(C)
            clipped += int(eigenvalues.min() < 1e-5)
            d = -Q @ ((Q.T @ g) / eigenvalues.clamp_min(1e-5))
        else:
            d = -M @ g
        if not torch.isfinite(d).all() or g @ d >= 0:
            M, d, resets = torch.eye(11), -g, resets + 1
        try:
            alpha, fn, gn, calls = wolfe(theta, f, g, d, c2)
        except RuntimeError:
            status = "line search failed"
            break
        evaluations += calls
        assert fn <= f + 1e-4 * alpha * float(g @ d)
        assert abs(float(gn @ d)) <= c2 * abs(float(g @ d))
        next_theta = theta + alpha * d
        s, z = next_theta - theta, gn - g
        if method in ("DFP", "BFGS"):
            if s @ z > 1e-12 * torch.linalg.vector_norm(s) * torch.linalg.vector_norm(z):
                M = update(M, s, z, method)
            else:
                skipped += 1
        theta, f, g = next_theta, fn, gn
        trace.append((time.perf_counter() - start, f, float(g.abs().max())))
    if g.abs().max() <= 1e-7:
        status = "converged"
    return dict(method=method, c2=c2, status=status, theta=theta.tolist(),
                iterations=len(trace) - 1, seconds=time.perf_counter() - start,
                nll=f, gradient=float(g.abs().max()), evaluations=evaluations,
                clipped=clipped, skipped=skipped, resets=resets, trace=trace)

Derivative checks

A small subset is sufficient to compare the analytic Hessian with second-order autodiff. We also check the individual gradients against direct differentiation of $q_i=\pi p_{i1}+(1-\pi)p_{i2}$, and check the two secant equations.

rng = torch.Generator().manual_seed(104729)
starts = [0.7 * torch.randn(11, generator=rng) for _ in range(3)]
A, b = X[:91], y[:91]
errors = []
for t in [truth, *starts]:
    H_auto = torch.func.hessian(lambda v: losses(v, A, b).mean())(t)
    errors.append(float((hessian(t, A, b) - H_auto).abs().max()))
    torch.testing.assert_close(hessian(t, A, b), H_auto, atol=1e-12, rtol=1e-10)
    p1, p2 = (A @ t[:10].reshape(2, 5).T).sigmoid().T
    w = t[-1].sigmoid()
    q = w * p1 + (1 - w) * p2
    dq = torch.cat((w * (p1 * (1 - p1))[:, None] * A,
                    (1 - w) * (p2 * (1 - p2))[:, None] * A,
                    (w * (1 - w) * (p1 - p2))[:, None]), dim=1)
    G = ((q - b) / (q * (1 - q)))[:, None] * dq
    torch.testing.assert_close(individual_grad(t, A, b), G)
    torch.testing.assert_close(opg(t, A, b), G.T @ G / len(b))
    for method in ("DFP", "BFGS"):
        s = t
        z = (G.T @ G + torch.eye(11)) @ s
        M = update(torch.eye(11), s, z, method)
        torch.testing.assert_close(M @ z, s)
        assert torch.linalg.eigvalsh(M).min() > 0
print(f"Largest Hessian discrepancy: {max(errors):.2e}; score and secant checks passed.")

Largest Hessian discrepancy: 1.39e-16; score and secant checks passed.

Million-observation comparison

Each method receives the same three starts, drawn independently of the true parameters. The main comparison uses $c_2=0.9$. A second comparison uses $c_2=0.1$ to examine DFP’s sensitivity to line-search accuracy; this sensitivity comparison was motivated by its slow convergence with $c_2=0.9$.

These measurements use the Apple M4 CPU, four Torch threads, and float64. Each time covers one complete fit, including trial gradient evaluations, curvature calculations, and line search. Data generation and the following warm-up precede timing. Each initial condition is timed once. The three starts are indexed 0, 1, and 2; method order rotates across them.

methods = ["Newton", "BHHH", "DFP", "BFGS"]
fg(starts[0])
hessian(starts[0], A, b)
opg(starts[0], A, b)
runs = []
for c2 in (0.9, 0.1):
    for k, initial in enumerate(starts):
        for method in methods[k:] + methods[:k]:
            runs.append(dict(start=k, **fit(method, initial, c2)))
results = pd.DataFrame(runs)

summary = results.groupby(["c2", "method"], sort=False).agg(
    converged=("status", lambda v: f"{sum(v == 'converged')}/{len(v)}"),
    median_seconds=("seconds", "median"),
    iterations=("iterations", lambda v: f"{v.min()}–{v.max()}"),
    gradient_evaluations=("evaluations", lambda v: f"{v.min()}–{v.max()}"))
display(summary.round(3))

Medians include all three starts, including time spent on iteration-limited fits. The per-fit table reports each stopping status and final gradient norm.

display(results[["c2", "method", "start", "status", "iterations", "seconds", "gradient"]]
        .style.hide(axis="index").format({"c2": "{:.1f}", "seconds": "{:.3f}",
                                          "gradient": "{:.2e}"}))

best_value = results.nll.min()
fig, axes = plt.subplots(2, 2, figsize=(11, 7), constrained_layout=True)
colors = dict(zip(methods, ["#2765ab", "#167d58", "#b94435", "#8a54a2"]))
for row, c2 in enumerate((0.9, 0.1)):
    for run in runs:
        if run["start"] != 0 or run["c2"] != c2:
            continue
        trace = np.array(run["trace"])
        gap = np.maximum(trace[:, 1] - best_value, 1e-14)
        for col, x in enumerate((np.arange(len(trace)), trace[:, 0])):
            axes[row, col].semilogy(x, gap, color=colors[run["method"]], label=run["method"])
    for col, ax in enumerate(axes[row]):
        ax.set(xlabel="Accepted iteration" if col == 0 else "Elapsed seconds",
               ylabel="Mean NLL − best found", title=f"Strong Wolfe: c₂ = {c2}")
        ax.grid(alpha=.2)
    axes[row, 0].legend()
plt.show()

Recovery at the fitted solution

Recovery errors coincide across the converged fits. The figure reports coefficient RMSE over the ten label-aligned coefficients and probability RMSE over the million observed covariate vectors, using the generating values as truth. These are single-sample checks; the three starts assess optimization sensitivity. The iteration-limited DFP fit remains visible.

estimates = torch.tensor(results.theta.tolist())
swapped = torch.cat((estimates[:, 5:10], estimates[:, :5], -estimates[:, -1:]), dim=1)
swap = ((swapped[:, :10] - truth[:10]).square().sum(1)
        < (estimates[:, :10] - truth[:10]).square().sum(1))
estimates = torch.where(swap[:, None], swapped, estimates)
results["beta_rmse"] = (estimates[:, :10] - truth[:10]).square().mean(1).sqrt().numpy()
results["pi_hat"] = estimates[:, -1].sigmoid().numpy()
results["pi_abs_error"] = abs(results.pi_hat - pi)

def probability(theta):
    w = theta[-1].sigmoid()
    return w * (X @ theta[:5]).sigmoid() + (1 - w) * (X @ theta[5:10]).sigmoid()

q_true = probability(truth)
results["q_rmse"] = [float((probability(t) - q_true).square().mean().sqrt())
                     for t in estimates]

Plot recovery across optimizers

fig, axes = plt.subplots(1, 2, figsize=(10, 3.4), layout="constrained")
for ax, metric, title in zip(axes, ("beta_rmse", "q_rmse"),
                           ("Coefficient RMSE", "Probability RMSE")):
    for c2, offset, marker, color in ((0.9, -.15, "o", "#2765ab"),
                                     (0.1, .15, "^", "#b65b16")):
        subset = results[results.c2 == c2]
        x = subset.method.map(dict(zip(methods, range(4)))) + offset + .07 * (subset.start - 1)
        faces = [color if status == "converged" else "none" for status in subset.status]
        ax.scatter(x, subset[metric], marker=marker, facecolors=faces,
                   edgecolors=color, s=42, label=f"c₂ = {c2}", zorder=3)
        for i in subset.index[subset.status != "converged"]:
            ax.annotate("Iteration limit", (x.loc[i], subset.loc[i, metric]),
                        xytext=(0, -15), textcoords="offset points", ha="center", fontsize=8)
    ax.set(xticks=range(4), xticklabels=methods, yscale="log", title=title)
    ax.margins(x=.13, y=.25)
    ax.grid(axis="y", alpha=.2)
    ax.legend(loc="upper left", frameon=False, fontsize=9)
plt.show()

Summarize the common recovery error

converged = results[results.status == "converged"]
display(Markdown(
    f"The **{len(converged)} converged fits** have coefficient RMSE about "
    f"**{converged.beta_rmse.median():.5f}** and probability RMSE about "
    f"**{converged.q_rmse.median():.6f}**. The main distinction is the work required "
    "to reach this common solution."))

The 23 converged fits have coefficient RMSE about 0.00854 and probability RMSE about 0.000953. The main distinction is the work required to reach this common solution.

Local curvature at the Newton fit

For a concrete curvature diagnostic, take safeguarded Newton with $c_2=0.9$, start 0. The following calculation evaluates both curvature matrices at that fit’s label-aligned estimate.

reference = results[(results.method == "Newton") & (results.c2 == 0.9)
                    & (results.start == 0)].iloc[0]
theta_hat = estimates[reference.name]
H, B = hessian(theta_hat, X, y), opg(theta_hat, X, y)
eigenvalues = torch.linalg.eigvalsh(H)

display(Markdown(
    f"For **Newton, c₂ = 0.9, start 0** ({reference.status}), "
    f"the smallest Hessian eigenvalue is **{eigenvalues[0]:.3g}**, "
    f"its condition number is **{eigenvalues[-1] / eigenvalues[0]:.1f}**, "
    f"and the relative Frobenius difference between OPG and Hessian is "
    f"**{100 * torch.linalg.matrix_norm(B - H) / torch.linalg.matrix_norm(H):.3f}%**."))

A positive-definite Hessian at a stationary point establishes a strict local minimum. The small OPG–Hessian discrepancy helps explain BHHH’s rapid local convergence on this sample. Agreement across starts gives evidence of numerical reproducibility; global optimality would require further analysis.

Reproduce the document

Source and replication downloads below are hosted in the private Tailscale library.

Download the Quarto document or the source bundle. The numerical methods appear beside their derivations; setup, checks, formatting, and exports are in expandable code blocks. Rendering this single document regenerates the sample, runs all 24 fits, checks derivatives, and creates every table and figure using only the code above and installed packages. It requires Python 3.13, Quarto, and the following packages:

uv venv --python 3.13
uv pip install torch==2.9.1 numpy==2.5.2 pandas==3.0.5 matplotlib==3.11.1 ipykernel
QUARTO_PYTHON=.venv/bin/python quarto render notes.qmd

The timing summary, per-fit recovery metrics, and full iteration histories are also saved as timing CSV, recovery CSV, and JSON:

recovery = results[["c2", "method", "start", "status", "beta_rmse", "pi_hat",
                    "pi_abs_error", "q_rmse"]]
summary.to_csv("summary.csv")
recovery.to_csv("recovery.csv", index=False)
results.to_json("fits.json", orient="records", indent=2, double_precision=15)