Algorithme de Kahan: comment obtenir la différence exacte des produits

Je suis récemment revenu à l'analyse des erreurs en virgule flottante pour affiner certains détails de la prochaine édition du livre sur le rendu physique . Les nombres à virgule flottante sont un domaine de calcul intéressant, plein de surprises (bonnes et mauvaises), ainsi que des astuces délicates pour se débarrasser des mauvaises surprises.

Dans le processus, je suis tombé sur ce post sur StackOverflow , à partir duquel j'ai appris un algorithme élégant pour le calcul exact $$$a \ fois b-c \ fois d$ .

Mais avant de poursuivre avec l'algorithme, vous devez comprendre ce qui est si rusé dans l'expression $$$a \ fois b-c \ fois d$ ? Prenez $$$a = 33962.035$ , $$$b = -30438,8$ , $$$c = 41563,4$ et $$$d = -24871,969$ . (Ce sont les vraies valeurs que j'ai obtenues lors du lancement de pbrt .) Pour les valeurs flottantes 32 bits, nous obtenons: $$$a \ fois b = -1.03376365 \ fois 10 ^ 9$ et $$$c \ fois d = -1.03376352 \ fois 10 ^ 9$ . Nous effectuons la soustraction, et nous obtenons $$$-128$ . Mais si vous effectuez les calculs avec une double précision, et à la fin les convertissez en float, vous obtenez $$$-75.1656$ . Qu'est-il arrivé?

Le problème est que la valeur de chaque œuvre peut aller bien au-delà du résultat net $$$-1 \ fois 10 ^ 9$ , où la distance entre les valeurs à virgule flottante représentables est très grande - 64. Autrement dit, lors de l'arrondi $$$a \ fois b$ et $$$c \ fois d$ individuellement, au flotteur représentable le plus proche, ils se transforment en nombres qui sont des multiples de 64. À leur tour, leur différence sera un multiple de 64, et il n'y aura aucun espoir qu'elle devienne $$$-75.1656$ plus proche que $$$-64$ . Dans notre cas, le résultat était encore plus loin en raison de l'arrondi des deux œuvres $$$-1 \ fois 10 ^ 9$ . Nous rencontrerons directement la bonne vieille réduction catastrophique ¹ .

Voici une meilleure solution que ² :

inline float DifferenceOfProducts(float a, float b, float c, float d) { float cd = c * d; float err = std::fma(-c, d, cd); float dop = std::fma(a, b, -cd); return dop + err; } 

DifferenceOfProducts() calcule $$$a \ fois b-c \ fois d$ d'une manière qui évite une contraction catastrophique. Cette technique a été décrite pour la première fois par le légendaire William Kahan dans l'article sur le coût du calcul en virgule flottante sans arithmétique extra-précise . Il convient de noter que le travail de Kahan est intéressant à lire dans son ensemble, ils ont de nombreux commentaires sur l'état actuel du monde des virgules flottantes, ainsi que des considérations mathématiques et techniques. Voici une de ses conclusions:

Ceux d'entre nous qui ont lutté avec les vicissitudes de l'arithmétique en virgule flottante et des "optimisations" mal pensées des compilateurs peuvent à juste titre être fiers de la victoire dans cette bataille. Mais si nous transmettons la poursuite de cette bataille aux générations futures, cela contredira toute l'essence de la civilisation. Notre expérience montre que les langages de programmation et les systèmes de développement sont des sources de trop de chaos auxquels nous devons faire face. Trop d'erreurs peuvent être supprimées, ainsi que quelques «optimisations» attrayantes qui sont sans danger pour les entiers, mais s'avèrent parfois fatales pour les nombres à virgule flottante.

Après avoir rendu hommage à son esprit, revenons à DifferenceOfProducts() : la base de la maîtrise de cette fonction est l'utilisation d'instructions FMA ³ multipliées-ajoutées fusionnées. D'un point de vue mathématique, FMA(a,b,c) est $$$a \ fois b + c$ Par conséquent, au début, il semble que cette opération ne soit utile que comme microoptimisation: une instruction au lieu de deux. Cependant fma
Il a une propriété spéciale - il ne se termine qu'une seule fois.

Comme d'habitude $$$a \ fois b + c$ premier calculé $$$a \ fois b$ , puis cette valeur, qui dans le cas général ne peut pas être représentée au format virgule flottante, est arrondie au flottant le plus proche. Ensuite, à cette valeur arrondie est ajoutée $$$c$ , et ce résultat est à nouveau arrondi au flottant le plus proche. FMA est implémenté de telle manière que l'arrondi n'est effectué qu'à la fin - une valeur intermédiaire $$$a \ fois b$ conserve une précision suffisante, donc après l'avoir ajouté $$$c$ le résultat final sera le plus proche de la vraie valeur $$$a \ fois b + c$ valeur du flotteur.

Après avoir traité avec FMA, nous reviendrons sur DifferenceOfProducts() . Encore une fois, je vais en montrer les deux premières lignes:

  float cd = c * d; float err = std::fma(-c, d, cd); 

Le premier calcule la valeur arrondie $$$c \ fois d$ et le second ... soustrait $$$c \ fois d$ de leur travail? Si vous ne savez pas comment fonctionne FMA, vous pourriez penser que l' err sera toujours nulle. Mais lorsque vous travaillez avec FMA, la deuxième ligne extrait en fait la valeur de l'erreur d'arrondi dans la valeur calculée $$$c \ fois d$ et l'enregistre pour se err . Après cela, la sortie est très simple:

  float dop = std::fma(a, b, -cd); return dop + err; 

Le deuxième FMA calcule la différence entre les œuvres utilisant le FMA, en n'effectuant l'arrondi qu'à la toute fin. Par conséquent, il résiste à une réduction catastrophique, mais il doit fonctionner avec une valeur arrondie. $$$c \ fois d$ . L' return "corrige" ce problème, en ajoutant l'erreur mise en évidence dans la deuxième ligne. Dans un article de Jeannenrod et al. il est démontré que le résultat est vrai jusqu'à 1,5 ulps (mesures de précision unitaire), ce qui est excellent: les opérations FMA et simples en virgule flottante sont précises jusqu'à 0,5 ulps, donc l'algorithme est presque parfait.

Nous utilisons un nouveau marteau

Lorsque vous commencez à chercher des moyens d'utiliser DifferenceOfProducts() , cela s'avère étonnamment utile. Calcul du discriminant de l'équation quadratique? Appelez DifferenceOfProducts(b, b, 4 * a, c) ⁴ . Calcul du déterminant d'une matrice 2x2? L'algorithme résoudra ce problème. Dans la prochaine version de pbrt, j'ai trouvé environ 80 utilisations. De tous, la fonction du produit vectoriel est la plus appréciée. Cela a toujours été une source de problèmes, à cause duquel vous avez dû lever la main et utiliser le double dans la mise en œuvre pour éviter une réduction catastrophique:

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { double v1x = v1.x, v1y = v1.y, v1z = v1.z; double v2x = v2.x, v2y = v2.y, v2z = v2.z; return Vector3f(v1y * v2z - v1z * v2y, v1z * v2x - v1x * v2z, v1x * v2y - v1y * v2x); } 

Et maintenant, nous pouvons continuer à travailler avec float et utiliser DifferenceOfProducts() .

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { return Vector3f(DifferenceOfProducts(v1.y, v2.z, v1.z, v2.y), DifferenceOfProducts(v1.z, v2.x, v1.x, v2.z), DifferenceOfProducts(v1.x, v2.y, v1.y, v2.x)); } 

Cet exemple astucieux du début de l'article fait en fait partie d'un travail vectoriel. À un certain stade, le code pbrt doit calculer le produit vectoriel des vecteurs $$$(33962.035, 41563.4, 7706.415)$ et $$$(- 24871.969, -30438.8, -5643.727)$ . Lorsque calculé à l'aide de float, nous obtiendrions un vecteur $$$(1552, -1248, -128)$ . (Règle générale: si dans les calculs en virgule flottante où de grands nombres sont impliqués, vous n'obtenez pas des valeurs entières si grandes, alors c'est un signe presque sûr qu'une réduction catastrophique s'est produite.)

Avec une double précision, le produit vectoriel est $$$(1556.0276, -1257.5151, -75.1656)$ . On voit qu'avec flotteur $$$x$ semble normal $$$y$ déjà assez mauvais aussi $$$z$ devient la catastrophe qui est devenue la motivation pour trouver une solution. Et quels résultats obtiendrons-nous avec DifferenceOfProducts() et les valeurs flottantes? $$$(1556.0276, -1257.5153, -75.1656)$ . Les valeurs $$$x$ et $$$z$ correspondent à une double précision, et $$$y$ légèrement décalé - d'où le ulp supplémentaire est venu.

Et la vitesse? DifferenceOfProducts() effectue deux FMA, ainsi que la multiplication et l'addition. Un algorithme naïf peut être implémenté avec un FMA et une multiplication, ce qui, semble-t-il, devrait prendre la moitié du temps. Mais en pratique, il s'avère qu'après avoir obtenu les valeurs des registres, cela n'a pas d'importance: dans le cas-test synthétique tenu sur mon ordinateur portable, DifferenceOfProducts() que 1,09 fois plus cher que l'algorithme naïf. Le fonctionnement en double précision était 2,98 fois plus lent.

Dès que vous apprenez la possibilité d'une réduction catastrophique, toutes sortes d'expressions d'apparence innocente dans le code commencent à sembler suspectes. DifferenceOfProducts() semble être un bon remède pour la plupart d'entre eux. Il est facile à utiliser et nous n'avons aucune raison particulière de ne pas l'utiliser.

Remarques

1. La réduction catastrophique n'est pas un problème lors de la soustraction de quantités avec des signes différents ou lors de l'ajout de valeurs avec le même signe. Cependant, cela peut devenir un problème lors de l'ajout de valeurs avec des signes différents. Autrement dit, les montants doivent être considérés avec le même soupçon que les différences.
2. Comme exercice pour le lecteur, je suggère d'écrire la fonction SumOfProducts() , qui offre une protection contre les contractions catastrophiques. Si vous voulez compliquer la tâche, expliquez pourquoi dans DifferenceOfProducts() , dop + err == dop , car les signes a*b et c*d sont différents.
3. L'instruction FMA est disponible sur le GPU depuis plus d'une décennie et sur la plupart des CPU depuis au moins cinq ans. Il peut être nécessaire d'ajouter des drapeaux de compilation au CPU pour les générer directement lors de l'utilisation de std::fma() ; dans gcc et clang, -march=native travaux -march=native .
4. Dans le format à virgule flottante IEEE, la multiplication par des puissances de deux est effectuée exactement, donc 4 * a ne provoque aucune erreur d'arrondi, sauf en cas de débordement.