Algorytm sumowania Kahana - Kahan summation algorithm

W analizie numerycznej The algorytm sumowanie Kahan , znany również jako wyrównaną sumowania , co znacznie zmniejsza błąd liczbową w całości otrzymanego przez dodanie sekwencji z finite- precyzyjnych liczb zmiennoprzecinkowych w porównaniu z oczywistym podejściem. Odbywa się to poprzez zachowanie oddzielnej kompensacji bieżącej (zmiennej do akumulacji małych błędów).

W szczególności, po prostu zsumowanie liczb w kolejności ma błąd najgorszy, który rośnie proporcjonalnie do , a średni kwadratowy błąd, który rośnie jak na przypadkowych wejść (na błędy zaokrąglenia tworzą przypadkowy spacer ). W przypadku sumowania skompensowanego granica najgorszego przypadku błędu jest skutecznie niezależna od , więc duża liczba wartości może być sumowana z błędem zależnym tylko od precyzji zmiennoprzecinkowej .

Algorytm przypisuje się William Kahan ; Wydaje się, że Ivo Babuška samodzielnie wymyślił podobny algorytm (stąd suma Kahana–Babuški ). Podobne, wcześniejsze techniki to na przykład algorytm liniowy Bresenhama , śledzący skumulowany błąd w operacjach na liczbach całkowitych (choć najpierw udokumentowany w tym samym czasie) oraz modulacja delta-sigma

Algorytm

W pseudokodzie algorytmem będzie;

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

Przykład pracy

Ten przykład zostanie podany w postaci dziesiętnej. Komputery zazwyczaj używają arytmetyki binarnej, ale ilustrowana zasada jest taka sama. Załóżmy, że używamy sześciocyfrowej dziesiętnej arytmetyki zmiennoprzecinkowej, sumosiągnęła wartość 10000,0, a kolejne dwie wartości input[i]to 3,14159 i 2,71828. Dokładny wynik to 10005.85987, który zaokrągla się do 10005,9. W przypadku zwykłego sumowania każda przychodząca wartość byłaby wyrównana z sum, a wiele cyfr niższego rzędu zostałoby utraconych (przez obcięcie lub zaokrąglenie). Pierwszy wynik po zaokrągleniu to 10003.1. Drugi wynik to 10005.81828 przed zaokrągleniem i 10005.8 po zaokrągleniu. To nie jest poprawne.

Jednak przy sumowaniu skompensowanym otrzymujemy poprawny zaokrąglony wynik 10005.9.

Załóżmy, że cma wartość początkową zero.

  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.

Suma jest tak duża, że ​​gromadzone są tylko cyfry wyższego rzędu liczb wejściowych. Ale w następnym kroku cpodaje błąd.

  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.

Tak więc sumowanie odbywa się za pomocą dwóch akumulatorów: sumprzechowuje sumę i cakumuluje części, które nie są zasymilowane w sum, aby przesunąć część niższego rzędu sumprzy następnym cyklu. W ten sposób sumowanie przebiega z „cyframi ochronnymi” w c, co jest lepsze niż ich brak, ale nie jest tak dobre, jak wykonywanie obliczeń z podwójną dokładnością danych wejściowych. Jednak samo zwiększenie precyzji obliczeń nie jest ogólnie praktyczne; jeśli inputjuż jest w podwójnej precyzji, niewiele systemów zapewnia poczwórną precyzję , a gdyby tak było, to inputmogłoby być w poczwórnej precyzji.

Precyzja

Konieczna jest dokładna analiza błędów w sumowaniu skompensowanym, aby docenić jego charakterystykę dokładności. Chociaż jest dokładniejsze niż sumowanie naiwne, nadal może dawać duże błędy względne dla źle uwarunkowanych sum.

Załóżmy, że sumujemy wartości , dla . Dokładna suma to

(obliczane z nieskończoną precyzją).

Przy sumowaniu skompensowanym otrzymujemy , gdzie błąd jest ograniczony przez

gdzie jest precyzja maszynowa stosowanej arytmetyki (np. dla standardu IEEE zmiennoprzecinkowego podwójnej precyzji ). Zwykle wielkość zainteresowania jest błędem względnym , który jest zatem ograniczony powyżej przez

