Warum ist np.dot ungenau? (n-dim Arrays)

15

Angenommen, wir nehmen np.dotzwei 'float32'2D-Arrays:

res = np.dot(a, b)   # see CASE 1
print(list(res[0]))  # list shows more digits
[-0.90448684, -1.1708503, 0.907136, 3.5594249, 1.1374011, -1.3826287]

Zahlen. Außer, sie können sich ändern:


Fall 1 : Scheibea

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(6, 6).astype('float32')

for i in range(1, len(a)):
    print(list(np.dot(a[:i], b)[0])) # full shape: (i, 6)
[-0.9044868,  -1.1708502, 0.90713596, 3.5594249, 1.1374012, -1.3826287]
[-0.90448684, -1.1708503, 0.9071359,  3.5594249, 1.1374011, -1.3826288]
[-0.90448684, -1.1708503, 0.9071359,  3.5594249, 1.1374011, -1.3826288]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]

Die Ergebnisse unterscheiden sich, obwohl die gedruckte Schicht aus genau den gleichen multiplizierten Zahlen stammt.


FALL 2 : Abflachen a, eine 1D-Version von nehmen bund dann in Scheiben schneiden a:

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(1, 6).astype('float32')

for i in range(1, len(a)):
    a_flat = np.expand_dims(a[:i].flatten(), -1) # keep 2D
    print(list(np.dot(a_flat, b)[0])) # full shape: (i*6, 6)
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]

Fall 3 : stärkere Kontrolle; Setzen Sie alle nicht beteiligten Entires auf Null : Fügen Sie a[1:] = 0den CASE 1-Code hinzu. Ergebnis: Diskrepanzen bleiben bestehen.


FALL 4 : Überprüfen Sie andere Indizes als [0]; Wie bei [0]beginnen die Ergebnisse, eine feste Anzahl von Array-Vergrößerungen von ihrem Erstellungspunkt aus zu stabilisieren. Ausgabe

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(6, 6).astype('float32')

for j in range(len(a) - 2):
    for i in range(1, len(a)):
        res = np.dot(a[:i], b)
        try:    print(list(res[j]))
        except: pass
    print()

Daher unterscheiden sich die Ergebnisse für den 2D * 2D-Fall - sind jedoch für 1D * 1D konsistent. Aus einigen meiner Messwerte geht hervor, dass dies auf 1D-1D mit einfacher Addition zurückzuführen ist, während 2D-2D eine „schickere“, leistungssteigernde Addition verwendet, die möglicherweise weniger genau ist (z. B. die paarweise Addition bewirkt das Gegenteil). Trotzdem kann ich nicht verstehen, warum Diskrepanzen verschwinden, wenn 1 einmal aüber einen festgelegten „Schwellenwert“ hinausgeschnitten wird. Je größer aund bdesto später scheint diese Schwelle zu liegen, aber sie existiert immer.

Alle sagten: Warum ist np.dotND-ND-Arrays ungenau (und inkonsistent)? Relevante Git


Zusätzliche Informationen :

  • Umgebung : Win-10 OS, Python 3.7.4, Spyder 3.3.6 IDE, Anaconda 3.0 2019/10
  • CPU : i7-7700HQ 2,8 GHz
  • Numpy v1.16.5

Mögliche Täterbibliothek: Numpy MKL - auch BLASS-Bibliotheken; Vielen Dank an Bi Rico für die Notiz


Stresstest-Code : Wie bereits erwähnt, verschärfen sich die Abweichungen in der Frequenz mit größeren Arrays. Wenn oben nicht reproduzierbar ist, sollte unten sein (wenn nicht, versuchen Sie es mit größeren Abmessungen). Meine Ausgabe

np.random.seed(1)
a = (0.01*np.random.randn(9, 9999)).astype('float32') # first multiply then type-cast
b = (0.01*np.random.randn(9999, 6)).astype('float32') # *0.01 to bound mults to < 1

for i in range(1, len(a)):
    print(list(np.dot(a[:i], b)[0]))

Schweregrad des Problems : Die angezeigten Abweichungen sind "gering", jedoch nicht mehr, wenn Sie in einem neuronalen Netzwerk arbeiten, in dem Milliarden von Zahlen über einige Sekunden und Billionen über die gesamte Laufzeit multipliziert werden. Die gemeldete Modellgenauigkeit unterscheidet sich pro Thread um ganze 10 Prozent .

