numerische Präzision von Kreuzprodukt in numpy

Question

Ich erlebe, was ich denke, ist seltsames Verhalten mit Cross-Produkt in numpy mit etwas großen Werten.

Z. B. scheint die folgende richtige:

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 

Aber wenn wir die Exponenten um 1 erhöhen, erhalten wir:

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 

Die Zahlen für größere Werte des Exponenten noch bizarrer bekommen. Weiß jemand, was hier vor sich geht? Vielen Dank!

Answer 1

Sie sehen Rundungsfehler. Standardmäßig gibt array() ein Objekt mit dtype=float64 zurück. Wenn Sie r größer und größer machen, gehen Ihnen die Mantissenräume aus, um die Array-Produkte genau darzustellen. Hier ist ein Weg, dies zu testen:

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) 

Mit dem Ergebnis:

 -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 

Hinweis sind die Dinge nicht "perfekt" auch für float64 mit r=5.99484e13. Dies zeigt, dass die Genauigkeit beginnt auseinander zu fallen, lange bevor Sie zu r=1e15 kommen, auch für float64. Wie erwartet, sind die Dinge mit dem weniger präzisen float32 viel schlechter.

Dem Vorschlag von OP folgend: die Mantissenfelder für 32- und 64-Bit Gleitkommadarstellung sind 24 bzw. 53 Bits (einschließlich des implizierten Bits). Wenn wir log10([2**24, 2**53]) betrachten, sehen wir, dass dies ungefähr 7 bzw. 16 Grßenordnungen entspricht. Dies entspricht den tabellarisierten Fehlern, die sich um r=4.6e7 für float32 und r=1e16 wie ursprünglich angegeben ergeben. Die Abrundung tritt auf, wenn das Skalarprodukt bewirkt, dass die zugrundeliegende Matrixberechnung große Zahlen subtrahiert, und die Differenzen können nicht als eine oder die andere Mantisse einer großen Zahl dargestellt werden.