Proximal Newton Adaptive Importance Sampler

• The paper introduces a novel adaptive importance sampling method that combines second-order local approximations with proximal operators to handle non-smooth, constrained targets.
• The methodology employs quadratic surrogate approximations and proximal gradient updates to iteratively adapt multivariate Gaussian proposals.
• Empirical results demonstrate orders-of-magnitude improvements in mean squared error and enhanced proposal concentration compared to traditional adaptive sampling techniques.

The Proximal Newton Adaptive Importance Sampler (PNAIS) is an adaptive importance sampling (AIS) methodology devised for efficiently estimating expectations under non-smooth and/or constrained target distributions, especially those encountered in modern statistics, signal processing, and machine learning. PNAIS combines local second-order information from a differentiable component of the target density with proximal operations for handling a non-smooth (convex) component, thus extending the applicability of adaptive proposal adaptation techniques to composite models beyond the reach of gradient-based methods (Elvira et al., 2024).

1. Problem Setting and Target Formulation

PNAIS addresses the estimation of expectations of the form

$$I = \mathbb{E}_{\pi}[h(X)] = \int_{\mathcal{X}} h(x)\, \pi(x)\,dx,$$

where the unnormalized target density $\pi$ on $\mathcal{X} \subseteq \mathbb{R}^d$ is specified by the composite exponential structure

$\pi(x) \propto \exp(-f(x)-g(x)),$

with $f : \mathbb{R}^d \rightarrow \mathbb{R}$ differentiable (possibly non-convex), and $g : \mathbb{R}^d \rightarrow (-\infty, +\infty]$ proper, lower-semicontinuous, and convex (e.g., indicator of a constraint or $\ell_1$ penalty). Standard adaptive IS methodologies that utilize gradient-based adaptation are inapplicable for non-smooth $g$, thus motivating the use of proximal methods for proposal adaptation.

2. Proximal Newton Proposal Adaptation

Each PNAIS proposal is a multivariate Gaussian

$$q_n(x; \mu_n, \Sigma_n) = \mathcal{N}(x|\mu_n,\Sigma_n)$$

where $n$ indexes the population of proposals. Adaptation is carried out iteratively over $\pi$0 rounds. At every step, the mean $\pi$1 is relocated by minimizing a local quadratic surrogate of $\pi$2 plus $\pi$3. Specifically, for proposal $\pi$4 at location $\pi$5, the quadratic approximation

$\pi$6

(with $\pi$7 a Hessian-based scaling) is combined with $\pi$8 to form the subproblem

$\pi$9

This admits the scaled proximal gradient update

$\mathcal{X} \subseteq \mathbb{R}^d$0

where

$\mathcal{X} \subseteq \mathbb{R}^d$1

for $\mathcal{X} \subseteq \mathbb{R}^d$2. The scaling $\mathcal{X} \subseteq \mathbb{R}^d$3, with $\mathcal{X} \subseteq \mathbb{R}^d$4 if positive definite and $\mathcal{X} \subseteq \mathbb{R}^d$5 otherwise, is stabilized via backtracking selection of $\mathcal{X} \subseteq \mathbb{R}^d$6.

When $\mathcal{X} \subseteq \mathbb{R}^d$7 allows, closed-form solutions for $\mathcal{X} \subseteq \mathbb{R}^d$8 are used (e.g., soft-thresholding for $\mathcal{X} \subseteq \mathbb{R}^d$9 penalties, Euclidean projection for constraints); otherwise, a numerical solver is invoked.

3. Sampling, Weighting, and Resampling Mechanisms

Each iteration consists of:

1. Sampling: For every proposal ($\pi(x) \propto \exp(-f(x)-g(x)),$0), generate $\pi(x) \propto \exp(-f(x)-g(x)),$1 samples $\pi(x) \propto \exp(-f(x)-g(x)),$2.
2. Weighting: Employ deterministic mixture multiple importance sampling (MIS) weights:

$\pi(x) \propto \exp(-f(x)-g(x)),$3

with self-normalization:

$\pi(x) \propto \exp(-f(x)-g(x)),$4

1. Resampling: At every $\pi(x) \propto \exp(-f(x)-g(x)),$5 iterations, perform global resampling by drawing $\pi(x) \propto \exp(-f(x)-g(x)),$6 new anchors from the aggregated weighted sample pool; otherwise, conduct local resampling within each proposal’s batch. Anchors provide the updated positions for adaptation.
2. Adaptation: Update proposal means via the proximal Newton step, set new covariances $\pi(x) \propto \exp(-f(x)-g(x)),$7, and iterate.

The output consists of all weighted samples $\pi(x) \propto \exp(-f(x)-g(x)),$8, enabling estimation

$\pi(x) \propto \exp(-f(x)-g(x)),$9

