Mean Shift Interacting Particle Systems

• Mean Shift Interacting Particle Systems are a framework that approximates target distributions with a finite weighted sum of Dirac masses by minimizing the maximum mean discrepancy.
• The approach employs gradient flows, fixed-point updates, and regularized kernel quadrature to optimize both particle locations and weights for accurate numerical integration.
• It extends classical mean shift and Lloyd-type algorithms, providing robust solutions for multi-modal, anisotropic, and high-dimensional Bayesian inference challenges.

Mean Shift Interacting Particle Systems (MSIP) are interacting-particle constructions for approximating a target probability law by a finitely supported weighted measure, typically with the objective of minimizing maximum mean discrepancy (MMD) and thereby producing a kernel quadrature rule or a weighted quantization. In the recent literature, MSIP is formulated as a multi-particle extension of classical mean shift: instead of locating a single mode of a kernel density estimate, one seeks a configuration of multiple interacting particles and weights that jointly minimize a kernel discrepancy to the target distribution. In Bayesian inference, the same framework is developed for continuous targets known only through an unnormalized density, and in amortized form it is learned as a map from observations and seed posterior samples to weighted quadrature nodes (Belhadji et al., 14 Feb 2025, Belhadji et al., 13 May 2026). The term should be distinguished from the broader literature on mean-field or McKean–Vlasov interacting particle systems, where each particle depends on the empirical law or its mean-field limit but no classical mean-shift mechanism is necessarily involved (Maestra et al., 2020).

1. Variational formulation and core objective

The basic MSIP problem is to approximate a target distribution $\pi \in \mathcal P(\mathcal X)$ by a finite weighted sum of Dirac masses

$\mu = \sum_{i=1}^M w_i \delta_{y_i},$

where $y_i$ are particle locations and $w_i$ are weights. The discrepancy is measured by the RKHS-induced maximum mean discrepancy,

$$\mathrm{MMD}(\mu,\nu) := \sup_{\|f\|_{\mathcal H}\le 1} \left| \int_{\mathcal X} f(x)\,d\mu(x) - \int_{\mathcal X} f(x)\,d\nu(x) \right|,$$

and the squared objective is

$F(\mu) := \frac12\,\mathrm{MMD}(\mu,\pi)^2.$

For a discrete weighted rule, minimizing MMD is equivalent to minimizing worst-case integration error over the RKHS unit ball, which gives MSIP an immediate quadrature interpretation (Belhadji et al., 14 Feb 2025).

In the continuous-target Bayesian setting, the same idea is written as

$$F^\pi(\mathbf w, Y) := \frac12 \mathrm{MMD}^2\!\left(\pi,\sum_{i=1}^M w_i\delta_{y_i}\right),$$

with $Y=\{y_i\}_{i=1}^M$. The target expectation

$$\mathbb{E}_\pi[f(Y)] = \int_{\mathcal X} f(y)\,d\pi(y)$$

is then approximated by

$$\int_{\mathcal X} f(y)\,d\pi(y)\approx \sum_{i=1}^M w_i f(y_i).$$

A central feature of this formulation is that the particle system is explicitly designed for nonuniform weighted quadrature, rather than uniform-weight sampling, because the accuracy criterion is numerical integration rather than sample generation alone (Belhadji et al., 13 May 2026).

The kernelized structure is built from the kernel matrix and kernel mean embeddings. With squared exponential kernel

$\mu = \sum_{i=1}^M w_i \delta_{y_i},$0

one introduces

$\mu = \sum_{i=1}^M w_i \delta_{y_i},$1

and $\mu = \sum_{i=1}^M w_i \delta_{y_i},$2 with entries $\mu = \sum_{i=1}^M w_i \delta_{y_i},$3. For fixed locations, the optimal weights satisfy

$\mu = \sum_{i=1}^M w_i \delta_{y_i},$4

which the literature identifies with the standard optimal kernel quadrature formula (Belhadji et al., 14 Feb 2025).

2. Particle dynamics, gradient flows, and fixed-point maps

One derivation of MSIP starts from a Wasserstein–Fisher–Rao (WFR) gradient flow for the MMD functional. In that geometry, the measure-valued evolution is

$\mu = \sum_{i=1}^M w_i \delta_{y_i},$5

where the reaction term $\mu = \sum_{i=1}^M w_i \delta_{y_i},$6 allows mass creation and destruction. This is important because the formulation does not force the measure to remain normalized at every step, which is compatible with weighted quantization by particles with variable weights (Belhadji et al., 14 Feb 2025).

Representing the evolving measure by finitely many Dirac masses,

