CAMA : ma31 x.T A x sur un maillage en Numpy

Cours du 17/05

Calculons

$x^{T}Ax$ avec Numpy

Nous voulons tracer la courbe, avec

$x\in {\mathbb{R}}^2$ :

$$J_A(x)=x^{T}Ax$$
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

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])

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^{T}Ax$, on commence par calculer

$x^{T}A$ pour tous points du maillage. On utilise la matrice identité pour vérifier nos calculs.

Cas avec

$A=2\ast Id$

A = 2 * np.diag([1, 1]) # Construit un array diagonal
MA = np.einsum("ijk, ka -> ija", M, A) # Notation d'Einstein

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 retrouve 2*M.

On peut vérifier sur un certain point:

M[0,1] @ A

array([ 0., -2.])

Notez qu'on a écrit

$xA$ et non

$x^{T}A$. 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

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 la matrice identité x2, on retrouve pour tout point sa norme au carré

Optimisons :np.tensordot

Un tensor est une matrice en N dimensions.

Comparons les temps de calcul.

%timeit np.einsum("ijk, ijk -> ij", np.einsum("ijk, ka -> ija", M, A), M)

135 µs ± 2.56 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%%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)   

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