Algoritmul de însumare al lui Kahan

De la Wikipedia, enciclopedia liberă.
Salt la navigare Salt la căutare

În analiza numerică , algoritmul de însumare Kahan (cunoscut și sub numele de însumare compensată [1] ) reduce semnificativ eroarea numerică din totalul obținut prin adăugarea unei secvențe de numere în virgulă mobilă cu precizie finită, în comparație cu procedura obișnuită. Acest lucru se realizează prin menținerea unei compensări progresive separate (o variabilă pentru acumularea unor erori mici).

Când reprezentăm un număr real generic în virgulă mobilă cu precizie finită, adică cu un număr finit de cifre semnificative , această reprezentare, în raport cu numărul real considerat, diferă de o anumită valoare, care corespunde erorii de rotunjire și este diferența dintre reprezentarea în virgulă mobilă și numărul în sine. Prin efectuarea unei însumări simple a mai multor numere reale, dar folosind reprezentările respective în virgulă mobilă, totalul obținut prezintă o anumită eroare dată de suma algebrică a erorilor de rotunjire unice, adică a singurelor aruncări și, în plus, un anumit pătrat deviația este obținută . medie , în engleză rădăcină înseamnă pătrat eroare , intenționată ca rădăcină pătrată a mediei aritmetice a pătratelor resturilor individuale.

În special, sumarea simplă a n numere în ordine prezintă o eroare care, în cel mai rău caz, crește proporțional cu n și o deviație standard care crește ca pentru completări aleatorii (erorile de rotunjire produc o plimbare aleatorie ). [2] Pe ​​de altă parte, cu însumarea compensată, cea mai gravă eroare posibilă este independentă de n , deci se poate adăuga un număr mare de valori cu o eroare care depinde doar de precizia reprezentării în virgulă mobilă. [2]

Algoritmul este atribuit lui William Kahan . [3] Tehnici anterioare similare sunt, de exemplu, algoritmul de linie Bresenham , care ține evidența erorii acumulate în operațiile întregi (chiar dacă este deja prezentă într-un articol publicat aproximativ în aceeași perioadă [4] ) și modulația Sigma-Delta [5] ] (care integrează și nu doar rezumă eroarea).

Algoritmul

În pseudocod , algoritmul este:

 funcția KahanSum (intrare)
    var sum = 0,0
    var c = 0,0 // O compensare în curs pentru biții cei mai puțin semnificativi pierduți.
    pentru i = 1 la input.length do
        var y = input [i] - c // Până acum, atât de bine: c este zero.
        var t = suma + y // Din păcate, suma este mare, y este mic, astfel încât mai puțin cifrele semnificative ale y sunt pierdute.
        c = (t - suma) - y // (t - suma) anulează cea mai semnificativă parte din y ; scăderea y recuperează partea mai puțin semnificativă din y schimbată în semn
        sum = t // Algebric, c ar trebui să fie întotdeauna zero. Feriți-vă de compilatoare cu optimizare excesiv de agresivă!
    următoarea i // Data viitoare, cea mai puțin semnificativă parte pierdută va fi adăugată la „y” într-o reîncercare. .
    returnează suma

Exemplu practic

Acest exemplu va fi dat în zecimal. Calculatoarele folosesc de obicei aritmetica binară, dar principiul ilustrat acum este același. Să presupunem că folosim aritmetică în virgulă zecimală cu șase cifre și că suma a atins valoarea 10000.0 și următoarele două valori ale intrărilor (i) sunt 3.14159 și 2.71828. Rezultatul exact este 10005.85987, care este rotunjit la 10005.9. Cu o sumă simplă, fiecare valoare introdusă ar fi aliniată la sumă , suma și multe cifre mai puțin semnificative s-ar pierde (prin tăiere sau rotunjire). Primul rezultat, după rotunjire, ar fi 10003,1. Al doilea rezultat ar fi 10005.81828 înainte de rotunjire și 10005.8 după. Acest lucru este incorect.

Cu toate acestea, cu însumarea compensată, obținem rezultatul rotunjit corectat de 10005,9.

