How to find location of Minimum in Interpolation function?

[CosmoProf]

Hi, when using an Interpolation object to approximate a function with values known over a few dozen discrete data points, how can I find where the (local) minimum of the function is? (Assuming that the minimum is not actually on one of the discrete data point locations, but somewhere in between two of them.)

Also, each of the commands below seem to produce the exact same thing (a connect-the-dots result with straight line segments between points); am I missing something? Thanks for any info!

plot(interpolate(data,BSpline(Constant())))
plot(interpolate(data,BSpline(Linear())))
plot(interpolate(data, BSpline(Quadratic(Reflect(OnCell())))))
plot(interpolate(data, BSpline(Cubic(Line(OnGrid())))))
(…etc…)

[anon92994695]

What you’re asking is basically an optimization question.

So when we have a discrete set or collection of numbers, ie: [1,2,3,4], we can instruct our computer to find the minimum value or index quite easily(min or argmin).

But when we have a “function”. We need to use numerical optimization methods more often then not. So there’s an awesome Julia package called Optim.jl(Optim.jl). Highly reccommend it!

Here’s an example…

using Interpolations, Optim
t= 1:0.1:2pi
y = sin.(t)
interp_cubic = CubicSplineInterpolation(t, y)
bounds = [1.0, (2pi) - 0.1]
t_interp = first(bounds):0.01:last(bounds) 
scatter([t], [y], label = "raw")
plot!([t_interp], [interp_cubic(t_interp)], label = "Spline")

opt = optimize(x -> interp_cubic(x), first(bounds), last(bounds))
t_min = Optim.minimizer(opt)
scatter!( [t_min], [interp_cubic(t_min)], color = :red,label = "estimated min")

Hope this helps.

[CosmoProf]

Hi ckneale, thanks a lot for your help & sample code.

I’ve been implementing it, and unfortunately the optimizer’s results are a bit poor & inconsistent.

For example, where y = y(x) is a simple function of x, and trying to minimize the distance z where z^2 = x^2 + y^2, I use cubic interpolation of z to find where z is smallest. x and y are (artificially) sampled as discrete data points.

It finds xMin, and zMin(xMin). Then using the distance formula, y = sqrt(z^2 - x^2), at xMin, I get y = 1.173. BUT, doing a cubic interpolation of the y data points themselves (and using the exact y(x)), yields the value 1.159 at xMin. That’s a pretty big difference!

Also, using Mathematica to find the exact minimum, gets xMin = 0.547 instead of Julia’s xMin = 0.554.

I understand some inexactness, but the inconsistency in y(xMin) is alarming. Is it supposed to be this off, or am I doing something wrong? I’ve included a copy of my sample code below… Thanks again for any feedback!

using Plots, Interpolations, Optim

k = 0.47
xSample = -12.0:1.0:12.0
y = Array{Float64,1}(undef, 25); z = Array{Float64,1}(undef, 25)
for i in 1:25
    y[i] = (1.42-(k*xSample[i]))
    z[i] = sqrt(xSample[i]^2+y[i]^2)
end
y_Interp_cubic = CubicSplineInterpolation(tSample, y)
z_Interp_cubic = CubicSplineInterpolation(tSample, z)
opt = optimize(Stuff -> z_Interp_cubic(Stuff), -12.0, 12.0)
xFORzMin = Optim.minimizer(opt)
yATmin = y_Interp_cubic(xFORzMin)
zMin = z_Interp_cubic(xFORzMin)
println("xFORzMin = ",xFORzMin,"   yATmin = ",yATmin,"   zMin = ",zMin)
InterpDisagreement = zMin-sqrt((xFORzMin^2)+(yATmin^2))
println("InterpDisagreement = ",InterpDisagreement)
yMinValDisagreement = yATmin - sqrt((zMin^2)-(xFORzMin^2))
println("yMinValDisagreement = ",yMinValDisagreement)
#
plot(xSample, abs.(xSample), xlim=(-2,3), ylim=(0,2.5)); plot!(xSample, abs.(y)); plot!(xSample, z)
scatter!(xSample, abs.(xSample)); scatter!(xSample, abs.(y)); scatter!(xSample, z)
xFine = -12.0:0.01:12.0
plot!(xFine,z_Interp_cubic(tFine), linewidth=3); scatter!([xFORzMin], [z_Interp_cubic(xFORzMin)], markersize=7)
scatter(xSample, y); plot!(xFine,y_Interp_cubic(tFine))

[anon92994695]

You should be able to change the tolerance of the optimization method to get a better approximation. That said your error will depend highly upon your interpolant, and grid spacing. The finer the grid the better your spline will fit. Similarly, you might obtain a closer match with a better interpolant function. For example, there’s no reason to expect your function will follow a cubic spline. Sometimes it’s good enough, but often times interpolation is more cosmetic/approximate than anything. We might want to ask ourselves, “should we be interpolating in the first place?”.

Some of the optimization methods benefit from things like derivatives, hessians, etc. I’d reccomend reading through the docs a bit and trying to see what you can come up with for your problem statement. Might get more lift that way. If your function is analytic, you can analytically determine your minima using calculus :D.

