Getting started

Installing

using Pkg
Pkg.add(url = "https://github.com/USER/TrustRegionRadius.jl")

For development, clone and Pkg.develop(path = ".").

A first solve

using TrustRegionRadius, ADNLPModels, NLPModels

rosen(x) = (1 - x[1])^2 + 100(x[2] - x[1]^2)^2
nlp = ADNLPModel(rosen, [-1.2, 1.0])

stats = tr_solve(nlp; rule = RDelta())
(stats.status, stats.solution, stats.iter)

(:first_order, [0.9999999999982704, 0.999999999994377], 32)

tr_solve returns a TRResult, which is JSO's GenericExecutionStats, so SolverBenchmark and the rest of the ecosystem accept it unchanged. Evaluation counts live on the model — neval_obj(nlp), neval_grad(nlp), neval_hprod(nlp) — not on the result.

Choosing the axes

stats = tr_solve(nlp;
    rule      = RGrad(μ = 1.0),      # how Δ moves
    model     = SR1Model(mem = 5),   # what curvature the model reports
    subsolver = SteihaugCG(),        # how the subproblem is solved
    params    = TRParams(η = 0.1, η1 = 0.1, η2 = 0.9, Δ0 = 1.0, tol = 1e-8))
(stats.status, stats.iter)

(:first_order, 106)

Defaults are RDelta(), ExactHessian(), SteihaugCG(). The threshold keywords are ASCII (η1, η2, Δ0); the subscript spellings η₁, η₂, Δ₀ are accepted as aliases and readable as properties either way.

The fourth axis, the sampling rule, belongs to the oracle rather than to tr_solve, because it decides what the oracle returns:

prob = LogisticRegression(K = 5, M = 2_000)

exact   = tr_solve(FullBatchNLP(prob);                       model = BHHHModel())
sampled = tr_solve(FiniteSumNLP(prob, RadiusProportional()); model = BHHHModel())
(exact.status, exact.iter, sampled.status, sampled.iter)

(:first_order, 7, :first_order, 22)

tr_solve picks the solver from the oracle: DeterministicTRSolver, FiniteSumTRSolver or ExpectationTRSolver. See Problem classes for what each one does differently and why the three are kept apart.

The three thresholds

η decides whether a step is accepted; η₁ and η₂ decide only how the radius is scaled, and are the only two a rule ever sees:

p1 = TRParams(η1 = 0.1, η2 = 0.9)           # η defaults to η1: the classical algorithm
p2 = TRParams(η = 0.0, η1 = 0.1, η2 = 0.9)  # accept every step with positive pred. reduction
(p1.η, p2.η)

(0.1, 0.0)

They must satisfy 0 ≤ η ≤ η1 ≤ η2 < 1, checked in the constructor. Setting η < η₁ opens a band where a step is accepted while the radius still contracts. Thresholds and factors covers the consequences, and the matching convention 0 < γ₁ ≤ γ₂ < 1 < γ₃ on the scaling factors.

Tracing

trace = true attaches per-iteration trajectories to stats.solver_specific:

stats = tr_solve(nlp; rule = RGrad(), trace = true)
ss = stats.solver_specific

# Δₖ, ‖gₖ‖, f(xₖ) are state trajectories, of length k+1;
# ρₖ, ‖sₖ‖ and the two Bool records are per-iteration, of length k.
[k => length(ss[k]) for k in (:delta_trajectory, :grad_trajectory, :obj_trajectory,
                              :ratio_trajectory, :step_trajectory,
                              :active_trajectory, :accepted_trajectory)]

7-element Vector{Pair{Symbol, Int64}}:
    :delta_trajectory => 38
     :grad_trajectory => 38
      :obj_trajectory => 38
    :ratio_trajectory => 37
     :step_trajectory => 37
   :active_trajectory => 37
 :accepted_trajectory => 37

The one-entry difference is the alignment convention, and it is set out in full on the Diagnostics page.

A sampled run adds :grad_sample_trajectory, :obj_sample_trajectory, :samples_total, :sample_cap_hits and — where a truth exists — :true_grad_trajectory. A finite-sum run adds :full_batch_trajectory, marking the iterations at which Nₖ = M and the step was therefore exact.

The first three are one entry longer than the rest, since they have a value before the first iteration; align on the tail when plotting them together.

Two are worth recording even when nothing else is. :active_trajectory decides whether the constraint eventually stops binding: two mechanisms can have identical iteration counts and identical first-order behaviour while one keeps the constraint permanently active and the other does not, and nothing else in a standard diagnostic distinguishes them. :accepted_trajectory cannot be reconstructed from :ratio_trajectory and η₁ once η < η₁, so it is the only record of which iterations belong to $\mathcal{S}$.

When the constraint stopped binding

The tail fraction answers "is it still binding at the end?" but not "when did it stop?". For that, count the active iterations still ahead of each index:

inactivity_index answers it directly. It returns the first iteration after which the constraint never binds again, counted from zero, or nothing if it was still binding at the end.

(inactivity_index(stats), active_fraction(stats; tail = 0.1))

(32, 0.0)

