完全向量化numpy polyfit

Question

概述

使用 polyfit 时，我遇到了性能问题，因为它似乎无法接受广播数组。我知道来自此帖子的信息，如果使用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参数，这样我就可以并行计算，但不幸的是没有。有人有什么建议吗？我很乐意听取意见。

- Greg

---

计算一阶数值导数而非使用polyfit应该更快、更准确。 - M4rtini

@M4rtini，你说得对，但我这样做是为了与实验者使用的一种方法保持一致。这在问题中已经提到了，但我知道有太多的文本让任何人都不想去读它。 - Greg

哦，我错过了问题的那一部分。我猜我的答案完全无关紧要了。 - M4rtini

Answer 1

polyfit 的工作原理是通过解决形式为最小二乘问题的方式实现的：

y = [X].a

其中y是您的依赖坐标，[X]是相应独立坐标的范德蒙矩阵，a是拟合系数的向量。

在您的情况下，您始终计算一次多项式逼近，并且实际上只对一次多项式系数感兴趣。这有一个已知的封闭形式解决方案，您可以在任何统计书籍中找到，或通过乘以以上方程的转置来创建2x2线性方程组自行生成。这一切都加起来，您要计算的值是:

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

此外，您可以在数据上运行滑动窗口，因此您可以使用类似于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倍。

- Jaime

Answer 2

非常抱歉我要回答自己的问题，但是经过20分钟的努力，我得到了以下解决方案：

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

< p>有些令人困惑的是，np.polyfit返回的系数以最高次幂优先。在np.polynomial.polynomial.polyfit中，最高阶次幂是最后的（因此使用-1而不是0索引）。

另一个令人困惑的是我们只使用了x的第一部分（x_array [0] ）。我认为这样做是可以的，因为使用的不是自变量x的绝对值，而是它们之间的差异。或者换句话说，就像更改参考x值一样。

如果有更好的方法，我还是很乐意听取建议！

- Greg

Answer 3

使用另一种计算变化率的方法可能是提高速度和准确性的解决方案。

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 

相对误差： 结果：

- M4rtini

---

正如问题评论中所讨论的那样，我想要解决的确切问题就是这些错误！虽然我非常感谢您在这里所做出的努力 - 实际上我之前并没有意识到这种方法有多糟糕！因此，我会将这个留在这里。 - Greg

我记得数值数学课上讲过，多项式本质上不适合用于拟合谐波函数。因此，对于其他类型的数据，相对误差可能不会这么大。尽管如此，速度应该仍然是相同的。我会把这个留下来供将来参考，即使它没有直接回答问题。 - M4rtini