Big picture, we are opening up lots of degrees of freedom for us to get things wrong (add error). Barring that, I think I need more information about what you are trying to do. For some reason I think you will want to use Kriging here… But again, if you could specify your exact needs we can probably help you get where you’re trying to go.

[CosmoProf]

Hello ckneale, thanks for the feedback so far, I’m diving deep into Pandora’s box here by necessity.

I’ve determined that interpolation error is definitely a limiting problem right now for my calc accuracy.
The cubic splines work best so far (better than BSpline(Constant()), BSpline(Linear()), or BSpline(Quadratic())), but not yet good enough. I don’t know which are the best boundary conditions to use for the cubic spline, but I’m using commands like:

scale(interpolate(DataArray, BSpline(Cubic(Line(OnCell())))), xStepRange)

A finer grid would improve performance, but it rapidly increases computer time and resources used, so there’s a hard limit I will soon reach for grid spacing. So assume I can’t make the sampling any finer.

Your Kriging idea is very intriguing, do you think that’s the best direction to move to if I go beyond cubic splines? I have installed the packages KrigingEstimators, Geostats, and Mads. Can you suggest good commands to use for correctly implementing Kriging? Is there some other approach I should try?

For background info, I am numerically solving equations like the wave equation from physics, Ftt = Fxx, in two dimensions (t, x). My code will be adapted soon to more complicated wave-like equations with no known solutions, so I need purely numerical algorithms and results. I am using the Method of Lines, with discrete arrays and finite-difference derivatives in x, but using ODE solvers in t for each x value. After I get the solutions, I am writing an algorithm to “walk” through the waves and calculate things. The ODE solutions in t have built-in interpolations, but to step through in the x-direction I need to interpolate somehow by hand. Hopefully this bigger picture helps. Thanks for any more recommendations!

[anon92994695]

So I’ve been in a position like yours before. Not a ton of fun.

Boundary conditions for splines depends on a few factors. most importantly, how far from the boundary you are when you are interpolating, or heaven forbid extrapolating… For numerical optimization in my experience, you never want to be extrapolating or introducing error in the beginning…

What I am not super seasoned with - is ODE solvers. Sure I’ve solved some ODE’s but the nuances for the method you’re dealing with are outside of my wheelhouse. I am surprised that there’s no ‘ready made’ way to poke through the X axis for what you want? The DiffEq crew have set up a literal world class ecosystem - I mean - it’s the best I’ve ever seen. It might be worth asking them if you are using a hammer when you need a drill.

Barring that - I have been in a situation like this back when I was doing CEM. What I found to get the most result per unit compute was an adaptive screening. Is there anyway to “know” if you are near an X that’s interesting? If so - can you adaptively add resolution only when you approach those areas? This is a common technique, it sounds hacky, but even for fluid mech. or FEM you’ll see people doing this from time to time. In my case this brought 180 days of simulations down to something more like 5-7 after 2 days of coding + debugging.

If you’re dealing with waves, maybe you are looking for a mode, or a dip/ripple? Little things like that can save months of compute.

[CosmoProf]

Luckily I’m not extrapolating or going near the edge of the data points. Usually the interpolation will cover ~7-15 or so actual data points, and the minimum location will be within the ~3 central points.

I spent a month figuring out the right way to do the 2D PDE integration, I’m confident this is the best way. But I need better interpolation over discrete x. Do you have any specific knowledge of how to use the Kriging interpolations? That’s a completely new package I haven’t seen before. Any command suggestions or links would be helpful. Thanks again!

[jkbest2]

Maybe Surrogates.jl would be helpful? Not my area of expertise, but it looks like it deals with a very similar problem.

[DNF]

There’s a risk that I’ve missed some nuances of the discussion. But I don’t think a general optimization routine is needed to find optima of a cubic spline. You’re just dealing with a bunch of 3rd degree polynomials, each of which should have an easily accessible analytic optimum for each interval.

Calculating optima of splines should be a well-known and solved problem. If the spline/interpolation package doesn’t provide a built-in max-method, can’t you find the solution on Wikipedia?

[anon92994695]

jkbest12 - nice find this looks to be exactly what he actually needs.

DNF - you are correct. Once you have a fitted spline you can easily find the minima between each knot analytically. I think the issue here is that - the spline isn’t matching the data “well” enough to predict the correct minima? Could be wrong.

[CosmoProf]

Hi guys, to clarify, the issue is that I’m training my code on a simpler problem with a known answer, and the results are weirder and off by more than I’d expect. Some of the problem could be with my Partial Diff Eq’n solutions, but I think it’s more likely with the Interpolation method applied to the solutions after I have them. (The PDE solutions are continuous in t but discrete in x, and the interpolation over the discrete points produces the odd results I think.)

For example, the attached picture, obtained with cubic spline interpolation, should’ve been a horizontal line with value = 1. (Each dot is basically the distance of one step taken through the PDE solution to the wave equation. The destination of the step is figured out by the Interpolator. The variation from step to step should be very small and random.)

[CosmoProf]

Also, I’d appreciate any suggestions on how to apply Kriging interpolation on an array of discrete values.
(Each data array representing values of my function at different x locations, for a given time slice.)