Nombres flottants et calcul stable

Pourquoi 0.1 + 0.2 diffère-t-il de 0.3 ? Pourquoi une autre écriture de la même formule mathématique peut-elle donner un résultat plus précis ? Il faut examiner la représentation des nombres, l’arrondi de chaque opération et la sensibilité du problème aux données. La vue d’ensemble de l’analyse numérique présente les principales catégories d’erreur ; cette note les développe à partir de calculs concrets.

Signe, exposant et significande

L’article de Goldberg sur l’arithmétique flottante, reproduit avec autorisation par Oracle décrit un flottant par son signe, sa significande et une puissance de la base. Pour un nombre normal au format IEEE 754 binary64,

$$x=(-1)^s(1+f)2^e,\qquad 0\le f<1.$$

Sur les 64 bits, un contient le signe, 11 l’exposant et 52 la partie fractionnaire. Le premier chiffre binaire d’un nombre normal vaut implicitement 1, ce qui donne 53 bits de précision. Si le champ d’exposant contient l’entier $E$, alors $e=E-1023$. Les valeurs $E=0$ et $E=2047$ sont réservées aux zéros, aux nombres sous-normaux, aux infinis et aux NaN : la formule des nombres normaux ne s’y applique pas.

Cette précision ne fixe pas un nombre de chiffres après la virgule. Dans l’intervalle $[2^e,2^{e+1})$, l’écart entre deux nombres normaux positifs consécutifs est

$$\Delta=2^{e-52}.$$

Quand l’ordre de grandeur double, l’espacement double aussi. L’écart vers le haut à partir de 1 est $2^{-52}$ ; à partir de 2, il vaut $2^{-51}$. Près de $10^{16}$, il atteint déjà 2 : ajouter 1 peut donc perdre cet incrément. Les formats flottants plus petits modifient la précision et l’étendue des valeurs. Les approximations de poids sur peu de bits décrites dans la quantification des modèles font aussi intervenir des échelles et des codages ; le nombre total de bits ne suffit pas à déterminer la précision du calcul.

Affichage décimal et valeur stockée

Le tutoriel Python sur les flottants explique pourquoi le développement binaire de 0.1 est périodique et infini, et pourquoi Python affiche généralement la plus courte représentation décimale qui permet de retrouver la même valeur stockée. Afficher 0.1 ne signifie pas stocker exactement $1/10$. Afficher davantage de chiffres ne modifie pas non plus la valeur.

Les exemples utilisent la bibliothèque standard de CPython 3.9 ou ultérieur (math.ulp est disponible depuis la version 3.9), avec un float binary64 et l’arrondi au plus proche, les égalités étant départagées par le chiffre pair. La première assertion vérifie la base binaire et la précision. Chaque bloc de code s’exécute indépendamment.

import math
import sys
from fractions import Fraction

assert sys.float_info.radix == 2 and sys.float_info.mant_dig == 53
x = 0.1
print(repr(x))
print(format(x, ".17g"))
print(x.as_integer_ratio())
print(0.1 + 0.2, (0.1 + 0.2) == 0.3)
error = Fraction.from_float(x) - Fraction(1, 10)
print(error)
print(float(error), float(error / Fraction(1, 10)))
print(error / Fraction.from_float(math.ulp(x)))
print(math.ulp(1.0), math.ulp(2.0), math.ulp(1e16))
print(1e16 + 1.0 == 1e16)

0.1
0.10000000000000001
(3602879701896397, 36028797018963968)
0.30000000000000004 False
1/180143985094819840
5.551115123125783e-18 5.551115123125783e-17
2/5
2.220446049250313e-16 4.440892098500626e-16 2.0
True

as_integer_ratio() donne la fraction exacte stockée. En lui soustrayant $1/10$, on obtient $1/180143985094819840$, soit environ $5.55\times10^{-18}$. Les opérations peuvent ajouter d’autres arrondis : 0.1 + 0.2 se trouve un pas représentable au-dessus de la valeur stockée pour 0.3.