$\mu = \sum_{i=1}^M w_i \delta_{y_i},$7

yields the interacting ODE system

$\mu = \sum_{i=1}^M w_i \delta_{y_i},$8

With $\mu = \sum_{i=1}^M w_i \delta_{y_i},$9, this becomes

$y_i$0

$y_i$1

Steady states satisfy the kernel quadrature optimality relation

$y_i$2

together with the corresponding gradient condition. Under the kernel assumption

$y_i$3

these steady-state equations can be rearranged into a mean-shift-like fixed-point map, called mean shift interacting particles (MSIP). For the squared exponential kernel, the resulting update simplifies to

$y_i$4

which is the cleanest form of the algorithm (Belhadji et al., 14 Feb 2025).

In the continuous-target Bayesian formulation, the same fixed-point logic is regularized through

$y_i$5

with weight solve

$y_i$6

and map

$y_i$7

The damped iteration is

$y_i$8

and the regularized reduced objective is

$y_i$9

A key theorem gives the explicit gradient

$w_i$0

so the fixed-point update is a preconditioned gradient descent on the reduced MMD objective (Belhadji et al., 13 May 2026).

3. Continuous-target Bayesian inference and normalization invariance

For Bayesian inference, MSIP is designed for targets known only through an unnormalized density $w_i$1, with

$w_i$2

The central algorithmic challenge is that the kernel means $w_i$3 and $w_i$4 are generally unavailable in closed form. The continuous-target formulation therefore constructs them through an inner quadrature rule

$w_i$5

where $w_i$6 is standard Gaussian. Using the Gaussian convolution identity,

$w_i$7

one obtains

$w_i$8

For the first moment,

$w_i$9

and also

$$\mathrm{MMD}(\mu,\nu) := \sup_{\|f\|_{\mathcal H}\le 1} \left| \int_{\mathcal X} f(x)\,d\mu(x) - \int_{\mathcal X} f(x)\,d\nu(x) \right|,$$0

which leads to two principal estimators: $$\mathrm{MMD}(\mu,\nu) := \sup_{\|f\|_{\mathcal H}\le 1} \left| \int_{\mathcal X} f(x)\,d\mu(x) - \int_{\mathcal X} f(x)\,d\nu(x) \right|,$$1 and the Stein or gradient-informed form

$$\mathrm{MMD}(\mu,\nu) := \sup_{\|f\|_{\mathcal H}\le 1} \left| \int_{\mathcal X} f(x)\,d\mu(x) - \int_{\mathcal X} f(x)\,d\nu(x) \right|,$$2

The one-point special case $$\mathrm{MMD}(\mu,\nu) := \sup_{\|f\|_{\mathcal H}\le 1} \left| \int_{\mathcal X} f(x)\,d\mu(x) - \int_{\mathcal X} f(x)\,d\nu(x) \right|,$$3, $$\mathrm{MMD}(\mu,\nu) := \sup_{\|f\|_{\mathcal H}\le 1} \left| \int_{\mathcal X} f(x)\,d\mu(x) - \int_{\mathcal X} f(x)\,d\nu(x) \right|,$$4, gives

$$\mathrm{MMD}(\mu,\nu) := \sup_{\|f\|_{\mathcal H}\le 1} \left| \int_{\mathcal X} f(x)\,d\mu(x) - \int_{\mathcal X} f(x)\,d\nu(x) \right|,$$5

This yields the named variants MSIP-F, MSIP-GI-1, MSIP-GI-10, and MSIP-GF; a hybrid estimator is also used (Belhadji et al., 13 May 2026).

A defining property of this construction is invariance to the unknown normalizing constant. If $$\mathrm{MMD}(\mu,\nu) := \sup_{\|f\|_{\mathcal H}\le 1} \left| \int_{\mathcal X} f(x)\,d\mu(x) - \int_{\mathcal X} f(x)\,d\nu(x) \right|,$$6 for $$\mathrm{MMD}(\mu,\nu) := \sup_{\|f\|_{\mathcal H}\le 1} \left| \int_{\mathcal X} f(x)\,d\mu(x) - \int_{\mathcal X} f(x)\,d\nu(x) \right|,$$7, then

$$\mathrm{MMD}(\mu,\nu) := \sup_{\|f\|_{\mathcal H}\le 1} \left| \int_{\mathcal X} f(x)\,d\mu(x) - \int_{\mathcal X} f(x)\,d\nu(x) \right|,$$8

so

$$\mathrm{MMD}(\mu,\nu) := \sup_{\|f\|_{\mathcal H}\le 1} \left| \int_{\mathcal X} f(x)\,d\mu(x) - \int_{\mathcal X} f(x)\,d\nu(x) \right|,$$9

