# Gaussian Mixture Model with unknown diagonal covariances via continuous sparse regularization

Romane Giard, Yohann de Castro, Clément Marteau

<figure><img src="image_1.png" /></figure>

## To cite this version:

Romane Giard, Yohann de Castro, Clément Marteau. Gaussian Mixture Model with unknown diagonal covariances via continuous sparse regularization. 2026. hal-05237411v3

HAL Id: hal-05237411

https://hal.science/hal-05237411v3

Preprint submitted on 8 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Gaussian Mixture Model with unknown diagonal covariances via continuous sparse regularization

Romane Giard$^{(1)}$, Yohann De Castro$^{(1,3)}$ and Clément Marteau$^{(2)}$

$^{1}$Centrale Lyon, CNRS UMR 5208, Institut Camille Jordan, Écully, France.

$^{2}$Univ. Claude Bernard, CNRS UMR 5208, Institut Camille Jordan, Villeurbanne, France.

$^{3}$Institut Universitaire de France (IUF)

version as of January 8, 2026

## Abstract

This paper addresses the statistical estimation of Gaussian Mixture Models (GMMs) with unknown diagonal covariances from independent and identically distributed samples. We employ the Beurling-LASSO (BLASSO), a convex optimization framework that promotes sparsity in the space of measures, to simultaneously estimate the number of components and their parameters.

Our main contribution extends the BLASSO methodology to multivariate GMMs with component-specific unknown diagonal covariance matrices. This setting is significantly more flexible than previous approaches, which required known and identical covariances. We establish non-asymptotic recovery guarantees with nearly parametric convergence rates for component means, diagonal covariances, and weights, as well as for density prediction.

A key theoretical contribution is the identification of an explicit separation condition on mixture components that enables the construction of non-degenerate dual certificates—essential tools for establishing statistical guarantees for the BLASSO. Our analysis leverages the Fisher-Rao geometry of the statistical model and introduces a novel semi-distance adapted to our framework, providing new insights into the interplay between component separation, parameter space geometry, and achievable statistical recovery.

Corresponding author: Romane Giard

E-mail address: romane.giard@ec-lyon.fr

## 1 Introduction

Gaussian Mixture Models (GMMs) are a cornerstone of statistical modeling, offering a flexible and powerful framework for representing complex data distributions. They are widely applied in diverse fields, including clustering and density estimation [McLachlan and Peel, 2000; Bouveyron et al., 2019], image processing [Houdard et al., 2018], bioinformatics and economics [Fruhwirth-Schnatter et al., 2019, Section 8.4]. Despite their ubiquity, the estimation of GMMs, particularly when the number of components and their covariance structures are unknown, presents significant statistical and computational challenges.

The predominant method for GMM estimation is the Expectation-Maximization (EM) algorithm, which iteratively maximizes the log-likelihood. While the EM algorithm guarantees a non-decreasing likelihood and converges to a stationary point under some conditions [McLachlan and Krishnan, 1997; Wu, 1983; Balakrishnan et al., 2017], it faces significant practical challenges. The log-likelihood function is non-convex and multi-modal, making the algorithm sensitive to initialization and prone to converging to local maxima rather than the global optimum. Crucially, it requires the number of mixture components to be specified in advance. We aim to overcome these limitations.

This paper introduces a novel approach to the estimation of multivariate GMMs with unknown diagonal covariances, leveraging the Beurling-LASSO (BLASSO) methodology [De Castro and Gamboa, 2012; Candès and Fernandez-Granda, 2014; Duval and Peyré, 2015; Boyer et al., 2017; Poon, 2019; De Castro et al., 2021a]. The BLASSO offers a convex optimization framework in the space of measures. It promotes sparsity, thereby offering a principled way to simultaneously estimate the number of components, their respective weights, means, and diagonal covariance matrices. It has not yet been applied to settings where each mixture component is allowed to have its own covariance structure. This extension makes it possible to handle significantly more realistic models. For instance, in clustering applications, assuming identical covariances across components yields Voronoi partitions (linear decision boundaries), as in k-means. In contrast, allowing distinct covariance matrices leads to boundaries defined by quadratic equations [Murphy, 2012, Section 4.2].
A pivotal contribution of our research is the identification of a separation condition for the mixture components. This condition is instrumental in the construction of so-called non-degenerate dual certificates [Candès and Fernandez-Granda, 2014; Duval and Peyré, 2015], which play a key role in establishing estimation guarantees. Our analysis is grounded in the Fisher-Rao geometry of the statistical model, providing theoretical insights into the intricate relationships between component separation, the underlying geometry of the parameter space, and the achievable statistical recovery.

## 1.1 Continuous sparse regression for Gaussian Mixture Models

We denote by $N(t, C)$ a multivariate Gaussian distribution in dimension $d \in \mathbb{N}^*$, with mean $t \in \mathbb{R}^d$ and covariance matrix $C \in \mathbb{S}_{++}^d$, where $\mathbb{S}_{++}^d$ denotes the space of positive-definite symmetric matrices of size $d \times d$. Given $u = (u_1, \dots, u_d) \in (\mathbb{R}_+^*)^d$ a vector with positive entries, we denote by $\text{diag}((u_1)^2, \dots, (u_d)^2)$ the $d \times d$ diagonal matrix with diagonal entries $u_k^2$, where each $u_k$ is interpreted as a marginal standard deviation. The density of $N(t, \text{diag}((u_1)^2, \dots, (u_d)^2))$ is denoted by $\varphi_{(t,u)}$.

In this paper, we observe a $n$-sample $X_1, \dots, X_n \in \mathbb{R}^d$ drawn from a Gaussian Mixture distribution with diagonal covariances, defined by:

$$X_i \stackrel{\text{i.i.d.}}{\sim} \sum_{j=1}^{s} a_j^0 \varphi_{(t_j^0, u_j^0)} := f^0 \quad \text{with} \quad \sum_{j=1}^{s} a_j^0 = 1 \quad \text{and} \quad a_j^0 > 0 \qquad (1)$$

where $t_j^0 \in \mathbb{R}^d$, $u_j^0 = (u_{j,k})_{k=1,\dots,d} \in (\mathbb{R}_+^*)^d$. Our goal is to estimate the number of components $s$ (also called the sparsity index), the weights $a_j^0$ and the location-scale parameters $(t_j^0, u_j^0)$ of each mixture component, indexed by $j \in \{1, \dots, s\}$, from the observations $X_1, \dots, X_n$.

**Remark 1.1.** As a matter of fact, the Gaussian Mixture Model with unknown number of components, weights, means and covariance matrices:

$$\left\{ \sum_{j=1}^{p} a_j N(t_j, C_j) : p \ge 1, a_j > 0, t_j \in \mathbb{R}^d, C_j \in \mathbb{S}_{++}^d, \forall i \ne j, (t_j, C_j) \ne (t_i, C_i), \sum_{j=1}^{p} a_j = 1 \right\}$$

is identifiable, i.e., if two distributions of this model are equal, then they have the same number of components, and the same components (up to a permutation) with the same associated weights. We refer to [Teicher, 1961] for the proof in dimension $d=1$, and [Yakowitz and Spragins, 1968] for the generalization in dimension $d \ge 1$. Under mild assumptions, an even stronger result holds for continuous mixtures of Gaussian distributions, defined by a density on the space of mean and covariance matrix [Bruni and Koch, 1985].

Our estimation strategy relies on the Beurling-LASSO (BLASSO). This framework, introduced in [De Castro and Gamboa, 2012; Candès and Fernandez-Granda, 2014], has been successfully applied to various statistical estimation problems, particularly in the context of sparse signal recovery [Duval and Peyré, 2015] and compressed sensing [Poon et al., 2023]. The BLASSO approach is characterized by its ability to promote sparsity in the space of measures, enabling the recovery of discrete measures from continuous data. This framework has been extended to various settings, including the estimation of Gaussian Mixture Models with known covariances [De Castro et al., 2021a; Poon et al., 2023]. Algorithmic implementations are discussed in [Chizat, 2022] and [De Castro et al., 2023]. One key modeling idea of this approach is to lift the parameter space into the space of measures via the embedding

$$(a_j, t_j, u_j)_{j=1}^p \mapsto \mu = \sum_{j=1}^p a_j \delta(t_j, u_j),$$

where $(t_j, u_j)_{j=1}^p$ are referred to as the particles of $\mu$. Any discrete probability measure on $\mathbb{R}^d \times (0, +\infty)^d$ with finite support, of size $p$ for any $p \ge 1$, describes a set of parameters. We stress that the parameter $p$, the number of components, is free and not prescribed. The law of the $n$-sample $(X_i)_{i=1}^n$ (see (1)) is unambiguously described by the target parameters $(a_j^0, t_j^0, u_j^0)_{j=1}^s$ (by identifiability of the model, see Remark 1.1), and hence unambiguously represented by the so-called target measure $\mu^0 := \sum_{j=1}^s a_j^0 \delta(t_j^0, u_j^0)$.

Our estimator will be defined as a solution of a minimization problem over the space of measures, constructed from the observations $X_1, \dots, X_n$. It is designed to estimate the target measure $\mu^0$. In particular, we want our estimator to put a mass close to $a_j^0$ around the particle $(t_j^0, u_j^0)$. We will consider a convex loss of the form

$$F_{n,\tau}(\mu) + \kappa R(\mu) \qquad (2)$$

where $F_{n,\tau}(\mu)$ is a data fidelity term, $R(\mu)$ is a regularization term enforcing sparsity, and $\kappa > 0$ a tuning parameter. The data fidelity term compares a predicted density encoded by $\mu$ with an empirical approximation
of the target density. This approximation is obtained from the empirical distribution of $(X_i)_{i=1}^n$ by convolution with a Gaussian kernel of covariance $\tau^2\text{Id}_d$, depending on a smoothing parameter $\tau$. This approach is standard in kernel density estimation [Tsybakov, 2008], but the main difference here is that $\tau$ will not be necessarily chosen to match some nonparametric rate. Its calibration will sometimes answer to an alternative purpose. The regularization term $R(\mu)$ is the total variation (TV) norm, analog of the $\ell^1$-norm for measures, and aims to concentrate the mass of our estimator in a few regions. We refer to Section 2 and $(P_\kappa)$ for a complete description of the BLASSO procedure.

## 1.2 Related works

Gaussian mixture estimation ranges from clustering to parameter estimation. Our focus is primarily on the latter. Prior work has investigated a variety of settings under different modeling assumptions.

**Gaussian location mixtures** When covariances are known, [Heinrich and Kahn, 2018] derived local asymptotic minimax rates for univariate mixtures. Generalizing to dimension $d$, [Doss et al., 2023] established minimax rates under $W_1$ and Hellinger distances. [Regev and Vijayaraghavan, 2017] further highlighted the gap between statistical and computational limits: while a mean separation of $\sqrt{\ln s}$ suffices for identifiability with a polynomial sample size, efficient algorithms typically require a separation of order $\sqrt{d}$. Other approaches include Bayesian methods with priors on the number of components [Ohn and Lin, 2023] or approximations by finite mixtures [Ma et al., 2025].

**Unknown covariances** Intermediate settings with shared unknown variance have been studied using the Denoised Method of Moments (DMM) [Wu and Yang, 2020] in the univariate case, or EM [Bing et al., 2025] in dimension $d$ under a separation condition on the components.

The general case of varying unknown covariances remains challenging. Restricted to two components, efficient estimators exist given sufficient separation [Hardt and Price, 2015] or specific weight constraints [Wu and Zhou, 2021].

When the number of mixture components $s$ is unknown, [Ho and Nguyen, 2016] analyzed the Maximum Likelihood Estimation (MLE) in over-parameterized settings, establishing slow polynomial convergence rates dependent on the excess number of components. [Guha et al., 2021] also studied this setting in a bayesian framework. In both works, the resulting convergence guarantees are pointwise in the sense that the bounds depend on the underlying target mixture. Alternatively, [Belkin and Sinha, 2015] proposed a method of moments (MoM) with polynomial sample complexity when $s$ is known; however, their grid-search approach remains exponentially hard in the dimension. Finally, [Chen and Zhang, 2024] derived minimax bounds for clustering risks with anisotropic covariances using Lloyd's algorithm.

**Positioning** Our work departs from these approaches by adopting the Beurling-LASSO (BLASSO) framework. Unlike the EM algorithm [Bing et al., 2025] or MoM [Belkin and Sinha, 2015], this convex optimization approach does not require specifying $s$. While BLASSO is established for fixed covariances [De Castro et al., 2021a], extending it to multivariate GMMs with component-specific unknown diagonal covariances presents a distinct theoretical challenge. We address this by introducing a geometry-adapted semi-distance and establishing non-asymptotic recovery guarantees for both parameter estimation and density prediction. Our results are not pointwise (except for the sparsity index estimation displayed in Section 6); rather, they hold uniformly over classes of mixtures satisfying a prescribed separation condition.

## 1.3 Contributions

Our primary contribution is the extension of the BLASSO framework to estimate GMMs with unknown diagonal covariances. This setting introduces a key challenge: the associated kernel (see Section 5) is not translation-invariant, departing from many standard BLASSO applications. We address this by:

* Introducing a reparametrization of the measures to work with a normalized kernel—an object essential for the theoretical analysis of the BLASSO, describing the correlation between 2 location-scale parameters $(t, u)$ and $(t', u')$.

* Establishing recovery guarantees through the construction of non-degenerate dual certificates. This involves proving a modified version of the *local positive curvature* assumption (LPC) from [Poon et al., 2023, Assumption 1] for our specific non translation-invariant kernel.

* Using a semi-distance naturally aligned with the kernel-induced geometry. This allows us to formulate more tractable conditions for the certificate construction, relaxing the reliance on the Fisher-Rao distance used in related works.
* Establishing conditions on the separation between the components of the mixture (i.e., the particles of $\mu^0$) with respect to the semi-distance. These conditions depend on bounds on the variance, the sparsity $s$, the dimension $d$, and the smoothing parameter $\tau$. Higher sparsity levels or looser bounds on the variances necessitate a wider minimal separation between the components of $\mu^0$. See (17) for a precise condition.

We also provide novel prediction guarantees for the target density, achieving nearly parametric rates of convergence in different regularization regimes.

**Informal results on error bounds** We define $\mu_{\omega}^0$ as a reparametrized version of the target measure $\mu^0$,

$$ \mu_{\omega}^{0} = \sum_{j=1}^{s} \omega_{j}^{0} \delta(t_{j}^{0}, u_{j}^{0}) \quad \text{with} \quad \omega_{j}^{0} = W(x_{j}^{0}) a_{j}^{0}, $$

where $W$ is a positive function that will be specified later. According to the procedure displayed in Section 2, we estimate $\mu_{\omega}^0$ rather than $\mu^0$. Our estimator is given by the argument minimum of a function $J_W$ of the form (2) over nonnegative measures with support restricted to a compact set $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, +\infty)^d$ (see $(\mathcal{P}_{\kappa})$ for a formalized version). Here, $u_{\min}$ denotes a lower bound on the diagonal elements of the covariance matrix.

Most of our results remain valid for approximate solutions, that is for any measure $\mu$ satisfying $J_W(\mu) \le J_W(\mu_{\omega}^0)$, and not only the exact argument minimum solution of $J_W$. Accordingly, we will use the notation $\mu_{n,\omega}^*$ when an exact solution is required, and $\hat{\mu}_{n,\omega}$ when an approximate solution suffices.

We evaluate the estimator by comparing its mass against the mass of $\mu_{\omega}^0$ within “near regions” $\mathcal{X}_j^{\text{near}}(r_e)$ centered on the true parameters $(t_j^0, u_j^0)$, and in the complementary “far region” within $\mathcal{X}$. Near regions correspond to balls of radius $r_e$ for a semi-distance $d$ (given by (11) below). Figure 1 illustrates this partitioning. The near region $\mathcal{X}_j^{\text{near}}(r_e)$ depends only on its radius $r_e$, $x_j^0$, the dimension $d$, and the smoothing parameter $\tau$.

<figure><img src="image_3.png" /><figcaption>Figure 1: Schematic representation for Gaussian mixture models in dimension $d=1$. Both parameters $u$ (standard deviation) and $t$ (mean) are in dimension 1, resulting in a 2-dimensional plot in location-scale space $(t, u)$. The discs represent the near regions, shown schematically: these regions correspond to balls defined with respect to a semi-distance, not the Euclidean distance. The hatched area corresponds to the far region.</figcaption></figure>

Our estimator satisfies the following properties, given in expected value over $X_1, \dots, X_n$.

**Theorem 1.1** (Recovery guarantees for the estimation of $\mu_{\omega}^0$, informal result). Assume that the particles $(t_j^0, u_j^0)$ of $\mu^0$ are sufficiently separated, where the minimal separation constraint only depends on the dimension $d$, the sparsity index $s$, bounds on the variance and choice of a smoothing parameter $\tau \le u_{\min}$. Choosing as regularization parameter $\kappa = \frac{\sqrt{2}}{(2\pi)^{d/4}\tau^{d/2}\sqrt{n}}$, for any $r_e$ such that $0 < r_e \le r$ with $r$ a fixed constant depending on $d$, it holds that, omitting the dependence on $d$,

$$ \mathbb{E} [ |\mu_{\omega}^{0}(\mathcal{X}_{j}^{\text{near}}(r_e)) - \hat{\mu}_{n,\omega}(\mathcal{X}_{j}^{\text{near}}(r_e)) | ] \lesssim \frac{s}{r_e^2 \sqrt{n\tau^{d/2}}} \quad \forall j = 1, \dots, s \quad (3) $$

and

$$ \mathbb{E} \left[ \left| \hat{\mu}_{n,\omega} \left( \mathcal{X} \setminus \bigcup_j \mathcal{X}_j^{\text{near}}(r) \right) \right| \right] \lesssim \frac{s}{\sqrt{n\tau^{d/2}}} . $$

We refer to Theorems 3.1 and 5.1, Proposition 3.1 with Lemma 5.4 for precise statements and related discussions. Note that the constant factor omitted in the bound (3) does not depend on the target mixture, only on the
dimension $d$. Our estimator allows us to locate the particles of the reparametrized target measure on fixed regions with a parametric convergence rate. The radius $r_e$ can be chosen as desired (sufficiently small), leading to degraded bounds (the rate depends on $r_e^{-2}$). We give in Corollary 3.1 a variation of the previous result, for a direct estimator of $\mu^0$—and with decreasing size of regions.

Remark 1.2. Note that in the above result, $\kappa$ does not depend on $s$ (unknown in practice). We call it the s-agnostic choice. Making the s-dependent choice $\kappa = \frac{\sqrt{2}}{(2\pi)^{d/4}\tau^{d/2}\sqrt{sn}}$, we obtain a rate of $\frac{\sqrt{s}}{r_e^2\sqrt{n\tau^{d/2}}}$ for the bound (3) (better dependence on $s$). See Remark 3.6.

We provide at the same time prediction error bounds. We achieve an almost parametric rate (up to a logarithmic factor) for the prediction of the target density $f^0$, in two distinct regimes: under small regularization (Proposition 4.1) and with larger regularization (Theorem 4.1) when the assumption of Theorem 1.1 is verified. We present an informal version of these results.

**Theorem 1.2** (Recovery guarantees for the prediction of $f^0$, informal results). Let $\tau = \frac{\sqrt{2}u_{\min}}{\sqrt{\ln n}}$. We construct an estimator $\hat{h}_n$ of $f^0$ from $\hat{\mu}_{n,\omega}$. We consider two regimes.

* Assume that the particles $(t_j^0, u_j^0)$ of $\mu^0$ are sufficiently separated. Choosing $\kappa = \frac{\sqrt{2}(\ln n)^{d/4}}{(2\pi)^{d/4}(2u_{\min}^2)^{d/4}\sqrt{n}}$, omitting the dependence on $d$ and on bounds on the variance,

$$ \mathbb{E}\left[\left\|\hat{h}_n - f^0\right\|^2_{L^2(\mathbb{R}^d)}\right] \lesssim \frac{s(\ln n)^{d/2}}{n}. $$

* Without separation assumptions on $(t_j^0, u_j^0)_j$, choosing $\kappa = \frac{4(\ln n)^{d/2}}{(2\pi)^{d/2}(2u_{\min}^2)^{d/2}n}$, omitting the dependence on $d$ and on bounds on the variance,

$$ \mathbb{E}\left[\left\|\hat{h}_n - f^0\right\|^2_{L^2(\mathbb{R}^d)}\right] \lesssim \frac{(\ln n)^{d/2}}{n}. $$

The regularization parameter $\kappa$ is allowed to depend on $n$. Its value has some importance on the performances of the procedure. It is calibrated according to the objective to be achieved: estimation, prediction or both. Under some conditions on $\mu^0$, large regularization provides good estimation and prediction. Small regularization allows the BLASSO to focus on minimizing the data fidelity term, yielding good prediction without any separation condition on the particles of $\mu^0$. For theoretical purposes, we will also briefly discuss an s-dependent choice of $\kappa$, which leads to improved bounds in the large regularization regime (the bound is of order $\frac{(\ln n)^{d/2}}{n}$ for $s$ reasonably small, see Remark 4.4 below).

We finally derive an alternative result for $\mu_{n,\omega}^*$, an exact solution to $(\mathcal{P}_\kappa)$, under some specific conditions. We provide below an informal result, formalized in Corollary 6.1.

**Theorem 1.3** (Sparsity of the estimator for a large sample size, informal result). Let $0 < \tau \le u_{\min}$ (fixed). Assume that the particles of $\mu^0$ are sufficiently separated, where the minimal separation constraint only depends on $d, s, \tau$ and bounds on the variance. Choosing $\kappa = \frac{\alpha\sqrt{\ln n}}{(2\pi)^{d/4}\tau^{d/2}\sqrt{n}}$ for any $\alpha > 0$, there exists $n_0 \in \mathbb{N}$ depending on $\mu^0, \mathcal{X}, \tau, \alpha$, such that if the sample size $n$ verifies $n \ge n_0$, then $\mu_{n,\omega}^*$ is s-sparse with probability greater than $1 - C_\Gamma n - \frac{\gamma_0^2\alpha^2}{C_\Gamma^2}$, with $C_\Gamma$ a universal positive constant and $\gamma_0$ depending on $\mu^0, \mathcal{X}, \tau$. Moreover, writing $\mu_{n,\omega}^* = \sum_{j=1}^s \omega_j^*\delta_{x_j^*}$, omitting constants depending on $\mu^0, \mathcal{X}, \tau$ we have for all $j = 1, \dots, s$

$$ |\omega_j^0 - \omega_j^*| \lesssim \alpha \sqrt{\frac{\ln n}{n}} \quad \text{and} \quad d((t_j^*, u_j^*), (t_j^0, u_j^0))^2 \lesssim \alpha \sqrt{\frac{\ln n}{n}}, \qquad (4) $$

where $d(\cdot, \cdot)$ is the semi-distance between location-scale parameters defined in (11) below.

The result displayed here provides a different flavor on the BLASSO performances. First, this bound holds provided the sample size is large enough (and under separation conditions that are slightly stronger than those required in Theorem 1.1). In such case, we first establish that the measure $\mu_{n,\omega}^*$ has, with high probability, exactly the same sparsity index $s$ than the target $\mu^0$. Moreover, the bound (4) provides theoretical guarantees on the estimation of the mixture parameters themselves instead of a control on far and near regions (Theorem 1.1).
**Outline** Section 2 introduces the statistical framework and the BLASSO estimators $\hat{\mu}_{n,\omega}$, $\mu_{n,\omega}^*$ used for recovering Gaussian Mixture Models. In Section 3, we establish recovery guarantees for the estimation of the target measure, relying on the existence of non-degenerate dual certificates. Section 4 focuses on prediction guarantees for the density $f^0$, providing rates of convergence under different regularization regimes. In Section 5, we construct non-degenerate dual certificates by analyzing the kernel properties and deriving sufficient conditions on $\mu^0$. Finally, Section 6 demonstrates that, for sufficiently large sample sizes and under sufficient separation between components, the estimator $\mu_{n,\omega}^*$ is, with high probability, a discrete measure. Section 7 concludes with a discussion of potential extensions, algorithmic considerations, and open problems. It also comments on the separation condition required in our analysis. The main results are summarized in Tables 1 and 2 with corresponding assumptions, choices of smoothing parameter and regularization. Table 3 summarizes notation used throughout the paper. All proofs and technical results are gathered in the appendix at the end of the paper. To facilitate verification, we also provide a notebook (the associated Zenodo repository can be found at [Giard, 2025]) implementing our calculations using a symbolic computation library.

# 2 Statistical modeling

Our approach operates in the space of Radon measures, targeting the recovery of a sparse (i.e. discrete) measure. Each particle of this discrete measure encodes the parameters of a Gaussian component, and its associated mass corresponds to the component's proportion in the mixture. In this section, we introduce the necessary notation, describe the statistical model, and present the BLASSO estimator.

## 2.1 Radon measures

We first start with some definitions allowing for a rigorous introduction of the BLASSO principle.

For $A \subset \mathbb{R}^p$ locally compact, we denote $C(A)$ the set of all continuous functions from $A$ to $\mathbb{R}$ and $C_0(A)$ the set of continuous functions that vanish at infinity, i.e.

$$C_0(A) := \{ f \in C(A) : \forall \varepsilon > 0, \exists C \text{ compact s.t. } |f| \le \varepsilon \text{ on } A \setminus C \}.$$

When $A$ is compact, $C_0(A) = C(A)$.

**Definition 2.1** (Radon measure on $A \subset \mathbb{R}^p$). Let $A \subset \mathbb{R}^p$ be locally compact. The space of (real-valued) Radon measures $M(A)$ is defined as the dual of $(C_0(A), \|\cdot\|_\infty)$. Its dual norm is the total variation norm:

$$\|\mu\|_{\text{TV}} = \sup_{\substack{\eta \in C_0(A) \\ \|\eta\|_{\infty} \le 1}} \int_A \eta(x) d\mu(x).$$

We denote $M(A)^+$ the set of nonnegative measures: $\mu \in M(A)^+$ if for all $\eta \in C_0(A)$ such that $\eta \ge 0$, $\int_A \eta d\mu \ge 0$. A discrete mixture measure, or a sparse measure, is a measure that can be expressed as a finite weighted sum of Dirac measures:

$$\mu = \sum_{j=1}^{s} a_j \delta_{x_j}$$

where $s \ge 1$, $a_1, \dots, a_s \in \mathbb{R}$, $x_1, \dots, x_s \in A$. Remark that $\|\mu\|_{\text{TV}} = \sum_{j=1}^{s} |a_j|$ when the $(x_j)_j$ are distinct. Other details about the functional framework can be found in Appendix A.

## 2.2 Model and estimator

**Model** Our aim is to recover the target measure encoding the Gaussian Mixture distribution:

$$\mu^0 = \sum_{j=1}^{s} a_j^0 \delta_{x_j^0} \quad \text{where} \quad a_1^0, \dots, a_s^0 > 0, \sum_{j=1}^{s} a_j^0 = 1.$$

We assume that the location-scale parameters $(x_1^0, \dots, x_s^0)$ are distinct points, each of the form $x_j^0 = (t_j^0, u_j^0)$ with $t_j^0 = (t_{j,k}^0)_{k=1}^d \in \mathbb{R}^d$ and $u_j^0 = (u_{j,k}^0)_{k=1}^d \in [u_{\min}, +\infty)^d$ where $u_{\min}$ is some positive lower bound on $(u_{j,k}^0)_k$. The parameters $t_j^0$ and $u_j^0$ represent respectively the mean and square root of the diagonal covariance of a Gaussian component with weight $a_j^0$.
We denote $\varphi$ the standard Gaussian density in $\mathbb{R}$, and $\sigma := \mathcal{F}[\varphi] = e^{-\frac{2}{2}}$ its associate Fourier transform. We can rewrite $f^0 = \Phi\mu^0$ where $\Phi : \mathcal{M}(\mathbb{R}^d \times [u_{\min}, +\infty)^d) \to L^2(\mathbb{R}^d)$ is the linear operator defined by

$$ \Phi\mu : z \in \mathbb{R}^d \mapsto \int_{\mathbb{R}^d \times [u_{\min}, +\infty)^d} \prod_{k=1}^{d} \frac{1}{u_k} \varphi \left( \frac{z_k - t_k}{u_k} \right) d\mu(t, u) \quad \forall \mu \in \mathcal{M}(\mathbb{R}^d \times [u_{\min}, +\infty)^d). $$

Our aim is to recover $\mu^0$ from a $n$-sample $X_1, \dots, X_n \stackrel{i.i.d.}{\sim} f^0$ (see (1)), considering that the weights $\{a_j^0\}_{j=1}^n$ the location-scale parameters $\{x_j^0\}_{j=1}^n$ and the sparsity index $s$ are unknown. The empirical density associated with our sample $X_1, \dots, X_n$ is defined as

$$ \hat{f}_n := \frac{1}{n} \sum_{i=1}^{n} \delta_{X_i}. $$

This empirical measure is smoothed, hence allowing for a comparison with any prediction $\Phi\mu$. We will use a Gaussian convolution with smoothing parameter $\tau > 0$, introducing

$$ \lambda : z \in \mathbb{R}^d \mapsto \frac{e^{-\frac{\|z\|^2}{2\tau^2}}}{(2\pi\tau^2)^{d/2}} \quad \text{and} \quad \Lambda := \mathcal{F}[\lambda] = e^{-\frac{\tau^2 \|z\|^2}{2}}. $$

Then, we set

$$ L \circ \hat{f}_n = \lambda * \hat{f}_n \quad \text{and} \quad L \circ f = \lambda * f \quad \forall f \in L^2(\mathbb{R}^d). $$

The term $L \circ \hat{f}_n$ is related to Kernel Density Estimation (KDE, [Tsybakov, 2008]). The choice of $\tau$ is important if we want to estimate $f^0$ in addition to $\mu^0$, and will be discussed in Section 4. Our setting differs in this respect from super-resolution [Candès and Fernandez-Granda, 2014], where the analogous parameter $\lambda_c$ is imposed by the experimental conditions (namely, the frequency cut-off in that paper).

**Hilbert space for the data fidelity term** The Hilbert space $\mathcal{L}$ in which we compare the observation and the prediction is the RKHS associated with $\lambda$. We define it using Mercer's theorem:

$$ \mathcal{L} := \left\{ f : \mathbb{R}^d \to \mathbb{R} : f \in L^2(\mathbb{R}^d) \text{ s.t. } \|f\|_{\mathcal{L}}^2 = \int_{\mathbb{R}^d} \frac{|\mathcal{F}[f](\xi)|^2}{\mathcal{F}[\lambda](\xi)} d\xi < +\infty \right\}, $$

with dot product

$$ \forall f, g \in \mathcal{L}, \quad \langle f, g \rangle_{\mathcal{L}} = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} \frac{\overline{\mathcal{F}[f](\xi)} \mathcal{F}[g](\xi)}{\mathcal{F}[\lambda](\xi)} d\xi. \quad (5) $$

**Problem** We will search the target measure over the space of nonnegative Radon measures $\mathcal{M}(\mathcal{X})^+$ where $\mathcal{X}$ is a compact set of $\mathbb{R}^d \times [u_{\min}, +\infty)^d$. We will restrict the possible values of the covariance even further later.

The BLASSO problem constructed in [De Castro et al., 2021a] can be generalized to our setting. We may consider the following optimization problem:

$$ \arg\min_{\mu \in \mathcal{M}(\mathcal{X})^+} J(\mu) \quad \text{where} \quad J(\mu) := \frac{1}{2} \| L \circ \hat{f}_n - L \circ \Phi \mu \|_{\mathcal{L}}^2 + \kappa \| \mu \|_{\text{TV}}, \quad \kappa > 0. \quad (\bar{P}) $$

The data fidelity term $\frac{1}{2} \| L \circ \hat{f}_n - L \circ \Phi \mu \|_{\mathcal{L}}^2$ compares smooth versions of the empirical density to any candidate for the predicted density function $\Phi\mu$. The regularization $\kappa \| \mu \|_{\text{TV}}$ promotes the sparsity of the solution. The choice of $\kappa$ will be discussed in Section 3: a balance must be found between the data fidelity term and the regularization term.

