SuperMix: Sparse Regularization for Mixtures

TL;DR

This paper introduces a grid-less convex method called Beurling-LASSO for estimating discrete mixing measures in kernel mixture models, providing theoretical support for support localization and stability, with practical algorithms and performance analysis.

Contribution

It develops a new regularization-based approach for mixture estimation, deriving bounds and stability results, and proposes efficient algorithms for implementation.

Findings

Supports localization bounds depend on support separation

Estimator recovers the correct number of mixture components with large enough samples

Algorithms based on Sliding Frank-Wolfe and Conic Particle Gradient Descent are effective

Abstract

This paper investigates the statistical estimation of a discrete mixing measure $\mu$0 involved in a kernel mixture model. Using some recent advances in l1-regularization over the space of measures, we introduce a "data fitting and regularization" convex program for estimating $\mu$0 in a grid-less manner from a sample of mixture law, this method is referred to as Beurling-LASSO. Our contribution is twofold: we derive a lower bound on the bandwidth of our data fitting term depending only on the support of $\mu$0 and its so-called "minimum separation" to ensure quantitative support localization error bounds; and under a so-called "non-degenerate source condition" we derive a non-asymptotic support stability property. This latter shows that for a sufficiently large sample size n, our estimator has exactly as many weighted Dirac masses as the target $\mu$0 , converging in amplitude and localization towards the true ones. Finally, we also introduce some tractable algorithms for solving this convex program based on "Sliding Frank-Wolfe" or "Conic Particle Gradient Descent". Statistical performances of this estimator are investigated designing a so-called "dual certificate", which is appropriate to our setting. Some classical situations, as e.g. mixtures of super-smooth distributions (e.g. Gaussian distributions) or ordinary-smooth distributions (e.g. Laplace distributions), are discussed at the end of the paper.

Figures

Equations

$${\mu }^0:=\sum _{k=1}^K{a}_k^0{\delta }_{t_k}$$

$$\sum _{k=1}^K{a}_k^0=1\text{and}\forall k\in [K]:=\{1,\dots ,K\}:a_k^0>0.$$

$$X_i\stackrel{\mathrm{iid}}{\sim }\sum _{k=1}^K{a}_k^0{F}_{t_k},\forall i\in [n]:=\{1,\dots ,n\},$$

$$f^0(x)=\sum _{k=1}^K{a}_k^0\phi (x-t_k),\forall x\in \text{mathds}{R}^d.$$

$$X_i\sim U^0+E,$$

$$\Delta :=\underset{k\ne \ell }{\min }\Vert t_k-t_{\ell }{\Vert }_2.$$

$$m\gtrsim \sqrt{K}{d}^{3/2}{\Delta }_{+}^{-1}\mathrm{where}{\Delta }_{+}=\min (\Delta ,1).$$

$$\rho_{n}=\mathcal{O}\Big{(}\sqrt{\frac{m^{d}}{n}}\,\Big{)},$$

$$\forall A\subset \text{mathds}{R}^d,\text{mathds}E[{\hat{\mu }}_n(A)]\gtrsim {\rho }_n{\mathcal{C}}_m(\phi )⟹\underset{k\in [K]}{\min }\underset{t\in A}{\inf }\Vert t-t_k{\Vert }_2^2\lesssim \frac{1}{m^2}.$$

$$\forall k\in [K]:\text{mathds}E\left[|a_k^0-{\hat{\mu }}_n(\text{mathds}{N}_k(ϵ))|\right]\lesssim {\rho }_n{\mathcal{C}}_m(\phi ),$$

$${\hat{\mu }}_n=\sum _{k=1}^{\hat{K}}{\hat{a}}_k{\delta }_{{\hat{t}}_k},$$

$$\Vert \mu {\Vert }_1:={\int }_{\text{mathds}{R}^d}\mathrm{d}|\mu |\forall \mu \in \mathcal{M}(\text{mathds}{R}^d,\text{mathds}R),$$

$$\Vert \mu {\Vert }_1=\sup \left\{{\int }_{\text{mathds}{R}^d}f\mathrm{d}\mu :f\text{is}\text{μ-measurable and }|f|\le 1\right\}.$$

$$\forall x\in \text{mathds}{R}^d,\forall f\in L^1(\text{mathds}{R}^d),\mathcal{F}[f](x):={\int }_{\text{mathds}{R}^d}{e}^{-ıx^{\top }\omega }f(\omega )\mathrm{d}\omega .$$

$$\forall x\in \text{mathds}{R}^d,\forall \mu \in \mathcal{M}(\text{mathds}{R}^d,\text{mathds}R),\mathcal{F}[\mu ](x):={\int }_{\text{mathds}{R}^d}{e}^{-ıx^{\top }\omega }\mathrm{d}\mu (\omega ).$$

$$\mu \mapsto \Phi (\mu ):=\phi \star \mu ={\int }_{\text{mathds}{R}^d}\phi (\cdot -x)\mathrm{d}\mu (x),\mu \in \mathcal{M}(\text{mathds}{R}^d,\text{mathds}R),$$

$$\forall \mu \in \mathcal{M}(\text{mathds}{R}^d,\text{mathds}R),\mathcal{F}[\Phi (\mu )]=\mathcal{F}[\phi ]\mathcal{F}[\mu ].$$

$$\text{The function}\phi \text{is a bounded continuous symmetric function of positive definite type.}$$

$$\sum _{i=1}^n\sum _{j=1}^n\phi (x_i-x_j)z_i{\bar{z}}_j\ge 0.$$

$$\phi ={\mathcal{F}}^{-1}[\sigma ]\text{ for some nonnegative }\sigma \in L^1(\text{mathds}{R}^d).$$

$$\mathrm{Supp}(\sigma )=\left\{\omega \in \text{mathds}{R}^d:\sigma (\omega )\ne 0\right\}.$$

$${\mu }^0={\mathcal{F}}^{-1}\left[\mathcal{F}(f^0){\sigma }^{-1}\right],$$

$$\underset{\mu \in \mathcal{M}(\text{mathds}{R}^d,\text{mathds}R)\text{:}\Phi (\mu )=f^0}{\min }\Vert \mu {\Vert }_1.$$

$$\mathcal{M}(f^0):=\{\mu \in \mathcal{M}(\text{mathds}{R}^d,\text{mathds}R)\text{:}\Phi (\mu )=f^0\}.$$

$$\phi ={\mathcal{F}}^{-1}[\sigma ]\text{, }\sigma (\omega )=\sigma (-\omega )\text{ a.e. with }\mathrm{Supp}(\sigma )=\text{mathds}{R}^d:\forall \omega \in \text{mathds}{R}^d\sigma (\omega )>0.$$

$$\mathcal{M}(f^0)=\{{\mu }^0\}.$$

$$\phi ={\mathcal{F}}^{-1}[\sigma ]\text{, }\sigma (\omega )=\sigma (-\omega )\text{ a.e. with }[-\eta ,\eta {]}^d\subset \mathrm{Supp}(\sigma ).$$

$${\hat{f}}_n:=\frac{1}{n}\sum _{i=1}^n{\delta }_{X_i},$$

$$\hat{\mu}_{n}:=\arg\min_{\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})}\big{\{}C(\Phi\mu,\hat{f}_{n})+\kappa\|\mu\|_{1}\big{\}},$$

$$⟨f,\ell (x,.){⟩}_{\text{mathds}L}=f(x)\forall f\in \text{mathds}L,\forall x\in \text{mathds}{R}^d.$$

Full text

Our contribution is two-fold: we derive a lower bound on the bandwidth of our data fitting term depending only on the support of $\mu^{0}$ and its so-called “minimum separation” to ensure quantitative support localization error bounds; and under a so-called “non-degenerate source condition” we derive a non-asymptotic support stability property. This latter shows that for a sufficiently large sample size $n$, our estimator has exactly as many weighted Dirac masses as the target $\mu^{0}$, converging in amplitude and localization towards the true ones. Finally, we also introduce some tractable algorithms for solving this convex program based on “Sliding Frank-Wolfe” or “Conic Particle Gradient Descent”.

Statistical performances of this estimator are investigated designing a so-called “dual certificate”, which is appropriate to our setting. Some classical situations, as e.g. mixtures of super-smooth distributions (e.g. Gaussian distributions) or ordinary-smooth distributions (e.g. Laplace distributions), are discussed at the end of the paper.

62G05,

90C25,

49M29,

Beurling Lasso; Mixture recovery; Dual certificate; Kernel approach; Super-resolution,

1 Introduction

1.1 Mixture problems

In this paper, we are interested in the estimation of a mixture distribution $\mu^{0}$ using some i.i.d. observations $\mathbf{X}:=(X_{1},\ldots,X_{n})\in(\mathds{R}^{d})^{n}$ with the help of some $\ell_{1}$-regularization methods. More precisely, we consider the specific situation of a discrete distribution $\mu^{0}$ that is given by a finite sum of $K$ components:

[TABLE]

where the set of positive weights $(a_{k}^{0})_{1\leq k\leq K}$ defines a discrete probability distribution, i.e. each $\delta_{t_{k}}$ is a Dirac mass at point $t_{k}\in\mathds{R}^{d}$ while

[TABLE]

We denote by $S^{0}:=\{t_{1},\ldots,t_{K}\}$ the support of the target distribution $\mu^{0}$. This distribution is indirectly observed: we assume that our set of observations $\mathbf{X}$ in $\mathds{R}^{d}$ satisfies

[TABLE]

where $(F_{t})_{t\in\mathds{R}^{d}}$ is a family of known distributions on $\mathds{R}^{d}$. Below, we consider the so-called location model where each distribution $F_{t}$ has a density with respect to the Lebesgue measure on $\mathds{R}^{d}$ given by the density function $\varphi(\cdot-t)$, where $\varphi$ denotes a known density function. In this case, the density function $f^{0}$ of the data $\mathbf{X}$ can be written as a convolution, namely

[TABLE]

Remark 1**.**

Equation (2) has a simple interpretation in the context considered here: the law of one observation $X_{i}$ is given by a sum of two independent random variables $U^{0}$ and $E$:

[TABLE]

where $U^{0}\in S^{0}$ is distributed according to $\mu^{0}$ $($i.e., the mixing law (1)$)$ and $E$ has a distribution of density $\varphi$ with respect to the Lebesgue measure on $\mathds{R}^{d}$. In this context, recovering the distribution of $U^{0}$ from the sample $\mathbf{X}$ appears to be an inverse (deconvolution) problem. The main difference with former contributions (see, e.g. [28] for a comprehensive introduction) is that the probability measure associated to $U^{0}$ is discrete, which avoids classical regularization approaches.

Equation (2) is known in the literature as a mixture model. A mixture model allows to describe some practical situations where a population of interest is composed of $K$ different sub-populations, each of them being associated to a proportion $a_{k}^{0}$ and to a location parameter $t_{k}$. Mixture models have been intensively investigated during the last decades and have been involved in several fields as biology, genetics, astronomy, among others. We refer to [20, 27] for a complete overview.

1.2 Previous works

The main goal of this paper is to provide an estimation of the discrete mixture law $\mu^{0}$ introduced in (1). When the component number $K$ is available, the maximum likelihood estimator (MLE) appears to be the most natural candidate. Although no analytic expression is available for the model (2), it can be numerically approximated. We mention for instance the well-known EM-algorithm and refer to [40], who established some of the most general convergence results known for the EM algorithm. However, the MLE (and the related EM-algorithm) does not always provide satisfactory results. First, the MLE suffers from several drawbacks (see, e.g., [25]) such as non-uniqueness of the solution, and second, obtaining theoretical guarantees for the EM-algorithm is still a difficult question (see, e.g., the recent contributions [3, 17]). Several alternative methods have been proposed in this context. Some contributions extensively use the MLE point of view to derive consistent properties in general semi-parametric models, including the Gaussian case (see e.g. [38]), whereas some other ones developed some contrast functions in a semi-parametric framework: with symmetry and number of component assumptions in [6, 8], or with a fixed number of component settings in [21] and a $L^{2}$ contrast. As a particular case, the Gaussian setting has attracted a lot of attention: a model selection strategy is developed in [26] and a specific analysis of the EM algorithm with two Gaussian components is provided in [41]. The article [3] provides a general theoretical framework to analyze the convergence of the EM updates in a neighborhood of the MLE, and derives some non-asymptotic bounds on the Euclidean error of sample-based EM iterates. Some of the aforementioned papers provide better results (for instance with parametric rates of convergence for the estimation of the weights $a^{0}_{k}$, see e.g. [29, 22]), but are obtained in more constrained settings: known fixed number of components (often $K=2$), univariate case, …

Our estimator will be any solution to a convex program and it does not require to know the number $K$ of components in the mixing law $\mu^{0}$. This estimator is based on ideas from super-resolution and “off-the-grid” methods [5, 9], where one aims at recovering a discrete measure from linear measurements. The so-called “sparse deconvolution” problem fits this framework since it concerns with estimating a target measure from the observation of some product of convolution between the target measure and known kernel as $f^{0}$ in (2). Note that in mixture models, we do not observe $f^{0}$ but rather a sample drawn from it, and standard strategies such that $(1.15)$ in [9] cannot be invoked here. However, remark that one of the main advances has been the construction of the so-called “dual certificate” in [9] which is the key to demonstrate the success of discreteness inducing norm regularization (see e.g. [13, 9, 16, 14]).

Recent works have addressed mixture models while assuming that the sampling law is known. For example, the authors of [31] study some dimension reduction techniques such as random “sketching” problems using “off-the-grid” minimization scheme. They prove convergence of random feature kernel towards the population kernel. We emphasize that the statistical estimation in terms of the sample size $n$ has not been considered in the super-resolution research field. To the best of our knowledge, this paper is the first that bridges the gap between the recent “off-the-grid” sparse regularization methods and a sharp statistical study of this estimation procedure in terms of the sample size and the bandwidth of the data fitting term.

1.3 Contribution

In this paper, we propose an estimator $\hat{\mu}_{n}$ of the measure $\mu^{0}$ (see Equation (1)) inspired by some recent results in $\ell_{1}$-regularization on the space of measures, sometimes referred to as super-resolution methods (see, e.g., [13, 9]). We investigate the statistical theoretical performances of $\hat{\mu}_{n}$. This estimator $\hat{\mu}_{n}$ is built according to the minimization of a criterion on the space of real measures on $\mathds{R}^{d}$ and does not require any grid for its computation. The stability result and the construction of the dual certificate given in [9] played a central role in our work to obtain the statistical recovery. However, these authors work on the torus and their construction provides periodic dual certificates which are not useful in our present framework. One important contribution of this paper is thus a novel dual certificate construction, interpolating phases/signs on $\mathds{R}^{d}$ (and not the $d$-dimensional torus as in [9]). We also investigate the stability with respect to sampling of our estimation strategy, i.e. the ability of our procedure to recover the mixture when we compute $\hat{\mu}_{n}$ up to some i.i.d. observations $(X_{i})_{i\in[n]}$ with $n\to+\infty$, which is a different problem from the stability issue studied in [9] that asserts the variation of the super-resolution solutions with respect to an $\ell_{1}$ norm control on the low-frequency data.

The minimized criterion requires to tune two parameters: a bandwidth parameter of the data fitting term denoted by $m\geq 1$ and an $\ell_{1}$-regularization tuning parameter denoted by $\kappa>0$ below. We prove that the bandwidth parameter $m$ depends only on the intrinsic hardness of estimating the support $S^{0}$ of the target $\mu^{0}$ through the so-called “minimum separation” $\Delta$ introduced in [9] that refers to the minimal distance between two spikes:

[TABLE]

We now assess briefly the performances of $\hat{\mu}_{n}$. We emphasize that a complete version is displayed in Theorem 10 (for points $i)$ and $ii)$) and Theorem 11 (for point $iii)$) later on.

Theorem 1**.**

Assume that the kernel $\varphi$ satisfies ($\mathcal{H}_{\eta}$) with $\eta=4m$ (see Section 2.3 for a definition) for a bandwidth $m$ verifying

[TABLE]

Then, some quantity $\mathcal{C}_{m}(\varphi)>0$ exists such that, setting

[TABLE]

our estimator $\hat{\mu}_{n}$ satisfies:

• $i)$

Spike detection property:

[TABLE]

• $ii)$

Weight reconstruction property:

[TABLE]

where $\mathds{N}_{k}(\epsilon)$ denotes a region that contains $t_{k}$ and $\epsilon=\epsilon_{n,m}(d)$ is made explicit later on.

• $iii)$

Support stability property: if $\varphi$ satisfies the Non-Degenerated Bandwidth condition (NDB) (see Section 4.4 for a definition), for $n$ large enough, with an overwhelming probability, $\hat{\mu}_{n}$ can be written as

