Algorithme de sommation de Kahan - Kahan summation algorithm

En analyse numérique , l' algorithme de sommation de Kahan , également connu sous le nom de sommation compensée , réduit considérablement l' erreur numérique dans le total obtenu en ajoutant une séquence de nombres à virgule flottante à précision finie , par rapport à l'approche évidente. Cela se fait en gardant une compensation de fonctionnement séparée (une variable pour accumuler de petites erreurs).

En particulier, la simple addition de nombres en séquence a une erreur dans le pire des cas qui croît proportionnellement à , et une erreur quadratique moyenne qui croît comme pour les entrées aléatoires (les erreurs d'arrondi forment une marche aléatoire ). Avec la sommation compensée, la limite d'erreur dans le pire des cas est effectivement indépendante de , de sorte qu'un grand nombre de valeurs peuvent être sommées avec une erreur qui ne dépend que de la précision à virgule flottante .

L' algorithme est attribué à William Kahan ; Ivo Babuška semble avoir proposé un algorithme similaire de manière indépendante (d'où la sommation de Kahan-Babuška ). Des techniques similaires et antérieures sont, par exemple, l' algorithme de ligne de Bresenham , gardant une trace de l'erreur accumulée dans les opérations entières (bien que documentées pour la première fois à la même époque) et la modulation delta-sigma

L'algorithme

En pseudocode , l'algorithme sera ;

function KahanSum(input)
    var sum = 0.0                    // Prepare the accumulator.
    var c = 0.0                      // A running compensation for lost low-order bits.

    for i = 1 to input.length do     // The array input has elements indexed input[1] to input[input.length].
        var y = input[i] - c         // c is zero the first time around.
        var t = sum + y              // Alas, sum is big, y small, so low-order digits of y are lost.
        c = (t - sum) - y            // (t - sum) cancels the high-order part of y; subtracting y recovers negative (low part of y)
        sum = t                      // Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
    next i                           // Next time around, the lost low part will be added to y in a fresh attempt.

    return sum

Exemple travaillé

Cet exemple sera donné en décimal. Les ordinateurs utilisent généralement l'arithmétique binaire, mais le principe illustré est le même. Supposons que nous sumutilisions une arithmétique décimale à virgule flottante à six chiffres, que nous ayons atteint la valeur 10000,0 et que les deux valeurs suivantes de input[i]soient 3,14159 et 2,71828. Le résultat exact est 10005.85987, qui est arrondi à 10005.9. Avec une sommation simple, chaque valeur entrante serait alignée avec sum, et de nombreux chiffres de poids faible seraient perdus (par troncature ou arrondi). Le premier résultat, après arrondi, serait 10003.1. Le deuxième résultat serait 10005.81828 avant arrondi et 10005.8 après arrondi. Ce n'est pas correct.

Cependant, avec une sommation compensée, nous obtenons le résultat arrondi correct de 10005.9.

Supposons que cla valeur initiale soit zéro.

  y = 3.14159 - 0.00000             y = input[i] - c
  t = 10000.0 + 3.14159
    = 10003.14159                   But only six digits are retained.
    = 10003.1                       Many digits have been lost!
  c = (10003.1 - 10000.0) - 3.14159 This must be evaluated as written! 
    = 3.10000 - 3.14159             The assimilated part of y recovered, vs. the original full y.
    = -0.0415900                    Trailing zeros shown because this is six-digit arithmetic.
sum = 10003.1                       Thus, few digits from input(i) met those of sum.

La somme est si grande que seuls les chiffres de poids fort des nombres d'entrée sont accumulés. Mais à l'étape suivante, cdonne l'erreur.

  y = 2.71828 - (-0.0415900)        The shortfall from the previous stage gets included.
    = 2.75987                       It is of a size similar to y: most digits meet.
  t = 10003.1 + 2.75987             But few meet the digits of sum.
    = 10005.85987                   And the result is rounded
    = 10005.9                       To six digits.
  c = (10005.9 - 10003.1) - 2.75987 This extracts whatever went in.
    = 2.80000 - 2.75987             In this case, too much.
    = 0.040130                      But no matter, the excess would be subtracted off next time.