However, we do not work with the loss $J$ defined by $(\bar{P})$. Indeed, the model we consider is more complex than that of [De Castro et al., 2021a]. Specifically, the addition of diagonal covariances in the parametrization of $\mu$ renders the problem $(\bar{P})$ unsuitable for the use of standard proof techniques from the BLASSO literature [Poon, 2019]. In fact, a key object is the associated kernel $\bar{K}$ describing the correlation between 2 features, i.e.

$$ \bar{K}(x, x') := \langle L \circ \Phi \delta_x, L \circ \Phi \delta_{x'} \rangle_{\mathcal{L}} $$

for $x, x' \in \mathcal{X}$. This kernel is not normalized, namely $\bar{K}(x, x)$ is not constant (it depends on $u$), which causes issues in the investigation of recovery guarantees for a solution of $(\bar{P})$.
**Reparametrization** To address this, we reparametrize the problem to work with a normalized kernel. This is achieved by introducing a positive weighting function $W(x)$. First notice that if $W: x \in \mathcal{X} \mapsto W(x) \in \mathbb{R}_+^*$ is continuous, then as $\mathcal{X}$ is compact, $\frac{1}{W}$ is continuous and bounded on $\mathcal{X}$. Hence we can reparametrize the weights of a measure by $\frac{1}{W}$. We will work with

$$W(x) := \prod_{k=1}^{d} (2\pi)^{-1/4} (2u_k^2 + \tau^2)^{-1/4} \quad \forall x = ((t_1, \dots, t_d), (u_1, \dots, u_d)) \in \mathbb{R}^d \times [u_{\min}, +\infty)^d. \quad (6)$$

This function allows us to renormalize the kernel (see Section 5). The BLASSO problem we consider from now on is therefore

$$\min_{\mu \in \mathcal{M}(\mathcal{X})^+} J_W(\mu) \quad \text{where} \quad J_W(\mu) := \frac{1}{2} \left\| L \circ \hat{f}_n - L \circ \Phi \left( \frac{\mu}{W} \right) \right\|_{\mathbb{L}}^2 + \kappa \|\mu\|_{\text{TV}}, \quad (\mathcal{P}_{\kappa})$$

with $\kappa > 0$ and where we define, for all $\mu \in \mathcal{M}(\mathbb{R}^d \times [u_{\min}, +\infty)^d)$, $\frac{\mu}{W}$ as the measure such that

$$\forall \eta \in C_0(\mathbb{R}^d \times [u_{\min}, +\infty)^d), \quad \int_{\mathbb{R}^d \times [u_{\min}, +\infty)^d} \eta \, d\left(\frac{\mu}{W}\right) = \int_{\mathbb{R}^d \times [u_{\min}, +\infty)^d} \frac{1}{W(x)} \eta(x) \, d\mu(x).$$

We can show that the optimization problem $(\mathcal{P}_{\kappa})$ has a solution (see the appendix, Proposition B.1).

**Estimator** According to the problem $(\mathcal{P}_{\kappa})$, we consider an estimator $\mu_{n,\omega}^*$ defined as

$$\mu_{n,\omega}^* \in \arg \min_{\mu \in \mathcal{M}(\mathcal{X})^+} J_W(\mu), \quad (7)$$

Most of our results also hold for an approximate solution, and we will use the notation $\hat{\mu}_{n,\omega}$ for our estimator when it suffices that

$$\hat{\mu}_{n,\omega} \in \{\mu \in \mathcal{M}(\mathcal{X})^+ : J_W(\mu) \le J_W(\mu_\omega^0)\}, \quad (8)$$

where $\mu_\omega^0$ is the weighted measure defined as

$$\mu_{\omega}^{0} = \sum_{j=1}^{s} \omega_{j}^{0} \delta_{x_{j}^{0}} \quad \text{with} \quad \omega_{j}^{0} = W(x_{j}^{0}) a_{j}^{0} \quad \forall j \in \{1, \dots, s\}. \quad (9)$$

Remark that $\mu_{n,\omega}^*$ satisfies the relaxed condition (8). This flexibility allows us to obtain the guarantees presented in Sections 3 and 4 without requiring the exact minimization of $J_W$—this is computationally advantageous. While algorithmic details are beyond the scope of this paper, they will be briefly discussed in Section 7 (see also Chizat, 2022; De Castro et al., 2023). In contrast to other sections, the results of Section 6 specifically require an exact solution $\mu_{n,\omega}^*$.

Remark 2.1. Historically, the BLASSO method has been applied to recover a sparse target measure that is not necessarily a probability measure. As $\mu^0$ is a probability measure in our setting, one may instead consider the constrained problem

$$\arg \min_{\mu \in \mathcal{M}(\mathcal{X})^+, \|\mu\|_{\text{TV}}=1} \frac{1}{2} \left\| L \circ \hat{f}_n - L \circ \Phi \mu \right\|_{\mathbb{L}}^2. \quad (\mathcal{P}_C)$$

Such an approach deserves some comments. First, the unconstrained programs $(\mathcal{P}_{\kappa})$ and $(\bar{\mathcal{P}})$ belong to a family of convex potentials studied in optimization on the space of measures. A line of work, e.g. [Chizat, 2022; De Castro et al., 2023], shows convergence of various gradient descent strategies (particle, Wasserstein gradient flows, stochastic gradient), while such results do not exist for the constrained program $(\mathcal{P}_C)$ to the best of our knowledge. Second, for the program that we deal with (see $(\mathcal{P}_{\kappa})$), the reparametrized target $\mu_\omega^0$ does not verify $\|\mu_\omega^0\|_{\text{TV}} = 1$, and enforcing the constraint $\|\frac{\hat{\mu}_{n,\omega}}{W}\|_{\text{TV}} = 1$ is an obstacle to establishing guarantees.

# 3 Estimation guarantees

In this section, we give guarantees for the recovery of $\mu_\omega^0$ and $\mu^0$ using the estimator $\hat{\mu}_{n,\omega}$ introduced in (8). We provide bounds on the difference in mass assigned by the estimator and the target measure over relevant regions of the space. Our guarantees hold under the existence of so-called non-degenerate dual certificates, which are the key objects for analyzing the properties of the BLASSO estimator [Candès and Fernandez-Granda, 2014; Duval and Peyré, 2015; De Castro et al., 2021a; Poon et al., 2023]. In a general context, they are functions interpolating the signs of the particles of the sparse target measure, with some prescribed smoothness and shape constraints. We investigate the construction of such objects in Section 5.
## 3.1 Non-degenerate dual certificates

In this section, we define non-degenerate certificates and introduce their properties, closely related to the optimization problem at hand. These certificates are linked to the feature map, which is defined as the linear operator

$$ \Psi: \mathcal{M}(\mathbb{R}^d \times [u_{\min}, +\infty)^d) \to \mathbb{L}, \mu \mapsto L \circ \Phi \frac{\mu}{W}. $$

Remark that the loss function of ($\mathcal{P}_\kappa$) can be rewritten as $J_W: \mu \in \mathcal{M}(\mathcal{X})^+ \mapsto \frac{1}{2} \|L \circ \hat{f}_n - \Psi \mu\|_{\mathbb{L}}^2 + \kappa \|\mu\|_{\text{TV}}$. The adjoint operator of $\Psi$ restricted to $\mathcal{M}(\mathcal{X})$ verifies

$$ \langle \Psi \mu, p \rangle_{\mathbb{L}} = \langle \Psi^* p, \mu \rangle_{\mathcal{C}(\mathcal{X}), \mathcal{M}(\mathcal{X})} = \int \Psi^* p \, d\mu \quad \forall p \in \mathbb{L} \qquad (10) $$

for all $\mu \in \mathcal{M}(\mathcal{X})$. In particular, for all $x \in \mathcal{X}$, $[\Psi^* p](x) = \langle p, \Psi \delta_x \rangle_{\mathbb{L}}$.

Dual certificates are continuous functions of the form $\Psi^*p$ for some $p \in \mathbb{L}$. Our estimation guarantees are based on controls of these certificates on regions defined by $\mu^0$ and some appropriate distance. In [De Castro et al., 2021a] the Euclidean distance is used. The latter is however not adapted to the geometry of our problem, as we deal with unknown covariances. In particular, the associated kernel is not translation-invariant, see Section 5. [Poon et al., 2023] work alternatively with the Fisher-Rao distance, which is however impractical to use in our context. In fact, interactions between the kernel and the Fisher-Rao appear to be quite intricate to manage. A key aspect of our contribution is to introduce greater flexibility in the control of certificates. We relax the assumptions from [Poon et al., 2023] which are challenging to verify in our context, by expressing controls using a specific distance on regions defined by a different semi-distance. More specifically, we use the semi-distance $d$ defined by

$$ d(x, x')^2 = \sum_{k=1}^{d} \left( \frac{(t_k - t'_k)^2}{u_k'^2 + u_k^2 + \tau^2} + \ln \left( \frac{u_k'^2 + u_k^2 + \tau^2}{\sqrt{2u_k^2 + \tau^2} \sqrt{2u_k'^2 + \tau^2}} \right) \right) \quad \forall x, x' \in \mathbb{R}^d \times [u_{\min}, +\infty)^d. \quad (11) $$

This semi-distance is symmetric, nonnegative, verifies $d(x, x') = 0 \iff x = x'$ for all $x, x' \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$, but the triangle inequality does not hold.

We use this semi-distance because it is closely tied to our problem and to the feature map. Section 5 makes this connection with the optimization problem ($\mathcal{P}_\kappa$) precise. Moreover, $d$ rescales the separation between means by the corresponding variances, which is statistically natural: two Gaussian components with identical means but large variances are effectively “closer” than those with small variances.

Controls of the certificates also depend on some distance $\mathfrak{d}_g$ on $\mathcal{X}$, left unspecified in this section. We will see in Section 5 that the Fisher-Rao distance is well-suited for this problem. To establish the existence of non-degenerate dual certificates, we ensure compatibility between the semi-distance and $\mathfrak{d}_g$ (see Section 5). Up to the semi-distance $d$, the elements displayed in this section are quite generic in the BLASSO literature and have been discussed in various contexts [Duval and Peyré, 2015; De Castro et al., 2021a; Poon, 2019].

First, we introduce the so-called far and near regions.

**Definition 3.1 (Near and far regions).** For $r > 0$, we define the near regions associated with each particle $x_j^0$, $j \in \{1, \dots, s\}$ as

$$ \mathcal{X}_j^{\text{near}}(r) := \{x \in \mathcal{X} : d(x, x_j^0) \le r\} $$

and the far region

$$ \mathcal{X}^{\text{far}}(r) := \mathcal{X} \setminus \mathcal{X}^{\text{near}}(r) \quad \text{where} \quad \mathcal{X}^{\text{near}}(r) := \bigcup_j \mathcal{X}_j^{\text{near}}(r). $$

The jth near region identifies the points in $\mathcal{X}$ that are close to $x_j^0$, the jth particle of $\mu^0$. To ensure accurate recovery, we require the existence of a (global) dual certificate $\eta$ that satisfies non-degeneracy conditions. These conditions are adapted below from [Poon et al., 2023, Definition 2], with modifications to account for the definition of near and far regions, all based on the semi-distance $d$ rather than the Fisher-Rao distance.

**Definition 3.2 (Global non-degenerate certificate).** Let $\varepsilon_0, \varepsilon_2, r \in \mathbb{R}_+^*$. A function $\eta \in \text{Im}(\Psi^*)$ is a $(\varepsilon_0, \varepsilon_2, r)$-non-degenerate certificate associated with the measure $\mu^0$ if

1. $\eta(x_j^0) = 1$ for all $j = 1, \dots, s$,

2. $|\eta(x)| \le 1 - \varepsilon_0$ for all $x \in \mathcal{X}^{\text{far}}(r)$,

3. $\eta(x) \le 1 - \varepsilon_2 \mathfrak{d}_g(x, x_j^0)^2$ for all $x \in \mathcal{X}_j^{\text{near}}(r)$, $j \in \{1, \dots, s\}$.
A global non-degenerate certificate is thus a function whose regularity is prescribed (it must be in Im($\Psi^*$)), which interpolates (item 1 of Definition 3.2) and localizes (items 2 and 3) the particles of the target $\mu^0$.

We will make use of additional non-degenerate certificates, localizing the jth particle of $\mu^0$.

**Definition 3.3** (Local non-degenerate certificates). Let $\tilde{\varepsilon}_0, \tilde{\varepsilon}_2, r \in \mathbb{R}_+^*$. For $j \in \{1, \dots, s\}, \eta_j \in \text{Im}(\Psi^*)$ is a $(\tilde{\varepsilon}_0, \tilde{\varepsilon}_2, r)$-non-degenerate certificate for the jth near region if

1. $\eta_j(x_j^0) = 1$ and $\eta_j(x_i^0) = 0$ for all $i \in \{1, \dots, s\}$ such that $i \neq j$,

2. $|\eta_j(x)| \le 1 - \tilde{\varepsilon}_0$ for all $x \in \mathcal{X}^{far}(r)$,

3. $|1 - \eta_j(x)| \le \tilde{\varepsilon}_2 \mathfrak{D}_g(x, x_j^0)^2$ for all $x \in \mathcal{X}_j^{near}(r)$,

4. $|\eta_j(x)| \le \tilde{\varepsilon}_2 \mathfrak{D}_g(x, x_i^0)^2$ for all $x \in \mathcal{X}_i^{near}(r)$, for all $i \in \{1, \dots, s\}$ such that $i \neq j$.

Note that the near and far regions along with $\Psi^*$ depend on the parameter $\tau$ appearing in $\lambda$. The distance $\mathfrak{D}_g$ will also depend on it. For the estimation of $\mu^0$, the choice of $\tau$ is related via $(x_1^0, \dots, x_s^0)$ to the possibility of constructing non-degenerate certificates (see Section 5), and influences the convergence rate (see Theorem 3.1). For the estimation of the density $f^0$, this choice is more critical: we suggest possible values for $\tau$ in Section 4.

## 3.2 Error bounds

We extend here the results of [De Castro et al., 2021a] to our setting, where both means and covariances of the mixture model are unknown. Providing estimation guarantees requires to bound the error we made by approximating $L \circ f^0$ from our random observations (see Lemma 3.1). We give results taken in expected value over $X_1, \dots, X_n \stackrel{\text{i.i.d.}}{\sim} f^0$.

**Lemma 3.1** (Control of the noise. [De Castro et al., 2021b, Lemma 3]). We define the so-called noise term $\Gamma_n$ as

$$ \Gamma_n := L \circ \hat{f}_n - L \circ \Phi \mu^0. \quad (12) $$

We have

$$ \mathbb{E} \left[ \| \Gamma_n \|_{\mathbb{L}}^2 \right] \le \frac{4 \int_{\mathbb{R}^d} \Lambda}{(2\pi)^d n} = \frac{4}{(2\pi)^{d/2} \tau^d n} =: \rho_n^2. $$

Moreover, universal constants $C_\Gamma, \tilde{C}_\Gamma > 0$ exist such that

$$ \mathbb{E} \left[ \| \Gamma_n \|_{\mathbb{L}}^4 \right] \le \tilde{C}_\Gamma \rho_n^4 $$

and

$$ \forall \rho > 0, \quad \|\Gamma_n\|_{\mathbb{L}}^2 \le \rho \frac{C_\Gamma^2 \int_{\mathbb{R}} \Lambda}{n(2\pi)^d} = \rho \frac{C_\Gamma^2}{n\tau^d(2\pi)^{d/2}} \quad \text{with probability greater than } 1 - C_\Gamma e^{-\rho}. $$

The proof is given in Appendix C. The constants $C_\Gamma, \tilde{C}_\Gamma$ in the previous lemma can be made explicit (see [Houdré and Reynaud-Bouret, 2003, Theorem 3.1]).

The control of any estimator $\hat{\mu}_{n,\omega}$ satisfying (8) also involves the non-degenerate dual certificates introduced in Definitions 3.2 and 3.3 above. Our recovery guarantees are based on the following assumption.

**Assumption 1.** There exists $\eta = \Psi^*p$ a global $(\varepsilon_0, \varepsilon_2, r)$-non-degenerate certificate associated with $\mu^0$, and local $(\tilde{\varepsilon}_0, \tilde{\varepsilon}_2, r)$-non-degenerate certificates $\eta_j = \Psi^*p_j$ ($j \in \{1, \dots, s\}$) for each near region. Moreover, there exists $c_p > 0$ such that $\|p\|_{\mathbb{L}}^2 \le c_p s$ and $\|p_j\|_{\mathbb{L}}^2 \le c_p$ for $j = 1, \dots, s$.

The dependence on $s$ for $\|p\|_{\mathbb{L}}^2$ is quite natural (related to our construction of certificates, see Proposition I.1 in the appendix) and already appears in the literature, c.f. [Poon, 2019]. In Section 5, we show that Assumption 1 holds under some conditions on $\mathcal{X}$ and $\mu^0$ by explicitly constructing certificates. In particular, we consider $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, u_{\max}]^d$ and we impose a minimal separation between the particles of $\mu^0$ depending on $s, u_{\min}, u_{\max}, d, \tau$. We provide fixed values for $r, \varepsilon_i, \tilde{\varepsilon}_i, c_p$, depending only on the dimension $d$. For a precise statement, see Theorem 5.1 below.

Now, we have all the ingredients to provide our first results concerning the performances of our estimator. We present a control of the mass of $\hat{\mu}_{n,\omega}$ on the near and far regions (Definition 3.1) in the next theorem, whose proof is displayed in Appendix D.

**Theorem 3.1** (Estimation error). Assume that Assumption 1 holds. Setting $\kappa = \frac{\rho_n}{\sqrt{c_p}}$ in $(\mathcal{P}_\kappa)$, we have the following controls on $\hat{\mu}_{n,\omega}$.
1. Control over the mass of the estimator on the far region:

$$ \mathbb{E}[\hat{\mu}_{n,\omega}(\mathcal{X}^{\text{far}}(r))] \leq \frac{\sqrt{c_p}}{2\varepsilon_0} \rho_n(1+\sqrt{s})^2. $$

2. Accuracy of the mass reconstruction: for any $j \in \{1, \dots, s\}$,

$$ \mathbb{E}[|\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r))|] \leq 2\sqrt{c_p}\rho_n(1+\sqrt{s}) + \max\left\{\frac{1-\tilde{\varepsilon}_0}{\varepsilon_0}, \frac{\tilde{\varepsilon}_2}{\varepsilon_2}\right\} \frac{\sqrt{c_p}}{2}\rho_n(1+\sqrt{s})^2. $$

3. Stability of the mass:

$$ -2\sqrt{c_p s} \rho_n (1 + \sqrt{s}) \leq \mathbb{E}[\|\hat{\mu}_{n,\omega}\|_{\text{TV}}] - \|\mu_{\omega}^0\|_{\text{TV}} \leq \frac{\sqrt{c_p}}{2} \rho_n. $$

We present below several remarks on this result and its proof.

Remark 3.1 (A Bregman divergence approach). The non-degenerate certificates from Assumption 1 allow us to derive recovery guarantees from the control of the so-called Bregman divergence, defined for $\eta \in \mathcal{C}(\mathcal{X})$ by

$$ D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^{0}) = \|\hat{\mu}_{n,\omega}\|_{\text{TV}} - \|\mu_{\omega}^{0}\|_{\text{TV}} - \int_{\mathcal{X}} \eta d(\hat{\mu}_{n,\omega} - \mu_{\omega}^{0}). \quad (13) $$

In particular, if $\eta$ is a global non-degenerate certificate, we get $\int_{\mathcal{X}} \eta d\mu_{\omega}^0 = \|\mu_{\omega}^0\|_{\text{TV}}$. In such a case, the Bregman divergence is nonnegative: $D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0) = \int_{\mathcal{X}} (1-\eta) d\hat{\mu}_{n,\omega} \geq 0$. The proof of Theorem 3.1 is based on lower and upper bounds on the Bregman divergence with $\eta$ from Assumption 1. First, using that $\eta = \Psi^*p$ with the Cauchy-Schwarz inequality, we get

$$ \left| \int_{\mathcal{X}} \eta d(\hat{\mu}_{n,\omega} - \mu_{\omega}^{0}) \right| = \left| \langle p, L \circ \Phi(\hat{\mu}_{n,\omega} - \mu_{\omega}^{0}) \rangle_{\mathbb{L}} \right| \leq \|p\|_{\mathbb{L}} \|L \circ \Phi(\hat{\mu}_{n,\omega} - \mu_{\omega}^{0})\|_{\mathbb{L}}. $$

We control this quantity by using the inequality $J_W(\hat{\mu}_{n,\omega}) \leq J_W(\mu_{\omega}^0)$ (recall (8)), along with the control on the noise term established in Lemma 3.1. This leads to the following upper bound on the Bregman divergence:

$$ \mathbb{E}[D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^{0})] \leq \frac{\sqrt{c_p}}{2} \rho_n(1 + \sqrt{s})^2. $$

At the same time, since $\eta$ satisfies Definition 3.2, we get

$$ D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^{0}) = \int (1-\eta) d\hat{\mu}_{n,\omega} \geq \varepsilon_0 \hat{\mu}_{n,\omega}(\mathcal{X}^{\text{far}}(r)) + \varepsilon_2 \sum_{j=1}^{s} \int_{\mathcal{X}_j^{\text{near}}(r)} \mathfrak{d}_g(x, x_j^0)^2 d\hat{\mu}_{n,\omega}(x) $$

from which we deduce the control of the estimator on the far region. The local dual certificates enable us to retrieve the control on the near regions, using again the lower bound on the Bregman divergence.

Remark 3.2 (Choice of $\kappa$, dependence on $s$). The regularization parameter $\kappa$ is chosen as $\frac{\rho_n}{\sqrt{c_p}}$ in Theorem 3.1. In particular, it does not depend on the sparsity index $s$ of $\mu^0$, unknown in practice—we call this choice $s$-agnostic. As a result, the user does not require prior knowledge of the target measure to select the regularization parameter: the proposed value can be directly used to solve the BLASSO numerically. The $s$-dependent choice $\kappa = \frac{\rho_n}{\sqrt{s\varepsilon_p}}$ results in better rates for the estimation (i.e., for the bounds presented in Theorem 3.1)—linear on $\sqrt{s}$ rather than on $s$ (see Equation (33) in Appendix). In all cases, the estimation error depends on $s$; recovering a mixture with a larger number of components incurs a higher estimation cost. To conclude this discussion, we stress that these possible choices for $\kappa$ are proposed according to theoretical considerations. For practical applications, this regularization parameter can be calibrated via cross-validation.

Remark 3.3 (Soft-thresholding effect). It is known that the LASSO estimator has a soft-thresholding effect, e.g. [Friedman et al., 2007]. It suggests that the $\ell_1$-regularization is biased in the sense that each of the $s$ weight components is under-estimated by an additive factor proportional to the regularization term $\kappa$. Having $s$ components, we expect this bias to be of the order of $s\kappa$. Our result in (3) (Theorem 3.1) is aligned with this comment since

$$ \left| \mathbb{E}[\|\hat{\mu}_{n,\omega}\|_{\text{TV}}] - \|\mu_{\omega}^0\|_{\text{TV}} \right| = \mathcal{O}(s\kappa), $$

up to a constant that may depend on $c_p$.
Remark 3.4 (Choice of $\tau$). Although this section focuses on the estimation of $\mu^0$, an additional objective that can be pursued in parallel is the estimation of the associated density $f^0$—a task we refer to as prediction. This prediction objective influences the choice of the smoothing parameter. As we will see in Section 4, when $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, +\infty)^d$ the choice $\tau = \sqrt{2}\frac{u_{\min}}{\sqrt{\ln n}}$ results in almost parametric convergence rates for both the estimation of $f^0$ (Theorem 4.1) and the estimation of $\mu^0$. In fact, under Assumption 1, Theorem 3.1 entails that with $\tau = \sqrt{2}\frac{u_{\min}}{\sqrt{\ln n}}$, setting $\kappa = \frac{\rho_n}{\sqrt{c_p}}$, keeping only the dependence on $n$ and $s$ we have

$$ \mathbb{E} [|\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r))|] \lesssim \frac{s(\ln n)^{d/4}}{\sqrt{n}} \quad \forall j = 1, \dots, s. $$

If we do not need to predict the density $f^0$, we can alternatively take $0 < \tau \le u_{\min}$ fixed. This leads to a rate of $\frac{s}{\sqrt{n}}$ for the estimation. Note that Assumption 1 is more difficult to check when $\tau$ is large, c.f. Theorem 5.1. Also note that the near and far regions depend on $\tau$ through the semi-distance $d$.

## 3.3 Effective near regions

Theorem 3.1 describes the performances of the BLASSO estimator in our setting. It provides control over the proximity between $\hat{\mu}_{n,\omega}$ and $\mu^0$ on associated far and near regions. However, our initial target is $\mu^0$ instead of $\mu^0$. Using Theorem 3.1 and (9), we can easily give a bound for $\mathbb{E} [|a_j^0 - \frac{\hat{\mu}_{n,\omega}}{W(x_j^0)}(\mathcal{X}_j^{\text{near}}(r))|]$, but we cannot access $W(x_j^0)$ without knowledge of $\mu^0$. In this context, a fair estimator of $\mu^0$ is the renormalized measure $\frac{\hat{\mu}_{n,\omega}}{W}$. Providing recovery guarantees for $\frac{\hat{\mu}_{n,\omega}}{W}$ requires to control $W$ on the near regions.

Another restrictive aspect of Theorem 3.1 is that we consider controls on regions with a fixed radius $r$ (which is directly related to the dual certificate). We would like to locate the mass of $\hat{\mu}_{n,\omega}$ more precisely.

However, we can overcome these initial limitations by providing controls of the estimator on

$$ \mathcal{X}_j^{\text{near}}(r_e) := \{x \in \mathcal{X} : d(x, x_j^0) \le r_e\} \quad (14) $$

for $r_e \le r$. Such regions are called effective near regions, and have been introduced by [De Castro et al., 2025]. To extend controls on $\mathcal{X}_j^{\text{near}}(r)$ to controls on $\mathcal{X}_j^{\text{near}}(r_e)$, the choice of $\mathfrak{d}_g$ plays an important role. It should be indeed compatible with our semi-distance $d$ in a sense which is made precise in Proposition 3.1 below.

**Proposition 3.1.** Let $r > 0$. Assume that Assumption 1 holds. Assume that there exists $\tilde{\varepsilon}_3 > 0$ such that for all $0 < r_e \le r$, for all $j \in \{1, \dots, s\}$, we have

$$ \mathfrak{d}_g(x_j^0, x)^2 \ge \frac{r_e^2}{\tilde{\varepsilon}_3} \quad \forall x \in \mathcal{X}_j^{\text{near}}(r) \setminus \mathcal{X}_j^{\text{near}}(r_e). \quad (15) $$

Then, for $\kappa = \frac{\rho_n}{\sqrt{c_p}}$ and for any $0 < r_e \le r$,

$$ \mathbb{E} [|\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))|] \le 2\sqrt{c_p}\rho_n(1 + \sqrt{s}) + \max \left\{ \frac{1 - \tilde{\varepsilon}_0}{\varepsilon_0}, \frac{1}{\varepsilon_2} \left( \frac{\tilde{\varepsilon}_3}{r_e^2} + \tilde{\varepsilon}_2 \right) \right\} \rho_n \frac{\sqrt{c_p}}{2} (1 + \sqrt{s})^2. \quad (16) $$

The proof is similar to that of Theorem 3.1, and is presented in Appendix F.1.

We check the assumption (15) in Lemma H.5 (in Appendix), for $\mathfrak{d}_g$ chosen as the Fisher-Rao distance. For the specific choice of $r$ which is made in Section 5, (depending on $d$), we show in Lemma 5.4 that $\tilde{\varepsilon}_3$ can be chosen as a constant not depending on $d$. Note that, as we decrease the size of near region $r_e$ towards zero, the bound (16) grows as $\frac{s\rho_n}{r_e^2}$, omitting the dependence on $c_p, \varepsilon_0, \varepsilon_2, \tilde{\varepsilon}_0, \tilde{\varepsilon}_2$ and $\tilde{\varepsilon}_3$.

Remark 3.5. We can make choices for $r_e$ that depend on $n$ in Proposition 3.1. We give the associated bound for $\mathbb{E} [|\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))|]$ (see (16)), that holds under the assumptions of Proposition 3.1 and with $n$ large enough such that $r_e \le r$. We only keep the dependence on $n, \tau$ and $s$. Let $\alpha > 0$.

* For $r_e = (\ln n)^{-\alpha}$,

$$ \mathbb{E} [|\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}((\ln n)^{-\alpha}))|] \lesssim \frac{s(\ln n)^{2\alpha}}{\sqrt{n\tau^{d/2}}}. $$

* For $r_e = n^{-\alpha}$, $\alpha < \frac{1}{4}$,

$$ \mathbb{E} [|\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(n^{-\alpha}))|] \lesssim \frac{s}{n^{\frac{1}{2}-2\alpha} \tau^{d/2}}. $$

These choices provide an overview of possible convergence rates, and show the trade-off between localization and control of the mass.
The above proposition provides, as Theorem 3.1, guarantees on the estimation on $\mu_{\omega}^0$ rather than on $\mu^0$. However, by combining Proposition 3.1 and a control of the function $W$ on the effective regions, we can derive guarantees for the estimation of $\mu^0$ by $\hat{\mu}_n := \frac{\hat{\mu}_{n,\omega}}{W}$, as displayed in the following corollary.

**Corollary 3.1.** We work under the same assumption as Proposition 3.1. We choose $\kappa = \frac{\rho_n}{\sqrt{c_p}}$ and $r_e = n^{-1/6}$. Assume that $0 < r_e \le r$. Omitting the dependence on $d, c_p, \varepsilon_i, \tilde{\varepsilon}_i, r$ we have

$$ \mathbb{E}\left[\left|a_j^0 - \frac{\hat{\mu}_{n,\omega}}{W}(\mathcal{X}_j^{\text{near}}(n^{-1/6}))\right|\right] \lesssim \left(s\tau^{-d/2}W(x_j^0)^{-1} + a_j^0\right)n^{-1/6}. $$

The proof can be found in Appendix F.2. The difficulty that appears here is that providing control on $\frac{\hat{\mu}_{n,\omega}}{W}$ requires local control of the function $W$, which slightly deteriorates the bound for the estimation of the weights $(a_j^0)_{j=1}^s$. The rate $n^{-1/6}$ provides the best compromise for the size of the effective near regions $r_e$ when dealing with this upper bound.

**Remark 3.6.** Making the $s$-dependent choice of regularization $\kappa = \frac{\rho_n}{\sqrt{c_p s}}$, a modified version of Proposition 3.1 yields, keeping only the dependence on $r_e, \tau, s, n$,

$$ \mathbb{E}\left[\left|\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))\right|\right] \lesssim \frac{\sqrt{s}}{\tau^{d/2}\sqrt{n}r_e^2}, $$

c.f. (38) in Appendix. With this choice of $\kappa$, the result of Corollary 3.1 becomes

$$ \mathbb{E}\left[\left|a_j^0 - \frac{\hat{\mu}_{n,\omega}}{W}(\mathcal{X}_j^{\text{near}}(n^{-1/6}))\right|\right] \lesssim \left(\sqrt{s}\tau^{-d/2}W(x_j^0)^{-1} + a_j^0\right)n^{-1/6}, $$

see (42) in Appendix.

Remark that all the results given in this section are controls in expected value. Controls with high probability could also be given, using Lemma 3.1. Also note that the target $\mu^0$ is fixed (does not vary with $n$).

# 4 Prediction

In Section 3, we provided controls for the estimator $\hat{\mu}_{n,\omega}$. In this section, we look at the prediction $\Phi\hat{\mu}_n = \Phi\frac{\hat{\mu}_{n,\omega}}{W}$ made by the BLASSO of the target density $\Phi\mu^0 = f^0$, and we derive bounds for the so-called prediction error $\|\Phi(\hat{\mu}_n - \mu^0)\|_{L^2}^2$. This is a question that has been quite overlooked in the BLASSO literature. We can however mention [Butucea et al., 2024] where the prediction error is investigated in a different context and for a slightly different problem.

