Module fd_hessian 

Sparse finite-difference Lagrangian Hessian, recovered by graph coloring from the analytic Jacobian (Curtis, Powell & Reid 1974; Coleman & Moré 1983).

Why this exists

The models this targets — a direct-collocation transcription built from an FMU or a CasADi DaeBuilder — supply analytic first derivatives and no second derivatives. Today that leaves only crate::hess::lim_mem_quasi_newton, and on benchmarks/large_scale laptime the difference is stark: the exact Hessian converges in 30 iterations where limited-memory takes 246, and one mesh refinement later limited-memory does not converge at all.

But an analytic Jacobian is already enough to build the Hessian. The Lagrangian gradient

    ∇ₓL(x, y) = ∇f(x) + J_c(x)ᵀ y_c + J_d(x)ᵀ y_d

is available in closed form, so its directional derivative

    ∇²ₓₓL · d ≈ [ ∇ₓL(x + d, y) − ∇ₓL(x, y) ] / h

costs one gradient and one Jacobian evaluation — and with a known sparsity pattern, one probe recovers a whole group of structurally orthogonal columns at once rather than a single one.

Why it is affordable

Measured on laptime at N = 160: a Jacobian evaluation costs 5.4 ms against a 92.6 ms iteration, and the Hessian pattern has rho_max = 15 with a mean row of 5.68 — unchanged at N = 320 (POUNCE_HESS_PATTERN_CENSUS). The row width is set by the per-stage stencil, not the horizon, so the number of probes does not grow with the mesh. That is the property that makes the whole scheme viable: the cost per Hessian is constant in problem size while the iteration count it buys is the exact path’s.

The partition

Columns are grouped so that no two columns in a group share a row of the pattern (Curtis-Powell-Reid structural orthogonality). Probing group g with d = Σ_{j∈g} h_j e_j then gives, for every row i,

    w_i = Σ_{j∈g} H_ij h_j = H_ij h_j   for the unique j ∈ g with H_ij ≠ 0

so each entry is read off directly with no linear solve.

A star colouring (fd_hessian_coloring=star) lets an entry be read from either endpoint’s probe and so needs fewer groups: 76 → 42 on the Jacobian-derived pattern here, 17 → 16 on the declared one. Its recovery is algebraically exact — overlapping_cliques_are_validated_not_assumed verifies that by recovering a known matrix through it.

And it is still the wrong choice on a dense pattern, which is why CPR is the default. On laptime, star colouring over the Jacobian-derived pattern takes 404 iterations to an objective of 65.368334 where CPR takes 38 to 65.371106.

The cause is not group size — the measurement rules that out, since declared/star packs the largest groups of the four (580 columns per probe against jacobian/cpr’s 122) and converges in 30 iterations with the exact objective:

pattern / colouring	groups	cols per group	result
declared / cpr	17	546	Optimal, 30 it
declared / star	16	580	Optimal, 30 it
jacobian / cpr	76	122	Optimal, 38 it
jacobian / star	42	221	Acceptable, 404 it, wrong objective

The cause is the finite-difference remainder. Direct-recovery theory assumes exact Hessian-vector products; a forward difference also carries ½ Σ_{m,p ∈ g} T_imp h_m h_p into row i, where T is the third derivative. T_imp ≠ 0 needs i, m and p in a common constraint’s support, hence H_im ≠ 0 and H_ip ≠ 0. CPR’s distance-2 property forbids two such columns in one group, so those cross terms vanish structurally. A star colouring only guarantees the single-neighbour property for the pair being recovered, so the cross terms survive — and they matter exactly when the pattern is dense (rho_max 59 here against the declared pattern’s 15).

So star colouring is safe on a sparse declared pattern and unsafe on a Jacobian-derived one, which is the mode most models need. It stays opt-in. Every colouring is additionally validated entry by entry before use, unconditionally — that check began as a debug_assert, which is compiled out in release, and a Hessian that is wrong but plausible is the failure this module is most exposed to.

The pattern

Two sources, selected by fd_hessian_pattern:

• declared — the TNLP’s own Hessian sparsity, when it declares one (every .nl does, through AMPL’s AD). Requires no values, only the structure call, so it is available to a model that cannot evaluate second derivatives.

• jacobian — derived as (O ⊗ O) ∪ ⋃_j supp(∇g_j) ⊗ supp(∇g_j), where O is the objective’s nonlinear variables. It needs nothing beyond the Jacobian pattern every TNLP must declare, plus the objective linearity the TNLP will state. This is a strict superset of the true pattern, which is safe (a superset costs extra groups, never a wrong answer) but not free: on laptime it is 146 267 nonzeros against the true 28 000.

  The O ⊗ O term is not optional. The constraint Jacobian says nothing about ∇²f, so without it a model whose objective couples two variables that never share a constraint row gets a pattern that is a subset of the truth, and the recovery drops that curvature silently.

A superset is always safe; a subset would silently drop curvature, so there is no fallback that guesses. That sentence was written before the code honoured it: the objective clique’s fallback read the first ∇f’s values, which for f = x₀x₁ at the origin is the zero vector, so the clique came back empty and the pattern was a subset after all — and a different subset from a different starting point. Every level of the fallback is structural now: stated objective linearity, else the nonlinear-variable set N, else all n. The last two are conservative and can cost a great many probes; FdStats::objective_clique_widened reports when one of them was taken. gh#823 review (@srikanth-gm).

Probe cost scales with the Hessian’s row width, not with the model’s size. laptime needs 17 groups because its per-stage stencil makes rho_max = 15; that is a property of that transcription and does not generalise. A model with rho_max = 176 needs ~181 groups and therefore ~180 gradient-plus-Jacobian evaluations per Hessian, which can be far more expensive per iteration than the limited-memory path it replaces — measured on a 60k-variable model in gh#823 review. Read rho_max from POUNCE_FD_HESSIAN_DEBUG before assuming this mode is affordable.

Structs

FdHessianUpdater

FdStats

Enums

FdColoring

How columns are grouped into probes.

FdPatternSource

Where the Hessian sparsity pattern comes from.