Outer-product Hessians

BHHHModel(; ridge = 0.0, dense_max = 5_000)

B = (1/N) Σₙ sₙ sₙᵀ, the average outer product of the per-observation scores.

Two properties, and the trade between them is the whole story:

• Cheap. The scores are already computed for the gradient, so B costs one rank-N product and no second derivatives.
• Positive semidefinite by construction. Every step is a descent direction, even where the true Hessian is indefinite.

The second is a guarantee and a blind spot at once: B ⪰ 0 always, so the model can never report negative curvature, and a run over it converges contentedly to a saddle while every first-order diagnostic looks healthy. Declared through reports_negative_curvature; the solver warns when it meets SecondOrder or tol_H.

Requires a LikelihoodProblem — see the file header. Use information_identity_error to measure how far the identity is from holding on the problem at hand rather than assuming it.

ridge > 0 matters when N < n: B is a sum of N rank-one terms, so a small batch makes it singular regardless of the problem, and under sampling the ridge and the batch size are not independent choices. It is also the right place for an L2 penalty on the problem, whose Hessian is λI — folding λθ into each score instead gives a rank-one λ²θθᵀ.

BHHH2Model(; ridge = 0.0, dense_max = 5_000)

W = (1/N) Σₙ (sₙ − ḡ)(sₙ − ḡ)ᵀ: the sample covariance of the scores rather than their average outer product.

B = W + ḡḡᵀ, so the two coincide exactly at a stationary point and differ by a rank-one term that is largest early in a run — precisely where the identity does not hold and neither matrix approximates the Hessian. Which is better there is an empirical question this package can answer, which is the reason to carry both.

Requires a LikelihoodProblem.

GaussNewtonModel(; ridge = 0.0, dense_max = 5_000)

B = (1/N) JᵀJ for f(x) = (1/2N) Σ rₙ(x)², whose exact Hessian is (1/N)(JᵀJ + Σ rₙ ∇²rₙ).

Gauss-Newton drops the second term. The approximation is good exactly when that term is small — small residuals at the solution, or nearly linear rₙ — and can be arbitrarily bad otherwise, which is the large-residual case where it loses its superlinear rate. gauss_newton_error measures the discarded term beside the RMS residual that controls it.

Requires an NLSProblem: the residual Jacobian is what it needs and nothing above that type has one.

OuterProductOperator(S, ḡ, centred, ridge)

Matrix-free B·v for the outer-product models: Bv = (1/N) S(Sᵀv) + ridge·v, two matrix–vector products against the n × N score matrix and no n × n array anywhere.

This is what makes the models usable at neural-network scale. Forming B densely costs n² — 25 450 parameters is 648 million entries — while S costs n·N and stays in memory for a batch of a few hundred. Pair it with SteihaugCG.

S is AbstractMatrix{T} and eltype is defined. Both matter: the field was Matrix{T}, so an Adjoint (which NLS scores used to be) failed to match the constructor, and Krylov.cg queries eltype, so KrylovCG over a large BHHHModel could not run at all.

information_identity_error(prob, x; batch = full_batch(prob)) -> NamedTuple

How badly V = −H fails at x, in relative Frobenius norm:

(B_err = ‖B − ∇²f‖ / ‖∇²f‖,
 W_err = ‖W − ∇²f‖ / ‖∇²f‖,
 grad_norm = ‖ḡ‖,
 BW_gap = ‖B − W‖ / ‖∇²f‖)     # equals ‖ḡḡᵀ‖/‖∇²f‖

For a correctly specified model at the true parameters these decay like 1/√N; anywhere else they do not decay at all. Belonging to LikelihoodProblem asserts that f is a negative log-likelihood; whether the model is correctly specified is a statement about the data that no type can carry, and this is how to check it.

scores(prob, x, batch) -> Matrix

The n × |batch| matrix whose column i is ∇f_i(x).

The column mean is the gradient, so this is strictly more information than a gradient evaluation and usually costs no more to produce. It is what the outer-product models consume, and the reason a likelihood admits a Hessian approximation a general objective does not.

Implementations must return a genuine Matrix, not a lazy Adjoint: OuterProductOperator stores what it is given.

loss_terms(prob, x, batch) -> Vector

The per-observation values f_i(x), whose mean is f(x). Kept separate from the sum because their sample variance is what the adaptive rules need.

score_matrix(nlp, x) -> Matrix

The n × N matrix of per-observation gradients at x, which the outer-product models consume.

Defined for FullBatchNLP (all terms) and for the two sampled oracles over a ScoredProblem (the current batch — so under sampling, BHHH is formed from the same realisations as the gradient, which is what keeps B and g consistent within an iteration).

LikelihoodNLP

Alias of FullBatchNLP, kept because the type was called that before the problem classes were separated. The new name says what it is — the full-batch, and therefore deterministic, view of a finite sum — rather than what it is usually used for.

residuals(prob, x, batch) -> Vector

