Sparse Grid Quadrature for Moderate-Dimensional Effects

Why Combine Lower-Dimensional Rules?

A GLMM marginal likelihood requires integrating over random effects. If a cluster has a $d$-dimensional random-effects vector, even a moderate number of one-dimensional nodes can create a large tensor product: using $m$ nodes per dimension requires $m^d$ evaluations. Sparse-grid quadrature reduces this cost by combining selected lower-dimensional rules rather than evaluating every combination at the highest one-dimensional order.

The reduction is useful when the integrand is reasonably smooth and its important variation can be represented by low-order interactions among coordinates. It is not a promise of equal accuracy at a lower cost: sparse grids omit some high-order coordinate combinations, and their combination weights can be negative. The rule, the integrand, and the requested accuracy must therefore be considered together.

Definition: Let $Q_\ell$ denote a one-dimensional quadrature rule at level $\ell$, with $Q_0=0$. Its hierarchical increment is $\Delta_\ell=Q_\ell-Q_{\ell-1}$. For dimension $d$ and sparse-grid level $L\geq 0$, the Smolyak rule is $$ S_{d,L}=\sum_{\substack{i_1,\ldots,i_d\geq 1\\ (i_1-1)+\cdots+(i_d-1)\leq L}} \Delta_{i_1}\otimes\cdots\otimes\Delta_{i_d}. $$ Here a tensor product applies one one-dimensional rule in each coordinate. The level constraint includes the base rule and selected increments, while excluding combinations whose total level is too high.

In a GLMM, one common target is a cluster contribution of the form $$ L_j=\int p(y_j\mid b)\,p(b)\,db, $$ where $y_j$ denotes the cluster's observed responses and $b$ its random effects. After expressing the random-effects integral in suitable standard-normal coordinates $z$, write its integrand as $f_j(z)$. Sparse-grid quadrature approximates the expectation by a weighted sum $S_{d,L}(f_j)$. If the marginal likelihood is formed from several cluster contributions, each contribution must be approximated adequately before taking its logarithm.

How the Combination Works

The simplest nontrivial sparse-grid level has $L=1$. Its index set consists of the all-level-1 rule and terms that raise one coordinate to level 2. For two dimensions, expanding the increments gives $$ S_{2,1}=Q_2\otimes Q_1+Q_1\otimes Q_2-Q_1\otimes Q_1. $$ The subtraction corrects for the base rule counted in both coordinate-specific terms. In higher dimensions, the same inclusion-and-correction idea combines rules that refine a limited number of coordinates at a time.

Key idea: A sparse grid allocates evaluations to selected lower-dimensional interactions. Its efficiency is strongest when those interactions carry much of the integrand's variation; it can be weaker when important behavior depends on many coordinates acting jointly.

The node count depends on the one-dimensional rules and whether their nodes are nested. With nested rules, higher levels reuse earlier nodes, so increment terms can often be evaluated using only newly introduced points. With nonnested rules, points generally do not coincide, and the cost must be counted from the actual component rules and their shared points. The formal Smolyak sum alone does not guarantee a particular savings.

Worked Example: Building a Two-Dimensional Rule

Use standard-normal quadrature. Let $Q_1$ use the single node $0$ with weight 1, and let $Q_2$ be the three-point rule with nodes $0,\pm\sqrt{3}$ and weights $2/3,1/6,1/6$, respectively. Apply the level-$1$ sparse rule to $f(x,y)=x^2+y^2$.

The two-dimensional combination is $Q_2\otimes Q_1+Q_1\otimes Q_2-Q_1\otimes Q_1$. Because the one-dimensional rules share the node 0, the combined points are $(0,0)$, $(\sqrt{3},0)$, $(-\sqrt{3},0)$, $(0,\sqrt{3})$, and $(0,-\sqrt{3})$. The central point's combined weight is $2/3+2/3-1=1/3$; each of the other four points has weight $1/6$.

At the central point $f(0,0)=0$. At each of the four axis points, $f=3$. Thus $$ S_{2,1}(f)=\frac{1}{3}(0)+4\left(\frac{1}{6}\right)(3)=2. $$ The rule uses 5 distinct points. A full tensor product of the same three-point rule uses $3^2=9$ points. Here the sparse rule is exact for this function and saves four evaluations; this example does not imply that it is exact for general integrands.

Node Savings and Their Conditions

A useful way to understand the savings is to count hierarchical increments. Suppose a nested one-dimensional sequence has $1,3,5$ total nodes at its first three levels. The first level adds 1 node, and each of the next two levels adds 2 new nodes. For $d=4$ and $L=2$, the admissible increments are the base term, one level-2 increment in any one coordinate, one level-3 increment in any one coordinate, and pairs of level-2 increments in two coordinates.

Worked Example: Counting a Four-Dimensional Sparse Grid

For the nested sequence just described, count the new hierarchical points required at $d=4,L=2$, and compare with a full tensor product at level 3. The base term contributes $1^4=1$ point. There are four choices of coordinate for a level-2 increment, each contributing $2(1)(1)(1)=2$ new points, for $4(2)=8$. The four level-3 increments contribute another $4(2)=8$. Finally, there are $\binom{4}{2}=6$ pairs of level-2 increments; each contributes $2\cdot2=4$ new points, for $6(4)=24$.

The sparse-grid count is therefore $1+8+8+24=41$ new points. A full tensor product with 5 nodes per coordinate uses $5^4=625$ points. This count illustrates savings for a nested sequence; it is not a universal node count for Gaussian quadrature. In particular, nonnested rules may require more evaluations because their component grids share fewer points.

