python - 计算 Vandermonde 矩阵的有效方法

Question

我正在计算 Vandermonde matrix对于相当大的一维阵列。做到这一点的自然而干净的方法是使用 np.vander() .但是，我发现这大约是。比基于列表理解的方法慢 2.5 倍。

In [43]: x = np.arange(5000)
In [44]: N = 4

In [45]: %timeit np.vander(x, N, increasing=True)
155 µs ± 205 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

# one of the listed approaches from the documentation
In [46]: %timeit np.flip(np.column_stack([x**(N-1-i) for i in range(N)]), axis=1)
65.3 µs ± 235 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [47]: np.all(np.vander(x, N, increasing=True) == np.flip(np.column_stack([x**(N-1-i) for i in range(N)]), axis=1))
Out[47]: True

我试图了解瓶颈在哪里以及执行原生 np.vander() 的原因是 ~ 2.5x 慢。

效率对我的实现很重要。因此，我们也欢迎使用更快的替代方案!

Answer 1

广播求幂怎么样？

%timeit (x ** np.arange(N)[:, None]).T
43 µs ± 348 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

完整性检查 -

np.all((x ** np.arange(N)[:, None]).T == np.vander(x, N, increasing=True))
True

这里需要注意的是，只有当您的输入数组 x 的 dtype 为 int 时，这种加速才有可能。正如@Warren Weckesser 在评论中指出的那样， float 组的广播求幂速度变慢了。

---

至于为什么np.vander慢，看source code -

x = asarray(x)
if x.ndim != 1:
    raise ValueError("x must be a one-dimensional array or sequence.")
if N is None:
    N = len(x)

v = empty((len(x), N), dtype=promote_types(x.dtype, int))
tmp = v[:, ::-1] if not increasing else v

if N > 0:
    tmp[:, 0] = 1
if N > 1:
    tmp[:, 1:] = x[:, None]
    multiply.accumulate(tmp[:, 1:], out=tmp[:, 1:], axis=1)

return v

该函数必须满足除您之外的更多用例，因此它使用更通用的计算方法，该方法可靠但速度较慢(我特别指出 multiply.accumulate) .

---

出于兴趣，我找到了另一种计算 Vandermonde 矩阵的方法，结果如下:

%timeit x[:, None] ** np.arange(N)
150 µs ± 230 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

它做同样的事情，但速度要慢得多。答案在于操作是广播的，但效率很低。

另一方面，对于 float 数组，这实际上最终表现最佳。