[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.

2

Nice! Is it conceivable that this construction will allow you to build nested multidimensional rules (ala Gauss–Kronrod in 1d) in order to have a “free-ish” error estimate?

3

Conceivable, yes. Easy, no.

First, Quadriceps.jl only delivers the rules: it does not compute them. The code to compute them is public at PositiveWeightQuadratureSolvers.jl but not in the registry to avoid confusion since computing them is not cheap.

As is, the best you can do is to take the tensor of a GK rule in a difficult dimension with one of the multivariate rules, but that too comes with caveats.