Kahanův součtový algoritmus - Kahan summation algorithm - Wikipedia

v numerická analýza, Kahanův součtový algoritmus, také známý jako kompenzovaný součet,[1] výrazně snižuje numerická chyba v součtu získaném přidáním a sekvence konečný-přesnost čísla s plovoucí desetinnou čárkou ve srovnání se zřejmým přístupem. Toho se dosáhne oddělením průběžná kompenzace (proměnná k hromadění malých chyb).

Zejména jednoduše shrnutí n čísla v pořadí mají nejhorší chybu, která roste úměrně na střední kvadratická chyba, která roste jako pro náhodné vstupy (chyby zaokrouhlení tvoří a náhodná procházka ).[2] S kompenzovaným součtem je nejhorší případ vázané chyby skutečně nezávislý na n, takže lze sčítat velké množství hodnot s chybou, která závisí pouze na pohyblivé řádové čárce přesnost.[2]

The algoritmus je přičítáno William Kahan.[3] Podobné dřívější techniky jsou například Bresenhamův liniový algoritmus, sledování nahromaděné chyby v celočíselných operacích (i když poprvé dokumentováno přibližně ve stejnou dobu[4]) a delta-sigma modulace[5]

Algoritmus

v pseudo kód, algoritmus bude;

funkce KahanSum (vstup) var součet = 0,0 // Připravte akumulátor. var c = 0,0 // Provozní kompenzace za ztracené bity nižšího řádu. pro i = 1 na délka vstupu dělat     // Pole vstup má prvky indexované vstup [1] na vstup [input.length]. var y = vstup [i] - c // C je poprvé nula. var t = součet + y // Bohužel, součet je velký, y malé číslice tak malého řádu y jsou ztraceny. c = (t - součet) - y // (t - součet) zruší část vysoké objednávky y; odečítání y získá negativní (nízká část y) sum = t // Algebraicky, C by měla být vždy nula. Dejte si pozor na příliš agresivní optimalizační kompilátory! další i // Příště bude ztracená spodní část přidána do y v novém pokusu. vrátit se součet

Pracoval příklad

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

S kompenzovaným součtem však získáme správný zaokrouhlený výsledek 10005,9.

Předpokládat, že C má počáteční hodnotu nula.

  y = 3,14159 - 0,00000 y = vstup [i] - c  t = 10000,0 + 3,14159 = 10003,14159 Je však zachováno pouze šest číslic. = 10003,1 Mnoho číslic bylo ztraceno! c = (10003,1 - 10 000,0) - 3,14159 Toto musí být vyhodnocen jako psaný! = 3,10000 - 3,14159 Asimilovaná část y obnoveno oproti původní plné y. = -0,0415900 Zobrazeny koncové nuly, protože se jedná o šestimístný aritmetický součet. = 10003.1 Tedy několik číslic od vstup (tj) se setkal s těmi z součet.

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

  y = 2,71828 - (-0,0415900) Zahrnuje se schodek z předchozí fáze. = 2,75987 Má podobnou velikost y: většina číslic se setká. t = 10003,1 + 2,75987 Ale jen málo z nich splňuje číslice součet. = 10005.85987 A výsledek je zaokrouhlený = 10005,9 na šest číslic. c = (10005,9 - 10003,1) - 2,75987 Toto extrahuje, co vstoupilo dovnitř. = 2,80000 - 2,75987 V tomto případě příliš mnoho. = 0,040130 Ale bez ohledu na to, přebytek by byl příště odečten. Součet = 10005,9 Přesný výsledek je 10005,85987, toto je správně zaokrouhleno na 6 číslic.

Součet se tedy provádí pomocí dvou akumulátorů: součet drží součet a C akumuluje části, které nejsou asimilovány součet, postrčit část nižšího řádu součet příště. Součet tedy pokračuje "strážnými číslicemi" C, což je lepší než žádný žádný, ale není to 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 obecně není praktické; -li vstup je již ve dvojité přesnosti, málo systémů dodává čtyřnásobná přesnost, a pokud ano, vstup by pak mohl být ve čtyřnásobné přesnosti.

