UKOnline

Calcul matriciel

Maintenant que l'on maitrise les vecteurs et les calculs qu'il est possible de faire avec ces derniers, on va s'intéresser aux opérations sur les matrices. La première section nous a permis de découvrir comment créer des matrices et on va maintenant voir comment les utiliser dans des calculs. Pour rappel, voici un exemple de code qui crée trois matrices :

u = np.array([1, 2, 3]) A = np.matrix('1 2 3;4 5 6') B = np.identity(2) C = np.diag(u)

La première matrice est créée directement en spécifiant ses valeurs, la deuxième est une matrice identité d'ordre $2$ et enfin la troisième est une matrice diagonale dont les éléments sont ceux du vecteur $\overrightarrow{u}$ :

$$A = \left( \begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \end{array} \right), \quad B = \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right), \quad C = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 3 \end{array} \right).$$

Les matrices $B$ et $C$ sont des matrices carrées, c'est-à-dire qu'elles ont le même nombre de lignes que de colonnes, et la matrice $A$ est une matrice rectangulaire. L'exécution du code produit exactement ces trois matrices, sans surprise :

[[1 2 3]
 [4 5 6]]
[[1. 0.]
 [0. 1.]]
[[1 0 0]
 [0 2 0]
 [0 0 3]]

Propriété

Une fois un objet matrix créé, on peut obtenir la valeur de propriétés spécifiques aux matrices, en plus de toutes les propriétés des ndarray, évidemment. Certaines propriétés sont directement accessibles par un attribut et d'autres sont obtenues en appelant une fonction. Enfin, pour certaines autres propriétés, comme tester si une matrice est carrée ou non, il faudra effectuer un calcul.

Dimension

La première information intéressante que l'on peut obtenir à propos d'une matrice est sa dimension, également appelée taille. On l'obtient simplement avec l'attribut shape. On l'utilise également pour tester si une matrice est carrée, en comparant les deux valeurs de sa forme, comme dans l'exemple ci-dessous :

dim = A.shape print(dim) print(dim[0] == dim[1])

La matrice $A$ étant une matrice de dimension $2 \times 3$, elle n'est pas carrée comme le confirme le résultat de l'exécution :

(2, 3)
False

Diagonale principale

Une deuxième propriété d'une matrice est sa diagonale principale, c'est-à-dire les éléments $(a)_{ij}$ tels que $i = j$ ou, dis autrement, ceux situés à l'intersection d'un même numéro de ligne et de colonne. On obtient cette diagonale principale avec la fonction diag, que l'on a déjà vue à la section 3.1.2. L'exemple suivant extrait deux diagonales d'une matrice :

print(np.diag(A)) print(np.diag(A, k=1))

La première instruction extrait la diagonale principale et la seconde extrait celle se trouvant une position plus haut, grâce au paramètre optionnel $k$ déjà expliqué précédemment :

[1 5]
[2 6]

Trace

La trace d'une matrice carrée est la somme des éléments de sa diagonale principale. On peut l'obtenir avec la fonction trace. L'exemple suivant calcule la trace d'une matrice identité d'ordre $2$ :

print(np.trace(B))

Une matrice identité n'ayant que des $1$ sur sa diagonale principale et la matrice $B$ étant d'ordre $2$, sa trace vaut $2$. Comme les éléments de la matrice sont des données de type float, le résultat de l'exécution est aussi un nombre flottant :

2.0

Rang

Le rang d'une matrice est une autre information intéressante qui s'obtient avec la fonction matrix_rank du module numpy.linalg. Il représente notamment le nombre de lignes linéairement indépendantes d'une matrice quelconque. L'exemple suivant calcule le rang de deux matrices :

D = np.matrix('1 2 3;2 4 6;3 6 9') print(matrix_rank(A)) print(matrix_rank(D))

Le rang de la matrice $A$ vaut $2$, car ses deux lignes ne sont pas liées entre elles. Par contre, celui de la matrice $D$ vaut $1$, même si elle a trois lignes, car peu importe la ligne que l'on sélectionne, on peut obtenir les deux autres à partir de celle choisie. Le résultat de l'exécution donne bien les valeurs attendues :

