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

105

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.

Yippie-Ki-Yay
fonte
17
Na verdade, isso tem tudo a ver com a programação do mundo real. No entanto, muitos aplicativos realmente não se importam com a melhor precisão absoluta do cálculo, desde que seja 'muito próximo'. Aplicações de engenharia? Extremamente importante. Aplicações médicas? Extremamente importante. Estatísticas em grande escala? Um pouco menos precisão é aceitável.
Zéychin
18
Não responda a menos que você realmente saiba e possa apontar para uma página que explique seu raciocínio em detalhes. Já existe tanta porcaria sobre números de ponto flutuante voando por aí que não queremos adicionar a isso. Se você acha que sabe. PARE. porque se você apenas pensa que sabe, então provavelmente está errado.
Martin York
4
@ Zéychin "Aplicações de engenharia? Extremamente importantes. Aplicações médicas? Extremamente importantes." ??? Acho que você ficaria surpreso se soubesse a verdade :)
BЈовић
3
@Zeychin O erro absoluto é irrelevante. O que é importante é o erro relativo. Se poucos centésimos de radiano são 0,001%, então quem se importa?
BЈовић
3
Eu realmente recomendo esta leitura: "o que todo cientista da computação precisa saber sobre ponto flutuante" perso.ens-lyon.fr/jean-michel.muller/goldberg.pdf
Mohammad Alaggan

Respostas:

108

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.

Steve Jessop
fonte
8
1 para explicação. Isso é um pouco contra-intuitivo, pois a adição é geralmente numericamente estável (ao contrário da subtração e divisão).
Konrad Rudolph
2
@Konrad, pode ser numericamente estável, mas não é preciso dadas as diferentes magnitudes de operandos :)
MSN
3
@ 6502: eles são classificados em ordem de magnitude, então -1 vem no final. Se o verdadeiro valor do total for de magnitude 1, tudo bem. Se você estiver somando três valores: 1 / bilhão, 1 e -1, então, você obterá 0, ponto no qual você tem que responder a pergunta prática interessante - você precisa de uma resposta que seja precisa na escala do soma verdadeira ou você só precisa de uma resposta precisa na escala dos maiores valores? Para algumas aplicações práticas, o último é bom o suficiente, mas quando não é, você precisa de uma abordagem mais sofisticada. A física quântica usa renormalização.
Steve Jessop
8
Se você vai ficar com este esquema simples, eu sempre adicionaria os dois números com a magnitude mais baixa e reinserirei a soma no conjunto. (Bem, provavelmente uma classificação por mesclagem funcionaria melhor aqui. Você poderia usar a parte da matriz que contém os números somados anteriormente como uma área de trabalho para as somas parciais.)
Neil
2
@Kevin Panko: A versão simples é que um float de precisão simples tem 24 dígitos binários, o maior dos quais é o maior bit definido no número. Então, se você somar dois números que diferem em magnitude em mais de 2 ^ 24, você sofre perda total do valor menor, e se eles diferem em magnitude em um grau menor, você perde um número correspondente de bits de precisão do menor número.
Steve Jessop
88

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
Daniel Pryden
fonte
3
1 adorável adição a este tópico. Qualquer compilador que "otimiza avidamente" essas instruções deve ser banido.
Chris A.
1
É um método simples para quase dobrar a precisão, usando duas variáveis ​​de soma sume cde magnitude diferente. Pode ser estendido trivialmente para N variáveis.
MSalters de
2
@ChrisA. bem, você pode controlar isso explicitamente em todos os compiladores que contam (por exemplo, -ffast-mathno GCC).
Konrad Rudolph
6
@Konrad Rudolph obrigado por apontar que esta é uma otimização possível com -ffast-math. O que aprendi com essa discussão e este link é que, se você se preocupa com a precisão numérica, provavelmente deve evitar o uso, -ffast-mathmas isso em muitas aplicações onde você pode estar limitado pela CPU, mas não se preocupa com cálculos numéricos precisos (programação de jogos, por exemplo ), -ffast-mathé razoável de usar. Assim, gostaria de emendar meu comentário "banido" com palavras fortes.
Chris A.
Usar variáveis ​​de precisão dupla para sum, c, t, yajudará. Você também precisa adicionar sum -= cantes return sum.
G. Cohen
34

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!

Andrew Stein
fonte
Meu valor "grande" é igual a 1, não 1e9. Sua segunda resposta, adicionada em ordem crescente de tamanho, é matematicamente correta (1 bilhão, mais um bilhão de bilionésimos, é 1 bilhão e 1), embora mais por sorte qualquer solidez geral do método :-) Observe que doublenão sofre mal perda de precisão na soma de um bilhão de bilionésimos, uma vez que possui 52 bits significativos, enquanto o IEEE floatpossui apenas 24 e teria.
Steve Jessop
@Steve, erro meu, desculpas. Eu atualizei o código de exemplo para o que você pretendia.
Andrew Stein
4
Kahan ainda tem precisão limitada, mas para construir um caso matador você precisa que a soma principal e o acumulador de erros ccontenham valores muito maiores do que a próxima soma. Isso significa que a soma é muito, muito menor do que a soma principal, então terá que haver uma quantidade enorme deles para somar muito. Especialmente com doublearitmética.
Steve Jessop
14

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 .

NPE
fonte
1
@Inverse: Ainda existem bibliotecas convencionais por aí. Como alternativa, comprar o PDF online custa US $ 5 a US $ 15 (dependendo se você é um membro ACM). Por último, a DeepDyve parece estar se oferecendo para emprestar o jornal por 24 horas por US $ 2,99 (se você é novo no DeepDyve, pode até conseguir obtê-lo gratuitamente como parte de seu teste gratuito): deepdyve.com/lp/acm /…
NPE
2

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

Quamrana
fonte
2

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.

fredoverflow
fonte
2

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.

rjmunro
fonte
0

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 ++.

Raymond Hettinger
fonte
0

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;
}
rcgldr
fonte
Isso soa um pouco como o método de Malcolm 1971 ou, mais ainda, sua variante que usa o expoente de Demmel e Hida ("Algoritmo 3"). Há outro algoritmo por aí que faz um loop baseado em carry como o seu, mas não consigo encontrar no momento.
ZachB de
@ZachB - o conceito é semelhante à ordenação de mesclagem ascendente para lista vinculada , que também usa um pequeno array, onde array [i] aponta para lista com 2 ^ i nós. Não sei até onde isso vai. No meu caso, foi a autodescoberta na década de 1970.
rcgldr
-1

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).

gnasher729
fonte
1
“Isso lhe dará mais precisão adicional do que qualquer outra técnica pode”. Você percebeu que sua resposta veio mais de um ano depois de uma resposta anterior anterior que descreveu como usar a soma exata?
Pascal Cuoq
O tipo "longo duplo" é horrível e você não deve usá-lo.
Jeff
-1

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.

KOAD
fonte