Analyse d'une erreur numérique dans la fonction C ++

Supposons que j'ai une fonction qui prend en entrée plusieurs valeurs à virgule flottante (simple ou double), effectue certains calculs et produit des valeurs à virgule flottante en sortie (également simples ou doubles). Je travaille principalement avec MSVC 2008, mais je prévois également de travailler avec MinGW / GCC. Je programme en C ++.

Quelle est la façon typique de mesurer par programme la quantité d'erreur dans les résultats? En supposant que j'ai besoin d'utiliser une bibliothèque de précision arbitraire: quelle est la meilleure bibliothèque de ce type si je ne me soucie pas de la vitesse?

Si vous recherchez une bonne limite à votre erreur d'arrondi, vous n'avez pas nécessairement besoin d'une bibliothèque de précision aribtraire. Vous pouvez plutôt utiliser l'analyse des erreurs en cours d'exécution.

Je n'ai pas pu trouver une bonne référence en ligne, mais tout est décrit dans la section 3.3 du livre de Nick Higham "Accuracy and Stability of Numerical Algorithms". L'idée est assez simple:

1. Re-factorisez votre code afin d'avoir une seule affectation d'une seule opération arithmétique sur chaque ligne.
2. Pour chaque variable, par exemple x, créez une variable x_errqui est initialisée à zéro lorsque xune constante lui est affectée.
3. Pour chaque opération, par exemple z = x * y, mettez à jour la variable en z_errutilisant le modèle standard d'arithmétique à virgule flottante zet les erreurs résultantes et en cours d'exécution x_erret y_err.
4. La valeur de retour de votre fonction doit alors également avoir une _errvaleur respective attachée à elle. Il s'agit d'une limite dépendante des données de votre erreur d'arrondi totale.

La partie délicate est l'étape 3. Pour les opérations arithmétiques les plus simples, vous pouvez utiliser les règles suivantes:

• z = x + y -> z_err = u*abs(z) + x_err + y_err
• z = x - y -> z_err = u*abs(z) + x_err + y_err
• z = x * y -> z_err = u*abs(z) + x_err*abs(y) + y_err*abs(x)
• z = x / y -> z_err = u*abs(z) + (x_err*abs(y) + y_err*abs(x))/y^2
• z = sqrt(x) -> z_err = u*abs(z) + x_err/(2*abs(z))

où u = eps/2est l'arrondi de l'unité. Oui, les règles +et -sont les mêmes. Les règles pour toute autre opération op(x)peuvent être facilement extraites en utilisant l'expansion de la série Taylor du résultat appliqué op(x + x_err). Ou vous pouvez essayer de googler. Ou en utilisant le livre de Nick Higham.

À titre d'exemple, considérons le code Matlab / Octave suivant qui évalue un polynôme dans les coefficients aà un point en xutilisant le schéma de Horner:

function s = horner ( a , x )
    s = a(end);
    for k=length(a)-1:-1:1
        s = a(k) + x*s;
    end

Pour la première étape, nous avons divisé les deux opérations en s = a(k) + x*s:

function s = horner ( a , x )
    s = a(end);
    for k=length(a)-1:-1:1
        z = x*s;
        s = a(k) + z;
    end

Nous introduisons ensuite les _errvariables. Notez que les entrées aet xsont supposées être exactes, mais nous pourrions tout aussi bien exiger que l'utilisateur passe des valeurs correspondantes pour a_erret x_err:

function [ s , s_err ] = horner ( a , x )
    s = a(end);
    s_err = 0;
    for k=length(a)-1:-1:1
        z = x*s;
        z_err = ...;
        s = a(k) + z;
        s_err = ...;
    end

Enfin, nous appliquons les règles décrites ci-dessus pour obtenir les termes d'erreur:

function [ s , s_err ] = horner ( a , x )
    u = eps/2;
    s = a(end);
    s_err = 0;
    for k=length(a)-1:-1:1
        z = x*s;
        z_err = u*abs(z) + s_err*abs(x);
        s = a(k) + z;
        s_err = u*abs(s) + z_err;
    end

Notez que puisque nous n'avons pas a_errou x_err, par exemple, ils sont supposés être nuls, les termes respectifs sont simplement ignorés dans les expressions d'erreur.

Et voilà! Nous avons maintenant un schéma de Horner qui renvoie une estimation d'erreur dépendante des données (remarque: il s'agit d'une limite supérieure de l'erreur) à côté du résultat.

En guise de remarque, puisque vous utilisez C ++, vous pouvez envisager de créer votre propre classe de valeurs à virgule flottante qui porte le _errterme et de surcharger toutes les opérations arithmétiques pour mettre à jour ces valeurs comme décrit ci-dessus. Pour les grands codes, il peut s'agir de l'itinéraire le plus facile, bien que moins efficace sur le plan des calculs. Cela dit, vous pourrez peut-être trouver un tel cours en ligne. Une recherche rapide sur Google m'a donné ce lien .

PS Notez que tout cela ne fonctionne que sur les machines respectant strictement IEEE-754, c'est-à-dire que toutes les opérations arithmétiques sont précises à . Cette analyse donne également une limite plus stricte et plus réaliste que l'utilisation de l'arithmétique d'intervalle car, par définition, vous ne pouvez pas représenter un nombre en virgule flottante, c'est-à-dire que votre intervalle serait juste arrondi au nombre lui-même.x ( 1 ± u $\pm u$$x(1 \pm u)$

Une belle bibliothèque portable et open source pour l'arithmétique à virgule flottante de précision arbitraire (et bien d'autres) est le NTL de Victor Shoup , qui est disponible sous forme de source C ++.

