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

Answer 1

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:

1. Use a 4th degree spline for interpolation, so that the roots of its derivative can be found easily.
2. 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