Module L10 · Partie G · Mathématiques de l’informatique

Algèbre linéaire numérique : le langage de l’IA et de la robotique.

Une image est une matrice. Un réseau de neurones est une suite de produits matriciels. La position d’un bras robot est une chaîne de matrices de rotation. Une régression est une projection orthogonale. Ce module apprend à calculer avec des vecteurs et des matrices — avec NumPy et en réimplémentant les algorithmes clés — pour que les modules d’IA et de robotique n’aient plus de mystère mathématique.

Durée : 3 séances · Prérequis : séances 14 et 19 (vecteurs, NumPy). Objectifs : produit matriciel et ses coûts, résolution de systèmes (Gauss, LU), moindres carrés, transformations géométriques 2D/3D, valeurs propres et SVD (intuition + usage), conditionnement et erreurs numériques.

Ce que vous saurez faire à la fin

Écrire l’élimination de Gauss et expliquer pourquoi NumPy fait mieux.
Résoudre un problème de moindres carrés (calibration de capteur, ajustement de droite/parabole).
Composer rotations et translations en matrices homogènes pour un bras robot.
Utiliser la SVD pour compresser une image et comprendre l’ACP.

Références : Linear Algebra and Its Applications (Strang) et ses cours MIT 18.06 en vidéo, 3Blue1Brown « Essence of linear algebra », programme de mathématiques MPSI/MP.

Fiche de cours · Définitions

Algèbre linéaire : définitions

Définition (espace vectoriel, combinaison, famille libre, base). ℝⁿ avec l’addition et la multiplication par un scalaire. Une combinaison linéaire est Σ λ_i v_i. La famille (v₁, …, v_k) est libre si Σ λ_i v_i = 0 ⇒ tous les λ_i = 0 ; génératrice de E si tout vecteur de E est combinaison ; une base est libre et génératrice ; toutes les bases ont le même cardinal, la dimension.

Définition (application linéaire, matrice, rang, noyau). f : ℝⁿ → ℝ^m est linéaire si f(u + λv) = f(u) + λf(v). Dans des bases fixées, f est représentée par une matrice A (m×n) : f(x) = Ax, la colonne j de A est l’image du j-ième vecteur de base. Im(A) = {Ax} (engendré par les colonnes), rang(A) = dim Im(A) ; Ker(A) = {x : Ax = 0}.

Définition (produit scalaire, norme, orthogonalité). ⟨x, y⟩ = xᵀy = Σ x_i y_i ; ‖x‖ = √⟨x, x⟩ ; x ⊥ y si ⟨x, y⟩ = 0. Une matrice Q est orthogonale si QᵀQ = I (colonnes orthonormées) : elle conserve les normes et les angles (rotations, réflexions).

Définition (déterminant, inverse). det(A) : volume signé de l’image du cube unité ; A est inversible ⇔ det(A) ≠ 0 ⇔ rang(A) = n ⇔ Ker(A) = {0}. det(AB) = det(A)det(B), det(Aᵀ) = det(A).

Définition (valeur propre, vecteur propre). Av = λv avec v ≠ 0. λ est racine du polynôme caractéristique det(A − λI). A est diagonalisable s’il existe une base de vecteurs propres : A = PDP⁻¹. Une matrice symétrique réelle est toujours diagonalisable en base orthonormée (théorème spectral).

Définition (SVD). Toute matrice A (m×n) s’écrit A = UΣVᵀ avec U, V orthogonales et Σ diagonale à coefficients σ₁ ≥ σ₂ ≥ … ≥ 0 (valeurs singulières). rang(A) = nombre de σ_i > 0 ; ‖A‖₂ = σ₁ ; conditionnement κ(A) = σ₁/σ_n.

Fiche de cours · Formules

Formules à connaître

(AB)_ij = Σ_k A_ik B_kj — coût O(mnp) pour (m×n)(n×p) ; (AB)ᵀ = BᵀAᵀ ; (AB)⁻¹ = B⁻¹A⁻¹

Théorème du rang : dim Ker(A) + rang(A) = n (nombre de colonnes)

Rotation plane R(θ) = [[cos θ, −sin θ], [sin θ, cos θ]] ; R(θ)R(φ) = R(θ+φ) ; R(θ)⁻¹ = R(θ)ᵀ = R(−θ) coordonnées homogènes : [[R, t], [0, 1]] enchaîne rotation et translation par un seul produit

Moindres carrés : x* = argmin ‖Ax − b‖² vérifie AᵀA x* = Aᵀb (équations normales) ; x* = (AᵀA)⁻¹Aᵀb si A est de rang n

Projection orthogonale sur Im(A) : P = A(AᵀA)⁻¹Aᵀ ; P² = P, Pᵀ = P ; sur un vecteur unitaire u : P = uuᵀ

2×2 : det = ad − bc ; A⁻¹ = (1/det)·[[d, −b], [−c, a]] ; λ² − tr(A)λ + det(A) = 0 ; tr(A) = Σ λ_i, det(A) = Π λ_i

Aᵏ = PDᵏP⁻¹ ; ‖Ax‖ ≤ ‖A‖‖x‖ ; erreur relative sur la solution de Ax = b ≤ κ(A) × erreur relative sur b

Fiche de cours · Théorèmes et démonstrations

Démonstrations à savoir refaire (1/2)

Théorème 1 (équations normales). x* minimise ‖Ax − b‖² si et seulement si AᵀA x* = Aᵀb. Le résidu r = b − Ax* est orthogonal aux colonnes de A.

