SimpleSolvers

SimpleSolvers.DEFAULT_WOLFE_c₁ — Constant
const DEFAULT_WOLFE_c₁

A constant $\epsilon$ on which a finite difference approximation of the derivative of the objective is computed. This is then used in the following stopping criterion:

\[\frac{f(\alpha) - f(\alpha_0)}{\epsilon} < \alpha\cdot{}f'(\alpha_0).\]

Extended help

source
SimpleSolvers.Backtracking — Type
Backtracking <: LinesearchMethod

The backtracking method.

Constructors

Backtracking()

Extended help

The backtracking algorithm starts by setting $y_0 \gets f(0)$ and $d_0 \gets \nabla_0f$.

The algorithm is executed by calling the functor of BacktrackingState.

The following is then repeated until the stopping criterion is satisfied or config.max_iterations (MAX_ITERATIONS by default) is reached:

if value!(obj, α) ≥ y₀ + ls.ϵ * α * d₀
    α *= ls.p
else
    break
end

The stopping criterion as an equation can be written as:

\[f(\alpha) < y_0 + \epsilon \alpha \nabla_0f = y_0 + \epsilon (\alpha - 0)\nabla_0f.\]

Note that if the stopping criterion is not reached, $\alpha$ is multiplied with $p$ and the process continues.

Sometimes the parameters $p$ and $\epsilon$ have different names such as $\tau$ and $c$.

source
SimpleSolvers.BacktrackingState — Type
BacktrackingState <: LinesearchState

Corresponding LinesearchState to Backtracking.

Keys

The keys are:

  • config::Options
  • α₀:
  • ϵ=$(DEFAULT_WOLFE_c₁): a default step size on whose basis we compute a finite difference approximation of the derivative of the objective. Also see DEFAULT_WOLFE_c₁.
  • p=$(DEFAULT_ARMIJO_p): a parameter with which $\alpha$ is decreased in every step until the stopping criterion is satisfied.

Functor

The functor is used the following way:

ls(obj, α = ls.α₀)

Implementation

The algorithm starts by setting:

\[x_0 \gets 0, y_0 \gets f(x_0), d_0 \gets f'(x_0), \alpha \gets \alpha_0,\]

where $f$ is the univariate objective (of type AbstractUnivariateObjective) and $\alpha_0$ is stored in ls. It then repeatedly does $\alpha \gets \alpha\cdot{}p$ until either (i) the maximum number of iterations is reached (the max_iterations keyword in Options) or (ii) the following holds:

\[ f(\alpha) < y_0 + \epsilon \cdot \alpha \cdot d_0,\]

where $\epsilon$ is stored in ls.

Info

The algorithm allocates an instance of SufficientDecreaseCondition by calling SufficientDecreaseCondition(ls.ϵ, x₀, y₀, d₀, one(α), obj), here we take the value one for the search direction $p$, this is because we already have the search direction encoded into the line search objective.

source
SimpleSolvers.Bisection — Type
Bisection <: LinesearchMethod

The bisection method.

Constructors

Bisection()

Extended help

The bisection algorithm starts with an interval and successively bisects it into smaller intervals until a root is found. See bisection.

source
SimpleSolvers.CurvatureCondition — Type
CurvatureCondition <: LinesearchCondition

The second of the Wolfe conditions [3]. The first one is the SufficientDecreaseCondition.

This encompasses the standard curvature condition and the strong curvature condition.

Constructor

CurvatureCondition(c, xₖ, gradₖ, pₖ, obj, grad; mode)

Here grad has to be a Gradient and obj an AbstractObjective. The other inputs are either arrays or numbers.

Implementation

For computational reasons CurvatureCondition also has a field gradₖ₊₁ in which the temporary new gradient is saved.

source
SimpleSolvers.GradientAutodiff — Type
GradientAutodiff <: Gradient

A struct that realizes Gradient by using ForwardDiff.

Keys

The struct stores:

  • F: a function that has to be differentiated.
  • ∇config: result of applying ForwardDiff.GradientConfig.

Constructors

GradientAutodiff(F, x::AbstractVector)
GradientAutodiff(F, nx::Integer)

Functor

The functor does:

grad(g, x) = ForwardDiff.gradient!(g, grad.F, x, grad.∇config)
source
SimpleSolvers.GradientFiniteDifferences — Type
GradientFiniteDifferences <: Gradient

A struct that realizes Gradient by using finite differences.

Keys

The struct stores:

  • F: a function that has to be differentiated.
  • ϵ: small constant on whose basis the finite differences are computed.
  • e: auxiliary vector used for computing finite differences. It's of the form $e_1 = \begin{bmatrix} 1 & 0 & \cdots & 0 \end{bmatrix}$.
  • tx: auxiliary vector used for computing finite differences. It stores the offset in the x vector.

Constructor(s)

GradientFiniteDifferences{T}(F, nx::Integer; ϵ)

By default for ϵ is DEFAULT_GRADIENT_ϵ.

Functor

The functor does:

for j in eachindex(x,g)
    ϵⱼ = grad.ϵ * x[j] + grad.ϵ
    fill!(grad.e, 0)
    grad.e[j] = 1
    grad.tx .= x .- ϵⱼ .* grad.e
    f1 = grad.F(grad.tx)
    grad.tx .= x .+ ϵⱼ .* grad.e
    f2 = grad.F(grad.tx)
    g[j] = (f2 - f1) / (2ϵⱼ)
end
source
SimpleSolvers.GradientFunction — Type
GradientFunction <: Gradient

A struct that realizes a Gradient by explicitly supplying a function.

Keys

The struct stores:

  • ∇F!: a function that can be applied in place.

Functor

The functor does:

grad(g, x) = grad.∇F!(g, x)
source
SimpleSolvers.Hessian — Type
Hessian

Abstract type. structs derived from this need an associated functor that computes the Hessian of a function (in-place).

Also see Gradient.

Implementation

When a custom Hessian is implemented, a functor is needed:

function (hessian::Hessian)(h::AbstractMatrix, x::AbstractVector) end

This functor can also be called with compute_hessian!.

Examples

Examples include:

These examples can also be called with e.g. Hessian(x; mode = :autodiff).

source
SimpleSolvers.HessianAutodiff — Type
HessianAutodiff <: Hessian

A struct that realizes Hessian by using ForwardDiff.

Keys

The struct stores:

  • F: a function that has to be differentiated.
  • H: a matrix in which the (updated) Hessian is stored.
  • Hconfig: result of applying ForwardDiff.HessianConfig.

Constructors

HessianAutodiff(F, x::AbstractVector)
HessianAutodiff(F, nx::Integer)

Functor

The functor does:

hes(g, x) = ForwardDiff.hessian!(hes.H, hes.F, x, grad.Hconfig)
source
SimpleSolvers.HessianFunction — Type
HessianFunction <: Hessian

A struct that realizes a Hessian by explicitly supplying a function.

Keys

The struct stores:

  • H!: a function that can be applied in place.

Functor

The functor does:

hes(H, x) = hes.H!(H, x)
source
SimpleSolvers.JacobianAutodiff — Type
JacobianAutodiff <: Jacobian

A struct that realizes Jacobian by using ForwardDiff.

Keys

The struct stores:

  • F: a function that has to be differentiated.
  • Jconfig: result of applying ForwardDiff.JacobianConfig.
  • ty: vector that is used for evaluating ForwardDiff.jacobian!

Constructors

JacobianAutodiff(F, y::AbstractVector)
JacobianAutodiff(F, nx::Integer)

Functor

The functor does:

jac(J, x) = ForwardDiff.jacobian!(J, jac.ty, x, grad.Jconfig)
source
SimpleSolvers.JacobianFiniteDifferences — Type
JacobianFiniteDifferences <: Jacobian

A struct that realizes Jacobian by using finite differences.

Keys

The struct stores:

  • F: a function that has to be differentiated.
  • ϵ: small constant on whose basis the finite differences are computed.
  • f1:
  • f2:
  • e1: auxiliary vector used for computing finite differences. It's of the form $e_1 = \begin{bmatrix} 1 & 0 & \cdots & 0 \end{bmatrix}$.
  • e2:
  • tx: auxiliary vector used for computing finite differences. It stores the offset in the x vector.

Constructor(s)

JacobianFiniteDifferences{T}(F, nx::Integer, ny::Integer; ϵ)

By default for ϵ is DEFAULT_JACOBIAN_ϵ.

Functor

The functor does:

for j in eachindex(x)
    ϵⱼ = jac.ϵ * x[j] + jac.ϵ
    fill!(jac.e, 0)
    jac.e[j] = 1
    jac.tx .= x .- ϵⱼ .* jac.e
    f(jac.f1, jac.tx)
    jac.tx .= x .+ ϵⱼ .* jac.e
    f(jac.f2, jac.tx)
    for i in eachindex(x)
        J[i,j] = (jac.f2[i] - jac.f1[i]) / (2ϵⱼ)
    end
end
source
SimpleSolvers.JacobianFunction — Type
JacobianFunction <: Jacobian

A struct that realizes a Jacobian by explicitly supplying a function.

Keys

The struct stores:

  • DF!: a function that can be applied in place.

Functor

The functor does:

jac(g, x) = jac.DF!(g, x)
source
SimpleSolvers.LUSolver — Type
struct LUSolver <: LinearSolver

A custom implementation of an LU solver.

Routines that use LUSolver include factorize! and ldiv!. In practice the LUSolver is used by calling its constructor together with ldiv! as shown in the Example section of this docstring.

Example

We use the LUSolver together with ldiv! to compute multiplication of a matrix inverse onto a vector (from the left):

A = [1. 2. 3.; 5. 7. 11.; 13. 17. 19.]

lu = LUSolver(A)
v = rand(3)
x = similar(v)

ldiv!(x, lu, v) ≈ inv(A) * v

# output

true

When calling LUSolver on an integer alone, a matrix with all zeros is allocated:

LUSolver{Float32}(2)

# output

LUSolver{Float32}(2, Float32[0.0 0.0; 0.0 0.0], [1, 2], [1, 2], 1)

Keys

  • n::Int
  • A::Matrix{T}
  • pivots::Vector{Int}
  • perms::Vector{Int}
  • info::Int
source
SimpleSolvers.LinearSolver — Type
LinearSolver <: AbstractSolver

A supertype that comprises e.g. LUSolver and LUSolverLAPACK.

Constructor

LinearSolver(x; linear_solver = :julia)

The convenience constructor allocates a specific struct derived from LinearSolver based on what is supplied to liner_solver. The default :julia calls the constructor for LUSolver. Another option would be :lapack which calls LUSolverLAPACK and uses the LinearAlgebra.BLAS package.

source
SimpleSolvers.LinesearchMethod — Type
LinesearchMethod

Examples include StaticState, Backtracking, Bisection and Quadratic. See these examples for specific information on linesearch algorithms.

Extended help

A LinesearchMethod always goes together with a LinesearchState and each of those LinesearchStates has a functor implemented:

ls(obj, α)

where obj is a AbstractUnivariateObjective and α is an initial step size. The output of this functor is then a final step size that is used for updating the parameters.

source
SimpleSolvers.LinesearchState — Type
LinesearchState

Abstract type.

Examples include StaticState, BacktrackingState, BisectionState and QuadraticState.

Implementation

A struct that is subtyped from LinesearchState needs to implement the functors:

ls(x; kwargs...)
ls(obj::AbstractUnivariateObjective, x; kwargs...)

Additionaly the following function needs to be extended:

LinesearchState(algorithm::LinesearchMethod; kwargs...)

Functors

The following functors are auxiliary helper functions:

ls(f::Callable; kwargs...) = ls(TemporaryUnivariateObjective(f, missing); kwargs...)
ls(f::Callable, x::Number; kwargs...) = ls(TemporaryUnivariateObjective(f, missing), x; kwargs...)
ls(f::Callable, g::Callable; kwargs...) = ls(TemporaryUnivariateObjective(f, g); kwargs...)
ls(f::Callable, g::Callable, x::Number; kwargs...) = ls(TemporaryUnivariateObjective(f, g), x; kwargs...)
source
SimpleSolvers.NewtonOptimizerCache — Type
NewtonOptimizerCache

Keys

  • x̄: the previous iterate,
  • x: current iterate (this stores the guess called by the functions generated with linesearch_objective),
  • δ: direction of optimization step (difference between x and x̄); this is obtained by multiplying rhs with the inverse of the Hessian,
  • g: gradient value (this stores the gradient associated with x called by the derivative part of linesearch_objective),
  • rhs: the right hand side used to compute the update.

To understand how these are used in practice see e.g. linesearch_objective.

source
SimpleSolvers.OptimizationAlgorithm — Type

An OptimizationAlgorithm is a data structure that is used to dispatch on different algorithms.

It needs to implement three methods,

initialize!(alg::OptimizationAlgorithm, ::AbstractVector)
update!(alg::OptimizationAlgorithm, ::AbstractVector)
solver_step!(::AbstractVector, alg::OptimizationAlgorithm)

that initialize and update the state of the algorithm and perform an actual optimization step.

Further the following convenience methods should be implemented,

objective(alg::OptimizationAlgorithm)
gradient(alg::OptimizationAlgorithm)
hessian(alg::OptimizationAlgorithm)
linesearch(alg::OptimizationAlgorithm)

which return the objective to optimize, its gradient and (approximate) Hessian as well as the linesearch algorithm used in conjunction with the optimization algorithm if any.

source
SimpleSolvers.Options — Type

Configurable options with defaults (values 0 and NaN indicate unlimited):

  • x_abstol = -Inf,
  • x_reltol = 4.440892098500626e-16,
  • x_suctol = 4.440892098500626e-16
  • f_abstol = 1.0e-50,
  • f_reltol = 4.440892098500626e-16,
  • f_suctol = 4.440892098500626e-16,
  • f_mindec = 0.0001,
  • g_restol = 1.4901161193847656e-8,
  • x_abstol_break = Inf,
  • x_reltol_break = Inf,
  • f_abstol_break = Inf,
  • f_reltol_break = Inf,
  • g_restol_break = Inf,
  • f_calls_limit = 0,
  • g_calls_limit = 0,
  • h_calls_limit = 0,
  • allow_f_increases = true,
  • min_iterations = 0,
  • max_iterations = 1000,
  • warn_iterations = 1000,
  • show_trace = false,
  • store_trace = false,
  • extended_trace = false,
  • show_every = 1,
  • verbosity = 1
source
SimpleSolvers.QuadraticState — Type
QuadraticState <: LinesearchState

Quadratic Polynomial line search.

Quadratic line search works by fitting a polynomial to a univariate objective (see AbstractUnivariateObjective) and then finding the minimum of that polynomial. Also compare this to BierlaireQuadraticState. The algorithm is taken from [1].

Keywords

source
SimpleSolvers.Static — Type
Static <: LinesearchMethod

The static method.

Constructors

Static(α)

Keys

Keys include: -α: equivalent to a step size. The default is 1.

Extended help

source
SimpleSolvers.SufficientDecreaseCondition — Type
SufficientDecreaseCondition <: LinesearchCondition

The condition that determines if $\alpha_k$ is big enough.

Constructor

SufficientDecreaseCondition(c₁, xₖ, fₖ, gradₖ, pₖ, obj)

Functors

sdc(xₖ₊₁, αₖ)
sdc(αₖ)

The second functor is shorthand for sdc(compute_new_iterate(sdc.xₖ, αₖ, sdc.pₖ), T(αₖ)), also see compute_new_iterate.

Extended help

We call the constant that pertains to the sufficient decrease condition $c$. This is typically called $c_1$ in the literature [3]. See DEFAULT_WOLFE_c₁ for the relevant constant

source
SimpleSolvers.UnivariateObjective — Type
UnivariateObjective <: AbstractUnivariateObjective

Keywords

It stores the following:

  • F: objective
  • D: derivative of objective
  • f: cache for function output
  • d: cache for derivative output
  • x_f: x used to evaluate F (stored in f)
  • x_d: x used to evaluate D (stored in d)
  • f_calls: number of times F has been called
  • d_calls: number of times D has been called

Constructor

There are several constructors, the most generic (besides the default one) is:

UnivariateObjective(F, D, x; f, d)

Where no keys are inferred, except x_f and x_d (via alloc_f and alloc_d). f_calls and d_calls are set to zero.

The most general constructor (i.e. the one the needs the least specification) is:

f(x::Number) = x ^ 2
UnivariateObjective(f, 1.)

# output

UnivariateObjective:

    f(x)              = NaN
    d(x)              = NaN
    x_f               = NaN
    x_d               = NaN
    number of f calls = 0
    number of d calls = 0

where ForwardDiff is used to generate the derivative of the (anonymous) function.

Functor

The functor calls value!.

source
GeometricBase.value — Method
value(obj::AbstractObjective, x)

Evaluates the objective value at x (i.e. computes obj.F(x)).

Examples

using SimpleSolvers

obj = UnivariateObjective(x::Number -> x^2, 1.)
value(obj, 2.)
obj.f_calls

# output

1

Note that the f_calls counter increased by one!

If value is called on obj (an AbstractObjective) without supplying x than the output of the last obj.F call is returned:

using SimpleSolvers

obj = UnivariateObjective(x::Number -> x^2, 1.)
value(obj)

# output

NaN

In this example this is NaN since the function hasn't been called yet.

source
SimpleSolvers.adjust_α — Method
adjust_α(ls, αₜ, α)

Check which conditions the new αₜ is in $[\sigma_0\alpha_0, \simga_1\alpha_0]$ and return the updated α accordingly (it is updated if it does not lie in the interval).

We first check the following:

\[ \alpha_t < \alpha_0\alpha,\]

where $\sigma_0$ is stored in ls (i.e. in an instance of QuadraticState). If this is not true we check:

\[ \alpha_t > \sigma_1\alpha,\]

where $\sigma_1$ is again stored in ls. If this second condition is also not true we simply return the unchanged $\alpha_t$. So if \alpha_t does not lie in the interval $(\sigma_0\alpha, \sigma_1\alpha)$ the interval is made bigger by either multiplying with $\sigma_0$ (default DEFAULT_ARMIJO_σ₀) or $\sigma_1$ (default DEFAULT_ARMIJO_σ₁).

source
SimpleSolvers.adjust_α — Method
adjust_α(αₜ, α)

Adjust αₜ based on the previous α. Also see adjust_α(::QuadraticState{T}, ::T, ::T) where {T}.

The check that $\alpha \in [\sigma_0\alpha_\mathrm{old}, \sigma_1\alpha_\mathrm{old}]$ should safeguard against stagnation in the iterates as well as checking that $\alpha$ decreases at least by a factor $\sigma_1$. The defaults for σ₀ and σ₁ are DEFAULT_ARMIJO_σ₀ and DEFAULT_ARMIJO_σ₁ respectively.

Implementation

Wee use defaults DEFAULT_ARMIJO_σ₀ and DEFAULT_ARMIJO_σ₁.

source
SimpleSolvers.alloc_j — Function
alloc_j(x, f)

Allocate NaNs of the size of the Jacobian of f (with respect to x).

This is used in combination with a MultivariateObjective.

Examples

x = rand(3)
fₓ = rand(2)
alloc_j(x, fₓ)

# output

2×3 Matrix{Float64}:
 NaN  NaN  NaN
 NaN  NaN  NaN
source
SimpleSolvers.bisection — Method
bisection(f, xmin, xmax; config)

Perform bisection of f in the interval [xmin, xmax] with Options config.

The algorithm is repeated until a root is found (up to tolerance config.f_abstol which is F_ABSTOL by default).

implementation

When calling bisection it first checks if $x_\mathrm{min} < x_\mathrm{max}$ and else flips the two entries.

Extended help

The bisection algorithm divides an interval into equal halves until a root is found (up to a desired accuracy).

We first initialize:

\[\begin{aligned} x_0 \gets & x_\mathrm{min}, x_1 \gets & x_\mathrm{max}, \end{aligned}\]

and then repeat:

\[\begin{aligned} & x \gets \frac{x_0 + x_1}{2}, \\ & \text{if $f(x_0)f(x) > 0$} \\ & \qquad x_0 \gets x, \\ & \text{else} \\ & \qquad x_1 \gets x, \\ & \text{end} \end{aligned}\]

So the algorithm checks in each step where the sign change occurred and moves the $x_0$ or $x_1$ accordingly. The loop is terminated (and errors) if config.max_iterations is reached (see MAX_ITERATIONS and the Options struct).

Warning

The obvious danger with using bisections is that the supplied interval can have multiple roots (or no roots). One should be careful to avoid this when fixing the interval.

source
SimpleSolvers.bracket_minimum — Method
bracket_minimum(f, x)

Move a bracket successively in the search direction (starting at x) and increase its size until a local minimum of f is found. This is used for performing Bisections when only one x is given (and not an entire interval). This bracketing algorithm is taken from [4]. Also compare it to bracket_minimum_with_fixed_point.

Keyword arguments

Extended help

For bracketing we need two constants $s$ and $k$ (see DEFAULT_BRACKETING_s and DEFAULT_BRACKETING_k).

Before we start the algorithm we initialize it, i.e. we check that we indeed have a descent direction:

\[\begin{aligned} & a \gets x, \\ & b \gets a + s, \\ & \mathrm{if} \quad f(b) > f(a)\\ & \qquad\text{Flip $a$ and $b$ and set $s\gets-s$.}\\ & \mathrm{end} \end{aligned}\]

The algorithm then successively computes:

\[c \gets b + s,\]

and then checks whether $f(c) > f(b)$. If this is true it returns $(a, c)$ or $(c, a)$, depending on whether $a<c$ or $c<a$ respectively. If this is not satisfied $a,$ $b$ and $s$ are updated:

\[\begin{aligned} a \gets & b, \\ b \gets & c, \\ s \gets & sk, \end{aligned}\]

and the algorithm is continued. If we have not found a sign chance after $n_\mathrm{max}$ iterations (see DEFAULT_BRACKETING_nmax) the algorithm is terminated and returns an error. The interval that is returned by bracket_minimum is then typically used as a starting point for bisection.

Info

The function bracket_root is equivalent to bracket_minimum with the only difference that the criterion we check for is:

\[f(c)f(b) < 0,\]

i.e. that a sign change in the function occurs.

See bracket_root.

source
SimpleSolvers.bracket_minimum_with_fixed_point — Method
bracket_minimum_with_fixed_point(f, x)

Find a bracket while keeping the left side (i.e. x) fixed. The algorithm is similar to bracket_minimum (also based on DEFAULT_BRACKETING_s and DEFAULT_BRACKETING_k) with the difference that for the latter the left side is also moving.

The function bracket_minimum_with_fixed_point is used as a starting point for Quadratic (taken from [1]), as the minimum of the polynomial approximation is:

\[p_2 = \frac{f(b) - f(a) - f'(0)b}{b^2},\]

where $b = \mathtt{bracket\_minimum\_with\_fixed\_point}(a)$. We check that $f(b) > f(a)$ in order to ensure that the curvature of the polynomial (i.e. $p_2$ is positive) and we have a minimum.

source
SimpleSolvers.check_gradient — Method
check_gradient(g)

Check norm, maximum value and minimum value of a vector.

Examples

using SimpleSolvers

g = [1., 1., 1., 2., 0.9, 3.]
SimpleSolvers.check_gradient(g; digits=3)

# output

norm(Gradient):               4.1
minimum(|Gradient|):          0.9
maximum(|Gradient|):          3.0
source
SimpleSolvers.check_hessian — Method
check_hessian(H)

Check the condition number, determinant, max and min value of the Hessian H.

using SimpleSolvers

H = [1. √2.; √2. 3.]
SimpleSolvers.check_hessian(H)

# output

Condition Number of Hessian: 13.9282
Determinant of Hessian:      1.0
minimum(|Hessian|):          1.0
maximum(|Hessian|):          3.0
source
SimpleSolvers.check_jacobian — Method
check_jacobian(J)

Check the condition number, determinant, max and min value of the Jacobian J.

using SimpleSolvers

J = [1. √2.; √2. 3.]
SimpleSolvers.check_jacobian(J)

# output

Condition Number of Jacobian: 13.9282
Determinant of Jacobian:      1.0
minimum(|Jacobian|):          1.0
maximum(|Jacobian|):          3.0
source
SimpleSolvers.compute_new_iterate — Method
compute_new_iterate(xₖ, αₖ, pₖ)

Compute xₖ₊₁ based on a direction pₖ and a step length αₖ.

Extended help

In the case of vector spaces this function simply does:

xₖ + αₖ * pₖ

For manifolds we instead perform a retraction [5].

source
SimpleSolvers.factorize! — Method
factorize!(lu, A)

Factorize the matrix A and store the result in lu.A.

Examples

A = [1. 2. 3.; 5. 7. 11.; 13. 17. 19.]
lu = LUSolver{Float64}(3, similar(A), zeros(Int, 3), zeros(Int, 3), 0)
factorize!(lu, A)
lu.A

# output

3×3 Matrix{Float64}:
 13.0        17.0       19.0
  0.0769231   0.692308   1.53846
  0.384615    0.666667   2.66667

Here lu.A stores the factorized result. If we want to save this factorized matrix in the same A to save memory we can write:

A = [1. 2. 3.; 5. 7. 11.; 13. 17. 19.]
lu = LUSolver{Float64}(3, A, zeros(Int, 3), zeros(Int, 3), 0)
factorize!(lu, A)
A

# output

3×3 Matrix{Float64}:
 13.0        17.0       19.0
  0.0769231   0.692308   1.53846
  0.384615    0.666667   2.66667
source
SimpleSolvers.increase_iteration_number! — Method
increase_iteration_number!(status)

Increase iteration number of status.

Examples

status = NonlinearSolverStatus{Float64}(5)
increase_iteration_number!(status)
status.i

# output

1
source
SimpleSolvers.linearsolver — Method
linearsolver(solver)

Return the linear part (i.e. a LinearSolver) of an AbstractNewtonSolver.

Examples

x = rand(3)
y = rand(3)
F(x) = tanh.(x)
s = NewtonSolver(x, y; F = F)
linearsolver(s)

# output

LUSolver{Float64}(3, [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0], [1, 2, 3], [1, 2, 3], 1)
source
SimpleSolvers.linesearch_objective — Method
linesearch_objective(objective, cache)

Create TemporaryUnivariateObjective for linesearch algorithm. The variable on which this objective depends is $\alpha$.

Example

x = [1, 0., 0.]
f = x -> sum(x .^ 3 / 6 + x .^ 2 / 2)
obj = MultivariateObjective(f, x)
gradient!(obj, x)
value!(obj, x)
cache = NewtonOptimizerCache(x)
hess = Hessian(obj, x; mode = :autodiff)
update!(hess, x)
update!(cache, x, obj.g, hess)
x₂ = [.9, 0., 0.]
gradient!(obj, x₂)
value!(obj, x₂)
update!(hess, x₂)
update!(cache, x₂, obj.g, hess)
ls_obj = linesearch_objective(obj, cache)
α = .1
(ls_obj.F(α), ls_obj.D(α))

# output

(0.4412947468016475, -0.8083161485821551)

In the example above we have to apply update! twice on the instance of NewtonOptimizerCache because it needs to store the current and the previous iterate.

Implementation

Calling the function and derivative stored in the TemporaryUnivariateObjective created with linesearch_objective does not allocate a new array, but uses the one stored in the instance of NewtonOptimizerCache.

source
SimpleSolvers.linesearch_objective — Method
linesearch_objective(objective!, jacobian!, cache)

Make a line search objective for a Newton solver (the cache here is an instance of NewtonSolverCache).

Implementation

Producing a single-valued output

Different from the linesearch_objective for NewtonOptimizerCaches, we apply l2norm to the output of objective!. This is because the solver operates on an objective with multiple outputs from which we have to find roots, whereas an optimizer operates on an objective with a single output of which we should find a minimum.

Also see linesearch_objective(::MultivariateObjective{T}, ::NewtonOptimizerCache{T}) where {T}.

source
SimpleSolvers.meets_stopping_criteria — Method
meets_stopping_criteria(status, config)

Check if the optimizer has converged.

Implementation

meets_stopping_criteria first calls assess_convergence! and then checks if one of the following is true:

  • converged (the output of assess_convergence!) is true and status.i $\geq$ config.min_iterations,
  • if config.allow_f_increases is false: status.f_increased is true,
  • status.i $\geq$ config.max_iterations,
  • status.rxₐ $>$ config.x_abstol_break
  • status.rxᵣ $>$ config.x_reltol_break
  • status.rfₐ $>$ config.f_abstol_break
  • status.rfᵣ $>$ config.f_reltol_break
  • status.rg $>$ config.g_restol_break
  • status.x_isnan
  • status.f_isnan
  • status.g_isnan
source
SimpleSolvers.print_status — Method
print_status(status, config)

Print the solver staus if:

  1. The following three are satisfied: (i) config.verbosity $\geq1$ (ii) assess_convergence!(status, config) is false (iii) status.i > config.max_iterations
  2. config.verbosity > 1.
source
SimpleSolvers.residual! — Method
residual!(status, x, x̄, f, f̄, g, ḡ)

Compute the residual based on previous iterates (x̄, f̄, ḡ) and current iterates (x, f, g).

source
SimpleSolvers.triple_point_finder — Method
triple_point_finder(f, x)

Find three points a > b > c s.t. f(a) > f(b) and f(c) > f(b). This is used for performing a quadratic line search (see QuadraticState).

Implementation

For δ we take DEFAULT_BRACKETING_s as default. For nmax we take [DEFAULTBRACKETINGnmax`](@ref) as default.

Extended help

The algorithm is taken from [2, Chapter 11.2.1].

source
SimpleSolvers.update! — Method
update!(cache::NewtonOptimizerCache, x, g, hes)

Update an instance of NewtonOptimizerCache based on x.

This sets:

\[\bar{x}^\mathtt{cache} \gets x, x^\mathtt{cache} \gets x, g^\mathtt{cache} \gets g, \mathrm{rhs}^\mathtt{cache} \gets -g, \delta^\mathtt{cache} \gets H^{-1}\mathrm{rhs}^\mathtt{cache},\]

where we wrote $H$ for the Hessian (i.e. the input argument hes).

Also see update!(::NewtonSolverCache, ::AbstractVector).

Implementation

The multiplication by the inverse of $H$ is done with LinearAlgebra.ldiv!.

source
SimpleSolvers.update! — Method
update!(status, x, f)

Update status based on x and the function f.

The new x and x̄ stored in status are used to compute δ. The new f and f̄ stored in status are used to compute γ.

source
SimpleSolvers.value!! — Method
value!!(obj::AbstractObjective, x)

Set obj.x_f to x and obj.f to value(obj, x) and return value(obj).

source