18

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?

Physics314
  • 341
  • 1
  • 2
  • 9

4 Answers4

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.

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)
James
  • 32,991
  • 4
  • 47
  • 70
8

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

inversefunction = interpolated.interp1d(data, xvariable, kind='cubic')
Ferenc Lengyel
  • 71
  • 1
  • 9
Elin
  • 81
  • 1
  • 1
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.

Eric
  • 81
  • 5