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

Question

我一直假设numpy 使用一种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

此行为是由于 numpy 在归约操作期间访问内存的方式（“添加”只是一种特殊情况）以提高缓存的利用率。

\n\n

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

\n\n\n\n

成对求和可以看作是“添加”操作的一种非常具体的优化，如果满足一些约束（稍后会详细介绍），就会完成这种优化。

\n\n

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

\n\n

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

\n\n

这就是为什么 numpy 在这种情况下按行求和：对第一行和第二行求和，然后将第三行添加到结果中，然后是第四行，依此类推。然而，成对求和仅针对一维求和实现，不能在这里使用。

\n\n

当满足以下条件时，执行成对求和：

\n\n

\n
• sum在一维 numpy 数组上调用
\n
• sum沿连续轴调用
\n

\n\n

numpy 还没有提供一种执行成对求和而不会对性能产生重大负面影响的方法。

\n\n

我的看法是：目标应该是沿连续轴执行求和，这不仅更精确，而且可能更快：

\n\n

A=np.ones((N,2), dtype=np.float32, order="C") #non-contiguous\n%timeit A.sum(axis=0)\n# 326 ms \xc2\xb1 9.17 ms\n\nB=np.ones((N,2), dtype=np.float32, order="F") # contiguous\n%timeit B.sum(axis=0)\n# 15.6 ms \xc2\xb1 898 \xc2\xb5s \n

\n\n

在这种特殊情况下，一行中只有 2 个元素，开销太大（另请参阅此处解释的类似行为）。

\n\n

它可以做得更好，例如通过仍然不精确einsum：

\n\n

%timeit np.einsum("i...->...", A)\n# 74.5 ms \xc2\xb1 1.47 ms \nnp.einsum("i...->...", A)\n# array([16777216.,  16777216.], dtype=float32)\n

\n\n

甚至：

\n\n

%timeit np.array([A[:,0].sum(), A[:,1].sum()], dtype=np.float32)\n# 17.8 ms \xc2\xb1 333 \xc2\xb5s \nnp.array([A[:,0].sum(), A[:,1].sum()], dtype=np.float32)\n# array([17000000., 17000000.], dtype=float32)\n

\n\n

它不仅几乎与连续版本一样快（加载内存两次的代价并不像加载内存 16 次那么高），而且也很精确，因为sum用于一维 numpy 数组。

\n\n

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

\n\n

B=np.ones((N,16), dtype=np.float32, order="F")\n%timeit B.sum(axis=0)\n# 121 ms \xc2\xb1 3.66 ms \n\nA=np.ones((N,16), dtype=np.float32, order="C")\n%timeit A.sum(axis=0)\n# 457 ms \xc2\xb1 12.1 ms \n\n%timeit np.einsum("i...->...", A)\n# 139 ms \xc2\xb1 651 \xc2\xb5s per loop\n

\n\n

但对于“精确”技巧来说，性能非常糟糕，可能是因为计算无法再隐藏延迟：

\n\n

def do(A):\n    N=A.shape[1]\n    res=np.zeros(N, dtype=np.float32)\n    for i in range(N):\n        res[i]=A[:,i].sum()\n    return res\n%timeit do(A)\n# 1.39 s \xc2\xb1 47.8 ms\n

\n\n\n\n

以下是 numpy 实现的详细细节。

\n\n

差异可以在代码中看到FLOAT_add可以在此处的定义：

\n\n

#define IS_BINARY_REDUCE ((args[0] == args[2])\\\n    && (steps[0] == steps[2])\\\n    && (steps[0] == 0))\n\n#define BINARY_REDUCE_LOOP(TYPE)\\\n   char *iop1 = args[0]; \\\n   TYPE io1 = *(TYPE *)iop1; \\\n\n/** (ip1, ip2) -> (op1) */\n#define BINARY_LOOP\\\n    char *ip1 = args[0], *ip2 = args[1], *op1 = args[2];\\\n    npy_intp is1 = steps[0], is2 = steps[1], os1 = steps[2];\\\n    npy_intp n = dimensions[0];\\\n    npy_intp i;\\\n    for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1)\n\n/**begin repeat\n* Float types\n*  #type = npy_float, npy_double, npy_longdouble#\n*  #TYPE = FLOAT, DOUBLE, LONGDOUBLE#\n*  #c = f, , l#\n*  #C = F, , L#\n*/\n\n/**begin repeat1\n * Arithmetic\n * # kind = add, subtract, multiply, divide#\n * # OP = +, -, *, /#\n * # PW = 1, 0, 0, 0#\n */\nNPY_NO_EXPORT void\n@TYPE@_@kind@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))\n{\n    if (IS_BINARY_REDUCE) {\n#if @PW@\n        @type@ * iop1 = (@type@ *)args[0];\n        npy_intp n = dimensions[0];\n\n        *iop1 @OP@= pairwise_sum_@TYPE@(args[1], n, steps[1]);\n#else\n        BINARY_REDUCE_LOOP(@type@) {\n            io1 @OP@= *(@type@ *)ip2;\n        }\n        *((@type@ *)iop1) = io1;\n#endif\n    }\n    else if (!run_binary_simd_@kind@_@TYPE@(args, dimensions, steps)) {\n        BINARY_LOOP {\n            const @type@ in1 = *(@type@ *)ip1;\n            const @type@ in2 = *(@type@ *)ip2;\n            *((@type@ *)op1) = in1 @OP@ in2;\n        }\n    }\n}\n

\n\n

生成后如下所示：

\n\n

NPY_NO_EXPORT void\nFLOAT_add(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))\n{\n    if (IS_BINARY_REDUCE) {\n#if 1\n        npy_float * iop1 = (npy_float *)args[0];\n        npy_intp n = dimensions[0];\n\n        *iop1 += pairwise_sum_FLOAT((npy_float *)args[1], n,\n                                        steps[1] / (npy_intp)sizeof(npy_float));\n#else\n        BINARY_REDUCE_LOOP(npy_float) {\n            io1 += *(npy_float *)ip2;\n        }\n        *((npy_float *)iop1) = io1;\n#endif\n    }\n    else if (!run_binary_simd_add_FLOAT(args, dimensions, steps)) {\n        BINARY_LOOP {\n            const npy_float in1 = *(npy_float *)ip1;\n            const npy_float in2 = *(npy_float *)ip2;\n            *((npy_float *)op1) = in1 + in2;\n        }\n    }\n}\n

\n\n

FLOAT_add可用于一维缩减，在本例中：

\n\n

\n
• args[0]是指向结果/初始值的指针（与args[2]）
\n
• args[1]是输入数组
\n
• steps[0]并且steps[2]是0，即指针指向标量。
\n

\n\n

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

\n\n

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

\n\n

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

\n\n

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

\n