Nachfolgend finden Sie eine gif von Arrays von der Fütterung auf ein Modell resultierenden was im Grunde ist a[0]w /, len(a)==1gegen len(a)==32:


ANDERE PLATTFORMEN Ergebnisse, und dank Pauls Tests:

Fall 1 (teilweise) reproduziert :

  • Google Colab VM - Intel Xeon 2,3 G-Hz - Jupyter - Python 3.6.8
  • Win-10 Pro Docker Desktop - Intel i7-8700K - Jupyter / Scipy-Notebook - Python 3.7.3
  • Ubuntu 18.04.2 LTS + Docker - AMD FX-8150 - Jupyter / Scipy-Notebook - Python 3.7.3

Hinweis : Diese ergeben einen viel geringeren Fehler als oben gezeigt. Zwei Einträge in der ersten Zeile sind in der niedrigstwertigen Ziffer von den entsprechenden Einträgen in den anderen Zeilen um 1 versetzt.

Fall 1 nicht reproduziert :

  • Ubuntu 18.04.3 LTS - Intel i7-8700K - IPython 5.5.0 - Python 2.7.15+ und 3.6.8 (2 Tests)
  • Ubuntu 18.04.3 LTS - Intel i5-3320M - IPython 5.5.0 - Python 2.7.15+
  • Ubuntu 18.04.2 LTS - AMD FX-8150 - IPython 5.5.0 - Python 2.7.15rc1

Anmerkungen :

  • Die verknüpften Colab-Notebook- und Jupyter-Umgebungen weisen eine weitaus geringere Diskrepanz auf (und nur für die ersten beiden Zeilen) als auf meinem System. Auch Fall 2 zeigte (noch) nie Ungenauigkeit.
  • In diesem sehr begrenzten Beispiel ist die aktuelle (Dockerized) Jupyter-Umgebung anfälliger als die IPython-Umgebung.
  • np.show_config()zu lang zum Posten, aber zusammenfassend: IPython-Envs basieren auf BLAS / LAPACK; Colab basiert auf OpenBLAS. In IPython Linux envs sind BLAS-Bibliotheken systeminstalliert - in Jupyter und Colab stammen sie aus / opt / conda / lib

UPDATE : Die akzeptierte Antwort ist korrekt, aber breit und unvollständig. Die Frage bleibt offen für jeden, der das Verhalten bei der erklären kann Code - Ebene - nämlich ein genauer Algorithmus durch np.dot, und wie es im Einklang Inkonsistenzen 'beobachtete in obigen Ergebnissen (siehe auch Anmerkungen) erläutert. Hier sind einige direkte Implementierungen, die über meine Entschlüsselung hinausgehen: sdot.c - arraytypes.c.src

OverLordGoldDragon
quelle
Kommentare sind nicht für eine ausführliche Diskussion gedacht. Dieses Gespräch wurde in den Chat verschoben .
Samuel Liew
Die allgemeinen Algorithmen ndarraysignorieren normalerweise den numerischen Präzisionsverlust. Da der Einfachheit halber sie reduce-sumentlang jeder Achse kann die Reihenfolge der Operationen möglicherweise nicht der optimale ... Beachten Sie, dass , wenn Sie Genauigkeit Fehler dagegen könnten Sie auch nutzenfloat64
Vitor SRG
Möglicherweise habe ich morgen keine Zeit für eine Überprüfung, daher vergebe ich jetzt das Kopfgeld.
Paul
@ Paul Es würde sowieso automatisch an die Antwort mit der höchsten Stimme vergeben - aber in Ordnung, danke für die Benachrichtigung
OverLordGoldDragon

Antworten:

7

Dies sieht nach unvermeidbarer numerischer Ungenauigkeit aus. Wie hier erläutert , verwendet NumPy eine hochoptimierte, sorgfältig abgestimmte BLAS-Methode zur Matrixmultiplikation . Dies bedeutet, dass sich wahrscheinlich die Reihenfolge der Operationen (Summe und Produkte), die zum Multiplizieren von 2 Matrizen befolgt werden, ändert, wenn sich die Größe der Matrix ändert.

