NumPy — algèbre linéaire et opérations avancées

NumPy — algèbre linéaire et opérations avancées#

Hide code cell source

 
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
import seaborn as sns

sns.set_theme(style="whitegrid", palette="muted", font_scale=1.1)

Manipulation de formes#

La capacité à remodeler des tableaux sans copier leurs données est l’une des forces de NumPy. La plupart des opérations de transformation de forme retournent des vues — des objets qui partagent le même bloc mémoire que le tableau d’origine. Modifier la vue modifie donc aussi le tableau source.

a = np.arange(24)
print("a :", a)
print("shape :", a.shape)

# reshape : changer la forme sans copier les données
b = a.reshape(4, 6)
print("\nreshape(4, 6) :\n", b)

c = a.reshape(2, 3, 4)
print("\nreshape(2, 3, 4) :\n", c)

# L'argument -1 laisse NumPy inférer la dimension manquante
d = a.reshape(6, -1)    # 6 lignes, nombre de colonnes inféré = 4
print("\nreshape(6, -1) :\n", d)

a : [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23]
shape : (24,)

reshape(4, 6) :
 [[ 0  1  2  3  4  5]
 [ 6  7  8  9 10 11]
 [12 13 14 15 16 17]
 [18 19 20 21 22 23]]

reshape(2, 3, 4) :
 [[[ 0  1  2  3]
  [ 4  5  6  7]
  [ 8  9 10 11]]

 [[12 13 14 15]
  [16 17 18 19]
  [20 21 22 23]]]

reshape(6, -1) :
 [[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]
 [16 17 18 19]
 [20 21 22 23]]

# ravel et flatten : mise à plat
a = np.array([[1, 2, 3], [4, 5, 6]])

vue = a.ravel()          # vue (si possible) — modification répercutée sur a
copie = a.flatten()      # toujours une copie

print("ravel()   :", vue)
print("flatten() :", copie)

# Modifier la vue modifie le tableau source
vue[0] = 99
print("a après modification de vue :", a)

ravel()   : [1 2 3 4 5 6]
flatten() : [1 2 3 4 5 6]
a après modification de vue : [[99  2  3]
 [ 4  5  6]]

# transpose et T
M = np.arange(12).reshape(3, 4)
print("M :\n", M)
print("M.T :\n", M.T)

# Pour un tableau 3D, transpose accepte les axes dans l'ordre voulu
T3 = np.ones((2, 3, 4))
print("\nT3.shape :", T3.shape)
print("T3.transpose(2, 0, 1).shape :", T3.transpose(2, 0, 1).shape)