Presupunem că c are valoarea inițială zero.

 y = 3,14159 - 0 y = intrare [i] - c
  t = 10000,0 + 3,14159
    = 10003.14159 Dar sunt păstrate doar șase cifre. 
    = 10003.1 Multe cifre s-au pierdut!
  c = (10003.1 - 10000.0) - 3.14159 Acest lucru trebuie evaluat așa cum este scris! 
    = 3.10000 - 3.14159 Partea asimilată a lui y , recuperată, spre deosebire de valoarea completă originală a lui y .
    = -.0415900 Se afișează zerouri finale pentru că este aritmetică din șase cifre.
sum = 10003.1 Astfel, câteva cifre de la intrarea (i ) sunt lăsate în cele de sumă .

Suma este atât de mare încât se acumulează doar cele mai semnificative cifre ale numerelor introduse. Dar, la pasul următor, c dă o eroare.

 y = 2.71828 - -.0415900 Deficitul obținut în etapa anterioară este inclus.
    = 2.75987 Are dimensiuni similare cu y : sunt lăsate cele mai semnificative cifre.
  t = 10003,1 + 2,75987 Dar puțini au rămas în sumă
    = 10005.85987 și rezultatul este rotunjit
    = 10005,9 șase cifre.
  c = (10005,9 - 10003,1) - 2,75987 Aceasta extrage orice a fost obținut.
    = 2.80000 - 2.75987 În acest caz, există un exces.
    = .040130 Dar nu contează, excesul ar fi scăzut data viitoare.
sum = 10005.9 Rezultatul exact este 10005.85987, este rotunjit corect la 6 cifre.

În acest fel, suma se acumulează folosind două variabile: suma care se referă la suma și c unde se acumulează părțile neasimilate în sumă , pentru a muta din nou partea mai puțin semnificativă din sumă data viitoare. Apoi însumarea continuă cu „cifre de gardă” în c, care este mai bine decât nimic, dar care nu este la fel de bun precum ar fi efectuarea calculelor cu precizie dublă a intrării. Cu toate acestea, pur și simplu creșterea preciziei calculelor nu este în general practică; dacă intrarea este deja cu precizie dublă, puține sisteme oferă precizie cvadruplă, iar dacă o fac, intrarea ar putea fi, prin urmare, cu precizie cvadruplă.

Precizie

Este necesară o analiză atentă a erorilor din însumarea compensată pentru a aprecia caracteristicile sale de precizie. Deși este mai precis decât o simplă însumare, poate da totuși erori relative mari în cazul sumelor în condiții deosebit de proaste.

Să presupunem că adăugăm n valori x i , pentru i = 1, ..., n . Suma exactă este:

(calculat cu o precizie infinită)

Pe de altă parte, cu însumarea compensată, obținem , unde eroarea este limitat de: [2]

unde ε este precizia mașinii aritmeticii utilizate (de exemplu, ε≈10 −16 pentru formatul cu virgulă mobilă de dublă precizie conform standardului IEEE). În general, valoarea dobânzii este eroarea relativă , care este, prin urmare, limitat în partea de sus de:

În expresia limitei relative de eroare, fracția Σ | x i | / | Σ x i | este condiționarea problemei însumării. În esență, condiționarea reprezintă sensibilitatea inerentă la erori ale problemei de însumare, indiferent de modul în care este calculată. [6] Limita relativă de eroare a fiecărei metode de sumare ( stabilă înapoi ) bazată pe un algoritm de precizie fixă ​​(adică nu cele care utilizează aritmetică de precizie arbitrară, nici algoritmi ale căror cerințe de memorie și timp se modifică pe baza datelor), este proporțională cu această condiționare . [2] O problemă de însumare slab condiționată este una în care acest raport este mare și, în acest caz, chiar și sumarea compensată poate avea o eroare relativă mare. De exemplu, dacă aditivele x i sunt numere aleatorii fără legătură cu medie zero, suma este o plimbare aleatorie și condiționarea va crește proporțional cu . Pe de altă parte, pentru intrări aleatorii medii diferite de zero, condiționarea tinde la o constantă finită. Dacă intrările sunt toate non-negative, atunci condiționarea este 1.

