19

I am trying to invert an interpolated function using scipy's interpolate function. Let's say I create an interpolated function,

import scipy.interpolate as interpolate
interpolatedfunction = interpolated.interp1d(xvariable,data,kind='cubic')

Is there some function that can find x when I specify a:

interpolatedfunction(x) == a

In other words, "I want my interpolated function to equal a; what is the value of xvariable such that my function is equal to a?"

I appreciate I can do this with some numerical scheme, but is there a more straightforward method? What if the interpolated function is multivalued in xvariable?

2
  • 1
    so you want x where interpolatedfunction(x) == a?
    – James
    Nov 14, 2017 at 0:56
  • That is correct. Is there a built-in function that does such a thing?
    – Physics314
    Nov 14, 2017 at 0:58

4 Answers 4

18

There are dedicated methods for finding roots of cubic splines. The simplest to use is the .roots() method of InterpolatedUnivariateSpline object:

spl = InterpolatedUnivariateSpline(x, y)
roots = spl.roots()

This finds all of the roots instead of just one, as generic solvers (fsolve, brentq, newton, bisect, etc) do.

x = np.arange(20)
y = np.cos(np.arange(20))
spl = InterpolatedUnivariateSpline(x, y)
print(spl.roots())

outputs array([ 1.56669456, 4.71145244, 7.85321627, 10.99554642, 14.13792756, 17.28271674])

However, you want to equate the spline to some arbitrary number a, rather than 0. One option is to rebuild the spline (you can't just subtract a from it):

solutions = InterpolatedUnivariateSpline(x, y - a).roots()

Note that none of this will work with the function returned by interp1d; it does not have roots method. For that function, using generic methods like fsolve is an option, but you will only get one root at a time from it. In any case, why use interp1d for cubic splines when there are more powerful ways to do the same kind of interpolation?

Non-object-oriented way

Instead of rebuilding the spline after subtracting a from data, one can directly subtract a from spline coefficients. This requires us to drop down to non-object-oriented interpolation methods. Specifically, sproot takes in a tck tuple prepared by splrep, as follows:

tck = splrep(x, y, k=3, s=0)
tck_mod = (tck[0], tck[1] - a, tck[2])
solutions = sproot(tck_mod)    

I'm not sure if messing with tck is worth the gain here, as it's possible that the bulk of computation time will be in root-finding anyway. But it's good to have alternatives.

1
  • 1
    Excellent alternative as .roots() finds all of the roots (if they exist).
    – Physics314
    Nov 14, 2017 at 11:17
16

After creating an interpolated function interp_fn, you can find the value of x where interp_fn(x) == a by the roots of the function

interp_fn2 = lambda x: interp_fn(x) - a

There are number of options to find the roots in scipy.optimize. For instance, to use Newton's method with the initial value starting at 10:

from scipy import optimize

optimize.newton(interp_fn2, 10)

Actual example

Create an interpolated function and then find the roots where fn(x) == 5

import numpy as np
from scipy import interpolate, optimize

x = np.arange(10)
y = 1 + 6*np.arange(10) - np.arange(10)**2
y2 = 5*np.ones_like(x)
plt.scatter(x,y)
plt.plot(x,y)
plt.plot(x,y2,'k-')
plt.show()

Plot of f(x) and y=5

# create the interpolated function, and then the offset
# function used to find the roots

interp_fn = interpolate.interp1d(x, y, 'quadratic')
interp_fn2 = lambda x: interp_fn(x)-5

# to find the roots, we need to supply a starting value
# because there are more than 1 root in our range, we need 
# to supply multiple starting values.  They should be 
# fairly close to the actual root

root1, root2 = optimize.newton(interp_fn2, 1), optimize.newton(interp_fn2, 5)

root1, root2
# returns:
(0.76393202250021064, 5.2360679774997898)
1
  • Very nice solution!
    – Pieter-Jan
    May 20, 2022 at 15:51
9

If your data are monotonic you might also try the following:

inversefunction = interpolated.interp1d(data, xvariable, kind='cubic')
1
  • 1
    This is brilliant. Especially useful for sampling from custom distributions.
    – marnix
    Feb 22, 2021 at 18:08
1

Mentioning another option because I found this page in a google search and the other option works for my simple use case. Hopefully it'll be of use to someone.

If the function you're interpolating is very simple and always has a 1:1 relationship between y and x, then you can simply take your data, swap x and y when you pass it into interp1d, and then call the interpolation function in that direction.

Adapting code from https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
x = np.arange(0, 10)
y = np.exp(-x/3.0)
f = interpolate.interp1d(x, y)
xnew = np.arange(0, 9, 0.1)
ynew = f(xnew)
plt.plot(x, y, 'o', xnew, ynew, '-')
plt.show()

enter image description here

When x and y have been swapped you can call swappedInterpolationFunction(a) to get the x value where that would occur.

f = interpolate.interp1d(y, x)
xnew = np.arange(np.exp(-9/3), np.exp(0), 0.01)
ynew = f(xnew)
plt.plot(y, x, 'o', xnew, ynew, '-')
plt.title("Inverted")
plt.show()

enter image description here

Of course, if the function ever has multiple x values for a given y value (like sine or a parabola) then this will not work because it will no longer be a 1:1 function from x to y, and the above answers are necessary. This is just a simplification in a limited use case.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.