[TABLE]

with $\hat{K}=K$. Furthermore, $(\hat{a}_{k},\hat{t}_{k})\to(a^{0}_{k},t_{k})$ for all $k\in[K]$, as $n$ tends to infinity.

Note that the constant $\mathcal{C}_{m}(\varphi)$ will depend on other quantities introduced later. It will be specified in Proposition 9.

These three results deserve several comments. First, $i)$ indicates that when a set $A$ has enough mass w.r.t. the estimated measure $\hat{\mu}_{n}$, it includes a true spike with an accuracy of the order $m^{-2}$. The second result $ii)$ provides some statistical guarantees on the mass set by $\hat{\mu}_{n}$ near a true spike $t_{k}$ that converges to $\mu^{0}(\{t_{k}\})=a_{k}^{0}$. Condition (NDB) is inspired from the so-called “non-degenerated source condition” (NDSC) introduced in [16] and allows to derive the support stability. The last result $iii)$ shows that, for large enough sample size, $\ell_{1}$-regularization successfully recovers the number of mixing components. The estimated weights on the Dirac masses then converge towards the true ones in amplitudes and localizations.

The bandwidth $m$ has to be adjusted to avoid over and under-fitting. Condition (3) ensures that the target point is admissible for our convex program and it may be seen as a condition to avoid a large bias term and under-fitting. Condition (4) ensures that the sample size is sufficiently large with respect to the model size $m$ and it might be seen as a condition to avoid over-fitting and therefore to upper-bound the variance of estimation.

Below, we will pay attention to the role of Fourier analysis of $\varphi$ and to the dimension $d$ of the ambient space. These results are applied to specific settings (super-smooth and ordinary-smooth mixtures).

1.4 Outline

This paper is organized as follows. Section 2 introduces some standard ingredients of $\ell_{1}$ regularization methods and gives a deterministic analysis of the exact recovery property of $\mu^{0}$ from $f^{0}$. Section 3 provides a description of the statistical estimator $\hat{\mu}_{n}$ derived from a deconvolution with a Beurling-LASSO strategy (BLASSO) (see e.g. [13]). Tractable algorithms solving BLASSO when the observation is a sample from a mixing law are introduced in Section 3.3. Section 4 focuses on the statistical performances of our estimator whereas Section 5 details the rates of convergence for specific mixture models. The main proofs are gathered in Section 6 whereas the most technical ones are deferred to the appendix.

2 Assumptions, notation and first results

This section gathers the main assumptions on the mixture model (2). Preliminary theoretical results in an “ideal” setting are stated in order to ease the understanding of the forthcoming paragraphs.

2.1 Functional framework

We introduce some notation used all along the paper.

Definition 1** (Set $(\mathcal{M}(\mathds{R}^{d},\mathds{R}),\|\cdot\|_{1})$).**

We denote by $(\mathcal{M}(\mathds{R}^{d},\mathds{R}),\|\cdot\|_{1})$ the space of real valued measures on $\mathds{R}^{d}$ equipped with the total variation norm $\|\cdot\|_{1}$, which is defined as

[TABLE]

where $|\mu|=\mu^{+}+\mu^{-}$ and $\mu=\mu^{+}-\mu^{-}$ is the Jordan decomposition associated to a given measure $\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$.

A standard argument proves that the total variation of $\mu$ is also described with the help of a variational relationship:

[TABLE]

Recall that $\varphi$ used in Equation (2) is a probability density function so that $\varphi\in L^{1}(\mathds{R}^{d})$.

Definition 2** (Fourier transform over $L^{1}(\mathds{R}^{d})$ and $\mathcal{M}(\mathds{R}^{d},\mathds{R})$).**

We denote by $\mathcal{F}$ the Fourier transform defined by:

[TABLE]

A standard approximation argument extends the Fourier transform to $\mathcal{M}(\mathds{R}^{d},\mathds{R})$ with:

[TABLE]

We denote by $\mathcal{C}_{0}(\mathds{R}^{d},\mathds{R})$ the space of continuous real valued functions vanishing at infinity on $\mathds{R}^{d}$ and recall that $\mathcal{F}\left(L^{1}(\mathds{R}^{d})\right)$ is a dense subset of $\mathcal{C}_{0}(\mathds{R}^{d},\mathds{R})$. We shall also introduce the convolution operator $\Phi$ as

[TABLE]

and it holds equivalently that (see e.g. [33, Section 9.14]):

[TABLE]

Concerning the density $\varphi$ involved in (2), we will do the following assumption.

[TABLE]

In particular, the positive definite type property involved in Assumption ($\mathcal{H}_{0}$) is equivalent to require that for any finite set of points $\{x_{1},\ldots,x_{n}\}\in\mathds{R}^{d}$ and for any $(z_{1},\ldots,z_{n})\in\mathds{C}^{n}$:

[TABLE]

In what follows, we consider $h:\mathds{R}^{d}\times\mathds{R}^{d}\longrightarrow\mathds{R}$ the function defined by $h(x,y)=\varphi(x-y)$ for all $x,y\in\mathds{R}^{d}$. In such a case, Assumption ($\mathcal{H}_{0}$) entails that $h(\cdot,\cdot)$ is a bounded continuous symmetric positive definite kernel. By Bochner’s theorem (see, e.g., [33, Theorem 11.32]), $\varphi$ is the inverse Fourier transform of a nonnegative measure $\Sigma$ referred to as the spectral measure. The Fourier inversion theorem states that $\Sigma$ has a nonnegative density $\sigma\geq 0$ with respect to the Lebesgue measure on $\mathds{R}^{d}$ such that $\sigma\in L^{1}(\mathds{R}^{d})$. Hence, it holds from the preceding discussion that

[TABLE]

Below, the set of points where the Fourier transform of a function does not vanish will play an important role. We will denote this support by $\mathrm{Supp}(\sigma)$:

[TABLE]

Some examples of densities $\varphi$ that satisfies ($\mathcal{H}_{0}$) will be given and discussed in the forthcoming sections. We emphasize that this assumption is not restrictive and concerns for instance Gaussian, Laplace or Cauchy distributions, this list being not exhaustive.

Additional notation. Given two real sequences $(a_{n})_{n\in\mathds{N}}$ and $(b_{n})_{n\in\mathds{N}}$, we write $a_{n}\lesssim b_{n}$ (resp. $a_{n}\gtrsim b_{n}$) if there exists a constant $C>0$ independent of $n$ such that $a_{n}\leq b_{n}$ (resp. $a_{n}\geq b_{n}$) for all $n\in\mathds{N}$. Similarly, we write $a_{n}\ll b_{n}$ if $a_{n}/b_{n}\rightarrow 0$ as $n\rightarrow+\infty$. The set $\mathds{N}^{*}$ stands for $\mathds{N}\setminus\{0\}$.

2.2 Exact Recovery of $\mu^{0}$ from $f^{0}$ - Case $\mathrm{Supp}(\sigma)=\mathds{R}^{d}$

In this paragraph, we are interested in an “ideal” problem where we are looking for $\mu^{0}$ not from a sample $X_{1},\dots,X_{n}$ distributed according to Equation (2), but from the population law $f^{0}$ itself. Of course, this situation does not occur in practice since in concrete situations, we do not observe $f^{0}$ but an empirical version of it and we will have to preliminary use an estimation of $f^{0}$ before solving the deconvolution inverse problem. Nevertheless, this toy problem already provides the first ingredients for a better understanding of the difficulties that arise in the context we consider.

We stress that $f^{0}:=\Phi(\mu^{0})$ where $\Phi$ is defined by (5). Hence, this paragraph concerns the recovery of $\mu^{0}$ from its convolution by the kernel $\varphi$. We thus face an inverse (deconvolution) problem. Several solutions could be provided and a standard method would rely on Fourier inversion

[TABLE]

where $\sigma$ is given by (7).

Here, we prove in a first step that this deconvolution problem can be efficiently solved using a $\ell_{1}$-regularization approach. We will be interested in the convex program (8) given by:

[TABLE]

In particular, we investigate under which conditions the solution set of (8) is the singleton $\{\mu^{0}\}$, that we referred to as the “Perfect Recovery” property. We introduce the set of admissible points to the program (8), denoted by $\mathcal{M}(f^{0})$ and defined as:

[TABLE]

In this context, some different assumptions on the kernel $\varphi$ shall be used in our forthcoming results.

A first reasonable situation is when the spectral density $\sigma=\mathcal{F}(\varphi)$ has its support equal to $\mathds{R}^{d}$ and in this case we denote $\sigma>0$. This requirement can be summarized in the next assumption on the function $\varphi$:

[TABLE]

Example 1**.**

It may be shown that the set of densities $\varphi$ that satisfy both Assumptions ($\mathcal{H}_{0}$) and ($\mathcal{H}_{\infty}$) include the Gaussian, Laplace, $B_{2\ell+1}$-spline, inverse multi-quadrics, Matérn class (see, e.g., [36, top of page 2397]) examples.

Under Assumptions ($\mathcal{H}_{0}$) and ($\mathcal{H}_{\infty}$), any target measure $\mu^{0}\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$ is the only admissible point of the program (8).

Theorem 2** (Perfect Recovery under ($\mathcal{H}_{0}$) and ($\mathcal{H}_{\infty}$)).**

Assume that the convolution kernel satisfies ($\mathcal{H}_{0}$) and ($\mathcal{H}_{\infty}$), then for any target $\mu^{0}$ the program (8) has $\mu^{0}$ as unique solution point:

[TABLE]

We emphasize that the previous result also holds for measures $\mu^{0}$ that are not necessarily discrete. The proof is given in Appendix C.1.

2.3 The Super-resolution phenomenon

Theorem 2 entails that the measure $\mu^{0}$ can be recovered as soon as the spectrum of $f^{0}$ is observed and as soon as its support is $\mathds{R}^{d}$. Surprisingly, this latter assumption can be relaxed and reconstruction can be obtained in some specific situations. Such a phenomenon is associated to the super-resolution theory and has been popularized by [9] among others.

Of course, this reconstruction is feasible at the expense of an assumption on the Fourier transform of $\varphi$. For the sake of simplicity, we assume that the spectral density $\sigma$ has a support that contains the hypercube $[-\eta,\eta]^{d}$ for some frequency threshold $\eta>0$:

[TABLE]

Remark 2**.**

The densities $\varphi$ that satisfy ($\mathcal{H}_{\eta}$) and for which $\mathrm{Supp}(\sigma)=[-\eta,\eta]^{d}$ act as “ low pass filters”. The convolution operator $\Phi$ described in (5) cancels all frequencies above $\eta$, see for instance (6). Of course, the larger $\eta$, the easier the inverse deconvolution problem.

Under ($\mathcal{H}_{0}$) and ($\mathcal{H}_{\eta}$), the target measure $\mu^{0}\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$ is not the only admissible point in $\mathcal{M}(f^{0})$ to the program (8). We will need to ensure the existence of a specific function, called in what follows a dual certificate, that will entail that $\mu^{0}$ is still the only solution of the program (8).

Theorem 3** (Dual Certificate for (8)).**

Assume that the density $\varphi$ satisfies ($\mathcal{H}_{0}$) and ($\mathcal{H}_{\eta}$) for some $\eta>0$. Assume that $\mu^{0}$ and $S^{0}=\{t_{1},\ldots,t_{K}\}$ are given by Equation (1) and that a function $\mathcal{P}_{\eta}$ exists such that it satisfies the interpolation conditions:

• •

$\forall t\in\{t_{1},\ldots,t_{K}\}\,:\mathcal{P}_{\eta}(t)=1\quad\text{and}\quad\forall t\notin\{t_{1},\ldots,t_{K}\}\,:|\mathcal{P}_{\eta}(t)|<1,$**

and the smoothness conditions:

• •

$\mathcal{P}_{\eta}\in\mathcal{C}_{0}(\mathds{R}^{d},\mathds{R})\cap L^{1}(\mathds{R}^{d})$,

• •

the support of the Fourier transform $\mathcal{F}[\mathcal{P}_{\eta}]$ satisfies $\mathrm{Supp}\left(\mathcal{F}[\mathcal{P}_{\eta}]\right)\subset[-\eta,\eta]^{d}$.

Then the program (8) has $\mu^{0}$ as unique solution point $($Perfect Recovery$)$.

The proof is given in Appendix C.2. A construction of such a certificate $\mathcal{P}_{\eta}$ is presented in Appendix Fwith some additional constraints. In particular, it will make it possible to address the more realistic statistical problem where only an empirical measure of the data is available.

Remark 3**.**

The previous theorem can be extended to the case where the convolution kernel is bounded, continuous and symmetric positive definite. The proof is the same substituting $[-\eta,\eta]^{d}$ by the support $\Omega$ of its spectral density. Remark that since $\sigma$ is nonzero, necessarily $\Omega$ has a nonempty interior.

3 Off-The-Grid estimation via the Beurling-LASSO (BLASSO)

In this section, we consider the statistical situation where the density $f^{0}$ is not available and we deal instead with a sample $\mathbf{X}=(X_{1},\dots,X_{n})$ of i.i.d. observations distributed with the density $f^{0}$. In this context, only the empirical measure

[TABLE]

is available, and our aim is to recover $\mu^{0}$ from $\hat{f}_{n}$. To this end, we use in this paper a BLASSO procedure (see e.g. [2]). Namely we deal with the following estimator $\hat{\mu}_{n}$ of the unknown discrete measure $\mu^{0}$ defined as:

[TABLE]

where $\kappa$ is a regularization parameter whose value will be made precise later on, and $C(\Phi\mu,\hat{f}_{n})$ is a data fidelity term that depends on the sample $\mathbf{X}$. The purpose of the data fidelity term is to measure the distance between the target $\mu^{0}$ and any candidate $\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$.

Some examples of possible cost functions $C:\mathds{H}\times\mathcal{M}(\mathds{R}^{d},\mathds{R})\to\mathds{R}$ are discussed in Section 3.1. Our goal is then to derive some theoretical results associated to this estimation procedure.

3.1 Kernel approach

3.1.1 RKHS functional structure

In order to design the data fidelity term, we need to define a space where we can compare the observations $\mathbf{X}=(X_{1},\dots,X_{n})$ and any model $f=\varphi\star\mu=\Phi\mu$ for $\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$. In this work, we focus our attention on a kernel approach.

Reminders on RKHS

The difficulty lies in the fact that the empirical law $\hat{f}_{n}$ introduced in (9) does not belong to $\mathcal{C}_{0}(\mathds{R}^{d},\mathds{R})$. To compare the prediction $\Phi\mu$ with $\hat{f}_{n}$, we need to embed these quantities in the same space. We consider here a Reproducing Kernel Hilbert Space (RKHS) structure, which provides a lot of interesting properties and has been at the core of several investigations and applications in approximation theory [39], as well as in the statistical and machine learning communities, (see [36] and the references therein). We briefly recall the definition of such a space.

Definition 3**.**

Let $(\mathds{L},\|.\|_{\mathds{L}})$ be a Hilbert space containing function from $\mathds{R}^{d}$ to $\mathds{R}$. The space $\mathds{L}$ is said to be a RKHS if $\delta_{x}:f\mapsto f(x)$ are continuous for all $x\in\mathds{R}^{d}$ from $(\mathds{L},\|.\|_{\mathds{L}})$ to $(\mathds{R},|.|)$.

The Riesz theorem leads to the existence of a function $\ell$ that satisfies the representation property:

[TABLE]

The function $\ell$ is called the reproducing kernel associated to $\mathds{L}$. Below, we consider a kernel $\ell$ such that $\ell(x,y)=\lambda(x-y)$ for all $x,y\in\mathds{R}^{d}$ where $\lambda$ satisfies ($\mathcal{H}_{0}$). Again, the Bochner theorem yields the existence of a nonnegative measure $\Lambda\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$ such that $\lambda$ is its inverse Fourier transform

[TABLE]

Moreover, since $\lambda$ is continuous, $\Lambda$ is then a bounded measure and the Mercer theorem (see e.g. [4]) proves that the RKHS $\mathds{L}$ is exactly characterized by

[TABLE]

with dot product

[TABLE]

Convolution in the RKHS

The RKHS structure associated to the kernel $\lambda$ entails a comparison between the empirical measure and any candidate $\Phi\mu$. Indeed, a convolution operator $L$ similar to the one defined in Equation (5) can be associated to the RKHS as pointed out by the next result.

Proposition 4**.**

For any $\nu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$, the convolution $L\nu=\lambda\star\nu$ belongs to $\mathds{L}$.

