# 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)’$$ 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.