Proximal Sampler for Non-smooth Distributions

Updated 3 July 2026

Proximal Sampler is a Markov chain Monte Carlo algorithm that employs proximal mappings in an alternating Gibbs sampling framework to tackle challenging distributions.
It utilizes a Restricted Gaussian Oracle to efficiently sample from non-smooth or weakly smooth potentials, ensuring convergence under logarithmic Sobolev or Poincaré conditions.
The method achieves state-of-the-art sampling complexity in convex and certain non-convex settings, outperforming traditional Langevin-based approaches.

A Proximal Sampler is a Markov chain Monte Carlo (MCMC) algorithm designed for sampling from probability distributions whose potential functions are non-smooth, weakly smooth, or even non-convex, and which may only be required to satisfy a logarithmic Sobolev or Poincaré inequality, rather than strong convexity or smoothness. The core mechanism is an alternating Gibbs sampler, closely analogous to the proximal point algorithm in convex optimization, but adapted to the sampling setting. Unlike discretizations of Langevin diffusions—which are known to degrade in the presence of non-smooth or weakly smooth targets—the proximal sampler leverages the regularization properties of proximal mappings, replacing gradient-based Markov transitions with randomized, regularized steps. This approach enables state-of-the-art sampling complexity in a broad range of settings, including strongly log-concave, log-concave, and certain classes of non-log-concave targets (Liang et al., 2022, Fan et al., 2023).

1. Algorithmic Definition and Alternating Sampling Framework

The proximal sampler operates on a joint augmented density

$\Pi(x, y) \propto \exp\left(-f(x) - \frac{1}{2\eta}\|x - y\|^2\right)$

for $x, y \in \mathbb{R}^d$ and stepsize $\eta>0$ . The Markov chain alternates the following steps:

Gaussian randomization: $y_{k} \sim \mathcal{N}(x_k, \eta I)$
Proximal update (restricted Gaussian oracle): $x_{k+1} \sim \pi^{X|Y}(x | y_k) \propto \exp\left(-f(x) - \frac{1}{2\eta}\|x-y_k\|^2\right)$

This chain has $\pi(x) \propto \exp(-f(x))$ as its stationary marginal. The key nontrivial operation is efficiently sampling from the conditional $\pi^{X|Y}(\cdot|y)$ , referred to as the Restricted Gaussian Oracle (RGO) (Liang et al., 2022, Chen et al., 2022, Fan et al., 2023).

2. Restricted Gaussian Oracle: Realizations and Rejection Sampling

Implementation of the RGO is central to the efficiency of the proximal sampler. When $f$ is convex (possibly non-smooth), the RGO can be realized via rejection sampling using a carefully constructed Gaussian proposal. The proposal is centered at a suitably chosen stationary point—typically, a (possibly approximate) minimizer of the regularized functional $f(x) + \frac{1}{2\eta}\|x-y\|^2$ . Specifically, a quadratic envelope (derived via a local quadratic or linearization argument based on the semi-smoothness of $f$ ) is used to sandwich the target, allowing control of acceptance probability. With careful selection of $x, y \in \mathbb{R}^d$ 0, the expected number of proposals is $x, y \in \mathbb{R}^d$ 1 per sample, even in high dimensions (Liang et al., 2022, Fan et al., 2023, Liang et al., 2021).

In the non-convex case, provided that the target measure $x, y \in \mathbb{R}^d$ 2 satisfies either a log-Sobolev (LSI) or Poincaré (PI) inequality and that the sum of semi-smoothness parameters and constants is appropriately controlled, rejection sampling remains efficient and maintains rigorous approximation guarantees (Liang et al., 2022, Liang et al., 2022).

3. Convergence Analysis and Complexity Guarantees

Assume that the target $x, y \in \mathbb{R}^d$ 3 satisfies LSI with constant $x, y \in \mathbb{R}^d$ 4, or PI with constant $x, y \in \mathbb{R}^d$ 5, and that $x, y \in \mathbb{R}^d$ 6 is a (possibly composite) sum of functions with semi-smoothness parameters $x, y \in \mathbb{R}^d$ 7. For suitable $x, y \in \mathbb{R}^d$ 8 (dependent on $x, y \in \mathbb{R}^d$ 9, $\eta>0$ 0, and target accuracy), the proximal sampler achieves:

