Model Hessians

The solver never forms a matrix: it needs only the operator action H·v, which is what hessian_op returns. This keeps large problems feasible and lets quasi-Newton models be supplied as LinearOperators without changing the solver loop.

model	curvature	can be indefinite?
`ExactHessian`	`∇²f(xₖ)`	yes
`LBFGSModel`	limited-memory BFGS	no — enforced positive definite
`SR1Model`	symmetric rank one	yes, and that is the point
`ScaledIdentity`	`c·I`, i.e. none	no
`SPDTarget`	constructed, minimiser pinned	no — positive definite by design
`BHHHModel`, `BHHH2Model`	outer product of scores	no — see Outer-product Hessians
`GaussNewtonModel`	`JᵀJ/N`	no — same

The last three are restricted to a problem class: required_problem declares the narrowest class they are defined for, and the solver constructor rejects the rest.

Which to use

ExactHessian for small problems and for any experiment where the model should not be a confounder. Note it does not satisfy the secant equation, so the Fan–Pan–Song route to ρ̃ → 1 does not apply to it; the Bastin route, which needs asymptotic second-order coherence instead, does.

LBFGSModel never reports negative curvature. On a problem whose true Hessian is indefinite along the trajectory this is invisible to every first-order diagnostic: ρ stays healthy, ‖g‖ decreases, and the limit can still fail second-order optimality.

SR1Model may be indefinite, which is what lets a subsolver that exploits negative curvature escape a saddle a positive-definite model would converge to. Pair it with SteihaugCG (which detects dᵀHd ≤ 0) or ExactMS (which handles the hard case).

ScaledIdentity is a diagnostic instrument, not a practical model: with it the step is -g/c truncated to the region, so the method is exactly gradient descent and any difference between radius rules is attributable to the rule alone.

`SPDTarget`

A two-dimensional construction whose unconstrained model minimiser is pinned to target while H ≻ 0. With d = target − x, u = d/‖d‖, v = u^⊥:

\[H = a\,uu^\top + b\,(uv^\top + vu^\top) + c\,vv^\top, \qquad a = -\frac{g^\top u}{\|d\|},\quad b = -\frac{g^\top v}{\|d\|},\quad c = \frac{b^2}{a} + \lambda_\perp .\]

This has minimiser target for every c, and is positive definite exactly when a > 0, equivalently when

\[\varphi(x) = g^\top(\text{target} - x) < 0 ,\]

that is, the target lies downhill. That condition is necessary for any positive definite model with that minimiser, so a DomainError is not an artefact of the construction: no such model exists at x. Query it directly with phi_target.

Point a SPDTarget at a saddle to obtain a run in which every hypothesis of the first-order theory holds, every ρ is successful, ‖g‖ → 0 — and the limit is a saddle. The failure is the model, and no radius rule can repair it.

using TrustRegionRadius, ADNLPModels

nlp = ADNLPModel(p -> p[1]^4 - p[1]^3 + (0.25 - p[1]/2)*p[2]^2 + p[2]^4/4, [-0.5, 0.6])
stats = tr_solve(nlp; rule = RDelta(), model = SPDTarget(target = [0.75, 0.0]))
stats.solution     # the saddle

2-element Vector{Float64}:
  0.75
 -2.7755575615628914e-17

Adding a model

struct MyModel <: ModelHessian end

TrustRegionRadius.hessian_op(::MyModel, nlp, x) = # anything supporting B * v
TrustRegionRadius.dense_hessian(::MyModel, nlp, x) = # Matrix, for ExactMS and diagnostics

Add reset_model! and update_model! if the model carries state; both default to no-ops. update_model! receives the secant pair (sₖ, yₖ = gₖ₊₁ − gₖ) on accepted steps only.

Two further traits are worth declaring:

TrustRegionRadius.reports_negative_curvature(::MyModel) = false  # if B ⪰ 0 by construction
TrustRegionRadius.required_problem(::MyModel) = LikelihoodProblem  # if it needs one

Not executed

The two blocks above are extension templates with elided right-hand sides. They are shown for reference and are not run when the documentation is built.

The first makes the solver warn when the model meets SecondOrder or tol_H, where it would be an expensive no-op reporting :second_order at a saddle. The second is checked at solver construction. Both default to the permissive answer.

API

TrustRegionRadius.ModelHessian — Type

ModelHessian

Abstract supertype for the model Hessian H_k used in the trust-region model

m_k(s) = f(x_k) + g_kᵀ s + ½ sᵀ H_k s.

See ExactHessian, LBFGSModel, SR1Model, ScaledIdentity, SPDTarget.

source

TrustRegionRadius.ExactHessian — Type

ExactHessian()

The true Hessian ∇²f(x_k), as an operator via NLPModels.hess_op.

Note that the exact Hessian does not satisfy the secant equation H_{k+1} s_k = y_k, so retrospective rules relying on that condition (the Fan–Pan–Song route to ρ̃ → 1) do not apply to it; the Bastin route, which needs asymptotic second-order coherence instead, does.

source

TrustRegionRadius.LBFGSModel — Type

LBFGSModel(; mem = 5, T = Float64)

Limited-memory BFGS model Hessian, backed by LinearOperators.LBFGSOperator.

