Wir haben Memoization am Beispiel von Last Fibonacci series betrachtet. Dieses Mal werden wir die Implementierung von Memoization und die Ausführungsgeschwindigkeit nach Sprache am Beispiel des Hermite-Polypolys betrachten. Die behandelten Sprachen sind Python, Wolfram (Mathematica) und Julia (in Vorbereitung).
Zuallererst ist es überwältigend besser, die Bewertung von super-grundlegenden orthogonalen Polynomsystemen in der Informatik wie Hermite-Polynomen zu verwenden, da sie als Paket bereitgestellt werden, ohne sie selbst zu implementieren.
In Python
In Wolfram-Sprache
Eigentlich verstehe ich in Julia (2020) wirklich nicht. Ich habe es nicht wirklich benutzt, aber soweit ich es überprüft habe
Es scheint, dass ein Paket, das sich auf orthogonale Polynome bezieht, in verwendet werden kann. In Bezug auf spezielle Funktionen
Scheint verfügbar zu sein. Ich verstehe nicht, ob Julia zum De-facto-Standard der Informatik in Python geworden ist, wie Numpy und Scipy, aber zumindest schaue ich auf Github und es ist ziemlich aktiv, also freue ich mich auf die Zukunft.
Letzte Fibonacci-Reihe, die Fibonacci-Reihe speichert die einmal ausgewerteten Begriffe zwischen (Memoisierung) und verbessert die Berechnungsgeschwindigkeit rekursiver Reihen, indem unnötige doppelte Auswertungen vermieden werden. tat. Betrachten Sie in ähnlicher Weise die Simulation der rekursiven Definitionsfunktion am Beispiel des Hermite-Polynoms.
Wie Sie in https://mathworld.wolfram.com/HermitePolynomial.html sehen können, können Hermite-Polynome auf verschiedene Arten definiert werden, aber hier Lassen Sie uns die folgende Definition bewerten, die relativ einfach als rekursive Funktion implementiert werden kann.
Bei symbolischer Auswertung bis zu n = 10 in der Wolfram-Sprache
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
Wie Sie sofort sehen können, ist Memoization wie Last Fibonacci series, da $ H_n (x) $ von $ x $ im Hermite-Polynom n-ter Ordnung abhängt. Dann hat es keinen Sinn.
Auch der Wert am Ursprung der geraden Bestellung ist
Die Bewertung eines einfachen Hermite-Polypolys als Referenz mit Python (numpy, mpmath) wäre wie folgt.
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
Da $ H_n (x) $ ein $ x $ -Polymorphismus vom Grad $ n $ ist, stellen wir fest, dass $ x $ nicht ausgewertet wurde und nur die Koeffizienten ausgewertet werden müssen. Wenn Sie den Lambda-Ausdruck verwenden, können Sie $ x $ zwischenspeichern, ohne ausgewertet zu werden, sodass die Berechnungsgeschwindigkeit verbessert wird. Dies entspricht jedoch der Effizienz der oben genannten einfachen Implementierung. Dies liegt daran, dass der Lambda-Ausdruck nur die Reihenfolge der Berechnungsoperationen zwischenspeichert und die Berechnungsergebnisse nicht bewertet werden (ein spezielles Beispiel wird später in mpmath.polynomial beschrieben).
Übrigens führt das Auswendiglernen des Lambda-Ausdrucks nicht zu einer Verbesserung der Berechnungsgeschwindigkeit. Daher wird nur der Koeffizient von $ H_n (x) $ ohne Verwendung der Lamda-Gleichung ausgewertet. Das heißt, ein Polymorphismus vom Grad $ m $
Ist gegeben. Eine Liste von Koeffizienten der Länge $ m $
Dann ist $ p_m (x) \ pm q_m (x) $
Multiplikator $ x \ times p_m $
Entspricht einer Liste mit der Länge $ m + 1 $. Durch Manipulieren dieser Operation kann Addition / Subtraktion (Division), Integration und Differenzierung von $ p_m (x) $ realisiert werden, dies reicht jedoch zur Bewertung des Hermite-Polynoms aus. Die Teilungsmethode steht in Klammern, da es mit Ausnahme des speziellen $ p_m $ im Allgemeinen nicht möglich ist, die Operation genau auszuwerten. Wenn die Koeffizientenliste L (p_m (x)) verwendet wird, kann der konkret ausgewertete Wert zwischengespeichert werden, so dass erwartet wird, dass eine Verbesserung der Berechnung durch Memoisierung erwartet werden kann.
Probieren Sie einige naive Implementierungen mit Python Numpy aus.
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)
Das Ausführungsergebnis ist
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
Wie erwartet verbessert Memoization die Geschwindigkeit, aber es gibt einige Dinge zu beachten. numpy.zeros gibt ein nullinitialisiertes Array mit der Länge $ n $ zurück, jedoch mit dtype Wenn nicht angegeben, wird der Float mit Null gefüllt. Bei der Auswertung des Hermit-Polynoms ist das Element der Koeffizientenliste int, daher ist es nutzlos. Daher ist das Ergebnis der Einstellung (Memoization) + dtype = np.int das zweite. Es ist gut, dass die Geschwindigkeit schneller ist, aber das Ergebnis von $ H_ {30} (0) $ ist anders. Dies kann durch vorherigen Artikel überfüllt sein.
Daher besteht die dritte Möglichkeit darin, Ganzzahlen mit mehreren Längen mit Memoization und Array zu verarbeiten, sodass bestätigt werden kann, dass das Ergebnis der Einstellung von dtype = object etwas langsamer ist als der Fall von int, aber denselben Wert wie die erste Definition zurückgibt. Das Ergebnis der Einstellung der Reihenfolge $ n = 100 $ ist wie folgt.
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
Wenn dtype auf np.int gesetzt ist, ist der Wert völlig zufällig, sodass man sich vorstellen kann, dass er immer noch überläuft. Da die Geschwindigkeit des hausgemachten Cache-Mechanismus bei lru_cache fast gleich ist, ist es besser, den vorhandenen lru_cache zu verwenden.
Das Ergebnis von $ n = 250 $ wird ebenfalls als nächstes zum Vergleich herangezogen (obwohl der falsche Teil 100 Stellen hat und daher um etwa 0,1 ... vom neuen Wert abweicht).
$ 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
In der einfachen Implementierung des vorherigen Kapitels wurde es von Hand implementiert, aber da es keine Vielseitigkeit aufweist und der Bedarf an algebraischer Arithmetik von Polynomen extrem groß ist, wäre es die Neuentwicklung von Rädern, die selbst implementiert werden müssen. Ich habe mich tatsächlich mit Google Teacher beraten
In Python
Erhielt die Offenbarung. Sie können erwarten, die naive Implementierung zu verbessern (aber es funktioniert nicht). Außerdem sympy und wolfram language, die symbolische Berechnungen ermöglichen guide / PolynomialAlgebra.html) kann einfacher implementiert werden (wird später erläutert).
numpy.polynomial
Auch wenn Sie das Dokument numpy.polynomial.polynomial.Polynomial lesen Ich bin mir nicht sicher,
>>> 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.])
Gibt eine polymorphe Klasse von $ P (x) = 1 + 2x + 3x ^ 2 $ zurück. [1,2,3] wurde als Ganzzahl für die Elemente der Koeffizientenliste zugewiesen, aber jedes Element hat einen Gleitkomma. Da dies einen potenziellen Überlauf hat, möchten wir ihn mit einer Ganzzahl mit mehreren Längen bewerten. Obwohl nicht im Dokument geschrieben
>>> f = Polynomial(np.array([1,2,3], dtype=object))
>>> f.coef
array([1, 2, 3], dtype=object)
Es sieht gut aus. Da die Memoisierung auch mit einem selbst erstellten Cache-Mechanismus dieselbe Geschwindigkeit aufweist, wird lru_cache verwendet.
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
Wenn int in dtype angegeben wird, wird es von float überschrieben, aber der Koeffizient von dtype = object bleibt objecto. Wenn Sie jedoch sorgfältig überlegen, sollte nur int im Koeffizienten erscheinen, sodass c.max () float. Der Dezimalpunkt ist, dass dtype ein Objekt ist, das tatsächliche Element jedoch auf float gerundet ist. In der Tat, wenn Sie es auf n = 260 erhöhen
$ 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
Daher tritt ein Gleitkommaüberlauf mit doppelter Genauigkeit auf. Dieses Problem kann vorerst mit mpmath gelöst werden
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)
Wenn ja, ist es etwas lang
0.30547070503234863 [sec] -1.7171591075700597e+283 4742832273990616849979871677859751938138722866700194098725540674413432950130755339448602832297348736154189959265033164596368647282052560604201151158911400057598325345216981770764974144951365648905162496309065906209718571075651658676878296726900331786598665420800000000000000000000000000000000.0 object
Es stellt sich heraus, dass Da der Koeffizient jedoch auch in Gleitkomma-Notation vorliegt, wahrscheinlich in ABCPolyBase, der übergeordneten Klasse von numpy.polynomial.polynomial. Ist es der Einfluss der internen Umsetzung von? Ich bin nicht so weit gefolgt.
Da die Berechnungsgeschwindigkeit offensichtlich höher ist, wenn Sie sie selbst implementieren, ist es besser, die Verwendung von numpy.polynomial.polynomial.Polynomial (beginnend mit einem Großbuchstaben !!) bei der mehrfachen Berechnung mit langer Genauigkeit zu vermeiden.
mpmath.polynomial
In erster Linie ist numpy eine c-Sprachbasis, sodass Sie es nicht übertreiben können. Da Mochi ein Mochi-Shop ist, werde ich versuchen, es in mpmath zu implementieren (aber dies hat auch ein Problem). In mpmath stellt mpmath (Orthogonales polygonales System) Hermite-Polynome bereit, sodass Sie sie verwenden können, aber in allgemeinen Fällen Betrachten Sie die Erweiterung in und prüfen Sie, ob sie mit mpmath.polynomial simuliert werden kann.
In umgekehrter Reihenfolge von für numpy die Koeffizientenliste Beachten Sie, dass wenn $ [c_0, c_1, c_2] $ angegeben ist, $ P (x) = c_0x ^ 2 + c_1x ^ 1 + c_2 $.
Lassen Sie uns mit 1000 Ziffern des falschen Teils auswerten.
import mpmath as mp
mp.mp.dps = 1000
import time
Bei naiver Umsetzung
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)
Ist es? Es scheint, dass ein Array (eine Liste) x in polyval nicht zugewiesen werden kann. Wenn Sie also einfach auswendig lernen
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
Wird es sein
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))
Ausführen als. Ergebnis ist
$ python mpmath-poly.py
17.37535572052002 [sec] 670442572800.0
17.98806405067444 [sec] 670442572800.0
Es kann nicht beschleunigt werden. Da der Lamda-Ausdrucks-Cache nur der Cache in der Reihenfolge der Auswertung ist, Dies liegt daran, dass das tatsächliche Bewertungsergebnis nicht zwischengespeichert wird. Ich habe mpmath.Polynomial nicht so oft verwendet, daher gibt es vielleicht einen besseren Weg, aber mpmath.Polynomial ist die Memoisierung. Ich konnte mir keinen Weg vorstellen, es zu tun. Ich bin gespannt auf die interne Implementierung von mp.hermite.
Im Fall der Hermite-Polypolyse muss nur der Koeffizient mehrmals ausgewertet werden, wenn er also ad hoc ist
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)))
Wenn du magst
$ python mpmath-poly.py
0.19260430335998535 [sec] -1.7171591075700597e+283
Es ist nicht so, dass man sich nichts merken kann.
sympy
Als nächstes versuchen wir, Sympy zu verwenden.
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))
Wenn Sie als laufen
$ python sympy-poly.py
1.612032413482666 0.1680293083190918 670442572800
1.5776398181915283 0.012787342071533203 670442572800
0.2140791416168213 0.07543253898620605 670442572800
Ich kann es mir merken, aber ich habe Sympy nicht so oft benutzt, deshalb bin ich mir nicht sicher, ob dies das Beste ist.
Mit Mathematica (Wolfram-Sprache) können Sie den allmählichen Ausdruck des Hermite-Polynoms sowie der Fibonacci-Reihe definieren. Ungefähr 14 Sekunden für eine rustikale Implementierung
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}
Wenn das Berechnungsergebnis zwischengespeichert wird
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}
Es scheint, dass Expand lange dauert
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}
Es ist sogar noch schneller, wenn Sie HermiteH verwenden, das in Mathematica definiert ist.
In[14]:= {First@Timing@HermiteH[30, x], "sec"}
Out[14]= {0.000216, sec}
Was um alles in der Welt machst du drinnen?
Julia
Operationen im Zusammenhang mit Polynomoperationen in Python können mithilfe von Polynomials.jl realisiert werden.
(Ich werde müde, also bereite ich das Programm vor)
Recommended Posts