Având în vedere o anumită valoare pentru condiționare, eroarea relativă a însumării compensate este efectiv independentă de n . În principiu, există O ( n ε 2 ) care crește liniar cu n , dar în practică acest termen este efectiv zero: deoarece rezultatul final este rotunjit la o precizie ε, termenul n ε 2 tinde la zero, cu excepția cazului în care n este aproximativ 1 / ε sau mai mare. [2] În dublă precizie, aceasta corespunde unui n de aproximativ 10 16 , mult mai mare decât majoritatea sumelor. Prin urmare, pentru o valoare fixă ​​a condiționării, erorile însumării compensate sunt efectiv O (ε), indiferent de n .

În comparație, limita relativă de eroare pentru însumarea simplă (pur și simplu adăugând numerele în ordine, rotunjind la fiecare pas) crește ca înmulțit cu condiționarea. [2] Această cea mai gravă eroare posibilă, însă, este rar observată în practică, deoarece apare doar dacă erorile de rotunjire sunt toate în aceeași direcție. În practică, erorile de rotunjire sunt mult mai susceptibile de a avea un semn aleatoriu, cu medie zero, în așa fel încât să producă o mers aleatoriu; în acest caz, suma simplă are o abatere standard, referindu-se la erorile relative, care crește ca înmulțit cu condiționarea. [7] Cu toate acestea, acest lucru este încă mult mai rău decât suma compensată. Rețineți, totuși, că dacă suma poate fi realizată cu precizie dublă, atunci ε este înlocuită cu ε 2 și suma simplă are o eroare mai gravă posibilă comparabilă cu termenul O ( n ε 2 ) în suma compensată cu precizia inițială.

Din același motiv, termenul Σ | x i | care apare mai sus, în , este cea mai proastă limită posibilă care apare numai dacă toate erorile de rotunjire au același semn (și sunt de mărimea maximă posibilă). [2] În practică, este mai probabil ca erorile să aibă un semn aleatoriu, caz în care termenii din Σ | x i | sunt înlocuite de o plimbare aleatorie; în acest caz, chiar și pentru zero intrări aleatorii medii, eroarea doar crește ca. (ignorând termenul n ε 2 ), suma crește cu aceeași viteză , ștergerea factorilor când se calculează eroarea relativă. Astfel, chiar și pentru sumele necondiționate asimptotic, eroarea relativă pentru însumarea compensată poate fi adesea mult mai mică decât ar sugera o analiză în cel mai rău caz.

Îmbunătățiri suplimentare

Neumaier [8] a introdus o ușoară modificare a algoritmului lui Kahan, care acoperă și cazul în care următorul termen care trebuie adăugat este mai mare în valoare absolută decât suma curentă, schimbând efectiv rolul a ceea ce este mare și a ceea ce este mic. În pseudocod , algoritmul este:

 funcția NeumaierSum (intrare)
    var sum = input [1]
    var c = 0,0 // O compensare în curs pentru biții cei mai puțin semnificativi pierduți.
    pentru i = 2 la input.length do
        var t = sumă + intrare [i]
        dacă | suma | > = | intrare [i] | do
            c + = (sumă - t) + intrare [i] // Dacă suma este mai mare, cele mai puțin semnificative cifre ale intrării [i] se pierd.
        altceva
            c + = (input [i] - t) + sum // În caz contrar, se pierd cifrele de sumă mai puțin semnificative.
        sumă = t
    return sum + c // Corecție aplicată o singură dată la final.

Pentru multe secvențe numerice, ambii algoritmi sunt de acord, dar un exemplu simplu al lui Peters [9] arată că pot diferi. Se adaugă în dublă precizie, algoritmul lui Kahan produce 0,0, în timp ce algoritmul lui Neumaier dă valoarea corectă 2.0.

Alternative

Deși cu algoritmul lui Kahan, prin adăugarea de n numere, o creștere a erorii de ordinul lui , doar o creștere ușor mai slabă în ordine se poate realiza printr-o însumare recursivă în perechi : setul de numere este împărțit recursiv în două jumătăți, adăugând fiecare jumătate și apoi se adaugă cele două sume. [2] Aceasta are avantajul de a solicita același număr de operații aritmetice ca sumarea simplă (spre deosebire de algoritmul lui Kahan, care, în comparație cu suma simplă, necesită un număr de operații aritmetice și are o latență de patru ori mai mare), iar calculul poate fi efectuat în paralel. Modelul de bază al recursivității ar putea fi, în principiu, suma doar a unui număr (sau a numerelor zero), dar pentru a conține cheltuielile generale ale recursivității, se folosește în mod normal un model de bază mai larg. Echivalentul însumării recursive perechi este utilizat în mulți algoritmi de transformare Fourier rapidă (FFT) și este responsabil pentru creșterea logaritmică a erorilor de rotunjire în astfel de FFT-uri. [10] În practică, cu erori de rotunjire ale semnului aleator, abaterile standard ale însumării recursive perechi cresc de fapt ca . [7]

