Quadrature Methods

A crucial part of transport-map applications is selecting a suitable quadrature method. Quadrature is used in map optimization, in particular when evaluating the Kullback–Leibler divergence

Here, denotes the variable in the reference space with density , is the target density, and is the transport map parameterized by .

In practice we approximate the integral by a quadrature sum:

The quadrature points and weights must be chosen so the sum approximates expectations with respect to the reference measure .

Especially in Bayesian inference, where evaluating the target density can be expensive, using efficient quadrature methods is important to reduce the number of target evaluations [4].

Reference Density and Quadrature Choice

Since quadrature is performed in the reference space , the quadrature method must match said space and reference density. An overview is provided in the table below.

Quadrature rule	Reference	Notes
MonteCarloWeights	any (default: )	Uses random sampling from the default reference
LatinHypercubeWeights	any (default: )	Quasi-random sampling from the default reference
TensorProductWeights	any	Quadrature by full tensor-product of one-dimensional quadrature rule
GaussHermiteWeights		Tensor-product Gauss-Hermite quadrature
GaussLegendreWeights	or	Tensor-product Gauss-Legendre quadrature
SparseSmolyakWeights	any (default: )	Sparse grid quadrature rule; fewer points than full tensor product

Sampling-based Quadrature

Monte Carlo

Monte Carlo estimates an expectation by sampling from the reference density (default Normal() for MonteCarloWeights(n, d)) and forming the sample average. It is simple and robust; weights are uniform. Convergence is stochastic () but independent of dimension.

mc = MonteCarloWeights(500, 2)
scatter(mc.points[:, 1], mc.points[:, 2], ms=3,
    label="MC samples", title="Monte Carlo (500 pts)", aspect_ratio=1)

Latin Hypercube Sampling

Latin Hypercube is a stratified sampling design that improves space-filling over plain Monte Carlo; weights remain uniform. It often reduces variance for low to moderate dimensions.

lhs = LatinHypercubeWeights(500, 2)
scatter(lhs.points[:, 1], lhs.points[:, 2], ms=3,
    label="LHS samples", title="Latin Hypercube (500 pts)", aspect_ratio=1)

Tensor-Product Quadrature

Tensor-product quadrature constructs multidimensional rules by combining one-dimensional quadrature rules across each coordinate. Given a 1D rule with points, the tensor product forms all   combinations of points in dimension . The number of points per dimension is given as    with the level .

General tensor-product quadrature is performed with TensorProductWeights:

quad = TensorProductWeights(level, dim, knots)

Different knot choices are available: GaussHermiteKnots(), GaussLegendreKnots(), ClenshawCurtisKnots()

Tensor-product Gauss–Hermite

Gauss–Hermite quadrature is optimized for Gaussian reference measures on   and achieves spectral accuracy for smooth integrands.

hermite = GaussHermiteWeights(3, 2)
# alias: TensorProductWeights(3, 2, GaussHermiteKnots())
scatter(hermite.points[:, 1], hermite.points[:, 2],
    ms=4, label="Gauss–Hermite", title="Tensor Gauss–Hermite (level 3)", aspect_ratio=1)

Tensor-product Gauss-Legendre

Gauss–Legendre quadrature is optimized for bounded uniform reference measures on   (or after rescaling).

legendre = GaussLegendreWeights(3, 2)
# alias: TensorProductWeights(3, 2, GaussLegendreKnots())
scatter(legendre.points[:, 1], legendre.points[:, 2],
    ms=4, label="Gauss–Legendre", title="Tensor Gauss–Legendre (level 3) on [-1, 1]", aspect_ratio=1)

legendre_scaled = GaussLegendreWeights(3, 2, [0, 1])
scatter(legendre_scaled.points[:, 1], legendre_scaled.points[:, 2],
    ms=4, label="Gauss–Legendre", title="Tensor Gauss–Legendre (level 3) on [0, 1]", aspect_ratio=1)

Sparse-Grid Quadrature

Sparse-grid methods reduce the number of quadrature points compared to full tensor products while maintaining accuracy for smooth functions.

Sparse Smolyak grids use carefully chosen coefficients to cancel redundant high-dimensional interactions, reducing points from to approximately . Negative weights are expected and necessary for high-order correction terms.

Different knot choices (GaussHermiteKnots(), GaussLegendreKnots(), ClenshawCurtisKnots()) optimize sparse grids for different reference measures, similar to their tensor-product counterparts:

sparse = SparseSmolyakWeights(3, 2) # knots = GaussHermiteKnots()
scatter(sparse.points[:, 1], sparse.points[:, 2], ms=6,
    label="Smolyak", title="Sparse Smolyak GaussHermiteKnots", aspect_ratio=1)

sparse_legendre = SparseSmolyakWeights(3, 2, GaussLegendreKnots())
scatter(sparse_legendre.points[:, 1], sparse_legendre.points[:, 2], ms=6,
    label="Smolyak", title="Sparse Smolyak GaussLegendreKnots", aspect_ratio=1)

sparse_cc= SparseSmolyakWeights(3, 2, ClenshawCurtisKnots())
scatter(sparse_cc.points[:, 1], sparse_cc.points[:, 2], ms=6,
    label="Smolyak", title="Sparse Smolyak ClenshawCurtisKnots", aspect_ratio=1)

---

This page was generated using Literate.jl.