Algèbre linéaire pour développeurs : vecteurs, matrices et NumPy en pratique

📍 Guide principal : Mathématiques appliquées à l’informatique : le socle pratique
Ce tutoriel approfondit la partie algèbre linéaire du guide principal — à parcourir d’abord pour situer pourquoi tous les frameworks d’apprentissage automatique reposent sur des opérations matricielles.

L’algèbre linéaire est passée du statut de matière mathématique théorique à celui de langue de travail des ingénieurs en données. Une image, un texte, un signal audio, une transaction client, un journal de logs : toutes ces choses se représentent comme des vecteurs ou des matrices, et toutes les transformations utiles (filtrage, projection, comparaison, classification) se réduisent à des opérations linéaires. Ce tutoriel pose, étape par étape, les outils pratiques pour manipuler vecteurs et matrices en NumPy, comprendre ce que font les opérations sous le capot, et reconnaître quand une boucle Python doit céder la place à un appel matriciel vectorisé.

Prérequis

• Python 3.10 ou plus récent
• NumPy 1.26 ou plus récent (toute version stable récente convient)
• Niveau attendu : intermédiaire — vous savez ce qu’est une liste, une boucle, un produit
• Temps estimé : 90 minutes

Étape 1 — Installer NumPy et créer ses premiers vecteurs

NumPy est la fondation de l’écosystème scientifique Python. Quasi toutes les bibliothèques de data science et d’apprentissage automatique (pandas, scikit-learn, scipy, matplotlib, PyTorch via torch.from_numpy) interagissent avec ses tableaux ndarray. L’installation est immédiate sur la plupart des plateformes — NumPy fournit des wheels pré-compilées qui embarquent les routines BLAS optimisées.

pip install numpy

Vérifiez la sortie Successfully installed numpy-X.Y.Z. Si l’installation déclenche une compilation depuis les sources (rare mais possible sur des architectures exotiques), vérifiez la présence de gcc ou utilisez la distribution Anaconda qui livre des binaires précompilés. Les versions récentes de NumPy supportent le standard du tableau Python 2024 et offrent des performances proches des bibliothèques C natives.

# vectors.py
import numpy as np

u = np.array([1.0, 2.0, 3.0])
v = np.array([4.0, 5.0, 6.0])

print("u + v =", u + v)
print("3u    =", 3 * u)
print("⟨u,v⟩ =", np.dot(u, v))
print("‖u‖   =", np.linalg.norm(u))
print("cos(u,v) =", np.dot(u, v) / (np.linalg.norm(u) * np.linalg.norm(v)))

L’exécution affiche la somme des deux vecteurs, le produit par un scalaire, le produit scalaire (32.0), la norme euclidienne de u (≈ 3.74) et le cosinus de l’angle entre les deux vecteurs (≈ 0.97, presque colinéaires). Le cosinus, sorte de produit scalaire normalisé, est l’une des métriques les plus utilisées en pratique : c’est elle qui mesure la similarité entre deux représentations vectorielles dans les moteurs de recommandation et dans la recherche sémantique.

Étape 2 — Manipuler des matrices et leurs dimensions

Une matrice est un tableau 2D. NumPy l’expose avec np.array sur une liste de listes ou avec des constructeurs comme np.zeros, np.ones, np.eye (identité), np.random.rand. La forme (shape) d’une matrice est cruciale : un produit A × B exige que le nombre de colonnes de A égale le nombre de lignes de B. Toute incohérence remonte sous forme de ValueError immédiate, ce qui est précieux pour repérer les bugs de dimension avant qu’ils ne se propagent.

# matrices.py
import numpy as np

A = np.array([[1, 2, 3], [4, 5, 6]])  # shape (2, 3)
B = np.array([[7, 8], [9, 10], [11, 12]])  # shape (3, 2)

print("A.shape  =", A.shape)
print("B.shape  =", B.shape)
print("A @ B    =\n", A @ B)
print("A.T      =\n", A.T)
print("eye(3)   =\n", np.eye(3))

Le script affiche le produit matriciel A @ B (matrice 2×2), la transposée de A (3×2), et la matrice identité 3×3. L’opérateur @, introduit en Python 3.5 spécifiquement pour le calcul matriciel (PEP 465), est plus lisible que np.matmul(A, B) ou np.dot(A, B). Pour le produit élément-par-élément (Hadamard), utilisez l’opérateur * standard, qui exige des shapes identiques ou compatibles via le broadcasting.

Étape 3 — Comprendre le broadcasting

Le broadcasting est le mécanisme qui permet à NumPy d’appliquer une opération entre tableaux de formes différentes en étendant virtuellement la dimension la plus petite. C’est ce qui rend l’écriture vectorisée concise : ajouter un vecteur ligne à chaque ligne d’une matrice, normaliser des colonnes, soustraire la moyenne — tous ces gestes deviennent des opérations en une ligne.

# broadcasting.py
import numpy as np

