Weighted Support Points from Random Measures: An Interpretable Alternative for Generative Modeling

TL;DR

This paper introduces a novel generative modeling approach using random weighted support points inspired by Bayesian methods, offering an interpretable, scalable, and efficient alternative to neural network-based models like GANs and DDPMs.

Contribution

It develops a new framework for generative modeling based on random weighted support points with a theoretical foundation and an efficient optimization algorithm, avoiding neural networks.

Findings

Produces high-quality, diverse samples on MNIST and CelebA-HQ

Achieves comparable results at lower computational cost

Generates interpolative samples that preserve data structure

Abstract

Support points summarize a large dataset through a smaller set of representative points that can be used for data operations, such as Monte Carlo integration, without requiring access to the full dataset. In this sense, support points offer a compact yet informative representation of the original data. We build on this idea to introduce a generative modeling framework based on random weighted support points, where the randomness arises from a weighting scheme inspired by the Dirichlet process and the Bayesian bootstrap. The proposed method generates diverse and interpretable sample sets from a fixed dataset, without relying on probabilistic modeling assumptions or neural network architectures. We present the theoretical formulation of the method and develop an efficient optimization algorithm based on the Convex--Concave Procedure (CCP). Empirical results on the MNIST and CelebA-HQ datasets show that our approach produces high-quality and diverse outputs at a fraction of the computational cost of black-box alternatives such as Generative Adversarial Networks (GANs) or Denoising Diffusion Probabilistic Models (DDPMs). These results suggest that random weighted support points offer a principled, scalable, and interpretable alternative for generative modeling. A key feature is their ability to produce genuinely interpolative samples that preserve underlying data structure.

Equations

$$\underset{G}{\min }\underset{D}{\max }V(D,G)$$

$$+{\mathbb{E}}_{\mathbf{z}\sim p_z}[\log (1-D(G(\mathbf{z})))].$$

$$D_G^{\ast }(\mathbf{x})=\frac{p_{\text{data}}(\mathbf{x})}{p_{\text{data}}(\mathbf{x})+p_g(\mathbf{x})},$$

$$\underset{G}{\min }V(D_G^{\ast },G)=-\log 4+2\cdot \text{JSD}(p_{\text{data}}\Vert p_g).$$

$$\underset{G}{\max }{\mathbb{E}}_{\mathbf{z}\sim p_z}[\log D(G(\mathbf{z}))],$$

$${\mathbf{x}}_t=\sqrt{1-{\beta }_t}{\mathbf{x}}_{t-1}+\sqrt{{\beta }_t}{\mathbf{ϵ}}_t,{\mathbf{ϵ}}_t\sim \mathcal{N}(0,\mathbf{I}),t=1,\dots ,T,$$

$${\mathbf{x}}_t=\sqrt{{\bar{\alpha }}_t}{\mathbf{x}}_0+\sqrt{1-{\bar{\alpha }}_t}\mathbf{ϵ},\mathbf{ϵ}\sim \mathcal{N}(0,\mathbf{I}),$$

$$p_{\theta }({\mathbf{x}}_{t-1}|{\mathbf{x}}_t)=\mathcal{N}({\mathbf{x}}_{t-1};{\mathbf{\mu }}_{\theta }({\mathbf{x}}_t,t),{\sigma }_t^2\mathbf{I}),$$

$$q({\mathbf{x}}_{t-1}|{\mathbf{x}}_t,{\mathbf{x}}_0)=\mathcal{N}({\mathbf{x}}_{t-1};{\tilde{\mathbf{\mu }}}_t,{\tilde{\beta }}_t\mathbf{I}),$$

$${\tilde{\mathbf{\mu }}}_t({\mathbf{x}}_t,{\mathbf{x}}_0)=\frac{\sqrt{{\bar{\alpha }}_{t-1}}{\beta }_t}{1-{\bar{\alpha }}_t}{\mathbf{x}}_0+\frac{\sqrt{{\alpha }_t}(1-{\bar{\alpha }}_{t-1})}{1-{\bar{\alpha }}_t}{\mathbf{x}}_t.$$

$${\mathbf{x}}_0=\frac{1}{\sqrt{{\bar{\alpha }}_t}}\left({\mathbf{x}}_t-\sqrt{1-{\bar{\alpha }}_t}\mathbf{ϵ}\right),$$

$${\mathbf{\mu }}_{\theta }({\mathbf{x}}_t,t)=\frac{1}{\sqrt{{\alpha }_t}}\left({\mathbf{x}}_t-\frac{{\beta }_t}{\sqrt{1-{\bar{\alpha }}_t}}{\mathbf{ϵ}}_{\theta }({\mathbf{x}}_t,t)\right).$$

$$L_{\text{simple}}(\theta )={\mathbb{E}}_{t,{\mathbf{x}}_0,\mathbf{ϵ}}\left[{\Vert \mathbf{ϵ}-{\mathbf{ϵ}}_{\theta }({\mathbf{x}}_t,t)\Vert }^2\right],$$

$${\mathbf{x}}_{t-1}={\mathbf{\mu }}_{\theta }({\mathbf{x}}_t,t)+{\sigma }_t\mathbf{z},\mathbf{z}\sim \mathcal{N}(0,\mathbf{I}).$$

$$E(F,G)=2\mathbb{E}\Vert \mathbf{X}-\mathbf{Y}{\Vert }_2-\mathbb{E}\Vert \mathbf{X}-{\mathbf{X}}^{\prime }{\Vert }_2-\mathbb{E}\Vert \mathbf{Y}-{\mathbf{Y}}^{\prime }{\Vert }_2,$$

$$\mathbf{S}=\underset{\mathbf{A}\in {\mathbb{R}}^{d\times n}}{\mathrm{argmin}}E(F,{\hat{F}}_n),$$

$${\hat{F}}_n(\mathbf{x})=\frac{1}{n}\sum _{i=1}^n\mathbb{I}\{{\mathbf{x}}_i\le \mathbf{x}\}.$$

$$E(F,{\hat{F}}_n)$$

$$-\frac{1}{n^2}\sum _{i=1}^n\sum _{j=1}^n\Vert {\mathbf{x}}_i-{\mathbf{x}}_j{\Vert }_2$$

$$-{\mathbb{E}}_{\mathbf{Y},{\mathbf{Y}}^{\prime }}\left[\Vert \mathbf{Y}-{\mathbf{Y}}^{\prime }{\Vert }_2\right].$$

$${\mathbb{E}}_{\mathbf{X}\sim {\hat{F}}_n}[h(\mathbf{X})]=\int h(\mathbf{x})d{\hat{F}}_n(\mathbf{x})=\frac{1}{n}\sum _{i=1}^{n}h({\mathbf{x}}_i).$$