Pour conserver le sens d’une entrée décimale, construire un Decimal à partir d’une chaîne ; pour des rationnels exacts, utiliser Fraction. Passer d’abord par float, puis convertir vers l’une de ces représentations, conserve la valeur déjà arrondie et ne restitue pas l’entrée d’origine.

Erreur absolue, erreur relative et ulp

Pour une valeur vraie $y$ et une valeur calculée $\hat y$,

$$E_{\mathrm{abs}}=|\hat y-y|,\qquad E_{\mathrm{rel}}=\frac{|\hat y-y|}{|y|}\quad(y\ne0).$$

L’erreur absolue a la même unité que le résultat ; l’erreur relative mesure la proportion perdue. Elle n’est pas définie quand la valeur vraie est nulle. Près de zéro, une petite erreur absolue peut correspondre à une grande erreur relative.

Une ulp, ou unité du dernier chiffre de la significande, permet de mesurer l’erreur à l’échelle de l’espacement des flottants. L’erreur de représentation de 0.1 ci-dessus vaut exactement $2/5$ de math.ulp(0.1). Pour un nombre positif fini, math.ulp donne l’écart vers le voisin supérieur, avec une règle particulière pour le plus grand nombre fini. Aux puissances de deux, les écarts vers le bas et vers le haut peuvent différer : il faut donc préciser la référence utilisée.

Dans le domaine des nombres normaux, un arrondi correct au plus proche introduit au plus une demi-ulp locale d’erreur. Pour un résultat normal non nul d’une opération élémentaire, sans débordement ni sous-dépassement, le modèle usuel est

$$\operatorname{fl}(a\circ b)=(a\circ b)(1+\delta), \qquad |\delta|\le u=2^{-53},$$

où $\circ$ désigne une opération arithmétique élémentaire. L’unité d’arrondi $u$ vaut la moitié de sys.float_info.epsilon, l’écart $2^{-52}$ vers le voisin supérieur de 1. Une borne pour une opération ne garantit pas que l’erreur relative de tout le programme reste inférieure à $u$.

Arrondi, débordement, sous-dépassement et valeurs spéciales

L’arrondi au plus proche choisit la valeur représentable la moins éloignée du résultat exact. À mi-chemin, il choisit celle dont le dernier chiffre de la significande est pair. 1e16 + 1 est précisément un tel cas : la trace précédente montre un retour à 1e16.

Situation	Sens en binary64	Conséquence pour le calcul
Débordement	Le résultat dépasse l’étendue finie ; la valeur finie maximale vaut environ $1.798\times10^{308}$	Avec l’arrondi au plus proche, un résultat suffisamment grand devient un infini signé ; l’interface du langage peut aussi lever une exception
Nombre sous-normal	Valeur vérifiant $0<\lvert x\rvert<2^{-1022}$ (valeur absolue inférieure au plus petit nombre normal positif), avec un premier chiffre implicite égal à 0	Le sous-dépassement progressif conserve de petites valeurs avec un espacement fixe $2^{-1074}$, mais leur précision relative diminue
Arrondi à zéro	Valeur trop petite pour rester un sous-normal non nul	La moitié du plus petit sous-normal positif s’arrondit à zéro avec la règle du chiffre pair
Infini	inf et -inf sont des valeurs spéciales	Ils peuvent représenter un dépassement de capacité ou certains résultats non bornés ; inf - inf produit NaN
NaN	« Not a number », pouvant résulter d’une opération invalide	N’est égal à aucune valeur, pas même à lui-même ; le tester avec math.isnan, et vérifier la finitude avec math.isfinite

Un nombre sous-normal n’est pas nécessairement inexact. Zéro possède deux codages, positif et négatif. L’article historique de Goldberg explique le sous-dépassement progressif et les valeurs spéciales ; les interfaces des langages ont leurs propres règles. Les conventions d’exception de math en Python prévoient des exceptions pour certaines opérations invalides et certains débordements, plutôt qu’un retour direct des valeurs spéciales IEEE :

import math
import sys

