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:

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\).

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

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.

## No comments:

Post a Comment