#python #numpy #row-major-order #column-major-order
#python #numpy #строка-основной порядок #столбец-основной порядок
Вопрос:
Мне любопытно, может ли кто-нибудь объяснить, что именно приводит к несоответствию в этой конкретной обработке упорядоченных массивов 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
Комментарии:
1. Это странно. Если я настрою ваш код на суммирование по другой оси, несоответствие исчезнет.
Ответ №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 возвращается к обычному суммированию.
Комментарии:
1. после обсуждения с одним из разработчиков numpy это правильное объяснение. Я был в плохой ситуации, когда я использовал 32-разрядные числа с плавающей запятой и суммировал миллионы из них с наивными суммами, что привело к очень большим ошибкам округления в моем случае. Размещение ссылки здесь, поскольку другие могут счесть это полезным, если они столкнутся с одним и тем же: github.com/numpy/numpy/issues/13240
Ответ №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
тогда как порядок Fortran использует 1e100 (-1e100) 1 0
. Проблема в том, что (1e100 1) == 1e100
поскольку значения с плавающей запятой недостаточно точны для представления этой небольшой разницы, 1
значение теряется.
В общем, не проводите проверку на равенство для чисел с плавающей запятой, вместо этого сравнивайте, используя маленький эпсилон ( if abs(float1 - float2) < 0.00001
или np.isclose
). Если вам нужна произвольная точность с плавающей запятой, используйте Decimal
библиотеку или представление с фиксированной точкой и int
s.
Комментарии:
1.
a.T
этоview
ofa
, один и тот же буфер данных, просто разные формы и шаги; эффективно меняется отC
кF
порядку. Дляsum()
я полагаю, что это просто суммирует databuffer, без учета порядка поиска.