Dans un article précédent, j'ai écrit sur l'accélération de NumPy et SciPy:
Accélération du calcul numérique avec NumPy: Basics Accélération du calcul numérique à l'aide de NumPy / SciPy: Application 1 Accélération du calcul numérique à l'aide de NumPy / SciPy: Application 2
Est-ce vraiment rapide? Découvrons-le.
Comparez avec l'implémentation oleore par Python. Ce n'est peut-être pas une évaluation valide car c'est une comparaison avec l'implémentation qui met l'accent sur la simplicité plutôt que sur la vitesse. Python ʻanaconda3, le % timeit` d'IPython pour la mesure du temps Je vais l'utiliser.
Par exemple, l'initialisation de la matrice.
before
N = 3000
a = [[i + j for i in range(N)] for j in range(N)]
Disons que vous l'avez écrit en inclusion comme ceci. Le temps d'exécution est ** 661ms </ font> **. D'un autre côté, dans le ʻarange` de NumPy
after
a = np.arange(0, N) + np.arange(0, N).reshape(N, 1)
** 15,5 ms </ font> **. Ce n'est pas 100 fois plus rapide, mais 40 fois plus rapide.
Il existe de nombreuses opérations algébriques linéaires, mais essayons un produit matriciel simple. Préparez deux matrices 200x200 et calculez le produit:
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)])
Si implémenté avec la définition du produit matriciel
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]
Le temps d'exécution est de ** 7,7s </ font> **.
Utilisation de la fonction universelle de NumPy
after
c = np.dot(a, b)
C'est bien d'être simple avant la vitesse. Le temps d'exécution est ** 202us </ font> **. C'est un moment qu'il est déraisonnable de dire combien de fois. C'est multi-cœur par BLAS de MKL La puissance du traitement. Même avec 1 000 × 1 000 ** 22,2 ms </ font> **. À ce stade, l'implémentation de la boucle for devient ingérable.
Différencions $ \ sin x $. Le nombre de divisions dans l'espace est de 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
Conservez la définition de la différenciation. Ne vous inquiétez pas de la précision pour le moment:
before
diff = []
for i in x:
diff.append((f(i + dx) - f(i)) / dx)
Le temps d'exécution est ** 91,9ms </ font> **. C'est peut-être la raison pour laquelle ʻappend` est lent. D'un autre côté, dans NumPy
after
diff = np.gradient(g(x), dx)
C'est simple et ne réduit pas le nombre d'éléments dans la valeur de retour (combien cela est reconnaissant pour les calculs numériques!). Le temps d'exécution est ** 8,25 ms </ font> **. Lumière C'est environ 10 fois plus rapide.
g = np.vectorize(f)
Vous pouvez également convertir l'argument et la valeur de retour dans la spécification ndarray en faisant.
Préparons une petite intégration de base:
\int_{-\infty}^{\infty} dx\int_{-\infty}^{x} dy\; e^{-(x^2 + y^2)}
Il s'agit d'une double intégration telle que $ x $ est inclus dans l'intervalle d'intégration de $ y $.
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
Si vous implémentez la double intégration telle que définie,
before
ans = 0
for index, X in enumerate(x):
for Y in y[:index+1]:
ans += dx * dy * h(X, Y)
Est-ce que c'est? Le temps d'exécution est ** 5,89 s </ font> **. Mais c'est trop inexact. Si vous l'écrivez vous-même, vous devriez utiliser l'intégration Simpson, mais dans ce cas Le code devient assez compliqué. De plus, même avec l'intégration Simpson, l'intégration au sens large ne peut correspondre qu'à "prendre un espace suffisamment large". Cependant, si la différence est affinée, le temps d'exécution augmente de l'ordre de $ O (N ^ 2) $. D'un autre côté, le dblquad de SciPy résout tout ce problème:
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]
L'intervalle d'intégration de $ x $ est $ [a, b] $ et l'intervalle d'intégration de $ y $ est $ [{\ rm gfun}, {\ rm hfun}] $. Avec l'intégration adaptative, vous pouvez également spécifier la plage d'erreur. La valeur de retour de dblquad est tuple, et il semble que l'erreur absolue soit également renvoyée. Le temps d'exécution est ** 51,1ms </ font> **. Plus rien à dire Non, génial.
(C'est un bonus. La solution numérique de l'équation des valeurs propres est un peu maniaque ...)
Résolvons l'équation de Schroedinger du système d'oscillateur harmonique. Est-ce la méthode Jacobi si nous l'implémentons nous-mêmes? Si vous êtes intéressé, jetez-y un œil:
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)
#Méthode Jacobi
flag = True
while(flag):
#Trouvez le maximum et l'indice des composants hors diagonale
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
#Jugement de convergence
if(maxValue < 1e-4):
flag = False
# print(maxValue)
#Préparation de la matrice de rotation
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)
#Fonctionnement de la matrice
H = np.array(np.matrix(J.T)*np.matrix(H)*np.matrix(J))
I = np.array(np.matrix(I)*np.matrix(J))
#Stockage des valeurs propres et des vecteurs propres
v, w = I.transpose(), []
for i in range(N):
w.append([H[i][i], i])
w.sort()
La convergence se produit lorsque la valeur maximale du terme hors diagonale devient suffisamment petite. Ce n'est pas très pratique car les valeurs propres ne sont pas dans l'ordre croissant. Le temps d'exécution est de ** 15,6 s </ font> **. Il est plus difficile d'augmenter le nombre de divisions. Avec NumPy
after
#Les paramètres du système
L, N = 10, 80
x, dx = np.linspace(-L/2, L/2, N), L / N
#Terme d'exercice 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)
#Terme potentiel
V = np.diag(np.linspace(-L/2, L/2, N)**2)
#Équation des valeurs propres de la matrice Hermeet
#w est une valeur unique,v est le vecteur propre
H = (K + V) / 2
w, v = np.linalg.eigh(H)
Le temps d'exécution est ** 1.03ms </ font> **. Quand il s'agit d'un algorithme aussi compliqué que l'équation des valeurs propres, il n'y a aucune chance de gagner à bien des égards. C'est aussi le déchaînement de BLAS. Si c'est C, vous devrez compter sur Lapack ou GSL, mais c'est aussi difficile.
Cela devrait suffire. ** Vous pouvez voir à quelle vitesse NumPy / SciPy fonctionne et à quel point il peut être écrit. ** Si C / C ++ effectue le calcul ci-dessus sans compter sur une bibliothèque externe Si vous essayez de l'écrire, le code sera tel qu'il est écrit au début de chaque chapitre.De plus, par expérience, il semble que de nombreuses bibliothèques externes de C / C ++ sont difficiles à toucher.
Pour référence seulement, la différence de vitesse est résumée ci-dessous:
Initialiser la liste: ** 661ms </ font> ** → ** 15,5ms </ font> ** Produit matriciel: ** 7.7s </ font> ** → ** 202us </ font> ** Différenciation: ** 91,9 ms </ font> ** → ** 8,25 ms </ font> ** (Double) Intégration: ** 5,89 s </ font> ** → ** 51,1 ms </ font> ** Équation de valeur intrinsèque: ** 15,6 s </ font> ** → ** 1,03 ms </ font> **
** Accablant. ** Je pense que la calculatrice vaut la peine d'utiliser Python.
Recommended Posts