# Bounded gradient descent on bounded interpolation searches [NaN, NaN]

I’m confused how`[NaN, NaN]` could ever come up, and am at a loss as to how go about debugging. The same code works for several other example data. It seems like the issue is that there are so many zeros but I would think that would lead to a `[0, 0]` gradient… I also feel like the answer is staring me in the face and I’ve been looking at it too long. Can anyone offer any input?

EDIT: copy-pasteable script in second post

I have data that looks like this:

``````data = 5×5 AxisArray{Float64,2}
• dim_1 - 14.0:24.0:110.0
• dim_2 - 14.0:24.0:110.0
14.0   38.0   62.0   86.0   110.0
14.0    0.2   0.16    0.0    0.0     0.0
38.0   0.24   0.16    0.0    0.0     0.0
62.0   0.08    0.0    0.0    0.0     0.0
86.0    0.0    0.0    0.0    0.0     0.0
110.0    0.0    0.0    0.0    0.0     0.0
``````

and I’m trying to find where the gradient’s norm is maximized. So I use the following code:

``````function Interpolations.interpolate(data::AbstractAxisArray, args...; kwargs...)
axs = axes_keys(data)
interpolate(axs, data, args...; kwargs...)
end
interpolation = interpolate(data, Gridded(Linear()))

xs = 14.0:24.0:110.0
ys = 14.0:24.0:110.0
center_pt = [62.0, 62.0]
[xs[end], ys[end]],
center_pt,
autodiff=:forward)
``````

where

``````map(neg_gradient_norm, Iterators.product(xs, ys)) = [
-0.002357022603955158 -0.001666666666666667 -0.006666666666666666 -0.0 -0.0;
-0.0037267799624996476 -0.0033333333333333322 -0.006666666666666666 -0.0 -0.0;
-0.007453559924999298 -0.007453559924999299 -0.0 -0.0 -0.0;
-0.003333333333333333 -0.0 -0.0 -0.0 -0.0; -0.0 -0.0 -0.0 -0.0 -0.0]
``````

but I get the error:

``````ERROR: LoadError: BoundsError: attempt to access 5×5 interpolate((14.0:24.0:110.0,14.0:24.0:110.0), ::Array{Float64,2}, Gridded(Interpolations.Linear())) with element type Float64 at index [NaN, NaN]
Stacktrace:
[1] throw_boundserror(::Interpolations.GriddedInterpolation{Float64,2,Float64,Interpolations.Gridded{Interpolations.Linear},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::Tuple{Float64,Float64}) at ./abstractarray.jl:541
{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::Float64, ::Float64) at /home/graham/.julia/packages/Interpolations/TBuvH/src/gridded/indexing.j
l:20
[4] optimize(::NLSolversBase.OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Optim.Fminbox{Optim.GradientDescent{
LineSearches.InitialPrevious{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Nothing,Optim.var"#13#15"},Float64,Optim.var"#43#45"}, ::Optim.Options{Float64,Nothing}) at /home/graham/.julia/packages/Optim/SFpsz/src/multivariate/solvers/constrained/fminbox.jl:318
[5] #optimize#59 at /home/graham/.julia/packages/Optim/SFpsz/src/multivariate/solvers/constrained/fminbox.jl:176 [inlined]
``````

I haven’t figured out any solution, but I’ve turned it into a copy-pasteable script that doesn’t pirate interpolate. Any input super appreciated.

``````using Interpolations, Optim, LinearAlgebra

A = [0.2   0.16    0.0    0.0     0.0;
0.24   0.16    0.0    0.0     0.0;
0.08    0.0    0.0    0.0     0.0;
0.0    0.0    0.0    0.0     0.0;
0.0    0.0    0.0    0.0     0.0]

xs = 14.0:24.0:110.0
ys = 14.0:24.0:110.0
center_pt = [62.0, 62.0]

interpolation = interpolate((xs, ys), A, Gridded(Linear()))

[xs[end], ys[end]],
center_pt,
autodiff=:forward)
``````

I guess it’s an issue with ForwardDiff through the interpolation. Just remove `autodiff=:forward`.

``````julia> using Interpolations, Optim, LinearAlgebra

julia> function main()
A = [0.2   0.16    0.0    0.0     0.0;
0.24   0.16    0.0    0.0     0.0;
0.08    0.0    0.0    0.0     0.0;
0.0    0.0    0.0    0.0     0.0;
0.0    0.0    0.0    0.0     0.0]

xs = 14.0:24.0:110.0
ys = 14.0:24.0:110.0
center_pt = [62.0, 62.0]

interpolation = interpolate((xs, ys), A, Gridded(Linear()))

end

result_optim = optimize(
[xs[begin], ys[begin]],
[xs[end], ys[end]],
center_pt,
# autodiff = :forward
)
return result_optim
end
main (generic function with 1 method)

julia>

julia> main()
* Status: success

* Candidate solution
Final objective value:     -7.453560e-03

* Found with

* Convergence measures
|x - x'|               = 0.00e+00 ≤ 0.0e+00
|x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
|g(x)|                 = 0.00e+00 ≤ 1.0e-08

* Work counters
Seconds run:   0  (vs limit Inf)
Iterations:    2
f(x) calls:    184
∇f(x) calls:   184
``````