$$MC(\mathbf{A};\mathbf{P})$$

$$-\underset{\text{repulsion term}}{\underbrace{\frac{1}{n^2}\sum _{i=1}^n\sum _{j=1}^n\Vert {\mathbf{x}}_i-{\mathbf{x}}_j{\Vert }_2}},$$

$$MC(\mathbf{A};\mathbf{P})=f_{\text{att}}(\mathbf{A};\mathbf{P})-f_{\text{rep}}(\mathbf{A}),$$

$${\mathbf{x}}_i^{(l+1)}\leftarrow M_i({\mathbf{A}}^{(l)};\mathbf{P})\text{(see equation eqrefeq1)}$$

$${\mathbf{x}}_i$$

$$=\frac{1}{q({\mathbf{x}}_i^{\prime },\mathbf{P})}\left(\sum _{m=1}^N\frac{{\mathbf{y}}_m}{\Vert {\mathbf{x}}_i^{\prime }-{\mathbf{y}}_m{\Vert }_2}+\frac{N}{n}\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n\frac{{\mathbf{x}}_i^{\prime }-{\mathbf{x}}_j^{\prime }}{\Vert {\mathbf{x}}_i^{\prime }-{\mathbf{x}}_j^{\prime }{\Vert }_2}\right),$$

$$q({\mathbf{x}}_i^{\prime },\mathbf{P})=\sum _{m=1}^N\frac{1}{\Vert {\mathbf{x}}_i^{\prime }-{\mathbf{y}}_m{\Vert }_2}.$$

$${\tilde{F}}_N(\mathbf{y})=\sum _{m=1}^N{w}_m\mathbb{I}\{{\mathbf{y}}_m\le \mathbf{y}\},$$

$$\mathbf{S}(\mathbf{w})=\underset{\mathbf{A}\in {\mathbb{R}}^{d\times n}}{\mathrm{argmin}}E({\tilde{F}}_N,{\hat{F}}_n),$$

$$MC(\mathbf{A};\mathbf{P},\mathbf{w})$$

$$-\frac{1}{n^2}\sum _{i=1}^n\sum _{j=1}^n\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert .$$

$$\frac{\mathrm{d}}{\mathrm{d}{\mathbf{x}}_i}\left(\frac{2}{n}\sum _{k=1}^n\sum _{m=1}^N{w}_m\Vert {\mathbf{x}}_k-{\mathbf{y}}_m\Vert \right)=\frac{2}{n}\sum _{m=1}^N{w}_m\frac{{\mathbf{x}}_i-{\mathbf{y}}_m}{\Vert {\mathbf{x}}_i-{\mathbf{y}}_m\Vert }.$$

Full text

Keywords: GANs, DDPMs, Energy Distance, Dirichlet Process, Truncated Measure.

1 Introduction

Generative models are the cornerstone of modern machine learning. Their goal is to capture the underlying probability distribution of observed data in order to generate new samples that resemble the original dataset (Ho et al.,, 2020; Kingma and Welling,, 2013). This ability to create realistic and diverse data, particularly in complex domains such as images, has opened up transformative applications in computer vision, digital art, data augmentation, and scientific simulation (Karras et al.,, 2019; Rombach et al.,, 2022; Dhariwal and Nichol,, 2021; Brock et al.,, 2019; Karras et al.,, 2020). Advances in deep learning have played a central role in this progress, enabling models to learn intricate data structures that were previously out of reach.

Two of the most prominent families of deep generative models are Generative Adversarial Networks (GANs) (Goodfellow et al.,, 2014; Radford et al.,, 2015) and Denoising Diffusion Probabilistic Models (DDPMs) (Ho et al.,, 2020; Sohl-Dickstein et al.,, 2015; Song and Ermon,, 2019). GANs are based on a game-theoretic setup, where a generator learns to produce realistic samples while a discriminator tries to distinguish them from real data (Goodfellow et al.,, 2014). In contrast, DDPMs are inspired by physical diffusion processes and rely on a two-stage procedure: a forward process that gradually corrupts the data by adding noise, and a learned reverse process that reconstructs samples from noise (Ho et al.,, 2020; Sohl-Dickstein et al.,, 2015). Both approaches have achieved remarkable results in image synthesis, producing high-resolution and visually compelling samples in a variety of applications (Ho et al.,, 2020; Karras et al.,, 2019).

Despite their impressive performance, established generative models like GANs and DDPMs face important limitations. At the same time, they share a reliance on deep neural networks to learn complex data distributions, while differing significantly in their architectures, training procedures, and computational trade-offs.

GANs consist of two neural networks, a generator and a discriminator, that are trained simultaneously in an adversarial setting. The generator learns to produce realistic samples, while the discriminator tries to distinguish real from synthetic data (Goodfellow et al.,, 2014). Although GANs can generate high-quality images efficiently once trained, the training process is notoriously unstable. It often suffers from issues such as mode collapse, where the generator produces a narrow set of outputs, and requires a delicate balance between the two networks to converge successfully (Goodfellow,, 2016; Salimans et al.,, 2016).

In contrast, DDPMs typically involve a single, large neural network, often based on a U-Net architecture with millions of parameters, trained to reverse a noise corruption process. Generating a single sample requires hundreds or even thousands of iterative denoising steps, each invoking the network. This makes both training and inference computationally intensive. Moreover, the multi-step generative process complicates interpretation and renders the model a black box, despite its strong empirical performance (Nichol and Dhariwal,, 2021; Song et al.,, 2020).

These limitations encourage exploring alternative generative frameworks that provide simpler training, reduced computational costs, or fresh theoretical perspectives for producing high-quality data.

This work explores the use of support points for generative modeling, proposing a mathematically grounded, nonparametric alternative to black-box neural methods such as GANs and diffusion models. Originally introduced by Mak and Joseph, (2018), support points are defined as a set of representative locations that minimize the energy distance between a target distribution $F$ and its discrete approximation $\widehat{F}_{n}$ for a chosen sample size $n$ (Székely and Rizzo,, 2004, 2013). The energy distance is a geometric divergence widely used in nonparametric hypothesis testing and distributional approximation. While support points have proven effective for tasks such as distribution compression, uncertainty quantification, and MCMC summarization (Mak and Joseph,, 2018; Riachy et al.,, 2019; Mak and Joseph,, 2017; Joseph et al.,, 2015, 2019), their potential for generative modeling, particularly in high-dimensional settings, remains un-explored.

