The dogleg method.

A simple and efficient heuristic for solving local trust region problems. Does not necessarily find the optimum, but improves on the Cauchy point.

struct ForwardDiffBuffer{TX, TF, TR, TC}

A buffer for wrapping a (nonlinear) mapping f to calculate the Jacobian with ForwardDiff.jacobian.

Generalized eigenvalue solver (Adachi et al 2017).

struct LocalModel{TF<:Real, TG<:(AbstractVector{T} where T), TB<:(AbstractMatrix{T} where T)}

A quadratic model for a minimization problem, relative to the origin, with the objective function

\[m(p) = f + p' g + 1/2 g' B g\]

SolverStoppingCriterion(residual_norm_tolerance)

Stopping criterion for trust region colver.

Arguments:

• residual_norm_tolerance: convergence is declared when the norm of the residual is below this value

TrustRegionParameters(η, Δ̄)

Trust region method parameters.

• η: trust reduction threshold
• Δ̄: initial trust region radius

struct TrustRegionResult{T<:Real, TX<:AbstractArray{T<:Real, 1}, TFX}

Fields

• Δ

  The final trust region radius.

• x

  The last value (the root only when converged).

• fx

  f(x) at the last x.

• residual_norm

  The Euclidean norm of the residual at x.

• converged

  A boolean indicating convergence.

• iterations

  Number of iterations (≈ number of function evaluations).

ForwardDiff_wrapper(f, n, jacobian_config)

Wrap an $ℝⁿ$ function f in a callable that can be used in trust_region_solver directly. Remaining parameters are passed on to ForwardDiff.JacobianConfig, and can be used to set eg chunk size.

Non-finite residuals will be treated as infeasible.

julia> f(x) = [1.0 2; 3 4] * x - ones(2)
f (generic function with 1 method)

julia> ff = ForwardDiff_wrapper(f, 2)
AD wrapper via ForwardDiff for ℝⁿ→ℝⁿ function, n = 2

julia> ff(ones(2))
(r = [2.0, 6.0], J = [1.0 2.0; 3.0 4.0])

julia> trust_region_solver(ff, [100.0, 100.0])
Nonlinear solver using trust region method converged after 9 steps
  with ‖x‖₂ = 3.97e-15, Δ = 128.0
  x = [-1.0, 1.0]
  r = [-1.78e-15, 3.55e-15]

_check_residual_Jacobian(x, fx)

Checks residual and Jacobian and throws an appropriate error if they are not both finite.

Internal, each evaluation of f should call this.

_factorize(J)

Return (F, is_valid_F) where

1. F is a form (usually a factorization) of the argument that supports F \ r for vectors r,

2. is_valid_F is a boolean indicating whether F can be used in this form (eg false for a singular factorization).

cauchy_point(Δ, model)

Find the Cauchy point (the minimum of the linear part, subject to $\| p \| ≤ Δ$) of the problem.

Return three values:

1. the Cauchy point vector pC,

2. the (Euclidean) norm of pC

3. a boolean indicating whether the constraint was binding.

Caller guarantees non-zero gradient.

dogleg_boundary(Δ, D, pC)
dogleg_boundary(Δ, D, pC, pC_norm)

Finds au such that

\[\| p_C + τ D \|_2 = Δ\]

pC_norm is \| p_C \|_2^2, and can be provided when available.

Caller guarantees that

1. pC_norm ≤ Δ,

2. norm(pC + D, 2) ≥ Δ.

3. dot(pC, D) ≥ 0.

These ensure a meaningful solution 0 ≤ τ ≤ 1. None of them are checked explicitly.

dogleg_implementation(Δ, model, pU, pU_norm)

Implements the branch of the dogleg method where it is already determined that the unconstrained optimum is lies outside the Δ ball. Returns the same kind of results as solve_model.

ellipsoidal_norm(x, _)

Ellipsoidal norm $\| x \|_B = x'Bx$.

ges_kernel(Δ, model, S)

Kernel for the generalized eigenvalue solver (methods specialize on the ellipsoidal norm matrix S). Solves the M̃(λ) pencil in Adachi et al (2017), eg Algorithm 5.1.

Returns

• λ, the largest eigenvalue,

• the gap to the next eigenvalue (see note below)

• y1 and y2, the two parts of the corresponding generalized eigenvalue,

All values are real, methods are type stable.

When theorerical assumptions are violated, gap will be non-finite and a debug statement is emitted. No other values should be used in this case.

local_residual_model(r, J)

Let r and J be residuals and their Jacobian at some point. We construct a local model as

\[m(p) = 1/2 \| J p - r \|^2_2 = \| r \|^2_2 + p' J ' r + 1/2 p' J' J p\]

model_reduction(model, p)

Reduction of model objective at p, compared to the origin.

reduction_ratio(model, p, p_objective)

Ratio between predicted (using model, at p) and actual reduction (using p_objective, the minimand at p).

Will return an arbitrary negative number for infeasible coordinates.

residual_minimand(r)

Convert a residual from a system of equations (in a rootfinding problem) to an objective for minimization.

solve_model(method, Δ, model)

Optimize the model with the given constraint Δ using method. Returns the following values:

• the optimum (a vector),
• the norm of the optimum (a scalar),
• a boolean indicating whether the constraint binds.

trust_region_solver(f, x; parameters, local_method, stopping_criterion, maximum_iterations, Δ, debug)

Solve f ≈ 0 using trust region methods, starting from x.

f is a callable (function) that

1. takes vectors real numbers of the same length as x,

2. returns a value with properties residual and Jacobian, containing a vector and a conformable matrix, with finite values for both or a non-finite residual (for infeasible points; Jacobian is ignored). A NamedTuple can be used for this, but any structure with these properties will be accepted. Importantly, this is treated as a single object and can be used to return extra information (a “payload”), which will be ignored by this function.

Returns a TrustRegionResult object.

Keyword arguments (with defaults)

• parameters = TrustRegionParameters(): parameters for the trust region method

• local_method = Dogleg(): the local method to use

• stopping_criterion = SolverStoppingCriterion(): the stopping criterion

• maximum_iterations = 500: the maximum number of iterations before declaring non-convergence,

• Δ = 1.0, the initial trust region radius

• debug = nothing: when ≢ nothing, a function that will be called with an object that has properties iterations, Δ, x, fx, converged, residual_norm.

trust_region_step(parameters, local_method, f, Δ, x, fx)

Take a trust region step using local_method.

f is the function that returns the residual and the Jacobian (see trust_region_solver).

Δ is the trust region radius, x is the position, fx = f(x). Caller ensures that the latter is feasible.

unconstrained_optimum(model)

Calculate the unconstrained optimum of the model, return its norm as the second value.

When the second value is infinite, the unconstrained optimum should not be used as this indicates a singular problem.