Kahan summeringsalgoritm

I numerisk analys minskar Kahan -summeringsalgoritmen , även känd som kompenserad summering , avsevärt det numeriska felet i den totala summan som erhålls genom att lägga till en sekvens av flyttalstal med ändlig precision , jämfört med den uppenbara metoden. Detta görs genom att behålla en separat löpande kompensation (en variabel för att ackumulera små fel), vilket i själva verket utökar summans precision med precisionen hos kompensationsvariabeln.

I synnerhet, att bara summera $n$ tal i följd har ett värsta tänkbara fel som växer proportionellt mot $n$ och ett rotmedelkvadratfel som växer som ${\sqrt {n} }$ för slumpmässiga inmatningar (avrundningsfelen bildar en slumpmässig gång ). Med kompenserad summering, med användning av en kompensationsvariabel med tillräckligt hög precision, är den värsta felgränsen effektivt oberoende av ${\displaystyle n},$ så ett stort antal värden kan summeras med ett fel som bara beror på flyttalsprecisionen av resultatet.

Algoritmen tillskrivs William Kahan ; Ivo Babuška verkar ha kommit på en liknande algoritm oberoende (därav Kahan–Babuška summering ) . Liknande tidigare tekniker är till exempel Bresenhams linjealgoritm , som håller reda på det ackumulerade felet i heltalsoperationer (även om det först dokumenterades ungefär samtidigt) och delta-sigma-modulationen .

Algoritmen

I pseudokod kommer algoritmen att vara:

 function  KahanSum(input)  var  summa = 0.0 // Förbered ackumulatorn.  var  c = 0.0 // En löpande kompensation för förlorade bitar av låg ordning.  for  i = 1  till  input.length  do  // Arrayingången  har element indexerad input[1] till input[input.length]  .  var  y = input[i] - c //  c  är noll första gången.  var  t = summa + y // Tyvärr,  summan  är stor,  y  liten, så siffror av låg ordning i  y  går förlorade. c = (t - summa) - y //   (t - summa)  avbryter högordningens del av  y  ; subtrahera   y  återvinner negativ (låg del av  y  ) summa = t // Algebraiskt bör   c  alltid vara noll. Akta dig för alltför aggressiva optimeringskompilatorer!   next  i // Nästa gång kommer den förlorade låga delen att läggas till  y  i ett nytt försök.  retursumma  _

Denna algoritm kan också skrivas om för att använda Fast2Sum- algoritmen:

 function  KahanSum2(input)  var  sum = 0.0 // Förbered ackumulatorn.  var  c = 0.0 // En löpande kompensation för förlorade bitar av låg ordning.  for  i = 1  till  input.length  do  // Arrayingången  har element indexerad input[1] till input[input.length]  .  var  y = input[i] + c //  c  är noll första gången. (summa,c) = Fast2Sum(summa,y) //   summa  +  c  är en approximation till den exakta summan.  next  i // Nästa gång kommer den förlorade låga delen att läggas till  y  i ett nytt försök.  retursumma  _

Arbetat exempel

Detta exempel kommer att ges med decimaler. Datorer använder vanligtvis binär aritmetik, men principen som illustreras är densamma. Anta att vi använder sexsiffrig decimal flyttalsaritmetik , summan har uppnått värdet 10000,0 och de följande två värdena för indata[i] är 3,14159 och 2,71828. Det exakta resultatet är 10005,85987, vilket avrundas till 10005,9. Med en vanlig summering skulle varje inkommande värde justeras med summa och många siffror av låg ordning skulle gå förlorade (genom trunkering eller avrundning). Det första resultatet, efter avrundning, skulle vara 10003,1. Det andra resultatet skulle vara 10005,81828 före avrundning och 10005,8 efter avrundning. Detta är inte korrekt.

Men med kompenserad summering får vi det korrekta avrundade resultatet på 10005,9.

Antag att c har startvärdet noll.

 y = 3,14159 - 0,00000  y = input[i] - c  t = 10000,0 + 3,14159 = 10003,14159 Men bara sex siffror behålls. = 10003.1 Många siffror har gått förlorade!  c = (10003,1 - 10000,0) - 3,14159 Detta   måste  utvärderas som skrivet! = 3,10000 - 3,14159 Den assimilerade delen av   y  återhämtade sig, kontra det ursprungliga fulla  y  . = -0,0415900 Efterföljande nollor visas eftersom detta är sexsiffrig aritmetik.  summa = 10003.1 Således, några siffror från   input(i  ) uppfyllde de för  summa  .