W wyrażeniu na względną granicę błędu ułamek jest numerem warunku problemu sumowania. Zasadniczo numer warunku reprezentuje wewnętrzną wrażliwość problemu sumowania na błędy, niezależnie od sposobu jego obliczenia. Względna granica błędu każdej ( stabilnej wstecz ) metody sumowania przez stały algorytm o stałej precyzji (tj. nie te, które używają arytmetyki o dowolnej precyzji , ani algorytmów, których wymagania dotyczące pamięci i czasu zmieniają się na podstawie danych), jest proporcjonalna do tego warunku liczba . Źle uwarunkowane problemem jest podsumowanie, w którym wskaźnik ten jest duży, aw tym przypadku nawet skompensowane podsumowanie może mieć duży błąd względny. Na przykład, jeśli summands są nieskorelowanymi liczbami losowymi o zerowej średniej, suma jest błądzeniem losowym , a liczba warunku będzie rosła proporcjonalnie do . Z drugiej strony, dla losowych danych wejściowych o wartości niezerowej oznaczają asymptoty liczby warunku do skończonej stałej jako . Jeśli wszystkie dane wejściowe są nieujemne , numer warunku wynosi 1.

Biorąc pod uwagę numer warunku, względny błąd skompensowanego sumowania jest skutecznie niezależny od . W zasadzie istnieje ten, który rośnie liniowo z , ale w praktyce ten termin jest faktycznie zero: ponieważ końcowy wynik jest zaokrąglany do dokładności , termin zaokrągla się do zera, chyba że jest w przybliżeniu lub większy. W podwójnej precyzji odpowiada to z grubsza , znacznie większe niż większość sum. Tak więc, dla ustalonej liczby warunku, błędy skompensowanego sumowania są efektywnie niezależne od .

Dla porównania, względny błąd związany z naiwnym sumowaniem (po prostu dodawanie liczb w kolejności, zaokrąglanie na każdym kroku) rośnie jako pomnożony przez liczbę warunku. Ten najgorszy przypadek jest jednak rzadko obserwowany w praktyce, ponieważ występuje tylko wtedy, gdy wszystkie błędy zaokrągleń są w tym samym kierunku. W praktyce znacznie bardziej prawdopodobne jest, że błędy zaokrągleń mają znak losowy, ze średnią zerową, tak że tworzą błądzenie losowe; w tym przypadku naiwne sumowanie ma pierwiastek średniokwadratowy względny błąd, który rośnie jako pomnożony przez liczbę warunku. Jest to jednak nadal znacznie gorsze niż sumowanie wyrównane. Jeśli jednak suma może być wykonana z podwójną precyzją, to jest zastępowana przez , a sumowanie naiwne ma błąd najgorszego przypadku porównywalny z terminem sumowania skompensowanego z pierwotną precyzją.

Z tego samego powodu, to, co pojawia się powyżej, jest ograniczeniem najgorszego przypadku, które występuje tylko wtedy, gdy wszystkie błędy zaokrągleń mają ten sam znak (i ​​mają maksymalną możliwą wielkość). W praktyce bardziej prawdopodobne jest, że błędy mają znak losowy, w którym to przypadku wyrazy w są zastępowane błądzeniem losowym, w którym to przypadku nawet dla losowych danych wejściowych o średniej zerowej błąd rośnie tylko wtedy, gdy (ignorując wyraz), w tym samym tempie suma rośnie, anulując czynniki, gdy obliczany jest błąd względny. Tak więc, nawet w przypadku sum asymptotycznie źle uwarunkowanych, względny błąd dla sumowania skompensowanego może często być znacznie mniejszy niż może to sugerować analiza najgorszego przypadku.

Dalsze ulepszenia

Neumaier wprowadził ulepszoną wersję algorytmu Kahana, który nazywa „ulepszonym algorytmem Kahana-Babuški”, który obejmuje również przypadek, gdy następny dodany termin ma wartość bezwzględną większą niż suma bieżąca, skutecznie zamieniając rolę tego, co jest duże i małe. W pseudokodzie algorytm to:

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.

W przypadku wielu ciągów liczb oba algorytmy są zgodne, ale prosty przykład Petersa pokazuje, jak mogą się różnić. W przypadku sumowania z podwójną precyzją algorytm Kahana daje 0,0, podczas gdy algorytm Neumaiera poprawną wartość 2,0.

Możliwe są również modyfikacje wyższego rzędu o lepszej dokładności. Na przykład wariant zaproponowany przez Kleina, który nazwał „iteracyjnym algorytmem Kahana-Babuški” drugiego rzędu. W pseudokodzie algorytm to:

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

Alternatywy