We establish in this section almost parametric rates for the bound on the prediction error, in two distinct regimes characterized by the value of the regularization parameter $\kappa$. With small regularization, no assumptions on $\mu^0$ are needed (Section 4.1). With stronger regularization, we show in Section 4.2 that the BLASSO achieves a good prediction when dual certificates exist (namely when Assumption 1 holds).

## 4.1 Prediction under small regularization

We can achieve good prediction results with small regularization, using control of the low frequencies of the estimated density (using Lemma 3.1) and of its high frequencies, resulting from an appropriate choice of $\tau$. These controls require an upper bound on the variance. We hence assume from now on that $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, u_{\max}]^d$.

**Proposition 4.1 (Prediction error under small regularization).** Assume that $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, u_{\max}]^d$. Choosing $\tau = \frac{\sqrt{2}u_{\min}}{\sqrt{\ln n}}$ and $\kappa = \rho_n^2$, keeping only the dependence on $n$ we have

$$ \mathbb{E}\left[\|\Phi(\hat{\mu}_n - \mu^0)\|_{L^2}^2\right] \lesssim \frac{(\ln n)^{d/2}}{n}. $$

We do not make any assumption on the existence of dual certificates in the above proposition. The proof is given in Appendix G.1. This result shows that in the small regularization regime corresponding to $\kappa = \rho_n^2$, we obtain a quasi-parametric rate (up to a log factor) for the estimation of the density $f^0$.

**Remark 4.1 (Choice of $\tau$, bounds on the variances).** Keeping $\tau$ fixed does not lead to a good prediction rate: to control the high frequencies of the predicted density with Lemma G.1 in Appendix, we need to lower bound the variances of $\mathcal{X}$ with a parameter not depending on $n$, while $\tau$ must decrease as $n$ grows. In Proposition 4.1, we also use an upper bound on the variances in $\mathcal{X}$ to control $\mathbb{E}[\|\hat{\mu}_n\|_{\text{TV}}^2]$.
**Remark 4.2 (Comparison with Kernel Density Estimation).** Note that $L \circ \hat{f}_n$ is itself a kernel density estimator of the density (see [Tsybakov, 2008]). However, the resulting rate of convergence deteriorates very quickly with the dimension: it can be shown that with $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, +\infty)^d$, setting $\tau = \frac{1}{\sqrt{\ln nn^{4+d/4}}}$, omitting the dependence on $d$ we have

$$ \mathbb{E} \left[ \left\| L \circ \hat{f}_n - \Phi \mu^0 \right\|_{L^2}^2 \right] \lesssim \frac{(\ln n)^{d/2}}{n^{\frac{4}{d+4}} u_{\min}^{d+4}}. $$

The proof is in Appendix G.2. This bound can be explained by the fact that the Gaussian kernel is not appropriate for the estimation of so-called super-smooth target densities. Nevertheless, as displayed in Proposition 4.1, an almost parametric rate is obtained for $\Phi \hat{\mu}_n$, which is based on a regularized version of $L \circ \hat{f}_n$.

**Choice of regularization** Proposition 4.1 does not take advantage on the fact that $\mu^0$ is a discrete measure. This is not surprising, as the considered regularization is very small. The choice $\kappa = \rho_n^2$ is not suited for the estimation of the target measure $\mu^0$ (see Section 3). In fact, the result displayed in Proposition 4.1 does not take into account the trade-off we would like to have between making a good prediction and guaranteeing a good estimation of the target measure $\mu^0$. In particular, choosing $\kappa = \rho_n^2$, we have no control over the proximity between $\|\mu_\omega^0\|_{TV}$ and $\|\hat{\mu}_{n,\omega}\|_{TV}$ as $n$ grows: we only know that $\mathbb{E}[\|\hat{\mu}_{n,\omega}\|_{TV}] \le \|\mu_\omega^0\|_{TV} + \frac{1}{2}$ (as $\kappa \mathbb{E}[\|\hat{\mu}_{n,\omega}\|_{TV}] \le \frac{\rho_n^2}{2} + \kappa \|\mu_\omega^0\|_{TV}$), which is enough to control the high frequencies of $\Phi \hat{\mu}_n$. For an estimation purpose, we need more regularization (i.e. a larger $\kappa$). This is exactly what has been done in Theorem 3.1 where $\kappa = \frac{\rho_n}{\sqrt{c_p}}$.

**Remark 4.3.** This property of the BLASSO estimator, i.e. obtaining a good prediction rate under small regularization, highlights a strong difference from the LASSO framework. In the LASSO setting [Tibshirani, 1996; Tibshirani, 2023], regularization plays a greater role in controlling the prediction error: to get a parametric rate, we need to use bounds on the estimation error.

## 4.2 Prediction under large regularization

Provided that non-degenerate dual certificates exist, one can obtain appropriate bounds for both estimation and prediction, with $\kappa$ chosen as in Theorem 3.1. The following theorem presents the prediction bound in this regime.

**Theorem 4.1.** Assume that $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, u_{\max}]^d$. Let $\tau^2 = \frac{2u_{\min}^2}{\ln n}$. If Assumption 1 holds, choosing $\kappa = \frac{\rho_n}{\sqrt{c_p}}$, keeping only the dependence on $s$ and $n$ we have

$$ \mathbb{E} \left[ \left\| \Phi \hat{\mu}_n - \Phi \mu^0 \right\|_{L^2(\mathbb{R}^d)}^2 \right] \lesssim \frac{s(\ln n)^{d/2}}{n}. $$

The proof is given in Appendix G.3. The upper bound on the prediction error displayed here is slightly larger than in Proposition 4.1. In particular, this bound depends linearly on the sparsity index $s$. We refer for instance to [Butucea et al., 2024] for a similar bound in a different context. We recall that this specific choice for $\kappa$ allows us to control the estimation performances of the estimator (see Theorem 3.1).

**Remark 4.4.** Similarly to the estimation task (see Remark 3.2), a choice of regularization depending on $s$ leads to a better rate. With $\kappa = \frac{\rho}{\sqrt{c_p}}$, we get $\mathbb{E} \left[ \left\| \Phi \hat{\mu}_n - \Phi \mu^0 \right\|_{L^2(\mathbb{R}^d)}^2 \right] \lesssim \frac{(\ln n)^{d/2}}{n}$ for $s = O(n(\ln n)^{d/2})$. We refer to Equation (45) in Appendix G.3.

# 5 Dual certificates

The proofs of Theorem 3.1 and Proposition 3.1 rely on the existence of non-degenerate dual certificates associated with $\mu^0$ (Assumption 1). In this section, we explain how to construct such objects and the assumptions needed.

**Connection with previous works** In the general framework of BLASSO, the objective is to recover a sparse target measure from the observed signal. A key analytical tool in this setting is the dual certificate—a smooth function associated with the underlying feature map of the problem at hand. These certificates identify the positions of the target particles by interpolating the signs of the target measure.

Dual certificates can be traced back to super-resolution [Candès and Fernandez-Granda, 2014] and minimal extrapolation [De Castro and Gamboa, 2012]. They are related to the dual solutions of the BLASSO when the observation is noiseless, c.f. [Duval and Peyré, 2015] (in our framework, the noiseless observation corresponds to $\mathbb{E}[L \circ \hat{f}_n] = L \circ f^0$). Our particular construction of dual certificates is inspired from the BLASSO literature (e.g., “vanishing derivatives pre-certificate” in [Duval and Peyré, 2015], “limit certificate” in [Poon et al., 2023]).
Although our setting differs—we observe a sample drawn from some mixture distribution—the construction of the dual certificates follows the same principles: they depend solely on the target measure $\mu^0$ and the feature map $\Psi$, and not on the observed data nor the regularization. Let us mention that other constructions exist for translation invariant kernels (e.g., *pivot certificates* in [De Castro et al., 2021a; De Castro et al., 2025]).

The existence of certificates requires assumptions on $\mu^0$, mainly on the separation between its particles. In [Poon et al., 2023], precise conditions are proposed, based on controls on the kernel. In this section, we adapt and check these assumptions.

**Contributions** Our main contribution is the construction of dual certificates for a BLASSO problem involving a kernel that is not translation invariant.

In Section 2, we have adapted the data fitting term to obtain a normalized kernel, and introduced reparametrized measures (see ($P_\kappa$)). It is however challenging to check that $K_{\text{norm}}$ satisfies the Local Positive Curvature assumption [Poon et al., 2023, Assumption 1] with the Fisher-Rao distance. We therefore consider regions for controls defined instead by the semi-distance $d$ introduced in (11). Even so, controls of the certificates are expressed with the Fisher-Rao distance, since they are carried out using Taylor expansions along the Fisher-Rao geodesics.

A technical contribution is the derivation of bounds and local controls for the Riemannian derivatives of the normalized Gaussian kernel. Additionally, we establish the compatibility between the Fisher-Rao metric and the semi-distance. It is essential to control the Fisher-Rao geodesics within the balls with respect to $d$ in order to make use of the local bounds on the kernel. Combining these elements, we derive sufficient conditions on the target $\mu^0$ that guarantee the existence of non-degenerate dual certificates. The next subsections are dedicated to the proof of the following theorem.

**Theorem 5.1** (Main result: existence of certificates under a minimal separation). Assume that $\mathcal{X}$ is a compact set of $\mathbb{R}^d \times [u_{\min}, u_{\max}]^d$. Let $s \ge 2$, $0 < \tau \le u_{\min}$ and $\{x_j^0\}_{j=1}^s \subset \mathcal{X}$. We set $\mathfrak{d}_g$ as the Fisher-Rao distance (associated with the metric $g$ defined by (20)). If

$$ \min_{i \neq j} d(x_i^0, x_j^0) \ge \max \left\{ \frac{\sqrt{u_{\max}^2 + \frac{1}{4} \frac{0.3025^2}{d}(2u_{\max}^2 + \tau^2)}}{u_{\min}} \left( \Delta + \frac{0.3025}{\sqrt{d}} \right), 2 \frac{u_{\max}}{u_{\min}} \Delta \right\} + \sqrt{d \ln \left( \frac{u_{\max}^2}{u_{\min}^2} \right)} \quad (17) $$

where $\Delta = 2\sqrt{11.9 + 3\ln(d+6.62) + \ln(s-1)}$, then Assumption 1 holds with $r = \frac{0.3025}{\sqrt{d}}$, $(\varepsilon_0, \varepsilon_2) = (\frac{0.03911}{d}, 0.06158)$, $(\tilde{\varepsilon}_0, \tilde{\varepsilon}_2) = (\frac{0.03911}{d}, \frac{\sqrt{4d^2+10d}}{2} + 0.004106)$ and $c_p = 2$.

This result stems from Theorems 5.2 and 5.3 below. The calculation of the parameters is detailed in [Giard, 2025, Section VII.3].

Theorem 5.1 ensures the existence of non-degenerate dual certificates when the particles $\{x_j^0\}_{j=1}^s$ are sufficiently separated. This separation, involving both means and covariances of the Gaussian components, is expressed via the semi-distance $d$, rather than the Fisher-Rao distance.

Note also that we give explicit constants, but a less constructive approach could be considered, based on the continuity of $K_{\text{norm}}$ and its derivatives.

## 5.1 Geometrical framework

**Kernel** Building on the work of [Duval and Peyré, 2015] and [Poon et al., 2023], we consider dual certificates of the form

$$ \eta_{\alpha, \beta} = \sum_{j=1}^{s} \alpha_j K_{\text{norm}}(x_j^0, \cdot) + \sum_{j=1}^{s} \beta_j^T \nabla_1 K_{\text{norm}}(x_j^0, \cdot) \quad (18) $$

where $\alpha_j \in \mathbb{R}$, $\beta_j \in \mathbb{R}^{2d}$ for all $j \in \{1, \dots, s\}$, and $K_{\text{norm}}$ is the real-valued kernel defined, for all $x, x' \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$, by

$$ K_{\text{norm}}(x, x') = \langle \Psi \delta_x, \Psi \delta_{x'} \rangle_{\mathbb{L}} = \prod_{k=1}^{d} (2u_k^2 + \tau^2)^{1/4} (2u_k'^2 + \tau^2)^{1/4} \frac{e^{-\frac{(t_k-t_k')^2}{2(u_k^2+u_k'^2+\tau^2)}}}{(u_k^2+u_k'^2+\tau^2)^{1/2}} \quad (19) $$

The kernel $K_{\text{norm}}$ is normalized according to the definition of $W$ (involved in $\Psi$). It is indeed easy to check that $K_{\text{norm}}(x, x) = 1$ for all $x \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$.

The gradient of $K_{\text{norm}}$ with respect to its first variable is written as $\nabla_1 K_{\text{norm}}$ ($\nabla_2 K_{\text{norm}}$ denotes the gradient with respect to the second variable).
**Remark 5.1** (Smoothness of the feature map). Note that $K_{\text{norm}}$ is $C^{\infty}$ on $(\mathbb{R}^d \times [u_{\text{min}}, +\infty)^d)^2$. By iteratively applying [Christmann and Steinwart, 2008, Lemma 4.34], it comes that

$$ (x \mapsto \Psi \delta_x) \in C^{\infty}(\mathbb{R}^d \times [u_{\text{min}}, +\infty)^d, \mathbb{L}) $$

and that $\langle \partial_1 \Psi \delta_x, \partial_2 \Psi \delta_{x'} \rangle_{\mathbb{L}} = \partial_1 \partial_2 K_{\text{norm}}(x, x')$, where $\partial_1$ (resp. $\partial_2$) denotes here any derivative w.r.t. to $x$ (resp. $x'$).

The next lemma shows the relevance of constructing $\eta$ of the form (18): such functions belong to $\text{Im}(\Psi^*)$.

**Lemma 5.1.** Let $\alpha \in \mathbb{R}^s, \beta \in \mathbb{R}^{2d \times s}$. The function $\eta_{\alpha, \beta}$ introduced in (18) verifies, for all $x \in \mathbb{R}^d \times [u_{\text{min}}, +\infty)^d$,

$$ \eta_{\alpha, \beta}(x) = \langle p_{\alpha, \beta}, \Psi \delta_x \rangle_{\mathbb{L}} \quad \text{with} \quad p_{\alpha, \beta} = \sum_{j=1}^{s} \alpha_j \Psi \delta_{x_j^0} + \sum_{j=1}^{s} \beta_j \nabla_x (\Psi \delta_{x_j^0}) \in \mathbb{L}. $$

Furthermore, $\eta_{\alpha, \beta} \in C^{\infty}(\mathbb{R}^d \times [u_{\text{min}}, +\infty)^d)$ and in particular, $\eta_{\alpha, \beta}|_{\chi} \in \text{Im}(\Psi^*)$.

The proof is an immediate consequence of Remark 5.1 and of the definition of $K_{\text{norm}}$ in (19).

The construction of a non-degenerate (global) dual certificate follows several steps: first we construct a function $\eta$ of the form (18), choosing $\alpha_j, \beta_j$ such that $\eta(x_j^0) = 1$ and $\nabla \eta(x_j^0) = 0$. This amounts to solve some linear system (see (63) in Appendix). Then we use Taylor expansions on geodesics associated with the distance $\mathfrak{d}_g$ to control $\eta$ on the near and far regions, in the spirit of [Poon et al., 2023]. The Fisher-Rao distance appears to be well suited for this purpose.

**Fisher-Rao metric** We work in a Riemannian geometry framework, and we use the Fisher-Rao metric induced by $K_{\text{norm}}$ (c.f. [Poon et al., 2023, Lemma 1]), defined for $x \in \mathbb{R}^d \times [u_{\text{min}}, +\infty)^d$ by

$$ \mathfrak{g}_x := \nabla_1 \nabla_2 K_{\text{norm}}(x, x) = \text{diag} \left( \frac{1}{2u_1^2 + \tau^2}, \dots, \frac{1}{2u_d^2 + \tau^2}, \frac{2u_1^2}{(2u_1^2 + \tau^2)^2}, \dots, \frac{2u_d^2}{(2u_d^2 + \tau^2)^2} \right). \quad (20) $$

We will use the associated norm, defined by

$$ \|v\|_x = \sqrt{v^T \mathfrak{g}_x v} \quad \forall v \in \mathbb{R}^{2d}, x \in \mathbb{R}^d \times [u_{\text{min}}, +\infty)^d. $$

Details about this metric are provided in Appendix H, together with a description of the associated geodesics. We recall that a geodesic for the metric $\mathfrak{g}$ between $x, x' \in \mathbb{R}^d \times [u_{\text{min}}, +\infty)^d$ is a piecewise continuously differentiable function $\gamma: [0, 1] \to \mathbb{R}^d \times [u_{\text{min}}, +\infty)^d$ such that $\gamma(0) = x, \gamma(1) = x'$, minimizing the quantity $\int_0^1 \|\dot{\gamma}(y)\|_{\gamma(y)} dy$. In dimension $d=1$, our Fisher-Rao geodesics share the same paths as those of the Poincaré half-plane model (c.f. Lemma H.2 in Appendix): they are portions of straight lines parallel to $\{t=0\}$ and semicircles whose origin is on $\{u=0\}$. The notation $\mathfrak{d}_g$ refers to the associated distance. It is defined by $\mathfrak{d}_g(x, x') = \int_0^1 \|\dot{\gamma}(y)\|_{\gamma(y)} dy$ with $\gamma$ the geodesic between $x, x'$. This Fisher-Rao distance is not practical to define near and far regions since it is not directly linked to the correlation between 2 features calculated with the kernel. This motivates the use of another distance on $\mathcal{X}$ to express the recovery results: we use the semi-distance defined by

$$ d(x, x') = \sqrt{-2 \ln(K_{\text{norm}}(x, x'))}, \quad (21) $$

whose expression is given by (11).

**Riemannian derivatives** Following the work of [Poon et al., 2023], we will show that the dual certificates we construct are non-degenerate (see Definitions 3.2 and 3.3) by controlling the Riemannian derivatives of $K_{\text{norm}}$. For details about this framework, see [Poon et al., 2023, p.263].

Riemannian derivatives involve the Christoffel symbols associated with $\mathfrak{g}$. We will use the notation $\Gamma^{t_k}, \Gamma^{u_k}$, and refer to (46) in Appendix H.1 for a precise definition.

**Definition 5.1.** Let $\psi \in C^2(\mathbb{R}^d \times [u_{\text{min}}, +\infty]^d)$. Let $x, x' \in \mathbb{R}^d \times [u_{\text{min}}, +\infty]^d$.

**Riemannian Hessian:** we define

$$ H^g \psi(x) = \nabla^2 \psi(x) - \sum_{k=1}^{d} \Gamma^{t_k} \partial_{t_k} \psi(x) - \sum_{k=1}^{d} \Gamma^{u_k} \partial_{u_k} \psi(x). $$

**Covariant derivatives:** Let $v, v' \in \mathbb{R}^{2d}$. We define

