Subproblem solvers

All solve

\[\min_{\|s\| \le \Delta} \; g^\top s + \tfrac12 s^\top H s\]

and return active::Bool, true when ‖s‖ = Δ. That flag is worth propagating: the fraction of iterations on which the constraint binds is the observable that separates radius mechanisms whose first-order behaviour is identical.

solver	needs	handles indefinite `H`?	size
`SteihaugCG`	`B * v` only	yes, stops on negative curvature	any
`ExactMS`	dense eigendecomposition	yes, including the hard case	`n ≤ nmax`
`KrylovCG`	`B * v`	no — assumes `H ≻ 0`	any
`KrylovCR`	`B * v`	yes, stops on negative curvature	any

`SteihaugCG` and the Cauchy point

The first CG direction is -g. A step that leaves the trust region on CG iteration 1 is therefore exactly the Cauchy point, and the model Hessian has had no influence on its direction whatsoever.

This is not a corner case. When the radius is small — a low μ_max in RGradCapped, a low ζ in RDFO — it happens at every iteration, and the method is gradient descent in disguise. Nothing in ρ, ‖g‖ or the radius trace reveals it. cg_step_info does:

using TrustRegionRadius, ADNLPModels, NLPModels

nlp = ADNLPModel(x -> (1 - x[1])^2 + 100(x[2] - x[1]^2)^2, [-1.2, 1.0])
x = [-1.2, 1.0]
g = grad(nlp, x)
info = cg_step_info(SteihaugCG(), ExactHessian(), nlp, x, g, 0.01)
(info.cg_iters, info.active, info.cos_cauchy)

(1, true, 1.0)

A radius small enough to truncate CG on its first iteration returns the Cauchy point: cg_iters == 1, the constraint active, and cos(s, -g) == 1, so the model Hessian played no part in choosing the direction.

cg_iters == 1 && active together with cos_cauchy ≈ 1 certify the degeneration.

`ExactMS` and the hard case

Uses the Moré–Sorensen characterisation: s* is optimal iff there is λ ≥ 0 with

\[(B + \lambda I)s^* = -g, \qquad B + \lambda I \succeq 0, \qquad \lambda(\Delta - \|s^*\|) = 0 .\]

‖s(λ)‖ decreases strictly in λ, so the root is bracketed and bisected.

The hard case — g orthogonal to the eigenspace of λ_min(B), where no λ attains ‖s(λ)‖ = Δ — is handled explicitly by adding the right multiple of that eigenvector. This matters with ExactHessian or SR1Model, where indefiniteness is exactly the information one wants to exploit; it cannot arise with a positive definite model, which is why solvers written for that case often omit it.

Requires an eigendecomposition, hence n ≤ nmax (default 200). Larger problems raise rather than silently allocating an n × n array.

Adding a solver

struct MySolver <: SubproblemSolver end

function TrustRegionRadius.solve_subproblem!(::MySolver, model::ModelHessian, nlp,
                                             x, g, Δ, s, Hs, ws::SubWorkspace;
                                             curv = nothing)
    B = hessian_op(model, nlp, x)
    # ... write the step into s, using ws.r, ws.d, ws.cand, ws.Hd as scratch ...
    return norm(s) ≥ (1 - 1e-8) * Δ     # active?
end

Not executed

An extension template with an elided body; shown for reference, not run.

Take the curvature from model, not from nlp: that is what lets one solver serve every model Hessian.

Three optional pieces:

ws is a caller-owned SubWorkspace, so nothing is allocated per iteration. An eight-argument convenience method allocates one for you, which is right at the REPL and wrong in the loop.
Hs is an output. Declare returns_hprod(::MySolver) = true if you leave B·s in it — CG forms it incrementally anyway, and the solver then skips its own model_hprod!, saving one Hessian-vector product per iteration.
curv carries (λ_min, v_min) when the solver has already computed it. Declare needs_eigenvector(::MySolver) = true to be given the eigenvector too, rather than recomputing an eigendecomposition the solver has just done.

API

TrustRegionRadius.SubproblemSolver — Type

SubproblemSolver

Abstract supertype for trust-region subproblem solvers. A concrete subtype implements

solve_subproblem!(sub, model, nlp, x, g, Δ, s, Hs, ws; curv = nothing) -> Bool

solving min_{‖s‖ ≤ Δ} gᵀs + ½ sᵀHs and writing the step into s.

source

TrustRegionRadius.SubWorkspace — Type

SubWorkspace(x)

Scratch vectors for the subproblem solvers, allocated once per TRSolver and parametric in the problem's vector type.

Four buffers: the CG residual r, the direction d, the trial iterate cand, and Hd = B·d. EigenPoint reuses cand and Hd after the inner solver has finished with them, and the trust-region loop reuses cand after solve_subproblem! has returned, to form the two residual diagnostics without allocating.

It also carries iters, the inner-iteration count of the last solve. That count is what distinguishes a step the model shaped from the Cauchy point: truncated CG returning after one iteration has taken the step along -g and the model Hessian has influenced nothing. 0 for solvers with no meaningful iteration count.

source

TrustRegionRadius.SteihaugCG — Type

