Tutoriel SciPy — Calcul scientifique en Python

SciPy (Scientific Python) est une bibliothèque open source qui s'appuie sur NumPy et ajoute une large collection d'algorithmes pour les mathématiques, les sciences et l'ingénierie. Là où NumPy vous offre le ndarray et les opérations de base, SciPy vous fournit les solveurs spécialisés : intégration numérique, optimisation, interpolation, traitement du signal, statistiques, algorithmes spatiaux, et bien plus encore.

Ce chapitre couvre :

L'installation de SciPy et la convention d'importation
L'algèbre linéaire (scipy.linalg)
L'intégration et la différentiation numériques (scipy.integrate)
L'optimisation — trouver des minima et des racines (scipy.optimize)
L'interpolation (scipy.interpolate)
Les statistiques (scipy.stats)
Le traitement d'images à N dimensions (scipy.ndimage)

Installation de SciPy

SciPy est inclus dans la distribution Anaconda. Pour l'installer avec pip :

pip install scipy

SciPy dépend de NumPy, que pip installe automatiquement s'il n'est pas déjà présent.

Convention d'importation

SciPy est organisé en sous-packages. N'importez que les sous-packages dont vous avez besoin plutôt que la bibliothèque entière :

import numpy as np
from scipy import linalg, integrate, optimize, interpolate, stats, ndimage

Vous pouvez également vérifier la version installée :

import scipy
print(scipy.__version__)  # e.g. 1.13.0

Algèbre linéaire

scipy.linalg étend les routines d'algèbre linéaire de NumPy avec des décompositions et des solveurs supplémentaires. Toutes les fonctions s'appliquent aux tableaux NumPy classiques.

Déterminant et inverse

import numpy as np
from scipy import linalg

a = np.array([[1, 2],
              [3, 4]])

det = linalg.det(a)
print(det)          # -2.0

inv = linalg.inv(a)
print(inv)
# [[-2.   1. ]
#  [ 1.5 -0.5]]

det([[1,2],[3,4]]) = 1×4 − 2×3 = −2. L'inverse satisfait a @ inv == I.

Valeurs propres et vecteurs propres

Les valeurs propres décrivent comment une matrice étire l'espace ; les vecteurs propres donnent les directions qui ne tournent pas.

import numpy as np
from scipy import linalg

a = np.array([[1, 2],
              [3, 4]])

eigenvalues, eigenvectors = linalg.eig(a)
print(eigenvalues)
# [-0.37228132+0.j  5.37228132+0.j]
print(eigenvectors)
# [[-0.82456484 -0.41597356]
#  [ 0.56576746 -0.90937671]]

Chaque colonne de eigenvectors correspond à la valeur propre correspondante. Les valeurs propres sont retournées sous forme de nombres complexes même lorsque la partie imaginaire est nulle.

Décomposition en valeurs singulières (SVD)

La SVD factorise une matrice A en trois matrices U, s, Vt telles que A = U @ diag(s) @ Vt. C'est le fondement de l'analyse en composantes principales (ACP) et de nombreuses techniques de réduction de dimensionnalité.

import numpy as np
from scipy import linalg

a = np.array([[1, 2],
              [3, 4]])

u, s, vt = linalg.svd(a)
print(u)
# [[-0.40455358 -0.9145143 ]
#  [-0.9145143   0.40455358]]
print(s)    # [5.4649857  0.36596619]
print(vt)
# [[-0.57604844 -0.81741556]
#  [ 0.81741556 -0.57604844]]

Résolution d'un système linéaire

linalg.solve est la méthode correcte pour résoudre Ax = b. Elle est plus rapide et plus stable numériquement que de calculer l'inverse puis de multiplier.

import numpy as np
from scipy import linalg

# Solve: 1x + 2y = 5
#        3x + 4y = 11
A = np.array([[1, 2],
              [3, 4]])
b = np.array([5, 11])

x = linalg.solve(A, b)
print(x)   # [1. 2.]

# Verify: A @ x should equal b
print(np.allclose(A @ x, b))  # True

Intégration numérique

scipy.integrate fournit des routines pour calculer des intégrales définies lorsqu'une solution analytique est impraticable.

Intégration à une variable avec `quad`

integrate.quad utilise une quadrature adaptative pour intégrer une fonction sur un intervalle. Elle retourne le résultat et une estimation de l'erreur absolue.

import numpy as np
from scipy import integrate

