Tutorials › Generalized Linear Mixed Models (GLMM) › Why Tensor Quadrature Explodes with Dimension

Laplace Approximation and Adaptive Quadrature · Tutorial 186 of 698

Why Tensor Quadrature Explodes with Dimension

Learn to calculate tensor-product quadrature workloads, estimate their growth as dimension or order increases, and judge when the resulting node counts may be impractical.

Advanced 11 min read

What You'll Learn

  • Calculate the number of tensor-product nodes from the per-coordinate quadrature orders.
  • Explain why equal-order tensor rules grow exponentially with random-effects dimension.
  • Compare node counts when quadrature order changes by one.
  • Translate a per-integral node count into a total evaluation workload across clusters.
  • Use a stated computational budget to identify feasible integer quadrature orders.
  • Distinguish node-count feasibility from quadrature accuracy and overall runtime.

Why Count the Nodes?

Tensor-product quadrature extends a one-dimensional rule to an integral over several random-effects coordinates by taking every combination of the one-dimensional nodes. That construction is systematic, but it can require many evaluations even when the number of coordinates is only moderate. Before choosing an order for numerical integration, it is therefore useful to estimate the workload implied by that choice.

The key issue is the interaction between two quantities: the number of nodes used in each coordinate and the number of random-effects coordinates being integrated. Increasing either quantity raises the total node count; increasing both can make the computation grow quickly. This count is a planning tool, not a guarantee of accuracy or a complete prediction of runtime.

Definition: Suppose a tensor-product rule uses \(q_j\) one-dimensional nodes in coordinate \(j\), for \(j=1,\ldots,d\). Its full tensor product has $$ N_{\mathrm{tensor}}=\prod_{j=1}^{d}q_j $$ multidimensional node combinations. If every coordinate uses the same order \(q\), then $$ N_{\mathrm{tensor}}=q^d. $$ Here, \(d\) is the dimension of the random-effects vector being integrated for one likelihood contribution, and \(q\) is the number of one-dimensional nodes per coordinate.

For a cluster \(i\), a generic marginal likelihood contribution has the form

$$ L_i=\int_{\mathbb{R}^d} L_i(\mathbf b_i) f(\mathbf b_i)\,d\mathbf b_i, $$

where \(\mathbf b_i\) is the \(d\)-dimensional random-effects vector, \(L_i(\mathbf b_i)\) is the conditional likelihood contribution, and \(f(\mathbf b_i)\) is its random-effects density. A tensor rule approximates this integral by evaluating its integrand at a Cartesian product of quadrature nodes and combining the evaluations with product weights. The count \(q^d\) is the number of such combinations when each coordinate uses \(q\) nodes.

This count assumes a full, unpruned tensor product: every combination is included, and each combination is evaluated once for the integral. It applies whether the one-dimensional nodes are nonadaptive or are shifted and scaled by an adaptive construction. Adaptation can change where nodes lie and what weights they receive; it does not, by itself, change the number of combinations. The calculation also assumes a common order \(q\). If coordinate \(j\) uses a different order \(q_j\), use the product \(\prod_j q_j\) instead.

The Growth Mechanism

The power \(q^d\) gives two useful comparisons. Holding dimension fixed, increasing order from \(q\) to \(q+1\) changes the node count by the factor \(((q+1)/q)^d\). Holding order fixed, adding one coordinate multiplies the count by \(q\). These are exact ratios, not approximations.

$$ \frac{(q+1)^d}{q^d}=\left(1+\frac{1}{q}\right)^d, \qquad \frac{q^{d+1}}{q^d}=q. $$

Thus, the cost of increasing the order by one becomes more pronounced in higher dimensions. Likewise, one additional random-effect coordinate multiplies the work by the full one-dimensional order. This is exponential growth in dimension when \(q>1\): each added coordinate multiplies the preceding workload by \(q\), rather than adding a fixed number of nodes.

Worked Example: Comparing Two and Six Random-Effects Coordinates

