tpapp / MultistartOptimization.jl

Multistart optimization methods in Julia.
Other
64 stars 10 forks source link

allow for differentiable objective functions and local gradient-based methods #27

Closed tbeason closed 2 years ago

tbeason commented 2 years ago

Although Guvenen et. al. really stress the performance and usefulness of TikTak for nondifferentiable problems, there is really no reason to limit ourselves to that. In particular, I am hoping to use it for a differentiable (but possibly nonconvex) problem. Luckily, the modifications necessary to make it work were super easy! I found a function called applicable which lets us test whether a method exists (I did not know about this before). So the relevant change is

# if a method objective(x,grad) exists, use it. otherwise assume objective is not differentiable (which was the previous behavior)
opt.min_objective = applicable(objective,x,x) ? objective : nloptwrapper(objective)

with

function nloptwrapper(fn)  
    function f̃(x,grad)              # wrapper for NLopt
        @argcheck isempty(grad)     # ensure no derivatives are asked for
        return fn(x)
    end
    return f̃
end

All tests continue to pass for existing nonsmooth local methods, yet now also work when gradients are computed via ForwardDiff.jl.

I updated to add a test/Project.toml (the only way I actually know how to deal with test dependencies) and bumped the package version.

This is really mostly for the NLopt backend. I don't think any functionality was changed regarding other backends.

Hoping that you like this PR and will accept/merge/tag new release.

tbeason commented 2 years ago

Obviously it is more performant to use this information when you know it is reliable. See below for example. (the histograms didn't copy well)


julia> using MultistartOptimization, NLopt, ForwardDiff, BenchmarkTools

julia> function autodiff(fn)
               # adapted from here https://github.com/JuliaOpt/NLopt.jl/issues/128
               function f(x)
                   return fn(x)
               end

               function f(x,∇f)
                           if !(∇f == nothing) && (length(∇f) != 0)
                       ForwardDiff.gradient!(∇f,fn,x)
                   end

                   fx = fn(x)
                   return fx
               end
               return f
           end
autodiff (generic function with 1 method)

julia> function rosenbrock(x)
           x1 = x[1:(end - 1)]
           x2 = x[2:end]
           sum(@. 100 * abs2(x2 - abs2(x1)) + abs2(x1 - 1))
       end
rosenbrock (generic function with 1 method)

julia> P = MinimizationProblem(rosenbrock, -30 .* ones(10), 30 .* ones(10))
MinimizationProblem{typeof(rosenbrock), Vector{Float64}}(rosenbrock, [-30.0, -30.0, -30.0, -30.0, -30.0, -30.0, -30.0, -30.0, -30.0, -30.0], [30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0])

julia> multistart_method = TikTak(100)
TikTak(100, 10, 0.1, 0.995, 0.5)

julia> local_method = NLoptLocalMethod(NLopt.LN_BOBYQA)
NLoptLocalMethod{Set{Symbol}}(NLopt.LN_BOBYQA, 1.0e-8, 1.0e-8, 0, 0.0, Set([:SUCCESS, :XTOL_REACHED, :FTOL_REACHED, :MAXTIME_REACHED, :STOPVAL_REACHED, :MAXEVAL_REACHED]))

julia> @benchmark multistart_minimization($multistart_method, $local_method, $P)
BenchmarkTools.Trial: 94 samples with 1 evaluation.
 Range (min … max):  49.958 ms … 74.214 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     51.535 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   53.246 ms ±  4.871 ms  ┊ GC (mean ± σ):  0.25% ± 0.73%

  ▁▆█▃ ▂
  ████▄██▇▅▃▁▅▁▅▄▁▁▃▁▄▁▁▁▃▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▃▃ ▁
  50 ms           Histogram: frequency by time          74 ms <

 Memory estimate: 4.97 MiB, allocs estimate: 66952.

julia> adrosenbrock = autodiff(rosenbrock)
(::var"#f#1"{typeof(rosenbrock)}) (generic function with 2 methods)

julia> P = MinimizationProblem(adrosenbrock, -30 .* ones(10), 30 .* ones(10))
MinimizationProblem{var"#f#1"{typeof(rosenbrock)}, Vector{Float64}}(var"#f#1"{typeof(rosenbrock)}(rosenbrock), [-30.0, -30.0, -30.0, -30.0, -30.0, -30.0, -30.0, -30.0, -30.0, -30.0], [30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0])

julia> local_method = NLoptLocalMethod(NLopt.LD_LBFGS)
NLoptLocalMethod{Set{Symbol}}(NLopt.LD_LBFGS, 1.0e-8, 1.0e-8, 0, 0.0, Set([:SUCCESS, :XTOL_REACHED, :FTOL_REACHED, :MAXTIME_REACHED, :STOPVAL_REACHED, :MAXEVAL_REACHED]))

julia> @benchmark multistart_minimization($multistart_method, $local_method, $P)
BenchmarkTools.Trial: 289 samples with 1 evaluation.
 Range (min … max):  15.104 ms …  25.720 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     17.305 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   17.298 ms ± 934.376 μs  ┊ GC (mean ± σ):  0.43% ± 1.31%

                             █▇▄▁     ▂▄
  ▃▄▁▃▄▃▄▃▄▁▃▄▃▄▃▄▃▃▃▄▁▃▃▄▃▃▆████▆▆▆████▇▇▇▇▆▄▄▅▃▄▄▄▄▃▃▄▃▃▃▁▃▄ ▃
  15.1 ms         Histogram: frequency by time           19 ms <

 Memory estimate: 3.99 MiB, allocs estimate: 11668.

julia> multistart_minimization(multistart_method, local_method, P)
(value = 4.24418266037787e-21, location = [1.0000000000001186, 1.0000000000002918, 1.0000000000007943, 1.0000000000010416, 0.9999999999995424, 0.9999999999975635, 0.9999999999917186, 0.9999999999815047, 0.9999999999631132, 0.999999999925814], ret = :SUCCESS)
tbeason commented 2 years ago

@tpapp any chance you could have a quick look? Like I said, I don't existing functionality sees any changes, it just expandsthe set of possible local NLopt solvers.

tpapp commented 2 years ago

@tbeason: apologies for not responding, and thanks for the reminder. LGTM, just two minor suggestions, then I will merge.

tbeason commented 2 years ago

@tpapp no worries on delay. I pushed those changes so this should be good to go.

thanks again for implementing this algorithm

jonasmac16 commented 2 years ago

On a side note you might also want to look at GalacticOptim it hooks into MultistartOptimization and allows to pass any GalacticOptim linked optimiser via the generic objective interface. https://galacticoptim.sciml.ai/dev/optimization_packages/multistartoptimization/