Go implementation of logaddexp

[Fumin Wang]

Does anyone know of a robust and correct Go implementation of logaddexp http://docs.scipy.org/doc/numpy/reference/generated/numpy.logaddexp.html ?

The mathematical definition of logaddexp is

logaddexp(x1, x2) = log(exp(x1) + exp(x2))

One might be tempted to simply implement it as math.Log(math.Exp(x1) + math.Exp(x2)), but then this easily runs into underflow issues on real world statistical/machine learning applications/problems.

[Roberto Zanotto]

Use this:

log(exp(a) + exp(b)) = a + log(1 + exp(b-a)) = a + math.Log1p(exp(b-a))

Source: http://lingpipe-blog.com/2009/06/25/log-sum-of-exponentials/ (see also the comment section)

[Tamás Gulácsi]

How shall we implement it?
What's the matter with the naïve method of log(exp(a)+exp(b)) ?

[Roberto Zanotto]

Doing more research it looks like my answer was a bit rushed. Log1p helps a bit but we are calculating log(1 + exp(x)), so doing exp and then log still has precision problems. Here it is explained how to compute log(1 + exp(x)) more accurately: http://sachinashanbhag.blogspot.it/2014/05/numerically-approximation-of-log-1-expy.html

The end result for logaddexp looks like this: 

func logaddexp(a, b float64) float64 {

return a + log1exp(b-a)

}

func log1exp(x float64) float64 {

if x > 35 {

return x

}

if x < -10 {

return math.Exp(x)

}

return math.Log1p(math.Exp(x))

}

Here you have it on the playground with some preliminary testing: https://play.golang.org/p/46FU2BWLP-

[Fumin Wang]

Hi Robert

For the implementation of log1exp, for my case of negative inputs, my observation is that the current implementation of math.Log1p is pretty close to math.Exp. Moreover, after consulting the implementation of numpy, I ended up the following implementation:

// logaddexp performs log(exp(x) + exp(y))

func logaddexp(x, y float64) float64 {

  tmp := x - y

  if tmp > 0 {

    return x + math.Log1p(math.Exp(-tmp))

  } else if tmp <= 0 {

    return y + math.Log1p(math.Exp(tmp))

  } else {

    // Nans, or infinities of the same sign involved

    log.Printf("logaddexp %f %f", x, y)

    return x + y

  }

}

Thanks for your answers and inspirations, though.