# numerically stable log1pexp calculation

What are good approximations for computing log1pexp for single precision and double precision floating point numbers?

Note: log1pexp(x) is log(1 + exp(x))

I have found few implementations of log1pexp for double precision but they don't provide an explanation on how they arrived at the approximations. Hence, I am not able to implement log1exp for single precision numbers (without converting to double precision intermediates of course).

Reference implementation: https://github.com/davisking/dlib/blob/master/dlib/dnn/utilities.h#L16-L29

• I don't know if this should be at stackoverflow or codereview. I choose CSSE because this had to do with a bit of theory and mathematics. – Yashas Jun 17 at 12:03
• It is on-topic on this site to ask for an algorithm or high level description to obtain an approximation, but not for an implementation. So I think here is indeed more appropriate than Stack Overflow or Code Review. – Discrete lizard Jun 17 at 12:20

Let $$0 < \varepsilon \lll 1$$ be the relative error bound of the floating-point system—$$2^{-53}$$ in IEEE 754 binary64 arithmetic.

First, the naive formula log1p(exp(x)) always gives a good approximation unless exp(x) overflows: If exp(x) computes $$(1 + \delta_1) e^x$$, and if log1p(y) computes $$(1 + \delta_2) \log(1 + y)$$, where $$|\delta_1|, |\delta_2| \leq \varepsilon$$, then log1p(exp(x)) yields

$$$$(1 + \delta_2) \log\bigl[1 + (1 + \delta_1) e^x\bigr] = (1 + \delta_2) (1 + \delta_1') \log(1 + e^x)$$$$

where $$|\delta_1'|$$ is bounded by a constant multiple of $$\delta_1$$ since $$(1 + \delta_1) e^x$$ is positive and $$y \mapsto \log(1 + y)$$ has small condition number on the positive real axis. Thus the relative error is bounded by $$|\delta_2| + |\delta_1'| + |\delta_2 \delta_1'|$$ which in turn is bounded by small polynomial function of $$\varepsilon$$.

Exercise: Compute a bound on $$|\delta_1'/\delta_1|$$.

But you can do better, with less computation (only one transcendental evaluation, not two), in some special cases:

1. If $$x \lll 0$$, then $$e^x$$ is near 0 so $$\log(1 + e^x) \approx e^x$$. Specifically, $$0 < e^x < 1$$, so:

• $$\log(1 + e^x) \geq e^x/2$$, since $$y \mapsto \log(1 + y)$$ is a concave function with derivative $$y \mapsto 1/(1 + y)$$ bounded below by $$1/2$$ for $$y < 1$$; and
• the Taylor expansion of $$y \mapsto \log(1 + y)$$ converges absolutely at $$e^x$$.

Thus the relative error of $$e^x$$ from $$\log(1 + e^x)$$ is:

\begin{align} \frac{|e^x - \log(1 + e^x)|}{\left|\log(1 + e^x)\right|} &\leq 2 e^{-x} |e^x - \log(1 + e^x)| \\ &= 2 e^{-x} |e^x - (e^x - e^{2x}/2 + e^{3x}/3 - \cdots)| \\ &= 2 e^{-x} |e^{2x}/2 - e^{3x}/3 + e^{4x}/4 - \cdots| \\ &= 2 |e^x/2 - e^{2x}/3 + e^{3x}/4 - \cdots| \\ &< 2 |e^x/2| \\ &= e^x. \end{align}

For $$x \leq \log \varepsilon$$, this relative error is bounded by $$\varepsilon$$, which justifies the formula exp(x).

In IEEE 754 binary64 arithmetic, $$\log \varepsilon \approx -37$$.

2. If $$x \gg 0$$, we can use the identity

$$$$\log(1 + e^x) = \log\bigl[e^x (1 + e^{-x})\bigr] = x + \log(1 + e^{-x})$$$$

to rewrite $$\log(1 + e^x)$$ in terms of $$\log(1 + y)$$ for $$0 < y < 1$$, with $$y = e^{-x}$$, as in case (1). The value of $$\log(1 + e^x)$$ lies between $$x + e^{-x}$$ and $$x + e^{-x} - e^{-2x}/2$$, so the relative error of $$x + e^{-x}$$ from $$\log(1 + e^x)$$ is below the relative error of $$x + e^{-x}$$ from $$x + e^{-x} - e^{-2x}/2$$, which is

$$$$\frac{|x + e^{-x} - (x + e^{-x} - e^{-2x}/2)|}{|x + e^{-x} - e^{-2x}/2|} = \frac{e^{-2x}/2}{|x + e^{-x} - e^{-2x}/2|} \leq e^{-2x}/2.$$$$

For $$x \geq \log(1/\sqrt{2\varepsilon})$$, this is bounded by $$\varepsilon$$, which justifies the formula x + exp(-x).

In IEEE 754 binary64 arithmetic, $$\log(1/\sqrt{2\varepsilon}) \approx 18$$.

3. If $$x \ggg 0$$, then $$x + \log(1 + e^{-x})$$ may simply be rounded to $$x$$. Specifically, the relative error of $$x$$ from $$\log(1 + e^x) = x + \log(1 + e^{-x})$$ is

$$$$\frac{\bigl|x - \bigl[x + \log(1 + e^{-x})\bigr]\bigr|}{\left|\log(1 + e^x)\right|} = \frac{\log(1 + e^{-x})}{\log(1 + e^x)} < \log(1 + e^{-x}),$$$$

as long as $$e^x \geq e - 1$$ so that the denominator $$\log(1 + e^x)$$ is at least $$1$$.

For $$x \geq \log[1/(e^\varepsilon - 1)]$$, then this error is below $$\varepsilon$$, which justifies the formula x.

In IEEE 754 binary64 arithmetic, $$\log[1/(e^\varepsilon - 1)] \approx 37$$.

Exercise: Find a nice formula for a bound on $$x$$ that corresponds to the cutoff $$33.3$$ in the source you cited, by taking the denominator $$\log(1 + e^x)$$ into account.