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})=\exp\left\{-\frac{1}{2} \boldsymbol{x}^T\boldsymbol{x} \right\}=\exp\left\{-\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 Lx\le u} \right\}$$ with
- the lower triangular matrix $\boldsymbol{L}\in \Bbb R^{n\times 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})d\boldsymbol{x}=\int_{\{\boldsymbol{l \le Lx\le u} \}}\exp\left\{-\frac{1}{2} \boldsymbol{x}^T\boldsymbol{x} \right\}d\boldsymbol{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})d\boldsymbol{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\} d\boldsymbol{y}$$ with $\mathcal{D}' = \{\boldsymbol{l' \le y \le u'}\}$ is a rectangular region.
Thank you in advance!