# Integrate f(x) = x^2 + 2x + 1 from 0 to 1
# Analytical result: [x^3/3 + x^2 + x] from 0 to 1 = 1/3 + 1 + 1 = 7/3
def f(x):
    return x**2 + 2*x + 1

result, error = integrate.quad(f, 0, 1)
print(result)   # 2.3333333333333335
print(error)    # ~2.6e-14  (absolute error estimate)

Différentiation numérique

La fonction approx_fprime de SciPy calcule un gradient par différences finies. Pour les fonctions scalaires, la derivative par différences centrales de scipy.misc est plus simple :

import numpy as np
from scipy.optimize import approx_fprime

# Derivative of sin(x) at x = 0 should be cos(0) = 1
result = approx_fprime([0.0], lambda x: np.sin(x[0]), 1e-8)
print(result[0])   # ~1.0

Optimisation

scipy.optimize trouve des minima, des maxima (en minimisant la valeur négative) et des racines de fonctions.

Minimisation d'une fonction multivariée

optimize.minimize prend en charge de nombreuses méthodes (Nelder-Mead, BFGS, L-BFGS-B, …). La méthode par défaut est choisie automatiquement.

from scipy import optimize

# Minimize f(x) = x^2 + 2x + 1 = (x + 1)^2  — minimum at x = -1
def f(x):
    return x**2 + 2*x + 1

result = optimize.minimize(f, x0=0)   # x0 is the starting guess
print(result.success)   # True
print(result.x)         # [-1.00000001]  (near -1)
print(result.fun)       # ~0.0           (minimum value)

Minimisation d'une fonction scalaire sur un intervalle borné

optimize.minimize_scalar est plus simple pour les problèmes à une seule variable. Fournissez toujours bounds avec method='bounded' lorsque la fonction n'a pas de minimum global (c'est-à-dire qu'elle est non bornée) :

from scipy import optimize

# Minimize h(x) = x^2 - 4x + 3 over [0, 4]  — minimum at x = 2, h(2) = -1
def h(x):
    return x**2 - 4*x + 3

result = optimize.minimize_scalar(h, bounds=(0, 4), method='bounded')
print(result.x)     # ~2.0
print(result.fun)   # -1.0

Recherche de racines

optimize.root_scalar trouve où une fonction s'annule :

from scipy.optimize import root_scalar

# Solve x^2 - 4 = 0 in the interval [0, 3]  — root at x = 2
res = root_scalar(lambda x: x**2 - 4, bracket=[0, 3])
print(res.root)     # 2.0

Interpolation

scipy.interpolate ajuste une courbe lisse à travers des points de données afin que vous puissiez estimer les valeurs entre eux.

Interpolation à 1 dimension

interp1d crée une fonction d'interpolation appelable à partir de paires discrètes (x, y). Le paramètre kind sélectionne la méthode : 'linear' (par défaut), 'quadratic' ou 'cubic'.

import numpy as np
from scipy.interpolate import interp1d

# Sample points from y = x^2
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 4, 9, 16])

f_linear = interp1d(x, y)           # piecewise linear
f_cubic  = interp1d(x, y, kind='cubic')  # cubic spline

# Estimate y at x = 2.5  (exact value: 2.5^2 = 6.25)
print(float(f_linear(2.5)))   # 6.5   (linear — slightly off)
print(float(f_cubic(2.5)))    # 6.25  (cubic — matches exactly for polynomials)

L'interpolation cubique retrouve le résultat exact parce que les données proviennent d'un polynôme quadratique, et la spline cubique est suffisamment flexible pour l'ajuster parfaitement.

Statistiques

scipy.stats contient plus de 80 distributions de probabilités continues et discrètes ainsi qu'une collection de tests statistiques.

Statistiques descriptives et distribution normale

import numpy as np
from scipy.stats import norm

# Standard normal distribution (mean=0, std=1)
print(norm.pdf(0))     # 0.3989422804014327  — probability density at x = 0
print(norm.cdf(1.96))  # 0.9750021048517795  — P(X <= 1.96)
print(norm.ppf(0.975)) # 1.9599639845400536  — inverse CDF (quantile function)

# Fit a normal distribution to data
data = np.array([2.1, 3.3, 2.8, 3.1, 2.5, 3.0, 2.7])
mu, sigma = norm.fit(data)
print(f'Fitted mean: {mu:.4f}, std: {sigma:.4f}')

Tests d'hypothèses

scipy.stats inclut des tests t, des tests du chi-deux, l'ANOVA, les tests de Kolmogorov-Smirnov, et bien d'autres.