2
1

Si on sélectionne la première ligne $L_1 = (1, 2, 3)$, on se rend directement compte que l'on a $L_2 = 2L_1$ et $L_3 = 3L_1$. Une seule ligne suffit donc à caractériser la matrice d'où son rang de $1$.

Déterminant

Une autre propriété, uniquement définie sur les matrices carrées cette fois-ci, peut également être obtenue à partir d'une fonction du module numpy.linalg. Il s'agit du déterminant, puissant outil intervenant notamment dans la résolution de systèmes d'équations linéaires. On l'obtient avec la fonction det, illustrée dans l'exemple suivant :

D = np.matrix('1 2 3;2 4 6;3 6 9') print(det(C)) print(det(D))

La matrice $C$ étant diagonale, son déterminant est simplement égal à la somme des éléments de sa diagonale, c'est-à-dire sa trace. Le déterminant de la matrice $D$ est, quant à lui, nul, car ses lignes sont liées entre elles, comme on l'a vu précédemment lorsque l'on a calculé son rang qui est strictement inférieur à son ordre.

Le résultat de l'exécution confirme les résultats attendus :

6.0
0.0

Opération de base

Tout comme pour les vecteurs, on peut facilement faire la somme et la différence de deux matrices de mêmes dimensions, et la multiplication par un scalaire, à l'aide des opérateurs arithmétiques habituels vus avec les ndarray. Voici un exemple de calcul matriciel qui combine ces différentes opérations :

print(2 * A[:,:2] + B)

On commence par multiplier par deux la sous-matrice de $A$ qui ne prend que ses deux premières colonnes, de sorte à avoir une matrice carrée d'ordre $2$. On ajoute ensuite la matrice $B$ au résultat du calcul précédent. Le calcul réalisé est donc le suivant :

$$2 \left( \begin{array}{cc} 1 & 2 \\ 4 & 5 \end{array} \right) + \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right) = \left( \begin{array}{cc} 3 & 4 \\ 8 & 11 \end{array} \right).$$

Comme la matrice $B$, créée avec la fonction identity, est de type float par défaut, le résultat final sera également de ce type, comme le confirme le résultat de l'exécution, qui donne d'ailleurs bien la réponse attendue :

[[ 3.  4.]
 [ 8. 11.]]

Transposée

Calculer la transposée d'une matrice revient à construire une nouvelle matrice dont les lignes sont les colonnes de la matrice originale, et vice-versa. Cette opération est tellement importante et utilisée en algèbre linéaire que l'on obtient la transposée d'une matrice directement avec l'attribut T. On peut, par exemple, écrire :

print(A.T)

Le résultat de l'exécution est donc une matrice avec 3 lignes et 2 colonnes, dont les lignes sont les colonnes de la matrice $A$, et inversement :

[[1 4]
 [2 5]
 [3 6]]

Cet attribut construit une vue de la matrice originale, simplement avec une autre forme. Si l'on veut modifier les éléments de la transposée d'une matrice, il faut donc au préalable en faire une copie.

Il faut également faire attention si on doit transposer un vecteur, c'est-à-dire un tableau à une dimension. Alors que la transposée d'une matrice-ligne devrait être une matrice-colonne, on va avoir un tout autre résultat, si on ne fait pas attention. Voyons ce qui se passe avec l'exemple suivant :

u = np.array([1, 2, 3]) print(u.shape, u.T.shape)

Le vecteur u a donc comme forme $(3,)$ et on pourrait s'attendre à ce que sa transposée soit de forme $(1, 3)$, mais ce n'est pas le cas comme le montre le résultat de l'exécution :

(3,) (3,)

Pour éviter ce problème, lorsque l'on travaille avec des vecteurs à transposer, il faut avant tout les convertir en matrices-ligne. On fait cela avec la fonction mat qui va interpréter son paramètre pour le convertir en une matrice, à deux dimensions donc :

U = np.mat(u) print(U.shape, U.T.shape)

