Efficient Multiple Importance Sampling Estimators

I Introduction

Importance sampling (IS) methods approximate statistical moments of a variable of interest by sets of samples, drawn from a proposal distribution different from the targeted one, and weights, assigned to the samples in order to measure their adequacy in approximating the target . In its standard form, IS uses one proposal distribution from which all samples are drawn. However, a more powerful strategy to reduce the variance of the estimators consists of using a set of different proposal distributions. This is the basis of multiple importance sampling (MIS) techniques . The samples drawn from the different proposals in MIS are then assigned weights proportional to the ratio between the target and the proposal densities evaluated at the sample values. Several strategies for calculating the weights have been considered, depending on the way in which the proposal evaluation in the denominator is interpreted. In the standard MIS, the evaluated proposal is exactly the one from which the sample was generated. This constitutes the simplest method in terms of computational complexity. A different approach, referred to as deterministic mixture (DM) MIS, interprets the set of generating proposals as a whole mixture distribution and calculates the weight of a given sample by considering the entire mixture as the global proposal . This method attains a variance reduction at the expense of an increase in the computational load .

In this work, we propose a novel MIS method that provides an efficient tradeoff in terms of computational complexity and variance reduction of the associated IS estimators. The approach is based on creating a partition of the set of available proposals and considering that each partitioned set constitutes a mixture distribution. A sample drawn from a mixand in one of the partitions (mixtures) is then assigned a weight that only accounts for that particular mixture, instead of the entire composite mixture as in the full DM-MIS. A remarkable reduction in computational complexity is achieved by this approach, while the variance of the associated partial DM-MIS estimator remains comparable to that of the full DM-MIS estimator. The method can not only applied to IS methods with static distributions (i.e., characterized by fixed parameters), but also to IS methods that adapt the parameters of the proposal distribution in an iterative way (i.e., adaptive IS (AIS) methods ). Computer simulations show the excellent performance of the proposed scheme in terms of variance reduction for a given computational load.

II Problem statement and background

where $f(\cdot)$ can be any integrable function of ${\bf x}$. In many practical scenarios, we cannot obtain an analytical solution of (2) and Monte Carlo methods are used to approximate it.

Let us consider $N$ samples (${\bf x}_{1},\ldots,{\bf x}_{N}$) drawn from a proposal pdf, $q({\bf x})$, with heavier tails than the target, $\pi({\bf x})$. The samples have associated importance weights given by

Using the samples and weights, the moment of interest can be approximated as

where $\hat{Z}=\frac{1}{N}\sum_{j=1}^{N}w_{j}$ is an unbiased estimator of $Z=\int_{\mathcal{D}}\pi({\bf x})d{\bf x}$ . Note that Eq. (4) always provides a consistent estimator of $I$, but its variance is directly related to the discrepancy between $\pi({\bf x})|f({\bf x})|$ and $q({\bf x})$ (for a specific choice of $f$) or to the mismatch between the target $\pi({\bf x})$ and the proposal $q({\bf x})$ (for a general and arbitrary $f$) .

II-B Multiple Importance Sampling (MIS)

Finding a good proposal pdf, $q({\bf x})$, is critical and can also be very challenging . An alternative and more robust strategy consists of using a population of different proposal pdfs. This scheme is often known in the literature as multiple importance sampling (MIS) . Let us consider a set of $N$ (normalized) proposal pdfs, $q_{1}({\bf x}),\ldots,q_{N}({\bf x})$, and let us assume that exactly one sample is drawn from each of them, i.e., ${\bf x}_{i}\sim q_{i}({\bf x})$, $i=1,...,N$. The importance weights associated to these samples can then be obtained according to one of the following strategies:

$\mbox{{\it Standard} MIS:}\;\;\;w_{i}=\frac{\pi({\bf x}_{i})}{{q_{i}({\bf x}_{i})}},\quad i=1,\ldots,N.$

where $\psi({\bf x})=\frac{1}{N}\sum_{j=1}^{N}q_{j}({\bf x})$ is the mixture pdf, composed of all the proposal pdfs. This approach interprets the complete set of samples, $\{{\bf x}_{i}\}_{i=1}^{N}$, as being distributed according to the mixture $\psi({\bf x})$, i.e., $\{{\bf x}_{1},\ldots,{\bf x}_{N}\}\sim\psi({\bf x})$. See Appendix A for further details.

In both cases, the consistency of the estimators is ensured. The main advantage of the DM-MIS weights is that they yield more stable and efficient estimators , i.e., with less variance (as proved in Appendix B). However, the DM-MIS estimator is computationally more expensive, as it requires $N$ evaluations of the proposal to obtain each weight instead of just one. Note that the number of evaluations of the target $\pi({\bf x})$ is the same regardless of whether the weights are calculated according to (3) or (5). In some practical scenarios, this additional load may be excessive (especially for large values of $N$) and alternative efficient solutions must be developed.

II-C Adaptive Importance Sampling (AIS)

III Partial Deterministic Mixture approach