The central contribution of this work is a framework for weighted support points driven by a finite random measure designed to balance representativeness and diversity. Instead of reweighting all $N$ observations under Dirichlet–process or Bayesian–bootstrap schemes (Ferguson,, 1973; Rubin,, 1981), we construct a reweighted empirical measure by first selecting a random subset of the reference data and then assigning exchangeable weights from a symmetric Dirichlet distribution. The resulting measure remains centered at the empirical distribution and has tunable dispersion, encouraging genuine run-to-run variability without collapsing mass onto a few atoms.

Given this random measure, we compute support points of size $n$ by minimizing the energy distance between it and the empirical distribution of the candidate set. This leads to a stochastic optimization whose solutions differ across realizations, yielding diverse yet data-aligned configurations from the same reference dataset. The subsetting rule and the calibration of weight dispersion are detailed in the body of the paper. Empirical results on MNIST and CelebA-HQ datasets, including comparisons with GANs and DDPMs, show that the proposed method produces high-quality, diverse outputs at a fraction of the computational cost, offering a principled, interpretable, and scalable alternative for generative modeling.

2 Background: Established Generative Models

This section provides a mathematical overview of two prominent generative modeling frameworks: Generative Adversarial Networks (GANs) and Denoising Diffusion Probabilistic Models (DDPMs). Understanding their core principles is essential for contextualizing and comparing the support points approach.

2.1 Generative Adversarial Networks (GANs)

GANs, introduced by Goodfellow et al., (2014), represent a milestone in generative modeling by casting the learning process as a two-player minimax game between a generator and a discriminator.

Training unfolds as a competitive interaction: the generator attempts to fool the discriminator, while the discriminator aims to correctly distinguish real from synthetic samples. The generator, denoted $G(\mathbf{z};\theta_{g})$, maps latent noise vectors $\mathbf{z}\sim p_{z}(\mathbf{z})$ (typically sampled from a standard normal or uniform distribution) to synthetic data samples $\mathbf{x}_{g}=G(\mathbf{z})$. The discriminator, $D(\mathbf{x};\theta_{d})$, receives a sample $\mathbf{x}$ and outputs the estimated probability that $\mathbf{x}$ originates from the real data distribution $p_{\text{data}}$, rather than from the generator’s distribution $p_{g}$.

The networks are trained by solving the following minimax optimization problem:

[TABLE]

This adversarial setup incentivizes the generator to produce increasingly realistic data. At equilibrium, the generator perfectly replicates the true data distribution ($p_{g}=p_{\text{data}}$), and the discriminator is maximally uncertain, yielding $D(\mathbf{x})=1/2$ for all $\mathbf{x}$.

For a fixed generator $G$, the optimal discriminator is

[TABLE]

which, when substituted into the value function $V(D,G)$, leads to an objective for the generator equivalent to minimizing the Jensen–Shannon divergence:

[TABLE]

In practice, training GANs can be challenging. The original loss in Equation (1) often results in vanishing gradients for the generator, particularly when $D(G(\mathbf{z}))\approx 0$. To address this, a widely adopted modification is to train the generator using the non-saturating heuristic:

[TABLE]

which yields stronger and more informative gradients during early stages of training.

GANs have been successfully applied in tasks such as image generation, domain translation, and data augmentation. Numerous extensions address training instability, including Deep Convolutional GANs (DCGAN) (Radford et al.,, 2015), Wasserstein GANs (WGAN) (Arjovsky et al.,, 2017), among others.

2.2 Denoising Diffusion Probabilistic Models (DDPMs)

DDPMs have become a leading generative modeling technique, achieving remarkable performance in image, audio, and multimodal generation tasks (Ho et al.,, 2020; Nichol and Dhariwal,, 2021; Dhariwal and Nichol,, 2021). Their central idea is conceptually simple but mathematically rich: progressively corrupt the training data by adding Gaussian noise, and then learn to reverse this process step by step to generate new data.

The diffusion process is structured as a pair of Markov chains: a fixed forward process that adds noise to data, and a learned reverse process that denoises this input back to the data distribution.

2.2.1 Forward Process

Let $\boldsymbol{x}_{0}\sim q(\boldsymbol{x}_{0})$ be a data point. The forward process generates a sequence $\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{T}$ by recursively adding Gaussian noise:

[TABLE]

where $\{\beta_{t}\}_{t=1}^{T}$ is a predefined variance schedule. Typically, $\beta_{t}$ increases linearly from small values (e.g., $10^{-4}$) to moderate values (e.g., $0.02$) (Ho et al.,, 2020).

A key result is that $\boldsymbol{x}_{t}$ can be sampled directly from $\boldsymbol{x}_{0}$ without computing intermediate steps:

[TABLE]

where $\alpha_{t}=1-\beta_{t}$ and $\bar{\alpha}_{t}=\prod_{s=1}^{t}\alpha_{s}$.

2.2.2 Reverse Process

At test time, generation begins with $\boldsymbol{x}_{T}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I})$, and a learned reverse process approximates the true posterior $q(\boldsymbol{x}_{t-1}|\boldsymbol{x}_{t})$. Each transition is modeled as a Gaussian:

[TABLE]

where $\boldsymbol{\mu}_{\theta}$ is predicted by a neural network, often a U-Net (Ronneberger et al.,, 2015; Ho et al.,, 2020). The variance $\sigma_{t}^{2}$ is fixed or learned, and commonly chosen as $\beta_{t}$ or $\tilde{\beta}_{t}=\frac{1-\bar{\alpha}_{t-1}}{1-\bar{\alpha}_{t}}\beta_{t}$ (Nichol and Dhariwal,, 2021).

Using Bayes’ rule, we can write in closed form:

[TABLE]

with

[TABLE]

Since $\boldsymbol{x}_{0}$ is not available at test time, the network learns to predict either $\boldsymbol{x}_{0}$ or, more commonly, the noise $\boldsymbol{\epsilon}$ used in the forward process. Using Eq. (3), we write:

[TABLE]

so the model can instead predict $\boldsymbol{\epsilon}_{\theta}(\boldsymbol{x}_{t},t)$ and compute:

[TABLE]

2.2.3 Training and Sampling

The model is trained by minimizing a simplified version of the variational lower bound:

[TABLE]

where $\boldsymbol{x}_{t}$ is sampled via Eq. (3). This trains the network to recover the added noise at any given timestep. At inference time, the generation algorithm applies Eq. (4) recursively, starting from $\boldsymbol{x}_{T}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ and sampling $\boldsymbol{x}_{t-1}$ from:

[TABLE]

Foundational papers and tutorials include (Ho et al.,, 2020; Sohl-Dickstein et al.,, 2015; Song and Ermon,, 2019; Nichol and Dhariwal,, 2021; Song et al.,, 2020; Dhariwal and Nichol,, 2021; Ronneberger et al.,, 2015; Kingma and Welling,, 2019; Turner et al.,, 2024; Weng,, 2021).

