GitHub

Gauss–Hermite cubature

Gauss–Hermite cubature approximates an expectation under a normal by a weighted sum of the function at fixed points. It suits expectations that have no closed form, in one or a few dimensions.

The method

For a scalar normal $N(m, v)$, the substitution $x = m + \sqrt{2v}\, ξ$ turns the expectation into an integral against $e^{-ξ^2}$, which the Gauss–Hermite rule with p nodes $ξ_i$ and weights $w_i$ evaluates:

\[\mathrm{E}[f(x)] = \frac{1}{\sqrt{π}} \int f\bigl(m + \sqrt{2v}\, ξ\bigr) e^{-ξ^2} \, \mathrm{d}ξ ≈ \frac{1}{\sqrt{π}} \sum_{i=1}^{p} w_i f\bigl(m + \sqrt{2v}\, ξ_i\bigr).\]

The nodes are the roots of the Hermite polynomial of degree p, and the sum is exact when f is a polynomial of degree up to 2p - 1. In d dimensions the points form a tensor-product grid of p^d points, mapped through a square root of the covariance, so the method stays practical for a handful of dimensions only.

getpoints and getweights return the mapped points and the normalised weights, $w_i / \sqrt{π}$, which sum to one. Three points integrate $x^2$ exactly, so they recover $\mathrm{E}[x^2] = m^2 + v$:

using MessagePassingRulesApproximations

gh = ghcubature(3)
m, v = 1.0, 4.0
points, weights = collect(getpoints(gh, m, v)), collect(getweights(gh, m, v))
sum(weights .* points .^ 2), m^2 + v
(5.0000000000000036, 5.0)

Moments of a reweighted normal

approximate_meancov returns the mean and variance of the density proportional to $g(x) N(x \mid m, v)$, for a non-negative $g$. A normal prior multiplied by a likelihood $g$ has this form, and its moments define the normal closest to the product. For a standard normal prior and a logistic likelihood of a positive label:

likelihood(x) = 1 / (1 + exp(-4x))
approximate_meancov(ghcubature(20), likelihood, 0.0, 1.0)
(0.7316388768995489, 0.46470455380917797)

The likelihood favours positive values, so the mean moves above zero and the variance shrinks below one.

API

MessagePassingRulesApproximations.GaussHermiteCubature — Type
GaussHermiteCubature(p::Int)

Gauss–Hermite cubature with p points per dimension, for expectations under a normal given by its mean and variance, or mean vector and covariance. It is exact for polynomials of degree up to 2p - 1 in each coordinate. In d dimensions it evaluates p^d points, a tensor product, so it suits a handful of dimensions only.

It provides getpoints and getweights, and through them approximate_meancov. ghcubature(p) builds the same.

Examples

julia> gh = GaussHermiteCubature(3);

julia> sum(w * x^2 for (w, x) in zip(getweights(gh, 1.0, 4.0), getpoints(gh, 1.0, 4.0))) ≈ 5.0  # E[x²] = m² + v
true
source
MessagePassingRulesApproximations.getpoints — Function
getpoints(method, mean::Real, variance::Real)
getpoints(method, mean::AbstractVector, covariance::AbstractMatrix)

The cubature points of method for a normal with these moments: a lazy generator of numbers, or of vectors in the multivariate case, in the order of getweights.

Warning

The multivariate generator yields one buffer, rewritten on every iteration, so that a sum over the points allocates nothing. collect therefore returns copies of the last point; use map(copy, getpoints(…)) to keep them. A consumer may overwrite the point it is given, as approximate_meancov does.

source
MessagePassingRulesApproximations.getweights — Function
getweights(method, mean::Real, variance::Real)
getweights(method, mean::AbstractVector, covariance::AbstractMatrix)

The cubature weights of method for a normal with these moments: a lazy generator of numbers summing to one, in the order of getpoints, so that E[f(x)] ≈ Σ wᵢ f(xᵢ).

source
MessagePassingRulesApproximations.approximate_meancov — Function
approximate_meancov(method, g, m::Real, v::Real) -> (mean, variance)
approximate_meancov(method, g, m::AbstractVector, P::AbstractMatrix) -> (mean, covariance)

The mean and variance (or covariance) of the density proportional to g(x) N(x | m, v), computed by the cubature method over the points of N(m, v). This is the projection of a normal reweighted by a likelihood onto a normal, the core of an expectation-propagation update.

Arguments

  • method: a cubature, such as GaussHermiteCubature(p); anything with getpoints and getweights;
  • g: a non-negative function of a point, a number for a scalar normal and a vector otherwise. It need not be normalised, but must not vanish at every point;
  • m, v or P: the mean and variance, or mean vector and covariance matrix, of the normal.

Examples

julia> m, v = approximate_meancov(ghcubature(21), x -> exp(-x^2 / 2), 0.0, 1.0);  # N(0, 1) × N(0, 1)

julia> isapprox(m, 0.0; atol = 1e-12) && v ≈ 0.5
true
source