À un niveau inférieur se trouve la bibliothèque Bignum GNU Multiple Precision (GMP) , également un package open source.

NTL peut être utilisé avec GMP si des performances plus rapides sont requises, mais NTL fournit ses propres routines de base qui sont certainement utilisables si vous "ne vous souciez pas de la vitesse". GMP prétend être la "bibliothèque bignum la plus rapide". GMP est largement écrit en C, mais possède une interface C ++.

Ajouté: Bien que l'arithmétique des intervalles puisse donner des limites supérieures et inférieures à la réponse exacte de manière automatisée, cela ne mesure pas précisément l'erreur dans un calcul de précision "standard" car la taille de l'intervalle augmente généralement avec chaque opération (soit dans un rapport sens d'erreur absolu).

La manière typique de trouver la taille d'erreur, soit pour les erreurs d'arrondi, soit pour les erreurs de discrétisation, etc., est de calculer une valeur de précision supplémentaire et de la comparer à la valeur de précision "standard". Seule une précision supplémentaire modeste est nécessaire pour déterminer la taille de l'erreur elle-même avec une précision raisonnable, car les erreurs d'arrondi seules sont sensiblement plus grandes en précision "standard" que dans le calcul de précision supplémentaire.

Le point peut être illustré en comparant les calculs de simple et double précision. Notez qu'en C ++, les expressions intermédiaires sont toujours calculées en (au moins) double précision, donc si nous voulons illustrer à quoi ressemblerait un calcul en simple précision "pure", nous devons stocker les valeurs intermédiaires en simple précision.

Extrait de code C

    float fa,fb;
    double da,db,err;
    fa = 4.0;
    fb = 3.0;
    fa = fa/fb;
    fa -= 1.0;

    da = 4.0;
    db = 3.0;
    da = da/db;
    da -= 1.0;

    err = fa - da;
    printf("Single precision error wrt double precision value\n");
    printf("Error in getting 1/3rd is %e\n",err);
    return 0;

La sortie d'en haut (chaîne d'outils Cygwin / MinGW32 GCC):

Single precision error wrt double precision value
Error in getting 1/3rd is 3.973643e-08

Ainsi, l'erreur concerne ce que l'on attend en arrondissant 1/3 à la précision simple. On ne se soucierait pas (je soupçonne) d'obtenir plus de quelques décimales dans l'erreur correcte, car la mesure de l'erreur est pour l'amplitude et non pour l'exactitude.

GMP (c'est-à-dire la bibliothèque GNU Multiple Precision) est la meilleure bibliothèque de précision arbitraire que je connaisse.

Je ne connais aucun moyen programmatique de mesurer l'erreur dans les résultats d'une fonction arbitraire en virgule flottante. Une chose que vous pourriez essayer est de calculer une extension d'intervalle d'une fonction en utilisant l' arithmétique d'intervalle . En C ++, vous devrez utiliser une sorte de bibliothèque pour calculer les extensions d'intervalle; une de ces bibliothèques est la bibliothèque arithmétique d'intervalle Boost. Fondamentalement, pour mesurer l'erreur, vous fourniriez comme arguments à vos intervalles de fonction qui ont une largeur de 2 fois l'unité arrondie (à peu près), centrée sur les valeurs d'intérêt, puis votre sortie serait une collection d'intervalles, la largeur de ce qui vous donnerait une estimation prudente de l'erreur. Une difficulté de cette approche est que l'arithmétique d'intervalle utilisée de cette manière peut surestimer l'erreur de manière significative, mais cette approche est la plus "programmatique" à laquelle je puisse penser.

Une estimation d'erreur rigoureuse et automatique peut être obtenue par analyse d'intervalle . Vous travaillez avec des intervalles plutôt qu'avec des nombres. Par exemple addition:

[a,b] + [c,d] = [min(a+c, a+d, b+c, b+d), max (a+c, a+d, b+c, b+d)] = [a+c, b+d]

L' arrondi peut également être géré de manière rigoureuse, voir Arithmétique des intervalles arrondis .

Tant que votre entrée se compose d'intervalles étroits, les estimations sont correctes et ne sont pas bon marché à calculer. Malheureusement, l'erreur est souvent surestimée, voir le problème de dépendance .

Je ne connais aucune bibliothèque arithmétique d'intervalle de précision arbitraire.

Cela dépend de votre problème actuel si l'arithmétique d'intervalle peut répondre à vos besoins ou non.

La bibliothèque GNU MPFR est une bibliothèque flottante de précision arbitraire qui a une grande précision (en particulier, un arrondi correct pour toutes les opérations, ce qui n'est pas aussi simple qu'il y paraît) comme l'un de leurs principaux points de focalisation. Il utilise GNU MP sous le capot. Il a une extension appelée MPFI qui fait de l'arithmétique d'intervalle, qui - comme le suggère la réponse de Geoff - pourrait être utile à des fins de vérification: continuez à augmenter la précision de travail jusqu'à ce que l'intervalle résultant tombe dans de petites limites.

Cela ne fonctionnera pas toujours, cependant; en particulier, ce n'est pas nécessairement efficace si vous faites quelque chose comme l'intégration numérique, où chaque étape comporte une "erreur" indépendante des problèmes d'arrondi. Dans ce cas, essayez un package spécialisé tel que COZY infinity qui le fait très bien en utilisant des algorithmes spécifiques pour limiter l'erreur d'intégration (et en utilisant des modèles dits de Taylor au lieu d'intervalles).

On me dit que MPIR est une bonne bibliothèque à utiliser si vous travaillez avec Visual Studio:

http://mpir.org/