print(sys.float_info.min)
tiny = math.ulp(0.0)
print(tiny, tiny / 2)
print(1e308 * 1e308)
nan = float("inf") - float("inf")
print(nan, nan == nan, math.isnan(nan))
for operation in (
    lambda: 1.0 / 0.0,
    lambda: math.sqrt(-1.0),
    lambda: math.exp(1000.0),
):
    try:
        operation()
    except (ZeroDivisionError, ValueError, OverflowError) as exc:
        print(type(exc).__name__)

2.2250738585072014e-308
5e-324 0.0
inf
nan False True
ZeroDivisionError
ValueError
OverflowError

5e-324 est l’affichage décimal court du plus petit sous-normal positif, pas sa valeur décimale exacte. Mettre les quantités intermédiaires à une échelle appropriée est souvent plus efficace que tenter de récupérer l’information après l’apparition de zéro ou d’un infini. Un résultat final fini ne prouve pas à lui seul sa précision.

Annulation : réécrire une formule équivalente

Soustraire deux valeurs proches peut faire dominer leurs erreurs préexistantes dans le résultat. La soustraction elle-même peut même être exacte ; l’information a pu disparaître pendant le calcul des deux opérandes. La section de Goldberg sur l’annulation distingue ce cas de la soustraction d’opérandes exacts.

Considérons $g(x)=\sqrt{1+x}-1$ près de zéro. Pour un réel $x\ge-1$, multiplier par le conjugué donne

$$g(x)=\frac{(\sqrt{1+x}-1)(\sqrt{1+x}+1)}{\sqrt{1+x}+1} =\frac{x}{\sqrt{1+x}+1}.$$

Avec $x=10^{-16}$, la formule directe arrondit d’abord 1 + x à 1, puis soustrait 1 et renvoie zéro. La formule réécrite conserve la petite quantité $x$ au numérateur. Le dénominateur reste proche de 2 : même arrondi à 2, il permet de garder un résultat d’environ $5\times10^{-17}$.

import math
from decimal import Decimal, localcontext

x = 1e-16
naive = math.sqrt(1.0 + x) - 1.0
stable = x / (math.sqrt(1.0 + x) + 1.0)
with localcontext() as ctx:
    ctx.prec = 80
    dx = Decimal.from_float(x)
    reference = dx / ((Decimal(1) + dx).sqrt() + Decimal(1))
    print(f"reference = {reference:.18e}")
    for name, value in [("naive", naive), ("stable", stable)]:
        absolute = abs(Decimal.from_float(value) - reference)
        if reference != 0:
            relative = absolute / abs(reference)
            print(f"{name} = {value:.18e}; relative error = {relative:.3e}")
        else:
            print(f"{name} = {value:.18e}; absolute error = {absolute:.3e}")

reference = 4.999999999999999770e-17
naive = 0.000000000000000000e+00; relative error = 1.000e+0
stable = 4.999999999999999895e-17; relative error = 2.500e-17

La référence est une approximation de haute précision, calculée avec 80 chiffres décimaux. Decimal.from_float(x) conserve la même entrée binaire ; la comparaison mesure donc l’erreur algorithmique sans y mêler l’écart entre entrée décimale et entrée binaire. La formule directe a une erreur relative de 1 ; celle de la formule réécrite est d’environ $2.5\times10^{-17}$. Indépendamment, le développement $g(x)=x/2-x^2/8+O(x^3)$ donne environ $5\times10^{-17}$, avec un deuxième terme d’environ $-1.25\times10^{-33}$.

Si x vaut 0, le résultat exact est lui aussi nul : le code indique alors l’erreur absolue. La formule réécrite évite l’annulation près de zéro, mais pas le sous-dépassement du résultat lui-même. Pour la plus petite entrée sous-normale positive, le résultat vaut environ la moitié de cette entrée et s’arrondit encore à zéro en binary64.

Des fonctions spécialisées suivent le même principe : log1p(x) calcule $\log(1+x)$ avec précision près de zéro ; expm1(x) évite l’annulation du calcul direct de $e^x-1$. Leurs domaines de définition réels restent à respecter. Ces réécritures modifient aussi les valeurs intermédiaires parcourues par la différentiation automatique, qui ne supprime pas elle-même les erreurs flottantes.

Ordre de sommation et sommation compensée