Summan är så stor att endast de höga siffrorna i de inmatade numren ackumuleras. Men i nästa steg c felet.

 y = 2,71828 - (-0,0415900) Underskottet från föregående steg inkluderas. = 2,75987 Den har en storlek som liknar   y  : de flesta siffror möts. t = 10003,1 + 2,75987 Men få uppfyller siffrorna för   summa  . = 10005,85987 Och resultatet är avrundat = 10005,9 till sex siffror.  c = (10005,9 - 10003,1) - 2,75987 Detta extraherar allt som gick in. = 2,80000 - 2,75987 I det här fallet, för mycket.  = 0,040130 Men oavsett, överskottet skulle dras av nästa gång.  summa = 10005,9 Exakt resultat är 10005,85987, detta är korrekt avrundat till 6 siffror.

Så summeringen utförs med två ackumulatorer: summan håller summan och c ackumulerar delarna som inte assimilerats till summan , för att knuffa den låga delen av summan nästa gång. Således fortsätter summeringen med "vaktsiffror" i c , vilket är bättre än att inte ha några, men är inte lika bra som att utföra beräkningarna med dubbel precision av inmatningen. Men att bara öka precisionen i beräkningarna är inte praktiskt i allmänhet; om inmatningen redan är i dubbel precision, är det få system som levererar fyrfaldig precision , och om de gjorde det, skulle inmatningen kunna vara i fyrfaldig precision.

Noggrannhet

En noggrann analys av felen i kompenserad summering behövs för att uppskatta dess noggrannhetsegenskaper. Även om det är mer korrekt än naiv summering, kan det ändå ge stora relativa fel för dåligt betingade summor.

Antag att man summerar $n$ värden $x_{i}$ , för $i=1,\,\ldots ,\,n$ . Den exakta summan är

S_{n}=\summa _{i=1}^{n}x_{i}

(beräknad med oändlig precision).

Med kompenserad summering får man istället $S_{n}+E_{n}$ , där felet $E_{n}$ begränsas av

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

där $\varepsilon$ är maskinprecisionen för den aritmetik som används (t.ex. $\varepsilon \approx 10^{-16}$ för IEEE-standard med dubbelprecision flytpunkt). Vanligtvis är kvantiteten av intresse det relativa felet ${\displaystyle |E_{n}|/|S_{n}|} ,$ som därför ovan begränsas av

{\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|} }.

I uttrycket för den relativa felgränsen, bråkdelen $\Sigma |x_{i}|/|\Sigma x_{i}|$ är villkorsnumret för summeringsproblemet. I huvudsak representerar villkorsnumret summationsproblemets inneboende känslighet för fel, oavsett hur det beräknas. Den relativa felgränsen för varje ( bakåtstabil ) summeringsmetod med en fast algoritm i fast precision (dvs inte de som använder aritmetik med godtycklig precision , inte heller algoritmer vars minnes- och tidskrav ändras baserat på data), är proportionell mot detta villkorsnummer . Ett dåligt betingat summeringsproblem är ett där detta förhållande är stort, och i detta fall kan även kompenserad summering ha ett stort relativt fel. Till exempel, om summorna $x_{i}$ är okorrelerade slumptal med noll medelvärde, är summan en slumpmässig promenad , och villkorstalet kommer att växa proportionellt mot ${\sqrt {n}}$ . Å andra sidan, för slumpmässiga indata med icke-noll betyder villkorstalet asymptoter till en finit konstant som $n\to \infty$ . Om alla ingångar är icke-negativa är villkorsnumret 1.

Givet ett villkorsnummer är det relativa felet för kompenserad summering i praktiken oberoende av $n$ . I princip finns det $O(n\varepsilon ^{2})$ som växer linjärt med $n$ , men i praktiken är denna term faktiskt noll: eftersom slutresultatet avrundas till en precision $\varepsilon$ , termen $n\varepsilon ^{2}$ avrundas till noll, om inte $n$ är ungefär $1/ \varepsilon$ eller större. I dubbel precision motsvarar detta ett $n$ på ungefär ${\displaystyle 10^{16}} ,$ mycket större än de flesta summor. Så, för ett fast villkorsnummer, är felen för kompenserad summering i praktiken ${\displaystyle O(\varepsilon )} ,$ oberoende av $n$ .

