Numpy sum (Numpy 1.15.4, связана с MKL) - PullRequest
Новая
       13

Numpy sum (Numpy 1.15.4, связана с MKL)

0 голосов
/ 05 января 2019

Я использую последнюю версию Python / Numpy (1.15.4) от Anaconda, связанную с MKL. Я использую следующий эпсилон: 

epsilon = 2**(-53)

Так что 1,0 + эпсилон равен 1,0. Затем я определяю следующий массив numpy, который заполнен эпсилоном, кроме первых 8 элементов, равных единице. 

import numpy as np
n = 1000000
a = np.full(n, epsilon)
a[0:8] = 1.0

Если вы вычислите сумму массива с классическим сокращением слева направо, вы должны получить ровно 8,0, так как все эпсилоны не изменяют сокращение по мере его продолжения. Но я удивлен, увидев, что 

print(np.sum(a))

возвращает вам 8.000000000111008, что не соответствует ожиданиям. Я попытался поместить их в середину и конец массива, чтобы проверить, не сложилась ли сумма задом наперед, но ничего не изменилось. Любая идея о том, как сумма сделана?

PS: Мне хорошо известен тот факт, что арифметика с плавающей запятой сложна и + неассоциативна с плавающей запятой, что делает результаты редукции зависимыми от порядка суммирования. Но я не могу выяснить, какой порядок суммирования здесь используется.

1 Ответ

0 голосов
/ 05 января 2019

Как объяснено в комментариях, numpy использует парное суммирование. Таким образом, суммирование (примерно) будет выглядеть так в конце рекурсивного попарного суммирования, когда стек вызовов начнет разрешаться: 

(1+1) + (1+1) + (1+1) + (1+1) + (epsilon + epsilon) + ... + (epsilon + epsilon)
(2+2) + (2+2) + (2*epsilon) + (2*epsilon) + ... + (2*epsilon)
(4+4) + (4*epsilon) + (4*epsilon) + ... + (4*epsilon)
8 + (8*epsilon) + (8*epsilon) + ... + (8*epsilon)
8 + (16*epsilon) + (16*epsilon) + ... + (16*epsilon)
...
8 + (999992*epsilon)

Вы правы, утверждая, что 1.0 + epsilon равно 1.0. Заманчиво думать, что x + epsilon == x для всех x. Это действительно, когда x является «большим», однако, это не удерживается, когда x == epsilon (то есть epsilon + epsilon != epsilon). Таким образом epsilon + epsilon условия начнут накапливаться: 

In [27]: epsilon = 2**(-53)

In [28]: 1.0 + epsilon == 1.0
Out[28]: True

In [29]: 2.0 + epsilon == 2.0
Out[29]: True

In [30]: epsilon + epsilon == epsilon
Out[30]: False

In [31]: epsilon
Out[31]: 1.1102230246251565e-16

In [32]: epsilon + epsilon
Out[32]: 2.220446049250313e-16

In [33]: 123*epsilon
Out[33]: 1.3655743202889425e-14

Я не могу получить numpy s ответ, но мы можем подойти очень близко: 

In [36]: 8 + (999992*epsilon)
Out[36]: 8.000000000111022

In [62]: def pairwise_sum(arr):
    ...:   if len(arr) <= 2:
    ...:     return sum(arr)
    ...:   midpoint = len(arr)//2
    ...:   first_half = arr[:midpoint]
    ...:   second_half = arr[midpoint:]
    ...:   return pairwise_sum(first_half) + pairwise_sum(second_half)
    ...:

In [63]: pairwise_sum(a)
Out[63]: 8.0000000001110205

Так что здесь явно есть некоторое другое различие (некоторые детали реализации), но, надеюсь, это убеждает вас, что numpy на самом деле использует парное суммирование.

НТН.

...