Algoritmo di sommatoria di Kahan

In analisi numerica, l'algoritmo di sommatoria di Kahan (conosciuto anche come sommatoria compensata) riduce significativamente l'errore numerico nel totale ottenuto sommando una sequenza di numeri in virgola mobile con precisione finita, se confrontato con il consueto procedimento. Ciò è ottenuto mantenendo una compensazione progressiva separata (una variabile per accumulare piccoli errori).

Quando rappresentiamo un generico numero reale in virgola mobile con precisione finita, ossia con un numero finito di cifre significative, tale rappresentazione, rispetto al numero reale considerato, differisce di un certo valore, il quale corrisponde all'errore di arrotondamento ed è lo scarto tra la rappresentazione in virgola mobile e il numero stesso. Eseguendo una semplice sommatoria di più numeri reali, utilizzando però le rispettive rappresentazioni in virgola mobile, il totale ottenuto presenta un certo errore dato dalla somma algebrica dei singoli errori di arrotondamento, ossia dei singoli scarti, ed, inoltre, si ottiene un certo scarto quadratico medio, in inglese root mean square error, inteso come la radice quadrata della media aritmetica dei quadrati dei singoli scarti.

In particolare, la sommatoria semplice di $n$ numeri in sequenza presenta un errore che, nel caso peggiore, cresce proporzionalmente ad $n$ ed uno scarto quadratico medio che cresce come ${\sqrt {n}}$ per addendi casuali (gli errori di arrotondamento producono una passeggiata aleatoria). Invece, con la sommatoria compensata, l'errore peggiore possibile è indipendente da $n$, dunque un gran numero di valori possono essere sommati con un errore che dipende solo dalla precisione della rappresentazione in virgola mobile.

L'algoritmo è attribuito a William Kahan. Tecniche precedenti simili sono, per esempio, l'algoritmo della linea di Bresenham, che mantiene traccia dell'errore accumulato nelle operazioni sugli interi (anche se già presente in un articolo pubblicato circa nello stesso periodo) e la modulazione Sigma-Delta (che integra e non soltanto somma l'errore).

L'algoritmo

In pseudocodice, l'algoritmo è:

function KahanSum(input)
    var sum = 0.0
    var c = 0.0                 // Una compensazione in esecuzione per i bit meno significativi persi.
    for i = 1 to input.length do
        var y = input[i] - c    // Fin qui, tutto bene: c è zero.
        var t = sum + y         // Purtroppo, sum è grande, y è piccola, dunque le cifre meno significative di y vengono perse.
        c = (t - sum) - y       // (t - sum) annulla la parte più significativa di y; sottraendo y si recupera la parte meno significativa di y cambiata di segno
        sum = t                 // Algebricamente, c dovrebbe sempre essere zero. Occorre prestare attenzione ai compilatori con ottimizzazione troppo aggressiva!
    next i                      // La volta successiva, la parte meno significativa persa verrà aggiunta a "y" in un nuovo tentativo. .
    return sum

Esempio pratico

Questo esempio sarà dato in decimale. I computer, tipicamente, utilizzano l'aritmetica binaria, ma il principio che ora viene illustrato è lo stesso. Supponiamo di stare utilizzando l'aritmetica in virgola mobile decimale a sei cifre e che sum, la somma, abbia raggiunto il valore 10000,0 e che i prossimi due valori degli input(i) siano 3,14159 e 2,71828. Il risultato esatto è 10005,85987, che si arrotonda a 10005,9. Con una sommatoria semplice ogni valore introdotto sarebbe allineato a sum, la somma, e molte cifre meno significative andrebbero perse (per troncamento o arrotondamento). Il primo risultato, dopo l'arrotondamento, sarebbe 10003,1. Il secondo risultato sarebbe 10005,81828 prima dell'arrotondamento e 10005,8 dopo. Ciò non è corretto.

Tuttavia, con la sommatoria compensata, otteniamo il risultato arrotondato corretto di 10005,9.

Assumiamo che c abbia il valore iniziale zero.

  y = 3.14159 - 0                   y = input[i] - c
  t = 10000.0 + 3.14159
    = 10003.14159                   Ma solo sei cifre vengono mantenute. 
    = 10003.1                       Molte cifre sono andate perse!
  c = (10003.1 - 10000.0) - 3.14159 Questo deve essere valutato come scritto! 
    = 3.10000 - 3.14159             La parte assimilata di y, recuperata, in contrapposizione al valore originale completo di y.
    = -.0415900                     Gli zeri finali vengono visualizzati perché si tratta di aritmetica a sei cifre.
