Numpy: Разность между точкой (а, б) и (а * б) .sum()

Question

В 1-D Numpy массивов, эти два выражения должны давать один и тот же результат (theorically):

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

последний использует dot() и быстрее. Но какой из них более точен? Зачем?

Ниже приведен контекст.

Я хотел вычислить взвешенную дисперсию образца с использованием numpy. Я нашел выражение dot() в another answer с комментарием о том, что он должен быть более точным. Однако никаких объяснений здесь нет.

Answer 1

Numpy dot - одна из подпрограмм, которая вызывает библиотеку BLAS, которую вы связываете при компиляции (или создает свои собственные). Важность этого заключается в том, что библиотека BLAS может использовать операции Multiply-accumulate (обычно Fused-Multiply Add), которые ограничивают количество округлений, которые выполняет вычисление.

Возьмите следующее:

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

Не точно, но достаточно близко.

>>> 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) будет более точным из двух, как это использовать приблизительно половину числа с плавающей точкой закруглений, что наивные (a*a).sum() делает.

Книга Nvidia имеет следующий пример для 4 цифр точности. rn обозначает 4 тура до ближайших 4-х цифр:

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 

чисел с плавающей точкой, конечно, не округляется до 16-го знака после запятой в базе 10, но вы получите идею.

Размещение np.dot(a,a) в обозначениях с некоторым дополнительным псевдокодом:

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

Хотя (a*a).sum() является:

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) 

Отсюда его легко видеть, что число округляется в два раза больше времени, используя (a*a).sum() по сравнению с np.dot(a,a). Эти небольшие различия, которые могут быть суммированы, могут изменить ответ. Дополнительные примеры можно найти: here.