5

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

enter image description here

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

3
  • Calculating the first order numerical derivative instead of using polyfit should be faster and more accurate
    – M4rtini
    Commented Jan 30, 2015 at 15:02
  • @M4rtini you are correct, but I'm doing it this way to be consistent with a method used by experimentalist. It is in the question but I appreciate there is way too much text for anyone to be bothered reading.
    – Greg
    Commented Jan 30, 2015 at 15:31
  • Ohh, missed that part of the question. I guess my answer is totally irrelevant then
    – M4rtini
    Commented Jan 30, 2015 at 15:35

3 Answers 3

8

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.

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!

1

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: Relative error

Results: Results-closeup

2
  • 1
    As discussed in the comments of the question, it's exactly these errors I want to account for! Although I really appreciate the effort you made here - I actually hadn't realised just how bad this method is! For that reason I would leave this here!
    – Greg
    Commented Jan 30, 2015 at 15:49
  • 1
    I seem to remember from numerical-math classes that polynomials are inherently bad for fitting harmonic functions. So the relative error may not be this bad for other types of data. The speed should still be the same though. I'll leave this here for future reference then, even though it doesent directly answer the question.
    – M4rtini
    Commented Jan 30, 2015 at 15:57

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.