stringtranslate.com

Algoritmo de suma de Kahan

En el 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 punto flotante de precisión finita , en comparación con el enfoque obvio. Esto se hace manteniendo una compensación de ejecución separada (una variable para acumular pequeños errores), lo que en efecto extiende la precisión de la suma por la precisión de la variable de compensación.

En particular, simplemente sumar números en secuencia tiene un error en el peor de los casos que crece proporcionalmente a , y un error cuadrático medio que crece como 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 del error en el peor de los casos es efectivamente independiente de , 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 realiza un seguimiento del error acumulado en operaciones con números enteros (aunque documentado por primera vez aproximadamente en la misma época [5] ) y la modulación delta-sigma . [6]

El algoritmo

En pseudocódigo , el algoritmo será:

función KahanSum(entrada) //Preparar el acumulador. var suma = 0.0 // Una compensación en ejecución para los bits de orden bajo perdidos. var c = 0.0 // La matriz input tiene elementos indexados input[1] a input[input.length]. for i = 1 a input.length do // c es cero la primera vez. var y = input[i] - c // Lamentablemente, la suma es grande, y pequeña, por lo que se pierden los dígitos de orden inferior de y . var t = suma + y // (t - suma) cancela la parte de orden superior de y ; // restando y recupera el negativo (parte baja de y ) c = (t - suma) - y // Algebraicamente, c siempre debería ser cero. ¡Cuidado! // ¡compiladores optimizadores demasiado agresivos! suma = t// La próxima vez, la parte baja perdida se agregará a y en un nuevo intento.  devolver suma

Este algoritmo también se puede reescribir para utilizar el algoritmo Fast2Sum : [7]

función KahanSum2(entrada) //Preparar el acumulador. var suma = 0.0 // Una compensación en ejecución para los bits de orden bajo perdidos. var c = 0.0 // La matriz input tiene elementos indexados para i = 1 en 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. A continuación, devuelvo la suma .

Ejemplo resuelto

El algoritmo no exige ninguna elección específica de base , solo que la aritmética "normalice las sumas de punto flotante antes de redondear o truncar". [3] Las computadoras suelen utilizar aritmética binaria, pero para que el ejemplo sea más fácil de leer, se dará en decimal. Supongamos que estamos utilizando aritmética de punto flotante decimal de seis dígitos , sumha alcanzado el valor 10000.0, y los dos siguientes valores de input[i]son 3.14159 y 2.71828. El resultado exacto es 10005.85987, que se redondea a 10005.9. Con una suma simple, cada valor entrante se alinearía con sum, y se perderían muchos dígitos de orden bajo (por truncamiento o redondeo). El primer resultado, después del redondeo, sería 10003.1. El segundo resultado sería 10005,81828 antes del redondeo y 10005,8 después del redondeo. Esto no es correcto.

Sin embargo, con la suma compensada, obtenemos el resultado correctamente redondeado de 10005,9.

Supongamos que cel valor inicial es cero. Los ceros finales se muestran cuando son significativos para el número de coma flotante de seis dígitos.

 y = 3,14159 - 0,00000 y = entrada[i] - c t = 10000,0 + 3,14159 t = suma + y = 10003.14159 Normalización realizada, siguiente redondeo a seis dígitos. = 10003.1 Algunos dígitos de input[i] coinciden con los de sum . ¡Se han perdido muchos dígitos! c = (10003,1 - 10000,0) - 3,14159 c = (t - suma) - y (Nota: ¡Los paréntesis deben evaluarse primero!) = 3,10000 - 3,14159 La parte asimilada de y menos la y original completa . = -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

