diff options
Diffstat (limited to 'PVCM/cama/fr/ma1 np05 Notation Einstein.ipynb')
| -rw-r--r-- | PVCM/cama/fr/ma1 np05 Notation Einstein.ipynb | 323 |
1 files changed, 323 insertions, 0 deletions
diff --git a/PVCM/cama/fr/ma1 np05 Notation Einstein.ipynb b/PVCM/cama/fr/ma1 np05 Notation Einstein.ipynb new file mode 100644 index 0000000..1fc6ed7 --- /dev/null +++ b/PVCM/cama/fr/ma1 np05 Notation Einstein.ipynb @@ -0,0 +1,323 @@ +{ + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2, + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "lang": "fr" + }, + "source": [ + "## Présentation de la notation d'Einstein\n", + "\n", + "Cette notation permet d'exprimer de facon concise un nombre extraordinaire d'opérations. Cette section n'est pas\n", + "nécessaire pour la suite mais c'est un exercice intellectuel intéressant.\n", + "\n", + "L'idée est base est de sommer les termes d'une équation lorsqu'on y trouve 2 fois un même indice qui n'est pas défini par ailleurs. Ainsi :\n", + "\n", + "$$ A_{i,i} \\quad \\textrm{veut dire} \\quad \\sum_{i=1}^N A_{i,i} \\; \\textrm{(la trace de la matrice)}$$\n", + "\n", + "Si on regarde la multiplication matricielle $A\\, B$ on a pour tout indice (i,j) du résultat :\n", + "\n", + "$$ C_{i,j} = A_{i,k} \\, B_{k,j} \\quad \\textrm{c.a.d} \\quad C_{i,j} = \\sum_{k=1}^N A_{i,k} \\, B_{k,j} $$ \n", + "\n", + "Le nom complet de la notation d'Einstein étant la convention de sommation d'Einstein, la fonction de Numpy est [`einsum`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html). Voici ce que cela donne pour nos 2 premiers exemples :" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Trace 'ii' : 12 \n", + "\n", + "Multiplication matricielle A A 'ij,jk->ik' :\n", + "[[ 15 18 21]\n", + " [ 42 54 66]\n", + " [ 69 90 111]]\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "\n", + "A = np.arange(9).reshape(3,3)\n", + "\n", + "print(\"Trace 'ii' : \", np.einsum('ii', A), '\\n') # 0 + 4 + 8 = 12\n", + "\n", + "print(\"Multiplication matricielle A A 'ij,jk->ik' :\")\n", + "print(np.einsum('ij,jk->ik', A, A)) # notez que j'ai nommé différement les indices" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 15, 18, 21],\n", + " [ 42, 54, 66],\n", + " [ 69, 90, 111]])" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A.dot(A) # on vérifie" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "lang": "fr" + }, + "source": [ + "On note que les arguments de `einsum` sont :\n", + "\n", + "* la règle de sommation dans une chaine de caractères avec une virgule pour séparer chaque composant\n", + "* les composants sur lesquels s'applique la règle\n", + "\n", + "On peut aller un peu plus loin avec Numpy. Voici l'ensemble des règles qu'utilise `einsum` :\n", + "\n", + "#### Règles de base et complémentaires\n", + "\n", + "1. un indice répété implique une sommation sur cet indice sauf si cet indice est cité dans le résultat<br/>(cf exemple diagonale de A ci-dessous pour l'exception)\n", + "2. les indices qui se répètent d'un composant à un autre impliquent que les éléments référencés seront multipliés entre eux<br/> (cf exemple du produit matriciel)\n", + "3. un lettre omise dans le résultat (après `->`) implique que l'on sommme les éléments sur cet indice<br/> (cf exemple de la somme des éléments d'un vecteur ci-dessous)\n", + "4. si on ne met pas la flèche, einsum la met pour vous avec à droite tous les indices qui ne sont pas doublés rangés par ordre alphabétique<br/> (cf exemple de la transposée ci-dessous)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "lang": "fr" + }, + "source": [ + "Voici une liste d'opérations prises sur le [blog du Dr Goulu](https://www.drgoulu.com/2016/01/17/einsum) :\n", + "\n", + "<table align=\"center\">\n", + "<tbody>\n", + "<tr>\n", + "<th>Signature de einsum</th>\n", + "<th>équivalent NumPy</th>\n", + "<th>Description</th>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('i->', v)</code></td>\n", + "<td><code>sum(v)</code></td>\n", + "<td>somme des valeurs du vecteur v</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('i,i->i', u, v)</code></td>\n", + "<td><code>u \\* v</code></td>\n", + "<td>multiplication des vecteurs u et v élément par élément</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('i,i', u, v)</code></td>\n", + "<td><code>inner(u, v)</code></td>\n", + "<td>produit scalaire de u et v</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('i,j', u, v)</code></td>\n", + "<td><code>outer(u, v)</code></td>\n", + "<td>produit dyadique de u et v</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('ij', A)</code></td>\n", + "<td><code>A</code></td>\n", + "<td>renvoie la matrice A</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('ji', A)</code></td>\n", + "<td><code>A.T</code></td>\n", + "<td>transposée de A</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('ii->i', A)</code></td>\n", + "<td><code>diag(A)</code></td>\n", + "<td>diagonale de A</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('ii', A)</code></td>\n", + "<td><code>trace(A)</code></td>\n", + "<td>somme de la diagonale de A</td>\n", + "</tr>\n", + "<tr>\n", + "<td class=\"tg-031e\"><code>('ij->', A)</code></td>\n", + "<td class=\"tg-031e\"><code>sum(A)</code></td>\n", + "<td class=\"tg-031e\">somme des valeurs de A</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('ij->j', A)</code></td>\n", + "<td><code>sum(A, axis=0)</code></td>\n", + "<td>somme des colonnes de A</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('ij->i', A)</code></td>\n", + "<td><code>sum(A, axis=1)</code></td>\n", + "<td>somme des lignes de A</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('ij,ij->ij', A, B)</code></td>\n", + "<td><code>A \\* B</code></td>\n", + "<td>multiplication des matrices A et B élément par élément</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('ij,ji->ij', A, B)</code></td>\n", + "<td><code>A \\* B.T</code></td>\n", + "<td>multiplication des matrices A et B transposée élément par élément</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('ij,jk', A, B)</code></td>\n", + "<td><code>dot(A, B)</code></td>\n", + "<td>produit scalaire de A et B</td> </tr>\n", + "<tr>\n", + "<td><code>('ij,jk->ij', A, B)</code></td>\n", + "<td><code>inner(A, B)</code></td>\n", + "<td>produit intérieur de A et B</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('ij,jk->ijk', A, B)</code></td>\n", + "<td><code>A[:, None] \\* B</code></td>\n", + "<td>chaque ligne de A multipliée par B</td>\n", + "</tr>\n", + "<tr>\n", + "<td><code>('ij,kl->ijkl', A, B)</code></td>\n", + "<td><code>A[:, :, None, None] \\* B</code></td>\n", + "<td>chaque valeur de A multipliée par B</td>\n", + "</tr>\n", + "</tbody>\n", + "</table>\n", + "\n", + "Le `None` dans les deux dernières lignes est une facon de redimensionner un tableau. Ainsi `np.arange(6)` est un tableau de dimension 1, `np.arange(6)[:]` est le même tableau de dimension 1, alors que `np.arange(6)[:,None]` est un tableau de dimension 2 à savoir `6 x 1` quand `np.arange(6)[None,:,None]` a 3 dimensions : `1 x 6 x 1`." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "lang": "fr" + }, + "source": [ + "## Mise en pratique\n", + "\n", + "On va comparer les performances d'`einsum` et des fonctions Numpy correspondante. Pour cela calculez\n", + "\n", + "* le cube de chaque élement d'un vecteur\n", + "* d'une matrice carrée, $A^3$,\n", + "\n", + "avec `einsum` et sans. Comparez la vitesse d'exécution dans tous les cas avec 10000 éléments." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# a vous de jouer\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Solution :" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "11.5 µs ± 112 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n", + "15.2 µs ± 131 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" + ] + } + ], + "source": [ + "u = np.random.random(10000)\n", + "\n", + "%timeit u*u*u\n", + "%timeit np.einsum('i,i,i->i', u,u,u)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "138 µs ± 9.48 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n", + "134 ms ± 1.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" + ] + } + ], + "source": [ + "A = u.reshape(100,100)\n", + "\n", + "%timeit A.dot(A).dot(A)\n", + "%timeit np.einsum('ij,jk,kl->il', A, A, A)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "lang": "fr" + }, + "source": [ + "On constate qu'`einsum` est plus lent (\\*). Mais en regardant sur le web on constate que cela n'a pas toujours été le cas et que c'est lié à la version de la bibliothèque Numpy. Conclusion si on veut de la performance, il faut tester son code avant pour choisir la méthode la plus rapide.\n", + "\n", + "(\\*) un peu plus lent pour le calcul vectoriel mais 1000 fois plus lent sur mon portable pour `A.dot(A)` (tous mes processeurs tournent à 100% alors que le calcul par `einsum` n'est fait que sur 1 seul processeur). La version `A.dot(A)` est nettement plus rapide grace à la bibliothèque MKL qu'utilise Numpy sur ma machine." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ] +}
\ No newline at end of file |