Stable Evaluation of Log Integrals with Log-Sum-Exp

Why Evaluate an Integral on the Log Scale?

A GLMM marginal likelihood often involves integrating a conditional likelihood over random effects. In a numerical calculation, that integral is replaced by a finite sum of weighted contributions. The contributions can be extremely small: multiplying many conditional probabilities can make an integrand’s value far below the range represented by ordinary floating-point arithmetic. Conversely, some intermediate exponentials can be so large that they overflow. Either failure can occur even when the final log of the integral is a perfectly useful finite number.

The remedy for a sum of positive terms is to calculate its logarithm without first calculating the terms themselves. The central identity is called the log-sum-exp identity. It is an algebraic rearrangement, not an approximation to the integral: the quadrature rule still determines which points and weights are used, while log-sum-exp makes the resulting positive sum numerically safer.

Definition: For a finite set of real numbers \(a_1,\ldots,a_K\), the log-sum-exp is \( \operatorname{LSE}(a_1,\ldots,a_K)=\log\left(\sum_{k=1}^{K}\exp(a_k)\right) \). If \(m=\max_k a_k\), its stable form is \( \operatorname{LSE}(a_1,\ldots,a_K)=m+\log\left(\sum_{k=1}^{K}\exp(a_k-m)\right) \).

From a GLMM Integral to a Positive Sum

For cluster \(j\), write its marginal likelihood contribution as an integral over its random effects \(b_j\):

L_j(\theta)=\int h_j(b_j;\theta)\,db_j, \qquad h_j(b_j;\theta)\geq 0.

Here \(\theta\) denotes the model parameters, and \(h_j\) includes the conditional likelihood and the random-effects density. A quadrature rule approximates this integral by evaluating the integrand at nodes \(b_{jk}\) and combining those values with weights \(\omega_{jk}\):

Q_j(\theta)=\sum_{k=1}^{K}\omega_{jk}h_j(b_{jk};\theta).

For the direct log-sum-exp calculation, the contributions must be nonnegative. In the standard setting here, the quadrature weights are positive and the integrand is nonnegative. Define the log contribution \(a_{jk}=\log(\omega_{jk})+\log h_j(b_{jk};\theta)\). Then the log of the quadrature approximation is

\log Q_j(\theta)=\operatorname{LSE}(a_{j1},\ldots,a_{jK}).

This formula means that neither the tiny integrand values nor their weighted products need to be formed on the ordinary scale. A zero weight or zero integrand contributes zero to the sum and can be represented by a log contribution of \(-\infty\), since \(\exp(-\infty)=0\). The formula assumes that at least one contribution is positive if a finite log result is required.

Why Maximum-Shifting Works

Let \(m=\max_k a_k\). Since \(a_k=m+(a_k-m)\), exponentiation gives \(\exp(a_k)=\exp(m)\exp(a_k-m)\). Summing and then taking the logarithm yields

\begin{aligned} \log\left(\sum_{k=1}^{K}\exp(a_k)\right) &=\log\left(\exp(m)\sum_{k=1}^{K}\exp(a_k-m)\right)\\ &=m+\log\left(\sum_{k=1}^{K}\exp(a_k-m)\right). \end{aligned}

Every shifted exponent satisfies \(a_k-m\leq 0\), so \(\exp(a_k-m)\leq 1\). The largest term becomes exactly 1, and no exponential of a large positive number is needed. Contributions far below the maximum may still round to zero after shifting, but then they are negligible relative to the maximum at the machine’s numerical precision. Shifting does not alter the mathematical sum.

Form log contributions.
For each positive term, add the log weight to the log integrand. Do not exponentiate the log integrand first.

Find the maximum.
Set \(m\) equal to the largest log contribution. This is a comparison on the log scale.

Sum shifted exponentials.
Calculate \(S=\sum_k\exp(a_k-m)\). The summands lie between zero and one.

Restore the scale in the logarithm.
Return \(m+\log(S)\). This is the log of the original positive sum.

Worked Example: Contributions Too Small to Exponentiate Directly

Suppose a hypothetical one-dimensional quadrature calculation has three log contributions \(a_1=-1000.2\), \(a_2=-1001.0\), and \(a_3=-1003.5\). Find the log of their sum.

The largest value is \(m=-1000.2\). The shifted values are \(0\), \(-0.8\), and \(-3.3\), so

\begin{aligned} \log\left(\sum_{k=1}^{3}\exp(a_k)\right) &=-1000.2+\log\left(1+\exp(-0.8)+\exp(-3.3)\right)\\ &=-1000.2+\log(1+0.4493+0.03688)\\ &=-1000.2+\log(1.4862)\\ &\approx -999.8038. \end{aligned}

The log of the sum is approximately \(-999.8038\). Directly computing \(\exp(-1000.2)\), \(\exp(-1001.0)\), and \(\exp(-1003.5)\) in standard floating-point arithmetic can underflow to zero, losing the sum. Maximum-shifting instead uses moderate exponentials and retains the finite log result.

Include Weights on the Log Scale

A frequent source of error is to stabilize the integrand values but forget that the quadrature weights also multiply them. The correct log contribution for a positive weight and positive integrand value is \(a_k=\log(\omega_k)+\log h(b_k;\theta)\). The weight and integrand logs are added because the logarithm turns multiplication into addition. The quadrature sum itself is still formed by log-sum-exp, not by adding the individual log contributions.

Worked Example: Combining Log Weights and Log Integrands

In a hypothetical quadrature calculation, three positive weights are \(0.25\), \(0.50\), and \(0.25\). The corresponding log integrand values are \(-990\), \(-1000\), and \(-995\). Calculate the log quadrature sum.