4. Theoretical Properties

For convex $f : \mathbb{R}^d \rightarrow \mathbb{R}$0 with $f : \mathbb{R}^d \rightarrow \mathbb{R}$1-Lipschitz gradient and convex $f : \mathbb{R}^d \rightarrow \mathbb{R}$2, the proximal Newton step guarantees monotonic decrease of the local quadratic surrogate and global convergence to minimizers of $f : \mathbb{R}^d \rightarrow \mathbb{R}$3 for step parameters $f : \mathbb{R}^d \rightarrow \mathbb{R}$4. For non-convex $f : \mathbb{R}^d \rightarrow \mathbb{R}$5, provided $f : \mathbb{R}^d \rightarrow \mathbb{R}$6 satisfies the Kurdyka–Łojasiewicz property, the sequence of iterates converges to a critical point. Regarding the AIS estimation itself, as $f : \mathbb{R}^d \rightarrow \mathbb{R}$7, the self-normalized IS estimator is asymptotically unbiased, and the adaptive mechanism reduces the variance of importance weights under regularity conditions (Elvira et al., 2024).

The per-iteration complexity is dominated by three operations:

Operation	Complexity	Notes
Sampling and weighting	$f : \mathbb{R}^d \rightarrow \mathbb{R}$8	For full Gaussian proposals
Hessian inversion	$f : \mathbb{R}^d \rightarrow \mathbb{R}$9	Reduced via diagonal/approximate Hessians
Proximal solves	Depends on $g : \mathbb{R}^d \rightarrow (-\infty, +\infty]$0, often $g : \mathbb{R}^d \rightarrow (-\infty, +\infty]$1 or $g : \mathbb{R}^d \rightarrow (-\infty, +\infty]$2	Closed form for some $g : \mathbb{R}^d \rightarrow (-\infty, +\infty]$3, iterative otherwise

A plausible implication is that PNAIS scales efficiently for moderate $g : \mathbb{R}^d \rightarrow (-\infty, +\infty]$4 and/or when proximal operators are computationally simple.

5. Numerical Performance and Empirical Validation

PNAIS is validated on two canonical target families:

• Convex-constraint mixture: An equal mixture of two Gaussians on $g : \mathbb{R}^d \rightarrow (-\infty, +\infty]$5, truncated to the probability simplex $g : \mathbb{R}^d \rightarrow (-\infty, +\infty]$6. Here, $g : \mathbb{R}^d \rightarrow (-\infty, +\infty]$7 is the negative log-mixture, and $g : \mathbb{R}^d \rightarrow (-\infty, +\infty]$8 is the indicator of the simplex constraint. Metrics include relative MSE of $g : \mathbb{R}^d \rightarrow (-\infty, +\infty]$9, $\ell_1$0, and the normalizing constant $\ell_1$1.
• Sparse-prior regression: Gaussian likelihood $\ell_1$2 with mean $\ell_1$3, covariance $\ell_1$4, and Laplace (L1) prior $\ell_1$5. Here, $\ell_1$6, $\ell_1$7 is $\ell_1$8, with similar metrics.

Competing algorithms include DM-PMC (no adaptation), PNAIS-grad (first-order proximal gradient with isotropic covariance), PNAIS-rcov (robust covariance), and full PNAIS. In both settings, PNAIS achieves orders-of-magnitude improvements in MSE and superior proposal concentration relative to ablated variants and standard approaches. This suggests that second-order adaptation and the use of proximal steps confer substantive empirical benefits for non-smooth and constrained targets (Elvira et al., 2024).

6. Connections to Related Adaptive Sampling and Proximal Methods

PNAIS generalizes the principle of leveraging local geometry for proposal design, traditionally restricted to smooth targets, by incorporating proximal operations for composite models. The approach is motivated by and related to advances in deterministic mixture multiple importance sampling, stochastic optimization with composite objectives, and recent adaptive IS frameworks using robust or second-order covariance adaptation. The scheme directly addresses settings where $\ell_1$9 encodes hard constraints or induces sparsity, which arise in modern Bayesian inference and high-dimensional statistics.

7. References and Implementation Considerations

The original development, comprehensive theoretical analysis, and extensive empirical study of the Proximal Newton Adaptive Importance Sampler are detailed in

• V. Elvira & E. Chouzenoux, “A Proximal Newton Adaptive Importance Sampler,” IEEE Trans. Signal Processing, 2024 (Elvira et al., 2024).

Further context and related variance analysis can be found in prior work (see Akyıldız & Elvira 2021 as cited in (Elvira et al., 2024)). Attention must be given to selection of proposal numbers $g$0, sample size $g$1, surrogate accuracy, and practical computation of scaled proximal operators (with available closed forms or efficient numerics for structured penalties) for maximal efficiency and robustness.