Posons J(x) = ‖Ax − b‖² = (Ax − b)ᵀ(Ax − b) = xᵀAᵀAx − 2bᵀAx + bᵀb. Son gradient est ∇J = 2AᵀAx − 2Aᵀb (dérivée de xᵀMx est 2Mx pour M symétrique, de cᵀx est c). J est convexe (AᵀA est semi-définie positive : xᵀAᵀAx = ‖Ax‖² ≥ 0), donc ∇J = 0 caractérise les minima : AᵀA x* = Aᵀb. Cela s’écrit Aᵀ(b − Ax*) = 0 : le résidu est orthogonal à chaque colonne de A — Ax* est la projection orthogonale de b sur Im(A). Interprétation géométrique : le point de Im(A) le plus proche de b est le pied de la perpendiculaire (Pythagore : pour tout y = Ax, ‖b − y‖² = ‖b − Ax*‖² + ‖Ax* − y‖² ≥ ‖b − Ax*‖²).

Théorème 2 (vecteurs propres de valeurs propres distinctes sont libres).

Récurrence sur k. k = 1 : v₁ ≠ 0. Hérédité : supposons Σ_i≤k+1 α_i v_i = 0. Appliquons A : Σ α_i λ_i v_i = 0. Multiplions la première relation par λ_k+1 et soustrayons : Σ_i≤k α_i(λ_i − λ_k+1) v_i = 0. Par hypothèse de récurrence, α_i(λ_i − λ_k+1) = 0, et λ_i ≠ λ_k+1 donne α_i = 0 pour i ≤ k ; puis α_k+1 v_k+1 = 0 donne α_k+1 = 0. Conséquence : n valeurs propres distinctes ⇒ diagonalisable.

Théorème 3 (une matrice symétrique a des vecteurs propres orthogonaux). Si A = Aᵀ, Av = λv, Aw = μw avec λ ≠ μ, alors v ⊥ w.

λ⟨v, w⟩ = (Av)ᵀw = vᵀAᵀw = vᵀAw = μ⟨v, w⟩, donc (λ − μ)⟨v, w⟩ = 0 et ⟨v, w⟩ = 0. (Le théorème spectral complet ajoute que les valeurs propres sont réelles et qu’il existe une base orthonormée de vecteurs propres — c’est ce qui fonde l’ACP : les axes principaux sont les vecteurs propres de la matrice de covariance, symétrique.)

Fiche de cours · Théorèmes et démonstrations

Démonstrations à savoir refaire (2/2)

Théorème 4 (élimination de Gauss : correction et coût). Les opérations élémentaires sur les lignes (échange, multiplication par un scalaire non nul, ajout d’un multiple d’une ligne à une autre) ne changent pas l’ensemble des solutions de Ax = b ; la réduction à une forme triangulaire coûte ~(2/3)n³ opérations.

Chaque opération est inversible (échange ↔ échange ; ×c ↔ ×1/c ; L_i += cL_j ↔ L_i −= cL_j) et une solution du système vérifie toute combinaison de ses équations : les deux systèmes ont les mêmes solutions. Coût : éliminer sous le pivot k met à jour (n − k) lignes de (n − k) coefficients, soit Σ_k 2(n − k)² ≈ 2·n³/3. La remontée coûte O(n²). Le choix du pivot de plus grand module (pivot partiel) évite de diviser par de petits nombres, source d’amplification des erreurs d’arrondi.

Théorème 5 (puissance itérée). Si A est diagonalisable avec |λ₁| > |λ₂| ≥ … et x₀ a une composante non nulle sur v₁, alors A^kx₀ / ‖A^kx₀‖ converge vers ±v₁, à la vitesse (|λ₂|/|λ₁|)^k.

Écrivons x₀ = Σ c_i v_i avec c₁ ≠ 0. Alors A^kx₀ = Σ c_i λ_i^k v_i = λ₁^k(c₁v₁ + Σ_i≥2 c_i(λ_i/λ₁)^k v_i). Les termes i ≥ 2 tendent vers 0 comme (|λ₂|/|λ₁|)^k. Après normalisation, la direction tend vers v₁. C’est l’algorithme de PageRank (vecteur propre dominant de la matrice de transition) et le point de départ de l’itération QR.

Théorème 6 (meilleure approximation de rang k, Eckart-Young). Si A = UΣVᵀ, alors A_k = Σ_i≤k σ_i u_i v_iᵀ minimise ‖A − B‖ parmi les matrices B de rang ≤ k, et ‖A − A_k‖₂ = σ_k+1.

(Idée pour la norme spectrale.) ‖A − A_k‖₂ = ‖Σ_i>k σ_iu_iv_iᵀ‖₂ = σ_k+1. Pour B de rang ≤ k, Ker(B) est de dimension ≥ n − k ; l’espace engendré par v₁, …, v_k+1 est de dimension k + 1 ; ils s’intersectent en un vecteur unitaire z. Alors ‖A − B‖₂ ≥ ‖(A − B)z‖ = ‖Az‖ = ‖Σ_i≤k+1 σ_i⟨v_i, z⟩u_i‖ ≥ σ_k+1‖z‖ = σ_k+1. C’est la justification de la compression d’image par SVD et de l’ACP.

Fiche de cours · Méthodes

Méthodes et pièges

Méthode — résoudre Ax = b numériquement. Jamais inv(A) @ b : utiliser np.linalg.solve (LU avec pivot, plus stable, 3× moins cher). Système surdéterminé : np.linalg.lstsq (QR/SVD, plus stable que les équations normales dont le conditionnement est κ(A)²). Vérifier le résidu ‖Ax − b‖ et le conditionnement np.linalg.cond.

Méthode — enchaîner des transformations. Écrire chaque transformation en homogène 3×3 (2D) ou 4×4 (3D) ; la transformation composée « d’abord T₁ puis T₂ » est T₂·T₁ (ordre inverse de l’écriture). Pour un bras : T_base→outil = T₁(θ₁)·T₂(θ₂)·… ; la position de l’outil est la dernière colonne.

Méthode — lire une SVD. Les colonnes de U : directions principales de l’espace d’arrivée ; de V : de l’espace de départ ; σ_i : gains. Rang numérique = nombre de σ_i > ε·σ₁. Énergie conservée par k composantes : Σ_i≤k σ_i² / Σ σ_i².

