¿En qué orden se deben agregar los flotadores para obtener el resultado más preciso?

Esta fue una pregunta que me hicieron en mi entrevista reciente y quiero saber (en realidad no recuerdo la teoría del análisis numérico, así que por favor ayúdenme :)

Si tenemos alguna función, que acumula números de punto flotante:

std::accumulate(v.begin(), v.end(), 0.0);

ves un std::vector<float>, por ejemplo.

• ¿Sería mejor ordenar estos números antes de acumularlos?

• ¿Qué orden daría la respuesta más precisa?

Sospecho que la clasificación de los números en orden ascendente en realidad lo haría el error numérico menos , pero por desgracia no puedo demostrar que yo mismo.

PD: Me doy cuenta de que esto probablemente no tiene nada que ver con la programación del mundo real, solo tengo curiosidad.

Su instinto es básicamente correcto, ordenar en orden ascendente (de magnitud) generalmente mejora un poco las cosas. Considere el caso en el que estamos agregando flotadores de precisión simple (32 bits), y hay mil millones de valores iguales a 1 / (mil millones) y un valor igual a 1. Si el 1 es primero, entonces la suma vendrá a 1, ya que 1 + (1/1 billón) es 1 debido a la pérdida de precisión. Cada adición no tiene ningún efecto en el total.

Si los valores pequeños vienen primero, al menos sumarán algo, aunque incluso entonces tengo 2 ^ 30 de ellos, mientras que después de 2 ^ 25 más o menos estoy de vuelta en la situación en la que cada uno individualmente no afecta el total nunca más. Así que todavía voy a necesitar más trucos.

Ese es un caso extremo, pero en general, sumar dos valores de magnitud similar es más exacto que sumar dos valores de magnitudes muy diferentes, ya que "descarta" menos bits de precisión en el valor más pequeño de esa manera. Al ordenar los números, agrupa valores de magnitud similar y, al sumarlos en orden ascendente, le da a los valores pequeños una "posibilidad" de alcanzar acumulativamente la magnitud de los números más grandes.

Aún así, si se trata de números negativos, es fácil "burlar" este enfoque. Considere tres valores Resumiendo, {1, -1, 1 billionth}. La suma aritméticamente correcta es 1 billionth, pero si mi primera adición involucra el valor minúsculo, entonces mi suma final será 0. De los 6 órdenes posibles, solo 2 son "correctos" - {1, -1, 1 billionth}y {-1, 1, 1 billionth}. Los 6 órdenes dan resultados que son precisos en la escala del valor de mayor magnitud en la entrada (0,0000001% fuera), pero para 4 de ellos el resultado es inexacto en la escala de la solución verdadera (100% fuera). El problema particular que está resolviendo le dirá si el primero es lo suficientemente bueno o no.

De hecho, puede jugar muchos más trucos que simplemente agregarlos en orden ordenado. Si tiene muchos valores muy pequeños, un número medio de valores medios y un número pequeño de valores grandes, entonces podría ser más exacto sumar primero todos los pequeños, luego sumar los medianos por separado, sumar esos dos totales juntos luego agregue los grandes. No es en absoluto trivial encontrar la combinación más precisa de adiciones de punto flotante, pero para hacer frente a casos realmente malos, puede mantener un conjunto completo de totales acumulados en diferentes magnitudes, agregar cada nuevo valor al total que mejor coincida con su magnitud, y cuando un total acumulado comience a ser demasiado grande para su magnitud, agréguelo al siguiente total y comience uno nuevo. Llevado a su extremo lógico, este proceso es equivalente a realizar la suma en un tipo de precisión arbitraria (por lo que ' haría eso). Pero dada la opción simplista de sumar en orden de magnitud ascendente o descendente, ascender es la mejor apuesta.

Tiene alguna relación con la programación del mundo real, ya que hay algunos casos en los que su cálculo puede salir muy mal si accidentalmente corta una cola "pesada" que consiste en una gran cantidad de valores, cada uno de los cuales es demasiado pequeño para afectar individualmente la suma, o si descarta demasiada precisión de una gran cantidad de valores pequeños que individualmente solo afectan a los últimos bits de la suma. En los casos en que la cola es insignificante de todos modos, probablemente no le importe. Por ejemplo, si solo está sumando una pequeña cantidad de valores en primer lugar y solo está usando algunas cifras significativas de la suma.

También hay un algoritmo diseñado para este tipo de operación de acumulación, llamado Kahan Summation , que probablemente debería conocer.

Según Wikipedia,

El algoritmo de suma de Kahan (también conocido como suma compensada ) reduce significativamente el error numérico en el total obtenido al agregar 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).

