Calcul numérique avec NumPy · Algèbre linéaire

Cette annexe fait le pont entre le cours et la pratique de l'ingénieur: chaque notion des douze chapitres y est mise en regard de la fonction NumPy ou SciPy qui la calcule, avec un exemple court dont la sortie est reproduite en commentaire. Tous les résultats affichés ont été obtenus avec NumPy 2.0 et SciPy 1.13; ils sont arrondis par l'option d'affichage indiquée ci-dessous. L'objectif n'est pas de remplacer le calcul à la main, indispensable pour comprendre, mais de vous donner les bons réflexes lorsque la matrice a $1 0^{4}$ lignes plutôt que $3$ .

Pourquoi NumPy et SciPy?

Python seul ne sait pas multiplier deux matrices. La bibliothèque NumPy fournit le type ndarray (tableau multidimensionnel homogène) et les opérations vectorisées; son sous-module numpy.linalg expose les routines d'algèbre linéaire. SciPy complète NumPy avec scipy.linalg (décompositions supplémentaires, noyau, image, exponentielle de matrice), scipy.sparse (matrices creuses) et scipy.integrate (équations différentielles).

Sous le capot, ces deux bibliothèques appellent LAPACK et BLAS, les bibliothèques Fortran/C qui font tourner le calcul scientifique depuis les années 1990: élimination de Gauss avec pivot partiel, algorithme QR pour les valeurs propres, bidiagonalisation pour la SVD. Vous bénéficiez donc d'implémentations optimisées, testées sur des décennies, et compilées pour votre processeur. Le prix à payer est que tout se fait en virgule flottante (float64, environ 16 chiffres significatifs): il n'y a ni fractions exactes ni symboles. Pour du calcul exact ou symbolique, on se tourne vers SymPy (sympy.Matrix), mentionné au chapitre 1.

Installation, dans un terminal:

pip install numpy scipy

Conventions d'importation utilisées dans toute l'annexe:

import numpy as np
from numpy import linalg as la
import scipy.linalg as sla

Tableaux, matrices et produit

Un vecteur est un tableau à une dimension, une matrice un tableau à deux dimensions. On les crée à partir de listes Python; les entiers restent des entiers (int64), ce qui peut surprendre plus loin: écrivez 1.0 ou passez dtype=float dès que des divisions interviennent.

A = np.array([[1, 2], [3, 4]])
B = np.array([[1.0, 2.0], [3.0, 4.0]])
print(A.shape, A.dtype, B.dtype)   # (2, 2) int64 float64

Le produit matriciel s'écrit avec l'opérateur @ (ou np.matmul, ou A.dot(B)). Il fonctionne aussi entre une matrice et un vecteur: A @ x calcule $A x$ , et x @ A calcule $x^{T} A$ . NumPy ne distingue pas vecteur ligne et vecteur colonne pour un tableau à une dimension; c'est presque toujours une simplification bienvenue.

L'étoile n'est pas le produit matriciel

L'étoile * est le produit terme à terme (produit de Hadamard), pas le produit des matrices du chapitre 2:

print(A * A)   # [[ 1  4]
               #  [ 9 16]]
print(A @ A)   # [[ 7 10]
               #  [15 22]]

L'ancienne classe np.matrix, pour laquelle * désignait le produit matriciel, est obsolète: n'utilisez que des ndarray et @.

Comparer des flottants avec np.allclose

En virgule flottante, $0, 1 + 0, 2 \neq = 0, 3$ et $(2)^{2} \neq = 2$ . Ne testez jamais une égalité exacte entre résultats calculés; utilisez np.isclose (scalaires) ou np.allclose (tableaux), qui acceptent une tolérance relative de $1 0^{- 5}$ et absolue de $1 0^{- 8}$ par défaut.

print(0.1 + 0.2 == 0.3)                       # False
print(np.allclose(0.1 + 0.2, 0.3))            # True
print(np.isclose(np.sqrt(2) ** 2, 2))         # True

La même remarque vaut pour « $det A = 0$ », « $A x = b$ » ou « $λ$ est entière»: on écrit np.isclose(la.det(A), 0) et np.allclose(A @ x, b).

Enfin, pour un affichage lisible, on limite le nombre de décimales et on supprime la notation scientifique pour les petits nombres. Toutes les sorties ci-dessous utilisent ce réglage:

np.set_printoptions(precision=4, suppress=True)

Chapitre 1 – Systèmes d'équations linéaires

La fonction centrale est la.solve(A, b), qui résout $A x = b$ pour $A$ carrée inversible par élimination de Gauss avec pivot partiel (routine LAPACK gesv). Elle n'affiche pas les étapes; elle renvoie directement la solution.

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])
x = la.solve(A, b)
print(x)                         # [ 2.  3. -1.]
print(np.allclose(A @ x, b))     # True
print(la.matrix_rank(A))         # 3

Si la matrice est singulière, solve lève LinAlgError: Singular matrix. Le rang se calcule avec la.matrix_rank (par SVD, avec un seuil relatif; voir le chapitre 12). Le critère de compatibilité du chapitre 1, $rg A = rg (A ∣ b)$ , se programme avec np.column_stack:

M = np.array([[1.0, 2.0, 3.0], [2.0, 4.0, 6.0]])
c = np.array([1.0, 2.0])
print(la.matrix_rank(M), la.matrix_rank(np.column_stack([M, c])))   # 1 1  -> compatible
c2 = np.array([1.0, 3.0])
print(la.matrix_rank(M), la.matrix_rank(np.column_stack([M, c2])))  # 1 2  -> incompatible