Přesnost

Je třeba pečlivě analyzovat chyby v kompenzovaném součtu, abychom ocenili jeho přesnost. I když je přesnější než naivní součet, pro špatně podmíněné součty může stále dávat velké relativní chyby.

Předpokládejme, že jeden sčítá n hodnoty Xi, pro i = 1, ... ,n. Přesná částka je

(počítáno s nekonečnou přesností).

S kompenzovaným součtem jeden místo toho získá , kde došlo k chybě je ohraničen[2]

kde ε je přesnost stroje použité aritmetiky (např. ε ≈ 10−16 pro standard IEEE dvojnásobná přesnost plovoucí bod). Obvykle je množství zájmu relativní chyba , který je tedy výše omezen

Ve výrazu pro relativní vázanou chybu zlomek Σ |Xi| / | ΣXi| je číslo podmínky součtové úlohy. Číslo podmínky v zásadě představuje vnitřní citlivost součtového problému na chyby, bez ohledu na to, jak je vypočítán.[6] Relativní chyba vázaná na každý (dozadu stabilní ) metoda součtu pevným algoritmem s pevnou přesností (tj. ne těmi, které používají.) libovolná přesnost aritmetika, ani algoritmy, jejichž paměťové a časové požadavky se mění na základě dat), jsou úměrné tomuto počtu podmínek.[2] An špatně podmíněný součtový problém je takový, ve kterém je tento poměr velký, a v tomto případě může mít i kompenzovaný součet velkou relativní chybu. Například pokud sčítá Xi jsou nekorelovaná náhodná čísla s nulovým průměrem, součet je a náhodná procházka a počet podmínek bude ú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.

Vzhledem k číslu podmínky je relativní chyba kompenzovaného součtu fakticky nezávislá na n. V zásadě existuje O (nε2), který roste lineárně s n, ale v praxi je tento termín v podstatě nulový: protože konečný výsledek je zaokrouhlen na přesnost ε, nε2 termín zaokrouhluje na nulu, pokud n je zhruba 1 /ε nebo větší.[2] Ve dvojité přesnosti to odpovídá n zhruba 1016, mnohem větší než většina částek. Takže u čísla s pevnou podmínkou jsou chyby kompenzovaného součtu efektivní Ó(ε), nezávislý nan.

Ve srovnání relativní chyba vázaná na naivní součet (jednoduché přidávání čísel v pořadí, zaokrouhlování v každém kroku) roste, jak vynásobeno číslem podmínky.[2] Tato chyba v nejhorším případě je v praxi zřídka pozorována, protože k ní dochází, pouze pokud jsou chyby zaokrouhlování ve stejném směru. V praxi je mnohem pravděpodobnější, že chyby zaokrouhlování mají náhodné znaménko s nulovým průměrem, takže tvoří náhodnou procházku; v tomto případě má naivní součet a střední kvadratická relativní chyba, která roste jako vynásobeno číslem podmínky.[7] To je však stále mnohem horší než kompenzovaný součet. Pokud však lze součet provést s dvojnásobnou přesností, pak ε je nahrazen ε2a naivní součet má nejhorší chybu srovnatelnou s O (nε2) výraz v kompenzovaném součtu s původní přesností.

Ze stejného důvodu je the |Xi| který se objeví v výše je vázán na nejhorší případ, ke kterému dochází pouze v případě, že všechny chyby zaokrouhlování mají stejné znaménko (a jsou maximální možné velikosti).[2] V praxi je pravděpodobnější, že chyby mají náhodné znaménko, v takovém případě výrazy v Σ |Xi| jsou nahrazeny náhodnou chůzí, v takovém případě, dokonce i pro náhodné vstupy s nulovým průměrem, chyba roste jen jako (ignoruje nε2 termín), stejná sazba součet roste, ruší faktory při výpočtu relativní chyby. Takže i pro asymptoticky špatně podmíněné součty může být relativní chyba kompenzovaného součtu často mnohem menší, než by mohla naznačovat analýza v nejhorším případě.

