En el análisis numérico , la suma por pares , también llamada suma en cascada , es una técnica para sumar una secuencia de números de punto flotante de precisión finita que reduce sustancialmente el error de redondeo acumulado en comparación con la acumulación ingenua de la suma en secuencia. [1] Aunque existen otras técnicas, como la suma de Kahan, que normalmente tienen errores de redondeo aún más pequeños, la suma por pares es casi tan buena (difiere solo por un factor logarítmico) mientras que tiene un costo computacional mucho menor; se puede implementar para tienen casi el mismo costo (y exactamente el mismo número de operaciones aritméticas) que la suma ingenua.
En particular, la suma por pares de una secuencia de n números x n funciona dividiendo recursivamente la secuencia en dos mitades, sumando cada mitad y sumando las dos sumas: un algoritmo de divide y vencerás . Sus errores de redondeo en el peor de los casos crecen asintóticamente como máximo O (ε log n ), donde ε es la precisión de la máquina (asumiendo un número de condición fijo , como se analiza a continuación). [1] En comparación, la técnica ingenua de acumular la suma en secuencia (sumando cada x i uno a la vez para i = 1, ..., n ) tiene errores de redondeo que, en el peor de los casos, crecen como O (ε n ). [1] La suma de Kahan tiene un error en el peor de los casos de aproximadamente O (ε), independiente de n , pero requiere varias veces más operaciones aritméticas. [1] Si los errores de redondeo son aleatorios y, en particular, tienen signos aleatorios, forman una caminata aleatoria y el crecimiento del error se reduce a un promedio depara la suma por pares. [2]
Una estructura recursiva de suma muy similar se encuentra en muchos algoritmos de transformada rápida de Fourier (FFT), y es responsable de la misma acumulación lenta de redondeo de esas FFT. [2] [3]
El algoritmo
En pseudocódigo , el algoritmo de suma por pares para una matriz x de longitud n > 0 se puede escribir:
s = por pares ( x [1… n ]) si n ≤ N caso base: suma ingenua para una matriz suficientemente pequeña s = x [1] para i = 2 an s = s + x [ i ] de lo contrario divide y vencerás: Sumar recursivamente dos mitades de la matriz m = piso ( n / 2) s = por pares ( x [1… m ]) + por pares ( x [ m + 1… n ]) final si
Para algunos N suficientemente pequeños , este algoritmo cambia a una suma ingenua basada en bucles como caso base , cuyo límite de error es O (Nε). [4] La suma completa tiene un error en el peor de los casos que crece asintóticamente como O (ε log n ) para n grande , para un número de condición dado (ver más abajo).
En un algoritmo de este tipo (como para los algoritmos de divide y vencerás en general [5] ), es deseable utilizar un caso base más grande para amortizar la sobrecarga de la recursividad. Si N = 1, entonces hay aproximadamente una llamada de subrutina recursiva para cada entrada, pero más generalmente hay una llamada recursiva para (más o menos) cada N / 2 entradas si la recursión se detiene en exactamente n = N . Haciendo N lo suficientemente grande, la sobrecarga de la recursividad puede hacerse insignificante (precisamente esta técnica de un caso base grande para la suma recursiva se emplea en implementaciones de FFT de alto rendimiento [3] ).
Independientemente de N , se realizan exactamente n -1 adiciones en total, lo mismo que para la suma ingenua, por lo que si la sobrecarga de recursividad se hace insignificante, la suma por pares tiene esencialmente el mismo costo computacional que para la suma ingenua.
Una variación de esta idea es dividir la suma en b bloques en cada etapa recursiva, sumando cada bloque de forma recursiva y luego sumando los resultados, que sus proponentes denominaron un algoritmo de "superbloque". [6] El algoritmo pairwise anterior corresponde a b = 2 para cada etapa, excepto para la última etapa, que es b = N .
Precisión
Suponga que uno suma n valores x i , para i = 1, ..., n . La suma exacta es:
(calculado con precisión infinita).
Con la suma por pares para un caso base N = 1, se obtiene, donde el error está delimitado por encima de: [1]
donde ε es la precisión de la máquina de la aritmética que se está empleando (por ejemplo, ε − 10 −16 para coma flotante de precisión doble estándar ). Por lo general, la cantidad de interés es el error relativo. , que, por lo tanto, está delimitado anteriormente por:
En la expresión para el límite de error relativo, la fracción (Σ | x i | / | Σ 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. [7] 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. [1] Un problema de suma mal condicionado es aquel en el que esta razón es grande y, en este caso, incluso la suma por pares puede tener un gran error relativo. Por ejemplo, si los sumandos 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. Por otro lado, para entradas aleatorias con distinto de cero significa las asíntotas del número de condición a una constante finita como. Si todas las entradas son no negativas , entonces el número de condición es 1.
Tenga en cuenta que el El denominador es efectivamente 1 en la práctica, ya que es mucho menor que 1 hasta que n pasa a ser de orden 2 1 / ε , que es aproximadamente 10 10 15 en doble precisión.
En comparación, el error relativo ligado a la suma ingenua (simplemente sumando los números en secuencia, redondeando en cada paso) crece a medida que multiplicado por el número de condición. [1] 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 de raíz cuadrada media que crece a medida que y la suma por pares tiene un error que crece a medida que de media. [2]
Implementaciones de software
La suma por pares es el algoritmo de suma predeterminado en NumPy [8] y el lenguaje técnico-informático de Julia , [9] donde en ambos casos se encontró que tenía una velocidad comparable a la suma ingenua (gracias al uso de un caso base grande).
Otras implementaciones de software incluyen la biblioteca HPCsharp [10] para el C agudo lenguaje y la biblioteca sumatoria estándar [11] en el D .
Referencias
- ^ a b c d e f g Higham, Nicholas J. (1993), "La precisión de la suma de coma flotante", SIAM Journal on Scientific Computing , 14 (4): 783–799, CiteSeerX 10.1.1.43.3535 , doi : 10.1137 / 0914050
- ^ a b c Manfred Tasche y Hansmartin Zeuner Manual de métodos analítico-computacionales en matemáticas aplicadas Boca Raton, FL: CRC Press, 2000).
- ^ a b S. G. Johnson y M. Frigo, " Implementación de FFT en la práctica , en Fast Fourier Transforms , editado por C. Sidney Burrus (2008).
- ^ Higham, Nicholas (2002). Precisión y estabilidad de algoritmos numéricos (2 ed) . SIAM. págs. 81–82.
- ^ Radu Rugina y Martin Rinard, " Recursion desenrollado para programas de divide y vencerás ", en Idiomas y compiladores para computación paralela , capítulo 3, págs. 34-48. Notas de la conferencia en Ciencias de la Computación vol. 2017 (Berlín: Springer, 2001).
- ^ Anthony M. Castaldo, R. Clint Whaley y Anthony T. Chronopoulos, "Reducción del error de punto flotante en el producto escalar utilizando la familia de algoritmos de superbloques", SIAM J. Sci. Computación. , vol. 32, págs. 1156-1174 (2008).
- ^ LN Trefethen y D. Bau, Álgebra lineal numérica (SIAM: Filadelfia, 1997).
- ^ ENH: implemente la suma por pares , github.com/numpy/numpy pull request # 3685 (septiembre de 2013).
- ^ RFC: use la suma por pares para sum, cumsum y cumprod , github.com/JuliaLang/julia pull request # 4039 (agosto de 2013).
- ^ https://github.com/DragonSpit/HPCsharp Paquete nuget HPCsharp de algoritmos C # de alto rendimiento
- ^ "std.algorithm.iteration - Lenguaje de programación D" . dlang.org . Consultado el 23 de abril de 2021 .