O altă alternativă este utilizarea aritmeticii de precizie arbitrară , care, în principiu, nu necesită deloc rotunjire cu costul unui efort de calcul mult mai mare. O modalitate de a face sume perfect rotunjite folosind precizie arbitrară este de a le extinde în mod adaptiv folosind mai multe componente în virgulă mobilă. Acest lucru va reduce costul de calcul în cazurile obișnuite în care nu este necesară o precizie ridicată. [11] [9] O altă metodă care folosește doar aritmetica întreagă, dar un acumulator mare a fost descrisă de Kirchner și Kulisch; [12] O implementare hardware a fost descrisă de Müller, Rüb și Rülling. [13]

Posibilă invalidare din optimizarea compilatorului

În principiu, o optimizare suficient de agresivă a compilatorului ar putea compromite eficacitatea însumării Kahan: de exemplu, dacă compilatorul are expresii simplificate în conformitate cu regulile de asociativitate ale aritmeticii reale, ar putea simplifica al doilea pas al secvenței din t = sum + y; c = (t - sum) - y; a ((sum + y) - sum) - y; apoi la c = 0; , eliminând compensarea erorilor. [14] În practică, mulți compilatori nu folosesc reguli de asociativitate (care, în aritmetică în virgulă mobilă, sunt doar aproximative) în simplificări, cu excepția cazului în care li se cere în mod explicit să facă acest lucru prin opțiuni de compilare care permit optimizări nesigure, [15] [16] [17 ] ] [18] deși compilatorul Intel C ++ este un exemplu care permite implicit transformări bazate pe asociativitate. [19] Versiunea originală K&R C a limbajului C a permis compilatorului să reordineze expresii în virgulă mobilă bazate pe regulile de asociativitate ale aritmeticii reale, dar ulterior standardul ANSI C nu a permis reordonarea pentru a face limbajul C mai potrivit pentru numerică aplicații (și mai mult ca Fortran , care nu permite reordonarea), [20] deși în practică opțiunile compilatorului pot reactiva reordonarea așa cum s-a menționat mai sus.

Suport cu biblioteci

În general, funcțiile de sumă incluse în limbajele de programare nu oferă de obicei nicio garanție că va fi utilizat un anumit algoritm de sumare, cum ar fi sumarea Kahan. [ Citație necesară ] BLAS standard (Basic Linear Algebra Subprograms) pentru subrutine de algebră liniară evită în mod explicit impunerea unei secvențe particulare de calcul pentru operațiuni din motive de performanță, [21] și implementările BLAS, de obicei, nu folosesc însumarea Kahan.

Biblioteca standard de limbaj Python specifică o funcție fsum pentru însumări exact rotunjite, utilizând algoritmul lui Shewchuk [9] pentru a trasa mai multe sume care rulează.

În limbajul Julia , implementarea implicită a funcției sum utilizează sumare recursivă pereche pentru o precizie ridicată cu performanțe bune [22], dar biblioteca standard are și o variantă de implementare sum_kbn numită sum_kbn pentru cazurile în care este necesară o precizie mai mare. [23]