En pseudocódigo, el algoritmo es:

function kahanSum(input)
 var sum = input[1]
 var c = 0.0          //A running compensation for lost low-order bits.
 for i = 2 to input.length
  y = input[i] - c    //So far, so good: c is zero.
  t = sum + y         //Alas, sum is big, y small, so low-order digits of y are lost.
  c = (t - sum) - y   //(t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
  sum = t             //Algebraically, c should always be zero. Beware eagerly optimising compilers!
 next i               //Next time around, the lost low part will be added to y in a fresh attempt.
return sum

Probé el ejemplo extremo en la respuesta proporcionada por Steve Jessop.

#include <iostream>
#include <iomanip>
#include <cmath>

int main()
{
    long billion = 1000000000;
    double big = 1.0;
    double small = 1e-9;
    double expected = 2.0;

    double sum = big;
    for (long i = 0; i < billion; ++i)
        sum += small;
    std::cout << std::scientific << std::setprecision(1) << big << " + " << billion << " * " << small << " = " <<
        std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    sum = 0;
    for (long i = 0; i < billion; ++i)
        sum += small;
    sum += big;
    std::cout  << std::scientific << std::setprecision(1) << billion << " * " << small << " + " << big << " = " <<
        std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    return 0;
}

Obtuve el siguiente resultado:

1.0e+00 + 1000000000 * 1.0e-09 = 2.000000082740371    (difference = 0.000000082740371)
1000000000 * 1.0e-09 + 1.0e+00 = 1.999999992539933    (difference = 0.000000007460067)

El error en la primera línea es más de diez veces mayor en la segunda.

Si cambio la doublesa floats en el código anterior, obtengo:

1.0e+00 + 1000000000 * 1.0e-09 = 1.000000000000000    (difference = 1.000000000000000)
1000000000 * 1.0e-09 + 1.0e+00 = 1.031250000000000    (difference = 0.968750000000000)

Ninguna de las respuestas se acerca siquiera a 2.0 (pero la segunda está un poco más cerca).

Usando la suma de Kahan (con doubles) como lo describe Daniel Pryden:

#include <iostream>
#include <iomanip>
#include <cmath>

int main()
{
    long billion = 1000000000;
    double big = 1.0;
    double small = 1e-9;
    double expected = 2.0;

    double sum = big;
    double c = 0.0;
    for (long i = 0; i < billion; ++i) {
        double y = small - c;
        double t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }

    std::cout << "Kahan sum  = " << std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    return 0;
}

Obtengo exactamente 2.0:

Kahan sum  = 2.000000000000000    (difference = 0.000000000000000)

E incluso si cambio las doubles por floats en el código anterior, obtengo:

Kahan sum  = 2.000000000000000    (difference = 0.000000000000000)

¡Parece que Kahan es el camino a seguir!

Existe una clase de algoritmos que resuelven este problema exacto, sin la necesidad de ordenar o reordenar los datos .

En otras palabras, la suma se puede realizar en una pasada sobre los datos. Esto también hace que dichos algoritmos sean aplicables en situaciones en las que el conjunto de datos no se conoce de antemano, por ejemplo, si los datos llegan en tiempo real y es necesario mantener la suma acumulada.

Aquí está el resumen de un artículo reciente:

Presentamos un novedoso algoritmo en línea para la suma exacta de un flujo de números de punto flotante. Por "en línea" queremos decir que el algoritmo necesita ver solo una entrada a la vez, y puede tomar un flujo de entrada de longitud arbitraria de tales entradas mientras requiere solo memoria constante. Por "exacta" queremos decir que la suma de la matriz interna de nuestro algoritmo es exactamente igual a la suma de todas las entradas, y el resultado devuelto es la suma correctamente redondeada. La prueba de corrección es válida para todas las entradas (incluidos los números no normalizados pero el desbordamiento intermedio del módulo) y es independiente del número de sumandos o del número de condición de la suma. El algoritmo necesita asintóticamente solo 5 FLOP por suma y, debido al paralelismo a nivel de instrucción, se ejecuta solo entre 2 y 3 veces más lento que lo obvio, Bucle de “suma recursiva ordinaria” rápido pero tonto cuando el número de sumandos es mayor que 10,000. Por lo tanto, hasta donde sabemos, es el más rápido, más preciso y más eficiente en memoria de todos los algoritmos conocidos. De hecho, es difícil ver cómo podría existir un algoritmo más rápido o uno que requiera significativamente menos FLOP sin mejoras de hardware. Se proporciona una aplicación para una gran cantidad de sumandos.

Fuente: Algoritmo 908: Suma exacta en línea de corrientes de punto flotante .

Sobre la base de la respuesta de Steve de ordenar primero los números en orden ascendente, presentaré dos ideas más:

1. Decida la diferencia en exponente de dos números por encima del cual podría decidir que perdería demasiada precisión.

2. Luego sume los números en orden hasta que el exponente del acumulador sea demasiado grande para el siguiente número, luego coloque el acumulador en una cola temporal y comience el acumulador con el siguiente número. Continúe hasta agotar la lista original.

Repite el proceso con la cola temporal (habiéndola ordenado) y con una diferencia posiblemente mayor en el exponente.

Creo que esto será bastante lento si tienes que calcular exponentes todo el tiempo.

Probé rápidamente un programa y el resultado fue 1.99903

Creo que puedes hacerlo mejor que ordenar los números antes de acumularlos, porque durante el proceso de acumulación, el acumulador se hace cada vez más grande. Si tiene una gran cantidad de números similares, comenzará a perder precisión rápidamente. Esto es lo que sugeriría en su lugar:

while the list has multiple elements
    remove the two smallest elements from the list
    add them and put the result back in
the single element in the list is the result

Por supuesto, este algoritmo será más eficiente con una cola de prioridad en lugar de una lista. Código C ++:

template <typename Queue>
void reduce(Queue& queue)
{
    typedef typename Queue::value_type vt;
    while (queue.size() > 1)
    {
        vt x = queue.top();
        queue.pop();
        vt y = queue.top();
        queue.pop();
        queue.push(x + y);
    }
}

conductor:

#include <iterator>
#include <queue>

template <typename Iterator>
typename std::iterator_traits<Iterator>::value_type
reduce(Iterator begin, Iterator end)
{
    typedef typename std::iterator_traits<Iterator>::value_type vt;
    std::priority_queue<vt> positive_queue;
    positive_queue.push(0);
    std::priority_queue<vt> negative_queue;
    negative_queue.push(0);
    for (; begin != end; ++begin)
    {
        vt x = *begin;
        if (x < 0)
        {
            negative_queue.push(x);
        }
        else
        {
            positive_queue.push(-x);
        }
    }
    reduce(positive_queue);
    reduce(negative_queue);
    return negative_queue.top() - positive_queue.top();
}

Los números en la cola son negativos porque topproduce el número más grande , pero queremos el más pequeño . Podría haber proporcionado más argumentos de plantilla a la cola, pero este enfoque parece más simple.

Esto no responde del todo a su pregunta, pero una cosa inteligente que puede hacer es ejecutar la suma dos veces, una con el modo de redondeo "redondear hacia arriba" y una vez con "redondear hacia abajo". Compare las dos respuestas y sabrá / cómo / inexactos son sus resultados, y si por lo tanto necesita utilizar una estrategia de suma más inteligente. Desafortunadamente, la mayoría de los lenguajes no hacen que cambiar el modo de redondeo de coma flotante sea tan fácil como debería ser, porque la gente no sabe que es realmente útil en los cálculos diarios.