A lower count matters only if the resulting approximation is accurate enough. Sparse grids are often attractive when $d$ is moderate, the integrand is smooth, and a limited set of coordinate interactions dominates. Their performance can deteriorate when the integrand has sharp changes, pronounced tail behavior, or important high-order interactions. The relevant property is the integrand's structure in the chosen coordinates, not dimension alone.

Numerical Accuracy for a Smooth Integrand

A second calculation shows both how to use the combination and why it needs checking. Consider the hypothetical smooth integrand $f(x,y)=\exp(-0.2x^2-0.1y^2)$, integrated with respect to two independent standard-normal variables. This is a numerical illustration, not a fitted GLMM.

Worked Example: Sparse Approximation Versus a Tensor Product

Use the one- and three-point rules from the earlier calculation. The one-dimensional three-point approximation for $\exp(-aZ^2)$, where $Z\sim N(0,1)$, is $$ Q_2\{\exp(-aZ^2)\}=\frac{2}{3}+\frac{1}{3}\exp(-3a). $$ For $a=0.2$, this is $2/3+(1/3)\exp(-0.6)\approx 0.8496$. For $a=0.1$, it is $2/3+(1/3)\exp(-0.3)\approx 0.9136$. The one-point rule gives 1 for either factor.

The level-$1$ sparse approximation is therefore $$ S_{2,1}(f)=0.8496+0.9136-1=0.7632. $$ The full tensor product gives $0.8496(0.9136)\approx 0.7762$ using 9 points. The exact expectation is the product $1/\sqrt{1.4}\times1/\sqrt{1.2}\approx0.7715$. Thus the sparse result is about $0.0083$ below the exact value, while this tensor result is about $0.0047$ above it. Both approximations have error, and the sparse approximation is not automatically the more accurate one.

Assumptions, Diagnostics, and Alternatives

Before using a sparse rule for a random-effects integral, check that the integral and coordinate system match the quadrature rule. For example, the standard-normal rules above apply to an expectation with respect to independent standard-normal coordinates. A transformation must include any required density or Jacobian factors in the integrand. The integrand should also be sufficiently regular for the selected quadrature sequence; discontinuities or concentrated tail contributions can undermine its accuracy. Correlated random-effects coordinates require appropriate handling before applying rules that assume independent standard normals.

In computation, increase the sparse-grid level and compare the resulting cluster contributions or marginal log likelihood. Also compare fitted estimates across refinements when quadrature is being used during estimation. Large changes signal that the current level may be inadequate. Inspect numerical signs as well: sparse-grid combination weights can be negative, and a weighted sum for a positive integrand can consequently be nonpositive. A nonpositive approximation cannot be used as a valid positive marginal likelihood contribution.

Approach	Main advantage	Main limitation
Tensor-product quadrature	Direct, systematic refinement across all coordinates	Node count grows exponentially with dimension
Sparse-grid quadrature	Can reduce evaluations when lower-order interactions dominate	May miss important high-order structure; weights may be negative
Laplace approximation	A local approximation can be much cheaper than many quadrature evaluations	Accuracy may suffer when the integrand is poorly represented by its local Gaussian approximation

These methods answer a computational trade-off, not a single universal ranking. Tensor products are a useful reference when dimension is small and reliable refinement is affordable. Sparse grids are an option when a full tensor product is too costly and the interaction structure is favorable. The Laplace approximation is another alternative, with its own local-approximation limitations. If a sparse-grid fit changes materially with refinement, a higher level, a tensor check in a smaller problem, or a different integration strategy may be needed.

Common pitfall: “The sparse grid has fewer nodes, so it must be at least as accurate as the tensor rule.” Fewer nodes establish a computational saving, not an accuracy guarantee. Accuracy depends on the integrand, the coordinate system, the rule levels, and the interactions omitted by the sparse construction. Check convergence rather than inferring it from node count.

When Sparse Grids Are Useful

Sparse grids are most compelling when moderate dimension makes a full tensor product expensive, but the integrand remains smooth enough that selectively combining lower-dimensional rules captures its important behavior. They are less attractive when accuracy depends on interactions spanning many random-effect coordinates, when quadrature weights create numerical instability, or when added levels fail to stabilize the likelihood. In a GLMM, approximation quality matters because errors in cluster integrals feed into the marginal likelihood used for estimation. Node savings are valuable only when the likelihood approximation remains stable for the intended analysis.

Check Your Understanding

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

For $d=2,L=1$, write the Smolyak rule in terms of $Q_1$ and $Q_2$.
Why can nested rules reduce the number of new evaluations compared with nonnested rules?
In the worked smooth-integrand calculation, which approximation was closer to the exact expectation: the sparse or tensor result?
Name two diagnostics that can reveal an inadequate sparse-grid approximation.

Solutions. (1) $S_{2,1}=Q_2\otimes Q_1+Q_1\otimes Q_2-Q_1\otimes Q_1$. (2) Nested rules reuse points from earlier levels, so shared function evaluations need not be repeated; nonnested rules generally share fewer points. (3) The tensor result, approximately $0.7762$, differs from the exact $0.7715$ by about $0.0047$, less than the sparse result's difference of about $0.0083$. (4) Examples include comparing likelihoods across successively higher sparse-grid levels and checking whether a positive-integrand approximation is nonpositive; changes in fitted estimates across refinements are another warning.