что вызывает различия в сумме массива вдоль оси для C и F упорядоченных массивов в NumPy

Question

Мне любопытно, если кто-нибудь может объяснить, что именно приводит к расхождению в этой конкретной обработке массивов C и Fortran в numpy. Смотрите код ниже:

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

Answer 1

Это почти наверняка является следствием того, что numpy иногда использует парное суммирование и иногда не .

Давайте построим диагностический массив:

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 возвращается к обычному суммированию.

Answer 2

Математика с плавающей точкой не обязательно ассоциативно , то есть (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() по-прежнему дает 0, как и aT = a.T; aT.sum(), поэтому я не уверен, как именно это реализовано вбэкэнд)

В порядке C используется последовательность операций (слева направо) 1e100 + 1 + (-1e100) + 0, тогда как в порядке Фортрана используется 1e100 + (-1e100) + 1 + 0.Проблема в том, что (1e100+1) == 1e100, потому что поплавки не обладают достаточной точностью, чтобы представить эту небольшую разницу, поэтому 1 теряется.

Как правило, не выполняйте тестирование на равенство для чисел с плавающей запятой, вместо этого сравните, используя маленький эпсилон (if abs(float1 - float2) < 0.00001 или np.isclose).Если вам нужна произвольная точность с плавающей точкой, используйте библиотеку Decimal или представление с фиксированной запятой и int s.