2.3 Computational Considerations

2.3.1 GANs

Despite their theoretical elegance, GANs are notoriously difficult to train in practice (Goodfellow,, 2016; Salimans et al.,, 2016). One central issue is training instability: since the generator and discriminator are simultaneously optimized in a minimax game, updates that improve one network can destabilize the other. This often leads to oscillatory dynamics or divergence, rather than convergence to a Nash equilibrium—especially in high-dimensional, non-convex settings (Arjovsky and Bottou,, 2017). Another common pathology is mode collapse, where the generator produces limited and repetitive outputs, capturing only a subset of the data distribution. This typically arises when the generator overfits to fooling the current state of the discriminator, rather than learning a general mapping.

A further challenge lies in maintaining a balance between the generator and the discriminator. If the discriminator becomes too accurate, the generator receives vanishing gradients; if the generator dominates, it may exploit weaknesses in the discriminator instead of improving output fidelity. Addressing this requires careful tuning of architectures, learning rates, and optimization strategies. Practical remedies include architectural heuristics, such as those introduced in DCGAN (Radford et al.,, 2015), and regularization techniques like one-sided label smoothing. Evaluation is also nontrivial: GANs do not optimize a tractable likelihood, so metrics like the Inception Score (IS) and Fréchet Inception Distance (FID) are commonly used, though each has known limitations (Salimans et al.,, 2016).

2.3.2 DDPMs

DDPMs are widely recognized for their stable training dynamics and their ability to generate diverse, high-quality samples. However, these benefits come at a significant computational cost. Training involves simulating millions of forward and reverse diffusion steps, and sampling is inherently sequential: each generated example requires $T$ full passes through the denoising network. As a result, DDPMs are considerably slower at inference time compared to single-step generative models such as GANs (Nichol and Dhariwal,, 2021; Song et al.,, 2020; Dhariwal and Nichol,, 2021).

These trade-offs, namely training efficiency versus sample quality and diversity, highlight the need for alternative generative modeling frameworks that offer a better balance between computational demands, interpretability, and generative performance.

2.3.3 GANs vs. DDPMs

The differences between GANs and DDPMs are substantial, both in architecture and in training philosophy. GANs learn a direct mapping from latent noise to data through an adversarial game between a generator and a discriminator. This setup allows for fast sample generation, but it often suffers from training instability, mode collapse, and sensitivity to hyperparameter tuning (Goodfellow et al.,, 2014; Salimans et al.,, 2016).

In contrast, DDPMs learn to generate data via an iterative denoising process, guided by a stable, likelihood-inspired objective. While this tends to produce samples with greater diversity and realism, the generation process is much slower due to the required sequence of denoising steps (Nichol and Dhariwal,, 2021; Ho et al.,, 2020).

Despite their conceptual differences, both frameworks rely heavily on deep neural networks as function approximators. In practice, the idealized objectives—such as the optimal discriminator in GANs or perfect noise prediction in DDPMs—are difficult to achieve. This mismatch between theoretical guarantees and empirical behavior contributes to the practical challenges each model faces (Salimans et al.,, 2016; Ho et al.,, 2020; Nichol and Dhariwal,, 2021).

3 Support Points: Theory and Algorithm

Support points offer a nonparametric approach to approximate a continuous probability distribution $F$ using a finite set of representative points. The resulting point set balances fidelity to the distribution with diversity among elements, making support points particularly attractive for tasks that require sample efficiency, interpretability, and structured coverage of the input space.

3.1 Mathematical Foundations of Support Points

The central idea behind support points is to construct a discrete probability measure that approximates a target distribution $F$ as closely as possible under a rigorous statistical criterion. This criterion relies on the energy distance, a metric between probability distributions, defined for two distributions $F$ and $G$ on a subset $\mathcal{X}\subseteq\mathbb{R}^{d}$, both with finite means, as (Székely and Rizzo,, 2005, 2013; Mak and Joseph,, 2018)

[TABLE]

where $\mathbf{X},\mathbf{X}^{\prime}\overset{\text{i.i.d.}}{\sim}G$, $\mathbf{Y},\mathbf{Y}^{\prime}\overset{\text{i.i.d.}}{\sim}F$, and $\|\cdot\|_{2}$ denotes the Euclidean norm. This metric can be interpreted as a generalized notion of potential energy between distributions and was originally proposed for multivariate, nonparametric goodness-of-fit testing (Székely and Rizzo,, 2004).

Importantly, $E(F,G)\geq 0$ for all $F,G$, and $E(F,G)=0$ if and only if $F=G$. These properties make the energy distance a valid statistical metric and a natural criterion for distributional approximation.

3.2 Definition of Support Points

Given a target distribution $F$ and a desired number of support points $n$, the goal is to construct a set of vectors $\boldsymbol{\xi}_{1},\ldots,\boldsymbol{\xi}_{n}\subset\mathcal{X}\subseteq\mathbb{R}^{d}$ that approximates $F$ as closely as possible. Let $\mathbf{S}\in\mathbb{R}^{d\times n}$ denote the matrix whose columns are the support points. These points are defined as the solution to the minimization problem

[TABLE]

where $\widehat{F}_{n}$ denotes the empirical distribution associated with the columns $\mathbf{x}_{1},\dots,\mathbf{x}_{n}$ of a candidate matrix $\mathbf{A}\in\mathbb{R}^{d\times n}$

[TABLE]

Using the definition of the energy distance between two distributions $F$ and $G$, and applying it to the case $G=\widehat{F}_{n}$, we obtain

[TABLE]

This expression follows from the fact that $\widehat{F}_{n}$ is a discrete probability measure that assigns equal mass $1/n$ to each point $\mathbf{x}_{i}$. Expectations under $\widehat{F}_{n}$ therefore take the form

[TABLE]

Let $\mathcal{P}_{n}(\mathbb{R}^{d})$ denote the space of probability measures supported on $n$ points. From this perspective, support points correspond to an optimal discrete measure $\widehat{F}_{n}\in\mathcal{P}_{n}(\mathbb{R}^{d})$ that minimizes the energy distance to a target measure $F\in\mathcal{P}(\mathbb{R}^{d})$.

This variational formulation highlights the interpretation of support points as the best approximation to $F$ within a finite-dimensional class of discrete measures, using the geometry induced by the energy distance.

3.2.1 Optimization: The Convex-Concave Procedure (CCP)

Direct minimization of the objective function (8) presents two main challenges: (i) the expectation term $\mathbb{E}_{\mathbf{Y}\sim F}\|\mathbf{x}_{i}-\mathbf{Y}\|_{2}$ is generally intractable in closed form, and (ii) the overall function is non-convex due to the interplay between attractive and repulsive components.

