Diagnostics

Quantities derived from a completed run. Two jobs: computing the thresholds the theory names but the solver cannot know, and doing the trace alignment once, here, rather than in every script that reads a trajectory.

Everything on this page needs trace = true.

The alignment convention

This is the one thing to get right, and it is easy to get wrong.

Entry j of every trajectory describes iteration j-1, counting from the initial point. The state trajectories carry one entry more than the per-iteration ones, because they also record the state the final iteration left behind:

length `k+1` (state)	length `k` (per-iteration)
`:delta_trajectory`, `:grad_trajectory`, `:obj_trajectory`	`:ratio_trajectory`, `:step_trajectory`
`:lambda_min_trajectory`, `:tau_trajectory`	`:active_trajectory`, `:accepted_trajectory`
`:dist_trajectory`, `:lambda_min_true_trajectory`	`:branch_trajectory`, `:cg_iters_trajectory`
`:hessian_norm_trajectory`	`:rho_tilde_trajectory`, `:gamma_trajectory`, `:xi_trajectory`, `:cos_cauchy_trajectory`

The two families agree at the head, so to pair a state with the iteration it produced, drop the last entry:

using TrustRegionRadius, ADNLPModels, NLPModels

nlp = ADNLPModel(x -> (1 - x[1])^2 + 100(x[2] - x[1]^2)^2, [-1.2, 1.0],
                 name = "Rosenbrock")
stats = tr_solve(nlp; rule = RGrad(μ = 1.0), model = ExactHessian(),
                 params = TRParams(tol = 1e-8), trace = true)

Δ = stats.solver_specific[:delta_trajectory]
s = stats.solver_specific[:step_trajectory]
(length(Δ), length(s))

(38, 37)

A radius can never be shorter than the step taken inside it, and that inequality is exactly what fails if the alignment is off by one:

all(s[j] <= Δ[j] * (1 + 1e-10) for j in eachindex(s))

true

Aligning on the tail instead pairs Δ_k with ‖s_{k-1}‖ and shifts every activity, ratio and step-versus-radius plot by one iteration. Prefer the helpers below, which do the alignment once and correctly.

Reading a trace

theta_trajectory returns θ_k = Δ_k/‖g_k‖, the ratio the local theory of Part II is about. For the gradient-scaled rules it is the multiplier μ_k.

θ = theta_trajectory(stats)
(length(θ), θ[1], θ[end])

(38, 1.0, NaN)

inactivity_index is the first iteration after which the constraint never binds again, counted from zero, or nothing if it was still binding at the end. active_fraction reports the fraction of a tail slice on which it did bind.

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

(32, 0.25)

branch_counts says which branch of the radius rule fired, and how often. On a rule with a half-step guard, :hold counts the iterations where the guard refused an expansion that the ratio alone would have allowed.

branch_counts(stats)

Dict{Symbol, Int64} with 4 entries:
  :expand   => 12
  :contract => 3
  :hold     => 19
  :shrink   => 3

Thresholds

kappa_bar is the inactivity threshold κ̄ of Condition C.Sg at a solution: the criticality-anchored rules reach eventual inactivity only if their parameter exceeds it. It is a property of the solution, so it cannot be checked before the run.

There are two conventions and they differ by a factor of two, which is the width of the interval most parameter sweeps resolve. Report which one produced a number.

(kappa_bar(nlp, [1.0, 1.0]),                                # 8/λ*, the default
 kappa_bar(nlp, [1.0, 1.0]; convention = :eigenvalue))      # 4/λ*

(20.03201278465082, 10.01600639232541)

kappa_bar_empirical is the largest realised ‖s_k‖/‖g_k‖ along the run, so comparing the two measures how loose the theoretical constant is:

kappa_bar_empirical(stats)

2.5039338811644374

Rates and hypotheses

observed_order fits the local convergence order by regressing log‖g_{k+1}‖ on log‖g_k‖, restricted by default to iterations beyond inactivity_index and to accepted steps. It returns NaN rather than a number when fewer than three usable points remain, with npts reporting the actual count, so "too short" stays distinguishable from "measured and small".

observed_order(stats)

(order = 1.1516582039782757, intercept = -2.7102041592046193, npts = 4, from = 32)

