Lien de la note Hackmd
Cours du 17/05
Calculons $x^TAx$ avec Numpy
Nous voulons tracer la courbe, avec $x\in\mathbb{R}^2$ : \(J_A(x) = x^TAx\) On prendra $x\in\mathbb{R}^n$ plus tard.
Pour une valeur de $x$ le calcul est x.T @ A @ x, mais on veut faire ce calcul pour un ensemble de $x$
On construit un maillage: un ensemble de point $x$ pour lesquels on calculera $J_A(x)$
On utilise np.meshgrid
1
2
3
4
5
6
x = np.linspace(-1,1,3) # retourne des nombres espaces egalement dans un intervalle specifie
y = np.linspace(-1,2,4)
mesh = np.meshgrid(x,y) # donne les x puis les y du maillage
M = np.array(mesh)
M = M.transpose([1,2,0])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
shape of M: (4, 3, 2)
M = 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.]]])
Pour calculer $x^TAx$, on commence par calculer $x^TA$ pour tous points du maillage. On utilise la matrice identité pour vérifier nos calculs.
Cas avec $A = 2*Id$
1
2
A = 2 * np.diag([1, 1]) # Construit un array diagonal
MA = np.einsum("ijk, ka -> ija", M, A) # Notation d'Einstein
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
MA = 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.]]])
On peut vérifier sur un certain point:
1
M[0,1] @ A
1
array([ 0., -2.])
Notez qu’on a écrit $xA$ et non $x^TA$. Lorsqu’on a un vecteur, Numpy privilégie le produit matrice vecteur qui donne un vecteur. Ainsi: m[0,1] @ A == m[0,1].T @ A
Si on veut une différence entre un vecteur vertical et horizontal, il faut utiliser des arrays 2D de taille 1*n ou n*1
1
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
1
2
3
4
array([[ 4., 2., 4.],
[ 2., 0., 2.],
[ 4., 2., 4.],
[10., 8., 10.]])
Optimisons :np.tensordot
Un tensor est une matrice en N dimensions.
Comparons les temps de calcul.
1
%timeit np.einsum("ijk, ijk -> ij", np.einsum("ijk, ka -> ija", M, A), M)
1
135 µs ± 2.56 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
1
2
3
4
%%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)
1
102 µs ± 1.2 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
On peut avoir un gain de temps 30% et plus ou un peu moins bien que einsum