To address these difficulties, Mak and Joseph, (2018) propose a two-step strategy. First, the distribution $F(\mathbf{y})$ is approximated by an empirical reference measure $\widehat{F}_{N}(\mathbf{y})$, supported on a sample $\{\mathbf{y}_{m}\}_{m=1}^{N}\subset\mathbb{R}^{d}$. Letting $\mathbf{P}=[\mathbf{y}_{1},\dots,\mathbf{y}_{N}]$, the objective becomes the (unnormalized) empirical energy distance between $\widehat{F}_{n}$ and $\widehat{F}_{N}$

[TABLE]

where $\mathbf{A}=[\mathbf{x}_{1},\dots,\mathbf{x}_{n}]\in\mathbb{R}^{d\times n}$ denotes the matrix of candidate points.

In the second step, the optimization is carried out via the Convex-Concave Procedure (CCP), an iterative method for minimizing functions expressed as the difference of two convex terms. The structure of (9) naturally fits this framework:

[TABLE]

where $f_{\text{att}}$ corresponds to the attraction term, encouraging each point $\mathbf{x}_{i}$ to move toward high-density regions of $\widehat{F}_{N}$, and $f_{\text{rep}}$ encodes a pairwise repulsion that spreads the points across the domain. Both components are convex, but their difference yields a non-convex objective.

At each CCP iteration, the concave part $f_{\text{rep}}$ is linearized around the current iterate, resulting in a convex surrogate problem. Solving this surrogate yields an updated configuration of points. Conceptually, the iteration corresponds to approximating the gradient of the full objective by freezing the non-convexity, and updating each $\mathbf{x}_{i}$ through an explicit expression involving the data and the remaining support points. This process is repeated from a suitable initialization to produce a sequence $\mathbf{A}^{(0)},\mathbf{A}^{(1)},\dots$ that converges to a local minimum.

This formulation effectively minimizes $E(\widehat{F}_{N},\widehat{F}_{n})$ up to an additive constant, yielding a principled and computationally tractable surrogate for the original objective $E(F,\widehat{F}_{n})$.

We now present Algorithm 1 which approximates the support points matrix $\mathbf{S}\in\mathbb{R}^{d\times n}$ that minimizes the objective in equation (9).

The mapping $M_{i}$, which updates the $i$-th column, is defined using a reference matrix $\mathbf{A}^{\prime}=[\mathbf{x}_{1}^{\prime},\dots,\mathbf{x}_{n}^{\prime}]\in\mathbb{R}^{d\times n}$, where each $\mathbf{x}_{j}^{\prime}$ is treated as fixed during the current iteration. The updated value of $\mathbf{x}_{i}$ is given by

[TABLE]

where the normalizing factor is:

[TABLE]

In the original formulation by Mak and Joseph, (2018), the initial matrix $\mathbf{A}^{(0)}\in\mathbb{R}^{d\times n}$ is constructed by randomly selecting $n$ columns from the reference matrix $\mathbf{P}$. This strategy is simple and ensures that the initial candidate points lie within the domain of the data. However, it may lead to numerical instabilities when duplicate or nearly identical vectors are chosen, particularly in high-dimensional settings.

To address this, our implementation allows for alternative initialization schemes. By default, $\mathbf{A}^{(0)}$ can be generated by drawing independent samples from a uniform distribution over the bounding box of $\mathbf{P}$. This method avoids exact duplications and introduces greater spatial diversity in the initial configuration. The initialization strategy is configurable, enabling users to experiment with different options and select the one that yields the best empirical performance for their task.

Unlike generative models such as GANs or diffusion-based methods, which produce samples by transforming latent noise through a trained neural network, support points offer a nonparametric alternative. They approximate the target distribution directly by selecting a finite set of representative points that minimize a statistical distance. This makes them particularly appealing in applications where interpretability, sample efficiency, and structured coverage are critical—especially in data-limited or resource-constrained environments.

4 Weighted Support Points

The standard minimization problem, expression (7), and its associated algorithm converge deterministically to a single set of support points. While this property is desirable in approximation tasks, where a stable and representative summary of the distribution is required, it limits the method’s applicability in generative settings where diversity is essential.

To address this limitation, we introduce a simple yet effective modification to the optimization objective, motivated by the idea of replacing the target distribution $F$ with a random probability measure $\widetilde{F}$. This leads to a randomized version of the energy distance objective, which in turn produces multiple, diverse configurations of support points. We refer to these as weighted support points, as each realization of $\widetilde{F}$ corresponds to a weighted empirical measure supported on the original dataset.

4.1 Random Measures

In the classical formulation, support points approximate a fixed probability measure $F\in\mathcal{P}(\mathbb{R}^{d})$ by minimizing the energy distance to a discrete uniform measure supported on $n$ points. While this is appropriate for summarization tasks, it limits the generation of diverse outputs, an essential feature in generative modeling.

To introduce stochasticity, we replace the fixed distribution $F$ with a random probability measure $\widetilde{F}$, supported on a finite reference set $\{\mathbf{y}_{m}\}_{m=1}^{N}\subset\mathbb{R}^{d}$. We consider random measures of the form

[TABLE]

where $\mathbf{w}=(w_{1},\dots,w_{N})$ is a vector of nonnegative weights summing to one. This general formulation encompasses several classical constructions

• •

The empirical distribution $\widehat{F}_{N}$, where $w_{m}=1/N$ for all $m$,

• •

The Bayesian bootstrap of Rubin, (1981), where $\mathbf{w}\sim\mathrm{Dirichlet}(1,\dots,1)$,

• •

Truncated Dirichlet processes (Ferguson,, 1973), where the weights $\mathbf{w}$ are generated via the stick-breaking construction of Sethuraman, (1994).

Substituting $F$ with a realization of $\widetilde{F}_{N}$ in the energy distance objective yields a randomized variant of the support points algorithm. Each draw of $\mathbf{w}$ defines a different optimization landscape, producing distinct configurations of representative points. We refer to these as weighted support points. This approach offers a principled way to introduce controlled randomness into the method, enabling the generation of diverse outputs while maintaining fidelity to the data-generating distribution.

4.2 Definition of Weighted Support Points

We define the weighted support points as the solution to the following stochastic optimization problem:

[TABLE]

where $\widetilde{F}_{N}$ is a discrete random probability measure supported on a reference dataset $\mathbf{P}=\{\mathbf{y}_{m}\}_{m=1}^{N}\subset\mathbb{R}^{d}$, with weights $\mathbf{w}=(w_{1},\dots,w_{N})$ randomly drawn and held fixed during each optimization run.

Ignoring the constant term that does not depend on $\mathbf{A}$, the problem reduces to minimizing the following weighted empirical energy distance