LSI (in KL):

$\eta>0$ 1

iterations to reach KL divergence $\eta>0$ 2 (Liang et al., 2022, Fan et al., 2023).

PI (in $\eta>0$ 3):

$\eta>0$ 4

for $\eta>0$ 5 divergence $\eta>0$ 6.

Log-concave (convex $\eta>0$ 7):

$\eta>0$ 8

to drive KL to $\eta>0$ 9 (Liang et al., 2022). No polynomial-in- $y_{k} \sim \mathcal{N}(x_k, \eta I)$ 0 dependence in the non-convex, semi-smooth regime is required.

In the strongly log-concave and smooth case, the proximal sampler achieves optimal dimension dependence (up to logs), improving over Langevin-based methods:

$y_{k} \sim \mathcal{N}(x_k, \eta I)$ 1

for TV error $y_{k} \sim \mathcal{N}(x_k, \eta I)$ 2 (Fan et al., 2023).

4. Theoretical Extensions: General Potentials, Heavy Tails, and Inexactness

Non-convex and Composite Potentials

The proximal sampler also applies to non-convex, composite, or sum-of-semi-smooth potentials, requiring only that $y_{k} \sim \mathcal{N}(x_k, \eta I)$ 3 satisfy LSI or PI. Key to this extension is the realization that the regularized functional is "nearly" locally convex, so the rejection step remains efficient at scale (Liang et al., 2022, Liang et al., 2022).

Heavy Tail Regimes

Proximal samplers using a Gaussian kernel (standard RGO) suffer a low-accuracy barrier for heavy-tailed targets, only achieving polynomial-in- $y_{k} \sim \mathcal{N}(x_k, \eta I)$ 4 convergence rates. However, recent advances show that substituting the Gaussian with an $y_{k} \sim \mathcal{N}(x_k, \eta I)$ 5-stable heat kernel (in both forward and backward Gibbs steps) overcomes this barrier, achieving $y_{k} \sim \mathcal{N}(x_k, \eta I)$ 6 mixing times when the target and kernel tail indices match (He et al., 2024).

Inexact RGO and Accumulated Error

Controlled inexactness in RGO can be tolerated. By new concentration inequalities for semi-smooth functions, the accumulated total variation error scales linearly with the number of steps, so with per-step tolerance $y_{k} \sim \mathcal{N}(x_k, \eta I)$ 7 one achieves overall TV error $y_{k} \sim \mathcal{N}(x_k, \eta I)$ 8 (Fan et al., 2023). Similar logic holds for $y_{k} \sim \mathcal{N}(x_k, \eta I)$ 9 and $x_{k+1} \sim \pi^{X|Y}(x | y_k) \propto \exp\left(-f(x) - \frac{1}{2\eta}\|x-y_k\|^2\right)$ 0 distances.

5. Practical Realizations and Comparison with Other Methods

Efficient practical realization of the RGO is demonstrated for numerous classes:

Non-smooth convex $x_{k+1} \sim \pi^{X|Y}(x | y_k) \propto \exp\left(-f(x) - \frac{1}{2\eta}\|x-y_k\|^2\right)$ 1: Proximal-bundle methods combined with rejection sampling facilitate efficient sampling even without explicit proximal maps (Liang et al., 2021).
Constrained or composite sampling: Epigraph lifting and cutting-plane methods implement RGOs using only minimal oracles (separation for the constraint, subgradient for $x_{k+1} \sim \pi^{X|Y}(x | y_k) \propto \exp\left(-f(x) - \frac{1}{2\eta}\|x-y_k\|^2\right)$ 2) (Dang et al., 16 Feb 2026).
Stochastic or high-dimensional regimes: Stochastic gradient variants achieve optimal (or near-optimal) dependence on dimension and error (Huang et al., 2024).
Non-Euclidean geometry or manifold support: Proximal samplers generalize via the log-Laplace transform or Riemannian constructs (Guan et al., 11 Feb 2025, Gopi et al., 2023).

Compared to Langevin and Metropolis-adjusted methods, the proximal sampler strictly improves complexity in the semi-smooth or non-smooth regime and does not require strong convexity, dissipativity, or smoothness (Liang et al., 2022, Fan et al., 2023, Liang et al., 2021).