Le résultat obtenu est cette fois-ci celui attendu :

(1, 3) (3, 1)

En plus de la transposée classique obtenue avec l'attribut T, on peut obtenir la matrice adjointe, c'est-à-dire la transposée de la matrice conjuguée, avec l'attribut H. Cette matrice adjointe s'avère utile pour certaines opérations dans le monde des nombres complexes. L'exemple suivant calcule cette matrice adjointe pour une matrice complexe :

E = np.matrix('3+1j 5;2-2j 1j') print(E.H)

On commence donc par calculer la matrice conjuguée (changement de signe des parties imaginaires), et puis on transpose le résultat :

$$\left( \begin{array}{cc} 3 + i & 5 \\ 2 - 2i & i \end{array} \right)^* = \left( \begin{array}{cc} 3 - i & 5 \\ 2 + 2i & -i \end{array} \right)^t = \left( \begin{array}{cc} 3 - i & 2 + 2i \\ 5 & -i \end{array} \right)$$

Le résultat de l'exécution produit bel et bien la même réponse, avec j au lieu de i pour représenter le nombre imaginaire :

[[3.-1.j 2.+2.j]
 [5.-0.j 0.-1.j]]

Produit

Tout comme pour les vecteurs, différentes possibilités existent pour le produit de matrices. Ces différentes opérations sont évidemment supportées par NumPy.

Produit matriciel

Étant donné deux matrices compatibles, c'est-à-dire que le nombre de colonnes de la première est égal au nombre de lignes de la seconde, on peut effectuer un produit matriciel entre ces deux matrices.

Pour rappel, le résultat du produit matriciel est une matrice dont chaque élément en position $(i, j)$ correspond au produit scalaire de la $i$^e ligne de la première matrice avec la $j$^e colonne de la seconde. Soient les deux matrices $A$ et $B$ suivantes, respectivement de dimensions $m \times n$ et $n \times p$, et donc compatibles :

$$A = \left( \begin{array}{ccc} a_{11} & \dots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{m1} & \dots & a_{mn} \end{array} \right), \qquad B = \left( \begin{array}{ccc} b_{11} & \dots & b_{1p} \\ \vdots & \ddots & \vdots \\ b_{n1} & \dots & b_{np} \end{array} \right),$$

leur produit $AB$ est une matrice $C$ de dimensions $m \times p$ dont la valeur des éléments est donnée par :

$$c_{ij} = a_{i\bullet} \cdot b_{\bullet j} = \left( \begin{array}{ccc} a_{i1} & \dots & a_{in} \end{array} \right) \left( \begin{array}{c} b_{1j} \\ \vdots \\ b_{nj} \end{array} \right) = \sum_{k = 0}^n a_{ik} b_{kj}.$$

En NumPy, le produit matriciel s'obtient avec l'opérateur *. On peut donc, par exemple, écrire les instructions suivantes pour multiplier entre elles deux matrices compatibles :

print(D * C[:,:2]) print(U * D)

La première instruction multiplie une matrice $3 \times 3$ par une $3 \times 2$ pour produire une matrice $3 \times 2$ comme résultat, tandis que la seconde multiplication se fait entre une matrice $1 \times 3$ et une $3 \times 3$ pour produire une matrice $1 \times 3$ comme résultat :

$$\left( \begin{array}{ccc} 1 & 2 & 3 \\ 2 & 4 & 6 \\ 3 & 6 & 9 \end{array} \right) \left( \begin{array}{cc} 1 & 0 \\ 0 & 2 \\ 0 & 0 \end{array} \right) = \left( \begin{array}{cc} 1 & 4 \\ 2 & 8 \\ 3 & 12 \end{array} \right)$$ $$\left( \begin{array}{ccc} 1 & 2 & 3 \end{array} \right) \left( \begin{array}{ccc} 1 & 2 & 3 \\ 2 & 4 & 6 \\ 3 & 6 & 9 \end{array} \right) = \left( \begin{array}{ccc} 14 & 28 & 42 \end{array} \right)$$

