Introduction à Python pour la Data Science¶

TP 1 — NumPy¶

Objectifs¶

À l’issue de ce TP, vous devrez être capable de :

  • créer et inspecter des arrays multidimensionnels ;
  • raisonner sur les attributs shape, ndim, size et dtype ;
  • utiliser l’indexation, le slicing et les boolean masks ;
  • distinguer une view d’une copie ;
  • vectoriser un calcul et exploiter le broadcasting ;
  • effectuer des agrégations selon un axis ;
  • mettre en œuvre des opérations d’algèbre linéaire ;
  • contrôler la reproductibilité d’une expérience numérique.

1. Introduction¶

En analyse de données, on manipule fréquemment des séries de valeurs numériques sur lesquelles on applique des opérations statistiques, matricielles ou algébriques. Les listes Python peuvent être utilisées pour des calculs simples, mais elles ne sont pas conçues pour le calcul scientifique intensif.

NumPy est la librairie de référence pour le calcul numérique en Python. Elle fournit notamment :

  • des tableaux multidimensionnels homogènes, appelés ndarray ;
  • des opérations vectorisées ;
  • des fonctions d’algèbre linéaire et de statistique ;
  • une gestion compacte et efficace de la mémoire ;
  • une forte interopérabilité avec pandas, SciPy, scikit-learn et les librairies de deep learning.

Documentation officielle : NumPy documentation

Dans ce TP, le terme array désigne un objet numpy.ndarray.

In [ ]:
import numpy as np

print("NumPy version:", np.__version__)

2. Pourquoi utiliser NumPy ?¶

La différence essentielle entre une liste Python et un array NumPy ne se limite pas à la syntaxe.

Un array NumPy :

  • contient en général des valeurs d’un même dtype ;
  • stocke ses données dans une zone mémoire compacte ;
  • délègue de nombreux calculs à du code compilé ;
  • permet d’exprimer un calcul sur tout un tableau sans boucle Python explicite.

Cette approche est appelée vectorisation.

In [ ]:
def mult_list(n: int):
    a = range(n)
    b = range(n)

    result = []
    for i in range(n):
        result.append(a[i] * b[i])

    return result


def mult_np(n: int):
    a = np.arange(n)
    b = np.arange(n)
    return a * b
In [ ]:
n = 100_000

resultat_liste = mult_list(n)
resultat_numpy = mult_np(n)

assert resultat_liste == resultat_numpy.tolist()
print("Les deux implémentations produisent le même résultat.")
In [ ]:
%%timeit -n 10
mult_list(n)
In [ ]:
%%timeit -n 10
mult_np(n)

Analyse¶

Le gain observé dépend de la machine et de la taille des données. Il ne faut donc pas retenir un facteur d’accélération universel.

La vectorisation améliore généralement :

  • les performances ;
  • la concision ;
  • la lisibilité ;
  • la possibilité de raisonner directement sur des vecteurs et des matrices.

Elle peut cependant créer des tableaux intermédiaires volumineux. La performance doit donc être analysée à la fois en temps de calcul et en mémoire.

3. Les arrays¶

L’objet central de NumPy est le ndarray. Il peut représenter :

  • un scalaire : ndim == 0 ;
  • un vecteur : ndim == 1 ;
  • une matrice : ndim == 2 ;
  • un tenseur d’ordre supérieur : ndim >= 3.

Les attributs suivants sont fondamentaux :

  • shape : dimensions de l’array ;
  • ndim : nombre d’axes ;
  • size : nombre total d’éléments ;
  • dtype : type numérique des éléments ;
  • itemsize : taille en octets d’un élément ;
  • nbytes : mémoire occupée par les données.
In [ ]:
x = np.array([[1, 2, 3], [4, 5, 6]], dtype=np.float64)

print("shape   :", x.shape)
print("ndim    :", x.ndim)
print("size    :", x.size)
print("dtype   :", x.dtype)
print("itemsize:", x.itemsize, "octets")
print("nbytes  :", x.nbytes, "octets")

3.1 Création d’arrays¶

La fonction np.array construit un array à partir d’une structure Python.