Další vylepšení

Neumaier[8] představil vylepšenou verzi Kahanova algoritmu, kterému říká „vylepšený Kahan – Babuška algoritmus“, který pokrývá také případ, kdy je další termín, který má být přidán, v absolutní hodnotě větší než průběžný součet, čímž se efektivně vymění role toho, co je velké a co je malé. v pseudo kód, algoritmus je:

funkce NeumaierSum (vstup) var součet = 0,0 var c = 0,0 // Provozní kompenzace za ztracené bity nižšího řádu. pro i = 1 na délka vstupu dělat        var t = součet + vstup [i] -li | součet | > = | vstup [i] | pak            c + = (součet - t) + vstup [i] // Pokud součet je větší číslice nízkého řádu vstup [i] jsou ztraceny. jiný            c + = (vstup [i] - t) + součet // Jiné číslice nízkého řádu součet jsou ztraceny. endif        součet = t další i vrátit se sum + c // Oprava byla použita pouze jednou na samém konci.

Pro mnoho sekvencí čísel se oba algoritmy shodují, ale jednoduchý příklad kvůli Petersovi[9] ukazuje, jak se mohou lišit. Pro sčítání s dvojitou přesností získá Kahanův algoritmus hodnotu 0,0, zatímco Neumaierův algoritmus získá správnou hodnotu 2,0.

Jsou také možné úpravy vyššího řádu s lepší přesností. Například varianta navržená Kleinem,[10] který nazval „druhého řádu“ „iterativní algoritmus Kahan – Babuška“. v pseudo kód, algoritmus je:

funkce KleinSum (vstup) var součet = 0,0 var cs = 0,0 var ccs = 0,0 pro i = 1 na délka vstupu dělat        var t = součet + vstup [i] -li | součet | > = | vstup [i] | pak            c = (součet - t) + vstup [i] jiný            c = (vstup [i] - t) + součet endif        součet = t t = cs + c -li | cs | > = | c | pak            cc = (cs - t) + c jiný            cc = (c - t) + cs endif        cs = t ccs = ccs + cc další i vrátit se součet + cs + ccs

Alternativy

Ačkoli Kahanův algoritmus dosahuje nárůst chyb při sčítání n čísla, jen o něco horší růstu lze dosáhnout párový součet: jeden rekurzivně rozděluje množinu čísel na dvě poloviny, sčítá každou polovinu a poté sčítá dvě součty.[2] 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ásobku jednoduchého součtu) a lze jej vypočítat paralelně. Základním případem rekurze by v zásadě mohl být součet pouze jednoho (nebo nulového) čísla, ale k amortizovat režie rekurze, jeden by normálně používal větší základní případ. Ekvivalent párového součtu se používá v mnoha rychlá Fourierova transformace (FFT) a odpovídá za logaritmický růst chyb zaokrouhlování v těchto FFT.[11] V praxi, s chybami zaokrouhlování náhodných znaků, chyby střední kvadratické hodnoty párového součtu ve skutečnosti rostou jako .[7]

Další alternativou je použití aritmetika s libovolnou přesností, které v zásadě nepotřebují žádné zaokrouhlování s náklady na mnohem větší výpočetní úsilí. Způsob provádění přesně zaokrouhlených součtů pomocí libovolné přesnosti je adaptivní rozšíření 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í nutná vysoká přesnost.[12][9] Další metodu, která používá pouze celočíselnou aritmetiku, ale velký akumulátor, popsali Kirchner a Kulisch;[13] hardwarovou implementaci popsali Müller, Rüb a Rülling.[14]

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

