Algoritmus souhrnu Kahana - Kahan summation algorithm

V numerické analýze je algoritmus Kahan sčítání , také známý jako kompenzovaným shrnutí , výrazně snižuje numerické chyby v součtu dosaženo přídavkem sekvence z finite- přesných čísel s plovoucí čárkou , ve srovnání se zřetelnými přístupu. To se provádí udržováním oddělené kompenzace běhu (proměnná pro akumulaci malých chyb).

Zejména jednoduché sečtení čísel v pořadí má chybu nejhoršího případu, která roste úměrně , a odmocninu střední chyby, která roste jako u náhodných vstupů (chyby zaokrouhlení tvoří náhodnou procházku ). Při kompenzovaném součtu je vazba chyby nejhoršího případu skutečně nezávislá na , takže velký počet hodnot lze sčítat s chybou, která závisí pouze na přesnosti s plovoucí desetinnou čárkou . ${\ displaystyle n}$${\ displaystyle n}$${\ displaystyle {\ sqrt {n}}}$${\ displaystyle n}$

Algoritmus je přičítán Williama Kahan ; Zdá se, že Ivo Babuška přišel s podobným algoritmem nezávisle (proto součet Kahan – Babuška ). Podobné, dřívější techniky jsou například Bresenhamův liniový algoritmus , sledování akumulované chyby v celočíselných operacích (i když nejprve zdokumentovány přibližně ve stejnou dobu) a modulace delta-sigma

Algoritmus

V pseudokódu bude algoritmus;

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

Fungoval příklad

Tento příklad bude uveden v desítkové soustavě. Počítače obvykle používají binární aritmetiku, ale princip, který je znázorněn, je stejný. Předpokládejme, že používáme šestimístnou desítkovou aritmetiku s plovoucí desetinnou čárkou, sumdosáhla hodnoty 10 000,0 a další dvě hodnoty input[i]jsou 3,14159 a 2,71828. Přesný výsledek je 10005,85987, který se zaokrouhlí na 10005,9. S prostým součtem by každá příchozí hodnota byla zarovnána suma mnoho číslic nízkého řádu by bylo ztraceno (zkrácením nebo zaokrouhlením). První výsledek po zaokrouhlení bude 10003,1. Druhý výsledek by byl 10005,81828 před zaokrouhlením a 10005,8 po zaokrouhlení. To není správné.

S kompenzovaným součtem však dostaneme správný zaokrouhlený výsledek 10005,9.

Předpokládejme, že cmá počáteční hodnotu nula.

  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.

Součet je tak velký, že se shromažďují pouze číslice vyšších řádů vstupních čísel. Ale v dalším kroku cdává chybu.

  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.

Sčítání se tedy provádí pomocí dvou akumulátorů: sumdrží součet a chromadí části, které do něj nejsou asimilovány sum, aby sumpříště posunul část nižšího řádu . Součet tedy pokračuje „ochrannými číslicemi“ c, což je lepší než nemít žádné, ale není tak dobré, jako provádět výpočty s dvojnásobnou přesností vstupu. Pouhé zvýšení přesnosti výpočtů však není obecně praktické; pokud inputje již v dvojnásobné přesnosti, několik systémů poskytuje čtyřnásobnou přesnost , a pokud ano, inputpak by mohly být ve čtyřnásobné přesnosti.

Přesnost

K zhodnocení charakteristik přesnosti je nutná pečlivá analýza chyb v kompenzovaném součtu. I když je přesnější než naivní součet, stále může u špatně podmíněných částek dávat velké relativní chyby.

Předpokládejme, že jedna sčítá hodnoty , pro . Přesná částka je ${\ displaystyle n}$${\ displaystyle x_ {i}}$${\ Displaystyle i = 1, \, \ ldots, \, n}$

${\ Displaystyle S_ {n} = \ sum _ {i = 1}^{n} x_ {i}}$ (počítáno s nekonečnou přesností).

