Bibliothèque Numpy

1Le tableau NumPy : creation, types, indexation

Les listes Python sont souples mais lentes pour le calcul numerique : chaque element est un objet independant, et chaque operation passe par l'interpreteur. NumPy fournit un type ndarray homogene, stocke de facon contigue, sur lequel les operations sont executees en C : le gain est d'un facteur 10 a 100.

import numpy as np

a = np.array([1, 2, 3])              # tableau 1D depuis une liste
b = np.array([[1, 2], [3, 4]])       # 2D
a.shape      # (3,)          forme
b.shape      # (2, 2)
b.ndim       # 2             nombre de dimensions
b.size       # 4             nombre total d'elements
b.dtype      # dtype('int64')  type des elements

Creation de tableaux

np.zeros(5)                    # [0. 0. 0. 0. 0.]        flottants
np.zeros((2, 3))               # matrice 2x3 de zeros
np.ones((2, 2))
np.full((2, 2), 7)             # rempli de 7
np.eye(3)                      # matrice identite 3x3
np.arange(0, 10, 2)            # [0 2 4 6 8]   comme range, mais tableau
np.linspace(0, 1, 5)           # [0. 0.25 0.5 0.75 1.]  n points, bornes INCLUSES
np.random.random((2, 2))       # uniformes dans [0,1[
np.random.randn(3)             # loi normale centree reduite
np.array([1, 2], dtype=float)  # type impose

np.arange(0, 1, 0.1)      # 10 valeurs, borne 1 EXCLUE, pas fixe
np.linspace(0, 1, 11)     # 11 valeurs, borne 1 INCLUSE, nombre de points fixe

Le type des elements

np.array([1, 2, 3]).dtype           # int64
np.array([1.0, 2, 3]).dtype         # float64  : un seul flottant suffit
np.array([1, "a"]).dtype            # <U21     : tout devient chaine !
np.zeros(3, dtype=int)              # [0 0 0]  entiers

t = np.array([1, 2, 3])
t[0] = 3.7                          # 3 : TRONQUE, car le dtype est entier
print(t)                            # [3 2 3]

t = np.array([2**62, 2**62], dtype=np.int64)
print(t.sum())      # -4611686018427387904 : DEBORDEMENT, aucun avertissement

n = 2**62 + 2**62   # 9223372036854775808 : correct, les int Python sont exacts

Indexation et tranches

A = np.arange(12).reshape(3, 4)
# [[ 0  1  2  3]
#  [ 4  5  6  7]
#  [ 8  9 10 11]]

A[1, 2]        # 6         ligne 1, colonne 2 (et non A[1][2], moins efficace)
A[0]           # [0 1 2 3]     premiere ligne
A[:, 1]        # [1 5 9]       deuxieme COLONNE : impossible aussi simplement en liste
A[0:2, 1:3]    # [[1 2], [5 6]]  sous-bloc
A[::2, ::2]    # une ligne sur deux, une colonne sur deux
A[-1, -1]      # 11
A.T            # transposee (sans copie !)

A = np.arange(5)
B = A[1:4]           # VUE sur la memoire de A
B[0] = 99
print(A)             # [ 0 99  2  3  4]  : A est modifiee !

C = A[1:4].copy()    # copie explicite
C[0] = -1
print(A)             # inchangee

Indexation booleenne et avancee

t = np.array([3, -1, 4, -5, 9])

t > 0                     # [True False True False True]  : tableau de booleens
t[t > 0]                  # [3 4 9]        selection
t[t < 0] = 0              # remplacement conditionnel
np.where(t > 2, t, 0)     # [3 0 4 0 9]    "si ... alors ... sinon"
np.nonzero(t)             # indices des elements non nuls

t[[0, 2, 4]]              # [3 4 9]  indexation par liste d'indices
t[np.array([4, 0])]       # [9 3]    dans l'ordre demande

Ces formes remplacent des boucles entieres : t[t > 0] tient en une expression la ou une liste demanderait une comprehension, et s'execute cent fois plus vite.

2Vectorisation, agregation et algebre lineaire

Vectorisation : le coeur de NumPy

t = np.array([1, 2, 3, 4])

t + 10          # [11 12 13 14]   element par element
t * 2           # [2 4 6 8]
t ** 2          # [1 4 9 16]
t + t           # [2 4 6 8]
t * t           # [1 4 9 16]      produit TERME A TERME, pas matriciel
np.sqrt(t)      # [1. 1.414 1.732 2.]
np.exp(t)       # exponentielle de chaque element
np.sin(t)