The proof of Proposition 4 is given in Appendix B.1.

3.1.2 Data fidelity term

For any $\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$, both $L\hat{f}_{n}$ and $L\circ\Phi\mu$ belong to $\mathds{L}$. Hence, one may use the following data fidelity term

[TABLE]

Example 2**.**

An important example is given by the sinus-cardinal kernel $\mathrm{sinc}$. Given a frequency “cut-off” $1/\tau>0$, one can consider the kernel

[TABLE]

Then, the spectral measure is given by

[TABLE]

In this particular case, we deduce that the convolution $L$ is a low-pass filter with a frequency cut-off $1/\tau$ and the RKHS (denoted by $\mathds{L}_{\tau}$) is given by:

[TABLE]

where $B_{\infty}(1/\tau)$ denotes the centered $\ell_{\infty}$ ball of radius $1/\tau$. The RKHS $\mathds{L}_{\tau}$ then corresponds to the band-limited functions in $L^{2}(\mathds{R}^{d})$ whose Fourier transform vanishes for a frequency larger than $1/\tau$. In this context, our criterion becomes

[TABLE]

and it may be checked that

[TABLE]

This loss focuses on the $L^{2}$-error of $\Phi\mu-\hat{f}_{n}$ for frequencies in the Fourier domain $[-1/\tau,1/\tau]^{d}$. In some sense, the kernel estimator $\lambda_{\tau}\star\hat{f}_{n}$ has a bandwidth $\tau$ that will prevent from over-fitting.

We stress that, as it is the case in the previous low-pass filter example, $C_{\lambda_{\tau}}(\Phi\mu,\hat{f}_{n})$ may depend on a tuning parameter (the bandwidth $\tau$ in Example 2). For the ease of presentation, this parameter is not taken into account in the notation. However, its value will be discussed in Section 5.

3.1.3 Data-dependent computation

The next proposition entails that the criterion $C_{\lambda}$ introduced in Equation (13) can be used in practice giving a useful expression to compute it.

Proposition 5**.**

For all $\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$, we have:

[TABLE]

We stress that $\|L\hat{f}_{n}\|_{\mathds{L}}^{2}$ does no depend on $\mu$ and can be removed from the criterion when it is used in the program (8). The proof of Proposition 5 is given in Appendix B.2.

3.2 Estimation by convex programming

Our estimator is defined as a solution of the following optimization program with the data-fidelity term $C_{\lambda}(\Phi\mu,\hat{f}_{n})$ introduced in (13). Hence, we consider the optimization problem:

[TABLE]

where $\|\cdot\|_{\mathds{L}}$ is the norm associated to the RKHS generated by $\lambda$ (see Section 3.1) and $\kappa>0$ is a tuning parameter whose value will be made precise later on. We emphasize that ($\mathbf{P}_{\kappa}$) is a convex programming optimization problem (convex function to be minimized on a convex constrained set). The estimator $\hat{\mu}_{n}$ is then any solution of

[TABLE]

Algorithmic issues related to the computation of (15) are sketched in Section 3.3 and discussed in depth in Appendix A.

Remark 4**.**

The tuning parameter $\kappa>0$ needs to be chosen carefully. First note that it may depend on the choice of the frequency cut-off $1/\tau$ in $\lambda_{\tau}$, which is the bandwidth feature map (see Remark 2 for a definition). Our analysis shows that $\tau=1/{4m}$ as in (28), and $m$ is a standard nonparametric bandwidth in mixture models for which rates are given in Section 5. The main message being that it depends only on the regularity of $\varphi$ and on the sample size $n$ for $n$ large enough. From a practical view point, it is not excluded to use a Cross-Validation scheme as it heuristically performs well for $L^{1}$-based methods such as LASSO. In this case, the score function can be chosen to be the data fidelity term $\|L\hat{f}_{n}-L\circ\Phi\mu_{\mathrm{cv}}\|_{\mathds{L}}^{2}$ evaluated on the validation set. From a theoretical view point, one may choose $\kappa$ as in Equation (24). Then, Equation (25) results in

[TABLE]

and these quantities depend only on the regularity of $\varphi$ and the sample size $n$ for $n$ large enough, as studied in Section 5.

Super-resolution is the ability to recover a discrete measure on the torus from some Fourier coefficients (recall that the Pontryagin’s dual of the torus is $\mathds{Z}^{d}$) while we want to recover a discrete measure on $\mathds{R}^{d}$ from some Fourier transform over $\mathds{R}^{d}$ (recall that the Pontryagin’s dual of $\mathds{R}^{d}$ is $\mathds{R}^{d}$). In particular the dual of ($\mathbf{P}_{\kappa}$) does not involve a set of fixed degree trigonometric polynomials as in super-resolution but inverse Fourier transform of some tempered distribution.

Hence, new theoretical guarantees are necessary in order to properly define the estimator $\hat{\mu}_{n}$. This is the aim of the next theorem. In this view, we consider primal variables $\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$ and $z\in\mathds{L}$ and introduce the dual variable $c\in\mathds{L}$ as well as the following Lagrangian:

[TABLE]

It is immediate to check that if $z\neq L\circ\Phi\mu$, then the supremum of $\mathcal{L}(\mu,z,c)$ over $c$ is $+\infty$. Therefore, the primal expression coincides with the supremum in the dual variables, namely

[TABLE]

In the meantime, the dual program of ($\mathbf{P}_{\kappa}$) is given by

[TABLE]

Theorem 6** (Primal-Dual programs, strong duality).**

The following statements are true.

• $i)$

The primal problem ($\mathbf{P}_{\kappa}$) has at least one solution and it holds that

[TABLE]

hence, they do not depend on the choice of the solution $\hat{\mu}_{n}$.

• $ii)$

The dual program of ($\mathbf{P}_{\kappa}$), given by ($\mathbf{P}^{*}_{\kappa}$) satisfies

[TABLE]

and there is no duality gap $($strong duality holds$)$. Furthermore, it has a unique solution

[TABLE]

• $iii)$

Any solution $\hat{\mu}_{n}$ to ($\mathbf{P}_{\kappa}$) satisfies

[TABLE]

where

[TABLE]

i.e. it is a sub-gradient of the total variation norm at point $\hat{\mu}_{n}$.

• $iv)$

If $d=1$ and if at least one of the spectral measures $\Lambda$ or $\sigma$ has a bounded support, then $\{x\in\mathds{R}\,:\,|\hat{\eta}_{n}|(x)=1\}$ is discrete with no accumulation point, any primal solution $\hat{\mu}_{n}$ has an $($at most countable$)$ discrete support $\hat{S}\subset\mathds{R}$ with no accumulation point:

[TABLE]

The proof of this result can be found in Appendix D.

It is generally numerically admitted, see for instance [9, Page 939], that the extrema of the dual polynomial $\hat{\eta}_{n}=\Phi\hat{c}_{n}$ are located in a discrete set, so that any solution to ($\mathbf{P}_{\kappa}$) has a discrete support by using $iii)$. However, this issue remains an open question. In practice, all solvers of ($\mathbf{P}_{\kappa}$) lead to discrete solutions: greedy methods are discrete by construction, and $L^{1}$-regularization methods empirically lead to discrete solutions, see e.g. [9]. Furthermore, as presented in Theorem 11, our theoretical result shows that for large enough $n$ and under the so-called (NDB) condition, the support stability property holds. In this case, the solution of ($\mathbf{P}_{\kappa}$) is discrete with $\hat{K}=K$ atoms.

Example 3**.**

Observe that the low-pass filter defined in Example 2 satisfies the requirements of $iv)$ in Theorem 6: we deduce that when $d=1$, all solutions $\hat{\mu}_{n}$ are of the form (18).

3.3 Tractable Algorithms for BLASSO Mixture Models

Available algorithms for solving (15) with “*​off-the-grid *” methodology can be roughly divided into two categories: greedy methods and Riemannian descent methods. We emphasize that if the BLASSO has been studied in the past decade, the formulation (15) has two new important features. First the observation is a sample from a mixing law. Second, the data fidelity term has been tuned to incorporate a low pass filter kernel $\lambda$. For both methods, we refer to Appendix A for further details and references.

Greedy method: the Sliding Frank-Wolfe algorithm (SFW)

The Frank-Wolfe algorithm is an interesting avenue for solving differentiable convex programs on weakly compact convex sets, see [15] and references therein for further details, which can be adapted to compute approximate solutions of the BLASSO Mixture Models (15) with a supplementary sliding step. For a measure $\mu_{a,t}$ that may be decomposed into a finite sum of Dirac masses, we define $F_{N}$ the data-fitting term:

[TABLE]

The SFW method is then described in Algorithm 1. It is a greedy method that recursively builds

[TABLE]

see Line 3 of Algorithm 1.

Conic Particle Gradient Descent (CPGD)

Conic Particle Gradient Descent [11] is an alternative promising avenue for solving BLASSO for Mixture Models (15). The idea is still to discretize a positive measure into a system of particles, i.e. a sum of $N$ Dirac masses following (47) with $a_{i}=r_{i}^{2}$ and use a mean-field approximation in the Wasserstein space jointly associated with a Riemannian gradient descent with the conic metric. We refer to [11] and the references therein for further details. This method may be shown to be rapid, with a $\log(\epsilon^{-1})$ cost instead of $\epsilon^{-1/2}$ for standard convex programs. Adapted to the BLASSO for Mixture Models, we derive in Algorithm 2 a version of the Conic Particle Gradient Descent of [11, Algorithm 1] and we implemented this algorithm for Mixture Models in Figure 1.

More precisely, Figure 1 is a proof-of-concept and CPGD for Mixture Models would be investigated in future work. One may see that this method uncovers the right number of targets Dirac masses and their locations as some particules cluster around three poles. Some of particules vanishes and do not detect the support. Notice that a soft-thresholding effect tends to zero the small amplitudes as it may standardly be shown in $L^{1}$ regularization.

4 Statistical recovery of $\mu^{0}$

This section provides some theoretical results for $\hat{\mu}_{n}$, built as the solution of ($\mathbf{P}_{\kappa}$). Contrary to $\ell_{1}$-regularization in high-dimensions, standard RIP or REC compatibility conditions do not hold in our context, and all the cornerstone results of high-dimensional statistics cannot be used here. In our situation, the statistical analysis is performed using a “dual certificate” $\mathcal{P}_{m}$ as in super-resolution, see [2, 5, 9, 16] for instance. The construction and the key properties satisfied by $\mathcal{P}_{m}$ are detailed in Section 4.1. However, our framework is quite different from super-resolution and we had to address two issues: build a dual certificate on $\mathds{R}^{d}$ and adapt its “frequency cut-off” (namely $4m$ in $iii)$ of Theorem 7) to the sample size $n$ and the tail of the kernel. This latter point is addressed in Section 5.

4.1 Strong dual certificate

Let $S^{0}=\{t_{1},\ldots,t_{K}\}$ be a fixed set of points in $\mathds{R}^{d}$ and define $\Delta:=\min_{k\neq\ell}\|t_{k}-t_{\ell}\|_{2}.$ For any $m\in\mathds{N}^{*}$, we consider the function $p_{m}^{\alpha,\beta}$ parameterized by a vector $\alpha$ and a matrix $\beta$ of coefficients, defined as:

[TABLE]

where $\alpha=(\alpha_{1},\ldots,\alpha_{K})^{T}$, $\beta=(\beta_{k}^{i})_{1\leq k\leq K,1\leq i\leq d}$ with

[TABLE]

One important feature of $\psi_{m}$ is its ability to interpolate $1$ at the origin, while being positive and compactly supported in the Fourier domain. We then state the next result, which is of primary importance for the statistical accuracy of our procedure.

Theorem 7** (Strong dual certificate).**

Let be given a set of $K$ points $S^{0}=\{t_{1},\ldots,t_{K}\}$ in $\mathds{R}^{d}$ with $\Delta:=\min_{k\neq\ell}\|t_{k}-t_{\ell}\|_{2}$ and $\Delta_{+}=\min(\Delta,1)$. Then, the following properties hold:

• •

$i)$* A function $\mathcal{P}_{m}$ defined by $\mathcal{P}_{m}(t)=[p_{m}^{\alpha,\beta}(t)]^{2}$ exists with $\displaystyle m\gtrsim\sqrt{K}d^{3/2}\Delta_{+}^{-1}$ such that*

[TABLE]

and

[TABLE]

• •

$ii)$* A universal pair $(\upsilon,\gamma)$ independent from $n,m$ and $d$ exists such that for

$\epsilon=\frac{\upsilon}{md}:$*

• –

Near region:* If we define*

[TABLE]

a positive constant $\mathcal{C}$ exists such that:

[TABLE]

• –

Far region:**

[TABLE]

• •

$iii)$* The support of the Fourier transform of $\mathcal{P}_{m}$ is growing linearly with $m$:*

[TABLE]

• •

$iv)$* If ($\mathcal{H}_{\eta}$) holds with $\eta=4m$, then an element $c_{0,m}\in\mathds{L}$ exists such that $\mathcal{P}_{m}=\Phi c_{0,m}$.*

The proof of this result is proposed in Appendix F. This construction is inspired from the one given in [9], which has been adapted to our specific setting. We emphasize that the size of the spectrum of $\mathcal{P}_{m}$ increases linearly with $m$, while the effect of the number of points $K$, the dimension $d$, and the spacing $\Delta$ between the location parameters $\{t_{1},\dots,t_{K}\}$ is translated in the initial constraint $m\gtrsim\sqrt{K}d^{3/2}\Delta_{+}^{-1}$.

We also state a complementary result, that will be useful for the proof of Theorem 10, $iii)$.

Corollary 8**.**

Let be given a set of $K$ points $S^{0}=\{t_{1},\ldots,t_{K}\}$ such that $\Delta:=\min_{k\neq\ell}\|t_{k}-t_{\ell}\|_{2}$. Let $m\gtrsim\sqrt{K}d^{3/2}\Delta_{+}^{-1}$. Then, for any $k\in[K]$, a function $\mathcal{Q}_{m}^{k}$ exists such that

[TABLE]

and a universal couple of constants $(\upsilon,\gamma)$ exists such that the function $\mathcal{Q}_{m}^{k}$ satisfies for $\epsilon=\frac{\upsilon}{md}$:

• $i)$

Near region $\mathds{N}_{k}(\epsilon)$: a positive constant $\widetilde{\mathcal{C}}$ exists such that:

[TABLE]

• $ii)$

Near region $\mathds{N}(\epsilon)\setminus\mathds{N}_{k}(\epsilon)$:

[TABLE]

• $iii)$

Far region $\mathds{F}(\epsilon)$:

[TABLE]

• $iv)$

A $c_{k,m}\in\mathds{L}$ exists such that $\mathcal{Q}_{m}^{k}=\Phi c_{k,m}$.

Proofs of $i),ii),iii)$ are similar to those of Theorem 7 and are omitted: the construction of $\mathcal{Q}_{m}^{k}$ obeys the same rules as the construction of $\mathcal{P}_{m}$ (the interpolation conditions only differ at points $t_{i},i\neq k$ and are switched from $1$ to [math]).

4.2 Bregman divergence $D_{\mathcal{P}_{m}}(\hat{\mu}_{n},\mu^{0})$

Below, the statistical loss between $\hat{\mu}_{n}$ and $\mu^{0}$ will be obtained in terms of the Bregman divergence associated to the dual certificate $\mathcal{P}_{m}$ obtained in Theorem 7. This divergence is defined by:

[TABLE]

We also introduce the term $\Gamma_{n}$ defined as

[TABLE]

which models the difference between the target $f^{0}=\Phi\mu^{0}$ and its empirical counterpart $\hat{f}_{n}$ in the RKHS. The next result provides a control between $\hat{\mu}_{n}$ and $\mu^{0}$ with the Bregman divergence.

Proposition 9**.**

Let $\mathcal{P}_{m}=\Phi c_{0,m}$ the dual certificate obtained in Theorem 7. Let $(\rho_{n})_{n\in\mathds{N}^{*}}$ be a sequence such that $\mathds{E}[\|\Gamma_{n}\|_{\mathds{L}}^{2}]\leq\rho_{n}^{2}$ for all $n\in\mathds{N}^{*}$. If $\kappa$ is chosen such that

[TABLE]

and if $\hat{\mu}_{n}$ is defined in ($\mathbf{P}_{\kappa}$), then:

• $i)$

For any integer $n$:

[TABLE]

• $ii)$

$c_{0,m}\in\mathds{L}$* satisfies*

[TABLE]

The proof of Proposition 9 is postponed to Section 6.1. The previous results indicate that the Bregman divergence between our estimator $\hat{\mu}_{n}$ and the target measure $\mu^{0}$ depends, up to some constants, on three main quantities:

