Optimization, from MATLAB to Julia

Specific Domains Optimization (Mathematical)
,
1

Hi there!

I am fairly new to Julia, I was converted after seeing some awesome performance compared to MATLAB, which I currently use. I tried rewriting in Julia an optimization problem that I had in MATLAB, and everything is much faster (as I hoped) except for the actual optimization. What am I doing wrong? To be clear, the optimization works, it’s just very slow compared to MATLAB. I attach both.
Also, is there in Julia something analogous to “Run and time” of MATLAB? Something showing the code, line-by-line, with the time spent in each line?

Julia:

using Optimizion, OptimizationOptimJL, Plots, Dierckx, OptimizationMOI, Ipopt, ForwardDiff, ModelingToolkit, Symbolics

# functions
function action(b::Vector{Float64},I::Vector{Float64})
    hh=[I[2:1:end];0]-I
    pop!(hh)
    bx=I[1:length(hh)]+hh./2
    bs=Spline1D(I, b; k=1)(bx)
    sx=diff(b)./hh
    Lag=lagrangian(bx,bs,sx)
    Int=hh.*Lag
    S=sum(Int)./(2*pi)
    return S
end
function lagrangian(x::Vector{Float64},s::Vector{Float64},sx::Vector{Float64})
	L=sqrt.(sx.^2 + s.^4)
	return L
end
function interval(N::Int)
    I=range(0,(L/2)^3,length=N)
    I=I.^(1/3)
    I=[-I[end:-1:1];I]
    filter!(e->hash(e)≠hash(-0.0),I)
    return I
end

# definitions
rmax = 1e2
L = 0.3
I = interval(40)
b0 = ones(length(I))
S(b, I) = action(b, I)

# optimization
cons(res, b, I) = (res .= [b[1], b[end]])
optprob = OptimizationFunction(S, Optimization.AutoFiniteDiff(), cons = cons)
prob = OptimizationProblem(optprob, b0, I, lcons = [rmax, rmax], ucons = [rmax, rmax], lb = zeros(length(I)), ub = rmax.*ones(length(I)))
sol = solve(prob, IPNewton());

MATLAB:

% parameters
h = 1e-2;
rmax = 1e2;
L = 0.3;

% creates interval, denser near the extrema
I = linspace(0,(L/2)^3,40);
I=I.^(1/3);
I=unique([-flip(I),I]);

b0 = ones(numel(I),1);
S = @(b) action(I,b);

% impose constraints
AA = zeros(length(b0)); AA(1) = 1; AA(end) = 1;
bb = zeros(size(b0)); bb(1) = rmax; bb(end) = rmax;
lb=zeros(size(b0))-h;
ub=rmax*ones(size(b0))+h;

% optimize
options = optimoptions('fmincon','Display','iter-detailed',...
    'OptimalityTolerance',1e-5,'MaxIterations',1e4,'ConstraintTolerance',1e-16,'MaxFunctionEvaluations',1e6);
[bopt,Sopt,exitflag,output]=fmincon(S,b0,[],[],AA,bb,lb,ub,[],options);

%% functions
function S = action(I,b)
    hh=circshift(I,-1)-I;
    hh(end)=[];
    bx=I(1:numel(hh))+hh./2; 
    bs=interp1(I,b,bx); 
    sx=diff(b)'./hh;
    L=lagrangian(bx,bs,sx);
    Int= hh .* L;
    S=sum(Int)./(2*pi);
end
function L=lagrangian(~,s,sx)
    L = sqrt(sx.^2 + s.^4);
end
1 Like
2

I would start by checking the tolerance that the Julia code is using and make sure that is comparable to the MATLAB code. Alternatively, you could check the number of iterations being done on each side.

2 Likes
3

If you are using a completely different optimization algorithm in Matlab, maybe you are comparing apples and oranges? e.g. are you also using finite-difference derivatives in Matlab (which are pretty inefficient in 40 dimensions)?

