Python pour le calcul scientifique/Manipulation de matrices

Les matrices sont un élément primordial du calcul scientifique sur ordinateur pour deux raisons :

L'algèbre linéaire est au cœur de nombreux calculs.
Les matrices sont l'élément de base du calcul vectorisé qui permet un gain de temps appréciable.

Pour pouvoir expliter les matrices, il faut charger le module NumPy ; nous utilisons également Matplotlib pour les graphiques. Ainsi, les programmes contiennent tous au début :

#!/usr/bin/python3

import numpy as np
import matplotlib.pyplot as plt

Rappels et complément sur le tranchage

Soit une liste L composée de n éléments.

Les éléments sont numérotés de 0 à n – 1.
Le premier élément s'obtient par L[0].
Le i-ème élément s'obtient par L[i - 1].
L'avant-dernier élément s'obtient par L[-2].
Le dernier élément s'obtient par L[-1].
La sous-liste composée des éléments contigus i à j s'obtient par L[i - 1:j].

Ainsi, L[0:1] va extraire le premier élément, L[1:-1] va extraire tous les éléments sauf le premier et le dernier.

NumPy fournit des fonctions permettant de manipuler les matrices :

np.append(A, B) : fusionne les vecteurs A et B ;
s'il s'agit de matrices ou de tenseurs, la fonction les « aplatit », les transforme en vecteur ;
si l'on veut intégrer B dans A, on utilise A = np.append(A, B) ;
np.append(A, B, axis = i) : fusionne les tenseurs selon l'indice i (0 pour le premier indice, 1 pour le deuxième…) ;
np.insert(A, i, m) : insère le vecteur m dans le vecteur A (ou la matrice A aplatie) à l'emplacement i ;
np.insert(A, i, M, axis = j) : insère le tenseur M dans le tenseur A à l'emplacement i de l'indice j ;
np.delete(A, I) : efface les éléments définis par le tranchage I du vecteur A (ou de la matrice A aplatie) ;
np.delete(A, I, axis = j) : efface les éléments définis par le tranchage I selon l'indice j du tenseur A.

Définir un tenseur

Un tenseur est similaire à une liste mais il est défini par la fonction np.array(). La définition et l'extraction de composante utilise la méthode du découpage en tranches (slicing).

Exemples

a = np.array([1, 3, 5, 7]) # vecteur
bc = np.array([[1], [2], [3], [4]]) # matrice 4 × 1 (matrice colonne)
bl = np.array([[1, 2, 3, 4]]) # matrice 1 × 4 (matrice ligne)
c = np.array([[1, 2, 3],
              [2, 3, 4],
              [3, 4, 5]]) # matrice 3 × 3
d = np.array([[[1,  2, 3],
               [2,  3, 4]],
              [[10, 9, 8],
               [ 7, 6, 5]]]) # tenseur d'ordre 3, de dimension 3 × 2 × 2

Notez que dans NumPy, un vecteur n'est pas la même chose qu'une matrice ligne ou colonne. Un vecteur de dimension n est un tenseur d'ordre 1 et de dimension n ; une matrice ligne ou colonne est un tenseur d'ordre 2 et de dimension 1 × n ou n × 1.

a = np.array([1, 2, 3])
b = np.array([[1, 2, 3]])
c = np.array([[1], [2], [3]])
print(a.size, b.size, c.size) # 3 3 3
print(a.ndim, b.ndim, c.ndim) # 1 2 2
print(a.shape, b.shape, c.shape) # (3,) (1, 3) (3, 1)

La fonction np.arange() est similaire à la fonction range() pour les liste ; elle génère un vecteur de réels. La fonction np.linspace() permet également de créer un vecteur de même type, mais on indique le dernier nombre alors que la règle du découpage en tranches fait que le nombre maximal indiqué à np.arange() est le premier nombre qui ne figure pas dans le vecteur.). La fonction np.zeros() génère une matrice nulle, np.zeros_like() une matrice nulle ayant les dimensions d'une matrice fournie comme modèle. De même, np.ones() et np.ones_like() crée des matrices, dont toutes les composantes sont à 1. La fonction np.eye() crée une matrice unité.