• •

The variance of the empirical measure through the operator $L$ quantified by $\rho_{n}$,

• •

The Fourier transform $\sigma$ of the convolution kernel $\varphi$ over the interval $[-4m;4m]^{d}$. This term measures the ill-posedness of the inverse problem, which is associated to the difficulty to recover $\mu^{0}$ with indirect observations (here $f^{0}=\Phi\mu^{0}$ and we need to invert $\Phi$),

• •

The structure of the RKHS used to smooth the problem identified through the kernel $\lambda$.

Remark 5**.**

By using similar arguments to prove item ii) of Proposition 9, we can complete item (iv) of Corollary 8 as follows: A $c_{k,m}\in\mathds{L}$ exists such that $\mathcal{Q}_{m}^{k}=\Phi c_{k,m}$ and

[TABLE]

Remark 6**.**

We will derive from Proposition 9 some explicit convergence rates in each specific situation, i.e. as soon as the quantities involved in Equation (25) are made precise on some concrete examples. These rates will depend on the tuning parameter $m$ for solving the optimization problem ($\mathbf{P}_{\kappa}$), and on the choice of the kernel $\lambda$. Some examples will be discussed in Section 5. Indeed, $\kappa$ is related to $m$ through the relationship $\kappa=\rho_{n}/\|c_{0,m}\|_{\mathds{L}}$. Similarly, we will see in Section 5 that the kernel $\lambda$ is also linked to $m$ in a transparent way. We stress that according to Proposition 9, $m\gtrsim\sqrt{K}d^{3/2}\Delta_{+}^{-1}$. Such a condition will be satisfied provided $m$ is allowed to go to infinity with $n$ whereas $K,\Delta,d$ are kept fixed.

Remark 7**.**

The upper bound proposed in Proposition 9 only uses items (iii) and (iv) of Theorem 7. An enhanced control on the performances of $\hat{\mu}_{n}$ is provided in the next section. Alternative features will be also proposed with the alternative certificate $\mathcal{Q}_{m}$ introduced in Corollary 8 .

4.3 Statistical recovery of far and near regions

The next result sheds light on the performance of the BLASSO estimator introduced in Equation (10). The goodness-of-fit reconstruction of the mixture distribution $\mu^{0}$ by $\hat{\mu}_{n}$ is translated by the statistical properties of the computed weights of $\hat{\mu}_{n}$ around the spikes of $\mu^{0}$ (the support points $S^{0}=\{t_{1},\ldots,t_{K}\}$), which will define a family of $K$ near regions, as well as the behaviour of $\hat{\mu}_{n}$ in the complementary set, the far region. The sets $\mathds{F}(\epsilon)$ and $\mathds{N}(\epsilon)$ have already been introduced in Theorem 7. Our result takes advantage on the previous bounds and on $i)$ and $ii)$ of Theorem 7.

Theorem 10**.**

Let $m\gtrsim\sqrt{K}d^{3/2}\Delta_{+}^{-1}$ and let $\mathcal{P}_{m}$ be a dual certificate given in Theorem 7. Assume that $\hat{\mu}_{n}$ is the BLASSO estimator given by ($\mathbf{P}_{\kappa}$) with $\kappa=\kappa_{n}$ chosen in Proposition 9. Let $\mathcal{C}_{m}(\varphi,\lambda)$ be the quantity introduced in Proposition 9, $\hat{\mu}_{n}=\hat{\mu}_{n}^{+}-\hat{\mu}_{n}^{-}$ the Jordan decomposition of $\hat{\mu}_{n}$. A universal couple of constants $(\gamma,\upsilon)$ exists such that, if

[TABLE]

• $i)$

Far region and negative part:

[TABLE]

• $ii)$

Near region (spike detection): a positive constant $\mathcal{C}$ exists such that

[TABLE]

• $iii)$

Near region (weight reconstruction): for any $k\in[K]$:

[TABLE]

The proof of this important result is deferred to Section 6.2.

Remark 8**.**

It can be shown in specific situations (see, e.g., $iv)$ of Theorem 6) that the solution of ($\mathbf{P}_{\kappa}$) is indeed a discrete measure that can be written as

[TABLE]

In such a case, the relevance of the locations $\hat{S}$ of the reconstructed spikes $\hat{a}_{t}$ can be derived from the results of Theorem 10. A discussion is given in some specific cases in Section 5.

4.4 Support stability for large sample size

We introduce $\mathcal{P}_{0}:=\Phi c_{0}$ the “minimal norm certificate” (see e.g. [16]), which is defined by:

[TABLE]

when it exists.

We say that the support $S^{0}=\{t_{1},\ldots,t_{K}\}$ of $\mu^{0}$ satisfies the Non-Degenerate Bandwidth condition (NDB) if there exists $0<q<1$, $r>0$ and $\rho>0$ such that:

[TABLE]

We then have the support stability result for large values of $n$.

Theorem 11**.**

Let the triple $\lambda,\varphi,\mu^{0}$ be such that (NDB) holds. Let $r_{\kappa}\in(0,\frac{1}{2})$ and set $\kappa_{n}=\sqrt{\lambda(0)}\,n^{-r_{\kappa}}$. Let $\hat{\mu}_{n}$ be the BLASSO estimator ($\mathbf{P}_{\kappa}$) with a tuning parameter $\kappa=\kappa_{n}$.

Then for $n$ large enough, and with probability at least $1-Ce^{-n^{\frac{1}{2}-r_{\kappa}}}$ for a universal constant $C>0$, it holds that $\hat{\mu}_{n}$ has $K$ spikes with exactly one spike $\hat{t}_{k}$ in each region $\mathds{N}_{k}(r)$. These spikes converge to the true ones, and so do the amplitudes $\hat{a}_{k}$, as $n$ tends to infinity.

The proof can be found in Appendix E. We emphasize that $C$ is independent from the dimension $d$, from the RKHS used $\mathds{L}$ or the location of the spikes for example.

Remark 9**.**

In Theorem 11, note that the data fidelity kernel $\lambda$ is fixed but in practice, the bandwidth of $\lambda$ often depends on the sample size $n$. Theorem 11 suggests the heuristics that the data fidelity kernel $\lambda=\lambda_{n}$ may depend on $n$ and it might be such that $\kappa_{n}=\sqrt{\lambda_{n}(0)}\,n^{-r_{\kappa}}$ vanishes as $n$ tends to infinity.

Remark 10**.**

Assume that the mixing kernel $\varphi$ is such that $\varphi\star\lambda=\psi_{m}\star\lambda$ where $\psi_{m}$ is defined by (84), assume that the data fidelity kernel is such that $\lambda=\lambda_{1/(4m)}$ and assume that $\displaystyle m\gtrsim\sqrt{K}d^{3/2}\Delta_{+}^{-1}$. Then our certificate $\mathcal{P}_{m}$ is called the vanishing derivatives pre-certificate by [16, Section 4, Page 1335]. According to Theorem 7, we know that $\|\mathcal{P}_{m}\|_{\infty}\leq 1$. In this case, vanishing derivatives pre-certificate and certificate of minimal norm coincide so that $\mathcal{P}_{m}$ is the minimal norm certificate $\mathcal{P}_{0}$ appearing in (NDB), and Theorem 7 shows that (NDB) holds.

5 Rates of convergence for some usual mixture models

5.1 Frequency cut-off and sinc kernel

In this section, we describe the consequences of Theorem 10 for some mixture models with classical densities $\varphi$. For this purpose, we will consider the sinus-cardinal kernel sinc with a frequency cut-off $1/\tau$, which is introduced in Example 2. As a band-limited function $\lambda_{\tau}$, we have that

[TABLE]

Hence, to obtain a tractable version of Theorem 10 with $\mathcal{C}_{m}(\varphi,\lambda)<+\infty$ (see Equation (25)) we are led to consider $\tau$ such that

[TABLE]

In that case, $\mathcal{F}[\lambda_{\tau}]$ is a constant function over its support and the term $\mathcal{C}_{m}(\varphi,\lambda_{\tau})$ involved in Proposition 9 and Theorem 10 appears to be equal to

[TABLE]

To make use of Theorem 10, we also need an explicit expression of $(\rho_{n})_{n\in\mathds{N}^{*}}$, which itself strongly depends on the kernel $\lambda_{\tau}$. In this context, some straightforward and standard computations yield

[TABLE]

This provides a natural choice for the sequence $(\rho_{n})_{n\in\mathds{N}^{*}}$ as

[TABLE]

Therefore, the statistical rate obtained in Theorem 10 satisfies

[TABLE]

We should understand the previous inequality as an upper bound that translates a tradeoff between the sharpness of the window where spikes are located (given by $\epsilon=\mathcal{O}(1/(md))$ in (27)) and the associated statistical ability to recover a such targeted accuracy (given by the bound $\rho_{n}\mathcal{C}_{m}(\varphi,\lambda_{\tau})$ on the Bregman divergence). A careful inspection of the previous tradeoff leads to the following conclusion: the window size $\epsilon$ is improved for large values of $m$ but the statistical variability is then degraded according to the decrease rate of the Fourier transform $\sigma$ of $\varphi$, which typically translates an inverse problem phenomenon.

Finally, we emphasize that the dimensionality effect is not only involved in the term $2^{3d/2}$ of Equation (29) but is also hidden in the constraint

[TABLE]

used to build our dual certificate in Theorem 7. By the way, we stress that at the end, the only tuning parameter involved in ($\mathbf{P}_{\kappa}$) appears to be $m$.

We now focus our attention to some specific and classical examples in mixture models:

• •

the case of severely ill-posed inverse problems with an exponential decrease of the Fourier transform for large frequencies, which corresponds to super-smooth distributions. We emphasize that this class contains the standard benchmark of the Gaussian case, which will be discussed in details.

• •

the case of mildly ill-posed inverse problems which encompasses multivariate Laplace distributions, Gamma distributions, double exponentials among others.

5.2 Super-smooth mixture deconvolution and Gaussian case

5.2.1 Description of the distributions

We consider in this paragraph the statistically hard situation of the general family of mixing distribution $\varphi$ with an exponential decrease of the Fourier transform. More precisely, we assume that the spectral density $\sigma$ of $\varphi$ satisfies:

[TABLE]

where for any $j\in\mathds{N}^{\star}$, $\|.\|_{j}$ denotes the $\ell^{j}$-norm. The assumption ($\mathcal{H}^{supersmooth}_{\alpha,\beta}$) includes obviously the Gaussian distribution but also many other distributions as suggested by the list of examples displayed below (among others).

$\bullet$ The multivariate Cauchy distribution

For a dispersion parameter $\alpha$, $\varphi$ is defined by:

[TABLE]

$\bullet$ The tensor product of univariate Cauchy distribution

An alternative example is:

[TABLE]

$\bullet$ The multivariate Gaussian distribution

A standard benchmark study of the Gaussian law:

[TABLE]

5.2.2 General recovery result

In the situations covered by assumption ($\mathcal{H}^{supersmooth}_{\alpha,\beta}$), we shall observe that $\|t\|_{j}\leq d^{1/j}\|t\|_{\infty}$ and we verify that:

[TABLE]

In that case, we obtain that

[TABLE]

A straightforward application of Theorem 10 leads to the following result.

Proposition 12**.**

Assume that $\varphi$ satisfies $(\mathcal{H}^{supersmooth}_{\alpha,\beta})$. Let $m\gtrsim\sqrt{K}d^{3/2}\Delta_{+}^{-1}$. Let $\hat{\mu}_{n}$ be the BLASSO estimator given by ($\mathbf{P}_{\kappa}$) with $\kappa=\kappa_{n}$ chosen as in Proposition 9, then up to some universal constants (independent from $n,d,K$ and $m$):

• $i)$

Far region and negative part: if $\epsilon=\mathcal{O}\left(\frac{1}{md}\right)$, then:

[TABLE]

• $ii)$

Near region (spike detection): a couple of constants $(c,\mathcal{C})$ exists such that

[TABLE]

• $iii)$

Near region (weight reconstruction): for any $k\in[K]$:

[TABLE]

According to the results displayed in Proposition 12, our estimation procedure $\hat{\mu}_{n}$ leads to a consistent estimation as soon as $m$ is chosen as

[TABLE]

In such a case,

[TABLE]

and every set $A$ such that $\mathds{E}[\hat{\mu}_{n}^{+}(A)]\gtrsim n^{-\frac{1}{2}+\delta}$ is at least at a logarithmic distance ($\mathcal{O}(m^{-2})$) of a true spike.

We observe that as it is commonly observed in severely-ill conditioned inverse problems, we can expect only logarithmic rates of convergence. This logarithmic limitation in the super-smooth situation has been intensively discussed in the literature and we refer among others to [19]. To make the situation more explicit, we illustrate it in the Gaussian mixture model.

5.2.3 Multivariate Gaussian mixtures

As a specific case of super-smooth distribution with $\beta=j=2$ and $\alpha=1/2$, Proposition 12 holds and we obtain that if $m\gtrsim\sqrt{K}d^{3/2}\Delta_{+}^{-1}$ and if $\epsilon=\mathcal{O}(\frac{1}{md})$, then the weights of the far region and of the negative parts are upper bounded by:

[TABLE]

Similarly, a couple of constants $(c,\mathcal{C})$ exists such that:

[TABLE]

whereas the weights recovery is ensured by the following inequality: for any $k\in[K]$:

[TABLE]

• •

Quantitative considerations When the dimension $d$ is kept fixed (as the number of components $K$ and the minimal value for the spacings between the spikes $\Delta$), the statistical ability of the BLASSO estimator $\hat{\mu}_{n}$ is driven by the term $e^{8dm^{2}}/\sqrt{n}$. In particular, this sequence converges to [math] provided that the following condition holds:

[TABLE]

In other words, the maximal admissible value for $m$ is $\sqrt{\frac{\log(n)}{16d}}$. In particular, if we consider $m=\sqrt{\frac{\delta}{16}\frac{\log(n)}{d}}$ for $\delta$ small enough, we observe that

[TABLE]

The counterpart of this admissible size for $m$ is a slow rate for $\epsilon_{n}$:

[TABLE]

Said differently, the size of the near regions recovered with an almost parametric rate $n^{-1/2}$ are of the order $(d\log(n))^{-1/2}$.

• •

Nature of the results Item $i)$ of Proposition 12 and Equation (30) both indicate that the mass set by $\hat{\mu}_{n}$ on the negative part and on the far region tends to [math] as the sample size $n$ grows under Condition (32). Our estimator is consistent: the mass allowed on the near region will be close to $1$ as soon as $n$ is large enough. At this step, we stress that the parameter $m$ plays the role of an accuracy index: if $m$ is constant, the mass of the near region converges to $1$ at a parametric rate… but this near region is in this case not really informative. On the opposite hand, if $m$ is close to the limit admissible value expressed in (32), Item $ii)$ of Proposition 12 and Equation (31) translate the fact that the near region is close to the support of the measure $\mu^{0}$ but the convergence of the associated mass will be quite slow.

• •

Case of dimension 1 and number of spikes detection According to Item $ii)$ of Proposition 12 and Equation (31), any set with a sufficiently large mass is close to a true spike $(a_{k}^{0},t_{k})$ for some $k\in[K]$. We stress that in the specific situation where $d=1$, $\hat{\mu}_{n}$ is a discrete measure (see Theorem 6), namely

[TABLE]

In such a case, we get from Proposition 12 that if a reconstructed spike $(\hat{a}_{\hat{t}},\hat{t})$ is large enough, it is in some sense close to a true spike. More formally, if $m=\mathcal{O}(\sqrt{\delta\log(n)})$ and $\hat{t}\in\hat{S}$, then

[TABLE]

In particular, the BLASSO estimator $\hat{\mu}_{n}$ provides a lower bound on the number of true spikes. Once again, the value of $m$ is critical in such a case. In particular, according to (32), we cannot expect more than a logarithmic precision.

• •

Importance of the mixture parameters It is also interesting to pay attention to the effect of the number of components $K$, the size of the minimal spacing $\Delta$ and of the dimension $d$ on the statistical accuracy of our method. In the Gaussian case, the rate is of the order $K^{2}C^{d}e^{8dm^{2}}n^{-1/2}$ but an important effect is hidden in the constraint brought by Theorem 7:

[TABLE]

In particular, the behavior of our estimator is seriously damaged in the Gaussian situation when $(\Delta^{-1}\vee K\vee d)\rightarrow+\infty$ since in that case, taking the minimal value of $m$ satisfying the previous contraint, we obtain a rate of the order

[TABLE]