Le résultat de l'exécution confirme que NumPy a bel et bien calculé des produits matriciels avec l'opérateur * :

[[ 1  4]
 [ 2  8]
 [ 3 12]]
[[14 28 42]]

Produit par élément

Un autre type de produit que l'on peut vouloir faire entre deux matrices de mêmes dimensions est le produit par élément, comme celui que l'on a déjà vu plus haut pour les vecteurs.

Ce produit est en fait celui qui se fait entre deux tableaux multidimensionnels avec l'opérateur *, en profitant éventuellement de l'opérateur de broadcasting. Néanmoins, lorsque l'on utilise cet opérateur entre deux matrix, on n'obtient pas un produit par élément, même si cette dernière est une sous-classe de ndarray.

Pour comprendre les différences entre ce qui se passe avec l'opérateur * entre deux ndarray ou entre deux matrix, examinons l'exemple suivant :

data = np.array([[1, 2], [3, 4]]) E = np.mat(data) print(data * data) print(E * E)

La variable data contient un ndarray et l'opération data * data est donc un produit par élément. Par contre, la variable E contient une matrice, obtenue avec la fonction mat par conversion du tableau data et, du coup, l'opération E * E est un produit matriciel :

[[ 1  4]
 [ 9 16]]
[[ 7 10]
 [15 22]]

Pour être sûrs de ce que l'on calcule, on peut utiliser des fonctions de NumPy. La fonction matmul permet de calculer un produit matriciel entre deux ndarray, tandis que la fonction multiply fait un produit par élément entre deux matrices. Pour compléter l'exemple précédent, et calculer les deux produits manquants, on peut donc écrire :

print(np.matmul(data, data)) print(np.multiply(E, E))

Cette fois-ci, on a fait un produit matriciel avec data et un produit par élément avec E, comme le confirme le résultat de l'exécution :

[[ 7 10]
 [15 22]]
[[ 1  4]
 [ 9 16]]

Comme détaillé plus loin, en section 3.4.4, il est recommandé de travailler avec des tableaux de type ndarray et de minimiser l'utilisation des matrix. Pour éviter de rendre le code illisible, avec des appels à la fonction matmul, NumPy propose l'opérateur @ comme raccourci de cette fonction. On aurait donc pu écrire l'instruction suivante :

print(data @ data)

Exponentiation

Enfin, une dernière opération qui intervient dans de nombreux calculs en algèbre linéaire est la puissance matricielle, c'est-à-dire la multiplication d'une matrice par elle-même, un certain nombre de fois. On peut réaliser cette opération avec l'opérateur **, le même qu'avec Python. Voici, par exemple, comment élever une matrice au carré :

print(E ** 2)

La variable E étant un objet matrix, l'élever au carré revient bien à la multiplier par elle-même suivant le produit matriciel, ce qui est correct et produit le résultat suivant :

[[ 7 10]
 [15 22]]

Puisque les dimensions doivent être compatibles pour pouvoir réaliser une multiplication matricielle, on en déduit que la puissance matricielle est une opération limitée aux matrices carrées. Pour calculer une telle puissance, on peut également passer par la fonction matrix_power du module numpy.linalg. Cette dernière prend en paramètres une matrice et la puissance à laquelle il faut l'élever. L'exemple précédent aurait donc aussi pu s'écrire comme suit :

print(matrix_power(E, 2))

Si E avait été un ndarray, l'opérateur ** aurait fait plusieurs produits par élément successifs, et donc pas une puissance matricielle. Par contre, utiliser la fonction matrix_power garanti que l'opération réalisée est bien une puissance matricielle, même avec des ndarray. En cas de doutes, mieux vaut donc utiliser cette fonction à la place de l'opérateur **.

Inverse

Une autre opération très importante en algèbre linéaire, notamment pour résoudre des systèmes d'équations linéaires, consiste à calculer l'inverse d'une matrice carrée. Si la matrice est inversible, alors on obtient son inverse avec l'attribut I. Repartons de la matrice $E$ de dimensions $2 \times 2$, vue un peu plus haut :