In this section we develop a partial DM-MIS scheme, which groups the proposal pdfs, $\{q_{i}({\bf x})\}_{i=1}^{N}$, into $P$ mixtures composed of $M$ mixands, with $PM=N$ (recall that ${\bf x}_{i}\sim q_{i}({\bf x}),i=1,\cdots,N$). For the sake of simplicity, we assume that all the mixtures contain the same number of proposal pdfs. However, any strategy that clusters the $N$ proposals into $P$ disjoint mixtures (regardless of their size) is valid. Namely, we define a partition of $\{1,\ldots,N\}$ into $P$ disjoint subsets of $M$ indices, $\mathcal{S}_{p}$ with $p=1,\ldots,P$, s.t.

The resulting partial DM-MIS (p-DM-MIS) estimator is then given by

which coincides with the expression in (4), but using the weights given by Eq. (7). Note that the particular cases $P=1$ and $P=N$ correspond to the full DM-MIS (f-DM-MIS) and the standard MIS (s-MIS) approaches, respectively. The number of evaluations of the proposal pdfs is then $PM^{2}$. Since $N\leq PM^{2}=NM\leq N^{2}$, the computational cost is larger than that of s-MIS approach ($M$ times larger), but lower than that of the f-DM-MIS approach (since $M\leq N$).

The good performance of the novel approach is ensured by Theorem 1 and Corollary 2 (see Appendix B) and can be summarized by the following expression:

which holds regardless of the choice of $P$ and the strategy followed to group the original proposals $\{q_{i}({\bf x})\}_{i=1}^{N}$ into mixtures. Therefore, there is a tradeoff between performance and computational cost: using a smaller number of mixtures ($P$) leads to a reduced variance, but at the expense of an increase in the number of evaluations of the proposal densities.

A simple strategy to choose the number of mixtures consists of starting with $P=N$ (which coincides with the standard MIS scheme), computing the corresponding estimator in Eq. (8), and iteratively reducing the number of mixtures $P$ (thus increasing $M=\frac{N}{P}$) while the estimation significantly changes w.r.t. the previous step. This iterative approach does not require a significant additional computational cost (since the proposal evaluations can be stored and re-used) and results in an efficient and automatic procedure to select $P$.

III-B Design of the PP mixtures

Developing an optimal strategy to cluster the proposals into $P$ mixtures is a difficult task, since the number of possible configurations is extremely large. Indeed, unless this clustering strategy is computationally inexpensive, the additional computational effort might be better invested in decreasing $P$ (thus increasing $M=\frac{N}{P}$ and reducing the variance of the estimators). Therefore, we propose applying a simple random clustering strategy, where $M=\frac{N}{P}$ different proposals are randomly assigned to each partition. This approach provides an excellent performance for large values of $M$ (see Table II), so there seems to be little room for improvement (except maybe for small values of $M$). Note that this result is not surprising: randomness is the key element behind compressive sampling , and many randomized clustering algorithms have been developed for applications such as data mining , image processing or blind channel identification .

III-C Application to AIS schemes

For the sake of simplicity, we have focused on p-DM-MIS for a static framework, where the parameters of the proposals are fixed. However, all the previous considerations can be easily applied in AIS schemes, where the proposals are iteratively updated. In methods like PMC or APIS , a population of $J$ proposal densities is adapted during $T$ iterations. At the $t$-th iteration, the $j$-th sample ${\bf x}_{j}^{(t)}$ is drawn from the $j$-th proposal, i.e., ${\bf x}_{j}^{(t)}\sim q_{j}^{(t)}({\bf x})$ for $j=1,\ldots,J$ and $t=1,\ldots,T$. Thus, after $T$ iterations we have $N=JT$ samples drawn from $N=JT$ different proposal pdfs.

Regardless of the adaptation procedure followed by each algorithm, different strategies can be used to design an efficient DM-MIS estimator when considering all the $N=JT$ proposal pdfs. Table I summarizes the weight calculation for three well-known AIS methods (PMC, AMIS and APIS), comparing them to the DM-MIS estimators. We also analyze the complexity in terms of the number of proposal evaluations and the performance in terms of the variance reduction. We can see that the novel p-DM-MIS approach provides a very good compromise between performance and computational cost. Finally, it is important to remark that Table I does not take into account the specific adaptive procedures of each algorithm, which can also have a large influence on the final performance.

IV Numerical example

We consider a bivariate multimodal target pdf, defined as a mixture of $5$ Gaussians:

We apply the MIS algorithm in a setup with $N=4096$ Gaussian proposal pdfs, $\{q_{i}({\bf x})=\mathcal{N}({\bf x};{\bm{\mu}}_{i},{\bf C}_{i})\}_{i=1}^{N}$, where ${\bm{\mu}}_{i}\sim\mathcal{U}(\times)$ is the randomly chosen location parameter and ${\bf C}_{i}=\sigma^{2}{\bf I}_{2}$, with $\sigma=5$, is the scale matrix. We proceed as follows. First, we draw a sample from each proposal. Then, we compute the corresponding weight according to (7). Finally, we build the estimator $\hat{I}_{\textrm{p-DM-MIS}}$ using (8) for different number of mixtures $P=2^{k}$ for $k=0,1,\ldots,12$ (i.e., $P=1,2,\ldots,2048,4096$). Since $PM=N$, the number of proposals per mixture is $M=2^{12-k}=4096,2048,\ldots,2,1$. Note that the case $P=1$ (i.e., $M=4096$) corresponds to the f-DM-MIS approach, while $P=4096$ (i.e., $M=1$) corresponds to the s-MIS. A random assignment of the proposals to the mixtures is performed in all cases.

Table II shows the mean square error (MSE) in the estimation of the mean of the target $E[{\bf X}]$ (averaged over both dimensions) and the normalizing constant $Z$. All results are averaged over $500$ independent runs. The last column of Table II shows the total number of proposal evaluations for each value of $M$. Figure 1 shows the MSE in the estimation of the mean of the target w.r.t. the total number of proposal evaluations. The results show that decreasing $P$ reduces the MSE (as expected) at the expense of an increase in the computational cost (measured by the number of proposal evaluations). However, note that decreasing the number of mixtures below $P=64$ does not improve the performance significantly. Indeed, the p-DM-MIS with $P=64$ obtains an MSE close to that of the f-DM-MIS estimator while performing $98.4\%$ less proposal evaluations, thus attaining an excellent performance-cost tradeoff.

V Conclusion

In this paper we propose a novel approach for the calculation of the weights in multiple importance sampling schemes that provides an efficient tradeoff between computational complexity and variance reduction of the associated estimator. The proposed scheme is based on constructing a random partition of the set of available proposals and then calculating the weight of each sample locally using only the corresponding subset of the partition. Computer simulations reveal a very good performance of the method, which is able to attain an excellent performance-computational cost tradeoff.

Appendix A Drawing samples from a mixture of pdfs

Let us consider a mixture of $N$ normalized pdfs with equal weights, i.e.,

The classical procedure for drawing $N$ samples from $\psi({\bf x})$ is (starting with $k=1$):

Draw $j^{*}\in\{1,\ \ldots,\ N\}$ with equal weights $\frac{1}{N}$.

Draw ${\bf x}_{k}\sim q_{j^{*}}({\bf x})$.

In this way, each sample ${\bf x}_{k}$ is distributed according to $\psi({\bf x})$ and, as a consequence, $\{{\bf x}_{1},\ \ldots,\ {\bf x}_{N}\}\sim\psi({\bf x})$. An alternative procedure, more deterministic than the previous one, consists of the following steps (starting with $i=1$):

Draw one sample from each $q_{i}({\bf x})$, i.e., ${\bf x}_{i}\sim q_{i}({\bf x})$.

In this case, we have ${\bf x}_{i}\sim q_{i}({\bf x})$ for $i=1,\dots,N$, but the joint set is again distributed as $\psi({\bf x})$, i.e., $\{{\bf x}_{1},{\bf x}_{2},\ \ldots,\ {\bf x}_{N}\}\sim\psi({\bf x})$. This is the underlying theoretical argument of the deterministic mixture (DM) weights approach used throughout this work. Furthermore, given $M$ indices $j_{r}\in\{1,\ldots,N\}$ with $r=1,\ldots,M$:

Appendix B Variance reduction of deterministic mixture weights

The variance of the f-DM-IS estimator is always lower or equal than the variance of the s-MIS estimator,

For the standard MIS, full DM-MIS and partial DM-MIS estimators the following inequality holds:

Proof: For the first inequality, note that the full DM-MIS estimator of (12b) can also be expressed as

is the proposal associated to the $p$-th mixture. Similarly, the partial DM-MIS estimator is given by

Hence, applying Theorem 1 with $\psi_{p}({\bf x}_{i})$ instead of $q_{i}({\bf x}_{i})$ and $\frac{1}{P}\sum_{p=1}^{P}{\psi_{p}({\bf x}_{i})}$ instead of $\frac{1}{N}\sum_{j=1}^{N}{q_{j}({\bf x}_{i})}$, we have $\textrm{Var}(\hat{I}_{\textrm{f-DM-MIS}})\leq\textrm{Var}(\hat{I}_{\textrm{p-DM-MIS}})$.

For the second inequality, following the same approach, the standard MIS estimator of (12a) can be rewritten as

Applying Theorem 1 with $\psi_{p}({\bf x}_{i})=\frac{1}{P}\sum_{p=1}^{P}{\psi_{p}({\bf x}_{i})}$ instead of $\frac{1}{N}\sum_{j=1}^{N}{q_{j}({\bf x}_{i})}$, we have $\textrm{Var}(\hat{I}_{\textrm{p-DM-MIS}})\leq\textrm{Var}(\hat{I}_{\textrm{s-MIS}})$. $\hfill\Box$

References