L’addition réelle est associative ; l’addition flottante ne l’est généralement pas. Prenons $10^{16}$, dix valeurs égales à 1 et $-10^{16}$. La somme exacte vaut 10. L’addition séquentielle de gauche à droite peut faire disparaître chaque 1 dans la grande somme partielle. Annuler d’abord les deux grandes valeurs conserve ces contributions.

L’analyse des erreurs de sommation de Goldberg explique la sommation compensée de Kahan : à côté de la somme partielle, elle garde l’information de poids faible perdue dans une addition pour corriger la suivante. La boucle séquentielle est écrite explicitement ci-dessous afin de ne pas dépendre de l’implémentation de la fonction intégrée sum dans une version donnée de Python.

import math

def sequential(values):
    total = 0.0
    for value in values:
        total += value
    return total

def kahan(values):
    total = correction = 0.0
    for value in values:
        adjusted = value - correction
        new_total = total + adjusted
        correction = (new_total - total) - adjusted
        total = new_total
    return total

values = [1e16] + [1.0] * 10 + [-1e16]
print(sequential(values))
print(sequential([1e16, -1e16] + [1.0] * 10))
print(kahan(values), math.fsum(values))
hard = [1e16, 1.0, -1e16]
print(kahan(hard), math.fsum(hard))

0.0
10.0
10.0 10.0
0.0 1.0

correction = (new_total - total) - adjusted estime l’écart entre cette addition et l’incrément voulu ; l’étape suivante soustrait cette correction. Elle récupère les contributions perdues dans l’exemple des dix valeurs égales à 1. La dernière ligne montre pourtant que le Kahan ordinaire perd encore 1 pour [1e16, 1.0, -1e16]. La compensation ne constitue pas une arithmétique exacte.

La fonction Python math.fsum conserve plusieurs sommes partielles intermédiaires et donne la somme exacte dans les deux exemples. Sa précision dépend des garanties de l’arithmétique IEEE et du mode habituel d’arrondi au plus proche avec égalités départagées par le chiffre pair. La documentation précise aussi que, sur certaines versions compilées hors Windows, une addition en précision étendue peut provoquer un double arrondi intermédiaire et affecter le dernier bit.

Pour $n\ge1$ entrées flottantes finies, considérons leurs valeurs stockées comme exactes et posons $S=\sum_i x_i$. La borne suivante pour la sommation séquentielle ne compte que les arrondis des additions, sans inclure l’erreur de conversion des entrées. On suppose qu’il n’y a ni débordement ni sous-dépassement, que les résultats non nuls respectent le modèle d’arrondi élémentaire précédent et que les additions exactement nulles ont $\delta=0$. Chaque entrée peut subir jusqu’à $n-1$ arrondis ; en développant leurs facteurs multiplicatifs, on obtient la borne d’erreur absolue

$$|\hat S-S|\le\gamma_{n-1}\sum_{i=1}^n|x_i|, \qquad \gamma_k=\frac{ku}{1-ku}\quad(ku<1).$$

Quand $ku$ est petit, $\gamma_k$ est approximativement égal à $ku$. Si $S\ne0$, la borne relative comporte aussi le facteur $\sum_i|x_i|/|S|$, qui peut devenir grand lorsque les termes positifs et négatifs s’annulent presque.

Pour de grands jeux de données, la sommation par paires équilibrée organise les additions en arbre et réduit la chaîne la plus longue de $n-1$ à $\lceil\log_2 n\rceil$. Sous les mêmes conditions d’arrondi, on remplace alors $\gamma_{n-1}$ dans la borne par le $\gamma_k$ correspondant à cette profondeur. Additionner d’abord les petites valeurs réduit souvent leur risque de disparaître dans une grande somme partielle. Avec des signes mélangés, aucun ordre simple ne garantit toutefois le meilleur résultat. Modifier l’ordre d’une réduction parallèle peut également modifier les derniers bits. Choisir un ordre, une compensation ou une précision supérieure dépend des ordres de grandeur et de l’erreur visée. La note sur NumPy présente les réductions et les comparaisons de tableaux.

Conditionnement, stabilité rétrograde et tolérances

