Em que ordem os flutuadores devem ser adicionados para obter o resultado mais preciso?

Esta foi uma pergunta que me fizeram em minha entrevista recente e eu quero saber (na verdade não me lembro da teoria da análise numérica, então por favor me ajude :)

Se tivermos alguma função, que acumula números de ponto flutuante:

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

vé um std::vector<float>, por exemplo.

• Seria melhor classificar esses números antes de acumulá-los?

• Qual ordem daria a resposta mais precisa?

Eu suspeito que classificar os números em ordem ascendente iria realmente fazer o erro numérico menos , mas infelizmente não posso provar isso a mim mesmo.

PS: Eu sei que isso provavelmente não tem nada a ver com a programação do mundo real, apenas estou curioso.

Seu instinto está basicamente certo, a classificação em ordem crescente (de magnitude) geralmente melhora as coisas um pouco. Considere o caso em que estamos adicionando flutuadores de precisão simples (32 bits) e há 1 bilhão de valores iguais a 1 / (1 bilhão) e um valor igual a 1. Se o 1 vier primeiro, a soma virá para 1, visto que 1 + (1/1 bilhão) é 1 devido à perda de precisão. Cada adição não tem efeito algum no total.

Se os valores pequenos vierem primeiro, eles pelo menos somarão alguma coisa, embora mesmo assim eu tenha 2 ^ 30 deles, enquanto depois de 2 ^ 25 ou mais estou de volta à situação em que cada um individualmente não está afetando o total não mais. Ainda vou precisar de mais truques.

Esse é um caso extremo, mas em geral adicionar dois valores de magnitude semelhante é mais preciso do que adicionar dois valores de magnitudes muito diferentes, já que você "descarta" menos bits de precisão no valor menor dessa forma. Classificando os números, você agrupa valores de magnitude semelhante e, ao adicioná-los em ordem crescente, dá aos valores pequenos uma "chance" de atingir cumulativamente a magnitude dos números maiores.

Ainda assim, se números negativos estiverem envolvidos, é fácil "enganar" essa abordagem. Considere três valores para somar {1, -1, 1 billionth},. A soma aritmeticamente correta é 1 billionth, mas se minha primeira adição envolver o valor minúsculo, minha soma final será 0. Das 6 ordens possíveis, apenas 2 são "corretas" - {1, -1, 1 billionth}e {-1, 1, 1 billionth}. Todas as 6 ordens fornecem resultados que são precisos na escala do valor de maior magnitude na entrada (0,0000001% de saída), mas para 4 delas o resultado é impreciso na escala da solução verdadeira (100% de saída). O problema específico que você está resolvendo dirá se o primeiro é bom o suficiente ou não.

Na verdade, você pode fazer muito mais truques do que apenas adicioná-los em ordem. Se você tiver muitos valores muito pequenos, um número médio de valores médios e um pequeno número de valores grandes, então pode ser mais preciso primeiro somar todos os pequenos e, em seguida, somar separadamente os médios, adicionar esses dois totais juntos, em seguida, adicione os grandes. Não é nada trivial encontrar a combinação mais precisa de adições de ponto flutuante, mas para lidar com casos realmente ruins, você pode manter toda uma matriz de totais em execução em diferentes magnitudes, adicionar cada novo valor ao total que melhor corresponda à sua magnitude, e quando um total corrente começar a ficar muito grande para sua magnitude, some-o ao próximo total e comece um novo. Levado ao extremo lógico, este processo é equivalente a realizar a soma em um tipo de precisão arbitrária (então você ' d fazer isso). Mas dada a escolha simplista de adicionar ordem de magnitude ascendente ou descendente, ascender é a melhor aposta.

Ele tem alguma relação com a programação do mundo real, já que há alguns casos em que seu cálculo pode dar muito errado se você acidentalmente cortar uma cauda "pesada" consistindo de um grande número de valores, cada um dos quais é muito pequeno para afetar individualmente a soma ou se você descartar precisão demais de muitos valores pequenos que individualmente afetam apenas os últimos bits da soma. Nos casos em que a cauda é insignificante, você provavelmente não se importa. Por exemplo, se você estiver apenas adicionando um pequeno número de valores em primeiro lugar e estiver usando apenas alguns algarismos significativos da soma.

Também existe um algoritmo projetado para esse tipo de operação de acumulação, chamado Kahan Summation , do qual você provavelmente deve estar ciente.

De acordo com a Wikipedia,

O algoritmo de soma de Kahan (também conhecido como soma compensada ) reduz significativamente o erro numérico no total obtido pela adição de uma sequência de números de ponto flutuante de precisão finita, em comparação com a abordagem óbvia. Isso é feito mantendo uma compensação de funcionamento separada (uma variável para acumular pequenos erros).

Em pseudocódigo, o algoritmo é:

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

Experimentei o exemplo extremo na resposta fornecida 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;
}

Obtive o seguinte 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)

O erro na primeira linha é mais de dez vezes maior na segunda.

Se eu alterar doubles para floats no código acima, obtenho:

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

Nenhuma das respostas está nem perto de 2.0 (mas a segunda está um pouco mais perto).

Usando o somatório Kahan (com doubles), conforme descrito por 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;
}

Eu recebo exatamente 2.0:

Kahan sum  = 2.000000000000000    (difference = 0.000000000000000)

E mesmo se eu mudar o doubles para floats no código acima, obtenho:

Kahan sum  = 2.000000000000000    (difference = 0.000000000000000)

Parece que Kahan é o caminho a percorrer!

Existe uma classe de algoritmos que resolve esse problema exato, sem a necessidade de classificar ou reordenar os dados .

