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

Specific Domains Optimization (Mathematical)
1

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()))
neg_gradient_norm(coord) = -norm(Interpolations.gradient(interpolation, coord...))
  
xs = 14.0:24.0:110.0                                                                                                                                                                           
ys = 14.0:24.0:110.0                      
center_pt = [62.0, 62.0]           
result_optim = optimize(neg_gradient_norm, [xs[begin], ys[begin]], 
                                               [xs[end], ys[end]], 
                                               center_pt, 
                                               Fminbox(GradientDescent()); 
                                               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
 [2] gradient(::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}}}}, ::Float64, ::Float64) at /home/graham/.julia/packages/Interpolations/TBuvH/src/gridded/indexing.j
l:20
 [3] (::TravelingWaveSimulationsPlotting.var"#neg_gradient_norm#68"{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}}}}})(::Array{Float64,1}) at /home/graham/git/TravelingWaveSimulationsPlotting/src/space_reduction.jl:74
 [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]
2

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]   

neg_gradient_norm(coord) = -norm(Interpolations.gradient(interpolation, coord...))
interpolation = interpolate((xs, ys), A, Gridded(Linear()))

result_optim = optimize(neg_gradient_norm, [xs[begin], ys[begin]], 
                                            [xs[end], ys[end]], 
                                            center_pt, 
                                            Fminbox(GradientDescent()); 
                                            autodiff=:forward)
3

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()))

           function neg_gradient_norm(coord)
               return -norm(Interpolations.gradient(interpolation, coord...))
           end

           result_optim = optimize(
               neg_gradient_norm, 
               [xs[begin], ys[begin]], 
               [xs[end], ys[end]], 
               center_pt, 
               Fminbox(GradientDescent()); 
               # 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
    Algorithm:     Fminbox with Gradient Descent

 * 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