Stable Sigmoid: Avoiding Overflow in Exponentials, Log-Sigmoid, and Losses
In exact real arithmetic, $\sigma(z)=1/(1+e^{-z})$ is always well-defined and strictly between zero and one. In floating-point arithmetic, the intermediate exponential $e^{-z}$ can overflow long before the final output should have been a harmless number near zero or one.
Where the naive formula fails
Every IEEE-754 floating-point format has a largest finite value, so $\exp(u)$ overflows to $+\infty$ above a format-dependent threshold: roughly $u>11.09$ in float16, $u>88.72$ in float32, and $u>709.78$ in float64.
When a logit is moderately negative—say $z=-100$ in float32 or $z=-1000$ in float64—evaluating $1/(1+e^{-z})$ asks for $e^{+100}$ or $e^{+1000}$, raising an overflow warning. Conversely, the algebraically equivalent fraction $e^z/(1+e^z)$ is safe for negative $z$, but overflows in both numerator and denominator for large positive $z$, producing $\infty/\infty=\mathrm{NaN}$.
Keep the exponent nonpositive
Multiply the numerator and denominator of $1/(1+e^{-z})$ by $e^z$ to get $e^z/(1+e^z)$. Choosing the first form when $z\geq0$ and the second when $z<0$ ensures that the exponential argument is never positive:
In both branches, the exponent is $-|z|\leq0$, so $e^{-|z|}$ lies in $(0,1]$ and never overflows in any floating-point precision. The denominator $1+e^{-|z|}$ stays between $1$ and $2$. Using $\tanh(z/2)$ expresses the same stabilization via the half-angle identity $\sigma(z)=\tfrac{1}{2}[1+\tanh(z/2)]$.
The derivative also benefits from writing the exponential in terms of $-|z|$:
When $z$ is positive and moderately large (for example $z=20$ in float32), $\sigma(z)$ rounds to $1.0$, so $\sigma(z)(1-\sigma(z))$ rounds directly to zero. Evaluating $e^{-|z|}/(1+e^{-|z|})^2$ preserves non-zero tail gradients until $e^{-|z|}$ itself underflows.
Why log-sigmoid needs its own primitive
Making $\sigma(z)$ overflow-free does not make $\log(\sigma(z))$ safe. At $z=-100$ in float32 (or $z=-1000$ in float64), even the stable branch $e^z/(1+e^z)$ underflows to $0.0$, and taking $\log(0.0)$ returns $-\infty$. Yet the true value of $\log\sigma(z)$ is smooth and nearly linear in $z$:
For $z=-100$, $-\max(0,100)=-100$ and $\log(1+e^{-100})\approx0$, so $\log\sigma(-100)\approx-100$ with no overflow or $\log(0)$. By symmetry, $1-\sigma(z)=\sigma(-z)$, so $\log(1-\sigma(z))=\log\sigma(-z)$.
| Quantity | −100 | −1 | 0 | 1 | 100 |
|---|---|---|---|---|---|
| $\sigma(z)$ | 0.0000 | 0.2689 | 0.5000 | 0.7311 | 1.0000 |
| $\log\sigma(z)$ | -100.0000 | -1.3133 | -0.6931 | -0.3133 | 0.0000 |
| BCE ($y=1$) | 100.0000 | 1.3133 | 0.6931 | 0.3133 | 0.0000 |
Combine sigmoid and binary cross-entropy
For a binary target $y\in[0,1]$ and logit $z$, binary cross-entropy is $-y\log\sigma(z)-(1-y)\log(1-\sigma(z))$. Substitute $\log\sigma(z)=-\log(1+e^{-z})$ and $\log(1-\sigma(z))=-z-\log(1+e^{-z})$ to obtain a single stable formula:
This is the exact identity inside PyTorch's BCEWithLogitsLoss. The exponential only ever sees $-|z|\leq0$, and $\operatorname{log1p}(e^{-|z|})$ stays accurate when $e^{-|z|}$ is tiny.
Advantages and disadvantages
| Formulation | Advantages | Disadvantages |
|---|---|---|
| Naive $1/(1+e^{-z})$ | Matches the textbook formula directly and is accurate for nonnegative $z$. | Overflows in $e^{-z}$ for large negative $z$; computing $\log(\sigma(z))$ afterward still fails when $\sigma(z)$ rounds to 0 or 1. |
Two-branch $\sigma(z)$ (expit) | Keeps the exponential argument nonpositive ($-|z|\leq0$), avoiding overflow across all inputs. | Requires a branch or mask; extreme negative inputs still underflow to $0.0$, so taking $\log$ afterward can still produce $-\infty$. |
| Log-space $\log\sigma(z)$ and fused BCE | Preserves linear tail behavior ($\approx z$ for $z\ll0$) and avoids both $e^{-z}$ overflow and $\log(0)$ underflow. | Must receive raw logits $z$ directly rather than pre-computed probabilities $p=\sigma(z)$. |
Reference implementation
import numpy as np
def stable_sigmoid(z):
z = np.asarray(z, dtype=float)
pos = z >= 0
out = np.empty_like(z)
out[pos] = 1.0 / (1.0 + np.exp(-z[pos]))
ez = np.exp(z[~pos])
out[~pos] = ez / (1.0 + ez)
return out
def log_sigmoid(z):
return -np.logaddexp(0.0, -np.asarray(z, dtype=float))
def bce_with_logits(z, y):
z = np.asarray(z, dtype=float)
return np.maximum(z, 0.0) - z * y + np.log1p(np.exp(-np.abs(z)))In SciPy, scipy.special.expit implements the stable sigmoid primitive used by SiLU and Swish and SwiGLU. Next, ReLU and its variants replace bounded exponential tails with piecewise linear branches.