Fully vectorise numpy polyfit

Question

Overview

I am running into issues with performance using polyfit because it doesn't appear able to accept broadcast arrays. I am aware from this post that the dependant data y can be multidimensional if you use numpy.polynomial.polynomial.polyfit. However, the x dimension cannot be multidimensional. Is there anyway around this?

Motivation

I need to compute the rate of change of some data. To match with an experiment I want to use the following method: take data y and x, for short sections of data fit a polynomial, then use the fitted coefficient as an estimate of the rate of change.

Illustration

import numpy as np
import matplotlib.pyplot as plt

n = 100
x = np.linspace(0, 10, n)
y = np.sin(x)

window_length = 10
ydot = [np.polyfit(x[j:j+window_length], y[j:j+window_length], 1)[0] 
                                  for j in range(n - window_length)]
x_mids = [x[j+window_length/2] for j in range(n - window_length)]

plt.plot(x, y)
plt.plot(x_mids, ydot)

plt.show()

The blue line is the original data (a sine curve), while the green is the first differential (a cosine curve).

The problem

To vectorise this I did the following:

window_length = 10
vert_idx_list = np.arange(0, len(x) - window_length, 1)
hori_idx_list = np.arange(window_length)
A, B = np.meshgrid(hori_idx_list, vert_idx_list)
idx_array = A + B 

x_array = x[idx_array]
y_array = y[idx_array]

This broadcasts the two 1D vectors to 2D vectors of shape (n-window_length, window_length). Now I was hoping that polyfit would have an axis argument so I could parallelise the calculation, but no such luck.

Does anyone have any suggestion for how to do this? I am open to

Answer 1

The way polyfit works is by solving a least-square problem of the form:

y = [X].a

where y are your dependent coordinates, [X] is the Vandermonde matrix of the corresponding independent coordinates, and a is the vector of fitted coefficients.

In your case you are always computing a 1st degree polynomial approximation, and are actually only interested in the coefficient of the 1st degree term. This has a well known closed form solution you can find in any statistics book, or produce your self by creating a 2x2 linear system of equation premultiplying both sides of the above equation by the transpose of [X]. This all adds up to the value you want to calculate being:

>>> n = 10
>>> x = np.random.random(n)
>>> y = np.random.random(n)
>>> np.polyfit(x, y, 1)[0]
-0.29207474654700277
>>> (n*(x*y).sum() - x.sum()*y.sum()) / (n*(x*x).sum() - x.sum()*x.sum())
-0.29207474654700216

On top of that you have a sliding window running over your data, so you can use something akin to a 1D summed area table as follows:

def sliding_fitted_slope(x, y, win):
    x = np.concatenate(([0], x))
    y = np.concatenate(([0], y))
    Sx = np.cumsum(x)
    Sy = np.cumsum(y)
    Sx2 = np.cumsum(x*x)
    Sxy = np.cumsum(x*y)

    Sx = Sx[win:] - Sx[:-win]
    Sy = Sy[win:] - Sy[:-win]
    Sx2 = Sx2[win:] - Sx2[:-win]
    Sxy = Sxy[win:] - Sxy[:-win]

    return (win*Sxy - Sx*Sy) / (win*Sx2 - Sx*Sx)

With this code you can easily check that (notice I extended the range by 1):

>>> np.allclose(sliding_fitted_slope(x, y, window_length),
                [np.polyfit(x[j:j+window_length], y[j:j+window_length], 1)[0]
                 for j in range(n - window_length + 1)])
True

And:

%timeit sliding_fitted_slope(x, y, window_length)
10000 loops, best of 3: 34.5 us per loop

%%timeit
[np.polyfit(x[j:j+window_length], y[j:j+window_length], 1)[0]
 for j in range(n - window_length + 1)]
100 loops, best of 3: 10.1 ms per loop

So it is about 300x faster for your sample data.

Answer 2

Using an alternative method for calculating the rate of change may be the solution for both speed and accuracy increase.

n = 1000
x = np.linspace(0, 10, n)
y = np.sin(x)


def timingPolyfit(x,y):
    window_length = 10
    vert_idx_list = np.arange(0, len(x) - window_length, 1)
    hori_idx_list = np.arange(window_length)
    A, B = np.meshgrid(hori_idx_list, vert_idx_list)
    idx_array = A + B 

    x_array = x[idx_array]
    y_array = y[idx_array]

    ydot = np.polynomial.polynomial.polyfit(x_array[0], y_array.T, 1)[-1]
    x_mids = [x[j+window_length/2] for j in range(n - window_length)]

    return ydot, x_mids

def timingSimple(x,y):
    dy = (y[2:] - y[:-2])/2
    dx =  x[1] - x[0]
    dydx = dy/dx
    return dydx, x[1:-1]

y1, x1 = timingPolyfit(x,y)
y2, x2 = timingSimple(x,y)

polyfitError = np.abs(y1 - np.cos(x1))
simpleError = np.abs(y2 - np.cos(x2))

print("polyfit average error: {:.2e}".format(np.average(polyfitError)))
print("simple average error: {:.2e}".format(np.average(simpleError)))

result = %timeit -o timingPolyfit(x,y)
result2 = %timeit -o timingSimple(x,y)

print("simple is {0} times faster".format(result.best / result2.best))

---

polyfit average error: 3.09e-03 
simple average error: 1.09e-05 
100 loops, best of 3: 3.2 ms per loop 
100000 loops, best of 3: 9.46 µs per loop 
simple is 337.995634151131 times faster 

Relative error:

Results:

Answer 3

Sorry for answering my own question, but 20 minutes more of trying to get to grips with it I have the following solution:

ydot = np.polynomial.polynomial.polyfit(x_array[0], y_array.T, 1)[-1]

One confusing part is that np.polyfit returns the coefficients with the highest power first. In np.polynomial.polynomial.polyfit the highest power is last (hence the -1 instead of 0 index).

Another confusion is that we use only the first slice of x (x_array[0]). I think that this is okay because it is not the absolute values of the independent vector x that are used, but the difference between them. Or alternatively it is like changing the reference x value.

If there is a better way to do this I am still happy to hear about it!