[TABLE]

Solving this optimization problem for multiple independent realizations of $\mathbf{w}$ yields a diverse collection of representative point sets, each adapted to a different randomized view of the underlying data distribution.

4.2.1 Optimization

To derive the update rule for weighted support points, we compute the gradient of the objective function in equation (12) with respect to each point $\mathbf{x}_{i}$, assuming the weights $\mathbf{w}=(w_{1},\dots,w_{N})$ are fixed during the optimization.

The gradient of the attraction term with respect to $\mathbf{x}_{i}$ is:

[TABLE]

The gradient of the repulsion term with respect to $\mathbf{x}_{i}$ is:

[TABLE]

Setting the total gradient to zero yields the first-order condition

[TABLE]

Rearranging this expression leads to an explicit update rule. Let $\mathbf{x}_{i}^{\prime}$ denote the current value of $\mathbf{x}_{i}$. Then the next iterate is given by

[TABLE]

with normalization constant

[TABLE]

This fixed-point iteration is repeated until convergence to obtain a locally optimal configuration of weighted support points for the current realization of the weights $\mathbf{w}$.

The optimization follows the same iterative structure as the original sp.ccp algorithm (Algorithm 1), with two key modifications: (i) the objective function is replaced by its weighted form (12), and (ii) the update rule $M_{i}$ is modified accordingly, as shown in equation (13). To efficiently evaluate the objective and update steps, we exploit symmetries and reuse intermediate computations using inner products and matrix caching. These optimizations significantly reduce cost in large-scale settings. Full details are provided in Appendix A.

4.3 Promoting Diversity via Random Subsetting and Symmetric Dirichlet Weights

In generative settings, we seek support points that are not only representative of the full dataset but also diverse across independent runs. When relying on a random measure $\widetilde{F}_{N}$ defined over a fixed reference set $\{\mathbf{y}_{m}\}_{m=1}^{N}$, there is an inherent trade-off between

• •

Structural representativeness: keeping the full dataset accessible so that all reference points can be informative;

• •

Sample diversity: inducing run-to-run variability by emphasizing different subsets of the data in each realization.

Although (11) is random, under common constructions (e.g., Bayesian bootstrap or finite Dirichlet-based schemes) it is centered at the empirical distribution $\widehat{F}_{N}$.

Two benchmarks illustrate opposite limitations. The Bayesian bootstrap assigns weights to all $N$ atoms with $\mathbf{w}\sim\mathrm{Dir}(1,\ldots,1)$, so $\mathbb{E}[w_{i}]=1/N$ and $\mathrm{Var}(w_{i})=\frac{N-1}{N^{2}(N+1)}=\mathcal{O}(N^{-2})$. For large $N$, the induced perturbation of empirical averages is small in absolute terms, and consecutive runs remain close to $\widehat{F}_{N}$, limiting diversity. By contrast, a Dirichlet–process draw supported on $\{\mathbf{y}_{1},\ldots,\mathbf{y}_{N}\}$ yields $\mathbf{w}\sim\mathrm{Dir}(\alpha/N,\ldots,\alpha/N)$, with relative dispersion $\mathrm{CV}\approx\sqrt{N/(\alpha+1)}$; for fixed $\alpha$ this produces spiky allocations where a few atoms dominate, harming coverage and stability.

To strike a balance, we adopt a finite random measure built by random subsetting and symmetric–Dirichlet weighting calibrated by a target coefficient of variation (CV). Let $I=\{i_{1},\ldots,i_{N_{0}}\}\subset\{1,\ldots,N\}$ be a random subset sampled uniformly without replacement (the rule for choosing $N_{0}$ is given below). Given $I$, draw $\mathbf{w}=(w_{1},\ldots,w_{N_{0}})\sim\mathrm{Dirichlet}(\alpha,\ldots,\alpha)$ with total concentration $\kappa=N_{0}\alpha$ calibrated via the target CV (see expression (15)). The resulting random probability measure is

[TABLE]

Because $\mathbb{E}[w_{j}\mid I]=1/N_{0}$ for a symmetric Dirichlet and each index has equal inclusion probability under uniform subsetting, it follows that

[TABLE]

as shown in Appendix C. Hence $\widetilde{F}_{N_{0}}$ is centered at $\widehat{F}_{N}$, while its dispersion is controlled by $\kappa$ (equivalently, by the chosen CV), which enhances diversity across runs without concentrating mass on a few atoms. The subsections below detail the subsetting rule and the CV–based calibration of $\kappa$.

4.3.1 Subset size and selection

Our aim is to balance coverage of the reference set with between-run variability in the active atoms. We randomize the retained size while enforcing a minimum retention rate. Specifically, draw $\theta\sim\mathrm{Unif}(0.7,0.9)$ and set

[TABLE]

The binomial mechanism corresponds to independent thinning with inclusion probability $\theta$. Conditional on $\theta$, $\mathbb{E}[S\mid\theta]=\theta N$ and $\mathrm{Var}(S\mid\theta)=N\theta(1-\theta)$. For large $N$ and $\theta\in[0.7,0.9]$, the lower bound $\lceil 0.6N\rceil$ is rarely active; it simply prevents overly aggressive pruning so the retained dictionary remains well anchored to the data geometry. Given $N_{0}$, we sample $I$ uniformly without replacement and work with $\{\mathbf{y}_{i_{1}},\ldots,\mathbf{y}_{i_{N_{0}}}\}$. Randomizing $\theta$ across runs induces additional variability in $N_{0}$, which helps decorrelate successive realizations and mitigates near-identical configurations upon reinitialization.

4.3.2 Symmetric Dirichlet weights and CV calibration

On the selected subset, draw $\mathbf{w}\sim\mathrm{Dirichlet}(\alpha,\ldots,\alpha)$ with $\kappa=N_{0}\alpha>0$. Equivalently, sample $h_{j}\sim\mathrm{Gamma}(\alpha,1)$ and set

[TABLE]

For any $\kappa>0$,

[TABLE]

Hence the dispersion relative to the mean $1/N_{0}$ is

[TABLE]

Given a target $\mathrm{CV}>0$, we invert to obtain

[TABLE]

which is feasible whenever $\mathrm{CV}<\sqrt{N_{0}-1}$. Larger $\kappa$ (equivalently, larger $\alpha$) yields more even weights; smaller $\kappa$ yields more irregular allocations while preserving exchangeability within the subset.

The complete strategy is described in Appendix C.

5 Examples

5.1 Example: MNIST Digits 6 and 9