NumPy ne fournit pas la forme échelonnée réduite: elle n'est pas stable numériquement (le choix «ce pivot est-il nul?» devient arbitraire en flottant). Pour l'obtenir exactement, avec des fractions, on utilise SymPy: sympy.Matrix(A).rref() renvoie la forme réduite et la liste des colonnes pivots. Pour un système non carré ou incompatible, la.lstsq(M, c, rcond=None) renvoie la solution de norme minimale au sens des moindres carrés (chapitre 7), avec le rang de $M$ en troisième position:

sol, res, rank, sv = la.lstsq(M, c, rcond=None)
print(sol, rank)                  # [0.0714 0.1429 0.2143] 1
print(np.allclose(M @ sol, c))    # True

Ici le système admet une infinité de solutions et lstsq choisit celle de plus petite norme, $\frac{1}{14} (1, 2, 3)$ ; il ne décrit pas l'ensemble des solutions, ce que fait la fonction null_space du chapitre 5.

Chapitre 2 – Calcul matriciel et matrices inversibles

Les objets de base: np.eye(n) pour $I_{n}$ , np.zeros((m, n)), np.diag(d) pour une matrice diagonale, A.T pour la transposée, np.trace(A). Le produit n'est pas commutatif, et NumPy le rappelle immédiatement:

A = np.array([[1.0, 2.0], [3.0, 4.0]])
B = np.array([[0.0, 1.0], [1.0, 0.0]])
print(A @ B)                                       # [[2. 1.]
                                                   #  [4. 3.]]
print(B @ A)                                       # [[3. 4.]
                                                   #  [1. 2.]]
print(np.allclose((A @ B).T, B.T @ A.T))           # True
print(la.inv(A))                                   # [[-2.   1. ]
                                                   #  [ 1.5 -0.5]]
print(la.matrix_power(A, 3))                       # [[ 37.  54.]
                                                   #  [ 81. 118.]]

la.matrix_power(A, k) calcule $A^{k}$ par élévations au carré successives, et accepte $k$ négatif pour $A^{- k} = (A^{- 1})^{k}$ . Une matrice élémentaire s'obtient en appliquant l'opération sur les lignes à np.eye(n); la multiplier à gauche reproduit l'opération $L_{3} \leftarrow L_{3} - 2 L_{1}$ :

E = np.eye(3)
E[2, :] = E[2, :] - 2 * E[0, :]
M = np.arange(1.0, 10.0).reshape(3, 3)      # lignes (1 2 3), (4 5 6), (7 8 9)
print(E @ M)                                # [[1. 2. 3.]
                                            #  [4. 5. 6.]
                                            #  [5. 4. 3.]]

Préférez solve à inv

Résoudre $A x = b$ par la.inv(A) @ b coûte trois fois plus d'opérations que la.solve(A, b) et est moins précis, car l'inverse est lui-même calculé avec des erreurs d'arrondi qui se propagent ensuite au produit. Sur la matrice de Hilbert $12 \times 12$ (très mal conditionnée), avec la solution exacte $x = (1, \dots, 1)$ :

H = sla.hilbert(12)
b = H @ np.ones(12)
print(la.norm(la.solve(H, b) - 1))      # 0.3567   (solve)
print(la.norm(la.inv(H) @ b - 1))       # 21.661   (inv)

N'appelez inv que lorsque vous avez réellement besoin des coefficients de $A^{- 1}$ ; pour résoudre plusieurs systèmes de même matrice, passez une matrice de seconds membres: la.solve(A, np.column_stack([b1, b2])).

Chapitre 3 – Déterminants

la.det(A) calcule le déterminant par décomposition LU: produit des pivots, avec le signe de la permutation. Le résultat est un flottant, même pour une matrice à coefficients entiers, et il porte les erreurs d'arrondi:

A = np.array([[2.0, 1.0, -1.0], [-3.0, -1.0, 2.0], [-2.0, 1.0, 2.0]])
print(la.det(A))                                                 # -0.9999999999999996
M = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])
print(la.det(M))                                                 # -9.516e-16
print(np.isclose(la.det(M), 0), la.matrix_rank(M))               # True 2

Le déterminant vaut $- 1$ dans le premier cas et $0$ dans le second: la lecture « $det M \approx 1 0^{- 15}$ donc $M$ est presque inversible» est fausse, $M$ est singulière. Pour tester l'inversibilité, préférez la.matrix_rank ou le conditionnement la.cond (chapitre 12). Un second piège est l'échelle: le déterminant est homogène de degré $n$ , donc $det (0, 1 I_{20}) = 1 0^{- 20}$ et $det (10 I_{20}) = 1 0^{20}$ , alors que ces deux matrices sont parfaitement inversibles. Pour de grandes matrices, la.slogdet(A) renvoie le signe et $ln ∣ det A ∣$ , ce qui évite les dépassements de capacité.

La règle de Cramer s'écrit en quelques lignes, mais elle n'est jamais utilisée en calcul numérique: $n + 1$ déterminants coûtent chacun autant qu'une élimination de Gauss complète, et la division par $det A$ amplifie les erreurs. Elle reste un outil théorique (formule explicite, dépendance continue en $b$ ).

b = np.array([8.0, -11.0, -3.0])
xc = np.empty(3)
for j in range(3):
    Aj = A.copy()
    Aj[:, j] = b                       # colonne j remplacée par b
    xc[j] = la.det(Aj) / la.det(A)
