计算 Vandermonde 矩阵的有效方法
Efficient way to compute the Vandermonde matrix
我正在计算 Vandermonde matrix for a fairly large 1D array. The natural and clean way to do this is using 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 的原因。
效率对我的实施很重要。因此，我们也欢迎更快的替代品！
广播求幂如何？
%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 在评论中指出的那样，浮点数组的广播求幂速度变慢了。
至于为什么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）。
出于兴趣，我找到了另一种计算范德蒙矩阵的方法，结果如下：
%timeit x[:, None] ** np.arange(N)
150 µs ± 230 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
它做同样的事情，但速度要慢得多。答案在于操作是广播的，但是效率低下。
另一方面，对于 float 个数组，这实际上最终表现最好。
这里有一些方法，其中一些方法（在我的电脑上）比目前发布的方法快很多。
我认为最重要的观察是，这在很大程度上取决于你想要多少学位。求幂（我认为这是小整数指数的特殊情况）仅对小指数范围有意义。指数越多，基于乘法的方法越好。
我想强调一个基于 multiply.accumulate 的方法 (ma)，它类似于 numpy 的内置方法但速度更快（并不是因为我在检查时略过 - nnc， numpy-no-checks 证明了这一点）。除了最小的指数范围外，它实际上对我来说是最快的。
由于我不明白的原因，据我所知，numpy 实现做了三件缓慢且不必要的事情：(1) 它制作了相当多的基本向量副本。 (2) 它使它们不连续。 (3) 它就地积累，我相信这会强制缓冲。
我想提的另一件事是，对于小范围的整数（out_e_1 本质上是 ma 的手动版本）最快的速度会慢两倍以上通过升级到更大的 dtype 的简单预防措施（safe_e_1 可以说有点用词不当）。
基于广播的方法称为 bc_*，其中 * 表示广播轴（b 表示基数，e 表示 exp）'cheat' 表示结果不连续。
计时（3 个中的最佳）：
rep=100 n_b=5000 n_e=4 b_tp=<class 'numpy.int32'> e_tp=<class 'numpy.int32'>
vander 0.16699657 ms
bc_b 0.09595204 ms
bc_e 0.07959786 ms
ma 0.10755240 ms
nnc 0.16459018 ms
out_e_1 0.02037535 ms
out_e_2 0.02656622 ms
safe_e_1 0.04652272 ms
safe_e_2 0.04081079 ms
cheat bc_e_cheat 0.04668466 ms
rep=100 n_b=5000 n_e=8 b_tp=<class 'numpy.int32'> e_tp=<class 'numpy.int32'>
vander 0.25086462 ms
bc_b apparently failed
bc_e apparently failed
ma 0.15843041 ms
nnc 0.24713077 ms
out_e_1 apparently failed
out_e_2 apparently failed
safe_e_1 0.15970622 ms
safe_e_2 0.19672418 ms
bc_e_cheat apparently failed
rep=100 n_b=5000 n_e=4 b_tp=<class 'float'> e_tp=<class 'numpy.int32'>
vander 0.16225773 ms
bc_b 0.53315020 ms
bc_e 0.56200830 ms
ma 0.07626799 ms
nnc 0.16059748 ms
out_e_1 0.03653416 ms
out_e_2 0.04043702 ms
safe_e_1 0.04060494 ms
safe_e_2 0.04104209 ms
cheat bc_e_cheat 0.52966076 ms
rep=100 n_b=5000 n_e=8 b_tp=<class 'float'> e_tp=<class 'numpy.int32'>
vander 0.24542852 ms
bc_b 2.03353578 ms
bc_e 2.04281270 ms
ma 0.11075758 ms
nnc 0.24212880 ms
out_e_1 0.14809043 ms
out_e_2 0.19261359 ms
safe_e_1 0.15206112 ms
safe_e_2 0.19308420 ms
cheat bc_e_cheat 1.99176601 ms
代码：
import numpy as np
import types
from timeit import repeat
prom={np.dtype(np.int32): np.dtype(np.int64), np.dtype(float): np.dtype(float)}
def RI(k, N, dt, top=100):
return np.random.randint(0, top if top else N, (k, N)).astype(dt)
def RA(k, N, dt, top=None):
return np.add.outer(np.zeros((k,), int), np.arange(N)%(top if top else N)).astype(dt)
def RU(k, N, dt, top=100):
return (np.random.random((k, N))*(top if top else N)).astype(dt)
def data(k, N_b, N_e, dt_b, dt_e, b_fun=RI, e_fun=RA):
b = list(b_fun(k, N_b, dt_b))
e = list(e_fun(k, N_e, dt_e))
return b, e
def f_vander(b, e):
return np.vander(b, len(e), increasing=True)
def f_bc_b(b, e):
return b[:, None]**e
def f_bc_e(b, e):
return np.ascontiguousarray((b**e[:, None]).T)
def f_ma(b, e):
out = np.empty((len(b), len(e)), prom[b.dtype])
out[:, 0] = 1
np.multiply.accumulate(np.broadcast_to(b, (len(e)-1, len(b))), axis=0, out=out[:, 1:].T)
return out
def f_nnc(b, e):
out = np.empty((len(b), len(e)), prom[b.dtype])
out[:, 0] = 1
out[:, 1:] = b[:, None]
np.multiply.accumulate(out[:, 1:], out=out[:, 1:], axis=1)
return out
def f_out_e_1(b, e):
out = np.empty((len(b), len(e)), b.dtype)
out[:, 0] = 1
out[:, 1] = b
out[:, 2] = c = b*b
for i in range(3, len(e)):
c*=b
out[:, i] = c
return out
def f_out_e_2(b, e):
out = np.empty((len(b), len(e)), b.dtype)
out[:, 0] = 1
out[:, 1] = b
out[:, 2] = b*b
for i in range(3, len(e)):
out[:, i] = out[:, i-1] * b
return out
def f_safe_e_1(b, e):
out = np.empty((len(b), len(e)), prom[b.dtype])
out[:, 0] = 1
out[:, 1] = b
out[:, 2] = c = (b*b).astype(prom[b.dtype])
for i in range(3, len(e)):
c*=b
out[:, i] = c
return out
def f_safe_e_2(b, e):
out = np.empty((len(b), len(e)), prom[b.dtype])
out[:, 0] = 1
out[:, 1] = b
out[:, 2] = b*b
for i in range(3, len(e)):
out[:, i] = out[:, i-1] * b
return out
def f_bc_e_cheat(b, e):
return (b**e[:, None]).T
for params in [(100, 5000, 4, np.int32, np.int32),
(100, 5000, 8, np.int32, np.int32),
(100, 5000, 4, float, np.int32),
(100, 5000, 8, float, np.int32)]:
k = params[0]
dat = data(*params)
ref = f_vander(dat[0][0], dat[1][0])
print('rep={} n_b={} n_e={} b_tp={} e_tp={}'.format(*params))
for name, func in list(globals().items()):
if not name.startswith('f_') or not isinstance(func, types.FunctionType):
continue
try:
assert np.allclose(ref, func(dat[0][0], dat[1][0]))
if not func(dat[0][0], dat[1][0]).flags.c_contiguous:
print('cheat', end=' ')
print("{:16s}{:16.8f} ms".format(name[2:], np.min(repeat(
'f(b.pop(), e.pop())', setup='b, e = data(*p)', globals={'f':func, 'data':data, 'p':params}, number=k)) * 1000 / k))
except:
print("{:16s} apparently failed".format(name[2:]))
我正在计算 Vandermonde matrix for a fairly large 1D array. The natural and clean way to do this is using 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 的原因。
效率对我的实施很重要。因此，我们也欢迎更快的替代品！
广播求幂如何？
%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 在评论中指出的那样，浮点数组的广播求幂速度变慢了。
至于为什么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）。
出于兴趣，我找到了另一种计算范德蒙矩阵的方法，结果如下：
%timeit x[:, None] ** np.arange(N)
150 µs ± 230 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
它做同样的事情，但速度要慢得多。答案在于操作是广播的，但是效率低下。
另一方面，对于 float 个数组，这实际上最终表现最好。
这里有一些方法，其中一些方法（在我的电脑上）比目前发布的方法快很多。
我认为最重要的观察是，这在很大程度上取决于你想要多少学位。求幂（我认为这是小整数指数的特殊情况）仅对小指数范围有意义。指数越多，基于乘法的方法越好。
我想强调一个基于 multiply.accumulate 的方法 (ma)，它类似于 numpy 的内置方法但速度更快（并不是因为我在检查时略过 - nnc， numpy-no-checks 证明了这一点）。除了最小的指数范围外，它实际上对我来说是最快的。
由于我不明白的原因，据我所知，numpy 实现做了三件缓慢且不必要的事情：(1) 它制作了相当多的基本向量副本。 (2) 它使它们不连续。 (3) 它就地积累，我相信这会强制缓冲。
我想提的另一件事是，对于小范围的整数（out_e_1 本质上是 ma 的手动版本）最快的速度会慢两倍以上通过升级到更大的 dtype 的简单预防措施（safe_e_1 可以说有点用词不当）。
基于广播的方法称为 bc_*，其中 * 表示广播轴（b 表示基数，e 表示 exp）'cheat' 表示结果不连续。
计时（3 个中的最佳）：
rep=100 n_b=5000 n_e=4 b_tp=<class 'numpy.int32'> e_tp=<class 'numpy.int32'>
vander 0.16699657 ms
bc_b 0.09595204 ms
bc_e 0.07959786 ms
ma 0.10755240 ms
nnc 0.16459018 ms
out_e_1 0.02037535 ms
out_e_2 0.02656622 ms
safe_e_1 0.04652272 ms
safe_e_2 0.04081079 ms
cheat bc_e_cheat 0.04668466 ms
rep=100 n_b=5000 n_e=8 b_tp=<class 'numpy.int32'> e_tp=<class 'numpy.int32'>
vander 0.25086462 ms
bc_b apparently failed
bc_e apparently failed
ma 0.15843041 ms
nnc 0.24713077 ms
out_e_1 apparently failed
out_e_2 apparently failed
safe_e_1 0.15970622 ms
safe_e_2 0.19672418 ms
bc_e_cheat apparently failed
rep=100 n_b=5000 n_e=4 b_tp=<class 'float'> e_tp=<class 'numpy.int32'>
vander 0.16225773 ms
bc_b 0.53315020 ms
bc_e 0.56200830 ms
ma 0.07626799 ms
nnc 0.16059748 ms
out_e_1 0.03653416 ms
out_e_2 0.04043702 ms
safe_e_1 0.04060494 ms
safe_e_2 0.04104209 ms
cheat bc_e_cheat 0.52966076 ms
rep=100 n_b=5000 n_e=8 b_tp=<class 'float'> e_tp=<class 'numpy.int32'>
vander 0.24542852 ms
bc_b 2.03353578 ms
bc_e 2.04281270 ms
ma 0.11075758 ms
nnc 0.24212880 ms
out_e_1 0.14809043 ms
out_e_2 0.19261359 ms
safe_e_1 0.15206112 ms
safe_e_2 0.19308420 ms
cheat bc_e_cheat 1.99176601 ms
代码：
import numpy as np
import types
from timeit import repeat
prom={np.dtype(np.int32): np.dtype(np.int64), np.dtype(float): np.dtype(float)}
def RI(k, N, dt, top=100):
return np.random.randint(0, top if top else N, (k, N)).astype(dt)
def RA(k, N, dt, top=None):
return np.add.outer(np.zeros((k,), int), np.arange(N)%(top if top else N)).astype(dt)
def RU(k, N, dt, top=100):
return (np.random.random((k, N))*(top if top else N)).astype(dt)
def data(k, N_b, N_e, dt_b, dt_e, b_fun=RI, e_fun=RA):
b = list(b_fun(k, N_b, dt_b))
e = list(e_fun(k, N_e, dt_e))
return b, e
def f_vander(b, e):
return np.vander(b, len(e), increasing=True)
def f_bc_b(b, e):
return b[:, None]**e
def f_bc_e(b, e):
return np.ascontiguousarray((b**e[:, None]).T)
def f_ma(b, e):
out = np.empty((len(b), len(e)), prom[b.dtype])
out[:, 0] = 1
np.multiply.accumulate(np.broadcast_to(b, (len(e)-1, len(b))), axis=0, out=out[:, 1:].T)
return out
def f_nnc(b, e):
out = np.empty((len(b), len(e)), prom[b.dtype])
out[:, 0] = 1
out[:, 1:] = b[:, None]
np.multiply.accumulate(out[:, 1:], out=out[:, 1:], axis=1)
return out
def f_out_e_1(b, e):
out = np.empty((len(b), len(e)), b.dtype)
out[:, 0] = 1
out[:, 1] = b
out[:, 2] = c = b*b
for i in range(3, len(e)):
c*=b
out[:, i] = c
return out
def f_out_e_2(b, e):
out = np.empty((len(b), len(e)), b.dtype)
out[:, 0] = 1
out[:, 1] = b
out[:, 2] = b*b
for i in range(3, len(e)):
out[:, i] = out[:, i-1] * b
return out
def f_safe_e_1(b, e):
out = np.empty((len(b), len(e)), prom[b.dtype])
out[:, 0] = 1
out[:, 1] = b
out[:, 2] = c = (b*b).astype(prom[b.dtype])
for i in range(3, len(e)):
c*=b
out[:, i] = c
return out
def f_safe_e_2(b, e):
out = np.empty((len(b), len(e)), prom[b.dtype])
out[:, 0] = 1
out[:, 1] = b
out[:, 2] = b*b
for i in range(3, len(e)):
out[:, i] = out[:, i-1] * b
return out
def f_bc_e_cheat(b, e):
return (b**e[:, None]).T
for params in [(100, 5000, 4, np.int32, np.int32),
(100, 5000, 8, np.int32, np.int32),
(100, 5000, 4, float, np.int32),
(100, 5000, 8, float, np.int32)]:
k = params[0]
dat = data(*params)
ref = f_vander(dat[0][0], dat[1][0])
print('rep={} n_b={} n_e={} b_tp={} e_tp={}'.format(*params))
for name, func in list(globals().items()):
if not name.startswith('f_') or not isinstance(func, types.FunctionType):
continue
try:
assert np.allclose(ref, func(dat[0][0], dat[1][0]))
if not func(dat[0][0], dat[1][0]).flags.c_contiguous:
print('cheat', end=' ')
print("{:16s}{:16.8f} ms".format(name[2:], np.min(repeat(
'f(b.pop(), e.pop())', setup='b, e = data(*p)', globals={'f':func, 'data':data, 'p':params}, number=k)) * 1000 / k))
except:
print("{:16s} apparently failed".format(name[2:]))