La suma es tan grande que solo se acumulan los dígitos de orden superior de los números de entrada. Pero en el siguiente paso, c, una aproximación del error de ejecución, contrarresta el problema.

 y = 2,71828 - (-0,0415900) La mayoría de los dígitos coinciden, ya que c es de un tamaño similar a y . = 2.75987 El déficit (dígitos de orden inferior perdidos) de la iteración anterior se restableció exitosamente. t = 10003,1 + 2,75987 Pero aún así sólo unos pocos cumplen 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 acelerar el redondeo.c = (10005,9 - 10003,1) - 2,75987 Estime el error acumulado, con base en la y ajustada . = 2,80000 - 2,75987 Como era de esperar, las partes de orden bajo se pueden conservar en c sin efectos de redondeo o con efectos de redondeo mínimos. = 0.0401300 En esta iteración, t fue un poco demasiado alto, el exceso se restará en la próxima iteración.suma = 10005.9 El resultado exacto es 10005.85987, la suma es correcta, redondeada a 6 dígitos.

El algoritmo realiza la suma con dos acumuladores: sumretiene la suma y cacumula las partes no asimiladas en sum, para empujar la parte de orden inferior de sumla siguiente vez. Por lo tanto, la suma procede con "dígitos de guarda" en c, lo cual es mejor que no tener ninguno, pero no es tan bueno como realizar los cálculos con el doble de precisión de la entrada. Sin embargo, simplemente aumentar la precisión de los cálculos no es práctico en general; si inputya está en doble precisión, pocos sistemas proporcionan precisión cuádruple , y si lo hicieran, inputpodrían estar en precisión cuádruple.

Exactitud

Es necesario realizar 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 simple, aún puede arrojar grandes errores relativos para sumas mal condicionadas.

Supongamos que se suman valores , para . La suma exacta es

(calculado con precisión infinita).

Con la suma compensada, en cambio, se obtiene , donde el error está limitado por [2]

donde es la precisión de la máquina de la aritmética que se está empleando (por ejemplo, para el punto flotante de doble precisión estándar IEEE ). Por lo general, la cantidad de interés es el error relativo , que por lo tanto está limitado por encima por

En la expresión para el límite de error relativo, la fracción 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 ) por un algoritmo fijo en precisión fija (es decir, no aquellos que usan aritmética de precisión arbitraria , ni algoritmos cuyos requisitos de memoria y tiempo cambian según los datos), es proporcional a este número de condición. [2] Un problema de suma mal condicionado es uno en el que esta relación es grande, y en este caso incluso la suma compensada puede tener un gran error relativo. Por ejemplo, si los sumandos son números aleatorios no correlacionados con media cero, la suma es un paseo aleatorio , y el número de condición crecerá proporcional a . Por otro lado, para entradas aleatorias con media distinta de cero, el número de condición asíntota a una constante finita como . Si las entradas son todas 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 . En principio, existe el que crece linealmente con , pero en la práctica este término es efectivamente cero: dado que el resultado final se redondea a una precisión , el término se redondea a cero, a menos que sea aproximadamente o mayor. [2] En precisión doble, esto corresponde a un de aproximadamente , mucho mayor 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 , independientes de .

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 al multiplicarse por el número de condición. [2] Sin embargo, este error del peor caso rara vez se observa en la práctica, porque solo ocurre si los errores de redondeo están todos 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 al multiplicarse por el número de condición. [9] Sin embargo, esto sigue siendo mucho peor que la suma compensada. Sin embargo, si la suma se puede realizar con el doble de precisión, entonces se reemplaza por , y la suma ingenua tiene un error del peor caso comparable al término en la suma compensada con la precisión original.

De la misma manera, lo que aparece en 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 se reemplazan por un paseo aleatorio, en cuyo caso, incluso para entradas aleatorias con media cero, el error crece solo como (ignorando el término), la misma tasa en que crece la suma , cancelando los factores cuando se calcula el error relativo. Entonces, incluso para sumas asintóticamente mal condicionadas, el error relativo para la suma compensada a menudo puede ser mucho menor de lo que un análisis del peor caso podría sugerir.

Mejoras adicionales

