Gauss-Hermite Integration: Examples in Python

In Gaussian filtering, we are interested in approximating Gaussian integrals of the form:

n-dimensional Gaussian Distribution times g(x)

Let’s break down this nasty looking integral.

Here, g(x) is any function of x (which can be a Multidimensional Vector). In the following examples, we will explore different values for this function. N(x | m,P) represents a Normal (Gaussian) distribution (which also can be multidimensional) with mean vector ‘m’ and covariance matrix ‘P’.

Other examples of this type of Integral include ‘First Moment’ (mean), ‘Second Moment’ (variance or Covariance Matrix), ‘Entropy’ of the Gaussian distribution. I will be covering these examples here.

Now, we already know the closed form solutions of the above mentioned examples. But for any arbitrary g(x), for example in case of filtering (or smoothing), this integral becomes arbitrarily complex and we need to rely on approximations. We can use Monte-Carlo estimates of the integrals, which I will do for some examples. However, with increasing dimensionality of the integral (which is the case in filtering), Monte-Carlo estimates can be computationally expensive. Especially when approximating this integral is just one step of many in an already complex algorithm like Gauss-Hermite Kalman Filter (GHKF) or Unscented Kalman Filter (UKF).

We start with a simple case. A 1D Standard Normal Gaussian Distribution (mean = 0, sd = 1). We can write the approximated integral as:

Gauss-Hermite Approximation

Wi are the ‘weights’ and xi’s are the so called ‘sigma-points’ at which we will evaluate the function values. It is essentially a weighted sum of the function at some special points (sigma points). In Gauss–Hermite integration, as in all Gaussian quadrature, the weights and sigma points are chosen such that with a polynomial integrand, the approximation becomes exact. It turns out that the polynomial order with a given number of points is maximized if we choose the sigma points to be roots of Hermite polynomials. The ‘p-th’ order Hermite Polynomial is given below:

p-th order Hermite Polynomial

The weights are given as:

i-th weight for the corresponding Hermite polynomial

Before you start getting nervous about coding it all up, I will relieve your pain by saying that there is a neat little Python library called NumPy, which you might have heard of, that comes to our rescue. We will get to that in a minute. But first, a caveat. We above expression is for a standard Normal distribution. But we want to generalise it for any Normal Distribution (with mean ‘μ’ and standard deviation of ‘σ’. There is a simple “reparameterization trick” that comes in handy here. We can scale a standard normal distribution by σ and shift it by μ to get a our desired Normal distribution

The Reparameterization Trick

With this, we just need to ‘reparameterize’ our sigma points according to the mean and standard deviation. The rest of the algorithm will remain the same. Now, for the examples, I have taken g(x) = x, which should give us the first moment (mean) of the Normal distribution. In the second example, g(x) = (x-μ)², which should give us the second moment or variance of the Normal Distribution.

import numpy as np
from numpy.polynomial.hermite_e import hermegauss

Getting the weights and sigma points (for a standard normal distribution) of a 3rd order Hermite polynomial is as easy as:

unit_sigma_points, hermite_weights = hermegauss(3) #3rd order
hermite_weights /= np.sum(hermite_weights) #need to scale the weights 

mu = 3  #mean
sd = 2  #standard deviation
s = unit_sigma_points*sd + mu #reparameterised sigma points

Now, g(x) = x. We should get the mean:

mu_est = np.sum(hermite_weights*s) #estimated mean ( = 3)

g(x) = (x-μ)². We should get the variance:

var_est = (np.sum(hermite_weights*(s-mu)**2) #estimated variance
sd_est = np.sqrt(var_est)  #estimated sd = 2

Suppose g(x) = -log(N(x | μ, σ), which is the Entropy of Normal distribution. Although a closed form expression for entropy does exist for a Normal Distribution, we will use TensorFlow Probability’s distribution.entropy() method.

import tensorflow as tf
import tensorflow_probability as tfp
dis = tfd.Normal(3,2)  #Normal Distribution
true_int = dis.entropy() # = -2.11
log_prob = -dis.log_prob(s) #calculates log_Normal @ sigma points
est_int = np.dot(log_prob, hermite_weights) # calculates weighted sum. 
#est_int = -2.11 (= true value)

Here I have looked at 1D Gaussians and a few examples. I have some more examples, albiet a bit messy, here at my GitHub:

[embed]Projects/Bayesian/Numerical Integration (MC, Gauss Quadrature).ipynb at main · sob-ANN/Projects Project Portfolio . Contribute to sob-ANN/Projects development by creating an account on GitHub.github.com

Hopefully, there is a Part 2, which looks at Multidimensional examples and some more complex 2-Dimensional Integrals. If you’re interested in that, here is a another messy code which does just that:

[embed]Projects/Bayesian/Probabilistic Kalman Filter.ipynb at main · sob-ANN/Projects Project Portfolio . Contribute to sob-ANN/Projects development by creating an account on GitHub.github.com

Note: The Jupyter Notebook linked contains Monte-Carlo estimates of the integrals as well.

Thanks for making it to the end.

Reference:

1. Simo Särkkä and Lennart Svensson (2023). Bayesian Filtering and Smoothing. Second Edition. Cambridge University Press.