Why the Quadrature Rule Moves During Optimization
In a generalized linear mixed model, integrating over random effects produces a marginal likelihood that is generally not available in closed form. Quadrature approximates that integral using a finite set of nodes. In adaptive quadrature, the nodes are positioned and scaled using the shape of the integrand near its conditional mode. That shape depends on the model parameters. Consequently, when an outer optimizer proposes new parameter values, the centers and scales used by the integration rule should be recomputed for those values.
This is not merely a way to choose better initial nodes. The mode and curvature define the particular approximation being evaluated at a parameter value. Reusing their old values after changing parameters evaluates a different, potentially poor quadrature approximation. The essential update is: find the new conditional mode, evaluate the new local curvature there, and construct the quadrature nodes from both.
The mode is conditional on the observed cluster data and on the current value of \(\theta\); it is not generally the same as the random-effects mean, which is often zero. When \(d=1\), write \(H_\theta=-h_{bb}(m_\theta;\theta)>0\) and \(s_\theta=H_\theta^{-1/2}\). A large curvature means a narrow integrand and a smaller scale; a small curvature means a wider integrand and a larger scale.
How the Center and Scale Define the Rule
For clarity, use a Gauss–Hermite rule whose nodes \(x_q\) and weights \(w_q\) approximate integrals weighted by \(\exp(-x^2)\). In one dimension, an order-\(Q\) adaptive rule uses locations \(b_q=m_\theta+\sqrt{2}s_\theta x_q\). The factor \(\sqrt{2}\) is part of the change of variables for this convention. The corresponding integral approximation is
The added \(x_q^2\) compensates for the \(\exp(-x^2)\) weight built into the Gauss–Hermite rule. The scale factor is the Jacobian of the transformation from \(x\) to \(b\). Both the node locations and the multiplier depend on the current mode and curvature.
For \(d\) random effects, let \(L_\theta L_\theta^\mathsf{T}=H_\theta^{-1}\), where \(L_\theta\) is a Cholesky factor. A tensor rule transforms each vector node \(x\) to \(b=m_\theta+\sqrt{2}L_\theta x\); its Jacobian factor is \(2^{d/2}|L_\theta|\). The dimension and node-set choice determine the computational burden, but do not change the update principle: the center and local scale come from the current integrand.
Why the Mode Must Be Recomputed
At an interior conditional mode, the score with respect to the random effects is zero: \(\nabla_b h(m_\theta;\theta)=0\). If \(\theta\) changes, this equation changes too. Under differentiability and nonsingular curvature, the implicit-function rule shows how the mode responds. For a scalar parameter \(\theta_j\), differentiate the zero-score equation:
Thus, except in special cases, even a fixed-effect update moves the conditional mode. The curvature also changes: it is the negative Hessian evaluated at the new mode and parameter value. Its change can reflect both the direct effect of changing \(\theta\) and the indirect effect of evaluating the Hessian at a different \(b\). A scale computed before the parameter update is therefore not generally the correct local scale afterward.
The outer optimization method supplies a new candidate \(\theta\).
At the candidate parameters, solve \(\nabla_bh(b;\theta)=0\) and verify that the solution is a local maximum.
Evaluate \(H_\theta=-\nabla_b^2h(m_\theta;\theta)\), then obtain the univariate scale or multivariate factor \(L_\theta\).
Transform the rule’s nodes using the new center and scale, and include the appropriate Jacobian and weight adjustment.
This procedure assumes that \(h\) is sufficiently smooth near the mode, that the relevant maximum is interior, and that its curvature is positive definite. In practice, a mode-finding algorithm may use the preceding iteration’s mode as a starting value; that is a useful computational shortcut. It must still converge to the mode for the new parameters before the new quadrature rule is evaluated.
Worked Example: Updating a Random-Intercept Center After a Fixed-Effect Step
Consider a hypothetical cluster with two binary responses, one success and one failure. Let \(b\sim N(0,1)\), and use a logit model with common fixed intercept \(\beta\). Ignoring constants that do not depend on \(b\), the cluster log integrand is \(h(b;\beta)=\beta+b-2\log\{1+\exp(\beta+b)\}-b^2/2\). Its first derivative is \(h_b=1-2p-b\), where \(p=\operatorname{logit}^{-1}(\beta+b)\), and its curvature is \(-h_{bb}=1+2p(1-p)\). The latter is positive for every finite \(b\), so the log integrand is strictly concave and its stationary point is the unique mode.
At \(\beta=0\), symmetry gives \(m_0=0\) and \(p=0.5\). Therefore \(H_0=1+2(0.5)(0.5)=1.5\), and \(s_0=1/\sqrt{1.5}\approx0.8165\). With the three-point physicists’ Gauss–Hermite nodes \(x=-\sqrt{3/2},0,\sqrt{3/2}\), the transformed locations are \(m_0+\sqrt{2}s_0x\), or approximately \(-1.4142,0,1.4142\).
Now suppose the outer optimizer proposes \(\beta=0.4\). The new mode solves \(1-2\operatorname{logit}^{-1}(0.4+b)-b=0\). Numerically, \(m_{0.4}\approx-0.1328\); at this value \(p\approx\operatorname{logit}^{-1}(0.2672)\approx0.5664\). Substitution gives \(H_{0.4}=1+2(0.5664)(0.4336)\approx1.4912\), so \(s_{0.4}=1/\sqrt{1.4912}\approx0.8189\).
The new outer-node displacement is \(\sqrt{3}s_{0.4}\approx1.4184\), giving updated locations approximately \(-1.5512,-0.1328,1.2856\). These are not the old locations shifted by an arbitrary optimizer step: the center comes from the new score equation, and the scale comes from curvature at that new center. For an order-three rule, the integral is then approximated by the stated weighted sum using these updated locations and \(h(b;0.4)\).
Worked Example: A Variance Update Changes the Scale Even When the Center Does Not Move
Return to the same hypothetical cluster, set \(\beta=0\), and let \(b\sim N(0,\tau^2)\). For finding the mode and curvature, use the log-kernel \(h(b;0,\tau)=b-2\log\{1+\exp(b)\}-b^2/(2\tau^2)\). The normalized random-effects density adds \(-\log\tau-\tfrac12\log(2\pi)\) to this log-kernel, and this term must be included when evaluating the marginal likelihood. At \(b=0\), the likelihood score is \(1-2(0.5)=0\), and the Gaussian prior score is also zero. Strict concavity therefore gives \(m_\tau=0\) for any \(\tau>0\).
Although the center is unchanged, the prior curvature changes. At the mode, \(H_\tau=1/\tau^2+2(0.5)(0.5)=1/\tau^2+0.5\). For \(\tau=1\), \(H_1=1.5\), \(s_1\approx0.8165\), and the outer three-point nodes are approximately \(\pm\sqrt{3}(0.8165)=\pm1.4142\), together with zero.
If the proposed standard deviation is \(\tau=2\), then \(H_2=1/4+0.5=0.75\), \(s_2=1/\sqrt{0.75}\approx1.1547\), and the outer nodes become \(\pm\sqrt{3}(1.1547)\approx\pm2.0000\), again with zero between them. This example isolates the scale update: a correct recomputation can change the quadrature locations even when the mode happens to remain fixed.
Numerical Checks and Failure Modes
The update is only as reliable as the mode and curvature calculations. At each candidate parameter value, check that the random-effects score at the computed center is small, that the Hessian is negative definite for \(h\) (equivalently, \(H_\theta\) is positive definite), and that the mode solver has converged. In one dimension, check \(h_b(m_\theta;\theta)\approx0\) and \(H_\theta>0\). In multiple dimensions, inspect the smallest eigenvalue of \(H_\theta\) or the success of its Cholesky factorization. A nearly singular curvature matrix implies an extremely broad local scale in at least one direction and can make the transformed rule numerically unstable.
Other warning signs include a mode that jumps substantially between nearby parameter values, a solver that reaches its iteration limit, or a quadrature objective that changes abruptly when the optimizer makes a small step. A jump can reflect a failed mode solve, but it can also indicate multiple local modes or a change in which mode is dominant. A local Gaussian scale around one mode does not by itself resolve a genuinely multimodal integrand. In such cases, do not interpret a failed factorization or erratic objective as ordinary optimizer noise; investigate the conditional integrand and the mode search.
A practical implementation may warm-start the mode search from the previous outer iteration, but should use the converged current mode, not the starting value, to construct the rule. Likewise, if numerical derivatives of the approximated objective are used, each perturbed parameter evaluation should recompute the mode, curvature, and nodes. Holding those quantities fixed computes derivatives of a frozen-node expression rather than of the adaptive approximation being optimized.
Alternatives and Limits
| Approach | Center and scale behavior | Main trade-off |
|---|---|---|
| Adaptive quadrature | Recomputed from the current mode and curvature | Tracks local integrand shape, but requires reliable mode and Hessian calculations at each parameter value |
| Fixed Gauss–Hermite quadrature | Uses a predetermined center and scale, commonly based on the random-effects distribution | Simpler setup, but can place too few nodes where the likelihood is concentrated |
| Laplace approximation | Uses the mode and curvature in a local second-order approximation | Efficient, but replaces the integral by a local approximation rather than summing a chosen quadrature rule |
| Simulation-based integration | Uses random draws rather than a deterministic mode-centered node set | Can handle complex integrands, but introduces simulation variability and its own computational demands |
Adaptive quadrature is useful when a deterministic approximation should concentrate nodes near the conditional integrand’s important region. It is not a guarantee of accuracy: the adequacy of a particular node order depends on the integrand and the rule, and a local mode and curvature may poorly describe a broad or multimodal shape. The key point here is narrower: whatever quadrature order is selected, its adaptive center and scale belong to the current parameter values. They must be updated as the outer optimizer moves.
Common Mistake
Wrong reasoning: “The optimizer changed only \(\beta\) a little, so the old center and scale are close enough; I can keep the same quadrature nodes throughout optimization.” This confuses a possible warm start with the final center used in the integration. Even a modest parameter change alters the mode equation, and the scale can change through both the likelihood and the random-effects distribution.
Correction: Use old quantities, at most, to initialize the new mode search. Then solve the current score equation, compute curvature at the new mode, and transform the quadrature nodes using that new scale. A rule held fixed across parameter updates may be a deliberate fixed-quadrature method, but it is not the adaptive rule described here.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- In the definition of adaptive quadrature, why is the scale based on \(H_\theta^{-1}\) rather than \(H_\theta\)?
- For the first worked example, what score equation must the updated mode satisfy when \(\beta=0.4\)?
- In the variance-update example, why do the nodes change even though the conditional mode remains zero?
- What numerical condition should hold for a multivariate curvature matrix before using its Cholesky factor to scale nodes?
- Is it valid to use the preceding iteration’s mode as a starting value? Explain the distinction between that use and reusing the preceding adaptive center.
Solutions. (1) \(H_\theta\) measures local curvature: greater curvature means a narrower peak. The local covariance-like scale is the inverse curvature \(H_\theta^{-1}\), and in one dimension its standard-deviation scale is \(H_\theta^{-1/2}\). (2) The equation is \(1-2\operatorname{logit}^{-1}(0.4+b)-b=0\), whose solution is approximately \(-0.1328\). (3) The prior contribution to curvature changes from \(1/\tau^2=1\) at \(\tau=1\) to \(1/4\) at \(\tau=2\), changing the total curvature from \(1.5\) to \(0.75\) and the scale from about \(0.8165\) to \(1.1547\). (4) \(H_\theta\) must be positive definite so that its inverse exists and a valid Cholesky factor of \(H_\theta^{-1}\) can be computed. (5) Yes. The old mode may initialize a search, but the search must converge to a mode at the new parameters; that newly computed mode and curvature define the adaptive center and scale.