提问人：ead 提问时间：4/4/2019 更新时间：4/8/2019 访问量：1085
如何避免具有多列的 numpy 数组的精确求和
How to avoid less precise sum for numpy-arrays with multiple columns
问：
我一直认为，numpy 使用一种成对求和，这也确保了 - 操作的高精度：float32
import numpy as np
N=17*10**6 # float32-precision no longer enough to hold the whole sum
print(np.ones((N,1),dtype=np.float32).sum(axis=0))
# [17000000.], kind of expected
但是，如果矩阵具有多个列，则似乎使用了不同的算法：
print(np.ones((N,2),dtype=np.float32).sum(axis=0))
# [16777216. 16777216.] the error is just to big
print(np.ones((2*N,2),dtype=np.float32).sum(axis=0))
# [16777216. 16777216.] error is bigger
可能只是天真地将所有值相加。一个指示是，例如：sum16777216.f+1.0f=16777216.f
one = np.array([1.], np.float32)
print(np.array([16777215.], np.float32)+one) # 16777216.
print(np.array([16777216.], np.float32)+one) # 16777216. as well
为什么 numpy 不对多列使用成对求和，是否可以强制 numpy 也对多列使用成对求和？
我的numpy版本是1.14.2，如果这起作用的话。
答：
2赞
Stop harming Monica
4/4/2019
#1
我真的没有解释，但它似乎与内存布局有关。使用 fortran 顺序而不是默认的 C 顺序，我得到了所需的输出。
>>> np.ones((N,2),dtype=np.float32, order='C').sum(axis=0)
array([16777216., 16777216.], dtype=float32)
>>> np.ones((N,2),dtype=np.float32, order='F').sum(axis=0)
array([17000000., 17000000.], dtype=float32)
评论
0赞
ead
4/4/2019
有趣的是，再次导致错误的结果......np.ones((2,N),dtype=np.float32, order='F').sum(axis=1)
0赞
Stop harming Monica
4/4/2019
@ead我所期望的那样，对连续数据求和似乎更准确。我不知道为什么会这样。
4赞
Paul Panzer
4/4/2019
这是因为沿连续轴使用成对求和，而不是沿跨步轴使用：github.com/numpy/numpy/issues/8786#issuecomment-286864207
3赞
ead
4/7/2019
#2
此行为是由于 numpy 在 reduce-operation 期间访问内存的方式（“add”只是一个特例）以提高缓存的利用率。
对于某些情况（如上所述），可以在不对性能产生重大影响的情况下强制执行成对求和。但一般来说，强制执行它会导致巨大的性能损失 - 使用双精度可能更容易，这在大多数情况下可以缓解上述问题。
成对求和可以看作是对“加法”运算的非常具体的优化，如果满足一些约束（稍后会详细介绍），就会完成。
求和（以及许多其他 reduce-operations）受内存带宽限制。如果我们沿着一个连续的轴求和，寿命是好的：提取到索引缓存中的内存将直接重用用于索引的计算，,...在使用之前，不会从缓存中逐出。ii+1i+2
当求和不是沿着连续的轴时，情况就不同了：要添加一个 float32 元素，16-float32 被提取到缓存中，但其中 15 个在可以使用之前就被逐出，并且必须再次提取 - 真是浪费。
这就是为什么 numpy 在这种情况下逐行求和的原因：将第一行和第二行相加，然后将第三行添加到结果中，然后是第四行，依此类推。但是，成对求和仅用于一维求和，不能在此处使用。
在以下情况下执行成对求和：
sum在一维 numpy-array 上调用sum沿连续轴调用
Numpy（还？）没有提供一种在不对性能产生重大负面影响的情况下强制执行成对求和的方法。
我从中得到的启示是：目标应该是沿连续轴执行求和，这不仅更精确，而且可能更快：
A=np.ones((N,2), dtype=np.float32, order="C") #non-contiguous
%timeit A.sum(axis=0)
# 326 ms ± 9.17 ms
B=np.ones((N,2), dtype=np.float32, order="F") # contiguous
%timeit B.sum(axis=0)
# 15.6 ms ± 898 µs
在这种特殊情况下，如果连续只有 2 个元素，开销就太大了（另请参阅此处解释的类似行为）。
它可以做得更好，例如通过仍然不精确：einsum
%timeit np.einsum("i...->...", A)
# 74.5 ms ± 1.47 ms
np.einsum("i...->...", A)
# array([16777216., 16777216.], dtype=float32)
甚至：
%timeit np.array([A[:,0].sum(), A[:,1].sum()], dtype=np.float32)
# 17.8 ms ± 333 µs
np.array([A[:,0].sum(), A[:,1].sum()], dtype=np.float32)
# array([17000000., 17000000.], dtype=float32)
这不仅几乎与连续版本一样快（加载内存两次的惩罚不如加载内存 16 次高），而且很精确，因为它用于一维 numpy-arrays。sum
对于更多的列，numpy 和 einsum-ways 与连续大小写的差异要小得多：
B=np.ones((N,16), dtype=np.float32, order="F")
%timeit B.sum(axis=0)
# 121 ms ± 3.66 ms
A=np.ones((N,16), dtype=np.float32, order="C")
%timeit A.sum(axis=0)
# 457 ms ± 12.1 ms
%timeit np.einsum("i...->...", A)
# 139 ms ± 651 µs per loop
但是对于“精确”技巧来说，性能非常糟糕，可能是因为延迟不能再通过计算来隐藏：
def do(A):
N=A.shape[1]
res=np.zeros(N, dtype=np.float32)
for i in range(N):
res[i]=A[:,i].sum()
return res
%timeit do(A)
# 1.39 s ± 47.8 ms
以下是numpy实现的血腥细节。
#define IS_BINARY_REDUCE ((args[0] == args[2])\
&& (steps[0] == steps[2])\
&& (steps[0] == 0))
#define BINARY_REDUCE_LOOP(TYPE)\
char *iop1 = args[0]; \
TYPE io1 = *(TYPE *)iop1; \
/** (ip1, ip2) -> (op1) */
#define BINARY_LOOP\
char *ip1 = args[0], *ip2 = args[1], *op1 = args[2];\
npy_intp is1 = steps[0], is2 = steps[1], os1 = steps[2];\
npy_intp n = dimensions[0];\
npy_intp i;\
for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1)
/**begin repeat
* Float types
* #type = npy_float, npy_double, npy_longdouble#
* #TYPE = FLOAT, DOUBLE, LONGDOUBLE#
* #c = f, , l#
* #C = F, , L#
*/
/**begin repeat1
* Arithmetic
* # kind = add, subtract, multiply, divide#
* # OP = +, -, *, /#
* # PW = 1, 0, 0, 0#
*/
NPY_NO_EXPORT void
@TYPE@_@kind@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
{
if (IS_BINARY_REDUCE) {
#if @PW@
@type@ * iop1 = (@type@ *)args[0];
npy_intp n = dimensions[0];
*iop1 @OP@= pairwise_sum_@TYPE@(args[1], n, steps[1]);
#else
BINARY_REDUCE_LOOP(@type@) {
io1 @OP@= *(@type@ *)ip2;
}
*((@type@ *)iop1) = io1;
#endif
}
else if (!run_binary_simd_@kind@_@TYPE@(args, dimensions, steps)) {
BINARY_LOOP {
const @type@ in1 = *(@type@ *)ip1;
const @type@ in2 = *(@type@ *)ip2;
*((@type@ *)op1) = in1 @OP@ in2;
}
}
}
生成后如下所示：
NPY_NO_EXPORT void
FLOAT_add(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
{
if (IS_BINARY_REDUCE) {
#if 1
npy_float * iop1 = (npy_float *)args[0];
npy_intp n = dimensions[0];
*iop1 += pairwise_sum_FLOAT((npy_float *)args[1], n,
steps[1] / (npy_intp)sizeof(npy_float));
#else
BINARY_REDUCE_LOOP(npy_float) {
io1 += *(npy_float *)ip2;
}
*((npy_float *)iop1) = io1;
#endif
}
else if (!run_binary_simd_add_FLOAT(args, dimensions, steps)) {
BINARY_LOOP {
const npy_float in1 = *(npy_float *)ip1;
const npy_float in2 = *(npy_float *)ip2;
*((npy_float *)op1) = in1 + in2;
}
}
}
FLOAT_add可用于一维还原，在本例中：
args[0]是指向结果/初始值的指针（与args[2])args[1]是输入数组steps[0]和 are ，即指针指向标量。steps[2]0
然后可以使用成对求和（用 勾选 ）。IS_BINARY_REDUCE
FLOAT_add可用于添加两个向量，在本例中：
args[0]第一个输入数组args[1]第二个输入数组args[2]输出数组steps- 对于上述数组，从数组中的一个元素步进到另一个元素。
参数仅用于求和 - 对于所有其他操作，不使用成对求和。@PW@1
评论
print(np.ones((N,2),dtype=np.float32).sum() / 2)17000000.0