Что вызывает различия в сумме массива вдоль оси для 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 библиотека или представление с фиксированной запятой и ints.