The residual vector r(x) of an NLSProblem, restricted to batch.

jacobian(prob, x, batch) -> Matrix

The |batch| × n Jacobian ∂r/∂x. Only first derivatives appear: dropping the Σ rₙ ∇²rₙ term is the whole of the Gauss-Newton approximation.

LogisticRegression(X, y)
LogisticRegression(; K = 5, M = 2_000, β_star = nothing, seed = 0,
                     intercept = true, signal = 1.0)

Binary logistic regression: minimise the average negative log-likelihood

\[f(\beta) = -\frac1M \sum_n \bigl[ y_n \ln \sigma(x_n^\top\beta) + (1-y_n)\ln(1-\sigma(x_n^\top\beta)) \bigr] .\]

The keyword form generates synthetic data with the model correctly specified: xₙ ~ N(0, I) (first coordinate 1 if intercept), and yₙ ~ Bernoulli(σ(xₙᵀβ*)) drawn from the very model being fitted. K sets the number of parameters, so the identity can be watched as the dimension grows.

The generating parameters are passed as β_star and retrieved with β_true; the keyword was previously β_true as well, colliding with the accessor of the same name.

Why the identity holds here, exactly

The score is sₙ = (yₙ − pₙ) xₙ and the per-observation Hessian is pₙ(1−pₙ) xₙxₙᵀ. At β*, E[(y−p)² | x] = Var(y | x) = p(1−p), so

\[V = E\bigl[(y-p)^2 x x^\top\bigr] = E\bigl[p(1-p) x x^\top\bigr] = -H\]

in the population — not approximately, identically. In a sample of size M the two agree to O(M^{-1/2}), which is what information_identity_error measures. Away from β* the first equality fails, because E[(y−p̂)²] = p(1−p) + (p − p̂)², and the gap (p − p̂)² does not shrink with M at all. That extra term is precisely why BHHH takes poor steps far from the optimum and good ones near it.

MultinomialLogit(X, y)
MultinomialLogit(; K = 10, M = 100_000, J = 5, β_star = nothing, seed = 0,
                   signal = 1.0, noise = 0.0)

Multinomial logit over J alternatives with utilities linear in the parameters.

Individual n faces J alternatives; alternative j carries the attribute vector x_{nj} ∈ R^K and the representative utility V_{nj} = β' x_{nj}. Under independent Gumbel errors the choice probabilities are

\[P_{nj} = \frac{e^{\beta^\top x_{nj}}}{\sum_{l} e^{\beta^\top x_{nl}}} ,\]

and the objective is the average negative log-likelihood

\[f(\beta) = -\frac1M \sum_{n=1}^{M} \sum_{j=1}^{J} y_{nj} \ln P_{nj} , \qquad y_{nj} = 1 \text{ iff } n \text{ chose } j .\]

X is M × J × K and y holds the chosen alternative in 1:J. There is no alternative-specific constant: the utility is exactly β' x_{nj}, with one generic coefficient vector shared across alternatives.

This is the genuine multiple-alternative case, not J binary fits. LogisticRegression is the binary problem in its own parameterisation and stays the cheaper choice when J = 2.

Scores and curvature

Writing p_n for the probability vector and x̄_n = Σ_j P_{nj} x_{nj} for the probability-weighted mean attribute, the per-observation score and Hessian are

\[s_n = \bar x_n - x_{n c_n} , \qquad H_n = \sum_j P_{nj} x_{nj} x_{nj}^\top - \bar x_n \bar x_n^\top ,\]

so H_n is the covariance of the attribute vector under the choice probabilities: positive semidefinite, and singular exactly when the attributes carry no variation across alternatives that the current β can see. The objective is therefore convex in β, which the binary case shares.

The information identity