Les fonctions np.arange et np.linspace créent des séquences numériques :

  • np.arange(start, stop, step) utilise un pas ;
  • np.linspace(start, stop, num) impose un nombre de valeurs, bornes incluses par défaut.

Pour les nombres flottants, np.linspace est souvent préférable lorsque l’on souhaite contrôler précisément le nombre de points.

In [ ]:
liste = [1, 2, 3, 4, 5]
a = np.array(liste)

print(a)
print(np.arange(0, 50, 10))
print(np.linspace(0, 1, 5))

On peut également initialiser des tableaux particuliers avec :

  • np.zeros ;
  • np.ones ;
  • np.full ;
  • np.eye pour une matrice identité.
In [ ]:
print(np.zeros(5, dtype=float))
print(np.ones((2, 3), dtype=int))
print(np.full((2, 2), 3.14))
print(np.eye(3))

3.2 dtype et conversion¶

Le dtype influence :

  • la précision numérique ;
  • la plage de valeurs représentables ;
  • l’occupation mémoire ;
  • le comportement de certaines opérations.

Une conversion explicite peut être effectuée avec astype. Elle crée généralement un nouvel array.

In [ ]:
a_int = np.array([1, 2, 3], dtype=np.int32)
a_float = a_int.astype(np.float64)

print(a_int, a_int.dtype)
print(a_float, a_float.dtype)

Attention aux dépassements de capacité
Les types entiers ont une plage de valeurs finie. Un calcul avec un type trop petit peut produire un overflow sans lever d’exception.

In [ ]:
x8 = np.array([120], dtype=np.int8)
print("Valeur initiale :", x8)
print("Après addition  :", x8 + np.int8(20))

3.3 Génération aléatoire et reproductibilité¶

L’API moderne repose sur un générateur créé avec np.random.default_rng.

Fixer une seed permet de reproduire les mêmes tirages, ce qui est indispensable pour :

  • débuguer ;
  • comparer deux méthodes ;
  • documenter une expérience ;
  • vérifier un résultat.
In [ ]:
rng = np.random.default_rng(seed=42)

uniforme = rng.uniform(low=0.0, high=1.0, size=(3, 3))
normale = rng.normal(loc=0.0, scale=1.0, size=(3, 3))

print("Uniforme :\n", uniforme)
print("\nNormale :\n", normale)

4. Indexation, slicing et boolean masks¶

Comme pour les listes Python, les indices commencent à 0. Pour un tableau 2D, on utilise généralement la syntaxe array[ligne, colonne].

In [ ]:
a = np.array([[1, 2, 3], [4, 5, 6]])

print(a[1, 2])
print(a[:, 1])
print(a[0, :])

Le slicing suit la notation :

array[start:stop:step]

La borne stop est exclue.

In [ ]:
arr = np.arange(10)

print(arr[2:6])
print(arr[:5])
print(arr[5:])
print(arr[::2])
print(arr[::-1])

mat = np.arange(1, 10).reshape(3, 3)
print("\nSous-matrice :\n", mat[:2, 1:])

4.1 boolean masks¶

Une comparaison appliquée à un array produit un tableau de booléens de même shape. Ce masque peut ensuite être utilisé pour filtrer les données.

In [ ]:
x = np.arange(10)
mask = (x > 3) & (x % 2 == 0)

print(mask)
print(x[mask])

Les opérateurs and, or et not ne doivent pas être utilisés directement entre des arrays. On utilise :

  • & pour le ET élément par élément ;
  • | pour le OU ;
  • ~ pour la négation.

Chaque condition doit être placée entre parenthèses.

4.2 view et copie¶

Un slicing simple produit généralement une view, c’est-à-dire un objet qui partage les mêmes données sous-jacentes que l’array d’origine.

Modifier la view peut donc modifier l’array original.

In [ ]:
x = np.arange(6)
view = x[1:4]

view[0] = 999

print("view :", view)
print("x    :", x)

Pour obtenir un tableau indépendant, il faut utiliser explicitement copy.

In [ ]:
x = np.arange(6)
copie = x[1:4].copy()

copie[0] = 999

print("copie:", copie)
print("x    :", x)

