# Optim does not optimize in several minutes and it takes 10 seconds for scipy

**URL:** <https://discourse.julialang.org/t/optim-does-not-optimize-in-several-minutes-and-it-takes-10-seconds-for-scipy/101756>\
**Category:** Optimization (Mathematical)\
**Tags:** performance, optim\
**Created:** [July 18, 2023, 10:58pm UTC](https://discourse.julialang.org/t/optim-does-not-optimize-in-several-minutes-and-it-takes-10-seconds-for-scipy/101756 "2023-07-18T22:58:57Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![kunzaatko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kunzaatko/32/25153_2.png) [@kunzaatko](https://discourse.julialang.org/u/kunzaatko)\
**Post date:** [July 18, 2023, 10:58pm UTC](https://discourse.julialang.org/t/optim-does-not-optimize-in-several-minutes-and-it-takes-10-seconds-for-scipy/101756/1 "2023-07-18T22:58:57Z")

</div>

I am optimizing cross correlation of two arrays taken over there transfer function intersection. My function looks like this:

```julia
otf_lr_1 = otf(tf, fft_lr_1, Δxy) # does not change with Δϕ
fxs = fftfreq(size(fft_lr_2, 1), 1 / Δxy) .* ones(size(fft_lr_2, 2))'
fys = ones(size(fft_lr_2, 1)) .* fftfreq(size(fft_lr_2, 2), 1 / Δxy)'

function xcorr(Δϕ)
    otf_lr_2 = otf.(tf,
        fxs .+ (Δϕ[1] / (Δxy * (size(fft_lr_2, 1)))),
        fys .+ (Δϕ[2] / (Δxy * (size(fft_lr_2, 2))))
    ) # non-shifted
    fft_lr_2_sΔϕ = fft(
        ifft(fft_lr_2) .* exp.(-2im * π * (Δϕ[1] * fftshift(fxs * Δxy) + Δϕ[2] * fftshift(fys * Δxy)))
    ) # shifted by Δϕ (fourier shift theorem)
    Δϕr = hypot.(fxs .* Δxy, fys .* Δxy) ./ splat(hypot)(Δϕ ./ size(fft_lr_2_sΔϕ))
    intersection = dropdims(all(isone, cat(
                otf_lr_1 .> alg.MTF_min,
                otf_lr_2 .> alg.MTF_min,
                Δϕr .> alg.Δϕr_min,
                Δϕr .< 1 - alg.Δϕr_min;
                dims=3); dims=3); dims=3)
    fft_lr_1_i = zeros(eltype(fft_lr_1), size(fft_lr_1)) # at intersection
    fft_lr_1_i[intersection] = fft_lr_1[intersection] ./ otf_lr_1[intersection] # transfer normalized
                                                                                                                                                                          
    fft_lr_2_i = zeros(eltype(fft_lr_2_sΔϕ), size(fft_lr_2_sΔϕ)) # at intersection
    fft_lr_2_i[intersection] = fft_lr_2_sΔϕ[intersection] ./ otf_lr_1[intersection] # same OTF... Imitating the image formation where `lr_img_1` is at the center
                                                                                                                                                                          
    lr_1_i = ifft(fft_lr_1_i) # spatial from intersection
    lr_2_i = ifft(fft_lr_2_i)
                                                                                                                                                                          
    return -abs(sum(lr_2_i .* conj(lr_1_i)) / sum(abs.(lr_1_i) .^ 2))
end

```

When using `Optim.jl` and after waiting for about 10 minutes, the optimization did not finish running.

```julia
res = Optim.optimize(xcorr,
    [(alg.Δϕ_0 .- box_size)...],
    [(alg.Δϕ_0 .+ box_size)...],
    [float.(alg.Δϕ_0)...], Fminbox(LBFGS(; m=4))
)
                                                 
if Optim.converged(res)
    return Optim.minimizer(res)
else
    @error "Optimization did not converge"
end

```

Then after changing to `scipy`, I got the result (correct results) in about 10 seconds with

```julia
using PyCall
optimize = pyimport("scipy.optimize")

res = optimize.minimize(
    xcorr,
    alg.Δϕ_0;
    bounds=((alg.Δϕ_0[1] .+ (-box_size, +box_size)),
        (alg.Δϕ_0[2] .+ (-box_size, +box_size)))
)
                                                              
if res["success"]
    return res["x"] |> Tuple
else
    @error "Optimization did not converge: $(res["message"])"
end

```

It seems to me that the parameters of the optimization is the same in both cases.

Here are the packages that I am using and the serialized data that is input

```julia
using FFTW
using Unitful
using TransferFunctions
using OffsetArrays
using OffsetArrays: centered
using TransferFunctions
using TransferFunctions: otf_support
using Serialization
using Optim

fft_lr_1 = deserialize("file_one.ser")
fft_lr_2 = deserialize("file_two.ser")

struct Dummy
    Δϕ_0
    Δϕr_min
    MTF_min
end
alg = Dummy((-176,-24), 0.15, 0.05)
tf = IdealOTFwithCurvature(488u"nm", 1.4, 1.0, 0.3)
Δxy = 61u"nm"
box_size = Inf
otf_lr_1 = otf(tf, fft_lr_1, Δxy) # does not change with Δϕ
fxs = fftfreq(size(fft_lr_2, 1), 1 / Δxy) .* ones(size(fft_lr_2, 2))'
fys = ones(size(fft_lr_2, 1)) .* fftfreq(size(fft_lr_2, 2), 1 / Δxy)'

```

The `TransferFunctions` package is of my creation and is here on [GH](https://github.com/kunzaatko/TransferFunctions.jl).  
I hope I provided everything that is necessary. Thank you very much for helping me untangle this!

[files.zip](https://github.com/kunzaatko/TransferFunctions.jl/files/12089557/files.zip)

---

<div class="post-metadata">

**Author:** ![tmigot](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tmigot/32/23914_2.png) [@tmigot](https://discourse.julialang.org/u/tmigot)\
**Post date:** [July 31, 2023, 4:46pm UTC](https://discourse.julialang.org/t/optim-does-not-optimize-in-several-minutes-and-it-takes-10-seconds-for-scipy/101756/2 "2023-07-31T16:46:35Z")

</div>

Hi @kunzaatko !

Not sure why Optim goes on for such a long time, I tried using alternatives, but there seem to be some issues when computing gradient:

```julia
x0 = [float.(alg.Δϕ_0)...]
ForwardDiff.gradient(xcorr, x0)

```

fails to compute the derivatives of `hypot` and `fft`. So, there might be something along these lines.