$$ D_0[\psi](x) = \psi(x), \quad D_1[\psi](x)[v] = v^T \nabla \psi(x), \quad D_2[\psi](x)[v, v'] = v^T H^g \psi(x)v'. $$
**Operator norms:** For $j \in \{0, 1, 2\}$, we define the operator norm

$$ \|D_j[\psi](x)\|_x := \sup_{\substack{V=[v_1, \dots, v_j] \in (\mathbb{R}^{2d})^j \\ \forall l=1, \dots, j, \|v_l\|_x \le 1}} D_j[\psi](x)[V]. $$

**Kernel derivatives and associated operators:** Let $i, j \in \{0, 1, 2\}$. We define the covariant derivative of the kernel of order $i$ with respect to the first variable $x$ and of order $j$ with respect to the second variable $x'$ by

$$ [Q] K_{\text{norm}}^{(ij)}(x, x')[V] = \langle D_i[\Psi](x)[Q], D_j[\Psi](x')[V] \rangle_{\mathbb{L}} \quad \forall Q \in (\mathbb{R}^{2d})^i, V \in (\mathbb{R}^{2d})^j. $$

The associated operator norm is

$$ \| K_{\text{norm}}^{(ij)}(x, x') \|_{x, x'} := \sup_{\substack{Q=[Q_1, \dots, Q_i] \in (\mathbb{R}^{2d})^i, V=[V_1, \dots, V_j] \in (\mathbb{R}^{2d})^j \\ \forall l \|Q_l\|_x, \|V_l\|_{x'} \le 1}} [Q] K_{\text{norm}}^{(ij)}(x, x') [V]. $$

Simplified expressions are given in Appendix (Lemma J.1).

## 5.2 Construction of dual certificates

In what follows, we give some results allowing us to adapt and check the hypothesis of [Poon et al., 2023, Theorem 2], which shows that the (global) non-degenerate certificate exists under some conditions on the kernel. The corresponding proof can be adapted (see for instance [Poon et al., 2023, Section 6.7]) to show the existence of additional certificates for the near regions, under the same conditions.

**Local positive curvature assumption** The following definition allows us to adapt [Poon et al., 2023, Assumption 1] to our framework. In particular, we require some smoothness and structural properties for the kernel that enable the construction of dual certificates, as stated below.

**Definition 5.2** (Kernel of local positive curvature with parameters $s, \Delta, r, \bar{\epsilon}_0$ and $\bar{\epsilon}_2$). Let $K$ be a real-valued normalized kernel of positive type, in $C^3((\mathbb{R}^d \times [u_{\min}, +\infty)^d)^2)$. It is said to satisfy the local positive curvature assumption (LPC) if the following holds:

1. For all $i, j \in \{0, 1, 2\}$ with $i + j \le 3$,

$$ \sup_{x, x' \in \mathbb{R}^d \times [u_{\min}, +\infty)^d} \| K^{(ij)}(x, x') \|_{x, x'} \le B_{ij} < +\infty. $$

For $i = 0, 1, 2$, we denote $B_i := 1 + B_{0i} + B_{1i}$.

2. There exists $r > 0$ such that $K$ has strictly positive curvature constants $\bar{\epsilon}_0$ and $\bar{\epsilon}_2$ with

$$ K(x, x') \le 1 - \bar{\epsilon}_0, \quad \forall x, x' \in \mathbb{R}^d \times [u_{\min}, +\infty)^d \text{ s.t. } d(x, x') \ge r, \\ -K^{(02)}(x, x')[v, v] \ge \bar{\epsilon}_2 \|v\|_{x'}^2, \quad \forall v \in \mathbb{R}^{2d}, \quad \forall x, x' \in \mathbb{R}^d \times [u_{\min}, +\infty)^d \text{ s.t. } d(x, x') < r. $$

3. There exist $s \ge 2$, $\Delta(s) > 0$ such that for all $\{x_l\}_{l=1}^s \in S_{\Delta(s)}$,

$$ \sum_{l=2}^{s} \| K^{(ij)}(x_1, x_l) \|_{x_1, x_l} \le \frac{1}{64} \min \left( \frac{\bar{\epsilon}_0(r)}{B_0}, \frac{\bar{\epsilon}_2(r)}{B_2} \right) \quad \forall (i, j) \in \{0, 1\} \times \{0, 1, 2\}, $$

$$ \text{where } S_{\Delta} := \{\{x_l\}_{l=1}^s \subset \mathbb{R}^d \times [u_{\min}, +\infty)^d : \min_{m \neq l} d(x_m, x_l) \ge \Delta\}. $$

Such kernel is said to verify the LPC with parameters $s, \Delta(s), r, \bar{\epsilon}_0(r)$ and $\bar{\epsilon}_2(r)$.

**Remark 5.2.** This definition differs on some points from [Poon et al., 2023, Assumption 1]. First, since we solve $(P_\kappa)$ for nonnegative measures, we do not have to control the negative part of the certificates, and thus do not require $r < B_{02}^{-1/2}$. Moreover, we use controls on regions defined by the semi-distance $d$ instead of $\partial_g$. We need to provide bounds on the kernel on $\mathbb{R}^d \times [u_{\min}, +\infty)^d$ and not $\mathcal{X}$, because a Fisher-Rao geodesic between 2 points of $\mathcal{X}$ could be outside $\mathcal{X}$. The use of a semi-distance, which better describes the correlation between 2 features (as expressed by the kernel), also leads to difficulties. We need to ensure that $d$ is compatible with the Fisher-Rao distance. This is detailed in the paragraph below.
**Compatibility between the Fisher-Rao metric and the semi-distance** For $x \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$ and $R > 0$, we denote

$$B_d(x, R) := \{x' \in \mathbb{R}^d \times [u_{\min}, +\infty)^d : d(x, x') \le R\}$$

and

$$\dot{B}_d(x, R) := \{x' \in \mathbb{R}^d \times [u_{\min}, +\infty)^d : d(x, x') < R\}.$$

For $A \subset \mathbb{R}^d \times [u_{\min}, +\infty)^d$ and $x_0 \in A$, we define the set containing all points from geodesics connecting $x_0$ with points of A,

$$\mathcal{G}_{x_0}(A) := \{\gamma(y) : \gamma \text{ is a Fisher-Rao geodesic}, y \in [0, 1], \gamma(0) = x_0, \gamma(1) \in A\}.$$

**Lemma 5.2** (Fisher-Rao geodesics remain within balls w.r.t. the semi-distance). Let $r > 0$ and $x_0 \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$. Then

$$\mathcal{G}_{x_0}(B_d(x_0, r)) \subset B_d(x_0, r).$$

The proof is detailed in Appendix H.3.1. This lemma allows us to use the controls we have on the balls $B_d(x_j^0, r)$ on the paths of the geodesics between $x_j^0$ and a point in $B_d(x_j^0, r)$. It entails that $d$ and $\mathfrak{d}_g$ are in some sense compatible.

To control the certificates, we also require that certain balls (with respect to the semi-distance) around the particles be disjoint. We restrict the possible values for the variances to establish this property: we need upper and lower bounds to establish a “pseudo-quasi triangle inequality” for $d$ (see Appendix H.3.2).

**Lemma 5.3** (Separation of the balls w.r.t. the semi-distance). Let $r, \Delta > 0$. Assume that $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, u_{\max}]^d$. If

$$\min_{i \neq j} d(x_i^0, x_j^0) \ge \max \left\{ \frac{\sqrt{u_{\max}^2 + \frac{1}{4}r^2(2u_{\max}^2 + \tau^2)}}{u_{\min}} (\Delta + r), 2\frac{u_{\max}}{u_{\min}}\Delta \right\} + \sqrt{d \ln \left( \frac{u_{\max}^2}{u_{\min}^2} \right)} =: \Delta_{\tau},$$

then the balls $\dot{B}_d(x_j^0, \Delta) \cap \mathcal{X}$ are disjoint, and for all $j \neq i$, $\mathcal{G}_{x_j^0}(B_d(x_j^0, r) \cap \mathcal{X})$ does not intersect $\dot{B}_d(x_i^0, \Delta)$.

The proof is given in Appendix H.3.2. We stress that this lemma leads to the separation condition in Theorem 5.1. It is crucial to manage the kernel on the far and near regions according to the local positive curvature assumption.

**Non-degeneration of certificates under a minimal separation** Using the compatibility between the semi-distance and the Fisher-Rao metric, and controls on $K_{\text{norm}}$, we can show the existence of non-degenerate certificates under some condition on the minimal separation between the particles of $\mu^0$. This is what the following theorem entails.

**Theorem 5.2.** Assume that $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, u_{\max}]^d$ and that $K_{\text{norm}}$ satisfies the LPC (Definition 5.2) with parameters $s, \Delta, r, \bar{\epsilon}_0, \bar{\epsilon}_2$.

If $\{x_j^0\}_{j=1}^s \subset \mathcal{X}$ satisfies $\min_{i \neq j} d(x_i^0, x_j^0) \ge \Delta_{\tau}$ (defined in Lemma 5.3), then there exists a $(\epsilon_0 = \frac{7}{8}\bar{\epsilon}_0, \epsilon_2 = \frac{15}{32}\bar{\epsilon}_2, r)$-global non-degenerate certificate $\eta$ of the form (18).

Under the same assumptions, for all $j = 1, \dots, s$ there exists a certificate $\eta_j$ for the $j$th near region of the form (18), of parameters $(\bar{\epsilon}_0 = \frac{7}{8}\bar{\epsilon}_0, \bar{\epsilon}_2 = \frac{B_{02}+\bar{\epsilon}_2/16}{2}, r)$.

$$\text{Moreover, } \eta|_{\mathcal{X}} = \Psi^*p \text{ and } \eta_j|_{\mathcal{X}} = \Psi^*p_j \text{ where } \|p\|_{\mathbb{L}} \le \sqrt{2s} \text{ and } \|p_j\|_{\mathbb{L}} \le \sqrt{2}.$$

The proof is based on an adaptation of [Poon et al., 2023, Theorem 2], and is given in Appendix I.

It remains to prove that $K_{\text{norm}}$ satisfies the LPC—this is the purpose of the following theorem, whose proof is provided in Appendix J. We give general bounds for $K_{\text{norm}}$ in dimension $d \ge 1$, and we provide a tighter constant for $\Delta(s)$ in the case $d = 1$.

**Theorem 5.3.** Let $s \ge 2$, $d \ge 1$. Assume that $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, u_{\max}]^d$ and that $\tau \le u_{\min}$. Then $K_{\text{norm}}$ satisfies the LPC with parameters $s, r = \frac{0.3025}{\sqrt{d}}, \bar{\epsilon}_2(r) = 0.13139, \bar{\epsilon}_0(r) = \frac{0.0894}{2d}$, and

$$\Delta(s) = 2\sqrt{11.9 + 3\ln(d + 6.62) + \ln(s - 1)}.$$

Moreover, we can take $B_{02} = \sqrt{4d^2 + 10d}$ (item 1 of Definition 5.2).

**Remark 5.3.** When $d = 1$, the minimal separation condition can be improved: $K_{\text{norm}}$ satisfies the LPC with the same parameters $r, \bar{\epsilon}_0, \bar{\epsilon}_2$ as in Theorem 5.3, with $\Delta(s) = 2\sqrt{13.88 + \ln(s - 1)}$ (Appendix J).
**Remark 5.4.** Note that the assumption on the separation between particles given in Theorem 5.2 depends on $\tau$ through the semi-distance. We write $d = d_\tau$ to emphasize this dependency. At first sight, the separation condition thus appears to depend on $n$ when choosing $\tau^2 = \frac{2u_{\min}^2}{\ln n}$. However, as $\tau \in \mathbb{R}^+ \mapsto d_\tau(x, x')$ is decreasing and $\tau \in \mathbb{R}^+ \mapsto \Delta_\tau$ is increasing, if $\{x_j^0\}_{j=1}^s$ satisfies $\min_{i \neq j} d_{\tau_1}(x_i^0, x_j^0) \ge \Delta_{\tau_1}$ for some $\tau_1 > 0$, then it also satisfies $\min_{i \neq j} d_{\tau_2}(x_i^0, x_j^0) \ge \Delta_{\tau_2}$ for $0 < \tau_2 \le \tau_1$. This motivates the introduction of the following assumption:

$$ \min_{x_i^0 \neq x_j^0} d_{\tau_{\max}}(x_i^0, x_j^0) \ge \Delta_{\tau_{\max}} \quad \text{and} \quad \forall j = 1, \dots, s, \ x_j^0 \in \mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, u_{\max}]^d, \ \tau_{\max} \le u_{\min}. \quad (22) $$

The recovery guarantees are not expressed using this assumption, but it is useful when thinking about the convergence of the solution towards $\mu_\omega^0$ as $n \to \infty$. The target is fixed, but for a good prediction $\tau$ must decrease with the number of observations. The constraint (22) provides in this context a condition on $\mu^0$ that does not depend on $n$.

**Controls on the effective near regions** Theorem 5.3 sets the size of near regions we consider for the certificates. We can lower-bound the Fisher-Rao distance with the semi-distance on these regions, allowing us to extend the control of the estimator to near regions of smaller size (Proposition 3.1). This result is key to provide recovery guarantees for the proximity of $\mu^0$ and the renormalized estimator $\frac{\hat{u}_{n,\omega}}{W}$ (c.f. Corollary 3.1). Its proof can be found in Appendix H.4.

**Lemma 5.4.** Assume that $x, x_0 \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$. If $r_e \le d(x, x_0) \le r = \frac{0.3025}{\sqrt{d}}$, then $\mathfrak{D}_g(x, x_0)^2 \ge \frac{r_e^2}{\tilde{\varepsilon}_3}$ with $\tilde{\varepsilon}_3 = 2.84$.

# 6 Sparsity of the solution under large sample sizes

The certificates constructed in Section 5.2 satisfy the requirements of Assumption 1. Handling such certificates allows us to obtain the non-asymptotic recovery guarantees presented in Section 3. These results are stated in terms of control of the mass in the far and near regions. Nevertheless, they say nothing about the sparsity of the estimator.

In this section, we present results of a different nature. Our estimator is here an exact solution $\mu_{n,\omega}^*$, verifying (7). Following [Duval and Peyré, 2015], we show that for large sample sizes and with high probability, $\mu_{n,\omega}^*$ is sparse and has exactly 1 particle in each near region. We stress that in the following, $0 < \tau \le u_{\min}$ is fixed (does not decrease as $n \to \infty$).

**Non Degenerate Source Condition** We present here the key tool of our analysis, based again on non-degenerate dual certificates (Definition 3.2), with additional properties. We show the existence of a global non-degenerate certificate $\eta_{NDSC}$ verifying $\eta_{NDSC} > -1$. It is said to satisfy the Non Degenerate Source Condition (NDSC).

**Lemma 6.1** (Non-degeneration of $\eta_{NDSC}$). Let $s \ge 2$ and assume that $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, u_{\max}]^d$. There exist $r_{NDSC} > 0$ depending only on $d$, and $\Delta_{NDSC} > 0$ depending only on $u_{\min}, u_{\max}, d, s, \tau$ such that if $\{x_j^0\}_{j=1}^s$ satisfies $\min_{x_i^0 \neq x_j^0} d(x_i^0, x_j^0) \ge \Delta_{NDSC}$, then there exists a $(\varepsilon_{0,NDSC}, \tilde{\varepsilon}_{2,NDSC}, r_{NDSC})$-non-degenerate certificate of the form (18) (of parameters only depending on $d$), that we denote by $\eta_{NDSC}$, verifying

$$ |\eta_{NDSC}| \le 1, \quad |\eta_{NDSC}(x)| = 1 \iff x \in \{x_j^0\}_{j=1}^s, \quad \nabla^2 \eta_{NDSC}(x_j^0) < 0 \quad \forall j = 1, \dots, s. \quad (\text{NDSC}) $$

It is said to satisfy the Non Degenerate Source Condition (NDSC).

Under the same assumptions, for all $j = 1, \dots, s$ there exists a certificate $\eta_{j,NDSC}$ for the $j$th near region of the form (18), of parameters $(\tilde{\varepsilon}_{0,NDSC}, \tilde{\varepsilon}_{2,NDSC}, r_{NDSC})$ only depending on $d$.

$$ \text{Moreover, } \eta_{NDSC}|_{\mathcal{X}} = \Psi^* p_{NDSC} \text{ and } \eta_{j,NDSC}|_{\mathcal{X}} = \Psi^* p_{j,NDSC} \text{ where } \|p_{NDSC}\|_{\mathbb{L}} \le \sqrt{2s} \text{ and } \|p_{j,NDSC}\|_{\mathbb{L}} \le \sqrt{2}. $$

The proof is in Appendix K.1. We construct $\eta_{NDSC}$ in the same way as in Theorem 5.2. The main difference with Theorem 5.1 is the condition $\eta_{NDSC} > -1$, which is verified under a stronger separation condition on $\{x_j^0\}_{j=1}^s$.

The certificate $\eta_{NDSC}$ from Lemma 6.1 is—as the global dual certificate constructed in Theorem 5.1—the vanishing derivative pre-certificate of [Duval and Peyré, 2015, Section 4]. We can then show (see Lemma K.2 in Appendix) that $\eta_{NDSC}$ corresponds to the minimal norm certificate [Duval and Peyré, 2015, Proposition 7]: we have $\eta_{NDSC}|_{\mathcal{X}} = \Psi^* p_{0,0}$ where

$$ p_{0,0} := \arg \min_{p \in \mathbb{L}} \{ \|p\|_{\mathbb{L}} : \Psi^* p \in \partial \|\mu^0\|_{\text{TV}} \}, \quad (p_{0,0}) $$

with $\partial \|\cdot\|_{\text{TV}}$ denoting the subdifferential of the TV-norm. This certificate allows us to control the sparsity index of $\mu_{n,\omega}^*$ for large sample sizes, leading to the following theorem.
**Theorem 6.1.** Under the assumptions of Lemma 6.1, there exists $\kappa_0 > 0$ and $\gamma_0 > 0$ (depending on $\mathcal{X}$, $\tau$ and $\mu^0$) such that for all $\kappa \le \kappa_0$ and if $\|\Gamma_n\|_{\mathbb{L}} \le \gamma_0\kappa$, then $\mu_{n,\omega}^*$ is s-sparse and has exactly 1 particle in each $\mathcal{X}_j^{\text{near}}(r_{\text{NDSC}})$.

The proof can be found in Appendix K.2.

This result should be seen in conjunction with Lemma 3.1. The latter indicates that the assumption $\|\Gamma_n\|_{\mathbb{L}} \le \gamma_0\kappa_0$ is verified with high probability for $n$ large enough. We detail this in the following corollary.

**Corollary 6.1.** Let $c_\kappa > 0$. Assume that the conditions given in Lemma 6.1 hold. Let $n \ge \frac{c_\kappa^2}{\kappa_0^2(2\pi)^{d/2}\tau^{d/2}}$ and $\kappa = \frac{c_\kappa}{(2\pi)^{d/4}\tau^{d/2}\sqrt{n}}$. With probability greater than $1 - C_\Gamma e^{-\left(\frac{\gamma_0\kappa_0}{c_\Gamma}\right)^2}$,

$$ \mu_{n,\omega}^* = \sum_{j=1}^{s} \omega_j^* \delta_{x_j^*} \quad \text{with} \quad \omega_j^* > 0, \quad x_j^* \in \mathcal{X}_j^{\text{near}}(r_{\text{NDSC}}) \quad \forall j = 1, \dots, s. $$

Moreover,

$$ |\omega_j^0 - \omega_j^*| \le c_0 \frac{c_\kappa}{\sqrt{n}} \qquad (23) $$

where $c_0 > 0$ depends on $\mathcal{X}$, $\tau$ and $\mu^0$. In addition, if $c_\kappa$ is chosen as $c_{\kappa,n} = o(\sqrt{n})$, there exists $n_0 \in \mathbb{N}$ depending on $(c_{\kappa,n})_n, \mathcal{X}, \tau, \mu^0$ such that for all $n \ge n_0$,

$$ d(x_j^*, x_j^0)^2 \le \tilde{c}_0 \frac{c_{\kappa,n}}{\sqrt{n}} $$

with probability at least $1 - C_\Gamma e^{-\left(\frac{\gamma_0\kappa_0}{c_\Gamma}\right)^2}$, where $\tilde{c}_0$ depends on $\mathcal{X}$, $\tau$ and $\mu^0$.

The proof is given in Appendix K.3. The result is established for a generic value of $c_\kappa$, but different specific choices can be considered. For instance, choosing $c_{\kappa,n} = \alpha\sqrt{\ln n}$ with $\alpha > 0$, there exists $n_{0,\alpha} \in \mathbb{N}$ depending

on $\mathcal{X}, \tau, \mu^0, \alpha$ such that for all $n \ge n_{0,\alpha}$ and with probability at least $1 - C_\Gamma n^{-\frac{\gamma_0^2\alpha^2}{c_\Gamma^2}}$, $\mu_{n,\omega}^* = \sum_{j=1}^s \omega_j^* \delta_{x_j^*}$ where for all $j = 1, \dots, s$,

$$ |\omega_j^0 - \omega_j^*| \le c_0 \alpha \frac{\sqrt{\ln n}}{\sqrt{n}} \quad \text{and} \quad d(x_j^*, x_j^0)^2 \le \tilde{c}_0 \alpha \frac{\sqrt{\ln n}}{\sqrt{n}}. \qquad (24) $$

Compared to Theorem 3.1, Inequality (24) yields more classical results in parametric estimation. We indeed obtain parametric rates of convergence (up to a logarithmic factor) for the estimation of both weight and location-scale parameters. Moreover, the estimator has a sparsity index exactly matching that of the target measure $\mu^0$. These bounds hold under slightly more restrictive conditions and for sufficiently large sample sizes. We also stress that these results are obtained under a specific tuning of the parameters $\kappa$ and $\tau$, which differs from that used in the previous sections.

# 7 Discussion and open problems

In this paper, we have described some theoretical properties of the BLASSO in the framework of Gaussian mixtures. In particular, we have considered the case where each underlying component has an unknown diagonal covariance. At this stage, further improvements and extensions are still possible and are discussed below.

**Algorithmic considerations** Recall that Theorems 3.1 and 4.1 apply with any $\hat{\mu}_{n,\omega}$ that verifies $J_W(\hat{\mu}_{n,\omega}) \le J_W(\mu^0)$: our estimator is not necessarily the solution of $(P_\kappa)$. So we can restrict $(P_\kappa)$ to any subset of $M(\mathcal{X})^+$ containing $\mu^0$.

This remark is interesting from an algorithmic point of view. We do not know how to solve $(P_\kappa)$ numerically on the entire measure space $M(\mathcal{X})^+$, due to the absence of a parametrization of this space. Hardy, 2023 studied a version of the BLASSO over the space of $K$-sparse measures (discrete measures with less than $K$ particles). This problem, although not convex, is closer to a realistic algorithmic framework, as we can parametrize the space of $K$-sparse measures. If $K \ge s$, $\mu^0$ belongs to this space.

Several algorithms have been proposed to solve the BLASSO on sparse measures, such as the sliding Frank-Wolfe algorithm [Denoyelle et al., 2019] or the Conic Particle Gradient Descent (CPGD) method [Chizat, 2022]. The latter provides a fairly general algorithmic framework together with convergence rate guarantees. CPGD has been previously studied in the context of location mixture models in De Castro et al., 2021a and De Castro et al., 2023. It is therefore natural to explore its applicability to the more challenging location-scale
framework considered in this work. This investigation is part of our ongoing research. We present here a preliminary illustration from this line of work, relying on so-called “natural gradient descent” combined with “conic retraction”. Figure 3 illustrates the convergence behavior of the algorithm when targeting a mixture of two bivariate Gaussian components (see Figure 2). We have access to a sample of size $n = 400$, and start from an initialization with three particles. The method successfully identifies the correct number of components, with one particle effectively disappearing as its associated weight converges to zero.

<figure><img src="image_4.png" /><figcaption>Figure 2: We have access to $n = 400$ samples of the target mixture. The 2 target components have parameters $(a_1^0, a_2^0) = (\frac{1}{2}, \frac{1}{2})$, $(t_1^0, t_2^0) = ((-1, 0), (1, 0))$, $(u_1^0, u_2^0) = ((1, 2), (3, 1))$. The corresponding separation with $\tau = 0.8$ is $d(x_1^0, x_2^0) \approx 1.1$ while $\Delta$ from Theorem 5.3 is around 8.6.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: We take $\tau = 0.8$, $\kappa = 0.01$ and initialize CPGD with 3 particles. From left to right: the evolution of the loss during training; the evolution of the weights $\hat{a}$; of the parameters $\hat{t}$; and $\hat{u}$. The target is represented by green full stars, the initialization by empty black stars, the last estimate by empty red stars. The other iterates are represented by small dots, of different color for each particle.</figcaption></figure>

**Non-diagonal covariance matrices** In this contribution, we are only dealing with the case where the covariances of the mixture are diagonal. Note that we can easily extend our model and results to the case where all covariances share the same (known) orientation, i.e. we can diagonalize them in the same basis. However, we are not yet able to process general covariances with varying or unknown orientations. Indeed, the Fisher-Rao metric seems impractical to work with, mainly because it might be tricky to prove the compatibility between the associated distance and the semi-distance, and obtain the desired controls on the kernel. A possible outcome could be to consider an alternative metric.

**Minimal separation condition** Our results are established under a separation condition for the mixture components (Theorem 5.1). In [De Castro et al., 2021a], the analysis is carried out for location Gaussian mixture models with a common, fixed covariance matrix across components. In this pure location setting, the minimal separation between components is allowed to shrink to zero as the sample size $n$ increases. This favorable behavior stems from the translation-invariant structure of the kernel arising in the pure location setting. Such kernels have been investigated in [De Castro et al., 2025]: this structure enables the construction of certificates based on a “pivot” kernel, whose decay can increase when dealing with closer particles.

In contrast, we consider a substantially more complex location-scale model, where component-specific covariances are unknown and allowed to vary across components. In this setting, the associated kernel loses
translation invariance and does not admit a practical pivot construction. The decay of $K_{\text{norm}}$ does not allow us to consider $\Delta \to 0$ (for $(t, u) \neq (t', u') \in \mathbb{R} \times [u_{\min}, +\infty)$, $e^{-\frac{(t-t')^2}{2(u^2+u'^2+\tau^2)}}$ cannot be as close to 0 as wanted by changing $\tau$).

**Technical side notes** We have used the semi-distance $d$ to define the near and far regions, and the Fisher-Rao metric to control the certificates with Taylor expansions. These are somewhat arbitrary or improvable choices. The semi-distance $d$ seems appropriate to control the kernel, because it measures precisely how 2 points are spaced for the kernel. But the downside is that it does not satisfy the triangle inequality. As a consequence, we need to take $\Delta_{\tau}$ much larger than $\Delta$. The Fisher-Rao metric allows us to retrieve quite easily global bounds for the kernel (i.e. evaluations of the bounds $B_{ij}$, see item 1 of Definition 5.2), but we could have used the Euclidean metric—although it leads to controls that depend on bounds on the variance.
# Summary of the main results

The target measure is $\mu^0 = \sum_{j=1}^s a_j^0 \delta_{(t_j^0, u_j^0)} \in \mathcal{M}(\mathcal{X})$. We assume in the following that $\mathcal{X}$ is a compact set of $\mathbb{R}^d \times [u_{\min}, u_{\max}]^d$, and that $n \ge 2$, $s \ge 2$, $0 < \tau \le u_{\min}$. The results displayed in the following table require more assumptions, that we explicit here. Keep in mind that $d$ depends on $\tau$.

**Assumption 2.**

$$ \min_{i \neq j} d(x_i^0, x_j^0) \geq \max \left\{ \frac{\sqrt{u_{\max}^2 + \frac{1}{4} \frac{0.3025^2}{d}(2u_{\max}^2 + \tau^2)}}{u_{\min}} \left( \Delta + \frac{0.3025}{\sqrt{d}} \right), 2 \frac{u_{\max}}{u_{\min}} \Delta \right\} + \sqrt{d \ln \left( \frac{u_{\max}^2}{u_{\min}^2} \right)} $$

$$ \text{where } \Delta = 2\sqrt{11.9 + 3\ln(d + 6.62) + \ln(s - 1)}. $$

**Assumption 3.**

$$ \min_{i \neq j} d(x_i^0, x_j^0) \geq \Delta_{NDSC} $$

(see Lemma 6.1—$\Delta_{NDSC}$ depends on $s, \tau, d, u_{\min}, u_{\max}$).

Some constants are omitted in the bounds, this is expressed using $\lesssim_c$ when the dependence on $c$ is not taken into account.

<table><caption>Table 1: Recovery guarantees, $\kappa$ not depending on $s$</caption><thead><tr><th></th><th>Parameters</th><th>Asm.</th><th>Bounds</th></tr></thead><tbody><tr><td>Estimation<br>(Prop. 3.1,<br>Thm. 5.1,<br>Lem. 5.4,<br>Cor. 3.1)</td><td>κ = <sup>√2</sup>⁄(2πτ<sup>2</sup>)<sup>d/4</sup>⁄√n</td><td>2</td><td>For 0 &lt; r<sub>e</sub> &le; <sup>0.3025</sup>⁄√d,<br>E [||μ<sub>n,ω</sub>(X<sub>j</sub><sup>near</sup>(r<sub>e</sub>)) - μ<sub>ω</sub><sup>0</sup>(X<sub>j</sub><sup>near</sup>(r<sub>e</sub>))||] &lesssim_d <sup>s</sup>⁄(r<sub>e</sub><sup>2</sup>√nτ<sup>d/2</sup>)<br>For n<sup>-1/6</sup> &le; <sup>0.3025</sup>⁄√d,<br>E [||a<sub>j</sub><sup>0</sup> - <sup>μ<sub>n,ω</sub></sup>⁄W (X<sub>j</sub><sup>near</sup>(n<sup>-1/6</sup>))||] &lesssim_d (<sup>s</sup>⁄(τ<sup>d/2</sup>W(x<sub>j</sub><sup>0</sup>)) + a<sub>j</sub><sup>0</sup>) n<sup>-1/6</sup></td></tr><tr><td>Sparsity<br>(Cor. 6.1)</td><td>c<sub>κ,n</sub> > 0,<br>c<sub>κ,n</sub> = o(√n)<br>κ = <sup>c<sub>κ,n</sub></sup>⁄(2πτ<sup>2</sup>)<sup>d/4</sup>⁄√n</td><td>3</td><td>For n &ge; n<sub>0</sub> (depends on μ<sup>0</sup>, X, τ, (c<sub>κ,n</sub>)<sub>n</sub>),<br>with prob. &ge; 1 - C<sub>Γ</sub>e<sup>-<sup>(γ<sub>0</sub>c<sub>κ,n</sub>)<sup>2</sup></sup>⁄C<sub>Γ</sub></sup> :<br>μ<sub>n,ω</sub><sup>*</sup> = &sum;<sub>j=1</sub><sup>s</sup> ω<sub>j</sub><sup>*</sup> δ<sub>x<sub>j</sub><sup>*</sup></sub> and<br>d(x<sub>j</sub><sup>*</sup>, x<sub>j</sub><sup>0</sup>)<sup>2</sup> &lesssim_μ<sup>0</sup>, X, τ <sup>c<sub>κ,n</sub></sup>⁄√n, |ω<sub>j</sub><sup>0</sup> - ω<sub>j</sub><sup>*</sup>| &lesssim_μ<sup>0</sup>, X, τ <sup>c<sub>κ,n</sub></sup>⁄√n</td></tr><tr><td>Prediction<br>with small<br>regulariza-<br>tion (Prop. 4.1)</td><td>τ = <sup>√2u<sub>min</sub></sup>⁄√(ln n)<br>κ = <sup>4(ln n)<sup>d/2</sup></sup>⁄(4πu<sub>min</sub><sup>2</sup>)<sup>d/2</sup>⁄n</td><td>none</td><td>E [||Φ <sup>μ<sub>n,ω</sub></sup>⁄W - Φμ<sup>0</sup>||<sup>2</sup>]<sub>L<sup>2</sup></sub> &lesssim_d u<sub>min</sub>, u<sub>max</sub> <sup>(ln n)<sup>d/2</sup></sup>⁄n</td></tr><tr><td>Prediction<br>with a good<br>estimator<br>(Thm. 4.1)</td><td>τ = <sup>√2u<sub>min</sub></sup>⁄√(ln n)<br>κ = <sup>√2(ln n)<sup>d/2</sup></sup>⁄(4πu<sub>min</sub><sup>2</sup>)<sup>d/2</sup>⁄√n</td><td>2</td><td>E [||Φ <sup>μ<sub>n,ω</sub></sup>⁄W - Φμ<sup>0</sup>||<sup>2</sup>]<sub>L<sup>2</sup></sub> &lesssim_d u<sub>min</sub>, u<sub>max</sub> <sup>s(ln n)<sup>d/2</sup></sup>⁄n</td></tr></tbody></table>
<table><caption>Table 2: Recovery guarantees, $κ$ depending on $s$</caption><thead><tr><th>Parameters</th><th>Asm.</th><th>Bounds</th></tr></thead><tbody><tr><td>Estimation, v2 (Rk. 3.6, Thm. 5.1, Lem. 5.4)</td><td>κ = <sup>√2</sup>&frasl;<sub>(2πτ<sup>2</sup>)<sup>d/4</sup>√sn</sub></td><td>For 0 &lt; r<sub>e</sub> ≤ <sup>0.3025</sup>&frasl;<sub>√d</sub>,<br>E [||μ<sub>n,ω</sub>(χ<sub>j</sub><sup>near</sup>(r<sub>e</sub>)) - μ<sub>ω</sub><sup>0</sup>(χ<sub>j</sub><sup>near</sup>(r<sub>e</sub>))||] &lesssim; d <sup>√s</sup>&frasl;<sub>r<sub>e</sub><sup>2</sup>√nτ<sup>d/2</sup></sub><br>For n<sup>-1/6</sup> ≤ <sup>0.3025</sup>&frasl;<sub>√d</sub>,<br>E [|a<sub>j</sub><sup>0</sup> - <sup>μ<sub>n,ω</sub></sup>&frasl;<sub>W</sub>(χ<sub>j</sub><sup>near</sup>(n<sup>-1/6</sup>))|] &lesssim; d (<sup>√s</sup>&frasl;<sub>τ<sup>d/2</sup>W(x<sub>j</sub><sup>0</sup>)</sub> + a<sub>j</sub><sup>0</sup>) n<sup>-1/6</sup></td></tr><tr><td>Prediction with a good estimator, v2 (Rk. 4.4)</td><td>τ = <sup>√2u<sub>min</sub></sup>&frasl;<sub>√ln n</sub><br>κ = <sup>√2</sup>&frasl;<sub>(2πτ<sup>2</sup>)<sup>d/4</sup>√sn</sub></td><td>If s = O(n(ln n)<sup>d/2</sup>),<br>E [||Φ<sup>μ<sub>n,ω</sub></sup>&frasl;<sub>W</sub> - Φμ<sup>0</sup>||<sub>L<sup>2</sup></sub><sup>2</sup>] &lesssim; d, u<sub>min</sub>, u<sub>max</sub> <sup>(ln n)<sup>d/2</sup></sup>&frasl;<sub>n</sub></td></tr></tbody></table>
Notation

<table><caption>Table 3: Table of notations</caption><thead><tr><th></th><th>Global notation</th></tr></thead><tbody><tr><td><i>X</i></td><td>compact of ℝ<sup>d</sup> × [u<sub>min</sub>, +∞)<sup>d</sup> with u<sub>min</sub>, u<sub>max</sub> > 0</td></tr><tr><td><i>M</i>(<i>X</i>)<sup>+</sup></td><td>nonnegative Radon measures on <i>X</i></td></tr><tr><td><i>W</i></td><td>reparametrization function. For <i>x</i> = ((t<sub>1</sub>, ..., t<sub>d</sub>), (u<sub>1</sub>, ..., u<sub>d</sub>)) ∈ ℝ<sup>d</sup> × [u<sub>min</sub>, +∞)<sup>d</sup>, <i>W</i>(<i>x</i>) = Π<sub><i>k</i>=1</sub><sup><i>d</i></sup>(2π)<sup>-1/4</sup>(2u<sub><i>k</i></sub><sup>2</sup> + τ<sup>2</sup>)<sup>-1/4</sup></td></tr><tr><td>μ<sup>0</sup>, μ<sub>ω</sub><sup>0</sup></td><td>resp. the target probability measure Σ<sub><i>j</i>=1</sub><sup><i>s</i></sup> a<sub><i>j</i></sub><sup>0</sup>δ<sub>x<sub><i>j</i></sub></sub> where x<sub><i>j</i></sub><sup>0</sup> = (t<sub><i>j</i></sub><sup>0</sup>, u<sub><i>j</i></sub><sup>0</sup>) ∈ <i>X</i>; its reparametrized version <i>W</i>μ<sup>0</sup> = Σ<sub><i>j</i>=1</sub><sup><i>s</i></sup> ω<sub><i>j</i></sub><sup>0</sup>δ<sub>x<sub><i>j</i></sub></sub></td></tr><tr><td>φ, σ, Φ</td><td>resp. the function <i>z</i> ∈ ℝ ↦ <sup>e<sup>-<i>z</i></sup></sup>/<sub>√2π</sub>; its Fourier transform; the operator μ ↦ (z ∈ ℝ<sup><i>d</i></sup> ↦ ∫<sub>ℝ<sup><i>d</i></sup>×[u<sub>min</sub>,+∞)<sup><i>d</i></sup></sub> <sup>d</sup> ∏<sub><i>k</i>=1</sub><sup><i>d</i></sup> <sup>1</sup>/<sub>u<sub><i>k</i></sub></sub>φ(<sup>z<sub><i>k</i></sub> - t<sub><i>k</i></sub></sup>/<sub>u<sub><i>k</i></sub></sub>) dμ(<i>t</i>, <i>u</i>)</td></tr><tr><td><i>X</i><sub>1</sub>, ..., <i>X</i><sub>n</sub></td><td>i.i.d. observations in ℝ<sup><i>d</i></sup>, drawn from <i>f</i><sup>0</sup> = Φμ<sup>0</sup></td></tr><tr><td><i>E</i></td><td>expected value w.r.t. <i>X</i><sub>1</sub>, ..., <i>X</i><sub>n</sub></td></tr><tr><td><i>hat</i><sub><i>n</i></sub></td><td>empirical density <sup>1</sup>/<sub>n</sub> Σ<sub><i>i</i>=1</sub><sup><i>n</i></sup> δ<sub><i>X</i><sub><i>i</i></sub></sub></td></tr><tr><td>κ</td><td>regularization constant</td></tr><tr><td>μ̂<sub><i>n</i>,ω</sub>, μ̂<sub><i>n</i></sub>, μ<sub>n,ω</sub><sup>*</sup></td><td>resp. a measure of <i>M</i>(<i>X</i>)<sup>+</sup> such that <i>J</i><sub><i>W</i></sub>(μ̂<sub><i>n</i>,ω</sub>) ≤ <i>J</i><sub><i>W</i></sub>(μ<sub>ω</sub><sup>0</sup>) (see (<i>P</i><sub>κ</sub>)); the measure μ̂<sub><i>n</i>,ω</sub>/<i>W</i>; an exact solution of (<i>P</i><sub>κ</sub>)</td></tr><tr><td>Γ<sub><i>n</i></sub>, ρ<sub><i>n</i></sub><sup>2</sup></td><td>resp. the noise <i>L</i> ◦ <i>hat</i><sub><i>n</i></sub> − <i>L</i> ◦ <i>f</i><sup>0</sup>; a bound on <i>E</i> [||Γ<sub><i>n</i></sub>||<sub>ℝ</sub><sup>2</sup>] equal to <sup>4</sup>/<sub>(2π)<sup>d/2</sup>τ<sup><i>d</i></sup><i>n</i></sub></td></tr><tr><td colspan="2"><strong>Kernel, differential geometry</strong></td></tr><tr><td>τ, λ, Λ, <i>L</i></td><td>resp. the smoothing parameter τ; the function <i>z</i> ∈ ℝ<sup><i>d</i></sup> ↦ <sup>e<sup>-<i>||z||</sup>2</sup>/<sub>2πτ<sup>2</sup></sub>; its Fourier transform; the operator <i>f</i> ↦ λ * <i>f</i></td></tr><tr><td><i>L</i></td><td>RKHS associated with λ (scalar product given by (5))</td></tr><tr><td><i>Ψ</i>, <i>K</i><sub>norm</sub></td><td>resp. the feature map μ ∈ <i>M</i>(ℝ<sup><i>d</i></sup> × [u<sub>min</sub>, +∞)<sup><i>d</i></sup>) ↦ <i>L</i> ◦ Φ<sub><i>W</i></sub><sup>μ</sup>; the normalized kernel (<i>x</i>, <i>x</i>′) ∈ ℝ<sup><i>d</i></sup> × [u<sub>min</sub>, +∞)<sup><i>d</i></sup> ↦ ⟨Ψδ<sub><i>x</i></sub>, Ψδ<sub><i>x</i>′</sub>⟩<sub><i>L</i></sub></td></tr><tr><td><i>d</i></td><td>semi-distance (<i>x</i>, <i>x</i>′) ↦ √-2 ln(<i>K</i><sub>norm</sub>(<i>x</i>, <i>x</i>′))) (expression given by (11))</td></tr><tr><td><i>X</i><sub><i>j</i></sub><sup>near</sup>(<i>r</i>), <i>X</i><sup>far</sup>(<i>r</i>)</td><td>resp. the <i>j</i>th near region of radius <i>r</i>, {<i>x</i> ∈ <i>X</i> : <i>d</i>(<i>x</i>, <i>x</i><sub><i>j</i></sub><sup>0</sup>) ≤ <i>r</i>}; the far region <i>X</i> \ Σ<sub><i>j</i>=1</sub><sup><i>s</i></sup> <i>X</i><sub><i>j</i></sub><sup>near</sup>(<i>r</i>)</td></tr><tr><td>η</td><td>dual certificate for μ<sup>0</sup>, of the form Ψ<sup>*</sup><i>p</i> with <i>p</i> ∈ ℝ</td></tr><tr><td><i>D</i><sub>η</sub>(μ̂<sub><i>n</i>,ω</sub>, μ<sub>ω</sub><sup>0</sup>)</td><td>Bregman divergence ||μ̂<sub><i>n</i>,ω</sub>||<sub>TV</sub> − ||μ<sub>ω</sub><sup>0</sup>||<sub>TV</sub> − ∫<sub><i>X</i></sub> η <i>d</i> (μ̂<sub><i>n</i>,ω</sub> − μ<sub>ω</sub><sup>0</sup>)</td></tr><tr><td><i>g</i>, <i>δ</i><sub><i>g</i></sub></td><td>resp. the Fisher-Rao metric (see (20)); the associated distance</td></tr><tr><td>γ, ˜γ</td><td>geodesic for the Fisher-Rao metric (parametrized on [0, 1] and by arc-length respectively)</td></tr></tbody></table>
References

Arcones M. A. and Giné E. (1993). "Limit Theorems for U-Processes". *The Annals of Probability* 21.3, pp. 1494-1542. DOI: 10.1214/aop/1176989128.

Balakrishnan S., Wainwright M. J., and Yu B. (2017). "Statistical guarantees for the EM algorithm: From population to sample-based analysis". *The Annals of Statistics* 45.1, pp. 77-120. DOI: 10.1214/16-AOS1435.

Beardon A. F. (1983). *The geometry of discrete groups*. Springer. DOI: 10.1007/978-1-4612-1146-4.

Belkin M. and Sinha K. (2015). "Polynomial Learning of Distribution Families". *SIAM Journal on Computing* 44.4, pp. 889-911. DOI: 10.1137/13090818X.

Bing X., Kong D., and Li B. (2025). "Convergence and Optimality of the EM Algorithm Under Multi-Component Gaussian Mixture Models". arXiv: 2509.08237 [math.ST].

Bouveyron C., Celeux G., Murphy T. B., and Raftery A. E. (2019). *Model-based clustering and classification for data science: with applications in R*. Cambridge series in statistical and probabilistic mathematics. Cambridge University Press. DOI: 10.1017/9781108644181.

Boyer C., De Castro Y., and Salmon J. (2017). "Adapting to unknown noise level in sparse deconvolution". *Information and Inference: A Journal of the IMA* 6.3, pp. 310-348. DOI: 10.1093/imaiai/iaw024.

Bruni C. and Koch G. (1985). "Identifiability of Continuous Mixtures of Unknown Gaussian Distributions". *The Annals of Probability* 13.4, pp. 1341-1357. DOI: 10.1214/aop/1176992817.

Butucea C., Delmas J.-F., Dutfoy A., and Hardy C. (2024). "Off-the-grid prediction and testing for linear combination of translated features". *Electronic Journal of Statistics* 18.2, pp. 2931-2973. DOI: 10.1214/24-EJS2267.

Candès E. J. and Fernandez-Granda C. (2014). "Towards a Mathematical Theory of Super-resolution". *Communications on Pure and Applied Mathematics* 67. DOI: 10.1002/cpa.21455.

Chen X. and Zhang A. Y. (2024). "Achieving optimal clustering in Gaussian mixture models with anisotropic covariance structures". *Advances in Neural Information Processing Systems* 37, pp. 113698-113741. DOI: 10.5220/079017-3610.

Chizat L. (2022). "Sparse optimization on measures with over-parameterized gradient descent". *Mathematical Programming* 194.1, pp. 487-532. DOI: 10.1007/s10107-021-01636-z.

Christmann A. and Steinwart I. (2008). *Support vector machines*. Springer. DOI: 10.1007/978-0-387-77242-4.

De Castro Y., Gadat S., Marteau C., and Maugis C. (2021a). "SuperMix: Sparse Regularization for Mixtures". *Annals of Statistics* 49.3, pp. 1779-1809. DOI: 10.1214/20-AOS2022.

De Castro Y., Gadat S., and Marteau C. (2023). *FastPart: Over-Parameterized Stochastic Gradient Descent for Sparse optimisation on Measures*. arXiv: 2312.05993 [math.OC].

De Castro Y., Gadat S., Marteau C., and Maugis-Rabusseau C. (2021b). "Supplement to "SuperMix: Sparse regularization for mixtures"". *Annals of Statistics* 49.3. DOI: 10.1214/20-AOS2022SUPP.

De Castro Y. and Gamboa F. (2012). "Exact reconstruction using Beurling minimal extrapolation". *Journal of Mathematical Analysis and Applications* 395.1, pp. 336-354. DOI: 10.1016/j.jmaa.2012.05.011.

De Castro Y., Gribonval R., and Jouvin N. (2025). *Effective regions and kernels in continuous sparse regularization, with application to sketched mixtures*. arXiv: 2507.08444 [math.ST].

Denoyelle Q., Duval V., Peyré G., and Soubies E. (2019). "The sliding Frank-Wolfe algorithm and its application to super-resolution microscopy". *Inverse Problems* 36.1, p. 014001. DOI: 10.1088/1361-6420/ab2a29.

Doss N., Wu Y., Yang P., and Zhou H. H. (2023). "Optimal estimation of high-dimensional Gaussian location mixtures". *The Annals of Statistics* 51.1, pp. 62-95. DOI: 10.1214/22-AOS2207.

Duval V. and Peyré G. (2015). "Exact Support Recovery for Sparse Spikes Deconvolution". *Foundations of Computational Mathematics* 15.5, pp. 1315-1355. DOI: 10.1007/s10208-014-9228-6.

Friedman J., Hastie T., Höfling H., and Tibshirani R. (2007). "Pathwise coordinate optimization". *The Annals of Applied Statistics* 1.2, pp. 302-332. DOI: 10.1214/07-AOS131.

Fruhwirth-Schnatter S., Celeux G., and Robert C. P. (2019). *Handbook of mixture analysis*. Chapman & Hall/CRC Press. DOI: 10.1201/9780429055911.
Giard R. (2025). *giardr/blasso_gmm_sympy: Symbolic computation for the BLASSO applied to GMMs with diagonal covariances*. Zenodo. Version v1.1. DOI: 10.5281/zenodo.16942355.

Guha A., Ho N., and Nguyen X. (2021). "On posterior contraction of parameters and interpretability in Bayesian mixture modeling". *Bernoulli* 27.4, pp. 2159–2188. DOI: 10.3150/20-BEJ1275.

Hardt M. and Price E. (2015). "Tight Bounds for Learning a Mixture of Two Gaussians". *Proceedings of the Forty-Seventh Annual ACM Symposium on Theory of Computing*. STOC '15. New York, NY, USA: Association for Computing Machinery, 753–760. ISBN: 9781450335362. DOI: 10.1145/2746539.2746579.

Hardy C. (2023). "Estimation and testing of mixtures of Hilbert-valued features issued from a continuous dictionary". PhD thesis. École des Ponts ParisTech.

Heinrich P. and Kahn J. (2018). "Strong identifiability and optimal minimax rates for finite mixture estimation". *The Annals of Statistics* 46.6A, pp. 2844–2870. DOI: 10.1214/17-AOS1641.

Ho N. and Nguyen X. (2016). "Convergence rates of parameter estimation for some weakly identifiable finite mixtures". *The Annals of Statistics* 44.6, pp. 2726–2755. DOI: 10.1214/16-AOS1444.

Houdard A., Bouveyron C., and Delon J. (2018). "High-Dimensional Mixture Models for Unsupervised Image Denoising (HDMI)". *SIAM Journal on Imaging Sciences* 11.4, pp. 2815–2846. DOI: 10.1137/17M1135694.

Houdré C. and Reynaud-Bouret P. (2003). "Exponential Inequalities, with Constants, for U-statistics of Order Two". *Stochastic Inequalities and Applications*. Basel: Birkhäuser Basel, pp. 55–69. DOI: 10.1007/978-3-0348-8069-5_5.

Ma Y., Wu Y., and Yang P. (2025). "On the Best Approximation by Finite Gaussian Mixtures". *IEEE Transactions on Information Theory* 71.7, pp. 5469–5492. DOI: 10.1109/TIT.2025.3558841.

McLachlan G. J. and Krishnan T. (1997). *The EM algorithm and extensions*. John Wiley & Sons. DOI: 10.1002/9780470191613.

McLachlan G. J. and Peel D. (2000). *Finite mixture models*. Vol. 299. Probability and Statistics – Applied Probability and Statistics Section. Wiley. DOI: 10.1002/0471721182.

Murphy K. P. (2012). *Machine learning : a probabilistic perspective*. MIT Press. ISBN: 9780262018029.

Ohn I. and Lin L. (2023). "Optimal Bayesian estimation of Gaussian mixtures with growing number of components". *Bernoulli* 29.2, pp. 1195–1218. DOI: 10.3150/22-BEJ1495.

Poon C. (2019). *An introduction to sparse spikes recovery via the BLASSO [Lecture Notes]*. URL: https://cmhsp2.github.io/files/teaching/sparsity/blasso.pdf.

Poon C., Keriven N., and Peyré G. (2023). "The geometry of off-the-grid compressed sensing". *Foundations of Computational Mathematics* 23, pp. 241–327. DOI: 10.1007/s10208-021-09545-5.

Regev O. and Vijayaraghavan A. (2017). "On Learning Mixtures of Well-Separated Gaussians". *2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)*, pp. 85–96. DOI: 10.1109/FOCS.2017.17.