sum = 10003.1                       Dunque, poche cifre da input(i) sono rimaste in quelle di sum.

La somma è così grande che solo le cifre più significative dei numeri immessi vengono accumulate. Ma, al passo successivo, c dà un errore.

  y = 2.71828 - -.0415900           Il deficit ottenuto al passo precedente viene incluso.
    = 2.75987                       È di dimensioni simili a  y : le cifre più significative sono rimaste.
  t = 10003.1 + 2.75987             Ma poche sono rimaste nelle cifre di sum
    = 10005.85987                   e il risultato è arrotondato
    = 10005.9                       a sei cifre.
  c = (10005.9 - 10003.1) - 2.75987 Ciò estrae qualunque cosa sia stata ottenuta.
    = 2.80000 - 2.75987             In tal caso, c'è un eccesso.
    = .040130                       Ma non importa, l'eccesso verrebbe sottratto alla volta successiva.
sum = 10005.9                       Il risultato esatto è 10005.85987, esso è correttamente arrotondato a 6 cifre.

In tal modo, la sommatoria è accumulata utilizzando due variabili: sum che si riferisce alla somma e c in cui si accumulano le parti non assimilate in sum, per spostare nuovamente la parte meno significativa di sum la volta successiva. Quindi la sommatoria procede con delle "cifre di guardia" in c che è meglio di niente ma che non è una cosa buona quanto sarebbe effettuare i calcoli con precisione doppia rispetto all'input. Tuttavia, aumentare semplicemente la precisione dei calcoli, in generale, non è pratico; se l'input è già in precisione doppia, pochi sistemi forniscono una precisione quadrupla e se lo fanno, l'input potrebbe quindi essere in precisione quadrupla.

Accuratezza

È necessaria un'attenta analisi degli errori nella sommatoria compensata per apprezzare le sue caratteristiche di accuratezza. Pur essendo più accurata di una sommatoria semplice, essa può ancora dare grandi errori relativi nel caso di somme in particolari cattive condizioni.

Supponiamo di stare sommando $n$ valori $x_{i}$, per $i=1,\,\ldots ,\,n$. La somma esatta è:

$S_{n}=\sum _{i=1}^{n}x_{i}$ (calcolata con precisione infinita)

Con la sommatoria compensata, invece, si ottiene $S_{n}+E_{n}$, dove l'errore $E_{n}$ è limitato da:

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

dove ε è la precisione di macchina dell'aritmetica utilizzata (per esempio, $\varepsilon \approx 10^{-16}$ per il formato a precisione doppia in virgola mobile secondo lo standard IEEE). Generalmente, la quantità di interesse è l'errore relativo ${\frac {|E_{n}|}{|S_{n}|}}$, che è quindi limitato superiormente da:

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

Nell'espressione per il limite dell'errore relativo, la frazione ${\frac {\sum |x_{i}|}{|\sum x_{i}|}}$ è il condizionamento del problema della sommatoria. Essenzialmente, il condizionamento rappresenta la sensibilità intrinseca agli errori del problema della sommatoria, indipendentemente da come è calcolato. Il limite dell'errore relativo di ogni metodo di sommatoria (stabile all'indietro) basato su un algoritmo fissato con precisione fissata (cioè non quelli che usano aritmetica di precisione arbitraria, né algoritmi i cui requisiti di memoria e tempo cambiano in base ai dati), è proporzionale a tale condizionamento. Un problema di sommatoria mal condizionato è uno in cui questo rapporto è grande e in tal caso anche la sommatoria compensata può avere un grande errore relativo. Per esempio, se gli addendi $x_{i}$ sono numeri casuali non correlati con media zero, la somma è una passeggiata aleatoria e il condizionamento crescerà proporzionalmente a ${\sqrt {n}}$. D'altra parte, per input casuali con media non nulla il condizionamento tende per $n\to \infty$ a una costante finita. Se gli input sono tutti non negativi, allora il condizionamento è $1$.

Dato un certo valore per il condizionamento, l'errore relativo della sommatoria compensata è effettivamente indipendente da $n$. In linea di principio, c'è $\mathrm {O} (n\varepsilon ^{2})$ che cresce linearmente con $n$, ma in pratica questo termine è effettivamente zero: poiché il risultato finale è arrotondato ad una precisione $\varepsilon$, il termine $n\varepsilon ^{2}$ tende a zero, tranne il caso in cui $n$ sia approssimativamente $1/\varepsilon$ o superiore. In precisione doppia, ciò corrisponde a un $n$ di approssimativamente $10^{16}$, molto più grande della maggior parte delle somme. Quindi, per un valore fissato del condizionamento, gli errori della sommatoria compensata sono effettivamente $\mathrm {O} (\varepsilon )$, indipendentemente da $n$.

