Tensor-Product Quadrature and Its Computational Cost

Why the Product Structure Matters

A subject-specific likelihood in a joint model is obtained by integrating over the subject’s latent random effects. With several random effects, a direct numerical rule must account for every coordinate of that vector. Tensor-product Gauss–Hermite quadrature does this by combining a one-dimensional rule across all coordinates. The construction is straightforward, but its cost can rise quickly: using \(m\) nodes in each of \(q\) dimensions requires \(m^q\) evaluations of the integrand for one subject.

The earlier material on marginal likelihoods and Gauss–Hermite quadrature motivates approximating these integrals, while the earlier material on Cholesky transformation handles correlated random effects. Here the focus is the product construction and its computational consequences.

Definition and Construction

Definition: Let \(q\) be the number of random effects, and suppose their covariance matrix \(D\) is positive definite. Write \(D=LL^{\mathsf T}\), where \(L\) is a Cholesky factor, and express the random-effects vector as \(b=Lz\), with \(z\sim N(0,I_q)\). Let \(z_1,\ldots,z_m\) and \(a_1,\ldots,a_m\) be nodes and weights for an \(m\)-node Gauss–Hermite rule for a standard normal expectation, with \(\sum_{r=1}^m a_r=1\). The tensor-product rule approximates an expectation of \(f(Lz)\) by summing over all \(q\)-tuples of univariate nodes, multiplying the weights across coordinates. It uses \(m^q\) node tuples.

For subject \(i\), let \(g_i(b;\theta)\) denote the conditional joint-data likelihood contribution given random effects \(b\), at parameter value \(\theta\). The marginal contribution is the Gaussian expectation of this conditional likelihood. After the transformation \(b=Lz\), the tensor approximation is

L_i(\theta) = \int_{\mathbb R^q} g_i(b;\theta)\,\phi_q(b;0,D)\,db = E\{g_i(Lz;\theta)\} \ \approx\ \sum_{r_1=1}^m\cdots\sum_{r_q=1}^m \left(\prod_{j=1}^q a_{r_j}\right) g_i\!\left(L \begin{bmatrix} z_{r_1}\\ \vdots\\ z_{r_q} \end{bmatrix};\theta\right).

The coordinates of \(z\) are independent standard normal variables. That independence is why the product weights are appropriate, even when the original random effects are correlated: their correlation is represented by \(L\). This construction assumes a valid positive-definite covariance matrix, a correctly specified one-dimensional standard-normal rule, and an integrable conditional likelihood under the Gaussian distribution. It approximates the integral; the formula does not itself guarantee a particular level of accuracy.

In practice, a rule is often written as a set of multi-indices \((r_1,\ldots,r_q)\), with each index ranging from 1 to \(m\). For each tuple, the implementation forms the corresponding standard-normal vector, transforms it by \(L\), evaluates the conditional likelihood, and multiplies by the tuple’s product weight. If an adaptive rule changes the node locations or scaling around a subject-specific posterior region, the same Cartesian-product counting applies when it still uses \(m\) nodes in each coordinate.

Worked Example: Forming a Two-Dimensional Product Rule

Consider a hypothetical two-dimensional standard-normal expectation and use the three-node rule with nodes \(-\sqrt{3},0,\sqrt{3}\) and respective weights \(1/6,2/3,1/6\). There are \(3^2=9\) node pairs. To check the construction, approximate the expectation of the nonnegative function \(f(z_1,z_2)=1+z_1^2+2z_2^2\).

The product-rule sum is

\sum_{r=1}^3\sum_{s=1}^3 a_r a_s \left(1+z_r^2+2z_s^2\right).

The univariate weights sum to \(1/6+2/3+1/6=1\), and their weighted second moment is \[ \frac{1}{6}(\sqrt{3})^2+\frac{2}{3}(0)^2+\frac{1}{6}(\sqrt{3})^2 =\frac{3}{6}+\frac{3}{6}=1. \] Therefore the nine-term sum equals \(1(1)(1)+1(1)+2(1)(1)=4\). The exact expectation is also \(1+E(z_1^2)+2E(z_2^2)=1+1+2=4\). This check confirms the arithmetic and illustrates how product weights sum over a grid; the nine evaluations are still the computational units, even though the sum can be simplified for this particular function.