Chociaż algorytm Kahana osiąga wzrost błędu przy sumowaniu n liczb, tylko nieco gorszy wzrost można osiągnąć przez sumowanie parami : jeden rekurencyjnie dzieli zbiór liczb na dwie połowy, sumuje każdą połowę, a następnie dodaje dwie sumy. Ma to tę zaletę, że wymaga takiej samej liczby operacji arytmetycznych jak sumowanie naiwne (w przeciwieństwie do algorytmu Kahana, który wymaga czterokrotności arytmetyki i ma opóźnienie czterokrotnie większe od prostego sumowania) i może być obliczane równolegle. Podstawowy przypadek rekurencji może w zasadzie być sumą tylko jednej (lub zerowej) liczby, ale aby zamortyzować narzut rekurencji, zwykle używa się większego przypadku podstawowego. Odpowiednik sumowania parami jest używany w wielu algorytmach szybkiej transformacji Fouriera (FFT) i jest odpowiedzialny za logarytmiczny wzrost błędów zaokrągleń w tych FFT. W praktyce, przy błędach zaokrąglania znaków losowych, pierwiastek błędów średniokwadratowych sumowania parami faktycznie rośnie jako .

Inną alternatywą jest użycie arytmetyki o dowolnej precyzji , która w zasadzie nie wymaga zaokrąglania, co wiąże się ze znacznie większym wysiłkiem obliczeniowym. Sposobem na wykonanie dokładnie zaokrąglonych sum przy użyciu dowolnej precyzji jest adaptacyjne rozszerzanie przy użyciu wielu składników zmiennoprzecinkowych. Zminimalizuje to koszty obliczeniowe w typowych przypadkach, w których nie jest potrzebna wysoka precyzja. Inna metoda, która wykorzystuje tylko arytmetykę liczb całkowitych, ale duży akumulator, została opisana przez Kirchnera i Kulischa ; implementację sprzętową opisali Müller, Rüb i Rülling.

Możliwe unieważnienie przez optymalizację kompilatora

W zasadzie wystarczająco agresywny kompilator optymalizujący mógłby zniszczyć skuteczność sumowania Kahana: na przykład, gdyby kompilator uprościł wyrażenia zgodnie z regułami asocjatywności prawdziwej arytmetyki, mógłby „uprościć” drugi krok w sekwencji

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

do

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

a potem do

c = 0;

w ten sposób eliminując kompensację błędów. W praktyce wiele kompilatorów nie używa reguł asocjacji (które są tylko przybliżone w arytmetyce zmiennoprzecinkowej) w uproszczeniu, chyba że wyraźnie nakazano to za pomocą opcji kompilatora umożliwiających „niebezpieczne” optymalizacje, chociaż kompilator Intel C++ jest jednym z przykładów umożliwiających asocjatywność przekształcenia oparte domyślnie. Oryginalna wersja K&R C języka programowania C pozwalała kompilatorowi na zmianę kolejności wyrażeń zmiennoprzecinkowych zgodnie z regułami asocjatywności rzeczywistej arytmetyki, ale późniejszy standard ANSI C zabronił zmiany kolejności w celu lepszego dopasowania C do aplikacji numerycznych ( i bardziej podobny do Fortran , który również zabrania ponownego zamawiania), chociaż w praktyce opcje kompilatora mogą ponownie włączyć ponowne zamawianie, jak wspomniano powyżej.

Przenośnym sposobem na lokalne powstrzymanie takich optymalizacji jest przełamanie jednej z linii w oryginalnym sformułowaniu na dwa zdania i uczynienie dwóch produktów pośrednich niestabilnymi :

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

Wsparcie przez biblioteki

Ogólnie rzecz biorąc, wbudowane funkcje „sumowania” w językach komputerowych zazwyczaj nie dają gwarancji, że zostanie zastosowany konkretny algorytm sumowania, a tym bardziej sumowanie Kahana. Standard BLAS dla podprogramów algebry liniowej wyraźnie unika nakazu jakiejkolwiek określonej kolejności obliczeniowej operacji ze względu na wydajność, a implementacje BLAS zwykle nie używają sumowania Kahana.

Standardowa biblioteka języka komputerowego Python określa funkcję fsum dla dokładnie zaokrąglonych sum, używając algorytmu Shewchuk do śledzenia wielu sum częściowych.

W języku Julia domyślna implementacja sumfunkcji wykonuje sumowanie parami w celu uzyskania wysokiej dokładności z dobrą wydajnością, ale zewnętrzna biblioteka zapewnia implementację wariantu Neumaiera nazwanego sum_kbndla przypadków, gdy potrzebna jest wyższa dokładność.

W języku C# pakiet nuget HPCsharp implementuje wariant Neumaiera i sumowanie parami : zarówno jako skalarne, równoległe do danych przy użyciu instrukcji procesora SIMD , jak i równoległe wielordzeniowe.

Zobacz też

Bibliografia

Zewnętrzne linki