Computationally Efficient Nonparametric Importance Sampling

The variance reduction established by importance sampling strongly depends on the choice of the importance sampling distribution. A good choice is often hard to achieve especially for high-dimensional integration problems. Nonparametric estimation of the optimal importance sampling distribution (known as nonparametric importance sampling) is a reasonable alternative to parametric approaches.In this article nonparametric variants of both the self-normalized and the unnormalized importance sampling estimator are proposed and investigated. A common critique on nonparametric importance sampling is the increased computational burden compared to parametric methods. We solve this problem to a large degree by utilizing the linear blend frequency polygon estimator instead of a kernel estimator. Mean square error convergence properties are investigated leading to recommendations for the efficient application of nonparametric importance sampling. Particularly, we show that nonparametric importance sampling asymptotically attains optimal importance sampling variance. The efficiency of nonparametric importance sampling algorithms heavily relies on the computational efficiency of the employed nonparametric estimator. The linear blend frequency polygon outperforms kernel estimators in terms of certain criteria such as efficient sampling and evaluation. Furthermore, it is compatible with the inversion method for sample generation. This allows to combine our algorithms with other variance reduction techniques such as stratified sampling. Empirical evidence for the usefulness of the suggested algorithms is obtained by means of three benchmark integration problems. As an application we estimate the distribution of the queue length of a spam filter queueing system based on real data.

💡 Research Summary

This paper addresses a central challenge in importance sampling (IS) for high‑dimensional integration: the difficulty of selecting an effective proposal distribution. While the optimal proposal is proportional to |ϕ| p (or |ϕ−Iϕ| p for self‑normalized IS), it is rarely available in closed form, and parametric approximations often fail to capture the true shape. The authors therefore develop non‑parametric importance sampling (NIS) algorithms that estimate the optimal proposal directly from data, but they replace the computationally heavy kernel density estimator with a linear‑blend frequency polygon (LBFP) estimator, also known as a multivariate histogram with linear interpolation.

The paper proposes two algorithms. Algorithm 1 (NIS) first draws M samples from a trial distribution q₀, computes importance weights w_j = |ϕ(x̃_j)| p(x̃_j)/q₀(x̃_j), and builds an LBFP estimate \hat q_{IS} of the optimal proposal. In a second stage, N−M samples are drawn from \hat q_{IS} and the standard IS estimator \hat I_{NIS} = (N−M)^{-1} Σ ϕ(x_i) p(x_i)/\hat q_{IS}(x_i) is evaluated. Algorithm 2 (NSIS) follows the same two‑stage structure but targets the self‑normalized estimator, using weights based on |ϕ−Î| p.

A rigorous asymptotic analysis is carried out under smoothness and boundedness assumptions on ϕ and p. For non‑negative (or non‑positive) integrands, the authors prove that with an optimally chosen bin width h* = C M^{-1/(d+4)} the mean‑square error (MSE) of the NIS estimator decays as O(N^{-(d+8)/(d+4)}), which is faster than the standard Monte Carlo rate O(N^{-1}) and matches the optimal variance achievable by the true importance distribution. They also derive the optimal proportion λ* = M/N of samples to allocate to the proposal‑estimation stage: λ* = 4/(d+8). This result quantifies the trade‑off between learning the proposal and using it for estimation, and shows that a larger fraction of the total budget must be spent on learning as the dimension d grows.

When ϕ takes both positive and negative values, the optimal variance cannot be attained directly because the optimal proposal no longer yields zero variance. The paper shows that, with the same optimal bandwidth, the MSE becomes (I_ϕ²−I_ϕ²)/(N−M) + o(N^{-1}), i.e., the estimator asymptotically reaches the optimal IS variance but does not improve the convergence rate. To recover the faster rate, the authors propose the NIS+/- strategy: decompose ϕ into its positive and negative parts (ϕ = ϕ⁺−ϕ⁻) and apply the NIS algorithm separately to each part. This yields the O(N^{-(d+8)/(d+4)}) rate again.

For the self‑normalized case, analogous results are obtained: the NSIS estimator is biased for finite N but asymptotically unbiased, and its MSE shares the same optimal bandwidth and λ* as the unnormalized version.

A key contribution is the computational advantage of the LBFP estimator. Unlike kernel estimators, LBFP allows O(1) sampling and density evaluation via the inversion method, and its memory footprint is comparable to a simple histogram. Consequently, the overall cost of the two‑stage NIS procedure is dramatically lower, making it feasible for dimensions up to at least d=10 in the experiments.

The empirical section validates the theory on three benchmark problems (multivariate normal density integration, rare‑event option pricing, and a complex nonlinear function) and on a real‑world application: estimating the queue‑length distribution of a spam‑filtering system from observed data. In all cases, the LBFP‑based NIS (and NIS+/- when needed) outperforms kernel‑based NIS, standard parametric IS, and plain Monte Carlo, often achieving several‑fold reductions in MSE for the same computational budget.

In summary, the paper delivers a practically efficient non‑parametric importance sampling framework. By leveraging the linear‑blend frequency polygon, it resolves the computational bottleneck of previous non‑parametric approaches, provides explicit optimal tuning rules (bandwidth and sample split), and demonstrates that the resulting estimator can asymptotically attain the optimal IS variance. The work opens the door to broader adoption of non‑parametric IS in high‑dimensional problems and suggests future extensions such as adaptive bandwidth selection, mixture proposals, and integration with modern machine‑learning based density estimators.