Rudin W. (1991). *Functional Analysis*. International series in pure and applied mathematics. McGraw-Hill. ISBN: 9780070542365.

Stahl S. (1993). *The Poincaré half-plane: A gateway to modern geometry*. Jones & Bartlett Learning. ISBN: 086720298X.

Teicher H. (1961). "Identifiability of Mixtures". *The Annals of Mathematical Statistics* 32.1, pp. 244–248. DOI: 10.1214/aoms/1177705155.

Tibshirani R. (1996). "Regression Shrinkage and Selection via the Lasso". *Journal of the Royal Statistical Society. Series B (Methodological)* 58.1, pp. 267–288. DOI: 10.1111/j.2517-6161.1996.tb02080.x.

Tibshirani R. (2023). *High-Dimensional Regression: Lasso [Lecture Notes]*. URL: https://www.stat.berkeley.edu/~ryantibs/statlearn-s23/lectures/lasso.pdf.

Tsybakov A. B. (2008). *Introduction to Nonparametric Estimation*. 1st. Springer Publishing Company, Incorporated. DOI: 10.1007/b13794.

Wu C. F. J. (1983). "On the Convergence Properties of the EM Algorithm". *The Annals of Statistics* 11.1, pp. 95–103. DOI: 10.1214/aos/1176346060.

Wu Y. and Yang P. (2020). "Optimal estimation of Gaussian mixtures via denoised method of moments". *The Annals of Statistics* 48.4, pp. 1981–2007. DOI: 10.1214/19-AOS1873.
Wu Y. and Zhou H. H. (2021). "Randomly initialized EM algorithm for two-component Gaussian mixture achieves near optimality in $O(\sqrt{n})$ iterations." *Mathematical Statistics & Learning* 4. DOI: 10.4171/MSL/29.

Yakowitz S. J. and Spragins J. D. (1968). "On the Identifiability of Finite Mixtures". *The Annals of Mathematical Statistics* 39.1, pp. 209-214. DOI: 10.1214/aoms/1177698520.
# Gaussian Mixture Model with unknown diagonal covariances via continuous sparse regularization

## Appendix

## A Functional framework

The BLASSO operates on the space of Radon measures. In this section, we provide definitions of the operators on measures used throughout the paper.

Let $A \subset \mathbb{R}^d$. In the following, $\mathcal{L}^\infty(A)$ is the set of bounded functions from $A$ to $\mathbb{R}$; $L^1(A)$ (resp. $L^2(A)$) the set of functions whose absolute value (resp. square) have a finite integral. The notions of convolution and Fourier transform can be extended to measures.

**Definition A.1 (Convolution $g * \mu$).** Let $A \subset \mathbb{R}^p$. The convolution between $g \in \mathcal{L}^\infty(A)$ and $\mu \in \mathcal{M}(A)$ is defined by

$$g * \mu := \int_A g(\cdot - x) d\mu(x).$$

**Definition A.2 (Fourier transform over $L^1(A)$ and $\mathcal{M}(A)$).** Let $A \subset \mathbb{R}^p$. We use the Fourier transform defined for $g \in L^1(A)$ by

$$\mathcal{F}[g](\xi) := \int_A e^{-i\langle\xi,x\rangle} g(x) dx.$$

We also use its extension to $\mathcal{M}(A)$:

$$\mathcal{F}[\mu](\xi) := \int_A e^{-i\langle\xi,x\rangle} d\mu(x) \quad \forall \mu \in \mathcal{M}(A).$$

The standard properties of convolution and Fourier transform apply, such as $\mathcal{F}[g * \mu] = \mathcal{F}[g]\mathcal{F}[\mu]$ for all $g \in L^1(A) \cap \mathcal{L}^\infty(A)$ and all $\mu \in \mathcal{M}(A)$ (see, for example, [Rudin, 1991, Part 2]).

## B Existence of a solution to the BLASSO

**Proposition B.1.** The problem $(\mathcal{P}_\kappa)$ has a solution.

*Proof.* $J_W$ is lower semi-continuous on $\mathcal{M}(\mathcal{X})^+$ for the weak* convergence: The TV norm is lower semi-continuous for the weak* convergence. It then suffices to show that $\mu \mapsto L \circ \Phi \frac{\mu}{W}$ is weak* to weak continuous. Let $\mu_j \stackrel{*}{\to} \mu$ in $\mathcal{M}(\mathcal{X})$. We want to show that for $f \in \mathbb{L}$,

$$\langle L \circ \Phi \frac{\mu_j}{W}, f \rangle_{\mathbb{L}} \to \langle L \circ \Phi \frac{\mu}{W}, f \rangle_{\mathbb{L}}.$$

Using [Christmann and Steinwart, 2008, Lemma 4.29], as

$$K_{\text{norm}} : (x, x') \in (\mathbb{R}^d \times [u_{\text{min}}, +\infty)^d)^2 \mapsto \left\langle L \circ \Phi \frac{\delta_x}{W}, L \circ \Phi \frac{\delta_{x'}}{W} \right\rangle_{\mathbb{L}}$$

is continuous (see (19)), we have

$$(x \mapsto L \circ \Phi \frac{\delta_x}{W}) \in \mathcal{C}(\mathcal{X}, \mathbb{L}).$$

Hence

$$(x \mapsto \langle L \circ \Phi \frac{\delta_x}{W}, f \rangle_{\mathbb{L}}) \in \mathcal{C}(\mathcal{X}).$$

As $\mathcal{M}(\mathcal{X}) = \mathcal{C}(\mathcal{X})^*$, it comes

$$\int \langle L \circ \Phi \frac{\delta_x}{W}, f \rangle_{\mathbb{L}} d\mu_j(x) \to \int \langle L \circ \Phi \frac{\delta_x}{W}, f \rangle_{\mathbb{L}} d\mu(x),$$
which concludes this part of the proof.

**Conclusion:** The result follows noticing that $J_W(\mu) \xrightarrow{\|\mu\|_{TV} \to \infty}$ we can restrict the problem to a closed ball of $M(\mathcal{X})^+$. By the Banach-Alaoglu theorem, this ball is weakly* compact (a ball of $M(\mathcal{X})$ is weakly* compact, and $M(\mathcal{X})^+$ is a weakly* closed subspace of $M(\mathcal{X})$). We deduce the existence of a minimizer of $J_W$ on $M(\mathcal{X})^+$. □

## C Proof of Lemma 3.1

Expected value of $\|\Gamma_n\|_{\mathbb{L}}^2$: Note that, according to the definition of $\Gamma_n$,

$$ \Gamma_n = \frac{1}{n} \sum_{i=1}^{n} (L \circ \delta_{X_i} - L \circ f^0). $$

For all $h \in \mathbb{L}$, for all $i \in \{1, \dots, n\}$, we have $\mathbb{E}[\langle L \circ \delta_{X_i}, h \rangle_{\mathbb{L}}] = \langle L \circ f^0, h \rangle_{\mathbb{L}}$. This entails that $\mathbb{E}[\langle \Gamma_n, h \rangle_{\mathbb{L}}] = 0$. Defining

$$ Z_i := L \circ \delta_{X_i} - L \circ f^0 $$

(where we recall that $f^0 = \Phi\mu^0$), observe that $Z_1, \dots, Z_n$ are i.i.d. Since for all $i \in \{1, \dots, n\}$,

$$ \|Z_i\|_{\mathbb{L}}^2 = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} \Lambda |\mathcal{F}[\delta_{X_i} - f^0]|^2 \le \frac{4 \int_{\mathbb{R}^d} \Lambda}{(2\pi)^d}, $$

we have

$$ \begin{aligned} \|\Gamma_n\|_{\mathbb{L}}^2 &= \frac{1}{n^2} \sum_{i=1}^{n} \|Z_i\|_{\mathbb{L}}^2 + \frac{1}{n^2} \sum_{i \ne j} \langle Z_i, Z_j \rangle_{\mathbb{L}}, \\ &\le \frac{4 \int_{\mathbb{R}^d} \Lambda}{(2\pi)^d n} + \frac{1}{n^2} \sum_{i \ne j} \langle Z_i, Z_j \rangle_{\mathbb{L}}. \end{aligned} $$

We deduce that $\mathbb{E}[\|\Gamma_n\|_{\mathbb{L}}^2] \le \frac{4 \int_{\mathbb{R}^d} \Lambda}{(2\pi)^d n}$ using that $\mathbb{E}[Z_i] = 0$ for all $i \in \{1, \dots, n\}$.

**Control in probability:** The control in probability of $\|\Gamma_n\|_{\mathbb{L}}^2$ comes from [De Castro et al., 2021b, Lemma 3], and stems from results on U-processes (see [Arcones and Giné, 1993, Proposition 2.3]). In particular, we have

$$ \forall \rho > 0, \quad \mathbb{P} \left( \|\Gamma_n\|_{\mathbb{L}}^2 > \rho \frac{C_{\Gamma}^2}{n\tau^d(2\pi)^{d/2}} \right) \le C_{\Gamma} e^{-\rho}, \qquad (25) $$

for some positive constant $C_{\Gamma} > 0$.

Expected value of $\|\Gamma_n\|_{\mathbb{L}}^4$: Using (25), it comes that

$$ \begin{aligned} \mathbb{E}[\|\Gamma_n\|_{\mathbb{L}}^4] &= \int_0^\infty \mathbb{P} (\|\Gamma_n\|_{\mathbb{L}}^4 > x) \, dx, \\ &= \int_0^\infty \mathbb{P} (\|\Gamma_n\|_{\mathbb{L}}^2 > \sqrt{x}) \, dx, \\ &\le \int_0^\infty C_\Gamma e^{-\sqrt{x} \frac{n\tau^d(2\pi)^{d/2}}{C_\Gamma^2}} \, dx, \end{aligned} $$

so

$$ \mathbb{E}[\|\Gamma_n\|_{\mathbb{L}}^4] \le C_{\Gamma} \frac{C_{\Gamma}^4}{n^2 \tau^{2d} (2\pi)^{2d}} \int_0^\infty e^{-\sqrt{x}} \, dx = \frac{2C_{\Gamma}^5}{n^2 \tau^{2d} (2\pi)^{2d}}. $$

This proves the desired result.

## D Proof of Theorem 3.1

The following proof is standard when dealing with the BLASSO procedure. Its main arguments can be found in [De Castro et al., 2021a]; we adapt their proof to our setting. It is based on bounds on the Bregman divergence introduced in (13). With W defined by (6), we denote in the following $\hat{\mu}_n = \frac{\hat{\mu}_{n,\omega}}{W}$. We also recall from (9) that
$\mu^0 = \frac{\mu_{\omega}^0}{W}$.

Upper bound on the Bregman divergence: Since $J_W(\hat{\mu}_{n,\omega}) \le J_W(\mu_{\omega}^0)$, we get

$$ \left\| L \circ \hat{f}_n - L \circ \Phi \hat{\mu}_n \right\|_{\mathbb{L}}^2 + 2\kappa \left\| \hat{\mu}_{n,\omega} \right\|_{\text{TV}} \le \left\| L \circ \hat{f}_n - L \circ \Phi \mu^0 \right\|_{\mathbb{L}}^2 + 2\kappa \left\| \mu_{\omega}^0 \right\|_{\text{TV}} \quad (26) $$

which can be rewritten as

$$ \left\| L \circ \hat{f}_n - L \circ \Phi \hat{\mu}_n \right\|_{\mathbb{L}}^2 + 2\kappa D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0) + 2\kappa \int_{\mathcal{X}} \eta d(\hat{\mu}_{n,\omega} - \mu_{\omega}^0) \le \|\Gamma_n\|_{\mathbb{L}}^2, \quad (27) $$

for any dual certificate $\eta$ satisfying the requirements of Assumption 1. As $\eta = \Psi^*p$, recalling (10) we have

$$ \int_{\mathcal{X}} \eta d(\hat{\mu}_{n,\omega} - \mu_{\omega}^0) = \langle p, L \circ \Phi(\hat{\mu}_n - \mu^0) \rangle_{\mathbb{L}}. $$

So using Cauchy-Schwarz inequality, (27) leads to

$$ \left\| L \hat{f}_n - L \circ \Phi \hat{\mu}_n \right\|_{\mathbb{L}}^2 + 2\kappa D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0) - 2\kappa \|p\|_{\mathbb{L}} \left\| L \circ \Phi \mu^0 - L \circ \Phi \hat{\mu}_n \right\|_{\mathbb{L}} \le \|\Gamma_n\|_{\mathbb{L}}^2. $$

The triangle inequality

$$ \|L \circ \Phi \mu^0 - L \circ \Phi \hat{\mu}_n\|_{\mathbb{L}} \le \|\Gamma_n\|_{\mathbb{L}} + \|L \hat{f}_n - L \circ \Phi \hat{\mu}_n\|_{\mathbb{L}} $$

then leads to

$$ \left( \left\| L \circ \hat{f}_n - L \circ \Phi \hat{\mu}_n \right\|_{\mathbb{L}} - \kappa \|p\|_{\mathbb{L}} \right)^2 + 2\kappa D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0) \le \left( \|\Gamma_n\|_{\mathbb{L}} + \kappa \|p\|_{\mathbb{L}} \right)^2. $$

As the Bregman divergence is positive, we deduce from the previous inequality that

$$ D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0) \le \frac{\|\Gamma_n\|_{\mathbb{L}}^2}{2\kappa} + \frac{\kappa}{2} \|p\|_{\mathbb{L}}^2 + \|\Gamma_n\|_{\mathbb{L}} \|p\|_{\mathbb{L}} \quad \text{and} \quad \left\| L \circ \hat{f}_n - L \circ \Phi \hat{\mu}_n \right\|_{\mathbb{L}} \le \|\Gamma_n\|_{\mathbb{L}} + 2\kappa \|p\|_{\mathbb{L}}. \quad (28) $$

The bound in expected value on the Bregman divergence follows by applying Lemma 3.1 and using $\|p\|_{\mathbb{L}} \le \sqrt{c_p s}$. We have

$$ \mathbb{E} [D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0)] \le \frac{\rho_n^2}{2\kappa} + \kappa \frac{c_p}{2}s + \rho_n \sqrt{c_p s}. \quad (29) $$

Taking $\kappa = \frac{\rho_n}{\sqrt{c_p}}$ gives

$$ \mathbb{E} [D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0)] \le \frac{\sqrt{c_p}}{2} \rho_n (1 + \sqrt{s})^2. \quad (30) $$

Lower bound on the Bregman divergence: We use the controls of the non-degenerate certificate $\eta$ on the near and far regions. Since $\hat{\mu}_{n,\omega}$ is nonnegative, according to Definition 3.2 we have

$$ D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0) = \int (1 - \eta) d\hat{\mu}_{n,\omega} \ge \varepsilon_0 \hat{\mu}_{n,\omega}(\mathcal{X}^{\text{far}}(r)) + \varepsilon_2 \sum_{j=1}^{s} \int_{\mathcal{X}_j^{\text{near}}(r)} \mathfrak{d}_g(x, x_j^0)^2 d\hat{\mu}_{n,\omega}(x). \quad (31) $$

Combining (30) and (31), we deduce the bound for the far region (item 1 of Theorem 3.1). To control the mass of the estimator on the jth near region, we make use of the local non-degenerate certificate $\eta_j$ (Definition 3.3). We have, for all $j = 1, \dots, s$,

$$ \begin{aligned} |\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r))| &= \left| \omega_j^0 - \int \eta_j \, d\hat{\mu}_{n,\omega} + \int \eta_j \, d\hat{\mu}_{n,\omega} - \int_{\mathcal{X}_j^{\text{near}}(r)} \, d\hat{\mu}_{n,\omega} \right|, \\ &\le \left| \int \eta_j \, d(\mu_{\omega}^0 - \hat{\mu}_{n,\omega}) \right| + \int_{\mathcal{X}_j^{\text{near}}(r) \setminus \mathcal{X}_j^{\text{near}}(r)} |\eta_j| \, d\hat{\mu}_{n,\omega} + \int_{\mathcal{X}_j^{\text{near}}(r)} |1 - \eta_j| \, d\hat{\mu}_{n,\omega} \\ &\quad + \int_{\mathcal{X}^{\text{far}}(r)} |\eta_j| \, d\hat{\mu}_{n,\omega}, \\ &\le |\langle p_j, L \circ \Phi \mu^0 - L \circ \Phi \hat{\mu}_n \rangle_{\mathbb{L}}| + \tilde{\varepsilon}_2 \sum_{l=1}^s \int_{\mathcal{X}_l^{\text{near}}(r)} \mathfrak{d}_g(x, x_l^0)^2 \, d\hat{\mu}_{n,\omega}(x) \\ &\quad + (1 - \tilde{\varepsilon}_0) \hat{\mu}_{n,\omega}(\mathcal{X}^{\text{far}}(r)), \\ &\le \|p_j\|_{\mathbb{L}} (2 \|\Gamma_n\|_{\mathbb{L}} + 2\kappa \|p\|_{\mathbb{L}}) + \max \left\{ \frac{1 - \tilde{\varepsilon}_0}{\varepsilon_0}, \frac{\tilde{\varepsilon}_2}{\varepsilon_2} \right\} D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0), \end{aligned} \quad (32) $$
where for the last inequality, we have used (31), and

$$|\langle p_j, L \circ \Phi \mu^0 - L \circ \Phi \hat{\mu}_n \rangle_{\mathbb{L}}| \le \|p_j\|_{\mathbb{L}} \left( \|\Gamma_n\|_{\mathbb{L}} + \|L \circ \hat{f}_n - L \circ \Phi \hat{\mu}_n\|_{\mathbb{L}} \right) \le \|p_j\|_{\mathbb{L}} (2 \|\Gamma_n\|_{\mathbb{L}} + 2\kappa \|p\|_{\mathbb{L}})$$

together with (28). As $\|p_j\|_{\mathbb{L}} \le \sqrt{c_p}$, $\|p\|_{\mathbb{L}} \le \sqrt{c_p s}$ and $\kappa = \frac{\rho_n}{\sqrt{c_p}}$, we finally have

$$|\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r))| \le 2\sqrt{c_p}(\|\Gamma_n\|_{\mathbb{L}} + \rho_n\sqrt{s}) + \max\left\{\frac{1 - \tilde{\varepsilon}_0}{\varepsilon_0}, \frac{\tilde{\varepsilon}_2}{\varepsilon_2}\right\} D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0).$$

The result in expected value follows from $\mathbb{E}[\|\Gamma_n\|_{\mathbb{L}}] \le \rho_n$ (Lemma 3.1 used with Jensen's inequality) and (30).

For the stability of the mass (item 3), we use a similar but simpler reasoning. As $\eta \le 1$ and $\hat{\mu}_{n,\omega}$ is nonnegative, we have

$$\begin{aligned} \|\mu_{\omega}^0\|_{\text{TV}} - \|\hat{\mu}_{n,\omega}\|_{\text{TV}} &= \int \eta (d\mu_{\omega}^0 - d\hat{\mu}_{n,\omega}) + \int (\eta - 1) d\hat{\mu}_{n,\omega}, \\ &\le \int \eta (d\mu_{\omega}^0 - d\hat{\mu}_{n,\omega}), \\ &= \langle p, L \circ \Phi \mu^0 - L \circ \Phi \hat{\mu}_n \rangle_{\mathbb{L}}, \\ &\le \|p\|_{\mathbb{L}} (2 \|\Gamma_n\|_{\mathbb{L}} + 2\kappa \|p\|_{\mathbb{L}}). \end{aligned}$$

Using (26), we also have $\|\hat{\mu}_{n,\omega}\|_{\text{TV}} \le \|\mu_{\omega}^0\|_{\text{TV}} + \frac{1}{2\kappa} \|\Gamma_n\|_{\mathbb{L}}^2$. We can conclude by taking the expectation in these inequalities and using Lemma 3.1.

With an s-dependent choice of regularization: Choosing $\kappa = \frac{\rho_n}{\sqrt{c_p s}}$, (29) gives $\mathbb{E}[D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0)] \le 2\rho_n \sqrt{c_p s}$. It comes that

$$\mathbb{E}[|\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r))|] \le 2\sqrt{c_p}\rho_n(1 + \sqrt{s}) + \max\left\{\frac{1 - \tilde{\varepsilon}_0}{\varepsilon_0}, \frac{\tilde{\varepsilon}_2}{\varepsilon_2}\right\} 2\rho_n \sqrt{c_p s}. \quad (33)$$

# E Basic inequalities

The following lemma gives inequalities useful when dealing with the semi-distance $d$. We will use them in various proofs.

**Lemma E.1 (Basic inequalities).** Let $a, b \in \mathbb{R}_+^*$, $c \ge 1$. Then

$$\frac{a^2 + b^2}{2ab} \le c \iff a \in [b(c - \sqrt{c^2 - 1}), b(c + \sqrt{c^2 - 1})] \quad (34)$$

and

$$\frac{a^2 + b^2}{2ab} \le c \implies \frac{|a^2 - b^2|}{a^2 + b^2} \le \sqrt{c^2 - 1} \quad (35)$$

along with

$$\frac{a^2 + b^2}{2ab} \le c \implies \frac{2a^2}{a^2 + b^2} \le c + \sqrt{c^2 - 1}. \quad (36)$$

*Proof.* Let $a, b \in \mathbb{R}_+^*$ and $c \ge 1$.

*Proof of (34):* We have

$$\frac{a^2 + b^2}{2ab} \le c \iff a^2 + b^2 - 2cab \le 0.$$

As $4c^2b^2 - 4b^2 = 4b^2(c^2 - 1) \ge 0$, the roots of this polynomial in $a$ are

$$\frac{2cb \pm 2b\sqrt{c^2 - 1}}{2} = b(c \pm \sqrt{c^2 - 1})$$

from which we deduce (34).

*Proof of (35):* Using (34) and $2ab \le a^2 + b^2$, we get

$$\begin{aligned} \frac{|a^2 - b^2|}{a^2 + b^2} &\le \frac{|a - b||a + b|}{2ab}, \\ &\le \frac{\min\{a, b\}(1 + c + \sqrt{c^2 - 1}) \max\{a, b\}(1 - (c - \sqrt{c^2 - 1}))}{2ab}, \\ &= \frac{(1 + c + \sqrt{c^2 - 1})(1 + \sqrt{c^2 - 1} - c)}{2}, \\ &= \sqrt{c^2 - 1}. \end{aligned}$$
Proof of (36): Using again (34) with $2ab \le a^2 + b^2$, we have

$$ \frac{2a^2}{a^2 + b^2} \le \frac{2ab(c + \sqrt{c^2 - 1})}{2ab}, \\ \le c + \sqrt{c^2 - 1}. $$

□

# F Proofs related to the control of the estimator on the effective near regions

## F.1 Proof of Proposition 3.1

The proof is similar to that of Theorem 3.1 in Section D. Let $j \in \{1, \dots, s\}$ and $\eta_j$ a corresponding local non-degenerate dual certificate satisfying the requirements of Assumption 1. According to Definition 3.3, for all $x \in \mathcal{X}_j^{\text{near}}(r)$, we have $|\eta_j| \le 1 + \tilde{\varepsilon}_2 \mathfrak{d}_g(x, x_j^0)^2$. We deduce that

$$ \begin{aligned} |\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))| &= \left| \omega_j^0 - \int \eta_j \, d\hat{\mu}_{n,\omega} + \int \eta_j \, d\hat{\mu}_{n,\omega} - \int_{\mathcal{X}_j^{\text{near}}(r_e)} \, d\hat{\mu}_{n,\omega} \right|, \\ &\le \left| \int \eta_j \, d(\mu_\omega^0 - \hat{\mu}_{n,\omega}) \right| + \int_{\mathcal{X}_j^{\text{near}}(r) \setminus \mathcal{X}_j^{\text{near}}(r_e)} |\eta_j| \, d\hat{\mu}_{n,\omega} + \int_{\mathcal{X}_j^{\text{near}}(r) \setminus \mathcal{X}_j^{\text{near}}(r_e)} |\eta_j| \, d\hat{\mu}_{n,\omega} \\ &\quad + \int_{\mathcal{X}_j^{\text{near}}(r_e)} |1 - \eta_j| \, d\hat{\mu}_{n,\omega} + \int_{\mathcal{X}_j^{\text{far}}(r)} |\eta_j| \, d\hat{\mu}_{n,\omega}, \\ &\le \|p_j\|_{\mathcal{L}} (2 \|\Gamma_n\|_{\mathcal{L}} + 2\kappa \|p\|_{\mathcal{L}}) + \tilde{\varepsilon}_2 \sum_{l \neq j} \int_{\mathcal{X}_l^{\text{near}}(r)} \mathfrak{d}_g(x, x_l^0)^2 \, d\hat{\mu}_{n,\omega} + (1 - \tilde{\varepsilon}_0) \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{far}}(r)) \\ &\quad + \tilde{\varepsilon}_2 \int_{\mathcal{X}_j^{\text{near}}(r_e)} \mathfrak{d}_g(x, x_j^0)^2 \, d\hat{\mu}_{n,\omega} + \int_{\mathcal{X}_j^{\text{near}}(r) \setminus \mathcal{X}_j^{\text{near}}(r_e)} (1 + \tilde{\varepsilon}_2 \mathfrak{d}_g(x, x_j^0)^2) \, d\hat{\mu}_{n,\omega}. \end{aligned} $$

From (15) we get $1 \le \frac{\tilde{\varepsilon}_3}{r_e^2} \mathfrak{d}_g(x_j^0, x)^2$ for all $x \in \mathcal{X}_j^{\text{near}}(r) \setminus \mathcal{X}_j^{\text{near}}(r_e)$, so $1 + \tilde{\varepsilon}_2 \mathfrak{d}_g(x_j^0, x)^2 \le (\frac{\tilde{\varepsilon}_3}{r_e^2} + \tilde{\varepsilon}_2) \mathfrak{d}_g(x_j^0, x)^2$. Using again (31), we deduce that

$$ |\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))| \le \|p_j\|_{\mathcal{L}} (2 \|\Gamma_n\|_{\mathcal{L}} + 2\kappa \|p\|_{\mathcal{L}}) + \max \left\{ \frac{1 - \tilde{\varepsilon}_0}{\varepsilon_0}, \frac{1}{\varepsilon_2} \left( \frac{\tilde{\varepsilon}_3}{r_e^2} + \tilde{\varepsilon}_2 \right) \right\} D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0). \quad (37) $$