Neumaier [10] introdujo una versión mejorada del algoritmo de Kahan, al que llama "algoritmo de Kahan–Babuška mejorado", que también cubre el caso en el que el siguiente término que se suma es mayor en valor absoluto que la suma corriente, intercambiando efectivamente el papel de lo que es grande y lo que es pequeño. En pseudocódigo , el algoritmo es:

función KahanBabushkaNeumaierSum(input) var sum = 0.0 var c = 0.0 // Una compensación en ejecución para los bits de orden bajo perdidos. para i = 1 a input.length hacer  var t = suma + input[i] si |sum| >= |input[i]| entonces c += (sum - t) + input[i] // Si suma es mayor, se pierden los dígitos de orden inferior de input[i] . de lo contrario c += (input[i] - t) + suma // De lo contrario, se pierden los dígitos de orden inferior de suma . endif suma = t el siguiente yo devuelve suma + c // La corrección solo se aplicó una vez al final.

Esta mejora es similar al reemplazo de Fast2Sum por 2Sum en el algoritmo de Kahan reescrito con Fast2Sum.

Para muchas secuencias de números, ambos algoritmos coinciden, pero un ejemplo sencillo de Peters [11] muestra cómo pueden diferir. Para la suma con precisión doble, el algoritmo de Kahan da como resultado 0,0, mientras que el algoritmo de Neumaier da como resultado el valor correcto 2,0.

También son posibles modificaciones de orden superior para lograr una mayor precisión. Por ejemplo, una variante sugerida por Klein [12] , a la que llamó "algoritmo iterativo de Kahan-Babuška" de segundo orden. En pseudocódigo , el algoritmo es:

función KahanBabushkaKleinSum(entrada) var suma = 0,0 var cs = 0,0 var ccs = 0,0 var c = 0,0 var cc = 0,0 para i = 1 a entrada.longitud hacer  var t = suma + entrada[i] si |suma| >= |entrada[i]| entonces c = (suma - t) + entrada[i] demás c = (entrada[i] - t) + suma fin si suma = t t = cs + c si |cs| >= |c| entonces c.c. = (c.s. - t.) + c. demás c.c. = (c.-t.) + c.s. fin si cs = t ccs = ccs + cc bucle final  devuelve suma + cs + ccs

Alternativas

Aunque el algoritmo de Kahan logra un crecimiento de error para sumar n números, solo se puede lograr un crecimiento ligeramente peor con la suma por pares : uno divide recursivamente el conjunto de números en dos mitades, suma cada mitad y luego suma 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 de cuatro veces una suma simple) y se puede calcular en paralelo. El caso base de la recursión podría, en principio, ser la suma de solo un número (o cero), pero para amortizar la sobrecarga de la recursión, normalmente se usaría un caso base más grande. 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. [13] En la práctica, con errores de redondeo de signos aleatorios, los errores cuadráticos medios de la suma por pares en realidad crecen como . [9]

Otra alternativa es utilizar aritmética de precisión arbitraria , que en principio no necesita redondeo alguno con un coste de esfuerzo computacional mucho mayor. Una forma de realizar sumas redondeadas correctamente utilizando precisión arbitraria es extender de forma adaptativa utilizando múltiples componentes de punto flotante. Esto minimizará el coste computacional en casos comunes en los que no se necesita alta precisión. [14] [11] Kirchner y Kulisch describieron otro método que utiliza solo aritmética de números enteros, pero un acumulador grande ; [15] Müller, Rüb y Rülling describieron una implementación de hardware. [16]

Posible invalidación por optimización del compilador

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. [17] En la práctica, muchos compiladores no utilizan reglas de asociatividad (que solo son aproximadas en aritmética de punto flotante) en simplificaciones, a menos que se indique explícitamente que lo hagan las opciones del compilador que habilitan optimizaciones "inseguras", [18] [19] [20] [21] aunque el compilador Intel C++ es un ejemplo que permite transformaciones basadas en asociatividad de forma predeterminada. [22] 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 prohibía la reordenación para hacer que C fuera más adecuado para aplicaciones numéricas (y más similar a Fortran , que también prohíbe la reordenación), [23] aunque en la práctica las opciones del compilador pueden volver a habilitar la reordenación, como se mencionó anteriormente.

