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

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

``````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 .+ (Δϕ / (Δxy * (size(fft_lr_2, 1)))),
fys .+ (Δϕ / (Δxy * (size(fft_lr_2, 2))))
) # non-shifted
fft_lr_2_sΔϕ = fft(
ifft(fft_lr_2) .* exp.(-2im * π * (Δϕ * fftshift(fxs * Δxy) + Δϕ * 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.

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

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

res = optimize.minimize(
xcorr,
alg.Δϕ_0;
bounds=((alg.Δϕ_0 .+ (-box_size, +box_size)),
(alg.Δϕ_0 .+ (-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

``````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.
I hope I provided everything that is necessary. Thank you very much for helping me untangle this!

files.zip

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:

``````x0 = [float.(alg.Δϕ_0)...]
fails to compute the derivatives of `hypot` and `fft`. So, there might be something along these lines.