Comparaison de la vitesse de calcul en implémentant python mpmath (calcul de précision arbitraire) (Note)
Contexte
Cela fait environ deux ans depuis le dernier Post. Nous avons implémenté différentes langues dans l'article précédent et comparé les vitesses de calcul. La motivation est la ligne cachée qui veut comparer la vitesse de calcul par un calcul de précision arbitraire. Pour un calcul de précision arbitraire, l'implémentation de GNU MP (GMP), MPFR (à la fois C et C ++) est majeure. ， Les wrappers et extensions de ces packages sont disponibles dans certaines langues. En python, il est disponible dans mpmath (je ne connais pas Decimal de la bibliothèque standard), et julia et Mathematica (langage wolfram) supportent par défaut des opérations de longueur multiple. Les deux semblent utiliser GMP et MPFR, donc le but ultime est de savoir combien de temps système et de codage la facilité de mise en œuvre est (pas atteint cette fois).
En passant, le calcul de précision 4x est pris en charge à partir de gcc 4.6, et le calcul de précision 4x devrait être possible en C, C ++ et Fortran, mais cela n'est pas traité ici. (En C, le type à virgule flottante de précision 4x nécessitait une déclaration comme __float128, mais qu'en est-il maintenant ...)
Exemples nécessitant une précision arbitraire
(Unidimensionnel) [carte du boulanger](https://ja.wikipedia.org/wiki/%E3%83%91%E3%82%A4%E3%81%93%E3%81%AD % E5% A4% 89% E6% 8F% 9B)
Pour le moment, écrivez (grossièrement) le code qui trace les itérations jusqu'au 3ème étage avec les conditions initiales appropriées.
import numpy as np
import matplotlib.pyplot as plt
import time
def B(x):
return 2*x if x <= 1/2 else (2*x - 1)
fig, ax = plt.subplots(1, 1, figsize=(5,5))
ax.set_xlabel("x")
ax.set_ylabel("B(x)")
x = np.linspace(0, 1, 300)
ax.set_ylim(-0.02,1.02)
ax.set_xlim(-0.02,1.02)
ax.plot(x,2*x,"-y")
ax.plot(x,2*x-1,"-y")
ax.plot(x,x,"--k")
x = (np.sqrt(5) - 1)/2
web = [np.array([x]),
np.array([0])]
traj = np.array([x])
nmax = 3
bboxprop = dict(facecolor='w', alpha=0.7)
ax.text(x+0.05,0, r"$x=%f$" % x, bbox=bboxprop)
a = time.time()
for n in range(1,nmax+1):
y = B(x)
web[0] = np.append(web[0], x)
web[1] = np.append(web[1], y)
traj = np.append(traj, y)
ax.text(x+0.02, y, r"$f^{(%d)}(x)$" % n, bbox=bboxprop)
x = y
if n != nmax:
ax.text(y+0.02, y, r"$x=f^{(%d)}(x)$" % n, bbox=bboxprop)
web[0] = np.append(web[0], x)
web[1] = np.append(web[1], y)
print(time.time() - a, "[sec.]")
ax.plot(web[0][:], web[1][:], "-o")
ax.plot(web[0][0], web[1][0], "or")
ax.plot(web[0][-1], web[1][-1], "or")
ax.grid()
print(traj)
plt.savefig("web.png ")
plt.show()
Si vous exécutez, vous obtiendrez la séquence orbitale de $ \ {x, B (x), B \ circ B (x), B \ circ B \ circ B (x) \} $.
$ python baker.py
0.0013127326965332031 [sec.]
[0.61803399 0.23606798 0.47213595 0.94427191]
Aussi, je pense que la figure (diagramme web ou portrait de phase) peut visuellement comprendre le sens de la synthèse itérative (le point de départ et le point final sont entourés de rouge).
Maintenant, que se passe-t-il si nous fixons le nombre d'itérations à 50?
nmax = 3
À
nmax = 50
Et courir
$ python baker-2.py
0.0012233257293701172 [sec.]
[0.61803399 0.23606798 0.47213595 0.94427191 0.88854382 0.77708764
0.55417528 0.10835056 0.21670112 0.43340224 0.86680448 0.73360896
0.46721792 0.93443584 0.86887168 0.73774336 0.47548671 0.95097343
0.90194685 0.8038937 0.60778741 0.21557482 0.43114964 0.86229928
0.72459856 0.44919711 0.89839423 0.79678845 0.59357691 0.18715382
0.37430763 0.74861526 0.49723053 0.99446106 0.98892212 0.97784424
0.95568848 0.91137695 0.82275391 0.64550781 0.29101562 0.58203125
0.1640625 0.328125 0.65625 0.3125 0.625 0.25
0.5 1. 1. ]
Le résultat de a été obtenu. Puisque $ B (1) = 1 $ par définition, une synthèse itérative suffisamment longue $ B ^ {(n)} (x) $ avec la condition initiale $ x = (\ sqrt {5} -1) / 2 $ est mise à 1. On attend du calcul numérique qu'il converge (le panoramique ne se mélange pas uniformément ...)
Cependant, cette prédiction basée sur un calcul numérique a été dérivée à tort en raison d'une précision de calcul insuffisante. Ce système se comporte de manière chaotique pour presque tous les points initiaux sauf pour les conditions initiales spéciales ($ x = 0,1 / 2,1 $, etc.), rendant $ x \ in [0,1] $ dense. remplir.
Pour confirmer cela, le même calcul est effectué en utilisant mpmath, qui est un package de précision arbitraire de python.
import numpy as np
import matplotlib.pyplot as plt
import mpmath
import time
mpmath.mp.dps = 100
def B(x):
return 2*x if x <= 1/2 else (2*x - 1)
fig, ax = plt.subplots(1, 1, figsize=(5,5))
ax.set_xlabel("x")
ax.set_ylabel("B(x)")
x = np.linspace(0, 1, 300)
ax.set_ylim(-0.02,1.02)
ax.set_xlim(-0.02,1.02)
ax.plot(x,2*x,"-y")
ax.plot(x,2*x-1,"-y")
ax.plot(x,x,"--k")
x = (np.sqrt(5) - 1)/2
web = [np.array([x]),
np.array([0])]
traj = np.array([x])
x = (mpmath.sqrt(5) - 1)/2
y = B(x)
web = [np.array([x]),
np.array([0])]
traj = np.array([x])
nmax = 50
a = time.time()
for n in range(1,nmax+1):
y = B(x)
web[0] = np.append(web[0], x)
web[1] = np.append(web[1], y)
traj = np.append(traj, y)
x = y
if n!=nmax:
web[0] = np.append(web[0], x)
web[1] = np.append(web[1], y)
print(time.time() - a, "[sec.]")
ax.plot(web[0], web[1], "-o")
ax.plot(web[0][0], web[1][0], "or")
ax.plot(web[0][-1], web[1][-1], "or")
ax.grid()
mpmath.nprint(traj.tolist())
plt.savefig("web-mp1.png ")
plt.show()
Le changement essentiel est
import mpmath
mpmath.mp.dps = 100
Et conditions initiales
x = (mpmath.sqrt(5) - 1)/2
Seulement. mpmath.mp.dps = 100 spécifie que l'évaluation doit être effectuée avec une précision de 100 chiffres dans la partie formelle. (Par défaut, la double précision (15) est spécifiée.) Si print (traj) est utilisé, 100 chiffres seront sortis tels quels, donc convertissez le traj de mpmath en liste une fois et utilisez nprint pour réduire le nombre de chiffres affichés. Le résultat de sortie est
$ python baker_mp.py
0.0019168853759765625 [sec.]
[0.618034, 0.236068, 0.472136, 0.944272, 0.888544, 0.777088, 0.554175, 0.108351, 0.216701, 0.433402, 0.866804, 0.733609, 0.467218, 0.934436, 0.868872, 0.737743, 0.475487, 0.950973, 0.901947, 0.803894, 0.607787, 0.215575, 0.43115, 0.862299, 0.724599, 0.449197, 0.898394, 0.796788, 0.593577, 0.187154, 0.374308, 0.748615, 0.49723, 0.994461, 0.988921, 0.977842, 0.955685, 0.911369, 0.822739, 0.645478, 0.290956, 0.581912, 0.163824, 0.327647, 0.655294, 0.310589, 0.621177, 0.242355, 0.48471, 0.96942, 0.93884]
Contrairement au résultat du calcul avec une double précision, il ne converge pas vers 1 en 50 pas. La vitesse de calcul est légèrement plus lente, mais elle se situe dans la plage autorisée. Développons le temps plus longtemps Développons le temps avec nmax = 350 La sortie devient également plus longue en fonction du nombre d'étapes, donc si vous affichez les 30 dernières étapes ou plus
mpmath.nprint(traj.tolist()[-30:])
$ python baker_mp.py
0.013526201248168945 [sec.]
[0.0256348, 0.0512695, 0.102539, 0.205078, 0.410156, 0.820313, 0.640625, 0.28125, 0.5625, 0.125, 0.25, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
Il devient. Même dans le calcul des 100 chiffres de la pièce impropre, il semble qu'elle se rapproche progressivement de 1 si le nombre de pas est suffisamment long. Et si nous améliorions encore la précision?
mpmath.mp.dps = 200
Exécutez en tant que.
$ python baker_mp.py
0.013622760772705078 [sec.]
[0.0256048, 0.0512095, 0.102419, 0.204838, 0.409676, 0.819353, 0.638705, 0.27741, 0.55482, 0.109641, 0.219282, 0.438564, 0.877127, 0.754255, 0.50851, 0.017019, 0.0340381, 0.0680761, 0.136152, 0.272304, 0.544609, 0.0892179, 0.178436, 0.356872, 0.713743, 0.427487, 0.854974, 0.709947, 0.419894, 0.839788]
Cette fois, il ne semble pas se rapprocher progressivement de 1. Afin d'obtenir avec précision un développement temporel à long terme, il est nécessaire d'augmenter la précision.
Comparaison de la vitesse de calcul par plusieurs implémentations
Puisque la mise en œuvre de l'exemple ci-dessus utilise np.append, la vitesse de calcul devient extrêmement lente. Dans mon environnement
mpmath.mp.dps = 100000
nmax = 10000
Lorsqu'il est exécuté avec
$ python baker_mp.py
1.895829439163208 [sec.]
[0.170384, 0.340767, 0.681535, 0.36307, 0.72614, 0.452279, 0.904559, 0.809118, 0.618235, 0.236471, 0.472941, 0.945882, 0.891765, 0.78353, 0.56706, 0.13412, 0.26824, 0.536479, 0.0729585, 0.145917, 0.291834, 0.583668, 0.167335, 0.334671, 0.669341, 0.338683, 0.677365, 0.35473, 0.70946, 0.418921]
Cela prend un peu moins de 2 secondes. Bien que l'ajout soit pratique, il est lent car il s'agit d'une opération qui répète la réacquisition d'un tableau (répétition de malloc et copie). Préparez les séquences nécessaires à l'avance et terminez avec juste l'opération d'affectation. Ça devrait être plus rapide.
Puisque la classe de la matrice est définie dans mpmath, utilisons-la. Cependant, puisqu'il s'agit d'une liste python qui correspond à un tableau unidimensionnel, comparez le cas de list avec le cas de numpy.array. Le code suivant montre uniquement les pièces modifiées
nmax = 10000
web = mpmath.zeros(2*nmax+1,2)
web[0,0] = x
web[0,1] = 0
traj = [mpmath.mpf("0")]*(nmax+1)
traj[0] = x
a = time.time()
for n in range(1,nmax+1):
y = B(x)
web[2*n - 1, 0] = x
web[2*n - 1, 1] = y
traj[n] = y
x = y
if n!=nmax:
web[2*n, 0] = x
web[2*n, 1] = y
print(time.time()-a, "[sec.]")
mpmath.nprint(traj[-30:])
$ python baker_mp1.py
0.355421781539917 [sec.]
[0.170384, 0.340767, 0.681535, 0.36307, 0.72614, 0.452279, 0.904559, 0.809118, 0.618235, 0.236471, 0.472941, 0.945882, 0.891765, 0.78353, 0.56706, 0.13412, 0.26824, 0.536479, 0.0729585, 0.145917, 0.291834, 0.583668, 0.167335, 0.334671, 0.669341, 0.338683, 0.677365, 0.35473, 0.70946, 0.418921]
Ensuite, essayez de couvrir traj avec numpy.array
traj = np.array([mpmath.mpf("0")]*(nmax+1))
dtype devient objet. Tous les chiffres sont affichés à moins que la sortie ne soit renvoyée une fois dans la liste.
mpmath.nprint(traj[-30:].tolist())
Le résultat de l'exécution semble être presque le même.
$ python baker_mp2.py
0.36023831367492676 [sec.]
[0.170384, 0.340767, 0.681535, 0.36307, 0.72614, 0.452279, 0.904559, 0.809118, 0.618235, 0.236471, 0.472941, 0.945882, 0.891765, 0.78353, 0.56706, 0.13412, 0.26824, 0.536479, 0.0729585, 0.145917, 0.291834, 0.583668, 0.167335, 0.334671, 0.669341, 0.338683, 0.677365, 0.35473, 0.70946, 0.418921]
Ensuite, c'est compliqué, mais mettez numpy.array dans la liste.
web = [np.array([mpmath.mpf(0)]*2*(nmax+1)),
np.array([mpmath.mpf(0)]*2*(nmax+1))]
traj = np.array([mpmath.mpf(0)]*(nmax+1))
traj[0] = x
a = time.time()
for n in range(1,nmax+1):
y = B(x)
web[0][2*n - 1] = x
web[1][2*n - 1] = y
traj[n] = y
x = y
if n != nmax:
web[0][2*n] = x
web[1][2*n] = y
print(time.time()-a, "[sec.]")
mpmath.nprint(traj[-30:].tolist())
Le résultat de l'exécution est un peu plus rapide.
$ python baker_mp3.py
0.2923922538757324 [sec.]
[0.170384, 0.340767, 0.681535, 0.36307, 0.72614, 0.452279, 0.904559, 0.809118, 0.618235, 0.236471, 0.472941, 0.945882, 0.891765, 0.78353, 0.56706, 0.13412, 0.26824, 0.536479, 0.0729585, 0.145917, 0.291834, 0.583668, 0.167335, 0.334671, 0.669341, 0.338683, 0.677365, 0.35473, 0.70946, 0.418921]
Avec numpy.array
web = np.zeros((2*nmax+1,2),dtype="object")
web[0,0] = x
web[0,1] = 0
traj = np.zeros(nmax+1, dtype="object")
traj[0] = x
a = time.time()
for n in range(1,nmax+1):
y = B(x)
web[2*n - 1, 0] = x
web[2*n - 1, 1] = y
traj[n] = y
x = y
if n!=nmax:
web[2*n, 0] = x
web[2*n, 1] = y
print(time.time()-a, "[sec.]")
mpmath.nprint(traj[-30:].tolist())
Même si cela ne change pas grand-chose.
$ python baker_mp4.py
0.2934293746948242 [sec.]
[0.170384, 0.340767, 0.681535, 0.36307, 0.72614, 0.452279, 0.904559, 0.809118, 0.618235, 0.236471, 0.472941, 0.945882, 0.891765, 0.78353, 0.56706, 0.13412, 0.26824, 0.536479, 0.0729585, 0.145917, 0.291834, 0.583668, 0.167335, 0.334671, 0.669341, 0.338683, 0.677365, 0.35473, 0.70946, 0.418921]
Sauvegardez les données avec une double précision.
web = np.zeros((2*nmax+1,2),dtype=np.float)
web[0,0] = float(x)
web[0,1] = 0
traj = np.zeros(nmax+1, dtype=np.float)
traj[0] = float(x)
a = time.time()
for n in range(1,nmax+1):
y = B(x)
web[2*n - 1, 0] = float(x)
web[2*n - 1, 1] = float(y)
traj[n] = float(y)
x = y
if n!=nmax:
web[2*n, 0] = float(x)
web[2*n, 1] = float(y)
print(time.time()-a, "[sec.]")
print(traj[-30:])
Il n'y a pratiquement aucun changement de vitesse.
$ python baker_mp5.py
0.3071897029876709 [sec.]
[0.17038373 0.34076746 0.68153493 0.36306985 0.72613971 0.45227941
0.90455883 0.80911766 0.61823531 0.23647062 0.47294124 0.94588249
0.89176498 0.78352995 0.5670599 0.13411981 0.26823961 0.53647923
0.07295846 0.14591691 0.29183383 0.58366766 0.16733532 0.33467064
0.66934128 0.33868256 0.67736512 0.35473024 0.70946048 0.41892095]
Conclusion
J'ai essayé diverses méthodes d'implémentation, mais j'ai conclu que python ne change pas de vitesse à moins que l'ajout ne soit utilisé.
Cette fois, une seule condition initiale a été ciblée, mais l'échantillonnage de la condition initiale est nécessaire lors du calcul des statistiques. Je devrais le faire dans ce cas aussi, mais cette fois je ferai de mon mieux.
référence
- mpmath
- Julia:BigFloats and BigInts
- À propos de l'implémentation interne de Mathematica (valeurs numériques et fonctions associées)
- Introduction à Chaos Mechanical Systems 2e édition Auteur original: Robert L. Devaney (Traducteur: Kenichi Goto, Nouvelle édition Traducteur: Hiroshi Kokubun, Yutaka Ishii, Shunsaku Arai, Masashi Kisaka) Kyoritsu Publishing
Recommended Posts