Analisando erro numérico na função C ++

Suponha que eu tenha uma função que tome como entrada vários valores de ponto flutuante (simples ou duplo), faça alguma computação e produza valores de ponto flutuante de saída (também únicos ou duplos). Estou trabalhando principalmente com o MSVC 2008, mas também planejo trabalhar com o MinGW / GCC. Estou programando em C ++.

Qual é a maneira típica de medir programaticamente quanto erro há nos resultados? Supondo que eu precise usar uma biblioteca de precisão arbitrária: qual é a melhor biblioteca se não me importo com velocidade?

Se você está procurando um bom limite para o seu erro de arredondamento, não precisa necessariamente de uma biblioteca de precisão arbitrária. Você pode usar a análise de erro em execução.

Não consegui encontrar uma boa referência on-line, mas está tudo descrito na Seção 3.3 do livro de Nick Higham "Exatidão e estabilidade de algoritmos numéricos". A ideia é bastante simples:

1. Re-fatore seu código para que você tenha uma única atribuição de uma única operação aritmética em cada linha.
2. Para cada variável, por exemplo x, crie uma variável x_errque seja inicializada com zero quando xfor atribuída uma constante.
3. Para cada operação, por exemplo z = x * y, atualize a variável z_errusando o modelo padrão da aritmética de ponto flutuante e os zerros resultantes e em execução x_erre y_err.
4. O valor de retorno da sua função também deve ter um _errvalor respectivo anexado. Este é um limite dependente de dados no seu erro total de arredondamento.

A parte complicada é a etapa 3. Para as operações aritméticas mais simples, você pode usar as seguintes regras:

• z = x + y -> z_err = u*abs(z) + x_err + y_err
• z = x - y -> z_err = u*abs(z) + x_err + y_err
• z = x * y -> z_err = u*abs(z) + x_err*abs(y) + y_err*abs(x)
• z = x / y -> z_err = u*abs(z) + (x_err*abs(y) + y_err*abs(x))/y^2
• z = sqrt(x) -> z_err = u*abs(z) + x_err/(2*abs(z))

onde u = eps/2é o arredondamento da unidade. Sim, as regras +e -são as mesmas. As regras para qualquer outra operação op(x)podem ser facilmente extraídas usando a expansão da série Taylor do resultado aplicado op(x + x_err). Ou você pode tentar pesquisar no Google. Ou usando o livro de Nick Higham.

Como exemplo, considere o seguinte código Matlab / Octave que avalia polinômios nos coeficientes aem um ponto xusando o esquema Horner:

function s = horner ( a , x )
    s = a(end);
    for k=length(a)-1:-1:1
        s = a(k) + x*s;
    end

Para a primeira etapa, dividimos as duas operações em s = a(k) + x*s:

function s = horner ( a , x )
    s = a(end);
    for k=length(a)-1:-1:1
        z = x*s;
        s = a(k) + z;
    end

Em seguida, apresentamos as _errvariáveis. Observe que as entradas ae xsão consideradas exatas, mas também podemos exigir que o usuário transmita valores correspondentes para a_erre x_err:

function [ s , s_err ] = horner ( a , x )
    s = a(end);
    s_err = 0;
    for k=length(a)-1:-1:1
        z = x*s;
        z_err = ...;
        s = a(k) + z;
        s_err = ...;
    end

Por fim, aplicamos as regras descritas acima para obter os termos do erro:

function [ s , s_err ] = horner ( a , x )
    u = eps/2;
    s = a(end);
    s_err = 0;
    for k=length(a)-1:-1:1
        z = x*s;
        z_err = u*abs(z) + s_err*abs(x);
        s = a(k) + z;
        s_err = u*abs(s) + z_err;
    end

Observe que, como não temos a_errou x_err, por exemplo, eles são assumidos como zero, os respectivos termos são simplesmente ignorados nas expressões de erro.

Et voilà! Agora temos um esquema Horner que retorna uma estimativa de erro dependente de dados (nota: esse é um limite superior do erro) ao lado do resultado.

Como observação lateral, como você está usando C ++, considere criar sua própria classe para valores de ponto flutuante, que carrega o _errtermo e sobrecarrega todas as operações aritméticas para atualizar esses valores, conforme descrito acima. Para códigos grandes, essa pode ser a rota mais fácil, embora computacionalmente menos eficiente. Dito isto, você poderá encontrar essa classe online. Uma rápida pesquisa no Google me deu esse link .

PS Observe que tudo isso funciona apenas em máquinas que cumprem estritamente a IEEE-754, ou seja, todas as operações aritméticas são precisas para . Essa análise também fornece um limite mais preciso e realista do que a aritmética do intervalo, pois, por definição, você não pode representar um número em ponto flutuante, ou seja, seu intervalo seria arredondado para o próprio número.$\pm u$$x(1 \pm u)$

Uma boa biblioteca portátil e de código aberto para aritmética de ponto flutuante de precisão arbitrária (e muito mais) é o NTL de Victor Shoup , disponível na forma de código-fonte C ++.

