Deriving Adaptive Nodes and Weights from a Local Gaussian Approximation

Why Adapt the Nodes?

Gauss-Hermite quadrature is most efficient when its nodes cover the region where an integrand contributes most to an integral. In a GLMM, the integrand for a cluster’s marginal likelihood can be sharply concentrated away from zero and can have a width quite different from the random effect’s prior standard deviation. A rule fixed at zero with a fixed scale may then place too few nodes in the important region. Adaptive Gauss-Hermite quadrature addresses this by placing nodes near the log integrand’s conditional mode and scaling their spacing according to its local curvature.

The key point is that adaptation changes more than the node locations. Once the integration variable is shifted and scaled, the integrand must also be expressed relative to the Gaussian weight used by the quadrature rule. That division produces a residual factor, which adjusts every quadrature weight. Omitting it generally changes the integral being approximated.

Definition: Let \(I=\int_{-\infty}^{\infty}\exp\{h(b)\}\,db\), where \(h(b)\) is the log integrand and the integral is finite. Let \(\hat b\) be an interior mode of \(h\), and define \(H=-h''(\hat b)\). When \(H>0\), the local Gaussian approximation has location \(\hat b\) and scale \(H^{-1/2}\). Adaptive Gauss-Hermite quadrature places standard nodes around \(\hat b\) on this scale and multiplies the standard weights by the residual between the true integrand and this local Gaussian.

In a marginal-likelihood calculation, \(h(b)\) includes all terms in the log integrand, including the random-effects density and the conditional likelihood. Constants can either be included in \(h\) or kept outside the integral, but they must be handled consistently. The integration target is the same throughout the derivation; the Gaussian factor is a device for constructing an efficient numerical rule, not a replacement for the integrand.

Deriving Location, Scale, and Weight Adjustment

