Integrating a barycentric monomial over a simplex

Question

Are there standard formulas for the integral over a simplex of a monomial in the barycentric coordinates? Can someone supply a reference? I think I have seen such formulas, but I am unable to find such material by searching on the web. This doesn't seem to be a difficult problem, but maybe there is a neat and memorable way of presenting the answer.

Answer 1

Not for the first time in my life, and almost certainly not the last, I have to eat humble pie. I thought that Igor had missed the point in his comment. In fact, it was I who missed the point. My answer has the advantage that it uses only very elementary methods. It needed a strong nudge and explanations from Igor before I realized the relevance of his different approach. What is great about Igor's answer is that it links a rather mundane question of mine to some very beautiful mathematics that was previously unfamiliar to me and probably to others. I have reversed my previous downvote (which would in any case have been uncalled for, even if his answer had been wrong) into an upvote. My own answer follows. See Igor's answer in his comments.

Thanks to Tim Goodman whose remarks eased my task of finding (a proof of) the formula. The following material is presumably all standard, but I have not been able to find a reference. Let the simplex $\sigma$ have dimension $n$, and let $x_0,\dots,x_n$ be barycentric coordinates. The expansion of $(x_0 + \dots x_n)^d=1^d=1$ is the sum of terms of the form $$\frac{d!}{i_0!\dots i_n!} x_0^{i_0}\dots x_n^{i_n}.$$ These terms form the Bernstein basis for the vector space of homogeneous polynomials of degree $d$ in the barycentric coordinates. The sum of the basis elements is 1, and there are $n+d \choose d$ of them (tip: Google "stars and bars"). It will turn out that the integral over $\sigma$ is independent of the decomposition $d=i_0+\dots+i_n$. It will follow that the integral of any Bernstein basis vector is $${\text{vol}(\sigma)}\bigg/{{n+d}\choose d}.$$

To prove the above assertion, we only need to prove it for a single $n$-simplex. The formula for other simplices will then follow by the change of variable formula for integration. First note that $$\frac{1}{i!j!}\int_{t=s}^{t=1} (t-s)^i (1-t)^j\,dt =\frac{(1-s)^{i+j+1}}{(i+j+1)!},$$ which is obvious when $i=0$, and, for $i>0$ follows by induction on $i$, using integration by parts to reduce $i$ and increase $j$. Taking $s=0$, we get the required result when $n=1$.

We have coordinates $0\le t_1\le\dots\le t_n\le 1$ for a certain standard $n$-simplex, with respect to which $x_i = t_{i+1}-t_i$ for $1\le i < n$, $x_0=t_1$ and $x_n=1-t_n$. We want to evaluate $$\frac{d!}{i_0!\dots i_n!}\int_{0\le t_1 \le \dots\le t_n\le 1} t_1^{i_0}(t_2-t_1)^{i_1}\dots (t_n-t_{n-1})^{i_{n-1}}(1-t_n)^{i_n} \,dt_1\dots dt_n .$$

From the integral with respect to $t$ above, we see that we can eliminate $t_n$ by integrating over the interval $[t_{n-1},1]$. We obtain the integral of another Bernstein basis element of degree $d+1$ over the standard simplex of dimension $n-1$. Explicitly, we obtain $$\frac{d!}{i_0!\dots i_{n-2}!(i_{n-1}+i_n+1)!}\int_{0\le t_1 \le \dots \le t_{n-1}\le 1} t_1^{i_0}(t_2-t_1)^{i_1}\dots (1-t_{n-1})^{i_{n-1}+i_n+1} \,dt_1\dots dt_{n-1} ,$$ which is $1/(d+1)$ times the integral of a Bernstein basis vector in dimension $n-1$ of degree $d+1$. This proves by induction on $n$ that the result is independent of the decomposition $d=i_0+\dots i_n$ and the proof is complete.

Answer 2

Yes, this is based on Brion's formula, discussed lucidly in Nicole Berline's notes here.

To address the OP's complaints:

1. For the standard simplex, defined by $\sum_{i=1}^n x_i = 1, \quad x_i \geq 0,$ the Euclidean coordinates and the baricentric coordinates are the same.