import numpy as np
from scipy.stats import ttest_1samp, ttest_ind

# One-sample t-test: is the sample mean significantly different from 5?
sample = np.array([4.8, 5.1, 4.9, 5.3, 5.2, 4.7, 5.0])
t_stat, p_value = ttest_1samp(sample, popmean=5.0)
print(f't = {t_stat:.4f}, p = {p_value:.4f}')
# p > 0.05 → no significant difference from 5

# Two-sample t-test: are these two groups different?
group_a = np.array([5.1, 5.3, 4.9, 5.2, 5.0])
group_b = np.array([6.1, 6.3, 5.8, 6.0, 5.9])
t_stat2, p_value2 = ttest_ind(group_a, group_b)
print(f't = {t_stat2:.4f}, p = {p_value2:.6f}')
# Very small p → groups are significantly different

Variables aléatoires et échantillonnage

Chaque distribution dans scipy.stats expose la même interface : pdf, cdf, ppf, rvs (variables aléatoires) et fit.

from scipy.stats import norm, poisson

# Draw 5 samples from a normal distribution with mean=10, std=2
samples = norm.rvs(loc=10, scale=2, size=5, random_state=42)
print(samples.round(2))   # [10.99  9.72 11.3  13.05  9.53]

# Poisson distribution: P(X = k) for mean lambda=3
for k in range(6):
    print(f'P(X={k}) = {poisson.pmf(k, mu=3):.4f}')

Traitement d'images à N dimensions

scipy.ndimage opère sur des tableaux de toute dimensionnalité (images, volumes, cubes de séries temporelles). Voici un exemple autonome qui utilise un tableau 2D synthétique afin qu'aucun fichier image externe ne soit nécessaire.

Flou gaussien et étiquetage des régions connexes

import numpy as np
from scipy import ndimage

# Create a synthetic 5x5 "image" with a bright spot in the centre
image = np.array([
    [0, 0, 0, 0, 0],
    [0, 1, 1, 1, 0],
    [0, 1, 5, 1, 0],
    [0, 1, 1, 1, 0],
    [0, 0, 0, 0, 0],
], dtype=float)

# Smooth the image with a Gaussian filter (sigma controls the blur radius)
blurred = ndimage.gaussian_filter(image, sigma=1)
print('Center value after blur:', round(blurred[2, 2], 4))
# 1.4166  — the bright peak is spread across neighbouring pixels

# Label connected non-zero regions (like counting distinct objects)
binary = image > 0
labeled, num_regions = ndimage.label(binary)
print('Number of connected regions:', num_regions)  # 1
print(labeled)
# [[0 0 0 0 0]
#  [0 1 1 1 0]
#  [0 1 1 1 0]
#  [0 1 1 1 0]
#  [0 0 0 0 0]]

Opérations courantes de `ndimage`

Fonction	Ce qu'elle fait
`gaussian_filter(a, sigma)`	Lissage avec un noyau gaussien
`sobel(a)`	Détection de contours (gradient de Sobel)
`label(a)`	Étiquetage des régions connexes dans un tableau binaire
`binary_dilation(a)`	Expansion des régions de premier plan
`zoom(a, factor)`	Redimensionnement d'un tableau
`rotate(a, angle)`	Rotation d'un tableau (en degrés)

Quand utiliser SciPy plutôt qu'une autre bibliothèque

Tâche	Outil recommandé
Création de tableaux et calculs de base	NumPy
DataFrames, séries temporelles, E/S	Pandas
Algorithmes scientifiques (intégration, optimisation)	SciPy
Apprentissage automatique	scikit-learn (s'appuie sur SciPy)
Visualisation	Matplotlib

SciPy ne remplace pas NumPy — il en dépend et l'étend. En pratique, vous importerez les deux.

Référence rapide

Sous-package	Fonctions clés
`scipy.linalg`	`det`, `inv`, `eig`, `svd`, `solve`
`scipy.integrate`	`quad`, `dblquad`, `solve_ivp`
`scipy.optimize`	`minimize`, `minimize_scalar`, `root_scalar`
`scipy.interpolate`	`interp1d`, `CubicSpline`, `griddata`
`scipy.stats`	`norm`, `ttest_1samp`, `ttest_ind`, `chi2_contingency`
`scipy.ndimage`	`gaussian_filter`, `label`, `sobel`, `zoom`
`scipy.signal`	`butter`, `lfilter`, `find_peaks`
`scipy.spatial`	`distance`, `KDTree`, `ConvexHull`