Una forma portátil de inhibir dichas optimizaciones localmente 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 a entrada.longitud hacer  var y = entrada[i] - c volátil var t = suma + y volátil var z = t - suma c = z - y suma = t el siguiente yo devolver suma

Apoyo de las bibliotecas

En general, las funciones de "suma" integradas en los lenguajes informáticos normalmente no ofrecen garantías de que se empleará un algoritmo de suma particular, mucho menos la suma de Kahan. [ cita requerida ] El estándar BLAS para subrutinas de álgebra lineal evita explícitamente exigir un orden computacional particular de operaciones por razones de rendimiento, [24] y las implementaciones de BLAS normalmente no utilizan la suma de Kahan.

La biblioteca estándar del lenguaje de programación Python especifica una función fsum para realizar sumas precisas. A partir de Python 3.12, la función "sum()" incorporada utiliza la suma de Neumaier. [25]

En el lenguaje Julia , la implementación predeterminada de la sumfunción realiza una suma por pares para lograr una alta precisión con un buen rendimiento, [26] pero una biblioteca externa proporciona una implementación de la variante de Neumaier llamada así sum_kbnpara los casos en los que se necesita una mayor precisión. [27]

En el lenguaje C# , el paquete nuget HPCsharp implementa la variante Neumaier y la suma por pares : tanto como escalar, como datos paralelos usando instrucciones de procesador SIMD y núcleos múltiples paralelos. [28]

Véase también

