Sparse-Grid Quadrature for Many Random Effects

Why Sparse Grids Are Used

A tensor-product quadrature rule assigns a one-dimensional rule to each random-effect coordinate and evaluates the integrand at every combination of those nodes. This is straightforward, but the number of evaluations grows rapidly with the number of random effects. Smolyak sparse-grid quadrature reduces that cost by combining a selected set of tensor-product rules. The combination uses both positive and negative coefficients, so it is not simply a smaller tensor grid with positive weights.

For an individual in an NLME model, write the random effects in standard-normal coordinates as \(b=\mu+Lz\), where \(z\) has independent standard-normal components and \(LL^\mathsf{T}\) is the random-effects covariance matrix. If \(g(z)\) is the individual’s conditional likelihood at the corresponding random effects, the marginal likelihood contribution is the expectation \(E\{g(Z)\}\). Quadrature approximates this expectation. The same integration idea applies if the integrand includes other factors, provided the chosen one-dimensional rules are constructed for the relevant measure or the measure has been transformed appropriately.

Definition: A Smolyak sparse-grid rule is a weighted combination of selected tensor-product quadrature rules at different one-dimensional levels. Its level selection limits how many coordinates can simultaneously use high-level rules. The combination coefficients can be negative; after coincident nodes are combined, individual node weights can also be negative.

Here, \(d\) is the number of integration dimensions, \(Q_\ell\) is a one-dimensional quadrature rule at level \(\ell\), and a multi-index \(j=(j_1,\ldots,j_d)\) has positive integer entries. Write \(|j|=j_1+\cdots+j_d\). The level \(q\) controls which tensor-product rules are included. The convention below uses \(q\geq d\), so the smallest possible multi-index has total level \(d\).

The Combination Formula and Its Coefficients

Let \(Q_j=Q_{j_1}\otimes\cdots\otimes Q_{j_d}\) denote the tensor-product rule at multi-index \(j\). The Smolyak combination formula is

A_{d,q} = \sum_{\substack{j_k\geq 1\\q-d+1\leq |j|\leq q}} (-1)^{q-|j|} \binom{d-1}{q-|j|}\,Q_j.

Only indices with \(0\leq q-|j|\leq d-1\) occur. Thus the coefficient depends on how far the index’s total level is below \(q\), not on the individual coordinate levels. The signed coefficient can be derived by counting the allowable coordinate-level choices with their signs. For a fixed \(j\), put \(s=q-|j|\). The contributions correspond to choosing \(r\) coordinates to shift up one level; there are \(\binom{d}{r}\) such choices and their sign is \((-1)^r\). The total-level constraint permits \(r=0,\ldots,s\), when \(0\leq s\leq d-1\). Therefore the coefficient is the signed sum

\sum_{r=0}^{s}(-1)^r\binom{d}{r} = (-1)^s\binom{d-1}{s}.

For completeness, the identity follows by substituting Pascal’s rule, \(\binom{d}{r}=\binom{d-1}{r}+\binom{d-1}{r-1}\). In the sum, terms from the second part cancel the preceding terms from the first part; the only uncancelled term is \((-1)^s\binom{d-1}{s}\). This is a signed sum, not the number of allowable choices. For example, with \(d=2\), \(q=3\), and \(j=(1,1)\), the choices are \(i=(1,1),(2,1),(1,2)\): one choice has sign \(+1\), and two have sign \(-1\), giving \(1-2=-1\). That equals \((-1)^1\binom{1}{1}\).

Worked Example: Expanding a Two-Dimensional Rule

Take two standard-normal Gauss-Hermite levels: \(Q_1 f=f(0)\) and \(Q_2 f=\{f(-1)+f(1)\}/2\). For \(d=2\) and \(q=3\), the allowed indices are \((1,1),(2,1),(1,2)\). Their coefficients are \(-1,+1,+1\), respectively, so

A_{2,3}=Q_{2}\otimes Q_{1}+Q_{1}\otimes Q_{2}-Q_{1}\otimes Q_{1}.

For \(g(x,y)=e^{x+y}\), \(Q_1(e^x)=1\) and \(Q_2(e^x)=\{e^{-1}+e\}/2=\cosh(1)\approx1.543\). Thus