print(xc)                              # [ 2.  3. -1.]   (mais utilisez la.solve!)

Les propriétés du chapitre 3 se vérifient numériquement: la.det(A.T) vaut -1.0, et la.det(A @ A) vaut 1.0 à $1 0^{- 15}$ près, comme $det (A)^{2}$ .

Chapitres 4 et 5 – Espaces vectoriels, noyau, image et rang

Une famille de vecteurs est libre si et seulement si la matrice dont ils sont les colonnes a pour rang le nombre de vecteurs. On assemble la matrice avec np.column_stack:

v1 = np.array([1.0, 2.0, 3.0]); v2 = np.array([2.0, 4.0, 6.0]); v3 = np.array([0.0, 1.0, 1.0])
V = np.column_stack([v1, v2, v3])
print(la.matrix_rank(V))          # 2   -> famille liée (v2 = 2 v1)

Les coordonnées d'un vecteur $x$ dans une base $B$ dont les vecteurs forment les colonnes de $B$ sont la solution de $B x^{'} = x$ , donc la.solve(B, x). Pour le noyau et l'image, SciPy fournit deux fonctions qui renvoient des bases orthonormées (calculées par SVD): sla.null_space(A) pour $ker A$ et sla.orth(A) pour $Im A$ . Les colonnes obtenues ne sont pas les vecteurs «à coefficients simples» que donne la résolution à la main, mais elles engendrent le même sous-espace.

A = np.array([[1.0, 2.0, 1.0, 3.0], [2.0, 4.0, 0.0, 2.0], [3.0, 6.0, 1.0, 5.0]])
print(la.matrix_rank(A))          # 2
N = sla.null_space(A)
print(N.shape)                    # (4, 2)   -> dim ker A = 2
print(np.allclose(A @ N, 0))      # True
Im = sla.orth(A)
print(Im.shape)                   # (3, 2)   -> dim Im A = 2
print(N)                          # [[-0.2092 -0.874 ]
                                  #  [-0.1126  0.467 ]
                                  #  [-0.8688  0.1199]
                                  #  [ 0.4344 -0.06  ]]

Le théorème du rang se lit sur les dimensions: la.matrix_rank(A) + N.shape[1] vaut 4, le nombre de colonnes. L'espace des lignes est sla.orth(A.T), et sla.null_space(A.T) donne le noyau de la transposée, orthogonal à l'image (chapitre 7). Pour tester si un vecteur $b$ appartient à $Im A$ , comparez la.matrix_rank(A) et la.matrix_rank(np.column_stack([A, b])), comme au chapitre 1. Pour déterminer l'intersection de deux sous-espaces $Vect (U)$ et $Vect (W)$ , on résout $U a = W c$ , c'est-à-dire qu'on calcule sla.null_space(np.column_stack([U, -W])).

Chapitre 6 – Applications linéaires et changements de base

La matrice d'une application linéaire $f$ dans les bases canoniques a pour colonnes les images des vecteurs de base. Si f est une fonction Python, on l'évalue sur les colonnes de np.eye(n):

def f(x):
    return np.array([x[0] + 2 * x[1], 3 * x[0] - x[1], x[1]])

F = np.column_stack([f(e) for e in np.eye(2)])
print(F)                          # [[ 1.  2.]
                                  #  [ 3. -1.]
                                  #  [ 0.  1.]]