Referencias

  1. ^ En sentido estricto, también existen otras variantes de la suma compensada: véase Higham, Nicholas (2002). Accuracy and Stability of Numerical Algorithms (2.ª edición) . SIAM. pp. 110–123. ISBN. 978-0-89871-521-7.
  2. ^ abcdefgh Higham, Nicholas J. (1993), "La precisión de la suma de coma flotante", SIAM Journal on Scientific Computing , 14 (4): 783–799, Bibcode :1993SJSC...14..783H, CiteSeerX 10.1.1.43.3535 , doi :10.1137/0914050, S2CID  14071038 .
  3. ^ ab Kahan, William (enero de 1965), "Observaciones adicionales sobre la reducción de errores de truncamiento" (PDF) , Communications of the ACM , 8 (1): 40, doi :10.1145/363707.363723, S2CID  22584810, archivado desde el original (PDF) el 9 de febrero de 2018.
  4. ^ Babuska, I.: Estabilidad numérica en el análisis matemático. Inf. Proc. ˇ 68, 11–23 (1969)
  5. ^ Bresenham, Jack E. (enero de 1965). "Algoritmo para el control informático de un trazador digital" (PDF) . IBM Systems Journal . 4 (1): 25–30. doi :10.1147/sj.41.0025. S2CID  41898371.
  6. ^ Inose, H.; Yasuda, Y.; Murakami, J. (septiembre de 1962). "Un sistema de telemetría mediante manipulación de código: modulación ΔΣ". Transacciones IRE sobre electrónica espacial y telemetría . SET-8: 204–209. doi :10.1109/IRET-SET.1962.5008839. S2CID  51647729.
  7. ^ Müller, Jean-Michel; Brunie, Nicolás; de Dinechin, Florent; Jeannerod, Claude-Pierre; Joldes, Mioara; Lefèvre, Vicente; Melquiond, Guillaume; Revol, Nathalie ; Torres, Serge (2018) [2010]. Manual de aritmética de coma flotante (2 ed.). Birkhäuser . pag. 179. doi :10.1007/978-3-319-76526-6. ISBN 978-3-319-76525-9. Número de serie LCCN  2018935254.
  8. ^ Trefethen, Lloyd N. ; Bau, David (1997). Álgebra lineal numérica . Filadelfia: SIAM. ISBN 978-0-89871-361-9.
  9. ^ por Manfred Tasche y Hansmartin Zeuner, Manual de métodos analíticos y computacionales en matemáticas aplicadas , Boca Raton, FL: CRC Press, 2000.
  10. ^ Neumaier, A. (1974). "Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen" [Análisis de errores de redondeo de algunos métodos para sumar sumas finitas] (PDF) . Zeitschrift für Angewandte Mathematik und Mechanik (en alemán). 54 (1): 39–51. Código bibliográfico : 1974ZaMM...54...39N. doi :10.1002/zamm.19740540106.
  11. ^ ab Hettinger, R. "Mejorar la precisión de la función sum() incorporada para entradas de coma flotante · Problema n.° 100425 · python/cpython". GitHub - Funciones añadidas de CPython v3.12 . Consultado el 7 de octubre de 2023 .
  12. ^ A., Klein (2006). "Un algoritmo de suma generalizado de Kahan–Babuška". Computing . 76 (3–4). Springer-Verlag: 279–293. doi :10.1007/s00607-005-0139-x. S2CID  4561254.
  13. ^ Johnson, SG; Frigo, MC Sidney Burns (ed.). "Transformadas rápidas de Fourier: implementación de FFT en la práctica". Archivado desde el original el 20 de diciembre de 2008.
  14. ^ Richard Shewchuk, Jonathan (octubre de 1997). "Aritmética de punto flotante de precisión adaptativa y predicados geométricos rápidos y robustos" (PDF) . Geometría discreta y computacional . 18 (3): 305–363. doi :10.1007/PL00009321. S2CID  189937041.
  15. ^ Kirchner, R.; Kulisch, U. (junio de 1988). "Aritmética precisa para procesadores vectoriales". Journal of Parallel and Distributed Computing . 5 (3): 250–270. doi :10.1016/0743-7315(88)90020-2.
  16. ^ Muller, M.; Rub, C.; Rulling, W. (1991). Acumulación exacta de números de punto flotante. Actas del 10.º Simposio IEEE sobre aritmética informática. págs. 64–69. doi :10.1109/ARITH.1991.145535.
  17. ^ Goldberg, David (marzo de 1991), "Lo que todo científico informático debería saber sobre aritmética de punto flotante" (PDF) , ACM Computing Surveys , 23 (1): 5–48, doi :10.1145/103162.103163, S2CID  222008826.
  18. ^ Manual de la colección del compilador GNU , versión 4.4.3: 3.10 Opciones que controlan la optimización, -fassociative-math (21 de enero de 2010).
  19. ^ Manual del usuario de Compaq Fortran para sistemas Tru64 UNIX y Linux Alpha Archivado el 7 de junio de 2011 en Wayback Machine , sección 5.9.7 Optimizaciones de reordenamiento aritmético (consultado en marzo de 2010).
  20. ^ Börje Lindh, Optimización del rendimiento de aplicaciones, Sun BluePrints OnLine (marzo de 2002).
  21. ^ Eric Fleegal, "Optimización de punto flotante de Microsoft Visual C++", Artículos técnicos de Microsoft Visual Studio (junio de 2004).
  22. ^ Martyn J. Corden, "Consistencia de resultados de punto flotante utilizando el compilador Intel", informe técnico de Intel (18 de septiembre de 2009).
  23. ^ MacDonald, Tom (1991). "C para computación numérica". Revista de supercomputación . 5 (1): 31–48. doi :10.1007/BF00155856. S2CID  27876900.
  24. ^ Foro técnico de BLAS, sección 2.7 (21 de agosto de 2001), archivado en Wayback Machine.
  25. ^ Novedades en Python 3.12.
  26. ^ RFC: utilizar suma por pares para suma, cumsum y cumprod, solicitud de extracción n.° 4039 de github.com/JuliaLang/julia (agosto de 2013).
  27. ^ Biblioteca KahanSummation en Julia.
  28. ^ Paquete nuget HPCsharp de algoritmos de alto rendimiento.

Enlaces externos