使用浮点 NumPy 数组进行比较和相关操作

Posted 2023-02-23

【中文标题】使用浮点 NumPy 数组进行比较和相关操作

【英文标题】：Working with floating point NumPy arrays for comparison and related operations

【发布时间】：2017-02-07 00:39:19

【问题描述】：

我有一个随机浮点数组，我需要将它与另一个具有不同顺序的相同值的数组进行比较。就此而言，我使用总和、乘积（以及其他组合，具体取决于表格的维度，因此需要的方程式数量）。

不过，当我根据值的顺序对数组执行求和（或乘积）时遇到了精度问题。

这是一个简单的独立示例来说明这个问题：

import numpy as np

n = 10
m = 4

tag = np.random.rand(n, m)

s1 = np.sum(tag, axis=1)
s2 = np.sum(tag[:, ::-1], axis=1)

# print the number of times s1 is not equal to s2 (should be 0)
print np.nonzero(s1 != s2)[0].shape[0]

如果您执行此代码，它有时会告诉您s1 和s2 不相等，并且差异取决于计算机精度。

问题是我需要在 np.in1d 这样的函数中使用那些我无法真正容忍的函数......

有没有办法避免这个问题？

【问题讨论】：

你永远不能期望浮点运算是精确的。您应该更改算法以适应一些错误。否则，statistics 模块中有一些更奇特的总和，但这不是现在的重点。特别是对于 numpy，向量化应该是一个基本工具，你永远不能依赖算术运算的顺序。

您如何/在哪里使用np.in1d？对于列出的代码，您可以使用np.isclose(s1,s2)。

@Divakar 我没有在示例中使用它，但在我的实际算法中，我会使用np.in1d(s1, s2) 和通过其他操作（如产品等）获得的其他等效数组...

@Andras Deak 我想我需要重新考虑我的标签的选择......正如你所指出的那样，随机浮动可能不是一个好主意，但它很方便，因为我正在这样做非常大的数组，我执行像tag**3 这样的操作，我担心如果使用整数会导致溢出......我愿意接受建议（我的标签数组必须有非重复值）......

不，不，解决方案通常不是限制您的数据来自浮点数:) 我的观点与 Divakar 的观点相同：不要使用精确测试，使用接近测试。但是您似乎已经意识到这一点，这就是为什么我只建议重构您的算法以避免精确测试。

【参考方案1】：

对于列出的代码，您可以使用np.isclose，也可以指定公差值。

使用提供的示例，让我们看看如何使用它 -

In [201]: n = 10
     ...: m = 4
     ...: 
     ...: tag = np.random.rand(n, m)
     ...: 
     ...: s1 = np.sum(tag, axis=1)
     ...: s2 = np.sum(tag[:, ::-1], axis=1)
     ...: 

In [202]: np.nonzero(s1 != s2)[0].shape[0]
Out[202]: 4

In [203]: (~np.isclose(s1,s2)).sum() # So, all matches!
Out[203]: 0

要在其他情况下使用公差值，我们需要根据具体情况进行工作。因此，假设对于像np.in1d 中涉及元素比较的实现，我们可以引入broadcasting 对第一个输入中的所有元素与第二个输入中的所有元素进行元素相等性检查。然后，我们使用np.abs 得到“接近因子”，最后与输入容差进行比较来决定匹配。根据需要模拟np.in1d，我们沿着其中一个轴进行任何操作。因此，np.in1d 使用 broadcasting 的容差可以像这样实现 -

def in1d_with_tolerance(A,B,tol=1e-05):
    return (np.abs(A[:,None] - B) < tol).any(1)

正如 OP 在 cmets 中所建议的，我们还可以在将浮点数按比例放大后对其进行舍入，这应该是内存效率高的，因为这是处理大型数组所需要的。所以，修改后的版本会是这样 -

def in1d_with_tolerance_v2(A,B,tol=1e-05):
    S = round(1/tol)
    return np.in1d(np.around(A*S).astype(int),np.around(B*S).astype(int))

示例运行 -

In [372]: A = np.random.rand(5)
     ...: B = np.random.rand(7)
     ...: B[3] = A[1] + 0.0000008
     ...: B[6] = A[4] - 0.0000007
     ...: 

In [373]: np.in1d(A,B) # Not the result we want!
Out[373]: array([False, False, False, False, False], dtype=bool)

In [374]: in1d_with_tolerance(A,B)
Out[374]: array([False,  True, False, False,  True], dtype=bool)

In [375]: in1d_with_tolerance_v2(A,B)
Out[375]: array([False,  True, False, False,  True], dtype=bool)

最后，关于如何使其适用于其他实现和用例 - 这将取决于实现本身。但在大多数情况下，np.isclose 和 broadcasting 应该会有所帮助。

【讨论】：

这种具有容差的 in1d 实现使用了更多内存（对于大数组）......所以我做了一个 np.around 我想要比较的数组与接近的小数位数浮点精度，它似乎可以完成这项工作。谢谢。

@thomleo 好点！将其添加到帖子中，希望没关系。