I jämförelse, det relativa felet gränsat för naiv summering (att helt enkelt addera talen i följd, avrundning vid varje steg) växer när $O(\varepsilon n)$ multiplicerat med villkorsnumret. Detta värsta tänkbara fel observeras dock sällan i praktiken, eftersom det bara inträffar om avrundningsfelen är alla i samma riktning. I praktiken är det mycket mer troligt att avrundningsfelen har ett slumpmässigt tecken, med noll medelvärde, så att de bildar en slumpmässig gång; i detta fall har naiv summering ett rotmedelkvadrat- relativt fel som växer som $O\left(\varepsilon {\sqrt {n}}\right)$ multiplicerat med villkorsnumret. Detta är dock fortfarande mycket värre än kompenserad summering. Men om summan kan utföras med dubbelt så hög precision, så ersätts ${\displaystyle \varepsilon } med$ $\varepsilon ^{2}$ , och naiv summering har ett värsta tänkbart fel jämförbart med $O(n\varepsilon ^{2})$ term i kompenserad summering vid den ursprungliga precisionen.

På samma sätt är $\Sigma |x_{i}|$ som visas i $E_{n}$ ovan är en gräns i värsta fall som uppstår endast om alla avrundningsfel har samma tecken (och är av maximalt möjliga magnitud). I praktiken är det mer troligt att felen har slumpmässiga tecken, i vilket fall termer i $\Sigma |x_{i}|$ ersätts av en slumpmässig gång, i vilket fall, även för slumpmässiga inmatningar med noll medelvärde, växer felet $E_{n}$ bara som $O\left(\varepsilon {\sqrt {n}}\right)$ (om man ignorerar termen $n\varepsilon ^{2}$ ), samma takt summan ${\ displaystyle S_{n}}$ växer, vilket tar bort de ${\displaystyle {\sqrt {n}}}-$ faktorerna när det relativa felet beräknas. Så även för asymptotiskt dåligt konditionerade summor kan det relativa felet för kompenserad summering ofta vara mycket mindre än en analys i värsta fall kan antyda.

Ytterligare förbättringar

Neumaier introducerade en förbättrad version av Kahan-algoritmen, som han kallar en "förbättrad Kahan–Babuška-algoritm", som också täcker fallet när nästa term som ska läggas till är större i absolut värde än den löpande summan, vilket effektivt byter ut rollen som vad som är stort och vad som är litet. I pseudokod är algoritmen:

         function  KahanBabushkaNeumaierSum(input)  var  sum = 0.0  var  c = 0.0 // En löpande kompensation för förlorade bitar av låg ordning.  för  i = 1  till  input.length  var  t = summa + input[i]  om  |  summa| >= |ingång[i]|   sedan  c += (summa - t) + input[i] // Om  summan är större,   försvinner  siffror av låg ordning i  input[i] .   annars  c += (input[i] - t) + summa // Annars går siffror av låg ordning för  summan  förlorade.  endif  summa = t  nästa  jag  returnerar  summa + c // Korrigering tillämpas endast en gång i slutet.

Denna förbättring liknar ersättningen av Fast2Sum med 2Sum i Kahans algoritm omskriven med Fast2Sum.

För många talsekvenser är båda algoritmerna överens, men ett enkelt exempel på grund av Peters visar hur de kan skilja sig åt. För att summera $[1.0,+10^{100},1.0,-10^{100}]$ i dubbel precision, ger Kahans algoritm 0.0, medan algoritmen ger rätt värde 2,0.

Modifieringar av högre ordning med bättre noggrannhet är också möjliga. Till exempel en variant som föreslagits av Klein, som han kallade en andra ordningens "iterativ Kahan–Babuška-algoritm". I pseudokod är algoritmen:

         

     funktion  KahanBabushkaKleinSum(input)  var  summa = 0,0  var  cs = 0,0  var  ccs = 0,0  var  c = 0,0  var  cc = 0,0  för  i = 1  till  input.length  do  var  t = summa + input[i]  om  |summa| >= |ingång[i]|   då  c = (summa - t) + input[i]  annars  c = (ingång[i] - t) + summa  endif  summa = t t = cs + c  om  |cs| >= |c|   sedan  cc = (cs - t) + c  annars  cc = (c - t) + cs  endif  cs = t ccs = ccs + cc  slut loop  retursumma  + cs + ccs

Alternativ

