# Optimization, from MATLAB to Julia

**URL:** <https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212>\
**Category:** Optimization (Mathematical)\
**Tags:** matlab, optimization\
**Created:** [March 5, 2024, 6:16pm UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212 "2024-03-05T18:16:55Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![cu2mauro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cu2mauro/32/207527_2.png) [@cu2mauro](https://discourse.julialang.org/u/cu2mauro)\
**Post date:** [March 5, 2024, 6:16pm UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212/1 "2024-03-05T18:16:55Z")

</div>

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:

```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:

```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

```

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [March 5, 2024, 6:31pm UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212/2 "2024-03-05T18:31:18Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 5, 2024, 6:52pm UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212/3 "2024-03-05T18:52:00Z")

</div>

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](https://julialang.org/blog/2017/01/moredots/#why_vectorized_code_is_not_as_fast_as_it_could_be) 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](https://docs.julialang.org/en/v1/manual/functions/#Argument-type-declarations) 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.

---

<div class="post-metadata">

**Author:** ![Zentrik](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zentrik/32/35409_2.png) [@Zentrik](https://discourse.julialang.org/u/Zentrik)\
**Post date:** [March 5, 2024, 7:26pm UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212/4 "2024-03-05T19:26:03Z")

</div>

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)](https://modernjuliaworkflows.github.io/optimizing/#profiling) for more details.

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [March 5, 2024, 10:10pm UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212/5 "2024-03-05T22:10:21Z")

</div>

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

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

```

The result, starting from a fresh Julia session, was

```julia
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

```julia
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](https://docs.julialang.org/en/v1/manual/arrays/#Generator-Expressions), 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

```julia
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

```julia
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:

```julia
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:

```julia
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

```julia
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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [March 5, 2024, 10:15pm UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212/6 "2024-03-05T22:15:33Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 5, 2024, 10:36pm UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212/7 "2024-03-05T22:36:09Z")

</div>

Here’s the same code in JuMP:

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [March 5, 2024, 10:50pm UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212/8 "2024-03-05T22:50:03Z")

</div>

> [@odow](#):
>
> Here’s the same code in JuMP:

Awesome!

---

<div class="post-metadata">

**Author:** ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)\
**Post date:** [March 5, 2024, 11:25pm UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212/9 "2024-03-05T23:25:09Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![cu2mauro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cu2mauro/32/207527_2.png) [@cu2mauro](https://discourse.julialang.org/u/cu2mauro)\
**Post date:** [April 13, 2024, 11:18am UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212/10 "2024-04-13T11:18:07Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)\
**Post date:** [April 13, 2024, 11:43am UTC](https://discourse.julialang.org/t/optimization-from-matlab-to-julia/111212/11 "2024-04-13T11:43:23Z")

</div>

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:

```julia
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

```