完全矢量化 numpy polyfit
Fully vectorise numpy polyfit
概览
我 运行 遇到使用 polyfit 的性能问题，因为它似乎无法接受广播数组。我知道 from this post 如果您使用 numpy.polynomial.polynomial.polyfit，依赖数据 y 可以是多维的。但是，x 维度不能是多维的。反正有这个吗？
动机
我需要计算一些数据的变化率。为了与实验相匹配，我想使用以下方法：取数据 y 和 x，对于短部分数据拟合多项式，然后使用拟合系数作为变化率的估计。
插图
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()
蓝线是原始数据（正弦曲线），而绿线是一阶微分（余弦曲线）。
问题
为了对其进行矢量化，我执行了以下操作：
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]
这会将两个一维向量广播到形状为 (n-window_length, window_length) 的二维向量。现在我希望 polyfit 会有一个 axis 参数，这样我就可以并行计算，但没有这样的运气。
有人对如何执行此操作有任何建议吗？我愿意
很抱歉回答我自己的问题，但我花了 20 分钟来尝试掌握它，我有以下解决方案：
ydot = np.polynomial.polynomial.polyfit(x_array[0], y_array.T, 1)[-1]
一个令人困惑的部分是 np.polyfit returns 具有最高幂的系数 首先 。在 np.polynomial.polynomial.polyfit 中，最高功率是 last（因此 -1 而不是 0 索引）。
另一个混淆是我们只使用了 x (x_array[0]) 的第一片。我认为这没问题，因为使用的不是独立向量 x 的绝对值，而是它们之间的差值。或者，它就像更改参考 x 值。
如果有更好的方法来做到这一点，我仍然很高兴听到！
使用替代方法计算变化率可能是提高速度和准确性的解决方案。
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
相对误差：
结果：
polyfit 的工作方式是解决以下形式的最小二乘问题：
y = [X].a
其中y是你的依赖坐标，[X]是对应独立坐标的Vandermonde matrix，a是拟合系数的向量。
在您的情况下，您总是在计算一次多项式近似值，并且实际上只对一次项的系数感兴趣。这有一个 well known closed form solution 你可以在任何统计书籍中找到，或者通过创建一个 2x2 线性方程组来产生你自己，将上述方程的两边乘以 [X] 的转置。这一切加起来就是您要计算的值：
>>> 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
最重要的是，您可以在数据上滑动 window 运行，因此您可以使用类似于 1D summed area table 的内容，如下所示：
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)
使用此代码，您可以轻松检查（注意我将范围扩大了 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
并且：
%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
所以您的示例数据的速度大约快 300 倍。
概览
我 运行 遇到使用 polyfit 的性能问题，因为它似乎无法接受广播数组。我知道 from this post 如果您使用 numpy.polynomial.polynomial.polyfit，依赖数据 y 可以是多维的。但是，x 维度不能是多维的。反正有这个吗？
动机
我需要计算一些数据的变化率。为了与实验相匹配，我想使用以下方法：取数据 y 和 x，对于短部分数据拟合多项式，然后使用拟合系数作为变化率的估计。
插图
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()
蓝线是原始数据（正弦曲线），而绿线是一阶微分（余弦曲线）。
问题
为了对其进行矢量化，我执行了以下操作：
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]
这会将两个一维向量广播到形状为 (n-window_length, window_length) 的二维向量。现在我希望 polyfit 会有一个 axis 参数，这样我就可以并行计算，但没有这样的运气。
有人对如何执行此操作有任何建议吗？我愿意
很抱歉回答我自己的问题，但我花了 20 分钟来尝试掌握它，我有以下解决方案：
ydot = np.polynomial.polynomial.polyfit(x_array[0], y_array.T, 1)[-1]
一个令人困惑的部分是 np.polyfit returns 具有最高幂的系数 首先 。在 np.polynomial.polynomial.polyfit 中，最高功率是 last（因此 -1 而不是 0 索引）。
另一个混淆是我们只使用了 x (x_array[0]) 的第一片。我认为这没问题，因为使用的不是独立向量 x 的绝对值，而是它们之间的差值。或者，它就像更改参考 x 值。
如果有更好的方法来做到这一点，我仍然很高兴听到！
使用替代方法计算变化率可能是提高速度和准确性的解决方案。
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
相对误差：
结果：
polyfit 的工作方式是解决以下形式的最小二乘问题：
y = [X].a
其中y是你的依赖坐标，[X]是对应独立坐标的Vandermonde matrix，a是拟合系数的向量。
在您的情况下，您总是在计算一次多项式近似值，并且实际上只对一次项的系数感兴趣。这有一个 well known closed form solution 你可以在任何统计书籍中找到，或者通过创建一个 2x2 线性方程组来产生你自己，将上述方程的两边乘以 [X] 的转置。这一切加起来就是您要计算的值：
>>> 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
最重要的是，您可以在数据上滑动 window 运行，因此您可以使用类似于 1D summed area table 的内容，如下所示：
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)
使用此代码，您可以轻松检查（注意我将范围扩大了 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
并且：
%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
所以您的示例数据的速度大约快 300 倍。