Suppose a hypothetical calculation uses a seven-node rule in every random-effects coordinate. Compare the node counts for a two-dimensional and a six-dimensional integral.

For \(d=2\), the count is $$ N_{\mathrm{tensor}}=7^2=7\times 7=49. $$ So one cluster’s integral requires 49 integrand evaluations under the full tensor rule.

For \(d=6\), the count is $$ N_{\mathrm{tensor}}=7^6=7\times7\times7\times7\times7\times7=117{,}649. $$ The six-dimensional rule therefore uses \(117{,}649/49=2{,}401\) times as many node combinations as the two-dimensional rule. This ratio also follows directly from adding four coordinates: \(7^4=2{,}401\).

The order is identical in both calculations, but the six-coordinate integral is much more expensive under this construction. These are counts, not evidence that either rule is accurate enough; accuracy must be assessed separately.

Using a Computational Budget

A budget can be expressed as a maximum number \(B\) of integrand evaluations per integral. Under a common order \(q\), feasibility requires

$$ q^d\leq B. $$

Because \(q\) must be a positive integer, the largest order that satisfies this bound is the greatest integer not exceeding \(B^{1/d}\). The root provides a boundary for planning; the integer order must still be checked by raising it to the \(d\)th power.

Worked Example: Finding the Largest Feasible Order

A hypothetical software run is allocated at most 40,000 integrand evaluations for each cluster’s integral. The random-effects vector has dimension \(d=6\), and a full tensor rule with common order \(q\) is planned. Find the largest feasible integer order.

The budget condition is $$ q^6\leq 40{,}000, \qquad q\leq (40{,}000)^{1/6}\approx 5.848. $$ Since the order is an integer, \(q\) cannot exceed 5. Check both the feasible order and the next order: $$ 5^6=15{,}625\leq40{,}000, \qquad 6^6=46{,}656>40{,}000. $$ Therefore, order 5 is the largest feasible integer order under this per-integral budget. The remaining budget is \(40{,}000-15{,}625=24{,}375\) evaluations, but it is not enough to accommodate order 6. The result says only that order 5 meets the stated node-count limit; it does not establish that order 5 gives adequate numerical accuracy.

The budget must be defined carefully. A per-integral limit differs from a limit on the full model calculation. If a model has \(m\) cluster-specific integrals and each uses the same tensor rule, the basic number of node-integral evaluations across clusters is \(m q^d\). This multiplication gives a useful first workload estimate, though it does not include all costs of fitting the model.

Worked Example: Scaling the Work Across Clusters

Consider a hypothetical GLMM with 120 clusters. For each cluster, suppose the random-effects vector has dimension 4 and an order-9 full tensor rule is used. Estimate the number of cluster-specific integrand evaluations in one pass through all clusters.

The node count per cluster is $$ 9^4=9\times9\times9\times9=6{,}561. $$ Across 120 clusters, the estimated total is $$ 120\times6{,}561=787{,}320 $$ cluster-specific evaluations in one pass.

If the order were raised to 11 with the same dimension and cluster count, the per-cluster count would be $$ 11^4=14{,}641, $$ and the total would be $$ 120\times14{,}641=1{,}756{,}920. $$ The ratio is \(14{,}641/6{,}561\approx2.231\), which equals \((11/9)^4\). Thus, increasing the order from 9 to 11 more than doubles the node evaluations in this four-dimensional example. The comparison assumes that every cluster uses the stated order and that every node entails one relevant integrand evaluation.

What the Count Does—and Does Not—Tell You

The node count is valuable because it exposes the structural cost of a full tensor product before a lengthy computation is attempted. It also helps distinguish the effect of dimension from the effect of sample size. Dimension changes the number of combinations per integral; the number of clusters multiplies that per-integral workload when cluster contributions are evaluated separately.