V zásadě dostatečně agresivní optimalizace kompilátoru by mohlo zničit účinnost Kahanova součtu: například pokud kompilátor zjednodušil výrazy podle asociativita Pravidla reálné aritmetiky by mohla „zjednodušit“ druhý krok v pořadí

t = součet + y;
c = (t - součet) - y;

na

c = ((součet + y) - součet) - y;

a pak do

c = 0;

čímž se eliminuje kompenzace chyb.[15] V praxi mnoho překladačů nepoužívá pravidla asociativity (která jsou pouze přibližná v aritmetice s plovoucí desetinnou čárkou) ve zjednodušení, pokud k tomu není výslovně nařízeno možnosti kompilátoru umožňující „nebezpečné“ optimalizace,[16][17][18][19] Ačkoliv Překladač Intel C ++ je jeden příklad, který ve výchozím nastavení umožňuje transformace založené na asociativitě.[20] Originál K&R C. verze Programovací jazyk C. dovolil kompilátoru přeuspořádat výrazy s plovoucí desetinnou čárkou podle pravidel reálné aritmetické asociativity, ale následující ANSI C. standardní zakázané opětovné objednávání, aby se C lépe hodilo pro numerické aplikace (a podobně jako Fortran, který rovněž zakazuje opětovné objednávání),[21] i když v praxi mohou možnosti kompilátoru znovu povolit opětovné řazení, jak je uvedeno výše.

Podpora ze strany knihoven

Obecně platí, že 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, mnohem méně Kahanův součet.[Citace je zapotřebí ] The BLAS standard pro lineární algebra podprogramy se výslovně vyhýbají tomu, aby bylo nutné z důvodu výkonu nařizovat jakékoli konkrétní výpočetní pořadí operací,[22] a implementace BLAS obvykle nepoužívají Kahanovu sumaci.

Standardní knihovna Krajta počítačový jazyk určuje fsum funkce pro přesně zaokrouhlený součet pomocí Shewchuk algoritmus[9] sledovat více dílčích součtů.

V Julie jazyk, výchozí implementace součet funkce ano párový součet pro vysokou přesnost s dobrým výkonem,[23] ale externí knihovna poskytuje implementaci pojmenované Neumaierovy varianty sum_kbn pro případy, kdy je potřeba vyšší přesnost.[24]

V C# Jazyk, Balíček nugetů HPCsharp implementuje Neumaierovu variantu a párový součet: oba jako skalární, datově paralelní použití SIMD instrukce procesoru a paralelní vícejádrový procesor.[25]

Viz také

