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

Question

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

6

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

然而，在对数组执行总和（或乘积）时，我遇到了精度问题，这取决于值的顺序。

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

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这样的函数中使用它们，但我无法提供容差...

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

- Thomas Leonard

你永远不能期望浮点数操作是精确的。你应该改变你的算法来适应一些误差。否则，在statistics模块中有一些更高级的求和方法，但现在这不是重点。特别是在numpy中，向量化应该是一种基本工具，你永远不能依赖算术运算的顺序。 - Andras Deak -- Слава Україні

你在哪里/如何使用np.in1d？对于所列代码，你可以使用np.isclose(s1,s2)。 - Divakar

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

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

不，不，解决方案通常不是限制您的数据为浮点数 :) 我的观点与Divakar的观点相同：不要使用精确测试，而要使用接近度测试。但是您似乎已经意识到了这一点，这就是为什么我建议重新构造您的算法以避免精确测试。 - Andras Deak -- Слава Україні

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Divakar · Accepted Answer

对于这段代码，您可以使用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，我们需要沿一个轴执行任意操作。因此，使用broadcasting的带有公差的np.in1d可以这样实现 -

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

正如OP在评论中建议的那样，我们可以在将浮点数放大后四舍五入，这样可以节省内存，特别适用于处理大型数组。因此，修改后的版本应该像这样 -

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应该有所帮助。