Exemples

e = np.arange(0, 2, 0.1) # vecteur [0, 0.1, 0.2…, 1.8, 1.9]
f = np.linspace(0, 2, 5) # 5 nombres entre 0 et 2 soit le vecteur [0, 0.5, 1, 1.5, 2]
g = np.zeros(3) # vecteur nul de dimension 3
h = np.zeros((3, 3)) # matrice nulle 3 × 3
k = np.ones_like(a) # matrice de 1 de même dimension que a 
u = np.eye(3) # matrice unité 3 × 3

Le paramètre dtype permet de forcer le type. Par exemple

a = np.array([1, 2, 3], dtype="complex")
k = np.ones_like(a, dtype="int")

La commande np.linspace() peut créer des matrices colonne : on donne la première et la dernière ligne ; par exemple : :

np.linspace([1, -1], [10, -10], 4)
# [[  1.  -1.]
#  [  4.  -4.]
#  [  7.  -7.]
#  [ 10. -10.]]

Les commandes np.geomspace() et np.geomspace() fonctionnent comme np.linspace(), mais avec une progression logarithmique. La commande np.geomspace(a, b, n) crée un vecteur (ou une matrice colonne) allant de a à b alors que np.logspace(a, b, n) crée un vecteur (ou une matrice colonne) allant de 10^a à 10^b.

La méthode .reshape() remet en forme une matrice. Par exemple, pour transformer un vecteur de dimension 9 en une matrice 3 × 3 :

a = np.arange(1, 10)
b = a.reshape(3, 3)
# ou bien directement
c = np.arange(1, 10).reshape(3, 3)

Avec la méthode .reshape(), on peut utiliser la valeur –1 pour une des dimensions ; sa valeur est alors automatiquement calculée en fonction du nombre d'éléments et de l'autre dimension. On peut aussi utiliser une dimension vide ; cela crée alors un vecteur. Par exemple, pour une matrice M quelconque :

M.reshape(-1, 1) # crée une matrice colonne
M.reshape(1, -1) # crée une matrice ligne
M.reshape(-1,) # crée un vecteur

La méthode .fill() remplit la matrice avec un scalaire :

b.fill(5) # remplace les valeurs de b par la valeur 5

Extraction matricielle

Si l'on veut extraire les colonnes 1 et 3 pour toutes les lignes d'une matrice M, on utilise :

M[:, [0, 2]]

Assemblage de matrices

Si l'on veut regrouper deux matrices, on utilise la commande np.concatenate(). On utilise le paramètre axis pour indiquer si l'on empile les matrices l'une au-dessus de l'autre (axis=0) ou bien l'une derrière l'autre (axis=1). Par exemple :

A = np.matrix([[1 ,2], [3, 4]])
B = np.matrix([[5 ,6], [7, 8]])
np.concatenate((A, B), axis = 0)
# [[1 2]
#  [3 4]
#  [5 6]
#  [7 8]]
np.concatenate((A, B), axis = 1)
# [[1 2 5 6]
#  [3 4 7 8]]

Si l'on veut transformer deux vecteurs en une matrice de deux colonnes, chaque vecteur occupant une colonne :

A = np.matrix([1 ,2])
B = np.matrix([3 ,4])
np.concatenate((A.reshape(-1, 1), B.reshape(-1, 1)), axis = 1)
# [[1 3]
#  [2 4]]

Opérations matricielles

Les quatre opérations classiques +, -, * et / ne fonctionnent qu'entre tenseurs de mêmes dimensions et sont des opérations élément par élément (elementwise operations) :

(A + B)[i, j, k] vaut A[i, j, k] + B[i, j, k] ;
(A - B)[i, j, k] vaut A[i, j, k] - B[i, j, k] ;
(A * B)[i, j, k] vaut A[i, j, k] * B[i, j, k] ;
(A / B)[i, j, k] vaut A[i, j, k] / B[i, j, k].

