Algoritmo de suma de Kahan

En análisis numérico , el algoritmo de suma de Kahan , también conocido como suma compensada , ^{[ 1 ]} reduce significativamente el error numérico en el total obtenido al sumar una secuencia de números de coma flotante de precisión finita , en comparación con el método ingenuo. Esto se logra manteniendo una compensación acumulativa independiente (una variable para acumular pequeños errores), extendiendo así la precisión de la suma en la precisión de la variable de compensación.

En particular, simplemente sumando${\displaystyle n}$Los números en secuencia tienen un error en el peor de los casos que crece proporcionalmente a${\displaystyle n}$y un error cuadrático medio que crece a medida que${\displaystyle {\sqrt {n}}}$para entradas aleatorias (los errores de redondeo forman un paseo aleatorio ). ^{[ 2 ]} Con la suma compensada, utilizando una variable de compensación con una precisión suficientemente alta, el límite de error en el peor de los casos es efectivamente independiente de${\displaystyle n}$, por lo que se puede sumar una gran cantidad de valores con un error que solo depende de la precisión de punto flotante del resultado. ^{[ 2 ]}

El algoritmo se atribuye a William Kahan ; ^{[ 3 ]} Ivo Babuška parece haber ideado un algoritmo similar de forma independiente (de ahí la suma de Kahan - Babuška ). ^{[ 4 ]} Técnicas similares anteriores son, por ejemplo, el algoritmo de línea de Bresenham , que registra el error acumulado en operaciones con enteros (aunque se documentó por primera vez aproximadamente al mismo tiempo ^{[ 5 ]} ) y la modulación delta-sigma . ^{[ 6 ]}

En pseudocódigo , el algoritmo será:

función KahanSum(entrada) // Preparar el acumulador. var suma = 0.0 // Una compensación continua por la pérdida de bits de orden bajo. var c = 0.0 // El array input tiene elementos indexados desde input[1] hasta input[input.length]. for i = 1 to input.length do // c es cero la primera vez. var y = input[i] - c // Por desgracia, la suma es grande, y es pequeña, por lo que se pierden los dígitos de menor orden de y . var t = suma + y // (t - suma) cancela la parte de orden superior de y ; // restar y recupera el valor negativo (parte baja de y ) c = (t - suma) - y // Algebraicamente, c siempre debe ser cero. Cuidado // ¡Compiladores de optimización excesivamente agresivos! suma = t // La próxima vez, la parte baja perdida se agregará a y en un nuevo intento. siguiente i suma de retorno

Este algoritmo también puede reescribirse para utilizar el algoritmo Fast2Sum : ^{[ 7 ]}

función KahanSum2(entrada) // Preparar el acumulador. var suma = 0.0 // Una compensación continua por la pérdida de bits de orden bajo. var c = 0.0 // El array input tiene elementos indexados desde input[1] hasta input[input.length]. for i = 1 to input.length do // c es cero la primera vez. var y = input[i] + c // suma + c es una aproximación a la suma exacta. (suma,c) = Fast2Sum(suma,y) // La próxima vez, la parte baja perdida se agregará a y en un nuevo intento. siguiente i suma de retorno

 y = 3.14159 - 0.00000 y = input[i] - c t = 10000.0 + 3.14159 t = suma + y = 10003.14159 Normalización realizada, siguiente redondeo a seis dígitos. = 10003.1 Pocos dígitos de input[i] coincidieron con los de sum . ¡Se han perdido muchos dígitos! c = (10003.1 - 10000.0) - 3.14159 c = (t - suma) - y (Nota: ¡Primero se deben evaluar los paréntesis!) = 3,10000 - 3,14159 La parte asimilada de y menos la y completa original . = -0,0415900 Debido a que c está cerca de cero, la normalización conserva muchos dígitos después del punto flotante. suma = 10003.1 suma = t

 y = 2,71828 - (-0,0415900) La mayoría de los dígitos coinciden, ya que c tiene un tamaño similar al de y . = 2,75987 El déficit (pérdida de dígitos de orden inferior) de la iteración anterior se ha restablecido correctamente. t = 10003.1 + 2.75987 Pero aún así, solo unos pocos cumplen con los dígitos de la suma . = 10005.85987 Normalización realizada, siguiente ronda a seis dígitos. = 10005.9 Nuevamente, se han perdido muchos dígitos, pero c ayudó a corregir el redondeo. c = (10005.9 - 10003.1) - 2.75987 Estime el error acumulado, basado en la y ajustada . = 2,80000 - 2,75987 Como era de esperar, las partes de orden bajo se pueden conservar en c con efectos de redondeo nulos o mínimos. = 0,0401300 En esta iteración, t fue un poco demasiado alto, el exceso se restará en la siguiente iteración. suma = 10005.9 El resultado exacto es 10005.85987, la suma es correcta, redondeada a 6 dígitos.

