Présentation de la notation d'Einstein

Cette notation permet d'exprimer de facon concise un nombre extraordinaire d'opérations. Cette section n'est pas nécessaire pour la suite mais c'est un exercice intellectuel intéressant.

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 :

$$ A_{i,i} \quad \textrm{veut dire} \quad \sum_{i=1}^N A_{i,i} \; \textrm{(la trace de la matrice)}$$

Si on regarde la multiplication matricielle $A\, B$ on a pour tout indice (i,j) du résultat :

$$ 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} $$

Le nom complet de la notation d'Einstein étant la convention de sommation d'Einstein, la fonction de Numpy est einsum. Voici ce que cela donne pour nos 2 premiers exemples :

import numpy as np

A = np.arange(9).reshape(3,3)

print("Trace 'ii' : ", np.einsum('ii', A), '\n')   # 0 + 4 + 8 = 12

print("Multiplication matricielle A A 'ij,jk->ik' :")
print(np.einsum('ij,jk->ik', A, A))  # notez que j'ai nommé différement les indices

Trace 'ii' :  12 

Multiplication matricielle A A 'ij,jk->ik' :
[[ 15  18  21]
 [ 42  54  66]
 [ 69  90 111]]

A.dot(A)  # on vérifie

array([[ 15,  18,  21],
       [ 42,  54,  66],
       [ 69,  90, 111]])

On note que les arguments de einsum sont :

• la règle de sommation dans une chaine de caractères avec une virgule pour séparer chaque composant
• les composants sur lesquels s'applique la règle

On peut aller un peu plus loin avec Numpy. Voici l'ensemble des règles qu'utilise einsum :

Règles de base et complémentaires

1. un indice répété implique une sommation sur cet indice sauf si cet indice est cité dans le résultat
   (cf exemple diagonale de A ci-dessous pour l'exception)
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
   (cf exemple du produit matriciel)
3. un lettre omise dans le résultat (après ->) implique que l'on sommme les éléments sur cet indice
   (cf exemple de la somme des éléments d'un vecteur ci-dessous)
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
   (cf exemple de la transposée ci-dessous)

Voici une liste d'opérations prises sur le blog du Dr Goulu :

Signature de einsum	équivalent NumPy	Description
('i->', v)	sum(v)	somme des valeurs du vecteur v
('i,i->i', u, v)	u \* v	multiplication des vecteurs u et v élément par élément
('i,i', u, v)	inner(u, v)	produit scalaire de u et v
('i,j', u, v)	outer(u, v)	produit dyadique de u et v
('ij', A)	A	renvoie la matrice A
('ji', A)	A.T	transposée de A
('ii->i', A)	diag(A)	diagonale de A
('ii', A)	trace(A)	somme de la diagonale de A
('ij->', A)	sum(A)	somme des valeurs de A
('ij->j', A)	sum(A, axis=0)	somme des colonnes de A
('ij->i', A)	sum(A, axis=1)	somme des lignes de A
('ij,ij->ij', A, B)	A \* B	multiplication des matrices A et B élément par élément
('ij,ji->ij', A, B)	A \* B.T	multiplication des matrices A et B transposée élément par élément
('ij,jk', A, B)	dot(A, B)	produit scalaire de A et B
('ij,jk->ij', A, B)	inner(A, B)	produit intérieur de A et B
('ij,jk->ijk', A, B)	A[:, None] \* B	chaque ligne de A multipliée par B
('ij,kl->ijkl', A, B)	A[:, :, None, None] \* B	chaque valeur de A multipliée par B

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.

Mise en pratique

On va comparer les performances d'einsum et des fonctions Numpy correspondante. Pour cela calculez

• le cube de chaque élement d'un vecteur
• d'une matrice carrée, $A^3$,

avec einsum et sans. Comparez la vitesse d'exécution dans tous les cas avec 10000 éléments.

# a vous de jouer

Solution :

u = np.random.random(10000)

%timeit u*u*u
%timeit np.einsum('i,i,i->i', u,u,u)

11.5 µs ± 112 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
15.2 µs ± 131 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

A = u.reshape(100,100)

%timeit A.dot(A).dot(A)
%timeit np.einsum('ij,jk,kl->il', A, A, A)

138 µs ± 9.48 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
134 ms ± 1.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

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.

(*) 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.