Accélérer les calculs numériques avec NumPy: principes de base

Cible

Il est écrit pour les débutants de Python et NumPy. Il peut être particulièrement utile pour ceux qui visent le calcul numérique de systèmes physiques correspondant à "Je peux utiliser le langage C mais j'ai récemment commencé Python" ou "Je ne comprends pas comment écrire comme Python". Peut être.

De plus, il peut y avoir des descriptions incorrectes en raison de mon manque d'étude. Veuillez me pardonner.

Résumé

Python possède une bibliothèque très riche de calculs numériques et est facile à utiliser. Beaucoup d'entre eux sont des wrappers d'héritage écrits en Fortran, mais ils sont beaucoup plus faciles que de les appeler depuis C / C ++ etc. Par exemple, si vous appelez LAPACK [^ 1] à partir de C / C ++ et essayez de calculer le problème de valeur propre

info = LAPACKE_zheevr(LAPACK_COL_MAJOR, 'V', 'A', 'L', H, (lapack_complex_double*)a, lda, vl, vu, il, iu, abstol, &mm, w, (lapack_complex_double*)z, ldz, isuppz);

Il y a trop d'arguments et c'est compliqué. Bien que plus il y a d'arguments, plus il y a de liberté de calcul, cela ne ressemble pas à un débutant. De plus, C / C ++ ne supporte pas nativement les complexes. Depuis, des types tels que lapack_complex_double sont définis à l'intérieur du wrapper, ce qui peut être déroutant. D'un autre côté, c'est très concis en Python:

w, v = numpy.linalg.eig(H)

Les arguments comme ceux de LAPACK sont également préparés correctement [^ 2], mais la valeur par défaut est définie, donc si vous voulez simplement le résoudre sans penser aux détails, vous pouvez le faire. Simple et bon À propos, divers solveurs semblent avoir des performances plus élevées dans SciPy que dans NumPy.

Lorsque cela se produit, il devient tentant de laisser divers calculs numériques à Python, mais cet enfant présente également des inconvénients. ** Les boucles For sont extrêmement lentes. ** Dans les opérations dites matricielles, des triples boucles se produisent fréquemment. Cependant, le calcul numérique est fondamentalement "finement différencié et bouclé", donc c'est assez fatal [^ 3]. Bien qu'il ne s'agisse pas d'un langage de compilation, il semble être plus lent que les autres langages LL Pour les opérations matricielles

for i in range(N):
	for j in range(N):
		for k in range(N):
			c[i][j] = a[i][k] * b[k][j]

Je l'ai écrit et c'est fini.

Fonctions et opérations intégrées de NumPy

Donc, ce qu'il faut faire est de tirer parti de la fonction universelle de NumPy. ** La fonction intégrée de NumPy fonctionne sur chaque élément du tableau NumPy ndarray et renvoie un nouveau ndarray. . ** Ceci est très important et a le potentiel de chasser les boucles pour les opérations tensorielles. Par exemple, math.sin (cmath.sin) ne peut transmettre que des scalaires, mais numpy.sin

>>> theta = numpy.array([numpy.pi/6, numpy.pi/4, numpy.pi])
>>> numpy.sin(theta) 
array([  5.00000000e-01,   7.07106781e-01,   1.22464680e-16])

Lancer un tableau NumPy fait agir numpy.sin sur chaque élément, et souvent ce type de traitement est aussi rapide que C / C ++. L'implémentation principale de NumPy est C. C'est parce que BLAS [^ 4] / LAPACK fait de son mieux dans les opérations linéaires.

De plus, les quatre règles de fonctionnement définies par «ndarray» sont également rapides et diverses.

Partie 1: vecteur