We observe that $d$, $K$ and $\Delta^{-1}$ cannot increase faster than a power of $\log(n)$: $d^{4}K\Delta^{-2}\ll\log(n)$.We will observe in Section 5.3 that a such hard constraint disappears in more favorable cases with smaller degrees of ill-posedness.

• •

Position of our result on Gaussian mixture models

To conclude this discussion, we would like to recall that the BLASSO estimator $\hat{\mu}_{n}$ depends on $m$. This parameter plays the role of a precision filter and only provides a quantification of the performances of our method. This is one of the main differences with the classical super-resolution theory where in general $m$ is fixed and constrained by the experiment. We should point out that many works have studied statistical estimation in Gaussian mixture models with a semi-parametric point of view (see, e.g. [38], [6]). These investigations are often reduced to the two-component case (K=2): we refer to [8], [22] or [21] among others. The general case ($K\in\mathds{N}^{*}$) has been for instance addressed in [26] using a model selection point of view: the selection of $K$ is achieved through the minimization of a criterion penalized by the number of components. We also refer to [7] where a Lasso-type estimator is built for mixture model using a discretization of the possible values of $t_{k}$. However, this last approach is limited by some constraints on the Gram matrix involved in the model that do not allow to consider situations where $\Delta$ is small: in [7], the minimal separation between two spikes has to satisfy $\Delta\geq\Delta_{0}>0$, i.e. has to be lower bounded by a positive constant $\Delta_{0}$, which depends on the mixing distribution $\varphi$. We emphasize that in our work, we only need an upper bound on $K$ and a lower bound on $\Delta$ or at least to assume that these quantities are fixed w.r.t. $n$. According to Proposition 12, our constraint expressed on these parameters already allows to cover a large number of interesting situations.

5.3 Ordinary smooth distributions

General result

Ordinary smooth distributions investigated in this section are described through a polynomial decrease of their Fourier transform. The corresponding deconvolution problem is then said to be mildly ill-posed. In this section, we assume that the density $\varphi$ satisfies

[TABLE]

We refer to [19] and the references therein for an extended description of the class of distributions involved by $(\mathcal{H}^{smooth}_{\beta})$ and some statistical consequences in the situation of standard non-parametric deconvolution (see also the end of this section for two examples). For our purpose, it is straightforward to verify that

[TABLE]

In that case, we obtain that

[TABLE]

We then deduce the following result (which is a direct application of Theorem 10).

Proposition 13**.**

Assume that $\varphi$ is ordinary smooth and satisfies $(\mathcal{H}^{smooth}_{\beta})$. Consider $m\gtrsim\sqrt{K}d^{3/2}\Delta_{+}^{-1}$. Let $\hat{\mu}_{n}$ be the BLASSO estimator given by ($\mathbf{P}_{\kappa}$) with $\kappa=\kappa_{n}$ chosen as in Proposition 9, then up to universal constants (independent from $n,d,K$ and $m$):

• $i)$

Far region and negative part: if $\epsilon=\mathcal{O}\left(\frac{1}{md}\right)$, then:

[TABLE]

• $ii)$

Near region (spike detection): a couple of constants $(c,\mathcal{C})$ exists such that

[TABLE]

• $iii)$

Near region (weight reconstruction): for any $k\in[K]$:

[TABLE]

The proof of this proposition is omitted, and we only comment on the consequences of this result for ordinary smooth mixtures. Provided $K,d$ and $\Delta$ are bounded (or fixed), we obtain a consistent estimation with the BLASSO estimator $\hat{\mu}_{n}$ when $m$ is chosen such that

[TABLE]

In such a case, $\epsilon_{n}=\mathcal{O}(d^{-1}n^{-\delta})$. Now, if $K\vee d\vee\Delta^{-1}$ is allowed to grow towards $+\infty$, setting $m\sim\sqrt{K}d^{3/2}\Delta_{+}^{-1}$ (the minimal value satisfying the constraint (4)) leads to a bound of order

[TABLE]

In particular, the maximal order for the dimension is $\mathcal{O}(\log(n))$ as $n\rightarrow+\infty$. In the same way, the minimal size of spacings to permit a consistent estimation should not be smaller than $n^{-1/(2\beta)}$. In particular, this indicates that a polynomial accuracy is possible (see e.g. Item $ii)$ of Proposition 13). This emphasized the strong role played by the mixture density $\varphi$ in our analysis. We present below two specific examples of ordinary smooth mixture density.

Multivariate Laplace distributions

In such a case:

[TABLE]

We obtain here an ordinary smooth density with $\beta=2$. The minimal spacing for a discoverable spike is therefore of the order $n^{-1/4}$ while the constraint on the dimension is not affected by the value of $\beta$. Concerning the number of components $K$, its value should not exceed $n^{1/6}$ and the smallest size of the window $\epsilon_{n}$ is $n^{-1/4}$.

Tensor product of Laplace distributions

Another interesting case is the situation where $\varphi$ is given by a tensor product of standard Laplace univariate distributions:

[TABLE]

In that case, $\beta=2d$ and the previous comments apply: the maximal value of $m$ is $n^{1/4d}$ with an optimal size of the window of the order $n^{-1/(4d)}$ whereas $d$ should be at least of order $\mathcal{O}(\log(n))$.

6 Proof of the Main Results

6.1 Analysis of the Bregman divergence

This paragraph is devoted to the statistical analysis of the Bregman divergence whose definition is recalled below:

[TABLE]

Proof of Proposition 9.

According to the definition of $\hat{\mu}_{n}$ as the minimum of our variational criterion (see Equation (15)), we know that:

[TABLE]

Proof of $i)$. With our notation $\Gamma_{n}=L\hat{f}_{n}-L\circ\Phi\mu^{0}$ introduced in Section 4.2, we deduce that:

[TABLE]

Using now $\mathcal{P}_{m}$ obtained in Theorem 7, we deduce that

[TABLE]

Hence, we deduce the following upper bound on the Bregman divergence:

[TABLE]

According to Theorem 7, $\mathcal{P}_{m}=\Phi c_{0,m}$ for some $c_{0,m}\in\mathds{L}$. In particular, we get

[TABLE]

where the last equality comes from the self-adjoint property of $\Phi$ in $L^{2}(\mathds{R}^{d})$. The reproducing kernel relationship yields:

[TABLE]

Gathering (34) and (35), we deduce that:

[TABLE]

Using now a straightforward computation with $\|.\|_{\mathds{L}}$, we conclude that:

[TABLE]

Since the first term of the left hand side is positive, the previous inequality leads to:

[TABLE]

where we have used $\|a+b\|_{\mathds{L}}^{2}\leq 1.5\|a\|_{\mathds{L}}^{2}+3\|b\|_{\mathds{L}}^{2}$ with $a=\Gamma_{n}$ and $b=-\kappa c_{0,m}/2$ for the right hand side. We now consider a sequence $(\rho_{n})_{n\in\mathds{N}^{*}}$ such that $\mathds{E}[\|\Gamma_{n}\|_{\mathds{L}}^{2}]\leq\rho_{n}^{2}$ for all $n\in\mathds{N}^{*}$ and we choose:

[TABLE]

Then we deduce from (36) that:

[TABLE]

Proof of $ii)$. We now derive an upper bound on $\|c_{0,m}\|_{\mathds{L}}$. Recall that according to ($\mathcal{H}_{0}$) and in particular (12) we have:

[TABLE]

Since $\varphi$ is symmetric and $\Phi^{\star}=\Phi$, we have according to Theorem 7 that:

[TABLE]

Indeed, $iii)$ of Theorem 7 entails that the support of the Fourier transform of $\mathcal{P}_{m}$ is contained in $[-4m,4m]^{d}$. This embedding, together with ($\mathcal{H}_{\infty}$) entails:

[TABLE]

which provides the last inequality. The inequality (38) can be rewritten as:

[TABLE]

We use (37), (39) and observe that $|\mathcal{F}[\varphi]|=\sigma$ to conclude the proof. ∎

6.2 Near and Far region estimations

In this paragraph, we provide the main result of the paper that establishes the statistical accuracy of our BLASSO estimation.

Proof of Theorem 10.

Proof of $i)$ In a first time, we provide a lower bound on the Bregman divergence. This bound takes advantage on the properties of the dual certificate associated to Theorem 7. First remark that

[TABLE]

since $\mathcal{P}_{m}(t_{k})=1$ for all $k$. This inequality yields the positiveness of the Bregman divergence:

[TABLE]

Now, using similar arguments and the Borel’s decomposition $\hat{\mu}_{n}=\hat{\mu}_{n}^{+}-\hat{\mu}_{n}^{-}$, we obtain

[TABLE]

Proposition 9 then implies that:

[TABLE]

Weight of the negative part. Since the dual certificate $\mathcal{P}_{m}$ is always positive, we have

[TABLE]

Moreover, according to item $ii)$ of Theorem 7,

[TABLE]

Therefore, we obtain that:

[TABLE]

Finally, the first part of $i)$ of Theorem 10 is a direct consequence of (40)-(42).

Weight of the far region. We consider $\gamma$ such that $d^{3}\geq\gamma\upsilon^{2}$ and we know that in the far region:

[TABLE]

Thus,

[TABLE]

We then conclude, using the previous expectation upper bound, that:

[TABLE]

Proof of $ii)$. Thanks to Theorem 7, we have:

[TABLE]

Then, for any subset $A\subset\mathds{R}^{d}$,

[TABLE]

Equations (40) and (43) lead to:

[TABLE]

Then,

[TABLE]

Proof of $iii)$. The idea of this proof is close to the one of [2, Theorem 2.1]. We consider the function $\mathcal{Q}_{m}^{k}$ given by Corollary 8 that interpolates $1$ at $t_{k}$ and [math] on the other points of the support of $\mu^{0}$. From the construction of $\mathcal{Q}_{m}^{k}$, we have that:

[TABLE]

We then use the decomposition:

[TABLE]

Study of $B+C+D$. On the set $\mathds{F}(\epsilon)$, we use that $\mathcal{Q}_{m}^{k}\leq 1-\gamma\frac{\upsilon^{2}}{d^{3}}$ so that:

[TABLE]

For the term $C$, we use the upper bound satisfied by $\mathcal{Q}_{m}^{k}$ in $\bigcup_{i\neq k}\mathds{N}_{i}(\epsilon)$ and obtain that:

[TABLE]

Finally, for $B$, we use that on the set $\mathds{N}_{k}(\epsilon)$, we have $|\mathcal{Q}_{m}^{k}-1|\leq\widetilde{\mathcal{C}}m^{2}\|t-t_{k}\|_{2}^{2}$. Therefore, we have:

[TABLE]

We then conclude that:

[TABLE]

Study of $A$. We use that $\mathcal{Q}_{m}^{k}$ may be written as:

[TABLE]

Since $\Phi$ is self-adjoint in $L^{2}$, we shall write that:

[TABLE]

where we used the Cauchy-Schwarz inequality and the triangle inequality in the last line. We then use (33) and obtain:

[TABLE]

Since we have obtained the positiveness of the Bregman divergence, we then conclude that:

[TABLE]

The Cauchy-Schwarz inequality yields:

[TABLE]

This inequality holds for any value of $\kappa$ and we choose:

[TABLE]

Using this value of $\kappa$, we then obtain:

[TABLE]

Now, we define $\Box_{n}=\|L\hat{f}_{n}-L\circ\Phi\hat{\mu}_{n}\|_{\mathds{L}}\|\Gamma_{n}\|_{\mathds{L}}^{-1}$ and remark that:

[TABLE]

This last inequality implies that $\Box_{n}\leq 1+\sqrt{3}$, which leads to:

[TABLE]

We then come back to $A$ and write that:

[TABLE]

Final bound. We use Equations (46) and (45) in the decomposition given in Equation (44) and obtain that:

[TABLE]

Finally, we conclude the proof using Equation (26) and $ii)$ of Proposition 9:

[TABLE]

∎

Appendix A Tractable Algorithms for BLASSO Mixture Models

We sketch three algorithms to compute approximate or exact solutions to (15) following the “*​off-the-grid *” methodology, e.g., [5, 37, 16, 2, 14]. This methodology searches in a gridless manner the location of the support points $t_{i}$ of the mixture distribution $\mu^{0}$. We will present the following methods:

• •

Greedy methods provide heuristic and theoretical results such as “Sliding Frank Wolfe ” [15] (also known as conditional gradient with “sliding” step) or “Continuous Orthogonal Matching Pursuit ” [24, 18]. We describe these methods in Section A.2.

• •

We discuss in Section A.3 on Conic Gradient Descent using particles. Here, $\hat{\mu}_{n}$ is approximated by a cloud of particles that is optimized all along a set of iterations.

We emphasize that if Beurling-LASSO has been studied in the past decade, the formulation (15) has two new important features. First the observation is a sample from a mixing law. Second, the data fidelity term has been tuned to incorporate a low pass filter kernel $\lambda$. The next paragraph carefully introduces these new features into the latter algorithms.

A.1 Notation for algorithm design solving BLASSO Mixture Models

We call that primal and dual convex programs of BLASSO for Mixture Models (15) are given by Theorem 6, and that strong duality holds, leading to Equation (17).

Gradient of the data fidelity term

The data fidelity term defined

[TABLE]

can be seen is related to the real-valued function $F$ on the space of measures $\mathcal{M}(\mathds{R}^{d},\mathds{R})$ endowed with the total-variation norm $\|\cdot\|_{1}$, namely:

[TABLE]

whose Fréchet differential at point $\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$ in the direction $\nu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$ is:

[TABLE]

Thanks to the convolution by $\varphi$ endowed in $\Phi$, the gradient $\nabla F(\mu)$ is given by:

[TABLE]

Dual functions

By (17), note that the dual function $\hat{\eta}_{n}$ is such that $\hat{\eta}_{n}=-{\nabla\mathrm{F}(\hat{\mu}_{n})}/{\kappa}$. Indeed, for a given $\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$, one may define its dual function by $\eta_{\mu}$:

[TABLE]

and we observe that $\eta_{\mu}$ corresponds to a residual , which involves the difference between $\hat{f}_{n}$ and $\Phi\mu$ smoothed by the convolution operator $L$. Its gradient is given by:

[TABLE]

Data fidelity term for measures with finite support

We pay a specific attention to discrete measures with finite support, namely a finite sum of Dirac masses. Given a number of atoms $N\geq 1$, of weights $a\in\mathds{R}^{N}$ and locations $t=(t_{1},\ldots,t_{N})\in(\mathds{R}^{d})^{N}$, we denote by:

[TABLE]

By Proposition 5, $a\longmapsto\mathrm{F}_{N}(a,t)$ is a positive semi-definite quadratic form and

[TABLE]

where

[TABLE]

A.2 Greedy methods: Sliding Frank-Wolfe / Continuous Orthogonal Matching Pursuit

Sliding Frank-Wolfe algorithm (SFW)

The Frank-Wolfe algorithm is an interesting avenue for solving differentiable convex programs on weakly compact convex sets, see [15] and references therein for further details. Stricto sensu, $(\mathbf{P}_{\kappa})$ is convex but not differentiable and the feasible set is convex but not weakly compact. Following [15, Lemma 4], note that $\hat{\mu}_{n}$ is a minimizer of $(\mathbf{P}_{\kappa})$ if and only if $(\|\hat{\mu}_{n}\|_{1},\hat{\mu}_{n})$ minimizes:

[TABLE]

and this latter program is a differentiable convex program on weakly compact convex set (for the weak-$\star$ topology).

Hence we can invoke the Frank-Wolfe scheme to compute approximate solutions to BLASSO Mixture Models (15). Unfortunately, the generated measures $\mu^{(k)}$ along this greedy algorithm are not very sparse compared to $\hat{\mu}_{n}$: each Dirac mass of $\hat{\mu}_{n}$ is approximated by a multitude of Dirac masses of $\mu^{(k)}$ with an inexact positions. This is why the improvement of sliding the Frank Wolfe algorithm is suggested in [15]: Equation (19) involved in the resolution of the BLASSO allows to move the Dirac masses. Algorithm 1 is the sliding Frank-Wolfe [15, Algorithm 2] adapted for the resolution of the BLASSO Mixture Models (15).

Following the analysis of [15], let us discuss the steps at lines $3$, $4$, $7$ and $8$ of Algorithm 1:

• •

Line $3$: This step is an optimal gradient step, with the notations in Section A.1:

[TABLE]

Note also that this step is the costly step of the algorithm since is relies on a black-box optimizer computing the global maximum of $|\eta^{(k)}|$. In general, this is done using a grid search and finding a local maxima by gradient descent.

• •