Reference

  1. ^ Striktně existují i ​​další varianty kompenzovaného součtu: viz Higham, Nicholas (2002). Přesnost a stabilita numerických algoritmů (2. vydání). SIAM. str. 110–123. ISBN  978-0-89871-521-7.
  2. ^ A b C d E F G h Higham, Nicholas J. (1993), "Přesnost sčítání s plovoucí desetinnou čárkou" (PDF), SIAM Journal on Scientific Computing, 14 (4): 783–799, CiteSeerX  10.1.1.43.3535, doi:10.1137/0914050, S2CID  14071038.
  3. ^ Kahan, William (leden 1965), „Další poznámky k omezení chyb zkrácení“ (PDF), Komunikace ACM, 8 (1): 40, doi:10.1145/363707.363723, S2CID  22584810, archivovány z originál (PDF) dne 9. února 2018.
  4. ^ Bresenham, Jack E. (leden 1965). "Algoritmus pro počítačové ovládání digitálního plotru" (PDF). IBM Systems Journal. 4 (1): 25–30. doi:10,1147 / sj.41.0025.
  5. ^ Inose, H .; Yasuda, Y .; Murakami, J. (září 1962). "Telemetrický systém podle manipulace s kódem - ΔΣ modulace". Transakce IRE v kosmické elektronice a telemetrii: 204–209. doi:10.1109 / IRET-SET.1962.5008839. S2CID  51647729.
  6. ^ Trefethen, Lloyd N.; Bau, David (1997). Numerická lineární algebra. Philadelphia: SIAM. ISBN  978-0-89871-361-9.
  7. ^ A b Manfred Tasche a Hansmartin Zeuner, Příručka analyticko-výpočetních metod v aplikované matematice, Boca Raton, FL: CRC Press, 2000.
  8. ^ Neumaier, A. (1974). „Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen“ [Analýza zaokrouhlovací chyby některých metod pro sčítání konečných součtů] (PDF). Zeitschrift für Angewandte Mathematik und Mechanik (v němčině). 54 (1): 39–51. Bibcode:1974ZaMM ... 54 ... 39N. doi:10,1002 / zamm.19740540106.
  9. ^ A b C Raymond Hettinger, Recept 393090: Binární sčítání s plovoucí desetinnou čárkou přesné s plnou přesností, Implementace algoritmu Python z článku Shewchuk (1997) (28. března 2005).
  10. ^ A., Klein (2006). „Zobecněný algoritmus Kahan – Babuška-Sumaration-Algorithm“. Výpočetní. Springer-Verlag. 76 (3–4): 279–293. doi:10.1007 / s00607-005-0139-x. S2CID  4561254.
  11. ^ S. G. Johnson a M. Frigo, “Implementace FFT v praxi, v Rychlé Fourierovy transformace, editoval C. Sidney Burrus (2008).
  12. ^ Jonathan R. Shewchuk, Adaptivní přesná aritmetika s pohyblivou řádovou čárkou a rychlé a robustní geometrické predikáty, Diskrétní a výpočetní geometrie, sv. 18, str. 305–363 (říjen 1997).
  13. ^ R. Kirchner, Ulrich W. Kulisch, Přesná aritmetika pro vektorové procesory, Journal of Parallel and Distributed Computing 5 (1988) 250–270.
  14. ^ M. Muller, C. Rub, W. Rulling [1], Přesná akumulace čísel s plovoucí desetinnou čárkouSborník 10. IEEE Symposium on Computer Arithmetic (Červen 1991), doi:10.1109 / ARITH.1991.145535.
  15. ^ Goldberg, David (březen 1991), „Co by měl každý počítačový vědec vědět o aritmetice s plovoucí desetinnou čárkou“ (PDF), ACM Computing Surveys, 23 (1): 5–48, doi:10.1145/103162.103163.
  16. ^ Sbírka překladačů GNU manuál, verze 4.4.3: 3.10 Možnosti, které řídí optimalizaci, -fassociativní matematika (21. ledna 2010).
  17. ^ Uživatelská příručka Compaq Fortran pro systémy Tru64 UNIX a Linux Alpha Archivováno 07.06.2011 na Wayback Machine, oddíl 5.9.7 Optimalizace aritmetického přeskupování (získaná v březnu 2010).
  18. ^ Börje Lindh, Optimalizace výkonu aplikace, Sun BluePrints OnLine (Březen 2002).
  19. ^ Eric Fleegal, “Microsoft Visual C ++ s plovoucí desetinnou čárkou ", Technické články k Microsoft Visual Studio (Červen 2004).
  20. ^ Martyn J. Corden, “Konzistence výsledků s plovoucí desetinnou čárkou pomocí kompilátoru Intel ", Technická zpráva Intel (18. září 2009).
  21. ^ MacDonald, Tom (1991). "C pro numerické výpočty". Journal of Supercomputing. 5 (1): 31–48. doi:10.1007 / BF00155856. S2CID  27876900.
  22. ^ Technické fórum BLAS, oddíl 2.7 (21. srpna 2001), Archivovány na Wayback Machine.
  23. ^ RFC: použijte párový součet pro součet, cumsum a cumprod, github.com/JuliaLang/julia pull request # 4039 (srpen 2013).
  24. ^ KahanSummation knihovna v Julii.
  25. ^ Balíček vysoce výkonných algoritmů HPCsharp nuget.

externí odkazy