CODEWITHSUNDEEP
pythonnumpy
Miss_Orchid
Miss_Orchid

Reputation: 374

numpy polyfit alternative for speed

I am utilizing numpy's polyfit a numerous amount of times to run calculations and get the slope between two datasets. However, the speed at which it performs these calculations are not fast enough for what would be desired.

Two things to note about the calculations:

  1. The value for x in the call numpy.polyfit(x,y,n) will always be the same slope value, and

  2. The value for n = 1. So it is a linear regression.

I know there are many different alternatives, including numpy.polynomial.polynomial.polyfit(x,y,n), but they seem to provide the same slow speed performance. I have had little luck getting np.linalg to work properly. Therefore, I am wondering what might be an alternative to speed up calculations?

Upvotes: 1

Views: 1206

Answers (1)

Nin17
Nin17

Reputation: 3492

As others have commented, this can be done using linear least squares. Using numpy.linalg.lstsq, this could look like:

import numpy as np

def lstsq(x, y):
    a = np.stack((x, np.ones_like(x)), axis=1)
    return np.linalg.lstsq(a, y)[0]

This offers a slight speed improvement over polyfit. To obtain a significant speed increase (at the expense of numerical stability - for a summary of methods see Numerical methods for linear least squares) you can instead solve the normal equations:

def normal(x, y):
    a = np.stack((x, np.ones_like(x)), axis=1)
    aT = a.T
    return np.linalg.solve(aT@a, aT@y)

As you say that x is constant, you can precompute a.T@a providing a further speed increase:

def normal2(aT, aTa,  y):
    return np.linalg.solve(aTa, aT@y)

Make up some test data and time:

rng = np.random.default_rng()
N = 1000

x = rng.random(N)
y = rng.random(N)

a = np.stack((x, np.ones_like(x)), axis=1)
aT = a.T
aTa = aT@a

assert np.allclose(lstsq(x, y), np.polyfit(x, y, 1))
assert np.allclose(normal(x, y), np.polyfit(x, y, 1))
assert np.allclose(normal2(aT, aTa, y), np.polyfit(x, y, 1))
%timeit np.polyfit(x, y, 1) 
%timeit lstsq(x, y)
%timeit normal(x, y)
%timeit normal2(aT, aTa, y)

Output:

256 µs ± 270 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
220 µs ± 1.87 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
20.2 µs ± 32.3 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
6.54 µs ± 13.5 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Upvotes: 4

Related Questions