Cancelamento catastrófico na soma do log

18

Estou tentando implementar a seguinte função no ponto flutuante de precisão dupla com baixo erro relativo :

logsum(x,y)=log(exp(x)+exp(y))

Isso é usado extensivamente em aplicativos estatísticos para adicionar probabilidades ou densidades de probabilidade representadas no espaço de log. Obviamente, exp(x) ou exp(y) poderiam facilmente estourar ou estourar, o que seria ruim porque o espaço de log é usado para evitar o estouramento em primeiro lugar. Esta é a solução típica:

logsum(x,y)=x+log1p(exp(yx))

O cancelamento de yx acontece, mas é mitigado pela exp . Pior ainda, de longe, é quando x e log1p(exp(yx)) encontram-se perto. Aqui está um gráfico de erro relativo:

insira a descrição da imagem aqui

A trama é cortada em para enfatizar a forma da curva de l o g s u m ( x , y ) = 0 , sobre a qual o cancelamento ocorrer. Eu vi erro até 10 - 11 e suspeito que ele fica muito pior. (FWIW, a função "verdade da terra" é implementada usando flutuadores de precisão arbitrária do MPFR com precisão de 128 bits.)1014logsum(x,y)=01011

Eu tentei outras reformulações, todas com o mesmo resultado. Com como a expressão exterior, o mesmo ocorre por erro tendo um registo de algo perto 1. Com l o g um p como a expressão exterior, cancelamento acontece na expressão interior.loglog1p

Agora, o absoluto de erro é muito pequena, de modo tem muito pequeno erro relativo (dentro de um epsilon). Alguém poderia argumentar que, porque um usuário de l o g s u m está realmente interessado em probabilidades (não probabilidades de log), este erro relativo terrível não é um problema. É provável que geralmente não seja, mas estou escrevendo uma função de biblioteca e gostaria que seus clientes pudessem contar com um erro relativo não muito pior do que um erro de arredondamento.exp(logsum(x,y))logsum

Parece que preciso de uma nova abordagem. O que poderia ser?

Neil Toronto
fonte
Eu não entendo o seu último parágrafo. "dentro de um epsilon" não significa nada para mim. Você quer dizer uma unidade em último lugar ? Quanto aos usuários interessados ​​em probabilidades, um pequeno erro de probabilidade de log resultará em um grande erro de probabilidade; portanto, esse não é o caso.
Aron Ahmadia 22/10/12
Por curiosidade, você já tentou tirar o melhor de seus dois métodos e traçar o erro disso? Então, tudo o que você precisa é a lógica certa para detectar em qual caso você está (espero que seja menos oneroso ou parte do custo exigido do algoritmo) e depois mude para o método apropriado.
Aron Ahmadia 22/10/12
axexp(x)exp(a)1xexp(x)
aa>1

Respostas:

12

logsum(x,y)=max(x,y)+log1p(exp(abs(xy))
logiexi=ξ+logiexiξ,   ξ=maxixi

Caso a soma do log esteja muito próxima de zero e você queira alta precisão relativa, provavelmente você pode usar usando um método preciso (isto é, mais do que precisão dupla) implementação de que é quase linear para pequeno .l e x p ( z ) : = log ( 1 + e - | z | ) z

logsum(x,y)=max(x,y)+lexp(xy)
lexp(z):=log(1+e|z|)
z
Arnold Neumaier
fonte
Em termos de erro absoluto, é. Em termos de erro relativo, é horrível quando a saída está próxima de zero.
Neil Toronto
@ NeilToronto: Por favor, dê um exemplo com duas entradas explícitas e , para que eu possa brincar com ela. yxy
Arnold Neumaier 22/10/12
Para x = -0,775 ey = -0,6175, recebo o erro 62271 ulps e o erro relativo 1.007e-11.
Neil Toronto
1
Calcule pontos de dados altamente precisos no intervalo de interesse - são necessários pelo menos dois intervalos diferentes devido ao comportamento assintótico. Pode-se usar a expressão definidora para z não próxima de zero. Para a faixa excepcional, ajuste uma função racional de grau suficientemente alto para obter a precisão desejada. Para estabilidade numérica, use polinômios de Bernstein ou polinômios de Tchebychev no numerador e denominador, adaptados ao intervalo de interesse. No final, expanda para uma fração contínua e descubra quanto é possível truncar os coeficientes sem medir a precisão.
Arnold Neumaier
1
Isso fornece Para obter faça o mesmo, mas aplicado à função lexp (z) -l (z). ml=l(z)m
Arnold Neumaier