At β* with noise = 0 the data were generated by the very model being fitted, so E[s s'] = -E[H] in the population and BHHH is justified. In a sample of M the two agree to O(M^{-1/2}); information_identity_error measures the gap.

Cost

X holds M·J·K entries: the 100_000 × 5 × 10 default is 40 MB. loss_terms and scores cost O(N·J·K) over a batch of N individuals, batch_hess costs O(N·J·K²).

MLPClassifier(X, y, n_class; hidden = 16, λ = 0.0)

Softmax cross-entropy over a one-hidden-layer tanh network, with per-observation scores so that BHHH applies.

X is M × d, y holds labels in 1:n_class. The parameter vector packs (W₁, b₁, W₂, b₂) in that order, so n = hidden·d + hidden + n_class·hidden + n_class. λ > 0 adds λ‖θ‖²/2.

B = (1/N) Σ (sᵢ + λθ)(sᵢ + λθ)ᵀ,

LeastSquares(t, y, φ, ∇φ!, n; x_star = nothing)

Nonlinear least squares with M observations and n parameters:

• t is M × d, row n being the design point tₙ;
• φ(tₙ, x) -> Real is the model;
• ∇φ!(tₙ, x, g) -> nothing writes its gradient in x into g.

The gradient callback is now in place. Previously it returned a fresh Vector per row, so jacobian allocated M vectors on every call, and jacobian is called once per iteration by GaussNewtonModel and again by every scores evaluation.

The per-observation score is rₙ ∇φₙ and the Jacobian row is ∇φₙᵀ, so GaussNewtonModel, BHHHModel and ExactHessian can all be run against the same problem — which is the comparison worth making, since Gauss–Newton and BHHH discard different terms and their failure modes separate on residual size.

Use linear_least_squares or exponential_fit for the two cases this package exercises.

linear_least_squares(A, b; x_star = nothing)
linear_least_squares(; n = 5, M = 2_000, noise = 0.1, seed = 0)

rₙ(x) = aₙᵀx − bₙ, so ∇²rₙ = 0 and Gauss–Newton is exact: JᵀJ/M = ∇²f identically, with no discarded term at all.

That makes it the control for the least-squares half of the comparison, exactly as a correctly specified logistic regression is the control for the likelihood half. gauss_newton_error should return zero to rounding here whatever the residual size, which is the property that distinguishes "Gauss–Newton is exact" from "Gauss–Newton is accurate because the residuals are small".

The keyword form generates A with independent normal entries and b = A x* + noise·ε.

exponential_fit(; n_terms_model = 2, M = 400, noise = 0.05, misfit = 0.0, seed = 0)

Sum-of-exponentials fitting, the classical nonlinear least-squares testbed:

\[\varphi(t, x) = \sum_{j=1}^{J} x_{2j-1}\, e^{-x_{2j} t}, \qquad n = 2J .\]

misfit > 0 adds a component to the data that the model cannot represent — a term misfit·sin(3t) — which raises the residual at the solution without changing the parameterisation. That knob is the point of this problem: it moves the fit continuously from small-residual to large-residual, and

• misfit = 0: residuals are noise-sized, Σ rₙ∇²rₙ is small, Gauss–Newton is a good approximation and converges superlinearly;
• misfit large: the residuals do not vanish, the discarded term is O(1), and Gauss–Newton degrades to linear convergence while still producing descent directions and healthy ratios.

So it exhibits the same pattern as BHHH under misspecification — a positive semidefinite model whose justification has quietly gone — arrived at from the other direction. gauss_newton_error measures it.

gauss_newton_error(p::LeastSquares, x; batch = full_batch(p)) -> NamedTuple

How much the Gauss–Newton approximation discards at x:

(GN_err   = ‖JᵀJ/N − ∇²f‖ / ‖∇²f‖,
 residual = ‖r‖ / √N,          # the RMS residual, which controls the gap
 rank_def = N < n)             # JᵀJ singular by construction when it is

The counterpart of information_identity_error for the other outer-product model. The discarded term is (1/N)Σ rₙ ∇²rₙ, so GN_err should track residual: small at a well-fitting solution, O(1) when the model cannot fit the data. That is a different failure mode from BHHH's — it depends on the fit, not on the specification — and it is worth measuring rather than assuming, for the same reason.

x_true(prob) -> Vector or nothing

The parameters that generated the data, when the problem was constructed synthetically. As with β_true, the least-squares estimate for a finite sample differs from them by the sampling error, so ‖x̂ − x*‖ does not vanish at the optimum.

init_params(p::MLPClassifier; seed = 0) -> Vector

Random initial weights, scaled by 1/√fan_in so the tanh units start unsaturated. Biases start at zero.

init_params(p::MultinomialLogit) -> Vector

zeros(K), the point at which every alternative is equally likely.

The PreDoc's logit runs start from x₀ = 0, and it is the natural origin here: at β = 0 all P_{nj} = 1/J, so f = log J whatever the data, and the gradient is the difference between the mean attribute and the mean chosen attribute.

accuracy(p::MLPClassifier, θ, batch = 1:p.M) -> Float64

Fraction of batch classified correctly at θ. Reported alongside the loss because a trust-region run reports objective values, and on a classification problem the objective is not the quantity anyone cares about.

accuracy(p::MultinomialLogit, β, batch = 1:p.M) -> Float64

Fraction of batch whose chosen alternative has the largest utility at β.

Reported alongside the loss because the loss is not the quantity a discrete-choice model is judged on. 1/J is the value at β = 0.

β_true(prob) -> Vector or nothing

The parameters that generated the synthetic data, when there are any.

A run should be scored against these rather than against its own stopping test — though note the maximum-likelihood estimate for a finite sample differs from them by the usual O(M^{-1/2}) sampling error, so ‖β̂ − β*‖ does not go to zero at the optimum.

β_true(p::MultinomialLogit) -> Vector or nothing

The generating coefficients. With noise > 0 these are the mean of the per-individual coefficients, not the minimiser of f: the plain-logit maximum likelihood estimate of a mixed-logit population is not β* even as M → ∞.