De même, les fonctions np.add(), np.subtract(), np.multiply() et np.divide() effectuent des opérations élément par élément sur des tenseurs de mêmes dimensions :

np.add(A, B)[i, j, k] vaut A[i, j, k] + B[i, j, k] ;
np.subtract(A, B)[i, j, k] vaut A[i, j, k] - B[i, j, k] ;
np.multiply(A, B)[i, j, k] vaut A[i, j, k] * B[i, j, k] ;
np.divide(A, B)[i, j, k] vaut A[i, j, k] / B[i, j, k].

La multiplication matricielle, au sens de l'algèbre linéaire, se fait avec les fonctions np.dot() ou np.matmul(), ou bien avec l'opérateur @.

a = np.matrix([[1, 2], [3, 4]])
b = np.matrix([[5, 6], [7, 8]])
print(a @ b) # np.matmul(a, b)
# [[19 22]
# [43 50]]

Pour un poduit scalaire :

A = np.matrix([1, 2, 3])
B = np.matrix([4, 5, 6])
print(np.dot(A, B.T))
# ［［</nowiki>32］］
print(np.dot(A, B.T)[0, 0])
# 32

Notez que :

pour le produit par un scalaire, les fonctions np.multiply() et l'opérateur * sont plus performants ;
la fonction np.dot() est plus performante pour le produit scalaire de deux vecteurs réels ;
lorsque l'on a des vecteurs complexes, la fonction np.vdot() fait le produit par le conjugué du premier membre (np.vdot(a, b) == np.dot(a.conj(), b)) ;
la fonction np.matmul() et l'opérateur @ (A @ B) sont plus performants pour un produit matriciel.

L'opérateur @= fait un produit matriciel en modifiant la variable elle-même (à l'image de *= pour les nombres).

Calcul vectorisé

Les fonctions de NumPy traitent en général les matrices en entier. Ainsi, il n'est pas nécessaire de créer une boucle pour faire défiler les indices un par un. Il en résulte un code clair et compact et surtout un plus grande rapidité d'exécution. Par exemple :

x = np.linspace(0, 2*np.pi, 50) # 50 points entre 0 et 2π
y = np.sin(x)
plt.plot(x, y)

La variable x est un vecteur de 50 valeurs et il est traité en une seule passe par la fonction sinus np.sin().

Outre le tranchage (slicing), on peut utiliser deux autres méthodes pour extraire certaines valeurs d'une matrice :

utiliser un vecteur ou une matrice d'indices, Python extrait alors les valeurs correspondant aux indices ;
utiliser un vecteur ou une matrice de booléens de même dimension que a matrice ; Python extrait alors les valeurs correspondant aux True, la matrice booléenne est un « masque » pour la matrice d'origine. Par exemple :

a = np.arange(0, 10)
b = np.array([1, 3, 5, 7, 9])
c = np.array([True, True, False, False, True, False, True, False, False, True])
print(a[b], "\n", a[c])

Si l'on veut inverser tous les éléments d'une matrice de bolléens, il faut utiliser la fonction np.logical_not()

Exercice

Écrire un programme Python mettant en œuvre le crible d'Érathostène pour trouver les nombres premiers inférieurs à une valeur donnée.

Solution

#!/usr/bin/env python
# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt

# ***************
# ***************
# ** Fonctions **
# ***************
# ***************

def eratosthene(limite):
    # Détermine la liste des nombres premiers entre 1 et N
    #     par le crible d'Ératosthène
    # Entrées : limite : nombre entier, N
    # Sorties : liste : vecteur d'entiers,liste des nombres premiers
    
    indices = (np.ones(limite) == 1) # vecteur de booléens tous à True
    # à la fin, indice(i-1) est True si i est premier, False sinon
    indices[0]=False # 1 n'est pas premier
    imax = int(limite)
    i = 2 # initialisation
    repete = (i <= imax)
    while repete:
        if indices[i-1]:
            jmax = int(limite/i)
            j = np.arange(1, jmax)+1
            indices[i*j-1]=False # élimination des multiples de i
            test = (i*jmax == imax)
        i = i + 1
        repete = i*i < limite # condition d'arrêt
    liste0 = np.arange(limite)
    liste = liste0[indices]+1
    return liste