Pièges : confondre A·B et B·A ; comparer des flottants avec == (utiliser allclose) ; inverser une matrice presque singulière (κ ≈ 10¹⁶ : résultat sans signification) ; oublier de centrer les données avant l’ACP ; degrés vs radians ; * (élément par élément) vs @ (produit matriciel) en NumPy.

Fiche de cours · Exercices corrigés

Exercices corrigés

Exercice 1. Un robot au point (2, 1) orienté à 90° avance de 3 puis tourne de −90°. Donner sa pose finale par produit de matrices homogènes.

Correction. Pose initiale T₀ = [[R(90°), (2, 1)ᵀ], [0, 1]] avec R(90°) = [[0, −1], [1, 0]]. Avancer de 3 dans le repère du robot : T_a = [[I, (3, 0)ᵀ], [0, 1]]. Tourner de −90° : T_r = [[R(−90°), 0], [0, 1]]. Pose finale T = T₀·T_a·T_r (les mouvements sont exprimés dans le repère courant, donc à droite). Translation : R(90°)·(3, 0)ᵀ + (2, 1)ᵀ = (0, 3) + (2, 1) = (2, 4). Orientation : 90° − 90° = 0°. Pose finale (2, 4, 0°). Vérification intuitive : orienté vers le haut, avancer de 3 monte de 3.

Exercice 2. Ajuster une droite y = ax + b aux points (0, 1), (1, 3), (2, 5), (3, 8) par moindres carrés, à la main.

Correction. A = [[0,1],[1,1],[2,1],[3,1]], x = (a, b)ᵀ, b = (1, 3, 5, 8)ᵀ. AᵀA = [[Σx², Σx], [Σx, n]] = [[14, 6], [6, 4]] ; Aᵀb = [Σxy, Σy] = [0+3+10+24, 17] = [37, 17]. Résoudre : det = 56 − 36 = 20 ; a = (4·37 − 6·17)/20 = (148 − 102)/20 = 2,3 ; b = (14·17 − 6·37)/20 = (238 − 222)/20 = 0,8. Droite y = 2,3x + 0,8. Résidus : (0,2 ; −0,1 ; −0,4 ; 0,3), somme nulle (car la colonne de 1 est dans A : le résidu lui est orthogonal).

Exercice 3. Soit A = [[2, 1], [1, 2]]. Diagonaliser A, calculer A¹⁰ et interpréter géométriquement.

Correction. Polynôme caractéristique λ² − 4λ + 3 = (λ − 1)(λ − 3) : λ = 1, 3. Vecteurs propres : (A − I)v = 0 ⇒ v₁ = (1, −1)ᵀ ; (A − 3I)v = 0 ⇒ v₂ = (1, 1)ᵀ, orthogonaux (A symétrique). P = [[1, 1], [−1, 1]], D = diag(1, 3), P⁻¹ = (1/2)[[1, −1], [1, 1]]. A¹⁰ = P·diag(1, 3¹⁰)·P⁻¹ = (1/2)[[1 + 3¹⁰, −1 + 3¹⁰], [−1 + 3¹⁰, 1 + 3¹⁰]] = [[29525, 29524], [29524, 29525]]. Géométriquement : A étire d’un facteur 3 le long de la diagonale (1, 1) et laisse l’anti-diagonale inchangée ; A¹⁰ écrase presque tout vecteur sur la diagonale (puissance itérée, Théorème 5).

01 / Vecteurs et matrices

Le produit matriciel : définition, coût, sens

Trois façons de lire A·x

1) Ligne par ligne : chaque composante du résultat est un produit scalaire (ligne de A) · x. 2) Colonne par colonne : A·x est une combinaison linéaire des colonnes de A, avec les coefficients x. 3) Application linéaire : A transforme l’espace (rotation, étirement, projection) et A·x est l’image de x. La deuxième lecture explique pourquoi « A·x = b a une solution ⇔ b est dans l’espace engendré par les colonnes ». NumPy appelle BLAS (code C/Fortran optimisé pour le cache et le SIMD) : 1000 fois plus rapide que la boucle Python, et c’est exactement ce que font les GPU en parallèle massif pour l’IA.

01 / Vecteurs et matrices

Transformations géométriques : rotations, homogènes, chaînes cinématiques

import numpy as np
def rot(theta):
    c, s = np.cos(theta), np.sin(theta)
    return np.array([[c, -s], [s, c]])

R = rot(np.pi / 2)
print(R @ np.array([1, 0]))                       # (1,0) tourné de 90° → (0,1)
print(np.round(R @ R.T, 6))                       # une rotation est orthogonale : R·Rᵀ = I, donc R⁻¹ = Rᵀ
print(np.linalg.det(R))                            # déterminant 1 : conserve les aires et l'orientation

# Coordonnées homogènes : une translation devient aussi une matrice → on peut tout composer
def T(dx, dy): return np.array([[1, 0, dx], [0, 1, dy], [0, 0, 1.]])
def Rh(theta): M = np.eye(3); M[:2, :2] = rot(theta); return M

# Bras à deux segments (module 3e « Le bras ») : L1 = 80, L2 = 60, angles a1, a2
def bras(a1, a2, L1=80, L2=60):
    base_vers_coude = Rh(a1) @ T(L1, 0)             # tourner de a1, puis avancer de L1
    coude_vers_main = Rh(a2) @ T(L2, 0)
    M = base_vers_coude @ coude_vers_main            # composition : de la base à la main
    coude = base_vers_coude @ np.array([0, 0, 1.])
    main = M @ np.array([0, 0, 1.])
    return coude[:2], main[:2]
