多个维度的 NumPy PolyFit 和 PolyVal？

Posted 2023-03-12

【中文标题】多个维度的 NumPy PolyFit 和 PolyVal？

【英文标题】：NumPy PolyFit and PolyVal in Multiple Dimensions?

【发布时间】：2013-12-10 17:57:09

【问题描述】：

假设一个 n 维的观察数组被重新整形为一个二维数组，每一行是一个观察集。使用这种重塑方法，np.polyfit 可以计算整个 ndarray（矢量化）的二阶拟合系数：

fit = np.polynomial.polynomialpolyfit(X, Y, 2)

其中 Y 是形状 (304000, 21)，X 是向量。这会产生一个 (304000,3) 的系数数组，fit。

使用迭代器可以为每一行调用np.polyval(fit, X)。当可能存在矢量化方法时，这是低效的。 fit 结果可以在不迭代的情况下应用于整个观察数组吗？如果有，怎么做？

这与this SO question 类似。

【问题讨论】：

仅供参考，不要在来自np.polynomial.polynomial.polyfit 的结果上调用np.polyval。使用np.polynomial.polynomial.polyval。见***.com/a/18767992/1730674。

一种更优雅但仍然很慢的方法是将polyfit与np.apply_along_axis...one example here一起使用

@askewchan 绝对！完整的通话路径让我很懒惰。 @Saullo Castro - 正如您所建议的那样，np.apply_along_axis 并不比 [i,j] 迭代器快。我想知道是否存在真正的矢量化（在 C 级别）方法。

【参考方案1】：

np.polynomial.polynomial.polyval 是一种非常精细（且方便）的有效评估多项式拟合的方法。

但是，如果您正在寻找“最快”，则只需构造多项式输入并使用基本的 numpy 矩阵乘法函数，计算速度就会稍快（大约快 4 倍）。

设置

使用与上述相同的设置，我们将创建 25 个不同的线配件。

>>> num_samples = 100000
>>> num_lines = 100
>>> x = np.random.randint(0,100,num_samples)
>>> y = np.random.randint(0,100,(num_samples, num_lines))
>>> fit = np.polyfit(x,y,deg=2)
>>> xx = np.random.randint(0,100,num_samples*10)

Numpy 的polyval 函数

res1 = np.polynomial.polynomial.polyval(xx, fit)

基本矩阵乘法

inputs = np.array([np.power(xx,d) for d in range(len(fit))])
res2 = fit.T.dot(inputs)

定时功能

使用上面相同的参数...

%timeit _ = np.polynomial.polynomial.polyval(xx, fit)
1 loop, best of 3: 247 ms per loop

%timeit inputs = np.array([np.power(xx, d) for d in range(len(fit))]);_ = fit.T.dot(inputs)
10 loops, best of 3: 72.8 ms per loop

打败一匹死马……

平均效率提升约 3.61 倍。速度波动可能来自后台的随机计算机进程。

【参考方案2】：

np.polynomial.polynomial.polyval 采用多维系数数组：

>>> x = np.random.rand(100)
>>> y = np.random.rand(100, 25)
>>> fit = np.polynomial.polynomial.polyfit(x, y, 2)
>>> fit.shape # 25 columns of 3 polynomial coefficients
(3L, 25L)
>>> xx = np.random.rand(50)
>>> interpol = np.polynomial.polynomial.polyval(xx, fit)
>>> interpol.shape # 25 rows, each with 50 evaluations of the polynomial
(25L, 50L)

当然还有：

>>> np.all([np.allclose(np.polynomial.polynomial.polyval(xx, fit[:, j]),
...                     interpol[j]) for j in range(25)])
True