sum = 10005.9                       Exact result is 10005.85987, this is correctly rounded to 6 digits.

La sommation est donc effectuée avec deux accumulateurs : sumdétient la somme, et caccumule les parties non assimilées dans sum, pour décaler la partie de poids faible de sumla prochaine fois. Ainsi, la sommation se déroule avec des "chiffres de garde" dans c, ce qui est mieux que de ne pas en avoir, mais n'est pas aussi bon que d'effectuer les calculs avec le double de la précision de l'entrée. Cependant, augmenter simplement la précision des calculs n'est pas pratique en général ; si inputest déjà en double précision, peu de systèmes fournissent une quadruple précision , et s'ils le faisaient, inputpourraient alors être en quadruple précision.

Précision

Une analyse minutieuse des erreurs dans la sommation compensée est nécessaire pour apprécier ses caractéristiques de précision. Bien qu'elle soit plus précise que la sommation naïve, elle peut toujours donner des erreurs relatives importantes pour des sommes mal conditionnées.

Supposons que l'on additionne des valeurs , pour . La somme exacte est

(calculé avec une précision infinie).

Avec la sommation compensée, on obtient à la place , où l'erreur est bornée par

où est la précision machine de l'arithmétique utilisée (par exemple pour la virgule flottante double précision standard IEEE ). Habituellement, la quantité d'intérêt est l' erreur relative , qui est donc bornée ci-dessus par

Dans l'expression de la borne d'erreur relative, la fraction est le numéro de condition du problème de sommation. Essentiellement, le numéro de condition représente la sensibilité intrinsèque du problème de sommation aux erreurs, quelle que soit la façon dont il est calculé. La borne d'erreur relative de chaque méthode de sommation ( stable vers l'arrière ) par un algorithme fixe à précision fixe (c'est-à-dire pas ceux qui utilisent l' arithmétique à précision arbitraire , ni les algorithmes dont les exigences en mémoire et en temps changent en fonction des données), est proportionnelle à ce nombre de condition . Un problème de sommation mal conditionné est un problème dans lequel ce rapport est grand, et dans ce cas même une sommation compensée peut avoir une erreur relative importante. Par exemple, si les sommes sont des nombres aléatoires non corrélés avec une moyenne nulle, la somme est une marche aléatoire et le nombre de conditions augmentera proportionnellement à . D'autre part, pour les entrées aléatoires avec une valeur non nulle, le nombre de conditions est asymptote à une constante finie telle que . Si les entrées sont toutes non négatives , alors le numéro de condition est 1.

Étant donné un numéro de condition, l'erreur relative de la sommation compensée est effectivement indépendante de . En principe, il y a le qui croît linéairement avec , mais en pratique ce terme est effectivement nul : puisque le résultat final est arrondi à une précision , le terme s'arrondit à zéro, à moins qu'il ne soit approximativement ou plus grand. En double précision, cela correspond à un d'environ , beaucoup plus grand que la plupart des sommes. Ainsi, pour un nombre de condition fixe, les erreurs de sommation compensée sont effectivement , indépendantes de .

En comparaison, l'erreur relative liée à la sommation naïve (simplement additionner les nombres en séquence, arrondir à chaque étape) augmente multipliée par le nombre de condition. Cette erreur du pire des cas est cependant rarement observée dans la pratique, car elle ne se produit que si les erreurs d'arrondi sont toutes dans le même sens. En pratique, il est beaucoup plus probable que les erreurs d'arrondi aient un signe aléatoire, de moyenne nulle, de sorte qu'elles forment une marche aléatoire ; dans ce cas, la somme naïve a une racine carrée moyenne erreur relative qui croît comme multiplié par le nombre de conditions. C'est encore bien pire que la sommation compensée, cependant. Cependant, si la somme peut être effectuée avec une précision deux fois plus élevée, alors est remplacé par , et la sommation naïve a une erreur dans le pire des cas comparable au terme de la sommation compensée à la précision d'origine.