c, m = bras(np.radians(30), np.radians(30))
print("coude", c.round(1), "main", m.round(1))
# Vérification avec la formule directe : x = L1 cos a1 + L2 cos(a1 + a2)
print(80 * np.cos(np.radians(30)) + 60 * np.cos(np.radians(60)))

L’ordre compte : Rh(a) @ T(L, 0) tourne puis translate dans le repère tourné (c’est ce qu’un bras fait) ; T @ Rh ferait l’inverse. En 3D, on ajoute une dimension (matrices 4×4) et les rotations se paramètrent par angles d’Euler ou quaternions : c’est ce que ROS appelle un « transform » (tf2). La convention de Denavit-Hartenberg standardise ces chaînes pour les robots industriels.

02 / Systèmes linéaires

Élimination de Gauss : résoudre A·x = b

import numpy as np

def gauss(A, b):
    """Résout A x = b par élimination avec pivot partiel. O(n³/3) multiplications."""
    A = A.astype(float).copy(); b = b.astype(float).copy(); n = len(b)
    for k in range(n):
        p = k + np.argmax(np.abs(A[k:, k]))          # pivot partiel : le plus grand en valeur absolue
        if abs(A[p, k]) < 1e-12: raise ValueError("matrice singulière")
        A[[k, p]], b[[k, p]] = A[[p, k]], b[[p, k]]  # échanger les lignes
        for i in range(k + 1, n):
            f = A[i, k] / A[k, k]
            A[i, k:] -= f * A[k, k:]; b[i] -= f * b[k]
    x = np.zeros(n)                                   # remontée : système triangulaire supérieur
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - A[i, i + 1:] @ x[i + 1:]) / A[i, i]
    return x

A = np.array([[2., 1., -1.], [-3., -1., 2.], [-2., 1., 2.]]); b = np.array([8., -11., -3.])
x = gauss(A, b); print(x, np.allclose(A @ x, b))
print(np.linalg.solve(A, b))

# Pourquoi le pivot : sans lui, un petit pivot amplifie les erreurs d'arrondi
eps = 1e-17
A2 = np.array([[eps, 1.], [1., 1.]]); b2 = np.array([1., 2.])
print("solve :", np.linalg.solve(A2, b2), "(vraie solution ≈ [1, 1])")
# Ne JAMAIS faire x = inv(A) @ b : 3× plus lent et moins précis que solve
print(np.allclose(np.linalg.inv(A) @ b, x))

LU, et pourquoi on ne calcule jamais l’inverse

Gauss factorise A = L·U (triangulaire inférieure × supérieure) ; résoudre A·x = b devient deux substitutions en O(n²). Si on a plusieurs b (plusieurs mesures, un filtre de Kalman à chaque pas), on factorise une fois en O(n³) et on résout chaque fois en O(n²). L’inverse explicite coûte 3 fois plus, remplit de zéros les matrices creuses, et amplifie les erreurs. Règle : « inverser une matrice » se traduit toujours par solve. Pour des matrices symétriques définies positives (covariances, moindres carrés), Cholesky est 2 fois plus rapide que LU.

02 / Systèmes linéaires

Moindres carrés : ajuster un modèle à des mesures bruitées

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)

# Calibration d'un capteur de distance infrarouge : tension → distance, mesures bruitées
d_vrai = np.linspace(10, 80, 25)
tension = 27 / d_vrai + 0.1 + np.random.normal(0, 0.03, 25)     # modèle physique ≈ a/d + b

# On veut d en fonction de la tension. Modèle linéaire dans les paramètres : d ≈ c0 + c1/v + c2/v²
X = np.column_stack([np.ones_like(tension), 1 / tension, 1 / tension**2])   # matrice de conception (25 × 3)
# Système surdéterminé X c ≈ d : pas de solution exacte. Minimiser ‖X c − d‖² ⇒ équations normales XᵀX c = Xᵀ d
c_normal = np.linalg.solve(X.T @ X, X.T @ d_vrai)
c_lstsq, residus, rang, sv = np.linalg.lstsq(X, d_vrai, rcond=None)    # méthode stable (QR/SVD)
print("coefficients :", c_lstsq.round(3), "| rang", rang)
print("écart entre les deux méthodes :", np.abs(c_normal - c_lstsq).max())
erreur = X @ c_lstsq - d_vrai
print(f"erreur RMS : {np.sqrt(np.mean(erreur**2)):.2f} cm")

v = np.linspace(tension.min(), tension.max(), 200)
plt.figure(figsize=(6, 3.5))
plt.plot(tension, d_vrai, "o", label="mesures")
plt.plot(v, c_lstsq[0] + c_lstsq[1] / v + c_lstsq[2] / v**2, label="modèle ajusté")
plt.xlabel("tension (V)"); plt.ylabel("distance (cm)"); plt.legend(); plt.show()

La géométrie : projeter

Les colonnes de X engendrent un sous-espace ; d n’y est pas (bruit). La meilleure approximation X·c est la projection orthogonale de d sur ce sous-espace, donc le résidu d − X·c est orthogonal aux colonnes : Xᵀ(d − X·c) = 0, d’où les équations normales. C’est la régression linéaire, celle du module L13, avec une solution exacte en une ligne. Les équations normales élèvent le conditionnement au carré ; lstsq (via QR ou SVD) est numériquement plus sûr — préférez-le toujours.

03 / Valeurs propres et SVD

Valeurs propres : les directions que la matrice ne fait qu’étirer

import numpy as np
A = np.array([[2., 1.], [1., 2.]])
val, vec = np.linalg.eig(A)
print("valeurs propres :", val, "\nvecteurs propres (colonnes) :\n", vec.round(3))
v = vec[:, 0]
print("A v = λ v ?", np.allclose(A @ v, val[0] * v))

# Méthode de la puissance itérée : trouver la plus grande valeur propre en itérant A·x
x = np.random.rand(2)
for _ in range(50):
    x = A @ x; x /= np.linalg.norm(x)
