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.