De même, le qui apparaît ci-dessus est une borne du pire des cas qui ne se produit que si toutes les erreurs d'arrondi ont le même signe (et sont d'une amplitude maximale possible). En pratique, il est plus probable que les erreurs aient un signe aléatoire, auquel cas les termes in sont remplacés par une marche aléatoire, auquel cas, même pour des entrées aléatoires avec une moyenne nulle, l'erreur ne croît que lorsque (en ignorant le terme), le même taux la somme augmente, annulant les facteurs lorsque l'erreur relative est calculée. Ainsi, même pour des sommes asymptotiquement mal conditionnées, l'erreur relative pour la sommation compensée peut souvent être beaucoup plus petite que ce qu'une analyse du pire des cas pourrait suggérer.

Autres améliorations

Neumaier a présenté une version améliorée de l'algorithme de Kahan, qu'il appelle un « algorithme de Kahan-Babuška amélioré », qui couvre également le cas où le prochain terme à ajouter est plus grand en valeur absolue que la somme cumulée, échangeant efficacement le rôle de ce qui est grand et ce qui est petit. En pseudocode , l'algorithme est :

function KahanBabushkaNeumaierSum(input)
    var sum = 0.0
    var c = 0.0                       // A running compensation for lost low-order bits.

    for i = 1 to input.length do
        var t = sum + input[i]
        if |sum| >= |input[i]| then
            c += (sum - t) + input[i] // If sum is bigger, low-order digits of input[i] are lost.
        else
            c += (input[i] - t) + sum // Else low-order digits of sum are lost.
        endif
        sum = t
    next i

    return sum + c                    // Correction only applied once in the very end.

Pour de nombreuses séquences de nombres, les deux algorithmes sont d'accord, mais un exemple simple dû à Peters montre en quoi ils peuvent différer. Pour la sommation en double précision, l'algorithme de Kahan donne 0,0, tandis que l'algorithme de Neumaier donne la valeur correcte 2,0.

Des modifications d'ordre supérieur d'une meilleure précision sont également possibles. Par exemple, une variante suggérée par Klein, qu'il a appelée un « algorithme itératif de Kahan-Babuška » de second ordre. En pseudocode , l'algorithme est :

function KahanBabushkaKleinSum(input)
    var sum = 0.0
    var cs = 0.0
    var ccs = 0.0
    var c
    var cc

    for i = 1 to input.length do
        var t = sum + input[i]
        if |sum| >= |input[i]| then
            c = (sum - t) + input[i]
        else
            c = (input[i] - t) + sum
        endif
        sum = t
        t = cs + c
        if |cs| >= |c| then
            cc = (cs - t) + c
        else
            cc = (c - t) + cs
        endif
        cs = t
        ccs = ccs + cc
    next i

    return sum + cs + ccs

Alternatives

Bien que l'algorithme de Kahan réalise une croissance d'erreur pour la somme de n nombres, seule une croissance légèrement pire peut être obtenue par sommation par paires : on divise récursivement l'ensemble de nombres en deux moitiés, additionne chaque moitié, puis additionne les deux sommes. Ceci a l'avantage de nécessiter le même nombre d'opérations arithmétiques que la sommation naïve (contrairement à l'algorithme de Kahan, qui nécessite quatre fois l'arithmétique et a une latence de quatre fois une sommation simple) et peut être calculé en parallèle. Le cas de base de la récursivité pourrait en principe être la somme d'un seul (ou de zéro) nombres, mais pour amortir le surcoût de la récursivité, on utiliserait normalement un cas de base plus grand. L'équivalent de la sommation par paires est utilisé dans de nombreux algorithmes de transformée de Fourier rapide (FFT) et est responsable de la croissance logarithmique des erreurs d'arrondi dans ces FFT. En pratique, avec les erreurs d'arrondi des signes aléatoires, les erreurs quadratiques moyennes de la sommation par paires augmentent en fait comme .

