численная точность кросс-продукта в numpy

Question

Я испытываю то, что, по моему мнению, странное поведение с использованием кросс-продукта в numpy с несколько большими значениями.

Е.Г., следующее кажется правильным:

r = 1e15 
a = array([1, 2, 3]) * r; 
b = array([-1, 2, 1]) * r; 

c = cross(a/norm(a), b/norm(b)); 

print(dot(c, a)) # outputs 0.0 

Но если поднять показатель на 1, мы получим:

r = 1e16 
a = array([1, 2, 3]) * r; 
b = array([-1, 2, 1]) * r; 

c = cross(a/norm(a), b/norm(b)); 

print(dot(c, a)) # outputs 2.0 

Числа получить еще более странным для больших значений показателя. Кто-нибудь знает, что здесь происходит? Благодаря!

Answer 1

Вы видите ошибку округления. По умолчанию array() возвращает объект с dtype=float64. По мере того как вы делаете r все больше и больше, у вас заканчивается пространство мантиссы, чтобы точно представлять массивные продукты. Вот способ проверить это:

def testcross(r, dt): 
    a = array([1, 2, 3], dtype=dt)*r 
    b = array([-1, 2, 1], dtype=dt)*r 
    c = cross(a/norm(a), b/norm(b)) 
    return dot(c, a) 

for rr in logspace(4, 15, 10): 
    print "%10.2f %10.2f %g" % (testcross(rr, float32), testcross(rr, float64) 

С результатом:

 -0.00  0.00 10000 
     0.00  -0.00 166810 
     0.00  0.00 2.78256e+06 
    -4.00  0.00 4.64159e+07 
    -64.00  0.00 7.74264e+08 
    1024.00  0.00 1.29155e+10 
     0.00  0.00 2.15443e+11 
-524288.00  0.00 3.59381e+12 
     0.00  -0.02 5.99484e+13 
-134217728.00  0.00 1e+15 

Примечание вещи не "идеальный" даже для float64 с r=5.99484e13. Это показывает, что точность начинает разваливаться задолго до того, как вы доберетесь до r=1e15, даже для float64. Как и ожидалось, дела обстоят намного хуже с менее точным float32.

Следующее предложение OP: поля мантиссы для 32 и 64-битного представления с плавающей запятой составляют 24 и 53 бит соответственно (включая подразумеваемый бит). Принимая log10([2**24, 2**53]), мы видим, что это соответствует примерно 7 и 16 порядков соответственно. Это соответствует табулированным ошибкам, появляющимся около r=4.6e7 для float32 и r=1e16, как первоначально отмечалось. Округление происходит, когда точечный продукт вызывает вычисление базовой матрицы для вычитания больших чисел, и различия не могут быть представлены одной или другой мантиссой большого числа.