We can conclude the proof using the controls on $\mathbb{E}[\|\Gamma_n\|_{\mathcal{L}}]$, $\|p_j\|_{\mathcal{L}}$, $\|p\|_{\mathcal{L}}$, $\mathbb{E}[D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0)]$ stemming from our assumptions along with Lemma 3.1 and (30). Choosing $\kappa = \frac{\rho_n}{\sqrt{c_p s}}$,

$$ |\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))| \le 2\sqrt{c_p}\rho_n(1 + \sqrt{s}) + \max \left\{ \frac{1 - \tilde{\varepsilon}_0}{\varepsilon_0}, \frac{1}{\varepsilon_2} \left( \frac{\tilde{\varepsilon}_3}{r_e^2} + \tilde{\varepsilon}_2 \right) \right\} \frac{\sqrt{c_p}}{2} \rho_n(1 + \sqrt{s})^2. $$

With an $s$-dependent choice of regularization: Choosing $\kappa = \frac{\rho_n}{\sqrt{c_p s}}$, using (37) and the control of $\mathbb{E}[D_{\eta}(\hat{\mu}_{n,\omega}, \mu_{\omega}^0)]$ stemming from (29) we get

$$ \begin{aligned} \mathbb{E}[|\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))|] &\le 4\sqrt{c_p}\rho_n + \max \left\{ \frac{1 - \tilde{\varepsilon}_0}{\varepsilon_0}, \frac{1}{\varepsilon_2} \left( \frac{\tilde{\varepsilon}_3}{r_e^2} + \tilde{\varepsilon}_2 \right) \right\} 2\rho_n \sqrt{c_p s}, \\ &\lesssim \frac{\sqrt{s}}{\tau^{d/2}\sqrt{n}r_e^2} \end{aligned} \quad (38) $$

keeping only the dependence on $s, \tau, r_e, n$.

## F.2 Proof of Corollary 3.1

Let $0 < r_e \le r$. We first prove that

$$ \left| \frac{\hat{\mu}_{n,\omega}}{W}(\mathcal{X}_j^{\text{near}}(r_e)) - a_j^0 \right| \le (1 + H(r)r_e)W(x_j^0)^{-1} \left| \omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e)) \right| + a_j^0 H(r)r_e \quad (39) $$

where

$$ H(r) = \frac{(e^{r^2} + \sqrt{e^{2r^2} - 1})^{d/2} - 1}{r}. \quad (40) $$
We use the triangle inequality

$$ \left| a_j^0 - \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W} \right| \le \underbrace{\left| a_j^0 - \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W(x_j^0)} \right|}_{=:A} + \underbrace{\left| \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W(x_j^0)} - \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W} \right|}_{=:B} . $$

**Control of B:** We recall the definitions of effective near regions (14) and of W (see (6)). We have

$$ \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W} = \int_{\mathcal{X}_j^{\text{near}}(r_e)} \frac{W(x_j^0)}{W(x)} d\left(\frac{\hat{\mu}_{n,\omega}}{W(x_j^0)}(x)\right) $$

so

$$ \left| \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W}(\mathcal{X}_j^{\text{near}}(r_e)) - \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W(x_j^0)} \right| \le \max \left\{ \sup_{x \in \mathcal{X}_j^{\text{near}}(r_e)} \frac{W(x_j^0)}{W(x)} - 1, 1 - \inf_{x \in \mathcal{X}_j^{\text{near}}(r_e)} \frac{W(x_j^0)}{W(x)} \right\} \times \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W(x_j^0)} . $$

For $x \in \mathcal{X}_j^{\text{near}}(r_e)$, as $\prod_{k=1}^d \frac{(u_{j,k}^0)^2+u_k^2+\tau^2}{\sqrt{2(u_{j,k}^0)^2+\tau^2}\sqrt{2u_k^2+\tau^2}} \le e^{r_e^2}$ and as each term of this product is greater than 1, we

have $\frac{(u_{j,k}^0)^2+u_k^2+\tau^2}{\sqrt{2(u_{j,k}^0)^2+\tau^2}\sqrt{2u_k^2+\tau^2}} \le e^{r_e^2}$ for all $k=1, \dots, d$. Using (34), this implies

$$ \sqrt{u_k^2 + \frac{\tau^2}{2}} \in \left[ \sqrt{(u_{j,k}^0)^2 + \frac{\tau^2}{2}(e^{r_e^2} - \sqrt{e^{2r_e^2} - 1})}, \sqrt{(u_{j,k}^0)^2 + \frac{\tau^2}{2}(e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})} \right] \quad \forall k = 1, \dots, d, $$

from which we deduce that

$$ \max \left\{ \sup_{x \in \mathcal{X}_j^{\text{near}}(r_e)} \frac{W(x_j^0)}{W(x)} - 1, 1 - \inf_{x \in \mathcal{X}_j^{\text{near}}(r_e)} \frac{W(x_j^0)}{W(x)} \right\} \le \max \left\{ (e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^{d/2} - 1, 1 - (e^{r_e^2} - \sqrt{e^{2r_e^2} - 1})^{d/2} \right\}, \\ = (e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^{d/2} - 1 $$

where we used that $(e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^{d/2} \ge 1$ to establish that

$$ 1 - (e^{r_e^2} - \sqrt{e^{2r_e^2} - 1})^{d/2} = 1 - (e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^{-d/2} \le (e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^{d/2} - 1. $$

Hence

$$ \left| \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W}(\mathcal{X}_j^{\text{near}}(r_e)) - \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W(x_j^0)} \right| \le \left( (e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^{d/2} - 1 \right) \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W(x_j^0)} . \quad (41) $$

**Control of A and proof of (39):** As

$$ \left| a_j^0 - \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W(x_j^0)} \right| = W(x_j^0)^{-1} \left| \omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e)) \right|, $$

from (41) we get

$$ \left| \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e))}{W}(\mathcal{X}_j^{\text{near}}(r_e)) - a_j^0 \right| \le (e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^{d/2} W(x_j^0)^{-1} \left| \omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(r_e)) \right| + a_j^0 \left( (e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^{d/2} - 1 \right). $$

We conclude the proof of (39) by noticing that $h: r_e \in \mathbb{R}^+ \mapsto (e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^{d/2}$ is convex (for all $d \in \mathbb{N}^*$), hence for $r_e \le r$ we have

$$ (e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^{d/2} \le h(0) + \frac{h(r) - h(0)}{r} r_e = 1 + \frac{h(r) - 1}{r} r_e. $$

**Conclusion:** Taking $r_e = n^{-\alpha}$ with $\alpha > 0$, (39) gives

$$ \mathbb{E} \left[ \left| a_j^0 - \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(n^{-\alpha}))}{W} \right| \right] \le W(x_j^0)^{-1} \mathbb{E} \left[ \left| \omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(n^{-\alpha})) \right| \right] (1 + H(r)n^{-\alpha}) + a_j^0 H(r)n^{-\alpha} $$

where $H(r)$ is defined by (40). Proposition 3.1 gives $\mathbb{E} [ |\omega_j^0 - \hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(n^{-\alpha}))| ] \lesssim \frac{s}{\tau^{d/2}\sqrt{n r_e^2}} = \frac{s}{\tau^{d/2}n^{1/2-2\alpha}}$.

We choose $r_e = n^{-1/6}$ to balance the terms.

With $\kappa = \frac{\rho_n}{\sqrt{c_p s}}$: Choosing $\kappa = \frac{\rho_n}{\sqrt{c_p s}}$, using (38) we get for $n^{-1/6} \le r$

$$ \mathbb{E} \left[ \left| a_j^0 - \frac{\hat{\mu}_{n,\omega}(\mathcal{X}_j^{\text{near}}(n^{-1/6}))}{W} \right| \right] \lesssim \left( W(x_j^0)^{-1} \sqrt{s} \tau^{-d/2} + a_j^0 \right) n^{-1/6}. \quad (42) $$
# G Proofs related to guarantees on the prediction

We will use the following lemma to go from controls of $\|L \circ \Phi(\hat{\mu}_n - \mu^0)\|_{\mathbb{L}}^2$ to controls on $\|\Phi(\hat{\mu}_n - \mu^0)\|_{L^2(\mathbb{R}^d)}^2$.

**Lemma G.1 (Control of the high frequencies).** Assume that $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, +\infty)^d$. Let $\tau > 0$. We work with $\Lambda(\xi) = e^{-\frac{1}{2}\tau^2}\|\xi\|_{\mathbb{L}^2}^2$. Let $\mu_1, \mu_2 \in \mathcal{M}(\mathcal{X})$. Then

$$ \|\Phi(\mu_1 - \mu_2)\|_{L^2}^2 \le e \|L \circ \Phi(\mu_1 - \mu_2)\|_{\mathbb{L}}^2 + \frac{2}{(2\pi)^d} (\|\mu_1\|_{\text{TV}}^2 + \|\mu_2\|_{\text{TV}}^2) \frac{\tau^d d^{d/2}}{2^{d/2} u_{\min}^{2d}} e^{-2 \frac{u_{\min}^2}{\tau^2}}. $$

*Proof.* Let $T > 0$. First write

$$ \|\Phi(\mu_1 - \mu_2)\|_{L^2}^2 = \frac{1}{(2\pi)^d} \int_{[-\frac{1}{T}, \frac{1}{T}]^d} |\mathcal{F}[\Phi(\mu_1 - \mu_2)]|^2 + \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d \setminus [-\frac{1}{T}, \frac{1}{T}]^d} |\mathcal{F}[\Phi(\mu_1 - \mu_2)]|^2. $$

Then, remark that

$$ \begin{aligned} \frac{1}{(2\pi)^d} \int_{[-\frac{1}{T}, \frac{1}{T}]^d} |\mathcal{F}[\Phi(\mu_1 - \mu_2)]|^2 &= \frac{1}{(2\pi)^d} \int_{[-\frac{1}{T}, \frac{1}{T}]^d} \frac{\Lambda}{\Lambda} |\mathcal{F}[\Phi(\mu_1 - \mu_2)]|^2, \\ &\le e^{\frac{d\tau^2}{2T^2}} \frac{1}{(2\pi)^d} \int_{[-\frac{1}{T}, \frac{1}{T}]^d} \Lambda |\mathcal{F}[\Phi(\mu_1 - \mu_2)]|^2, \\ &\le e^{\frac{d\tau^2}{2T^2}} \|L \circ \Phi(\mu_1 - \mu_2)\|_{\mathbb{L}}^2. \end{aligned} $$

Concerning the high frequencies of $\Phi(\mu_1 - \mu_2)$, recall that $u_1, \dots, u_d \ge u_{\min}$ for all $((t_1, \dots, t_d), (u_1, \dots, u_d)) \in \mathcal{X}$. Hence

$$ \begin{aligned} \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d \setminus [-\frac{1}{T}, \frac{1}{T}]^d} |\mathcal{F}[\Phi(\mu_1 - \mu_2)]|^2 &\le \frac{2}{(2\pi)^d} (\|\mu_1\|_{\text{TV}}^2 + \|\mu_2\|_{\text{TV}}^2) \int_{\mathbb{R}^d \setminus [-\frac{1}{T}, \frac{1}{T}]^d} e^{-u_{\min}^2 \|\xi\|_{\mathbb{L}^2}^2} d\xi, \\ &\le \frac{2}{(2\pi)^d} (\|\mu_1\|_{\text{TV}}^2 + \|\mu_2\|_{\text{TV}}^2) \left( \frac{T}{u_{\min}^2} e^{-\frac{u_{\min}^2}{T^2}} \right)^d \end{aligned} $$

using

$$ \int_{\mathbb{R}^d \setminus [-\frac{1}{T}, \frac{1}{T}]} e^{-u_{\min}^2 z^2} dz = 2 \frac{1}{T} \int_{[1, +\infty)} e^{-u_{\min}^2 \frac{z^2}{T^2}} dz \le 2 \frac{1}{T} \int_{[1, +\infty)} z e^{-u_{\min}^2 \frac{z^2}{T^2}} dz = \frac{T}{u_{\min}^2} e^{-\frac{u_{\min}^2}{T^2}}. \quad (43) $$

So

$$ \|\Phi(\mu_1 - \mu_2)\|_{L^2(\mathbb{R}^d)}^2 \le e^{\frac{d\tau^2}{2T^2}} \|L \circ \Phi(\mu_1 - \mu_2)\|_{\mathbb{L}}^2 + \frac{2}{(2\pi)^d} (\|\mu_1\|_{\text{TV}}^2 + \|\mu_2\|_{\text{TV}}^2) \left( \frac{T}{u_{\min}^2} e^{-\frac{u_{\min}^2}{T^2}} \right)^d. $$

Taking $T = \frac{\tau\sqrt{d}}{\sqrt{2}}$, we get

$$ \|\Phi(\mu_1 - \mu_2)\|_{L^2}^2 \le e \|L \circ \Phi(\mu_1 - \mu_2)\|_{\mathbb{L}}^2 + \frac{2}{(2\pi)^d} (\|\mu_1\|_{\text{TV}}^2 + \|\mu_2\|_{\text{TV}}^2) \frac{\tau^d d^{d/2}}{2^{d/2} u_{\min}^{2d}} e^{-2 \frac{u_{\min}^2}{\tau^2}}. $$

□

## G.1 Proof of Proposition 4.1

We do not make any assumption on the existence of dual certificates in this proof. From $J_W(\hat{\mu}_{n,\omega}) \le J_W(\mu_{\omega}^0)$, we have

$$ \begin{aligned} \|L \circ \Phi(\hat{\mu}_n - \mu^0)\|_{\mathbb{L}}^2 &\le 2 \|L \circ \hat{f}_n - L \circ \Phi \hat{\mu}_n\|_{\mathbb{L}}^2 + 2 \|L \circ \hat{f}_n - L \circ \Phi \mu^0\|_{\mathbb{L}}^2, \\ &\le 4 \|\Gamma_n\|_{\mathbb{L}}^2 + 4\kappa \|\mu_{\omega}^0\|_{\text{TV}}. \end{aligned} $$

Combining this inequality with Lemma G.1, we get

$$ \|\Phi(\hat{\mu}_n - \mu^0)\|_{L^2}^2 \le e \left( 4 \|\Gamma_n\|_{\mathbb{L}}^2 + 4\kappa \|\mu_{\omega}^0\|_{\text{TV}} \right) + \frac{2}{(2\pi)^d} (\|\mu_1\|_{\text{TV}}^2 + \|\mu_2\|_{\text{TV}}^2) \frac{\tau^d d^{d/2}}{2^{d/2} u_{\min}^{2d}} e^{-2 \frac{u_{\min}^2}{\tau^2}}. $$
Remark that for $\mu \in \mathcal{M}(\mathcal{X})$ we have $\|\frac{\mu}{W}\|_{\text{TV}} \le \|\mu\|_{\text{TV}} \sup_{\mathcal{X}} \frac{1}{W} \le (2\pi)^{d/4}(2u_{\text{max}}^2 + \tau^2)^{d/4} \|\mu\|_{\text{TV}}$, and in the same way $\|\mu\|_{\text{TV}} \le (2\pi)^{-d/4}(2u_{\text{min}}^2 + \tau^2)^{-d/4} \|\frac{\mu}{W}\|_{\text{TV}}$. As $J_W(\hat{\mu}_{n,\omega}) \le J_W(\mu_{\omega}^0)$ implies that $\|\hat{\mu}_{n,\omega}\|_{\text{TV}} \le \frac{1}{2\kappa} \|\Gamma_n\|_{\text{L}}^2 + \|\mu_{\omega}^0\|_{\text{TV}}$, we deduce that

$$ 
\begin{aligned}
\frac{2}{(2\pi)^d} (\|\hat{\mu}_n\|_{\text{TV}}^2 + \|\mu^0\|_{\text{TV}}^2) &\le \frac{2}{(2\pi)^d} \left( (2\pi)^{d/2}(2u_{\text{max}}^2 + \tau^2)^{d/2} \left( \frac{1}{2\kappa^2} \|\Gamma_n\|_{\text{L}}^4 + 2 \|\mu_{\omega}^0\|_{\text{TV}}^2 \right) + \|\mu^0\|_{\text{TV}}^2 \right), \\
&\le \frac{(2\pi)^{-d/2}(2u_{\text{max}}^2 + \tau^2)^{d/2}}{\kappa^2} \|\Gamma_n\|_{\text{L}}^4 + \frac{2}{(2\pi)^d} \left( 2 \left( \frac{2u_{\text{max}}^2 + \tau^2}{2u_{\text{min}}^2 + \tau^2} \right)^{d/2} + 1 \right) \|\mu^0\|_{\text{TV}}^2, \\
&\le \frac{(2\pi)^{-d/2}(2u_{\text{max}}^2 + \tau^2)^{d/2}}{\kappa^2} \|\Gamma_n\|_{\text{L}}^4 + \frac{2}{(2\pi)^d} \left( 2 \left( \frac{u_{\text{max}}}{u_{\text{min}}} \right)^d + 1 \right) \|\mu^0\|_{\text{TV}}^2
\end{aligned} 
\quad (44) $$

using that $\tau \in \mathbb{R}^+ \mapsto \frac{2u_{\text{max}}^2 + \tau^2}{2u_{\text{min}}^2 + \tau^2}$ is decreasing. Hence

$$ 
\begin{aligned}
\|\Phi(\hat{\mu}_n - \mu^0)\|_{L^2}^2 &\le e \left( 4 \|\Gamma_n\|_{\text{L}}^2 + 4\kappa(2\pi)^{-d/4}(2u_{\text{min}}^2 + \tau^2)^{-d/4} \|\mu^0\|_{\text{TV}} \right) \\
&\quad + \left( \frac{(2\pi)^{-d/2}(2u_{\text{max}}^2 + \tau^2)^{d/2}}{\kappa^2} \|\Gamma_n\|_{\text{L}}^4 + \frac{2}{(2\pi)^d} \left( 2 \left( \frac{u_{\text{max}}}{u_{\text{min}}} \right)^d + 1 \right) \|\mu^0\|_{\text{TV}}^2 \right) \frac{\tau^d d^{d/2}}{2^{d/2} u_{\text{min}}^{2d}} e^{-2 \frac{u_{\text{min}}^2}{\tau^2}}.
\end{aligned} $$

With Lemma 3.1, choosing $\tau = \frac{\sqrt{2}u_{\text{min}}}{\sqrt{\ln n}}$ and $\kappa = \rho_n^2$ it comes

$$ 
\begin{aligned}
\mathbb{E} \left[ \|\Phi(\hat{\mu}_n - \mu^0)\|_{L^2}^2 \right] &\le 4e \rho_n^2 \left( 1 + (2\pi)^{-d/4} (\sqrt{2}u_{\text{min}})^{-d/2} \left(1 + \frac{1}{\ln n}\right)^{-d/4} \|\mu^0\|_{\text{TV}} \right) \\
&\quad + \left( \tilde{C}_{\Gamma} \pi^{-d/2} \left( u_{\text{max}}^2 + \frac{u_{\text{min}}^2}{\ln n} \right)^{d/2} + 2(2\pi)^{-d} \|\mu^0\|_{\text{TV}}^2 \left( 2 \left( \frac{u_{\text{max}}}{u_{\text{min}}} \right)^d + 1 \right) \right) \frac{d^{d/2}}{(\ln n)^{d/2} u_{\text{min}}^d n} \\
&\lesssim \frac{(\ln n)^{d/2}}{n}.
\end{aligned} $$

## G.2 Prediction with Kernel Density Estimation

**Lemma G.2.** With $\mathcal{X} \subset \mathbb{R}^d \times [u_{\text{min}}, +\infty)^d$, setting $\tau = \frac{1}{\sqrt{\ln n} \frac{u_{\text{min}}}{4^d}}$, omitting the dependence on $d$ we have

$$ \mathbb{E} \left[ \left\| L \circ \hat{f}_n - \Phi \mu^0 \right\|_{L^2}^2 \right] \lesssim \frac{(\ln n)^{d/2}}{n^{\frac{4}{d+4}} u_{\text{min}}^{d+4}}. $$

*Proof.* Control of $\|L \circ \Phi \mu^0 - \Phi \mu^0\|_{L^2}^2$: As $|\mathcal{F}[\Phi \mu^0](\xi)|^2 \le \|\mu^0\|_{\text{TV}}^2 e^{-u_{\text{min}}^2 \|\xi\|_2^2}$ for $\xi \in \mathbb{R}^d$ and using (43), it comes that for $T > 0$,

$$ 
\begin{aligned}
\|L \circ \Phi \mu^0 - \Phi \mu^0\|_{L^2}^2 &= \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |\Lambda(\xi) - 1|^2 |\mathcal{F}[\Phi \mu^0](\xi)|^2 d\xi, \\
&\le \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d \setminus [-\frac{1}{T}, \frac{1}{T}]} |\mathcal{F}[\Phi \mu^0](\xi)|^2 d\xi + \frac{1}{(2\pi)^d} \frac{d^2 \tau^4}{4T^4} \int_{[-\frac{1}{T}, \frac{1}{T}]} |\mathcal{F}[\Phi \mu^0](\xi)|^2 d\xi, \\
&\le \frac{\|\mu^0\|_{\text{TV}}^2}{(2\pi)^d} \left( \frac{T^d}{u_{\text{min}}^{2d}} e^{-\frac{du_{\text{min}}^2}{T^2}} + \frac{2^d \tau^4 d^2}{4T^4} \right),
\end{aligned} $$

where we used that for $\xi \in \mathbb{R}^d$, $|\Lambda(\xi) - 1|^2 = |e^{-\frac{\tau^2}{2}} \|\xi\|_2^2 - 1| \le 1$ and

$$ |\Lambda(\xi) - 1|^2 \le |e^{-\frac{\tau^2}{2T^2}} - 1| \le \frac{d^2 \tau^4}{4T^4} \quad \forall \xi \in [-\frac{1}{T}, \frac{1}{T}]^d. $$

Control of $\mathbb{E} \left[ \left\| L \circ \hat{f}_n - \Phi \mu^0 \right\|_{L^2}^2 \right]$ : Using Lemma 3.1, we have

$$ \mathbb{E} \left[ \left\| L \circ \hat{f}_n - \Phi \mu^0 \right\|_{L^2}^2 \right] \le 2\rho_n^2 + 2 \left\| L \circ \Phi \mu^0 - \Phi \mu^0 \right\|_{L^2}^2 \le \frac{8}{(2\pi)^{d/2} \tau^d n} + \frac{2 \|\mu^0\|_{\text{TV}}^2}{(2\pi)^d} \left( \frac{T^d}{u_{\text{min}}^{2d}} e^{-\frac{du_{\text{min}}^2}{T^2}} + \frac{2^d \tau^4 d^2}{4T^4} \right). $$
To balance these terms, we choose $T = \frac{u_{\min}\sqrt{d}}{\sqrt{\ln n}}$ and $\tau = \frac{1}{\sqrt{\ln n}n^{\frac{4}{d+4}}}$. It comes

$$ \mathbb{E} \left[ \left\| L \circ \hat{f}_n - \Phi \mu^0 \right\|_{L^2}^2 \right] \le \frac{8(\ln n)^{d/2}}{(2\pi)^{d/2} n^{\frac{4}{d+4}}} + \frac{2 \|\mu^0\|_{\text{TV}}^2}{(2\pi)^d} \left( \frac{d^{d/2}}{u_{\min}^d (\ln n)^{d/2} n} + \frac{2^d (\ln n)^{d/2}}{4d^{d/2} u_{\min}^{d+4} n^{\frac{4}{d+4}}} \right). $$

□

## G.3 Proof of Theorem 4.1

Equation (28) gives

$$ \|L \circ \Phi(\hat{\mu}_n - \mu^0)\|_{\mathbb{L}} \le 2 \|\Gamma_n\|_{\mathbb{L}} + 2\kappa \|p\|_{\mathbb{L}} $$

from which we deduce, using Lemma 3.1 and $\|p\|_{\mathbb{L}} \le \sqrt{c_p s}$, that

$$ \mathbb{E} \left[ \left\| L \circ \Phi(\hat{\mu}_n - \mu^0) \right\|_{\mathbb{L}}^2 \right] \le 4(\rho_n + \kappa \sqrt{c_p s})^2. $$

To go from a control of $\mathbb{E} \left[ \left\| L \circ \Phi(\hat{\mu}_n - \mu^0) \right\|_{\mathbb{L}}^2 \right]$ to a bound on $\mathbb{E} \left[ \left\| \Phi \hat{\mu}_n - \Phi \mu^0 \right\|_{L^2(\mathbb{R}^d)}^2 \right]$, we use (44) with Lemmas G.1 and 3.1. We get

$$ \begin{aligned} \mathbb{E} \left[ \left\| \Phi \hat{\mu}_n - \Phi \mu^0 \right\|_{L^2(\mathbb{R}^d)}^2 \right] & \le e \mathbb{E} \left[ \left\| L \circ \Phi(\hat{\mu}_n - \mu^0) \right\|_{\mathbb{L}}^2 \right] + \frac{2}{(2\pi)^d} \left( \mathbb{E} \left[ \left\| \hat{\mu}_n \right\|_{\text{TV}}^2 \right] + \left\| \mu^0 \right\|_{\text{TV}}^2 \right) \frac{\tau^d d^{d/2}}{2^{d/2} u_{\min}^{2d}} e^{-2 \frac{u_{\min}^2}{\tau^2}}, \\ & \le 4e(\rho_n + \kappa \sqrt{c_p s})^2 \\ & \quad + \left( \frac{(2\pi)^{d/2}(2u_{\max}^2 + \tau^2)^{d/2}}{\kappa^2} \tilde{C}_{\Gamma} \rho_n^4 + 2 \left( 2 \left( \frac{u_{\max}}{u_{\min}} \right)^d + 1 \right) \left\| \mu^0 \right\|_{\text{TV}}^2 \right) \frac{(2\pi)^{-d} \tau^d d^{d/2}}{2^{d/2} u_{\min}^{2d}} e^{-2 \frac{u_{\min}^2}{\tau^2}}, \\ & \le 4e(\rho_n + \kappa \sqrt{c_p s})^2 \\ & \quad + \left( 4 \left( \frac{u_{\max}^2}{u_{\min}^2} + \frac{1}{\ln n} \right)^{d/2} \tilde{C}_{\Gamma} \frac{(\ln n)^{d/2}}{n \kappa^2} + 2 \left( 2 \left( \frac{u_{\max}}{u_{\min}} \right)^d + 1 \right) \left\| \mu^0 \right\|_{\text{TV}}^2 \right) \frac{(2\pi)^{-d} d^{d/2}}{u_{\min}^d (\ln n)^{d/2} n}. \end{aligned} $$

With the choice $\kappa = \frac{\rho_n}{\sqrt{c_p s}}$, it comes

$$ \begin{aligned} \mathbb{E} \left[ \left\| \Phi \hat{\mu}_n - \Phi \mu^0 \right\|_{L^2(\mathbb{R}^d)}^2 \right] & \le 4e \rho_n^2 (1 + \sqrt{s})^2 \\ & \quad + \left( 4 \left( \frac{u_{\max}^2}{u_{\min}^2} + \frac{1}{\ln n} \right)^{d/2} \tilde{C}_{\Gamma} c_p \frac{(\ln n)^{d/2}}{n} + 2 \|\mu^0\|_{\text{TV}}^2 \left( 2 \left( \frac{u_{\max}}{u_{\min}} \right)^d + 1 \right) \right) \frac{(2\pi)^{-d} d^{d/2}}{u_{\min}^d (\ln n)^{d/2} n}. \end{aligned} $$

This concludes the proof of Theorem 4.1. With the choice $\kappa = \frac{\rho_n}{\sqrt{c_p s}}$, we get

$$ \begin{aligned} \mathbb{E} \left[ \left\| \Phi \hat{\mu}_n - \Phi \mu^0 \right\|_{L^2(\mathbb{R}^d)}^2 \right] & \le 16e \rho_n^2 + \left( 4 \left( \frac{u_{\max}^2}{u_{\min}^2} + \frac{1}{\ln n} \right)^{d/2} \tilde{C}_{\Gamma} c_p \frac{s(\ln n)^{d/2}}{n} + 2 \|\mu^0\|_{\text{TV}}^2 \left( 2 \left( \frac{u_{\max}}{u_{\min}} \right)^d + 1 \right) \right) \frac{(2\pi)^{-d} d^{d/2}}{u_{\min}^d (\ln n)^{d/2} n}, \\ & \lesssim \left( \frac{s}{n(\ln n)^{d/2}} + 1 \right) \frac{(\ln n)^{d/2}}{n} \end{aligned} \quad (45) $$

keeping only the dependence on $n$ and $s$.

# H Properties of the Fisher-Rao metric

We present properties associated with the Fisher-Rao metric $g$, defined at point $x \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$ by $g_x = \nabla_1 \nabla_2 K_{\text{norm}}(x, x)$ (see also (20)). Note that this metric depends on the smoothing parameter $\tau > 0$ through $K_{\text{norm}}$. We recall the definition of the Riemannian norm: for $v \in \mathbb{R}^{2d}$ and $x \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$, we define $\|v\|_x = \sqrt{v^T g_x v}$.

## H.1 Christoffel symbols

The non-zero Christoffel symbols associated with $g$ are

$$ \begin{aligned} \Gamma^{t_k}_{u_k t_k} &= \Gamma^{t_k}_{t_k u_k} = \frac{-2u_k}{2u_k^2 + \tau^2}, \\ \Gamma^{u_k}_{t_k t_k} &= \frac{1}{u_k}, \\ \Gamma^{u_k}_{u_k u_k} &= \frac{\tau^2 - 2u_k^2}{u_k(2u_k^2 + \tau^2)}. \end{aligned} $$
with $k = 1, \dots, d$ (see [Giard, 2025, Section I.1]). We define

$$ \Gamma^{t_k} = \begin{pmatrix} (\Gamma^{t_k}_{t_l t_m})_{1 \le l, m \le d} & (\Gamma^{t_k}_{t_l u_m})_{1 \le l, m \le d} \\ (\Gamma^{t_k}_{u_l t_m})_{1 \le l, m \le d} & (\Gamma^{t_k}_{u_l u_m})_{1 \le l, m \le d} \end{pmatrix} \quad \text{and} \quad \Gamma^{u_k} = \begin{pmatrix} (\Gamma^{u_k}_{t_l t_m})_{1 \le l, m \le d} & (\Gamma^{u_k}_{t_l u_m})_{1 \le l, m \le d} \\ (\Gamma^{u_k}_{u_l t_m})_{1 \le l, m \le d} & (\Gamma^{u_k}_{u_l u_m})_{1 \le l, m \le d} \end{pmatrix}. \quad (46) $$

## H.2 Geodesics and geodesic distance

The next lemmas provide the parametrization of the geodesics associated with the Fisher-Rao metric. We denote $\tilde{\gamma}$ a geodesic parametrized by arc length connecting the points $x = \tilde{\gamma}(0)$, $x' = \tilde{\gamma}(l) \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$. The parameter $l$ is the geodesic distance between $x$ and $x'$, denoted by $\mathfrak{d}_g(x, x')$. We also denote $\gamma : y \in [0, 1] \mapsto \tilde{\gamma}(ly)$ (it is the geodesic such that $x = \gamma(0)$, $x' = \gamma(1)$) and $\dot{\gamma}$ its derivative.

We do not use the formula of the geodesic distance $\mathfrak{d}_g$ in this paper, but we give it in the next lemmas for information.

**Lemma H.1 (Geodesics of the Poincaré half-plane model).** The Poincaré half-plane is $\{x = (t, u) \in \mathbb{R} \times \mathbb{R}_+^*\}$, on which we consider the metric defined by $h_x = \begin{pmatrix} \frac{1}{u^2} & 0 \\ 0 & \frac{1}{u^2} \end{pmatrix}$ for all $x = (t, u) \in \mathbb{R} \times \mathbb{R}_+^*$. The associated norm is defined by $\|v\|_x = \sqrt{v^T h_x v}$ for $v \in \mathbb{R}^2$.

