Quando log1p e expm1 devem ser usados?

30

Eu tenho uma pergunta simples que é realmente difícil para o Google (além do artigo canônico O que todo cientista da computação deve saber sobre aritmética de ponto flutuante ).

Quando funções como log1pou expm1devem ser usadas em vez de loge exp? Quando eles não devem ser usados? Como as diferentes implementações dessas funções diferem em termos de uso?

Tim
fonte
2
Bem-vindo ao Scicomp.SE! Essa é uma pergunta muito razoável, mas seria mais fácil responder se você explicasse um pouco a que log1p você está se referindo (especialmente como ela é implementada, para que não tenhamos que adivinhar).
Christian Clason
4
Para argumentos com valor real, log1p e expm1 devem ser usados ​​quando for pequeno, por exemplo, quando na precisão do ponto flutuante. Consulte, por exemplo, docs.scipy.org/doc/numpy/reference/generated/numpy.expm1.html e docs.scipy.org/doc/numpy/reference/generated/numpy.log1p.html . (x)(x)x1 1+x=1 1
GoHokies
@ChristianClason obrigado, estou me referindo principalmente ao C ++ std ou R, mas, como você pergunta, começo a pensar que aprender sobre as diferenças nas implementações também seria muito interessante.
Tim
Aqui está um exemplo: scicomp.stackexchange.com/questions/8371/…
Juan M. Bello-Rivas
11
@ user2186862 "quando é pequeno" está correto, mas não apenas "quando 1 + x = 1 com precisão de ponto flutuante" (o que acontece para x 10 - 16 na aritmética usual de dupla precisão). As páginas de documentação que você vinculou mostram que elas já são úteis para x 10 - 10 , por exemplo. x1 1+x=1 1x10-16x10-10
Federico Poloni

Respostas:

25

Nós todos sabemos que

exp(x)=n=0 0xnn!=1 1+x+1 12x2+...
implica que para|x|1 1, temosexp(x)1 1+x. Isso significa que, se tivermos que avaliar no ponto flutuanteexp(x)-1 1, para|x|1 1cancelamento catastrófico pode ocorrer.

Isso pode ser facilmente demonstrado em python:

>>> from math import (exp, expm1)

>>> x = 1e-8
>>> exp(x) - 1
9.99999993922529e-09
>>> expm1(x)
1.0000000050000001e-08

>>> x = 1e-22
>>> exp(x) - 1
0.0
>>> expm1(x)
1e-22

Os valores exatos são

exp(10-8)-1 1=0.000000010000000050000000166666667083333334166666668...exp(10-22)-1 1=0.0000000000000000000001000000000000000000000000005000000...

Em geral, uma implementação "precisa" expe expm1deve estar correta para não mais que 1ULP (isto é, uma unidade do último local). No entanto, como atingir essa precisão resulta em código "lento", algumas vezes uma implementação rápida e menos precisa está disponível. Por exemplo, na CUDA, temos expfe expm1f, onde fsignifica rápido. De acordo com o guia de programação CUDA C, app. D o expftem um erro de 2ULP.

Se você não se importa com erros da ordem de poucos ULPS, geralmente implementações diferentes da função exponencial são equivalentes, mas cuidado para que os erros possam estar ocultos em algum lugar ... (Lembra do bug do Pentium FDIV ?)

Portanto, é bem claro que expm1deve ser usado para calcular exp(x)-1 1 para x pequeno . Usá-lo para o x geral não é prejudicial, pois expm1é esperado que seja preciso em toda a sua gama:

>>> exp(200)-1 == exp(200) == expm1(200)
True

(No exemplo acima 1 1 está bem abaixo de 1ULP de exp(200) , portanto, todas as três expressões retornam exatamente o mesmo número de ponto flutuante.)

Uma discussão semelhante vale para as funções inversas loge log1pdesde registro(1 1+x)x para |x|1 1 .

Stefano M
fonte
11
Essa resposta já estava contida nos comentários da pergunta do OP. No entanto, me senti útil em dar uma explicação mais longa (embora básica) apenas para maior clareza, na esperança de que seja útil para alguns leitores.
Stefano M
OK, mas, em seguida, pode-se simplesmente concluir "para que eu possa usar sempre expm1 vez de exp" ...
Tim
11
@tim sua conclusão está errada: você sempre pode usar em expm1(x)vez de exp(x)-1. Claro que exp(x) == exp(x) - 1não se aplica em geral.
Stefano M
OK, está claro. E existem critérios claros para ? x1 1
Tim
11
@ Tim não existe um limite definido, e a resposta depende da precisão da implementação do ponto flutuante (e do problema que está sendo resolvido). Embora expm1(x)deva ter precisão de 1ULP em todo o intervalo , perde progressivamente a precisão de poucos ULPs quando x 1 para uma quebra completa quando x < ϵ , onde ϵ é máquina-épsilon. 0 0x1 1exp(x) - 1x1 1x<ϵϵ
Stefano M
1

Para expandir a diferença entre loge log1ppode ser útil recuperar o gráfico se o logaritmo:

Logaritmo

logx0 0em(x)-x0 0em(x)em(1 1e)=-1 1em(1 1e10)=-10

x0 0em(x+1 1)0 0em(1 1+1 1e)0,31em(1 1+1 1e10)0,000045log1p

1 1log0 01 1log1p

sfmiller940
fonte