There is an amazing line of work on implementing exact and efficient integrators for multivariate polynomials over simplices. This integration is useful in probabilistic programming because the inference task in the continuous setting boils down to integration.

Formally, let \(\Delta\) be a simplex, a subset of \(\mathbb{R}^n\) of the form \(\{ x \in \mathbb{R}^n \mid Ax \le b\}\) for some set of linear constraints \(A\) and constants \(b\). The canonical simplex of dimension \(n\) is the simplex of the form \(\{\mathbf{x} \in R^n \mid \mathbf{e}^T \mathbf{x} \le 1\}\), where \(\mathbf{e} = (1, \cdots, 1) \in \mathbb{R}^n\).

Let \(f \in R[\mathbf{x}]\) be a polynomial over variables \(\mathbf{x}\). Then, our goal is to solve integrals of the form:

$$\int_\Delta f(\mathbf{x})~\mathrm{d}\mathbf{x}.$$

A naive approach to solving this integral to use iterated integration and Fourier-Motzkin elimination. This is exponential in the number of dimensions of the simplex.

In general, we would really like to avoid iterated integration, because there are issues with error compounding.

Lasserre Integration

Lasserre, Jean B. "Simple formula for integration of polynomials on a simplex." BIT Numerical Mathematics 61.2 (2021): 523-533.

Incredibly, there is a way of integrating any polynomial that is efficient (in the degree of the polynomial and simplex, for natural notions of degree), easy-to-implement, and does not involve solving any integrals at all. I do not yet appreciate the proof of why it works, but I can spell out the details of what you need to do to implement it.

This is a weird object, best illustrated with example. Suppose \(f(x, y) = 3x^3 + y + 2xy^2\). Then, \(\hat{f}(x,y)=(3!)3x^3 + y + (2!)2xy^2\).

As further notation, we will need to refer to the homogeneous forms of a polynomial. For a polynomial \(f\), let \(f_j\) denote the sum of all terms of degree \(j\). For instance, for the above polynomial, \(f_3(x, y) = 3x^3 + 2xy^2\), and \(f_1(x,y) = y\).

Prepare to have your mind blown. Lasserre integration lets us compute an integral \(\int_\Delta f(\mathbf{x})\mathrm{d}\mathbf{x}\) by evaluating the Bombieri form in a particular way on a finite number of evaluation points:

OK, now we are ready to state the main theorem of the paper that lets us integrate polynomials:

There is a straightforward way of generalizing the above to arbitrary simplexes using linear transformations, but let's see an example of it in action on our running example. For \(f(x, y) = 3x^3 + y + 2xy^2\) and \(\hat{f}(x,y)=(3!)3x^3 + y + (2!)2xy^2\), we have that:

$$\hat{f_0}(x,y) = 0 \quad \hat{f}_1(x,y) = y \quad \hat{f}_2(x,y) = 0 \quad \hat{f}_3(x,y) = 18x^3 + 4xy^2$$

Then, we can compute each evaluation point \(\xi_i\):

• \(\xi_1 = \frac{\mathbf{e}}{(2+1)^1} = \frac{\mathbf{e}}{3}\)
• \(\xi_2 = \frac{\mathbf{e}}{((2+1)(2+2))^{(1/2)}} = \frac{\mathbf{e}}{\sqrt{12}}\)
• \(\xi_3 = \frac{\mathbf{e}}{((2+1)(2+2)(2+3))^{(1/3)}} = \frac{\mathbf{e}}{120^{(1/3)}}\)

Now, integrating the polynomial is as simple as evaluating each of the terms above on the correct evaluation points:

• \(\hat{f}_1(\xi_1) = 1/3.\)
• \(\hat{f}_3(\xi_3) = \frac{18}{120} + \frac{4}{120}\)

Then, by the above theorem we have that \(\int_\Delta f(\mathbf{x}) \mathrm{d}\mathbf{x} = \frac{1}{2}(0 + \frac{1}{3} + \frac{18}{120} + \frac{4}{120}) = \frac{7}{20}\). You can check that this is in fact the correct answer to this integral.

Arbitrary simplex

So far we've only discussed how to integrate over the canonical simplex. What about a general simplex given as a system of constraints \(\{x \mid Ax \le b, x \ge 0\}\)? This is fairly easy: the transformation \(A^{-1}(y - b)\) maps all points in \(\Delta\) to the standard simplex, so we can this as a change of variables. Let \(\Delta_e\) be the standard simplex. Let \(g(y) = f(A^{-1}(y - b))\); note that \(g\) is a polynomial in the same degree as \(f\). Then, by the usual change of variables formula, we have (eq. 2.11 in the paper):

$$\int_\Delta f(x)dx = \frac{1}{\mathrm{det}(A^{-1})} \int_{\Delta_e} g(y)~dy.$$

The multivariate case is a straightforward generalization of this.

Resources

There is quite a long line of work on this topic, and the LattE page has a lot of good places to start. To my knowledge, the method presented here of Lasserre integration is significantly easier to implement, and looks like it would be much more effective in practice than existing methods.

This blog post has an implementation of Lasserre integration in several languages. For posterity here is the Julia implementation (which is totally undocumented, but appears to compute correct answers) test inline code omg:

using TypedPolynomials
using LinearAlgebra
function integratePolynomialOnSimplex(P, S)
    gens = variables(P)
    n = length(gens)
    v = S[end]
    B = Array{Float64}(undef, n, 0)
    for i in 1:n
        B = hcat(B, S[i] - v)
    end
    Q = P(gens => v + B * vec(gens))
    s = 0.0
    for t in terms(Q)
        coef = TypedPolynomials.coefficient(t)
        powers = TypedPolynomials.exponents(t)
        j = sum(powers)
        if j == 0
            s = s + coef
            continue
        end
        coef = coef * prod(factorial.(powers))
        s = s + coef / prod((n+1):(n+j))
    end
    return abs(LinearAlgebra.det(B)) / factorial(n) * s
end

To use this code:

using TypedPolynomials
@polyvar x y z
P = x^4 + y + 2*x*y^2 - 3*z

# simplex vertices
v1 = [1.0, 1.0, 1.0]
v2 = [2.0, 2.0, 3.0]
v3 = [3.0, 4.0, 5.0]

v4 = [3.0, 2.0, 1.0]
# simplex
S = [v1, v2, v3, v4]

# integrate
integratePolynomialOnSimplex(P, S)