Em outras palavras, o somatório pode ser feito em uma passagem pelos dados. Isso também torna esses algoritmos aplicáveis ​​em situações em que o conjunto de dados não é conhecido com antecedência, por exemplo, se os dados chegam em tempo real e a soma corrente precisa ser mantida.

Aqui está o resumo de um artigo recente:

Apresentamos um novo algoritmo online para a soma exata de um fluxo de números de ponto flutuante. Por "online" queremos dizer que o algoritmo precisa ver apenas uma entrada de cada vez e pode receber um fluxo de entrada de comprimento arbitrário de tais entradas enquanto requer apenas memória constante. Por “exato” queremos dizer que a soma da matriz interna de nosso algoritmo é exatamente igual à soma de todas as entradas, e o resultado retornado é a soma arredondada corretamente. A prova de exatidão é válida para todas as entradas (incluindo números não normalizados, mas estouro intermediário do módulo) e é independente do número de somas ou do número de condição da soma. O algoritmo precisa assintoticamente de apenas 5 FLOPs por soma e, devido ao paralelismo em nível de instrução, é executado apenas cerca de 2-3 vezes mais lento do que o óbvio, loop de "soma recursiva ordinária" rápido, mas burro, quando o número de somas é maior que 10.000. Portanto, até onde sabemos, é o mais rápido, mais preciso e mais eficiente em termos de memória entre os algoritmos conhecidos. Na verdade, é difícil ver como um algoritmo mais rápido ou requerendo significativamente menos FLOPs poderia existir sem melhorias de hardware. Um aplicativo para um grande número de summands é fornecido.

Fonte: Algoritmo 908: Soma Exata Online de Fluxos de Ponto Flutuante .

Com base na resposta de Steve de primeiro classificar os números em ordem crescente, eu apresentaria mais duas ideias:

1. Decida a diferença no expoente de dois números acima da qual você pode decidir que perderá muita precisão.

2. Em seguida, some os números em ordem até que o expoente do acumulador seja muito grande para o próximo número, em seguida, coloque o acumulador em uma fila temporária e inicie o acumulador com o próximo número. Continue até esgotar a lista original.

Você repete o processo com a fila temporária (tendo-a classificado) e com uma diferença possivelmente maior no expoente.

Acho que isso será bem lento se você tiver que calcular expoentes o tempo todo.

Tive uma experiência rápida com um programa e o resultado foi 1.99903

Acho que você pode fazer melhor do que ordenar os números antes de acumulá-los, porque durante o processo de acumulação, o acumulador fica cada vez maior. Se você tiver uma grande quantidade de números semelhantes, começará a perder a precisão rapidamente. Aqui está o que eu sugeriria:

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

É claro que esse algoritmo será mais eficiente com uma fila de prioridade em vez de uma 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);
    }
}

motorista:

#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();
}

Os números na fila são negativos porque topproduz o maior número, mas queremos o menor . Eu poderia ter fornecido mais argumentos de modelo para a fila, mas essa abordagem parece mais simples.

Isso não responde exatamente à sua pergunta, mas uma coisa inteligente a fazer é calcular a soma duas vezes, uma com o modo de arredondamento "arredondar para cima" e outra com "arredondar para baixo". Compare as duas respostas e você sabe / como / imprecisos são seus resultados e, portanto, precisa usar uma estratégia de soma mais inteligente. Infelizmente, a maioria das linguagens não torna a alteração do modo de arredondamento de ponto flutuante tão fácil quanto deveria ser, porque as pessoas não sabem que ele é realmente útil nos cálculos diários.

Dê uma olhada na aritmética de intervalo, onde você faz todas as contas assim, mantendo os valores mais altos e mais baixos conforme você avança. Isso leva a alguns resultados e otimizações interessantes.

A classificação mais simples que melhora a precisão é classificar pelo valor absoluto crescente. Isso permite que os menores valores de magnitude tenham a chance de se acumular ou cancelar antes de interagir com os valores de magnitude maiores que poderiam causar uma perda de precisão.

Dito isso, você pode fazer melhor rastreando várias somas parciais não sobrepostas. Aqui está um artigo que descreve a técnica e apresenta uma prova de precisão: www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps

Esse algoritmo e outras abordagens para soma exata de ponto flutuante são implementados em Python simples em: http://code.activestate.com/recipes/393090/ Pelo menos dois deles podem ser convertidos trivialmente para C ++.

Para IEEE 754 de precisão simples ou dupla ou números de formato conhecido, outra alternativa é usar uma matriz de números (passada pelo chamador, ou em uma classe para C ++) indexada pelo expoente. Ao adicionar números ao array, apenas números com o mesmo expoente são adicionados (até que um slot vazio seja encontrado e o número armazenado). Quando uma soma é solicitada, a matriz é somada do menor ao maior para minimizar o truncamento. Exemplo de precisão única:

/* 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;
}

exemplo de precisão dupla:

/* 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;
}

Seus flutuadores devem ser adicionados com precisão dupla. Isso lhe dará mais precisão adicional do que qualquer outra técnica pode. Para um pouco mais de precisão e significativamente mais velocidade, você pode criar, digamos, quatro somas e adicioná-las ao final.

Se você estiver adicionando números de precisão dupla, use long double para a soma - no entanto, isso só terá um efeito positivo em implementações onde long double realmente tem mais precisão do que double (normalmente x86, PowerPC dependendo das configurações do compilador).

Quanto à classificação, parece-me que, se você espera o cancelamento, os números devem ser somados em ordem decrescente de magnitude, não crescente. Por exemplo:

((-1 + 1) + 1e-20) resultará em 1e-20

mas

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

Na primeira equação, dois números grandes são cancelados, enquanto na segunda o termo 1e-20 se perde quando adicionado a 1, pois não há precisão suficiente para retê-lo.

Além disso, a soma de pares é bastante decente para somar muitos números.