import numpy as np
import numpy.linalg as lin
import matplotlib.pylab as plt
import plotly.offline as py
import plotly.graph_objects as go
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
np.set_printoptions(precision=3, linewidth=150, suppress=True)
plt.style.use(['seaborn-whitegrid','data/cours.mplstyle'])
Calculons ${\bf x}^T \, A \, {\bf x} $ avec Numpy¶
Le but est de tracer la courbe $J_A({\bf x})$ avec ${\bf x} \in ℝ^2$ pour que les pauvres petits humains que nous sommes puissent voir quelque chose.
$$ J_A({\bf x}) = {\bf x}^T \, A \, {\bf x} $$
Plus tard on prendra ${\bf x} \in ℝ^n$.
Pour une valeur x ce calcule est simple : x.T @ A @ x mais on veut faire ce calcul pour un ensemble de
x.
On commence par construire un maillage c.a.d. un ensemble de points ${\bf x}$ pour lesquels on calculera
$J_A({\bf x})$. Pour cela on utilise np.meshgrid. Regardons ce que cela donne.
x = np.linspace(-1,1,3)
y = np.linspace(-1,2,4)
mesh = np.meshgrid(x,y) # donne les x puis les y du maillage
M = np.array(mesh)
print("shape of M:", M.shape)
print(M)
M = M.transpose([1,2,0]) # je préfère avoir un tableau 4x3 de points
print("\nshape of M:", M.shape)
M
shape of M: (2, 4, 3) [[[-1. 0. 1.] [-1. 0. 1.] [-1. 0. 1.] [-1. 0. 1.]] [[-1. -1. -1.] [ 0. 0. 0.] [ 1. 1. 1.] [ 2. 2. 2.]]] shape of M: (4, 3, 2)
array([[[-1., -1.],
[ 0., -1.],
[ 1., -1.]],
[[-1., 0.],
[ 0., 0.],
[ 1., 0.]],
[[-1., 1.],
[ 0., 1.],
[ 1., 1.]],
[[-1., 2.],
[ 0., 2.],
[ 1., 2.]]])
Je veux appliquer ${\bf x}^T \, A \, {\bf x}$ à chaque point de mon maillage. Je commence avec la matrice identité pour que je puisse facilement vérifier les calculs.
En premier je calcule ${\bf x}^T \, A$ pour tous les points du maillage.
Cas test avec A = 2 Id¶
Je regarde sur un tout petit maillage comment écrire les calculs. Je prends pour A deux fois la matrice identité pour pouvoir vérifier les calculs de tête.
A = 2 * np.diag([1, 1])
MA = np.einsum("ijk, ka -> ija", M, A) # voici un cas où la notation d'Einstein simplifie bien la vie
MA # k est en double et pas dans le résultat, donc on somme sur lui
array([[[-2., -2.],
[ 0., -2.],
[ 2., -2.]],
[[-2., 0.],
[ 0., 0.],
[ 2., 0.]],
[[-2., 2.],
[ 0., 2.],
[ 2., 2.]],
[[-2., 4.],
[ 0., 4.],
[ 2., 4.]]])
Je retrouve 2 M ce qui est cohérent.
Dans un autre cas on pourrait vérifier certains points :
# vérification sur un point particulier
M[0,1] @ A
array([ 0., -2.])
Notez que j'ai écrit ${\bf x} \, A$ et non ${\bf x}^T \, A$ mais lorsqu'on a un vecteur Numpy s'adapte et privilégie le produit matrice vecteur qui donne un vecteur. Ainsi m[0,1] @ A == m[0,1].T @ A car, comme
dit la doc, "Transposing a 1-D array returns an unchanged view of the original array".
Si on veux avoir une différence entre un vecteur vertical et un vecteur horizontal il faut utiliser
des array 2D de taille 1xn ou nx1.
np.einsum("ijk, ijk -> ij", MA, M) # comme k n'est pas dans le résultat, c'est sur lui qu'on fait la somme
array([[ 4., 2., 4.],
[ 2., 0., 2.],
[ 4., 2., 4.],
[10., 8., 10.]])
Comme A est 2 fois l'identité on doit trouver pour tout point 2 fois sa norme au carré. C'est bon.
Un vrai cas¶
On prend un maillage de 100x100 et une matrice A aléatoire :
np.random.seed(2)
A = np.random.randint(-10,10, size=(2,2)) / 100
print(A)
x = np.linspace(-20,20,100)
y = np.linspace(-20,20,100)
mx, my = np.meshgrid(x,y)
M = np.array(np.meshgrid(x,y)).transpose([1,2,0])
MA = np.einsum("ijk, ka -> ija", M, A)
mz = np.einsum("ijk, ijk -> ij", MA, M)
[[-0.02 0.05] [ 0.03 -0.02]]
trace = go.Surface(x=mx, y=my, z=mz)
fig = go.Figure(data=[trace])
py.iplot(fig, show_link=False)
On a un superbe point selle.
Bien sûr avec une autre graine aléatoire nous aurions autre chose. Imaginez ce que vous pourriez avoir si toutes les valeurs de A sont positives.
def J(x):
return x @ A @ x
mz2 = np.array([J(x) for line in M for x in line]).reshape([M.shape[0], M.shape[1]])
np.max(np.square(mz2 - mz)) # pour vérifier qu'on obtient la même chose
0.0
# comparons les temps calculs :
%timeit np.einsum("ijk, ijk -> ij", np.einsum("ijk, ka -> ija", M, A), M)
143 µs ± 1.99 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit np.array([J(x) for line in M for x in line]).reshape([M.shape[0], M.shape[1]])
19.7 ms ± 494 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
On est plus de 100 fois plus rapide avec la sommation d'Einstein.
Notons qu'avec apply_along_axis c'est pire :
%timeit np.array([np.apply_along_axis(J, 1, line) for line in M]).reshape([M.shape[0], M.shape[1]])
34.6 ms ± 989 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Utiliser np.tensordot¶
Un tensor est une matrice en N dimension donc voici une fonction qui devrait pouvoir servir puisque notre ensemble de point est en 3 dimensions.
%%timeit # pour calculer le temps d'execution de toute la cellule %%
MA = np.tensordot(M, A, axes=(2,1)) # on somme sur l'axe 2 de M (les points) et l'axe 1 de A (les colonnes)
np.einsum("ijk, ijk -> ij", MA, M) # je n'arrive pas à faire cette opération autrement
100 µs ± 1.15 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
On a un gain de 30 % par rapport à einsum.
Pour avoir de bonne mesures il faudrait utiliser des cas plus gros que 100x100 et lancer les tests sur une machine qui ne fait que cela, donc en mode single user sous Linux.
Conclusion¶
Grosso modo on a :
| Méthode | Temps |
|---|---|
einsum |
140 µs |
tensordot |
100 µs |
boucle for |
20 000 µs |