Many models in machine learning and statistics multiply positive quantities like likelihood terms, transition probabilities, priors, and partition functions. In finite-precision arithmetic, these products can quickly overflow or underflow. The standard remedy is to work in log-space.
For positive values \(a,b\), \[ \log(ab)=\log a+\log b, \qquad \log\left(\frac{a}{b}\right)=\log a-\log b, \qquad \log(a^r)=r\log a. \] So multiplication, division, and powers become numerically stable additions and subtractions.
But one operation remains nontrivial: \[ \log(u+v)=\log\left(e^x+e^y\right) \quad\text{for}\quad x=\log u,\ y=\log v. \] This expression is the log-sum-exp operation.
Definition
For a vector \(\mathbf{x}=(x_1,\ldots,x_n)\), define
\[ \boxed{ \operatorname{LSE}(\mathbf{x}) := \log\left(\sum_{i=1}^{n}e^{x_i}\right) }. \]
If \(x_i=\log u_i\) with \(u_i>0\), then \[ \operatorname{LSE}(x_1,\ldots,x_n) = \log(u_1+\cdots+u_n). \] So LSE is addition in the logarithmic domain.
Why Naive Evaluation Fails
The normal formula \[ \operatorname{LSE}(\mathbf{x})=\log\left(\sum_i e^{x_i}\right) \] suggests to exponentiate, sum and then take the logarithm. In floating-point arithmetic thisin can fail in both directions, overflow and underflow.
Overflow
With \(\mathbf{x}=(1000,999,998)\), the true value is finite: \[ \log\left(e^{1000}+e^{999}+e^{998}\right)\approx 1000.4076059644444. \] But \(e^{1000}\) overflows in double precision, so the naive pipeline returns \(+\infty\).
Underflow
With \(\mathbf{x}=(-1000,-1001)\), both exponentials underflow to zero in ordinary double arithmetic, producing \[ \log(0+0)=-\infty, \] although the correct value is \[ \operatorname{LSE}(-1000,-1001)\approx -999.6867383124818. \]
Derivation of the Trick
Let \[ y=\log\left(\sum_{i=1}^{n}e^{x_i}\right). \] For any finite shift \(a\), write \(e^{x_i}=e^a e^{x_i-a}\). Then \[ \sum_i e^{x_i}=e^a\sum_i e^{x_i-a}, \] hence \[ \boxed{ \operatorname{LSE}(\mathbf{x}) = a+\log\left(\sum_{i=1}^{n}e^{x_i-a}\right) }. \] This identity is exact.
Proof: \[ \begin{array}{rrl} & y = & \log\sum_i e^{x_i} \\ \Leftrightarrow& e^y = & \sum_i e^{x_i} \\ \Leftrightarrow& e^{-a}e^{y} = & e^{-a}\sum_i e^{x_i} \\ \Leftrightarrow& e^{y-a} = & \sum_i e^{x_i-a} \\ \Leftrightarrow& y-a = & \log\sum_i e^{x_i-a} \\ \Leftrightarrow& y = & a+\log\sum_i e^{x_i-a}. \end{array} \]
Finding a Numerically Optimal Shift
Choose \[ m=\max_i x_i. \] Then each shifted exponent satisfies \[ x_i-m\le 0, \] so \[ 0<e^{x_i-m}\le 1, \] and at least one term is \(e^0=1\).
Therefore the stable formula is \[ \boxed{ \operatorname{LSE}(\mathbf{x}) = m+\log\left(\sum_{i=1}^{n}e^{x_i-m}\right), \qquad m=\max_i x_i }. \]
The shifted sum is bounded: \[ 1\le\sum_i e^{x_i-m}\le n. \] So it cannot overflow, and it cannot underflow to zero because one term is exactly 1.
Stable Overflow Example
For \((1000,999,998)\), \(m=1000\), so \[ \operatorname{LSE} =1000+\log\left(1+e^{-1}+e^{-2}\right) \approx 1000.4076059644444. \]
Stable Underflow Example
For \((-1000,-1001)\), \(m=-1000\), so \[ \operatorname{LSE} =-1000+\log\left(1+e^{-1}\right) \approx -999.6867383124818. \]
LSE as Maximum Plus Correction
The shifted form of the LSE is not only a stability trick. It also decomposes LSE into a dominant scale term and a residual that measures how much other entries still contribute beyond the maximum. This view is useful in optimization, probability normalization, and "soft argmax" intuition.
The stable form can be read as \[ \operatorname{LSE}(\mathbf{x}) = m + \underbrace{\log\left(\sum_i e^{x_i-m}\right)}_{\text{correction term}}. \] Since \(1\le\sum_i e^{x_i-m}\le n\), \[ 0\le \text{correction}\le \log n, \] so \[ \boxed{m\le\operatorname{LSE}(\mathbf{x})\le m+\log n}. \]
This is the key interpretation: \(\operatorname{LSE}\) is a differentiable (smooth) surrogate for \(\max\). In practice, you can treat \(\operatorname{LSE}(\mathbf{x})\) as "almost the maximum" plus a bounded nonnegative offset.
Two limiting cases make this more clear: \[ \text{if } x_k \gg x_j\ (j\neq k),\quad \operatorname{LSE}(\mathbf{x})\approx x_k=\max_i x_i, \] because the correction term is close to \(0\). On the other hand, if all values are equal, \(x_i=c\), then \[ \operatorname{LSE}(\mathbf{x})=c+\log n, \] which is the maximal possible gap to \(\max\) under the bound above.
Binary Form (logaddexp)
For two values \(u,v\), \[ \operatorname{LSE}(u,v)=\log(e^u+e^v). \] Let \(m=\max(u,v)\). Then \[ \operatorname{LSE}(u,v) = m+\log\left(e^{u-m}+e^{v-m}\right). \] Since one of \(u-m\) or \(v-m\) is \(0\) (which can also be seen as \(\max(u,v)-\max(u,v)=0\)) and the other is \[ \min(u,v)-\max(u,v)=-|u-v|, \] we get \[ e^{u-m}+e^{v-m}=1+e^{-|u-v|}. \] A symmetric stable formula is therefore \[ \boxed{ \operatorname{LSE}(u,v) = \max(u,v)+\log\left(1+e^{-\lvert u-v\rvert}\right) }. \]
Note: Numerical libraries use \(\operatorname{log1p}(t)=\log(1+t)\) for better accuracy for near zero values.
Softmax and Log-Softmax
For logits \(x_i\), the softmax produces normalized probabilities \[ p_i=\frac{e^{x_i}}{\sum_j e^{x_j}}, \qquad \log p_i=x_i-\operatorname{LSE}(\mathbf{x}). \]
With \(m=\max_j x_j\), the stable forms are \[ p_i=\frac{e^{x_i-m}}{\sum_j e^{x_j-m}}, \qquad \log p_i=(x_i-m)-\log\sum_j e^{x_j-m}. \] The common shift cancels in softmax, so subtracting \(m\) does not change the probabilities.
Translation Invariance
For any constant \(c\), \[ \operatorname{LSE}(\mathbf{x}+c\mathbf{1})=c+\operatorname{LSE}(\mathbf{x}), \] and \[ \operatorname{softmax}(\mathbf{x}+c\mathbf{1})=\operatorname{softmax}(\mathbf{x}). \]
Gradient, Hessian, Convexity
Let \( f(\mathbf{x})=\operatorname{LSE}(\mathbf{x}) \). Then \[ \frac{\partial f}{\partial x_i}= \frac{e^{x_i}}{\sum_j e^{x_j}} =\operatorname{softmax}(\mathbf{x})_i. \] Hence \[ \boxed{\nabla\operatorname{LSE}(\mathbf{x})=\operatorname{softmax}(\mathbf{x})}. \]
Writing \(\mathbf{p}=\operatorname{softmax}(\mathbf{x})\), \[ \nabla^2\operatorname{LSE}(\mathbf{x})=\operatorname{diag}(\mathbf{p})-\mathbf{p}\mathbf{p}^{\mathsf T}. \] For any \(\mathbf{v}\), \[ \mathbf{v}^{\mathsf T}\nabla^2\operatorname{LSE}(\mathbf{x})\mathbf{v} = \operatorname{Var}_{i\sim \mathbf{p}}(v_i) \ge 0, \] so LSE is convex.
Softplus as a Special Case
Softplus is the two-term LSE where one logit is fixed at zero. This makes it a smooth replacement for ReLU and a canonical example of how LSE turns hard piecewise maxima into differentiable transitions.
\[ \operatorname{softplus}(x)=\log(1+e^x)=\operatorname{LSE}(0,x). \] Stable form: \[ \operatorname{softplus}(x) = \max(0,x)+\log\left(1+e^{-\lvert x\rvert}\right). \]
Temperature-Scaled Smooth Maximum
A temperature parameter controls how sharply LSE tracks the hard maximum. Large temperatures average more broadly; small temperatures concentrate mass near the largest component. This scaling is central in calibration, sampling, and entropy-regularized objectives.
\[ \operatorname{LSE}_{\tau}(\mathbf{x}) := \tau\log\left(\sum_i e^{x_i/\tau}\right), \qquad \tau>0. \] With \(m=\max_i x_i\): \[ m\le \operatorname{LSE}_{\tau}(\mathbf{x})\le m+\tau\log n. \] As \(\tau\to 0^+\), \(\operatorname{LSE}_{\tau}(\mathbf{x})\to\max_i x_i\).
Log-Sum-Exp Variants
Two common variants appear frequently in practice: log-mean-exp for normalized aggregation and weighted LSE for mixture-like sums or prior-weighted evidence terms.
Log-Mean-Exp
If \(x_i=\log u_i\), then \[ \log\left(\frac1n\sum_i u_i\right) = \operatorname{LSE}(\mathbf{x})-\log n. \] This is often called log-mean-exp.
Weighted LSE
For weights \(w_i\ge 0\), \[ \log\left(\sum_i w_i e^{x_i}\right) = \operatorname{LSE}(x_1+\log w_1,\ldots,x_n+\log w_n), \] with \(w_i=0\Rightarrow \log w_i=-\infty\).
Applications
Naive Bayes Posterior Normalization
For a Naive Bayes classifier with classes \(C_1,\ldots,C_K\) and observed features \(\mathbf{x}=(x_1,\ldots,x_d)\), the unnormalized posterior score for class \(C_k\) is \[ P(C_k)\prod_{r=1}^{d}P(x_r\mid C_k). \] Multiplying many probabilities can underflow, so we move to log-space and define \[ a_k=\log P(C_k)+\sum_r \log P(x_r\mid C_k). \] Hence \[ e^{a_k}=P(C_k)\prod_{r=1}^{d}P(x_r\mid C_k). \] Bayes normalization becomes \[ P(C_k\mid\mathbf{x})=\frac{e^{a_k}}{\sum_{j=1}^{K}e^{a_j}}. \] In log-space: \[ \log P(C_k\mid\mathbf{x})=a_k-\operatorname{LSE}(a_1,\ldots,a_K). \]
A useful practical distinction:
- For MAP class prediction only, normalization is unnecessary: \(\widehat{k}=\arg\max_k a_k\).
- For calibrated posterior probabilities, the denominator is required and should be computed via LSE.
Multinomial Logistic Regression and Cross-Entropy
Let a classifier produce logits \(z_1,\ldots,z_K\). Softmax gives class probabilities \[ p_k=\frac{e^{z_k}}{\sum_{j=1}^{K}e^{z_j}}, \qquad \log p_k=z_k-\operatorname{LSE}(\mathbf{z}). \] If the correct class is \(y\), then the negative log-likelihood (multiclass cross-entropy) is \[ \mathcal{L}=-\log p_y=\operatorname{LSE}(\mathbf{z})-z_y. \] Stable evaluation uses \(m=\max_j z_j\): \[ \mathcal{L}=m-z_y+\log\left(\sum_{j=1}^{K}e^{z_j-m}\right). \] So there is no reason to evaluate potentially enormous terms \(e^{z_j}\) directly.
Differentiation gives the standard gradient \[ \frac{\partial\mathcal{L}}{\partial z_k}=p_k-\mathbb{1}[k=y]. \]
Bayesian Filtering and Hidden-State Models
Consider discrete hidden states \(h_t\) and observations \(v_t\). The forward quantity for state \(j\) is \[ \alpha_t(j)=P(v_{1:t},h_t=j). \] The standard recursion is \[ \alpha_t(j)=P(v_t\mid h_t=j)\sum_i\alpha_{t-1}(i)P(h_t=j\mid h_{t-1}=i). \] The outer multiplication is easy in log-space, but the predecessor-state sum requires log-sum-exp.
Define \(\ell_t(j)=\log\alpha_t(j)\). Then \[ \ell_t(j) = \log P(v_t\mid h_t=j) + \log\left(\sum_i \exp\left(\ell_{t-1}(i)+\log P(h_t=j\mid h_{t-1}=i)\right)\right). \] Therefore \[ \ell_t(j) = \log P(v_t\mid h_t=j) + \operatorname{LSE}_i\left(\ell_{t-1}(i)+\log P(h_t=j\mid h_{t-1}=i)\right). \] This is the log-domain forward algorithm form. The same pattern appears in HMMs, discrete Bayesian filters, switching or mixture Kalman filters, mixture filters, particle-weight normalization, and probabilistic sequence models.
Mixture Models
For a mixture distribution \[ p(x)=\sum_{k=1}^{K}\pi_k p_k(x), \qquad \pi_k\ge 0, \qquad \sum_k\pi_k=1, \] the log-likelihood is \[ \log p(x)=\log\left(\sum_{k=1}^{K}\pi_k p_k(x)\right). \]
Define component log-joint scores \[ a_k=\log\pi_k+\log p_k(x). \] Then \[ \log p(x)=\operatorname{LSE}(a_1,\ldots,a_K), \qquad \log P(k\mid x)=a_k-\operatorname{LSE}(\mathbf{a}), \qquad P(k\mid x)=\operatorname{softmax}(\mathbf{a})_k. \] This is why LSE repeatedly appears in Gaussian mixtures and EM-style responsibilities.
Attention Normalization
In scaled dot-product attention, one query produces a score vector \( \mathbf{s}=(s_1,\ldots,s_T) \) over \(T\) keys, typically with \( s_j=\frac{q^\top k_j}{\sqrt{d_k}} \). The attention weights are a softmax over these logits: \[ p_i=\frac{e^{s_i}}{\sum_{j=1}^{T}e^{s_j}}. \]
In log-space this is exactly an LSE normalization: \[ \log p_i = s_i-\operatorname{LSE}(\mathbf{s}). \] Directly computing \(e^{s_j}\) can overflow when logits are large, so we use the max-shift form with \(m=\max_j s_j\): \[ p_i = \frac{e^{s_i-m}}{\sum_{j=1}^{T}e^{s_j-m}}, \qquad \log p_i = (s_i-m)-\log\sum_{j=1}^{T}e^{s_j-m}. \]
This is the same LSE trick as in softmax classifiers, now applied row-wise to every query position. It keeps attention probabilities finite and stable even when score magnitudes vary strongly across sequence length, masking patterns, or precision formats.
Online and Blockwise LSE
Suppose a block \(A\) is represented by \( m_A=\max_{i\in A} x_i \) and \( \ell_A=\sum_{i\in A} e^{x_i-m_A} \), and similarly for block \(B\). For \(m=\max(m_A,m_B)\), merged statistics are \[ \ell=e^{m_A-m}\ell_A+e^{m_B-m}\ell_B, \qquad \operatorname{LSE}(A\cup B)=m+\log\ell. \]
The same merge principle extends to normalized weighted sums (for example in tiled attention), where numerator and denominator are accumulated in the same shifted scale.
Important Edge Cases
- All inputs \(-\infty\): \(\operatorname{LSE}=-\infty\) (since sum is zero).
- Any input \(+\infty\): \(\operatorname{LSE}=+\infty\).
- Empty set: by convention \(\operatorname{LSE}(\varnothing)=\log 0=-\infty\).
- Zero probabilities map to \(-\infty\), which is compatible because \(e^{-\infty}=0\).