Bon usage des nombres à virgule flottante : exemples

Calculer une variance peut paraître simple, mais il y a de très mauvaises formules qui circulent.

Considérons une variance portant sur des nombres énormes avec une très faible variance. Un cas fréquent, ce sont les dates ou paires dates/heures ! Ces valeurs sont typiquement représentées par un nombre de jours, secondes ou millisecondes depuis une certaine époque (1er janvier 1970, 1er janvier 1900 ou 1er janvier 1601). Le plus extrême est obtenu par les dates de création/modification des fichiers Windows qui sont enregistrés comme le nombre d’intervalles de 100 nanosecondes depuis le 1er janvier 1601. Ainsi, pour les dates récentes, les valeurs dépassent les 10 puissance 17.

Même si les dates sont très grandes, on peut vouloir les manipuler sur un intervalle de temps très court, tel que les quelques secondes d’une expérience de marche, comme j’ai pu voir dans un projet de recherche.

Considérons donc le code R suivant (exécuté sur R 4.0.3 sous Windows 10), basé sur deux dates/temps, exprimées en nombre de secondes depuis le 1er janvier 1970:

> x1=strptime("2021-01-01 00:00:01", format="%Y-%m-%d %H:%M:%S")
> x2=strptime("2021-01-01 00:00:02", format="%Y-%m-%d %H:%M:%S")
> var(as.numeric(c(x1,x2)))
[1] 0.5

Tout semble fonctionner bien avec la fonction var de base. Cependant si on réimplemente cette fonction avec la formule que l’on retrouve dans tous les livres, on obtient un résultat aberrant:

> badvar=function(x) (sum(x^2)-sum(x)^2/length(x))/(length(x)-1)
> badvar(as.numeric(c(x1,x2)))
[1] 0

Cela est explicable par le fait qu’on fait la soustraction de deux nombres très grand (sum(x^2) et sum(x)^2), dont la différence devrait être petite.

> sum(x^2) # représentation inexacte (dépasse 2^53, soit 9.0e15)
[1] 5.180695e+18
> sum(x)^2/(length(x))
[1] 5.180695e+18
> sum(x^2) - sum(x)^2/(length(x))
[1] 0

Une meilleure fonction commence par enlever la moyenne à toutes les valeurs, en utilisant la formule la plus naïve de la variance:

> bettervar=function(x) {sum((x-mean(x))^2)/(length(x)-1)}
> bettervar(x)
[1] 0.5

Sur de très grands vecteurs cette fonction peut quand même avoir une perte de précision parce que sa manière de calculer la moyenne des carrés, passe par une somme très grande.

Pour mieux comprendre ça, comparons l’implémentation de R de la fonction mean, à une implémentation naïve:

> badmean=function(v) {sum(v)/length(v)}
> v=rep(c(0,2^50 + 100), 2^24)
> mean(v[1:2])-2^49 # correct
[1] 50
> mean(v)-2^49 # correct aussi
[1] 50
> badmean(v)-2^49 # incorrect
[1] 7.5
> rowMeans(matrix(nrow=1,v))-2^49 # incorrect
[1] 7.5

En effet, la fonction sum() va conduire à un nombre aux alentours de 2^74, avec une représentation intermédiaire inexacte. On peut corriger ça, en rectifiant le biais, par deux appels à badmean.

> bettermean=function(v) {m = badmean(v); bias = badmean(v - m) ; m+bias}
> bettermean(v)-2^49
[1] 50

En effet, on commence par calculer une moyenne approximative avec badmean, puis on calcule la moyenne des écarts à cette moyenne, qui devrait correspondre à l’opposé du biais de la moyenne. Ce biais est estimé avec des sommes qui oscillent autour de zéro et ne vont donc pas dans des ordres de grandeur trop grands.

Sauf qu’il y a quand même un problème… Ces oscillations de la somme autour de zéro disparaissent si les valeurs sont plus ou moins ordonnées:

> w=sort(v)
> mean(w)-2^49 # incorrect
[1] 20.375
> bettermean(w)-2^49 # incorrect
[1] 20.375
> mean(w) - mean(v) # ouille
[1] -29.625

Ainsi, la fonction bettermean, comme la fonction mean de R, dépendent de l’ordre des valeurs dans le vecteur ! Heureusement sur plateforme PC, les FP80 utilisés par sum() fournissent une marge de sécurité durant le processus de sommation des nombres positifs et négatifs, de telle sorte que ce problème n’arrive que pour les grands vecteurs.

Par ailleurs la fonction badmean() dépend aussi de l’ordre des valeurs, comme dans l’exemple ci-dessous :

> v=rep(c(-2^52-1, 2^52), 1e6)
> badmean(v)
[1] -0.5
> badmean(sort(v))
[1] -0.002048

On peut alors implémenter encore une variance plus robuste en utilisant la fonction mean() deux fois:

> badvar=function(x) (sum(x^2) - sum(x)^2/length(x))/(length(x)-1)
> bettervar=function(x) {sum((x-mean(x))^2)/(length(x)-1)}
> bestvar = function(x) {mean((x-mean(x))^2)*(length(x)/(length(x)-1))}
> cst=2^25+10
> off=2^49
> v=rep(off+c(-cst, cst), each=1e7)
> vartheo = cst^2*(length(v)/(length(v)-1))
> bestvar(v) - vartheo # perfect
[1] 0
> bettervar(v) - vartheo # small error
[1] -87.5
> var(v) - vartheo # small error
[1] -87.5
> badvar(v) - vartheo # huge error
[1] -1.480078e+15
> badvar(v) # negative variance!
[1] -3.541775e+14