To discretize continually: Mean shift interacting particle systems for Bayesian inference

Abstract: Integration against a probability distribution given its unnormalized density is a central task in Bayesian inference and other fields. We introduce new methods for approximating such expectations with a small set of weighted samples -- i.e., a quadrature rule -- constructed via an interacting particle system that minimizes maximum mean discrepancy (MMD) to the target distribution. These methods extend the classical mean shift algorithm, as well as recent algorithms for optimal quantization of empirical distributions, to the case of continuous distributions. Crucially, our approach creates dynamics for MMD minimization that are invariant to the unknown normalizing constant; they also admit both gradient-free and gradient-informed implementations. The resulting mean shift interacting particle systems converge quickly, capture anisotropy and multi-modality, avoid mode collapse, and scale to high dimensions. We demonstrate their performance on a wide range of benchmark sampling problems, including multi-modal mixtures, Bayesian hierarchical models, PDE-constrained inverse problems, and beyond.

• The paper introduces MSIP, a novel algorithm that constructs weighted quadrature rules for unnormalized Bayesian inference using kernel embeddings and mean shift dynamics.
• It establishes a normalization-invariant framework that supports both gradient-informed and gradient-free variants, efficiently addressing multi-modal and high-dimensional challenges.
• Empirical results demonstrate that MSIP achieves lower kernel Stein discrepancy and improved mode coverage compared to traditional particle samplers like SVGD.

Mean Shift Interacting Particle Systems for Bayesian Inference: An Expert Overview

Motivation and Context

Integration against probability distributions with unnormalized densities is central to Bayesian inference, statistical physics, and numerous computational fields. Conventional particle-based and MCMC methods often suffer from slow convergence, mode collapse, or inefficiency, especially in high-dimensional or multi-modal distributions. The paper "To discretize continually: Mean shift interacting particle systems for Bayesian inference" (2605.14142) proposes a novel class of algorithms linking optimal quantization, mean-shift methods, and kernel discrepancy minimization for building weighted quadrature rules directly from an unnormalized target.

Theoretical Framework

The authors generalize the construction of quadrature rules via interacting particles so as to minimize the Maximum Mean Discrepancy (MMD) between the empirical measure and the target distribution, using both gradient-informed and gradient-free variants. Unlike classical approaches relying on uniform weights, MSIP provides a principled mechanism for assigning particle weights via optimal kernel interpolation, yielding increased flexibility and accuracy in numerical integration.

The key insight is that MSIP dynamics are invariant to the unknown normalization constant of the target density—a critical feature enabling their application to arbitrary unnormalized distributions encountered in Bayesian inference. The method formalizes particle movements as preconditioned gradient descent in the space of quadrature rules, leveraging kernel mean embeddings to capture nonlocal properties of the target.

A major theoretical contribution is the demonstration that regularization techniques required for numerical stability in kernel inversion can be reinterpreted as penalization in MMD optimization, preserving the invariance and robustness of the method.

Algorithmic Architecture and Implementation

Implementation of MSIP comprises the following workflow:

1. Kernel Embedding Estimators: For each particle configuration, kernel mean ($v_0$) and first moment ($v_1$) embeddings are estimated. Since the target density is unnormalized, these are computed using inner quadrature rules—sampling from a Gaussian centered at each particle, either via Monte Carlo, spherical quadrature, or one-point rules.
2. Weight Assignment: Optimal particle weights are determined by solving a regularized linear system involving the kernel matrix and the mean embeddings. This step is invariant to normalization.
3. Particle Update: The particle configuration is evolved via a mean-shift map, alternating between weight recomputation and an update towards kernel-weighted modes.
4. Gradient-Informed and Gradient-Free Variants: The method supports both score-informed (Stein estimator) approaches and fully gradient-free variants, broadening its applicability to scenarios where score functions are unavailable or intractable.

The authors prove that the dynamics retain invariance to normalization, and regularization preserves convergence guarantees.

Figure 1: Visualization of MSIP, SVGD, and gradient-free approaches, demonstrating modal exploration and weight allocation in a multi-modal density.

Comparisons with Existing Particle Samplers

MSIP is benchmarked against SVGD, ALDI, and consensus-based sampling (CBS) on a variety of canonical sampling problems, including highly anisotropic Gaussian mixtures, funnel distributions, and Bayesian neural network posteriors.

Figure 2: Qualitative comparison between MSIP and other samplers on a mixture of five anisotropic Gaussians, highlighting mode coverage and anisotropy adaptation.

Key empirical findings include:

• MSIP variants achieve consistently lower kernel Stein discrepancy (KSD) values across multi-modal and high-dimensional targets.
• SVGD and ALDI are susceptible to mode collapse or fail to cover modes in challenging geometries, even with extensive tuning.
• Gradient-free MSIP maintains competitive performance, particularly in domains lacking explicit score functions.
• Quantitative results exhibit strong robustness across geometries and problem dimensions.

Case Study: Bayesian Neural Network Posteriors

An additional experiment shows rapid convergence of MSIP (both Fredholm and gradient-informed variants) on the log-likelihood of Bayesian neural network ensembles, significantly outperforming SVGD.

Figure 3: Convergence of ensemble log-likelihood for TwoMoons classification, demonstrating MSIP’s efficiency and stability.

Modal Coverage and Anisotropy

Visualizations on canonical synthetic distributions further evidence the modal coverage and anisotropy adaptation capabilities of MSIP compared to standard samplers.

Figure 4: MSIP particle configuration on a 2D Gaussian mixture: all modes are covered and anisotropy is respected.

Figure 5: MSIP adaptation to non-Gaussian (Joker) distributions.

Figure 6: Coverage by MSIP in strongly anisotropic funnel distributions.

Figure 7: MSIP recovery of modes in the Himmelblau distribution.

Practical and Theoretical Implications

MSIP provides a scalable, normalization-invariant quadrature construction for arbitrary unnormalized densities, capable of robust inference on complex posteriors. For practical Bayesian computation, MSIP offers an attractive alternative to classical MCMC and particle-based samplers, especially in settings where efficiency, modality, or shape accuracy are crucial.

Theoretical implications include:

• MMD minimization for unnormalized densities, previously considered infeasible, is now tractable via kernel mean embedding estimation.
• Robustness to high-dimensionality and non-Gaussianity, achieved by preconditioned particle interactions.
• Potential for adaptive kernel and inner quadrature selection, offering avenues for further algorithmic improvement.

Outlook and Future Directions

The paper outlines several avenues for future research:

• Rigorous convergence and approximation analysis for the nonconvex optimization at the heart of MSIP.
• Adaptive selection strategies for inner quadrature rules and kernel hyperparameters.
• Application of MSIP dynamics to broader generative modeling and uncertainty quantification domains.

Conclusion

"To discretize continually: Mean shift interacting particle systems for Bayesian inference" (2605.14142) introduces a mathematically principled algorithmic framework for weighted quadrature construction via interacting particles, addressing fundamental obstacles in Bayesian integration for unnormalized densities. Both theoretical and empirical evidence demonstrates substantially increased robustness, systematic modal coverage, and scalability compared to established kernel and consensus-based methods. MSIP opens new directions for measure transport, kernel-based inference, and numerical integration in high-dimensional and complex geometries.