5. shape, reshape et ajout d’axes¶

reshape modifie l’organisation logique d’un tableau sans nécessairement recopier les données. Le nombre total d’éléments doit rester inchangé.

In [ ]:
x = np.arange(12)

matrice = x.reshape(3, 4)
print(matrice)

print("\nRetour en 1D :", matrice.reshape(-1))

La valeur -1 demande à NumPy d’inférer automatiquement la dimension correspondante.

Pour ajouter un axe, on peut utiliser :

  • np.newaxis ;
  • np.expand_dims.

Ces opérations sont très utiles pour préparer un broadcasting.

In [ ]:
x = np.array([10, 20, 30])

colonne = x[:, np.newaxis]
ligne = x[np.newaxis, :]

print("colonne shape:", colonne.shape)
print("ligne shape  :", ligne.shape)

6. Opérations vectorisées et broadcasting¶

Les opérateurs +, -, *, / et ** sont appliqués terme à terme.

Le produit matriciel utilise l’opérateur @.

In [ ]:
a = np.array([[2, 2, 1], [2, 5, 1]])
b = np.array([[3, 4, 1], [4, 7, 3]])

print("Somme terme à terme :\n", a + b)
print("\nProduit terme à terme :\n", a * b)

6.1 Règles du broadcasting¶

NumPy compare les dimensions depuis le dernier axe vers le premier. Deux dimensions sont compatibles si :

  1. elles sont égales ;
  2. l’une des deux vaut 1 ;
  3. l’une des dimensions est absente.

Le broadcasting n’effectue pas nécessairement une duplication physique des données. Il fournit une vue logique compatible avec l’opération.

In [ ]:
mat = np.arange(12).reshape(3, 4)
offset_colonnes = np.array([10, 20, 30, 40])
offset_lignes = np.array([100, 200, 300])[:, np.newaxis]

print("Ajout par colonne :\n", mat + offset_colonnes)
print("\nAjout par ligne :\n", mat + offset_lignes)

Question de compréhension¶

Prédire la shape du résultat avant d’exécuter chaque opération :

A = np.zeros((8, 1, 6, 1))
B = np.zeros((7, 1, 5))
A + B

Justifier la réponse axe par axe.

In [ ]:
A = np.zeros((8, 1, 6, 1))
B = np.zeros((7, 1, 5))

# Décommenter après avoir formulé votre prédiction.
# print((A + B).shape)

7. Agrégations et paramètre axis¶

Sans paramètre axis, une agrégation porte sur tous les éléments.

Avec axis=k, l’axe k est supprimé du résultat, sauf si keepdims=True.

In [ ]:
a = np.array([[1, 2, 3], [4, 5, 6]])

print("Somme totale       :", np.sum(a))
print("Somme selon axis=0 :", np.sum(a, axis=0))
print("Somme selon axis=1 :", np.sum(a, axis=1))
print("Avec keepdims      :\n", np.sum(a, axis=1, keepdims=True))

Pour une matrice de shape (n, p) :

  • axis=0 agrège les lignes et retourne une valeur par colonne ;
  • axis=1 agrège les colonnes et retourne une valeur par ligne.

L’expression « travailler sur les lignes » peut être ambiguë. Il est préférable de raisonner en termes d’axe supprimé et de shape du résultat.

7.1 Standardisation manuelle¶

La standardisation d’une matrice de données consiste à centrer puis réduire chaque colonne :

[ z_{ij} = \frac{x_{ij} - \mu_j}{\sigma_j} ]

où (\mu_j) et (\sigma_j) sont calculés sur les observations de la colonne (j).

In [ ]:
rng = np.random.default_rng(42)
X = rng.normal(size=(100, 4))

moyennes = X.mean(axis=0, keepdims=True)
ecarts_types = X.std(axis=0, keepdims=True)

X_standardise = (X - moyennes) / ecarts_types

print("Moyennes :", X_standardise.mean(axis=0))
print("Écarts-types :", X_standardise.std(axis=0))

8. Algèbre linéaire¶

NumPy fournit le sous-module np.linalg.