print("puissance itérée :", x.round(3), "λ ≈", (x @ A @ x).round(3))

# Application : une chaîne de Markov (module 3e « prévision météo ») ; l'état stationnaire est le vecteur propre pour λ = 1
P = np.array([[0.8, 0.2],      # beau → beau 0.8, beau → pluie 0.2
              [0.4, 0.6]])     # pluie → beau 0.4
etat = np.array([1., 0.])
for _ in range(30): etat = etat @ P
print("distribution à long terme :", etat.round(4), "(2/3, 1/3)")
val, vec = np.linalg.eig(P.T)
stationnaire = np.real(vec[:, np.argmin(abs(val - 1))]); print(stationnaire / stationnaire.sum())

# PageRank de Google : le vecteur propre principal de la matrice du web
liens = {0: [1, 2], 1: [2], 2: [0], 3: [2]}          # page → pages pointées
n = 4; M = np.zeros((n, n))
for i, js in liens.items():
    for j in js: M[j, i] = 1 / len(js)
G = 0.85 * M + 0.15 / n                              # amortissement : saut aléatoire
r = np.ones(n) / n
for _ in range(100): r = G @ r
print("PageRank :", r.round(3), "→ page la plus importante :", r.argmax())

Une matrice symétrique a des valeurs propres réelles et des vecteurs propres orthogonaux (théorème spectral) : c’est le cas des covariances, des laplaciens de graphes, des hessiennes — d’où son omniprésence en IA (ACP, optimisation) et en physique (modes de vibration d’un robot).

03 / Valeurs propres et SVD

SVD : compresser une image, comprendre l’ACP

import numpy as np
import matplotlib.pyplot as plt

# Une « image » synthétique 64×64 : disque + bandes + bruit
y, x = np.mgrid[0:64, 0:64]
img = ((x - 32)**2 + (y - 32)**2 < 300).astype(float) + 0.4 * np.sin(x / 4) + 0.1 * np.random.rand(64, 64)

U, s, Vt = np.linalg.svd(img)                         # img = U · diag(s) · Vt
print("valeurs singulières (premières) :", s[:6].round(1), "… énergie dans les 8 premières :",
      f"{(s[:8]**2).sum() / (s**2).sum():.1%}")

fig, axes = plt.subplots(1, 4, figsize=(9, 2.5))
for ax, k in zip(axes, [1, 4, 12, 64]):
    approx = (U[:, :k] * s[:k]) @ Vt[:k]              # somme des k premières composantes de rang 1
    ax.imshow(approx, cmap="gray"); ax.set_title(f"rang {k} : {k * (64 + 64 + 1)} nombres"); ax.axis("off")
plt.show()

# ACP (analyse en composantes principales) = SVD des données centrées
np.random.seed(1)
donnees = np.random.multivariate_normal([0, 0], [[3, 2.5], [2.5, 3]], 300)   # nuage allongé
centre = donnees - donnees.mean(axis=0)
U, s, Vt = np.linalg.svd(centre, full_matrices=False)
print("axe principal :", Vt[0].round(3), "| variance expliquée :", (s**2 / (s**2).sum()).round(3))
plt.figure(figsize=(4, 4)); plt.scatter(*centre.T, s=6)
for k in range(2): plt.plot([0, 3 * s[k] / 10 * Vt[k, 0]], [0, 3 * s[k] / 10 * Vt[k, 1]], lw=3)
plt.axis("equal"); plt.title("ACP : les deux axes principaux"); plt.show()

Ce que dit la SVD

Toute matrice A (même rectangulaire) s’écrit U·Σ·Vᵀ avec U, V orthogonales et Σ diagonale positive : A envoie une base orthonormée (colonnes de V) sur une base orthonormée (colonnes de U) étirée par les valeurs singulières σ. Tronquer aux k premières donne la meilleure approximation de rang k (théorème d’Eckart-Young) : compression, débruitage, recommandation (Netflix), réduction de dimension avant apprentissage (ACP), pseudo-inverse pour les moindres carrés, et le « conditionnement » σ_max/σ_min qui mesure la sensibilité d’un système aux erreurs.

04 / Précision

Conditionnement et arithmétique flottante : quand le calcul ment

import numpy as np
# Matrice de Hilbert : célèbre pour être mal conditionnée
def hilbert(n): return np.array([[1 / (i + j + 1) for j in range(n)] for i in range(n)])
for n in (4, 8, 12):
    H = hilbert(n); x_vrai = np.ones(n); b = H @ x_vrai
    x = np.linalg.solve(H, b)
    print(f"n={n:2d}  conditionnement {np.linalg.cond(H):.1e}  erreur sur x : {np.abs(x - x_vrai).max():.1e}")
# Règle : on perd log10(cond) chiffres significatifs. Avec cond = 1e16 et 16 chiffres de précision : plus rien.

# Les flottants : 53 bits de mantisse ⇒ ~16 chiffres, et des surprises
print(0.1 + 0.2 == 0.3, np.isclose(0.1 + 0.2, 0.3))
print(1e16 + 1 - 1e16, "(1 est perdu : absorption)")
a = 1e8; print((a + 1e-8) - a, "(cancellation catastrophique)")
# Sommer des petits nombres : l'ordre compte
petits = np.full(10**6, 1e-10, dtype=np.float32)
print(np.sum(petits), petits.astype(np.float64).sum(), "(float32 : 24 bits de mantisse)")
# Kahan : somme compensée
def kahan(t):
    s, c = 0.0, 0.0
    for x in t:
        y = x - c; t2 = s + y; c = (t2 - s) - y; s = t2
    return s
import math
print(sum([0.1] * 10) == 1.0, kahan([0.1] * 10) == 1.0, math.fsum([0.1] * 10) == 1.0, "(somme naïve, Kahan, fsum)")

