Use Importance sampling for multimodal and multivariate distribution draws, how to choose proposal distribution?

Question

I'm in trouble trying to generate samples following a particular distribution which is not numerically known perfectly. Let us consider a $R^n$ space provided with an orthonormal base $( e_{1},...,e_{n} )$ and N points $(X^{(1)},...,X^{(N)})$ in this $R^n$ space (the family of point is growing during the process). I want to draw uniform samples in the following subset of R^n : $$\begin{aligned}SubSet = R^n & - \bigcup_{i=1}^{N}\left (\{\zeta \in R^{n}, \zeta_{j} \leq X^{(i)}_{j}, \forall j \in (1,\dots,n) \} \right )\\& - \bigcup_{i=1}^{N}\left (\{\zeta \in R^{n}, \zeta_{j} \geq X^{(i)}_{j}, \forall j \in (1,\dots,n) \} \right ), \end{aligned} $$ where $\cdot_{j}$ refers to the $j^{th}$ component. When the space dimension is taken as dimension $2$, one can obtain the following figure : in red the SubSet, in green the $X^{(i)}$ points, in blue the areas to discard

As the $(X^{(i)})$ family is growing during the process, my SubSet becomes more and more small during the process. In dimension 2, the red areas are easy to compute and symbolic calculation is OK for a classic computer. Thus it is easy to draw uniform samples. But when the dimension arises, it becomes really complicated. I took some time looking at importance sampling in order to draw samples in the red areas without calculating them but the proposal distribution is really hard to define since it is both multimodal and multivariate distribution.

Do you think importance sampling is a good candidate that deserve some exploration ? How would you define the proposal distribution ? Do you have any other idea of how to do such thing ?

P.S. : From any point in the $R^n$ space, it is numerically not costly to know if it is in blue area or not since the number of points $X^{(i)}$ is rather low.

Answer 1

Depending on the shape of the domain $\mathcal{X}$ and how $N$ (the number of points) scales with $n$ (the dimension), plain-old rejection sampling might suffice. Specifically, sample uniformly from $\mathcal{X}$ and check if the sample lies in $SubSet$ (which is efficient to do); if not, draw another sample and try again. This method is efficient if (1) it's easy to sample from $\mathcal{X}$, and (2) $SubSet$ (i.e. the "red region") is a fairly large fraction of $\mathcal{X}$.

Since $\mathcal{X}$ isn't specified, it's not clear if that is true. With the hope of helping the original poster, let's consider one special case. Suppose that $\mathcal{X}$ is the unit hypercube $H_n=[0,1]^n$. Let us suppose that the $X_j$ lie in the interior of $H_n$, with a little extra margin: $$ \epsilon < X_j^{(i)} < 1-\epsilon $$ for some $0<\epsilon<1/2$.

In the language of the original poster, each point $X_j$ produces two forbidden "blue regions". However, as the dimension $n$ increases, the fraction of space taken up by these forbidden regions represents an exponentially smaller fraction of $\mathcal{X}$. Specifically, for any $i$, $$ Volume\left[ \{\zeta \in H_n, \zeta_{j} \leq X^{(i)}_{j}, \forall j \in (1,\dots,n) \} \right] < (1-\epsilon)^n $$ and $$ Volume\left[ \{\zeta \in H_n, \zeta_{j} \geq X^{(i)}_{j}, \forall j \in (1,\dots,n) \} \right] < (1-\epsilon)^n $$ Therefore, since $Volume(H_n)=1$, the probability $p$ of an arbitrary point in the $H_n$ lying within $SubSet$ is at least $$ p \geq 1 - 2N(1-\epsilon)^n $$ That means that if we use rejection sampling, the probability of a single sample of $H_n$ being accepted is $p$. For a fixed $\epsilon$ and assuming $N$ is scaling at most polynomially with $n$ (i.e., $N$ doesn't grow too fast), then $$ \lim_{n\rightarrow \infty} p = 1 $$ convergence is exponentially fast, and rejection sampling is very efficient.

Answer 2

Bill, Thank you for your post and your detailed developments. Rejection method is actually not really relevant in here as my code is giving me more and more $X^{(i)}$. Thus the subset domain size is decreasing to the point where the volume of the SubSet becomes really small compared to the whole domain $\mathcal{X}$ considered.

For the first iterations, rejection method could work... but after a few iteration, uniform draws into the $\mathcal{X}$ space will results with only a little fraction of draws into the SubSet of interest.

I edited my initial post to give some precision about the fact that the family $(X^{(1)},\dots,X^{(N)})$ is growing.