Pourquoi np.dot est-il imprécis? (tableaux n-dim)

Supposons que nous prenions np.dotdeux 'float32'tableaux 2D:

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]

Nombres. Sauf, ils peuvent changer:

CAS 1 : tranchea

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]

Les résultats diffèrent, même si la tranche imprimée provient des mêmes nombres multipliés.

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]

CAS 3 : contrôle renforcé; mettre tous les entiers non impliqués à zéro : ajouter a[1:] = 0au code CASE 1. Résultat: des écarts persistent.

CAS 4 : vérifier les indices autres que [0]; comme pour [0], les résultats commencent à stabiliser un nombre fixe d'agrandissements de réseau à partir de leur point de création. Production

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()

Par conséquent, pour le cas 2D * 2D, les résultats diffèrent - mais sont cohérents pour 1D * 1D. D'après certaines de mes lectures, cela semble provenir de 1D-1D utilisant un ajout simple, tandis que 2D-2D utilise un ajout plus sophistiqué et plus performant qui peut être moins précis (par exemple, l'addition par paire fait le contraire). Néanmoins, je n'arrive pas à comprendre pourquoi les écarts disparaissent dans le cas où 1 aest une fois dépassé un «seuil» défini; plus grand aet plus btard ce seuil semble se situer, mais il existe toujours.

Tout cela dit: pourquoi est np.dotimprécis (et incohérent) pour les matrices ND-ND? Git pertinent

Informations supplémentaires :

• Environnement : 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

Bibliothèque coupable possible : Numpy MKL - également bibliothèques BLASS; merci à Bi Rico d' avoir noté

Code de test de résistance : comme indiqué, les écarts aggravent la fréquence avec des réseaux plus grands; si ci-dessus n'est pas reproductible, ci-dessous devrait être (sinon, essayez des dims plus grands). Ma sortie

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]))

Gravité du problème : les écarts indiqués sont «faibles», mais plus lorsqu'ils fonctionnent sur un réseau de neurones avec des milliards de chiffres multipliés en quelques secondes et des milliards sur toute la durée d'exécution; la précision du modèle rapporté diffère de 10 pour cent entiers, pour ce fil .

Ci-dessous est un gif de tableaux résultant de l'alimentation d'un modèle ce qui est fondamentalement a[0], w / len(a)==1vs len(a)==32:

AUTRES PLATEFORMES résultats, selon et grâce aux tests de Paul :

Le cas 1 reproduit (en partie) :

• 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

Remarque : ceux-ci produisent une erreur beaucoup plus faible que celle indiquée ci-dessus; deux entrées sur la première ligne sont décalées de 1 dans le chiffre le moins significatif des entrées correspondantes dans les autres lignes.

Cas 1 non reproduit :

• Ubuntu 18.04.3 LTS - Intel i7-8700K - IPython 5.5.0 - Python 2.7.15+ et 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

Remarques :

• Les environnements de bloc-notes et de jupyter Colab liés présentent une différence bien moindre (et uniquement pour les deux premières lignes) que celle observée sur mon système. De plus, le cas 2 n'a jamais (encore) fait preuve d'imprécision.
• Dans cet échantillon très limité, l'environnement Jupyter actuel (Dockerisé) est plus sensible que l'environnement IPython.
• np.show_config()trop long pour poster, mais en résumé: les envs IPython sont basés sur BLAS / LAPACK; Colab est basé sur OpenBLAS. Dans les environnements IPython Linux, les bibliothèques BLAS sont installées par le système - dans Jupyter et Colab, elles proviennent de / opt / conda / lib

MISE À JOUR : la réponse acceptée est exacte, mais large et incomplète. La question reste ouverte à quiconque peut expliquer le comportement au niveau du code - à savoir, un algorithme exact utilisé par np.dot, et comment il explique les `` incohérences cohérentes '' observées dans les résultats ci-dessus (voir également les commentaires). Voici quelques implémentations directes au-delà de mon déchiffrement: sdot.c - arraytypes.c.src

Cela ressemble à une imprécision numérique inévitable. Comme expliqué ici , NumPy utilise une méthode BLAS hautement optimisée et soigneusement réglée pour la multiplication matricielle . Cela signifie que probablement la séquence d'opérations (somme et produits) suivie pour multiplier 2 matrices, change lorsque la taille de la matrice change.

En essayant d'être plus clair, nous savons que, mathématiquement , chaque élément de la matrice résultante peut être calculé comme le produit scalaire de deux vecteurs (séquences de nombres de longueur égale). Mais ce n'est pas ainsi que NumPy calcule un élément de la matrice résultante. En fait, il existe des algorithmes plus efficaces mais complexes, comme l' algorithme Strassen , qui obtiennent le même résultat sans calculer directement le produit scalaire ligne-colonne.

Lorsque vous utilisez de tels algorithmes, même si l'élément C ij d'une matrice résultante C = AB est défini mathématiquement comme le produit scalaire de la i-ème ligne de A avec la j-ème colonne de B , si vous multipliez une matrice A2 ayant la même i-ème ligne que A avec une matrice B2 ayant la même j-ème colonne que B , l'élément C2 ij sera effectivement calculé à la suite d'une séquence d'opérations différente (cela dépend de l'ensemble A2 et B2 matrices), pouvant conduire à différentes erreurs numériques.

C'est pourquoi, même si mathématiquement C ij = C2 ij (comme dans votre CAS 1), la séquence d'opérations différente suivie par l'algorithme dans les calculs (en raison du changement de taille de la matrice) conduit à des erreurs numériques différentes. L'erreur numérique explique également les résultats légèrement différents selon l'environnement et le fait que, dans certains cas, pour certains environnements, l'erreur numérique peut être absente.