Logarithm of A Sum of Exponentials: Part 3

This post is a continuation of the Logarithm of A Sum of Exponentials series. Here are part 1 and part 2 of the series. I also talked about handling \(f(y) = \log (1 + \exp(y))\) in an exception-safe fashion, here.

Suppose we are given \(x_1\) and \(x_2\). Without any loss of generality let us suppose \(x_1 > x_2\).

We can write \(F = \log(e^{x_1} + e^{x_2})\) as \[F =  \log\left(e^{x_1} (1 + e^{x_2-x_1}) \right) = x_1 + \log\left(1 + e^{x_2-x_1}\right).\] Let \(y = x_2 - x_1 < 0\). Therefore \(e^y\) is a number smaller than 1, and we can use our numerical evaluation of the function \(f(y) = \log (1 + \exp(y)\) to enumerate it.

In Octave, one could accomplish this with:

%
% Compute F = log( exp(x1) + exp(x2) )
%
function F = LogSumExp(x1, x2)

  m = max(x1,x2);
  x = min(x1,x2) - m;
  F = m + ApproxLogP1(x);

endfunction

This yields the following algorithm for solving the general problem, \[F = \log \left(e^{x_1} + e^{x_2} + ... + e^{x_N} \right),\] once one recognizes \(F = \log e^F\).

%
% Method: Protected Sum of Logs
%
F = x(1);
for i = 2:n
  F = LogSumExp(F, x(i));
endfor

When we evaluate the specific problem  \[F = \log \left(e^{-1000} + e^{-1001} + ... + e^{-1100} \right),\] we get F = -999.54 as expected.