I have a nonlinear optimization problem that I’ve tried solving with scipy. For small scale it (the ‘TNC’ option in scipy minimize, specifically) works well, but it’s too slow when I get to thousands of variables. I’m learning Julia to see if I can find something better here.
Here is an example of a small scale version of the problem:
(Using ComponentArrays for convenience)
k_unr_out = [ 0., 3., 1., 5., 2., 4., 5., 0., 0., 2., 0., 4., 5.,
6., 1., 22., 0., 1., 19., 0., 0., 7., 0., 4., 0., 6.]
k_unr_in = [ 7., 4., 9., 3., 2., 3., 4., 11., 4., 4., 5., 4., 3.,
3., 3., 0., 2., 6., 0., 3., 0., 1., 8., 2., 3., 3.]
k_recip = [ 1., 8., 2., 8., 3., 7., 6., 1., 1., 7., 4., 5., 4.,
2., 9., 3., 2., 6., 6., 2., 25., 6., 4., 3., 3., 6.]
N = length(k_unr_out)
v0 = ComponentArray(x=repeat([0.1], N), y=repeat([0.1], N), z=repeat([0.1], N))
function neg_llhood(v, k_unr_out, k_unr_in, k_recip)
llhood = sum(k_unr_out .* log.(v.x))
llhood += sum(k_unr_in .* log.(v.y))
llhood += sum(k_recip .* log.(v.z))
xy = v.x .* v.y'
zz = v.z .* v.z'
Q = log.(1 .+ xy + xy' + zz)
llhood -= sum(tril(Q, -1))
-llhood
end
f = v -> neg_llhood(v, k_unr_out, k_unr_in, k_recip)
Since logarithms are taken for everything, my only constraints are x_i, y_i, z_i > 0 for all i. I’d like to use Optim to solve the above but I’m not sure exactly how to do that.
In scipy I was using JAX to do auto differentiation so I’m trying to do that here as well (I could compute the gradient analytically if I need to tho). For box constraints the docs for Optim give this example,
lower = [1.25, -2.1]
upper = [Inf, Inf]
initial_x = [2.0, 2.0]
inner_optimizer = GradientDescent()
results = optimize(f, g!, lower, upper, initial_x, Fminbox(inner_optimizer))
I’m confused by the docs. What is Fminbox doing? What is an inner optimizer? How do I specify that I want to use AD?
I tried
lower = repeat([floatmin()], 3*N)
upper = repeat([Inf], 3*N)
inner_optimizer = GradientDescent()
results = optimize(f, lower, upper, v0, Fminbox(inner_optimizer); autodiff=:forward)
And I get
MethodError: no method matching ldiv!(::ComponentVector{Float64}, ::Optim.InverseDiagonal, ::ComponentVector{Float64})
Closest candidates are:
...
I’m not sure if JuMP + a solver is preferred for my given problem. I gave it a try
model = Model(Ipopt.Optimizer)
@variable(model, x[1:N] >= floatmin(), start=0.1)
@variable(model, y[1:N] >= floatmin(), start=0.1)
@variable(model, z[1:N] >= floatmin(), start=0.1)
@NLobjective(model, Max,
sum(k_unr_out[i] * log(x[i]) + k_unr_in[i] * log(y[i]) + k_recip[i] * log(z[i]) for i in 1:N)
- sum(log(1 + x[i] * y[j] + x[j] * y[i] + z[i] * z[j]) for i in 1:N for j in i+1:N))
optimize!(model)
This causes my Jupyter kernel to crash.
I’m a bit lost in how to solve this problem most efficiently. Any help appreciated.