A_{2,3}g = \cosh(1)\cdot 1+1\cdot\cosh(1)-1\cdot1 = 2(1.543)-1 = 2.086.

The exact standard-normal expectation is \(E\{e^{Z_1+Z_2}\}=e^{(1+1)/2}=e\approx2.718\). The sparse rule is an approximation, not an exact answer in this case. The calculation also shows why a low-level sparse rule can be materially inaccurate for a rapidly changing integrand.

The combination form is useful because each \(Q_j\) is a familiar tensor rule. The same construction can also be understood through one-dimensional hierarchical differences, \(Q_\ell-Q_{\ell-1}\), with \(Q_0=0\). That view helps explain why nested rules can share many nodes. It does not change the signed coefficients in the combination formula.

Nodes, Weights, and Computational Cost

Suppose a one-dimensional rule \(Q_\ell\) has nodes \(x_{\ell,a}\) and weights \(w_{\ell,a}\). A tensor rule \(Q_j\) evaluates the integrand at every \(d\)-tuple of its one-dimensional nodes, with product weights. In a Smolyak combination, each tensor rule’s weights are multiplied by its coefficient in the formula. If a node occurs in more than one tensor rule, its contributions are added. The resulting sparse-grid weight at that node can be negative, even if every original one-dimensional weight is positive.

Worked Example: A Signed Weight at the Centre

Use the two-dimensional rule from the previous example. \(Q_2\otimes Q_1\) places weight \(1/2\) at \((-1,0)\) and \((1,0)\); \(Q_1\otimes Q_2\) places weight \(1/2\) at \((0,-1)\) and \((0,1)\). The term \(-Q_1\otimes Q_1\) places weight \(-1\) at \((0,0)\). These five nodes are distinct, so the combined weights are four weights of \(1/2\) and one weight of \(-1\).

For the constant integrand \(g(x,y)=1\), the rule returns

\frac12+\frac12+\frac12+\frac12-1=1.

The negative centre weight is part of the combination, not an error. It helps make the rule exact for constants while cancelling contributions that would otherwise be counted repeatedly. For a nonconstant integrand, however, negative weights can make cancellation numerically important.

Worked Example: Counting Nodes in Five Dimensions

Consider \(d=5\), \(q=7\), and nested one-dimensional rules with \(m_1=1\), \(m_2=3\), and \(m_3=5\) nodes at their first three levels. The sparse combination includes indices with \(3\leq |j|\leq7\). Since every one of the five coordinates has level at least 1, the smallest possible total is 5; the largest level any one coordinate can reach is \(7-4=3\).

With nested rules, the new nodes introduced at levels 1, 2, and 3 number \(1,2,2\), respectively. Counting the distinct nodes by their hierarchical level indices gives the following totals. At total level 5, all coordinates are at level 1, contributing \(1\) node. At total level 6, one coordinate is at level 2, giving \(5\cdot2=10\) nodes. At total level 7, either one coordinate is at level 3, giving \(5\cdot2=10\), or two are at level 2, giving \(\binom{5}{2}\cdot2^2=40\). Hence the sparse grid has \(1+10+10+40=61\) distinct nodes.

The full tensor rule using level 3 in every dimension has \(5^5=3{,}125\) nodes. This comparison illustrates the potential computational saving when the relevant one-dimensional rules are nested. It does not establish equal accuracy: the full tensor rule resolves all coordinates at level 3 simultaneously, whereas this sparse grid does not.

Worked Example: A Three-Dimensional Combination

For \(d=3\) and \(q=4\), the indices with total level 4 have coefficient \(+1\), and the index \((1,1,1)\), with total level 3, has coefficient \(-\binom{2}{1}=-2\). The three total-level-4 indices are \((2,1,1),(1,2,1),(1,1,2)\). Therefore

A_{3,4} = Q_2\otimes Q_1\otimes Q_1 + Q_1\otimes Q_2\otimes Q_1 + Q_1\otimes Q_1\otimes Q_2 - 2Q_1\otimes Q_1\otimes Q_1.

Using the same \(Q_1,Q_2\) as above with \(g(x,y,z)=e^{x+y+z}\), each of the first three terms is \(\cosh(1)\), and the last tensor rule returns 1. The approximation is

