Efficiently summing log quantities

Working in C++, I'd like to find the sum of some quantities, and then take the log of the sum:

log(a_1 + a_2 + a_3 + ... + a_n)

However, I do not have the quantities themselves, I only have their log'd values:

l_1 = log(a_1), l_2 = log(a_2), ... , l_n = log(a_n)

Is there any efficient way to get the log the sum of开发者_如何学编程 the a_i's? I'd like to avoid

log(s) = log(exp(l_1) + exp(l_2) + ... + exp(l_n))

if possible - exp becomes a bottleneck as the calculation is done many times.

How large is n?

This quantity is known as log-sum-exp and Lieven Vandenberghe talks about it on page 72 of his book. He also has an optimization package which uses this operation, and from a brief look, it seems as if he doesn't do anything special there, just exponentiates and adds. Perhaps exponentiation is not a serious bottleneck when n is small enough for the vector to fit into memory.

This operation often comes up in modeling and the bottleneck there is the sheer number of terms. Magnitude of n=2^100 is common, where the terms are implicitly represented. In such cases, there are various tricks to approximating this quantity relying on convexity of log-sum-exp. The simplest trick -- approximate log(s) as max(l1,l2,....,ln)

I don't know any way, because, in general, there is no way to calculate

e^x + e^y

using addition and only one exponentiation, which is equivalent to what you're asking.

As was mentioned in Frédéric Hamidi's comment above, even if you do sum the exponents, you have another problem to worry about: overflow. The link he gave gives a pretty good solution (following Fortran code copied from that link):

function log_sum_exp(v) result(e)
  real, dimension(:), intent(in) :: v   ! Input vector
  real                           :: e   ! Result is log(sum(exp(v)))
  real                           :: c   ! Shift constant

  ! Choose c to be the element of v that is largest in absolute value.
  if ( maxval(abs(v)) > maxval(v) ) then
     c = minval(v)
  else
     c = maxval(v)
  end if

  e = log(sum(exp(v-c))) + c
end function log_sum_exp

You could use the following identity:

log( a + b ) = log(a) + log( 1 + (b/a) )

It's not very elegant, but you could try the following.

• Obtain lg a_i from log a_i (divide by log 2).
• Let lg a_i = k + q where k is integer and q is real and 0 >= q >= 1
• Get a_i with 2^kpow(2,q) (use bit shifting for 2^k = 1 << k).
• You can speed up pow(2,q) with precomputed tables of limited precision for [0,1]

So the whole idea is to leverage a fast power-of-2 function. Hope it helps!

If s_k := sum(a_1 + ... + a_k) then s_{k+1} == s_k + f(l_{k+1} - s_k), where

f(x) := log(1+exp(x))

This function f can probably be computed with a Taylor series or similar with a speed comparable to exp, and possibly even inlined.

That only saves about two mathematical functions, but it might be a useful starting point.