Quelques opérations fréquentes :

  • produit matriciel : A @ B ;
  • norme : np.linalg.norm ;
  • résolution de système : np.linalg.solve ;
  • valeurs et vecteurs propres : np.linalg.eigh ou np.linalg.eig ;
  • décomposition en valeurs singulières : np.linalg.svd.
In [ ]:
A = np.array([[3.0, 1.0], [1.0, 2.0]])
b = np.array([9.0, 8.0])

solution = np.linalg.solve(A, b)

print("Solution :", solution)
print("Vérification A @ x :", A @ solution)

Pour résoudre (Ax=b), il faut préférer np.linalg.solve(A, b) à np.linalg.inv(A) @ b.
Le calcul explicite de l’inverse est généralement plus coûteux et moins stable numériquement.

In [ ]:
v = np.array([3.0, 4.0])

print("Norme euclidienne :", np.linalg.norm(v))
print("Produit scalaire  :", np.dot(v, v))

9. Fonctions universelles et stabilité numérique¶

Les fonctions comme np.exp, np.log, np.sqrt, np.sin ou np.maximum sont appelées des ufuncs. Elles sont vectorisées et gèrent le broadcasting.

Les calculs flottants ont une précision finie. Deux valeurs obtenues par calcul ne doivent donc pas toujours être comparées avec ==.

In [ ]:
a = 0.1 + 0.2
b = 0.3

print(a == b)
print(np.isclose(a, b))

Pour comparer deux arrays, on utilise généralement :

np.allclose(A, B)

avec éventuellement des tolérances adaptées au problème.

10. Synthèse¶

Avant d’écrire un calcul NumPy, il est utile de vérifier systématiquement :

  1. la shape des entrées ;
  2. le dtype ;
  3. l’axe selon lequel le calcul doit être appliqué ;
  4. la compatibilité du broadcasting ;
  5. la création éventuelle d’une view ou d’une copie ;
  6. le coût mémoire des tableaux intermédiaires ;
  7. la reproductibilité des tirages aléatoires.

Exercices¶

Exercice 1 — Loi des gaz parfaits¶

On considère la relation :

[ PV = nRT ]

avec :

  • (P) : pression en Pa ;
  • (V) : volume en m³ ;
  • (n) : quantité de matière en mol ;
  • (T) : température absolue en K ;
  • (R = 8{,}314\ \mathrm{J,mol^{-1},K^{-1}}).

Créer une fonction vectorisée gazparfait retournant (n) à partir de P, V et T.

Calculer les quantités de matière associées aux volumes suivants, exprimés en litres :

[15, 89, 56, 78, 152, 66, 48, 77, 2, 96]

On prendra :

  • (P = 101325\ \mathrm{Pa}) ;
  • (T = 20^\circ\mathrm{C}).

Contraintes¶

  • convertir les litres en m³ ;
  • convertir les degrés Celsius en kelvins ;
  • ne pas utiliser de boucle Python ;
  • vérifier la shape et le dtype du résultat ;
  • identifier le volume correspondant à la quantité de matière maximale.
In [ ]:
# Votre code

Exercice 2 — Statistiques descriptives et reproductibilité¶

À l’aide de np.random.default_rng :

  1. générer X, contenant 10 000 réalisations d’une loi uniforme (U(0,1)) ;
  2. générer Y, contenant 10 000 réalisations d’une loi normale de moyenne 0 et d’écart-type 2 ;
  3. comparer les moyennes et variances empiriques aux valeurs théoriques ;
  4. écrire une fonction statdesc retournant, dans un dictionnaire :
    • la moyenne ;
    • la médiane ;
    • l’écart-type ;
    • le minimum ;
    • le maximum ;
    • les quantiles à 25 % et 75 % ;
  5. relancer les tirages avec la même seed, puis avec une autre seed, et commenter.
In [ ]:
# Votre code

Exercice 3 — Manipulation d’un array¶

À partir de la liste suivante :