Likewise,

$F(\mu) := \frac12\,\mathrm{MMD}(\mu,\pi)^2.$0

and therefore

$F(\mu) := \frac12\,\mathrm{MMD}(\mu,\pi)^2.$1

This invariance is one of the central reasons MSIP is adapted to Bayesian posteriors specified only up to evidence (Belhadji et al., 13 May 2026).

4. Amortized MSIP and sample-only posterior access

A later development replaces per-observation MSIP optimization by a learned, amortized construction. In that formulation, the quadrature measure is

$F(\mu) := \frac12\,\mathrm{MMD}(\mu,\pi)^2.$2

and the discrepancy to the posterior $F(\mu) := \frac12\,\mathrm{MMD}(\mu,\pi)^2.$3 is again

$F(\mu) := \frac12\,\mathrm{MMD}(\mu,\pi)^2.$4

The classical mean-shift interacting-particle update is characterized by

$F(\mu) := \frac12\,\mathrm{MMD}(\mu,\pi)^2.$5

but the expensive step is that this optimization must be solved separately for each new observation. In the score-based version, it also requires $F(\mu) := \frac12\,\mathrm{MMD}(\mu,\pi)^2.$6, which can be prohibitive in inverse problems with PDE forward models (Siahkoohi, 14 Jun 2026).

The amortized method learns a set-equivariant map $F(\mu) := \frac12\,\mathrm{MMD}(\mu,\pi)^2.$7 that predicts particle displacements from an observation, a seed set of posterior samples, and the requested node budget: $F(\mu) := \frac12\,\mathrm{MMD}(\mu,\pi)^2.$8 The network therefore learns node locations rather than weights or a density. Inputs are an observation $F(\mu) := \frac12\,\mathrm{MMD}(\mu,\pi)^2.$9, seed samples $$F^\pi(\mathbf w, Y) := \frac12 \mathrm{MMD}^2\!\left(\pi,\sum_{i=1}^M w_i\delta_{y_i}\right),$$0, and $$F^\pi(\mathbf w, Y) := \frac12 \mathrm{MMD}^2\!\left(\pi,\sum_{i=1}^M w_i\delta_{y_i}\right),$$1; outputs are the displacements $$F^\pi(\mathbf w, Y) := \frac12 \mathrm{MMD}^2\!\left(\pi,\sum_{i=1}^M w_i\delta_{y_i}\right),$$2. The weights are still solved in closed form at the emitted nodes: $$F^\pi(\mathbf w, Y) := \frac12 \mathrm{MMD}^2\!\left(\pi,\sum_{i=1}^M w_i\delta_{y_i}\right),$$3 with $$F^\pi(\mathbf w, Y) := \frac12 \mathrm{MMD}^2\!\left(\pi,\sum_{i=1}^M w_i\delta_{y_i}\right),$$4 chosen so that $$F^\pi(\mathbf w, Y) := \frac12 \mathrm{MMD}^2\!\left(\pi,\sum_{i=1}^M w_i\delta_{y_i}\right),$$5. The training loss is MMD-based and uses only samples from a reference posterior; the method does not evaluate density or score at inference (Siahkoohi, 14 Jun 2026).

A central theorem states that for fixed nodes, optimal reweighting cannot be worse than equal Monte Carlo weights: $$F^\pi(\mathbf w, Y) := \frac12 \mathrm{MMD}^2\!\left(\pi,\sum_{i=1}^M w_i\delta_{y_i}\right),$$6 This establishes the “provably no worse than equal weights” property for reweighting. The larger empirical gain is attributed to moving the particles, since reweighting alone cannot place nodes on modes or ridges absent from the seed set (Siahkoohi, 14 Jun 2026).

To address the high-dimensional wall, the amortized formulation uses a posterior-whitened Mahalanobis kernel

$$F^\pi(\mathbf w, Y) := \frac12 \mathrm{MMD}^2\!\left(\pi,\sum_{i=1}^M w_i\delta_{y_i}\right),$$7

with dimension-aware bandwidth

$$F^\pi(\mathbf w, Y) := \frac12 \mathrm{MMD}^2\!\left(\pi,\sum_{i=1}^M w_i\delta_{y_i}\right),$$8

In whitened coordinates this is the usual squared-exponential kernel, and the stated purpose is to preserve useful Gram-matrix geometry in high dimension (Siahkoohi, 14 Jun 2026).

5. Relation to classical mean shift, Lloyd-type quantization, and broader mean-field systems