Exactness and Its Justification

An \(m\)-node Gauss–Hermite rule for a standard normal expectation integrates any univariate polynomial of degree at most \(2m-1\) exactly. The tensor rule inherits a coordinatewise version of this property. It exactly integrates a multivariate polynomial if, in every coordinate separately, the polynomial’s degree is no greater than \(2m-1\). This is a statement about the rule’s exactness for polynomials, not a general claim that joint-model likelihood integrands are polynomials.

To justify the result, first consider a monomial \(z_1^{d_1}\cdots z_q^{d_q}\) with \(d_j\leq 2m-1\) for every \(j\). Independence of the standard-normal coordinates gives

E\left(\prod_{j=1}^q z_j^{d_j}\right) = \prod_{j=1}^q E(z_j^{d_j}).

The product-rule sum for the same monomial also separates:

\sum_{r_1=1}^m\cdots\sum_{r_q=1}^m \prod_{j=1}^q \left(a_{r_j}z_{r_j}^{d_j}\right) = \prod_{j=1}^q \left(\sum_{r=1}^m a_r z_r^{d_j}\right).

Each univariate sum on the right equals its corresponding normal moment by the one-dimensional exactness property. Thus the product sum equals the multivariate expectation for this monomial. Linearity extends the result to any polynomial satisfying the stated coordinatewise degree limits. The polynomial’s total degree can exceed \(2m-1\); for example, a product of several coordinate powers may still be integrated exactly if each individual coordinate degree remains within the limit.

For a general conditional likelihood, exactness does not follow from this argument. The integrand may be nonpolynomial, sharply concentrated, or vary substantially over the region important to the integral. Increasing \(m\) may improve an approximation, but the number of terms grows at the same time.

Worked Example: Comparing Orders as Dimension Increases

Suppose a hypothetical joint model has three random effects: a longitudinal intercept, a longitudinal slope, and a shared effect entering the hazard. Compare a five-node rule with a nine-node rule in each coordinate. The five-node tensor rule requires \(5^3=125\) evaluations per subject. The nine-node rule requires \(9^3=729\).

The increase is \(729-125=604\) evaluations per subject, and the ratio is \[ \frac{729}{125}=\left(\frac{9}{5}\right)^3=5.832. \] Thus, holding the cost of one conditional-likelihood evaluation fixed, the nine-node rule uses about 5.832 times as many evaluations per subject. This calculation describes computational work, not how much more accurate the nine-node approximation will be.

The dimension also matters even when the order stays fixed. At five nodes per coordinate, a two-dimensional rule has \(5^2=25\) tuples, while a four-dimensional rule has \(5^4=625\). Adding two random effects multiplies the node count by \(5^2=25\). The general count is \(m^q\), not \(mq\).

From Node Counts to Likelihood Cost

If there are \(N\) subjects and the same \(m\)-by-\(q\) product rule is applied to each subject, one full evaluation of all subject contributions requires \(N m^q\) conditional-likelihood evaluations. If evaluating subject \(i\)’s conditional likelihood at a node has cost \(C_i\), the work is proportional to \(m^q\sum_{i=1}^N C_i\). The cost of forming the transformed random-effects vectors and summing weighted contributions is additional, but the nodewise likelihood evaluations are often the dominant component.

This count applies to one evaluation of the sample likelihood at a fixed parameter value. An iterative optimizer or other estimation procedure will evaluate the likelihood repeatedly, so the total work also depends on the number of evaluations required. In joint models, conditional contributions may involve several longitudinal measurements and survival terms; subjects with more expensive contributions can therefore require more time per node even though the number of nodes is the same.

Worked Example: Cost for a Sample of Subjects

Suppose a hypothetical analysis has \(N=180\) subjects and four random effects. Compare nine and eleven nodes per coordinate. With nine nodes, the rule has \(9^4=6{,}561\) tuples per subject, giving \(180(6{,}561)=1{,}180{,}980\) conditional-likelihood evaluations for one full sample-likelihood evaluation. With eleven nodes, it has \(11^4=14{,}641\) tuples per subject, giving \(180(14{,}641)=2{,}635{,}380\) evaluations.