In confronto, il limite dell'errore relativo per la sommatoria semplice (addizionando semplicemente i numeri in sequenza, arrotondamento ad ogni passaggio) cresce come $\mathrm {O} (\varepsilon n)$ moltiplicato per il condizionamento. Questo errore peggiore possibile, tuttavia, raramente è osservato nella pratica, poiché si verifica solo se gli errori di arrotondamento sono tutti nella stessa direzione. In pratica, è molto più probabile che gli errori di arrotondamento abbiano un segno casuale, con media zero, in modo tale da produrre una passeggiata aleatoria; in questo caso, la sommatoria semplice ha uno scarto quadratico medio, riferito agli errori relativi, che cresce come $\mathrm {O} (\varepsilon {\sqrt {n}})$ moltiplicato per il condizionamento. Comunque, questo è ancora molto peggio della sommatoria compensata. Notare, tuttavia, che se la somma può essere eseguita con una precisione doppia, allora $\varepsilon$ viene sostituito da $\varepsilon ^{2}$ e la sommatoria semplice ha un errore peggiore possibile comparabile con il termine $\mathrm {O} (n\varepsilon ^{2})$ nella sommatoria compensata alla precisione originale.

Per lo stesso motivo, il termine $\sum |x_{i}|$ che appare sopra, in $E_{n}$, è il limite peggiore possibile che si verifica solo se tutti gli errori di arrotondamento hanno lo stesso segno (e sono della grandezza massima possibile). In pratica, è più probabile che gli errori abbiano un segno casuale, nel qual caso i termini in $\sum |x_{i}|$ sono sostituiti da una passeggiata aleatoria; in questo caso, anche per input casuali a media zero, l'errore $E_{n}$ cresce solo come $\mathrm {O} (\varepsilon {\sqrt {n}})$ (ignorando il termine $n\varepsilon ^{2}$), con la stessa velocità cresce la somma $S_{n}$, cancellando i fattori ${\sqrt {n}}$ quando viene calcolato l'errore relativo. Pertanto, anche per somme asintoticamente mal condizionate, l'errore relativo per la sommatoria compensata spesso può essere molto più piccolo di quanto un'analisi del caso peggiore possa suggerire.

Ulteriori miglioramenti

Neumaier ha introdotto una leggera modifica dell'algoritmo di Kahan che copre anche il caso in cui il termine successivo da aggiungere è maggiore in valore assoluto rispetto alla somma corrente, scambiando efficacemente il ruolo di ciò che è grande e ciò che è piccolo. In pseudocodice, l'algoritmo è:

 function NeumaierSum(input)
    var sum = input[1]
    var c = 0.0                 // Una compensazione in corso per i bit meno significativi persi.
    for i = 2 to input.length do
        var t = sum + input[i]
        if |sum| >= |input[i]| do
            c += (sum - t) + input[i] // Se sum è più grande, le cifre meno significative di input[i] vanno perse.
        else
            c += (input[i] - t) + sum // Altrimenti, le cifre meno significative di sum vanno perse.
        sum = t
    return sum + c              // Correzione applicata solo una volta alla fine.

Per molte sequenze di numeri, entrambi gli algoritmi concordano ma un simplice esempio dovuto a Peters mostra che essi possono differire. Sommando $[1.0,+10^{100},1.0,-10^{100}]$ in doppia precisione, l'algoritmo di Kahan produce $0.0$ mentre l'algoritmo di Neumaier fornisce il valore corretto $2.0$.

Alternative

