Блог

Мне любопытно, может ли кто-нибудь объяснить, что именно приводит к несоответствию в этой конкретной обработке упорядоченных массивов 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
  

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

Математика с плавающей запятой не обязательно ассоциативна, т.е. (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 тогда как порядок Fortran использует 1e100 (-1e100) 1 0 . Проблема в том, что (1e100 1) == 1e100 поскольку значения с плавающей запятой недостаточно точны для представления этой небольшой разницы, 1 значение теряется.

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