SteihaugCG(; χ = 0.1, θ = 0.5, max_iters = 100)

Steihaug–Toint truncated conjugate gradient. Stops on negative curvature, on reaching the boundary, or when the residual satisfies the forcing condition ‖r‖ ≤ min(χ, ‖g‖^θ)‖g‖.

Needs only B * v, so it scales to large problems with any operator-backed model. Accumulates B·s alongside s, so the predicted reduction costs no extra product.

The Cauchy point

The first CG direction is -g. A step that leaves the region on CG iteration 1 is therefore exactly the Cauchy point, and the model Hessian has had no influence on its direction. When the radius is small — a low μ_max in RGradCapped, a low ζ in RDFO — this happens at every iteration and the method is gradient descent in disguise. cg_iters in the return of cg_step_info reveals it.

source

TrustRegionRadius.ExactMS — Type

ExactMS(; nmax = 200, tol = 1e-12, maxit = 100)

Exact step by the Moré–Sorensen characterisation: s* is optimal iff there is λ ≥ 0 with

(B + λI)s* = −g,   B + λI ⪰ 0,   λ(Δ − ‖s*‖) = 0.

‖s(λ)‖ decreases strictly in λ, so the root is bracketed and bisected.

Handles indefinite B, including the hard case (g orthogonal to the eigenspace of λ_min(B), where no λ attains ‖s(λ)‖ = Δ), by adding the right multiple of a single eigenvector from that eigenspace. The previous implementation added τ to every index within tolerance of λ_min, so with a degenerate leftmost eigenvalue of multiplicity m it returned a step of norm √(‖s̄‖² + mτ²) ≠ Δ — and degenerate λ_min is common on symmetric test problems, which is precisely where the hard case arises.

Requires a dense eigendecomposition, hence n ≤ nmax; larger problems raise rather than silently allocating an n × n array.

source

TrustRegionRadius.KrylovCG — Type

KrylovCG(; atol = 1e-6, rtol = 1e-6, itmax = 0)

Krylov.cg with the radius keyword (Steihaug-style truncation). Assumes positive definite B; use KrylovCR otherwise.

source

TrustRegionRadius.KrylovCR — Type

KrylovCR(; atol = 1e-6, rtol = 1e-6, itmax = 0)

Krylov.cr with the radius keyword; handles indefinite B, stepping to the boundary on negative curvature exactly as SteihaugCG does.

This was `KrylovCGLanczos`, and that wrapper never ran

The type was previously generated from Krylov.cg_lanczos, which has no radius keyword — Base.kwarg_decl on it lists check_curvature and no radius in any released version — so every call raised a MethodError before doing any work. The docstring's claim that it solved the constrained subproblem was therefore never exercised.

Krylov.cr is the routine in the same package that does accept radius and handles indefinite B, so it is what the type always meant. The alternative, keeping cg_lanczos and rescaling its unconstrained iterate down to ‖s‖ = Δ, was rejected: a scaled Newton-ish direction is not the trust-region step, and on indefinite B the iterate being rescaled can already have diverged.

KrylovCGLanczos remains as a deprecated alias.

source

TrustRegionRadius.solve_subproblem! — Function

solve_subproblem!(sub, model, nlp, x, g, Δ, s, Hs, ws; curv = nothing) -> Bool

Solve the trust-region subproblem, writing the step into s.

sub: the solver
model: ModelHessian supplying the curvature
nlp: the problem, for hessian_op
x, g: current iterate and gradient
Δ: radius
s: output buffer for the step
Hs: output buffer; holds B·s on return iff returns_hprod(sub)
ws: SubWorkspace scratch
curv: optional (λ_min, v_min) already computed by the caller

Returns true if ‖s‖ = Δ.

source

TrustRegionRadius.cg_step_info — Function

cg_step_info(sub, model, nlp, x, g, Δ) -> (; cg_iters, active, cos_cauchy)

Run SteihaugCG and report how the step was produced:

cg_iters: CG iterations taken
active: whether the step ended on the boundary
cos_cauchy: cos(s, −g), which is 1 exactly when the step is the Cauchy point

cg_iters == 1 && active and cos_cauchy ≈ 1 together certify that the model Hessian played no part in choosing the direction.

Now a thin wrapper over _steihaug! rather than a second copy of the recurrence. Allocating, and intended for diagnostics rather than the solver loop.

source

TrustRegionRadius.returns_hprod — Function

returns_hprod(sub) -> Bool

Whether solve_subproblem! leaves B·s in the Hs buffer. When true the solver skips its own model_hprod! for the predicted reduction, saving one Hessian-vector product per iteration.

source

TrustRegionRadius.needs_eigenvector — Function

needs_eigenvector(sub) -> Bool

Whether the subsolver wants the leftmost eigenvector in curv. The solver then requests lambda_min_estimate(...; vector = true) when it is computing the curvature anyway, so the estimate is paid for once rather than twice.

source

Subproblem solvers

SteihaugCG and the Cauchy point

ExactMS and the hard case

Adding a solver

API

`SteihaugCG` and the Cauchy point

`ExactMS` and the hard case