Benchmarking

performance_profile(T; ntau = 500, taumax = nothing) -> (τ, prof)

Performance profile of Dolan & Moré (2002).

T[p, s] is the cost of solver s on problem p; use Inf or NaN for a failure. Returns the abscissa τ and the matrix

prof[i, s] = fraction of problems with r_{p,s} ≤ τ[i],
r_{p,s} = T[p,s] / min_{s'} T[p,s'].

prof[1, s] is the fraction of problems on which s is fastest (efficiency); the right-hand asymptote is the fraction it solves at all (reliability).

Example

τ, prof = performance_profile(T)
open(io -> profile_to_pgfplots(io, τ, prof, labels), "prof.tex", "w")

data_profile(N, dims; nkappa = 500, kmax = nothing) -> (κ, d)

Data profile of Moré & Wild (2009).

N[p, s] is the number of function evaluations solver s needed on problem p (Inf if unsolved) and dims[p] the dimension of problem p. The abscissa is the budget measured in units of n+1 evaluations — one simplex gradient — which is what makes problems of different size comparable.

Unlike the performance profile, the data profile is not normalised by the best solver, so it does not suffer the set-dependence noted above; it answers "what fraction is solved within this budget?" directly.

profile_to_pgfplots(io, x, y, labels; npoints = 60, styles = nothing, logx = true)

Write \addplot blocks ready to \input into a pgfplots axis.

The curves are subsampled onto a logarithmic grid (or a linear one if logx = false, appropriate for data profiles) so the emitted file stays small: a profile computed on a 2000-point grid becomes ~60 coordinates per curve with no visible difference.

styles, if given, is a vector of pgfplots option strings, one per curve.

Example

τ, prof = performance_profile(T)
open("prof_mech.tex", "w") do io
    profile_to_pgfplots(io, τ, prof, ["R-delta", "R-step", "R-DFO", "R-grad"];
                        styles = ["black,solid", "black,dashed",
                                  "gray,solid", "gray,dashdotted"])
end

then in LaTeX:

\begin{axis}[xmode=log, xlabel={$\tau$}, ylabel={$\pi_s(\tau)$}]
  \input{prof_mech.tex}
\end{axis}

run_matrix(problems, configs; cost = :iter, verbose = true)
    -> (T, stats_matrix)

Run every configuration on every problem.

Arguments

• problems: vector of thunks, each returning a fresh AbstractNLPModel. Thunks rather than models so that CUTEst handles are opened and finalised around a single run.
• configs: vector of TRConfig.

Keyword arguments

• cost: which quantity fills the matrix.

  cost	measures	valid on
  :iter	stats.iter	any
  :obj, :grad, :hprod	neval_*(nlp)	any
  :time	stats.elapsed_time	any
  :samples	cumulative term evaluations	sampled oracles only

  :iter is proportional to work only under FixedSample. Under any adaptive sampling rule, N_k varies per iteration and the same runs order differently by the two measures — which is the whole point of the sampling axis. Use :samples whenever the columns differ in sampling rule; it raises on a deterministic problem rather than returning a meaningless zero.

• verbose: print one line per problem as it completes.

Returns

• T[p, c]: the cost, or Inf where the run did not reach :first_order.
• stats_matrix[p, c]: the GenericExecutionStats, or nothing on error.

Runs terminating :first_order or :second_order count as solved; :max_iter, :stalled, :user, :exception and any thrown error give Inf. That convention keeps the profile honest: a solver that stops early with a large gradient has not solved the problem, whatever its iteration count.

A configuration rejected by the class checks — BHHHModel on a problem that is not a likelihood, a rule carrying N_max on a finite sum — raises from the solver constructor and is recorded as Inf with a warning, so one illegal cell does not abort the matrix.

Example

problems = [() -> ADNLPModel(x -> sum(x.^2), ones(5), name = "sphere")]
configs  = [TRConfig("R-delta"; rule = RDelta()),
            TRConfig("R-grad";  rule = RGrad())]
T, S = run_matrix(problems, configs)
summarise(T, [c.label for c in configs])

summarise(T, labels; io = stdout)

Print solved counts, median and mean cost per configuration.

The median is over solved problems only, so it is not comparable across columns with very different reliability — read it together with the solved count, or use a performance profile, which handles the two jointly.

TRConfig(label; rule, model, subsolver, params)

One cell of the experiment matrix: a named triple of radius rule, model Hessian, and subproblem solver, together with the solver parameters.

rule and model carry mutable state (μ for RGrad, the operator for L-BFGS/SR1), so run_matrix deep-copies them per run; a single TRConfig can therefore be reused across problems without contamination.

The solver constructor deep-copies its three axes as well, so the copy here is now belt-and-braces rather than load-bearing. It is kept deliberately: it makes this file correct on its own terms rather than by relying on an invariant maintained two modules away, and one deepcopy of a rule is free beside a solve.

The sampling rule is not a field here. It lives on the oracle, because it decides what the oracle returns — so a sweep over sampling rules varies the problem thunks, not the configurations:

problems = [() -> FiniteSumNLP(prob, r) for r in (FixedSample(64),
                                                  RadiusProportional(),
                                                  NormTest())]
configs  = [TRConfig("R-delta"; rule = RDelta())]
T, _ = run_matrix(problems, configs; cost = :samples)

Note that this puts the sampling rules down the rows, where performance_profile treats them as different problems rather than different solvers. To profile them as solvers, transpose T — or build one thunk per problem and one config per rule and accept that the sampling rule is then fixed across the matrix.

sweep_configs(name, values, mkrule; model, subsolver, params) -> Vector{TRConfig}

Build one TRConfig per parameter value, labelled name = value.

mkrule(v) must return a fresh rule for the value v.

configs = sweep_configs("ζ", [0.1, 1.0, 10.0], ζ -> RDFO(ζ = ζ))
T, _    = run_matrix(problems, configs)

A sweep of this kind is the right way to study ζ or μ_max: the survey's thresholds involve κ̄ = 8/λ*_min, which is a property of the solution and so cannot be chosen a priori, and only a sweep reveals where the transition sits on a given family of problems.

mkrule may return any RadiusRule, including a SecondOrder wrapper — in which case pair it with params = TRParams(tol_H = 1e-6) and a model that can report negative curvature, or the sweep measures an expensive no-op. Runs reaching a certified second-order point report :second_order, which run_matrix counts as solved.