Es necesario un análisis cuidadoso de los errores en la suma compensada para apreciar sus características de precisión. Si bien es más precisa que la suma ingenua, aún puede generar grandes errores relativos para sumas mal condicionadas.

Supongamos que uno está sumando${\displaystyle n}$valores${\displaystyle x_{i}}$, para${\displaystyle i=1,\,\ldots ,\,n}$La suma exacta es

${\displaystyle S_{n}=\sum _{i=1}^{n}x_{i}}$(calculado con precisión infinita).

Con la suma compensada, se obtiene en cambio${\displaystyle S_{n}+E_{n}}$donde el error${\displaystyle E_{n}}$está delimitado por ^{[ 2 ]}

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

dónde${\displaystyle \varepsilon }$es la precisión de la máquina de la aritmética que se emplea (por ejemplo${\displaystyle \varepsilon \approx 10^{-16}}$(para el estándar IEEE de punto flotante de doble precisión ). Por lo general, la cantidad de interés es el error relativo.${\displaystyle |E_{n}|/|S_{n}|}$, que por lo tanto está delimitado superiormente por

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

En la expresión para el límite de error relativo, la fracción${\displaystyle \Sigma |x_{i}|/|\Sigma x_{i}|}$es el número de condición del problema de suma. Esencialmente, el número de condición representa la sensibilidad intrínseca del problema de suma a los errores, independientemente de cómo se calcule. ^{[ 8 ]} El límite de error relativo de cada método de suma ( estable hacia atrás ) mediante un algoritmo fijo en precisión fija (es decir, no aquellos que utilizan aritmética de precisión arbitraria , ni algoritmos cuyos requisitos de memoria y tiempo cambian en función de los datos), es proporcional a este número de condición. ^{[ 2 ]} Un problema de suma mal condicionado es aquel en el que esta razón es grande, y en este caso incluso la suma compensada puede tener un gran error relativo. Por ejemplo, si los sumandos${\displaystyle x_{i}}$son números aleatorios no correlacionados con media cero, la suma es un paseo aleatorio y el número de condición crecerá proporcionalmente a${\displaystyle {\sqrt {n}}}$Por otro lado, para entradas aleatorias con media distinta de cero, el número de condición tiende asintóticamente a una constante finita como${\displaystyle n\to \infty }$Si todas las entradas son no negativas , entonces el número de condición es 1.

Dado un número de condición, el error relativo de la suma compensada es efectivamente independiente de${\displaystyle n}$. En principio, existe el${\displaystyle O(n\varepsilon ^{2})}$que crece linealmente con${\displaystyle n}$, pero en la práctica este término es efectivamente cero: puesto que el resultado final se redondea a una precisión${\displaystyle \varepsilon }$, el${\displaystyle n\varepsilon ^{2}}$El término se redondea a cero, a menos que${\displaystyle n}$es aproximadamente${\displaystyle 1/\varepsilon }$o mayor. ^{[ 2 ]} En doble precisión, esto corresponde a un${\displaystyle n}$de aproximadamente${\displaystyle 10^{16}}$, mucho mayores que la mayoría de las sumas. Por lo tanto, para un número de condición fijo, los errores de la suma compensada son efectivamente${\displaystyle O(\varepsilon)}$, independientemente de${\displaystyle n}$.

En comparación, el límite de error relativo para la suma ingenua (simplemente sumar los números en secuencia, redondeando en cada paso) crece a medida que${\displaystyle O(\varepsilon n)}$multiplicado por el número de condición. ^{[ 2 ]} Este error en el peor de los casos rara vez se observa en la práctica, sin embargo, porque solo ocurre si todos los errores de redondeo están en la misma dirección. En la práctica, es mucho más probable que los errores de redondeo tengan un signo aleatorio, con media cero, de modo que formen un paseo aleatorio; en este caso, la suma ingenua tiene un error relativo cuadrático medio que crece como${\displaystyle O\left(\varepsilon {\sqrt {n}}\right)}$multiplicado por el número de condición. ^{[ 9 ]} Sin embargo, esto sigue siendo mucho peor que la suma compensada. No obstante, si la suma se puede realizar con el doble de precisión, entonces${\displaystyle \varepsilon }$es reemplazado por${\displaystyle \varepsilon ^{2}}$y la suma ingenua tiene un error en el peor de los casos comparable al${\displaystyle O(n\varepsilon ^{2})}$término en la suma compensada con la precisión original.

Del mismo modo, el${\displaystyle \Sigma |x_{i}|}$que aparece en${\displaystyle E_{n}}$arriba es un límite del peor caso que ocurre solo si todos los errores de redondeo tienen el mismo signo (y son de la máxima magnitud posible). ^{[ 2 ]} En la práctica, es más probable que los errores tengan un signo aleatorio, en cuyo caso los términos en${\displaystyle \Sigma |x_{i}|}$se reemplazan por un paseo aleatorio, en cuyo caso, incluso para entradas aleatorias con media cero, el error${\displaystyle E_{n}}$crece solo como${\displaystyle O\left(\varepsilon {\sqrt {n}}\right)}$(ignorando el${\displaystyle n\varepsilon ^{2}}$término), la misma tasa la suma${\displaystyle S_{n}}$crece, cancelando el${\displaystyle {\sqrt {n}}}$factores al calcular el error relativo. Por lo tanto, incluso para sumas asintóticamente mal condicionadas, el error relativo para la suma compensada a menudo puede ser mucho menor de lo que podría sugerir un análisis del peor caso.

Aunque el algoritmo de Kahan logra${\displaystyle O(1)}$El crecimiento del error al sumar n números es solo ligeramente peor.${\displaystyle O(\log n)}$El crecimiento se puede lograr mediante la suma por pares : se divide recursivamente el conjunto de números en dos mitades, se suma cada mitad y luego se suman las dos sumas. ^{[ 2 ]} Esto tiene la ventaja de requerir el mismo número de operaciones aritméticas que la suma ingenua (a diferencia del algoritmo de Kahan, que requiere cuatro veces la aritmética y tiene una latencia cuatro veces mayor que una suma simple) y se puede calcular en paralelo. El caso base de la recursión podría ser, en principio, la suma de solo uno (o cero) números, pero para amortizar la sobrecarga de la recursión, normalmente se usaría un caso base mayor. El equivalente de la suma por pares se usa en muchos algoritmos de transformada rápida de Fourier (FFT) y es responsable del crecimiento logarítmico de los errores de redondeo en esas FFT. ^{[ 15 ]} En la práctica, con errores de redondeo de signos aleatorios, los errores cuadráticos medios de la suma por pares crecen como${\displaystyle O\left({\sqrt {\log n}}\right)}$. ^{[ 9 ]} La suma por pares es utilizada por NumPy y Julia, entre otros paquetes numéricos notables.

Otra alternativa es utilizar aritmética de precisión arbitraria , que en principio no necesita ningún redondeo, aunque a un coste computacional mucho mayor.

• Una forma de realizar sumas correctamente redondeadas con precisión arbitraria es extender de forma adaptativa utilizando múltiples componentes de punto flotante (método de Shewchuk). Esto minimizará el costo computacional en casos comunes donde no se requiere alta precisión. ^{[ 16 ]} Este estilo de suma es utilizado por la función de Python math.fsum. ^{[ 11 ]}
• Otro método que utiliza únicamente aritmética de enteros, pero un acumulador grande, fue descrito por Kirchner y Kulisch ; ^{[ 17 ]} una implementación de hardware fue descrita por Müller, Rüb y Rülling. ^{[ 18 ]}

En principio, un compilador optimizador suficientemente agresivo podría destruir la efectividad de la suma de Kahan: por ejemplo, si el compilador simplificara las expresiones según las reglas de asociatividad de la aritmética real, podría "simplificar" el segundo paso de la secuencia.

t = sum + y;

c = (t - sum) - y;

a

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

y luego a

c = 0;

eliminando así la compensación de errores. ^{[ 19 ]} En la práctica, muchos compiladores no utilizan reglas de asociatividad (que son solo aproximadas en aritmética de punto flotante) en simplificaciones, a menos que se les indique explícitamente que lo hagan mediante opciones del compilador que habilitan optimizaciones "inseguras", ^{[ 20 ] [ 21 ] [ 22 ] [ 23 ]} aunque el compilador Intel C++ es un ejemplo que permite transformaciones basadas en la asociatividad por defecto. ^{[ 24 ]} La versión original K&R C del lenguaje de programación C permitía al compilador reordenar expresiones de punto flotante de acuerdo con reglas de asociatividad de aritmética real, pero el estándar ANSI C posterior prohibió el reordenamiento para que C fuera más adecuado para aplicaciones numéricas (y más similar a Fortran , que también prohíbe el reordenamiento), ^{[ 25 ]} aunque en la práctica las opciones del compilador pueden volver a habilitar el reordenamiento, como se mencionó anteriormente.

Una forma portátil de inhibir localmente dichas optimizaciones es dividir una de las líneas de la formulación original en dos declaraciones y hacer que dos de los productos intermedios sean volátiles :

función KahanSum(entrada) var suma = 0,0 var c = 0,0 para i = 1 hasta input.length hacer var y = input[i] - c volatile var t = sum + y volatile var z = t - sum c = z - y suma = t a continuación yo suma de retorno

En general, las funciones de suma integradas en los lenguajes de programación no suelen garantizar que se emplee un algoritmo de suma específico, y mucho menos la suma de Kahan. El estándar BLAS para subrutinas de álgebra lineal evita explícitamente imponer un orden de operaciones computacional determinado por motivos de rendimiento ^{[ 26 ]} , y las implementaciones de BLAS no suelen utilizar la suma de Kahan.

La biblioteca estándar del lenguaje de programación Python utiliza la suma de Neumaier en la función integrada "sum()" a partir de CPython 3.12. ^{[ 27 ]} NumPy no garantiza el orden de la suma, aunque la suma parcial por pares se utiliza "a menudo". ^{[ 28 ]}

En el lenguaje Julia , la implementación predeterminada de la sumfunción realiza una suma por pares para obtener una alta precisión con un buen rendimiento, ^{[ 29 ]} pero una biblioteca externa proporciona una implementación de la variante de Neumaier denominada sum_kbnpara los casos en los que se necesita mayor precisión. ^{[ 30 ]}

En el lenguaje C# , el paquete NuGet HPCsharp implementa la variante de Neumaier y la suma por pares : ambas como escalares, en paralelo de datos usando instrucciones de procesador SIMD y en paralelo multinúcleo. ^{[ 31 ]}

• Algoritmos para calcular la varianza , que incluye la suma estable.

• Resumen en coma flotante, Dr. Dobb's Journal, septiembre de 1996