Nous avons considéré la mémorisation en utilisant Last Fibonacci series comme exemple. Cette fois, nous allons considérer l'implémentation de la mémorisation et la vitesse d'exécution par langage en utilisant le polypole Hermite comme exemple. Les langages manipulés sont python, wolfram language (Mathematica) et julia (en préparation).
Tout d'abord, il est extrêmement préférable d'utiliser l'évaluation des systèmes polynomiaux orthogonaux super-basiques en informatique tels que les polynômes Hermite, car ils sont fournis sous forme de package sans les implémenter par eux-mêmes.
En python
En langue Wolfram
En fait, je ne comprends vraiment pas dans Julia (2020). Je ne l'ai pas vraiment utilisé, mais d'après ce que j'ai vérifié
Il semble qu'un package lié aux polynômes orthogonaux puisse être utilisé dans. Concernant les fonctions spéciales
Semble être disponible. Je ne comprends pas si Julia est devenue la norme de facto de l'informatique en python comme numpy et scipy, mais au moins je regarde Github et il est assez actif, donc j'ai hâte à l'avenir.
Dernière série de Fibonacci, la série de Fibonacci met en cache (mémorisation) les termes évalués une fois, et améliore la vitesse de calcul des séries récursives en éliminant les évaluations en double inutiles. fait. De même, considérons la simulation de la fonction de définition récursive en utilisant le polynôme Hermite comme exemple.
Comme vous pouvez le voir dans https://mathworld.wolfram.com/HermitePolynomial.html, les polynômes Hermite peuvent être définis de différentes manières, mais ici Evaluons la définition suivante, qui peut être implémentée relativement facilement en tant que fonction récursive.
Lorsqu'elle est évaluée symboliquement jusqu'à n = 10 dans le langage wolfram
H[0, x_] = 1;
H[1, x_] = 2x;
H[n_, x_] := H[n, x] = 2x*H[n-1, x] - 2(n-1)H[n-2, x];
res = Table["H_"<>ToString[n]<>"(x)="<>ToString@TeXForm@Expand@H[n, x], {n, 0, 10}]//MatrixForm;
Print@res
Comme vous pouvez le voir immédiatement, dans le polynôme Hermite de degré n, $ H_n (x) $ dépend de $ x
L'évaluation d'un polypoly Hermite simple comme référence avec python (numpy, mpmath) serait la suivante.
import numpy as np
import time
def H(n,x):
if n==0:
return np.array([1]*len(x))
elif n==1:
return 2*x
else:
return 2*x*H(n-1, x) - 2*(n-1)* H(n-2,x)
time0 = time.time()
x = np.linspace(-1, 1, 100)
n=30
y = H(n,x)
y0 = y[len(x)//2]
print("numpy:", time.time()-time0, "sec", "H_%d(0)=%.18e" % (n, y0))
n=20
time0 = time.time()
x = np.linspace(-1, 1, 100)
y = H(n,x)
y0 =y[len(x)//2]
print("numpy:", time.time()-time0, "sec", "H_%d(0)=%.18e" % (n, y0))
import mpmath
mpmath.mp.dps = 100
time0 = time.time()
x = np.array(mpmath.linspace(-1, 1, 100), dtype="object")
y = H(n, x)
y0 = float(y[len(x)//2])
print("mpmath:", time.time()-time0, "sec", "H_%d(0)=%.18e" % (n,y0))
$ python hermite-poly.py
numpy: 10.97638750076294 sec H_30(0)=-2.022226152780266537e+20
numpy: 0.09095430374145508 sec H_20(0)=6.690748809755917969e+11
mpmath: 6.780889987945557 sec H_20(0)=6.690748809755917969e+11
Puisque $ H_n (x) $ est un polymorphisme $ x $ de degré $ n $, nous constatons que $ x $ n'a pas été évalué et que seuls les coefficients doivent être évalués. Si vous utilisez l'expression lambda, vous pouvez mettre en cache $ x $ sans être évalué, donc la vitesse de calcul sera améliorée ... mais ce n'est pas différent de l'implémentation simple mentionnée ci-dessus. Cela est dû au fait que l'expression lambda ne met en cache que l'ordre des opérations de calcul et que les résultats du calcul restent non évalués (un exemple spécifique sera décrit plus loin dans mpmath.polynomial).
D'ailleurs, la mémorisation de l'expression lambda n'améliore pas la vitesse de calcul. Par conséquent, seul le coefficient de $ H_n (x) $ est évalué sans utiliser l'équation lamda. Autrement dit, un polymorphisme de degré $ m $
Est donnée. Une liste de coefficients de longueur $ m $
Alors $ p_m (x) \ pm q_m (x) $ est
Multiplicateur $ x \ fois p_m $
Équivaut à faire une liste de longueur $ m + 1 $. En manipulant cette opération, l'addition / soustraction (division), l'intégration et la différenciation de $ p_m (x) $ peuvent être réalisées, mais cela est suffisant pour évaluer le polynôme Hermite. La méthode de division est entre parenthèses car, à l'exception du spécial $ p_m $, il n'est généralement pas possible d'évaluer exactement les opérations. Si la liste de coefficients L (p_m (x)) est utilisée, la valeur évaluée concrètement peut être mise en cache, on s'attend donc à une amélioration du calcul par mémorisation.
Essayez des implémentations naïves en utilisant python numpy.
import numpy as np
import time
from functools import lru_cache
def HermiteCoeff(n):
if n == 0:
return np.array([1])
elif n == 1:
return np.array([0, 2])
else:
coeff = np.zeros(n+1)
coeff[:n-1] += -2*(n-1)*HermiteCoeff(n-2)
coeff[1:] += 2*HermiteCoeff(n-1)
return coeff
known0={0:np.array([1], dtype=np.int64),
1:np.array([0, 2], dtype=np.int64)}
def HermiteCoeffMemo0(n):
if n in known0:
return known0[n]
else:
coeff = np.zeros(n+1, dtype=np.int64)
coeff[:n-1] += -2*(n-1)*HermiteCoeffMemo0(n-2)
coeff[1:] += 2*HermiteCoeffMemo0(n-1)
known0[n] = coeff
return coeff
known1={0:np.array([1], dtype=o),
1:np.array([0, 2], dtype=object)}
def HermiteCoeffMemo1(n):
if n in known1:
return known1[n]
else:
coeff = np.zeros(n+1, dtype=object)
coeff[:n-1] += -2*(n-1)*HermiteCoeffMemo1(n-2)
coeff[1:] += 2*HermiteCoeffMemo1(n-1)
known1[n] = coeff
return coeff
@lru_cache()
def HermiteCoeffMemo2(n):
if n == 0:
return np.array([1], dtype=object)
elif n == 1:
return np.array([0, 2], dtype=object)
else:
coeff = np.zeros(n+1, dtype=object)
coeff[:n-1] += -2*(n-1)*HermiteCoeffMemo2(n-2)
coeff[1:] += 2*HermiteCoeffMemo2(n-1)
return coeff
def Hermite(coeff):
return lambda x:sum(c*x**i for i, c in enumerate(coeff))
import time
funcs = [HermiteCoeff, HermiteCoeffMemo0, HermiteCoeffMemo1, HermiteCoeffMemo2]
x = np.linspace(-1, 1, 100)
n = 30
for f in funcs:
time0 = time.time()
coeff = f(n)
H = Hermite(coeff)
y = H(x)
y0 = y[len(x)//2]
print("eval. time:%10f[sec]" % (time.time()-time0), "H_%d(0)=%5e" % (n, y0), "by", f.__name__, coeff.dtype)
Le résultat de l'exécution est
eval. time: 11.646824[sec] H_30(0)=-2.022226e+20 by HermiteCoeff float64
eval. time: 0.000543[sec] H_30(0)=7.076253e+16 by HermiteCoeffMemo0 int64
eval. time: 0.000932[sec] H_30(0)=-2.022226e+20 by HermiteCoeffMemo1 object
eval. time: 0.000933[sec] H_30(0)=-2.022226e+20 by HermiteCoeffMemo2 object
Comme prévu, la mémorisation améliore la vitesse, mais il y a quelques points à garder à l'esprit. numpy.zeros renvoie un tableau initialisé à zéro de longueur $ n $, mais avec dtype S'il n'est pas spécifié, le flotteur est rempli de zéro. Dans l'évaluation du polynôme Hermit, l'élément de la liste de coefficients est int, ce qui est inutile. Par conséquent, le résultat de la définition (Memoization) + dtype = np.int est le second. C'est bien que la vitesse soit plus rapide, mais le résultat de $ H_ {30} (0) $ est différent. Cela peut être débordé par article précédent.
Par conséquent, le troisième est de permettre de gérer plusieurs entiers de longueur avec mémorisation et tableau, et à la suite de la définition de dtype = object, il peut être confirmé que la même valeur que la première définition est renvoyée, bien qu'elle soit légèrement plus lente que le cas de int. Le résultat de la définition de la commande $ n = 100 $ est le suivant.
eval. time: 0.001842[sec] H_100(0)=-2.410424e+18 by HermiteCoeffMemo0 int64
eval. time: 0.003694[sec] H_100(0)=3.037263e+93 by HermiteCoeffMemo1 object
eval. time: 0.003563[sec] H_100(0)=3.037263e+93 by HermiteCoeffMemo2 object
Lorsque dtype est défini sur np.int, la valeur est complètement aléatoire, on peut donc imaginer qu'elle déborde toujours. De plus, puisque la vitesse du mécanisme de cache maison est presque la même avec lru_cache, il est préférable d'utiliser le lru_cache existant.
Le résultat de $ n = 250 $ est également inclus pour la comparaison suivante (bien que la partie incorrecte soit de 100 chiffres, elle s'écarte donc de la nouvelle valeur d'environ 0,1 ...).
$ python hermite-poly-coeff.py
eval. time: 0.005363[sec] H_250(0)=0.000000e+00 by HermiteCoeffMemo0 int64
eval. time: 0.012851[sec] H_250(0)=-1.673543e+283 by HermiteCoeffMemo1 object
eval. time: 0.013394[sec] H_250(0)=-1.673543e+283 by HermiteCoeffMemo2 object
polynomial manipulation packages
Dans l'implémentation simple du chapitre précédent, il a été implémenté à la main, mais comme il n'a pas de polyvalence et que la demande d'arithmétique algébrique des polynômes est extrêmement grande, ce serait le redéveloppement des roues qu'il implémenterait par lui-même. J'ai en fait consulté un enseignant Google
En python
Reçu la révélation. Vous pouvez vous attendre à améliorer l'implémentation naïve (mais cela ne fonctionne pas). De plus, sympy et wolfram language qui peuvent effectuer des calculs symboliques guide / PolynomialAlgebra.html) peut être implémenté plus simplement (discuté plus tard).
numpy.polynomial
Même si vous lisez le document numpy.polynomial.polynomial.Polynomial Je ne suis pas sûr,
>>> import numpy as np
>>> from numpy.polynomial.polynomial import Polynomial
>>> f = Polynomial([1,2,3])
Polynomial([1., 2., 3.], domain=[-1, 1], window=[-1, 1])
>>> f.coef
array([1., 2., 3.])
Renvoie une classe polymorphe de $ P (x) = 1 + 2x + 3x ^ 2 $. [1, 2, 3] a été assigné comme un entier pour les éléments de la liste de coefficients, mais chaque élément a une virgule flottante. Comme cela a un débordement potentiel, nous aimerions l'évaluer avec un entier de plusieurs longueurs. Bien que non écrit dans le document
>>> f = Polynomial(np.array([1,2,3], dtype=object))
>>> f.coef
array([1, 2, 3], dtype=object)
Ça à l'air bon. Puisque la mémorisation a la même vitesse même avec un mécanisme de cache fait maison, lru_cache est utilisé.
import numpy as np
from numpy.polynomial.polynomial import Polynomial
from functools import lru_cache
import time
import mpmath
mpmath.mp.dps = 100
@lru_cache(maxsize=1000)
def HermitePoly(n, dtype=None):
if n == 0:
return Polynomial(np.array([1], dtype=dtype))
elif n == 1:
return Polynomial(np.array([0,2], dtype=dtype))
else:
xc = Polynomial(np.array([0,1], dtype=dtype) )
return 2*xc*HermitePoly(n-1) -2*(n-1)*HermitePoly(n-2)
mpmath.mp.dps = 100
x = np.array(mpmath.linspace(1, 1, 100))
n = 250
for dtype in [np.int, np.float, object]:
time0 = time.time()
H = HermitePoly(n, dtype=dtype)
y = H(x)
c = H.coef
print(time.time()-time0, "[sec]", H(0), c.max(), c.dtype)
HermitePoly.cache_clear()
$ python numpy-poly.py
0.14997124671936035 [sec] 5.922810347657652e+296 2.46775722331116e+305 float64
0.15068888664245605 [sec] 5.922810347657652e+296 2.46775722331116e+305 float64
0.14925169944763184 [sec] 5.922810347657652e+296 2.46775722331116e+305 object
Lorsque int est spécifié dans dtype, il est remplacé par float, mais le coefficient de dtype = object reste objecto, mais si vous réfléchissez bien, seul int doit apparaître dans le coefficient, donc c.max () flotte. Le point décimal est que dtype est un objet, mais l'élément réel est arrondi pour flotter. En fait, si vous le portez à n = 260
$ python numpy-poly.py
/home/hanada/Applications/anaconda3/lib/python3.7/site-packages/numpy/polynomial/polynomial.py:788: RuntimeWarning: invalid value encountered in double_scalars
c0 = c[-i] + c0*x
0.14511513710021973 [sec] nan nan float64
0.14671945571899414 [sec] nan nan float64
0.14644980430603027 [sec] nan inf object
Par conséquent, un débordement en virgule flottante double précision se produit. Pour le moment, ce problème peut être résolu en utilisant mpmath
mpmath.mp.dps = 10000
@lru_cache(maxsize=1000)
def HermitePoly1(n, dtype=None):
if n == 0:
return Polynomial(np.array([mpmath.mpf("1")]))
elif n == 1:
return Polynomial(np.array([mpmath.mpf("0"), mpmath.mpf("2")]))
else:
xc = Polynomial(np.array([mpmath.mpf("0"), mpmath.mpf("1")]))
return 2*xc*HermitePoly1(n-1) -2*(n-1)*HermitePoly1(n-2)
x = np.array(mpmath.linspace(1, 1, 100))
n = 250
time0 = time.time()
H = HermitePoly1(n, dtype=dtype)
y = H(x)
c = H.coef
print(time.time()-time0, "[sec]", float(H(0)), c.max(), c.dtype)
Si oui, c'est un peu long
0.30547070503234863 [sec] -1.7171591075700597e+283 4742832273990616849979871677859751938138722866700194098725540674413432950130755339448602832297348736154189959265033164596368647282052560604201151158911400057598325345216981770764974144951365648905162496309065906209718571075651658676878296726900331786598665420800000000000000000000000000000000.0 object
Il se trouve que Cependant, puisque le coefficient est également en notation à virgule flottante, probablement dans ABCPolyBase qui est la classe parente de numpy.polynomial.polynomial. Est-ce l'influence de la mise en œuvre interne de? Je n'ai pas suivi aussi loin.
De plus, étant donné que la vitesse de calcul est évidemment plus rapide si vous l'implémentez vous-même, il vaut mieux éviter d'utiliser numpy.polynomial.polynomial.Polynomial (en commençant par une majuscule !!) dans le calcul de précision longue et multiple.
mpmath.polynomial
En premier lieu, numpy est une base de langage C, vous ne pouvez donc pas en faire trop. Puisque mochi est un magasin de mochi, je vais essayer de l'implémenter dans mpmath (mais cela pose également un problème). Dans mpmath, mpmath (Orthogonal Polygonal System) fournit des polynômes Hermite, vous pouvez donc les utiliser, mais dans des cas généraux Considérez l'extension dans, et voyez si elle peut être simulée en utilisant mpmath.polynomial.
Dans l'ordre inverse de pour numpy, la liste des coefficients Notez que si $ [c_0, c_1, c_2] $ est donné, $ P (x) = c_0x ^ 2 + c_1x ^ 1 + c_2 $.
Évaluons avec 1000 chiffres de la pièce incorrecte.
import mpmath as mp
mp.mp.dps = 1000
import time
Lorsqu'il est mis en œuvre naïvement
def HermitePoly1(n,x):
if n == 0:
return mp.polyval([1], x)
elif n==1:
return mp.polyval([2,0], x)
else:
xc = mp.polyval([1,0], x)
return 2*xc* HermitePoly1(n-1, x) - 2*(n-1)*HermitePoly1(n-2,x)
C'est ça? Il semble qu'un tableau (liste) ne peut pas être assigné à x dans polyval, donc si vous mémorisez simplement
known={0:lambda x: mp.polyval([1], x),
1:lambda x: mp.polyval([2,0], x)}
def HermitePoly2(n):
if n in known:
return known[n]
else:
H = lambda x: 2*mp.polyval([1,0], x) * HermitePoly2(n-1)(x) - 2*(n-1)*HermitePoly2(n-2)(x)
known[n] = H
return H
Qu'il sera
x = mp.linspace(-5,5,100,endpoint=False)
n = 20
time0 = time.time()
y = [HermitePoly1(n, xx) for xx in x]
print(time.time()-time0, "[sec]", HermitePoly1(n, 0))
time0 = time.time()
H = HermitePoly2(n)
y = [H(xx) for xx in x]
print(time.time()-time0, "[sec]", H(0))
Exécutez en tant que. Le résultat est
$ python mpmath-poly.py
17.37535572052002 [sec] 670442572800.0
17.98806405067444 [sec] 670442572800.0
Cela ne peut pas être accéléré. Parce que le cache d'expression lamda n'est que le cache dans l'ordre d'évaluation, Cela est dû au fait que le résultat réel de l'évaluation n'est pas mis en cache. Je n'ai pas tellement utilisé mpmath.Polynomial, donc il y a peut-être un meilleur moyen, mais mpmath.Polynomial est la mémorisation. Je ne pouvais pas penser à un moyen de le faire. Je suis curieux de connaître l'implémentation interne de mp.hermite.
Dans le cas du polynôme Hermite, seul le coefficient doit être évalué plusieurs fois, donc s'il est ad hoc,
from functools import lru_cache
import numpy as np
@lru_cache()
def HermiteCoeff(n):
if n == 0:
return np.array([1], dtype=object )
elif n == 1:
return np.array([0, 2], dtype=object)
else:
coeff = np.zeros(n+1, dtype=object)
#2xH(n-1) - 2*(n-1)*H(n-2)
coeff[:n-1] += -2*(n-1)*HermiteCoeff(n-2)
coeff[1:] += 2*HermiteCoeff(n-1)
return coeff
def Herm(coeff):
return lambda x:sum(c*x**i for i, c in enumerate(coeff))
time0 = time.time()
coeff = HermiteCoeff(250)
H = Herm(coeff)
coeff = coeff.tolist()[::-1]
y = [mp.polyval(coeff, xx) for xx in x]
print(time.time()-time0, "[sec]", float(mp.polyval(coeff, 0)))
Si tu aimes
$ python mpmath-poly.py
0.19260430335998535 [sec] -1.7171591075700597e+283
Ce n'est pas que vous ne pouvez pas mémoriser.
sympy
Ensuite, essayons d'utiliser sympy.
import sympy, mpmath
import time
x = sympy.symbols('x')
def HermitePoly1(n):
if n==0:
return 1
elif n==1:
return x
else:
return 2*x*HermitePoly1(n-1) - 2*(n-1)*HermitePoly1(n-2)
known={0:1, 1:x}
def HermitePoly2(n):
if n in known:
return known[n]
else:
Hn = 2*x*HermitePoly2(n-1) - 2*(n-1)*HermitePoly2(n-2)
known[n] = Hn
return Hn
known1={0:1, 1:x}
def HermitePoly3(n):
if n in known1:
return known1[n]
else:
Hn = 2*x*HermitePoly3(n-1) - 2*(n-1)*HermitePoly3(n-2)
known1[n] = sympy.expand(Hn)
return known1[n]
mpmath.mp.dps = 100
x1 = mpmath.linspace(0, 1, 100)
n = 20
time0 = time.time()
H = HermitePoly1(n)
time1 = time.time()
y = [H.subs(x, i) for i in x1]
print(time.time() - time1, time1 - time0, H.subs(x, 0))
time0 = time.time()
H = HermitePoly2(n)
time1 = time.time()
y = [H.subs(x, i) for i in x1]
print(time.time() - time1, time1 - time0, H.subs(x, 0))
time0 = time.time()
H = HermitePoly3(n)
time1 = time.time()
y = [H.subs(x, i) for i in x1]
print(time.time() - time1, time1 - time0, H.subs(x, 0))
Si vous exécutez en tant que
$ python sympy-poly.py
1.612032413482666 0.1680293083190918 670442572800
1.5776398181915283 0.012787342071533203 670442572800
0.2140791416168213 0.07543253898620605 670442572800
Je peux le mémoriser, mais je n'ai pas tellement utilisé sympy, donc je ne sais pas si c'est le meilleur.
Avec Mathematica (langage wolfram), vous pouvez définir l'expression progressive du polynôme Hermite ainsi que de la série de Fibonacci. Environ 14 secondes pour une mise en œuvre rustique
He[n_, x_] := 2x*He[n - 1, x] - 2 (n - 1)*He[n - 2, x]
He[0, x_] = 1;
He[1, x_] = 2x;
{First@Timing@Expand@He[30, x], "sec"}
Out[4]= {14.214, sec}
Lorsque le résultat du calcul est mis en cache
H[n_, x_] := H[n, x] = 2 x*H[n - 1, x] - 2 (n - 1)*H[n - 2, x]
H[0, x_] = 1;
H[1, x_] = 2x;
{First@Timing@Expand@H[30, x], "sec"}
Out[8]= {7.89031, sec}
Il semble que Expand prenne beaucoup de temps
h[n_, x_] := h[n, x] = 2x*h[n - 1, x] - 2 (n - 1)*h[n - 2, x] // Expand
h[0, x_] = 1;
h[1, x_] = 2x;
{First@Timing@Expand@h[30, x], "sec"}
Out[12]= {0.007158, sec}
C'est encore plus rapide si vous utilisez HermiteH défini dans Mathematica.
In[14]:= {First@Timing@HermiteH[30, x], "sec"}
Out[14]= {0.000216, sec}
Que diable faites-vous à l'intérieur?
Julia
Les opérations liées aux opérations polynomiales en python peuvent être réalisées en utilisant Polynomials.jl.
(Je me fatigue, donc je prépare le programme)
Recommended Posts