Une autre alternative consiste à utiliser l'arithmétique de précision arbitraire , qui en principe n'a pas besoin d'être arrondie du tout avec un coût d'effort de calcul beaucoup plus important. Une façon d'effectuer des sommes exactement arrondies en utilisant une précision arbitraire consiste à étendre de manière adaptative à l'aide de plusieurs composants à virgule flottante. Cela minimisera les coûts de calcul dans les cas courants où une haute précision n'est pas nécessaire. Une autre méthode qui utilise uniquement l'arithmétique entière, mais un grand accumulateur, a été décrite par Kirchner et Kulisch ; une implémentation matérielle a été décrite par Müller, Rüb et Rülling.

Invalidation possible par optimisation du compilateur

En principe, un compilateur d'optimisation suffisamment agressif pourrait détruire l'efficacité de la sommation de Kahan : par exemple, si le compilateur simplifiait les expressions selon les règles d' associativité de l'arithmétique réelle, il pourrait "simplifier" la deuxième étape de la séquence

t = sum + y;
c = (t - sum) - y;

à

c = ((sum + y) - sum) - y;

et ensuite à

c = 0;

éliminant ainsi la compensation d'erreur. En pratique, de nombreux compilateurs n'utilisent pas de règles d'associativité (qui ne sont qu'approximatives en arithmétique à virgule flottante) dans les simplifications, à moins que cela ne soit explicitement indiqué par les options du compilateur permettant des optimisations "non sûres", bien que le compilateur Intel C++ soit un exemple qui permet l'associativité -basées sur les transformations par défaut. La version originale K&R C du langage de programmation C permettait au compilateur de réorganiser les expressions à virgule flottante selon les règles d'associativité arithmétique réelle, mais la norme ANSI C qui a suivi interdisait la réorganisation afin de rendre le C mieux adapté aux applications numériques ( et plus similaire à Fortran , qui interdit également la réorganisation), bien qu'en pratique les options du compilateur puissent réactiver la réorganisation, comme mentionné ci-dessus.

Un moyen portable d'inhiber localement de telles optimisations consiste à diviser l'une des lignes de la formulation d'origine en deux instructions et à rendre deux des produits intermédiaires volatils :

function KahanSum(input)
    var sum = 0.0
    var c = 0.0

    for i = 1 to input.length do
        var y = input[i] - c
        volatile var t = sum + y
        volatile var z = t - sum
        c = z - y
        sum = t
    next i

    return sum

Prise en charge par les bibliothèques

En général, les fonctions de "somme" intégrées dans les langages informatiques ne fournissent généralement aucune garantie qu'un algorithme de sommation particulier sera utilisé, encore moins la sommation de Kahan. La norme BLAS pour les sous-routines d' algèbre linéaire évite explicitement d'imposer un ordre de calcul particulier des opérations pour des raisons de performances, et les implémentations BLAS n'utilisent généralement pas la sommation de Kahan.

La bibliothèque standard du langage informatique Python spécifie une fonction fsum pour une sommation exactement arrondie, en utilisant l' algorithme Shewchuk pour suivre plusieurs sommes partielles.

Dans le langage Julia , l'implémentation par défaut de la sumfonction effectue une sommation par paire pour une précision élevée avec de bonnes performances, mais une bibliothèque externe fournit une implémentation de la variante de Neumaier nommée sum_kbnpour les cas où une précision plus élevée est nécessaire.

Dans le langage C# , le package nuget HPCsharp implémente la variante de Neumaier et la sommation par paires : à la fois en tant que scalaire, parallèle aux données à l'aide des instructions du processeur SIMD et multicœur parallèle.

Voir également

Les références

Liens externes