The difference is \(2{,}635{,}380-1{,}180{,}980=1{,}454{,}400\) evaluations. The ratio is \[ \frac{14{,}641}{6{,}561} =\left(\frac{11}{9}\right)^4 \approx 2.231. \] If, purely for illustration, one conditional evaluation takes \(0.0004\) seconds on the computing setup in use, the two totals correspond to \(1{,}180{,}980(0.0004)=472.392\) seconds and \(2{,}635{,}380(0.0004)=1{,}054.152\) seconds, respectively. These are hypothetical timings for one likelihood evaluation, not predictions for a complete model fit.

Implementation Checks, Alternatives, and Limitations

Count the coordinates.
Identify the dimension \(q\) of the random-effects vector being integrated for a subject.

Choose a univariate rule.
Use \(m\) standard-normal nodes and weights in each coordinate, with weights normalized to sum to one.

Generate every tuple.
Enumerate all \(m^q\) multi-indices, multiply their coordinate weights, and transform each node vector using \(L\).

Evaluate and sum.
Compute the conditional subject contribution at each transformed vector and sum the weighted values.

Check behavior as the rule changes.
Compare approximations at more than one order and inspect whether the resulting likelihood calculations are numerically stable.

Useful implementation checks include verifying that the one-dimensional weights sum to one, that the full product weights sum to one, and that the rule gives the known expectation for simple test functions such as \(1\) and \(z_j^2\). For a diagonal covariance matrix, the Cholesky transformation should not introduce cross-coordinate mixing. More generally, confirm that the transformed nodes use the intended factor \(L\), since using a factor with the wrong orientation changes the covariance represented by the nodes.

A common mistake is to reason that “four random effects and eleven nodes” requires \(4(11)=44\) evaluations per subject. The rule is a Cartesian product: every choice in one coordinate is combined with every choice in the others. The correct count is \(11^4=14{,}641\), not 44. This distinction is consequential for both runtime planning and memory use if all tuples are stored at once.

Tensor-product quadrature is systematic and convenient at low dimension, and it allows direct refinement by changing the univariate order. Adaptive Gauss–Hermite quadrature can place nodes more effectively for a concentrated subject-specific integrand, but if it still uses \(m\) nodes in each of \(q\) dimensions, it retains the \(m^q\) count. Monte Carlo or quasi-Monte Carlo integration can avoid a full Cartesian grid and can be attractive when dimension is high, but introduces sampling or sequence-based approximation error and its own choices about replication and stability. Sparse-grid quadrature is another approach that aims to reduce the cost of a full product construction in higher dimensions; its construction is a distinct method.

The main limitation is therefore structural: even modest increases in \(q\) can make the full grid impractical. A smaller grid may be computationally manageable but inaccurate if the integrand is concentrated or has features that the chosen nodes do not resolve. Conversely, increasing \(m\) is not automatically a practical remedy because it raises cost exponentially in \(q\). Compare successive likelihood approximations, check for stable computations, and interpret changes in light of the model’s fitting needs rather than treating a large node count as proof of accuracy.

Key takeaway: A tensor-product rule with \(m\) nodes per coordinate and \(q\) random effects has exactly \(m^q\) node tuples per subject. Its polynomial exactness follows coordinate by coordinate from one-dimensional Gauss–Hermite exactness, but accuracy for a general likelihood integrand must be assessed numerically.

Check Your Understanding

Answer each question, then compare with the solutions below.

How many tuples are in a tensor rule with six nodes per coordinate for three random effects?
A two-dimensional rule uses standard-normal weights summing to one in each coordinate. What is the sum of all product weights?
For \(m=4\), what is the coordinatewise polynomial degree limit for exactness? Does the result impose a total-degree limit of seven?
A sample has 75 subjects, two random effects, and eight nodes per coordinate. How many conditional-likelihood evaluations are used for one full likelihood evaluation?

Solutions. (1) The count is \(6^3=216\) tuples. (2) The product weights sum to \((\sum_r a_r)(\sum_s a_s)=1(1)=1\). (3) Each coordinate degree can be at most \(2m-1=7\); there is no requirement that total degree be at most seven, provided every coordinate separately satisfies the limit. (4) There are \(8^2=64\) tuples per subject, so the total is \(75(64)=4{,}800\) conditional-likelihood evaluations.