ComponentArrays.jl with NLopt.jl

Specific Domains Optimization (Mathematical)
,
1

These two packages seem to be incompatible. Without fully understanding the internals of NLopt.jl it seems that the vector of optimization variables ends up getting cast as an Array{Float64,1}. This basic snippet works:

using NLopt
function rosenbrock(x)
    (1 - x[1])^2 + 100*(x[2] - x[1]^2)^2
end

function rosenbrock(x, grad)
    if length(grad) > 0
        grad[1] = -400 * x[1] * (x[2] - x[1]^2) - 2 * (1 - x[1])
        grad[2] = 200 * (x[2] - x[1]^2)
    end
    rosenbrock(x)
end

opt = NLopt.Opt(:LD_SLSQP, 2);
opt.lower_bounds = [0., -0.5];
opt.upper_bounds = [1., 2.];
opt.min_objective = rosenbrock;

guess = [0.5, 0.];
(minf,minx,ret) = NLopt.optimize(opt, guess);
println("evals = $(opt.numevals)");
println("got $minf at $minx with constraints (returned $ret)");

But when I specify the rosenbrock function to work on ComponentArrays, it does not work (I get :FORCED_STOP):

using ComponentArrays
function rosenbrock(x)
    (1 - x.arg1)^2 + 100*(x.arg2 - x.arg1^2)^2
end
...
...
guess = ComponentArray(arg1=.5, arg2=0.)

Even when I pass this ComponentArray in as the initial vector, the NLopt.optimize function always returns Array{Float64,1}.

It would be great if it were possible to utilize ComponentArrays.jl within NLopt.jl. I’m solving optimization problems with 10’s or 100’s of optimization variables. Each specific problem I work on thus requires slicing the x vector with carefully constructed ranges in order to make sense of anything. With ComponentArrays.jl, I can store everything in a ComponentArray which appears as a flat array but can also be accessed via problem-specific field names.

2

You can write a function that changes the vector to a component vector and call this function first thing in your objective before proceeding with the rest of your objective definition.

4 Likes
3

I believe NLopt is a C/C++ library, hence they have no idea about fancy array types and only work with Array{Float64} and perhaps Float32. If you call C-code, you will probably have to do the conversion yourself like Mohamed suggests.

2 Likes
4

Doing it the way Mohamed suggests is a good way to go. The getaxes function from ComponentArrays will help here. This is what I usually do in situations like that:

guess = ComponentArray(arg1=.5, arg2=0.)
ax = getaxes(guess)

rosenbrock(x) = rosenbrock(ComponentArray(x, ax))
rosenbrock(x::ComponentArray) = (1 - x.arg1)^2 + 100*(x.arg2 - x.arg1^2)^2

Another option is to use a native Julia solver from Optim.jl, which will let ComponentArrays pass through just fine. I like using GalacticOptim.jl for a standard interface to all solvers, so here’s how I’d do it in GalacticOptim with an Optim solver:

using ComponentArrays, GalacticOptim

function rosenbrock(x, p)
    (p.a - x.arg1)^2 + p.b*(x.arg2 - x.arg1^2)^2
end

guess = ComponentArray(arg1=0.5, arg2=0)
lb = ComponentArray(arg1=0, arg2=-0.5)
ub = ComponentArray(arg1=1.0, arg2=2)

params = (a=1, b=100)

opt_fun = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
prob = OptimizationProblem(opt_fun, guess, params; lb=lb, ub=ub)

sol = solve(prob, Fminbox(BFGS()))
3 Likes