Python: 求大数组的乘积时，如何减少浮点误差？

Question

假设我有一个大数组，里面有一堆浮点数，我需要找到乘积，同时尽可能少地因浮点错误而损失精度：

import numpy as np
randoms = np.random.uniform(0.5, 1.61, 10000)
print(randoms[0:10])

array([ 1.01422339,  0.65581167,  0.8154046 ,  1.49519379,  0.96114304,
    1.20167417,  0.93667198,  0.66899907,  1.26731008,  1.59689486])

一个可能不好的方法是遍历数组并迭代相乘。这显然会在每次乘法时产生错误，因此应尽可能避免：

product_1 = 1
for i in randoms:
    product_1 = product_1 * i
print(product_1)

64355009.758539267

下一个方法是使用numpy内置的prod函数，但是这个returns和上面的值一模一样，说明是这样的prod 实际上是在计算它：

product_2 = np.prod(randoms)
print(product_2)

64355009.758539267

print(product_1 == product_2)

True

第三种方法是计算每一项的对数，将它们相加，最后取幂。每个对数都是单独计算的，因此不会出现相同的误差复合，但对数过程和求幂过程本身都会引入一些误差。无论如何它都会产生不同的答案：

product_3 = np.exp(np.sum(np.log(randoms)))
print(product_3)

64355009.758538999

print(product_3 == product_1)

False

我知道在这个例子中我并没有失去那么多精度，但是对于我实际需要做的事情，复合错误最终会造成麻烦，足以让我考虑使用一个可以做到的包符号/任意精度计算。那么，哪种方法最好呢？还有其他我没有考虑过的方式吗？

Answer 1

我尝试了一些实验。代码在下面，但首先是一些评论。

可以通过将值转换为精确的有理数，精确计算乘积，然后执行最终转换为浮点数来精确计算结果。它可以通过 Python 中包含的 fractions 模块来完成，但它最终会变得非常慢。我使用 gmpy2 模块来进行更快的有理算术。

用于显示的二进制浮点值的格式有一些细微差别。 Python return 的最新版本将产生原始值的最短十进制字符串。 numpy 浮点数具有不同的格式。 gmpy2.mpfr 类型也是如此。 Decimal 显然使用了不同的格式规则。所以我总是将计算结果转换为 Python float.

除了Decimal类型的自定义十进制精度外，我还使用了gmpy2.mpfr，因为它支持自定义二进制精度。

程序输出几个值：

使用 53 位精度的顺序乘法（IEEE 64 位格式）。
使用有理算术的精确值。
使用精度为 28 位小数的小数。
使用具有用户指定精度的十进制。
使用具有用户指定精度的 mpfr。
使用递归乘法的方法来减少乘法的次数。

这是代码。您可以修改 Decimal 和 mpfr 精度并测试准确性。

import numpy as np
from gmpy2 import mpq, mpfr, get_context, round2
from decimal import Decimal, getcontext

randoms = np.random.uniform(0.5, 1.61, 10000)

# Sequential multiplication using 53-bit binary precision.

product_1 = 1
for i in randoms:
    product_1 = product_1 * i
print("53-bit binary:     ", float(product_1))

# Exact value by converting all floats to fractions and then a final
# conversion to float. Uses gmpy2 for speed.

product_2 = 1
for i in randoms:
    product_2 = product_2 * mpq(i)
print("exact using mpq:   ", float(mpfr(product_2, precision=53)))

# Decimal math with 28 decimal digits (~93 bits of precision.)

product_3 = 1
for i in randoms:
    product_3 = product_3 * Decimal(i)
print("Decimal(prec=28):  ", float(product_3))

# Choose your own decimal precision.

getcontext().prec=18
product_4 = 1
for i in randoms:
    product_4 = product_4 * Decimal(i)
print("Decimal(prec=%s):   %s" % (getcontext().prec, float(product_4)))

# Choose your own binary precision.

get_context().precision = 60
product_5 = 1
for i in randoms:
    product_5 = product_5 * mpfr(i)
print("mpfr(precision=%s): %s" % (get_context().precision, float(product_5)))

# Recursively multiply pairs of numbers together.

def rmult(d):
    if len(d) == 1:
        return d[0]
    # If the length is odd, extend with 1.
    if len(d) & 1:
        d.append(1)
    temp = []
    for i in range(len(d)//2):
        temp.append(d[2*i] * d[2*i+1])
    return rmult(temp)

print("recursive 53-bit:  ", float(rmult(list(randoms))))

作为粗略的准则，随着乘法次数的增加，中间精度将需要提高。有理算术将有效地为您提供无限的中间精度。

结果 100% 准确有多重要？

Python: 求大数组的乘积时，如何减少浮点误差？

Python: When finding the product of a large array, how best to reduce floating point error?

python

precision

numpy

floating-accuracy