Při kompenzovaném součtu se místo toho získá , kde je chyba ohraničena ${\ displaystyle S_ {n}+E_ {n}}$${\ displaystyle E_ {n}}$

${\ Displaystyle | E_ {n} | \ leq {\ big [} 2 \ varepsilon +O (n \ varepsilon ^{2}) {\ big]} \ sum _ {i = 1} ^{n} | x_ { i} |,}$

kde je strojní přesnost použité aritmetiky (např. u standardu IEEE s plovoucí desetinnou čárkou s dvojitou přesností ). Obvykle je úroková částka relativní chybou , která je proto ohraničena výše ${\ displaystyle \ varepsilon}$${\ Displaystyle \ varepsilon \ cca 10^{-16}}$ ${\ displaystyle | E_ {n} |/| S_ {n} |}$

${\ displaystyle {\ frac {| E_ {n} |} {| S_ {n} |}} \ leq {\ big [} 2 \ varepsilon +O (n \ varepsilon ^{2}) {\ big]} { \ frac {\ sum \ limits _ {i = 1}^{n} | x_ {i} |} {\ left | \ sum \ limits _ {i = 1}^{n} x_ {i} \ right |} }.}$

Ve výrazu pro relativní chyby vázané frakce je číslo stav problému součtového. V zásadě číslo podmínky představuje vnitřní citlivost problému součtu na chyby, bez ohledu na to, jak je vypočítán. Relativní chyba svázaná každou ( zpětně stabilní ) metodou součtu fixním algoritmem s pevnou přesností (tj. Ne těmi, které používají aritmetiku s libovolnou přesností , ani algoritmy, jejichž požadavky na paměť a čas se mění na základě dat), je úměrná tomuto počtu podmínek . Špatně podmíněné součet Problém je ten, ve kterém je tento poměr je velký, a v tomto případě i kompenzován součtem mohou mít velký relativní chybu. Pokud jsou například součty nekorelovaná náhodná čísla s nulovým průměrem, součet je náhodná procházka a číslo podmínky poroste úměrně . Na druhou stranu pro náhodné vstupy s nenulovou hodnotou znamená počet podmínek asymptoty na konečnou konstantu jako . Pokud jsou všechny vstupy nezáporné , pak je číslo podmínky 1. ${\ displaystyle \ Sigma | x_ {i} |/| \ Sigma x_ {i} |}$${\ displaystyle x_ {i}}$${\ displaystyle {\ sqrt {n}}}$${\ displaystyle n \ to \ infty}$

Vzhledem k počtu podmínek je relativní chyba kompenzovaného součtu skutečně nezávislá na . V zásadě existuje lineárně rostoucí , ale v praxi je tento termín ve skutečnosti nulový: protože konečný výsledek je zaokrouhlen na přesnost , termín se zaokrouhlí na nulu, pokud není zhruba nebo větší. Ve dvojité přesnosti to odpovídá zhruba , mnohem větším částkám. Takže pro číslo pevné podmínky jsou chyby kompenzovaného součtu efektivně nezávislé na . ${\ displaystyle n}$${\ Displaystyle O (n \ varepsilon ^{2})}$${\ displaystyle n}$${\ displaystyle \ varepsilon}$${\ Displaystyle n \ varepsilon ^{2}}$${\ displaystyle n}$${\ Displaystyle 1/\ varepsilon}$${\ displaystyle n}$${\ displaystyle 10^{16}}$${\ Displaystyle O (\ varepsilon)}$${\ displaystyle n}$

