Rebinning algorithm in VEGAS

Question

I am trying to understand the rebinning algorithm of the VEGAS (original publication (preprint from LKlevin) and implementation notes) Monte Carlo integration. I will try to explain first what I think I understood and then pose my questions.

For simplicity, let's assume we have a 1-dimensional function $f(x)$ that is positive over the whole interval $[0,1]$. This interval is seperated into, let's say, $n$ bins. These bins are initially equally-sized. The bin sizes $\Delta x_i$ define a probability density

$$\rho(x) = \begin{cases} 0 \le x < \Delta x_1 : \frac{1}{n \Delta x_1} \\ \vdots \\ \Delta x_{n-1} \le x < \Delta x_n : \frac{1}{n \Delta x_n} \end{cases} .$$

The bin sizes must add up to the length of the interval in order to make $\rho(x)$ properly normalized:

$$\sum_{i=1}^n \Delta x_i = 1 \quad \Rightarrow \quad \int\limits_0^1 \mathrm{d}x \, \rho(x) = 1.$$

Using $N$ randomly chosen numbers $\{ x_i \}$ from the distribution $\rho(x)$ VEGAS computes an approximation $S^{(1)}$ of the integral:

$$S^{(1)} = \frac{1}{N} \sum_{\{x_i\}} \frac{f(x_i)}{\rho(x_i)} \approx \int\limits_0^1 \mathrm{d}x \, f(x)$$

Up to now this is just importance sampling (I am not interested in stratified sampling) using a variably-sized grid. The interesting piece in VEGAS is the rebinning algorithm, i.e. the algorithm which recomputes the bin sizes depending on the function values accumulated in the iteration before:

• For each bin the squared function values (?) are summed up (in the original publication the absolute values are summed up).
• There is also a dampening function applied to each value, to "avoid rapid, destablizing changes".
• After that, each value is smoothed with the neighboring bins. I guess this also adds some stability to the rebinning algorithm for certain functions (but I can not explain why). Let's call the final values $b_i$.
• The bin sizes are now set such that every new bin contains approximately the average:

$$\overline{b} = \frac{1}{n} \sum_{i=1}^n b_i$$

This algorithm makes the bins grow where the function is "small" and shrink where the function is "big". Small and big are understood in relation to each other, e.g. the function's maxima are considered "big" and everything else would be considered "smaller". Since the probability for a point $x$ to end up in any bin is equal (VEGAS is written to exhibit this behavior), the function is sampled most where it is biggest and thereby the error is reduced.

Why is it written that way? Why does one not tackle the problem more directly, for example by taking the averaged and binned function as a probability density for the next iteration?

Answer 1

Disclaimer: I'm very familiar with Monte Carlo algorithms in general, but not with the VEGAS algorithm specifically.

The VEGAS algorithm is based on the fact that if we knew the exact probability density that minimizes the variance of our Monte Carlo integration, we could get our answer using exactly 1 evaluation of the function f.

This is given by

$$p(\mathbf{x}) = \frac {|f(\mathbf{x})|} {\int_\Omega |f(\mathbf{x})| d\mathbf{x}}$$

We don't know the probability density, and estimating it exactly is not feasible for a high-dimensional function as it would quickly take up a massive amount of memory.

Instead the VEGAS algorithm approximates it as an M-step piecewise constant function.

I'm guessing you don't have access to the full article? The original article does NOT use the squared function, but the absolute one (a pre-print version can be found here).

I hope this helps answer your question. The original article answers most of this, so it might be worthwhile getting