# *************************
# *************************
# ** Programme principal **
# *************************
# *************************

print("***** Recherche de nombres premiers par le crible d'Ératosthène *****\n")

nmax = eval(input("Entrez la valeur maximale : "))
resultat = eratosthene(nmax)

print("\n", resultat.shape[0], "nombres premiers entre 1 et", nmax, ":\n")
print(resultat)

plt.plot(resultat, np.zeros_like(resultat), "|")

L'extraction par un vecteur d'indice intervient dans l'instruction :

indices[i*j-1]=False

qui élimine en une seule passe tous les multiples de i. L'extraction par un vecteur booléen intervient dans l'instruction :

liste = liste0[indices]+1

qui permet d'extraire toutes la valeurs conservées en une seule passe.

Attributs

La classe ndarray, qui définit les matrices, possède un certain nombre d'attributs :

.shape : dimensions de la matrice ;
.ndim : ordre du tenseur ;
.size : nombre d'éléments ;
.dtype : type des éléments.

a = np.linspace(1, 9, 9)
print("a", a, "\n",
    " ; shape :", a.shape,
    " ; dim : ", a.ndim,
    " ; size : ", a.size,
    " ; dtype : ", a.dtype, "\n")

.real, .imag : parties réelle et imaginaire de la matrice ;
.flat : liste des éléments de la matrice ; les éléments sont réorganisés en une liste ;
.T : transposée.

a = np.arange(0, 9).reshape(3, 3)
print(a)
# [[0 1 2]
#  [3 4 5]
#  [6 7 8]] 
print(a.T)
# [[0 3 6]
#  [1 4 7]
#  [2 5 8]] 
print(a.flat[:])
# array([0, 1, 2, 3, 4, 5, 6, 7, 8])

Ressources: Section « Attribute », « numpy.ndarray », sur Numpy and Scipy Documentation (consulté le 16 mars 2019)

Fonctions et méthodes de base

La classe ndarray possède un certain nombre de méthodes :

.min() et max() : valeurs respectivement minimale et maximale ;
.ptp() : amplitude « max – min » (peak to peak) ;
.argmin() et argmax() : indice où se trouvent les valeurs respectivement minimale et maximale ;
.sum(), prod() : somme et produit de tous les éléments de la matrice ;
.cumsum(), cumprod() : somme et produit cumulés.

a = np.linspace(1, 9, 9)
print("min : ", a.min(),
    "; max : ", a.max(), "\n",
    "sum : ", a.sum(),
    "; cumsum : ", a.cumsum(), "\n",
    "prod : ", a.prod(),
    "; cumprod : ", a.cumprod(), "\n")

Méthodes statistiques :

.mean() : moyenne ;
.std() : écart type (standard deviation).

Extraction de données :

.diagonal() : vecteur contenant les éléments de la diagonale ;
.flatten() : vecteur contenant les éléments réorganisés en liste ; par rapport à l'attribut .flat, on peut choisir le sens de linéarisation (par lignes, .flatten(C), ou par colonnes, .flatten(F)) mais cela crée une copie, on ne peut pas par exemple s'en servir pour modifier la matrice ;
.tofile() : crée un fichier texte contenant les valeurs de la matrice ; par exemple, pour une matrice a et pour séparer les valeurs par un point-virgule :

a.tofile("matriceA.txt", sep=" ; ")

.astype(type) : copie la matrice en convertissant le type de données :

a = a.astype(float) # pour avoir une matrice de réels en virgule flottante
a = a.astype(str) # pour avoir une matrice de chaînes de caractères

Tri :