[17, 38, 10, 25, 72]
  1. créer un array ;
  2. le trier sans modifier l’original ;
  3. ajouter la valeur 12 ;
  4. afficher l’array dans l’ordre inverse ;
  5. retrouver l’indice de la valeur 17 ;
  6. extraire les deuxième et troisième éléments ;
  7. extraire les deux premiers éléments ;
  8. extraire les éléments à partir du troisième ;
  9. afficher le dernier élément avec un indice négatif ;
  10. remplacer toutes les valeurs supérieures à 30 par 30 à l’aide d’un boolean mask.
In [ ]:
# Votre code

Exercice 4 — Distances pairwise et broadcasting¶

Créer une matrice X de shape (10, 2) contenant des points aléatoires dans le plan.

  1. fixer une seed ;
  2. représenter les points avec plt.scatter ;
  3. construire, sans boucle, une matrice D de shape (10, 10) telle que :

$$ D_{ij} = \lVert X_i - X_j \rVert_2 $$

Pour cela, utiliser :

X[:, np.newaxis, :]
X[np.newaxis, :, :]
  1. vérifier que :
    • D est symétrique ;
    • sa diagonale est nulle ;
    • toutes ses valeurs sont positives ou nulles ;
  2. déterminer les deux plus proches voisins de chaque point ;
  3. représenter les arêtes correspondantes ;
  4. comparer le résultat avec une implémentation utilisant deux boucles Python.
In [ ]:
# Votre code
Indice pour la représentation des voisins 

plt.scatter(X[:, 0], X[:, 1], s=100)


K = 2


for i in range(X.shape[0]):
for j in nearest_partition[i, :K]:
plt.plot(*zip(X[j], X[i]), color="black")

Exercice 5 — Standardisation avec NumPy¶

Générer une matrice X de shape (500, 5) dont les colonnes ont des moyennes et des échelles différentes.

  1. calculer les moyennes et écarts-types par colonne ;
  2. standardiser les données sans boucle ;
  3. vérifier numériquement que chaque colonne standardisée possède une moyenne proche de 0 et un écart-type proche de 1 ;
  4. conserver les dimensions lors des agrégations avec keepdims=True ;
  5. écrire une fonction standardize(X) ;
  6. gérer explicitement le cas d’une colonne constante.
In [ ]:
# Votre code

Exercice 6 — Produit matriciel et système linéaire¶

On considère :

[ A = \begin{pmatrix} 4 & 1 & 2 \\ 1 & 5 & 1 \\ 2 & 1 & 3 \end{pmatrix}, \qquad b = \begin{pmatrix} 7 \\ 8 \\ 5 \end{pmatrix} ]

  1. vérifier que les dimensions sont compatibles ;
  2. résoudre (Ax=b) avec np.linalg.solve ;
  3. vérifier la solution avec A @ x ;
  4. calculer la norme du résidu (Ax-b) ;
  5. comparer la solution avec np.linalg.inv(A) @ b ;
  6. expliquer pourquoi solve est préférable.
In [ ]:
# Votre code

Exercice 7 — Vectorisation d’un calcul métier¶

Une entreprise dispose d’une matrice consommation de shape (365, 24) représentant la consommation électrique horaire d’un bâtiment pendant un an.

Générer un jeu de données synthétique positif, puis calculer sans boucle :

  1. la consommation totale de chaque journée ;
  2. la consommation moyenne pour chaque heure de la journée ;
  3. le jour de consommation maximale ;
  4. les heures dépassant de plus de deux écarts-types la consommation horaire moyenne ;
  5. une version normalisée de chaque journée, où la somme des 24 heures vaut 1 ;
  6. les éventuelles divisions par zéro doivent être gérées explicitement.
In [ ]:
# Votre code

Exercice 8 — Question d’approfondissement¶

Soient deux matrices :

A.shape == (n, d)
B.shape == (m, d)

Construire une matrice D de shape (n, m) contenant toutes les distances euclidiennes entre les lignes de A et celles de B.

Proposer deux solutions :

  1. une solution directe par broadcasting ;
  2. une solution utilisant l’identité :

$$ \lVert a-b\rVert^2 = \lVert a\rVert^2 + \lVert b\rVert^2 - 2a^\top b $$

Comparer les besoins mémoire des deux approches lorsque n, m et d deviennent grands.

In [ ]:
# Votre code