Sur un microcontrôleur sans unité flottante (Arduino Uno), chaque opération en float coûte des dizaines de cycles, et double n’existe pas (c’est un float 32 bits déguisé) : on travaille en virgule fixe (entiers avec facteur d’échelle). Module L18.

Cours

Cours 1 — Espaces vectoriels, bases, rang : le vocabulaire qui rend NumPy lisible

Notion	Définition opérationnelle	En NumPy
Combinaison linéaire	Σ λ_i v_i	`V @ lam` (colonnes de V pondérées)
Espace engendré (span)	Ensemble des combinaisons linéaires	« b ∈ span(colonnes de A) » ⇔ Ax = b a une solution
Indépendance linéaire	Aucune combinaison non triviale ne donne 0	`np.linalg.matrix_rank(V) == V.shape[1]`
Base, dimension	Famille libre et génératrice ; son cardinal	—
Rang	Dimension de l’espace des colonnes (= des lignes)	`matrix_rank` (via SVD : nombre de σ > ε)
Noyau	{x : Ax = 0}	Derniers vecteurs de Vᵀ dans la SVD (σ = 0)
Inversible	Carrée et rang plein ⇔ det ≠ 0 ⇔ noyau = {0}	Ne pas tester `det != 0` en flottant : regarder le conditionnement
Orthogonalité	u·v = 0 ; base orthonormée : Qᵀ Q = I	`np.linalg.qr`

Théorème du rang. Pour A de taille m×n : rang(A) + dim(noyau(A)) = n. Conséquence pour les systèmes : Ax = b a une solution unique si rang = n = m ; une infinité si rang < n (noyau non nul) ; aucune ou une infinité si m > n selon que b est dans l’espace des colonnes — d’où les moindres carrés quand il n’y en a pas.

Cours

Cours 2 — Exemple travaillé : ajuster un modèle de moteur par moindres carrés, avec incertitudes

import numpy as np
rng = np.random.default_rng(0)
# Un moteur : vitesse mesurée ω (rad/s) en fonction de la tension u (V) et du couple de charge c (N·m)
# Modèle physique linéaire : ω = k1·u − k2·c − k0  (k0 : frottement sec). On mesure 40 points bruités.
n = 40; u = rng.uniform(3, 12, n); c = rng.uniform(0, 0.5, n)
omega = 12.0 * u - 30.0 * c - 8.0 + rng.normal(0, 3.0, n)
X = np.column_stack([u, -c, -np.ones(n)])                             # colonnes : coefficients de k1, k2, k0
k, residus, rang, sv = np.linalg.lstsq(X, omega, rcond=None)
print("k1, k2, k0 estimés :", k.round(2), "(vrais : 12, 30, 8)")
# Incertitude sur les coefficients : covariance = σ² (XᵀX)⁻¹, σ² estimé par les résidus
sigma2 = ((omega - X @ k) ** 2).sum() / (n - 3)
cov = sigma2 * np.linalg.inv(X.T @ X)
print("écarts-types :", np.sqrt(np.diag(cov)).round(2), "| σ résiduel :", np.sqrt(sigma2).round(2), "(vrai 3.0)")
print("conditionnement de X :", np.linalg.cond(X).round(1), "(< 1000 : correct)")
# Prédiction avec intervalle : à u = 9 V, c = 0.2
x0 = np.array([9, -0.2, -1]); pred = x0 @ k; err = np.sqrt(x0 @ cov @ x0)
print(f"ω(9 V, 0.2 N·m) ≈ {pred:.1f} ± {2 * err:.1f} rad/s (95 %)")

Lecture. C’est l’identification de système (module L21 : le modèle du LQR vient de là). Le conditionnement dit si les colonnes sont « presque dépendantes » — si u et c étaient corrélés dans les mesures, k1 et k2 seraient mal séparés et leurs écarts-types énormes : il faut alors concevoir l’expérience (varier u et c indépendamment). La covariance (XᵀX)⁻¹ est aussi ce que manipule le filtre de Kalman (L19).

Cours

Cours 3 — Diagonalisation, puissances de matrices, systèmes dynamiques

Théorème. Si A (n×n) possède n vecteurs propres indépendants (colonnes de P) avec valeurs propres λ_i, alors A = P D P⁻¹ avec D = diag(λ_i), et A^k = P D^k P⁻¹. Donc l’évolution x_{k+1} = A x_k a pour solution x_k = Σ c_i λ_i^k v_i : chaque mode croît ou décroît comme λ_i^k.

|λ_i| < 1 pour tout i ⇒ x_k → 0 : stable (L21 : pôles).
Une valeur propre λ = 1 avec les autres < 1 ⇒ convergence vers le vecteur propre associé : état stationnaire (Markov, PageRank).
|λ| > 1 ⇒ divergence ; λ complexes ⇒ oscillations (rotation dans le plan propre).
Fibonacci : A = [[1,1],[1,0]], λ = φ et −1/φ ; F_n ≈ φⁿ/√5.

Toutes les matrices ne sont pas diagonalisables (bloc de Jordan [[1,1],[0,1]]), mais les symétriques le sont toujours, en base orthonormée (théorème spectral) — et la SVD diagonalise « au sens large » n’importe quelle matrice, même rectangulaire. En pratique : eigh pour les symétriques (rapide, stable), eig sinon, svd pour le rang et les moindres carrés.

TP guidé

TP — Identification d’un robot et compression d’images (sur PC, 2 h 30)

