GALAHAD UGO package#
purpose#
The ugo package aims to find the global minimizer of a univariate
twice-continuously differentiable function \(f(x)\) of a single
variable over the finite interval \(x^l <= x <= x^u\). Function
and derivative values are provided via a subroutine call.
Second derivatives may be used to advantage if they are available.
See Section 4 of $GALAHAD/doc/ugo.pdf for additional details.
method#
The algorithm starts by splitting the interval \([x^l,x^u]\) into a specified number of subintervals \([x_i,x_{i+1}]\) of equal length, and evaluating \(f\) and its derivatives at each \(x_i\). A surrogate (approximating) lower bound function is constructed on each subinterval using the function and derivative values at each end, and an estimate of the first- and second-derivative Lipschitz constant. This surrogate is minimized, the true objective evaluated at the best predicted point, and the corresponding interval split again at this point. Any interval whose surrogate lower bound value exceeds an evaluated actual value is discarded. The method continues until only one interval of a maximum permitted width remains.
reference#
Many ingredients in the algorithm are based on the paper
D. Lera and Ya. D. Sergeyev, “Acceleration of univariate global optimization algorithms working with Lipschitz functions and Lipschitz first derivatives” SIAM J. Optimization 23(1), (2013) 508–529
but adapted to use second derivatives.
introduction to function calls#
To solve a given problem, functions from the ugo package must be called in the following order:
ugo_initialize - provide default control parameters and set up initial data structures
ugo_read_specfile (optional) - override control values by reading replacement values from a file
ugo_import - set up problem data structures and fixed values
ugo_reset_control (optional) - possibly change control parameters if a sequence of problems are being solved
solve the problem by calling one of
ugo_solve_direct - solve using function calls to evaluate function and derivative values, or
ugo_solve_reverse - solve returning to the calling program to obtain function and derivative values
ugo_information (optional) - recover information about the solution and solution process
ugo_terminate - deallocate data structures
See the examples section for illustrations of use.
parametric real type T#
Below, the symbol T refers to a parametric real type that may be Float32 (single precision), Float64 (double precision) or, if supported, Float128 (quadruple precision).
callable functions#
function ugo_initialize(T, data, control, status)
Set default control values and initialize private data
Parameters:
data |
holds private internal data |
control |
is a structure containing control information (see ugo_control_type) |
status |
is a scalar variable of type Int32, that gives the exit status from the package. Possible values are (currently):
|
function ugo_read_specfile(T, control, specfile)
Read the content of a specification file, and assign values associated with given keywords to the corresponding control parameters. An in-depth discussion of specification files is available, and a detailed list of keywords with associated default values is provided in $GALAHAD/src/ugo/UGO.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/ugo.pdf for a list of how these keywords relate to the components of the control structure.
Parameters:
control |
is a structure containing control information (see ugo_control_type) |
specfile |
is a one-dimensional array of type Vararg{Cchar} that must give the name of the specification file |
function ugo_import(T, control, data, status, x_l, x_u)
Import problem data into internal storage prior to solution.
Parameters:
control |
|
data |
holds private internal data |
status |
is a scalar variable of type Int32, that gives the exit status from the package. Possible values are:
|
x_l |
is a scalar variable of type T, that holds the value \(x^l\) of the lower bound on the optimization variable \(x\). |
x_u |
is a scalar variable of type T, that holds the value \(x^u\) of the upper bound on the optimization variable \(x\). |
function ugo_reset_control(T, control, data, status)
Reset control parameters after import if required.
Parameters:
control |
is a structure whose members provide control parameters for the remaining procedures (see ugo_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type Int32, that gives the exit status from the package. Possible values are:
|
function ugo_solve_direct(T, data, userdata, status, x, f, g, h, eval_fgh)
Find an approximation to the global minimizer of a given univariate function with a Lipschitz gradient in an interval.
This version is for the case where all function/derivative information is available by function calls.
Parameters:
data |
holds private internal data |
userdata |
is a structure that allows data to be passed into the function and derivative evaluation programs (see below). |
status |
is a scalar variable of type Int32, that gives the entry and exit status from the package. On initial entry, status must be set to 1. Possible exit values are:
|
x |
is a scalar variable of type T, that holds the value of the approximate global minimizer \(x\) after a successful (status = 0) call. |
f |
is a scalar variable of type T, that holds the the value of the objective function \(f(x)\) at the approximate global minimizer \(x\) after a successful (status = 0) call. |
g |
is a scalar variable of type T, that holds the the value of the gradient of the objective function \(f^{\prime}(x)\) at the approximate global minimizer \(x\) after a successful (status = 0) call. |
h |
is a scalar variable of type T, that holds the the value of the second derivative of the objective function \(f^{\prime\prime}(x)\) at the approximate global minimizer \(x\) after a successful (status = 0) call. |
eval_fgh |
is a user-provided function that must have the following signature: function eval_fgh(x, f, g, h, userdata) The value of the objective function \(f(x)\) and its first derivative \(f^{\prime}(x)\) evaluated at x= \(x\) must be assigned to f and g respectively, and the function return value set to 0. In addition, if control.second_derivatives_available has been set to true, when calling ugo_import, the user must also assign the value of the second derivative \(f^{\prime\prime}(x)\) in h; it need not be assigned otherwise. If the evaluation is impossible at x, return should be set to a nonzero value. |
function ugo_solve_reverse(T, data, status, eval_status, x, f, g, h)
Find an approximation to the global minimizer of a given univariate function with a Lipschitz gradient in an interval.
This version is for the case where function/derivative information is only available by returning to the calling procedure.
Parameters:
data |
holds private internal data |
status |
is a scalar variable of type Int32, that gives the entry and exit status from the package. On initial entry, status must be set to 1. Possible exit values are:
|
eval_status |
is a scalar variable of type Int32, that is used to indicate if objective function and its derivatives can be provided (see above). |
x |
is a scalar variable of type T, that holds the next value of \(x\) at which the user is required to evaluate the objective (and its derivatives) when status > 0, or the value of the approximate global minimizer when status = 0 |
f |
is a scalar variable of type T, that must be set by the user to hold the value of \(f(x)\) if required by status > 0 (see above), and will return the value of the approximate global minimum when status = 0 |
g |
is a scalar variable of type T, that must be set by the user to hold the value of \(f^{\prime}(x)\) if required by status > 0 (see above), and will return the value of the first derivative of \(f\) at the approximate global minimizer when status = 0 |
h |
is a scalar variable of type T, that must be set by the user to hold the value of \(f^{\prime\prime}(x)\) if required by status > 0 (see above), and will return the value of the second derivative of \(f\) at the approximate global minimizer when status = 0 |
function ugo_information(T, data, inform, status)
Provides output information
Parameters:
data |
holds private internal data |
inform |
is a structure containing output information (see ugo_inform_type) |
status |
is a scalar variable of type Int32, that gives the exit status from the package. Possible values are (currently):
|
function ugo_terminate(T, data, control, inform)
Deallocate all internal private storage
Parameters:
data |
holds private internal data |
control |
is a structure containing control information (see ugo_control_type) |
inform |
is a structure containing output information (see ugo_inform_type) |
available structures#
ugo_control_type structure#
struct ugo_control_type{T} error::Int32 out::Int32 print_level::Int32 start_print::Int32 stop_print::Int32 print_gap::Int32 maxit::Int32 initial_points::Int32 storage_increment::Int32 buffer::Int32 lipschitz_estimate_used::Int32 next_interval_selection::Int32 refine_with_newton::Int32 alive_unit::Int32 alive_file::NTuple{31,Cchar} stop_length::T small_g_for_newton::T small_g::T obj_sufficient::T global_lipschitz_constant::T reliability_parameter::T lipschitz_lower_bound::T cpu_time_limit::T clock_time_limit::T second_derivative_available::Bool space_critical::Bool deallocate_error_fatal::Bool prefix::NTuple{31,Cchar}
detailed documentation#
control derived type as a Julia structure
components#
Int32 error
error and warning diagnostics occur on stream error
Int32 out
general output occurs on stream out
Int32 print_level
the level of output required. Possible values are:
\(\leq\) 0 no output,
1 a one-line summary for every improvement
2 a summary of each iteration
\(\geq\) 3 increasingly verbose (debugging) output
Int32 start_print
any printing will start on this iteration
Int32 stop_print
any printing will stop on this iteration
Int32 print_gap
the number of iterations between printing
Int32 maxit
the maximum number of iterations allowed
Int32 initial_points
the number of initial (uniformly-spaced) evaluation points (<2 reset to 2)
Int32 storage_increment
incremenets of storage allocated (less that 1000 will be reset to 1000)
Int32 buffer
unit for any out-of-core writing when expanding arrays
Int32 lipschitz_estimate_used
what sort of Lipschitz constant estimate will be used:
1 = global contant provided
2 = global contant estimated
3 = local costants estimated
Int32 next_interval_selection
how is the next interval for examination chosen:
1 = traditional
2 = local_improvement
Int32 refine_with_newton
try refine_with_newton Newton steps from the vacinity of the global minimizer to try to improve the estimate
Int32 alive_unit
removal of the file alive_file from unit alive_unit terminates execution
NTuple{31,Cchar} alive_file
see alive_unit
T stop_length
overall convergence tolerances. The iteration will terminate when the step is less than .stop_length
T small_g_for_newton
if the absolute value of the gradient is smaller than small_g_for_newton, the next evaluation point may be at a Newton estimate of a local minimizer
T small_g
if the absolute value of the gradient at the end of the interval search is smaller than small_g, no Newton serach is necessary
T obj_sufficient
stop if the objective function is smaller than a specified value
T global_lipschitz_constant
the global Lipschitz constant for the gradient (-ve means unknown)
T reliability_parameter
the reliability parameter that is used to boost insufficiently large estimates of the Lipschitz constant (-ve means that default values will be chosen depending on whether second derivatives are provided or not)
T lipschitz_lower_bound
a lower bound on the Lipscitz constant for the gradient (not zero unless the function is constant)
T cpu_time_limit
the maximum CPU time allowed (-ve means infinite)
T clock_time_limit
the maximum elapsed clock time allowed (-ve means infinite)
Bool second_derivative_available
if .second_derivative_available is true, the user must provide them when requested. The package is generally more effective if second derivatives are available.
Bool space_critical
if .space_critical is true, every effort will be made to use as little space as possible. This may result in longer computation time
Bool deallocate_error_fatal
if .deallocate_error_fatal is true, any array/pointer deallocation error will terminate execution. Otherwise, computation will continue
NTuple{31,Cchar} prefix
all output lines will be prefixed by .prefix(2:LEN(TRIM(.prefix))-1) where .prefix contains the required string enclosed in quotes, e.g. “string” or ‘string’
ugo_time_type structure#
struct ugo_time_type{T} total::Float32 clock_total::T
detailed documentation#
time derived type as a Julia structure
components#
Float32 total
the total CPU time spent in the package
T clock_total
the total clock time spent in the package
ugo_inform_type structure#
struct ugo_inform_type{T} status::Int32 eval_status::Int32 alloc_status::Int32 bad_alloc::NTuple{81,Cchar} iter::Int32 f_eval::Int32 g_eval::Int32 h_eval::Int32 time::ugo_time_type{T}
detailed documentation#
inform derived type as a Julia structure
components#
Int32 status
return status. See UGO_solve for details
Int32 eval_status
evaluation status for reverse communication interface
Int32 alloc_status
the status of the last attempted allocation/deallocation
NTuple{81,Cchar) bad_alloc
the name of the array for which an allocation/deallocation error occurred
Int32 iter
the total number of iterations performed
Int32 f_eval
the total number of evaluations of the objective function
Int32 g_eval
the total number of evaluations of the gradient of the objective function
Int32 h_eval
the total number of evaluations of the Hessian of the objective function
ugo_time_type{T} ugo_time_type time
timings (see above)
example calls#
This is an example of how to use the package to minimize a univariate function; the code is available in $GALAHAD/src/ugo/Julia/test_ugo.jl .
# test_ugo.jl
# Simple code to test the Julia interface to UGO
using GALAHAD
using Test
using Printf
using Accessors
using Quadmath
function test_ugo(::Type{T}) where T
# Test problem objective
function objf(x::T)
a = 10.0
res = x * x * cos(a * x)
return res
end
# Test problem first derivative
function gradf(x::T)
a = 10.0
res = -a * x * x * sin(a * x) + 2.0 * x * cos(a * x)
return res
end
# Test problem second derivative
function hessf(x::T)
a = 10.0
res = -a * a * x * x * cos(a * x) - 4.0 * a * x * sin(a * x) + 2.0 * cos(a * x)
return res
end
# Derived types
data = Ref{Ptr{Cvoid}}()
control = Ref{ugo_control_type{T}}()
inform = Ref{ugo_inform_type{T}}()
# Initialize UGO
status = Ref{Cint}()
eval_status = Ref{Cint}()
ugo_initialize(T, data, control, status)
# Set user-defined control options
@reset control[].print_level = Cint(1)
# Test problem bounds
x_l = Ref{T}(-1.0)
x_u = Ref{T}(2.0)
# Test problem objective, gradient, Hessian values
x = Ref{T}(0.0)
f = Ref{T}(objf(x[]))
g = Ref{T}(gradf(x[]))
h = Ref{T}(hessf(x[]))
# import problem data
ugo_import(T, control, data, status, x_l, x_u)
# Set for initial entry
status[] = 1
# Solve the problem: min f(x), x_l ≤ x ≤ x_u
terminated = false
while !terminated
# Call UGO_solve
ugo_solve_reverse(T, data, status, eval_status, x, f, g, h)
# Evaluate f(x) and its derivatives as required
if (status[] ≥ 2) # need objective
f[] = objf(x[])
if (status[] ≥ 3) # need first derivative
g[] = gradf(x[])
if (status[] ≥ 4) # need second derivative
h[] = hessf(x[])
end
end
else # the solution has been found (or an error has occured)
terminated = true
end
end
# Record solution information
ugo_information(T, data, inform, status)
if inform[].status == 0
@printf("%i evaluations. Optimal objective value = %5.2f status = %1i\n",
inform[].f_eval, f[], inform[].status)
else
@printf("UGO_solve exit status = %1i\n", inform[].status)
end
# Delete internal workspace
ugo_terminate(T, data, control, inform)
return 0
end
@testset "UGO" begin
@test test_ugo(Float32) == 0
@test test_ugo(Float64) == 0
@test test_ugo(Float128) == 0
end