A nothing there means the run never entered the regime where the local rate is available. If you want the staircase itself for a plot, count the active iterations still ahead of each index:

a = ss[:active_trajectory]
K = length(a)
R = [count(@view a[k+1:end]) for k in 0:K-1]     # Rₖ, a staircase down to zero
(first(R), last(R))

(14, 0)

R starts at the total number of active iterations, steps down by one at each active iteration, and is flat on inactive ones. It reaches zero at the last active iteration k*, after which the step is the unconstrained model minimiser. The range stops at K-1 on purpose: R_K = 0 for every run, so including it would make every curve reach zero and erase the distinction.

Because the run stops at ‖g‖ ≤ tol rather than at infinity, a positive tail is evidence of inactivity over the iterations observed, not a proof of eventual inactivity — and a tail of zero may mean only that the budget ran out first. Experiment 8 in benchmark/ plots this against Δ, ‖g‖ and ρ on a shared axis.

Statuses

:stalled deserves a word. Once the radius collapses far enough, the subtraction f(xₖ) − f(xₖ + sₖ) is pure rounding, ρ becomes noise, every step is rejected, and the radius contracts forever. This is a genuine numerical limitation of the criticality-anchored rules in floating point rather than a bug, and it is reported separately so that it is not silently folded into "failure".

Adding your own rule

One struct and two methods:

struct MyRule <: RadiusRule
    γ1::Float64
    γ2::Float64
    γ3::Float64
    function MyRule(; γ1 = 0.25, γ2 = 0.5, γ3 = 2.0)
        check_factors(:MyRule; γ1 = γ1, γ2 = γ2, γ3 = γ3)
        new(γ1, γ2, γ3)
    end
end

function TrustRegionRadius.update_radius!(r::MyRule, Δ::Float64, ρ::Float64, ::Bool,
                                          η1::Float64, η2::Float64,
                                          ::Float64, ::Float64, ::Float64)
    ρ < η1  && return r.γ1 * Δ        # unsuccessful: the result MUST be < Δ
    ρ >= η2 && return r.γ3 * Δ
    return Δ
end

TrustRegionRadius.asymptotic_regime(::MyRule) = :bounded_below

tr_solve(nlp; rule = MyRule(), params = TRParams(tol = 1e-8)).status

:first_order

Two things in that template are easy to get wrong, and both were wrong in an earlier version of this page.

check_factors takes ASCII keywords, γ1, γ2, γ3. The Unicode subscripts γ₁, γ₂, γ₃ are accepted as aliases on TRParams (η₁, η₂, Δ₀) but not by the rule validators, so check_factors(:MyRule; γ₁ = 0.25) raises MethodError: no method matching check_factors(::Symbol; γ₁::Float64).

There is deliberately no initial_radius method above. The fallback initial_radius(::RadiusRule, Δ0::Float64, ::Float64) already returns Δ0, and adding an unannotated initial_radius(::MyRule, Δ0, g_norm) = Δ0 makes the two ambiguous rather than overriding one:

MethodError: initial_radius(::MyRule, ::Float64, ::Float64) is ambiguous.

Define it only when your rule ignores Δ0, as the μ-scaled rules do, and annotate the argument types when you do.

Then tr_solve(nlp; rule = MyRule()). update_radius! is exported, so extending it needs the qualified name; using brings it into scope without making it extensible. Arguments the rule ignores are written as bare type annotations.

Add reset_rule! if the struct carries mutable state, and needs_retrospective(::MyRule) = true if it should be driven by ρ̃ rather than ρ — the solver then computes the retrospective ratio, at the cost of one extra Hessian-vector product per accepted iteration. Add a validate_thresholds method if the rule cannot accept every admissible (η, η₁, η₂).

The three obligations

Contract on unsuccessful iterations. If ρ < η₁ can leave Δ unchanged, the rule falls outside weak admissibility, the convergence theorem of Part I no longer covers it, and — since a rejected step changes neither the model nor the iterate — the solver re-solves an identical subproblem until the budget runs out. test/test_thresholds.jl checks this for every rule in the package at five values of ρ below η; add new rules to it.

Return a strictly positive radius. A rule anchored to ‖sₖ‖ can collapse to zero. That is what a Δmin field is for.

Declare the regime. asymptotic_regime and, where they apply, needs_retrospective and is_criticality_anchored are what the experiment scripts group by.

The arguments in order

update_radius!(rule, Δ, ρ, accepted, η1, η2, s_norm, g_norm_old, g_norm_new)

accepted is third. It is passed rather than recomputed because a rule cannot derive it: acceptance is decided by η, which the rule never receives, and a retrospective rule is handed ρ̃ in the ρ slot and so never sees ρ at all.

Note the two gradient-norm arguments: g_norm_old is ‖gₖ‖ before the accept/reject decision and g_norm_new is ‖gₖ₊₁‖ after it. Rules anchored to criticality at the current iterate (RDFO) use the former; rules that set the next radius from the next gradient (RGrad and relatives) use the latter. Getting this backwards is the easiest mistake to make here, and it shifts every threshold silently.