如何避免具有多列的 numpy 数组的总和精度较低

Question

我一直假设 numpy uses a kind of pairwise-summation，这也确保了 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

可能 sum 只是天真地对所有值求和。一个指示是 16777216.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，如果这个起作用的话。

Answer 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)

Answer 2

此行为是由于 numpy 在减少操作期间访问内存的方式（"add" 只是一种特殊情况）以提高缓存的利用率。

对于某些情况（如上面的情况），可以在对性能没有太大影响的情况下强制执行成对求和。但总的来说，强制执行它会导致大量性能损失——使用双精度可能更容易，这在大多数情况下可以缓解上述问题。

可以看出，成对求和是对 "add" 运算的一种非常具体的优化，如果满足某些约束条件（稍后会详细介绍），则会完成此操作。

求和（以及许多其他减少操作）受内存带宽限制。如果我们沿着一个连续的轴求和，生活会很好：为索引 i 提取到缓存中的内存将直接重新用于索引 i+1, i+2,... 的计算而无需在使用之前从缓存中逐出。

情况不同，当求和不是沿着连续的轴时：要添加一个 float32 元素，16 个 float32 被提取到缓存中，但其中 15 个在使用之前被驱逐，并且必须再次被取回 - 真是浪费。

这就是 numpy 在这种情况下按行求和的原因：对第一行和第二行求和，然后将第三行添加到结果中，然后是第四行，依此类推。但是pairwise summation只是针对一维求和实现的，这里不能用

进行两两求和，当：

sum 在一维 numpy 数组上调用
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 个元素，开销太大（另请参阅解释的类似行为 here）。

它可以做得更好，例如通过仍然不精确 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 次那么高），而且精确，因为 sum 用于一次-dimensional numpy-arrays.

对于更多列，对于 numpy 和 einsum 方式，与连续大小写的差异要小得多：

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

但是 "precise" 技巧的性能非常糟糕，可能是因为延迟无法再通过计算隐藏：

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 实现的血淋淋的细节。

区别可以看FLOAT_add with defines from here的代码:

#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]是指向result/initial值的指针（与args[2]相同）
args[1]是输入数组
steps[0] 和 steps[2] 是 0，即指针指向标量。

然后可以使用成对求和（用IS_BINARY_REDUCE检查）。

FLOAT_add可用于将两个向量相加，在本例中：

args[0] 第一个输入数组
args[1] 第二个输入数组
args[2]输出数组
steps - 对于上述数组，从数组中的一个元素步进到另一个元素。

参数 @PW@ 1 仅用于求和 - 对于所有其他操作，不使用成对求和。

如何避免具有多列的 numpy 数组的总和精度较低

How to avoid less precise sum for numpy-arrays with multiple columns

python

numpy

floating-accuracy

ieee-754