$$E = \left( \begin{array}{cc} 1 & 2 \\ 3 & 4 \end{array} \right).$$

Pour rappel, l'inverse d'une matrice $E$, notée $E^{-1}$, est la matrice par laquelle il faut multiplier $E$ pour obtenir la matrice identité comme résultat, autrement dit, c'est la matrice $E^{-1}$ telle que $EE^{-1} = I$. Pour notre exemple, on peut écrire l'égalité suivante :

$$EE^{-1} = \left( \begin{array}{cc} 1 & 2 \\ 3 & 4 \end{array} \right) \left( \begin{array}{cc} -2 & 1 \\ 3/2 & -1/2 \end{array} \right) = \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right) = I.$$

Vérifions maintenant cette propriété en utilisant l'attribut I pour trouver l'inverse de la matrice E, puis en calculant le produit matriciel de E avec son inverse obtenu avec E.I :

print(E.I) print(E * E.I)

Le résultat obtenu après exécution est un peu différent de celui que l'on aurait pu espérer obtenir, c'est-à-dire np.identity(2) :

[[-2.   1. ]
 [ 1.5 -0.5]]
[[1.0000000e+00 0.0000000e+00]
 [8.8817842e-16 1.0000000e+00]]

On a presqu'une matrice identité puisque l'élément situé à la première colonne de la deuxième ligne ne vaut pas exactement zéro, mais il est par contre très petit ($8,\!88 \cdot 10^{16}$). Le problème vient du calcul en nombres flottants qui n'est, par nature, pas précis. On pourrait essayer de comparer le résultat obtenu avec la matrice identité d'ordre $2$ :

print(E * E.I == np.identity(2))

Le résultat obtenu est assez surprenant, deux des quatre valeurs ne sont pas égales à celles de la matrice identité ! On voit que l'élément à la deuxième colonne de la deuxième ligne n'est pas exactement égal à $1$ :

[[ True  True]
 [False False]]

Pour s'en sortir lorsque l'on fait des calculs en nombres flottants, il ne faut jamais comparer deux valeurs avec l'opérateur d'égalité. Une solution pour comparer deux matrices de float consiste à utiliser la fonction allclose qui permet de faire une égalité approximative :

print(np.allclose(E * E.I, np.identity(2)))

La fonction allclose teste si les éléments aux mêmes positions de deux matrices sont proches, c'est-à-dire qu'ils seraient égaux si l'on n'avait pas les imprécisions dues au calcul en nombres flottants. Cette fois-ci, la propriété que l'on voulait vérifier est bien satisfaite, E * E.I donne bien une matrice identité :

True

La matrice inverse peut se calculer de deux autres manières. La première consiste à utiliser la fonction d'exponentiation matrix_power pour élever la matrice à la puissance $-1$. On peut aussi utiliser la fonction inv du module numpy.linalg. Les deux instructions suivantes sont donc équivalentes et calculent l'inverse de $E$ :

print(matrix_power(E, -1)) print(inv(E))

L'avantage de passer par la fonction matrix_power est que si l'on désire élever à une certaine puissance l'inverse d'une matrice, on peut directement combiner les deux opérations. Par exemple, si on veut élever au cube l'inverse de E, il suffit d'appeler matrix_power(E, -3).

Concernant la seconde possibilité, l'intérêt de passer par numpy.linalg.inv est que cette fonction permet en fait d'inverser plusieurs matrices en même temps, en les empilant les unes au-dessus des autres. Cette possibilité est détaillée plus loin, à la section 3.4.4.

Autre opération

Le module scipy.linalg définit plusieurs fonctions utiles pour la manipulation de matrices, comme on a déjà pu s'en rendre compte. Ce module contient également des fonctions qui permettent de réaliser des opérations par élément sur des matrices.

On peut notamment calculer l'exponentielle (expm, expm2 ou expm3), le logarithme (logm), les fonctions trigonométriques (cosm, sinm et tanm), les fonctions trigonométriques hyperboliques (coshm, sinhm et tanhn), la racine carrée (sqrtm) et enfin la fonction de signe (signm).

• Précédent
• Suivant