2013-08-07 31 views
16

Do 1-D tablic NumPy, te dwa wyrażenia należy takiego samego wyniku (theorically)Numpy Różnica pomiędzy punktu (a, b) i (a * b) .sum()

(a*b).sum()/a.sum() 
dot(a, b)/a.sum() 

ostatnio używa dot() i jest szybszy. Ale który z nich jest dokładniejszy? Czemu?

Następuje kontekst.

Chciałem obliczyć ważoną wariancję próbki za pomocą numpy. Znalazłem wyrażenie dot() w another answer, z komentarzem stwierdzającym, że powinno być dokładniejsze. Jednak nie podano tam żadnego wyjaśnienia.

+1

Jest to średnia ważona, prawda? Możesz po prostu użyć ['np.average'] (http://docs.scipy.org/doc/numpy/reference/generated/numpy.average.html). – user2357112

+0

Myślę, że część dotycząca "dokładności liczbowej" odnosiła się do odejmowania średniej od wartości, zamiast używania "kropki". – user2357112

Odpowiedz

9

Numpy dot to jedna z procedur, które wywołują bibliotekę BLAS, którą łączymy podczas kompilacji (lub tworzymy ją samodzielnie). Ważność tego jest taka, że ​​biblioteka BLAS może korzystać z operacji kumulowania wielokrotności (zazwyczaj Fused-Multiply Add), które ograniczają liczbę zaokrągleń wykonywanych przez obliczenia.

podjąć następujące:

>>> a=np.ones(1000,dtype=np.float128)+1E-14 
>>> (a*a).sum() 
1000.0000000000199948 
>>> np.dot(a,a) 
1000.0000000000199948 

Nie dokładnie, ale wystarczająco blisko.

>>> a=np.ones(1000,dtype=np.float64)+1E-14 
>>> np.dot(a,a) 
1000.0000000000176 #off by 2.3948e-12 
>>> (a*a).sum() 
1000.0000000000059 #off by 1.40948e-11 

The np.dot(a, a) będzie dokładniejsza z nich korzystać, ponieważ około połowa liczby pływających zaokrągleń punktowych że naiwny (a*a).sum() robi.

Książka Nvidii ma następujący przykład dla 4 cyfr precyzji. rn oznacza 4 rundzie z dokładnością do 4 cyfr:

x = 1.0008 
x2 = 1.00160064     # true value 
rn(x2 − 1) = 1.6006 × 10−4   # fused multiply-add 
rn(rn(x2) − 1) = 1.6000 × 10−4  # multiply, then add 

oczywiście liczb zmiennoprzecinkowych nie są zaokrąglane do 16. miejsca po przecinku w bazie 10, ale masz pomysł.

Umieszczenie np.dot(a,a) w powyższym zapisie z jakimś dodatkowym kodem pseudo:

out=0 
for x in a: 
    out=rn(x*x+out) #Fused multiply add 

Podczas (a*a).sum() jest:

arr=np.zeros(a.shape[0]) 
for x in range(len(arr)): 
    arr[x]=rn(a[x]*a[x]) 

out=0 
for x in arr: 
    out=rn(x+out) 

Z tego łatwo jest zauważyć, że liczba jest zaokrąglana dwa razy tyle razy, używając (a*a).sum() w porównaniu z np.dot(a,a). Te niewielkie różnice sumowane mogą nieznacznie zmienić odpowiedź. Dodatkowe przykłady można znaleźć here.

+5

Jeśli numpy używa zoptymalizowanych błon dla maszyny użytkownika i ma procesor z fma. Nie należy przyjmować zbyt wielu założeń w oparciu o "the blas * can * doing that ..." –

+0

Tak, nawet dla Intel Ivy Bridge 'a + b * c' kompiluje się do' mulss', po którym następuje 'addss'. –

+0

Czy dodano odpowiednie opcje kompilatora? (jak -march = core-avx2 -mavx -mfma ... z gcc) –