MSIP is explicitly presented as an extension of classical mean shift. Classical mean shift identifies critical points of a kernel density estimate through the fixed-point relation

$$F^\pi(\mathbf w, Y) := \frac12 \mathrm{MMD}^2\!\left(\pi,\sum_{i=1}^M w_i\delta_{y_i}\right),$$9

The MSIP literature shows that when $Y=\{y_i\}_{i=1}^M$0, the interacting-particle construction reduces to this classical mean-shift relation. The difference is that MSIP uses multiple particles, optimized weights, and an MMD objective; classical mean shift is therefore the one-particle, no-interaction special case, whereas MSIP seeks a weighted multi-particle configuration minimizing MMD to a target law (Belhadji et al., 14 Feb 2025).

The same papers also position MSIP relative to Lloyd’s algorithm. In $Y=\{y_i\}_{i=1}^M$1-quantization, Lloyd’s method is a fixed-point or preconditioned gradient descent on a Voronoi-based objective. MSIP is described as a relaxation of Lloyd’s algorithm for clustering, but with weighted particles, kernel interactions, and an MMD objective rather than Voronoi tessellations. This suggests a conceptual bridge between clustering, kernel quadrature, and weighted quantization (Belhadji et al., 14 Feb 2025).

A recurring source of confusion is the phrase “mean-field interacting particle systems.” In the broader McKean–Vlasov literature, one studies systems such as

$Y=\{y_i\}_{i=1}^M$2

whose large-$Y=\{y_i\}_{i=1}^M$3 limit satisfies a nonlinear Fokker–Planck equation. That literature concerns interacting particle systems and mean-field limits, and it may include Vlasov-type convolution drifts $Y=\{y_i\}_{i=1}^M$4, but it does not address the classical mean-shift algorithm or MMD-based weighted quadrature as such (Maestra et al., 2020). Related ensemble-based samplers, including ALDI, CBS, and stochastic SVGD, also use particle interactions derived from ensemble statistics or kernelized drifts, yet their goal is sampling with correct invariant law rather than MSIP-style deterministic quadrature; Metropolization is introduced there to remove discretization and finite-ensemble bias (Sprungk et al., 2023).

A plausible implication is that “mean shift interacting particle systems” should be interpreted narrowly when referring to the MMD-minimizing quantization and quadrature methods, and only broadly when used as a loose label for mean-field particle dynamics more generally.

6. Reported behavior, application domains, and technical themes

The continuous-target MSIP paper reports experiments on synthetic multimodal and anisotropic targets, including 2D and higher-dimensional Gaussian mixtures, a “Joker” distribution, funnel distributions, and the Himmelblau density, as well as Bayesian benchmarks including Schools, COVID, CovType hierarchical logistic regression, ABPDE, and a TwoMoons Bayesian neural network with 201 parameters. The reported findings are that MSIP variants consistently achieve the lowest or near-lowest KSD across many problems; on five-component anisotropic Gaussian mixtures, MSIP-Fredholm captures all modes while SVGD and ALDI exhibit mode collapse; and MSIP is especially strong when the target is multi-modal, anisotropic, or high-dimensional. The same work also states that MSIP-Fredholm and MSIP-GI-1 have roughly the same density-evaluation cost as SVGD, while larger inner quadrature rules cost more but remain parallelizable (Belhadji et al., 13 May 2026).

The amortized formulation is reported across closed-form, sampled, learned, and physics-based posteriors, including a thousand-coefficient groundwater field. Its stated result is that the method integrates more accurately than the same number of samples at every budget, while remaining sample-based at inference and generalizing to unseen observations and integrands. The posterior-whitened, dimension-aware kernel is introduced specifically to remove the high-dimensional wall, and the method is described as a Pareto improvement on Monte-Carlo integration rather than a competitor to drawing more samples (Siahkoohi, 14 Jun 2026).

Several technical themes recur across the literature. First, kernel-matrix conditioning motivates explicit regularization through $Y=\{y_i\}_{i=1}^M$5. Second, access to $Y=\{y_i\}_{i=1}^M$6 and $Y=\{y_i\}_{i=1}^M$7 is the main interface between the abstract MMD objective and implementable algorithms, whether by Gaussian inner quadrature, score-based Stein identities, or posterior samples alone. Third, the distinction between reweighting and relocation is structurally important: the former is accompanied by a theorem, while the latter is presented as the principal empirical source of improvement. Finally, the literature consistently emphasizes that MSIP is not merely a mode-seeking method. Its stated objective is coverage of the full distribution through a weighted particle configuration that minimizes a global RKHS discrepancy, which explains its connection to kernel quadrature as much as to mean shift.