Le conditionnement mesure la sensibilité du problème mathématique aux perturbations des entrées. Pour $y=a-b\ne0$, si la perturbation relative de chaque entrée est au plus $\varepsilon$, l’inégalité triangulaire donne

$$|\Delta y|\le\varepsilon(|a|+|b|),\qquad \frac{|\Delta y|}{|y|}\le\kappa\varepsilon, \qquad \kappa=\frac{|a|+|b|}{|a-b|}.$$

En considérant $a=1.000001$ et $b=1.000000$ comme des entrées réelles exactes, la différence vaut $10^{-6}$ et $\kappa=2000001$. Une incertitude relative de $10^{-8}$ sur les entrées donne une borne d’incertitude absolue de $2.000001\times10^{-8}$ sur la sortie, soit environ 2 % en relatif. Même un calcul exact ne peut éliminer cette incertitude des données.

La stabilité rétrograde examine si le résultat calculé est la solution exacte d’un problème aux entrées voisines. Si le calcul de $f(x)$ donne $\hat y=f(x+\Delta x)$ et si la perturbation des entrées est petite dans une norme ou à une échelle par composante précisée, l’erreur rétrograde est petite. Un algorithme rétrogradement stable doit fournir des bornes de perturbation proportionnées à la précision d’arrondi sur son domaine d’application ; trouver une entrée voisine quelconque pour un seul résultat ne suffit pas. L’erreur directe compare, elle, $\hat y$ à $f(x)$. Pour de petites perturbations, avec $f$ différentiable et $f(x)\ne0$, la relation au premier ordre est approximativement « erreur directe relative $\lesssim$ nombre de conditionnement × erreur rétrograde relative ». Un problème mal conditionné peut donc présenter une erreur directe importante même avec un algorithme rétrogradement stable.

L’exemple de la racine carrée a une autre cause. Pour $x>0$, son nombre de conditionnement relatif est

$$\left|\frac{xg'(x)}{g(x)}\right| =\frac{\sqrt{1+x}+1}{2\sqrt{1+x}} \longrightarrow1\quad(x\to0^+).$$

Le problème est bien conditionné ici, mais la formule directe perd le résultat. Cette distinction permet de décider s’il faut améliorer les données ou réécrire l’algorithme.

Les tolérances doivent découler d’un budget d’erreur : incertitude des entrées, erreur de discrétisation ou d’approximation, et arrondis de l’algorithme. La tolérance absolue porte l’unité du résultat ; la tolérance relative s’appuie sur une échelle non nulle pertinente. La précision machine décrit la représentation et ne remplace pas ce budget.

Pour des valeurs finies, math.isclose utilise

$$|a-b|\le\max\bigl(\text{rel\_tol}\max(|a|,|b|),\text{abs\_tol}\bigr).$$

import math

print(math.isclose(0.1 + 0.2, 0.3, rel_tol=0.0, abs_tol=math.ulp(0.3)))
print(math.isclose(1e-12, 0.0, rel_tol=1e-9))
print(math.isclose(1e-12, 0.0, rel_tol=1e-9, abs_tol=1e-11))

True
False
True

La première ligne accepte un écart d’une ulp prise à la valeur stockée de 0.3, ce qui convient à ce court calcul déjà analysé. Les deux suivantes illustrent la nécessité d’une tolérance absolue pertinente près de zéro. 1e-11 sert de seuil de démonstration, pas de valeur par défaut universelle. NaN n’est proche d’aucune valeur ; un infini n’est proche que de l’infini de même signe. Pour les tableaux, choisir les tolérances selon la définition de l’interface NumPy plutôt que transposer la formule de math.isclose à une autre API.

En pratique, déterminer les échelles des entrées et de la sortie, repérer les valeurs intermédiaires non finies et les annulations, choisir une sommation fiable ou des fonctions spécialisées, puis vérifier l’erreur à l’aide d’une solution exacte, d’une référence de haute précision ou d’une borne démontrée. Deux algorithmes qui donnent des résultats proches sont cohérents entre eux ; leur précision demande une justification indépendante.

Mes lectures · Les notes ouvertes sont mémorisées dans ce navigateur.