Logsumexp implementation in C?

Question

Does anyone know the open source C library that provides the logsumexp function?

The logsumexp(a) function computes the sum of the exponentials of log (e ^ {a_1} + ... e ^ {a_n}) components of 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 allows you to effectively divide all arguments into the largest ones, and then add it back to the end to avoid overflow, so it behaved well to add a large number of equally scaled values ​​with errors creeping in if some arguments are many orders of magnitude larger than others.

If you want it to run smoothly when you set 0 arguments, you need to add a case to it :)