Pro srovnání relativní chyba vázaná na naivní součet (prosté sečtení čísel za sebou, zaokrouhlení v každém kroku) roste vynásobením číslem podmínky. Tato chyba nejhoršího případu je však v praxi jen zřídka pozorována, protože k ní dochází pouze v případě, že jsou chyby zaokrouhlení všechny ve stejném směru. V praxi je mnohem pravděpodobnější, že chyby zaokrouhlení mají náhodné znaménko s nulovým průměrem, takže tvoří náhodnou chůzi; v tomto případě má naivní součet kořenovou střední kvadratickou relativní chybu, která roste vynásobením číslem podmínky. To je však stále mnohem horší než kompenzované součty. Pokud však součet lze provést s dvojnásobnou přesností, pak je nahrazen a naivní součet má nejhorší možnou chybu srovnatelnou s výrazem v kompenzovaném součtu při původní přesnosti. ${\ Displaystyle O (\ varepsilon n)}$${\ Displaystyle O \ left (\ varepsilon {\ sqrt {n}} \ right)}$${\ displaystyle \ varepsilon}$${\ displaystyle \ varepsilon ^{2}}$${\ Displaystyle O (n \ varepsilon ^{2})}$

Ze stejného důvodu platí, že to, co se objeví výše, je vazba pro nejhorší případy, ke které dochází pouze v případě, že všechny chyby zaokrouhlení mají stejné znaménko (a mají maximální možnou velikost). V praxi je pravděpodobnější, že chyby mají náhodné znaménko, přičemž v tomto případě jsou výrazy v nahrazeny náhodnou chůzí, v takovém případě i u náhodných vstupů s nulovým průměrem chyba roste pouze tak, že (ignorování výrazu) součtem roste součet , což ruší faktory při výpočtu relativní chyby. Takže i pro asymptoticky špatně podmíněné částky může být relativní chyba pro kompenzované součty často mnohem menší, než by mohla naznačovat analýza nejhoršího případu. ${\ displaystyle \ Sigma | x_ {i} |}$${\ displaystyle E_ {n}}$${\ displaystyle \ Sigma | x_ {i} |}$${\ displaystyle E_ {n}}$${\ Displaystyle O \ left (\ varepsilon {\ sqrt {n}} \ right)}$${\ Displaystyle n \ varepsilon ^{2}}$${\ displaystyle S_ {n}}$${\ displaystyle {\ sqrt {n}}}$

Další vylepšení

Neumaier představil vylepšenou verzi Kahanova algoritmu, který nazývá „vylepšený Kahan – Babuškův algoritmus“, který také pokrývá případ, kdy je další přidávaný výraz větší v absolutní hodnotě než je běžná částka, čímž se efektivně vymění role toho, co je velký a co je malý. V pseudokódu je algoritmus:

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.

U mnoha sekvencí čísel oba algoritmy souhlasí, ale jednoduchý příklad kvůli Petersovi ukazuje, jak se mohou lišit. Pro sčítání s dvojitou přesností přináší Kahanův algoritmus 0,0, zatímco Neumaierův algoritmus dává správnou hodnotu 2,0. ${\ displaystyle [1,0,+10^{100}, 1,0, -10^{100}]}$

Jsou možné i úpravy vyššího řádu s lepší přesností. Například varianta navržená Kleinem, kterou nazval „iteračním algoritmem Kahan – Babuška druhého řádu“. V pseudokódu je algoritmus:

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

Alternativy

Ačkoli Kahanův algoritmus dosahuje nárůstu chyb při součtu n čísel, jen o málo horšího růstu lze dosáhnout párovým součtem : jeden rekurzivně rozdělí sadu čísel na dvě poloviny, každou polovinu sečte a poté sčítá dva součty. To má tu výhodu, že vyžaduje stejný počet aritmetických operací jako naivní součet (na rozdíl od Kahanova algoritmu, který vyžaduje čtyřnásobek aritmetiky a má latenci čtyřnásobek jednoduchého součtu) a lze jej vypočítat souběžně. Základní případ rekurze by v zásadě mohl být součtem pouze jednoho (nebo nula) čísel, ale k amortizaci režie rekurze by se normálně použil větší základní případ. Ekvivalent párového součtu se používá v mnoha rychlých algoritmech Fourierovy transformace (FFT) a je zodpovědný za logaritmický růst chyb zaokrouhlení v těchto FFT. V praxi s chybami zaokrouhlení náhodných znaků skutečně rostou střední kvadratické chyby párového součtu jako . ${\ displaystyle O (1)}$${\ displaystyle O (\ log n)}$${\ Displaystyle O \ left ({\ sqrt {\ log n}} \ right)}$