hypotheses_report measures the standing hypotheses of the local theory on the run rather than assuming them. A run whose rate looks wrong should be read here first.

hypotheses_report(stats)

(xi_final = 4.582473538140181e-14, xi_max_tail = NaN, gamma_final = 7.508712419958438e-12, gamma_max_tail = NaN, rho_final = 1.0000000001602591, cauchy_fraction = 0.3783783783783784)

radius_sums returns the three radius series of Part II. The third, Σ Δ_k²/M_k with M_k = L + max_{i≤k}‖H_i‖, is NaN unless the run was traced with hessian_norm = true, because there is then no recorded Hessian norm to divide by. Absence is reported as NaN rather than as zero, which would read as a converged sum.

radius_sums(stats)

(sum_delta = 282.49042284964963, sum_delta2 = 55214.0236444967, sum_delta2_over_M = NaN, n = 38)

withnorm = tr_solve(nlp; rule = RGrad(μ = 1.0), model = ExactHessian(),
                    params = TRParams(tol = 1e-8), hessian_norm = true, trace = true)
radius_sums(withnorm)

(sum_delta = 282.49042284964963, sum_delta2 = 55214.0236444967, sum_delta2_over_M = 36.65376654805906, n = 38)

n is reported beside the sums so that a partial sum is never mistaken for a converged one: on a run of forty iterations neither convergence nor divergence is established, and the honest reading is the shape of the partial sums rather than their value.

API

TrustRegionRadius.theta_trajectory — Function

theta_trajectory(stats) -> Vector{Float64}

The radius-to-criticality ratio θ_k = Δ_k/‖g_k‖, correctly aligned.

θ_k > κ̄ is the sufficient condition for the trust-region constraint to be inactive, and for the gradient-scaled rules θ_k is the multiplier μ_k itself, so this is the sequence that the local theory of Part II is about. Length k+1, one entry per iterate; pair it with the per-iteration arrays through θ[1:end-1].

source

TrustRegionRadius.inactivity_index — Function

inactivity_index(stats) -> Union{Int, Nothing}

The first iteration after which the trust-region constraint never binds again, counted from zero, or nothing if it was still binding at the end.

This is k_⋆ of Part II's inactivity theorem, and it is the quantity the mechanism is supposed to decide. Every statement about the rate belongs after it: a rate estimated across k_⋆ mixes the linear phase with the asymptotic one and lands between the two.

source

TrustRegionRadius.active_fraction — Function

active_fraction(stats; tail = 0.1) -> Float64

Fraction of the last tail of the run on which the constraint was active.

Zero for a configuration that reaches eventual inactivity, and one for a configuration trapped below its threshold. Part II predicts no intermediate value in the limit, so an intermediate value is a finding.

source

TrustRegionRadius.branch_counts — Function

branch_counts(stats) -> Dict{Symbol, Int}

How many times each branch of the radius rule fired. See last_branch for the vocabulary.

:expand counts the climb a gradient-scaled multiplier makes towards κ̄; :expand_capped counts the iterations on which the user's cap refused that climb, which is the mechanism by which a cap below κ̄ traps the run.

source

TrustRegionRadius.kappa_bar — Function

kappa_bar(λ_min::Real; convention = :neighbourhood) -> Float64
kappa_bar(nlp, x_star; convention = :neighbourhood, nmax, lanczos_k) -> Float64

The inactivity threshold κ̄ of Condition C.Sg at a solution.

κ̄ is the constant for which ‖s_k‖ ≤ κ̄‖g_k‖ near x*, and the criticality- anchored rules reach eventual inactivity only if their parameter exceeds it: ζ > κ̄ for RDFO and μ_max > κ̄ for RGradCapped. It is a property of the solution, so it cannot be checked before the run; this function computes it afterwards, which is what a threshold experiment needs in order to say where the observed transition sits relative to the predicted one.

The two conventions

The local lemma gives κ̄ = 4/m where m lower-bounds ∇²f on a neighbourhood of x*, and there are two ways to quote it.