M :
 [[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
M.T :
 [[ 0  4  8]
 [ 1  5  9]
 [ 2  6 10]
 [ 3  7 11]]

T3.shape : (2, 3, 4)
T3.transpose(2, 0, 1).shape : (4, 2, 3)

# np.newaxis : ajouter une dimension de taille 1
v = np.array([1, 2, 3])          # forme (3,)
print("v.shape         :", v.shape)
print("v[:, np.newaxis] :", v[:, np.newaxis].shape)  # (3, 1) — vecteur colonne
print("v[np.newaxis, :] :", v[np.newaxis, :].shape)  # (1, 3) — vecteur ligne

# squeeze : supprimer les dimensions de taille 1
a = np.zeros((1, 3, 1, 4))
print("\nBefore squeeze :", a.shape)
print("After squeeze  :", np.squeeze(a).shape)    # (3, 4)
print("squeeze(axis=2):", np.squeeze(a, axis=2).shape)  # (1, 3, 4)

v.shape         : (3,)
v[:, np.newaxis] : (3, 1)
v[np.newaxis, :] : (1, 3)

Before squeeze : (1, 3, 1, 4)
After squeeze  : (3, 4)
squeeze(axis=2): (1, 3, 4)

Note

np.newaxis est simplement None — les deux notations sont interchangeables. Il est cependant fortement conseillé d’utiliser np.newaxis pour la lisibilité, car l’intention est immédiatement claire. L’ajout de dimensions avec np.newaxis est un idiome fréquent pour rendre deux tableaux compatibles au broadcasting : transformer un vecteur de forme (n,) en (n, 1) ou (1, n) avant une opération.

Empilement et découpage#

Assembler plusieurs tableaux en un seul, ou au contraire diviser un tableau en morceaux, sont des opérations très fréquentes lors du prétraitement de données.

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

# Concaténation selon un axe existant
print("concatenate axis=0 :\n", np.concatenate([a, b], axis=0))  # empiler verticalement
print("concatenate axis=1 :\n", np.concatenate([a, b], axis=1))  # côte à côte

concatenate axis=0 :
 [[1 2]
 [3 4]
 [5 6]
 [7 8]]
concatenate axis=1 :
 [[1 2 5 6]
 [3 4 7 8]]

# vstack, hstack, dstack — raccourcis expressifs
v = np.array([1, 2, 3])
w = np.array([4, 5, 6])

print("vstack :\n", np.vstack([v, w]))   # (2, 3)
print("hstack :", np.hstack([v, w]))     # (6,)
print("column_stack :\n", np.column_stack([v, w]))  # (3, 2)

vstack :
 [[1 2 3]
 [4 5 6]]
hstack : [1 2 3 4 5 6]
column_stack :
 [[1 4]
 [2 5]
 [3 6]]

# np.stack : créer une NOUVELLE dimension
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

print("stack axis=0 :\n", np.stack([a, b], axis=0))  # (2, 3)
print("stack axis=1 :\n", np.stack([a, b], axis=1))  # (3, 2)

stack axis=0 :
 [[1 2 3]
 [4 5 6]]
stack axis=1 :
 [[1 4]
 [2 5]
 [3 6]]

# np.split : découper un tableau
M = np.arange(24).reshape(4, 6)
parties = np.split(M, 2, axis=0)   # 2 parties égales selon les lignes
print("split axis=0 :")
for p in parties:
    print(p)

colonnes = np.hsplit(M, 3)         # 3 parties égales selon les colonnes
print("\nhsplit en 3 :")
for c in colonnes:
    print(c)

split axis=0 :
[[ 0  1  2  3  4  5]
 [ 6  7  8  9 10 11]]
[[12 13 14 15 16 17]
 [18 19 20 21 22 23]]

hsplit en 3 :
[[ 0  1]
 [ 6  7]
 [12 13]
 [18 19]]
[[ 2  3]
 [ 8  9]
 [14 15]
 [20 21]]
[[ 4  5]
 [10 11]
 [16 17]
 [22 23]]

Algèbre linéaire avec np.linalg#

L’algèbre linéaire est au cœur de la quasi-totalité des algorithmes de machine learning et de deep learning. NumPy fournit un sous-module np.linalg complet, qui délègue les calculs aux bibliothèques BLAS et LAPACK — des implémentations hautement optimisées, souvent accélérées matériellement.

A = np.array([[2.0, 1.0, -1.0],
              [-3.0, -1.0, 2.0],
              [-2.0, 1.0, 2.0]])
b = np.array([8.0, -11.0, -3.0])

# Résolution du système linéaire Ax = b
x = np.linalg.solve(A, b)
print("Solution Ax = b :", x)
print("Vérification A @ x :", A @ x)

Solution Ax = b : [ 2.  3. -1.]
Vérification A @ x : [  8. -11.  -3.]

# Produit matriciel : l'opérateur @ (Python 3.5+)
A = np.random.default_rng(0).standard_normal((3, 4))
B = np.random.default_rng(1).standard_normal((4, 5))
C = A @ B
print("A @ B shape :", C.shape)

# Autres opérations linalg
M = np.array([[4.0, 2.0], [1.0, 3.0]])
print("\nDéterminant :", np.linalg.det(M))
print("Inverse :\n", np.linalg.inv(M))
print("Norme (Frobenius) :", np.linalg.norm(M))
print("Rang :", np.linalg.matrix_rank(M))

A @ B shape : (3, 5)

Déterminant : 10.000000000000002
Inverse :
 [[ 0.3 -0.2]
 [-0.1  0.4]]
Norme (Frobenius) : 5.477225575051661
Rang : 2

# Valeurs propres et vecteurs propres
M = np.array([[3.0, 1.0], [1.0, 3.0]])
valeurs, vecteurs = np.linalg.eig(M)
print("Valeurs propres :", valeurs)
print("Vecteurs propres (en colonnes) :\n", vecteurs)

# Vérification : M @ v = lambda * v
for i in range(len(valeurs)):
    v = vecteurs[:, i]
    print(f"M @ v{i} = {M @ v},  λ{i} * v{i} = {valeurs[i] * v}")

Valeurs propres : [4. 2.]
Vecteurs propres (en colonnes) :
 [[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
M @ v0 = [2.82842712 2.82842712],  λ0 * v0 = [2.82842712 2.82842712]
M @ v1 = [-1.41421356  1.41421356],  λ1 * v1 = [-1.41421356  1.41421356]

Décomposition en valeurs singulières (SVD)

La décomposition en valeurs singulières (Singular Value Decomposition, SVD) factorise toute matrice M de taille (m, n) en trois matrices : M = U Σ Vᵀ, où U est orthogonale (m, m), Σ est diagonale (m, n) avec des valeurs non négatives décroissantes appelées valeurs singulières, et Vᵀ est orthogonale (n, n). La SVD est fondamentale en data science : elle sous-tend l’ACP (Analyse en Composantes Principales), la recommandation collaborative, la compression d’images et la pseudo-inverse.

# SVD complète
M = np.random.default_rng(42).standard_normal((4, 6))
U, sigma, Vt = np.linalg.svd(M, full_matrices=True)
print("M shape :", M.shape)
print("U shape :", U.shape)
print("sigma   :", sigma)
print("Vt shape:", Vt.shape)

# Reconstruction exacte
Sigma = np.zeros_like(M)
np.fill_diagonal(Sigma, sigma)
M_rec = U @ Sigma @ Vt
print("\nErreur de reconstruction :", np.max(np.abs(M - M_rec)))

M shape : (4, 6)
U shape : (4, 4)
sigma   : [3.30733024 1.83214124 1.30524693 0.65280266]
Vt shape: (6, 6)

Erreur de reconstruction : 8.881784197001252e-16

Hide code cell source

 
# Illustration : compression d'image par SVD (approximation de rang k)
rng = np.random.default_rng(7)

# Créer une "image" synthétique (gradient + bruit structuré)
x = np.linspace(0, 1, 64)
y = np.linspace(0, 1, 64)
X, Y = np.meshgrid(x, y)
image = np.sin(4 * np.pi * X) * np.cos(3 * np.pi * Y) + 0.3 * rng.standard_normal((64, 64))

U, sigma, Vt = np.linalg.svd(image, full_matrices=False)

rangs = [1, 3, 8, 20, 64]
fig, axes = plt.subplots(1, len(rangs), figsize=(15, 3.5))

for ax, k in zip(axes, rangs):
    approx = U[:, :k] @ np.diag(sigma[:k]) @ Vt[:k, :]
    ax.imshow(approx, cmap='RdBu_r', aspect='auto',
              vmin=image.min(), vmax=image.max())
    ax.set_title(f"Rang k = {k}", fontsize=11, fontweight='bold')
    ax.set_xticks([])
    ax.set_yticks([])

fig.suptitle("Décomposition SVD — reconstruction progressive d'une image synthétique",
             fontsize=13, fontweight='bold', y=1.04)
plt.tight_layout()
plt.show()
_images/1a93df3830e2c3767a796380e3f452ccb22073f86ecb4355ca709f47eb84d74f.png

Statistiques#

NumPy fournit un ensemble complet de fonctions statistiques descriptives. La plupart acceptent un argument axis pour calculer la statistique le long d’une dimension donnée plutôt que sur l’ensemble du tableau.

rng = np.random.default_rng(0)
data = rng.standard_normal((5, 4))   # 5 observations, 4 variables
print("Données :\n", np.round(data, 3))
print()

# Statistiques globales
print("Moyenne globale :", np.round(data.mean(), 4))
print("Écart-type global:", np.round(data.std(), 4))
print()

# Statistiques par axe
print("Moyenne par colonne (axis=0) :", np.round(data.mean(axis=0), 4))
print("Moyenne par ligne  (axis=1) :", np.round(data.mean(axis=1), 4))
print()

# Autres statistiques
print("Médiane :", np.round(np.median(data), 4))
print("Percentiles [25, 50, 75] :", np.round(np.percentile(data, [25, 50, 75]), 4))
print("Minimum :", data.min(), "| Maximum :", data.max())
print("Argmin  :", data.argmin(), "| Argmax  :", data.argmax())

Données :
 [[ 0.126 -0.132  0.64   0.105]
 [-0.536  0.362  1.304  0.947]
 [-0.704 -1.265 -0.623  0.041]
 [-2.325 -0.219 -1.246 -0.732]
 [-0.544 -0.316  0.412  1.043]]

Moyenne globale : -0.1832
Écart-type global: 0.8509

Moyenne par colonne (axis=0) : [-0.7966 -0.3142  0.0974  0.2807]
Moyenne par ligne  (axis=1) : [ 0.1847  0.5193 -0.6378 -1.1305  0.1484]

Médiane : -0.1754
Percentiles [25, 50, 75] : [-0.6434 -0.1754  0.3741]
Minimum : -2.3250307746388343 | Maximum : 1.3040000451301372
Argmin  : 12 | Argmax  : 6

# Corrélation et histogramme
x = rng.standard_normal(100)
y = 0.7 * x + 0.3 * rng.standard_normal(100)

corrcoef = np.corrcoef(x, y)
print("Matrice de corrélation :\n", np.round(corrcoef, 4))

counts, bins = np.histogram(x, bins=10)
print("\nHistogramme de x :")
print("Comptages :", counts)
print("Bords des bacs :", np.round(bins, 2))

Matrice de corrélation :
 [[1.     0.9195]
 [0.9195 1.    ]]

Histogramme de x :
Comptages : [ 2  1 12 14 13 13 18  9 11  7]
Bords des bacs : [-2.25 -1.82 -1.4  -0.97 -0.55 -0.12  0.3   0.73  1.15  1.58  2.  ]

Entrées et sorties#

NumPy propose plusieurs fonctions pour persister des tableaux sur disque. Le choix du format dépend du contexte : interopérabilité avec d’autres outils, taille des données, nécessité de compression.

import tempfile
import os

# Créer un tableau à sauvegarder
a = np.arange(12).reshape(3, 4)

with tempfile.TemporaryDirectory() as tmpdir:
    # Format binaire .npy (un seul tableau)
    path_npy = os.path.join(tmpdir, "tableau.npy")
    np.save(path_npy, a)
    a_chargé = np.load(path_npy)
    print("Chargé depuis .npy :\n", a_chargé)

    # Format .npz (plusieurs tableaux, avec compression optionnelle)
    b = np.linspace(0, 1, 6)
    path_npz = os.path.join(tmpdir, "multi.npz")
    np.savez(path_npz, matrice=a, vecteur=b)
    données = np.load(path_npz)
    print("\nClés dans .npz :", list(données.keys()))
    print("matrice :\n", données["matrice"])
    print("vecteur :", données["vecteur"])

    # Format texte .csv
    path_csv = os.path.join(tmpdir, "matrice.csv")
    np.savetxt(path_csv, a, delimiter=",", fmt="%d", header="c0,c1,c2,c3")
    a_txt = np.loadtxt(path_csv, delimiter=",", skiprows=1)
    print("\nChargé depuis .csv :\n", a_txt)

Chargé depuis .npy :
 [[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]

Clés dans .npz : ['matrice', 'vecteur']
matrice :
 [[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
vecteur : [0.  0.2 0.4 0.6 0.8 1. ]

Chargé depuis .csv :
 [[ 0.  1.  2.  3.]
 [ 4.  5.  6.  7.]
 [ 8.  9. 10. 11.]]

Note

Le format .npy est le plus rapide à lire et écrire car il ne fait aucune conversion : les octets du tableau sont écrits directement sur disque avec un en-tête minimal décrivant le dtype et la shape. Pour des données volumineuses, préférer np.savez_compressed (format .npz avec compression zlib) ou les formats de haut niveau comme HDF5 (via h5py) ou Zarr, qui supportent l’accès partiel aux données et la lecture en parallèle.

Performance : vues, contiguïté et np.einsum#

Maîtriser les mécanismes de performance de NumPy permet d’écrire du code qui tire pleinement parti du matériel, notamment pour des calculs intensifs sur de grands tableaux.

Vues vs copies. Une vue partage la même zone mémoire que le tableau source. Les slices de base créent des vues ; le fancy indexing crée des copies. On peut vérifier si deux tableaux partagent la même mémoire avec np.shares_memory(a, b).

a = np.arange(12).reshape(3, 4)

vue = a[1:, ::2]        # slice → vue
copie = a[[0, 2], :]    # fancy indexing → copie

print("Vue partage la mémoire :", np.shares_memory(a, vue))
print("Copie partage la mémoire :", np.shares_memory(a, copie))

Vue partage la mémoire : True
Copie partage la mémoire : False

Contiguïté mémoire. Un tableau est dit C-contigu si les éléments sont ordonnés ligne par ligne (ordre C), et Fortran-contigu s’ils sont ordonnés colonne par colonne. NumPy crée par défaut des tableaux C-contigus. Les opérations qui traversent le tableau dans l’ordre de sa contiguïté sont plus rapides, car elles génèrent moins de défauts de cache.

a = np.arange(12).reshape(3, 4)
print("C-contigu :", a.flags['C_CONTIGUOUS'])
print("F-contigu :", a.flags['F_CONTIGUOUS'])

a_f = np.asfortranarray(a)
print("\nAprès asfortranarray :")
print("C-contigu :", a_f.flags['C_CONTIGUOUS'])
print("F-contigu :", a_f.flags['F_CONTIGUOUS'])

C-contigu : True
F-contigu : False

Après asfortranarray :
C-contigu : False
F-contigu : True

np.einsum — contraction tensorielle. einsum est une interface générale pour exprimer des sommes de produits sur des tableaux multidimensionnels via la notation d’Einstein. Elle permet d’exprimer en une seule ligne des opérations aussi variées que la transposition, le produit matriciel, la trace, les produits externes ou les contractions tensorielles.

A = np.random.default_rng(0).standard_normal((3, 4))
B = np.random.default_rng(1).standard_normal((4, 5))

# Produit matriciel : A_{ik} B_{kj} → C_{ij}
C_einsum = np.einsum('ik,kj->ij', A, B)
C_matmul = A @ B
print("Erreur max A@B vs einsum :", np.max(np.abs(C_einsum - C_matmul)))

# Trace d'une matrice
M = np.random.default_rng(2).standard_normal((4, 4))
print("Trace via einsum :", np.einsum('ii', M))
print("Trace via linalg :", np.trace(M))

# Produit externe de deux vecteurs
u = np.array([1.0, 2.0, 3.0])
v = np.array([4.0, 5.0])
print("Produit externe u⊗v :\n", np.einsum('i,j->ij', u, v))

Erreur max A@B vs einsum : 1.1102230246251565e-16
Trace via einsum : 2.2115886532394096
Trace via linalg : 2.2115886532394096
Produit externe u⊗v :
 [[ 4.  5.]
 [ 8. 10.]
 [12. 15.]]

Quand utiliser einsum ?

np.einsum est particulièrement utile lorsque l’opération souhaitée n’a pas d’équivalent direct dans NumPy, ou lorsque l’enchaînement de plusieurs opérations (matmul, sum, transpose) produirait des tableaux intermédiaires inutiles. Par exemple, le calcul de la trace du produit AB, np.trace(A @ B), alloue la matrice entière AB avant de sommer la diagonale. L’équivalent np.einsum('ij,ji->', A, B) effectue le calcul sans allouer AB. Sur de grands tableaux, l’argument optimize=True permet à NumPy de chercher le chemin de contraction optimal.

Résumé#

Ce chapitre a approfondi NumPy au-delà des bases, en couvrant les opérations avancées qui constituent le quotidien du travail en data science :

  • La manipulation des formes — reshape, ravel, flatten, transpose, np.newaxis, squeeze — permet de réorganiser les données sans copie, en exploitant le mécanisme des vues et des strides.

  • L”empilement et le découpage — concatenate, stack, vstack, hstack, split — permettent d’assembler et de désassembler des tableaux pour constituer des ensembles de données ou traiter des sous-blocs.

  • Le module np.linalg fournit toute l’algèbre linéaire nécessaire : résolution de systèmes (solve), déterminant, inverse, valeurs propres (eig) et décomposition SVD. L’opérateur @ est la notation recommandée pour le produit matriciel.

  • Les fonctions statistiques — mean, std, var, median, percentile, corrcoef, histogram — acceptent toutes un argument axis pour opérer sur une dimension précise.

  • Les formats .npy et .npz sont les formats natifs de NumPy, rapides et fiables pour persister des tableaux. np.savetxt et np.loadtxt servent pour l’interopérabilité avec les outils texte.

  • Comprendre la distinction vue vs copie et la contiguïté mémoire est essentiel pour éviter des bugs subtils et maximiser les performances. np.einsum offre une syntaxe élégante pour les contractions tensorielles complexes.

Dans le chapitre suivant, nous passerons à la visualisation avec Matplotlib et Seaborn, pour donner une forme graphique aux données que NumPy nous permet de manipuler si efficacement.

Contenu