Une matrice de passage $P$ a pour colonnes les vecteurs de la nouvelle base exprimés dans l'ancienne. Les formules du chapitre 6, $x = P x^{'}$ et $A^{'} = P^{- 1} A P$ , se traduisent sans jamais former $P^{- 1}$ :

P = np.array([[1.0, 1.0], [1.0, -1.0]])      # nouvelle base: (1,1) et (1,-1)
x = np.array([3.0, 1.0])
x_new = la.solve(P, x)                       # coordonnées dans la nouvelle base
print(x_new)                                 # [2. 1.]
A = np.array([[2.0, 1.0], [1.0, 2.0]])
A_new = la.solve(P, A @ P)                   # P^{-1} A P
print(A_new)                                 # [[ 3.  0.]
                                             #  [-0.  1.]]

La matrice $A$ devient diagonale dans cette base: les colonnes de $P$ étaient des vecteurs propres (chapitre 9). Le -0. affiché est un zéro négatif de la virgule flottante, sans signification. La composition $g \circ f$ correspond au produit G @ F, dans cet ordre; pour une rotation d'angle $θ$ dans le plan, np.array([[np.cos(t), -np.sin(t)], [np.sin(t), np.cos(t)]]) avec t = np.pi / 6 envoie $(1, 0)$ sur [0.866 0.5].

Chapitre 7 – Produit scalaire, orthogonalité et projections

Produit scalaire, norme et angle:

u = np.array([1.0, 2.0, 2.0]); v = np.array([2.0, 0.0, 1.0])
print(u @ v, np.dot(u, v))                 # 4.0 4.0
print(la.norm(u), la.norm(v))              # 3.0 2.2361
cos_t = u @ v / (la.norm(u) * la.norm(v))
print(np.degrees(np.arccos(cos_t)))        # 53.3957
print(la.norm(u, 1), la.norm(u, np.inf))   # 5.0 2.0   (normes 1 et infini)

Avant np.arccos, encadrez le cosinus par np.clip(cos_t, -1, 1): pour deux vecteurs colinéaires, l'arrondi peut donner $1, 0000000000000002$ et arccos renverrait nan.

Le procédé de Gram–Schmidt est implémenté, sous une forme numériquement plus robuste (réflexions de Householder), par la décomposition QR: Q, R = la.qr(A) avec $Q$ à colonnes orthonormées engendrant le même espace que celles de $A$ , et $R$ triangulaire supérieure. Les signes des colonnes de $Q$ peuvent différer de ceux du calcul à la main.

A = np.array([[1.0, 1.0], [1.0, 0.0], [0.0, 1.0]])
Q, R = la.qr(A)
print(Q)                                   # [[-0.7071  0.4082]
                                           #  [-0.7071 -0.4082]
                                           #  [-0.      0.8165]]
print(np.allclose(Q.T @ Q, np.eye(2)))     # True
Pw = Q @ Q.T                               # matrice de la projection orthogonale sur Im A
b = np.array([1.0, 2.0, 3.0])
p = Pw @ b
print(p)                                   # [2.3333 0.6667 1.6667]
print(np.allclose(A.T @ (b - p), 0))       # True   (b - p orthogonal aux colonnes de A)

La formule de projection du cours, $p = A (A^{T} A)^{- 1} A^{T} b$ , s'écrit A @ la.solve(A.T @ A, A.T @ b) et donne le même vecteur; le complément orthogonal $W^{⊥}$ d'un sous-espace engendré par les colonnes de $W$ est sla.null_space(W.T).

Les moindres carrés sont le pain quotidien de l'ingénieur: ajuster une droite $y = a + b t$ sur des mesures revient à résoudre $X β \approx y$ avec $X$ de colonnes $(1, \dots, 1)$ et $(t_{i})$ . Trois écritures équivalentes:

t = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
y = np.array([1.1, 2.9, 5.2, 6.8, 9.1])
X = np.column_stack([np.ones_like(t), t])
coef, res, rank, sv = la.lstsq(X, y, rcond=None)
print(coef, res)                              # [1.04 1.99] [0.107]   (résidu = somme des carrés)
print(la.solve(X.T @ X, X.T @ y))             # [1.04 1.99]   (équations normales)
print(np.polyfit(t, y, 1))                    # [1.99 1.04]   (coefficients par degré décroissant!)

np.polyfit(t, y, d) ajuste un polynôme de degré $d$ et renvoie les coefficients du plus haut degré au plus bas; np.polyval(coef, t) l'évalue. Les équations normales sont acceptables pour un petit problème bien conditionné, mais lstsq (par SVD) est plus stable: le conditionnement de $X^{T} X$ est le carré de celui de $X$ .

Chapitre 8 – Géométrie vectorielle de l'espace

Le produit vectoriel est np.cross, et le produit mixte $[u, v, w] = det (u, v, w)$ se calcule indifféremment par np.cross(u, v) @ w ou par le déterminant de la matrice des trois colonnes.

u = np.array([1.0, 2.0, 0.0]); v = np.array([0.0, 1.0, 3.0]); w = np.array([2.0, 0.0, 1.0])
c = np.cross(u, v)
print(c, c @ u, c @ v)                          # [ 6. -3.  1.] 0.0 0.0
print(la.norm(c))                               # 6.7823   (aire du parallélogramme)
print(la.det(np.column_stack([u, v, w])))       # 13.0     (volume du parallélépipède)
print(abs(np.cross(u, v) @ w) / 6)              # 2.1667   (volume du tétraèdre)

Les formules de distance du chapitre 8 s'écrivent directement. Plan par trois points $A$ , $B$ , $C$ , puis distance du point $M = (1, 1, 1)$ à ce plan, distance à une droite, et distance entre deux droites gauches:

A = np.array([1.0, 0.0, 0.0]); B = np.array([0.0, 2.0, 0.0]); C = np.array([0.0, 0.0, 3.0])
n = np.cross(B - A, C - A); d = n @ A
print(n, d)                                       # [6. 3. 2.] 6.0   -> plan 6x + 3y + 2z = 6
M = np.array([1.0, 1.0, 1.0])
print(abs(n @ M - d) / la.norm(n))                # 0.7143   (distance point–plan)
P0 = np.zeros(3); dv = np.array([1.0, 1.0, 0.0])
print(la.norm(np.cross(M - P0, dv)) / la.norm(dv))       # 1.0   (distance point–droite)
P1 = np.zeros(3); d1 = np.array([1.0, 0.0, 0.0])
P2 = np.array([0.0, 1.0, 1.0]); d2 = np.array([0.0, 1.0, 0.0])
nn = np.cross(d1, d2)
print(abs((P2 - P1) @ nn) / la.norm(nn))          # 1.0   (distance entre droites gauches)

En mécanique, le moment d'une force $F$ appliquée en $r$ est np.cross(r, F): pour r = [0, 2, 0] et F = [0, 0, -10] on obtient [-20. 0. 0.], un moment de $20 N\cdotm$ autour de l'axe $x$ dans le sens négatif. Attention: np.cross n'est défini que pour des vecteurs de dimension $3$ (et, par extension, $2$ ); il n'existe pas en dimension supérieure.

Chapitre 9 – Valeurs propres et vecteurs propres

la.eig(A) renvoie le tableau des valeurs propres et la matrice dont les colonnes sont les vecteurs propres associés, normés à $1$ . la.eigvals(A) ne renvoie que les valeurs propres. L'algorithme sous-jacent (QR implicite, LAPACK geev) ne passe pas par le polynôme caractéristique, qui est numériquement très mal conditionné; np.poly(A) le fournit néanmoins, par degré décroissant, si vous voulez le comparer à votre calcul.

A = np.array([[4.0, 1.0], [2.0, 3.0]])
w, V = la.eig(A)
print(w)                                        # [5. 2.]
print(V)                                        # [[ 0.7071 -0.4472]
                                                #  [ 0.7071  0.8944]]
print(np.allclose(A @ V[:, 0], w[0] * V[:, 0])) # True
print(np.poly(A))                               # [ 1. -7. 10.]   -> lambda^2 - 7 lambda + 10
print(np.trace(A), w.sum(), la.det(A), w.prod())# 7.0 7.0 10.0 10.0

Ordre des valeurs propres et vecteurs propres en colonnes

Les valeurs propres de la.eig sont renvoyées dans un ordre quelconque, ni croissant ni décroissant: pour la matrice diagonale np.diag([2, 5, 1]) on obtient [2. 5. 1.]. Le $i$ -ème vecteur propre est la colonne V[:, i], pas la ligne V[i]; c'est l'erreur la plus fréquente. Pour trier par valeur propre décroissante:

idx = np.argsort(w)[::-1]
w, V = w[idx], V[:, idx]

Les vecteurs propres sont normés, donc un vecteur propre $(1, 1)$ apparaît comme [0.7071, 0.7071], éventuellement avec le signe opposé.

Une matrice réelle peut avoir des valeurs propres complexes; eig renvoie alors des tableaux de type complex128, et les vecteurs propres sont complexes aussi:

R = np.array([[0.0, -1.0], [1.0, 0.0]])       # rotation d'un quart de tour
w, V = la.eig(R)
print(w)                                      # [0.+1.j 0.-1.j]

Le sous-espace propre associé à $λ$ est le noyau de $A - λ I$ : sla.null_space(A - lam * np.eye(n)), dont le nombre de colonnes est la multiplicité géométrique. Lorsque la multiplicité géométrique est inférieure à la multiplicité algébrique, eig renvoie tout de même $n$ colonnes, mais elles sont liées:

J = np.array([[2.0, 1.0], [0.0, 2.0]])
w, V = la.eig(J)
print(w)                                      # [2. 2.]
print(la.matrix_rank(V))                      # 1   -> non diagonalisable

Retenez enfin qu'une valeur propre double est fragile: une perturbation de $1 0^{- 16}$ de $J$ peut la scinder en deux valeurs distantes de $1 0^{- 8}$ . Regroupez les valeurs propres avec np.isclose, pas avec ==.

Chapitre 10 – Diagonalisation et applications

La diagonalisation $A = P D P^{- 1}$ se vérifie et s'exploite avec np.diag. La matrice est diagonalisable si les vecteurs propres renvoyés sont indépendants, ce qu'on teste par le rang de $P$ :

A = np.array([[4.0, 1.0], [2.0, 3.0]])
w, P = la.eig(A)
D = np.diag(w)
print(np.allclose(P @ D @ la.inv(P), A))              # True
print(la.matrix_rank(P) == A.shape[0])                # True   (diagonalisable)
print(P @ np.diag(w ** 10) @ la.inv(P))               # [[6510758. 3254867.]
                                                      #  [6509734. 3255891.]]
print(la.matrix_power(A, 10))                         # même résultat

Pour les suites récurrentes, la.matrix_power suffit. Piège d'entiers: avec un tableau int64, $F^{93}$ pour la matrice de Fibonacci dépasse $2^{63}$ et NumPy renvoie des nombres négatifs sans avertissement; en float64 le résultat est approché mais de bon ordre de grandeur. L'exponentielle de matrice $e^{A t}$ du cours est sla.expm; ce n'est pas np.exp, qui agit terme à terme:

A = np.array([[0.0, 1.0], [-1.0, 0.0]])
print(sla.expm(A * np.pi / 2))     # [[ 0.  1.]
                                   #  [-1.  0.]]       (rotation d'angle pi/2)
print(sla.expm(np.diag([1.0, 2.0])))   # diag(e, e^2) = [[2.7183 0.], [0. 7.3891]]
print(np.exp(np.diag([1.0, 2.0])))     # [[2.7183 1.], [1. 7.3891]]   FAUX pour e^A

Chaîne de Markov. Avec la convention du cours (colonnes stochastiques, $x_{k + 1} = P x_{k}$ ), on itère par produits successifs et on obtient le vecteur stationnaire comme vecteur propre de $P$ pour $λ = 1$ , normalisé pour que la somme vaille $1$ :

P = np.array([[0.9, 0.2], [0.1, 0.8]])
print(P.sum(axis=0))                            # [1. 1.]   (colonnes stochastiques)
x = np.array([1.0, 0.0])
for k in range(3):
    x = P @ x
    print(x)                                    # [0.9 0.1], [0.83 0.17], [0.781 0.219]
print(la.matrix_power(P, 50) @ [1.0, 0.0])      # [0.6667 0.3333]
w, V = la.eig(P)
k = np.argmin(abs(w - 1))                       # indice de la valeur propre 1
pi = np.real(V[:, k]); pi = pi / pi.sum()
print(pi, np.allclose(P @ pi, pi))              # [0.6667 0.3333] True

Si votre matrice est à lignes stochastiques (convention fréquente en anglais, $π_{k + 1} = π_{k} P$ ), le vecteur stationnaire est un vecteur propre de P.T: appelez la.eig(P.T).

Système différentiel $x^{'} = A x$ . La solution exacte est $x (t) = e^{A t} x_{0}$ , soit sla.expm(A * t) @ x0, ou, via la diagonalisation, P @ (c * np.exp(w * t)) avec c = la.solve(P, x0). Pour un système non linéaire, ou pour obtenir la trajectoire sur une grille de temps, on utilise l'intégrateur de SciPy:

from scipy.integrate import solve_ivp
A = np.array([[0.0, 1.0], [-2.0, -3.0]])                # oscillateur amorti, valeurs propres -1 et -2
sol = solve_ivp(lambda t, x: A @ x, (0.0, 5.0), [1.0, 0.0],
                t_eval=[0.0, 1.0, 2.0, 5.0], rtol=1e-8)
print(sol.y.T)                                          # [[ 1.      0.    ]
                                                        #  [ 0.6004 -0.4651]
                                                        #  [ 0.2524 -0.234 ]
                                                        #  [ 0.0134 -0.0134]]
print(sla.expm(A * 1.0) @ [1.0, 0.0])                   # [ 0.6004 -0.4651]

La première composante en $t = 1$ vaut $2 e^{- 1} - e^{- 2} = 0, 6004$ , conformément à la résolution par valeurs propres du chapitre 10.

Chapitre 11 – Matrices symétriques et formes quadratiques

Pour une matrice symétrique, utilisez la.eigh plutôt que la.eig: l'algorithme est deux fois plus rapide, garantit des valeurs propres réelles renvoyées en ordre croissant, et des vecteurs propres orthonormés, c'est-à-dire directement la matrice orthogonale $Q$ du théorème spectral. la.eigvalsh ne renvoie que les valeurs propres.

S = np.array([[2.0, 1.0, 0.0], [1.0, 2.0, 0.0], [0.0, 0.0, 3.0]])
w, Q = la.eigh(S)
print(w)                                          # [1. 3. 3.]
print(Q)                                          # [[-0.7071  0.7071  0.    ]
                                                  #  [ 0.7071  0.7071  0.    ]
                                                  #  [ 0.      0.      1.    ]]
print(np.allclose(Q.T @ Q, np.eye(3)))            # True
print(np.allclose(Q @ np.diag(w) @ Q.T, S))       # True

eigh ne lit que le triangle inférieur (paramètre UPLO) et ne vérifie pas la symétrie: appliquée à une matrice non symétrique, elle renvoie un résultat faux sans avertissement. Vérifiez np.allclose(S, S.T) au préalable.

La forme quadratique $q (x) = x^{T} S x$ s'écrit x @ S @ x; sa nature (définie positive, négative, indéfinie) se lit sur la.eigvalsh(S). La factorisation de Cholesky $S = L L^{T}$ n'existe que pour les matrices symétriques définies positives, ce qui en fait le test le plus rapide (et celui utilisé en pratique, par exemple pour vérifier une matrice de rigidité):

x = np.array([1.0, -1.0, 2.0])
print(x @ S @ x)                                  # 14.0
Aq = np.array([[1.0, 2.0], [2.0, -2.0]])          # q = x^2 + 4xy - 2y^2
print(la.eigvalsh(Aq))                            # [-3.  2.]   -> forme indéfinie

def is_positive_definite(M):
    try:
        la.cholesky(M)
        return True
    except la.LinAlgError:                        # "Matrix is not positive definite"
        return False

print(is_positive_definite(S), is_positive_definite(Aq))   # True False
print(la.cholesky(S))                             # [[1.4142 0.     0.    ]
                                                  #  [0.7071 1.2247 0.    ]
                                                  #  [0.     0.     1.7321]]

Pour classer une conique $5 x^{2} - 4 x y + 2 y^{2} = 1$ , on diagonalise la matrice symétrique associée: la.eigh([[5, -2], [-2, 2]]) donne les valeurs propres [1. 6.], toutes deux positives, donc une ellipse de demi-axes $1$ et $1/ 6$ ; l'angle des axes est np.degrees(np.arctan2(Q[1, 0], Q[0, 0])). Si $det Q = - 1$ , échangez ou changez le signe d'une colonne pour obtenir une rotation.

Analyse en composantes principales. On centre les données (une ligne par observation), on forme la matrice de covariance avec np.cov(X, rowvar=False) et on la diagonalise avec eigh; comme les valeurs propres sortent en ordre croissant, on les retourne:

rng = np.random.default_rng(0)
t = rng.normal(size=200)
X = np.column_stack([2 * t, t]) + 0.1 * rng.normal(size=(200, 2))   # nuage allongé selon (2, 1)
Xc = X - X.mean(axis=0)
C = np.cov(Xc, rowvar=False)
w, V = la.eigh(C)
w, V = w[::-1], V[:, ::-1]                    # ordre décroissant
print(w)                                      # [4.6357 0.0103]
print(V[:, 0])                                # [-0.8909 -0.4542]   ~ -(2, 1)/sqrt(5)
print(w / w.sum())                            # [0.9978 0.0022]     (variance expliquée)
scores = Xc @ V[:, 0]                         # coordonnées sur le premier axe

Le premier axe principal retrouve la direction $(2, 1) / 5 \approx (0, 894, 0, 447)$ , au signe près, et explique $99, 8 %$ de la variance. Pour de grandes matrices de données, on préfère la SVD de Xc directement (chapitre 12), qui évite de former $X^{T} X$ .

Chapitre 12 – Décompositions matricielles et calcul numérique

LU. sla.lu(A) renvoie $P$ , $L$ , $U$ avec $A = P LU$ (attention: la permutation est à gauche, contrairement à la convention $P A = LU$ de certains cours). Pour résoudre plusieurs systèmes de même matrice, on factorise une fois avec sla.lu_factor et on résout avec sla.lu_solve:

A = np.array([[2.0, 1.0, 1.0], [4.0, 3.0, 3.0], [8.0, 7.0, 9.0]])
P, L, U = sla.lu(A)
print(L)                                      # [[1.     0.     0.    ]
                                              #  [0.25   1.     0.    ]
                                              #  [0.5    0.6667 1.    ]]
print(U)                                      # [[ 8.      7.      9.    ]
                                              #  [ 0.     -0.75   -1.25  ]
                                              #  [ 0.      0.     -0.6667]]
print(np.allclose(P @ L @ U, A))              # True
lu, piv = sla.lu_factor(A)
print(sla.lu_solve((lu, piv), [4.0, 10.0, 24.0]))    # [1. 1. 1.]

Le pivot partiel a placé la ligne $(8, 7, 9)$ en tête: c'est pourquoi $U$ ne ressemble pas à l'élimination «à la main» sans échange.

QR (la.qr, chapitre 7) et SVD: U, s, Vt = la.svd(A) renvoie les valeurs singulières s en ordre décroissant et la matrice $V^{T}$ , pas $V$ . Les valeurs singulières sont les racines carrées des valeurs propres de $A^{T} A$ , et elles donnent la norme spectrale, la norme de Frobenius, le rang et le conditionnement:

U, s, Vt = la.svd(A)
print(s)                                              # [15.2486  1.1963  0.2193]
print(np.allclose(U @ np.diag(s) @ Vt, A))            # True
print(np.allclose(np.sqrt(la.eigvalsh(A.T @ A))[::-1], s))   # True
print(la.norm(A, 2), la.norm(A, "fro"))               # 15.2486 15.2971
print(la.cond(A), s[0] / s[-1])                       # 69.539 69.539
print(la.cond(sla.hilbert(8)))                        # 15257576253.48   ~ 1.5e10

Le conditionnement $κ (A) = σ_{1} / σ_{n}$ borne l'amplification des erreurs relatives: avec la matrice de Hilbert $8 \times 8$ , une perturbation relative de $2, 4 \cdot 1 0^{- 11}$ sur $b$ produit une erreur relative de $1, 5 \cdot 1 0^{- 5}$ sur $x$ , soit six ordres de grandeur, comme le prévoit $κ \approx 1 0^{10}$ . Règle pratique: avec $κ = 1 0^{k}$ , on perd $k$ chiffres significatifs sur les $16$ disponibles. la.matrix_rank compte les valeurs singulières supérieures à s.max() * max(m, n) * eps; pour des données bruitées, fixez vous-même le seuil avec tol=.

La pseudo-inverse de Moore–Penrose la.pinv(A) est définie pour toute matrice, même rectangulaire ou de rang déficient; la.pinv(A) @ b est la solution des moindres carrés de norme minimale, identique à la.lstsq. Pour $A$ de rang plein en colonnes, pinv coïncide avec $(A^{T} A)^{- 1} A^{T}$ .

M = np.array([[1.0, 2.0], [2.0, 4.0], [3.0, 6.0]])   # rang 1
print(la.pinv(M) @ [1.0, 2.0, 3.0])                  # [0.2 0.4]
print(la.lstsq(M, [1.0, 2.0, 3.0], rcond=None)[0])   # [0.2 0.4]

Compression d'image par SVD tronquée. Une image en niveaux de gris est une matrice; conserver ses $k$ premières composantes singulières, $A_{k} = \sum_{i = 1}^{k} σ_{i} u_{i} v_{i}^{T}$ , donne la meilleure approximation de rang $k$ (théorème d'Eckart–Young), avec une erreur $∥ A - A_{k} ∥_{2} = σ_{k + 1}$ . Sur une image synthétique $64 \times 64$ (un disque, un coin clair et un dégradé):

n = 64
xs = np.linspace(-1.0, 1.0, n)
X, Y = np.meshgrid(xs, xs)
img = 1.0 * ((X - 0.2) ** 2 + (Y + 0.1) ** 2 < 0.3) + 0.5 * (X + Y > 0.9) + 0.3 * X
U, s, Vt = la.svd(img)
print(s[:6])                                  # [35.1046 10.7433  6.4351  3.9129  2.9341  2.5245]
for k in (1, 2, 5, 10, 20, 40):
    A_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
    err = la.norm(img - A_k) / la.norm(img)   # erreur relative (Frobenius)
    print(k, round(err, 4), k * (2 * n + 1))  # k, erreur, nombres stockés (sur 4096)
# 1 0.3889 129
# 2 0.2679 258
# 5 0.1637 645
# 10 0.1068 1290
# 20 0.0527 2580
# 40 0.012 5160
print(la.norm(img - U[:, :5] @ np.diag(s[:5]) @ Vt[:5, :], 2), s[5])   # 2.5245 2.5245

Avec $k = 10$ , on stocke $10 \times (64 + 64 + 1) = 1290$ nombres au lieu de $4096$ (un facteur $3$ ) pour une erreur relative de $11 %$ ; la ligne finale vérifie numériquement Eckart–Young. Sur une vraie photographie, chargée avec matplotlib.image.imread et convertie en gris, le spectre décroît plus lentement, mais le principe est le même: c'est aussi celui de l'ACP.

Matrices creuses. Les matrices issues de la discrétisation d'équations différentielles ou de réseaux (éléments finis, PageRank) ont une écrasante majorité de zéros. scipy.sparse ne stocke que les coefficients non nuls; scipy.sparse.linalg.spsolve et eigs remplacent la.solve et la.eig. La matrice tridiagonale du laplacien discret $1000 \times 1000$ tient en $2998$ coefficients au lieu de $1 0^{6}$ :

import scipy.sparse as sp
from scipy.sparse.linalg import spsolve
T = sp.diags([-1, 2, -1], [-1, 0, 1], shape=(1000, 1000), format="csr")
print(T.nnz)                                  # 2998
x = spsolve(T, np.ones(1000))
print(x[:3])                                  # [ 500.  999. 1497.]

Ne convertissez jamais une grande matrice creuse en tableau dense avec .toarray() «pour voir»: c'est précisément ce que le format évite.

Tableau récapitulatif

Notion du cours	NumPy / SciPy
Matrice, vecteur	`np.array([[…], […]])`, `np.array([…])`
Identité, zéros, diagonale	`np.eye(n)`, `np.zeros((m, n))`, `np.diag(d)`
Produit matriciel $A B$ , $A x$	`A @ B`, `A @ x`
Produit terme à terme	`A * B` (pas le produit matriciel)
Transposée, trace	`A.T`, `np.trace(A)`
Résolution de $A x = b$	`la.solve(A, b)`
Inverse $A^{- 1}$	`la.inv(A)` (préférer `solve`)
Puissance $A^{k}$	`la.matrix_power(A, k)`
Rang	`la.matrix_rank(A)`
Forme échelonnée réduite exacte	`sympy.Matrix(A).rref()`
Déterminant, $ln ∣ det A ∣$	`la.det(A)`, `la.slogdet(A)`
Test d'indépendance linéaire	`la.matrix_rank(np.column_stack([…]))`
Coordonnées dans une base $B$	`la.solve(B, x)`
Base de $ker A$	`sla.null_space(A)`
Base de $Im A$	`sla.orth(A)`
Complément orthogonal $W^{⊥}$	`sla.null_space(W.T)`
Matrice d'une application $f$	`np.column_stack([f(e) for e in np.eye(n)])`
Changement de base $P^{- 1} A P$	`la.solve(P, A @ P)`
Produit scalaire, norme	`u @ v`, `la.norm(u)`
Angle entre vecteurs	`np.arccos(np.clip(u @ v / (la.norm(u) * la.norm(v)), -1, 1))`
Gram–Schmidt, base orthonormée	`Q, R = la.qr(A)`
Projection orthogonale sur $Im A$	`Q @ Q.T @ b` ou `A @ la.solve(A.T @ A, A.T @ b)`
Moindres carrés	`la.lstsq(X, y, rcond=None)`, `np.polyfit(t, y, d)`
Produit vectoriel, produit mixte	`np.cross(u, v)`, `np.cross(u, v) @ w` ou `la.det(…)`
Valeurs et vecteurs propres	`w, V = la.eig(A)` (colonnes de `V`)
Valeurs propres seules	`la.eigvals(A)`
Polynôme caractéristique	`np.poly(A)` (coefficients par degré décroissant)
Sous-espace propre $E_{λ}$	`sla.null_space(A - lam * np.eye(n))`
Diagonalisation $P D P^{- 1}$	`P @ np.diag(w) @ la.inv(P)`
Exponentielle $e^{A t}$	`sla.expm(A * t)`
Vecteur stationnaire de Markov	`la.eig(P)`, colonne pour $λ = 1$ , normalisée
Système $x^{'} = A x$	`sla.expm(A * t) @ x0`, `scipy.integrate.solve_ivp`
Matrice symétrique (théorème spectral)	`w, Q = la.eigh(S)` (valeurs propres croissantes)
Forme quadratique	`x @ S @ x`, signature via `la.eigvalsh(S)`
Test définie positive	`la.cholesky(S)` (lève `LinAlgError` sinon)
Covariance, ACP	`np.cov(X, rowvar=False)` puis `la.eigh`
Décomposition LU	`sla.lu(A)`, `sla.lu_factor` / `sla.lu_solve`
Décomposition en valeurs singulières	`U, s, Vt = la.svd(A)`
Conditionnement	`la.cond(A)`
Pseudo-inverse	`la.pinv(A)`
Approximation de rang $k$	`U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]`
Matrices creuses	`scipy.sparse`, `scipy.sparse.linalg.spsolve`
Comparaison de flottants	`np.isclose`, `np.allclose`

Références

NumPy Developers, NumPy Reference, section «Linear algebra (numpy.linalg)», numpy.org/doc.
SciPy Developers, SciPy Reference Guide, sections scipy.linalg, scipy.sparse et scipy.integrate, docs.scipy.org.
L. N. Trefethen et D. Bau, Numerical Linear Algebra, SIAM, 1997 (25th anniversary edition, 2022): la référence sur la stabilité, le conditionnement et les algorithmes QR, SVD et de valeurs propres.
R. Johansson, Numerical Python: Scientific Computing and Data Science Applications with NumPy, SciPy and Matplotlib, 3e édition, Apress, 2024.
G. Strang, Linear Algebra and Learning from Data, Wellesley-Cambridge Press, 2019: SVD, moindres carrés et ACP du point de vue des applications.
E. Anderson et al., LAPACK Users' Guide, 3e édition, SIAM, 1999, netlib.org/lapack: la description des routines appelées par NumPy et SciPy.