A_{3,4}g = 3(1.543)-2(1) = 2.629.

The exact expectation is \(e^{3/2}\approx4.482\). The example verifies both the \(-2\) coefficient and the practical limitation: a sparse rule’s low node count alone does not ensure adequate integration accuracy.

When the Cost Reduction Is Useful

Sparse grids are most attractive when the integrand is sufficiently smooth in each coordinate and its important variation can be represented through relatively low-order interactions among coordinates. In a random-effects integral, this can occur when many effects contribute smoothly but the likelihood is not dominated by a narrow, irregular feature requiring high resolution in every coordinate at once. Nested one-dimensional rules are especially useful because the same nodes can be reused across levels. For non-nested rules, tensor terms may share few or no nodes, so the number of distinct evaluations can be much larger than the nested-grid count.

The comparison with tensor-product quadrature is therefore conditional. A tensor product is simple and can be highly effective in a small number of dimensions, particularly when high accuracy is needed across all coordinates. A sparse grid limits simultaneous high levels and can sharply reduce the number of evaluations as \(d\) grows, but it may be less effective when the integrand has strong high-order interactions or sharp localized changes. Monte Carlo or importance sampling can be alternatives when dimension is high or the integrand is difficult to resolve deterministically; those methods have different sources of error and cost, so node counts alone do not decide between them.

For an NLME likelihood, quadrature is often repeated for many individuals and at many parameter values during estimation. Reducing evaluations per integral can therefore save substantial computation. Yet the signed combination means that a computed marginal likelihood is a sum with cancellation. In severe cases, rounding error can be large relative to the final sum, or a numerical approximation can even become nonpositive although the exact likelihood is positive. Stable summation and checking that the resulting likelihood is finite and positive are basic safeguards; a suspicious result should not be treated as a valid log-likelihood.

Key takeaway: Smolyak quadrature replaces one large tensor product with a signed combination of selected tensor rules. Its cost advantage depends on the level structure, node reuse, dimension, and integrand; it is not a universal accuracy improvement or a guarantee of fewer evaluations.

Common Mistakes and Limitations

Counting choices instead of summing their signs. For \(d=2,q=3,j=(1,1)\), there are three allowable choices, but the coefficient is \(1-2=-1\), not 3.
Assuming all sparse-grid weights are positive. Combination coefficients can be negative, and node weights after combining terms can be negative.
Claiming equal accuracy from fewer nodes. Sparse and full tensor rules allocate resolution differently. Their error depends on the integrand, levels, and one-dimensional rules.
Ignoring nesting and cancellation. Non-nested rules can require many distinct evaluations, while signed weights can amplify numerical error when large positive and negative contributions nearly cancel.

Sparse-grid quadrature is a deterministic integration strategy, not a substitute for examining whether its assumptions suit the likelihood being integrated. Smoothness and coordinate structure support its rationale; irregularity, strong interactions, and unstable cancellation limit it. A smaller evaluation count is valuable only if the resulting integral is sufficiently reliable for the estimation task.

Check Your Understanding

Work through these questions, then compare your reasoning with the solutions.

For \(d=4\), \(q=6\), and \(|j|=5\), calculate the coefficient of \(Q_j\).
In the two-dimensional \(q=3\) rule, why is the coefficient on \(Q_1\otimes Q_1\) negative even though the rule returns 1 for a constant integrand?
For the five-dimensional nested example, calculate the number of newly introduced nodes at total hierarchical level 7.
Name two conditions that make a sparse-grid cost reduction more plausible, and one reason it may fail to improve accuracy.

Solutions. (1) Here \(s=q-|j|=1\), so the coefficient is \((-1)^1\binom{3}{1}=-3\). (2) The coefficient corrects overlap among the tensor rules: the four positive weights sum to 2, and the central weight \(-1\) makes the constant integral equal to \(2-1=1\). (3) A level-3 coordinate contributes \(5\cdot2=10\) nodes; two level-2 coordinates contribute \(\binom{5}{2}\cdot2^2=40\). Together there are \(10+40=50\) new nodes at that total level. (4) Nested one-dimensional rules and a sufficiently smooth integrand with limited high-order coordinate interactions support savings. Strong interactions or sharp localized changes can make the sparse rule inaccurate at its chosen level.