Python:在求大数组的乘积时，如何最好地减少浮点错误？

Question

假设我有一个大数组，里面有一堆 float ，我需要找到乘积，同时尽可能减少浮点错误的精度:

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 函数，但是这会返回与上面完全相同的值，表明这就是 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

我尝试了一些实验。代码如下，但首先是一些注释。

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

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

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

程序输出几个值:

• 使用 53 位精度的顺序乘法(IEEE 64 位格式)。
• 使用有理算术得出精确值。
• 使用精度为 28 位小数的 Decimal。
• 使用具有用户指定精度的 Decimal。
• 按照用户指定的精度使用 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% 准确有多重要？