The log contributions, including the weights, are

\begin{aligned} a_1&=\log(0.25)-990=-991.3863,\\ a_2&=\log(0.50)-1000=-1000.6931,\\ a_3&=\log(0.25)-995=-996.3863. \end{aligned}

The maximum is \(m=-991.3863\). The differences from it are \(0\), \(-9.3068\), and \(-5.0000\). Therefore,

\begin{aligned} \log Q &=-991.3863+\log\left(1+\exp(-9.3068)+\exp(-5)\right)\\ &\approx -991.3863+\log(1+0.0000908+0.006738)\\ &\approx -991.3863+0.006806\\ &\approx -991.3795. \end{aligned}

Thus the log of the weighted quadrature sum is approximately \(-991.3795\). The calculation has included all three weights, and no value such as \(\exp(-1000)\) was needed. Omitting the log weights would calculate a different sum and generally a different answer.

Keep Cluster Integrals Separate

When clusters are conditionally independent given the model parameters, the marginal likelihood is the product of their marginal likelihood contributions. Consequently, the log likelihood is the sum of the cluster log likelihoods. Each cluster’s quadrature sum must first be evaluated separately. Applying log-sum-exp to every node contribution across all clusters would instead calculate a sum of contributions, not the product of cluster integrals.

Worked Example: Combining Two Cluster Log Integrals

Consider hypothetical quadrature calculations for two independent clusters. Cluster A has log contributions \(-100\) and \(-101\); cluster B has log contributions \(-2\) and \(-4\). Find the combined log likelihood contribution.

For cluster A, shift by its largest log contribution, \(-100\):

\log Q_A=-100+\log(1+\exp(-1)) =-100+0.3133=-99.6867.

For cluster B, shift by \(-2\):

\log Q_B=-2+\log(1+\exp(-2)) =-2+0.1269=-1.8731.

The combined log likelihood contribution is the sum of these cluster log integrals:

\log Q_A+\log Q_B=-99.6867-1.8731=-101.5598.

This calculation preserves the product structure across clusters. In contrast, taking one log-sum-exp over all four node contributions would give the logarithm of a sum of four terms; that is not the log of the product \(Q_AQ_B\).

Diagnostics, Alternatives, and Limits

Log-sum-exp addresses a specific computational problem: representing the sum of positive terms without forming potentially extreme exponentials. It does not guarantee that the quadrature rule approximates the integral accurately, that the model is well specified, or that the log contributions were constructed correctly. In particular, the numerical stability of a sum should not be confused with the accuracy of the integral approximation.

Approach	What it does	Main concern
Direct-scale summation	Exponentiates each log contribution and adds the results.	Terms may underflow or overflow before summation.
Maximum-shifted log-sum-exp	Subtracts the largest log contribution before exponentiating, then restores it in the final log.	Stabilizes a positive sum but cannot correct inaccurate nodes, weights, or integrand values.
Laplace approximation	Approximates an integral using local behavior near a mode.	It is an integration approximation, not an algebraic method for summing a given set of positive quadrature contributions.

Several checks help catch implementation errors:

Check signs and weights. The stated identity applies to a sum of nonnegative terms. If a rule uses negative weights, taking logs of the terms is not valid; signed-sum methods are needed instead.
Check the maximum shift. Each finite shifted value \(a_k-m\) must be at most zero, and at least one must equal zero.
Check zero and invalid terms separately. A genuine zero contribution has log value \(-\infty\). A NaN or an unexpected infinite value is not interchangeable with a legitimate zero and should trigger investigation.
Check the scale of the result. Since the sum includes the largest positive term, its logarithm cannot be less than \(m\). A computed value below \(m\) signals an arithmetic or coding error.
Check integration accuracy independently. A stable result can still be a poor approximation if important regions of the integrand are not represented by the quadrature points. Stability does not resolve multimodality or inadequate coverage.

A common mistake is to calculate \(\log(\sum_k a_k)\) when \(a_k\) are already log contributions, or to use \(\sum_k\log(a_k)\). Neither expression is the log of the desired sum. The correct expression is \(\log(\sum_k\exp(a_k))\), evaluated by shifting the largest \(a_k\). Another mistake is to treat a very small shifted exponential that rounds to zero as evidence that the original term was exactly zero. It may simply be negligible at the available precision; the distinction can matter when diagnosing the model or the computation.

Key takeaway: For positive quadrature contributions, add log weights to log integrand values, subtract the largest resulting log contribution before exponentiating, and then restore that maximum. This prevents range failures in the sum while leaving quadrature accuracy as a separate question.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

For log contributions \(-20\), \(-21\), and \(-23\), identify the shift and write the stable log-sum-exp expression.
Why does subtracting the largest log contribution avoid overflow in the exponential calculations?
A quadrature term has weight \(0.2\) and log integrand value \(-50\). Write its log contribution.
Two independent clusters have quadrature log integrals \(-3.2\) and \(-4.1\). What is their combined log likelihood contribution, and why is it not obtained by applying log-sum-exp to the two cluster log integrals?

Solutions. (1) The maximum is \(m=-20\), so the result is \(-20+\log(1+\exp(-1)+\exp(-3))\). (2) Every shifted value is at most zero, so each exponential is at most one; a large positive exponent is never formed. (3) The log contribution is \(\log(0.2)-50\approx -51.6094\). (4) The combined contribution is \(-3.2+(-4.1)=-7.3\), because independent cluster likelihoods multiply and their logarithms add. Log-sum-exp would instead represent the log of a sum, not the log of that product.