logsumexp implementation in C?

Question

Does anybody know of an open source numerical C library that provides the logsumexp-function?

The logsumexp(a) function computes the sum of exponentials log(e^{a_1}+...e^{a_n}) of the components of the array a, avoiding numerical overflow.

Answer 1

Here's a very simple implementation from scratch (tested, at least minimally):

double logsumexp(double nums[], size_t ct) {
  double max_exp = nums[0], sum = 0.0;
  size_t i;

  for (i = 1 ; i < ct ; i++)
    if (nums[i] > max_exp)
      max_exp = nums[i];

  for (i = 0; i < ct ; i++)
    sum += exp(nums[i] - max_exp);

  return log(sum) + max_exp;
}

This does the trick of effectively dividing all of the arguments by the largest, then adding its log back in at the end to avoid overflow, so it's well-behaved for adding a large number of similarly-scaled values, with errors creeping in if some arguments are many orders of magnitude larger than others.

If you want it to run without crashing when given 0 arguments, you'll have to add a case for that :)