Search code examples
pythonpython-2.7numpyscipyscientific-computing

How to improve the performance when 2d interpolating/smoothing lines using scipy?


I have a moderate size data set, namely 20000 x 2 floats in a two column matrix. The first column is the the x column which represents the distance to the original point along a trajectory, another column is the y column which represents the work has done to the object. This data set is obtained from lab operations, so it's fairly arbitrary. I've already turned this structure into numpy array. I want to plot y vs x in a figure with a smooth curve. So I hope the following code could help me:

x_smooth = np.linspace(x.min(),x.max(), 20000)
y_smooth = spline(x, y, x_smooth)
plt.plot(x_smooth, y_smooth)
plt.show()

However, when my program execute the line y_smooth = spline(x,y,x_smooth), it takes a very long time,say 10 min, and even sometimes it will blow my memory that I have to restart my machine. I tried to reduce the chunk number to 200 and 2000 and none of them works. Then I checked the official scipy reference: scipy.interpolate.spline here. And they said that spline is deprecated in v 0.19, but I'm not using the new version. If spline is deprecated for quite a bit of the time, how to use the equivalent Bspline now? If spline is still functioning, then what causes the slow performance

One portion of my data could look like this:

13.202      0.0
13.234738      -0.051354643759
12.999116      0.144464320836
12.86252      0.07396528119
13.1157      0.10019738758
13.357109      -0.30288563381
13.234004      -0.045792536285
12.836279      0.0362257166275
12.851597      0.0542649286915
13.110691      0.105297378401
13.220619      -0.0182963209185
13.092143      0.116647353635
12.545676      -0.641112204849
12.728248      -0.147460703493
12.874176      0.0755861585235
12.746764      -0.111583725833
13.024995      0.148079528382
13.106033      0.119481137144
13.327233      -0.197666132456
13.142423      0.0901867159545

Solution

  • Several issues here. First and foremost, spline fitting you're trying to use is global. This means that you're solving a system of linear equations of the size 20000 at the construction time (evaluations are weakly sensitive to the dataset size though). This explains why the spline construction is slow.

    scipy.interpolate.spline, furthermore, does linear algebra with full matrices --- hence memory consumption. This is precisely why it's deprecated from scipy 0.19.0 on.

    The recommended replacement, available in scipy 0.19.0, is the BSpline/ make_interp_spline combo:

    >>> spl = make_interp_spline(x, y, k=3)    # returns a BSpline object
    >>> y_new = spl(x_new)                     # evaluate 
    

    Notice it is not BSpline(x, y, k): BSpline objects do not know anything about the data or fitting or interpolation.

    If you are using older scipy versions, your options are:

    • CubicSpline(x, y) for cubic splines
    • splrep(x, y, s=0) / splev combo.

    However, you may want to think if you really need twice continuously differentiable functions. If only once differentiable functions are smooth enough for your purposes, then you can use local spline interpolations, e.g. Akima1DInterpolator or PchipInterpolator:

    In [1]: import numpy as np
    
    In [2]: from scipy.interpolate import pchip, splmake
    
    In [3]: x = np.arange(1000)
    
    In [4]: y = x**2
    
    In [5]: %timeit pchip(x, y)
    10 loops, best of 3: 58.9 ms per loop
    
    In [6]: %timeit splmake(x, y)    
    1 loop, best of 3: 5.01 s per loop
    

    Here splmake is what spline uses under the hood, and it's also deprecated.