[ANN] Quadriceps.jl: Gauss-Hermite and Gauss-Legendre accuracy in several dimensions, at a fraction of the nodes, with positive weights

Package Announcements
1

If you integrate against a multivariate normal or uniform in two to five dimensions, you are probably using one of three things, and each has a problem.

  • A product of one-dimensional Gauss rules. Exact, positive weights, but qᵈ nodes: a degree-9 rule in five dimensions costs 3125 evaluations, and that sits inside every likelihood evaluation, every optimization step, every bootstrap draw.
  • A sparse grid. Far fewer nodes, but some weights are negative. A rule with negative weights is not a probability distribution: a positive integrand can integrate to a negative number, a log-likelihood can become the log of something negative, and rounding error is amplified instead of averaged away.
  • Monte Carlo or quasi-Monte Carlo. Positive weights, any dimension, but the error falls slowly, and the integrand is evaluated at hundreds or thousands of points to get a few digits.

Quadriceps.jl offers a fourth: rules that integrate every polynomial up to a given degree exactly, like the product rule, whose weights are all positive, like the product rule, and which use a fraction of the nodes. The degree-9 rule in five dimensions has 244 nodes instead of 3125. For the cube at degree 21 in five dimensions, 10984 instead of 161051. They are the smallest positive-weight rules known to me; 116 of the 142 are new.

Because the weights are positive and sum to one, a rule is a discrete distribution, and anything you would do with a Gauss-Hermite grid you can do with it unchanged: mixed logit, random coefficients, Bayesian posteriors, expected utility, moments of a function of a normal vector.

The two functions follow gausshermite(q) and gausslegendre(q) from FastGaussQuadrature.jl, with the dimension in front. q is the size of the one-dimensional rule being replaced, so the degree is 2q − 1.

using Quadriceps

X, w = ghpos(3, 4)              # standard normal on ℝ³, degree 7: 27 nodes instead of 64
f(x) = x[1]^2 * x[2]^4
sum(w[i] * f(X[i, :]) for i in eachindex(w))        # E[Z₁² Z₂⁴] = 3.0

X, w = lepos(2, 5)              # uniform on [0,1]², degree 9: 17 nodes instead of 25
sum(w .* X[:, 1] .^ 3 .* X[:, 2] .^ 2)              # 1/4 · 1/3

ghpos(3; p = 7) == ghpos(3, 4)  # true: ask for a degree instead of a q

X is n × d, one node per row; w has n positive weights. The number type is an optional first argument: ghpos(Float128, 3, 4) gives the rule in quadruple precision and ghpos(BigFloat, 3, 4) the underlying 80-digit rule, with no extra dependency. Rules are stored up to degree 33 for the normal and 77 for the cube in two dimensions, with lower ceilings in higher dimensions; pragmatic = true serves anything beyond that as the cheapest tensor product of stored rules.

Install with Pkg.add("Quadriceps"); Julia 1.10 or later, and the only dependency is FastGaussQuadrature.jl.

The rules and how they were found are in Joris Pinkse, Positive weight Hermite and Legendre quadrature rules (2026), arXiv:2609.26840. The rules that descend from published ones are marked as such by Quadriceps.ruleinfo, so you can cite the source.

Requests for rules I have not stored are welcome.