(As @mkitti says, you should also look at the tolerances … it could also be that in Matlab you are converging to a much lower tolerance.)

PS. Your Julia code for action is also very Matlab-ish, allocating lots of unnecessary intermediate arrays to do a sequence of separate “vectorized” steps. This leaves a lot of performance on the table compared to simply writing a single loop that accumulates your objective function in a single pass over the array (with no array allocations), which looks easy to write at first glance at your code. But this is independent of how quickly (in iterations / function evaluations) your optimization algorithm converges.

PPS. Argument-type declarations like ::Vector{Float64} do nothing for the performance of your code, but they make it impossible to use certain automatic-differentiation tools. Ultimately, you should try very hard to get an analytical gradient of your objective function, either manually or using AD, since that enables much more efficient algorithms for local optimization.

3 Likes
4

If you’re using VSCode , you can use @profview in the REPL to profile. Else you can use something like ProfileCanvas or ProfileView. See Optimizing your code (modernjuliaworkflows.github.io) for more details.

3 Likes
5

I ran your posted Julia code after replacing the last line with

@time sol = solve(prob, IPNewton());

The result, starting from a fresh Julia session, was

19.150450 seconds (35.07 M allocations: 19.272 GiB, 8.71% gc time, 37.81% compilation time)

As you can see, a substantial fraction of the time is due to compilation, so I reran the last line with the result

11.569715 seconds (22.47 M allocations: 18.469 GiB, 12.24% gc time)

Better, but still very slow, and a huge amount of memory allocations. I then tried to make your code more Julian, in part by eliminating unnecessary allocations using generators, and other. Also, your call to Dierckx.Spline1D prevents auto-differentiation, because this library is written in Fortran and accepts only Float64 quantities. Fortunately, you are only asking for linear interpolation at the midpoints of the sampled values, i.e. their averaged values, which can be computed directly. With these changes the code becomes

using Optimization, OptimizationOptimJL, OptimizationMOI, Ipopt, ForwardDiff

# functions

lagrangian(s, sx) = sqrt(sx^2 + s^4) # For scalar arguments

function action(b, ivl)
    hhg = (ivl[i+1] - ivl[i] for i in 1:length(ivl)-1)
    bsg = ((b[i] + b[i+1])/2 for i in 1:length(b)-1)
    sxg = ((b[i+1] - b[i])/(ivl[i+1] - ivl[i]) for i in 1:length(ivl)-1)
    s = sum(lagrangian(bs, sx) * hh for (bs, sx, hh) in zip(bsg, sxg, hhg))
    s /= (2*pi)
    return s
end

function interval(n, L)
    ivl = Vector{Float64}(undef, 2n-1)
    for (i, x) in enumerate(range(0, (L/2)^3, length=n))
        y = x^(1/3)
        im1 = i - 1
        ivl[n - im1] = -y
        ivl[n + im1] = y
    end
    return ivl
end

# definitions

rmax = 1e2
L = 0.3
n = 40
ivl = interval(n, L)
b0 = ones(2n-1)
S(b, ivl) = action(b, ivl)
lb = zeros(2n-1)
ub = fill(rmax, 2n-1)

# optimization

cons(res, b, ivl) = (res .= [b[1], b[end]])
optprob = OptimizationFunction(S, Optimization.AutoFiniteDiff(), cons = cons)
prob = OptimizationProblem(optprob, b0, ivl; lcons = [rmax, rmax], ucons = [rmax, rmax], lb, ub)
@time sol = solve(prob, IPNewton());

and the output after the first call to solve is

7.857151 seconds (15.09 M allocations: 1.025 GiB, 6.60% gc time, 92.28% compilation time)

After re-submitting the final line in the same Julia session:

0.593737 seconds (2.43 M allocations: 221.809 MiB, 4.67% gc time)

There is a drastic reduction in both execution time and allocations.

