integration – Approximation of integral of gaussian function over a parallelepiped

Remark: I posted this question in math stackexchange here and computer science stackexchange https://cs.stackexchange.com/ few weeks ago but obtain no answer.

Given a multi-dimensional gaussian function, defined by
$$f(boldsymbol{x})=expleft{-frac{1}{2} boldsymbol{x}^Tboldsymbol{x} right}=expleft{-frac{1}{2} sum_{i=1}^nx_i^2 right}$$
And an integration region as the form of a $n$-dimensional parallelepiped, defined by
$$mathcal{D} = left{boldsymbol{l le Lxle u} right}$$
with

  • the lower triangular matrix $boldsymbol{L}in Bbb R^{ntimes n}$ where all lower elements and diagonal equal to $1$, all upper elements equal to $0$
    $$
    boldsymbol{L} = left(
    begin{matrix}
    1&0&0&ldots&0\
    1&1&0&ldots&0\
    1&1&1&ldots&0\
    vdots&vdots&vdots&ddots&0\
    1&1&1&ldots&1\
    end{matrix}
    right)
    $$
  • the vectors $boldsymbol{l,u} in Bbb R^n$: $boldsymbol{l} = (l_1,…,l_n)’$ and $boldsymbol{u} = (u_1,…,u_n)’$

The integration region

Are there any methods/algorithms that we can use to approximate the integral of $f(boldsymbol{x})$ over $mathcal{D}$
$$int_{mathcal{D}}f(boldsymbol{x})dboldsymbol{x}=int_{{boldsymbol{l le Lxle u} }}expleft{-frac{1}{2} boldsymbol{x}^Tboldsymbol{x} right}dboldsymbol{x}$$
satisfying

  • Fast computation (because later I must compute many integrals with different values of $boldsymbol{l,u}$)
  • The accuracy doesn’t need to be high (absolute error less than $10^{-3}$ is sufficient)

My attempt: we may use Monte Carlo simulation to approximate this integral but given the very specific form of the integration region and also the integrand, I hope there may exist a faster numerical method/algorithm/closed-form approximation.

Besides, we notices that the inversion matrix $boldsymbol{L}^{-1}$ is an upper bi-diagonal matrix
$$
boldsymbol{L}^{-1} = left(
begin{matrix}
1&-1&0&ldots&0\
0&1&-1&ldots&0\
0&ddots&ddots&ddots&0\
vdots&ddots&ddots&ddots&-1\
0&ldots&ddots&ddots&1\
end{matrix}
right)
$$

So, by making a change of variable $boldsymbol{y = Lx}$, we can transform the integral into $$int_{mathcal{D}}f(boldsymbol{x})dboldsymbol{x}= int_{{boldsymbol{l’ le y le u’}}} exp left{-frac{1}{2} left(y_n^2+sum_{i=1}^{n-1} (y_i-y_{i+1})^2 right) right} dboldsymbol{y}$$
with $mathcal{D}’ = {boldsymbol{l’ le y le u’}}$ is a rectangular region.

Thank you in advance!