The Poincaré geodesics are circular arcs whose origin is on the axis $\{u = 0\}$ and straight vertical lines (parallel to $\{t = 0\}$).

A Poincaré geodesic parametrized by arc-length, denoted by $\tilde{h} = (\tilde{h}_t, \tilde{h}_u)$, is of the form

$$ \tilde{h} : y \in [0, l] \mapsto \left( \frac{\tanh(C_2 + y)}{C_1} + C_3, \frac{1}{\cosh(C_2 + y)|C_1|} \right) $$

(semicircle) or

$$ \tilde{h} : y \in [0, l] \mapsto (C_3, |C_1|e^y) \quad \text{or} \quad \tilde{h} : y \in [0, l] \mapsto (C_3, |C_1|e^{-y}) $$

(straight line), where $C_1 \in \mathbb{R}^*$, $C_2, C_3 \in \mathbb{R}$, $l \in \mathbb{R}^+$.

Moreover, writing $\mathfrak{d}_h$ the Poincaré distance,

$$ \mathfrak{d}_h(x, x') = \ln \left( \frac{\sqrt{(t-t')^2 + (u+u')^2} + \sqrt{(t-t')^2 + (u-u')^2}}{\sqrt{(t-t')^2 + (u+u')^2} - \sqrt{(t-t')^2 + (u-u')^2}} \right) \quad \forall x, x' \in \mathbb{R} \times \mathbb{R}_+^* $$

**Proof.** The fact that the Poincaré geodesics are semicircles whose origin is on the axis $\{u = 0\}$ and straight lines parallel to $\{t = 0\}$ is well-known (see for instance [Stahl, 1993, Theorem 4.2.1]). The formula for $\mathfrak{d}_h$ can be found in [Beardon, 1983, Theorem 7.2.1].

We can check that the parametrizations given for the geodesics verify the geodesic equations

$$ \left\{ \begin{aligned} \frac{(\dot{h}_t)^2 + (\dot{h}_u)^2}{\dot{h}_u^2} &= 1 \\ \ddot{h}_t - 2 \frac{\dot{h}_t \dot{h}_u}{\dot{h}_u} &= 0 \\ \ddot{h}_u - \frac{(\dot{h}_u)^2}{\dot{h}_u} + \frac{(\dot{h}_t)^2}{\dot{h}_u} &= 0 \end{aligned} \right. , $$

as done in [Giard, 2025, Section I.2]. We found all the geodesics, because all the portions of Poincaré semicircles and straight lines can be obtained with appropriate choices of $C_1, C_2, C_3, l$. $\square$

**Lemma H.2 (Geodesics for $d = 1$).** Let $x, x' \in \mathbb{R} \times [u_{\min}, +\infty)$.

* If $t = t'$, the geodesic is of the form

$$ \tilde{\gamma}(y) = \left( c_3, \sqrt{\frac{c_1^2}{2} e^{\sqrt{8}y} - \frac{\tau^2}{2}} \right) \quad \forall y \in [0, l] \quad \text{or} \quad \tilde{\gamma}(y) = \left( c_3, \sqrt{\frac{c_1^2}{2} e^{-\sqrt{8}y} - \frac{\tau^2}{2}} \right) \quad \forall y \in [0, l] \quad (47) $$

where $c_1, c_3 \in \mathbb{R}$. It is a portion of a straight line parallel to $\{t = 0\}$.

* If $t \neq t'$, the geodesic is of the form

$$ \tilde{\gamma}(y) = \left( c_3 + \frac{\sqrt{2} \tanh\left(\frac{c_2}{2} + \sqrt{2}y\right)}{2c_1}, \sqrt{-\frac{\tau^2}{2} + \frac{1 - \tanh^2\left(\frac{c_2}{2} + \sqrt{2}y\right)}{2c_1^2}} \right) \quad \forall y \in [0, l] \quad (48) $$
where $c_1 \neq 0$, $c_2, c_3 \in \mathbb{R}$. It is a portion of a semicircle with center $(c_3, 0)$ and radius

$$ \sqrt{\frac{1}{2c_1^2} - \frac{\tau^2}{2}} = \sqrt{\frac{1}{4} \left( \frac{-(t-t')^2 + u^2 - u'^2}{t' - t} \right)^2 + u^2}. \quad (49) $$

Moreover, the Fisher-Rao distance between $x$ and $x'$ is

$$ \mathcal{D}_g(x, x') = \sqrt{2} \ln \left( \frac{\sqrt{(t-t')^2 + \left(\sqrt{u^2 + \frac{\tau^2}{2}} - \sqrt{u'^2 + \frac{\tau^2}{2}}\right)^2} + \sqrt{(t-t')^2 + \left(\sqrt{u^2 + \frac{\tau^2}{2}} + \sqrt{u'^2 + \frac{\tau^2}{2}}\right)^2}}{\sqrt{2}(2u^2 + \tau^2)^{1/4}(2u'^2 + \tau^2)^{1/4}} \right). $$

**Proof.** Link with the Poincaré half-plane model: Recall that the variational formulation of a (Fisher-Rao) geodesic $\gamma = (\gamma_t, \gamma_u)$ connecting $x, x'$ is

$$ \inf_{\gamma(0)=x, \gamma(1)=x'} \int_0^1 \sqrt{\dot{\gamma}_t(y)^2 \frac{1}{2\gamma_u(y)^2 + \tau^2} + \dot{\gamma}_u(y)^2 \frac{2\gamma_u(y)^2}{(2\gamma_u(y)^2 + \tau^2)^2}} dy. $$

We use the change of variable $h = (h_t, h_u) = (\gamma_t, \sqrt{\gamma_u^2 + \frac{\tau^2}{2}})$. Noticing that $(\dot{h}_u)^2 = \frac{\dot{\gamma}_u(y)^2 \gamma_u(y)^2}{\gamma_u^2 + \frac{\tau^2}{2}}$, it comes that the variational formulation is equivalent to the problem

$$ \inf \left\{ \frac{1}{\sqrt{2}} \int_0^1 \sqrt{(\dot{h}_t(y))^2 \frac{1}{h_u(y)^2} + (\dot{h}_u(y))^2 \frac{1}{h_u(y)^2}} dy : h(0) = \left(t, \sqrt{u^2 + \frac{\tau^2}{2}}\right), h(1) = \left(t', \sqrt{u'^2 + \frac{\tau^2}{2}}\right) \right\}, $$

and we recognize the Poincaré metric tensor $(\frac{dt^2+du^2}{u^2})$ in this formulation. So $h$ is the geodesic for the Poincaré half-plane metric connecting $(t, \sqrt{u^2 + \frac{\tau^2}{2}})$ and $(t', \sqrt{u'^2 + \frac{\tau^2}{2}})$.

**Geodesic distance:** In particular, using the formula for $\mathcal{D}_h$ in Lemma H.1, we have

$$ \begin{aligned} \mathcal{D}_g(x, x') &= \frac{1}{\sqrt{2}} \mathcal{D}_h \left( \left( t, \sqrt{u^2 + \frac{\tau^2}{2}} \right), \left( t', \sqrt{u'^2 + \frac{\tau^2}{2}} \right) \right), \\ &= \frac{1}{\sqrt{2}} \ln \left( \frac{\sqrt{(t-t')^2 + \left(\sqrt{u^2 + \frac{\tau^2}{2}} + \sqrt{u'^2 + \frac{\tau^2}{2}}\right)^2} + \sqrt{(t-t')^2 + \left(\sqrt{u^2 + \frac{\tau^2}{2}} - \sqrt{u'^2 + \frac{\tau^2}{2}}\right)^2}}{\sqrt{(t-t')^2 + \left(\sqrt{u^2 + \frac{\tau^2}{2}} + \sqrt{u'^2 + \frac{\tau^2}{2}}\right)^2} - \sqrt{(t-t')^2 + \left(\sqrt{u^2 + \frac{\tau^2}{2}} - \sqrt{u'^2 + \frac{\tau^2}{2}}\right)^2}} \right), \\ &= \sqrt{2} \ln \left( \frac{\sqrt{(t-t')^2 + \left(\sqrt{u^2 + \frac{\tau^2}{2}} - \sqrt{u'^2 + \frac{\tau^2}{2}}\right)^2} + \sqrt{(t-t')^2 + \left(\sqrt{u^2 + \frac{\tau^2}{2}} + \sqrt{u'^2 + \frac{\tau^2}{2}}\right)^2}}{\sqrt{2}(2u^2 + \tau^2)^{1/4}(2u'^2 + \tau^2)^{1/4}} \right). \end{aligned} $$

**Parametrization by arc-length:** We saw that $(\gamma_t, \sqrt{\gamma_u^2 + \frac{\tau^2}{2}})$ is a geodesic for the Poincaré half-plane model. Lemma H.1 gives its parametrization by arc-length, $\tilde{h}$. For the $h$-norm, $\|\dot{\tilde{h}}_t(y), \dot{\tilde{h}}_u(y)\|_{\tilde{h}(y)} = 1$. We deduce that for the $g$-norm, writing $g(y) = (\tilde{h}_t(y), \sqrt{\tilde{h}_u(y)^2 - \frac{\tau^2}{2}})$, we have $\|\dot{g}_t(y), \dot{g}_u(y)\|_{g(y)} = \frac{1}{\sqrt{2}}$. So defining $\tilde{\gamma}(y) = g(\sqrt{2}y)$, it comes that $\|\dot{\tilde{\gamma}}_t(y), \dot{\tilde{\gamma}}_u(y)\|_{\gamma(y)} = 1$ for the $g$-norm: $\tilde{\gamma}$ is a geodesic parametrized by arc-length. The geodesics we provide are of this form.

**Radius of the semicircle connecting $x$ and $x'$:** The parametrization of the semicircle gives

$$ (\tilde{\gamma}_t - c_3)^2 + \tilde{\gamma}_u^2 = \frac{1}{2c_1^2} - \frac{\tau^2}{2} $$

which is the square of the radius of the semicircle. We also have

$$ \frac{1}{2c_1^2} - \frac{\tau^2}{2} = (t - c_3)^2 + u^2 = (t' - c_3)^2 + u'^2, $$
from which we deduce that $c_3 = \frac{t'^2+u'^2-(t^2+u^2)}{2(t'-t)}$ along with $\frac{1}{2c_1^2} - \frac{\tau^2}{2} = \frac{1}{4} \left( \frac{-(t-t')^2+u^2-u'^2}{t'-t} \right)^2 + u^2$. $\square$

We can extend this result to higher dimensions. By abuse of notation, $\mathfrak{D}_g(x_k, x'_k)$ will refer to the Fisher-Rao distance in dimension 1 between $x_k$ and $x'_k$, for all $k \in \{1, \dots, d\}$. The notation $\mathfrak{g}_{x_k}$ follows the same principle.

**Lemma H.3 (Geodesics in dimension $d \ge 1$).** Let $x, x' \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$. For $k = 1, \dots, d$, we denote $\tilde{\gamma}_k$ the geodesic in dimension 1, parametrized by arc length connecting $x_k = (t_k, u_k)$ and $x'_k = (t'_k, u'_k)$. The geodesic connecting $x$ and $x'$ is of the form

$$ \tilde{\gamma} = (\tilde{\gamma}_{t_1}, \dots, \tilde{\gamma}_{t_d}, \tilde{\gamma}_{u_1}, \dots, \tilde{\gamma}_{u_d}) \quad \text{where} \quad (\tilde{\gamma}_{t_k}(y), \tilde{\gamma}_{u_k}(y)) = \tilde{\gamma}_k(\sqrt{g_k}y) \quad \text{with} \quad g_k \ge 0 \quad \text{and} \quad \sum_{i=1}^d g_i = 1. $$

Moreover,

$$ \mathfrak{D}_g(x, x') = \sqrt{\sum_{k=1}^{d} \mathfrak{D}_g(x_k, x'_k)^2}. $$

*Proof.* For all $k \in \{1, \dots, d\}$, we denote $\mathfrak{D}_k = \mathfrak{D}_g(x_k, x'_k)^2$ and $\tilde{\gamma}_k(y) = \tilde{\gamma}_k \left( \frac{\sqrt{\mathfrak{D}_k}}{\sqrt{\sum_j \mathfrak{D}_j}} y \right)$. Then, for all $v \in \mathbb{R}^{2d}$ and $x \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$, we have

$$ \|v\|_x^2 = v^T \mathfrak{D}_x v = \sum_{k=1}^{d} (v_k, v_{k+d}) \mathfrak{D}_{x_k} (v_k, v_{k+d})^T, $$

where $(v_k, v_{k+d})$ represents the $k$-th component of $v$ in $\mathbb{R}^d \times [u_{\min}, +\infty)^d$. Using this, the geodesic $\tilde{\gamma}$ connecting $x$ and $x'$ can be expressed as:

$$ \tilde{\gamma} = (\tilde{\gamma}_{t_1}, \dots, \tilde{\gamma}_{t_d}, \tilde{\gamma}_{u_1}, \dots, \tilde{\gamma}_{u_d}), $$

where each $(\tilde{\gamma}_{t_k}, \tilde{\gamma}_{u_k})$ corresponds to the geodesic $\tilde{\gamma}_k$ in dimension 1 connecting $x_k = (t_k, u_k)$ and $x'_k = (t'_k, u'_k)$. This ensures that $\tilde{\gamma}$ is the geodesic connecting $x$ and $x'$ in $\mathbb{R}^d \times [u_{\min}, +\infty)^d$. $\square$

# H.3 Compatibility with the semi-distance

## H.3.1 Proof of Lemma 5.2

Let $r > 0$ and $x_0 \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$. Let $\gamma$ be a geodesic between $x_0 = \gamma(0)$ and $x = \gamma(1) \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$ for the metric $\mathfrak{g}$.

**Lower bound on the variance** We write $\gamma = (\gamma_{t_1}, \dots, \gamma_{t_d}, \gamma_{u_1}, \dots, \gamma_{u_d})$. Recall that $(\gamma_{t_k}, \gamma_{u_k})$ is a portion of a straight line parallel to $\{t_k = 0\}$ or of a semicircle with center on $\{u_k = 0\}$. As $x_0, x_1 \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$, we deduce that $u_{\min} \le \gamma_{u_k}(y)$ for all $k = 1, \dots, d, y \in [0, 1]$.

**The function $y \in [0, 1] \mapsto d(x_0, \gamma(y))$ is non-decreasing** We can reduce the problem to the case $d = 1$. Indeed, for each $k \in \{1, \dots, d\}$, $\gamma_k := (\gamma_{t_k}, \gamma_{u_k})$ is the geodesic in dimension 1 connecting $x_{0,k} = \gamma_k(0)$ and $x_k = \gamma_k(1)$ (see Lemma H.3). Furthermore, $d(x_0, \gamma(y))^2 = \sum_{k=1}^d d(x_{0,k}, \gamma_k(y))^2$, where, by abuse of notation, $d(x_{0,k}, \gamma_k(y))$ denotes the semi-distance in dimension 1 between $x_{0,k}$ and $\gamma_k(y)$. Hence, if the function $y \in [0, 1] \mapsto d(x_{0,k}, \gamma_k(y))$ is increasing for all $k$, then $y \in [0, 1] \mapsto d(x_0, \gamma(y))$ is also increasing.

Until the end of the proof, we therefore concentrate our attention on the case $d = 1$. Let $x_0, x \in \mathbb{R} \times [u_{\min}, +\infty)$, and let $\gamma$ and $\tilde{\gamma}$ be the geodesic connecting $x_0$ and $x$, parametrized by $[0, 1]$ and by arc-length, respectively. Proving that $y \in [0, 1] \mapsto d(x_0, \gamma(y))$ is increasing is equivalent to proving that $h: y \in [0, l] \mapsto d(x_0, \tilde{\gamma}(y))^2$ is increasing, where $l = \mathfrak{D}_g(x_0, x)$. We will consider alternatively the cases where the geodesic $\gamma$ is a straight line and a semicircle. It suffices to show that $h' \ge 0$ in both cases.

*Proof for straight lines, $d = 1$:* See [Giard, 2025, Section II.2]. Using the form of the geodesic given in Lemma H.2, we can show that if $\tilde{\gamma}$ is a straight line,

$$ h(y) = \ln(\cosh(\sqrt{2}y)). $$

We obtain this formula by deducing from $\tilde{\gamma}(0) = (t_0, u_0)$ in (47) that $c_1^2 = 2u_0^2 + \tau^2$ (the two formulas of (47) give the same result). This function is non-decreasing on $\mathbb{R}^+$.

*Proof for semicircles, $d = 1$:* See [Giard, 2025, Section II.1]. If $\tilde{\gamma}$ is a semicircle, then

$$ h(y) = \ln \left( \frac{\cosh^2\left(\frac{c_2}{2}\right) + \cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right)}{2 \cosh\left(\frac{c_2}{2}\right) \cosh\left(\frac{c_2}{2} + \sqrt{2}y\right)} \right) + \frac{\sinh^2(\sqrt{2}y)}{\cosh^2\left(\frac{c_2}{2}\right) + \cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right)}. $$
We obtain this formula by deducing from $\tilde{\gamma}(0) = (t_0, u_0)$ in (48) that $c_3 = t_0 - \frac{\sqrt{2} \tanh(\frac{c_2}{2})}{2c_1}$ and $c_1 = \pm \frac{1}{\sqrt{\tau^2 + 2u_0^2 \cosh(\frac{c_2}{2})}}$. Let $y \ge 0$. We have $h(y) = A(y) + B(y)$ with

$$A(y) = \frac{\sinh^2(\sqrt{2}y)}{\cosh^2(\frac{c_2}{2}) + \cosh^2(\frac{c_2}{2} + \sqrt{2}y)}.$$

and

$$B(y) = \ln \left( \frac{\cosh^2(\frac{c_2}{2}) + \cosh^2(\frac{c_2}{2} + \sqrt{2}y)}{2 \cosh(\frac{c_2}{2}) \cosh(\frac{c_2}{2} + \sqrt{2}y)} \right).$$

Then,

$$A'(y) = \frac{\sqrt{2} \sinh(\sqrt{2}y) (\cosh(c_2 + \sqrt{2}y) + \cosh(\sqrt{2}y)(2 + \cosh(c_2)))}{(\cosh^2(\frac{c_2}{2}) + \cosh^2(\frac{c_2}{2} + \sqrt{2}y))^2} \ge 0$$

and

$$B'(y) = \frac{\sqrt{2} \sinh\left(\frac{c_2}{2} + \sqrt{2}y\right) \left(\cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right) - \cosh^2\left(\frac{c_2}{2}\right)\right)}{\cosh\left(\frac{c_2}{2} + \sqrt{2}y\right) \left(\cosh^2\left(\frac{c_2}{2}\right) + \cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right)\right)}.$$

If $c_2 \ge 0$, then as $y \ge 0$, $\cosh^2(\frac{c_2}{2} + \sqrt{2}y) \ge \cosh^2(\frac{c_2}{2})$ so $B' \ge 0$, hence $d(x_0, \tilde{\gamma}(y))^2$ is increasing on $\mathbb{R}^+$. We now deal with the case $c_2 < 0$. Remark that $B'(y) > 0$ if and only if $y \notin [\frac{-c_2}{2\sqrt{2}}, \frac{-c_2}{\sqrt{2}}]$. We want to show that $A' + B' \ge 0$ on the interval $[\frac{-c_2}{2\sqrt{2}}, \frac{-c_2}{\sqrt{2}}]$. First, using that for $y \in [\frac{-c_2}{2\sqrt{2}}, \frac{-c_2}{\sqrt{2}}]$, we have

$$\cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right) \le \cosh^2\left(\frac{c_2}{2}\right), \quad \sinh(\sqrt{2}y) \ge \sinh\left(\frac{-c_2}{2}\right), \\ \cosh(\sqrt{2}y) \ge \cosh\left(\frac{c_2}{2}\right) \quad \text{and} \quad 2 + \cosh(c_2) \ge 2 \cosh^2\left(\frac{c_2}{2}\right),$$

we obtain a lower bound for $A'$ on $[\frac{-c_2}{2\sqrt{2}}, \frac{-c_2}{\sqrt{2}}]$:

$$\begin{aligned} A'(y) &\ge \frac{\sqrt{2} \sinh\left(\frac{-c_2}{2}\right) \left(\cosh(c_2 + \sqrt{2}y) + \cosh(\sqrt{2}y)(2 + \cosh(c_2))\right)}{2 \cosh^2\left(\frac{c_2}{2}\right) \left(\cosh^2\left(\frac{c_2}{2}\right) + \cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right)\right)}, \\ &\ge \frac{\sqrt{2} \sinh\left(\frac{-c_2}{2}\right) \cosh\left(\frac{c_2}{2}\right) (2 + \cosh(c_2))}{2 \cosh^2\left(\frac{c_2}{2}\right) \left(\cosh^2\left(\frac{c_2}{2}\right) + \cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right)\right)}, \\ &\ge \frac{\sqrt{2} \sinh\left(\frac{-c_2}{2}\right) \cosh\left(\frac{c_2}{2}\right)}{\cosh^2\left(\frac{c_2}{2}\right) + \cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right)}. \end{aligned}$$

It follows that

$$A'(y) + B'(y) \ge \frac{\sqrt{2} \cosh\left(\frac{c_2}{2} + \sqrt{2}y\right) \sinh\left(\frac{-c_2}{2}\right) \cosh\left(\frac{c_2}{2}\right) + \sqrt{2} \sinh\left(\frac{c_2}{2} + \sqrt{2}y\right) \left(\cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right) - \cosh^2\left(\frac{c_2}{2}\right)\right)}{\left(\cosh^2\left(\frac{c_2}{2}\right) + \cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right)\right) \cosh\left(\frac{c_2}{2} + \sqrt{2}y\right)}.$$

The positivity of $A' + B'$ follows, since

$$\begin{aligned} & \cosh\left(\frac{c_2}{2} + \sqrt{2}y\right) \sinh\left(\frac{-c_2}{2}\right) \cosh\left(\frac{c_2}{2}\right) + \sinh\left(\frac{c_2}{2} + \sqrt{2}y\right) \left(\cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right) - \cosh^2\left(\frac{c_2}{2}\right)\right) \\ &= \frac{1}{4} \left(\sinh\left(\frac{3c_2}{2} + 3\sqrt{2}y\right) - 2\sinh\left(\frac{3c_2}{2} + \sqrt{2}y\right) - \sinh\left(\frac{c_2}{2} + \sqrt{2}y\right)\right), \\ & \ge \frac{1}{2} \left(\sinh\left(\frac{c_2}{2} + \sqrt{2}y\right) - \sinh\left(\frac{3c_2}{2} + \sqrt{2}y\right)\right), \\ & \ge 0, \end{aligned}$$

where we used

$$\begin{aligned} \sinh\left(\frac{3c_2}{2} + 3\sqrt{2}y\right) &= \sinh\left(\frac{c_2}{2} + \sqrt{2}y\right) \cosh(c_2 + 2\sqrt{2}y) + \cosh\left(\frac{c_2}{2} + \sqrt{2}y\right) \sinh(c_2 + 2\sqrt{2}y), \\ &\ge \sinh\left(\frac{c_2}{2} + \sqrt{2}y\right) + \sinh(c_2 + 2\sqrt{2}y), \end{aligned}$$
$$ \sinh(c_2 + 2\sqrt{2y}) = 2\sinh\left(\frac{c_2}{2} + \sqrt{2y}\right) \cosh\left(\frac{c_2}{2} + \sqrt{2y}\right) \geq 2\sinh\left(\frac{c_2}{2} + \sqrt{2y}\right) $$

and $c_2 < 0$ together with the fact that sinh is increasing.

We have just proved the following statement: if $d(x_0, x) \le r$, then

$$ \forall y \in [0, 1], \quad \gamma(y) \in \{x' \in \mathbb{R}^d \times [u_{\min}, +\infty)^d : d(x_0, x') \le r\}. $$

This concludes the proof of Lemma 5.2.

### H.3.2 Proof of Lemma 5.3

Let $r, \Delta > 0$. Assume that $\mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, u_{\max}]^d$. The semi-distance $d$ is defined by (11).

Control of the variance of the geodesics: Let $x_j^0 \in \mathcal{X}$. For $x \in B_d(x_j^0, r) \cap \mathbb{R}^d \times [u_{\min}, u_{\max}]^d$, we have $|t_k - t_{j,k}^0| \le r\sqrt{2u_{\max}^2 + \tau^2}$ for all $k \in \{1, \dots, d\}$. We deduce that

$$ B_d(x_j^0, r) \cap (\mathbb{R}^d \times [u_{\min}, u_{\max}]^d) \subset \left( \bigtimes_{k=1}^d [t_{j,k}^0 \pm r\sqrt{2u_{\max}^2 + \tau^2}] \right) \times [u_{\min}, u_{\max}]^d $$

where $\times_{k=1}^d$ denotes the $d$-ary Cartesian product.

For any $k \in \{1, \dots, d\}$, the semicircle $\gamma_k = (\gamma_{t_k}, \gamma_{u_k})$ (geodesic in dimension 1) connecting the points $(t_{j,k}^0 - r\sqrt{2u_{\max}^2 + \tau^2}, u_{\max})$ and $(t_{j,k}^0, u_{\max})$ satisfies $\gamma_{u_k}(y)^2 \le u_{\max}^2 + \frac{1}{4}r^2(2u_{\max}^2 + \tau^2)$. This follows because the square of the radius of the semicircle is $u_{\max}^2 + \frac{1}{4}r^2(2u_{\max}^2 + \tau^2)$ (see Lemma H.2). The same holds for the semicircle connecting $(t_{j,k}^0 + r\sqrt{2u_{\max}^2 + \tau^2}, u_{\max})$ and $(t_{j,k}^0, u_{\max})$.

Any geodesic connecting $(t_{j,k}^0, u_{j,k}^0)$ with a point of $[t_{j,k}^0 - r\sqrt{2u_{\max}^2 + \tau^2}, t_{j,k}^0 + r\sqrt{2u_{\max}^2 + \tau^2}] \times [u_{\min}, u_{\max}]$ is below these geodesics.

To see why it is the case, one can note that 2 geodesics in our model have at most 1 intersection point, or their union is a geodesic. So a geodesic between points of $[t_{j,k}^0 - r\sqrt{2u_{\max}^2 + \tau^2}, t_{j,k}^0] \times [u_{\min}, u_{\max}]$ cannot intersect twice the semicircle connecting $(t_{j,k}^0 - r\sqrt{2u_{\max}^2 + \tau^2}, u_{\max})$ and $(t_{j,k}^0, u_{\max})$. An analogous statement holds when considering $[t_{j,k}^0, t_{j,k}^0 + r\sqrt{2u_{\max}^2 + \tau^2}] \times [u_{\min}, u_{\max}]$.

Using Lemma H.3, we deduce that

$$ \mathcal{G}_{x_j^0}(B_d(x_j^0, r) \cap (\mathbb{R}^d \times [u_{\min}, u_{\max}]^d)) \subset \mathbb{R}^d \times \left[ u_{\min}, \sqrt{u_{\max}^2 + \frac{1}{4}r^2(2u_{\max}^2 + \tau^2)} \right]^d. $$

Pseudo-quasi triangle inequality: Let $x, x', \tilde{x} \in \mathbb{R}^d \times [u_{\min}, \tilde{u}_{\max}]^d$, where $\tilde{u}_{\max}$ is not necessarily equal to $u_{\max}$. We have

$$ \frac{u_k^2 + u_k'^2 + \tau^2}{\sqrt{2u_k^2 + \tau^2}\sqrt{2u_k'^2 + \tau^2}} \le \frac{2\tilde{u}_{\max}^2 + \tau^2}{2u_{\min}^2 + \tau^2} \le \frac{\tilde{u}_{\max}^2}{u_{\min}^2} $$

using that $\tau \in \mathbb{R}^+ \mapsto \frac{2\tilde{u}_{\max}^2 + \tau^2}{2u_{\min}^2 + \tau^2}$ is decreasing (because $u_{\min} \le \tilde{u}_{\max}$). Moreover, using the triangle inequality for the $l^2$-norm in dimension $d$ and that $u_k^2 + \tilde{u}_k^2 + \tau^2$ and $u_k'^2 + \tilde{u}_k^2 + \tau^2$ are smaller than $\frac{\tilde{u}_{\max}^2}{u_{\min}^2}(u_k^2 + u_k'^2 + \tau^2)$, we get

$$ \begin{aligned} \sqrt{\sum_{k=1}^d \frac{(t_k - t_k')^2}{u_k^2 + u_k'^2 + \tau^2}} &\le \sqrt{\sum_{k=1}^d \frac{(t_k - \tilde{t}_k)^2}{u_k^2 + \tilde{u}_k'^2 + \tau^2}} + \sqrt{\sum_{k=1}^d \frac{(t_k' - \tilde{t}_k)^2}{u_k^2 + \tilde{u}_k'^2 + \tau^2}}, \\ &\le \frac{\tilde{u}_{\max}}{u_{\min}} \left( \sqrt{\sum_{k=1}^d \frac{(t_k - \tilde{t}_k)^2}{u_k^2 + \tilde{u}_k^2 + \tau^2}} + \sqrt{\sum_{k=1}^d \frac{(t_k' - \tilde{t}_k)^2}{u_k'^2 + \tilde{u}_k^2 + \tau^2}} \right), \\ &\le \frac{\tilde{u}_{\max}}{u_{\min}} (d(x, \tilde{x}) + d(x', \tilde{x})). \end{aligned} $$

Hence

$$ d(x, x') \le \frac{\tilde{u}_{\max}}{u_{\min}} (d(x, \tilde{x}) + d(x', \tilde{x})) + \sqrt{d \ln \left( \frac{\tilde{u}_{\max}^2}{u_{\min}^2} \right)}. $$

$$ \text{So } d(x, \tilde{x}) \ge \frac{u_{\min}}{\tilde{u}_{\max}} \left( d(x, x') - \sqrt{d \ln \left( \frac{\tilde{u}_{\max}^2}{u_{\min}^2} \right)} - d(x', \tilde{x}) \right). $$

Conclusion: Let $x_j^0, x_i^0 \in \mathcal{X}$. Firstly, replacing $\tilde{u}_{\max}$ by $u_{\max}$ in the pseudo-quasi triangle inequality, for $\tilde{x} \in$
$B_d(x_j^0, \Delta) \cap \mathcal{X}$, if $\frac{u_{\min}}{u_{\max}} \left( d(x_i^0, x_j^0) - \sqrt{d \ln\left(\frac{u_{\max}^2}{u_{\min}^2}\right)} \right) - d(\tilde{x}, x_j^0) > \Delta$ we have $d(x_i^0, \tilde{x}) > \Delta$. Note also that $d(\tilde{x}, x_j^0) < \Delta$, hence it suffices that $d(x_j^0, x_i^0) \ge 2\frac{u_{\max}}{u_{\min}}\Delta + \sqrt{d \ln\left(\frac{u_{\max}^2}{u_{\min}^2}\right)}$ for the open balls $\tilde{B}_d(x_j^0, \Delta) \cap \mathcal{X}$ to be disjoint.

Secondly, for $\tilde{x} \in \mathcal{G}_{x_j^0}(B_d(x_j^0, r) \cap \mathcal{X})$, taking $\tilde{u}_{\max} = \sqrt{u_{\max}^2 + \frac{1}{4}r^2(2u_{\max}^2 + \tau^2)}$ we have $d(x_i^0, \tilde{x}) \ge \Delta$ as soon as $\frac{u_{\min}}{\tilde{u}_{\max}} \left( d(x_i^0, x_j^0) - \sqrt{d \ln\left(\frac{\tilde{u}_{\max}^2}{u_{\min}^2}\right)} \right) \ge \Delta + r$. We used that $\mathcal{G}_{x_j^0}(B_d(x_j^0, r)) \subset B_d(x_j^0, r)$, so $d(\tilde{x}, x_j^0) \le r$ (see Lemma 5.2). This condition can be rewritten as $d(x_i^0, x_j^0) \ge \frac{\tilde{u}_{\max}}{u_{\min}}(\Delta + r) + \sqrt{d \ln\left(\frac{u_{\max}^2}{u_{\min}^2}\right)}$. We proved that if