Další alternativou je použít aritmetiku s libovolnou přesností , která v zásadě nevyžaduje žádné zaokrouhlování s náklady na mnohem větší výpočetní úsilí. Způsob, jak provádět přesně zaokrouhlené částky pomocí libovolné přesnosti, je adaptivně rozšiřovat pomocí více komponent s plovoucí desetinnou čárkou. Tím se minimalizují výpočetní náklady v běžných případech, kdy není potřeba vysoká přesnost. Další metodu, která používá pouze celočíselnou aritmetiku, ale velký akumulátor, popsali Kirchner a Kulisch ; hardwarovou implementaci popsali Müller, Rüb a Rülling.

Možné zneplatnění optimalizací kompilátoru

V zásadě by dostatečně agresivní optimalizující kompilátor mohl zničit účinnost Kahanova součtu: například pokud kompilátor zjednoduší výrazy podle pravidel asociativity skutečné aritmetiky, může to „zjednodušit“ druhý krok v sekvenci

t = sum + y;

c = (t - sum) - y;

na

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

a potom do

c = 0;

čímž se eliminuje kompenzace chyb. V praxi mnoho kompilátorů nepoužívá při zjednodušování pravidla asociativity (která jsou v aritmetice s plovoucí desetinnou čárkou pouze přibližná), pokud to výslovně nenavrhují možnosti kompilátoru umožňující „nebezpečné“ optimalizace, ačkoli kompilátor Intel C ++ je jedním příkladem, který asociativitu umožňuje -ve výchozím nastavení transformace. Původní verze programovacího jazyka C pro K&R C umožnila kompilátoru změnit pořadí výrazů s plovoucí desetinnou čárkou podle pravidel reálné aritmetické asociativity, ale následná norma ANSI C zakázala přeuspořádání, aby se C lépe hodil pro numerické aplikace ( a více podobné Fortranu , který také zakazuje přeuspořádání), ačkoli v praxi možnosti kompilátoru mohou znovu povolit opětovné objednání, jak je uvedeno výše.

Přenosný způsob, jak lokálně inhibovat takovéto optimalizace, je rozbít jednu z řádků původní formulace na dvě prohlášení a dva z meziproduktů nastavit jako těkavé :

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

Podpora ze strany knihoven

Obecně vestavěné funkce „součtu“ v počítačových jazycích obvykle neposkytují žádné záruky, že bude použit konkrétní algoritmus součtu, natož pak součet Kahan. Standard BLAS pro podprogramy lineární algebry se výslovně vyhýbá nařizování jakéhokoli konkrétního výpočetního pořadí operací z důvodů výkonu a implementace BLAS obvykle nepoužívají součet Kahan.

Standardní knihovna počítačového jazyka Python specifikuje funkci fsum pro přesně zaokrouhlené součty pomocí algoritmu Shewchuk ke sledování více dílčích součtů.

V jazyce Julia výchozí implementace sumfunkce provádí párové součty pro vysokou přesnost s dobrým výkonem, ale externí knihovna poskytuje implementaci Neumaierovy varianty pojmenované sum_kbnpro případy, kdy je zapotřebí vyšší přesnosti.

V jazyce C# balíček nugetů HPCsharp implementuje variantu Neumaier a párové součty : oba jako skalární, datově paralelní pomocí instrukcí procesoru SIMD a paralelní vícejádrové.

Viz také

• Algoritmy pro výpočet rozptylu , které zahrnují stabilní součet

Reference

externí odkazy

• Souhrn s pohyblivou řádovou čárkou, Dr. Dobb's Journal, září 1996