numpy 中 C 与 F 有序数组的数组总和沿轴的不同原因
what causes different in array sum along axis for C versus F ordered arrays in numpy
我很好奇是否有人能解释 numpy
中 C 与 Fortran 有序数组的这种特殊处理方式究竟是什么导致了差异。请看下面的代码:
system:
Ubuntu 18.10
Miniconda python 3.7.1
numpy 1.15.4
def test_array_sum_function(arr):
idx=0
val1 = arr[idx, :].sum()
val2 = arr.sum(axis=(1))[idx]
print('axis sums:', val1)
print(' ', val2)
print(' equal:', val1 == val2)
print('total sum:', arr.sum())
n = 2_000_000
np.random.seed(42)
rnd = np.random.random(n)
print('Fortran order:')
arrF = np.zeros((2, n), order='F')
arrF[0, :] = rnd
test_array_sum_function(arrF)
print('\nC order:')
arrC = np.zeros((2, n), order='C')
arrC[0, :] = rnd
test_array_sum_function(arrC)
打印:
Fortran order:
axis sums: 999813.1414744433
999813.1414744079
equal: False
total sum: 999813.1414744424
C order:
axis sums: 999813.1414744433
999813.1414744433
equal: True
total sum: 999813.1414744433
浮点数学不一定是 associative,即 (a+b)+c != a+(b+c)
。
由于您是沿不同的轴添加,因此操作顺序不同,这会影响最终结果。作为一个简单的例子,考虑总和为 1 的矩阵。
a = np.array([[1e100, 1], [-1e100, 0]])
print(a.sum()) # returns 0, the incorrect result
af = np.asfortranarray(a)
print(af.sum()) # prints 1
(有趣的是,a.T.sum()
和 aT = a.T; aT.sum()
一样仍然给出 0,所以我不确定这在后端是如何实现的)
C 命令使用操作序列(从左到右)1e100 + 1 + (-1e100) + 0
,而 Fortran 命令使用 1e100 + (-1e100) + 1 + 0
。问题是 (1e100+1) == 1e100
因为浮点数没有足够的精度来表示那个小差异,所以 1
丢失了。
一般来说,不要对浮点数进行相等性测试,而是使用小的 epsilon(if abs(float1 - float2) < 0.00001
或 np.isclose
)进行比较。如果您需要任意浮点精度,请使用 Decimal
库或定点表示和 int
。
这几乎可以肯定是 numpy 有时使用 pairwise summation and sometimes not 的结果。
让我们构建一个诊断阵列:
eps = (np.nextafter(1.0, 2)-1.0) / 2
1+eps+eps+eps
# 1.0
(1+eps)+(eps+eps)
# 1.0000000000000002
X = np.full((32, 32), eps)
X[0, 0] = 1
X.sum(0)[0]
# 1.0
X.sum(1)[0]
# 1.000000000000003
X[:, 0].sum()
# 1.000000000000003
这强烈表明一维数组和连续轴使用成对求和,而多维数组中的跨轴不使用。
请注意,要看到效果,数组必须足够大,否则 numpy 会退回到普通求和。
我很好奇是否有人能解释 numpy
中 C 与 Fortran 有序数组的这种特殊处理方式究竟是什么导致了差异。请看下面的代码:
system:
Ubuntu 18.10
Miniconda python 3.7.1
numpy 1.15.4
def test_array_sum_function(arr):
idx=0
val1 = arr[idx, :].sum()
val2 = arr.sum(axis=(1))[idx]
print('axis sums:', val1)
print(' ', val2)
print(' equal:', val1 == val2)
print('total sum:', arr.sum())
n = 2_000_000
np.random.seed(42)
rnd = np.random.random(n)
print('Fortran order:')
arrF = np.zeros((2, n), order='F')
arrF[0, :] = rnd
test_array_sum_function(arrF)
print('\nC order:')
arrC = np.zeros((2, n), order='C')
arrC[0, :] = rnd
test_array_sum_function(arrC)
打印:
Fortran order:
axis sums: 999813.1414744433
999813.1414744079
equal: False
total sum: 999813.1414744424
C order:
axis sums: 999813.1414744433
999813.1414744433
equal: True
total sum: 999813.1414744433
浮点数学不一定是 associative,即 (a+b)+c != a+(b+c)
。
由于您是沿不同的轴添加,因此操作顺序不同,这会影响最终结果。作为一个简单的例子,考虑总和为 1 的矩阵。
a = np.array([[1e100, 1], [-1e100, 0]])
print(a.sum()) # returns 0, the incorrect result
af = np.asfortranarray(a)
print(af.sum()) # prints 1
(有趣的是,a.T.sum()
和 aT = a.T; aT.sum()
一样仍然给出 0,所以我不确定这在后端是如何实现的)
C 命令使用操作序列(从左到右)1e100 + 1 + (-1e100) + 0
,而 Fortran 命令使用 1e100 + (-1e100) + 1 + 0
。问题是 (1e100+1) == 1e100
因为浮点数没有足够的精度来表示那个小差异,所以 1
丢失了。
一般来说,不要对浮点数进行相等性测试,而是使用小的 epsilon(if abs(float1 - float2) < 0.00001
或 np.isclose
)进行比较。如果您需要任意浮点精度,请使用 Decimal
库或定点表示和 int
。
这几乎可以肯定是 numpy 有时使用 pairwise summation and sometimes not 的结果。
让我们构建一个诊断阵列:
eps = (np.nextafter(1.0, 2)-1.0) / 2
1+eps+eps+eps
# 1.0
(1+eps)+(eps+eps)
# 1.0000000000000002
X = np.full((32, 32), eps)
X[0, 0] = 1
X.sum(0)[0]
# 1.0
X.sum(1)[0]
# 1.000000000000003
X[:, 0].sum()
# 1.000000000000003
这强烈表明一维数组和连续轴使用成对求和,而多维数组中的跨轴不使用。
请注意,要看到效果,数组必须足够大,否则 numpy 会退回到普通求和。