np.sort(M, i) : crée une copie et trie la matrice selon l'axe i (0 pour le premier indice, 1 pour le deuxième… la valeur par défaut est –1 pour le dernier indice, la valeur none aplatit la matrice) ;
M.sort(i) : cette méthode trie la matrice en elle-même ;
np.argsort(M, i) : crée une matrice contenant le nouvel indice, selon l'axe i, de chaque élément.

Pour trier les lignes selon la deuxième colonne d'une matrice, on peut faire comme suit :

M = np.array([[1, 5, 7, 5], [2, 9, 4, 3], [9, 1, 7, 5]])
ind = np.argsort(M[:, 1]) # indices de tri selon la 2e colonne
N = M[ind, :] # réarrangement des lignes

Pour faire un tri lexicographique, on utilise la fonction np.lexsort() en indiquant les différentes colonnes :

M = np.array([[1, 5, 7, 5], [2, 9, 4, 3], [9, 1, 7, 5], [2, 9, 8, 5]])
ind = np.lexsort((M[:, 3], M[:, 2], M[:, 1], M[:, 0])) # indices de tri selon la 1re, puis la 2e, puis la 3e, puis la 4e colonne
N = M[ind, :] # réarrangement des lignes

Algèbre linéaire :

a.dot(b) : produit matriciel a⋅b ; on peut aussi écrire a@b ;
.trace() : trace de la matrice (somme des éléménts diagonaux) ;
.transpose() : transpose la matrice, résultat similaire à l'attribut .T ;
np.cross() : produit vectoriel dans ℝ³.

Matrices de booléens :

.all() : applique un « et » logique à toutes les valeurs de la matrice ;
.any() : applique un « ou » logique à toutes les valeurs de la matrice.

Autre méthodes :

.conj() : conjugué des valeurs complexes ;
.nonzero() : n-uplet contenant les indices des valeurs non-nulles ;
.round(n) : arrondit les valeurs à la n-ième décimale.

Ressources: Section « Method », « numpy.ndarray », sur Numpy and Scipy Documentation (consulté le 16 mars 2019)

Propagation

Le terme « propagation » (broadcasting) désigne la manière dont Python complète les matrice lorsque des dimensions manquent.

Supposons que l'on veuille additionner deux matrices M₁ et M₂ de dimensions m₁ × n₁ et m₂ × n₂ différentes. Alors :

le résultat a pour dimension max(m₁, m₂) × max(n₁, n₂) ;
si une des dimensions vaut 1, alors les valeurs de l'autre dimension sont dupliquées ;
sinon, les dimensions manquantes pour chaque matrice sont complétées par des 1.

Par exemple :

\mathrm {A} ={\begin{pmatrix}1&2\\3&4\end{pmatrix}}

5+\mathrm {A} ={\begin{pmatrix}5&5\\5&5\end{pmatrix}}+{\begin{pmatrix}1&2\\3&4\end{pmatrix}}

La matrice (5) est de dimension 1 × 1, la valeur « 5 » est donc répétée dans les deux dimensions

A = np.array([[1, 2], [3, 4]])
print(5 + A)
# [[6 7]
#  [7 8]]

(5,4)+\mathrm {A} ={\begin{pmatrix}5&4\\5&4\end{pmatrix}}+{\begin{pmatrix}1&2\\3&4\end{pmatrix}}

print(np.array([[5, 4]]) + A)
# [[6 6]
#  [6 8]]

{\begin{pmatrix}5\\4\end{pmatrix}}+\mathrm {A} ={\begin{pmatrix}5&5\\4&4\end{pmatrix}}+{\begin{pmatrix}1&2\\3&4\end{pmatrix}}

print(np.array([[5], [4]]) + A)
# [[6 7]
#  [7 8]]

Fonctions « universelles »

Les fonctions universelles (ufunc) sont les fonctions s'appliquant aux matrices, des fonctions vectorisées.

Ressource: « Available ufuncs », sur SciPy documentation (consulté le 21 mars 2019)

Algèbre linéaire

Pour plus de détails voir : Python pour le calcul scientifique/Algèbre linéaire.

Notes et références

Graphiques < ↑ > Polynômes