numpy的power和**在某些值上的区别

Question

我有一个numpy数组，其中的 f**2 和 f[i]**2 的值不一样，但这只发生在某个特定的值上。

import numpy as np
np.set_printoptions(precision = 16)

f = np.array([ -40709.6555510835, -40708.6555510835, -33467.081758611654, -27653.379955714125])

f2 = f**2
# f2 = np.power(f,2)

print("outside loop", np.abs(f[1]**2 - f2[1]), np.abs(1.0 - f2[1] / f[1]**2), f.dtype, f[1].dtype, f2.dtype, f2[1].dtype)

for i, val in enumerate(f):        
    print("inside loop", i, np.abs(val**2 - f2[i]), np.abs(1.0 - f2[i] / val**2), val.dtype, f2.dtype, f2[i].dtype)

输出结果是：

outside loop 2.384185791015625e-07 2.220446049250313e-16 float64 float64 float64 float64
inside loop 0 0.0 0.0 float64 float64 float64
inside loop 1 2.384185791015625e-07 2.220446049250313e-16 float64 float64 float64
inside loop 2 0.0 0.0 float64 float64 float64
inside loop 3 0.0 0.0 float64 float64 float64

我注意到这是一个相对误差，大小大约是epsilon。

当我在定义 f2 时使用 np.power 而不是 ** 时，这个问题就消失了。不过，我还是想知道为什么 f[i]**2 和 f**2 的第 i 个值不一样（即使这只在某些 f 的值上发生）。

我使用的是python 3.10.6和最新的numpy 1.26.4。

编辑：

这个根本问题可以在以下内容中找到：

import numpy as np
f = np.array([-40708.6555510835])
print((f[0])**2 - (f**2)[0])

它显示的值是

-2.384185791015625e-07

我想知道为什么那个特定的数字会有这个特定的问题。如果你想确认，或者尝试不同的 f 值，可以查看这个 演示。

Answer 1

我之前遇到过非常类似的问题，原因在于Python在x86架构上使用的是传统的80位浮点寄存器，而NumPy则使用了SIMD，这种方式是64位浮点。

这通常表现为NumPy中的浮点运算精度较低，但这并不是说它们是错误的。传统的x86（技术上讲是x87）因为有更宽的寄存器，所以可以提供更高的精度，但在使用SIMD指令时，这种更高的精度就无法实现了，很多其他架构也是如此。

Answer 2

这可能是numpy的一个问题，可能是在这个链接中引入的 https://github.com/numpy/numpy/pull/7496.

我是在查找numpy的更新日志时，搜索“pow”这个词发现的。看起来python和numpy都使用标准库中的pow函数来处理浮点数。所以我在寻找numpy在处理标量值和数组时对**运算符的不同处理方式。

我不知道这是否是导致这种行为的真正原因，但这是我目前找到的唯一一个numpy没有使用标准库中的pow()函数的地方。这让我发现了一个测试，似乎能重现你所看到的行为。

我在我的i7 iMac和M1 MBP上看到的结果完全相同。

>>> x = -40708.6555510835
>>> x**2
1657194636.7767615
>>> x * x
1657194636.7767618
>>> x**2 - x*x
-2.384185791015625e-07

所以如果我没错的话，f[0]**2是使用数学库的pow函数，而(f**2)[0]则跳过了pow，直接把值相乘。

Answer 3

结果不同是因为 f**2 调用了 numpy.square，而 f[0]**2 和 numpy.power(f, 2) 则调用了 numpy.power。

---

numpy.ndarray.__pow__ 是用 C 语言写的。它的样子是这样的：

static PyObject *
array_power(PyObject *a1, PyObject *o2, PyObject *modulo)
{
    PyObject *value = NULL;

    if (modulo != Py_None) {
        /* modular exponentiation is not implemented (gh-8804) */
        Py_INCREF(Py_NotImplemented);
        return Py_NotImplemented;
    }

    BINOP_GIVE_UP_IF_NEEDED(a1, o2, nb_power, array_power);
    if (fast_scalar_power(a1, o2, 0, &value) != 0) {
        value = PyArray_GenericBinaryFunction(a1, o2, n_ops.power);
    }
    return value;
}

这行 value = PyArray_GenericBinaryFunction(a1, o2, n_ops.power); 是在调用一个 Python 函数，指向 numpy.power 这个特殊的函数对象，但它首先会尝试一个叫 fast_scalar_power 的函数。这个函数会尝试优化一些常见的幂运算，比如说 2 的幂。

对于 f**2 这个操作，fast_scalar_power 检测到指数是 2，然后 把这个操作转交给 numpy.square：

else if (exponent ==  2.0) {
    fastop = n_ops.square;
}

而对于 numpy.power(f, 2)，这显然是直接调用 numpy.power。numpy.power 不会经过 fast_scalar_power，也没有对指数为 2 的情况进行特别处理。（不过，具体的实现可能会对 2 有特别的处理，这取决于它使用的底层实现。）

对于标量，我认为 numpy.float64.__pow__ 实际上只是调用了 array_power：

static PyObject *
gentype_power(PyObject *m1, PyObject *m2, PyObject *modulo)
{
    if (modulo != Py_None) {
        /* modular exponentiation is not implemented (gh-8804) */
        Py_INCREF(Py_NotImplemented);
        return Py_NotImplemented;
    }

    BINOP_GIVE_UP_IF_NEEDED(m1, m2, nb_power, gentype_power);
    return PyArray_Type.tp_as_number->nb_power(m1, m2, Py_None);
}

所以它会进入 fast_scalar_power，但是 在 fast_scalar_power 中的第一个检查之一是

if (PyArray_Check(o1) &&

一个 numpy.float64 的实例并没有通过 PyArray_Check，这个检查是用来确认对象的类型是否正好是 numpy.ndarray。因此，这个标量会走一般的 numpy.power 的处理路径。