:neighbourhood (default) takes the concrete choice m = λ*_min/2, giving κ̄ = 8/λ*_min. It is what a fixed neighbourhood costs, and the extra factor of 2 is removable. This is the default because it is the convention Part III of the survey is written in throughout, and a package whose default disagrees with the paper it accompanies is a trap.
:eigenvalue takes m = λ_min(∇²f(x*)), giving κ̄ = 4/λ*_min. The neighbourhood can be shrunk until m is as close to λ*_min as one likes, and every statement in which κ̄ is compared with a parameter is asymptotic, so this is the constant that is actually attained.

Report which one a number was produced under. A threshold quoted in one convention and compared against a parameter chosen in the other is off by a factor of two, which is the width of the interval most sweeps resolve.

κ = kappa_bar(nlp, x_star)            # 8/λ*_min
tr_solve(nlp; rule = RDFO(ζ = 2κ))    # comfortably above the threshold

source

TrustRegionRadius.kappa_bar_empirical — Function

kappa_bar_empirical(stats; inactive_only = true, tail = 1.0) -> Float64

The largest realised ‖s_k‖/‖g_k‖ along a run: the constant of Condition C.Sg that the run actually exhibited.

inactive_only = true restricts to iterations on which the trust-region constraint was inactive, where the step solves the model equation and the ratio is the one the sharp bound 1/λ_min(H_k) describes. On active iterations the ratio is Δ_k/‖g_k‖, which says more about the rule than about the geometry.

tail keeps only the final fraction of the run, so that the constant is measured where the local hypotheses hold rather than during the global phase.

Comparing this with kappa_bar is how one measures the looseness of the theoretical constant, which Part II gives as a factor of 4 asymptotically.

source

TrustRegionRadius.observed_order — Function

observed_order(stats; after_inactivity = true, last = 20, use = :grad)
    -> (; order, intercept, npts, from)

Local convergence order, estimated by regressing log‖g_{k+1}‖ on log‖g_k‖ (or on ‖x_k − x_ref‖ with use = :dist, when the trace carries it).

after_inactivity = true restricts the regression to iterations beyond inactivity_index and to accepted steps. This is not a refinement: a regression that includes binding iterations mixes the linear phase with the asymptotic one and returns an order between the two, which is an artefact.

Returns NaN order when fewer than three usable points remain, which is the honest outcome on a run that reaches machine precision in nine iterations.

The claim under test is that this number depends on the model and the solver but not on the radius mechanism, the mechanism deciding only from.

source

TrustRegionRadius.hypotheses_report — Function

hypotheses_report(stats) -> NamedTuple

The standing hypotheses of Part II §7, as measured on this run rather than assumed:

xi_final, xi_max_tail: the realised inexactness ‖H_ks_k + g_k‖/‖g_k‖, which Condition C.Inx requires to tend to zero;
gamma_final, gamma_max_tail: the Dennis–Moré residual, which the superlinear rate requires to tend to zero;
rho_final: the acceptance ratio, which the local theory says tends to one;
cauchy_fraction: the fraction of iterations on which the step was the Cauchy point to within 1e-12, i.e. on which the model Hessian did not influence the direction.

A run whose rate looks wrong should be read here first: an order estimate is only about the mechanism if these say the hypotheses held.

source

TrustRegionRadius.radius_sums — Function

radius_sums(stats; L = 0.0) -> (; sum_delta, sum_delta2, sum_delta2_over_M, n)

The three radius series of Part II, on one run.

sum_delta = Σ_k Δ_k, which the step-driven rules make finite in the local regime;
sum_delta2 = Σ_k Δ_k²;
sum_delta2_over_M = Σ_k Δ_k²/M_k with M_k = L + max_{i≤k}‖H_i‖, which is the quantity the criticality-anchored propositions actually prove finite, and is NaN unless the run was traced with hessian_norm = true.

The third is not a refinement of the second. Σ Δ_k² is equivalent to it only while {‖H_k‖} is bounded, and that boundedness is exactly the hypothesis the M_k form was introduced to drop.

L is the Lipschitz constant of the gradient. It is rarely known, and leaving it at zero changes the constant, not whether the series converges; say which was used when reporting a number.

Sums are over the whole trajectory including the final radius, and are reported with n so that a partial sum is never mistaken for a converged one: on a run of forty iterations neither convergence nor divergence is established, and the honest reading is the shape of the partial sums, not their value.

source