X = np.array([[1.0, 2.0, 3.0],
              [4.0, 5.0, 6.0],
              [7.0, 8.0, 9.0]])

mean_per_column = X.mean(axis=0)
X_centered = X - mean_per_column

std_per_column = X.std(axis=0)
X_standardized = X_centered / std_per_column

print("Moyenne par colonne :", mean_per_column)
print("X centré :\n", X_centered)
print("X standardisé :\n", X_standardized)
print("Vérif moyenne après centrage :", X_centered.mean(axis=0))

Le script centre chaque colonne (moyenne 0) puis la standardise (variance 1) — une étape de préparation systématique pour la régression linéaire et les réseaux de neurones. La ligne X - mean_per_column applique la soustraction d’un vecteur de taille 3 à chaque ligne d’une matrice 3×3 sans qu’on ait à écrire une boucle. Les règles du broadcasting sont précises (les shapes sont alignées de droite à gauche, les dimensions de taille 1 sont étendues, le reste doit correspondre), mais en pratique, le réflexe est de vérifier les shapes avant l’opération avec print(X.shape, mean_per_column.shape).

Étape 4 — Vectorisation : la différence qui change tout

La promesse principale de NumPy est la vectorisation : remplacer une boucle Python par un appel à une routine BLAS écrite en C ou Fortran qui exploite les instructions SIMD du processeur (SSE, AVX). Le gain typique sur un produit matriciel est d’un facteur 50 à 200, parfois davantage sur des architectures récentes. La règle d’or : si une opération se voit naturellement comme une opération sur tableau, la version vectorisée existe presque toujours et il faut la chercher.

# vectorize.py
import numpy as np
import time

n = 1_000_000
a = np.random.rand(n)
b = np.random.rand(n)

t0 = time.perf_counter()
result_loop = sum(a[i] * b[i] for i in range(n))
t1 = time.perf_counter()

t2 = time.perf_counter()
result_vec = np.dot(a, b)
t3 = time.perf_counter()

print(f"Boucle Python : {t1-t0:.3f}s")
print(f"NumPy dot     : {t3-t2:.5f}s")
print(f"Accélération  : {(t1-t0)/(t3-t2):.0f}x")

L’exécution affiche typiquement un facteur d’accélération de 100 à 500 en faveur de la version NumPy sur un million d’éléments. Cette différence n’est pas un détail : elle sépare un script qui finit en quelques millisecondes d’un script qui dure plusieurs minutes. Sur des matrices, l’écart s’amplifie encore parce que NumPy utilise des routines BLAS hautement optimisées (OpenBLAS, MKL) qui exploitent la mémoire cache et le parallélisme du processeur.

Étape 5 — Résoudre un système linéaire

Beaucoup de problèmes pratiques se ramènent à résoudre A × x = b où A et b sont connus et x est l’inconnu : régression linéaire, équilibrage de réactions chimiques, calcul de courants dans un circuit, équilibrage statique d’un système mécanique. NumPy expose numpy.linalg.solve qui choisit en interne la décomposition la plus efficace selon la structure de A.

# solve.py
import numpy as np

# Système : 2x + y = 5 ; x - 3y = -8
A = np.array([[2.0, 1.0], [1.0, -3.0]])
b = np.array([5.0, -8.0])
x = np.linalg.solve(A, b)
print("Solution x =", x)

# Vérification
print("A @ x  =", A @ x, "  (devrait être", b, ")")
print("Résidu :", np.linalg.norm(A @ x - b))

La solution affichée doit être [1.0, 3.0] et le résidu, distance entre A @ x et b, doit être de l’ordre de 10⁻¹⁵ — la précision de la double précision flottante (IEEE 754 binary64). N’inversez jamais A explicitement avec np.linalg.inv(A) @ b : c’est plus lent, plus fragile numériquement, et inutile dans la quasi-totalité des cas. Préférez toujours solve, qui factorise A une fois et propage la solution.

Étape 6 — Décomposition en valeurs singulières (SVD)

La SVD est la décomposition matricielle la plus fondamentale en data science. Toute matrice M de taille m×n se décompose en U × Σ × VT où U et V sont orthogonales et Σ est diagonale de valeurs positives décroissantes (les valeurs singulières). Les valeurs singulières quantifient l’importance de chaque direction principale ; en n’en gardant que les k plus grandes, on obtient la meilleure approximation de M de rang k au sens des moindres carrés (théorème d’Eckart-Young, 1936).

# svd.py
import numpy as np

# Matrice utilisateurs × films, scores de 1 à 5
M = np.array([
    [5, 4, 1, 1],
    [4, 5, 1, 2],
    [1, 1, 5, 4],
    [1, 2, 4, 5],
    [2, 1, 4, 5],
], dtype=float)

U, S, Vt = np.linalg.svd(M, full_matrices=False)
print("Valeurs singulières :", np.round(S, 3))