Données. Enregistrez (ou simulez avec bruit si vous n’avez pas de robot) 200 lignes pwm_gauche, pwm_droite, vitesse_lineaire, vitesse_angulaire à partir des encodeurs (séance 27). CSV → NumPy (np.loadtxt).
Identification. Modèle linéaire (v, ω) = M · (pwm_g, pwm_d) + biais. Estimez M et le biais par lstsq ; affichez les résidus (histogramme : gaussiens ? un biais systématique révèle une non-linéarité, par exemple la zone morte des moteurs). Ajoutez des termes quadratiques et comparez les résidus. Conditionnement ?
Inverser le modèle. Pour une consigne (v, ω) voulue, calculez les PWM par solve (M est 2×2). Vérifiez sur le robot : erreur en boucle ouverte. Ce modèle inverse est votre « cinématique inverse » calibrée (L21).
Compression par SVD. Chargez une image en niveaux de gris (PIL.Image.open(...).convert("L"), np.asarray). SVD ; reconstruction aux rangs 5, 20, 50, 100 ; erreur relative ‖A − A_k‖_F / ‖A‖_F et taux de compression ; figure 2×2. Trouvez le rang pour une erreur < 5 %.
ACP sur des données réelles. from sklearn.datasets import load_digits (1797 chiffres 8×8, 64 dimensions). Centrer, SVD, projeter sur 2 composantes, nuage coloré par chiffre. Variance expliquée cumulée : combien de composantes pour 90 % ? Reconstruisez un chiffre avec 10 composantes.
Livrable. Notebook ou script + figures + RESULTATS.md : matrice M avec incertitudes, choix de modèle justifié par les résidus, rang de compression, nombre de composantes ACP.

Exercices

Exercices auto-corrigés — matrices et transformations

Exercice 1 — Transformations homogènes 2D

Écrivez rotation(theta), translation(dx, dy), mise_a_echelle(sx, sy) (3×3 homogènes) et appliquer(M, points) (points : tableau n×2 → n×2). Puis rotation_autour(theta, cx, cy) : rotation autour d’un centre quelconque, par composition.

Correction

def rotation(t): c, s = np.cos(t), np.sin(t); return np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])
def translation(dx, dy): return np.array([[1, 0, dx], [0, 1, dy], [0, 0, 1.]])
def mise_a_echelle(sx, sy): return np.diag([sx, sy, 1.])
def appliquer(M, P): H = np.column_stack([P, np.ones(len(P))]); return (H @ M.T)[:, :2]
def rotation_autour(t, cx, cy): return translation(cx, cy) @ rotation(t) @ translation(-cx, -cy)

Exercice 2 — Gauss-Jordan et inverse

Implémentez inverser(A) par Gauss-Jordan avec pivot partiel sur la matrice augmentée [A | I] (sans np.linalg), qui lève ValueError si A est singulière. Vérifiez A·A⁻¹ = I.

Correction

def inverser(A):
    n = len(A); M = np.hstack([A.astype(float), np.eye(n)])
    for k in range(n):
        p = k + np.argmax(np.abs(M[k:, k]))
        if abs(M[p, k]) < 1e-12: raise ValueError("matrice singulière")
        M[[k, p]] = M[[p, k]]; M[k] /= M[k, k]
        for i in range(n):
            if i != k: M[i] -= M[i, k] * M[k]
    return M[:, n:]

Exercices

Exercices auto-corrigés — moindres carrés et valeurs propres

Exercice 3 — Ajuster un cercle par moindres carrés

Des points bruités sur un cercle (un obstacle rond vu par un lidar). L’équation x² + y² + a x + b y + c = 0 est linéaire en (a, b, c). Écrivez ajuster_cercle(points) → (cx, cy, r) avec lstsq.

Correction

def ajuster_cercle(P):
    x, y = P[:, 0], P[:, 1]
    A = np.column_stack([x, y, np.ones(len(P))]); b = -(x**2 + y**2)
    a, bb, c = np.linalg.lstsq(A, b, rcond=None)[0]
    cx, cy = -a / 2, -bb / 2; return cx, cy, np.sqrt(cx**2 + cy**2 - c)

Exercice 4 — Puissance itérée et état stationnaire

a) puissance_iteree(A, iters) renvoie (λ_max, vecteur propre normalisé). b) stationnaire(P) renvoie la distribution stationnaire d’une matrice de transition (lignes sommant à 1) par puissance itérée sur Pᵀ. c) Vérifiez sur la chaîne météo du module.

Correction

def puissance_iteree(A, iters=500):
    x = np.ones(len(A))
    for _ in range(iters): x = A @ x; x /= np.linalg.norm(x)
    return float(x @ A @ x), x
def stationnaire(P):
    pi = np.ones(len(P)) / len(P)
    for _ in range(1000): pi = pi @ P
    return pi / pi.sum()

05 / Défis

Défi ★ — Solveur et vérificateur

Consigne

1) Générez 200 systèmes aléatoires 5×5 et comparez votre gauss à np.linalg.solve (écart max). 2) Trouvez une matrice pour laquelle Gauss sans pivot échoue (division par zéro) alors que le système a une solution. 3) Ajustez une parabole y = a + bx + cx² par moindres carrés sur 30 points bruités et tracez.

Correction (extraits)

np.random.seed(0); ecart = 0
for _ in range(200):
    A = np.random.rand(5, 5) + np.eye(5); b = np.random.rand(5)
    ecart = max(ecart, np.abs(gauss(A, b) - np.linalg.solve(A, b)).max())
print(ecart)
A = np.array([[0., 1.], [1., 0.]]); b = np.array([1., 2.])     # pivot nul en (0,0) mais solution (2, 1)
x = np.linspace(-2, 2, 30); y = 1 + 0.5 * x - 0.8 * x**2 + np.random.normal(0, 0.2, 30)
X = np.column_stack([np.ones(30), x, x**2]); c = np.linalg.lstsq(X, y, rcond=None)[0]; print(c.round(2))

05 / Défis

Défi ★★ — Cinématique inverse par jacobienne

Consigne

Pour le bras à 2 segments, la position de la main p(a1, a2) est non linéaire. Calculez la jacobienne J = ∂p/∂(a1, a2) (2×2, par différences finies ou à la main), puis itérez la méthode de Newton : Δa = J⁻¹ (cible − p). Atteignez 5 cibles aléatoires accessibles en moins de 20 itérations chacune ; tracez la trajectoire de la main. Que se passe-t-il quand le bras est tendu (a2 = 0) ? Regardez le déterminant de J.