Since the code is now auto-differentiable, I tried replacing Optimization.AutoFiniteDiff() with Optimization.AutoForwardDiff() and got this output on the first run:

20.042682 seconds (9.70 M allocations: 736.020 MiB, 0.94% gc time, 99.18% compilation time)

which is almost all compilation time. Rerunning it produced

0.166772 seconds (30.91 k allocations: 81.451 MiB, 11.55% gc time)

I don’t have much experience solving differentiable optimization problems, but I think the takeaway is that Optimization.jl is much faster when you need to solve the same or similar problems multiple times. The startup costs seem to be very high. I’m guessing that much of this startup cost is from the large wrapper that is Optimization.jl. Optimization.jl is convenient because it allows you to access many different solvers using the same API. Maybe once you’ve settled on an acceptable solver, it would incur less startup time to directly access that solver without going through the Optimization.jl wrapper. Perhaps others with more experience can comment on this.

12 Likes
6

This is a minor thing, but replacing x^(1/3) with cbrt(x) will likely be a small speedup since ^ needs to work for all powers, but cbrt knows exactly what operation you are doing.

4 Likes
7

Here’s the same code in JuMP:

julia> using JuMP

julia> import Ipopt

julia> lagrangian(s, sx) = sqrt(sx^2 + s^4)
lagrangian (generic function with 1 method)

julia> function action(b, I)
           hh = diff(I)
           bs = (b[1:end-1] .+ b[2:end]) ./ 2
           sx = diff(b) ./ hh
           return sum(hh .* lagrangian.(bs, sx)) / (2π)
       end
action (generic function with 1 method)

julia> function interval(N::Int, L)
           I = range(0, (L / 2)^3; length = N).^(1 / 3)
           return filter!(x -> x !== -0.0, vcat(-reverse(I), I))
       end
interval (generic function with 1 method)

julia> function main()
           I = interval(40, 0.4)
           N = length(I)
           model = Model(Ipopt.Optimizer)
           set_silent(model)
           @variable(model, 0 <= b[1:N] <= 100, start = 1)
           fix(b[1], 100; force = true)
           fix(b[end], 100; force = true)
           @objective(model, Min, action(b, I))
           optimize!(model)
           @assert is_solved_and_feasible(model)
           return value.(b)
       end
main (generic function with 1 method)

julia> @time main();

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

  1.063289 seconds (936.18 k allocations: 63.424 MiB, 2.34% gc time, 98.60% compilation time: 8% of which was recompilation)

julia> @time main();
  0.011167 seconds (26.37 k allocations: 1.655 MiB)
10 Likes
8

Awesome!

1 Like
9

For the Optimization.jl approach, I’d also recommend switching from forward-mode to reverse-mode AD, which will be much faster for taking gradients of high-dimensional problems (say, more than about 10 parameters). I tried replacing AutoForwardDiff() with AutoReverseDiff(true) and the post-compilation time for the solve call went from 0.5 s to 0.03 s.

Note the true option to AutoReverseDiff tells it to compile the “tape” of gradient instructions, which provides the best performance but can produce incorrect results if the code has branches, e.g. if or while statements.

7 Likes
10

When trying the AutoReverseDiff(true), it tells me that
LoadError: AutoReverseDiff does not currently support constraints.
How did you manage to run solve? Did you remove the constraints?

11

What versions of the packages are you using? If I replace ForwardDiff with ReverseDiff in @PeterSimon’s code above it runs fine for me in a fresh environment, even with AutoReverseDiff(true).

I’ve got the following package versions:

Status `/tmp/jl_IPYiyE/Project.toml`
  [f6369f11] ForwardDiff v0.10.36
  [b6b21f68] Ipopt v1.6.2
  [7f7a1694] Optimization v3.24.3
  [fd9f6733] OptimizationMOI v0.4.2
  [36348300] OptimizationOptimJL v0.2.3
  [37e2e3b7] ReverseDiff v1.15.1
1 Like