Um klarer zu sein, wissen wir, dass mathematisch jedes Element der resultierenden Matrix als Punktprodukt zweier Vektoren (gleichlange Zahlenfolgen) berechnet werden kann . Aber das ist nicht , wie NumPy ein Element der resultierenden Matrix berechnet. Tatsächlich gibt es effizientere, aber komplexere Algorithmen wie den Strassen-Algorithmus , die das gleiche Ergebnis erzielen , ohne das Zeilen-Spalten-Punktprodukt direkt zu berechnen.

Selbst wenn das Element C ij einer resultierenden Matrix C = AB mathematisch als das Punktprodukt der i-ten Zeile von A mit der j-ten Spalte von B definiert ist , multiplizieren Sie bei Verwendung solcher Algorithmen eine Matrix A2 mit der Dieselbe i-te Zeile wie A mit einer Matrix B2 mit derselben j-ten Spalte wie B , das Element C2 ij wird tatsächlich nach einer anderen Folge von Operationen berechnet (die von der gesamten A2 und B2 abhängt Matrizen), was möglicherweise zu unterschiedlichen numerischen Fehlern führt.

Selbst wenn mathematisch C ij = C2 ij ist (wie in Ihrem FALL 1), führt die unterschiedliche Abfolge von Operationen, denen der Algorithmus in den Berechnungen folgt (aufgrund der Änderung der Matrixgröße), zu unterschiedlichen numerischen Fehlern. Der numerische Fehler erklärt auch die geringfügig unterschiedlichen Ergebnisse in Abhängigkeit von der Umgebung und die Tatsache, dass in einigen Fällen in einigen Umgebungen der numerische Fehler möglicherweise nicht vorhanden ist.

mmj
quelle
2
Vielen Dank für den Link, er scheint relevante Informationen zu enthalten - Ihre Antwort könnte jedoch detaillierter sein, da es sich bisher um eine Umschreibung von Kommentaren unterhalb der Frage handelt. Das verknüpfte SO zeigt beispielsweise direkten CCode an und bietet Erklärungen auf Algorithmenebene, sodass es in die richtige Richtung weist .
OverLordGoldDragon
Es ist auch nicht "unvermeidlich", wie am Ende der Frage gezeigt - und das Ausmaß der Ungenauigkeit variiert je nach Umgebung, was ungeklärt bleibt
OverLordGoldDragon
1
@OverLordGoldDragon: (1) Ein triviales Beispiel mit Zusatz: Nummer nnehmen, Nummer kso nehmen, dass sie unter der Genauigkeit der kletzten Mantissenziffer liegt. Für Pythons native schwebt n = 1.0und k = 1e-16funktioniert. Nun lass ks = [k] * 100. Sehen Sie, dass sum([n] + ks) == n, während sum(ks + [n]) > ndas heißt, die Summierungsreihenfolge wichtig war. (2) Moderne CPUs haben mehrere Einheiten, um Gleitkommaoperationen (FP) parallel auszuführen, und die Reihenfolge, in der a + b + c + dauf einer CPU berechnet wird, ist nicht definiert, selbst wenn der Befehl im Maschinencode a + bvorher kommt c + d.
9000
1
@OverLordGoldDragon Sie sollten sich bewusst sein, dass die meisten Zahlen, mit denen Ihr Programm umgehen soll, nicht genau durch einen Gleitkomma dargestellt werden können. Versuchen Sie es format(0.01, '.30f'). Wenn selbst eine einfache Zahl wie 0.01nicht genau durch ein NumPy-Gleitkomma dargestellt werden kann, müssen die tiefen Details des NumPy-Matrixmultiplikationsalgorithmus nicht bekannt sein, um den Punkt meiner Antwort zu verstehen. Das heißt, unterschiedliche Startmatrizen führen zu unterschiedlichen Abfolgen von Operationen , so dass sich mathematisch gleiche Ergebnisse aufgrund numerischer Fehler geringfügig unterscheiden können.
mmj
2
@OverLordGoldDragon re: schwarze Magie. Es gibt ein Papier, das bei einigen CS MOOCs gelesen werden muss. Mein Rückruf ist nicht so toll, aber ich denke, das ist es: itu.dk/~sestoft/bachelor/IEEE754_article.pdf
Paul