However, \(m q^d\) is not generally a complete runtime estimate. Evaluating the integrand at one node can have different costs across models. For example, the computation may involve many observations within a cluster, and the cost of each evaluation can depend on cluster size and model structure. In an estimation procedure, the integrations may also be repeated at many parameter values. Derivatives, optimization steps, and implementation details can add further work. Consequently, report the count as a node-evaluation estimate, not as an exact elapsed time.

Nor does a large node count prove that tensor quadrature is unusable. Feasibility depends on available computing resources, implementation efficiency, the number of repeated evaluations, and the accuracy required. Conversely, a count that fits a budget does not prove that the chosen order is sufficient. The numerical convergence checks introduced in the earlier material on quadrature order remain essential: compare results at increasing orders where computation allows, and assess whether the quantities relevant to the analysis are stable enough for the intended use.

Alternatives and Failure Modes

The main contrast is with a Laplace approximation, which approximates an integral using local information around a mode rather than evaluating a full Cartesian grid. It does not incur a \(q^d\) tensor-node count, although it has its own approximation error and computational requirements. Tensor quadrature instead offers a direct numerical integration strategy whose workload is explicit from the orders and dimension. Neither method is automatically preferable in every model.

Sparse-grid quadrature is another response to the cost of full tensor products; it is designed to reduce the number of nodes in appropriate moderate-dimensional settings. Its construction and accuracy considerations are distinct from the full tensor count derived here. A node-count calculation can signal when it is worth considering a different integration strategy, but does not determine which alternative will be accurate for a particular integrand.

Several practical checks help detect a planning or implementation problem:

  • Verify the dimension. Count the coordinates being integrated for each contribution, not the number of clusters or the total number of observations.
  • Verify the rule’s order convention. In this tutorial, \(q\) means the number of one-dimensional nodes per coordinate. A software setting described as “order” should be checked to confirm that it has this meaning.
  • Check the product independently. For a common order, compare a direct calculation of \(q^d\) with the result from repeated multiplication; for unequal orders, multiply the coordinate-specific counts.
  • Keep per-integral and total counts separate. Multiplying by the number of clusters is appropriate for a pass over separate cluster contributions, not for the node count of one integral.
  • Do not infer accuracy from affordability. A feasible order may still be too low for a stable numerical approximation; an expensive order may be unnecessary if lower orders already give sufficiently stable results.

A common mistake is to say that adding one random-effect coordinate adds \(q\) nodes. It does not: it multiplies the existing count by \(q\). For example, at order 7, moving from two to three coordinates changes the count from \(7^2=49\) to \(7^3=343\), an increase of 294 nodes and a multiplication by 7. Another mistake is to treat \(q^d\) as a direct accuracy measure. It is only a count of node combinations; accuracy depends on how well the quadrature rule represents the particular integrand.

Key takeaway: A full tensor rule with \(q\) nodes in each of \(d\) coordinates requires \(q^d\) node combinations per integral. Use this exponential growth to screen computational feasibility, then assess numerical accuracy separately.

Check Your Understanding

Use the full tensor-product rule described above. The solutions follow the questions so that each node-count calculation can be checked.

  1. How many node combinations are required for an order-4 rule in dimension 5?
  2. At order 8, by what factor does the node count change when dimension increases from 3 to 4?
  3. A per-integral budget is 1,000 evaluations in dimension 3. What is the largest feasible integer order?
  4. For 60 clusters, dimension 2, and order 10, calculate the total cluster-specific node evaluations in one pass. State one reason this may not equal the full fitting runtime.

Solutions. (1) \(4^5=1{,}024\) node combinations. (2) Adding one coordinate multiplies the count by \(8\), since \(8^4/8^3=8\). (3) The condition is \(q^3\leq1{,}000\). Since \(10^3=1{,}000\), the largest feasible integer order is 10. (4) Each cluster uses \(10^2=100\) nodes, so \(60\times100=6{,}000\) cluster-specific evaluations in one pass. Full fitting can require repeated evaluations during parameter optimization, and the cost per node can also vary with the model and cluster data.

Clinical Biostats