Em um nível inferior, está a Biblioteca Bignum da GNU Multiple Precision (GMP) , também um pacote de código aberto.

O NTL pode ser usado com o GMP, se for necessário um desempenho mais rápido, mas o NTL fornece suas próprias rotinas básicas que certamente são utilizáveis ​​se você "não se importar com a velocidade". O GMP afirma ser a "biblioteca de bignum mais rápida". O GMP é amplamente escrito em C, mas possui uma interface C ++.

Adicionado: embora a aritmética do intervalo possa fornecer limites superiores e inferiores da resposta exata de maneira automatizada, isso não mede com precisão o erro em um cálculo de precisão "padrão" porque o tamanho do intervalo geralmente cresce a cada operação (seja em senso absoluto de erro).

A maneira típica de encontrar o tamanho do erro, seja para erros de arredondamento ou de discretização, etc., é calcular um valor de precisão extra e comparar com o valor de precisão "padrão". Somente uma precisão extra modesta é necessária para determinar o tamanho do erro com precisão razoável, uma vez que os erros de arredondamento são substancialmente maiores na precisão "padrão" do que no cálculo da precisão extra.

O ponto pode ser ilustrado através da comparação de cálculos de precisão simples e dupla. Observe que, em C ++, as expressões intermediárias são sempre computadas em (pelo menos) precisão dupla, portanto, se quisermos ilustrar como seria um cálculo em precisão única "pura", precisamos armazenar os valores intermediários em precisão única.

Snippet de código C

    float fa,fb;
    double da,db,err;
    fa = 4.0;
    fb = 3.0;
    fa = fa/fb;
    fa -= 1.0;

    da = 4.0;
    db = 3.0;
    da = da/db;
    da -= 1.0;

    err = fa - da;
    printf("Single precision error wrt double precision value\n");
    printf("Error in getting 1/3rd is %e\n",err);
    return 0;

A saída de cima (cadeia de ferramentas Cygwin / MinGW32 GCC):

Single precision error wrt double precision value
Error in getting 1/3rd is 3.973643e-08

Portanto, o erro é sobre o que se espera no arredondamento de 1/3 para a precisão única. (Eu suspeitaria) que não se importaria em obter mais de duas casas decimais no erro correto, pois a medição do erro é por magnitude e não por exatidão.

GMP (ou seja, a biblioteca GNU Multiple Precision) é a melhor biblioteca de precisão arbitrária que eu conheço.

Não conheço nenhuma maneira programática de medir o erro nos resultados de uma função arbitrária de ponto flutuante. Uma coisa que você pode tentar é calcular uma extensão de intervalo de uma função usando aritmética de intervalo . No C ++, você precisaria usar algum tipo de biblioteca para calcular as extensões de intervalo; uma dessas bibliotecas é a Biblioteca Aritmética Boost Interval. Basicamente, para medir o erro, você forneceria como argumentos para os intervalos de suas funções uma largura de 2 vezes o arredondamento da unidade (aproximadamente), centralizada nos valores de interesse e, em seguida, sua saída seria uma coleção de intervalos, a largura de o que forneceria uma estimativa conservadora do erro. Uma dificuldade dessa abordagem é que a aritmética de intervalo usada dessa maneira pode superestimar o erro em quantidades significativas, mas essa abordagem é a mais "programática" em que consigo pensar.

A estimativa rigorosa e automática de erros pode ser obtida por análise de intervalo . Você trabalha com intervalos em vez de números. Por exemplo, adição:

[a,b] + [c,d] = [min(a+c, a+d, b+c, b+d), max (a+c, a+d, b+c, b+d)] = [a+c, b+d]

O arredondamento também pode ser tratado com rigor, consulte Aritmética de intervalos arredondados .

Contanto que sua entrada consista em intervalos estreitos, as estimativas são aceitáveis ​​e são baratas de calcular. Infelizmente, o erro geralmente é superestimado; veja o problema da dependência .

Não conheço nenhuma biblioteca aritmética de intervalo de precisão arbitrária.

Depende do seu problema em mãos se a aritmética de intervalo pode atender às suas necessidades ou não.

A biblioteca GNU MPFR é uma biblioteca flutuante de precisão arbitrária que possui alta precisão (em particular, arredondamento correto para todas as operações, o que não é tão fácil quanto parece) como um de seus principais pontos de foco. Ele usa o GNU MP sob o capô. Ele possui uma extensão chamada MPFI que faz a aritmética de intervalos, que - como sugere a resposta de Geoff - pode ser útil para fins de verificação: continue aumentando a precisão de trabalho até que o intervalo resultante caia dentro de limites pequenos.

Isso nem sempre funciona, no entanto; em particular, não é necessariamente eficaz se você estiver fazendo algo como integração numérica, em que cada etapa carrega um "erro" independente de problemas de arredondamento. Nesse caso, tente um pacote especializado como o COZY infinity, que faz isso muito bem usando algoritmos específicos para limitar o erro de integração (e usando os chamados modelos de Taylor em vez de intervalos).

Disseram-me que o MPIR é uma boa biblioteca para usar se você estiver trabalhando com o Visual Studio:

http://mpir.org/