Notă

  1. ^ Strict, există și alte variații ale însumării compensate: vezi Nicholas Higham, Accuracy and Stability of Numerical Algorithms (2 ed.) , SIAM, 2002, pp. 110–123.
  2. ^ a b c d e f g h Nicholas J. Higham, Precizia sumării în virgulă mobilă , în SIAM Journal on Scientific Computing , vol. 14, n. 4, 1993, pp. 783–799, DOI : 10.1137 / 0914050 .
  3. ^ William Kahan, Remarci suplimentare despre reducerea erorilor de trunchiere , în Comunicări ale ACM , vol. 8, nr. 1, ianuarie 1965, p. 40, DOI : 10.1145 / 363707.363723 .
  4. ^ Jack E. Bresenham, „Algorithm for computer control of a digital plotter” ( PDF ), în IBM Systems Journal , vol. 4, nr. 1, ianuarie 1965, pp. 25-30. Adus pe 29 aprilie 2019 (depus de „url original 28 mai 2008).
  5. ^ H. Inose, Y. Yasuda, J. Murakami, "A Telemetering System by Code Manipulation - ΔΣ Modulation", IRE Trans on Space Electronics and Telemetry, sept. 1962, pp. 204-209.
  6. ^ LN Trefethen și D. Bau, Algebra liniară numerică (SIAM: Philadelphia, 1997).
  7. ^ a b Manfred Tasche și Hansmartin Zeuner Handbook of Analytic-Computational Methods in Applied Mathematics Boca Raton, FL: CRC Press, 2000).
  8. ^ ( DE ) A. Neumaier, Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen , în Zeitschrift für Angewandte Mathematik und Mechanik , vol. 54, nr. 1, 1974, pp. 39-51.
  9. ^ a b c Raymond Hettinger, Rețeta 393090: Sumare în virgulă mobilă binară precisă până la precizie maximă , implementare Python a algoritmului din lucrarea Shewchuk (1997) (28 martie 2005).
  10. ^ SG Johnson și M. Frigo, Implementing FFTs in practice , în C. Sidney Burrus (ed.), Fast Fourier Transforms , 2008.
  11. ^ Jonathan R. Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates , Discrete and Computational Geometry , vol. 18, pp. 305–363 (octombrie 1997).
  12. ^ R. Kirchner, UW Kulisch, Accurate arithmetic for vector processors , Journal of Parallel and Distributed Computing 5 (1988) 250-270
  13. ^ M. Muller, C. Rub, W. Rulling [1] , Acumularea exactă a numerelor în virgulă mobilă , Proceedings 10th IEEE Symposium on Computer Arithmetic (iunie 1991), doi 10.1109 / ARITH.1991.145535
  14. ^ David Goldberg, Ce ar trebui să știe fiecare informatician despre aritmetica în virgulă mobilă ( PDF ), în ACM Computing Surveys , vol. 23, n. 1, martie 1991, pp. 5–48, DOI : 10.1145 / 103162.103163 .
  15. ^ GNU Compiler Collection manual, versiunea 4.4.3: 3.10 Opțiuni care controlează optimizarea , -fassociative-math (21 ianuarie 2010).
  16. ^ Compaq Fortran Manual de utilizare pentru Tru64 UNIX și Linux Alpha Systems Arhivat 7 iunie 2011 la Internet Archive . , secțiunea 5.9.7 Optimizări de reordonare aritmetică (recuperat în martie 2010).
  17. ^ Börje Lindh, Application Performance Optimization. Arhivat 2 iunie 2010 la Internet Archive ., Sun BluePrints OnLine (martie 2002).
  18. ^ Eric Fleegal, „ Microsoft Visual C ++ Floating-Point Optimization, articolele tehnice Microsoft Visual Studio (iunie 2004).
  19. ^ Martyn J. Corden, „ Coerența rezultatelor în virgulă mobilă utilizând compilatorul Intel ”, raport tehnic Intel (18 septembrie 2009).
  20. ^ Tom Macdonald, "C for Numerical Computing", Journal of Supercomputing vol. 5, pp. 31-48 (1991).
  21. ^ BLAS Technical Forum , secțiunea 2.7 (21 august 2001), arhivat pe Wayback Machine .
  22. ^ https://docs.julialang.org/en/stable/stdlib/collections/?highlight=sum#Base.sum
  23. ^ https://docs.julialang.org/en/stable/stdlib/arrays/?highlight=sum_kbn#Base.sum_kbn

Bibliografie

  • Luciano M. Barone, Enzo Marinari, Giovanni Organtini și Federico Ricci-Terseng, Programare științifică. Limbajul C, algoritmi și modele în știință , Pearson Education Italia, 2006, p. 108, paragraful 4.5: O problemă de rotunjire , ISBN 978-88-7192-242-3 .

Elemente conexe

linkuri externe