1

I would like to compute log( exp(A1) + exp(A2) ).
The formula below 

 log(exp(A1) + exp(A2) ) = log[exp(A1)(1 + exp(A2)/exp(A1))] = A1 + log(1+exp(A2-A1))

is useful when A1 and A2 are large and numerically exp(A1)=Inf (or exp(A2)=Inf). (this formula is discussed in this thread -> How to calculate log(sum of terms) from its component log-terms). The formula is true when the role of A1 and A2 are replaced.

My concern of this formula is when A1 and A2 are very small. For example, when A1 and A2 are:

 A1 <- -40000
 A2 <- -45000

then the direct calculation of log(exp(A1) + exp(A2) ) is:

 log(exp(A1) + exp(A2))
 [1] -Inf

Using the formula above gives:

 A1 + log(1 + exp(A2-A1))
 [1] -40000

which is the value of A1. Ising the formula above with flipped role of A1 and A2 gives:

A2 + log(1 + exp(A1-A2))
[1] Inf

Which of the three values are the closest to the true value of log(exp(A1) + exp(A2))? Is there robust way to compute log(exp(A1) + exp(A2)) that can be used both when A1, A2 are small and A1, A2 are large.

Thank you in advance

Community
  • 1
  • 1
FairyOnIce
  • 2,526
  • 7
  • 27
  • 48

2 Answers2

3

You should use something with more accuracy to do the direct calculation.

It’s not “useful when [they’re] large”. It’s useful when the difference is very negative.

When x is near 0, then log(1+x) is approximately x. So if A1>A2, we can take your first formula:

log(exp(A1) + exp(A2)) = A1 + log(1+exp(A2-A1))

and approximate it by A1 + exp(A2-A1) (and the approximation will get better as A2-A1 is more negative). Since A2-A1=-5000, this is more than negative enough to make the approximation sufficient.

Regardless, if y is too far from zero (either way) exp(y) will (over|under)flow a double and result in 0 or infinity (this is a double, right? what language are you using?). This explains your answers. But since exp(A2-A1)=exp(-5000) is close to zero, your answer is approximately -40000+exp(-5000), which is indistinguishable from -40000, so that one is correct.

Teepeemm
  • 4,331
  • 5
  • 35
  • 58
1

in such huge exponent differences the safest you can do without arbitrary precision is

  • chose the biggest exponent let it be Am = max(A1,A2)
  • so: log(exp(A1)+exp(A2)) -> log(exp(Am)) = Am
  • that is the closest you can get for such case
  • so in your example the result is -40000+delta
  • where delta is something very small

If you want to use the second formula then all breaks down to computing log(1+exp(A))

  • if A is positive then the result is far from the real thing
  • if A is negative then it will truncate to log(1)=0 so you get the same result as in above

[Notes]

  • your exponent difference is base^500
  • single precision 32bit float can store numbers up to (+/-)2^(+/-128)
  • double precision 64bit float can store numbers up to (+/-)2^(+/-1024)
  • so when your base is 10 or e then this is nowhere near enough what you need
  • if you have quadruple precision that should be enough but when you start changing the exp difference again yo will quickly get to the same point as now

[PS] if you need more precision without arbitrary precision

  • you can try to create own number class
  • with internal store of numbers like number=a^b
  • where a,b are floats
  • but for that you would need to code all basic functions
  • *,/ is easy
  • +,- is a nightmare but there could be some approaches/algorithms out there even for this
Spektre
  • 49,595
  • 11
  • 110
  • 380