The operator is stored in a typed field, Union{Nothing, LBFGSOperator{T}}, so hessian_op is type-stable and every mul! in the CG recurrence resolves at compile time. The previous op::Any silently defeated the M <: ModelHessian type parameter on TRSolver, whose whole purpose is that "every dispatch in the loop is resolved at compile time"; the comment justifying Any cited avoiding a hard dependency on LinearOperators, which the module already usings.

Pass T to match the problem's element type; TRSolver checks and raises on a mismatch.

The operator enforces positive definiteness, so the model never reports negative curvature — see reports_negative_curvature. On a problem whose true Hessian is indefinite along the trajectory this is invisible to every first-order diagnostic: ρ stays healthy, ‖g‖ decreases, and the limit can still fail second-order optimality.

source

TrustRegionRadius.SR1Model — Type

SR1Model(; mem = 5, T = Float64)

Symmetric rank-one model Hessian, backed by LinearOperators.LSR1Operator.

Unlike L-BFGS, SR1 may be indefinite — which is the point: it can represent negative curvature, so a subsolver that exploits it (truncated CG detecting dᵀHd ≤ 0, or an exact solver handling the hard case) can escape a saddle that a positive-definite model would converge to.

source

TrustRegionRadius.ScaledIdentity — Type

ScaledIdentity(; c = 1.0)

H_k = c·I. Carries no curvature information whatsoever.

Diagnostic instrument rather than a practical model: with this model the trust-region step is -g/c truncated to the region, so the method is exactly gradient descent and any difference between radius rules is attributable to the radius rule alone. Useful for isolating the Cauchy-point degeneration described in the survey.

source

TrustRegionRadius.SPDTarget — Type

SPDTarget(; target, λ⊥ = 1.0)

Two-dimensional construction whose unconstrained model minimiser is pinned to target, while H_k ≻ 0.

With d = target - x, u = d/‖d‖, v = u^⊥:

H = a·uuᵀ + b·(uvᵀ + vuᵀ) + c·vvᵀ,
a = -gᵀu/‖d‖,   b = -gᵀv/‖d‖,   c = b²/a + λ⊥.

This has minimiser target for every c, and is positive definite exactly when a > 0, equivalently when

φ(x) = gᵀ(target - x) < 0,

i.e. the target lies downhill. That condition is necessary for any SPD model with that minimiser, so a DomainError here is not an artefact of the construction: no such model exists at x.

Point a SPDTarget at a saddle to obtain a run in which every hypothesis of the first-order theory holds, every ρ is successful, ‖g‖ → 0, and the limit is nonetheless a saddle — the failure being the model, which no radius rule can repair.

This model evaluates the gradient

dense_hessian calls grad(nlp, x), so neval_grad counts gradients the algorithm never asked for. Since evaluation counts are the benchmark's cost measure, subtract model_grad_evals(solver) before reporting, or run this model only in the diagnostic experiments it exists for.

source

TrustRegionRadius.hessian_op — Function

hessian_op(model, nlp, x)

Return an object B supporting B * v and mul!(y, B, v), representing H_k at x.

source

TrustRegionRadius.dense_hessian — Function

dense_hessian(model, nlp, x) -> Matrix

Dense representation, for diagnostics on small problems only.

source

TrustRegionRadius.model_hprod! — Function

model_hprod!(model, nlp, x, v, Hv)

Write H_k · v into Hv, where H_k is the model Hessian.

Used for the predicted reduction -gᵀs - ½ sᵀH_k s. Using ∇²f here instead would make ρ measure agreement between the true function and a model the algorithm never minimised, so ρ would stop being the quantity the acceptance test and the radius rules are stated in terms of.

Subsolvers that already form B·s while computing the step declare returns_hprod, and the solver then skips this call entirely — see Subproblem/subproblem.jl.

source

TrustRegionRadius.update_model! — Function

update_model!(model, s, y)

Absorb the secant pair (s_k, y_k) with y_k = g_{k+1} - g_k, on accepted steps only. Default: no-op.

source

TrustRegionRadius.reset_model! — Function

reset_model!(model, n)

Clear any mutable state, for a problem of dimension n. Called at the start of every solve. Default: no-op.

source

TrustRegionRadius.phi_target — Function

phi_target(model, nlp, x) -> Float64

The descent functional φ(x) = ∇f(x)ᵀ(target − x). The SPDTarget construction exists at x if and only if φ(x) < 0.

source

TrustRegionRadius.reports_negative_curvature — Function

reports_negative_curvature(model) -> Bool

Whether the model is capable of reporting λ_min(B) < 0.

false for every model that is positive (semi)definite by construction: LBFGSModel, ScaledIdentity, SPDTarget, and the outer-product models BHHHModel, BHHH2Model, GaussNewtonModel.

This trait exists because four separate docstrings previously asserted the consequence — that SecondOrder over such a model gives τ ≡ ‖g‖ and a :second_order status certifying nothing — and one of them claimed the solver warned about it, while nothing in the solver did. TRSolver now consults this and warns once at construction. The failure mode is silent and the diagnosis is a one-line trait, so it should not be prose.

source

TrustRegionRadius.model_eltype — Function

model_eltype(model) -> Type or nothing

The element type a stateful model was constructed for, or nothing when the model carries no typed state. TRSolver checks it against the problem's element type and raises rather than letting a Float64 quasi-Newton operator be applied to Float32 vectors deep inside the CG recurrence.

source