$$ \min_{i \neq j} d(x_i^0, x_j^0) \geq \max \left\{ \frac{\sqrt{u_{\max}^2 + \frac{1}{4}r^2(2u_{\max}^2 + \tau^2)}}{u_{\min}} (\Delta + r), 2\frac{u_{\max}}{u_{\min}}\Delta \right\} + \sqrt{d \ln\left(\frac{u_{\max}^2}{u_{\min}^2}\right)} = \Delta_{\tau}, $$

then the open balls $\tilde{B}_d(x_j^0, \Delta) \cap \mathcal{X}$ are disjoint, and for all $j \neq i$, the ball $\mathcal{G}_{x_j^0}(B_d(x_j^0, r) \cap \mathcal{X})$ does not intersect $\tilde{B}_d(x_i^0, \Delta)$.

### H.3.3 Local control with the semi-distance

**Lemma H.4 (Local control of $\mathfrak{d}_g$ with the semi-distance, lower bound, $d=1$).** Let $x_0, x \in \mathbb{R} \times [u_{\min}, +\infty)$. If $d(x_0, x) \le r$, then

$$ \mathfrak{d}_g(x_0, x)^2 \geq \frac{d(x_0, x)^2}{\hat{\varepsilon}_3} $$

for $\hat{\varepsilon}_3 \ge 1 + \frac{1}{R(r)}$ with $R(r)$ defined by (57) below.

*Proof.* Let $x_0 = (t_0, u_0)$, $x = (t, u) \in \mathbb{R} \times [u_{\min}, +\infty)$ such that $r_e := d(x_0, x) \le r$. Without loss of generality, we assume that $u_0 \le u$.

We examine the case where $u$ is “far” from $u_0$, and then the case where $t$ is far from $t_0$ and $u, u_0$ are close. Bounds on $B_d(x_0, r)$, exhibition of 2 cases: Since $d(x_0, x)^2 = \frac{(t_0-t)^2}{u^2+u_0^2+\tau^2} + \ln\left(\frac{u^2+u_0^2+\tau^2}{\sqrt{2u^2+\tau^2}\sqrt{2u_0^2+\tau^2}}\right) = r_e^2$, we have that for all $0 < w < 1$,

$$ \frac{(t_0 - t)^2}{u^2 + u_0^2 + \tau^2} \geq w r_e^2 \qquad (50) $$

(first case) or

$$ \frac{u^2 + u_0^2 + \tau^2}{\sqrt{2u^2 + \tau^2}\sqrt{2u_0^2 + \tau^2}} \geq e^{(1-w)r_e^2} $$

(second case). The second case can be rewritten as

$$ \sqrt{u^2 + \frac{\tau^2}{2}} \geq \sqrt{u_0^2 + \frac{\tau^2}{2}}(e^{(1-w)r_e^2} + \sqrt{e^{2(1-w)r_e^2} - 1}) \qquad (51) $$

(see (34)).

Assuming the first case holds and the second case does not: Provided that (50) holds, then $t_0 \neq t$ and the geodesic connecting $x$ and $x_0$ is a semicircle (see Lemma H.2). We suppose additionally that we are not in the second case, namely that

$$ \sqrt{u^2 + \frac{\tau^2}{2}} \leq \sqrt{u_0^2 + \frac{\tau^2}{2}}(e^{(1-w)r_e^2} + \sqrt{e^{2(1-w)r_e^2} - 1}). $$

In particular, using (35) we have

$$ \frac{|u_0^2 - u^2|}{u_0^2 + u^2 + \tau^2} \leq \sqrt{e^{2(1-w)r_e^2} - 1}. \qquad (52) $$

Denoting $l$ the arc length (i.e. $l = \mathfrak{d}_g(x, x_0)$), using (48) and (49) we have

$$ \frac{|\tanh(\frac{c_2}{2} + \sqrt{2l}) - \tanh(\frac{c_2}{2})|}{\sqrt{2}} = |c_1||t - t_0| \qquad (53) $$
for some $c_2$ and $c_1$ verifying

$$ \frac{1}{2c_1^2} - \frac{\tau^2}{2} = \frac{1}{4} \left( \frac{-(t-t_0)^2 + u_0^2 - u^2}{t-t_0} \right)^2 + u_0^2. $$

According to (50) and (52),

$$ \begin{aligned} \frac{1}{|c_1||t-t_0|} &= \sqrt{\frac{1}{2}\left(-1 + \frac{u_0^2 - u^2}{(t-t_0)^2}\right)^2 + \frac{2u_0^2 + \tau^2}{(t-t_0)^2}}, \\ &= \sqrt{\frac{1}{2} + \frac{1}{2}\frac{(u_0^2 - u^2)^2}{(t-t_0)^4} + \frac{u_0^2 + u^2 + \tau^2}{(t-t_0)^2}}, \quad (54) \\ &\le \sqrt{\frac{1}{2} + \frac{e^{2(1-w)r_e^2} - 1}{2w^2r_e^4} + \frac{1}{wr_e^2}}, \\ &\le \frac{1}{\sqrt{wr_e}}\sqrt{\frac{r^2}{2} + \frac{1-w}{w}\frac{e^{2r^2}-1}{2r^2} + 1} \end{aligned} $$

where we used that $y \mapsto e^{2y} - 1$ is convex along with $(1-w)r_e^2 \le r^2$, and $wr_e^2 \le r^2$.

So as tanh is 1-Lipschitz,

$$ \mathcal{D}_g(x, x_0) = l \ge |c_1||t - t_0|, \\ \ge \frac{r_e \sqrt{w}}{\sqrt{1 + \frac{r^2}{2} + \frac{1-w}{w} \frac{e^{2r^2}-1}{2r^2}}}. \quad (55) $$

Assuming the second case holds: If (51) holds, we must consider the 2 possible geodesics (c.f. Lemma H.2).

If $t = t_0$, the geodesic is of the form $\left(c_3, \sqrt{-\frac{\tau^2}{2} + \frac{c_2^2}{2}e^{\sqrt{8}y}}\right)$. As $\tilde{\gamma}(0) = x_0$ we have $c_1^2 = 2u_0^2 + \tau^2$. Since $\tilde{\gamma}_u(l) = u$, we have $e^{\sqrt{8}l} = \frac{2u_0^2 + \tau^2}{2u_0^2 + \tau^2}$ from which we deduce that

$$ \mathcal{D}_g(x_0, x) = l = \frac{1}{\sqrt{2}} \ln \left( \frac{\sqrt{u^2 + \frac{\tau^2}{2}}}{\sqrt{u_0^2 + \frac{\tau^2}{2}}} \right) \ge \frac{1}{\sqrt{2}} \ln \left( e^{(1-w)r_e^2} + \sqrt{e^{2(1-w)r_e^2} - 1} \right). $$

If $t \ne t_0$, $\tilde{\gamma}_u$ is of the form $\sqrt{-\frac{\tau^2}{2} + \frac{1}{2\cosh^2\left(\frac{c_2}{2} + \sqrt{2}y\right)c_1^2}}$. Using that $\tilde{\gamma}_u(0) = u_0$ and $\tilde{\gamma}_u(l) = u$, we get $\frac{\cosh^2\left(\frac{c_2}{2}\right)}{\cosh^2\left(\frac{c_2}{2} + \sqrt{2}l\right)} = \frac{2u_0^2 + \tau^2}{2u_0^2 + \tau^2}$, from which we deduce using (51) that

$$ \frac{\cosh\left(\frac{c_2}{2}\right)}{\cosh\left(\frac{c_2}{2} + \sqrt{2}l\right)} \ge e^{(1-w)r_e^2} + \sqrt{e^{2(1-w)r_e^2} - 1}. $$

As

$$ \frac{\cosh\left(\frac{c_2}{2}\right)}{\cosh\left(\frac{c_2}{2} + \sqrt{2}l\right)} = \cosh(\sqrt{2}l) - \frac{\sinh\left(\frac{c_2}{2}\right)\sinh(\sqrt{2}l)}{\cosh\left(\frac{c_2}{2}\right)} \le \cosh(\sqrt{2}l) + |\sinh(\sqrt{2}l)| = e^{\sqrt{2}l}|, $$

it comes that $e^{\sqrt{2}l} \ge e^{(1-w)r_e^2} + \sqrt{e^{2(1-w)r_e^2} - 1}$ leading again to

$$ l \ge \frac{1}{\sqrt{2}} \ln \left( e^{(1-w)r_e^2} + \sqrt{e^{2(1-w)r_e^2} - 1} \right). $$

As $\frac{1}{\sqrt{2}} \ln(e^{y^2} + \sqrt{e^{2y^2} - 1}) \ge y$ for all $y \in \mathbb{R}$, we get for this second case

$$ \mathcal{D}_g(x, x_0) \ge \sqrt{1 - wr_e}. \quad (56) $$

**Conclusion:** We now choose $0 < w < 1$ depending only on $r$ to balance the bounds of the two cases (namely, (55) and (56)). Writing $X = \frac{1-w}{w}$, we want to pick $w$ such that

$$ \begin{aligned} \sqrt{1-w} &= \frac{\sqrt{w}}{\sqrt{1 + \frac{r^2}{2} + \frac{1-w}{w} \frac{e^{2r^2}-1}{2r^2}}} \\ &\iff \frac{e^{2r^2}-1}{2r^2} X^2 + \left(1 + \frac{r^2}{2}\right) X - 1 = 0. \end{aligned} $$
This polynomial has exactly 1 positive root

$$R(r) = \frac{-(1 + \frac{r^2}{2}) + \sqrt{(1 + \frac{r^2}{2})^2 + 2 \frac{e^{2r^2} - 1}{r^2}}}{\frac{e^{2r^2} - 1}{r^2}}. \qquad (57)$$

We can take $w_r = \frac{1}{R(r)+1}$. Then

$$\mathfrak{D}_{\mathfrak{g}}(x_0, x) \geq \sqrt{1 - w_r} r_e = \sqrt{\frac{R(r)}{1 + R(r)}} r_e.$$

**Lemma H.5** (Local control of $\mathfrak{D}_g$ with the semi-distance, lower bound, $d \ge 1$). Let $x_0, x \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$. If $d(x_0, x) \le r$, then

$$\mathfrak{D}_{\mathfrak{g}}(x_0, x)^2 \geq \frac{d(x_0, x)^2}{\tilde{\varepsilon}_3}$$

for $\tilde{\varepsilon}_3 \ge 1 + \frac{1}{R(r)}$ with $R(r)$ defined by (57).

*Proof.* Let $x_0, x \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$ such that $d(x_0, x) \le r$. For all $k \in \{1, \dots, d\}$, $d(x_{0,k}, x_k) \le r$. Lemma H.4 gives $\mathfrak{D}_{\mathfrak{g}}(x_k, x_{0,k})^2 \ge \frac{d(x_{0,k}, x_k)^2}{1 + \frac{1}{R(r)}}$. We conclude using that

$$\mathfrak{D}_{\mathfrak{g}}(x, x_0)^2 = \sum_{k=1}^{d} \mathfrak{D}_{\mathfrak{g}}(x_k, x_{0,k})^2 \ge \sum_{k=1}^{d} \frac{d(x_{0,k}, x_k)^2}{1 + \frac{1}{R(r)}} = \frac{d(x_0, x)^2}{1 + \frac{1}{R(r)}}.$$

**Lemma H.6** (Local control of $\mathfrak{D}_g$ with the semi-distance, upper bound, $d = 1$). Let $x, x_0 \in \mathbb{R} \times [u_{\min}, +\infty)$. If $d(x, x_0) < \sqrt{2}$, then $\mathfrak{D}_g(x, x_0)^2 \le F(d(x, x_0))$ (F defined by (62) below).

*Proof.* Let $x, x_0 \in \mathbb{R} \times [u_{\min}, +\infty)$. Without loss of generality, we can assume that $u_0 \le u$. We denote $l = \mathfrak{D}_g(x, x_0)$ and $r_e = d(x, x_0)$. The proof is close to that of Lemma H.4.

**First case:** If $t \ne t_0$, the geodesic connecting $x$ and $x_0$ is a semicircle. Recalling (53) and (54), we have

$$\frac{|\tanh(\frac{c_2}{2} + \sqrt{2}l) - \tanh(\frac{c_2}{2})|}{\sqrt{2}|c_1|} = |t - t_0|$$

for some $c_2$ and $c_1$ verifying

$$\frac{1}{|c_1| |t - t_0|} = \sqrt{\frac{1}{2} + \frac{1}{2} \frac{(u_0^2 - u^2)^2}{(t - t_0)^4} + \frac{u_0^2 + u^2 + \tau^2}{(t - t_0)^2}}, \\ \ge \frac{1}{r_e}.$$

Hence

$$|\tanh\left(\frac{c_2}{2} + \sqrt{2}l\right) - \tanh\left(\frac{c_2}{2}\right)| \le \sqrt{2}r_e. \qquad (58)$$

Furthermore, using $\frac{u_0^2 + u^2 + \tau^2}{\sqrt{2u^2 + \tau^2}\sqrt{2u_0^2 + \tau^2}} \le e^{r_e^2}$ together with (34) and recalling that $u_0 \le u$, we have

$$1 \le \frac{\cosh\left(\frac{c_2}{2}\right)}{\cosh\left(\frac{c_2}{2} + \sqrt{2}l\right)} = \sqrt{\frac{2u^2 + \tau^2}{2u_0^2 + \tau^2}} \le e^{r_e^2} + \sqrt{e^{2r_e^2} - 1}. \qquad (59)$$

Using (59) along with the equality $\frac{\cosh(A+B)}{\cosh(A)} = \cosh(B) + \tanh(A)\sinh(B)$, we deduce that

$$\cosh(\sqrt{2}l) + \tanh(\frac{c_2}{2}) \sinh(\sqrt{2}l) \le 1$$

and that

$$\cosh(\sqrt{2}l) - \tanh\left(\frac{c_2}{2} + \sqrt{2}l\right) \sinh(\sqrt{2}l) \le e^{r_e^2} + \sqrt{e^{2r_e^2} - 1}.$$

It comes

$$2 \cosh(\sqrt{2}l) + \left(\tanh\left(\frac{c_2}{2}\right) - \tanh\left(\frac{c_2}{2} + \sqrt{2}l\right)\right) \sinh(\sqrt{2}l) \le e^{r_e^2} + \sqrt{e^{2r_e^2} - 1} + 1,$$
so

$$2 \cosh(\sqrt{2}l) - \sinh(\sqrt{2}l) \sqrt{2}r_e \le e^{r_e^2} + \sqrt{e^{2r_e^2} - 1} + 1$$

where we have used (58). As $\cosh \ge 0$ and $\frac{\sinh(\sqrt{2}l)}{\cosh(\sqrt{2}l)} \le 1$, we get

$$\cosh(\sqrt{2}l) \left(1 - \frac{1}{\sqrt{2}}r_e\right) \le \frac{1}{2} \left(e^{r_e^2} + \sqrt{e^{2r_e^2} - 1} + 1\right).$$

For $r_e < \sqrt{2}$ we get

$$l^2 \le \frac{1}{2} \operatorname{arcosh} \left( \frac{\frac{1}{2} \left( e^{r_e^2} + \sqrt{e^{2r_e^2} - 1} + 1 \right)}{1 - \frac{1}{\sqrt{2}} r_e} \right)^2 . \quad (60)$$

**Second case:** If $t_0 = t$, using (47) we get $e^{\sqrt{8}l} = \frac{2u^2+\tau^2}{2u_0^2+\tau^2} \le (e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^2$ so

$$l^2 \le \frac{1}{2} \ln(e^{r_e^2} + \sqrt{e^{2r_e^2} - 1})^2 = \frac{1}{2} \operatorname{arcosh}(e^{r_e^2})^2. \quad (61)$$

**Comparison of the bounds found, conclusion:** Let $0 \le r_e < \sqrt{2}$. As $1 - \frac{1}{\sqrt{2}}r_e \le 1$ and $\sqrt{e^{2r_e^2} - 1} \ge e^{r_e^2} - 1$, the bound found for the first case (see (60)) is larger than the one found for the second case (see (61)).

So if $d(x, x_0) < \sqrt{2}$, $\mathfrak{D}_g(x, x_0)^2 \le F(d(x, x_0))$ with

$$F: y \in \mathbb{R}^+ \mapsto \frac{1}{2} \operatorname{arcosh} \left( \frac{\frac{1}{2} \left( e^{y^2} + \sqrt{e^{2y^2} - 1} + 1 \right)}{1 - \frac{1}{\sqrt{2}} y} \right)^2 . \quad (62)$$

**Lemma H.7** (Local control of $\mathfrak{D}_g$ with the semi-distance, upper bound, $d \ge 1$). Let $x, x_0 \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$. If $d(x, x_0) < \sqrt{2}$, then $\mathfrak{D}_g(x, x_0)^2 \le dF(d(x, x_0))$ ($F$ defined by (62)).

*Proof.* Let $x, x_0 \in \mathbb{R}^d \times [u_{\min}, +\infty)^d$ such that $d(x, x_0) = r_e < \sqrt{2}$. As for all $k = 1, \dots, d$, $d(x_k, x_{0,k}) \le d(x, x_0)$, by Lemma H.6 it comes that

$$\begin{aligned} \mathfrak{D}_g(x, x_0)^2 &= \sum_{k=1}^{d} \mathfrak{D}_g(x_k, x_{0,k})^2 \\ &\le \sum_{k=1}^{d} F(d(x_k, x_{0,k})), \\ &\le dF(d(x, x_0)) \end{aligned}$$

where we used that $F$ is non-decreasing.

## H.4 Proof of Lemma 5.4

We use Lemma H.5. It comes that we can choose $\tilde{\epsilon}_3 = 1 + \frac{1}{R(\frac{0.3025}{\sqrt{d}})}$ with $R$ defined by (57) in the appendix. We aim to find a more interpretable parameter. Remark that $r \in \mathbb{R}^+ \mapsto \frac{e^{2r^2}-1}{r^2}$ is increasing, along with $\frac{e^{2r^2}-1}{r^2} \ge 2$. We deduce that for all $0 < r \le 0.3025$, $R(r) \ge \frac{-(1+0.3025^2)+\sqrt{5}}{\frac{e^2 \times 0.3025^2-1}{0.3025^2}}$. So $1 + \frac{1}{R(\frac{0.3025}{\sqrt{d}})} \le 2.84$ (see [Giard, 2025, Section III]).

# I Proof of Theorem 5.2

**Construction of the certificates by solving a linear system** We give explicit formulas for $\eta, \eta_j$ of Theorem 5.2. These certificates are of the form (18). In the following, we write $\eta = \eta_{\alpha,\beta}$ and $\eta_j = \eta_{\alpha^j, \beta^j}$. Recalling Definitions 3.2 and 3.3, we want that for all $j = 1, \dots, s$, $\eta_{\alpha,\beta}(x_j^0) = 1$ and $\nabla\eta_{\alpha,\beta}(x_j^0) = 0_{2d}$. We also want $\eta_{\alpha^j, \beta^j}(x_j^0) = 1$, $\eta_{\alpha^j, \beta^j}(x_l^0) = 0$ for all $l \ne j$, and $\nabla\eta_{\alpha^j, \beta^j}(x_l^0) = 0_{2d}$.
These constraints translate to linear systems. Writing $e_j$ the vector of size $s$ containing a 1 at position $j$, and zeros elsewhere, along with $1_s = \sum_{j=1}^s e_j$, we want to solve

$$ \Upsilon \begin{pmatrix} \alpha \\ \beta \end{pmatrix} = \begin{pmatrix} 1_s \\ 0_{s2d} \end{pmatrix} := \mathbf{u}_s \quad \text{and} \quad \Upsilon \begin{pmatrix} \alpha^j \\ \beta^j \end{pmatrix} = \begin{pmatrix} e_j \\ 0_{s2d} \end{pmatrix} := \mathbf{u}_s^j \quad \forall j = 1, \dots, s \tag{63} $$

$$ \text{where } \Upsilon = \begin{pmatrix} (K_{\text{norm}}(x_i^0, x_j^0))_{i,j=1,\dots,s} & (\nabla_1 K_{\text{norm}}(x_i^0, x_j^0))_{i,j=1,\dots,s}^T \\ (\nabla_2 K_{\text{norm}}(x_i^0, x_j^0))_{i,j=1,\dots,s} & (\nabla_1 \nabla_2 K_{\text{norm}}(x_i^0, x_j^0))_{i,j=1,\dots,s} \end{pmatrix} \in \mathbb{R}^{s(1+2d) \times s(1+2d)}. $$

The following lemma entails that these linear systems admit a solution provided that a minimal separation condition holds.

**Lemma I.1.** Let {$x_j^0$}$_{j=1}^s \subset \mathbb{R}^d \times [u_{\min}, +\infty)^d$. Assume that $K_{\text{norm}}$ verifies Definition 5.2. If $\min_{i \neq j} d(x_i^0, x_j^0) \ge \Delta$, then $\Upsilon$ (see (63)) is invertible.

$$ \text{Moreover, } \begin{pmatrix} \alpha \\ \beta \end{pmatrix} = \Upsilon^{-1} \mathbf{u}_s \quad \text{and} \quad \begin{pmatrix} \alpha^j \\ \beta^j \end{pmatrix} = \Upsilon^{-1} \mathbf{u}_s^j \text{ are well-defined and verify, for all } j = 1, \dots, s, $$

$$ \begin{gathered} \|\alpha\|_{\infty} \vee \|\alpha^j\|_{\infty} \le \frac{1}{1-2h}, \\ \max_{l=1,\dots,s} \|\beta_l\|_{x_l^0} \vee \max_{l=1,\dots,s} \|\beta_l^j\|_{x_l^0} \le 4h \end{gathered} \tag{64} $$

$$ \text{and } \|\alpha - 1_s\|_{\infty} \vee |\alpha_j^j - 1| \vee \max_{l \neq j} |\alpha_l^j| \le \frac{2h}{1-2h}, $$

$$ \text{where } h = \frac{1}{64} \min \left( \frac{\bar{\varepsilon}_0(r)}{B_0}, \frac{\bar{\varepsilon}_2(r)}{B_2} \right). $$

The proof can be found in [Poon et al., 2023, pp. 269-270]. This proof explicitly requires to handle normalized kernels. This motivates the introduction of $W$ in $(\mathcal{P}_{\kappa})$.

**Controlling the certificates via the kernel** We provide controls on the certificates we constructed in Lemma I.1 on the far and near regions, in order to show that they are non-degenerate (Definitions 3.2 and 3.3). To do so, we use Taylor expansions on Fisher-Rao geodesics along with bounds on the kernel stemming from the LPC (Definition 5.2).

**Lemma I.2.** Let $r > 0$ and $x_0 \in \mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, +\infty)^d$. We define

$$ \mathcal{X}_0^{\text{near}}(r) := \{x' \in \mathcal{X} : d(x', x_0) \le r\}. $$

Assume that there exists $\bar{\varepsilon}_2 > 0$ such that, for all $x \in \mathcal{G}_{x_0}(\mathcal{X}_0^{\text{near}}(r))$ and $v \in \mathbb{R}^{2d}$,

$$ -K_{\text{norm}}^{(02)}(x_0, x)[v, v] \ge \bar{\varepsilon}_2 \|v\|_x^2 \quad \text{and} \quad \|K_{\text{norm}}^{(02)}(x_0, x)\|_x \le B_{02}. $$

Let $\eta: \mathbb{R}^d \times [u_{\min}, +\infty)^d \to \mathbb{R}$ be a $\mathcal{C}^2$ function. Then the following holds:

(i) If $\eta(x_0) = 0$, $\nabla\eta(x_0) = 0$ and $\|D_2[\eta](x)\|_x \le \delta$ for all $x \in \mathcal{G}_{x_0}(\mathcal{X}_0^{\text{near}}(r))$, then $|\eta(x)| \le \delta \mathfrak{D}_g(x, x_0)^2$ for all $x \in \mathcal{X}_0^{\text{near}}(r)$.

(ii) Let $a \in \{-1, 1\}$. If $\eta(x_0) = a$, $\nabla\eta(x_0) = 0$ and $\|aD_2[\eta](x) - K_{\text{norm}}^{(02)}(x_0, x)\|_x \le \delta$ for all $x \in \mathcal{G}_{x_0}(\mathcal{X}_0^{\text{near}}(r))$, for some $\delta < \bar{\varepsilon}_2$, then

$$ 1 - \frac{B_{02} + \delta}{2} \mathfrak{D}_g(x, x_0)^2 \le a\eta(x) \le 1 - \frac{\bar{\varepsilon}_2 - \delta}{2} \mathfrak{D}_g(x, x_0)^2 \quad \forall x \in \mathcal{X}_0^{\text{near}}(r). $$

*Proof.* This proof is based on [Poon et al., 2023, Lemma 2]. We only show (ii), as the proof for (i) is similar. Let $x \in \mathcal{X}_0^{\text{near}}(r)$. We denote by $\gamma$ the geodesic for the metric $g$ between $x_0$ and $x$, parametrized between 0 and 1. We refer to Section H.2 for a description of the geodesic properties and related notations. By definition, $\gamma(y) \in \mathcal{G}_{x_0}(\mathcal{X}_0^{\text{near}}(r))$ for all $y \in [0, 1]$. Hence, for all $y \in [0, 1]$,

$$ \|aD_2[\eta](\gamma(y)) - K_{\text{norm}}^{(02)}(x_0, \gamma(y))\|_{\gamma(y)} \le \delta \quad \text{and} \quad -K_{\text{norm}}^{(02)}(x_0, \gamma(y))[\dot{\gamma}(y), \dot{\gamma}(y)] \ge \bar{\varepsilon}_2 \|\dot{\gamma}(y)\|_{\gamma(y)}^2, $$
from which we deduce that

$$aD_2[\eta](\gamma(y))[\dot{\gamma}(y), \dot{\gamma}(y)] \le (-\bar{\epsilon}_2 + \delta) \|\dot{\gamma}(y)\|_{\gamma(y)}^2.$$

By a Taylor expansion, we get

$$
\begin{align*}
a\eta(x) &= a\eta(x_0) + a\nabla\eta(x_0)^T\dot{\gamma}(0) + \int_0^1 (1-y)aD_2[\eta](\gamma(y))[\dot{\gamma}(y), \dot{\gamma}(y)] dy, \\
&= 1 + \int_0^1 (1-y)aD_2[\eta](\gamma(y))[\dot{\gamma}(y), \dot{\gamma}(y)] dy, \\
&\le 1 - (\bar{\epsilon}_2 - \delta) \int_0^1 (1-y) \|\dot{\gamma}(y)\|_{\gamma(y)}^2 dy, \\
&= 1 - \frac{\bar{\epsilon}_2 - \delta}{2} \mathfrak{D}_g(x, x_0)^2.
\end{align*}
$$

Using the same reasoning, we also have $a\eta(x) \ge 1 - \frac{B_{02}+\delta}{2}\mathfrak{D}_g(x, x_0)^2$.

Note that we do not have (unlike in [Poon et al., 2023, Lemma 2]) the bound $a\eta(x) \ge -1 + \frac{\bar{\epsilon}_2-\delta}{2}\mathfrak{D}_g(x, x_0)^2$. This is not a problem in our framework: as we work with nonnegative measures, we do not need to control the negative part of the certificate. □

**Theorem I.1.** Let $\{x_j^0\}_{j=1}^s \subset \mathcal{X} \subset \mathbb{R}^d \times [u_{\min}, u_{\max}]^d$. Assume that $K_{\text{norm}}$ verifies Definition 5.2. If $\min_{i \neq j} d(x_i^0, x_j^0) \ge \Delta_{\tau}$ (defined in Lemma 5.3), then the certificates constructed in Lemma I.1 are non-degenerate. The global certificate $\eta_{\alpha,\beta}$ is $(\frac{7}{8}\bar{\epsilon}_0, \frac{15}{32}\bar{\epsilon}_2, r)$-non-degenerate and the local certificate $\eta_{\alpha_j,\beta j}$ is $(\frac{7}{8}\bar{\epsilon}_0, \frac{B_{02}+\bar{\epsilon}_2}{2}, r)$-non-degenerate.

*Proof.* This proof is largely based on [Poon et al., 2023, pp. 270-241]. First remark that Lemmas I.1 and 5.3 hold under these assumptions. We denote $h = \frac{1}{64} \min \left( \frac{\bar{\epsilon}_0(r)}{B_0}, \frac{\bar{\epsilon}_2(r)}{B_2} \right)$.

Control on the far region: Let $x \in \mathcal{X}^{far}(r)$. As the open balls $\tilde{B}_d(x_i^0, \Delta) \cap \mathcal{X}$ are disjoint (Lemma 5.3), there exists at most one index $j$ such that $d(x, x_j^0) < \Delta$ and for all $i \neq j$, $d(x, x_i^0) \ge \Delta$. So using bounds displayed in Lemma I.1 (see (64)),

$$
\begin{align*}
|\eta_{\alpha,\beta}(x)| &= \left| \alpha_j K_{\text{norm}}(x_j^0, x) + \sum_{j \neq i} \alpha_i K_{\text{norm}}(x_i^0, x) + \beta_j^T K_{\text{norm}}^{(10)}(x_j^0, x) + \sum_{j \neq i} \beta_i^T K_{\text{norm}}^{(10)}(x_i^0, x) \right|, \\
&\le \|\alpha\|_{\infty} \left( |K_{\text{norm}}(x_j^0, x)| + \sum_{j \neq i} |K_{\text{norm}}(x_i^0, x)| \right) + \max_i \|\beta_i\|_{x_i^0} \left( \left\| K_{\text{norm}}^{(10)}(x_j^0, x) \right\|_{x_j^0} + \sum_{j \neq i} \left\| K_{\text{norm}}^{(10)}(x_i^0, x) \right\|_{x_i^0} \right), \\
&\le \frac{1}{1-2h} (1-\bar{\epsilon}_0+h) + 4h(B_{10}+h), \\
&\le 1 - \frac{\bar{\epsilon}_0-3h}{1-2h} + 4h(B_{10}+h), \\
&\le 1 - \bar{\epsilon}_0 + 3h + 4h(B_{10}+h), \\
&\le 1 - \bar{\epsilon}_0 + 3\frac{\bar{\epsilon}_0}{64} + 4\frac{\bar{\epsilon}_0}{64} + 4\frac{\bar{\epsilon}_0}{64^2}, \\
&\le 1 - \frac{7}{8}\bar{\epsilon}_0.
\end{align*}
$$

We can apply the same reasoning to show that $|\eta_{\alpha^j, \beta^j}(x)| \le 1 - \frac{7}{8}\bar{\epsilon}_0$.

Controls on the near regions: Let $x \in \mathcal{G}_{x_j^0}(\mathcal{X}^{\text{near}}(r))$. We have

$$D_2[\eta_{\alpha,\beta}](x) = K_{\text{norm}}^{(02)}(x_j^0, x) + (\alpha_j - 1)K_{\text{norm}}^{(02)}(x_j^0, x) + \sum_{i \neq j} \alpha_i K_{\text{norm}}^{(02)}(x_i^0, x) + [\beta_j]K_{\text{norm}}^{(12)}(x_j^0, x) + \sum_{i \neq j} [\beta_i]K_{\text{norm}}^{(12)}(x_i^0, x).$$

As $\mathcal{G}_{x_j^0}(B_d(x_j^0, r) \cap \mathcal{X})$ is disjoint from $\tilde{B}_d(x_i^0, \Delta)$ for $i \neq j$ (Lemma 5.3), we have $d(x, x_i^0) \ge \Delta$ for all $i \neq j$. Using Lemma I.1, it comes

$$
\begin{align*}
\| D_2[\eta_{\alpha,\beta}](x) - K_{\text{norm}}^{(02)}(x_j^0, x) \|_{x} &\le |\alpha_j - 1| B_{02} + \|\alpha\|_{\infty} h + \max_i \|\beta_i\|_{x_i^0} (B_{12} + h), \\
&\le \frac{2h}{1-2h} B_{02} + \frac{h}{1-2h} + 4h(B_{12} + h), \\
&\le \frac{\bar{\epsilon}_2}{16}.
\end{align*}
$$