Very fast but inaccurate estimations of multivariate Gaussian integral over a hypercube

110 Views Asked by At

$\def\Z{\mathbb{Z}}\def\R{\mathbb{R}}\def\A{\mathcal{A}}\def\N{\mathcal{N}}$ I'm working on 4D positive real values, i.e. $\R^4_{\geq 0}$, where it is gridized with hypercubes of side length $a > 0$. So, an index vector $I = [i ~~ j ~~ k ~~ \ell]^T$ uniquely represents the region

$$\A_I := \{ x \in \R^4_{\geq 0} ~~|~~ a I_n \leq x_n < a (I_n+1),~~ n = 1,2,3,4 \}$$

where $x_n$ and $I_n$ represents the $n$th entity of $x$ and $I$ respectively.

I have a multivariate Gaussian distribution in $\R^4$, which is $N(m,P)$, i.e. with mean $m$ and covariance matrix $P$ where $m \in \R^4$ and $P\in\R^{4\times 4}$ is symmetric. Let $m \in \A_I$ for some $I$ and $\N_I := \{J \in \Z^4 ~~|~~ I_n-1 \leq J_n \leq I_n+1\}$ be the index set that represents the immediate neighbours of $\A_I$ (including itself). For each $J \in \N_I$ I want to estimate the Gaussian integral over $\A_J$, i.e.

$$\alpha_J := |2 \pi P|^{-1/2} \int_{\A_J} \exp\left(-\frac{1}{2} (x-m)^T P^{-1} (x-m)\right) dx $$

where $|\cdot|$ represents the determinant.

So, I need to estimate $3^4=81$ numbers very fast, e.g. under 1 milliseconds on a typical mid-range laptop if possible. The estimations do not need to be very accurate though, just rough values.

  1. One method that comes to mind is to calculate the integrand for the vertices and midpoint of regions, and use trapezoidal integration. This would require the calculation of the integrand $4^4+3^4=337$ times, which might be too much.
  2. Another method could be to calculate the integrand for only the midpoints and use them as estimations or maybe use some sort of interpolation and averaging centered around $m$.

Although both methods can run very fast, I'm not sure they would create relevant estimations especially when $m$ is close to some boundary between regions.

Is there a name for these kind of problems? Is there very fast and robust algorithms (produces relevant results) suited to this problem?

Note: The values of $m$ and $P$ are not known beforehand. So, a precomputation table cannot be used (I guess). But if the values in a table for a specific $P$ can be used by adjusting the values for the incoming $P$, that would be ok.