Även om Kahans algoritm uppnår ${\displaystyle O(1)}-$ feltillväxt för att summera n tal, kan endast något sämre ${\displaystyle O(\log n)}-$ tillväxt uppnås genom parvis summering : man delar rekursivt upp siffrorna i två halvor, summerar varje halva och adderar sedan de två summorna. Detta har fördelen av att kräva samma antal aritmetiska operationer som den naiva summeringen (till skillnad från Kahans algoritm, som kräver fyra gånger aritmetiken och har en latens på fyra gånger en enkel summering) och kan beräknas parallellt. Basfallet för rekursionen skulle i princip kunna vara summan av endast ett (eller noll) tal, men för att amortera omkostnaden för rekursionen skulle man normalt använda ett större basfall. Motsvarigheten till parvis summering används i många snabba Fourier-transformations (FFT)-algoritmer och är ansvarig för den logaritmiska tillväxten av avrundningsfel i dessa FFT. I praktiken, med avrundningsfel av slumpmässiga tecken, växer rotmedelkvadratfelen för parvis summering faktiskt som $O\left({\sqrt {\log n}}\right)$ .

Ett annat alternativ är att använda aritmetik med godtycklig precision , som i princip inte behöver någon avrundning alls med en kostnad av mycket större beräkningsansträngning. Ett sätt att utföra exakt avrundade summor med godtycklig precision är att förlänga adaptivt med hjälp av flera flyttalskomponenter. Detta kommer att minimera beräkningskostnaderna i vanliga fall där hög precision inte behövs. En annan metod som endast använder heltalsaritmetik, men en stor ackumulator, beskrevs av Kirchner och Kulisch ; en hårdvaruimplementering beskrevs av Müller, Rüb och Rülling.

Möjlig ogiltigförklaring genom kompilatoroptimering

I princip kan en tillräckligt aggressiv optimerande kompilator förstöra effektiviteten av Kahan-summering: till exempel, om kompilatorn förenklade uttryck enligt associativitetsreglerna för verklig aritmetik, kan det "förenkla" det andra steget i sekvensen

t = summa + y;

c = (t - summa) - y;

till

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

och sedan till

c = 0;

eliminerar därmed felkompensationen. I praktiken använder många kompilatorer inte associativitetsregler (som endast är ungefärliga i flyttalsaritmetik) i förenklingar, såvida de inte uttryckligen uppmanas att göra det av kompilatoralternativ som möjliggör "osäkra" optimeringar, även om Intel C++-kompilatorn är ett exempel som tillåter associativitet -baserade transformationer som standard. Den ursprungliga K&R C- versionen av programmeringsspråket C gjorde det möjligt för kompilatorn att ordna om flyttalsuttryck enligt realaritmetiska associativitetsregler, men den efterföljande ANSI C- standarden förbjöd omordning för att göra C bättre lämpad för numeriska applikationer ( och mer liknande Fortran , som också förbjuder ombeställning), även om kompilatoralternativ i praktiken kan återaktivera ombeställning, som nämnts ovan.

Ett bärbart sätt att hämma sådana optimeringar lokalt är att dela upp en av raderna i den ursprungliga formuleringen i två påståenden och göra två av mellanprodukterna flyktiga :

         funktion  KahanSum(ingång)  var  summa = 0,0  var  c = 0,0  för  i = 1  till  input.length  do  var  y = input[i] - c  volatile var  t = summa + y  volatile var  z = t - summa c = z - y summa = t  nästa  jag  returnerar  summa

Stöd av bibliotek

I allmänhet ger inbyggda "summa"-funktioner i datorspråk vanligtvis inga garantier för att en viss summeringsalgoritm kommer att användas, än mindre Kahan-summering. ^{[ citat behövs ]} BLAS - standarden för subrutiner för linjär algebra undviker uttryckligen att kräva att någon särskild beräkningsordning av operationer av prestandaskäl, och BLAS-implementeringar använder vanligtvis inte Kahan-summering.

Standardbiblioteket för datorspråket Python specificerar en fsumfunktion för exakt avrundad summering, med hjälp av Shewchuk -algoritmen för att spåra flera delsummor.

I Julia -språket gör standardimplementeringen av summafunktionen parvis summering för hög noggrannhet med bra prestanda, men ett externt bibliotek tillhandahåller en implementering av Neumaiers variant som heter sum_kbn för de fall då högre noggrannhet behövs.

I C# -språket implementerar HPCsharp nuget-paketet Neumaier-varianten och parvis summering : både som skalär, dataparallell med SIMD- processorinstruktioner och parallella flerkärniga.

Se även

Algoritmer för beräkning av varians, som inkluderar stabil summering

externa länkar

Floating-point Summation, Dr. Dobb's Journal september 1996