Find maximum/minimum of a 1d interpolated function

Question

I have a set of data which I am interpolating with kind = 'cubic'.

I would like to find the maximum of this cubic interpolation function.

Currently what I am doing is just find the maximum value in the array of interpolated data, but I was wondering whether the interpolated function, as an object, can be differentiated to find its extrema?

Code:

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

x_axis = np.array([ 2.14414414,  2.15270826,  2.16127238,  2.1698365 ,  2.17840062, 2.18696474,  2.19552886,  2.20409298,  2.2126571 ,  2.22122122])
y_axis = np.array([ 0.67958442,  0.89628424,  0.78904004,  3.93404167,  6.46422317, 6.40459954,  3.80216674,  0.69641825,  0.89675386,  0.64274198])

f = interp1d(x_axis, y_axis, kind = 'cubic')

x_new = np.linspace(x_axis[0], x_axis[-1],100)

fig = plt.subplots()
plt.plot(x_new, f(x_new))

a cubic function `ax**3+bx**2+cx+d` has its (local) [extrema](http://mathworld.wolfram.com/Extremum.html) on `-b+-sqrt(b**2-3ac)/3a` — Ma0, May 16 '18 at 12:47
but I don't know what the a, b,c and d parameters are for my cubic function? It's not a global cubic function though, so I canno just use global paramters/ — SuperCiocia, May 16 '18 at 12:48
not fast or anything but you can evaluate `f` repeatedly on different values incrementing by `np.finfo(float).eps` or some multiple thereof. — Mohammad Athar, May 16 '18 at 12:57
I once wrote my own cubic interpolation function, it's not too difficult. I adapted it over time, but at some point posted the code as an [answer](https://stackoverflow.com/a/49736163/6012085). Perhaps you could use it. — Peter9192, May 16 '18 at 13:01
In case you want an approximate solution, you could perform [parabolic interpolation](https://en.wikipedia.org/wiki/Successive_parabolic_interpolation) from the maximum and its two neighbors. — heltonbiker, May 16 '18 at 13:15

score 11 · Accepted Answer · edited Aug 11 '20 at 23:07

The derivative of a cubic spline is a quadratic spline. SciPy only has a built-in method to find the roots of a cubic spline. So there are two approaches:

Use a 4th degree spline for interpolation, so that the roots of its derivative can be found easily.
Use a cubic spline (which is often preferable), and write a custom function for the roots of its derivative.

I describe both solutions below.

4th degree spline

Use InterpolatedUnivariateSpline.which had .derivative method returning a cubic spline, to which .roots method can be applied.

from scipy.interpolate import InterpolatedUnivariateSpline
f = InterpolatedUnivariateSpline(x_axis, y_axis, k=4)
cr_pts = f.derivative().roots()
cr_pts = np.append(cr_pts, (x_axis[0], x_axis[-1]))  # also check the endpoints of the interval
cr_vals = f(cr_pts)
min_index = np.argmin(cr_vals)
max_index = np.argmax(cr_vals)
print("Maximum value {} at {}\nMinimum value {} at {}".format(cr_vals[max_index], cr_pts[max_index], cr_vals[min_index], cr_pts[min_index]))

Output:

Maximum value 6.779687224066201 at 2.1824928509277037
Minimum value 0.34588448400295346 at 2.2075868177297036

Cubic spline

We need a custom function for the roots of a quadratic spline. Here it is (explained below).

def quadratic_spline_roots(spl):
    roots = []
    knots = spl.get_knots()
    for a, b in zip(knots[:-1], knots[1:]):
        u, v, w = spl(a), spl((a+b)/2), spl(b)
        t = np.roots([u+w-2*v, w-u, 2*v])
        t = t[np.isreal(t) & (np.abs(t) <= 1)]
        roots.extend(t*(b-a)/2 + (b+a)/2)
    return np.array(roots)

Now proceed exactly as above, except for using the custom solver.

from scipy.interpolate import InterpolatedUnivariateSpline
f = InterpolatedUnivariateSpline(x_axis, y_axis, k=3)
cr_pts = quadratic_spline_roots(f.derivative())
cr_pts = np.append(cr_pts, (x_axis[0], x_axis[-1]))  # also check the endpoints of the interval
cr_vals = f(cr_pts)
min_index = np.argmin(cr_vals)
max_index = np.argmax(cr_vals)
print("Maximum value {} at {}\nMinimum value {} at {}".format(cr_vals[max_index], cr_pts[max_index], cr_vals[min_index], cr_pts[min_index]))

Output:

Maximum value 6.782781181150518 at 2.1824928579767167
Minimum value 0.45017143148176136 at 2.2070746522580795

The slight discrepancy with the output in the first method is not a bug; the 4th degree spline and 3rd degree spline are a little different.

Explanation of `quadratic_spline_roots`

Suppose we know the values of a quadratic polynomial at -1, 0, 1 are u, v, w. What are its roots on the interval [-1, 1]? With some algebra we can find that the polynomial is

((u+w-2*v) * x**2 + (w-u) * x + 2*v) / 2

Now the quadratic formula can be used, but it's better to use np.roots because it will also handle the case of leading coefficient being zero. The roots are then filtered to real numbers between -1 to 1. Finally, if the interval is some [a, b] instead of [-1, 1], a linear transformation is made.

Bonus: the width of cubic spline at midrange

Suppose we want to find where the spline takes the value equal to the average of its maximum and minimum (i.e., its midrange). Then we should definitely use the cubic spline for interpolation, because the roots method will now be needed for it. One can't just do (f - mid_range).roots(), as the addition of a constant to a spline is not supported in SciPy. Instead, build a shifted down spline from y_axis - mid_range.

mid_range = (cr_vals[max_index] + cr_vals[min_index])/2
f_shifted = InterpolatedUnivariateSpline(x_axis, y_axis - mid_range, k=3)
roots = f_shifted.roots()
print("Mid-range attained from {} to {}".format(roots.min(), roots.max()))

Mid-range attained from 2.169076230034363 to 2.195974299834667

Hey thanks so much this works beautifully. So what if now I want to find the full width at half maximum? Can I redefine f by f = f - cr_vals[max_index]/2 (based on the first method) and then solve it to find the x-axis positions? — SuperCiocia, May 17 '18 at 18:50
That won't work for two reasons: "spline - number" is not supported; and `roots` only works for cubic splines. I revised the second part of the answer (using cubic interpolation): this is what you'll need. — , May 17 '18 at 19:57
wouldn't all the maxima and minima be all the spots on the interpolated curve where the slope equals 0? is there an easy way to find where a specific slope exists on a curve? — Drafter250, Sep 01 '20 at 19:18

Find maximum/minimum of a 1d interpolated function

1 Answers1

4th degree spline

Cubic spline

Explanation of quadratic_spline_roots

Bonus: the width of cubic spline at midrange

Explanation of `quadratic_spline_roots`