Line $4$: the stopping condition implies that $\hat{\mu}^{(k)}$ is an exact solution and hence, $\eta^{(k)}$ is the dual function such that Equation (17) holds. In this sense, we may say that SFW iteratively construct a dual function such that (17) holds.

• •

Line $7$ (LASSO Step): note that the support is fixed and we are optimizing on the amplitudes $a$. It amounts in solving a standard LASSO, which can be efficiently done using proximal forward-backward schemes such as FISTA for instance.

• •

Line $8$: it requires solving a non-convex optimization program (19). As mentioned in [15], one does not need to exactly solve this program and their main result (namely finite convergence of the algorithm, see below) pertains if $(a^{(k+1)},t^{(k+1)})$, obtained by a gradient step initialized in $(a^{(k+\frac{1}{2})},t^{(k+\frac{1}{2})})$, diminishes the objective function. This is done by a bounded Broyden-Fletcher-Goldfarb-Shanno (BFGS) method in [15], which is numerically shown to be rapid with few iterations needed.

An important feature of SFW is that it can actually lead to exact solutions. Under a Non-Degeneracy Condition alike to (NDB), the result in [15, Theorem 3] proves that Algorithm 1 recovers exactly $\hat{\mu}_{n}$ in a finite number of steps. They also show [15, Proposition 5] that the generated measure sequence $(\mu^{(k)})_{k}$ converges towards $\hat{\mu}_{n}$ for the weak-$\star$ topology.

As mentioned in [15, Remark 8], Algorithm 1 can be adapted to build a positive measure as follows

• •

the stopping condition $|\eta^{(k)}(t_{\star}^{(k)})|\leq 1$ becomes $\eta^{(k)}(t_{\star}^{(k)})\leq 1$;

• •

the LASSO is solved on $a\in\mathds{R}_{+}^{N^{(k)}+1}$;

• •

the step (19) is solved on $\mathds{R}_{+}^{N^{(k)}+1}\times(\mathds{R}^{d})^{N^{(k)}+1}$.

Continuous Orthogonal Matching Pursuit (COMP)

Continuous Orthogonal Matching Pursuit (COMP) [18] is another greedy approach that is the Orthogonal Matching Pursuit approximation algorithm [30] adapted in the context of continuous parametric dictionaries. This framework fits ours and COMP can be applied to Mixture Models estimation. Continuous Orthogonal Matching Pursuit (COMP) is an iterative algorithm that add a Dirac mass one at the time, building a sequence of measures $\hat{\mu}^{(k)}$, but it does not solve the BLASSO Mixture problem (15) per se. Nevertheless it builds a sequence of dual functions $\eta^{(k)}=\frac{1}{\kappa}\Phi(L\hat{f}_{n}-L\circ\Phi\hat{\mu}^{(k)})$, referred to as the residual in the framework of COMP. The Dirac mass added to the model is defined as in Line 3 of Algorithm 1 (SFW for BLASSO Mixture Models) but the weights $a$ are updated differently (we referred to [24, Algorithm 1]): alternating between hard-thresholding [24, Step 3 in Algorithm 1] and even some sliding-flavour step [24, Step 5 in Algorithm 1].

A.3 Conic Particle Gradient Descent (CPGD)

Conic Particle Gradient Descent [11] is an alternative promising avenue for solving BLASSO for Mixture Models (15). The idea is to discretize a positive measure into a system of particles, i.e. a sum of $N$ Dirac masses, by:

[TABLE]

with $a_{i}=r_{i}^{2}$. We observe that the objective function involved in the minimization of Equation (15) is given by:

[TABLE]

where $F_{N}$ is defined in Equation (47). As already emphasized in [11], Equation (15) is a convex program in $\mu$ whereas the parametrization given in $H$ translates this minimization into a non-convex differentiable problem in terms of $r$ and $t$. This function $H$ can be seen as an instance of the BLASSO Equation (15) for the measure $\mu_{a,t}$, namely a convex program that does not depends on the number of Particles $N$. All the more, it is possible to run a gradient descent on positions $t_{i}\in\mathds{R}^{d}$ and weights $r_{i}>0$ of the $N$ particles system. The crucial ingredient is then to implement a gradient descent on the lifted problem in the Wasserstein space approximating the Wasserstein gradient flow. For two step-sizes $\alpha>0$ and $\beta>0$, and for any position $(r,t)$, we define the Riemannian inner product by: $\forall(\delta r_{1},\delta r_{2})\in\mathds{R}_{+}^{2}\quad\forall(\delta t_{1},\delta t_{2})\in\{\mathds{R}^{d}\}^{2}$:

[TABLE]

The gradient w.r.t. this conic metric is given by:

[TABLE]

and the Wasserstein gradient [11, Section 2.2] is $g_{\mu}(r,t)=(-2\alpha r\kappa(\eta_{\mu}(t)-1),-\beta\kappa\nabla\eta_{\mu}(t))$ for a.e. point $(r,t)\in\Omega$.

Using standard mean-field limits of gradient flows in Wasserstein space (see e.g. [34]), it is possible to prove that the approximate gradient flows of the $N$-particles system converge towards the gradient flow on the Wasserstein space when $N\longrightarrow+\infty$ (see Theorem 2.6 and Theorem 3.5 of [12]). Hence, for a large enough number of particles $N$, it then implies the convergence towards the global minimizer of $\mu_{a,t}\longmapsto F(\mu_{a,t})$ itself, despite the lack of convexity of the function $(r,t)\longmapsto F_{N}(r,t)$. We refer to Theorem 3.9 of [11] that establishes the convergence of the particle gradient descent with a constant step-size under some non-degeneracy assumptions, i.e. the convergence of the CPGD toward $\hat{\mu}_{n}$ (15) in Hellinger-Kantorovich metric, and hence in the weak-$\star$ sense. Furthermore, Theorem 3.9 provides an exponential convergence rate: CPGD has a complexity scaling as $\log(1/\varepsilon)$ in the desired accuracy $\varepsilon$, instead of $\varepsilon^{-1/2}$ for general accelerated convex methods.

Appendix B Proofs related to the kernel construction

B.1 Convolution in the RKHS

Proof of Proposition 4.

Consider $\mathcal{A}:f\longmapsto x\mapsto\int\ell(x-y)f(y)dy$, $\mathcal{A}$ is a self-adjoint operator. We denote by $(w_{k})_{k\geq 1}$ the non-negative eigenvalues of $\mathcal{A}$ and $(\psi_{k})_{k\geq 1}$ the associated eigenvectors. We shall remark that the following equality holds:

[TABLE]

while $\mathds{L}$ corresponds to the next Hilbert space

[TABLE]

We now consider a non-negative measure $\nu$ and we remark that

[TABLE]

We observe that the coefficients of $L\nu$ are $c_{k}(L\nu)=w_{k}\int\psi_{k}(y)\nu(y)dy$. We shall remark that

[TABLE]

The Jensen inequality yields

[TABLE]

where the last equality comes from the Tonelli Theorem. We then observe that

[TABLE]

giving the result. ∎

B.2 Computation of the data-fidelity terms

Proof of Proposition 5.

Recall that $\Phi\mu\in\mathcal{C}_{0}(\mathds{R}^{d},\mathds{R})\cap L^{1}(\mathds{R}^{d})$. Now, given $f\in\mathcal{C}_{0}(\mathds{R}^{d},\mathds{R})\cap L^{1}(\mathds{R}^{d})$, one can consider the measure $\mu$ with signed density function $f$ and we may define:

[TABLE]

The embedding $L$ allows to compare $\hat{f}_{n}$ with $\Phi\mu$ in $\mathcal{C}_{0}(\mathds{R}^{d},\mathds{R})\cap L^{1}(\mathds{R}^{d})$. One has:

[TABLE]

Replacing $f^{0}$, which is unknown, by the empirical measure $\hat{f}_{n}$ in the previous equation leads to the following criterion:

[TABLE]

In particular, for all $\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$ note that $\Phi\mu\in\mathcal{C}_{0}(\mathds{R}^{d},\mathds{R})\cap L^{1}(\mathds{R}^{d})$ and introduce the criterion:

[TABLE]

which will be investigated in this paper. Note that it holds

[TABLE]

as claimed. ∎

Appendix C Perfect recovery properties - Theorem 2 and Theorem 3

C.1 Perfect recovery

Proof of Theorem 2.

Remark first that under ($\mathcal{H}_{0}$) and ($\mathcal{H}_{\infty}$), the RKHS, denoted by $\mathds{H}$, generated by the kernel $h(.,.)=\varphi(.-.)$ is dense in $\mathcal{C}_{0}(\mathds{R}^{d},\mathds{R})$ with respect to the uniform norm, see [10, Proposition 5.6] for instance. Furthermore, using [36, Proposition 2], we can show that ($\mathcal{H}_{\infty}$) implies that the embedding $\Phi$ is injective onto $\mathds{H}$. This means that we have identifiability of $\mu$ from the knowledge of $\Phi\mu$. More precisely, denote $f^{0}:=\Phi\mu^{0}$, we deduce that if it holds $\|f^{0}-\Phi\mu\|_{\mathds{H}}^{2}=0$ then one has $\mu=\mu^{0}$. ∎

C.2 Perfect recovery with a dual certificate

Proof of Theorem  3.

Let

[TABLE]

Step 1: Support inclusion. We observe that both $\hat{\mu}$ and $\mu^{0}$ belong to $\mathcal{M}(f^{0})$ so that $\Phi\hat{\mu}=\Phi\mu^{0}$. Hence, considering the Fourier transform on both sides and using $\mathcal{F}[\Phi(\mu)]=\mathcal{F}[\varphi]\mathcal{F}[\mu]$, we get that $\sigma\mathcal{F}(\hat{\mu})=\sigma\mathcal{F}(\mu^{0})$ which is equivalent to $\mathcal{F}(\hat{\mu})=\mathcal{F}(\mu^{0})$ on the support of $\sigma$. Now, Assumption ($\mathcal{H}_{\eta}$) yields:

[TABLE]

Denote by $q_{\eta}:=\mathcal{F}(\mathcal{P}_{\eta})$ the Fourier transform of $\mathcal{P}_{\eta}$. By assumption, the support of $q_{\eta}$ is included in $[-\eta,\eta]^{d}$ and from (48) we get that:

[TABLE]

Since $\mathcal{P}_{\eta}\in L^{1}(\mathds{R}^{d})$, the Riemann-Lebesgue lemma shows that $q_{\eta}$ is continuous. Recall also that $q_{\eta}$ has a compact support so we deduce that $q_{\eta}\in L^{1}(\mathds{R}^{d})$. By Fourier inversion theorem, we have

[TABLE]

namely

[TABLE]

Remark that $\mathcal{P}_{\eta}$ satisfies

[TABLE]

and the Hölder inequality leads to

[TABLE]

From the definition of $\hat{\mu}$, one also has $\|\hat{\mu}\|_{1}\leq\|\mu^{0}\|_{1}$. Putting everything together, we deduce that

[TABLE]

Since $\mathcal{P}_{\eta}$ is continuous and strictly lower than one outside of the support of $\mu^{0}$, we deduce from the above equality that the support of $\hat{\mu}$ is included in the support of $\mu^{0}$:

[TABLE]

Step 2: Identifiability and conclusion. We prove that $\{\varphi(\cdot-t_{1}),\ldots,\varphi(\cdot-t_{K})\}$ spans a vector subspace of $\mathcal{C}_{0}(\mathds{R}^{d},\mathds{R})\cap L^{1}(\mathds{R}^{d})$ of dimension $K$. This proof is standard and relies on a Vandermonde argument. We assume first that $K$ coefficients $x_{1},\ldots,x_{K}\in\mathds{R}$ exist such that:

[TABLE]

Applying the Fourier transform and using $\varphi=\mathcal{F}^{-1}[\sigma]$, we deduce that:

[TABLE]

Since $\sigma$ is nonzero, there exists an open set $\Omega\subseteq\mathds{R}^{d}$ such that $\sigma>0$ on $\Omega$. We deduce that:

[TABLE]

Now, we can choose some points $u_{j}$ in $\Omega$ so that the Fourier matrix with entries $(e^{\imath u_{j}^{\top}t_{k}^{0}})_{kj}$ is invertible. It implies that $x_{k}=0$ and $\{\varphi(\cdot-t_{1}),\ldots,\varphi(\cdot-t_{K})\}$ spans a subspace of dimension $K$.

We now conclude the proof. We know from Step 1 that

[TABLE]

Since $\hat{\mu}$ and $\mu^{0}$ belong to $\mathcal{M}(f^{0})$, then

[TABLE]

which in turn implies that $x_{k}=a_{k}^{0}$ for all $k\in\{1,\ldots,K\}$, namely $\hat{\mu}=\mu^{0}$. ∎

Appendix D Primal-Dual problems and duality gap

Proof of Theorem 6.

We consider some primal variables $\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$ and $z\in\mathds{L}$ and introduce the dual variable $c\in\mathds{L}$. The Lagrangian is given in Equation (49)

[TABLE]

and we consider the dual problem ($\mathbf{P}^{*}_{\kappa}$).

Proof of $i)$. The existence of some solutions $\hat{\mu}_{n}$ to the primal problem ($\mathbf{P}_{\kappa}$) is obtained with the help of a standard argument: we prove that the primal objective function is a proper lower semi-continuous (for the weak-* topology) convex function on the Banach space $\mathcal{M}(\mathds{R}^{d},\mathds{R})$.

We now consider the “invariant property” related to the solutions of ($\mathbf{P}_{\kappa}$). The norm $\|\cdot\|_{\mathds{L}}$ satisfies

[TABLE]

Now consider two primal solutions $\mu_{1}$ and $\mu_{2}$ and define $\tilde{\mu}=(\mu_{1}+\mu_{2})/2$. Using (50) and the triangle inequality for $\|.\|_{1}$, one has:

[TABLE]

But, remind that:

[TABLE]

We then conclude that $\tilde{\mu}$ is also a solution to the primal problem and that $L\circ\Phi\mu_{1}=L\circ\Phi\mu_{2}$. We can repeat this argument for any pair of primal solutions so that the quantity $\hat{z}_{n}:=L\circ\Phi\hat{\mu}_{n}$ is uniquely defined and does not depend on the choice of the primal solution point $\hat{\mu}_{n}$. It also implies that $\hat{m}_{n}:=\|\hat{\mu}_{n}\|_{1}$ is uniquely defined (does not depend on the choice of the primal solution point).

Proof of $ii)$. We shall write the dual program ($\mathbf{P}^{*}_{\kappa}$) as follows: consider dual variable $c$ and write:

[TABLE]

and the previous infimum appears to be splitted in terms of the influence of $z$ and $\mu$. Optimizing in $z$ the first term \raisebox{-.9pt} {1}⃝ leads to $z=L\hat{f}_{n}-c$ so that:

[TABLE]

The second term \raisebox{-.9pt} {2}⃝ is more intricate. Observe that:

[TABLE]

and the Fubini yields:

[TABLE]

We deduce that:

[TABLE]

We use the $L^{1}-L^{\infty}$ Hölder inequality, namely $\int_{\mathds{R}^{d}}\Phi c\,\mathrm{d}\mu\leq\|\Phi c\|_{\infty}\|\mu\|_{1}$, which yields:

[TABLE]

Hence, we conclude that:

[TABLE]

where $\boldsymbol{I}_{\{\|\Phi c\|_{\infty}\leq\kappa\}}(c)$ is the constraint $\|\Phi c\|_{\infty}\leq\kappa$, namely it is [math] if $c$ such that $\|\Phi c\|_{\infty}\leq\kappa$ and $\infty$ otherwise. Finally, we obtain that for a fixed dual variable $c$:

[TABLE]

Hence, the dual problem ($\mathbf{P}^{*}_{\kappa}$) shall be written as

[TABLE]

Here again the dual objective function is lower semi-continuous and coercive on the Hilbert space $\mathds{L}$ so a minimizer $\hat{c}$ exists. Again, Inequality (50) implies the uniqueness of $\hat{c}$.

To prove that there is no duality gap, we use the Slater condition: we remark that a feasible point $(c^{\circ})$ exists in the interior of the constrained set $\{\|\Phi c\|_{\infty}\leq\kappa\}$. Now, the generalized Slater condition shall be used (see e.g. [32]). Indeed, given any nonzero $c\in\mathds{L}\subseteq\mathcal{C}_{0}(\mathds{R}^{d},\mathds{R})$, note that the convolution operator satisfies $\|\Phi c\|_{\infty}\leq\|c\|_{\infty}$. Hence, we set $c^{\circ}=\kappa c/(2\|c\|_{\infty})$ and these points are in the interior of the constrained set. The generalized Slater condition implies that strong duality holds, and there is no duality gap (i.e., strong duality holds):