Écrivons une opération linéaire liée à un vecteur sans boucle for. Comme on s'attend à ce que des fonctions telles que le produit interne et le produit externe soient correctement préparées, comment la somme, le produit, le quotient, etc. pour ndarray sont-ils principalement décrits ci-dessous? Voyons s'il est défini. ** L'opération sur ndarray est complètement différente de l'opération sur liste construite en Python ** Après cela, les alias sont définis comme ʻimport numpy as np`.

Vecteur et scalaire

Un "ndarray" unidimensionnel peut être considéré comme une sorte de vecteur. Scalar (complex, float, int) fonctionne sur tous les éléments de ndarray:

>>> a = np.array([1, 2, 3])
>>> a + 2
array([3, 4, 5])
>>> a * (3 + 4j)
array([ 3. +4.j,  6. +8.j,  9.+12.j])

Vecteurs

ndarray fonctionne avec des éléments avec des index correspondants:

>>> a = np.array([1, 2, 3])
>>> b = -np.array([1, 2, 3])
>>> a + b
array([0, 0, 0])
>>> a * b
array([-1, -4, -9])
>>> a / b
array([-1., -1., -1.])

Il est intéressant de noter qu'il ne s'agit ni d'un produit intérieur ni d'un produit extérieur. Il est possible de définir une division paire. Il est souvent utilisé lors de l'exécution d'une opération différentielle sur une certaine fonction différenciée. La somme / le produit de «ndarray» avec des dimensions différentes n'est pas défini.

Ndarrays remodelés

ndarray a une méthode appelée reshape, par exemple, vous pouvez réorganiser un vecteur 3x1 en 1x3. C'est comme un vecteur horizontal pour un vecteur vertical, mais un produit (*). Il n'y a pas de correspondance simple car la définition de est différente de celle de l'espace vectoriel. Le produit avec ce vecteur reshape est très intéressant:

>>> c = a.reshape(3, 1)
>>> c
array([[1],
       [2],
       [3]])
>>> c + b
array([[ 0, -1, -2],
       [ 1,  0, -1],
       [ 2,  1,  0]])
>>> c * b
array([[-1, -2, -3],
       [-2, -4, -6],
       [-3, -6, -9]])

Tout en ayant des règles comme (3 × 1) + ou * (1 × 3) = (3 × 3), mais comme l'opération précédente, ** cette opération est réalisable **. Cela semble étrange à première vue. C'est possible, mais si vous regardez, vous pouvez voir que cela peut être expliqué en utilisant les règles de la section précédente. Avec cela, il est possible que même l'initialisation de la matrice ** n'ait pas besoin d'utiliser la boucle for. Masu **. Il est très utile pour des tâches telles que l'initialisation d'une énorme matrice à chaque fois que les paramètres sont modifiés et la poursuite du calcul.

Personnellement, j'aime les spécifications très flexibles qui ne sont pas liées à ces types [^ 5].

Initialisation vectorielle

Ce n'est pas une opération, mais il est inutile d'utiliser une boucle for pour préparer un vecteur en premier lieu. Il existe de nombreuses fonctions qui créent ndarray, mais voici quelques-unes de celles que j'utilise le plus souvent:

>>> L = 1
>>> N = 10
>>> np.linspace(-L/2, L/2, N)
array([-0.5       , -0.38888889, -0.27777778, -0.16666667, -0.05555556,
        0.05555556,  0.16666667,  0.27777778,  0.38888889,  0.5       ])

>>> dL = 0.2
>>> np.arange(-L/2, L/2, dL)
array([-0.5, -0.3, -0.1,  0.1,  0.3])

>>> np.logspace(0, 1 ,10, base=np.e)
array([ 1.        ,  1.11751907,  1.24884887,  1.39561243,  1.5596235 ,
        1.742909  ,  1.94773404,  2.17662993,  2.43242545,  2.71828183])

Si vous combinez des opérations avec des scalaires et des vecteurs, cela semble tout à fait le cas même sans boucle for.

Partie 2: Matrice

Le déroulement de l'histoire est presque le même que celui du vecteur. Puisque le produit matriciel a une fonction dédiée, ce qui suit est aussi l'histoire de l'opération. Je pense que le résultat de l'opération peut être prédit.

Matrice et scalaire

>>> a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
>>> a
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
>>> a + (2 - 3j)
array([[  3.-3.j,   4.-3.j,   5.-3.j],
       [  6.-3.j,   7.-3.j,   8.-3.j],
       [  9.-3.j,  10.-3.j,  11.-3.j]])
>>> a * 0.2
array([[ 0.2,  0.4,  0.6],
       [ 0.8,  1. ,  1.2],
       [ 1.4,  1.6,  1.8]])

Matrice et vecteur

>>> b = np.array([1, 2, 3])
>>> a + b
array([[ 2,  4,  6],
       [ 5,  7,  9],
       [ 8, 10, 12]])
>>> a * b
array([[ 1,  4,  9],
       [ 4, 10, 18],
       [ 7, 16, 27]])
>>> a / b
array([[ 1. ,  1. ,  1. ],
       [ 4. ,  2.5,  2. ],
       [ 7. ,  4. ,  3. ]])

Matrice et matrice

>>> b = np.array([[-1, 2, -3], [4, -5, 6], [-7, 8, -9]])
>>> a + b
array([[ 0,  4,  0],
       [ 8,  0, 12],
       [ 0, 16,  0]])
>>> a * b
array([[ -1,   4,  -9],
       [ 16, -25,  36],
       [-49,  64, -81]])

Initialisation de la matrice

Il existe également des fonctions utiles pour les matrices. En voici quelques-unes:

>>> np.identity(4)
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  1.]])

>>> np.eye(4, 3)
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.],
       [ 0.,  0.,  0.]])

>>> c = np.array([1, 2, 3, 4])
>>> np.diag(c)
array([[1, 0, 0, 0],
       [0, 2, 0, 0],
       [0, 0, 3, 0],
       [0, 0, 0, 4]])

en conclusion

Je pense que d'autres personnes ont écrit de meilleurs articles sur les fonctions pratiques et les règles arithmétiques détaillées de NumPy. Si vous regardez cet article et que vous l'écrivez comme NumPy, vous pouvez faire divers calculs sans utiliser de boucles for. Nous espérons que vous vous sentirez capable de le faire. De nombreux algorithmes manuels sont déjà disponibles dans SciPy et autres.

Cela fait un peu longtemps, je voudrais donc résumer l'application spécifique au calcul numérique dans un autre article.

** Addendum (2016/11/25): Suite ci-dessous ** 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

[^ 1]: LAPACK (Linear Algebra PACKage) est une bibliothèque d'opérations d'algèbre linéaire écrites en Fortran 90. Certains wrappers pour C / C ++ sont appelés Lapacke. Parce qu'il a une interface compatible MKL, Le manuel MKL vous sera utile. [^ 2]: Pour plus d'informations, voir NumPy et [SciPy](https://docs.scipy.org/doc/scipy-0.18. Veuillez vous référer au manuel 1 / reference /). [^ 3]: A côté de ça, j'ai l'impression que la surcharge de l'appel de fonction est également importante. Dans le cas de Procon, il arrive souvent que "l'écriture en méthode de planification dynamique passe, mais en récurrence mémorielle TLE". [^ 4]: Une bibliothèque qui contient un ensemble d'opérations linéaires plus basique que LAPACK. Dans LAPACK, nous appelons BLAS. [^ 5]: Certaines personnes n'aiment pas "aucun type". Contrairement aux programmes conçus pour répondre à "certaines spécifications", les calculs numériques en physique sont si faciles à vérifier si la sortie est correcte. Ce n'est pas le cas. Et il est encore plus difficile de déboguer si "le type de variable change implicitement" est autorisé pendant l'exécution du programme. Puisque le montant observable est fixé à float, le typage statique crée des bogues. Je comprends l'affirmation selon laquelle c'est difficile, mais j'aime Python.

Recommended Posts

Accélérer les calculs numériques avec NumPy: principes de base
Accélération du calcul numérique avec NumPy / SciPy: Application 2
Accélération du calcul numérique avec NumPy / SciPy: Application 1
Accélération du calcul numérique avec NumPy / SciPy: ramasser les oreilles tombées
Découvrez la puissance de l'accélération avec NumPy / SciPy
Principes de base de NumPy
Les bases de #Python (#Numpy 1/2)
Les bases de #Python (#Numpy 2/2)
Principes de base de Python #Numpy
Résultats lors de l'accélération des calculs numériques avec Python et Numba
Remarques sur l'accélération du code Python avec Numba
Chargement / affichage et accélération de gif avec python [OpenCV]
Test numpy Python Basic 8
Moyenne mobile avec numpy
Premiers pas avec Numpy
Apprenez avec Chemo Informatics NumPy
Concaténation de matrices avec Numpy
Code de bourdonnement avec numpy
Effectuer une analyse de régression avec NumPy
Étendre NumPy avec Rust
Calcul numérique avec Python