Correction

def jacobienne(a, h=1e-6):
    J = np.zeros((2, 2))
    for k in range(2):
        d = np.zeros(2); d[k] = h
        J[:, k] = (p(a + d) - p(a - d)) / (2 * h)
    return J

def newton(cible, a=np.array([0.5, 0.5]), tol=1e-6):
    for it in range(50):
        e = cible - p(a)
        if np.linalg.norm(e) < tol: return a, it
        a = a + np.linalg.solve(jacobienne(a), e)          # pas de inv !
    return a, 50

np.random.seed(3)
for _ in range(5):
    r = np.random.uniform(30, 130); t = np.random.uniform(0, np.pi)
    cible = r * np.array([np.cos(t), np.sin(t)])
    a, it = newton(cible); print(cible.round(1), "→ angles", np.degrees(a).round(1), "en", it, "itérations")
print("det J bras tendu :", np.linalg.det(jacobienne(np.array([0.3, 0.0]))).round(6))

Quand det J → 0 (bras tendu ou replié), la jacobienne est singulière : une singularité cinématique. Le bras ne peut pas bouger dans certaines directions et Newton diverge ; les contrôleurs industriels utilisent une pseudo-inverse amortie (Levenberg-Marquardt). Ce défi contient tout un chapitre de robotique (module L21).

05 / Défis

Défi ★★★ — Esprit prépa : itération QR et laplacien de graphe

Consigne

1) Implémentez l’algorithme QR pour les valeurs propres d’une matrice symétrique : A₀ = A ; Aₖ = Q·R ⇒ Aₖ₊₁ = R·Q (avec np.linalg.qr). Montrez que la diagonale converge vers les valeurs propres et comparez à eig. 2) Construisez le laplacien L = D − A d’un graphe à deux communautés faiblement reliées (module L05). Calculez le vecteur propre de la 2e plus petite valeur propre (vecteur de Fiedler) : son signe partitionne le graphe. Vérifiez. 3) Expliquez pourquoi le nombre de valeurs propres nulles de L est le nombre de composantes connexes.

Correction

def qr_valeurs_propres(A, iters=200):
    A = A.copy()
    for _ in range(iters):
        Q, R = np.linalg.qr(A); A = R @ Q
    return np.sort(np.diag(A))
S = np.random.rand(6, 6); S = S + S.T
print(qr_valeurs_propres(S).round(4)); print(np.sort(np.linalg.eigvalsh(S)).round(4))

# Deux communautés de 6 sommets, 2 arêtes entre elles
n = 12; A = np.zeros((n, n)); rng = np.random.default_rng(0)
for grp in (range(0, 6), range(6, 12)):
    for i in grp:
        for j in grp:
            if i < j and rng.random() < 0.7: A[i, j] = A[j, i] = 1
A[2, 8] = A[8, 2] = 1; A[5, 6] = A[6, 5] = 1
L = np.diag(A.sum(1)) - A
val, vec = np.linalg.eigh(L)
print("valeurs propres :", val.round(3))
fiedler = vec[:, 1]
print("partition :", (fiedler > 0).astype(int))

xᵀLx = Σ_{(i,j) arête} (x_i − x_j)² : L mesure la « rugosité » d’une fonction sur le graphe. Le vecteur constant donne 0 ; sur k composantes, k vecteurs indicateurs indépendants donnent 0, et réciproquement une valeur nulle impose x constant sur chaque composante. Le vecteur de Fiedler minimise la rugosité sous contrainte d’orthogonalité à la constante : il « coupe » le graphe là où il y a le moins d’arêtes. C’est le spectral clustering, utilisé en vision, en bio-informatique et en analyse de réseaux.

06 / Vérification

Pour résoudre A·x = b, on écrit :

Deux questions supplémentaires

1. Pourquoi une rotation a-t-elle R⁻¹ = Rᵀ ? Ses colonnes sont orthonormées : RᵀR = I.

2. Combien de nombres pour stocker une image 1000×1000 en rang 20 ? 20 × (1000 + 1000 + 1) ≈ 40 000 au lieu d’un million : 25 fois moins.

Référence

Les mots à retenir

Mot	Définition
Produit matriciel	C[i,j] = Σ A[i,k]·B[k,j] ; O(n³) ; non commutatif.
Orthogonale	Matrice dont l’inverse est la transposée (rotations, réflexions).
Coordonnées homogènes	(x, y, 1) : translations et rotations en une seule matrice.
Gauss / LU	Élimination avec pivot ; factorisation pour résoudre vite plusieurs b.
Moindres carrés	Minimiser ‖Xc − d‖² ; projection orthogonale ; `lstsq`.
Valeur propre	λ tel que A·v = λ·v.
SVD	A = UΣVᵀ ; meilleure approximation de rang k.
ACP	SVD des données centrées : axes de variance maximale.
Conditionnement	σ_max/σ_min : amplification des erreurs.
Jacobienne	Matrice des dérivées partielles ; linéarisation locale.
Laplacien	D − A ; ses vecteurs propres révèlent la structure d’un graphe.

Pour continuer

Vous calculez avec des matrices

Module suivant : probabilités et statistiques computationnelles — simuler, estimer, tester, inférer : les fondations de l’apprentissage automatique et des filtres de capteurs.

À faire chez soi

Regarder les 15 vidéos de « Essence of linear algebra » (3Blue1Brown) et refaire chaque exemple en NumPy.
Implémenter la factorisation QR par Gram-Schmidt et la comparer à np.linalg.qr sur une matrice mal conditionnée.
Faire les exercices 1-10 du chapitre « moindres carrés » de Strang.

← L09 Sommaire L11 : Probabilités et statistiques →