[TABLE]

Furthermore, note that $\hat{z}_{n}:=L\circ\Phi\hat{\mu}_{n}$ and $z_{n}=L\hat{f}_{n}-\hat{c}$ from (51) and strong duality, we deduce that

[TABLE]

Proof of $iii)$. We consider the unique $\hat{c}$ solution to

[TABLE]

and the strong duality implies that:

[TABLE]

Since $\Phi\hat{c}$ is continuous, we verify, using the argument of Lemma A.1 in [13], that:

[TABLE]

where we recall that $\Phi\hat{c}\in L^{\infty}(\mathds{R}^{d})$ is such that its supremum norm is less than $\kappa$.

Proof of $iv)$. The last point is a consequence of the Schwartz-Paley-Wiener Theorem (see e.g. Theorem XVI, chapter VII in [35, Page 272]). Indeed, note that $\Phi\hat{c}$ is a continuous function whose inverse Fourier transform has a support included in the support of $\sigma\times\Lambda$. By assumption, this latter is bounded and one may apply the Schwartz-Paley-Wiener Theorem: we deduce that $\Phi\hat{c}$ can be extended to complex values $\mathds{C}^{d}$ into an analytic entire function of exponential type. In particular, $\Phi\hat{c}\pm\kappa$ has isolated zeros one the real line, which concludes the proof. ∎

Appendix E Support stability

Proof of Theorem 11.

We follow the ideas of [16] for the proof of Theorem 11. Consider the convex program

[TABLE]

whose Lagrangian expression is, for all $\mu\in\mathcal{M}(\mathds{R}^{d},\mathds{R})$, $c\in\mathds{L}$,

[TABLE]

using (52) in the last equation. Now, Equation (53) yields that the dual program is:

[TABLE]

Note also that the objective function of the dual program satisfies:

[TABLE]

which gives the following equivalent formulation of the dual:

[TABLE]

Note that the dual certificate $\mathcal{P}_{m}$ exists, then we know that $\mu^{0}$ is a solution to ($\mathbf{P}_{0}$) by Theorem 3. As in Section D, we use the Slater condition to prove that there is no duality gap: we remark that a feasible point $c$ exists in the interior of the set $\{\|\Phi c\|_{\infty}\leq 1\}$. Now, the generalized Slater condition shall be used (see e.g. [32]). We get that any solution $c$ to ($\mathbf{D}_{0}$) satisfies that $\Phi c$ is a sub-gradient of the total variation norm at point $\mu^{0}$. We recall the definition of condition (NDB).

[TABLE]

where $\rho>0$. Under condition (NDB), we know that $\mathcal{P}_{0}:=\Phi c_{0}$ is a solution to ($\mathbf{D}_{0}$).

Consider also the following convex program:

[TABLE]

which is the same as the one used in Section D and Theorem 6, exchanging $L\hat{f}_{n}$ by $L\circ\Phi\mu^{0}$. Following the arguments used in Section D, one may prove that there is no duality gap and the dual program is given by:

[TABLE]

We denote by $c_{\kappa}$ the solution to ($\mathbf{D}_{\kappa}(\Phi\mu^{0})$) (unicity can be proven by (50)) and by $\mathcal{P}_{\kappa}:=\Phi c_{\kappa}$ the dual polynomial. Its gradient is denoted by $\nabla\mathcal{P}_{\kappa}$, and its Hessian is denoted by $\nabla^{2}\mathcal{P}_{\kappa}$. We first state the next lemma.

Lemma 14**.**

If $c_{0}$ exists, then $\|c_{\kappa}-c_{0}\|_{\mathds{L}}\to 0$, $\nabla\mathcal{P}_{\kappa}\to\nabla\mathcal{P}_{0}$ uniformly, and $\nabla^{2}\mathcal{P}_{\kappa}\to\nabla^{2}\mathcal{P}_{0}$ uniformly as $\kappa\to 0$.

Proof.

Since $c_{\kappa}$ is a solution to ($\mathbf{D}_{\kappa}(\Phi\mu^{0})$), it holds that:

[TABLE]

leading to:

[TABLE]

and $c_{0}$ being a solution to ($\mathbf{D}_{0}$) implies that:

[TABLE]

We deduce that $\|c_{\kappa}\|_{\mathds{L}}\leq\|c_{0}\|_{\mathds{L}}$. Closed unit balls of Hilbert spaces being weakly sequentially compact we deduce that given $\kappa_{n}\to 0$, one shall extract a subsequence such that $c_{\kappa_{n}}$ weakly converges to some $c^{*}\in\mathds{L}$. Taking the limit as $\kappa\to 0$ in (54), we obtain that:

[TABLE]

Note that $\Phi c_{\kappa_{n}}$ converges weakly to $\Phi c^{*}$ so that:

[TABLE]

We deduce that $c^{*}\in\mathds{L}$ is a solution to ($\mathbf{D}_{0}$) and hence:

[TABLE]

Furthermore, $c^{*}$ is the solution of minimal norm since:

[TABLE]

The solution of minimal norm is unique by strict coercivity of the norm $\|\cdot\|_{\mathds{L}}$, see (50). We deduce that $c^{*}=c_{0}$, $\|c_{\kappa_{n}}\|_{\mathds{L}}\to\|c_{0}\|_{\mathds{L}}$, and $c_{\kappa_{n}}\to c_{0}$ strongly in ${\mathds{L}}$. Note that it implies that $\lim_{\kappa\to 0}\|c_{\kappa}-c_{0}\|_{\mathds{L}}=0$, since otherwise one can extract a subsequence $c_{\kappa_{n}}$ such that $\|c_{\kappa_{n}}-c_{0}\|_{\mathds{L}}>\varepsilon$, and by the above argument, one can extract a sequence such that $c_{\kappa_{n}}\to c_{0}$.

Now, the Cauchy-Schwarz inequality yields:

[TABLE]

which proves the uniform convergence. The same computation gives the uniform convergence of the functions and their gradients. ∎

We denote by $c_{\kappa,n}$ the dual solution of ($\mathbf{P}_{\kappa}$) , namely:

[TABLE]

and $\mathcal{P}_{\kappa,n}=\Phi c_{\kappa,n}$. The primal solution is denoted by $\hat{\mu}_{n}$.

Lemma 15**.**

Let $0<t\leq r$ and assume (NDB). If $\kappa$ and $\|\Gamma_{n}\|_{\mathds{L}}/\kappa$ are sufficiently small, any solution $\hat{\mu}_{n}$ has support of size $\hat{K}=K$ with one and only one spike in each near region $\mathds{N}_{k}(t)$ for $k\in[K]$.

Proof.

Note that ($\mathbf{D}_{\kappa}(\hat{f}_{n})$) and ($\mathbf{D}_{\kappa}(\Phi\mu^{0})$) are projection onto a closed convex set. We deduce that

[TABLE]

and that $\|\nabla^{2}\mathcal{P}_{\kappa}-\nabla^{2}\mathcal{P}_{\kappa,n}\|_{\infty}=\mathcal{O}(\frac{\|\Gamma_{n}\|_{\mathds{L}}}{\kappa})$ (the same result holds for the functions and their gradients). Under (NDB), we know that there exists $0<q<1$, $r>0$ and $\rho>0$ such that $\nabla^{2}\mathcal{P}_{0}\prec-\rho\mathrm{Id}_{d}$ on $\mathds{N}(r)$ and $|\mathcal{P}_{0}|<1-q$ on $\mathds{F}(r)$. We deduce that, for sufficiently small $\kappa$ and ${\|\Gamma_{n}\|_{\mathds{L}}}/\kappa$, $\mathcal{P}_{\kappa,n}$ is such that $\nabla^{2}\mathcal{P}_{\kappa,n}\prec-(\rho/2)\mathrm{Id}_{d}$ on $\mathds{N}(r)$ and $|\mathcal{P}_{\kappa,n}|<1-q/2$ on $\mathds{F}(r)$. We deduce that at most $1$ point in each $\mathds{N}_{k}(r)$ is such that $\mathcal{P}_{\kappa,n}(\hat{t}_{k})=1$.

But, since $\mu^{0}$ is the unique solution of ($\mathbf{P}_{0}$) (see Theorem 3), we deduce that $\hat{\mu}_{n}$ converges to $\mu^{0}$ in the weak-*topology as $\kappa$ and ${\|\Gamma_{n}\|_{\mathds{L}}}/\kappa$ go to zero. Hence, it holds that $\hat{\mu}_{n}(\mathds{N}_{k}(r))\to\mu^{0}(\mathds{N}_{k}(r))=a_{k}^{0}$. In particular, $\hat{\mu}_{n}$ has one spike in $\mathds{N}_{k}(r)$.

Now, by Taylor’s theorem, observe that if (NDB) with neighborhood size $r$ holds then it holds with neighborhood size $t$. ∎

It remains to bound $\|\Gamma_{n}\|_{\mathds{L}}^{2}$, which is the purpose of the next lemma.

Lemma 16**.**

A large enough universal constant $C>0$ exists such that for any RKHS $\mathds{L}$ associated to a nonnegative measure $\Lambda$:

[TABLE]

with probability at least $1-\alpha$. Or equivalently

[TABLE]

where $v_{n}:=\frac{C^{2}\Lambda(\mathds{R}^{d})}{n}=\frac{C^{2}\lambda(0)}{n}$.

Proof.

Let $X$ be a random variable with density $f^{0}$, we observe that $\mathds{E}_{X}L\delta_{X}=Lf^{0}$, and denote by $(Z_{i})_{i\in[n]}$ the i.i.d. random variables:

[TABLE]

which are i.i.d. centered random variables with values in $\mathds{L}$. Observe that $\|L\delta_{X_{i}}\|_{\mathds{L}}^{2}=\lambda(0)=\Lambda(\mathds{R}^{d})$ by the representation property of RKHS and the definition of its spectral measure $\Lambda$. We deduce that

[TABLE]

Using this inequality it holds that

[TABLE]

Now, consider the kernel $h(X_{i},X_{j})=\langle Z_{i},Z_{j}\rangle_{\mathbb{L}}$ and observe that the latter right hand side is a $U$-process. First, the Cauchy–Schwarz inequality and (55) lead to $\|h\|_{\infty}\leq 2\Lambda(\mathds{R}^{d})$. Second, check that this kernel is $\sigma$-canonical, namely:

[TABLE]

By Proposition 2.3 of [1], it follows that there exists two universal constants $C_{1},C_{2}>0$ such that

[TABLE]

with probability at least $1-\alpha$. ∎

Let $\delta_{\kappa}>0$ be arbitrarily small. Set $\kappa=\kappa_{n}=\sqrt{\Lambda(\mathds{R}^{d})}\,n^{-\frac{1}{2}+\delta_{\kappa}}$ so that

[TABLE]

with probability greater than $1-e_{n}:=1-Ce^{-n^{\delta_{\kappa}}}$. In this case, with an overwhelming probability, the requirements of the aforementioned Lemma 2 are met: $\kappa_{n}$ and $\|\Gamma_{n}\|_{\mathds{L}}\kappa_{n}^{-1}$ are small enough. Hence, a sequence of probability events $(e_{n})_{n\geq 1}$ exists such that $\lim_{n\rightarrow+\infty}e_{n}=0$ (almost) exponentially fast and for which the desired result holds (with $\delta_{\kappa}=1/2-r_{\kappa}$). It ends the proof of Theorem 11. ∎

Appendix F Construction of a dual certificate (proof of Theorem 7)

For a given set of points $S^{0}=\{t_{1},\ldots,t_{K}\}$, we recall that $\Delta:=\min_{k\neq\ell}\|t_{k}-t_{\ell}\|_{2}$. For any $\alpha\in\mathds{R}^{K}$ and $\beta\in\mathds{R}^{Kd}$, we consider the function

[TABLE]

For the sake of convenience, we omit the dependency in $\alpha$ and $\beta$ of the previous function and simply denote it by $p_{m}$. We are interested in the existence of a set of coefficients $(\alpha,\beta)$ such that $p_{m}$ defined in (56) satisfies the two interpolation conditions:

[TABLE]

The following proposition establishes the control of $\alpha$ and $\beta$ due to the conditions (57).

Proposition 17**.**

If $m$ is chosen such that $m\geq\frac{K^{1/4}d^{3/4}}{\mathcal{C}\Delta}$ for some positive constant $\mathcal{C}$ small enough, then $(\alpha,\beta)$ exists such that (57) holds and:

• •

$i)$* The supremum norms are upper bounded by:*

[TABLE]

• •

$ii)$* The Euclidean norm is upper bounded by:*

[TABLE]

Even though not directly usable in our framework, we emphasize that the stability result and the construction given in [9] played a central role in our work to prove Proposition 17.

Proof.

The proofs of $i)$ and $ii)$ are divided into four steps.

Step 1: Matricial formulation of (57).

The certificate $p_{m}$ should satisfy the following properties:

[TABLE]

We can organize the above equations to obtain a linear system of $K(d+1)$ equations with $K(d+1)$ parameters. In the following, we denote these parameters by $\alpha=(\alpha_{1},\ldots,\alpha_{K})^{T}\in\mathds{R}^{d}$ and $\beta=(\beta_{1}^{1},\ldots,\beta_{1}^{d},\beta_{2}^{1},\ldots,\beta_{2}^{d},\ldots,\beta_{K}^{1},\ldots,\beta_{K}^{d})^{T}\in\mathds{R}^{Kd}$. The above equations can be rewritten as:

[TABLE]

where $\mathbf{1}_{K}$ denotes the vector of size $K$ having all its entries equal to $1$ (similar definition for $\mathbf{0}_{Kd}$). The matrix $A_{m}\in\mathds{R}^{K\times K}$ acts on the coefficients $\alpha$ as:

[TABLE]

while $D_{m}\in\mathds{R}^{K\times Kd}$ describes the effect of the partial derivatives of $\psi_{m}$ on $\alpha$ and $\beta$ as:

[TABLE]

Finally, the squared matrix $B_{m}\in\mathds{R}^{Kd\times Kd}$ is given by:

[TABLE]

Step 2: Inversion of the system (60)

According to linear algebra results (see e.g. [23]), the system (60) is invertible if and only if

[TABLE]

and its Schur complement

[TABLE]

are both invertible. To prove this assertion, we remember that a symmetric matrix $M$ is invertible if $\|I-M\|_{\infty}<1$, where $\|.\|_{\infty}$ is the subordinate matrix infinity norm ($\|M\|_{\infty}=\underset{i}{\max}\sum_{j}|M_{ij}|$). In such a case $\|M^{-1}\|_{\infty}\leq\frac{1}{1-\|I-M\|_{\infty}}$.

Moreover, we will use in the sequel that $\|M\|_{1}=\underset{j}{\max}\sum_{i}|M_{ij}|=\|M^{T}\|_{\infty}$.

$\bullet$ Invertibility of $G_{m}$ and computation of $\|G_{m}^{-1}\|_{\infty}$:

For all $i,k\in[K]$ and $u,v\in[d]$,

[TABLE]

according to the definition of $B_{m}$. Setting $\tilde{G}_{m}=\frac{-3}{4m^{2}}G_{m}$, we get

[TABLE]

with

[TABLE]

according to Lemma 19 and $iii)$ of Lemma 20. Thus, if there exists a positive constant $\mathcal{C}$ small enough such that $\frac{Kd^{3}}{m^{4}\Delta^{4}}\leq\mathcal{C}$, then $\|I_{Kd}-\tilde{G}_{m}\|_{\infty}<1/2$ and the matrix $G_{m}$ is invertible. Moreover,

[TABLE]

$\bullet$ Invertibility of $H_{m}$ and computation of $\|H_{m}^{-1}\|_{\infty}$:

In the same way, we want to prove that $\|I_{K}-H_{m}\|_{\infty}<1$. According to the properties of the $\infty$-norm,

[TABLE]

In a first time, we provide an upper bound on $\|A_{m}\|_{\infty}$. Remark that for any $i\in[K]$

[TABLE]

Applying $i)$ of Lemma 20, we hence obtain

[TABLE]

Now, recall that

[TABLE]

Applying $ii)$ of Lemma 20, we deduce that:

[TABLE]

Following the same ideas, for any pair $(i,u)$ with $i\in[K]$ and $u\in[d]$, we have:

[TABLE]

Again, $ii)$ of Lemma 20 yields:

[TABLE]

Gathering (61), (63), (64) and (65) in (62), we deduce that:

[TABLE]

