缩放在numpy中在3D阵列上广播操作的时间
max*_*max 6 python numpy broadcasting numpy-ufunc
我试图在两个3D阵列上播放">"的简单操作.一个具有另一个维度(m,1,n)(1,m,n).如果我改变第三维(n)的值,我会天真地期望计算的速度将缩放为n.
然而,当我尝试明确地测量它时,我发现当将n从1增加到2时,计算时间增加约10倍,之后缩放是线性的.
当从n = 1到n = 2时,为什么计算时间会急剧增加？我假设它是numpy中的内存管理工件,但我正在寻找更具体的内容.
代码附在下面,附带结果图.
import numpy as np
import time
import matplotlib.pyplot as plt
def compute_time(n):
x, y = (np.random.uniform(size=(1, 1000, n)),
np.random.uniform(size=(1000, 1, n)))
t = time.time()
x > y
return time.time() - t
a = [
[
n, np.asarray([compute_time(n)
for _ in range(100)]).mean()
]
for n in range(1, 30, 1)
]
a = np.asarray(a)
plt.plot(a[:, 0], a[:, 1])
plt.xlabel('n')
plt.ylabel('time(ms)')
plt.show()
广播操作的时间图
我无法证明这一点,但我很确定这是由于一个简单的优化只能在n == 1时获得.
目前,numpy ufunc实现基于最内层循环的计算机生成代码,该代码映射到简单的C循环.封闭循环需要使用完全成熟的迭代器对象,该对象取决于有效负载,即最内层循环的大小和原子操作的成本可能是显着的开销.
现在,在n == 1时,问题本质上是2D(numpy足够聪明,可以检测到),最内层循环的大小为1000,因此迭代器对象的步长为1000.从n == 2向上,最里面的循环的大小为n,我们有100,000步迭代器对象,它们考虑了你正在观察的跳跃.
正如我所说,我不能证明它,但我可以使它看起来似乎合理:如果我们将变量维度移动到前面,那么最内层循环的常量大小为1000,外部循环在1000个迭代步骤中线性增长.事实上,这使得跳跃消失了.
码:
import numpy as np
import time
import matplotlib.pyplot as plt
def compute_time(n, axis=2):
xs, ys = [1, 10], [10, 1]
xs.insert(axis, n)
ys.insert(axis, n)
x, y = (np.random.uniform(size=xs),
np.random.uniform(size=ys))
t = time.perf_counter()
x > y
return time.perf_counter() - t
a = [
[
n,
np.asarray([compute_time(n) for _ in range(100)]).mean(),
np.asarray([compute_time(n, 0) for _ in range(100)]).mean()
]
for n in range(0, 10, 1)
]
a = np.asarray(a)
plt.plot(a[:, 0], a[:, 1:])
plt.xlabel('n')
plt.ylabel('time(ms)')
plt.show()
相关:https://stackoverflow.com/a/48257213/7207392
@保罗的理论是对的.在这个答案中,我使用perf和调试器深入研究以支持这一理论.
首先,让我们看一下运行时间的花费(参见run.py bellow的列表以获取确切的代码).
对于n=1我们看到以下内容:
Event count (approx.): 3388750000
Overhead Command Shared Object Symbol
34,04% python umath.cpython-36m-x86_64-linux-gnu.so [.] DOUBLE_less
32,71% python multiarray.cpython-36m-x86_64-linux-gnu.so [.] _aligned_strided_to_contig_size8_srcstride0
28,16% python libc-2.23.so [.] __memmove_ssse3_back
1,46% python multiarray.cpython-36m-x86_64-linux-gnu.so [.] PyArray_TransferNDimToStrided
相比n=2:
Event count (approx.): 28954250000
Overhead Command Shared Object Symbol
40,85% python libc-2.23.so [.] __memmove_ssse3_back
40,16% python multiarray.cpython-36m-x86_64-linux-gnu.so [.] PyArray_TransferNDimToStrided
8,61% python umath.cpython-36m-x86_64-linux-gnu.so [.] DOUBLE_less
8,41% python multiarray.cpython-36m-x86_64-linux-gnu.so [.] _contig_to_contig
对于n = 2,计算的事件多8.5倍,但仅为数据的两倍,因此我们需要解释减速因子4.
另一个重要的观察:运行时间由内存操作的主导n=2和(不太明显)也n=1(_aligned_strided_to_contig_size8_srcstride0是所有关于复制数据),他们超重比较成本- DOUBLE_less.
显然,PyArray_TransferNDimtoStrided这两种尺寸都需要,那为什么它在运行时间上的份额差别如此之大？
显示的自身时间PyArray_TransferNDimtoStrided不是复制所需的时间,而是开销:调整指针,以便在最后一个维度中可以通过以下方式复制stransfer:
PyArray_TransferNDimToStrided(npy_intp ndim,
....
/* A loop for dimensions 0 and 1 */
for (i = 0; i < shape1; ++i) {
if (shape0 >= count) {
stransfer(dst, dst_stride, src, src_stride0,
count, src_itemsize, data);
return 0;
}
else {
stransfer(dst, dst_stride, src, src_stride0,
shape0, src_itemsize, data);
}
count -= shape0;
src += src_stride1;
dst += shape0*dst_stride;
}
...
这些转移函数是_aligned_strided_to_contig_size8_srcstride0(参见下面列表中的生成代码)和_contig_to_contig:
_contig_to_contig用于n=2和传输2双(最后维度有2个值)的情况下,调整指针的开销非常高!_aligned_strided_to_contig_size8_srcstride0用于n=1并且每次调用传输1000个双精度(正如@Paul指出的那样,我们很快就会看到,numpy足够聪明地丢弃尺寸,这是1个元素的长度),调整指针的开销可以忽略不计.
顺便说一句,为了使用现代CPU的矢量化,使用这些函数而不是简单的for循环:在编译时已知步幅,编译器能够对代码进行矢量化(编译器通常无法对仅在运行时),因此numpy分析访问模式并分派给不同的预编译函数.
还有一个问题:如果我们的观察结果显示numpy的大小为1,那么numpy是否真的丢弃了最后一个维度？
使用debbuger很容易验证:
- ufunc通过迭代器访问数据,迭代器是在
iterator_loopvia中创建的NpyIter_AdvancedNew - 在
NpyIter_AdvancedNew,时,分析(和重新解释)维度npyiter_coalesce_axes
至于速度因子4被"丢失"比较时n=2到n=1:它没有特殊的意义,是对我的maschine只是一个随机值:从10 ^ 3至10 ^ 4更改矩阵的尺寸将进一步转向的优点(更少的开销)甚至更进一步 - n=1这导致我的机器失去了速度因子12.
run.py
import sys
import numpy as np
n=int(sys.argv[1])
x, y = (np.random.uniform(size=(1, 1000, n)),
np.random.uniform(size=(1000, 1, n)))
for _ in range(10000):
y<x
然后:
perf record python run.py 1
perf report
....
perf record python run.py 2
perf report
生成的来源_aligned_strided_to_contig_size8_srcstride0:
/*
* specialized copy and swap for source stride 0,
* interestingly unrolling here is like above is only marginally profitable for
* small types and detrimental for >= 8byte moves on x86
* but it profits from vectorization enabled with -O3
*/
#if (0 == 0) && 1
static NPY_GCC_OPT_3 void
_aligned_strided_to_contig_size8_srcstride0(char *dst,
npy_intp dst_stride,
char *src, npy_intp NPY_UNUSED(src_stride),
npy_intp N, npy_intp NPY_UNUSED(src_itemsize),
NpyAuxData *NPY_UNUSED(data))
{
#if 8 != 16
# if !(8 == 1 && 1)
npy_uint64 temp;
# endif
#else
npy_uint64 temp0, temp1;
#endif
if (N == 0) {
return;
}
#if 1 && 8 != 16
/* sanity check */
assert(npy_is_aligned(dst, _ALIGN(npy_uint64)));
assert(npy_is_aligned(src, _ALIGN(npy_uint64)));
#endif
#if 8 == 1 && 1
memset(dst, *src, N);
#else
# if 8 != 16
temp = _NPY_NOP8(*((npy_uint64 *)src));
# else
# if 0 == 0
temp0 = (*((npy_uint64 *)src));
temp1 = (*((npy_uint64 *)src + 1));
# elif 0 == 1
temp0 = _NPY_SWAP8(*((npy_uint64 *)src + 1));
temp1 = _NPY_SWAP8(*((npy_uint64 *)src));
# elif 0 == 2
temp0 = _NPY_SWAP8(*((npy_uint64 *)src));
temp1 = _NPY_SWAP8(*((npy_uint64 *)src + 1));
# endif
# endif
while (N > 0) {
# if 8 != 16
*((npy_uint64 *)dst) = temp;
# else
*((npy_uint64 *)dst) = temp0;
*((npy_uint64 *)dst + 1) = temp1;
# endif
# if 1
dst += 8;
# else
dst += dst_stride;
# endif
--N;
}
#endif/* @elsize == 1 && 1 -- else */
}
#endif/* (0 == 0) && 1 */