import numpy as np, time

n = 10**7
liste = list(range(n))
tableau = np.arange(n)

t0 = time.perf_counter()
r1 = [x*x for x in liste]
t1 = time.perf_counter()
r2 = tableau ** 2
t2 = time.perf_counter()

print(f"liste  : {t1-t0:.3f} s")     # ~ 1.2 s
print(f"numpy  : {t2-t1:.3f} s")     # ~ 0.02 s
print(f"gain   : {(t1-t0)/(t2-t1):.0f}x")

# AVANT : boucle Python
resultat = []
for x in donnees:
    if x > 0:
        resultat.append(sqrt(x))

# APRES : vectorise
donnees = np.array(donnees)
resultat = np.sqrt(donnees[donnees > 0])

Fonctions d'agregation

A = np.array([[1, 2, 3], [4, 5, 6]])

A.sum()              # 21       tout le tableau
A.sum(axis=0)        # [5 7 9]        somme des COLONNES (on ecrase les lignes)
A.sum(axis=1)        # [6 15]         somme des LIGNES
A.mean(), A.std()    # moyenne, ecart-type
A.min(), A.max()     # extremums
A.argmin(), A.argmax()   # INDICE de l'extremum (dans le tableau aplati)
A.cumsum()           # [1 3 6 10 15 21]   sommes cumulees
np.median(A)         # 3.5
(A > 3).sum()        # 3   nombre d'elements > 3 (True vaut 1)
(A > 3).any(), (A > 3).all()      # (True, False)

Algebre lineaire

A = np.array([[1, 2], [3, 4]])
B = np.array([[0, 1], [1, 0]])
v = np.array([1, 1])

A * B            # produit TERME A TERME : [[0 2], [3 0]]
A @ B            # produit MATRICIEL     : [[2 1], [4 3]]
A.dot(B)         # equivalent
A @ v            # produit matrice-vecteur : [3 7]

np.linalg.det(A)          # -2.0000000000000004
np.linalg.inv(A)          # matrice inverse
np.linalg.solve(A, v)     # resout A x = v  (JAMAIS via inv : moins precis et plus lent)
np.linalg.eig(A)          # valeurs et vecteurs propres
np.linalg.matrix_rank(A)  # rang
np.linalg.norm(v)         # norme euclidienne
A.T                       # transposee
np.trace(A)               # trace

x = np.linalg.inv(A) @ b       # MAUVAIS : plus lent et moins precis
x = np.linalg.solve(A, b)      # BON : factorisation LU avec pivot

Diffusion (broadcasting)

A = np.ones((3, 4))
A + 5                       # (3,4) + scalaire   -> (3,4)
A + np.array([1, 2, 3, 4])  # (3,4) + (4,)       -> ajoute a chaque LIGNE
A + np.array([[1], [2], [3]])   # (3,4) + (3,1)  -> ajoute a chaque COLONNE

x = np.linspace(0, 1, 3)          # (3,)
y = np.linspace(0, 2, 4)          # (4,)
D = x[:, None] - y[None, :]       # (3,1) - (1,4) -> (3,4) : toutes les differences

# application : matrice des distances entre n points, sans aucune boucle
P = np.random.random((5, 2))
diff = P[:, None, :] - P[None, :, :]        # (5,5,2)
D = np.sqrt((diff**2).sum(axis=2))          # (5,5) : toutes les distances
assert np.allclose(D, D.T) and np.allclose(np.diag(D), 0)

Le broadcasting evite des boucles imbriquees et l'allocation de tableaux intermediaires : c'est l'outil qui rend NumPy reellement expressif.

3Visualisation, comparaison et applications

Liste ou tableau NumPy ?

Application — resolution d'un systeme lineaire

import numpy as np

A = np.array([[2, 1, -1],
              [-3, -1, 2],
              [-2, 1, 2]], dtype=float)
b = np.array([8, -11, -3], dtype=float)

x = np.linalg.solve(A, b)
print(x)                                # [ 2.  3. -1.]
print(np.allclose(A @ x, b))            # True : verification systematique
print(f"conditionnement : {np.linalg.cond(A):.2f}")

H = np.array([[1/(i+j+1) for j in range(10)] for i in range(10)])
print(f"{np.linalg.cond(H):.2e}")     # 1.60e+13 : tres mal conditionnee