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

Posted 2023-02-23

【中文标题】Python：求大数组的乘积时，如何最好地减少浮点错误？

【英文标题】：Python: When finding the product of a large array, how best to reduce floating point error?

【发布时间】：2018-08-03 20:39:31

【问题描述】：

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

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

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

【问题讨论】：

您还没有考虑过的其他方式之一是NumPy 的cumprod，它代表累积产品。这就是你基本上正在做的事情。只需将最后一个元素[-1] 设为product_1 = np.cumprod(randoms)[-1]。你可以比较你的答案

看看 Python 的小数和分数库

我删除了我的答案，因为print(product_3 == product_1) 给了我False。 print(product_4 == product_2) 和 print(product_4 == product_1) 产生 True，其中 product_4 是使用 np.cumprod() 的结果。似乎 log 和 exp 是导致精度变化的罪魁祸首。

即使您使用decimal 并执行from decimal import * 然后getcontext().prec = 20，它仍然会导致False。现在取决于您想要什么精度。

是的，但我可以以不同的精度重复计算，并询问输出以查看它们开始发散的数字。

【参考方案1】：

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

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

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

除了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% 准确有多重要？