2. The notes begin by defining Maple functions to compute the integral of monomials over the simplex.

3. Brion's formula is easy in this case (not that I understand why this would be relevant - nowhere in the question does it say to derive the formula starting from Peano axioms). In fact, the authors of the notes give the proof as an exercise.

4. For this particular simplex, writing down the formula for the monomial from its Laplace transform is very easy (you don't need Maple).

Answer 3

For what it worth, here is a not-inductive calculation. Let $\alpha_1,\ldots,\alpha_n$ be arbitrary complex numbers with $\mathrm{Re}\,\alpha_k>-1$ for all $k=1,2,\ldots,n$. We compute the integral over the simplex $S:=\{(x_1,\ldots,x_n)\in \mathbb{R}^n\colon x_i\geqslant 0,\sum x_i=1\}$ of the monomial $\prod_k x_k^{\alpha_k}$. Denote the value of this integral by $I$ and consider $$ J:=\int_{x_1\geqslant 0,\ldots,x_n\geqslant 0} e^{-\sum x_k}\prod x_k^{\alpha_k}dx_1\ldots dx_n. $$ From one side, $J=\prod_k \Gamma(\alpha_k+1)$, since the integrand factorizes as a product of functions $x_k^{\alpha_k}e^{-x_k}$ depending on single variables. From the other side, we denote $(x_1,\ldots,x_n)=r(y_1,\ldots,y_n)$, where $y=(y_1,\ldots,y_n)\in S$ and make a change of variables with new coordinates being $(y,r)\in S\times (0,\infty)$. Computing the Jacobian of this change of variables we get $r^{n-1}/\sqrt{n}$, thus $$ J=n^{-1/2}\int_0^\infty \int_S e^{-r}r^{n-1+\sum \alpha_k}\prod y_k^{\alpha_k}dr dy=n^{-1/2} \Gamma(n+\sum \alpha_k)\cdot I, $$ thus $$I=\sqrt{n}\frac{\Gamma(n+\sum \alpha_k)}{\prod \Gamma(\alpha_k+1)}.$$

Answer 4

There are several methods. Here is a recent one. I implemented it in Python for $n=3$ (and for a polynomial, not only a monomial):

from math import factorial
from sympy import Poly
from sympy.abc import x, y, z
import numpy as np

def term(Q, monom):
    coef = Q.coeff_monomial(monom)
    powers = list(monom)
    j = sum(powers)
    if j == 0:
        return coef
    coef = coef * np.prod(list(map(factorial, powers)))
    return coef / np.prod(list(range(4, 4+j)))

def integral(P, v1, v2, v3, v4):
    X = Poly(x, x, y, z, domain="RR")
    Y = Poly(y, x, y, z, domain="RR")
    Z = Poly(z, x, y, z, domain="RR")
    B = np.column_stack((v1-v4, v2-v4, v3-v4))
    newx = B[0,0]*X + B[0,1]*Y + B[0,2]*Z + v4[0]
    newy = B[1,0]*X + B[1,1]*Y + B[1,2]*Z + v4[1]
    newz = B[2,0]*X + B[2,1]*Y + B[2,2]*Z + v4[2]
    let = {x: newx.as_expr(), y: newy.as_expr(), z: newz.as_expr()}
    Q = P.subs(let, simultaneous=True).as_expr().as_poly(x, y, z)
    monoms = Q.monoms()
    s = 0.0
    for monom in monoms:
        s = s + term(Q, monom)
    return np.abs(np.linalg.det(B)) / 6 * s

###############################################################################
# tetrahedron vertices
v1 = np.array([1.0, 1.0, 1.0]) 
v2 = np.array([2.0, 2.0, 3.0]) 
v3 = np.array([3.0, 4.0, 5.0]) 
v4 = np.array([3.0, 2.0, 1.0])

# polynomial to integrate    
P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR")

# integral
val = integral(P, v1, v2, v3, v4)
print(val)

---

EDIT

For an arbitrary $n$:

from math import factorial
from sympy import Poly
import numpy as np

def term(Q, monom):
    coef = Q.coeff_monomial(monom)
    powers = list(monom)
    j = sum(powers)
    if j == 0:
        return coef
    coef = coef * np.prod(list(map(factorial, powers)))
    n = len(monom)
    return coef / np.prod(list(range(n+1, n+j+1)))

def integral(P, S):
    gens = P.gens
    n = len(gens)
    dico = {}
    v = S[n,:]
    columns = []
    for i in range(n):
        columns.append(S[i,:] - v)    
    B = np.column_stack(tuple(columns))
    for i in range(n):
        newvar = v4[i]
        for j in range(n):
            newvar = newvar + B[i,j]*Poly(gens[j], gens, domain="RR")
        dico[gens[i]] = newvar.as_expr()
    Q = P.subs(dico, simultaneous=True).as_expr().as_poly(gens)
    monoms = Q.monoms()
    s = 0.0
    for monom in monoms:
        s = s + term(Q, monom)
    return np.abs(np.linalg.det(B)) / factorial(n) * s

###############################################################################
# tetrahedron vertices
v1 = np.array([1.0, 1.0, 1.0]) 
v2 = np.array([2.0, 2.0, 3.0]) 
v3 = np.array([3.0, 4.0, 5.0]) 
v4 = np.array([3.0, 2.0, 1.0])
S = np.array([v1, v2, v3, v4])

# polynomial to integrate
from sympy.abc import x, y, z
P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR")

# integral
val = integral(P, S)
print(val)

---

Julia implementation:

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


#### --- EXAMPLE --- ####

# define the polynomial to be integrated
@polyvar x y z
P = x^4 + y + 2*x*y^2 - 3*z

#= be careful
if your polynomial does not involve one of the variables, 
e.g. P(x,y,z) = x⁴ + 2xy², it must be defined as a 
polynomial in x, y, and z; you can do:
@polyvar x y z
P = x^4 + 2*x*y^2 + 0.0*z
then check that variables(P) returns (x, y, 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]

# integral
integratePolynomialOnSimplex(P, S)

---

R implementation:

library(spray)

integratePolynomialonSimplex <- function(P, S) {
  n <- ncol(S)
  v <- S[n+1L, ]
  B <- t(S[1L:n, ]) - v
  gens <- lapply(1L:n, function(i) lone(i, n))
  newvars <- vector("list", n)
  for(i in 1L:n) {
    newvar <- v[i]
    Bi <- B[i, ]
    for(j in 1L:n) {
      newvar <- newvar + Bi[j] * gens[[j]]
    }
    newvars[[i]] <- newvar
  }
  Q <- 0
  exponents <- P[["index"]]
  coeffs    <- P[["value"]] 
  for(i in 1L:nrow(exponents)) {
    coef <- coeffs[i]
    powers <- exponents[i, ]
    term <- 1
    for(j in 1L:n) {
      term <- term * newvars[[j]]^powers[j] 
    }
    Q <- Q + coef * term
  }
  s <- 0
  exponents <- Q[["index"]]
  coeffs    <- Q[["value"]] 
  for(i in 1L:nrow(exponents)) {
    coef <- coeffs[i]
    powers <- exponents[i, ]
    j <- sum(powers)
    if(j == 0L) {
      s <- s + coef
      next
    }
    coef <- coef * prod(factorial(powers))
    s <- s + coef / prod((n+1L):(n+j))
  }
  abs(det(B)) / factorial(n) * s
}

############################# --- EXAMPLE --- ##################################
library(spray)

# variables
x <- lone(1, 3)
y <- lone(2, 3)
z <- lone(3, 3)
# polynomial
P <- x^4 + y + 2*x*y^2 - 3*z

# simplex (tetrahedron) vertices
v1 <- c(1, 1, 1)
v2 <- c(2, 2, 3)
v3 <- c(3, 4, 5)
v4 <- c(3, 2, 1)
# simplex
S <- rbind(v1, v2, v3, v4)

# integral
integratePolynomialonSimplex(P, S)
```