Comparação relativa de números de ponto flutuante

10

Eu tenho uma função numérica f(x, y)retornando um número de ponto flutuante duplo que implementa alguma fórmula e quero verificar se ele está correto em relação às expressões analíticas para todas as combinações de parâmetros xe no yqual estou interessado. Qual é a maneira correta de comparar os valores calculados e números de ponto flutuante analítico?

Digamos que os dois números são ae b. Até agora, tenho me assegurado de que os erros absolutos ( abs(a-b) < eps) e relativos ( abs(a-b)/max(abs(a), abs(b)) < eps) sejam menores que eps. Dessa forma, ele detectará imprecisões numéricas, mesmo que os números sejam, digamos, em torno de 1e-20.

No entanto, hoje eu descobri um problema, o valor numérico ae o valor analítico beram:

In [47]: a                                                                     
Out[47]: 5.9781943146790832e-322

In [48]: b                                                                     
Out[48]: 6.0276008792632078e-322

In [50]: abs(a-b)                                                              
Out[50]: 4.9406564584124654e-324

In [52]: abs(a-b) / max(a, b)                                                  
Out[52]: 0.0081967213114754103

Portanto, o erro absoluto [50] é (obviamente) pequeno, mas o erro relativo [52] é grande. Então, eu pensei que eu tinha um bug no meu programa. Ao depurar, percebi, que esses números são anormais . Como tal, escrevi a seguinte rotina para fazer a comparação relativa adequada:

real(dp) elemental function rel_error(a, b) result(r)
real(dp), intent(in) :: a, b
real(dp) :: m, d
d = abs(a-b)
m = max(abs(a), abs(b))
if (d < tiny(1._dp)) then
    r = 0
else
    r = d / m
end if
end function

Where tiny(1._dp)retorna 2.22507385850720138E-308 no meu computador. Agora tudo funciona e eu simplesmente recebo 0 como o erro relativo e está tudo ok. Em particular, o erro relativo acima [52] está errado, é simplesmente causado pela precisão insuficiente dos números desnormais. Minha implementação da rel_errorfunção está correta? Devo apenas verificar se abs(a-b)é menor que minúsculo (= denormal) e retornar 0? Ou devo verificar alguma outra combinação, como max(abs(a), abs(b))?

Gostaria apenas de saber qual é o caminho "adequado".

floating-point Ondřej Čertík
fonte

Respostas:

11

Você pode verificar diretamente se há denormals usando a isnormal()partir de math.h(C99 ou posterior, POSIX.1 ou posterior). No Fortran, se o módulo ieee_arithmeticestiver disponível, você poderá usá-lo ieee_is_normal(). Para ser mais preciso sobre a igualdade difusa, você deve considerar a representação em ponto flutuante dos denormais e decidir o que você quer dizer com resultados suficientemente bons.

Mais precisamente, para acreditar que um dos resultados está correto, você deve ter certeza de que não perdeu muitos dígitos em uma etapa intermediária. Computar com denormals geralmente não é confiável e deve ser evitado com a redimensionamento interna do algoritmo. Para garantir que seu dimensionamento interno seja bem-sucedido, recomendo ativar as exceções de ponto flutuante usando o feenableexcept()C99 ou o ieee_arithmeticmódulo no Fortran.

Embora você possa fazer com que seu aplicativo capte o sinal gerado em exceções de ponto flutuante, todos os kernels que eu tentei redefinir o sinalizador de hardware para fetestexcept()que não retorne um resultado útil. Quando executados -fp_trap, os programas PETSc (por padrão) imprimem um rastreamento de pilha quando um erro de ponto flutuante é gerado, mas não identificam a linha incorreta. Se você executar em um depurador, o depurador preservará o sinalizador de hardware e interromperá a expressão ofensiva. Você pode verificar o motivo exato chamando fetestexceptdo depurador onde o resultado é um OR bit a bit dos seguintes sinalizadores (os valores podem variar de acordo com a máquina, consulte fenv.h; esses valores são para x86-64 com glibc).

  • FE_INVALID = 0x1
  • FE_DIVBYZERO = 0x4
  • FE_OVERFLOW = 0x8
  • FE_UNDERFLOW = 0x10
  • FE_INEXACT = 0x20
Jed Brown
fonte
Obrigado pela excelente resposta. A expressão analítica que eu comparo no regime assintótico é exp(log_gamma(m+0.5_dp) - (m+0.5_dp)*log(t)) / 2para m = 234, t = 2000. Ele chega a zero rapidamente à medida que aumenta m. Tudo o que quero ter certeza de que minha rotina numérica retorne números "corretos" (retornar zero também é perfeitamente adequado) para pelo menos 12 dígitos significativos. Portanto, se o cálculo retornar um número não normal, então é simplesmente zero e não deverá haver problema. Portanto, apenas a rotina de comparação precisa ser robusta contra isso.
Ondřej Čertík
5

Donald Knuth tem uma proposta para um algoritmo de comparação de ponto flutuante no volume 2 "Algoritmos seminuméricos" de "A arte da programação de computadores". Foi implementado em C por Th. Belding (consulte o pacote fcmp ) e está disponível no GSL .

GertVdE
fonte
2
Aqui está minha implementação do Fortran: gist.github.com/3776847 , observe que eu preciso lidar explicitamente com números desnormais de qualquer maneira. Caso contrário, eu acho que é praticamente equivalente ao erro relativo, a única diferença é que, em vez de fazer abs(a-b)/max(a, b) < eps, nós fazemos abs(a-b)/2**exponent(max(a, b)) < eps, que praticamente deixa cair a mantissa na max(a, b), então, na minha opinião, a diferença é insignificante.
Ondřej Čertík
5

Números desnormalizados idealmente arredondados podem realmente ter um erro relativo alto. (Liberar isso para zero enquanto ainda o chama de erro relativo é enganoso.)

Mas quase zero, calcular erros relativos não faz sentido.

Portanto, mesmo antes de atingir números não normalizados, você provavelmente deve mudar para a precisão absoluta (a saber, o que você deseja garantir neste caso).

yx|y−x|≤absacc+relacc∗max(|x|,|y|)

Em seguida, os usuários do seu código sabem exatamente quanta precisão eles realmente têm.

Arnold Neumaier
fonte
Você tem certeza de que não faz sentido calcular erros relativos próximos de zero? Eu acho que só faz sentido se houver perda de precisão (por qualquer motivo). Se, por exemplo, houver perda de precisão para x <1e-150 devido a alguns problemas numéricos (como subtrair dois números grandes), você está certo. No meu caso, no entanto, os números parecem ser precisos até zero, exceto quando atinge os números desnormais. Portanto, no meu caso, absacc = 1e-320 ou mais, e posso verificar abs(a-b) < tiny(1._dp)como faço acima.
Ondřej Čertík
@ OndřejČertík: Nesse caso, substitua o 1e-150 por 1e-300 ou qualquer outro limite que você possa verificar. De qualquer forma, muito próximo de zero, você comete um erro absoluto e sua reivindicação de erro deve refletir isso em vez de declarar que o erro relativo é zero.
Arnold Neumaier
Eu vejo. Posso verificar que tudo funciona para números maiores que tiny(1._dp)=2.22507385850720138E-308(cometi um erro no meu comentário anterior, é 2e-308, não 1e-320). Então este é o meu erro absoluto. Então eu preciso comparar o erro relativo. Entendo o seu ponto, acho que você está certo. Obrigado!
Ondřej Čertík
11
@ OndřejČertík: Para encontrar o erro relativo adicional dado absacc, monitore o máximo de . |y−x|−absaccmax(|x|,|y|)
Arnold Neumaier 25/09/12