In einem früheren Artikel habe ich über die Beschleunigung von NumPy und SciPy geschrieben:
Beschleunigung der numerischen Berechnung mit NumPy: Basics Beschleunigung der numerischen Berechnung mit NumPy / SciPy: Anwendung 1 Beschleunigung der numerischen Berechnung mit NumPy / SciPy: Anwendung 2
Ist es wirklich schnell? Lass es uns herausfinden.
Vergleichen Sie mit der Oleore-Implementierung von Python. Dies ist möglicherweise keine gültige Bewertung, da es sich um einen Vergleich mit der Implementierung handelt, bei dem eher Einfachheit als Geschwindigkeit im Vordergrund steht. Python ist "anaconda3", IPythons "% timeit" für die Zeitmessung Ich werde es benutzen.
Zum Beispiel Matrixinitialisierung.
before
N = 3000
a = [[i + j for i in range(N)] for j in range(N)]
Angenommen, Sie haben es in einer inklusiven Notation wie dieser geschrieben. Die Ausführungszeit beträgt ** 661 ms </ font> **. Andererseits in NumPys "Bereich"
after
a = np.arange(0, N) + np.arange(0, N).reshape(N, 1)
** 15,5 ms </ font> **. Es ist nicht 100-mal schneller, aber 40-mal schneller.
Es gibt viele lineare algebraische Operationen, aber versuchen wir es mit einem einfachen Matrixprodukt. Bereiten Sie zwei 200x200-Matrizen vor und berechnen Sie das Produkt:
setting
import numpy as np
from numpy.random import rand
N = 200
a = np.array(rand(N, N))
b = np.array(rand(N, N))
c = np.array([[0] * N for _ in range(N)])
Wenn mit der Definition des Matrixprodukts implementiert
befor
for i in range(N):
for j in range(N):
for k in range(N):
c[i][j] = a[i][k] * b[k][j]
Die Ausführungszeit beträgt ** 7,7 s </ font> **.
Verwendung der universellen Funktion von NumPy
after
c = np.dot(a, b)
Es ist gut, vor der Geschwindigkeit einfach zu sein. Die Ausführungszeit beträgt ** 202us </ font> **. Es ist ein Moment, in dem es unvernünftig ist zu sagen, wie oft. Dies ist Multi-Core von BLAS von MKL Die Verarbeitungsleistung. Selbst mit 1000 × 1000 ** 22,2 ms </ font> **. Zu diesem Zeitpunkt ist die Implementierung der for-Schleife nicht mehr verwaltbar.
Unterscheiden wir $ \ sin x $. Die Anzahl der Teilungen im Raum beträgt 100000:
setting
import math as m
def f(x):
return m.sin(x)
def g(x):
return np.sin(x)
N, L = 100000, 2 * m.pi
x, dx = np.linspace(0, L, N), L / N
Behalten Sie die Definition der Differenzierung bei. Machen Sie sich zu diesem Zeitpunkt keine Sorgen um die Genauigkeit:
before
diff = []
for i in x:
diff.append((f(i + dx) - f(i)) / dx)
Die Ausführungszeit beträgt ** 91,9 ms </ font> **. Vielleicht ist dies der Grund, warum "Anhängen" langsam ist. Andererseits in NumPy
after
diff = np.gradient(g(x), dx)
Es ist einfach und reduziert nicht die Anzahl der Elemente im Rückgabewert (wie dankbar dies für numerische Berechnungen ist!). Die Ausführungszeit beträgt ** 8,25 ms </ font> **. Licht Es ist ungefähr zehnmal schneller. Übrigens
g = np.vectorize(f)
Sie können auch die Argumente konvertieren und Werte in die ndarray-Spezifikation zurückgeben.
Bereiten wir eine kleine Kernintegration vor:
\int_{-\infty}^{\infty} dx\int_{-\infty}^{x} dy\; e^{-(x^2 + y^2)}
Es ist eine doppelte Integration, so dass $ x $ im Integrationsintervall von $ y $ enthalten ist. Der Integrand und die Anzahl der Teilungen im Raum
setting
def h(x, y):
return np.exp(-(x**2 + y**2))
N, L = 2000, 100
x, dx = np.linspace(-L/2, L/2, N), L / N
y, dy = np.linspace(-L/2, L/2, N), L / N
Wenn Sie die Doppelintegration wie definiert implementieren,
before
ans = 0
for index, X in enumerate(x):
for Y in y[:index+1]:
ans += dx * dy * h(X, Y)
Die Ausführungszeit beträgt ** 5.89s </ font> **. Dies ist jedoch zu ungenau. Wenn Sie es selbst schreiben, sollten Sie die Simpson-Integration verwenden, aber in diesem Fall Der Code wird ziemlich kompliziert. Selbst bei der Simpson-Integration kann die Wide-Sense-Integration nur "ausreichend großen Raum einnehmen" entsprechen. Wenn der Unterschied jedoch fein gemacht wird, erhöht sich die Ausführungszeit in der Größenordnung von $ O (N ^ 2) $ Auf der anderen Seite löst SciPys "dblquad" all dieses Problem:
after
from scipy.integrate import dblquad
ans = dblquad(h, a = -np.inf, b = np.inf, gfun = lambda x : -np.inf, hfun = lambda x : x)[0]
Das Integrationsintervall von $ x $ ist $ [a, b] $ und das Integrationsintervall von $ y $ ist $ [{\ rm gfun}, {\ rm hfun}] $. Bei der adaptiven Integration können Sie auch den Fehlerbereich angeben. Der Rückgabewert von "dblquad" ist Tupel, und es scheint, dass auch der absolute Fehler zurückgegeben wird. Die Ausführungszeit beträgt ** 51,1 ms </ font> **. Mehr nicht zu sagen Großartig.
(Es ist ein Bonus. Die numerische Lösung der Eigenwertgleichung ist ein wenig verrückt ...)
Lösen wir die Schrödinger-Gleichung des harmonischen Oszillatorsystems. Ist es die Jacobi-Methode, wenn wir sie selbst implementieren? Wenn Sie interessiert sind, schauen Sie sich Folgendes an:
before
I = np.eye(N)
H = [[0 for i in range(N)] for j in range(N)]
for i in range(N):
H[i][i] = 2.0/(2*dx**2) + 0.5*(-L/2+dx*i)**2
if(0 <= i+1 < N):
H[i][i+1] = -1.0/(2*dx**2)
if(0 <= i-1 < N):
H[i][i-1] = -1.0/(2*dx**2)
H = np.array(H)
#Jacobi-Methode
flag = True
while(flag):
#Finden Sie das Maximum und den Index von nicht diagonalen Komponenten
maxValue = 0
cI, rI = None, None
for j in range(N):
for i in range(j):
if(maxValue < abs(H[i][j])):
maxValue = abs(H[i][j])
rI, cI = i, j
#Konvergenzurteil
if(maxValue < 1e-4):
flag = False
# print(maxValue)
#Vorbereitung der Rotationsmatrix
theta = None
if(H[cI][cI] == H[rI][rI]):
theta = m.pi/4
else:
theta = 0.5*m.atan(2.0*H[rI][cI]/(H[cI][cI]-H[rI][rI]))
J = np.eye(N)
J[rI][rI] = m.cos(theta)
J[cI][cI] = m.cos(theta)
J[rI][cI] = m.sin(theta)
J[cI][rI] = -m.sin(theta)
#Matrixbetrieb
H = np.array(np.matrix(J.T)*np.matrix(H)*np.matrix(J))
I = np.array(np.matrix(I)*np.matrix(J))
#Speicherung von Eigenwerten und Eigenvektoren
v, w = I.transpose(), []
for i in range(N):
w.append([H[i][i], i])
w.sort()
Es konvergiert, wenn der Maximalwert des nicht diagonalen Terms ausreichend klein wird. Dies ist nicht sehr praktisch, da die Eigenwerte nicht in aufsteigender Reihenfolge vorliegen. Die Ausführungszeit beträgt ** 15,6 s </ font> ** Es ist schwierig, die Anzahl der Abteilungen nicht mehr zu erhöhen. Mit NumPy
after
#Systemeinstellungen
L, N = 10, 80
x, dx = np.linspace(-L/2, L/2, N), L / N
#Übungsbegriff K.
K = np.eye(N, N)
K_sub = np.vstack((K[1:], np.array([0] * N)))
K = dx**-2 * (2 * K - K_sub - K_sub.T)
#Mögliche Laufzeit
V = np.diag(np.linspace(-L/2, L/2, N)**2)
#Hermeet-Matrix-Eigenwertgleichung
#w ist ein eindeutiger Wert,v ist der Eigenvektor
H = (K + V) / 2
w, v = np.linalg.eigh(H)
Die Ausführungszeit beträgt ** 1,03 ms </ font> **. Wenn es um einen Algorithmus geht, der so kompliziert ist wie die Eigenwertgleichung, gibt es in vielerlei Hinsicht keine Gewinnchance. Wenn es C ist, müssen Sie sich auf Lapack oder GSL verlassen, aber das ist auch schwierig.
Dies sollte ausreichen. ** Sie können sehen, wie schnell NumPy / SciPy funktioniert und wie einfach es geschrieben werden kann. ** Wenn C / C ++ die obige Berechnung durchführt, ohne auf eine externe Bibliothek angewiesen zu sein Wenn Sie versuchen, es einzuschreiben, entspricht der Code dem am Anfang jedes Kapitels beschriebenen Code. Erfahrungsgemäß sind viele externe Bibliotheken von C / C ++ schwer zu berühren.
Nur als Referenz ist die Geschwindigkeitsdifferenz unten zusammengefasst:
Liste initialisieren: ** 661 ms </ font> ** → ** 15,5 ms </ font> ** Matrixprodukt: ** 7.7s </ font> ** → ** 202us </ font> ** Unterscheidung: ** 91,9 ms </ font> ** → ** 8,25 ms </ font> ** (Doppelte) Integration: ** 5,89s </ font> ** → ** 51,1 ms </ font> ** Eigenwertgleichung: ** 15,6 s </ font> ** → ** 1,03 ms </ font> **
** Überwältigend. ** Ich denke, der Taschenrechner ist es wert, Python zu verwenden.