Assume \(h\) is twice differentiable near \(\hat b\), and that \(\hat b\) is an interior mode. Then \(h'(\hat b)=0\). A second-order Taylor expansion gives the local approximation

h(b)\approx h(\hat b)-\frac{H}{2}(b-\hat b)^2.

The Gaussian-shaped term in this expression suggests centering at \(\hat b\) and measuring distance in units of \(H^{-1/2}\). To obtain the quadrature formula, factor this term out without approximating the remainder:

\exp\{h(b)\} = \exp\{h(\hat b)\} \exp\left\{-\frac{H}{2}(b-\hat b)^2\right\} \exp\left\{h(b)-h(\hat b)+\frac{H}{2}(b-\hat b)^2\right\}.

The last exponential is the residual correction. It equals one at the mode, but away from the mode it records how the actual log integrand differs from its quadratic approximation. Define \(u=\sqrt{H/2}(b-\hat b)\), so that \(b=\hat b+\sqrt{2/H}\,u\) and \(db=\sqrt{2/H}\,du\). The Gaussian factor becomes \(\exp(-u^2)\), which is exactly the weight used in the standard Gauss-Hermite rule.

I= \exp\{h(\hat b)\}\sqrt{\frac{2}{H}} \int_{-\infty}^{\infty}e^{-u^2} \exp\left\{h\left(\hat b+\sqrt{\frac{2}{H}}u\right)-h(\hat b)+u^2\right\}\,du.

Let \(x_j\) and \(w_j\), for \(j=1,\ldots,m\), be the nodes and weights for the convention \(\int e^{-x^2}f(x)\,dx\approx\sum_{j=1}^m w_j f(x_j)\). Substitution gives the adaptive nodes \(b_j\) and the corresponding approximation:

b_j=\hat b+\sqrt{\frac{2}{H}}x_j,\qquad I\approx \exp\{h(\hat b)\}\sqrt{\frac{2}{H}} \sum_{j=1}^m w_j \exp\left\{h(b_j)-h(\hat b)+x_j^2\right\}.

Thus the adjusted weight at node \(j\) is proportional to \(w_j\exp\{h(b_j)-h(\hat b)+x_j^2\}\), with the common prefactor \(\exp\{h(\hat b)\}\sqrt{2/H}\). The \(+x_j^2\) term is essential: it divides out the Gaussian factor already represented by the Gauss-Hermite rule. The identity before applying the finite sum is exact; the approximation enters when the remaining weighted integral is replaced by a finite quadrature sum.

Curvature condition: At a differentiable interior maximum, \(h''(\hat b)\leq 0\), so \(H\geq0\). A maximum need not have \(H>0\): for example, \(h(b)=-b^4\) has a strict maximum at zero but \(H=0\) there. The stated Gaussian scale requires \(H>0\). In multiple dimensions, a maximum implies that the negative Hessian is positive semidefinite, while the usual local Gaussian scale requires it to be positive definite.

In \(d\) dimensions, let \(\hat{\boldsymbol b}\) be the mode and \(H\) the positive-definite negative Hessian there. Using a symmetric inverse square root \(H^{-1/2}\), the local transformation is \(\boldsymbol b=\hat{\boldsymbol b}+\sqrt{2}H^{-1/2}\boldsymbol x\). Its Jacobian is \(2^{d/2}|H|^{-1/2}\), and the Gaussian exponent becomes \(-\boldsymbol x^\mathsf{T}\boldsymbol x\). The residual correction is \(\exp\{h(\boldsymbol b)-h(\hat{\boldsymbol b})+\boldsymbol x^\mathsf{T}\boldsymbol x\}\). This displays the same location, scale, and weight adjustment in vector form; the one-dimensional derivation above is the direct calculation used in the examples.

Worked Calculations

Worked Example: An Exactly Gaussian Integrand

Suppose \(h(b)=-3-2(b-1)^2\). This hypothetical integrand has mode \(\hat b=1\), with \(h(\hat b)=-3\). Since \(h''(b)=-4\), its negative curvature is \(H=4\), so the adaptive scale is \(H^{-1/2}=1/2\).

At a standard node \(x_j\), the transformed node is \(b_j=1+\sqrt{2/4}\,x_j=1+x_j/\sqrt{2}\). The log residual at any such node is \(h(b_j)-h(\hat b)+x_j^2=-2(b_j-1)^2+x_j^2\). Because \((b_j-1)^2=x_j^2/2\), this residual is zero. Every adjusted weight is therefore just its standard weight, apart from the common prefactor.

For any Gauss-Hermite order that integrates a constant exactly, \(\sum_jw_j=\sqrt{\pi}\). The estimate is

I\approx e^{-3}\sqrt{\frac{2}{4}}\sqrt{\pi} =e^{-3}\sqrt{\frac{\pi}{2}} \approx 0.06240.

This is also the exact integral, since the integrand is Gaussian. The calculation confirms that the adaptive scale and Jacobian preserve the integral’s magnitude; the nodes alone would not do so.

Worked Example: A Quartic Residual Changes the Weights

Consider \(h(b)=-b^2/2-0.1b^4\). The mode is \(\hat b=0\): \(h'(b)=-b-0.4b^3\) vanishes there, and \(h(b)<h(0)=0\) for every \(b\ne0\). Also \(h''(0)=-1\), so \(H=1\). The local Gaussian scale is one, but the quartic term makes the integrand non-Gaussian away from the mode.

Use the two-node standard rule with \(x_1=-1/\sqrt{2}\), \(x_2=1/\sqrt{2}\), and \(w_1=w_2=\sqrt{\pi}/2\). The adaptive nodes are \(b_1=-1\) and \(b_2=1\). At either node, the log residual is \(h(b_j)-h(0)+x_j^2=(-1/2-0.1)+1/2=-0.1\). Thus each standard weight is multiplied by \(e^{-0.1}\). The resulting estimate is

I\approx \sqrt{2}\left(\frac{\sqrt{\pi}}{2}e^{-0.1} +\frac{\sqrt{\pi}}{2}e^{-0.1}\right) =\sqrt{2\pi}\,e^{-0.1} \approx 2.2681.

Without the residual adjustment, the same two nodes would give \(\sqrt{2\pi}\approx2.5066\), the local Gaussian integral. The quartic residual appropriately reduces the estimate at these nodes. This calculation does not establish how accurate a particular finite rule is for every non-Gaussian integrand; it shows exactly where the correction enters.

Worked Example: Shifted Mode and Narrower Scale

Let \(h(b)=-1.2-4.5(b-2)^2-0.1(b-2)^4\), a hypothetical log integrand in a dimensionless random-effect coordinate. It has its maximum at \(\hat b=2\). Differentiating twice gives \(h''(2)=-9\), so \(H=9\) and the local Gaussian scale is \(1/3\). The two standard nodes above map to \(b_j=2+\sqrt{2/9}\,x_j=2\pm\sqrt{2}/3\).

For either node, \((b_j-2)^2=2/9\) and \((b_j-2)^4=4/81\). The quadratic contribution cancels against \(x_j^2\) in the residual, leaving \(h(b_j)-h(2)+x_j^2=-0.1(4/81)=-0.00493827\). Each standard weight is multiplied by \(e^{-0.00493827}\). The two-node estimate is

I\approx e^{-1.2}\sqrt{\frac{2}{9}}\sqrt{\pi}\,e^{-0.00493827} =e^{-1.2}\sqrt{\frac{2\pi}{9}}\,e^{-0.00493827} \approx 0.25042.

If the quartic term were absent, the residual would equal one and the corresponding Gaussian integral would be \(e^{-1.2}\sqrt{2\pi/9}\approx0.25166\). Here the mode shift and narrower curvature determine where the nodes are placed, while the quartic term supplies a small additional weight correction.

Assumptions, Diagnostics, and Alternatives

The construction requires a finite target integral and a mode at which the log integrand can be evaluated and differentiated. For the curvature-based Gaussian scale, the mode must be interior and \(H>0\). These conditions make the local quadratic term a proper Gaussian kernel. A single-mode adaptive rule is most compelling when that mode captures the dominant mass and the integrand is reasonably represented by its local shape over the region covered by the nodes. A local maximum alone does not guarantee either condition.

Several checks help identify problems before trusting an approximation:

Check the mode and curvature. Verify that the fitted point is stationary and that the negative second derivative is positive. A zero or negative \(H\) makes the stated scale undefined or unsuitable.
Inspect the transformed nodes. Confirm that the nodes lie in a region where \(h(b_j)\) can be evaluated reliably and where the model’s integrand is meaningful. For bounded integration domains, nodes may extend outside the domain; the transformation or integration treatment must then respect that domain.
Examine residual log corrections. The values \(h(b_j)-h(\hat b)+x_j^2\) indicate how much the integrand differs from the local Gaussian at the sampled locations. Large or rapidly varying corrections warn that the local shape alone is not a good description over those nodes.
Look for competing modes or long tails. A mode-centered rule can miss substantial mass in another mode or in a tail that is poorly covered by the transformed nodes. A well-behaved local curvature does not rule out those possibilities.
Compute on a stable scale. When the integrand is extremely small or large, evaluate log contributions and combine positive terms with a stable log-sum calculation rather than exponentiating large intermediate values unnecessarily.

Laplace approximation uses the same mode and curvature but replaces the residual factor by one, giving \(I\approx \exp\{h(\hat b)\}\sqrt{2\pi/H}\). It is particularly transparent and exact for a Gaussian integrand, but it does not correct for non-Gaussian shape. Adaptive Gauss-Hermite quadrature retains the residual and can therefore represent more of that shape at the selected nodes. In contrast, nonadaptive Gauss-Hermite quadrature centered at zero uses a fixed location and scale; it may be inefficient when the integrand’s mass is displaced or concentrated at a different width. Each approach remains a numerical approximation except in special cases, and adaptation around one mode cannot by itself resolve a genuinely multimodal target.

Common pitfall: It is incorrect to move the nodes to \(\hat b+\sqrt{2/H}\,x_j\) and then use only the original weights and Jacobian. That treats the local Gaussian as if it were the actual integrand. The correction \(\exp\{h(b_j)-h(\hat b)+x_j^2\}\) is needed to integrate the original integrand rather than only its quadratic approximation.

Practice and Solutions

Check Your Understanding

Use the standard convention \(\int e^{-x^2}f(x)\,dx\approx\sum_jw_jf(x_j)\). The solutions are included after the questions.

For \(h(b)=-5-8(b-1)^2\), find the mode, \(H\), adaptive scale, and transformed node corresponding to a standard node \(x_j=1\).
Why is the condition \(H>0\) stronger than the statement that \(\hat b\) is a strict interior maximum?
For \(h(b)=-b^2/2-0.2b^4\), find the residual multiplier at the two nodes of the two-point rule used above.
State which term would be missing if one shifted and scaled the nodes but did not adjust the weights, and explain the consequence.

Solutions. (1) The mode is \(\hat b=1\), and \(h''(b)=-16\), so \(H=16\) and the local scale is \(1/\sqrt{16}=0.25\). The transformed node is \(1+\sqrt{2/16}(1)=1+\sqrt{2}/4\approx1.354\). (2) A strict maximum can have zero second derivative, as in \(h(b)=-b^4\) at zero. The scale \(H^{-1/2}\) requires a strictly positive \(H\), not merely \(H\geq0\). (3) The mode and curvature are again zero and one. At either node \(b_j=\pm1\), the residual exponent is \(-0.2b_j^4=-0.2\), so the multiplier is \(e^{-0.2}\approx0.8187\). (4) The missing term is the residual multiplier \(\exp\{h(b_j)-h(\hat b)+x_j^2\}\). Without it, the sum represents the local Gaussian factor rather than the original integrand at the transformed nodes.