To illustrate the proposed method, we applied the weighted support point algorithm to a filtered subset of the MNIST dataset containing only images of the digits “6” and “9”. Each grayscale image was flattened into a column vector of dimension $784=28\times 28$, resulting in a reference matrix $\mathbf{P}\in\mathbb{R}^{784\times 11867}$.

We generated $n=10$ weighted support points in each of five independent runs (Figure 2), and also generated $n=50$ weighted support points in a single run (Figure 2). In each case, the procedure was as follows:

• •

The initial candidate matrix $\mathbf{A}^{(0)}\in\mathbb{R}^{784\times n}$ was sampled independently from a uniform distribution over $[0,1]$.

• •

A sparse random measure was generated using Algorithm

 gen.rmeasure (Appendix B) with $\mbox{CV}=0.4$, yielding a truncated set of atoms $\widetilde{\mathbf{P}}\subset\mathbf{P}$ and normalized weights $\widetilde{\mathbf{w}}$.

• •

The weighted support point algorithm esp.ccp-w (Appendix A) was run with a maximum of 1,000 iterations and a convergence tolerance of $10^{-5}$. Column-wise updates were parallelized across 16 CPU cores. Each run with $n=10$ completed in approximately 1 minute, while the run with $n=50$ took around 4 minutes.

• •

Final pixel values were clipped to $[0,255]$ and rounded to the nearest integer.

Each row of Figure 2 illustrates the diversity across independent runs, while Figure 2 shows no discernible pattern, indicating that different representative images are consistently selected. This demonstrates that the stochasticity introduced by the random measure effectively promotes variability while preserving structure and legibility.

5.1.1 Empirical comparisons

To enable an empirical comparison, we generated samples of digits “6” and “9” using two neural-network-based generative approaches: a GANs (Figure 4) and a DDPMs (Figure 4). Both models were trained on the same dataset used for the weighted support points strategy.

For the GANs, we used a lightweight fully connected architecture with spectral normalization in the discriminator. The generator maps a 100-dimensional noise vector to a $28\times 28$ grayscale image through three linear layers with ReLU activations and a final Tanh output. The discriminator processes flattened $28\times 28$ images through two linear layers with LeakyReLU activations and a final linear output for real/fake logits. Training used the Adam optimizer ($\mathrm{lr}=2\times 10^{-4}$, $\beta_{1}=0.5$, $\beta_{2}=0.999$), binary cross-entropy loss, batch size 128, label smoothing for real labels (0.9), $k=2$ discriminator updates per generator update, and up to 100 epochs. Once trained, the model generates images in milliseconds, producing recognizable digits 6 and 9 with relatively clean shapes.

For the DDPMs, we used a lightweight U-Net adapted to $32\times 32$ grayscale images as the noise predictor at each diffusion step. Training employed the Adam optimizer ($\mathrm{lr}=2\times 10^{-5}$), batch size 64, and 5,000 epochs, minimizing the mean squared error (MSE) between true and predicted noise under a linear noise schedule with $T=1,000$ steps. Compared to the GANs, training and image generation are considerably slower, but the resulting samples show smoother digit shapes and fewer visual artifacts, demonstrating the advantages of the diffusion-based approach even with a simple architecture.

Next, we compare these results with those obtained using our proposed weighted support points approach. As shown in Figures 2 and 2, the method selects clean and aesthetically pleasing representative digits “6” and “9”. Although slightly blurred due to the use of distance-based optimization rather than pixel-level generative modeling, the outputs remain legible and structurally coherent. In contrast, the GAN (Figure 4) produces recognizable digits almost immediately after training, but with some speckled artifacts. The DDPMs (Figure 4) yields smoother results with fewer artifacts, though it occasionally generates distorted shapes, still clearly identifiable as “6” or “9”. Overall, weighted support points offer a principled, interpretable, and computationally efficient alternative to neural generative models, producing diverse and stable samples without the need for heavy training.

A key methodological difference lies in how the data are represented and processed. Neural-network-based models (GANs and DDPMs) operate on images as multidimensional arrays, or tensors, which encode batches of images along with their channels and spatial dimensions. These models also require dataset-specific transformations tightly coupled to their architectures and training procedures. In contrast, the weighted support points approach works directly on vectorized image representations, avoiding architecture-dependent preprocessing and enabling distance-based criteria to be applied in their native form.

It is worth noting that both the GANs and DDPMs can, in principle, achieve arbitrarily high performance by increasing the complexity of their architectures and training procedures. Such enhancements typically yield sharper and more realistic digits but at a substantial computational cost. In contrast, our weighted support points approach deliberately avoids neural architectures, focusing instead on interpretable and efficient sample generation.

5.2 Example: CelebA-HQ Faces

As a second example, we apply our algorithm to a subset of the CelebA-HQ dataset, which contains 27,000 high-quality color images of celebrity faces captured under unconstrained conditions. Each image has a resolution of $256\times 256$ pixels and three color channels (RGB), yielding $256\times 256\times 3=196{,}608$ pixel values per image. Flattening all images into column vectors and stacking them produces a matrix

[TABLE]

with over 5.3 billion entries. Storing this matrix in double precision (8 bytes per entry) would require

[TABLE]

which is infeasible on standard hardware.

To reduce memory requirements, we downscale each image to $144\times 144$ pixels while preserving the RGB channels. Each image then becomes a vector in $\mathbb{R}^{62{,}208}$, and a random subset of 10,000 images is selected to form the reduced matrix $\mathbf{P}\in\mathbb{R}^{62{,}208\times 10{,}000}$. A sample of ten resized images is shown in Figure 5.

We now apply the same methodology as for the MNIST Digits “6” and “9” to the reduced CelebA-HQ dataset. In a single run, we computed 25 weighted support points using a random measure $\widetilde{F}_{N_{0}}$ constructed setting $\mbox{CV}=0.4$ and via algorithm gen.rmeasure in Appendix B. The optimization algorithm esp.ccp-w was initialized with a uniformly sampled matrix. The execution required approximately 3 hours to converge under tolerance $\epsilon=10^{-4}$ and a maximum of 1,000 iterations. Figure 6 displays the resulting weighted support points.

Each row corresponds to a different run, and while some of the generated faces appear visually similar across runs, closer inspection reveals subtle but meaningful differences, such as variations in facial orientation or lighting. This suggests that the stochastic mechanism introduced by the random measure promotes diversity without compromising structural consistency. The observed variability confirms the method’s ability to explore different representative configurations while remaining faithful to the overall data distribution.

6 Conclusion

We introduced a generative modeling framework based on weighted support points driven by a finite random measure built from random subsetting and symmetric–Dirichlet weights calibrated by a target coefficient of variation (CV) (see Section 4.3). This construction yields two complementary sources of stochasticity, subset selection and weight randomness, while remaining centered at the empirical distribution (expression (14) and Appendix C). Calibrating the total concentration via the CV identity (15) provides explicit control over dispersion, preventing mass from collapsing onto a few atoms and promoting genuine run–to–run diversity.