provided $\frac{Kd^{3}}{m^{4}\Delta^{4}}\leq\mathcal{C}$ for some constant $\mathcal{C}$ small enough. This implies that under such a condition, $\|I_{K}-H_{m}\|_{\infty}<1/2$. Moreover, the Schur complement $H_{m}=(I_{K}+A_{m})-D_{m}G_{m}^{-1}D_{m}^{T}$ is then invertible and

[TABLE]

for some positive constant $C$, provided the constraint $\frac{Kd^{3}}{m^{4}\Delta^{4}}\leq\mathcal{C}$ is satisfied.

To conclude this second step, the system (60) is invertible if $\frac{Kd^{3}}{\Delta^{4}m^{4}}\leq\mathcal{C}$ for some constant $\mathcal{C}$ small enough. In such a case

[TABLE]

**Step 3: Proof of i) ** In the sequel, we assume that $\frac{Kd^{3}}{m^{4}\Delta^{4}}\leq\mathcal{C}$ for some positive constant $\mathcal{C}$ small enough. First, according to (67) and (68), we obtain that:

[TABLE]

Moreover

[TABLE]

with $\tilde{H}_{m}=A_{m}-D_{m}G_{m}^{-1}D_{m}^{T}$. Hence, since for $\tilde{H}_{m}$ small enough (i.e. for a sufficiently small norm) we have

[TABLE]

Hence, for $\tilde{H}_{m}$ small enough,

[TABLE]

According to (66),

[TABLE]

and we can choose the constant $\mathcal{C}$ small enough in the constraint $\frac{Kd^{3}}{m^{4}\Delta^{4}}\leq\mathcal{C}$ such that $\|\tilde{H}_{m}\|_{\infty}\leq 1/2$. We conclude that

[TABLE]

In a second time, gathering (61), (65) and (67), we deduce that:

[TABLE]

**Step 4: Proof of ii) ** According to (68),

[TABLE]

Using (71), $\|G_{m}^{-1}D_{m}^{T}H_{m}^{-1}\|_{\infty}\leq\|G_{m}^{-1}\|_{\infty}\|D_{m}\|_{1}\|H_{m}^{-1}\|_{\infty}\lesssim\frac{\mathcal{C}}{md}$. For the second term, we use the dual relationship between $\|.\|_{\infty}$ and $\|.\|_{1}$ and that the matrices $G_{m}$ and $H_{m}$ are symmetric. Gathering (61), (64) and (67), we obtain that:

[TABLE]

We then deduce that

[TABLE]

∎

Thanks to the previous proposition, we are now ready to prove Theorem 7. Our strategy is inspired from the one of [9].

Proof of Theorem 7.

We define an integer $m$ that will be chosen large enough below and consider $\mathcal{P}_{m}=p_{m}^{2}$.

Proof of $i)$ and $ii)$: From Proposition 17, we know that if $m$ satisfies $m\geq\mathcal{C}\frac{K^{1/4}d^{3/4}}{\Delta},$ for a constant $\mathcal{C}$ large enough independent from $K$, $\Delta$ and $d$, then $\mathcal{P}_{m}$ satisfies the interpolation properties:

[TABLE]

Our strategy relies on a study of the variations of $\mathcal{P}_{m}$ near each support points $\{t_{1},\ldots,t_{K}\}$, whose union defines the near region, and far from these support points, which is then the far region.

Near region

Let $\epsilon\in\left(0,\frac{\Delta}{2}\right)$ a parameter whose value will be made precise later on. The near-region $\mathds{N}(\epsilon)$ is the union of $K$ sets that are defined by:

[TABLE]

The far region is therefore given by:

[TABLE]

Let $i\in\{1,\dots,K\}$ be fixed, the function $p_{m}$ involves a sum over $k\in[K]$ and we consider two cases:

• •

If $k\not=i$, then, for all $t\in\mathds{N}_{i}(\epsilon)$, $\xi_{t,i,k}$ exists such that

[TABLE]

with $\|\xi_{t,i,k}-t_{i}\|_{2}\leq\|t-t_{i}\|_{2}$. Moreover, for any $u\in[d]$, a $\tilde{\xi}_{t,i,k}^{u}$ exists such that:

[TABLE]

with $\|\tilde{\xi}_{t,i,k}^{u}-t_{i}\|_{2}\leq\|t-t_{i}\|_{2}$.

• •

If $k=i$, since $\nabla\psi_{m}(0)=0$ and $D^{3}(\psi_{m})(0)=0$, for all $t\in\mathds{N}_{i}(\epsilon)$, a $\xi_{t,i,i}$ exists such that:

[TABLE]

with $\|\xi_{t,i,i}-t_{i}\|_{2}\leq\|t-t_{i}\|_{2}$. We also have that for any $u\in\{1,\dots,d\}$, the existence of $\tilde{\xi}_{t,i,i}^{u}$ such that:

[TABLE]

with $\|\tilde{\xi}_{t,i,i}^{u}-t_{i}\|_{2}\leq\|t-t_{i}\|_{2}$.

Hence, for all $t\in\mathds{N}_{i}(\epsilon)$, we can use the previous Taylor formulas and obtain that:

[TABLE]

These terms can be re-arranged as follows:

[TABLE]

Of course, the construction of Proposition 17 yields

[TABLE]

and

[TABLE]

thanks to the constraints expressed on the function $p_{m}$. Hence, for all $t\in\mathds{N}_{i}(\epsilon)$ we have

[TABLE]

In the following, we prove that $C_{2}$ is a negative matrix and bounded from below. Thanks to Lemma 20, we can compute the first term of $C_{2}$ and we have

[TABLE]

which entails

[TABLE]

The second term of $C_{2}$ may be upper bounded with the help of the spectral radius of $D^{2}(\psi_{m})(\xi_{t,i,k}-t_{k})$: (denoted by $\rho(M)$ for any squared symmetric matrix $M$):

[TABLE]

To handle this last term, we use the fact that in the near region $\mathds{N}_{i}(\epsilon)$, $\|\xi_{t,i,k}-t_{k}\|_{2}$ is far from [math]. Using the triangle inequality, since $\epsilon<\frac{\Delta}{2}$, we have for any $k\in[K]$ with $k\not=i$

[TABLE]

Using the inequality $\rho(M)\leq\|M\|_{\infty}$ for any symmetric matrix, Proposition 17 and $iii)$ of Lemma 20, we obtain that:

[TABLE]

The third term of $C_{2}$ is described by the matrix:

[TABLE]

Using that $\|\text{sinc}^{\prime}\|_{\infty}\vee\|\text{sinc}^{(2)}\|_{\infty}\vee\|\text{sinc}^{(3)}\|_{\infty}\vee\|\text{sinc}^{(4)}\|_{\infty}\leq 1/2$, we obtain that $\partial_{u,v,p,q}\psi_{m}(\xi_{t,i,i}-t_{i})\lesssim m^{4}$. Therefore, for any $(u,v)\in[d]^{2}$, we have:

[TABLE]

where the last line comes from the Cauchy-Schwarz inequality. Again, the inequality $\rho(M)\leq\|M\|_{\infty}$ and Proposition 17 yield:

[TABLE]

The last term of $C_{2}$ is studied into two steps. We first consider the situation when $k\neq i$: the triangle inequality, $iv)$ of Lemma 20 and the inequality $\rho(M)\leq\|M\|_{\infty}$ yield:

[TABLE]

Hence, we deduce from Proposition 17 that:

[TABLE]

Now, we consider the situation where $k=i$. For any pair $(u,v)\in[d]^{2}$:

[TABLE]

where we used $iv)$ of Lemma 19, $\nabla\psi_{m}(0)=0$, $D^{3}\psi_{m}(0)=0$ and $\|m(\tilde{\xi}_{t,i,i}^{u}-t_{i})\|_{2}\leq m\epsilon$ and $i)$ of Proposition 17.

Using the previous bounds, we then conclude that

[TABLE]

provided $\frac{Kd^{3}}{m^{4}\Delta^{4}}\leq\mathcal{C}$ for a constant $\mathcal{C}$ small enough.

We now plug Equations (74), (75), (76), (77) and (78) in $C_{2}(t)$ and deduce that a constant $\square$ exists such that

[TABLE]

Then, we choose $\epsilon$ and $m$ such that $A_{\epsilon,m}\leq-\frac{\alpha_{i}}{3}$. A careful inspection of the above terms prove that a sufficiently small $\upsilon$ and large enough $C$ (both independent of $d$, $K$ and $\Delta$) exist such that

[TABLE]

Far region $\mathds{F}(\epsilon)$

The relationship between $\epsilon,m$ and $d$ being established in (79), we are looking for a value of $\eta>0$ such that

[TABLE]

The definition of $p_{m}$ and the Cauchy-Schwarz inequality yield

[TABLE]

We consider the second term of the right hand side with the help of Lemma 19 and Proposition 17:

[TABLE]

In particular, there exists a constant $\check{C}$ such that

[TABLE]

Let $t\in\mathds{F}(\epsilon)$ and $t_{i^{\star}}$ the closest point of $t$ in the set $\{t_{1},\dots,t_{K}\}$, the triangle inequality shows that $\forall k\not=i^{\star}$, we have $\|t-t_{k}\|_{2}>\frac{\Delta}{2}.$ Hence, since $\|\alpha\|_{\infty}$ is upper bounded by a universal constant (see Proposition 17), we deduce from $i)$ Lemma 20 that

[TABLE]

In the same time, the last term that involves $i^{\star}$ is upper bounded by

[TABLE]

where $C_{0}$ is a large enough universal constant. Using that

[TABLE]

and the fact that when $\|y\|_{2}\geq\upsilon d^{-1}$, then the absolute value of one of the coordinate of $y$ is greater than $\upsilon d^{-3/2}$, we deduce that

[TABLE]

where $\eta\asymp\upsilon^{2}d^{-3}$. This entails the desired result as soon as $m$ is chosen such that

[TABLE]

for some positive constant $\mathcal{C}$ small enough. It is easy to check that in this case, a small enough $\upsilon$ exists (independent of $d$, $K$, $m$ and $\Delta$) such that:

[TABLE]

Conclusion of the interpolation

To accomodate with conditions (79) and (81), we consider an integer $m$ such that $m\gtrsim K^{1/2}d^{3/2}\Delta^{-1}$ and $\epsilon=\upsilon m^{-1}d^{-1}$. We deduce that $p_{m}$ satisfies in the far region $\mathds{F}(\epsilon)$:

[TABLE]

while in the near region we have:

[TABLE]

We then set $\mathcal{P}_{m}=p_{m}^{2}$. This function satisfies both the constraints and the interpolation conditions in the statement of Theorem 7. We then obtain $i)$ and $ii)$.

Proof of $iii)$:

Remark first that $p_{m}$ is a linear combination of shifted sinus cardinal functions and derivatives of sinus cardinal functions up to the power $4$ used in $\psi$. Moreover, it is straightforward to check that

[TABLE]

Therefore, the Fourier transform of $\psi^{4}$ has a compact support of size $[-2,2]^{d}$ since the Fourier transform of the sinus cardinal is the rectangular indicator function of $[-1/2,1/2]$. Using the effect on the Fourier transform of scaling and shifting a function we deduce that the Fourier transform of $p_{m}$ has a compact support, which size varies linearly with $m$:

[TABLE]

Since $\mathcal{P}_{m}=p_{m}^{2}$, we have $\mathcal{F}[\mathcal{P}_{m}]=\mathcal{F}[p_{m}]\star\mathcal{F}[p_{m}]$ so that

[TABLE]

We now compute an upper bound of $\|\mathcal{P}_{m}\|_{2}$: the isometry property entails the several inequalities:

[TABLE]

where we used the standard inequality $\|g\star h\|_{2}\leq\|g\|_{2}\|h\|_{1}$.

Now, the triangle inequality yields

[TABLE]

where the last line comes from the Cauchy-Schwarz inequality.

We then deduce that

[TABLE]

where $|\mathcal{F}[\nabla\psi_{m}]|_{2}$ refers to the Euclidean norm of the $d$-dimensional vector $\mathcal{F}[\nabla\psi_{m}]$. Now, remark that a dilatation by a ratio $m$ yields on $L^{2}$ norms:

[TABLE]

Hence

[TABLE]

We use a similar argument and obtain that

[TABLE]

In the meantime, the effect of this dilatation on the $L^{1}$ norms is managed by:

[TABLE]

and with a same argument we obtain that: $\||\mathcal{F}[\nabla\psi_{m}]|_{2}\|_{1}\lesssim d.$ Hence

[TABLE]

We then obtain that

[TABLE]

provided

[TABLE]

for some constant $\mathcal{C}$ small enough.

Proof of $iv)$: The last point is a simple consequence of the convolution kernel induced by $\Phi$. Since $\varphi$ satisfies $(\mathcal{H}_{4m})$, then $\forall\xi\in[-4m,4m]^{d}$, we have $\sigma(\xi)\neq 0$. Hence, we can define $c_{0,m}$ through its Fourier transform:

[TABLE]

Moreover, the Fourier transform of $c_{0,m}$ is naturally compact, which entails that $c_{0,m}\in\mathds{L}$.

Conclusion: The constraints (80) and (82) together with $\epsilon\sim\frac{1}{md}$ and $\epsilon\leq\Delta/2$ leads to the condition

[TABLE]

Provided $\Delta$ is small or bounded, this condition reduces to $m\gtrsim\frac{K^{1/2}d^{3/2}}{\Delta}$ and $m\gtrsim\sqrt{Kd}$ when $\Delta$ is large. ∎

Some useful properties of the sinus cardinal function are detailed in the following basic lemma.

Lemma 18**.**

If $g(x)=\text{sinc}(x)$, then for any $x\in\mathds{R}$:

• $i)$

[TABLE]

• $ii)$

[TABLE]

• $iii)$

[TABLE]

• $iii)$

[TABLE]

Some additional ingredients on $\psi_{m}$ are detailed below where $\psi_{m}$ is defined by:

[TABLE]

In the sequel, we will use the shortcut of notation $\partial_{u}$ instead of $\partial^{|u|}_{u}\psi$ for any multi-index $u$.

Lemma 19**.**

Let $\psi_{m}$ be the function defined in (84). Then

• •

$i)$* $\psi_{m}(0)=g^{4}(0)^{d}=1$.*

• •

$ii)$* $\nabla\psi_{m}(0)=0$ and*

[TABLE]

• •

$iii)$* $D^{2}\psi_{m}(0)=-\frac{4}{3}m^{2}I_{d}$ and*

[TABLE]

• •

$iv)$* $(D^{3}\psi_{m})(0)=0$ and*

[TABLE]

• •

$v)$* Finally*

[TABLE]

with

[TABLE]

[TABLE]

and

[TABLE]

Several bounds on the successive derivatives of $\psi_{m}$ are given in the following lemma.

Lemma 20**.**

For any pair $(i,j)$ such that $i\neq j$:

• •

$i)$* $|\psi_{m}(t_{i}-t_{j})|\lesssim\frac{d^{2}}{m^{4}\Delta^{4}}.$*

• •

$ii)$* $|\partial_{u}\psi_{m}(t_{i}-t_{j})|\lesssim\frac{d^{2}}{m^{3}\Delta^{4}}$*

• •

$iii)$* $|\partial^{2}_{u,v}\psi_{m}(t_{i}-t_{j})|\lesssim\frac{d^{2}}{m^{2}\Delta^{4}}.$*

• •

$iv)$* $|\partial^{3}_{u,v,w}\psi_{m}(t_{i}-t_{j})|\lesssim\frac{d^{2}}{m\Delta^{4}}.$*

Proof.

In what follows, we deliberately choose to omit the multiplicative constants since the rest of the paragraph will be managed in the same way.

Point $i)$: we use $|\text{sinc}(x)|\leq|x|^{-1}$ and remark that $\|t_{i}-t_{j}\|_{2}\geq\Delta$ so that

[TABLE]

We then deduce that

[TABLE]

because one coordinate $\ell_{0}$ exists such that $|t_{i}^{\ell_{0}}-t_{j}^{\ell_{0}}|^{2}\geq\Delta^{2}d^{-1}$.

Point $ii)$: we use Lemma 18, Lemma 19 and

[TABLE]

associated with $|g(x)|\vee|g^{\prime}(x)|\lesssim\frac{1}{|x|}$. It yields

[TABLE]

We then obtain $ii)$.

Point $iii)$: we still use Lemma 18 and Lemma 19, the fact that

[TABLE]

and $|g(x)|\vee|g^{\prime}(x)|\vee|g^{\prime\prime}(x)|\lesssim\frac{1}{|x|}$. It leads to

[TABLE]

Point $iv)$: the proof follows the same lines with the help of the previous lemmas, we check that

[TABLE]

∎