Sebbene con l'algoritmo di Kahan, sommando $n$ numeri, si realizzi una crescita dell'errore dell'ordine di $\mathrm {O} (1)$, solo una crescita leggermente peggiore dell'ordine $\mathrm {O} (\log n)$ può essere realizzata da una sommatoria ricorsiva a coppie: si divide in modo ricorsivo l'insieme di numeri in due metà, sommando ogni metà, e quindi si addizionano le due somme. Questo presenta il vantaggio di richiedere lo stesso numero di operazioni aritmetiche della sommatoria semplice (a differenza dell'algoritmo di Kahan, che, rispetto a una sommatoria semplice, richiede un numero di operazioni aritmetiche e ha una latenza quattro volte maggiore) e il calcolo può essere eseguito in parallelo. Il modello di base della ricorsione potrebbe essere, in linea di principio, la somma solo di un numero (o zero numeri), ma per contenere il sovraccarico della ricorsione, normalmente, si usa un modello di base più ampio. L'equivalente della sommatoria ricorsiva a coppie è usato in molti algoritmi di trasformata di Fourier veloce (FFT), ed è responsabile della crescita logaritmica degli errori di arrotondamento in tali FFT. In pratica, con errori di arrotondamento di segno casuale, gli scarti quadratici medi della sommatoria ricorsiva a coppie, effettivamente, crescono come $\mathrm {O} ({\sqrt {\log n}})$.

Un'altra alternativa è usare l'aritmetica a precisione arbitraria, che in linea di principio non richiede affatto arrotondamenti con il costo di uno sforzo computazionale molto maggiore. Un modo per eseguire somme perfettamente arrotondate usando precisione arbitraria consiste nell'estenderle in modo adattivo usando componenti multipli in virgola mobile. Questo minimizzerà il costo computationale nei casi comuni in cui non è necessaria un'elevata precisione. Un altro metodo che utilizza solo aritmetica intera, ma un grande accumulatore è stato descritto da Kirchner e Kulisch; un'implementazione hardware è stata descritta da Müller, Rüb e Rülling.

Possibile invalidazione dall'ottimizzazione del compilatore

In linea di principio, un'ottimizzazione del compilatore sufficientemente aggressiva potrebbe compromettere l'efficacia della sommatoria di Kahan: per esempio, se il compilatore ha semplificato delle espressioni in base alle regole di associatività dell'aritmetica reale, potrebbe semplificare il secondo passo della sequenza da t = sum + y; c = (t - sum) - y; a ((sum + y) - sum) - y; poi a c = 0;, eliminando la compensazione dell'errore. In pratica, molti compilatori non usano le regole di associatività (che, nell'aritmetica in virgola mobile, sono solo approssimate) nelle semplificazioni, salvo che esplicitamente comandati a farlo dalle opzioni del compilatore che consentano ottimizzazioni non sicure, sebbene il compilatore Intel C++ sia un esempio che consente trasformazioni basate sull'associatività per impostazione predefinita. La versione originale K&R C del linguaggio C consentiva al compilatore di riordinare le espressioni in virgola mobile in base alle regole di associatività dell'aritmetica reale, ma il successivo standard ANSI C non consentiva il riordino allo scopo di rendere il linguaggio C più adatto per le applicazioni numeriche (e più simile al Fortran, che non permette il riordino), anche se, in pratica, le opzioni del compilatore possono abilitare nuovamente il riordino come menzionato sopra.

Supporto con le librerie

In generale, le funzioni di somma incluse nei linguaggi di programmazione tipicamente non forniscono alcuna garanzia che venga utilizzato un particolare algoritmo di sommatoria, come la sommatoria di Kahan. Lo standard BLAS (Basic Linear Algebra Subprograms) per le subroutines di algebra lineare evita esplicitamente di imporre un particolare ordine di calcolo delle operazioni per motivi di prestazioni, e le implementazioni BLAS, tipicamente, non utilizzano la sommatoria di Kahan.

La libreria standard del linguaggio Python specifica una funzione fsum. per sommatorie arrotondate esattamente, usando l'algoritmo di Shewchuk per tracciare più somme parziali.

Nel linguaggio Julia, l'implementazione predefinita della funzione sum utilizza la sommatoria ricorsiva a coppie per un'elevata accuratezza con buone prestazioni, ma la libreria standard ha anche un'implementazione della variante di Neumaier's chiamata sum_kbn per i casi in cui è richiesta un'accuratezza più elevata.

Bibliografia

• Luciano M. Barone, Enzo Marinari, Giovanni Organtini e Federico Ricci-Terseng, Programmazione Scientifica. Linguaggio C, algoritmi e modelli nella scienza, Pearson Education Italia, 2006, p. 108, paragrafo 4.5: Un problema di arrotondamento, ISBN 978-88-7192-242-3.

Voci correlate

• Algoritmi per il calcolo della varianza
• Sommatoria ricorsiva a coppie
• Condizionamento (matematica)
• Epsilon di macchina
• Formato a precisione doppia in virgola mobile
• Formato a precisione quadrupla in virgola mobile
• Aritmetica a precisione arbitraria

Collegamenti esterni

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