Empirically, the method produces high–quality and diverse representatives on MNIST and CelebA–HQ at a fraction of the computational cost of neural approaches. The optimization is geometry–aware, the random measure is simple and transparent, and the overall pipeline is memory–efficient because it operates on a subset of atoms of size $N_{0}$. In contrast to standard centering schemes, such as the Bayesian bootstrap, which perturbs all $N$ atoms only mildly in absolute terms for large $N$, or Dirichlet–process draws on a fixed support, which often induce overly spiky allocations for fixed $\alpha$ our CV-calibrated subsetting strategy strikes a principled balance between structure and diversity.

Beyond the present results, the framework invites several extensions: stratified or geometry–aware subsetting (for example, $k$–center or determinantal selection), data–driven calibration of the CV, and finite–sample guarantees linking CV and diversity of the resulting support–point configurations. Overall, weighted support points with CV–calibrated symmetric–Dirichlet random measures constitute a practical, interpretable, and scalable alternative for generative modeling when stability, transparency, or limited computational resources are paramount.

Supplementary Material and Data Availability

All code and datasets used to implement the algorithms and generate the experimental results presented in this paper is available as supplementary material. This includes the R implementation of the weighted support points algorithm, the sampling routine for truncated Dirichlet processes, and scripts for replicating the MNIST and CelebA-HQ experiments. Additionally, Python scripts for training and generating samples with both a simple GANs and a DDPMs on MNIST digits “6” and “9” are provided to facilitate reproducibility of the neural network–based comparison.

Funding

The authors did not receive any specific funding for this work.

Conflict of Interest

The authors declare that there is no conflict of interest.

Appendix A Appendix A: Efficient Evaluation for Weighted Support Points

The two main computational challenges in our method are evaluating the weighted cost function $MC(\mathbf{A}^{(l)};\mathbf{P},\mathbf{w})$ and updating the support points. For high-dimensional datasets, such as when $\mathbf{P}$ contains images, naive implementations become computationally expensive. This section presents an efficient strategy that leverages algebraic identities, matrix operations, and symmetry to reuse intermediate computations and accelerate both cost evaluation and optimization updates.

In the repulsion term, we exploit symmetry and the identity $\|\mathbf{x}_{i}^{(l)}-\mathbf{x}_{j}^{(l)}\|_{2}=\|\mathbf{x}_{j}^{(l)}-\mathbf{x}_{i}^{(l)}\|_{2}$. This allows us to compute each distance only once and scale appropriately. The weighted cost function, where the weights are placed on the reference atoms $\mathbf{y}_{m}$, takes the form

[TABLE]

To compute the pairwise distances efficiently, we rely on algebraic expansions of squared Euclidean norms. These allow us to express distances in terms of squared norms and inner products:

[TABLE]

These expressions can be evaluated efficiently using matrix operations. To structure the computation, we distinguish between static and dynamic quantities: static quantities are computed once and reused across all iterations, while dynamic quantities depend on the current value of $\mathbf{A}^{(l)}$ and are updated in each iteration.

The static component consists of the squared norms of the data points:

[TABLE]

At iteration $l$, the dynamic cache $\mathcal{C}^{(l)}=\texttt{Cache}(\mathbf{A}^{(l)},\mathbf{P},\mathbf{r}_{p})$ includes:

[TABLE]

where $r_{x,i}^{(l)}=\|\mathbf{x}_{i}^{(l)}\|_{2}^{2}$, and the inner product matrices are

[TABLE]

To improve numerical stability in the evaluation of $D_{xp}^{(l)}$ and $D_{xx}^{(l)}$, we replace the standard Euclidean norm with a smooth regularized form

[TABLE]

where $\varepsilon$ is scaled to the problem as \varepsilon=c\sqrt{\texttt{.Machine\double.eps}}\cdot s$. Here,$s$denotes a characteristic norm of the current configuration and$c$ is a small constant. This scale-aware regularization mitigates instabilities when points coincide or are extremely close, ensures gradients remain well-behaved, and avoids the inefficiencies that can arise from overly conservative numerical adjustments, while preserving the efficiency of the underlying vectorized matrix-based computation.

These cached quantities are used both to evaluate the cost function and to compute the update steps efficiently. The final form of the cost function becomes

[TABLE]

The optimization step is also accelerated using the cached matrices. Each candidate point is updated as:

[TABLE]

where the normalization constant is

[TABLE]

The key matrix operations, such as computing $G_{xp}^{(l)}=\mathbf{A}^{(l)\top}\mathbf{P}$ and $G_{xx}^{(l)}=\mathbf{A}^{(l)\top}\mathbf{A}^{(l)}$, are highly optimized in standard linear algebra libraries (e.g., BLAS level-3).

Appendix B Appendix B: Sampling a Random Measure via Random Subsetting and Symmetric Dirichlet Weights

This appendix (see Algorithm 3 below) describes the random subsetting and symmetric Dirichlet approach described in Section 4.3. The procedure draws a subset of atoms to introduce structured sparsity and then assigns Dirichlet–symmetric weights calibrated by a target coefficient of variation (CV). The resulting pair $(\widetilde{\mathbf{P}},\widetilde{\mathbf{w}})$ defines a discrete random measure that is centered at the empirical distribution while allowing explicit control of weight dispersion and improved run–to–run diversity.

Appendix C Appendix C: Centered

Random Measure

Define the random finite measure

[TABLE]

where $I=\{i_{1},\ldots,i_{N_{0}}\}$ is sampled uniformly without replacement from $\{1,\ldots,N\}$, and $\mathbf{w}=(w_{1},\ldots,w_{N_{0}})\sim\mathrm{Dirichlet}(\alpha,\ldots,\alpha)$ (possibly with $\alpha$ depending on $N_{0}$ through the CV calibration). Then

[TABLE]

Since the Dirichlet is symmetric, $\mathbb{E}[w_{j}\,|\,I,N_{0}]=1/N_{0}$ for all $j$, hence

[TABLE]

Conditionally on $N_{0}$, uniform sampling without replacement gives $\mathbb{P}(m\in I\,|\,N_{0})=N_{0}/N$, so

[TABLE]

Remarks. (i) The result holds for any distribution of $N_{0}$ provided that, conditional on $N_{0}$, the subset $I$ is uniformly distributed. (ii) The Dirichlet concentration (or the CV-derived $\kappa$) may depend on $N_{0}$; symmetry ensures $\mathbb{E}[w_{j}\,|\,I,N_{0}]=1/N_{0}$, so centering is unaffected.