# Approximation rang 2 (compression / dénoising)
k = 2
M_approx = U[:, :k] @ np.diag(S[:k]) @ Vt[:k, :]
print("Erreur Frobenius :", np.linalg.norm(M - M_approx))
print("M approchée (rang 2) :\n", np.round(M_approx, 2))

Le script extrait les valeurs singulières puis reconstruit la matrice avec seulement les deux plus grandes. La matrice reconstruite reste très proche de l’originale parce que la matrice de notes a effectivement un rang faible (deux groupes d’utilisateurs, deux types de films). C’est exactement le mécanisme derrière les premiers moteurs de recommandation collaborative et derrière l’analyse en composantes principales (PCA), qui n’est rien d’autre qu’une SVD appliquée à une matrice centrée.

Étape 7 — Construire un mini moteur de similarité par cosinus

Pour conclure, mettons les outils en pratique sur un cas concret : trouver le document le plus proche d’une requête dans un corpus représenté par des vecteurs (typiquement des embeddings produits par un modèle de plongement). C’est le cœur de la recherche sémantique et des moteurs de recommandation par contenu.

# cosine_search.py
import numpy as np

# Corpus de 5 documents en dimension 4 (vecteurs simulés)
corpus = np.random.RandomState(42).rand(5, 4)
norms = np.linalg.norm(corpus, axis=1, keepdims=True)
corpus_normalized = corpus / norms  # vecteurs unitaires

# Requête
query = np.random.RandomState(7).rand(4)
query_normalized = query / np.linalg.norm(query)

# Similarités = produit scalaire avec chaque ligne (vectorisé)
similarities = corpus_normalized @ query_normalized
ranking = np.argsort(-similarities)
print("Similarités :", np.round(similarities, 3))
print("Top 3 documents :", ranking[:3].tolist())

Le script normalise chaque ligne du corpus, normalise la requête, calcule en une seule opération vectorisée le cosinus de la requête avec chacun des 5 documents, et trie pour produire le top-K. Cette même structure se généralise à des corpus de millions de documents en dimension 768 ou 1536 (taille typique d’un embedding moderne) : il suffit de remplacer corpus par une matrice de la bonne taille. Pour les très grands corpus, des bibliothèques spécialisées comme FAISS de Meta indexent les vecteurs pour ramener la recherche à O(log N) au lieu de O(N).

Erreurs fréquentes

Tutoriels associés

• Notation big-O et analyse de complexité : mesurer le coût d’un algorithme
• Théorie des graphes en pratique : modéliser, BFS, DFS et Dijkstra

Pour aller plus loin

• 🔝 Retour au guide principal : Mathématiques appliquées à l’informatique : le socle pratique
• Documentation officielle NumPy — module linalg : numpy.org/doc/stable/reference/routines.linalg.html
• Gilbert Strang — Introduction to Linear Algebra (MIT) : cours et livre de référence
• Cours MIT 18.06 — Linear Algebra : ocw.mit.edu/courses/18-06
• Trefethen, Bau — Numerical Linear Algebra, SIAM, 1997

FAQ

Pourquoi NumPy est-il si rapide par rapport à du Python pur ?

Trois raisons combinées : les tableaux ndarray sont stockés contigüment en mémoire (cache-friendly), les opérations sont déléguées à des routines BLAS écrites en C ou Fortran avec instructions SIMD, et le coût des appels Python est amorti sur des millions d’éléments. Le gain typique est d’un facteur 50 à 500 selon l’opération.

Quelle bibliothèque pour le GPU ?

PyTorch et JAX permettent d’écrire du code très proche de NumPy qui s’exécute sur GPU (CUDA, ROCm, Metal). CuPy est l’alternative la plus directement compatible avec NumPy pour le GPU NVIDIA. Pour des matrices de plusieurs centaines de Mo et des opérations massives, le GPU peut donner un facteur 10 à 100 supplémentaire.

Comment diagnostiquer une instabilité numérique ?

Le nombre de conditionnement np.linalg.cond(A) mesure la sensibilité de la solution aux erreurs sur l’entrée. Un conditionnement supérieur à 10¹⁰ annonce des problèmes de précision sur des matrices flottantes en double précision. La régularisation (Tikhonov, ridge) ou un changement de représentation peut souvent restaurer la stabilité.

Quel format pour stocker des vecteurs en dimension élevée ?

Pour quelques millions de vecteurs en dimension 1024, numpy.savez ou numpy.memmap conviennent. Au-delà, des formats spécialisés comme HDF5 (via h5py), Zarr ou Parquet offrent compression et accès partiel. Pour la recherche par similarité à grande échelle, FAISS, Milvus ou Qdrant indexent les vecteurs en mémoire ou sur disque.

Quand passer à une bibliothèque spécialisée comme JAX ?

Quand on a besoin de différentiation automatique, de compilation just-in-time (jit) ou d’exécution sur TPU. JAX expose une API très proche de NumPy mais ajoute grad, jit, vmap et pmap qui permettent de transformer du code numérique en code dérivable, compilé et parallèle.