Eche un vistazo a la aritmética de intervalos, donde hace todas las matemáticas como esta, manteniendo los valores más altos y más bajos a medida que avanza. Conduce a algunas optimizaciones y resultados interesantes.

El más simple tipo que mejora la precisión es para ordenar por el valor absoluto ascendente. Eso permite que los valores de magnitud más pequeños tengan la oportunidad de acumularse o cancelarse antes de interactuar con valores de magnitud más grandes que provocarían una pérdida de precisión.

Dicho esto, puede hacerlo mejor si realiza un seguimiento de varias sumas parciales que no se superponen. Aquí hay un artículo que describe la técnica y presenta una prueba de precisión: www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps

Ese algoritmo y otros enfoques para la suma exacta de punto flotante se implementan en Python simple en: http://code.activestate.com/recipes/393090/ Al menos dos de ellos se pueden convertir trivialmente a C ++.

Para los números de formato conocido o de precisión simple o doble IEEE 754, otra alternativa es usar una matriz de números (pasados ​​por el llamador o en una clase para C ++) indexados por el exponente. Al agregar números a la matriz, solo se agregan números con el mismo exponente (hasta que se encuentra un espacio vacío y se almacena el número). Cuando se solicita una suma, la matriz se suma de menor a mayor para minimizar el truncamiento. Ejemplo de precisión simple:

/* clear array */
void clearsum(float asum[256])
{
size_t i;
    for(i = 0; i < 256; i++)
        asum[i] = 0.f;
}

/* add a number into array */
void addtosum(float f, float asum[256])
{
size_t i;
    while(1){
        /* i = exponent of f */
        i = ((size_t)((*(unsigned int *)&f)>>23))&0xff;
        if(i == 0xff){          /* max exponent, could be overflow */
            asum[i] += f;
            return;
        }
        if(asum[i] == 0.f){     /* if empty slot store f */
            asum[i] = f;
            return;
        }
        f += asum[i];           /* else add slot to f, clear slot */
        asum[i] = 0.f;          /* and continue until empty slot */
    }
}

/* return sum from array */
float returnsum(float asum[256])
{
float sum = 0.f;
size_t i;
    for(i = 0; i < 256; i++)
        sum += asum[i];
    return sum;
}

ejemplo de doble precisión:

/* clear array */
void clearsum(double asum[2048])
{
size_t i;
    for(i = 0; i < 2048; i++)
        asum[i] = 0.;
}

/* add a number into array */
void addtosum(double d, double asum[2048])
{
size_t i;
    while(1){
        /* i = exponent of d */
        i = ((size_t)((*(unsigned long long *)&d)>>52))&0x7ff;
        if(i == 0x7ff){         /* max exponent, could be overflow */
            asum[i] += d;
            return;
        }
        if(asum[i] == 0.){      /* if empty slot store d */
            asum[i] = d;
            return;
        }
        d += asum[i];           /* else add slot to d, clear slot */
        asum[i] = 0.;           /* and continue until empty slot */
    }
}

/* return sum from array */
double returnsum(double asum[2048])
{
double sum = 0.;
size_t i;
    for(i = 0; i < 2048; i++)
        sum += asum[i];
    return sum;
}

Sus flotadores deben agregarse con doble precisión. Eso le dará más precisión adicional que cualquier otra técnica. Para un poco más de precisión y significativamente más velocidad, puede crear, digamos, cuatro sumas y sumarlas al final.

Si está agregando números de doble precisión, use long double para la suma; sin embargo, esto solo tendrá un efecto positivo en implementaciones donde long double en realidad tiene más precisión que double (típicamente x86, PowerPC dependiendo de la configuración del compilador).

Con respecto a la clasificación, me parece que si espera una cancelación, los números deben agregarse en orden descendente de magnitud, no ascendente. Por ejemplo:

((-1 + 1) + 1e-20) dará 1e-20

pero

((1e-20 + 1) - 1) dará 0

En la primera ecuación se anulan dos números grandes, mientras que en la segunda el término 1e-20 se pierde cuando se suma a 1, ya que no hay suficiente precisión para retenerlo.

Además, la suma por pares es bastante decente para sumar muchos números.