Dealing with constant parameters
In many applications, there may be factors that are relevant to the function evaluations, but are fixed throughout the optimization. An obvious example is using data in a likelihood function, but it could also be parameters we wish to hold constant.
Consider a squared error loss function that depends on some data x and y, and parameters betas. As far as the solver is concerned, there should only be one input argument to the function we want to minimize, call it sqerror.
The problem is that we want to optimize a function sqerror that really depends on three inputs, and two of them are constant throughout the optimization procedure. To do this, we need to define the variables x and y
x = [1.0, 2.0, 3.0]
y = 1.0 + 2.0 * x + [-0.3, 0.3, -0.1]
We then simply define a function in three variables
function sqerror(betas, X, Y)
err = 0.0
for i in 1:length(X)
pred_i = betas[1] + betas[2] * X[i]
err += (Y[i] - pred_i)^2
end
return err
end
and then optimize the following anonymous function
res = optimize(b -> sqerror(b, x, y), [0.0, 0.0])
Alternatively, we can define a closure sqerror(betas) that is aware of the variables we just defined
function sqerror(betas)
err = 0.0
for i in 1:length(x)
pred_i = betas[1] + betas[2] * x[i]
err += (y[i] - pred_i)^2
end
return err
end
We can then optimize the sqerror function just like any other function
res = optimize(sqerror, [0.0, 0.0])
Avoid repeating computations
Say you are optimizing a function
f(x) = x[1]^2+x[2]^2
g!(storage, x) = copyto!(storage, [2x[1], 2x[2]])
In this situation, no calculations from f could be reused in g!. However, sometimes there is a substantial similarity between the objective function, and gradient, and some calculations can be reused.
To avoid repeating calculations, define functions fg! or fgh! that compute the objective function, the gradient and the Hessian (if needed) simultaneously. These functions internally can be written to avoid repeating common calculations.
For example, here we define a function fg! to compute the objective function and the gradient, as required:
function fg!(F,G,x)
# do common computations here
# ...
if G != nothing
# code to compute gradient here
# writing the result to the vector G
end
if F != nothing
# value = ... code to compute objective function
return value
end
end
Optim will only call this function with an argument G that is nothing (if the gradient is not required) or a Vector that should be filled (in-place) with the gradient. This flexibility is convenient for algorithms that only use the gradient in some iterations but not in others.
Now we call optimize with the following syntax:
Optim.optimize(Optim.only_fg!(fg!), [0., 0.], Optim.LBFGS())
Similarly, for a computation that requires the Hessian, we can write:
function fgh!(F,G,H,x)
G == nothing || # compute gradient and store in G
H == nothing || # compute Hessian and store in H
F == nothing || return f(x)
nothing
end
Optim.optimize(Optim.only_fgh!(fgh!), [0., 0.], Optim.Newton())
Provide gradients
As mentioned in the general introduction, passing analytical gradients can have an impact on performance. To show an example of this, consider the separable extension of the Rosenbrock function in dimension 5000, see SROSENBR in CUTEst.
Below, we use the gradients and objective functions from mastsif through CUTEst.jl. We only show the first five iterations of an attempt to minimize the function using Gradient Descent.
julia> @time optimize(f, initial_x, GradientDescent(),
Optim.Options(show_trace=true, iterations = 5))
Iter Function value Gradient norm
0 4.850000e+04 2.116000e+02
1 1.018734e+03 2.704951e+01
2 3.468449e+00 5.721261e-01
3 2.966899e+00 2.638790e-02
4 2.511859e+00 5.237768e-01
5 2.107853e+00 1.020287e-01
21.731129 seconds (1.61 M allocations: 63.434 MB, 0.03% gc time)
Results of Optimization Algorithm
* Algorithm: Gradient Descent
* Starting Point: [1.2,1.0, ...]
* Minimizer: [1.0287767703731154,1.058769439356144, ...]
* Minimum: 2.107853e+00
* Iterations: 5
* Convergence: false
* |x - x'| < 0.0: false
* |f(x) - f(x')| / |f(x)| < 0.0: false
* |g(x)| < 1.0e-08: false
* Reached Maximum Number of Iterations: true
* Objective Function Calls: 23
* Gradient Calls: 23
julia> @time optimize(f, g!, initial_x, GradientDescent(),
Optim.Options(show_trace=true, iterations = 5))
Iter Function value Gradient norm
0 4.850000e+04 2.116000e+02
1 1.018769e+03 2.704998e+01
2 3.468488e+00 5.721481e-01
3 2.966900e+00 2.638792e-02
4 2.511828e+00 5.237919e-01
5 2.107802e+00 1.020415e-01
0.009889 seconds (915 allocations: 270.266 KB)
Results of Optimization Algorithm
* Algorithm: Gradient Descent
* Starting Point: [1.2,1.0, ...]
* Minimizer: [1.0287763814102757,1.05876866832087, ...]
* Minimum: 2.107802e+00
* Iterations: 5
* Convergence: false
* |x - x'| < 0.0: false
* |f(x) - f(x')| / |f(x)| < 0.0: false
* |g(x)| < 1.0e-08: false
* Reached Maximum Number of Iterations: true
* Objective Function Calls: 23
* Gradient Calls: 23
The objective has obtained a value that is very similar between the two runs, but the run with the analytical gradient is way faster. It is possible that the finite differences code can be improved, but generally the optimization will be slowed down by all the function evaluations required to do the central finite differences calculations.
Separating time spent in Optim's code and user provided functions
Consider the Rosenbrock problem.
using Optim, OptimTestProblems
prob = UnconstrainedProblems.examples["Rosenbrock"];
Say we optimize this function, and look at the total run time of optimize using the Newton Trust Region method, and we are surprised that it takes a long time to run. We then wonder if time is spent in Optim's own code (solving the sub-problem for example) or in evaluating the objective, gradient or hessian that we provided. Then it can be very useful to use the TimerOutputs.jl package. This package allows us to run an over-all timer for optimize, and add individual timers for f, g!, and h!. Consider the example below, that is due to the author of the package (Kristoffer Carlsson).
using TimerOutputs
const to = TimerOutput()
f(x ) = @timeit to "f" prob.f(x)
g!(x, g) = @timeit to "g!" prob.g!(x, g)
h!(x, h) = @timeit to "h!" prob.h!(x, h)
begin
reset_timer!(to)
@timeit to "Trust Region" begin
res = Optim.optimize(f, g!, h!, prob.initial_x, NewtonTrustRegion())
end
show(to; allocations = false)
end
We see that the time is actually not spent in our provided functions, but most of the time is spent in the code for the trust region method.
Early stopping
Sometimes it might be of interest to stop the optimizer early. The simplest way to do this is to set the iterations keyword in Optim.Options to some number. This will prevent the iteration counter exceeding some limit, with the standard value being 1000. Alternatively, it is possible to put a soft limit on the run time of the optimization procedure by setting the time_limit keyword in the Optim.Options constructor.
using Optim, OptimTestProblems
problem = UnconstrainedProblems.examples["Rosenbrock"]
f = problem.f
initial_x = problem.initial_x
function slow(x)
sleep(0.1)
f(x)
end
start_time = time()
optimize(slow, zeros(2), NelderMead(), Optim.Options(time_limit = 3.0))
This will stop after about three seconds. If it is more important that we stop before the limit is reached, it is possible to use a callback with a simple model for predicting how much time will have passed when the next iteration is over. Consider the following code
using Optim, OptimTestProblems
problem = UnconstrainedProblems.examples["Rosenbrock"]
f = problem.f
initial_x = problem.initial_x
function very_slow(x)
sleep(.5)
f(x)
end
start_time = time()
time_to_setup = zeros(1)
function advanced_time_control(x)
println(" * Iteration: ", x.iteration)
so_far = time()-start_time
println(" * Time so far: ", so_far)
if x.iteration == 0
time_to_setup[:] = time()-start_time
else
expected_next_time = so_far + (time()-start_time-time_to_setup[1])/(x.iteration)
println(" * Next iteration ≈ ", expected_next_time)
println()
return expected_next_time < 13 ? false : true
end
println()
false
end
optimize(very_slow, zeros(2), NelderMead(), Optim.Options(callback = advanced_time_control))
It will try to predict the elapsed time after the next iteration is over, and stop now if it is expected to exceed the limit of 13 seconds. Running it, we get something like the following output
julia> optimize(very_slow, zeros(2), NelderMead(), Optim.Options(callback = advanced_time_control))
* Iteration: 0
* Time so far: 2.219298839569092
* Iteration: 1
* Time so far: 3.4006409645080566
* Next iteration ≈ 4.5429909229278564
* Iteration: 2
* Time so far: 4.403923988342285
* Next iteration ≈ 5.476739525794983
* Iteration: 3
* Time so far: 5.407265901565552
* Next iteration ≈ 6.4569235642751055
* Iteration: 4
* Time so far: 5.909044027328491
* Next iteration ≈ 6.821732044219971
* Iteration: 5
* Time so far: 6.912338972091675
* Next iteration ≈ 7.843148183822632
* Iteration: 6
* Time so far: 7.9156060218811035
* Next iteration ≈ 8.85849153995514
* Iteration: 7
* Time so far: 8.918903827667236
* Next iteration ≈ 9.870419979095459
* Iteration: 8
* Time so far: 9.922197818756104
* Next iteration ≈ 10.880185931921005
* Iteration: 9
* Time so far: 10.925468921661377
* Next iteration ≈ 11.888488478130764
* Iteration: 10
* Time so far: 11.92870283126831
* Next iteration ≈ 12.895747828483582
* Iteration: 11
* Time so far: 12.932114839553833
* Next iteration ≈ 13.902462200684981
Results of Optimization Algorithm
* Algorithm: Nelder-Mead
* Starting Point: [0.0,0.0]
* Minimizer: [0.23359374999999996,0.042187499999999996, ...]
* Minimum: 6.291677e-01
* Iterations: 11
* Convergence: false
* √(Σ(yᵢ-ȳ)²)/n < 1.0e-08: false
* Reached Maximum Number of Iterations: false
* Objective Function Calls: 24