Histoire
John Napier (1550 – 1617) est un mathématicien, physicien et astronome écossais qui inventa le logarithme. La première table de logarithme publiée par Napier était conçue pour effectuer des multiplications par le sinus d’angles et n’utilisait pas la fonction logarithme telle qu’elle est connue de nos jours (référence : Pochon, Luc-Olivier. « A propos d’une table de logarithmes. »).
La fonction logarithme fut extrêmement utile pour simplifier les calculs de multiplication et divisions en de simples additions et soustractions, en s’aidant de tables de logarithme à une époque où les calculatrices n’existaient pas. Ces tables sont restées une manière très efficace et précise de faire des multiplications de nombre décimaux jusqu’au XXième siècle. Les règles à calcul, basées elles aussi sur l’échelle logarithmique permettaient de faire des multiplications et divisions approximatives encore plus rapidement. Ces outils devirent progressivement désuets durant la deuxième moitié du XXième siècle, avec la concurrence des calculatrices mécaniques capables de faire des multiplications, puis dans les années 1970 et 1980, les calculatrices électroniques dont les prix, rapidement décroissants permirent à tous les foyers de s’équiper. Le logarithme et l’exponentielle restent des outils mathématiques indispensables, utilisés au quotidien par le biostatisticien, par exemple pour estimer les paramètres d’une régression logistique ou calculer une variance de Greenwood.
Problématique
Vous avez peut-être déjà atteint les limites des nombres à virgule flottante 64 bits que l’on trouve dans tous les logiciels statistiques. Un nombre inférieur en valeur absolue, à 10^-323 ou supérieur à 10^308, sera représenté comme 0 ou +Inf. Ces limites peuvent paraître, au premier abord, difficiles à atteindre, mais en pratique, on les atteint très vite dès qu’on cumule de nombreuses multiplications ou des exponentiations. Par exemple, dès qu’on calcule la vraisemblance d’un échantillon, on passe rapidement en dessous de 10^-323. Sur un échantillon de taille 1000, avec un pourcentage à 90%, on a une probabilité de 10^-1000 d’observer zéro événement sur les 1000 observations indépendantes, selon la loi binomiale. La probabilité d’observer une valeur au-delà de +100 écarts-types d’une loi normale est d’environ 3×10^-2218. Ce genre de probabilité devra être calculé dans de nombreuses procédures statistiques, telles que les estimateurs du maximum de vraisemblance. C’est une des raisons qui justifie le recours à la log-vraisemblance plutôt que la vraisemblance brute. Grâce à cette transformation logarithmique, les multiplications et les puissances sont respectivement remplacées par des additions et multiplications. On peut alors représenter des vraisemblances de l’ordre de 10^(-10^308).
Conceptualisation de l’outil
Trois points de vue peuvent être posés sur ces transformations logarithmiques :
- Il s’agit d’un outil mathématique simplifiant parfois les algorithmes (p.e. dérivation souvent plus simple sur la log-vraisemblance que sur la vraisemblance) ;
- Il s’agit d’un outil informatique permettant de dépasser les limitations des nombres à virgule flottante, mais obligeant à réécrire les algorithmes, ce qui est parfois délicat pour les novices ;
- Il s’agit d’un détail de représentation interne de l’ordinateur, permettant de gérer des nombres extrêmement grands et petits. Il n’y a pas de mise à jour notable d’un algorithme pour utiliser ces grands nombres.
Les points numéros 1 et 2 correspondent à la vision la plus classique de l’usage du log. Le point numéro 3 est un peu plus original et correspond à un package R que je suis en train de développer.
L’idée, c’est de représenter un nombre quelconque dans l’espace des nombres réels, par deux valeurs:
- Une valeur logique (TRUE/FALSE) représentant le signe positif ou négatif du nombre
- Une valeur numérique flottante (64 bits) représentant le logarithme du nombre
Une telle représentation sera nommée un flog, mot valise issu de float et log. Ainsi, des nombres aussi petits ou grands que exp(-10^308) et exp(10^308) peuvent être représentés par un flog. Les fonctions arithmétiques de base sont assez facilement réécrites afin de permettre un usage tout à fait normal des flogs, comme s’il s’agissait de valeurs numériques habituelles : additions, soustractions, divisions, multiplications, exponentiations, comparaisons peuvent être réalisées sur les flogs de manière transparente, comme si on manipulait les nombres non transformés. Grâce aux capacités de surcharge des opérateurs de R, un algorithme qui jadis, était programmé pour des nombres à virgule flottante normaux (allant de 10^-323 à 10^308) peut être réutilisé pour les flogs sans modification, gérant ainsi des nombres allant de 10^(-10^308) à 10^(10^308).
Apparement, une recherche Google permet de découvrir que je ne suis pas le premier à avoir eu cette idée, puisqu’un package Haskell existe pour ça (https://hackage.haskell.org/package/logfloat-0.13.4/docs/Data-Number-LogFloat.html).
Outre la simplification donnée par les flogs par rapport à une procédure plus classique de transformation logarithmique, ceux-là ont l’avantage d’être capables de représenter les nombres négatifs comme les nombres positifs.
Rapprochement d’un ancien concept : le float
Après tout, les nombres à virgule flottante (float) ont une représentation interne avec un signe, une mantisse et un exposant (qui se rapproche du log du nombre) de manière transparente. On les manipule comme s’il s’agissait de nombre réels, en ignorant les détails de représentation interne, et notamment, le fait qu’ils sont représentés en base 2. Les flogs poussent le concept un cran au-dessus.
En langage C++, on pourrait écrire un template flog<T>, paramétré par un type numérique, tel que float ou double, qui transformerait ce type numérique en un flog. Comme le flog lui-même est un type numérique, on pourrait créer un flog<flog<double> > qui serait une représentation numérique de très très grand nombres dont le logarithme est lui-même représentable comme un flog. C’est flog<flog<double> > pourraient ainsi représenter des nombres allant jusqu’à 10^(10^(10^308)). On peut appeler ça un flog carré. Bien sûr, on peut aller au-delà, avec les flog<flog<flog<double> > > que l’on appelerait les flogs cubes, et d’une manière plus générale, les flogs itérés à l’ordre k.
Les flogs itérés peuvent paraître superflus, mais mon expérience personnelle m’a montré que l’on pouvait malheureusement dépasser les capacités des flogs simples avec la famille des lois de Weibull. Si on veut représenter correctement la probabilité d’une survie dépassant un certain seuil. Par exemple, un flog est incapable de représenter le taux de survie au temps t=200 d’une distribution de Weibull de paramètre de forme égal à 200 et de paramètre d’échelle égal à 1. Un flog carré est capable de représenter le taux de survie au temps t=10^300 d’une distribution de Weibull de paramètre de forme égal à 10^300. On peut aussi utiliser un cflolog dont je parlerais plus tard.
Limitations générales des floats et flogs
Les limitations des nombres à virgule flottante (floats) sont exacerbés sur les flogs. Tout nombre réel n’est pas représentable avec un float ou un flog. Au maximum 2^64 valeurs distinctes sont représentables. Ainsi, on peut définir la précision d’un nombre float ou flog par la différence qui existe entre les deux nombres successifs représentables sur le float ou flog, sans que ne soit représentable le moindre autre nombre entre les deux. Ainsi, la précision de la valeur 1024 pour un float 64 bits est de 2.3×10^-13, parce que le plus petit nombre strictement supérieur à 1024 représentable par un float 64 bits est égal à (1024+2.3×10^-13). Il s’agit de la précision absolue. La précision relative est obtenue en divisant la précision absolue par le nombre considéré (1024 dans l’exemple). C’est ainsi que le nombre 1024 a une précision relative de 2.2×10^-16.
Contourner les problèmes de précision sur les floats
Les floats ont une précision relative presque indépendante de la valeur; cette précision relative varie entre 1.1×10^-16 et 2.2×10^-16 (sauf pour les dénormaux, mais c’est une autre histoire). On peut donc considérer qu’un float 64 bits a 15 à 16 chiffres numériquement significatifs, peu importe sa valeur. La précision absolue sera très mauvaise pour des valeurs énormes, comme 10^300 qui aura une précision absolue d’environ 10^284, mais sera excellente pour des valeurs petites, comme 10^-280 qui aura une précision absolue d’environ 10^-296. En conséquence, il faut être très délicat avec l’ordre dans lequel on effectue les opérations afin de ne pas perdre en précision inutilement. Imaginons que l’on veuille calculer z=y+x-y/(1+x^2) en sachant que x est petit (p.e. 10^-9) alors que y est grand (p.e. 10^9). Si on calcule cela directement avec des float 64 bits, on obtiendra, à tort, la valeur zéro comme résultat, parce que y+x = y et 1+x^2 = 1 en raison des erreurs d’arrondi. Si on change l’ordre du calcul à z=y-y/(1+x^2)+x, alors on obtiendra, pour résultat, 10^-9, parce qu’on aura supprimé l’erreur d’arrondi correspondant à l’addition de x. Cependant, on n’aura pas supprimé l’erreur d’arrondi conduisant à 1+x^2 = 1 alors qu’en réalité, 1+x^2 = 1+10^-18. Pour corriger cette deuxième erreur d’arrondi, il va falloir réécrire y-y/(1+x^2)+x = y*(1-1/(1+x^2))+x puis se concentrer sur l’estimation de 1-1/(1+x^2). Considérant que x^2 est très proche de zéro, on peut approximer la fonction inverse à son développement limité à l’ordre 1, et donc, approximer 1/(1+x^2) à 1-x^2; d’où une approximation de 1-1/(1+x^2) à x^2. On réécrit alors z ≃ y×x^2+x qui se calcule alors de manière précise à 2×10^-9. Paradoxalement, en utilisant un algorithme approximatif (développement limité), on obtient un résultat plus précis qu’avec un algorithme « exact » à ceci-près que ce dernier ignore le problème d’imprécision des nombres à virgule flottante. En pratique, si on veut implémenter cela de manière générique et propre, on utilisera un développement limité à un plus grand ordre (p.e. 1-1/(1+x^2) ≃ x^2 – x^4 +x^6 – x^8) et on définira un seuil sur la valeur de x (p.e. 10^-2) en-dessous duquel on utilise le développement limité, et au-dessus duquel, on utilise le calcul direct 1-1/(1+x^2). Ce seuil doit être choisi de telle sorte que la précision des deux méthodes soit satisfaisante à ce seuil. Ici, l’erreur absolue du développement limité est environ égal à x^10, c’est-à-dire, 10^-20 au niveau du seuil, soit une erreur relative d’environ 10^-16. C’est donc une précision satisfaisante. Pour le rapport 1-1/(1+x^2) calculé directement, l’erreur absolue de 1+x^2 est d’environ 10^-16, ce qui conduit à une erreur absolue du même ordre de grandeur pour son inverse ainsi que pour 1-1/(1+x^2) et donc, une erreur relative d’environ 10^-12. La précision est moins bonne, mais reste satisfaisante. Les développements limités s’avèrent très précieux dans la manipulation des nombres à virgule flottante. Pour cela, il ne faut pas hésiter à s’aider du logiciel libre Xcas qui calcule assez bien les développements limités de manière formelle. Pour certaines fonctions fréquemment utilisées, il existe des procédures dans la plupart des logiciels et langages de programmation : expm1 pour exp(x)-1, log1p pour log(1+x). De manière un peu plus rare, on trouve log1pmx(x)=log(1+x)-x et lgamma1p(x)=log(gamma(1+x)). Ces fonctions sont faciles à implémenter avec des développements limités si elles ne sont pas nativement disponibles dans le logiciel.
Limitations supplémentaires des flogs
La précision absolue des floats est décroissante avec leur taille alors que la précision relative est à peu près constante, aux alentours de 2×10^-16 pour les floats 64 bits. Avec les flogs, la précision relative devient décroissante avec la taille des nombres. Ainsi, la table suivante donne des points de repères :
| Valeur | Précision relative d’un float 64 bits | Précision relative d’un flog 64 bits |
| 10^-(10^100) | Non représentable* | 10^(3.9×10^84) |
| 10^-(10^16) | Non représentable* | 54 |
| 10^-(10^15) | Non représentable* | 0.65 |
| 10^-(10^10) | Non représentable* | 3.8×10^-6 |
| 10^-1000 | Non représentable* | 4.5×10^-13 |
| 10^-100 | ~ 2.2×10^-16 | 2.8×10^-14 |
| 10^-15 | ~ 2.2×10^-16 | 7.1×10^-15 |
| 1 | ~ 2.2×10^-16 | 4.9×10^-324 |
| 10^15 | ~ 2.2×10^-16 | 7.1×10^-15 |
| 10^100 | ~ 2.2×10^-16 | 2.8×10^-14 |
| 10^1000 | Non représentable** | 4.5×10^-13 |
| 10^(10^10) | Non représentable** | 3.8×10^-6 |
| 10^(10^15) | Non représentable** | 0.65 |
| 10^(10^16) | Non représentable** | 54 |
| 10^(10^100) | Non représentable** | 1.94×10^84 |
* Non représentable ou représenté comme strictement égal à 0, et donc, infiniement imprécis, selon le point de vue.
** Non représentable ou représenté comme égal à +Inf, et donc, infiniement imprécis, selon le point de vue.
La précision relative des flogs est optimale autour de 1, bien meilleure que celle des floats. La précision relative des flogs reste excellente jusqu’à 10^1000 et est tolérable jusqu’à 10^(10^10). À partir de 10^(10^16), la précision relative descent à 54, ce qui veut dire que le rapport entre deux flogs successifs est égal à 55. En conséquence, en considérant un mode de calcul qui arrondit au nombre représentable le plus proche, si on additionne un nombre avec lui-même, il est inchangé (x+x = x), et si on multiplie un nombre de l’ordre de 10^(10^16), par un facteur inférieur à 27.5, il est inchangé (p.e. x×20 = x). En première impression, l’addition n’a plus aucun sens sur des valeurs aussi élevées, puisque l’addition de deux nombres de cet ordre de grandeur, sera toujours égal au plus grand des deux. Pourtant, il est m’est arrivé d’additionner et soustraire des flogs bien plus grands, sans que cela pose problème. Lorsqu’on veut savoir si 10^(10^270)+ 10^(10^300) est supérieur ou inférieur à 10^(10^271)+10^(10^301), on comprend que les « petits » nombres 10^(10^270) et 10^(10^271) sont totalement négligeables dans la comparaison et peuvent être ignorés.
Pour 10^(10^100), on voit que la précision relative est très mauvaise, puisque deux nombres représentables successifs ont un rapport 1.94×10^84 l’un avec l’autre. Si on multiplie un nombre de cette taille par des milliards de milliards, il reste inchangé (x×10^18 = x). Pourtant, certaines opérations restent signifiantes à cette échelle : un nombre élevé au carré sera toujours estimé avec une bonne précision. En fait, pour passer d’un nombre représentable au plus petit nombre représentable strictement supérieur, il faut élever le premier à la puissance (1+2.2×10^-16) environ. L’opérateur d’exponentiation (^) garde alors pleinement son sens.
Exemple d’implémentation de l’addition de flogs
Comment peut-on calculer a+b à partir de la=log(a) et de lb=log(b) et le représenter sur une échelle logarithmique ? Si on calcule directement log(exp(la)+exp(lb)), les représentations intermédiaires de exp(la) et exp(lb) dépassent les capacités des floats. Il existe pourtant une stratégie simple et efficace de calcul, basé sur le fait que pour toute constante c > 0, exp(la)+exp(lb) = (exp(la)/c + exp(lb)/c)×c et donc log(exp(la)+exp(lb)) = log(exp(la – log(c)) + exp(lb – log(c))) + log(c). Il faut ensuite choisir une constante lc = log(c) qui garantisse que l’on ne dépasse pas les capacités des nombre à virgule flottante. Pour cela, la constante lc=max(la, lb) est idéale. Ainsi, la-lc et lb-lc seront tous deux inférieurs ou égaux à 0 et on ne dépassera pas les capacités des nombre à virgule flottante. Si on additionne des nombres qui n’ont pas le même ordre de grandeur, tels que e^10000 et e^10100, le plus petit des deux sera négligeable face au plus grand, et le résultat final sera égal au plus grand des deux, en raison des arrondis flottant.
Place des flogs et floats
Les floats sont tout à fait adaptés à des algorithmes ayant un recours massif aux additions et comprenant bon nombre de multiplications, mais on arrive à leurs limites dès que l’on commence à utiliser des puissances ou que le recours aux multiplications est trop massif. C’est là que les flogs ont leur place.
Bon usage des flogs
Les limites des flogs sont semblables à celles de floats, mais il est encore plus nécessaire d’en avoir conscience car ils sont conçus pour manipuler des valeurs extrêmes. C’est ainsi que contrairement à ce que j’ai décrit en introduction, il faut avoir conscience de la représentation des nombres pour écrire des algorithmes corrects. Par exemple, si x est potentiellement immensément grand, positif ou négatif, il faudra absolument éviter d’écrire des expressions telles que log(1+exp(x)), puisque cette expression dépassera les capacités des nombres pour un grand x positif et perdra toute précision pour un grand x négatif. Pour x positif et supérieur à 37 environ, on pourra approximer log(1+exp(x)) à x, sans perte de précision. Pour x négatif et inférieur à 37 environ, on pourra approximer log(1+exp(x)) à exp(x). Pour x entre les deux, on pourra s’aider de la fonction log1p pour écrire log1p(exp(x)).
Pour aller plus loin
La transformation logarithmique n’est pas la seule qui soit possible pour créer un nouvel espace de nombres à partir des floats. On peut notamment, pour des valeurs comprises entre 0 et 1 utiliser la transformation cloglog qui à x associe log(-log(x)) pour construire ainsi les cflologs. Seuls les nombres compris entre 0 et 1 sont représentables ainsi.
Là où les floats sont incapables de représenter des nombres compris entre 1-1.1×10^-16 et 1, les flogs sont capables de représenter 1-5×10^-325 et les cflologs sont capables de représenter 1-10^(-1.8×10^308). Là où le plus petit nombre strictement positif représentable par les floats est 5×10^-325, les flogs peuvent représenter 10^(-1.8×10^308) et les cflologs peuvent représenter 10^(-10^(1.8×10^308)). On peut en effet dépasser les capacités des flogs avec la distribution de Weibull, parce que celle-ci dépend un exponentiation dans une exponentielle, la fonction cumulative de distribution de Weibull étant F(x; forme, échelle) = 1-exp(-(x/échelle)^forme).
Le problème de cette transformation, c’est que l’espace de représentations numériques ainsi créé est limité à l’intervalle [0, 1], interdisant alors bon nombre d’opérations. Si on veut rendre ces représentations polyvalentes, il faut enrichir la représentation. Par exemple, on pourrait proposer deux floats x et y, représentant la valeur exp(-exp(x))×y.
Conclusion
Les flogs sont un nouveau point de vue, sur un outil ancien. Utilisez-les à bon escient !