Tensor-Product Quadrature and Its Dimensional Cost

Why Tensor-Product Quadrature Becomes Expensive

A subject’s marginal likelihood in a nonlinear mixed-effects model integrates over all of that subject’s random effects. One-dimensional Gauss-Hermite quadrature provides a deterministic way to approximate an integral over a Gaussian random effect. With several effects, the direct extension is to take every combination of nodes from a one-dimensional rule. This construction is simple and useful, but its cost grows rapidly: using \(m\) nodes for each of \(q\) effects requires \(m^q\) evaluations for each subject.

The key trade-off is therefore between the accuracy gained by increasing the quadrature order and the cost of evaluating the integrand at all combinations of nodes. The number of evaluations is not merely proportional to the number of random effects. It is exponential in that dimension when the order is held fixed.

From a Gaussian Random-Effects Integral to a Product Rule

Let \(b_i\) be the \(q\)-dimensional random-effects vector for subject \(i\), with \(b_i\sim N(0,\Omega)\), where \(\Omega\) is positive definite. Let \(f_i(b_i;\theta)\) denote the subject-specific conditional likelihood contribution given the random effects and population parameters \(\theta\). The marginal likelihood contribution is the integral of that conditional contribution over the random-effects distribution.

Here \(\phi_q(\,\cdot\,;0,\Omega)\) is the multivariate normal density. To use a standard-normal quadrature rule, factor the covariance matrix as \(\Omega=LL^\mathsf{T}\), for example by a Cholesky decomposition, and set \(b_i=Lz\), where \(z\sim N(0,I_q)\). The transformed integral is

The coordinates of \(z\) are independent standard normal variables. This independence is a property of the whitened coordinates; the original random effects \(b_i\) may be correlated. The transformation incorporates both the covariance structure and the change of variables, so no additional determinant factor is needed in this standard-normal expectation form.

With this notation, the approximation to the subject contribution is

There are exactly \(m^q\) node tuples in this sum. The weights sum to one because \(\sum_{j=1}^{m}a_j=1\), so the rule integrates the constant function exactly. For \(N\) subjects, applying the same rule separately to each subject requires \(Nm^q\) evaluations of the subject-specific integrand, before accounting for other model-fitting work.

Exactness: What the Product Rule Guarantees

A one-dimensional \(m\)-node Gauss-Hermite rule for a standard normal expectation is exact for every polynomial of degree at most \(2m-1\). The tensor rule consequently has a precise multivariate exactness property: it integrates exactly any polynomial whose degree in each individual coordinate is at most \(2m-1\). This is a coordinate-by-coordinate bound; it is not a claim that the total degree across all coordinates may be as large as \(2m-1\) only.

To justify the multivariate claim, write a polynomial \(P(z_1,\ldots,z_q)\) satisfying these coordinatewise bounds as a finite sum of monomials. Each monomial has the form \(c_\alpha\prod_{k=1}^{q}z_k^{\alpha_k}\), where every \(\alpha_k\leq 2m-1\). Under the independent standard-normal density, its integral factors into a product of one-dimensional moments. The tensor rule for that monomial also factors into a product of one-dimensional quadrature sums. Each one-dimensional sum equals the corresponding moment by the one-dimensional exactness result. Thus the product rule gives the exact integral for every monomial, and by linearity it gives the exact integral for their sum \(P\). For a general nonlinear likelihood integrand, this property does not imply exact integration; the rule is then an approximation.

Counting the Computational Cost

For a fixed order \(m\), increasing the random-effects dimension from \(q\) to \(q+1\) multiplies the number of node evaluations by \(m\). For fixed dimension \(q\), increasing the order from \(m\) to \(m+1\) multiplies it by \(((m+1)/m)^q\). In the latter case even a one-node increase can be costly when \(q\) is large.

If different coordinates use different one-dimensional orders \(m_1,\ldots,m_q\), the full tensor grid contains \(\prod_{k=1}^{q}m_k\) nodes. This can be useful when the desired orders differ by coordinate, but it does not remove the multiplicative growth. Counts should be stated per subject and, when relevant, multiplied by the number of subjects.

Correlated Effects and a Worked Transformation

A tensor rule is applied in independent standard-normal coordinates even when the random effects are correlated. Before using the rule, the covariance transformation must be applied to every quadrature node. The following example also shows how a tensor sum can reproduce a known multivariate-normal expectation.

Implementation Checks, Alternatives, and Limitations

A correct node count does not guarantee a correct integral. In implementation, each node tuple must be evaluated after transforming it to the original random-effects scale, and its weight must be the product of the corresponding one-dimensional weights. The one-dimensional weights must use the same standard-normal convention as the transformed integral. If weights from a different Gauss-Hermite convention are used without the required rescaling, both the nodes and the integral can be wrong.

Useful basic checks include confirming that the one-dimensional weights sum to one, the product weights sum to one, the rule integrates a constant to one, and the covariance factor satisfies \(LL^\mathsf{T}=\Omega\). Check counts independently: with equal order \(m\), the number of tuples must be \(m^q\), not \(mq\). For a study with \(N\) subjects, multiply the per-subject count by \(N\). A discrepancy in any of these checks can create a numerical error that may be mistaken for a problem with the statistical model.

Tensor-product Gauss-Hermite quadrature is deterministic, and its polynomial exactness is a useful theoretical guarantee. Its main limitation is dimensional cost: even moderate \(m\) and \(q\) can make a full grid impractical. Adaptive quadrature around individual-effect modes can place nodes in a more relevant region of the integrand, though a full tensor grid still incurs multiplicative node growth. Importance sampling instead estimates the integral using weighted draws and avoids enumerating every combination, but introduces Monte Carlo error and requires attention to weight stability. These approaches address computational challenges differently; none makes it safe to assume that an approximation is accurate without checking it.

The tensor rule is most attractive when the number of random effects is modest and a deterministic integration rule is affordable. When the dimension is larger, the node-count calculation should be done before fitting, not after an unexpectedly expensive computation begins. The choice of method depends on the dimension, likelihood shape, computational budget, and the accuracy required for the estimation task.