Blanche Buet

Antonin Chambolle

Simon Masnou

Gabriel Peyré

Keywords: inverse problems, total variation, sparsity. C

À Vincent qui a eu bien du courage de s'embarquer dans cette galère.

À Yohann pour m'avoir suivi depuis Lyon (et pour m'avoir poussé dans le bureau de Vincent).

À Marc Dambrine et Otmar Scherzer pour leurs précieuses relectures.

R

On s'intéresse dans cette thèse à une famille de problèmes inverses, qui consistent à reconstruire une image à partir de mesures partielles et possiblement bruitées. Si la résolution de problèmes inverses mal posés est une thématique classique en traitement du signal, une série de travaux menés depuis une trentaine d'années a particulièrement marqué le domaine. Ces derniers concernent la reconstruction de signaux dits parcimonieux à l'aide de méthodes variationnelles, qui reposent sur l'utilisation d'une fonctionnelle, appelée régulariseur, dont la minimisation produit des solutions de même structure que le signal inconnu.

La principale contribution de cette thèse est l'étude d'un choix particulier de régulariseur, la variation totale (du gradient). Cette fonctionnelle est utilisée en imagerie depuis les travaux de Rudin Osher et Fatemi, menés en 1992, notamment pour sa capacité à pénaliser les oscillations tout en préservant les discontinuités. Alors qu'il est bien connu que sa minimisation produit des images constantes par morceaux, présentant une forme de parcimonie (elles sont composées d'un petit nombre de formes simples), ce point de vue n'a à notre connaissance pas été privilégié pour analyser les performances de ce régulariseur.

Dans cette thèse, on se propose de mener cette étude. Dans la partie II, on considère les reconstructions obtenues par minimisation de la variation totale dans un régime de faible bruit, et on étudie leur proximité avec l'image inconnue. Cette dernière étant supposée parcimonieuse (c'est à dire composée d'un petit nombre de formes simples), on s'intéresse particulièrement à la structure de la reconstruction: est-elle elle-même parcimonieuse, est-elle composée du même nombre de formes, et ces formes sont-elles proches de celles présentes dans l'image inconnue. Dans la partie III, on propose une méthode numérique pour résoudre les problèmes variationnels associés à ce régulariseur. On adapte des travaux récents sur la reconstruction de mesures discrètes, a n d'obtenir un algorithme ne reposant pas sur l'introduction d'une discrétisation spatiale xe. Ceci a l'avantage, contrairement aux techniques existantes, de n'introduire ni ou ni anisotropie dans les images reconstruites, et d'en produire une représentation parcimonieuse.

Les contributions de cette thèse peuvent aussi être considérées du point de vue du calcul des variations. On s'intéresse en e et à une problème variationnel, qu'on analyse théoriquement et pour lequel on propose une méthode numérique de résolution. Dans la partie II, on s'intéresse à la convergence des solutions de ce problème lorsque l'un des paramètres le dé nissant tend vers zéro. Le type de convergence qu'on obtient est bien plus fort que ceux considérés usuellement en calcul des variations. Plus précisément, on montre la convergence des ensembles de niveau des solutions en termes de déformations normales, ainsi que celle du nombre d'ensemble de niveaux non triviaux (et de leur nombre de composantes connexes). Dans la partie III, on résout ce problème numériquement sans se ramener à un problème de dimension nie via l'introduction d'une discrétisation spatiale xe. On propose une méthode itérative qui construit des solutions approximées ayant la même structure que certaines solutions d'intérêt.

Mots-clés: problèmes inverses, variation totale, parcimonie.

C

Part I. We begin this manuscript with a general introduction to inverse problems and variational regularization techniques. We then introduce the total (gradient) variation, which is the speci c regularizer this thesis is dedicated to. We review existing theoretical results on total variation regularization, and discuss numerical methods for solving the associated variational problems. We argue that these subjects have not extensively been studied from a sparse recovery viewpoint, and state that the main aim of this thesis is to conduct such an analysis. We also stress that our contributions are attempts to adapt recently introduced ideas for sparse spikes recovery in a continuous domain.

Part II. This part is devoted to the noise robustness analysis of total variation regularization, and more speci cally to the proof of an exact support recovery result.

• Piecewise constant functions. The main purpose of this section is to introduce the class of sparse objects we aim at recovering, which are piecewise constant functions with a few speci c properties. We call them M -simple functions. We rst begin by collecting results on the faces of the total variation unit ball, which were rst proved in [Duval, 2022]. As we focus on exposed faces (instead of general linearly closed faces in the above-mentioned work), we are able to provide slightly di erent (and often simpler) proofs of some results. It turns out that the elements of nite dimensional exposed faces are M -simple functions. This analysis in particular allows to decompose M -simple functions belonging to the same face using a common set of atoms, which is crucial to prove our support recovery result. As a byproduct, we are also able to state an abstract su cient identi ability condition, which can be seen as a strengthened source condition.

• The prescribed curvature problem. We then focus on a geometric variational problem called the prescribed curvature problem, which naturally appears in the analysis conducted in the previous section. Our aim is to study the behaviour of its solutions under variations of the associated curvature functional. We prove that, for two su ciently close curvature functionals, solutions of the two problems are close in terms of C 2 normal deformations.

Assuming some solution of one problem is strictly stable, we show that it has a neighborhood (in terms of C 2 normal deformations) which contains at most one solution of the other problem.

• Exact support recovery. Using the results of the rst two sections, we are nally able to prove our support recovery result. To achieve this, we rst show that the dimension of the exposed faces of the total variation unit ball is stable in some sense. We then de ne a so-called non-degenerate source condition, under which exact support recovery is guaranteed.

7

We stress that the term exact refers to the estimation of some measure of sparsity of the unknown image. The reconstructed image is made of the same number of shapes, each a C 2 normal deformation of a shape appearing in the sougt-after image.

Part III. This third part is dedicated to the numerical resolution of variational problems associated to total variation regularization. We propose a grid-free approach based on the Frank-Wolfe (or conditional gradient) algorithm. A complete implementation of the proposed method is available at https://github.com/rpetit/pycheeger and https://github.com/rpetit/tvsfw.

• Frank-Wolfe approach. In this rst section, we describe the theoretical algorithm we propose to solve our variational problem of interest. We explain how the linear minimization step, which is at the core of Frank-Wolfe algorithm, is closely linked to a shape optimization problem called the (generalized) Cheeger problem. We then turn to convergence results, before discussing the interest of a special nal update, called the sliding step.

• Polygonal approximation of generalized Cheeger sets. This section is devoted to the numerical approximation of solutions of the (generalized) Cheeger problem. We propose to optimize its objective over simple polgons with a xed number of sides, and prove the existence of optimizers for this new problem. We then discuss a two-step numerical method to approach them. Finally, we discuss how polygonal maximizers compare to their continuous counterpart in the speci c case of radial weight functions.

• Numerical results. We begin this last section by discussing the implementation of the sliding step. We describe the iterative method we propose to locally optimize the objective, and also mention the issue of topology changes, which might appear over the course of this evolution. We conclude this part by assessing the performance of our algorithm on a few recovery examples. We provide comparisons with two standard xed-grid methods, and discuss parameter choices.

Part IV. We conclude this manuscript by discussing interesting avenues for future research.

Publications. This thesis gave rise to two articles, which are listed below.

• [De Castro et al., in preparation] Exact recovery of the support of piecewise constant images via total variation regularization, Y. De Castro, V. Duval and R. Petit, in preparation.

• [De Castro et al., 2022] Towards o -the-grid algorithms for total variation regularized inverse problems, Y. De Castro, V. Duval and R. Petit, Journal of Mathematical Imaging and Vision, 2022

A short version of the second article appeared in the proceedings of the 8th international conference on scale space and variational methods in computer vision (SSVM 2021) [De Castro et al., 2021].

(Ω, R n ) (C k b (Ω) if n = 1) the set of functions u ∈ C k (Ω, R n ) that satisfy u C k (Ω) < +∞ with u C k (Ω) def. = k i=1 D i u ∞ ,
where D i u denotes the i-th derivative of u. If a function u is bounded and uniformly continuous on Ω, it has a unique continuous extension to Ω. We denote by C k (Ω) the space of functions in C k (Ω) with bounded and uniformly continuous derivatives up to order k.

Hausdor measure. We denote by H 1 the 1-dimensional Hausdor measure on R 2 , and, for every Borel set A ⊂ R 2 , by H 1 A the measure H 1 restricted to A, i.e. such that for every Borel set E we have

H 1 A (E) = H 1 (A ∩ E) .

I

Broadly speaking, an inverse problem consists in reconstructing an unknown object from partial and often indirect observations of it. Many problems fall in this category. To name only a few, let us mention the reconstruction of a 3D object from a few 2D photographs, the location of an earthquake's epicenter from measurements of the resulting seismic waves, or the determination of the composition of a material from the way it interacts with light (or the way other kinds of waves travel through it). This thesis is dedicated to a particular family of inverse problems, which belongs to the class of imaging inverse problems, in which the unknown object to reconstruct is an image. In this section, we rst provide a broad introduction to this domain1 , before presenting the problem we are speci cally interested in. Then, we introduce a family of methods for solving inverse problems, namely variational regularization techniques, that are at the core of our work. Finally, we present the main aim of this thesis, which is to adapt recent advances in the eld to our speci c setting.

I

As mentioned above, an imaging inverse problem consists in reconstructing an image (in the general sense) from a set of measurements. Such problems are ubiquitous in biological, medical, and astronomical imaging, as well as in computational photography. Before discussing these applications in greater detail, let us introduce the mathematical framework commonly used to model these problems.

Framework and assumptions. The unknown image is modeled as an element f 0 of some vector space E. We assume that the measurement operator Φ, which accounts for the measurements acquisition process, is a linear 2 mapping from E to some Hilbert space H (a common scenario is to have access to m real-valued measurements, in which case H = R m ). Finally, the noiseless and noisy observations are respectively de ned by y 0 = Φf 0 and y = y 0 + w where w ∈ H is an additive noise. Solving the inverse problem means (approximately) recovering f from the knowledge of y 0 (noiseless case) or y (noisy case).

Modeling images. The choice of the vector space E in which the sought-after image lies depends on the application. When dealing with digital images (say grayscale images for simplicity), which are arrays of real values, it is natural to take E = R m×n for some pair of integers m, n ∈ N * , or simply E = R n for some n ∈ N * . On the other hand, if one wishes to model "continuous images", choosing E to be an in nite dimensional space of functions or measures is relevant. The main examples of interest for our purpose are E = M(Ω), the set of bounded Radon measures on Ω, and E = L 2 (Ω), where Ω is the image domain.

Structural assumptions on the sought-after image. In most cases, a major obstacle towards the resolution of an inverse problem is the non-injectivity of the measurement operator Φ. To put it another way, given a set of observations y 0 , multiple images f in E could satisfy Φf = y 0 , and the di culty is to select the right one among them. Even in situations where Φ is injective, the linear system Φf = y 0 is often ill-conditioned. Inverting Φ hence yields reconstructions that are highly sensitive to changes of y 0 , and which are therefore not robust to noise. The main ingredient usually used to overcome these issues is to make structural assumptions on the sought-after image f 0 . If we know a priori that f 0 has a special structure, then this information can be exploited to recover it among all images f satisfying Φf = y 0 . In the following subsection, we give examples of interesting structural assumptions on f 0 , which are all closely linked to the concept of sparsity.

Sparse signals

A class of signals that has been extensively studied since the 90's is the class of signals having some kind of sparsity. Given a set of atoms called a dictionary, a signal is said to be sparse if it can be decomposed as the weighted sum of a few atoms.

Sparse vectors. The rst example of sparse objects that were extensively studied are sparse vectors. The sparsity of a vector x of R n is de ned as its number of non-zero coordinates, and a vector is said to be sparse if its sparsity is signi cantly smaller than n (it can hence be written as a weighted sum of a few 1-sparse vectors). The problem of recovering a sparse vector from a set of linear measurements has a tremendous number of applications. It has been studied empirically from the 70's in di erent domains, including seismic imaging [START_REF] Claerbout | Robust modeling with erratic data[END_REF]. Its mathematical analysis started in the 90's with the work of David Donoho and his co-authors (see for instance [Donoho, 1992, Chen et al., 1998]). In some situations, it may happen than one is interested in recovering images directly represented by sparse vectors (for example, if one wishes to reconstruct a superposition of a few point light sources). However, it often happens that the sought-after image is not directly sparse, but admits a sparse decomposition in some dictionary. In other terms, there exists a sparse vector x 0 and a dictionary Ψ such that f 0 = Ψx 0 . In this case, in order to exploit the sparsity of x 0 for solving the considered inverse problem, one can de ne a new measurement operator A = ΦΨ, and solve the inverse problem associated to x 0 and A1 . This ultimately allows to obtain an approximation of f 0 . We refer the reader to [Mallat, 2008, Chapter 12 and 13] for more details on this subject.

Discrete measures. If one wishes to reconstruct a superposition of point light sources, it is often interesting, as explained in Section 1.3, to model the sought-after image as the weighted sum of a few Dirac masses. In this case, one choses E = M(X) the space of bounded Radon measures on a domain X (usual choices of X include the n-dimension torus T n or an open subset of R n ), and assumes that the unknown signal is a measure µ 0 of the form

µ 0 = s i=1 a i δ x i ,
with a ∈ R s and x ∈ X s . The recovery of such signals has been the subject of many works in recent years (see for instance [De Castro and Gamboa, 2012, Bredies and Pikkarainen, 2013, Candès and Fernandez-Granda, 2014]). We refer the reader to Section 1.3.2 below and the review article [START_REF] Laville | O -The-Grid Variational Sparse Spike Recovery: Methods and Algorithms[END_REF] for more details. A part of the work presented in this thesis is an attempt to adapt these results to the reconstruction of piecewise constant functions.

Time-dependant discrete measures. One could also consider the dynamic inverse problem which consists in recovering the evolution of a superposition of point light sources from a set of measurements. In this case the unknown signal is a time-dependent measure of the form

µ 0 : t → s i=1 a i (t)δ x i (t)
where a and x take values in R s and X s . This task was notably studied in the xed mass setting in [START_REF] Alberti | Dynamic Spike Superresolution and Applications to Ultrafast Ultrasound Imaging[END_REF], Bredies and Fanzon, 2020, Bredies et al., 2021], and then in the varying mass setting in [Bredies et al., 2022a] .

Measures supported on 1D curves. Another class of signals presenting some kind of sparsity are measures supported on a few 1D curves. Such measures could be relevant to model structures of interest in biological imaging. We provide in Figure 1 two images of microtubules to support this claim. In this context, one would chose E = M(X) but assume the sought-after signal is a measure µ 0 of the form

µ 0 = s i=1 a i H 1 Γ i ,
with a ∈ R s and Γ i a connected 1D curve for all i ∈ {1, ..., s}. Let us stress that a related task is considered in [START_REF] Laville | O -the-grid curve reconstruction through divergence regularisation: An extreme point result[END_REF]. This work however deals with vector-valued measures, which are hence associated (just as time-dependant measures) to oriented curves.

Piecewise constant functions. In this thesis, the signals we focus on are models for piecewise constant images (also called "cartoon images"). The recovery of such objects could have potential applications in e.g. cell imaging. We provide in Figure 2 a set of images which could be considered piecewise constant as a rst approximation. In this setting, one could choose E = L 2 (R 2 ) and assume the sought-after signal is an element u 0 of E of the form

u 0 = s i=1 a i 1 E i ,
with a ∈ R s and (E i ) s i=1 a collection of su ciently regular subsets of the plane. We give a precise de nition of the class of piecewise constant functions we deal with in Section 1 of Part 2.

Examples of inverse problems

We now give a few examples of measurement operators of interest.

Denoising. The choice of Φ = Id corresponds to the denoising setting, where one wishes to recover an image corrupted by some (additive) noise.

Deconvolution. Another setting of interest is when one has access to a blurry version of the sought-after image. This is modeled by de ning Φf = ϕ f with ϕ a convolution kernel.

The " " operation denotes a discrete or continuous convolution depending on the nature of f . The kernel ϕ could for example model the point spread function of an imaging device, or a motion blur.

Fourier transform sampling. The reconstruction of an image from a sampling of its Fourier transform has many practical applications. We describe two of them below. The rst is the detection of point light sources using radio interferometry in astronomical imaging. The second is medical imaging, and more precisely X-ray tomography and magnetic resonance imaging.

• Radio interferometry: if one wishes to observe point light sources in the sky, a natural model for the sought-after image is a measure µ 0 that is a combination of a few Dirac masses. In radio interferometry, one uses an array of antennas to record the light emitted by celestial bodies. Computing the cross-correlation between these recordings gives access to a sampling of the Fourier transform of µ 0 (see e.g. [START_REF] Pan | Towards Generalized FRI Sampling With an Application to Source Resolution in Radioastronomy[END_REF]). The measurement operator can hence be modeled as

Φ : M(R n ) → C |Ω| µ → ˆRd exp(-i ω, x ) dµ(x) ω∈Ω ,
where Ω ⊂ R n is a set of observed frequencies, and d is typically equal to two (e.g. in the context of astronomical imaging).

• X-ray tomography, magnetic resonance imaging: other situations in which one wishes to recover an image from a sampling of its Fourier transform include X-ray tomographgy and magnetic resonance imaging (MRI). Here, the aim is to to reconstruct a general digi-

tal (E = R n ) or continuous (e.g. E = L 2 (R 2 )) image.
The measurement operator is given by

Φf = f (ω) ω∈Ω ,
where f denotes the (discrete or continuous) Fourier transform of f and Ω is a set of frequencies. The main di erence between these two imaging techniques is the type of Fourier samples they give access to. For X-ray tomography Ω is a set of points located on radial lines in the Fourier plane. In the case of MRI, these points are located on smooth curves.

Laplace transform sampling ( uorescence microscopy). A major challenge in modern biology is to understand mechanisms at the sub-cellular level. This requires the development of so-called super-resolution microscopy techniques, able to overcome the di raction limit (see Figure 2 for an example of images obtained using such techniques). A popular method to achieve this is single-molecule localization microscopy (SMLM) [START_REF] Betzig | Imaging Intracellular Fluorescent Proteins at Nanometer Resolution[END_REF], Rust et al., 2006, Shro et al., 2008]. It consists in introducing uorophores in a biological sample, and recording the light they emit after illuminating a random subset of them. In this context, acquiring depth information on the sample is a major challenge that can be adressed using several techniques [START_REF] Huang | Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy[END_REF], Pavani et al., 2009]. Among them, the multi-angle total internal re ection uorescence method (MA-TIRF) [START_REF] Boulanger | Fast high-resolution 3D total internal re ection uorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF], Santos et al., 2014, Soubies et al., 2014] consists in illumating the biological sample along di erent directions. The resulting measurements correspond to a sampling of the Laplace transform (in depth) of the sought-after image, which is modeled as a discrete measure. In other words, neglecting the horizontal blur, this corresponds to the choice of E = M(R + ) and

Φ : M(R + ) → R K µ → ˆR+ exp(-s k x) dµ(x) k=1,...,K
, where s k ∈ R + for every k ∈ {1, ..., K}. Let us stress that this model is a gross simpli cation of the real models used in uorescence microscopy, and we refer the reader to [Denoyelle, 2018, Chapter 5] for more information on this subject.

V 1.2.1. Presentation

A fruitful idea for leveraging a priori knowledge on the sought-after signal f 0 is that of variational regularization. It consists, in the case of noiseless observations, in nding an approximation of f 0 by solving inf

f ∈E R(f ) s.t. Φf = y 0 , (P 0 (y 0 ))
where R : E → R ∪ {+∞} is called the regularizer. This mapping is chosen in order to promote solutions having the same structure as f 0 . Loosely speaking, the value of R(f ) should be small if f has the desired structure, and large otherwise. The choice of R given a priori assumptions on the sought-after signal is discussed in Section 1.2.2. If one has access to noisy observations, i.e. to y = Φf 0 + w, then minimizing R among functions f satisfying Φf = y makes little sense (even f 0 would not belong to the admissible set). One hence solves instead inf

f ∈E R(f ) + 1 2λ Φf -y 2 H , (P λ (y))
with λ a positive real number called the regularization parameter. Typically, λ should be chosen according to an estimate of w H . Having de ned these two variational problems, it is natural to investigate if solving them allows to solve the considered inverse problems. Three questions therefore arise.

1. Is f 0 the unique solution of (P 0 (y 0 ))?

2. If w is small and λ well chosen, are solution of (P λ (y 0 + w)) close to solutions of (P 0 (y 0 ))?

3. Can one numerically solve (P 0 (y 0 )) and (P λ (y))?

1.2.2. Choice of the regularizer and a priori information Discrete 1 regularization. The most important (and rst historical) example of regularizer is arguably the 1 -norm in the Euclidean space. From the 70's, 1 regularization has been used empirically in various elds (notably in seismic imaging) to recover sparse vectors. It was later studied mathematically by David Donoho and its co-authors (among others) from the 90's (see e.g. [Donoho, 1992, Chen et al., 1998] and also [Tibshirani, 1996] for applications in statistics).

Convex regularizers. In order to obtain convex problems (P 0 (y 0 )) and (P λ (y)), it is natural to choose R among convex functions. Although non-convex regularizers could also be considered, being able to leverage the theory of convex optimization is really useful for analyzing variation regulatization techniques.

Representer theorems and convex regularization. The following series of recent works [START_REF] Chandrasekaran | The Convex Geometry of Linear Inverse Problems[END_REF], Boyer et al., 2019, Bredies and Carioni, 2019] has recently shed a new light on convex regularization. In these articles, the link between the choice of a regularizer and the structural properties of some solutions of (P 0 (y 0 )) and (P λ (y)) is investigated. Their main nding is that the atoms promoted by a convex regularizer are the extreme points of its associated unit ball: loosely speaking, if H has nite dimension, then some solutions of (P 0 (y 0 )) and (P λ (y)) are linear combinations of extreme points of {R ≤ 1}. We provide below a few examples of regularizers and describe the atoms they promote.

• If E = R n and R = . 1 is the 1 -norm, then the extreme points of {R ≤ 1} are the vectors ±e i with i ∈ {1, ..., d}, where (e i ) i=1,...,n is the canonical basis of R n .

• If E = R n×m and R = . * is the nuclear norm of matrices then the extreme points of {R ≤ 1} are the matrices

u T v with u ∈ R n , v ∈ R m and u 2 = v 2 = 1. • If E = M(X) with X and R = |.|(X)
is the total variation (of measures) then the extreme points of {R ≤ 1} are the measures ±δ x with x ∈ X.

•

If E = L 2 (R 2 ) and R = |D.|(R 2
) is the total (gradient) variation then the extreme points of {R ≤ 1} are the functions ±1 E /P (E) with E a simple set with positive nite measure1 .

Besides extreme points, the results of [START_REF] Boyer | On Representer Theorems and Convex Regularization[END_REF] show that the structure of the solution set of (P 0 (y 0 )) or (P λ (y)) is linked to the structure of the faces of {R ≤ 1}. The importance of closely studying these faces is highlighted in Part 2. If one has access to m observations (i.e. H = R m ), then some solutions lie in a face of dimension at most m -1 (or m in degenerate cases). In particular, these solutions can be written as the sum of at most m (or m + 1) extreme points of {R ≤ 1}.

Identi ability

If the rst question stated in Section 1.2.1 has a positive answer, i.e. if f 0 is the unique solution of (P 0 (y 0 )), then it is said to be identi able. Even if we rather focus on noise robustness results in this thesis, let us brie y discuss this question.

Discrete 1 regularization. We begin by considering the recovery of sparse vectors using the 1 norm as a regularizer (i.e. E = R n , R = • 1 and H = R m ). The rst historical example of identi ability condition is the nullspace property. This name was rst used in [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF], although the condition also appeared in [START_REF] Gribonval | Sparse representations in unions of bases[END_REF] and implicitly in previous works. If Φ is the measurement matrix, then every s-sparse vector is identi able if and only if the nullspace property of order s is satis ed. This condition is di cult to check in practice. To circumvent this issue, the stronger restricted isometry property was introduced in [START_REF] Candes | Decoding by linear programming[END_REF]. Even if computing the restricted isoperimetry constants of a matrix is known to be NP-hard, large classes of random matrices satisfy the property with high probability (see for instance [START_REF] Candes | Near-Optimal Signal Recovery From Random Projections: Universal Encoding Strategies[END_REF]). We refer the reader to the reference book [START_REF] Foucart | A Mathematical Introduction to Compressive Sensing[END_REF] for more details.

Extensions. Su cient identi ability conditions also exist for low-rank matrix recovery, i.e. when E = R n×m and R = • * (see for instance [Foucart and Rauhut, 2013, Section 4.6]). Likewise, perfect recovery results were proven in the sparse spikes setting, i.e. with E = M(X) and R = | • |(X). We provide further details on this case in Section 1.3.2.

Noise robustness

As mentioned above, the second question stated in Section 1.2.1 is closely related to noise robustness. Indeed, if a signal f 0 is identi able and one proves that solutions of (P λ (y)) converge to a solution of (P 0 (y 0 )) when w is small and λ well chosen, then solving (P λ (y)) allows to approximately recover f 0 in a low noise regime.

Stability and model stability. Noise robustness results typically fall into two categories. The rst ones are concerned with the convergence of a solution f λ,w of (P λ (y 0 + w)) towards f 0 with respect to some distance on the signal space E. This distance is typically the

2 distance if E = R n , the L 2 distance if E = L 2 (R 2 ), or an optimal transport-based distance if E = M(X).
This notion of stability is however quite weak. Indeed, if f 0 is a combination of a few atoms, a convergence result of this kind does not tell anything about the structure of f λ,w . To be more precise, it does not allow to answer the following question: is f λ,w made of the same number of atoms as f 0 , and are these atoms related to those appearing in the decomposition of f 0 ? This shows the necessity of a second category of noise robustness results, allowing to answer these questions. They are sometimes called model stability results, and are the ones we focus on in this thesis.

Discrete 1 regularization. Given a sparse vector, several works investigated whether solving (P λ (y 0 + w)) (with E = R n and R = . 1 ) in a low noise regime allows to exactly recover its support. To name only a few, let us cite [Fuchs, 2004, Dossal and Mallat, 2005, Tropp, 2006], in which su cient conditions are derived. Loosely speaking and using a modern terminology, exact recovery is possible if some minimal norm dual certi cate interpolates the sign of the signal coordinates.

Extensions. In [START_REF] Vaiter | Model selection with low complexity priors[END_REF], Vaiter et al., 2018], the analysis of 1 regularization was extended to the family of so-called partly smooth regularizers, which promote signals having some notion of low complexity. Numerous popular regularizers fall in this category. To name only a few, let us mention the 1 norm, the 1 -2 norm (used to enforce group sparsity), and the nuclear norm of matrices. An extension to so-called mirror-strati able regularizers was also studied in [START_REF] Fadili | Sensitivity Analysis for Mirror-Strati able Convex Functions[END_REF], Fadili et al., 2019]. This last series of works has similarities with the dual-based approach of [START_REF] Duval | Exact Support Recovery for Sparse Spikes Deconvolution[END_REF], in which an exact support recovery result is proved for discrete measures. We discuss this topic in greater details in Section 1.3.2

G

At rst glance, solving (P λ (y)) or (P 0 (y 0 )) with E an in nite dimensional space seems intractable. Considering such problems hence only seems to be of theoretical interest. However, a recent line of works demonstrated that it is in fact possible to solve (P λ (y)) and (P 0 (y 0 )) e ciently when E = M(X) and R = |.|(X) is the total variation of measures. Doing so even has numerous advantages, both from a numerical point of view and with respect to recovery guarantees. This section provides a review of these advances, and situates the subject of this thesis within this context.

Sparse spikes super-resolution on thin grids

Sparse spikes recovery and the BLASSO. As explained in Sections 1.1.1 and 1.1.2, the recovery of a discrete measure µ 0 = s i=1 a 0,i δ x 0,i from noisy linear measurements y = Φµ 0 + w has numerous applications. We also mentioned that solving (P λ (y 0 + w)) where y 0 = Φµ 0 with E = M(X) and R = |.|(X) is a good way to approximate µ 0 . The resulting problem, which writes inf

µ∈M(X) |µ|(X) + 1 2λ Φµ -y 2 H (P λ (y))
and is known as the BLASSO1 , is however in nite-dimensional. How to numerically solve it is hence unclear.

Restriction to a nite grid. A rst approach one can consider is to look for an approximate solution of (P λ (y)) whose support lies in some discrete grid G = {x i } n-1 i=0 ⊂ X, and therefore to solve (P λ (y)) with the additional constraint Supp(µ) ⊂ G. This yields the nite dimensional problem

inf a∈R n a 1 + 1 2λ Φ G a -y 2 H , (P G λ (y))
which is the celebrated basis pursuit, also known as the LASSO (see for instance [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], Tibshirani, 1996]), where Φ G is de ned by

Φ G : R n → H a → Φ n i=1 a i δ x i .
Thin grids. Since no information on the support of µ 0 is known a priori, it is natural to consider increasingly "thinner" grids G, in order to cover the whole domain X. However, the analysis conducted in [START_REF] Duval | Sparse regularization on thin grids I: The Lasso[END_REF] shows this approach yields reconstructed measures whose structure is di erent from that of µ 0 . To be more precise, each spike in µ 0 is estimated by a pair of neighbooring spikes. This demonstrates the impossibility to exactly recover the support of µ 0 by relying on 1 regularization techniques on (even arbitrarily thin) nite grids.

Grid-free sparse spikes recovery

Considering the negative result presented above, it is natural to wonder if, despite its in nitedimensional nature, numerically solving the BLASSO problem (P λ (y)) is possible. Over the last ten years, numerous methods have in fact been developed to achieve this. We brie y review them here, before discussing available results on noise robustness. We refer the reader to the lecture notes [Poon, 2019] for a more detailed introduction to sparse spikes recovery using the BLASSO.

Numerical resolution of the BLASSO. The rst method allowing to solve the BLASSO was introduced in [Tang et al., 2013, Candès andFernandez-Granda, 2014]. It relies on the equivalence, in the speci c case of 1D Fourier measurements, of (P λ (y)) with a nite dimensional semi-de nite program, for which e cient solvers exist. Extensions to higher dimensions rely on the so-called Lasserre hierarchy (see e.g. [START_REF] Catala | A Low-Rank Approach to O -the-Grid Sparse Superresolution[END_REF], [Catala, 2020, Chapter 2] and also [De Castro et al., 2017] for the case of polynomial measurements). For general measurements, approaches based on the conditional gradient or Frank-Wolfe algorithm, which does not require any Hilbertian structure (and is hence particularly well suited for optimization on the space of measures), were developed in [START_REF] Bredies | Inverse problems in spaces of measures[END_REF], Boyd et al., 2017, Denoyelle et al., 2019]. Finally, let us mention the work of [Chizat, 2022], which is based on conic particle gradient descent.

Identi ability. The groundbreaking work [START_REF] Candès | Towards a Mathematical Theory of Super-resolution[END_REF] shows that a discrete measure can be exactly recovered from 2f c + 1 Fourier samples, provided the minimal distance ∆ between its atoms satis es ∆ ≥ 2/f c . In [START_REF] Tang | Compressed Sensing O the Grid[END_REF], it has also been shown that this result remains valid with high probability if a small number of Fourier samples are randomly selected. Let us stress that, for positive measures, the minimal separation condition can be lifted (see for instance [De Castro andGamboa, 2012, Schiebinger et al., 2018]).

Noise robustness. The convergence of solutions µ λ,w of (P λ (y 0 + w)) in a low noise regime has been investigated in several works. The weak-* convergence of µ λ,w towards a solution µ 0 of the noiseless problem is proved in [START_REF] Bredies | Inverse problems in spaces of measures[END_REF]. In a similar spirit, [START_REF] Candès | Super-Resolution from Noisy Data[END_REF], Fernandez-Granda, 2013, Azaïs et al., 2015] provide error bounds (based on [START_REF] Burger | Convergence rates of convex variational regularization[END_REF]) which show that µ λ,w concentrates around the support of µ 0 . However, these results do not address the question of support recovery: if µ 0 is discrete and identi able, is µ λ,w discrete, and does it have the same number of spikes as µ 0 ?

Exact support recovery. In [START_REF] Duval | Exact Support Recovery for Sparse Spikes Deconvolution[END_REF]], an answer to this last question is provided. The authors prove that, if a so-called non-degenerate source condition holds, µ λ,w has exactly the same number of spikes as µ 0 in a low noise regime, and the locations of these spikes (as well as the associated amplitudes) smoothly converge to the limit ones. This is in sharp contrast with the negative result presented for nite grids in Section 1.3.1.

Perspectives

In Section 1.3.2 above, we have seen that, in the context of sparse spikes recovery, grid-free approaches both have theoretical and numerical bene ts. In this subsection, we mention recent attempts to develop similar methods for recovering other kinds of signals. We nally introduce the problem this thesis is dedicated to, and advocate for studying noise robustness properties and developing numerical methods using the fruitful ideas presented in Section 1.3.2.

Dynamic sparse spikes recovery. As mentioned in Section 1.1.1, the recovery of timedependant discrete measures has recently been investigated. Several regularizers were introduced to tackle the xed and varying mass settings, and extreme points results were derived. Although we are not aware of any identi ability or noise robustness result, numerical methods received a lot of attention, and approaches based on the Frank-Wolfe (or conditional gradient) algorithm were very recently proposed in [Bredies et al., 2022b, Duval andTovey, 2021].

Curve recovery. We also discussed in Section 1.1.1, the reconstruction of images which are the superposition of a few 1D curves. To our knowledge, the only work in which this task is investigated is the recent work [START_REF] Laville | O -the-grid curve reconstruction through divergence regularisation: An extreme point result[END_REF]. Its authors prove an extreme point result which we brie y state below. Let E = M(X) 2 be the dual of C 0 (X, R 2 ) (equipped with the supremum norm), and let R be de ned by

R : M(X) 2 → R ∪ {+∞} µ → |µ|(X) + |div(µ)|(X) ,
with div(µ) the divergence of µ (de ned in the sense of distributions). Then the extreme points of {R ≤ 1} are the measures µ γ /R(µ γ ) where µ γ is the measure canonically associated to a simple oriented Lipschitz curve γ (see [Laville et al., 2022, De nition 7] for a precise de nition). An interesting avenue for future research is to investigate the performance of this regularizer for solving inverse problems (i.e. to derive identi ability and noise robustness results), and also to develop numerical methods allowing to solve (P 0 (y 0 )) and (P λ (y)) e ciently. Let us also stress that the measures this work deals with are associated to oriented curves, and that they also carry tangential information. Recovering measures supported on unoriented curves has, to our knowledge, not be considered, although they seem a natural model for images such as those of Figure 1.

Reconstruction of piecewise constant images. In this thesis, we are concerned with the recovery of piecewise constant images, i.e. which are the superposition of a few simple shapes. To recover such objects, the natural choice of regularizer is the total (gradient) variation. The rst appearance of this functional for imaging applications dates back to the pioneering work of Rudin Osher and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], which focuses on the denoising problem. Since then, the performance of this regularizer has been the subject of a large number of works. It has notably been shown in [START_REF] Chambolle | Geometric properties of solutions to the total variation denoising problem[END_REF], Iglesias et al., 2018] that, in a small noise regime, solutions of (P λ (y 0 + w)) converge to a solution of (P 0 (y 0 )) (which is unique and equal to the image to recover in the case of denoising). The type of convergence proved in this result is already quite strong. Namely, it is shown that the level sets of solutions of (P λ (y 0 + w)) converge to the level sets of solutions of (P 0 (y 0 )) (see Proposition 1.12 for a precise statement). This convergence result for level sets is of particular interest in the context of imaging. Even if the analysis presented in Section 1.2.2 suggests that the total variation is particularly well suited to recover piecewise constant functions, the above mentioned works do not focus on unknown images having this structure. If we make this assumption, we argue that the emphasis should be on structurepreserving convergence guarantees, in the spirit of [START_REF] Duval | Exact Support Recovery for Sparse Spikes Deconvolution[END_REF]. To be more speci c, one should study whether the reconstructed image is composed of the same number of shapes than the sought-after image, and whether these shapes are close to the unknown ones. Finally, let us mention that, in this setting, the problem of deriving su cient identi ability conditions is to our knowledge mostly open. Still, let us point out that, in [START_REF] Bredies | A perfect reconstruction property for PDE-constrained total-variation minimization with application in Quantitative Susceptibility Mapping[END_REF], a perfect reconstruction result is obtained in the case where the measurement associated to a function is its image under a linear partial di erential operator with unknown boundary conditions.

The main aim of this thesis is to study total variation regularization from a sparse recovery viewpoint.

T

The use of the total variation as a regularizer started with the pioneering work of Rudin Osher and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] in the context of image denoising. This functional penalizes oscillations while allowing discontinuities. Reconstruction methods based on its minimization therefore have the striking feature of preserving image edges (see Figure 3 for a few examples).

Although state of the art algorithms now exhibit much better performance, the total variation is still routinely used in speci c applications, and remains an import baseline. In this section, we provide an introduction to total variation regularization. We review existing theoretical results and numerical methods, and also mention the open problems that we attempt to address in this thesis.

Figure 3 -Reconstructions obtained with [Getreuer, 2012a[START_REF] Getreuer ; Getreuer | Total Variation Inpainting using Split Bregman[END_REF] in the case of deconvolution ( rst and third rows) and inpainting (second and fourth rows). Left: observations y = Φu 0 + w, middle: unknown image u 0 , right: approximate solution of (P λ (y 0 + w)) (source images: bazzier, Le Pays Malouin).

2.1. P

The total variation

We begin this section by de ning the total variation functional. Then, we introduce a few notions to be able to state the extreme point result mentioned in Section 1.2.2, which characterizes the atoms promoted by this regularizer. Finally, we give several results about the subdi erential of the total variation, which are useful for the analysis of the variational problems (P 0 (y 0 )) and (P λ (y)).

De nition. If u ∈ L 1 loc (R 2 )
, its total variation is de ned by

TV(u) def. = sup - ˆR2 u divφ φ ∈ C ∞ c (R 2 , R 2 ), φ ∞ ≤ 1 . (1)
div φ = ´∂E φ • ν E dH 1 , which yields TV(1 E ) = H 1 (∂E).
Further details about the total variation can be found in [Ambrosio et al., 2000, Chapter 3] and [Evans and Gariepy, 2015, Chapter 5]. Let us also mention [Chambolle et al., 2010a] for more information on its use in imaging.

Remark 1.1

In this thesis, we extensively work with square integrable functions with nite total variation. These do not belong to the more classical space BV(R 2 ) of integrable functions with nite total variation. The main interest of this choice is that, in R 2 , the total variation dominates the L 2 norm (and more generally the L N/(N -1) norm in R N ), as shown by the isoperimetric inequality (3) presented below. As a consequence, the total variation unit ball is weakly compact in L 2 (R 2 ), which is hence the natural space in which we can obtain existence of solutions for variational problems with a total variation term.

In all the following we consider TV as a mapping from L 2 (R 2 ) to R ∪ {+∞}. This function is convex, proper and lower semi-continuous.

Sets of nite perimeter. If a measurable set E ⊂ R 2 is such that P (E)

def.

= TV(1 E ) is nite, it is said to have nite perimeter. Such sets have similar properties to those of smooth sets, in a measure-theoretic sense. In particular, they satisfy a generalized Gauss-Green formula, and have an approximate normal and an approximate tangent space. They are central in the study of geometric variational problems, mainly thanks to useful compactness results. We refer the reader to [Maggi, 2012, Part II] for further details on this subject. We collect notations and a few properties of sets of nite perimeter in Appendix A.

The isoperimetric inequality

We now state the useful isoperimetric inequality, and take this opportunity to make a digression on its quantitative version. This allows us to introduce several topics of interest for our purpose, namely shape calculus and stability in geometric variational problems.

For every set of nite perimeter E ⊂ R 2 , the isoperimetric inequality states that

min(|E|, |E c |) ≤ c 2 P (E) , (2)
with equality if and only if E is a ball, and where c 2 def.

= 1/ √ 4π is the isoperimetric constant (see e.g. [Maggi, 2012, Chapter 14]). In particular, if E is a set of nite perimeter, either E or E c has nite measure. The functional version1 of (2) is:

∀u ∈ L 2 (R 2 ), u L 2 (R 2 ) ≤ c 2 TV(u) . (3)
A sharp quantitative version of (2) was rst proved in [START_REF] Fusco | The sharp quantitative isoperimetric inequality[END_REF] using symmetrization techniques (see also the review article [Fusco, 2015]). Denoting B = B(0, 1) the unit ball, for every measurable set E ⊂ R 2 with |E| = |B|, it writes:

min x∈R 2 |(x + E) B| 2 ≤ C[P (E) -P (B)] . (4)
Another proof, relying on optimal transport, was given in [START_REF] Figalli | A mass transportation approach to quantitative isoperimetric inequalities[END_REF]. A third approach, which is the most relevant for our purpose, was proposed in [START_REF] Cicalese | A Selection Principle for the Sharp Quantitative Isoperimetric Inequality[END_REF]. It relies on two main ingredients. The rst is a socalled selection principle, which allows to reduce the proof to nearly spherical sets, i.e. sets whose boundary is the normal graph of a su ciently small Lipschitz function on the sphere. The second is the proof of (4) for such sets, which dates back to [Fuglede, 1989].

The proof of Fuglede relies on shape calculus. Since nearly spherical sets are parametrized by a function ϕ on the sphere, one can look at the mapping ϕ → P (B ϕ ), where B ϕ is the nearly spherical set associated to ϕ, and use calculus to study the stability of B = B 0 as a minimizer. To be more precise, one can show that the second order derivative of this mapping at zero, which is called the second order shape derivative of the perimeter, is coercive, which ultimately yields the result.

The derivation of stability results using second order shape derivatives has been generalized to other geometric functionals in several works (see [START_REF] Dambrine | Stability in shape optimization with second variation[END_REF] and the references therein). The noise robustness analysis we conduct in Part 2 heavily relies on these ideas, which we apply to the prescribed curvature problem presented below.

Indecomposable and simple sets. Now, we introduce the notion of indecomposable and simple sets, which are the measure-theoretic analogues of connected and simply connected sets (we refer the reader to [START_REF] Ambrosio | Connected components of sets of nite perimeter and applications to image processing[END_REF] for more details). A set of nite perimeter E ⊂ R 2 is said to be decomposable if there exists a partition of E in two sets of positive measure A and B with P (E) = P (A) + P (B). We say that E is indecomposable if it is not decomposable

1 . An indecomposable set E is called simple if E = R 2 or |E| < +∞ and R 2 \ E is indecomposable.
Extreme points of the total variation unit ball. We are now able to state the extreme point result mentioned Section 1.2.2. Proposition 1.2 ( [Fleming, 1957, Ambrosio et al., 2001])

The extreme points of the convex set

{TV ≤ 1} def. = u ∈ L 2 (R 2 ) TV(u) ≤ 1
are the functions of the form ±1 E /P (E), where E is a simple set with 0 < |E| < +∞.

This result, in conjunction with the representer theorems mentioned in Section 1.2.2, justies the a rmation that using the total variation as a regularizer promotes piecewise constant functions (i.e. that are the weighted sum of a few indicator functions of simple sets).

Subdi erential. Let us now collect several results on the subdi erential of TV, which are useful to derive and analyze the dual problems of (P 0 (y 0 )) and (P λ (y)). Since TV is the support function of the convex set

C def. = div z z ∈ C ∞ c (R 2 , R 2 ), z ∞ ≤ 1 , its subdi erential at 0 is the closure of C in L 2 (R 2 ), that is ∂TV(0) = C = div z z ∈ L ∞ (R 2 , R 2 ), divz ∈ L 2 (R 2 ), ||z|| ∞ ≤ 1 .
(5)

We also have the following useful identity:

∂TV(0) = η ∈ L 2 (R 2 ) ∀u ∈ L 2 (R 2 ), ˆR2 η u ≤ TV(u) .
Finally, the sudi erential of TV at some u ∈ L 2 (R 2 ) is given by:

∂TV(u) = η ∈ ∂TV(0) ˆR2 η u = TV(u) .
Hence, if η ∈ ∂TV(u), then η is an element of ∂TV(0) = C for which the supremum in the de nition of the total variation is attained.

From the subdi erential to the level sets. Now, there is an important su cient and necessary condition for η ∈ ∂TV(u) to hold. In the rest of the article, given u : R 2 → R and t ∈ R, we use the following notation:

U (t) def. = {x ∈ R 2 | u(x) ≥ t} if t ≥ 0 , {x ∈ R 2 | u(x) ≤ t} otherwise.
It is worth noting that, if u ∈ L 2 (R 2 ), then |U (t) | < +∞ for all t = 0. With this notation, we may state the following result.

Proposition 1.3 (see e.g. [START_REF] Chambolle | Geometric properties of solutions to the total variation denoising problem[END_REF]) Let u ∈ L 2 (R 2 ) be such that TV(u) < +∞, and let η ∈ L 2 (R 2 ). Then the following conditions are equivalent.

(i) η ∈ ∂TV(u) .

(ii) η ∈ ∂TV(0) and the level sets of u satisfy ∀t > 0, P (U (t) ) = ˆU(t) η and ∀t < 0, P (U (t) ) = -ˆU(t) η .

(iii) The level sets of u satisfy

∀t > 0, U (t) ∈ Argmin E⊂R 2 , |E|<+∞ P (E) -ˆE η , ∀t < 0, U (t) ∈ Argmin E⊂R 2 , |E|<+∞ P (E) + ˆE η .
The geometric variational problem appearing in (iii), which is

inf E⊂R 2 , |E|<+∞ P (E) -ˆE η , (6)
is called the prescribed curvature problem associated to η. This terminology stems from the fact that, if η is su ciently regular, every solution of (6) has a (scalar) distributional curvature (de ned in Appendix A) equal to η. This problem plays a crucial role in the analysis of total variation regularization, as explained below. If η ∈ ∂TV(0), one should note that a set of nite measure E solves (6) if and only if P (E) = ´E η, which is why (ii) and (iii) are equivalent.

Primal and dual problems, dual certi cates

We can now specialize the variational problems introduced in Section 1.2.1 to our setting, which corresponds to the choice of E = L 2 (R 2 ) and R = TV.

Assumptions. In all the following, we assume that Φ :

L 2 (R 2 ) → H is a continuous linear operator, which is equivalent to the existence of ϕ ∈ L 2 (R 2 , H) such that Φ : L 2 (R 2 ) → H u → ˆR2 u(x) ϕ(x) dx .
We also assume that y 0 = Φu 0 with TV(u 0 ) < +∞.

Primal problems. As mentioned above, in our setting, the two variational problems introduced in Section 1.2.1 are inf

u∈L 2 (R 2 )
TV(u) s.t. Φu = y 0 , (P 0 (y 0 ))

inf u∈L 2 (R 2 ) TV(u) + 1 2λ Φu -y 2 H . (P λ (y))
Existence of solutions for (P 0 (y 0 )) and (P λ (y)) can be shown using the direct method of the calculus of variations and the isoperimetric inequality (3).

Dual problems. In the remaining of this section, we collect useful results following from standard duality arguments. We refer the reader to [Iglesias et al., 2018, Section 2] for their proof (which had been established in the particular case of denoising in [START_REF] Chambolle | Geometric properties of solutions to the total variation denoising problem[END_REF]).

The Fenchel-Rockafellar dual problems of (P 0 (y 0 )) and (P λ (y)) are:

sup p∈H p, y 0 H s.t. Φ * p ∈ ∂TV(0) , (D 0 (y 0 )) sup p∈H p, y H - λ 2 ||p|| 2 H s.t. Φ * p ∈ ∂TV(0) . (D λ (y))
The existence of a solution of (D 0 (y 0 )) is not always guaranteed. Requiring that there indeed exists a maximizer is called the source condition in the literature. On the contrary, (D λ (y)) can be reformulated as the projection problem of y/λ on the closed convex set {p ∈ H | Φ * p ∈ ∂TV(0)}.

It hence has a unique solution. Strong duality holds between (P 0 (y 0 )) and (D 0 (y 0 )) (respectively (P λ (y)) and (D λ (y))), as stated in the following two propositions.

Proposition 1.4 Strong duality holds between (P 0 (y 0 )) and (D 0 (y 0 )). Moreover, if there exists a solution p to (D 0 (y 0 )), then for every solution u of (P 0 (y 0 )) we have

Φ * p ∈ ∂TV(u) . (7)
Conversely, if (u, p) ∈ L 2 (R 2 ) × H with Φu = y 0 and (7) holds, then u and p respectively solve (P 0 (y 0 )) and (D 0 (y 0 )).

Proposition 1.5 Strong duality holds between (P λ (y)) and (D λ (y)). Moreover, denoting p the unique solution of (D λ (y)), for every solution u of (P λ (y)) we have

Φu = y -λ p , Φ * p ∈ ∂TV(u) . (8)
Conversely, if (8) holds, then u and p respectively solve (P λ (y)) and (D λ (y)).

Remark 1.6 Although there might not be a unique solution to (P λ (y)), Proposition 1.5 shows that all of them have the same image by Φ and the same total variation.

Dual certi cates. If η = Φ * p and η ∈ ∂TV(u), we call η a dual certi cate for u in (P 0 (y 0 )), as its existence proves the optimality of u for (P 0 (y 0 )) (assuming u is admissible). Similarly, if η = -Φ * (Φuy)/λ and η ∈ ∂TV(u), we call η a dual certi cate for u in (P λ (y)). There could be multiple dual certi cates associated to (P 0 (y 0 )). One of them, the minimal norm certi cate, plays a crucial role in our analysis. A quick look at the objective of (D λ (y)) indeed suggests that, as λ goes to 0, its solution converges to the solution of the limit problem (D 0 (y 0 )) with minimal norm. This is Proposition 1.8 below.

De nition 1.7

If there exists a solution to (D 0 (y 0 )), the minimal norm dual certi cate associated to (P 0 (y 0 )), denoted η 0 , is de ned as

η 0 = Φ * p 0 with p 0 = argmin p H s.t. p solves (D 0 (y 0 )) .
If λ > 0, we denote p λ,w the unique solution of (D λ (y 0 + w)), and η λ,w = Φ * p λ,w the associated dual certi cate. Noise robustness results extensively rely on the behaviour of η λ,w as λ and w go to zero. This behaviour is described by the following results.

Proposition 1.8

If there exists a solution to (D 0 (y 0 )), then p λ,0 converges strongly to p 0 as λ → 0.

Since p λ,w is the projection of (y 0 +w)/λ onto the closed convex set {p ∈ H | Φ * p ∈ ∂TV(0)}, the non-expansiveness of the projection mapping yields

∀(λ, w) ∈ R * + × H, p λ,w -p λ,0 H ≤ w H λ , (9)
and hence

∀(λ, w) ∈ R * + × H, η λ,w -η λ,0 L 2 (R 2 ) ≤ Φ * w H λ .
As a result, if λ → 0 and w H /λ → 0, the dual certi cate η λ,w converges in L 2 (R 2 ) to the minimal norm certi cate η 0 .

N

In this section, we review existing results concerning noise robustness (i.e. question 2 in Section 1.2.1). The main references on this topic are [START_REF] Chambolle | Geometric properties of solutions to the total variation denoising problem[END_REF], Iglesias et al., 2018].

Setting. In all the following we consider two sequences (w n ) n≥0 ∈ H N and (λ n ) n≥0 of noises and regularization parameters, and we assume that (u n ) n≥0 is such that u n solves (P λn (y 0 + w n )) for every n ≥ 0.

Set convergence. Let us de ne a useful notion of set convergence, known as Painlevé-Kuratowski set convergence. We refer the reader to [Rockafellar and Wets, 1998, Chapter 4] for more details.

De nition 1.9

If (E n ) n≥0 is a sequence of sets in R d , its inner and outer limits are de ned by

lim inf n→+∞ E n def. = x ∈ R d lim sup n→+∞ dist(x, E n ) = 0 , lim sup n→+∞ E n def. = x ∈ R d lim inf n→+∞ dist(x, E n ) = 0 .
If there exists a set E such that

lim inf n→+∞ E n = lim sup n→+∞ E n = E ,
we say that (E n ) n≥0 converges to E and write lim n→+∞

E n = E.
If a sequence of closed sets (E n ) n≥0 is bounded (i.e. there exists R > 0 such that E n ⊂ B(0, R) for every n ≥ 0), then the above notion of set convergence is equivalent to convergence in the Hausdor sense, that is to

lim n→+∞ dist(•, E n ) -dist(•, E) ∞ = 0 .
First convergence result. General convergence results given in [START_REF] Hofmann | A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators[END_REF] apply to our setting. This yields the following proposition.

Proposition 1.10 ( [Hofmann et al., 2007, Theorem 3.5]) If λ n → 0 and w n 2 H = o(λ n ) then, up to the extraction of a subsequence (not relabeled), we have that (u n ) n≥0 converges strongly in L1 loc (R 2 ) and weakly in L 2 (R 2 ) to a solution u * of (P 0 (y 0 )). Moreover, we have TV(u n ) → TV(u * ).

We stress that the strong L 1 loc convergence of (u n ) n≥0 towards u * and the fact TV(u n ) → TV(u * ) imply Du n * Du * , which in turn implies 1

Supp(Du * ) ⊆ lim inf n→+∞ Supp(Du n ) . ( 10
)
Properties of the level sets in the low noise regime. Stronger convergence results can be obtained by having a closer look at the level sets of u n . Indeed, from Propositions 1.3 and 1.5, we know they are solutions of the prescribed curvature problem associated to η λn,wn . This information can be exploited to obtain uniform properties of the level sets in the low noise regime.

The main ingredient to achieve this is the following lemma, that we use in di erent parts of this manuscript.

Lemma 1.11 ([Chambolle et al., 2016, Section 5]) Let (η n ) n≥0 ⊂ ∂TV(0) be a sequence converging strongly in L 2 (R 2 ) to η ∞ , and let E be de ned by

E def. = E ⊂ R 2 , 0 < |E| < +∞ ∃n ∈ N ∪ {∞}, P (E) = ˆE η n .
Then the following holds:

1. inf E∈E P (E) > 0 and sup E∈E P (E) < +∞, 2. inf E∈E |E| > 0 and sup E∈E |E| < +∞, 3. there exists R > 0 such that, for every E ∈ E, it holds E ⊂ B(0, R),
4. there exists r 0 > 0 and C ∈ (0, 1/2) such that for every r ∈ (0, r 0 ] and E ∈ E:

∀x ∈ ∂E, C ≤ |E ∩ B(x, r)| |B(x, r)| ≤ 1 -C .
Point 4 in Lemma 1.11 is a weak regularity property. In particular, it ensures that E does not have cusps.

Convergence of level sets. Before stating Proposition 1.12 below, Proposition 1.12 ( [START_REF] Chambolle | Geometric properties of solutions to the total variation denoising problem[END_REF], Iglesias et al., 2018]) If (D 0 (y 0 )) has a solution, λ n → 0 and

w n H λ n ≤ 1 4 c 2 Φ * ,
then (Supp(u n )) n≥0 is bounded and, up to the extraction of a subsequence (not relabeled), we have that (u n ) n≥0 converges strictly in BV(R 2 ) to a solution u * of (P 0 (y 0 )). Moreover, for almost every t ∈ R, we have:

U (t) n U (t) * -→ 0 and ∂U (t) n -→ ∂U (t) * .
Remark 1.13

The result of Proposition 1.12 slightly di ers from the one proved in [START_REF] Chambolle | Geometric properties of solutions to the total variation denoising problem[END_REF], Iglesias et al., 2018]. Indeed, the latter guarantees the convergence of level sets up to the extraction of a further subsequence. We argue that this can be avoided by using Proposition D.1, that we prove in Appendix D.

The extended support. The notion of extended support associated to u 01 was introduced in [Chambolle et al., 2016, Section 6]. If (D 0 (y 0 )) has a solution and η 0 is the minimal norm certi cate de ned in Appendix B, this set, denoted Ext(Du 0 ), is de ned as follows:

Ext(Du 0 ) def. = Supp(Du) η 0 ∈ ∂TV(u) .
From Proposition 1.3, the extended support is also characterized by

Ext(Du 0 ) = ∂ * E |E| < +∞ and P (E) = ˆE η 0 ,
where ∂ * E denotes the reduced boundary of E (see Appendix A for a de nition). In other words, Ext(Du 0 ) contains the reduced boundary of all solutions of the prescribed curvature problems associated to the minimum norm dual certi cate and its opposite. In the above-mentioned work, the following result was proved in the case of denoising (i.e. Φ = Id). One may check that it still holds in the case of general measurement operators.

Proposition 1.14

If (D 0 (y 0 )) has a solution, λ n → 0 and w n = o(λ n ) then lim sup n→+∞ Supp(Du n ) ⊆ Ext(Du 0 ) ,
or equivalently, for every r > 0 there exists n 0 ∈ N such that

∀n ≥ n 0 , Supp(Du n ) ⊆ {x ∈ R 2 | dist(x, Ext(Du 0 )) ≤ r} .
Combining Proposition 1.14 and (10) yields

Supp(Du * ) ⊆ lim inf n→+∞ Supp(Du n ) ⊆ lim sup n→+∞ Supp(Du n ) ⊆ Ext(Du 0 ) .
Structure-preserving convergence results. Considering the exact support recovery results mentioned for other reconstruction problems, one could wonder if stronger convergence guarantees can be derived. Namely, if u 0 is identi able and the sum of a few indicator functions of simple sets, are solutions of (P λ (y 0 + w)) made of the same number of atoms? And if so, are these atoms related to those appearing in the decomposition of u 0 ? In Part 2, we provide an answer to these questions by adapting the tools introduced in [START_REF] Duval | Exact Support Recovery for Sparse Spikes Deconvolution[END_REF].

N

In this section, we provide a quick overview of existing numerical methods for solving (P λ (y))2 . We refer the reader to the recent survey [Chambolle and Pock, 2021a] for a complete review of the topic. In the following, as we are interested in the practical resolution of (P λ (y)), we only consider the case of nitely many measurements, i.e. H = R m from some m ∈ N * .

The main idea behind most numerical methods is to introduce a xed spatial discretization and a discrete version of the total variation. To be more precise, one de nes a discretization parameter h, a nite dimensional space E h , and discrete versions of the measurement operator and the total variation Φ h : E h → H and TV h : E h → R, and looks for a solution of

inf u∈E h 1 2λ Φ h u -y 2 + TV h (u) . (11)
In most situations ( 11) is a non-smooth convex problem and can be solved e ciently using standard algorithms. The main di culty usually lies in the choice of a relevant discretization TV h of the total variation. In the remaining of this section, we brie y present a few standard choices for TV h . Our aim is to give an idea of why discretizing the total variation is challenging, and to underline the drawbacks which are common to most existing approaches.

Spatial discretization. As all solutions of (P λ (y)) have their support included in some common ball1 , one can equivalently solve (P λ (y)) in [-R, R] 2 (with Dirichlet boundary conditions) for a su ciently large R > 0. The most standard choice of spatial discretization is then to take a positive integer N , de ne h def.

= 2R/N , and choose E h to be the space of N by N matrices. Every element u = (u i,j ) (i,j)∈[1,N ] 2 of E h encodes the values of a piecewise constant function on the partition of [-R, R] 2 composed of squares of equal size. The resulting approximate solution of (P λ (y)) is hence piecewise constant on this pixel grid.

Finite di erences discretizations

In order to simplify the exposition, we only review the two most basic discretizations, and refer to [Chambolle and Pock, 2021a, Section 2] for more details.

For every matrix u = (u i,j

) (i,j)∈[1,N ] 2 ∈ E h we de ne ∂ h x u i,j def. = u i+1,j -u i,j , ∂ h y u i,j def. = u i,j+1 -u i,j ,
for all (i, j) ∈ [0, N ] 2 , with the convention u i,j = 0 if either i or j is in {0, N + 1}. We also de ne the discrete gradient

∇ h u i,j def. = ∂ h x u i,j , ∂ h y u i,j
Anisotropic total variation. The most basic idea to discretize the total variation is to dene TV h as follows:

TV h (u) def. = h N i=0 N j=0 |∇ h u i,j | 1 = h ∇ h u 1,1 .
With this choice TV h (u) = TV(ũ) with ũ the piecewise constant function naturally associated to u. The main drawback of this choice is that the functional TV h is strongly "anisotropic": it over-estimates the total variation of oblique discontinuities, and its minimization hence favors horizontal and vertical discontinuities. Figure 4 shows how this anisotropy induces a strong bias on the reconstructed images. In fact, one can prove this discretization is not even consistent, as it Γ-converges towards the perimeter associated with a crystalline total variation (see e.g. [START_REF] Chambolle | Continuous limits of discrete perimeters[END_REF]).

Original image Anistropic TV Isotropic TV

Figure 4 -Solutions of the discrete Rudin-Osher-Fatemi denoising problem using the anisotropic and isotropic total variations ( gure taken from [START_REF] Tabti | Symmetric Upwind Scheme for Discrete Weighted Total Variation[END_REF]).

Isotropic total variation. To avoid the anisotropy of the above-mentioned discretization, it is standard to consider the so-called isotropic discrete total variation, de ned by

TV h (u) def. = h N i=0 N j=0 |∇ h u i,j | 2 = h ∇ h u 2,1 .
It is arguably the most widely used discretization. As explained in [Chambolle and Pock, 2021a, Sections 2.1 and 2.2], it still has some anisotropy, since it is not invariation by a rotation of 90°of the input image. Reducing this anisotropy is however possible with slightly more sophisticated discretizations. Its main drawback (which is also shared by more complex discretizations) is that its minimization produces blurry images (see Figure 4 for an illustration of this phenomenon), while the (continuous) total variation is precisely used to recover sharp edges.

Dual-based discretizations

Another idea to de ne TV h is to consider the dual formulation of TV, given in (1), and to "discretize" the constraint φ ∞ ≤ 1 imposed on the dual variable. To be more precise, this amounts to de ning TV h as follows:

TV h (u) def. = sup h p, ∇ h u Lp 2,∞ ≤ 1 , ( 12
)
where L is a linear operator encoding the constraints to be imposed on the dual variable p.

Hintermüller, Rautenberg, Hahn and Condat's discretization. A discretization of the form ( 12) was introduced in two independant works [START_REF] Hintermüller | Functionalanalytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF], Condat, 2017]. The operator L is chosen in this case to enforce the constraint |p| ≤ 1 both at pixel centers and edges. The performance of this discrete total variation is acknowledged in [Chambolle and Pock, 2021a, Section 4.2]. The experiments presented in [Condat, 2017] indeed show reconstructions that are signi cantly less blurry than those obtained with the isotropic discretization, although recovered edges are never exactly sharp (i.e. the resulting images are not binary).

Learned dual discretization. In [Chambolle and Pock, 2021b], a learning approach is proposed to nd the best choice of operator L on a particular task. The consistency of the resulting discrete total variations is guaranteed: they all Γ-converge towards the continuous total variation1 . Interestingly, discretizations learned on di erent tasks signi cantly di er, which, according to the authors, show that it is hopeless to look for a universal best discrete total variation.

Motivation for grid-free approaches

We now motivate the numerical approach introduced in Part 3.

Common drawbacks of xed grid approaches. As mentioned above, most existing numerical methods for solving (P λ (y)) rely on the introduction of a xed spatial discretization. This often yield reconstruction artifacts, such as anisotropy or blur. Most importantly, they often fail to preserve the structure exhibited by (some) solutions of (P λ (y)), which are piecewise constant.

Adaptative discretizations. To circumvent these issues, mesh adaptation techniques were introduced (see e.g. [START_REF] Viola | A unifying resolutionindependent formulation for early vision[END_REF], Bartels et al., 2021]). The spatial discretization is hence adapted to the reconstructed image during the reconstruction process. The re nement rules these works propose are, however, either heuristic or too restrictive to faithfully recover edges. In any case, they still rely on a discretization of the whole domain, and hence do not provide a compact representation of the reconstruced image, which is highly desirable when working with simple images (i.e. that are the superposition of a few simple shapes).

A rst grid-free approach. In [START_REF] Ongie | O -the-Grid Recovery of Piecewise Constant Images from Few Fourier Samples[END_REF], a method for recovering piecewise constant images from few Fourier samples is introduced. Its orginiality is to produce a continuous domain representation of the image, assuming its edge set is a trigonometric curve. However, this approach heavily relies on relations satis ed by the Fourier coe cients of the image. As such, it does not seem possible to adapt it to handle other types of measurements. In some sense, this method can be seen as the counterpart of Prony-type methods used for the recovery of sparse spikes (see e.g. [Catala, 2020, Chapter 1]), which are speci c to certain type of measurements.

On the other hand, greedy approaches (that we try to adapt to our setting in Part 3) are agnostic to the choice of measurement operator.

Towards general grid-free numerical methods. Our goal is to design an algorithm which does not su er from some grid bias, while providing a continuous domain representation of solutions. To this aim, we propose to construct an approximate solution built from the atoms promoted by the total variation, namely indicator functions of simple sets. To obtain numerically tractable algorithms, we choose the simple sets we deal with to be simple polygons, although other choices could be considered. A schematic comparison of this approach and grid-based approaches is provided in Figure 5. We begin by collecting results about the total variation unit ball, and more speci cally about the structure of its (exposed) faces. The importance of such an analysis, hinted in Section 1.2.2, is showcased at then end, where we state that elements of a given face can be decomposed using a common set of atoms. Then, we turn to the analysis of the prescribed curvature problem. Its solutions are closely linked with elements of exposed faces, and our noise robustness analysis is intimately related to the stability of its minimizers. Finally, we introduce a so-called nondegenerate source condition, under which exact support recovery is obtained in a low noise regime. We end this part by a discussion on this condition, investigating whether it holds in a simple case.

Conclusion

P

Piecewise constant functions are the sparse objects associated to the total variation. They play the same role as sparse vectors for the 1 norm, and as discrete measures for the total variation of measures. However, these two kinds of objects are arguably easier to deal with, mainly because there is a canonical way to write them as the sum of a few atoms. Finding a good way to decompose piecewise constant functions is somewhat more intricate. This increased complexity is closely linked to properties of the total variation unit ball, whose faces have a non-trivial structure. In this section, we discuss these issues in detail, and precisely de ne the class of piecewise constant functions we aim at recovering, which we call M -simple.

E

Motivation. As a result of Proposition 1.5 (and recalling Remark 1.6), given (λ, w) ∈ R * + × H, we have that every nonzero solution u of (P λ (y

0 + w)) satis es ˆR2 η λ,w u TV(u) = 1 ,
with η λ,w the dual certi cate associated to this problem (de ned in Section 2.1.2). The solution set of (P λ (y 0 + w)) is hence included (up to a normalization) in the solution set of

sup u∈L 2 (R 2 ) ˆR2 η λ,w u s.t. TV(u) ≤ 1 , (13)
which is in fact the face of {TV ≤ 1} exposed by η λ,w (see [Rockafellar, 1970, Section 18] for more details on the notion of exposed face). This observation plays a crucial role in our analysis. In the following, we describe the structure of the exposed faces of the total variation unit ball, and derive a useful result allowing to represent its elements using a common set of atoms. This representation is at the core of the proof of our support recovery result. The following is based on the habilitation thesis of Vincent Duval [Duval, 2022] and on a forthcoming work [Boyer et al., in preparation] 1 . The results we present here are a subset of those given in [Duval, 2022]. In this last work, general linearly closed faces of the total variation unit ball are considered. We only treat the case of exposed faces, which is enough for our purpose. This allows slightly di erent (and often shorter) proofs.

Setting. Let us x η ∈ L 2 (R 2 ), and denote F the face of {TV ≤ 1} exposed by η, which is the solution set of sup

u∈L 2 (R 2 ) ˆR2 ηu s.t. TV(u) ≤ 1 . ( 14
)
The weak compactness of {TV ≤ 1} yields that F is non-empty (i.e. the supremum in ( 14) is attained). Moreover, if η is not equal to zero almost everywhere, it has a Lebesgue point x such that |η(x)| > 0. There hence exists r > 0 such that | ´B(x,r) η| > 0, which shows the value of ( 14) is strictly positive, say 1/α for some α ∈ R * + . Hence, we have that

∀u ∈ L 2 (R 2 ), TV(u) - ˆR2 (αη)u ≥ 0 ,
with equality if and only if u solves ( 14). For convenience reasons, we assume in the following that α = 1 (all the results below remain valid if α = 1 by replacing η with η def.

= αη), and consequently obtain

∀u ∈ L 2 (R 2 ), TV(u) - ˆR2 ηu ≥ 0 . ( 15
)
With this assumption we have the following characterization of F:

F = u ∈ L 2 (R 2 ) TV(u) ≤ 1 and ˆR2 η u = 1 ,
and every u ∈ F satis es ´R2 η u = TV(u).

Remark 2.1

Let us stress that Problem 14 has value 1 (i.e. α = 1) if and only if their exists a nonzero function u ∈ L 2 (R 2 ) which satis es ´R2 η u = TV(u) (i.e. η ∈ ∂TV(u)). Dual certi cates introduced above all fall in this category, unless the unique solution of (P 0 (y 0 )) (respectively (P λ (y))) is zero.

Indicator functions

From Proposition 1.2, we know that the extreme points of {TV ≤ 1} are the functions of the form ±1 E /P (E), with E a simple set such that 0 < |E| < +∞. Elements of F that are proportional to indicator functions hence play a special role in the analysis of its structure. Let us therefore de ne

E def. = E + ∪ E -∪ {∅, R 2 } , E + def. = E ⊂ R 2 |E| < +∞, 0 < P (E) < +∞, 1 E P (E) ∈ F , E -def. = E ⊂ R 2 |E c | < +∞, 0 < P (E c ) < +∞, -1 E c P (E c ) ∈ F . ( 16
)
Remark 2.2

We stress that an extreme point of {TV ≤ 1} is not necessarily exposed. In fact, E ∈ E + if and only if

P (E) = ´E η (and E ∈ E -if and only if P (E c ) = -´Ec η). This is equivalent to require that E (respectively E c
) is a solution of the prescribed curvature problem associated to η (respectively -η) studied in Section 2. In particular, this implies that E (respectively E c ) is bounded, that the number of Jordan curves in the decomposition of ∂ M E is nite, and that the following weak regularity property holds:

∃r 0 > 0, ∀r ∈ (0, r 0 ], ∀x ∈ ∂E, 1 16 ≤ |E ∩ B(x, r)| |B(x, r)| ≤ 1 - 1 16 .
We refer the reader to [Chambolle et al., 2016, Section 5] for more details.

The set E has the remarkable property of being closed under several operations that we describe below. We rst consider the intersection and union operations, before turning to M -connected components and holes.

Proposition 2.3

The set E is closed under countable unions and countable intersections.

Proof : If E ∈ E + and F ∈ E + the submodularity of the perimeter (see e.g. [START_REF] Ambrosio | Connected components of sets of nite perimeter and applications to image processing[END_REF], Proposition 1]) yields:

P (E ∩ F ) + P (E ∪ F ) ≤ P (E) + P (F ) = ˆE η + ˆF η = ˆE∩F η + ˆE∪F η .
We hence obtain:

P (E ∩ F ) - ˆE∩F η + P (E ∪ F ) - ˆE∪F η ≤ 0 .
Using (15), we get that the two terms above are nonnegative, which yields

E ∩ F ∈ E + (un- less |E ∩ F | = 0) and E ∪ F ∈ E + .
The same reasoning shows that the result also holds when replacing E + by E -.

If E ∈ E + and F ∈ F -we have:

P (E ∩ F ) + P ((E ∪ F ) c ) = P (E ∩ F ) + P (E ∪ F ) ≤ P (E) + P (F ) = P (E) + P (F c ) = ˆE η - ˆF c η = ˆE∩F η - ˆ(E∪F ) c η
Reasoning as above, we obtain

E ∩ F ∈ E + (unless |E ∩ F | = 0) and E ∪ F ∈ E -(un- less |(E ∪ F ) c | = 0), which shows E ∩ F ∈ E.
Now, we have proved that E is stable under nite unions and nite intersections. The countable case follows from the fact exposed faces are closed in

L 2 (R 2 ). Proposition 2.4 Let E, F ∈ E with F ⊂ E. If there exists C 1 , C 2 ⊂ R 2 such that 1. |C 1 | > 0 and |C 2 | > 0 2. E \ F = C 1 ∪ C 2 3. P (E \ F ) = P (C 1 ) + P (C 2 ) then F ∪ C 1 ∈ E and F ∪ C 2 ∈ E.
Proof : We rst note that C 1 R 2 and C 2 R 2 (otherwise we would have E = R 2 and F = ∅, which would contradict the indecomposability of R 2 ). Using 2. and 3., from [Ambrosio et al., 2001, Proposition 3] we hence obtain that H 1 (∂ * C 1 ∩ ∂ * C 2 ) = 0. Now, using Proposition 1 of the same reference, we obtain:

P (E) + P (F ) = 2 P (F ) + P (E \ F ) -2 H 1 (∂ * F ∩ ∂ * (E \ F )) = 2 P (F ) + P (E \ F ) -2 H 1 (∂ * F ∩ ∂ * (C 1 ∪ C 2 )) = (P (F ) + P (C 1 ) -2 H 1 (∂ * F ∩ ∂ * C 1 )) + (P (F ) + P (C 2 ) -2 H 1 (∂ * F ∩ ∂ * C 2 )) = P (G 1 ) + P (G 2 ) , with G 1 def. = F ∪ C 1 and G 2 def.
= F ∪ C 2 . We can then conclude as in the proof of Proposition 2.3.

As a consequence of Proposition 2.4, we obtain the following result. We remind the reader that M -connected components and holes are de ned in Appendix A.

Corollary 2.5 If E ∈ E, any M -connected component of E belongs to E. Moreover, for any hole Y of E, we have that E ∪ Y ∈ E.

Chains and maximal chains

Considering exposed faces (instead of general linearly closed faces) does not allow for shorter proofs of the results given in this section. We hence do not include them and refer the reader to [Duval, 2022, Section 2.3].

We say that a collection C of subsets of R 2 is a chain if for every E, F ∈ C we have E ⊆ F or F ⊆ E, i.e. if C is totally ordered for the inclusion relation. We call its cardinal (denoted |C|) the length of C. Denoting dim(F) the dimension of the a ne hull of F, we may state the following result.

Proposition 2.6

If C ⊆ E \ {∅, R 2 } is a chain and dim(F) = d, then C has length at most d + 1.

We say that a chain

C in E is maximal if there is no E ∈ E \ C such that C ∪ {E} is a chain. Such a chain must therefore contain ∅ and R 2 . As a result, if C = {E i } m i=0 is a maximal chain, then the collection of its increments {E i \E i-1 } m i=1 is a partition of R 2 .
The following proposition shows that every set in E can be written as the union of elements of this partition.

Proposition 2.7

Assume

C = {E i } m i=0 is a maximal chain in E with ∅ = E 0 ⊂ E 1 ⊂ ... ⊂ E m = R 2 .
Then the following holds.

1. For every i ∈ {1, ..., m}, the set

E i \ E i-1 is indecomposable. 2. For every F ∈ E, there exists I ⊆ {1, ..., m} such that F = i∈I (E i \ E i-1 ).
Even if there could be several maximal chains in F, they all have the same length and the same collection of increments, as the following proposition shows. In particular, we see that each exposed face F with nite dimension is naturally associated to a partition of R2 (which is the collection of increments of any maximal chain in F).

Proposition 2.8

If C = {E i } m i=0 and C = {E i } m i=0
are two maximal chains in E, we have:

1. |C| = |C | = dim(F) + 3 , 2. {E i \ E i-1 } m i=1 = {E i \ E i-1 } m i=1 .

C 1.2.1. De nition and properties

Let u ∈ L 2 (R 2 ) be a function with nite total variation, taking a nite number of values1

t 1 > ... > t i 0 -1 > t i 0 = 0 > t i 0 +1 > ... > t m . Let E i def. = {u ≥ t i } if 1 ≤ i ≤ m and E 0 def. = ∅. We hence have ∅ = E 0 ⊂ E 1 ⊂ ... ⊂ E m = R 2 , i.e. {E i } m i=0 is a chain. Now if F i def. = E i \ E i-1 = {u = t i } for every 1 ≤ i ≤ m, we have CC M (F i ) = {F i,j } j∈J i with J i at most countable. De ning J def. = {(i, j) | 1 ≤ i ≤ m, 1 ≤ j ≤ |J i |}, for every (i, j) ∈ J we introduce E i,j def. = E i-1   1≤k≤j F i,j   .
With these notations, the lexicographic order on J yields a total order on {E i,j } (i,j)∈J , which is hence a chain. We say that it is a chain associated to the function u (see Figure 6 for an illustration). By construction, u is constant on each of its increments (F i,j ) (i,j)∈J , which are all indecomposable.

Remark 2.9

Since the order we choose on CC M ({u = t i }) is arbitrary, multiple chains can be associated to u (see Figure 6 for an illustration). These chains however have the same collection of increments I, which is given by:

I = m i=1 CC M ({u = t i }) .
Considering the construction above, we now de ne the sparse objects naturally associated to the total variation, which we call M -simple functions. This is the class of piecewise constant functions we work with in the rest of this manuscript.

E 1 E 2 Y 2 Y 1 Figure 6 -A function u = 1 E 1 -1 E 2 . ∅ ∅ ∅ E1 E1 E1 E1 ∪ Y1 E c 2 \ (Y1 ∪ Y2) E c 2 \ (Y1 ∪ Y2) E1 ∪ Y1 ∪ Y2 E c 2 \ Y1 E c 2 \ Y2 E c 2 E c 2 E c 2 R 2 R 2 R 2
Figure 7 -Three chains associated to u (see Figure 6). Each column corresponds to a chain, whose elements are ordered from top to bottom.

De nition 2.10 (M -simple functions)

We say that a function u : R 2 → R is M -simple if there exists a nite set I, a collection {E i } i∈I of indecomposable sets of nite measure and t ∈ R I such that

u = i∈I t i 1 E i .
We stress that, with this de nition, any M -simple function u belongs to L 2 (R 2 ) and has nite total variation. We also have that a function u ∈ L 2 (R 2 ) with nite total variation is M -simple if and only if its associated chains have nite length. Moreover, the length of these chains is exactly the sum, over all values t taken by u, of the number of bounded M -connected components of {u = t}.

S

. The "right" measure of sparsity for a piecewise constant function u is the sum, over all values t taken by u, of the number of bounded M -connected components of {u = t}.

We conclude this section by stating the following useful result, whose importance is explained below.

Proposition 2.11

Assume that dim(F) = d and u ∈ F. Then u is M -simple. Moreover, any chain associated to u is a chain in E, and hence has at most d + 1 elements (not counting ∅ and R 2 ). Finally, u is constant on the increments of any maximal chain in E.

Proof : Using Carathéodory's theorem, we have that u can be written as a convex combination of d + 1 extreme points of F, which are of the form 1 E /P (E) with ∈ {-1, 1} and E a simple set such that 0 < |E| < +∞. As a consequence, u takes a nite number of values. Using the above-de ned notations we have:

u = i0-1 i=1 θ i 1 Ei + m-1 i=i0 θ i -1 E c i and ˆR2 η u = i0-1 i=1 θ i ˆEi η - m-1 i=i0 θ i ˆEc i η ,
with θ i > 0. By the coarea formula (see e.g. [Maggi, 2012, Chapter 13]), we also have:

i0-1 i=1 θ i P (E i ) + m-1 i=i0 θ i P (E c i ) = TV(u) .
Using ´R2 η u = TV(u) we hence obtain

P (E i ) = ´Ei η for 1 ≤ i ≤ i 0 -1 and P (E c i ) = -´Ec i η for i 0 ≤ i ≤ m -1, which shows {E i } 1≤i≤m-1 is a chain in E.
We conclude that {E i,j } (i,j)∈J is a chain in E by using Proposition 2.4. The last claim is a straightforward consequence of Proposition 2.7, since the level sets of u all belong to E.

The main interest of this proposition is that, denoting {F j } 1≤j≤m the collection of increments of any maximal chain in F, every function u ∈ F is of the form

u = m j=1 t j 1 F j ,
with t ∈ R m . To put it another way, the elements of F can be decomposed using a common set of atoms, which is the collection {1 F j } 1≤j≤m .

A su cient condition for identi ability

In this subsection, we derive a su cient identi ability condition using the above results. To this aim, we rst recall a condition, introduced in [START_REF] Burger | Convergence rates of convex variational regularization[END_REF] as the source condition, which ensures that a given function u 0 solves (P 0 (y 0 )).

De nition 2.12

A function u 0 satis es the source condition if there exists η ∈ Im Φ * such that η ∈ ∂TV(u 0 ).

If u 0 is an M -simple function, the source condition can be strengthened to ensure it is the unique solution of (P 0 (y 0 )) (i.e. u 0 is identi able). To this aim, for any M -simple function u and any chain C associated to u, whose ordered increments are denoted (F j ) 1≤j≤m , we de ne

Φ C : R m → H t → Φ   m j=1 t j 1 F j   .
We stress that this mapping only mildly depends on the choice of a chain C, in the sense that if C is another chain associated to u, then C and C have the same collection of increments. There hence exists a permutation of the coordinates σ such that

Φ C = Φ C • σ.
We may now state the following result, which is a natural analog of [Duval and Peyré, 2015, Proposition 5].

Proposition 2.13

Let u 0 be an M -simple function. Assume there exists η ∈ Im Φ * ∩ ∂TV(u 0 ). If a chain C 0 associated to u 0 is maximal in the face exposed by η, and if Φ C 0 is injective, then u 0 is the unique solution of (P 0 (y 0 )).

Proof : Since η ∈ Im Φ * and η ∈ ∂TV(u 0 ), by Proposition 1.4 we have that every solution u of (P 0 (y 0 )) satis es η ∈ ∂TV(u), and hence that u/TV(u) is in the face exposed by η. Since C 0 is a maximal chain in this face, denoting (F j ) 1≤j≤m its increments, we obtain that u/TV(u) is constant on each F j , which yields the existence of t ∈ R m such that

u = m j=1 t j 1 Fj .
The injectivity of Φ C0 nally allows to conclude.

Remark 2.14

We stress that a chain associated to an M -simple function u is maximal in an exposed face of {TV ≤ 1} if and only if all chains associated to u are maximal in this face. This directly follows from the fact all these chains have the same length. As a consequence, the result stated in Proposition 2.13 does not depend on the choice of a particular chain.

The main question which remains unanswered is the following: if u 0 is M -simple and identi able, with w and λ small enough, are solutions of (P λ (y 0 + w)) M -simple, what is the length of their associated chains, and how are their elements related to those of C 0 ? To answer an analog question for the sparse spikes problem, a non-degenerate version of the source condition was introduced in [Duval and Peyré, 2015, De nition 5].

From Proposition 1.3 and ( 8), the elements of the chains associated to solutions of (P λ (y 0 +w)) are all solutions of the prescribed curvature problem associated to η λ,w . In Section 2.1.2, we have also seen that, under a few assumptions, η λ,w converges to the minimal norm certi cate η 0 when w and λ go to zero. It is therefore natural to investigate how solutions of the prescribed curvature problem behave under variations of the curvature functional, which is the topic of the following section. The results therein allow us to state a natural analog of the non-degenerate source condition in Section 3, and to nally answer the above question.

T

As mentioned above, this section is dedicated to the study of the prescribed curvature problem associated to some function η ∈ ∂TV(0):

inf E⊂R 2 , |E|<+∞ J(E) def. = P (E) -ˆE η .
(PC(η))

Our main aim is to investigate how the solution set of (PC(η)) behaves when η varies. To be more speci c, given two su ciently close curvature functionals η and η we are interested in answering the following two questions.

(i) Are solutions of (PC(η )) close to solutions of (PC(η))?

(ii) How many solutions of (PC(η )) are there in a neighborhood of a given solution of (PC(η))?

In the following subsection, we answer the rst question using the notion of quasi-minimizers of the perimeter, as well as rst order optimality conditions for (PC(η)). Then, under a strict stability assumption on solutions of (PC(η)), we answer the second question using the notion of second order shape derivatives.

Tools. As discussed in the following subsection, if η is su ciently regular, the solutions of (PC(η)) also enjoy some regularity. In this section, we hence mainly work with smooth sets and extensively use related notions. Relevant de nitions and properties are collected in Appendix B. We also use the notion of curvature (de ned in Appendix A) and of second fundamental form (see e.g. [Maggi, 2012, Section 17.6]). Given a su ciently smooth set E, we often consider sets that are normal deformations of E, i.e. whose boundary is a normal graph over ∂E. Such sets are parametrized by real-valued functions on ∂E, which leads us to use the notion of tangential gradient, tangential Jacobian, and the spaces C k (∂E), L p (∂E) and H 1 (∂E) (along with their associated norms). We refer to the reader to [Henrot and Pierre, 2018, Sections 5.4.1, 5.4.3 and 5.9.1] for a precise de nition of these objects.

G

Existence of minimizers. We rst stress that existence of solutions is guaranteed for (PC(η)). Indeed, since η ∈ ∂TV(0), the objective J is always nonnegative, and equal to zero when evaluated at the empty set, which is admissible. From Proposition 1.3, we also know that it has a non trivial solution as soon as there exists a nonzero function u ∈ L 2 (R 2 ) such that η ∈ ∂TV(u).

Boundedness. As already mentioned, by Lemma 1.11 all solutions of (PC(η)) are included in some common ball, i.e. there exists R > 0 such that, for every solution E of (PC(η)), we have E ⊂ B(0, R).

Regularity of the solutions. Let us now discuss the regularity of solutions of (PC(η)). If we have η ∈ L ∞ (R 2 ), then any solution of (PC(η)) is a strong quasi-minimizer of the perimeter, and, consequently, is of class C 1,1 (see [Ambrosio, 2010, De nition 4.7.3 and Theorem 4.7.4]). If η is moreover continuous, then the boundary of any solution is locally the graph of a function u which solves (in the sense of distributions) the Euler-Lagrange equation associated to (PC(η)), which is (up to a translation and a rotation):

u 1 + u 2 (z) = u (z) 1 + u (z) 2 3/2 = η(z, u(z)) . (17)
This in turn implies that

u is C 2 (C k+2,α if η ∈ C k,α (R 2 ))
and solves (17) in the classical sense. As a result, assuming η is bounded and of class C 1 , every solution of (PC(η)) is of class C 3 .

Convergence result. Now, we state a result that, loosely speaking, tells that any neighborhood (in terms of C 2 -normal deformations) of the solution set of (PC(η 0 )) contains the solution set of (PC(η)) provided η is su ciently close to η 0 in C 1 (R 2 ) and L 2 (R 2 ). This provides an answer to question (i).

Proposition 2.15

Let η 0 ∈ ∂TV(0) ∩ C 1 b (R 2 ).
For every > 0 there exists r > 0 such that for every η ∈ ∂TV(0)

with η -η 0 L 2 (R 2 ) + η -η 0 C 1 (R 2 ) ≤ r, the following holds: every solution F of (PC(η)) is a C 2 -normal deformation of size at most of a solution E of (PC(η 0 )) (i.e., with the notation of Proposition B.5, F = E ϕ with ϕ C 2 (∂E) ≤ ).
Proof : We argue by contradiction and assume the existence of two sequences

(η n ) n∈N * and (F n ) n∈N * such that • for all n ∈ N * , η n ∈ ∂TV(0) ∩ C 1 b (R 2 ), • the sequence (η n ) n∈N * converges to η 0 in L 2 (R 2 ) and C 1 b (R 2 ),
• for all n ∈ N * , the set F n solves (PC(η n )) and cannot be written as a C 2 -normal deformation of size at most of a solution of (PC(η 0 )) .

We hence have that (F n ) n∈N * is bounded and that F n is a strong (Λ, r 0 )-quasi minimizer of the perimeter (in short form F n ∈ QM(Λ, r 0 ), see [Maggi, 2012, Section 21] and [Ambrosio, 2010, De nition 4.7.3]) for all n ∈ N * , with Λ = sup { η n ∞ , n ∈ N * } and r 0 any positive real number.

Taking r 0 small enough to have Λ r 0 ≤ 1, from [Maggi, 2012, Propositions 21.13 and 21.14] we get that, up to the extraction of a subsequence (not relabeled), (F n ) n≥0 converges in measure to a bounded set E ∈ QM(Λ, r 0 ), and that (∂F n ) n≥0 converges to ∂E. From |F n E| → 0 we obtain that E is a solution of (PC(η 0 )), and the convergence of (∂F n ) n≥0 towards E yields

∀r > 0, ∃n 0 ∈ N, ∀n ≥ n 0 , ∂F n ⊂ x∈∂E C(x, r, ν E (x)) ,
where C(x, r, ν E (x)) denotes the square of axis ν E (x) and side r centered at x, de ned in (60).

From [Ambrosio, 2010, 4.7.4], and arguing as in the proof of [Maggi, 2012, Theorem 26.6], for every x ∈ ∂E we obtain the existence of r > 0,

of n 0 ∈ N, of u ∈ C 1,1 ([-r, r]) and of a sequence (u n ) n≥n0 which is uniformly bounded in C 1,1 ([-r, r]), such that, in C(x, r, ν E (x))
, the set E is the hypograph of u and, for every n ≥ n 0 , the set F n is the hypograph of u n . Moreover, we have that u nu C 1 ([-r,r]) → 0. Now, we also have that u and u n (for n ≥ n 0 ) respectively solve (in the sense of distributions) the following equations in (-r, r):

u (z) (1 + u (z) 2 ) 3/2 = H(z, u(z)) , with H(z, t) def. = η 0 (x + R ν E (x) (z, t)) , u n (z) (1 + u n (z) 2 ) 3/2 = H n (z, u n (z)) , with H n (z, t) def. = η n (x + R ν E (x) (z, t)) . (18)
We hence immediately obtain that u and u n belong to C 2 ([-r, r]). Moreover, for every z ∈ (-r, r) we have:

|u n (z) -u (z)| = H n (z, u n (z)) 1 + u n (z) 2 3/2 -H(z, u(z)) 1 + u (z) 2 3/2 ≤ ( H n -H ∞ + |H(z, u n (z)) -H(z, u(z))|) 1 + u n (z) 2 3/2 + H ∞ 1 + u n (z) 2 3/2 -1 + u (z) 2 3/2 ,
from which we obtain that u nu ∞ → 0.

Using these new results in combination with (18), we get that u and u n belong to C 3 ([-r, r]). Di erentiating (18), we obtain, for every z ∈ (-r, r):

u (3) (z) = ∂ 1 H(z, u(z)) + u (z) ∂ 2 H(z, u(z)) (1 + u (z) 2 ) 3/2 + 3 H(z, u(z)) u (z) u (z) (1 + u (z) 2 ) 3/2 , u (3) n (z) = ∂ 1 H n (z, u n (z)) + u n (z) ∂ 2 H n (z, u n (z)) (1 + u n (z) 2 ) 3/2 + 3 H n (z, u n (z)) u n (z) u n (z) (1 + u n (z) 2 ) 3/2 ,
from which we can nally show u

(3) nu (3) ∞ → 0. Finally, using the compactness of ∂E, we obtain that (F n ) n≥0 converges in C 3 towards E, and Proposition B.3 allows to eventually write F n as a C 2 -normal deformation of E, whose norm converges to zero. This yields a contradiction.

S

Question (ii) is closely linked to the stability of minimizers of (PC(η)), that is to the behaviour of the objective J in a neighborhood of a solution. To analyze this behaviour, we use the general framework presented in [START_REF] Dambrine | Stability in shape optimization with second variation[END_REF], which relies on the notion of second order shape derivative.

Approach. The natural path to obtain our main stability result, which is Proposition 2.23, is to prove that J is in some sense of class C 2 , i.e. that its second order shape derivative is continuous at zero (see Proposition 2.21 for a precise statement). We could not nd this result in the literature. A large part of this section is hence dedicated to its proof. To the best of our knowledge, deriving such a result is necessary to obtain Proposition 2.23. In particular, we had to use a stronger condition than the "improved continuity condition" (IC H 1 ,W 2,∞ ) of [START_REF] Dambrine | Stability in shape optimization with second variation[END_REF], which is satis ed by our functional. This condition only requires some uniform control of second order directional derivatives at zero, which is weaker than the result of Proposition 2.21.

Structure of shape derivatives.

Given an open solution E of (PC(η)), we introduce the following mapping, where E ϕ denotes the normal deformation of E associated to ϕ, de ned in Proposition B.5:

j E : C 1 (∂E) → R ϕ → J(E ϕ ) .
With this notation, the following result holds.

Proposition 2.16 (See e.g. [Henrot and Pierre, 2018, Chapter 5])

If E is a bounded open set of class C 2 and η ∈ C 1 (R 2 )
, then j E is twice Fréchet di erentiable at 0 and, for every ψ ∈ C 1 (∂E), we have:

j E (0).(ψ) = ˆ∂E [η -H] ψ dH 1 j E (0).(ψ, ψ) = ˆ∂E |∇ τ ψ| 2 -H η + ∂η ∂ν ψ 2 dH 1
where H denotes the curvature of E and ∇ τ ψ def.

= ∇ψ -(∇ψ • ν) ν is the tangential gradient of ψ with respect to E.

From the expression of j E (0) and j E (0) given above, we immediately notice that j E (0) can be extended to a continuous linear form on L 1 (∂E), and j E (0) to a continuous bilinear form on H 1 (∂E).

Strict stability. Following [START_REF] Dambrine | Stability in shape optimization with second variation[END_REF], we say that a solution E of (PC(η)) is strictly stable if j E (0) is coercive in H 1 (∂E), i.e. if the following property holds:

∃α > 0, ∀ψ ∈ H 1 (∂E), j E (0).(ψ, ψ) ≥ α ψ 2 H 1 (∂E) .
Under a few assumptions (that are satis ed by our functional), this strict stability condition ensures that E is a strict local minimizer of J (see Theorem 1.1 in the above-mentioned reference), and is hence the only minimizer (modulo Lebesgue negligible sets) among sets E ϕ with ϕ in a neighborhood of 0. It plays a crucial role in our answer to question (ii).

Continuity results. Now, we prove a few results concerning the convergence of j E towards j 0,E and the continuity of ϕ → j E (ϕ), where j E and j 0,E are the functionals respectively associated to η and η 0 . To achieve this, we need to compute j E (ϕ) for ϕ ∈ C1 (∂E) in a neighborhood of 0. This may be done using Lemma 2.17 below. To state it, given a bounded open set E of class C 2 and ϕ in a neighborhood of 0 in C 1 (∂E), we introduce the mapping f ϕ = Id + ξ ϕ , with ξ ϕ de ned as in Lemma B.4. If ϕ C 1 (∂E) is su ciently small then f ϕ is a C 1 -di eomorphism, and we denote its inverse by g ϕ .

Lemma 2.17

Let E be a bounded open set of class C 2 . Then for every ϕ in a neighborhood of 0 in C 1 (∂E), and for every ψ ∈ H 1 (∂E), we have:

j E (ϕ).(ψ, ψ) = j Eϕ (0).(ξ ψ • g ϕ • ν ϕ , ξ ψ • g ϕ • ν ϕ ) + j Eϕ (0).(Z ϕ,ψ ) (19)
where ν ϕ is the unit outward normal to E ϕ and

Z ϕ,ψ = B ϕ ((ξ ψ • g ϕ ) τϕ , (ξ ψ • g ϕ ) τϕ ) -2(∇ τϕ (ξ ψ • g ϕ • ν ϕ )) • (ξ ψ • g ϕ ) τϕ ,
with ζ τϕ and ∇ τϕ ζ the tangential part and the tangential gradient of ζ with respect to E ϕ , and B ϕ the second fundamental form of E ϕ .

Proof : To prove this result, we need to introduce J E 1 de ned by

J E : C 1 b (R 2 , R 2 ) → R ξ → J((Id + ξ)(E)) .
We denote ν the outward unit normal to E and B its second fundamental form. We also denote ζ τ and ∇ τ ζ the tangential part and the tangential gradient of ζ with respect to E. The structure theorem (see e.g. [Henrot and Pierre, 2018, Theorem 5.9.2] or [Dambrine and Lamboley, 2019, Theorem 2.1]) then yields, for every su ciently smooth vector eld ζ:

J E (0).(ζ) = j E (0).(ζ ∂E • ν) , J E (0).(ζ, ζ) = j E (0).(ζ ∂E • ν, ζ ∂E • ν) + j E (0).(Z ζ ) ,
where

Z ζ def. = B(ζ τ , ζ τ ) -2 (∇ τ (ζ • ν)) • ζ τ .
Now, we rst notice that, for every pair of vector elds ξ, ζ such that Id + ξ is invertible, we have:

(Id + ξ + ζ)(E) = (Id + ζ • (Id + ξ) -1 )((Id + ξ)(E)) .
De ning

F def. = (Id + ξ)(E) we hence obtain J E (ξ + ζ) = J F (ζ • (Id + ξ) -1 , ζ • (Id + ξ) -1 ), which yields J E (ξ).(ζ, ζ) = J F (0).(ζ • (Id + ξ) -1 ) .
Using this with ξ = ξ ϕ and ζ = ξ ψ , we get:

j E (ϕ).(ψ, ψ) = J E (ξ ϕ )(ξ ψ , ξ ψ ) = J Eϕ (0).(ξ ψ • g ϕ , ξ ψ • g ϕ ) ,
and we nally obtain (19) by applying the structure theorem.

Most of the results below rely on the following technical lemma, whose rst part is contained in [Dambrine and Lamboley, 2019, Lemma 4.7].

Lemma 2.18

Let E be a bounded C 2 set. If ϕ C 1 (∂E) → 0 we have:

(i) f ϕ -Id C 1 (∂E) → 0 , ν ϕ • f ϕ -ν C 0 (∂E) → 0 , (iii) (ii) g ϕ -Id C 1 (∂Eϕ) → 0 , Jac τ f ϕ -1 C 0 (∂E) → 0 . (iv)
If ϕ C 2 (∂E) → 0 then we also have:

(v) H ϕ • f ϕ -H C 0 (∂E) → 0 , B ϕ • f ϕ -B C 0 (∂E) → 0 . (vi)
Moreover, the following holds:

(a) lim

ϕ C 1 (∂E) →0 sup ψ∈L 2 (∂E)\{0} (ξ ψ • g ϕ ) τϕ L 2 (∂Eϕ) ψ L 2 (∂E) = 0 , (b) lim ϕ C 1 (∂E) →0 sup ψ∈H 1 (∂E)\{0} ∇ τϕ (ξ ψ • g ϕ • ν ϕ ) L 2 (∂Eϕ) -∇ τ ψ L 2 (∂E) ψ H 1 (∂E) = 0 , (c) lim ϕ C 2 (∂E) →0 sup ψ∈H 1 (∂E)\{0} Z ϕ,ψ L 1 (∂Eϕ) ψ 2 H 1 (∂E)
= 0 .

Proof : See [Dambrine and Lamboley, 2019, Lemma 4.7] for a proof of the results stated in the rst part of the lemma. To prove (a) we use the fact that

(ξ ψ • g ϕ ) τϕ 2 L 2 (∂Eϕ) = ˆ∂E (ν • g ϕ ) 2 τϕ • f ϕ Jac τ f ϕ ψ 2 dH 1 ≤ Jac τ f ϕ C 0 (∂E) (ν • g ϕ ) τϕ • f ϕ 2 C 0 (∂E) ψ 2 L 2 (∂E) = Jac τ f ϕ C 0 (∂E) ν -(ν • ν ϕ • f ϕ ) ν ϕ • f ϕ 2 C 0 (∂E) ψ 2 L 2 (∂E) .
which, using (i), (iii) and (iv), yields the result.

To prove (b), we notice that:

∇ τϕ (ξ ψ • g ϕ • ν ϕ ) = c 1 ϕ ψ • g ϕ + c 2 ϕ • ∇ τ ψ • g ϕ τ ϕ , with τ = ν ⊥ , τ ϕ = ν ⊥ ϕ 1 and c 1 ϕ def. = τ • g ϕ • ν ϕ (J gϕ τ ϕ ) • τ • g ϕ + τ ϕ • ν • g ϕ , c 2 ϕ def. = ν • g ϕ • ν ϕ (J gϕ τ ϕ ) .
We hence obtain

|∇ τϕ (ξ ψ • g ϕ • ν ϕ ) • f ϕ Jac τ f ϕ -∇ τ ψ| ≤ c ϕ (|ψ| + |∇ τ ψ|)
with c ϕ independant of ψ. Moreover, using (ii) and (iii), we have:

lim ϕ C 1 (∂E) →0 c ϕ C 0 (∂E) → 0 .
1 These two vectors are de ned as the application of the rotation of angle π/2 to ν and νϕ.

Denoting A def. = ∇ τϕ (ξ ψ • g ϕ • ν ϕ ) L 2 (∂Eϕ) -∇ τ ψ L 2 (∂E) , this nally yields A ≤ ∇ τϕ (ξ ψ • g ϕ • ν ϕ ) • f ϕ Jac τ f ϕ -∇ τ ψ L 2 (∂E) ≤ √ 2 c ϕ C 0 (∂E) ψ H 1 (∂E) .
We now prove (c). Since

B ϕ ((ξ ψ • g ϕ ) τϕ , (ξ ψ • g ϕ ) τϕ ) L 1 (∂Eϕ) ≤ B ϕ C 0 (∂Eϕ) (ξ ψ • g ϕ ) τϕ 2 L 2 (∂Eϕ) and B ≤ ∇ τϕ (ξ ψ • g ϕ • ν ϕ )) L 2 (∂Eϕ) (ξ ψ • g ϕ ) τϕ L 2 (∂Eϕ) with B def. = (∇ τϕ (ξ ψ • g ϕ • ν ϕ )) • (ξ ψ • g ϕ ) τϕ L 1 (∂Eϕ)
, we get the result.

Using the above result, we now prove the continuity of ϕ → j E (ϕ) by proving the continuity of the two terms appearing in its expression. In all the following, if E is a (real) vector space, we denote Q(E) the set of quadratic forms over E.

Proposition 2.19

If E is a bounded C 2 set and p E : ϕ → P (E ϕ ), the mapping

p E : C 2 (∂E) → Q(H 1 (∂E)) ϕ → p E (ϕ)
is continuous at 0.

Proof : Using Lemma 2.18, for every ϕ ∈ C 2 (∂E) in a neighborhood of 0 and ψ ∈ H 1 (∂E), we obtain:

p E (ϕ).(ψ, ψ) -p E (0).(ψ, ψ) = A + p Eϕ (0).(Z ϕ,ψ ) , with A def. = p Eϕ (0).((ξ ψ • g ϕ ) • ν ϕ , (ξ ψ • g ϕ ) • ν ϕ ) -p E (0).(ψ, ψ)
. Now, we also have:

A = ∇ τϕ (ξ ψ • g ϕ • ν ϕ ) 2 L 2 (∂Eϕ) -∇ τ ψ 2 L 2 (∂E) ,
and using Lemma 2.18 we obtain

lim ϕ C 2 (∂E) →0 sup ψ∈H 1 (∂E)\{0} p Eϕ (0).((ξ ψ • g ϕ ) • ν Eϕ , (ξ ψ • g ϕ ) • ν Eϕ ) -p E (0).(ψ, ψ) ψ 2 H 1 (∂E) = 0 . Moreover |p E (0).(Z ϕ,ψ )| ≤ H ϕ L ∞ (∂Eϕ) Z ϕ,ψ L 1 (∂Eϕ) ,
and Lemma 2.18 allows to conclude.

Proposition 2.20

If E is a bounded C 2 set, η ∈ C 1 b (R 2 )
and g E : ϕ → ´Eϕ η, the mapping

g E : C 2 (∂E) → Q(H 1 (∂E)) ϕ → g E (ϕ)
is continuous at 0.

Proof : We proceed as in Proposition 2.19. De ning

A def. = g Eϕ (0).((ξ ψ • g ϕ ) • ν Eϕ , (ξ ψ • g ϕ ) • ν Eϕ )
we have:

A = ˆ∂Eϕ H ϕ η + ∂η ∂ν ϕ ((ψ ν) • g ϕ • ν Eϕ ) 2 dH 1 = ˆ∂E H ϕ η + ∂η ∂ν ϕ • f ϕ (ν • ν ϕ • f ϕ ) 2 Jac τ f ϕ ψ 2 dH 1 .
This yields:

g Eϕ (0).((ξ ψ • g ϕ ) • ν ϕ , (ξ ψ • g ϕ ) • ν ϕ ) -g E (0).(ψ, ψ) ψ 2 L 2 (∂E) ≤ c ϕ , with c ϕ def. = H ϕ η + ∂η ∂ν ϕ • f ϕ (ν • ν ϕ • f ϕ ) 2 Jac τ f ϕ -H η + ∂η ∂ν ∞ .
Using Lemma 2.18 we obtain

lim ϕ C 2 (∂E)→0 sup ψ∈H 1 (∂E)\{0} g Eϕ (0).((ξ ψ • g ϕ ) • ν ϕ , (ξ ψ • g ϕ ) • ν ϕ ) -g E (0).(ψ, ψ) ψ 2 H 1 (∂E) = 0 . Moreover |g Eϕ (0).(Z ϕ,ψ )| ≤ η ∞ Z ϕ,ψ L 1 (∂Eϕ) ,
and using again Lemma 2.18 we nally obtain the result.

As a consequence of the two propositions above, we obtain:

Proposition 2.21 If E is a bounded C 2 set and η ∈ C 1 b (R 2
), the mapping

j E : C 2 (∂E) → Q(H 1 (∂E)) ϕ → j E (ϕ)
is continuous at 0. Now, we prove that for ϕ in a neighborhood of

0 in C 1 (∂E), the mapping j E (ϕ) is uniformly close to j 0,E (ϕ) provided η -η 0 C 1 (R 2 ) is small engouh. Proposition 2.22 Let E be a bounded C 2 set and η 0 ∈ C 1 b (R 2 ). There exists > 0 such that lim η-η 0 C 1 (R 2 ) →0 sup ϕ C 2 (∂E) ≤ j E (ϕ) -j 0,E (ϕ) Q(H 1 (∂E)) = 0 . Proof : Since |(j E -j 0,E ) (ϕ).(ψ, ψ)| ≤ c 1 ϕ + c 2 ϕ with c 1 ϕ def. = ˆ∂Eϕ H ϕ (η -η 0 ) + ∂ (η -η 0 ) ∂ν ϕ (ξ ψ • g ϕ • ν ϕ ) 2 dH 1 , c 2 ϕ def. = ˆ∂Eϕ (η -η 0 ) Z ϕ,ψ dH 1 ,
the result readily follows from Lemma 2.18.

Stability result. We are now able to state the nal result of this section, which, loosely speaking, states that if E is a strictly stable solution of (PC(η 0 )), there is at most one ϕ in a neighborhood of 0 such that E ϕ is a solution of (PC(η)), provided η -

η 0 C 1 (R 2 ) is small engouh. Proposition 2.23 Let η 0 ∈ ∂TV(0) ∩ C 1 b (R 2
) and E be a strictly stable solution of (PC(η 0 )). Then there exists > 0 and r > 0 such that for every η ∈ ∂TV(0) with ηη 0 C 1 (R 2 ) ≤ r there is at most one ϕ ∈ C 2 (∂E) such that ϕ C 2 (∂E) ≤ and E ϕ solves (PC(η)).

Proof :

The fact E is a strictly stable solution of (PC(η 0 )) and Propositions 2.21 and 2.22 give the existence of > 0, r > 0 and α > 0 such that, for every

(ϕ, η) ∈ C 2 (∂E) × C 1 b (R 2 ) with ϕ C 2 (∂E) ≤ and η -η 0 C 1 (R 2 ) ≤ r, we have: sup ψ∈H 1 (∂E)\{0} j E (ϕ).(ψ, ψ) ψ 2 H 1 (∂E) ≥ α
As a result, j E (ϕ) is coercive (and hence positive de nite) for every ϕ such that ϕ C 2 (∂E) ≤ .

We therefore obtain that j E is strictly convex on this set and the result follows.

Summary. Combining the results of Propositions 2.15 and 2.23, we have proved that, provided η is su ciently close to η 0 in C 1 (R 2 ) and L 2 (R 2 ), every solution of (PC(η)) belongs to a neighborhood (in terms of C 2 -normal deformations) of a solution of (PC(η 0 )), and that, under a strict stability assumption, each of these neighborhoods contain at most one solution of (PC(η)).

In Theorem 2.27 below, we prove (under suitable assumptions) that, if η = η λ,w is the dual certicate associated to (P λ (y 0 + w)) and η 0 the minimal norml dual certi cate associated to (P 0 (y 0 )), then each neighborhood of a solution of (PC(η 0 )) contains exactly one solution of (PC(η λ,w )).

E

In this section, we provide an answer to the following question: if u 0 is M -simple and identi able, are solutions of (P λ (y 0 + w)) M -simple, what is the length of their associated chains, and how are they related to the chains associated to u 0 ? To answer it, we rst use the results proved in Section 2 to show that, under a few assumptions, the dimension of the faces of the total variation unit ball is in some sense stable. Then, we introduce a non-degenerate version of the source condition, and nally prove our exact support recovery result.

S

If η ∈ ∂TV(0) ∩ C 1 b (R 2 ) and F is the face exposed by η, de ning E as in ( 16), we say that E ⊂ R 2 is a strictly stable element of E if E is a strictly stable solution of (PC(η)) or E c is a strictly stable solution of (PC(-η)).

Theorem 2.24

Let

η 0 ∈ ∂TV(0) ∩ C 1 b (R 2 ).
Assume the face F 0 exposed by η 0 has a maximal chain C 0 whose elements are all strictly stable. Then for every > 0 there exists r > 0 such that, for every η ∈ ∂TV(0

) ∩ C 1 b (R 2 ) with η -η 0 L 2 (R 2 ) + η -η 0 C 1 (R 2 ) ≤ r ,
the face F exposed by η has a maximal chain whose elements are C 2 -normal deformations of size at most of elements of C 0 . In particular dim(F) ≤ dim(F 0 ).

Remark 2.25

To be more precise, Theorem 2.24 states that, if

C 0 = {E i } m i=0 , then F has a maximal chain C of length n ≤ m, say C = {F j } n j=0 , with ∀j ∈ {1, ..., n -1}, F j = E θ(j) ϕ j ,
where θ : {1, ..., n -1} → {1, ..., m -1} is a strictly increasing function and ∀j ∈ {1, ..., n -1}, ϕ j C 2 (∂E θ(j) ) ≤ .

Proof : We argue by contradiction and assume the existence of a sequence (η k ) k∈N * converging in L 2 (R 2 ) and C 1 (R 2 ) to η 0 . We denote by E 0 the set associated to F 0 de ned as in ( 16) (and similarly by E k the set associated to the face F k exposed by η k ). For simplicity, we rst treat the case where, for every i ∈ {1, ..., m -1}, the set E i belongs to E + 0 (i.e. E i is a solution of (PC(η 0 ))). The general case, where E i belongs to E - 0 (i.e. E c i is a solution of (PC(-η 0 ))) for every i ≥ i 0 + 1 with i 0 possibly strictly smaller than m -1, is discussed at the end of the proof.

From Propositions 2.15 and 2.23 we know that there exists > 0 and k 1 ∈ N * such that, for every k ≥ k 1 , every solution of (PC(η k )) is in a neighborhood of size (in terms of C 2 -normal deformations) of a solution of (PC(η 0 )), and that each of these neighborhoods contains at most one of these solutions. We hence have that, for every i ∈ {1, ..., m -1}, there exists in nitely many k such that the neighborhood of size of E j contains exactly one solution of (PC(η k )), or in nitely many k such that this neighborhood does not contain any solution of (PC(η k )). We therefore obtain the existence of a partition (I, J) of {1, ..., m -1} such that, up to extraction and possibly increasing k 1 , for every k ≥ k 1 the following holds:

• if i ∈ I then there is no solution of (PC(η k )) in the neighborhood of E i of size ,

• if i ∈ J there is exactly one solution of (PC(η k )) in the neighborhood of E i of size .

There hence exists n = |J| + 1 ≤ m and an increasing bijection θ : {1, ..., n -1} → J with the following property: for every k ≥ k 1 and j ∈ {1, ..., n -1}, there exists ϕ k,j such that F k,j

def. = E θ(j) ϕ k,j
is the unique solution of (PC(η k )) in the neighborhood of E θ(j) of size .

Now, let us show that there exists k 2 ∈ N such that, for every k ≥ k 2 , the collection

C k def. = {F k,j } 0≤j≤n with F k,0 def. = ∅ and F k,n def.
= R 2 is a chain. For every j ∈ {1, ..., n -2}, we have that (F k,j ∩ F k,j+1 ) k≥k1 solves (PC(η k )) and converges in measure to E θ(j) ∩ E θ(j+1) = E θ(j) . Arguing as in the proof of Proposition 2.15, we obtain the existence of k 2,j such that F k,j ∩ F k,j+1 is in the neighborhood of E θ(j) of size for every k ≥ k 2,j . Using Proposition 2.23 then yields F k,j ∩ F k,j+1 = F k,j for every k ≥ k 2,j .

Repeatedly applying this argument we obtain the existence of k 2,1 , ..., k 2,n-2 and de ning k 2 as the maximum of these integers we get that C k is a chain for every k ≥ k 2 . Now, let us show the existence of k 3 ∈ N * such that, for every k ≥ k 3 , the chain C k is maximal in the face F k exposed by η k . Arguing by contradiction, we assume the existence of

(G k ) k≥k2 such that G k ∈ E k \ C k and C k ∪ {G k } is a chain for every k ≥ k 2 1 .
Arguing as in the proof of Proposition 2.15, we get the existence of G ∈ E 0 such that, up to the extraction of a subsequence that we do not relabel,

G k = G ψ k for k large enough, with ψ k C 2 (∂G) → 0. Now, since C k ∪ {G k } is a chain for every k ≥ k 2 ,
we have that, for every j ∈ {1, n -1}, there exists in nitely many k such that F k,j ⊂ G k or in nitely many k such that G k ⊂ F k,j . Up to another extraction, we hence get the existence of j ∈ {0, ..., n -1} and

k 3 ∈ N * such that F k,j ⊂ G k ⊂ F k,j+1 for every k ≥ k 3 . As a result, we obtain      ∅ ⊂ G ⊂ E θ(1) if j = 0 , E θ(n-1) ⊂ G ⊂ R 2 if j = n -1 , E θ(j) ⊂ G ⊂ E θ(j+1) otherwise.
Using the maximality of C 0 , we get

G = E i with i ∈ I. Since G k = G ψ k with ψ k C 2 (∂G) → 0
and, for every k ≥ k 1 , su ciently small neighborhoods of E i do not contain any solution of (PC(η k )), we get a contradiction.

Finally, let us comment on the general case, where there exists i 0 ∈ {0, ..., m -1} such that, for every i ∈ {1, ..., m -1},

E i ∈ E + 0 if i ≤ i 0 and E i ∈ E - 0 if i ≥ i 0 + 1.
The same arguments can be applied to obtain the result, bearing in mind that if i ≥ i 0 + 1, then E c i is a solution of (PC(-η 0 )) (instead of E i being a solution of (PC(η 0 ))).

M

We are now able to introduce a non-degenerate version of the source condition, which ultimately allows us to state our main result.

De nition 2.26 (Non-degenerate source condition)

Let u 0 be an M -simple function and C 0 a chain associated to u 0 . We say that u 0 satis es the non-degenerate source condition if

1. Φ C 0 has full rank 2. η 0 ∈ ∂TV(u 0 ),
3. C 0 is maximal in the face exposed by η 0 , 4. every element of C 0 is strictly stable, with η 0 the minimal norm certi cate. In that case, we say that η 0 is non-degenerate.

Let us stress that, as pointed out in Remark 2.14, De nition 2.26 does not depend on the choice of a chain C 0 associated to u 0 , as one is maximal in an exposed face if and only if they all are.

Theorem 2.27 Let u 0 be an M -simple function and

C 0 = {E i } m i=0 an associated chain. Assume Φ * is contin- uous from H to C 1 b (R 2
) and u 0 satis es the non-degenerate source condition. Then there exists constants α, λ 0 ∈ R * + such that for every (λ, w) ∈ R * + × H with λ ≤ λ 0 and w H /λ ≤ α, any solution u λ,w of (P λ (y)) is M -simple and its associated chains have the same length as C 0 .

Moreover, writing

u 0 = m i=1 t i 1 E i \E i-1 ,
we have

u λ,w = m i=1 t λ,w i 1 E λ,w i \E λ,w i-1 , ( 20
)
with E λ,w

0 def. = ∅, E λ,w m def.
= R 2 , and

∀i ∈ {1, ..., m -1}, E λ,w i = (E i ) ϕ λ,w i with ϕ λ,w i ∈ C 2 (∂E i ) . (21)
Finally, if λ = w H /α, we have:

∀i ∈ {1, ..., m}, lim w→0 t λ,w i = t i , ∀i ∈ {1, ..., m -1}, lim w→0 ϕ λ,w i C 2 (∂E i ) = 0 .
Proof : Let and r be such that the results of Theorem 2.24 hold. Since Φ * is continuous from H to C 1 b (R 2 ), using Proposition 1.8 and ( 9) we obtain the existence of α and λ 0 such that for every (λ, w) ∈ R * + × H with λ ≤ λ 0 and w H /λ ≤ α 0 , we have:

η λ,w -η 0 L 2 (R 2 ) + η λ,w -η 0 C 1 (R 2 ) ≤ r .
Since C 0 is a maximal chain in the face exposed by η 0 , Theorem 2.24 ensures the existence of a maximal chain C λ,w in the face exposed by η λ,w , whose length is not greater than the length of C 0 . As a result, u λ,w is M -simple. We also know the elements of C λ,w are C 2 -normal deformations of elements of C 0 . Using the fact u λ,w is constant on each increment of C λ,w and converges in L 1 (R 2 ) towards u 0 , we obtain that C 0 and C λ,w have equal length, and that (20) and ( 21) hold.

A B u 0 A λ,w B λ,w u λ,w
Finally, the convergence of ϕ λ,w i towards 0 in C 2 (∂E i ) follows from Theorem 2.24, and the convergence of t λ,w i from the convergence of u λ,w towards u 0 in L 1 (R 2 ).

A

The core assumption in Theorem 2.27 is the strict stability of the elements of the chains associated to u 0 , as minimizers of the prescribed curvature functional associated to the minimal norm certi cate η 0 . In this section, we study this strict stability assumption and derive a natural su cient condition for it to hold. We then discuss to what extent this condition is necessary.

Setting. We x η ∈ ∂TV(0)∩C 1 b (R 2
) and E a (non-trivial) solution of the prescribed curvature problem (PC(η)) associated to η. We recall that E is a strictly stable solution of (PC(η)

) if j E (0) is coercive in H 1 (∂E), with ∀ψ ∈ H 1 (∂E), j E (0).(ψ, ψ) = ˆ∂E |∇ τ E ψ| 2 -H E η + ∂η ∂ν E ψ 2 dH 1 .
Equivalence of coercivity and positive de niteness. As explained (in a more general context) in [START_REF] Dambrine | Stability in shape optimization with second variation[END_REF], under a few assumptions, the bilinear form j E (0) is in fact coercive if and only if it is positive de nite. Instead of proving that our functional J ts the assumptions of Lemma 3.1 in the above reference, we use their arguments to provide a direct proof in our speci c setting.

Lemma 2.28 ([Dambrine and Lamboley, 2019, Lemma 3.1])

Let η ∈ C 1 (R 2
) and E be a bounded open set of class C 2 . Then the following propositions are equivalent:

(i) j E (0) is positive de nite, i.e. ∀ψ ∈ H 1 (∂E) \ {0}, j E (0).(ψ, ψ) > 0 , (ii) j E (0) is coercive, i.e. ∃α > 0, ∀ψ ∈ H 1 (∂E), j E (0).(ψ, ψ) ≥ α ψ 2 H 1 (∂E)
.

Proof : The fact (ii) implies (i) is trivial. We assume (i) and let (ψ n ) n≥0 be a minimizing sequence for inf

ψ∈H 1 (∂E) j E (0).(ψ, ψ) s.t. ψ H 1 (∂E) = 1 . (22)
Up to the extraction of a subsequence (not relabeled), we have that (ψ n ) n≥0 converges weakly in H 1 (∂E) and (by the compactness of the embedding H 1 (∂E) ⊂⊂ L 2 (∂E)) strongly in L 2 (∂E) to some ψ ∈ H 1 (∂E). We hence obtain

lim n→+∞ ˆ∂E H E η + ∂η ∂ν E ψ 2 n dH 1 = ˆ∂E H E η + ∂η ∂ν E ψ 2 dH 1 .
Let us now distinguish two cases.

1. If ψ = 0, then we use

ˆ∂E |∇ τ E ψ| 2 dH 1 ≤ lim inf n→+∞ ˆ∂E |∇ τ E ψ n | 2 dH 1 , which, using (i), yields lim inf n→+∞ j E (0).(ψ n , ψ n ) ≥ j E (0).(ψ, ψ) > 0.
2. If ψ = 0, using that ψ n H 1 (∂E) = 1 for all n, we obtain lim inf

n→+∞ j E (0).(ψ n , ψ n ) = lim inf n→+∞ ˆ∂E |∇ τ E ψ n | 2 dH 1 = 1 > 0 .
As a result, we obtain that, in both cases, the in mum in ( 22) is strictly positive, which shows that j E (0) is coercive.

A su cient condition for coercivity. Using the expression of j E (0) and the equivalence between coercivity and positive de niteness, the following result can be directly obtained.

Proposition 2.29 Let η ∈ C 1 (R 2 ) and E be a bounded open set of class C 2 such that j E (0) = 0. If sup x∈∂E H E (x) 2 + ∂η ∂ν E (x) < 0 , then j E (0) is coercive.
Necessity of the condition? A natural question to investigate is whether the condition given in Proposition 2.29 is necessary. A rst result in this direction is given in Proposition 2.30 below. Before stating it, let us de ne λ 1 (Γ), the rst Dirichlet eigenvalue of the Laplacian associated to some simple C 1 curve Γ with nite length (see for instance [START_REF] Kuttler | Eigenvalues of the Laplacian in Two Dimensions[END_REF] for the more classical case of open bounded sets):

λ 1 (Γ) def.
= inf

ψ∈H 1 0 (Γ)\{0} ∇ τ Γ ψ 2 L 2 (Γ) ψ 2 L 2 (Γ) . ( 23
)
The in mum in ( 23) is attained and is actually equal to the Dirichlet eigenvalue of the interval

I = (0, H 1 (Γ)) ⊂ R, which is 1/(πH 1 (Γ) 2 ).
To see this, one simply needs to consider an arc-length parameterization γ of Γ, and notice that, for every ψ ∈ H 1 0 (Γ) \ {0}, we have:

∇ τ Γ ψ 2 L 2 (Γ) ψ 2 L 2 (Γ) = ´I (ψ • γ) (t) 2 dt ´I (ψ • γ)(t) 2 dt .
We can now state the following proposition.

Proposition 2.30

If there exists α > 0 such that

H 2 E + ∂η ∂ν E ≥ α on a connected subset Γ of ∂E with α H 1 (Γ) 2 ≥ 1/π ,
then j E (0) is not coercive.

Proof : Since the in mum in the de nition of λ 1 (Γ) is attained, we have the existence of a nonzero

function ϕ ∈ H 1 0 (Γ) such that ∇ τΓ ϕ 2 L 2 (Γ) ϕ 2 L 2 (Γ) = λ 1 (Γ) = 1 π H 1 (Γ) 2 ≤ α . We hence obtain ˆΓ |∇ τΓ ϕ| 2 -H 2 E + ∂η ∂ν E ϕ 2 dH 1 ≤ ˆΓ |∇ τΓ ϕ| 2 -α ϕ 2 dH 1 ≤ 0 .
We can then extend ϕ to ψ ∈ H 1 (∂E) whose support is compactly included in Γ, which yields

ˆ∂E |∇ τ E ψ| 2 -H 2 E + ∂η ∂ν E ψ 2 dH 1 ≤ 0 .
We can therefore conclude that j E (0) is not coercive.

We stress that Proposition 2.30 does not show that the condition in Proposition 2.29 is necessary, since it does not cover the case where H 2 E + ∂η ∂ν E is greater than a "small" positive constant on a "small" portion of the boundary.

C

In this subsection, we describe two ways to construct candidates for the minimal norm dual certi cate de ned in De nition 1.7. The rst one relies on an attempt to adapt the notion of calibrability, which is exploited for the denoising case in [START_REF] Chambolle | Geometric properties of solutions to the total variation denoising problem[END_REF], to general measurement operators. The second is based on vanishing derivatives pre-certi cates, introduced in [Duval and Peyré, 2015, Section 4].

As computing the minimal norm certi cate in full generality is challenging, we provide a further study of a particular setting, where Φ is assumed to be a convolution operator with radial kernel h, and the unknown image u 0 is the indicator of a ball of radius R 0 > 0. Our aim is to investigate, in this simpli ed scenario, when the non-degenerate source condition is satis ed.

Φ-calibrability

In [Chambolle et al., 2016, Section 4.3], the minimal norm certi cate is computed when Φ is the identity (and (P λ (y)) is hence the Rudin Osher Fatemi denoising problem) for a special class of sets called calibrable sets (see e.g. [START_REF] Bellettini | The Total Variation Flow in RN[END_REF], Alter et al., 2005]). It is hence tempting to adapt this strategy to our setting. An extension of the notion of calibrability for a set E with respect to Φ is to require that 1 E is a singular vector of the total variation as de ned in [Benning and Burger, 2013, De nition 4]. This leads to the following de nition.

De nition 2.31

A set E ⊂ R 2 is said to be Φ-calibrable if λ E Φ * Φ1 E ∈ ∂TV(1 E ) with λ E def. = P (E) Φ1 E 2 .
Remark 2.32 If λΦ * Φ1 E ∈ ∂TV(1 E ) for some λ ∈ R, taking the inner product with 1 E shows λ = λ E .

In [Chambolle et al., 2016, Proposition 6], it is proven that, if Φ = Id and E is calibrable, then λ E 1 E is the minimal norm certi cate. With De nition 2.31, the following analog of this result is straightforward to obtain.

Proposition 2.33 If E is Φ-calibrable, then the minimal norm dual certi cate is λ E Φ * Φ1 E . Proof : Let p ∈ H such that Φ * p ∈ ∂TV(1 E )
. Then, we have:

P (E) = Φ * p, 1 E = p, Φ1 E ≤ p H Φ1 E H . Since λ E Φ1 E = P (E)/ Φ1 E H , we obtain that p H ≤ λ E Φ1 E .
Now, from Proposition 1.3, we can readily obtain the following useful characterization of Φcalibrable sets. Proposition 2.34

A set E is Φ-calibrable if and only if E maximizes F → K(E, F ) def. = Φ 1 E P (E) , Φ 1 F P (F )
among sets of nite perimeter with positive nite measure.

When Φ = Id, a complete characterization of calibrable sets is available [START_REF] Alter | A characterization of convex calibrable sets in Rn[END_REF]. We are not aware of any result outside this speci c setting. In fact, we provide numerical evidence below suggesting that even disks are not Φ-calibrable for natural choices of Φ (namely the convolution with the Gaussian and ideal low-pass kernels).

Deconvolution of a disk. Let us de ne

α(R 0 ) def. = sup K(R, R 0 ) K(R 0 , R 0 ) , R > 0 with K(R, R ) = K 1 B(0,R) , 1 B(0,R ) .
With this notation we have that, if 1 B(0,R 0 ) is Φ-calibrable, then α(R 0 ) ≤ 1. Figure 9 shows the graph of R 0 → α(R 0 ) when Φ is the convolution with the Gaussian or ideal low-pass kernel, for several values of the variance and cut-o frequency. In every situation, these plots seem to indicate that α(R 0 ) > 1 for every R 0 > 0, and hence that 1 B(0,R 0 ) is never Φ-calibrable. This questions the relevance of our de nition of Φ-calibrability. In fact, it could even be that there are no Φ-calibrable sets, even for reasonable choices of measurement operator Φ.

Pre-certi cates

In [Duval and Peyré, 2015, Section 4], the notion of dual pre-certi cate is introduced. Its authors de ne it as any "good candidate" p ∈ H for solving Φ * p ∈ ∂TV(u), where u is an admissible function for (P 0 (y 0 )) whose optimality is to be proven. In this subsection, we de ne an analog of the so-called vanishing derivative pre-certi cate, in order to study the optimality of 1 E with E a simple set. We assume as in Theorem 2.27 that Φ * is continuous from H to C 1 b (R 2 ). We also assume that E is of class C 3 , as any set E satisfying Φ * p ∈ ∂TV(1 E ) for some p ∈ H has this regularity.

Necessary optimality conditions. We know that Φ * p is a dual certi cate associated to 1 E if and only if Φ * p ∈ ∂TV(0) and E minimizes J(E)

def.

= P (E) -´E Φ * p. Using the expression of the rst order shape derivative of J, we obtain the following necessary optimality conditions:

ˆE Φ * p = P (E) , (24)
Φ * p ∂E = H E . (25)
Now, provided E is simple, (25) in particular implies that

ˆ∂E Φ * p = ˆ∂E H E = 2π . ( 26
)
One could then de ne the vanishing derivative pre-certi cate as the solution of ( 24) and ( 26) with minimal norm.

De nition 2.35

We call vanishing derivative pre-certi cate associated to a simple set E the function η v = Φ * p v with p v the unique solution of

min p∈H p H s.t. p, ˆE ϕ = P (E) and p, ˆ∂E ϕ = 2π . (27)
If the source condition holds (i.e. there exists η ∈ ImΦ * ∩ ∂TV(0) such that η ∈ ∂TV(1 E ), then p v is well-de ned as Problem ( 27) is feasible. Since any dual certi cate satis es ( 24) and ( 26), we have the following result.

Proposition 2.36 If η v ∈ ∂TV(0), then η v is the minimal norm dual certi cate.

The constraints in ( 27) can be rewritten as follows:

p, f E = 1 , p, g E = 1 , with        f E def. = 1 P (E) ˆE ϕ , g E def. = 1 2π ˆ∂E ϕ . (28)
We hence have that p v satis es 27) for E = 1 B(0,1) and Φ the convolution with the Gaussian kernel with variance σ = 0.2.

p v = αf E + βg E with α def. = g E 2 -f E , g E f E 2 g E 2 -f E , g E 2 and β def. = f E 2 -f E , g E f E 2 g E 2 -f E , g E 2 . x -2 -1 0 1 2 y -2 -1 0

Remark 2.37

In some sense, the function λ E Φ * Φ1 E considered above can be seen as the pre-certi cate associated to the single constraint (24). Imposing the second constraint (26) allows to build a pre-certi cate which is a combination of the two functions f E and g E (instead of f E only for λ E Φ * Φ1 E ). Other constraints could also be considered, as every certi cate satis es (25), which is only exploited "in average" in (26). Considering [Duval and Peyré, 2015, Propositions 7 and 8], it would be interesting to investigate whether the non-degenerate source condition holds if and only if η v = η 0 and η 0 is non-degenerate, and also whether η v = η 0 is a necessary condition for support recovery.

Deconvolution of the unit disk. We now focus on the case where E = 1 B(0,1) and Φ is the convolution with the Gaussian kernel with variance σ. We provide in Figure 10 a plot of Φ * f E and Φ * g E , which are the two "basis functions" from which η v is built. From (5), we know that a way to show

η v ∈ ∂TV(0) is to nd a vector eld z ∈ L ∞ (R 2 , R 2 ) such that div z = η v .
Since f E and g E are radial, so is η v . It is hence natural to look for a radial vector eld z (i.e. such that there exists z r : R + → R with z(x) = z r ( x ) x/ x for almost every x ∈ R 2 ). In this case we have div z = η v if and only if, for every r > 0:

ηv (r) = 1 r ∂ ∂r (r z r )(r) ⇐⇒ r ηv (r) = ∂ ∂r (r z r )(r) ⇐⇒ z r (r) = 1 r ˆr 0 ηv (s) s ds .
This shows one only needs to ensure the mapping f v de ned by satis es f v ∞ ≤ 1 to show η v ∈ ∂TV(0). Figure 11 shows the graph of f v for several values of σ. This suggests that there exists σ 0 > 0 such that η v is a dual certi cate (and hence the one with minimal norm) for every σ < σ 0 . It even seems that σ 0 ≥ 0.75.

f v : R + → R r → 1 r ˆr 0 ηv (s) s ds (29)

Remark 2.38

One could wonder whether looking for a radial vector eld is not restrictive. In fact, if a vector eld z is suitable, then so is the radial vector eld z de ned by z(x)

def.

= zr ( x ) x x with zr (r)

def. = 1 2π ˆ2π 0 z r (r, θ) dθ ,
where z r denotes the radial component of z. Indeed, we have |z(r)| ≤ 1 for all r with equality if and only if z r (r, θ) = 1 for almost every θ or z r (r, θ) = -1 for almost every θ. Moreover

η v (r) = 1 2π ˆ2π 0 η v (r) dr = 1 r ∂ ∂r r 1 2π ˆ2π 0 z r (r, θ) dθ + 1 r 1 2π ˆ2π 0 ∂z θ ∂θ (r, θ) dθ = 1 r ∂ ∂r (r zr ) = div z.
Finally, we can investigate if the non-degenerate source condition holds. As explained in Section 3.3, it is su cient to show that

sup x∈∂E H 2 E (x) + ∂η 0 ∂ν E (x) < 0 .
In our case H E is constant equal to one, and, since η 0 is radial, ∂η 0 ∂ν E is constant on ∂E. Proving that ∂η v ∂r (1) < -1 is hence su cient. In Figure 12, we numerically compute this quantity and notice it is always the case, even when η v / ∈ ∂TV(0). This suggests that there exists σ 0 > 0 such that, for every σ > 0, the non-degenerate source condition holds (and, from our experiments, it seems that σ 0 ≥ 0.75). Surprisingly, σ → ∂ηv ∂r (1) does not seem to be monotonous, even on [0, σ 0 ). Beyond the radial case. A way to ensure η v ∈ ∂TV(0) in the general case is to solve the (generalized) Cheeger problem associated to η v . This can be done using the numerical method presented in Section 2. Now, to ensure the non-degenerate source condition holds, one must also show that E is the unique (non-trivial) solution of the prescribed curvature problem associated to η v . In [START_REF] Buttazzo | On the selection of maximal Cheeger sets[END_REF], Carlier et al., 2009], a method is proposed to compute a "maximal" solution of the Cheeger problem, which is closely related to the prescribed curvature problem. It relies on the introduction of a small strictly convex penalization. We argue that these ideas could be adapted to numerically compute a "maximal" element of the face exposed by η v (i.e. a function whose chains are maximal in this face). If this element is proportional to 1 E , then 1 E is identi able and η v is the minimal norm certi cate.

C

Summary. In this part, we collected properties of the total variation unit ball and its (exposed) faces. The main interest of this analysis is a useful result, stating that elements of a given face can be decomposed using a common set of atoms. We also de ned the class of M -simple functions, which are the sparse objects naturally associated to the total variation. Then, we investigated the behaviour, under variations of the associated curvature functional, of solutions of the precribed curvature problem. This allowed us to prove that the dimension of the faces of the total variation unit ball is in some sense stable. Using this last result, we nally managed to prove that, under a so-called non-degenerate source condition, the jump set of an unknown M -simple function can be exactly recovered in a low noise regime.

Towards grid-free numerical methods. Our support recovery result is yet another motivation for designing grid-free numerical methods allowing to solve (P λ (y)). Indeed, it shows that, under suitable assumptions, if y = Φu 0 + w with u 0 an M -simple function, then solutions of (P λ (y)) are themselves M -simple. Being able to compute approximate solutions with this structure is hence highly desirable. As explained in Section 2.3, traditional grid-based method are not designed with that goal in mind. They do not produce a continuous domain representation of the reconstructed images, which are often blurry. The aim of Part 3 is to propose a numerical solver which does not su er from these drawbacks, and accurately estimates the jump set of solutions.

P 3 G

In this part, we focus on the numerical resolution of (P λ (y)). We present a so-called grid-free algorithm, which does not rely on the introduction of a xed spatial discretization. It iteratively constructs an approximate solution built from the atoms promoted by the total variation, namely indicator functions of simple sets. We choose to numerically represent these sets using simple polygons. We provide a complete implementation of the proposed method in the following online repositories.

• https://github.com/rpetit/pycheeger

• https://github.com/rpetit/tvsfw Its performance is investigated in the last section, and comparisons with grid-based approaches are provided. This part is based on the following publication.

Conclusion

F W

In [START_REF] Bredies | Inverse problems in spaces of measures[END_REF], Boyd et al., 2017, Denoyelle et al., 2019], variants of the conditional gradient algorithm, also known as Frank-Wolfe algorithm, were introduced to perform continuous domain sparse spikes recovery. In this section, we adapt this fruitful approach to our setting, and derive a modi ed Frank-Wolfe algorithm allowing to solve (P λ (y)) in a grid-free manner.

In all the following, as we are interested in the numerical resolution of (P λ (y)), we only consider the case of nitely many measurements, i.e. H = R m from some m ∈ N * .

F W

Frank-Wolfe algorithm (see Algorithm 1) allows to minimize a function f over a subset C of a Banach space. The objective f is assumed to be convex and di erentiable, and the admissible set C convex and weakly compact. Each step of the algorithm consists in minimizing the rst order expansion of f on C, and building the next iterate as a convex combination of the obtained point and the current iterate.

Algorithm 1: Frank-Wolfe algorithm k+1] ) // nal update 9 end 10 end Sparse greedy updates. Our main interest for Frank-Wolfe algorithm lies in the following classical observation: since the linear minimization step (Line 2) consists in minimizing a linear form over the weakly compact convex set C, at least one of its solutions is an extreme point of C. Choosing such a point at every iteration, we obtain that x [k] is a convex combination of at most k extreme points of C, provided x [0] = 0. This is particularly interesting in situations where the extreme points of C are "atoms" and where one is looking for a solution which is a sparse combination of those.

Data: objective f , domain C, starting point x [0] ∈ C Result: point x * 1 while true do 2 nd s [k] ∈ Argmin s∈C f (x [k] ) + df (x [k] ).(s -x [k] ) 3 if df (x [k] ).(s [k] -x [k] ) = 0 then 4 output x * ← x [k] , which is optimal 5 else 6 γ [k] ← Argmin γ∈[0,1] f (x [k] + γ(s [k] -x [k] )) // line search 7 x[k+1] ← x [k] + γ [k] (s [k] -x [k] ) // tentative update 8 choose any x [k+1] such that f (x [k+1] ) ≤ f (x [
Choice of the nal update. An important feature of the algorithm is that, while the classical update (Line 8) is to take x [k+1] to be equal to x[k+1] , all convergence guarantees are preserved if one chooses instead any x [k+1] ∈ C such that f (x [k+1] ) ≤ f (x [k+1] ). As suggested by [START_REF] Bredies | Inverse problems in spaces of measures[END_REF], Boyd et al., 2017, Denoyelle et al., 2019] in the context of sparse spikes recovery, on can extensively make use of this exibility to produce iterates that are as sparse as possible.

P

Frank-Wolfe algorithm can not be directly applied to (P λ (y)), as its objective

F (u) def. = TV(u) + 1 2λ Φu -y 2
is not di erentiable and its admissible set is unbounded. In this subsection, we derive an equivalent formulation of (P λ (y)) which ts into this framework, and then describe the resulting algorithm.

Epigraphical lift. To obtain a problem equivalent to (P λ (y)) with a di erentiable objective and a weakly compact admissible set, we perform an epigraphical lift. This strategy is used in [START_REF] Denoyelle | The Sliding Frank-Wolfe Algorithm and its Application to Super-Resolution Microscopy[END_REF] (following an idea of [START_REF] Harchaoui | Conditional gradient algorithms for norm-regularized smooth convex optimization[END_REF], Section 2, Paragraph "penalized norm minimization"]).

Proposition 3.1

Problem (P λ (y)) is equivalent to inf (u,t)∈L 2 (R 2 )×R G(u, t) def. = 1 2λ Φu -y 2 + t s.t. TV(u) ≤ t ≤ 1 2λ y 2 , (Q λ (y))
i.e. both problems have the same value and:

1. if (u, t) is a solution of (Q λ (y)), then t = TV(u) and (u, v) is a solution of (P λ (y)),

2. if u is a solution of (P λ (y)), then (u, TV(u)) is a solution of (Q λ (y)).

Proof : This is a straightforward consequence of the following fact: if u is a solution of (P λ (y)), then we have that F (u) ≤ F (0), which yields TV(u) ≤ 1 2λ y 2 . The admissible set of (P λ (y)) can hence be restricted to functions u satisfying TV(u) ≤ 1 2λ y 2 without modifying its solutions.

The admissible set of (Q λ (y)) is now weakly compact and convex. Moreover, its objective G is di erentiable with, for all (u, t) ∈ L 2 (R 2 ) × R:

dG(u, t) : L 2 (R 2 ) × R → R (v, s) → ˆR2 1 λ Φ * (Φu -y) v + s .
Linear minimization step. Applying Algorithm 1 to (Q λ (y)), the linear minimization step at iteration k (Line 2) writes inf

(u,t)∈L 2 (R 2 )×R ˆR2 1 λ Φ * (Φu [k] -y) u + t s.t. TV(u) ≤ t ≤ 1 2λ y 2 . ( 30
)
As explained above, there always exists a solution of (30) which is an extreme point of the admissible set. These points are described by the following lemma, whose proof is given in Appendix E.

Lemma 3.2

If t * > 0, the extreme points of C def. = {u ∈ L 2 (R 2 ) | TV(u) ≤ t ≤ t * } are:
• (0, 0) ,

• ( t * 1 E /P (E), t * ) with ∈ {-1, 1}, E ⊂ R 2 simple and 0 < |E| < +∞ .
With this result, we obtain that nding a solution of (30) among the extreme points of the admissible set is equivalent to solving sup

E⊂R 2 ∈{-1,1} P (E) ˆE η [k] s.t. E simple, 0 < |E| < +∞ , (31)
with η [k] def. = (-1/λ)Φ * (Φu [k] y). Indeed, a simple computation shows that, if the value of ( 31) is smaller than 1, then (0, 0) is optimal. Otherwise, de ning M def.

= y 2 /(2λ), for any solution (E, ) of ( 31), we get that ( M 1 E /P (E), M ) is optimal. Properties of (31), which is reminiscent of the Cheeger problem (see the surveys [Parini, 2011, Leonardi, 2015]), are discussed in Section 2.1.

Form of the iterates. Applying Algorithm 1 to (Q λ (y)) with (30) solved as described above, one obtains a sequence of iterates (u [k] , t [k] ) k≥0 whose form is described by the following lemma.

Lemma 3.3

If (u [0] , t [0] ) = (0, 0), for every k ∈ N, we either have that (u [k] , t [k] ) = (0, 0) or that there exists

N [k] ∈ N * , [k] ∈ {-1, 1} N [k] , a [k] ∈ (R * + ) N [k] and a collection E [k] = E [k] 1 , ..., E [k] N [k]
of simple sets of positive nite measure such that

(u [k] , t [k] ) =   N [k] i=1 [k] i a [k] i 1 E [k] i , N [k] i=1 a [k] i P (E [k] i )   . ( 32
)
Proof : We argue by induction. The result is already known for k = 0. Let us x k ∈ N and assume that (u [k] , t [k] ) has the right form.

1. If (u [k] , t [k] ) = (0, 0) then:

• if (0, 0) solves ( 30) then (u [k+1] , t [k+1] ) = (0, 0),

• otherwise there exists ∈ {-1, 1} and a simple set E with 0 < |E| < +∞ such that the couple ( M 1 E /P (E), M ) solves (30). We therefore obtain that (u [k] , t [k] ) is as in (32) with andE [k+1] = (E) .

N [k+1] = 1, [k+1] = ( ), a [k+1] = (γ [k] M/P (E))
2. If (u [k] , t [k] ) is as in (32) then:

• if (0, 0) solves ( 30), we get that (u [k+1] , t [k+1] ) is as in (32) with

N [k+1] = N [k] , [k+1] = [k] , a [k+1] = (1 -γ [k]
) a [k] and

E [k+1] = E [k] .
• Otherwise as above denoting ( M 1 E /P (E), M ) a solution of ( 30) we obtain the result with

N [k+1] = N [k] + 1, [k+1] = ( [k] , ), a [k+1] = ((1 -γ [k] ) a [k] , γ [k] M/P (E))
and

E [k+1] = (E [k] , E).
"Fully-corrective" variant. Lemma 3.3 allows us to use a so-called fully corrective variant of Frank-Wolfe, meaning that instead of obtaining the next iterate as a convex combination of the new atom and the previous iterate (as in Line 7 of Algorithm 1), we nd a [k+1] by minimizing G N ( , a, E) de ned by k+1] , E [k+1] ) xed. This minimization amounts to solving a LASSO problem with a positivity constraint and a weighted 1 penalty (the weights being the perimeters of the sets in E [k+1] ). Indeed, given N ∈ N * , a collection of sets E 1 , ..., E N , and two vectors a ∈ (R + ) N , ∈ {-1, 1} N we have:

G N ( , a, E) def. = G N i=1 i a i 1 E i , N i=1 a i P (E i ) with N = N [k] + 1 and ( , E) = ( [
G N ( , a, E) = 1 2λ Φ E a -y 2 + N i=1 P (E i ) a i with Φ E def. =   i ˆEi ϕ j 1≤i≤N 1≤j≤m   T ∈ R m×N .
The main interest of Lemma 3.3 is that, considering (32), this choice of a [k+1] yields a value of the objective no higher than that obtained with the standard update, and hence does not break convergence guarantees.

Modi ed Frank-Wolfe algorithm applied to (Q λ (y)). Dropping the dependence on the auxiliary variable t ∈ R we obtain Algorithm 2. As a result of the equivalence between (P λ (y)) and (Q λ (y)), this algorithm is a valid application of Algorithm 1 to (P λ (y)), in the sense that Proposition 3.5 below holds. Before studying convergence results, let us make a comment on the stopping condition.

Algorithm 2: modi ed Frank-Wolfe algorithm applied to (P λ (y)) Data: measurement operator Φ, observations y, regularization parameter λ Result:

function u * 1 u [0] ← 0 2 N [0] ← 0 3 while true do 4 η [k] ← -1 λ Φ * (Φu [k] -y) 5 (E * , * ) ← Argmax E⊂R 2 ∈{-1,1} ´E η [k] /P (E) s.t. E simple, 0 < |E| < +∞ 6 if ´E * η [k] ≤ P (E * ) then 7 output u * ← u [k]
, which is optimal

8 else 9 E [k+1] ← (E [k] 1 , ..., E [k] N [k] , E * ) 10 [k+1] ← ( [k] 1 , ..., [k] N [k] , * ) 11 a [k+1] ← Argmin a∈(R + ) N [k] +1 1 2λ Φ [k+1] E [k+1] a -y 2 + N [k] +1 i=1 P (E [k+1] i ) a i 12
remove atoms with zero amplitude

13

N [k+1] ← number of atoms in E [k+1] 14 u [k+1] ← N [k+1] i=1 a [k+1] i [k+1] i 1 E [k+1] i 15 end 16 end
Stopping condition. The stopping condition (Line 3 of Algorithm 1) is replaced in Algorithm 2 by

sup P (E) ˆE η [k] , E ⊂ R 2 simple , 0 < |E| < +∞, ∈ {-1, 1} ≤ 1 , (33)
which is equivalent to η [k] ∈ ∂TV(0). Since the optimality of a [k] at Line 11 always ensures

ˆR2 η [k] u [k] = N [k] i=1 a [k] i P (E [k] i ) ≥ TV(u [k] ) ,
we obtain that u [k] solves (P λ (y)) as soon as (33) holds.

C 1.3.1. Convergence in objective value

Curvature constant. As pointed out in [Jaggi, 2013], in the convergence analysis of Frank-Wolfe algorithm applied to a function f , a measure of the "nonlinearity" of f , called the curvature constant, naturally appears. Let f be a convex di erentiable function and C a weakly compact convex set. The curvature constant C f of f with respect to C is de ned by:

C f def. = sup x,s∈C, γ∈[0,1], y=x+γ(s-x) 2 γ 2 [f (y) -f (x) -df (x).(y -x)] .

Proposition 3.4

The curvature constant C G of G with respect to the admissible set C of (Q λ (y)) satis es:

C G ≤ 1 λ c 2 Φ y 2 λ 2 ,
where c 2 = 1/ √ 4π is the isoperimetric constant.

Proof : For any u, v ∈ L 2 (R 2 ) and t, s ∈ R, we have:

G(u, t) -G(v, s) -dG(u, t).(v, s) = 1 2λ Φ(u -v) 2 .
We hence obtain that

C G = sup u,v∈L 2 (R 2 ) 1 λ Φ(u -v) 2 s.t. TV(u) ≤ 1 2λ y 2 , TV(v) ≤ 1 2λ y 2 . ( 34
)
Now, if u, v are admissible for (34), using the isoperimetric inequality (3), we nally get:

Φ(u -v) ≤ Φ u -v L 2 (R 2 ) ≤ Φ c 2 TV(u -v) ≤ c 2 Φ y 2 λ .
Convergence in objective value. As already mentioned, Algorithm 2 is a valid application of Algorithm 1 to (P λ (y)), in the sense that the following property holds (see [Jaggi, 2013]):

Proposition 3.5 Let (u [k] ) k≥0 be a sequence produced by Algorithm 2. Then for any solution u * of (P λ (y)), we have:

∀k ∈ N, F (u [k] ) -F (u * ) ≤ 2 C G k + 2 .
Approximate linear minimization step. As discussed in [Jaggi, 2013], the linear minimization step, which consists in solving (30) or equivalently (31), can be carried out approximately.

In fact if there exists δ > 0 such that, for every k, the couple (E * , * ) computed at Line 5 of Algorithm 2 is an [k] -minimizer of ( 31) with [k] = δ C G /(k + 2), then

∀k ∈ N * , F (u [k] ) -F (u * ) ≤ 2 C G k + 2 (1 + δ) .

Convergence of the iterates

Convergence of minimizing sequences. Now, let us state a general property of minimizing sequences of (P λ (y)), which hence applies to the sequence of iterates produced by Algorithm 2.

Proposition 3.6 Let (u [k] ) k≥0 be a minimizing sequence for (P λ (y)). Then there exists a subsequence which converges strongly in L 1 loc (R 2 ) and weakly in L 2 (R 2 ) to a solution u * of (P λ (y)). Moreover, we have TV(u [k] ) → TV(u * ).

Proof : We rst notice that since (u [k] ) k≥0 is a minimizing sequence then (TV(u [k] )) k≥0 is bounded, and (from the isoperimetric inequality) (u [k] ) k≥0 is therefore bounded in L 2 (R 2 ). There hence exists a subsequence (not relabeled) which converges weakly in L 2 (R 2 ) to u * ∈ L 2 (R 2 ). Since this subsequence satis es

sup k∈N u [k] L 2 (R 2 ) + TV(u [k] ) < +∞ ,
from [Ambrosio et al., 2000, Theorem 3.23] we obtain that (up to the extraction of a further subsequence, still not relabeled) (u [k] ) k≥0 converges strongly in L 1 loc (R 2 ) (to a limit which is necessarily u * ). The fact u * solves (P λ (y)) then follows from the lower semicontinuity of the data delity term with respect to the weak L 2 (R 2 ) topology, and the lower semicontinuity of the total variation with respect to the strong L 1 loc (R 2 ) topology. Finally, since

1 2λ Φu * -y 2 ≤ lim inf k→+∞ 1 2λ Φu [k] -y 2 , TV(u * ) ≤ lim inf k→+∞ TV(u [k] ) , 1 2λ Φu * -y 2 + TV(u * ) = lim k→+∞ 1 2λ Φu [k] -y 2 + TV(u [k] ) ,
we obtain TV(u [k] ) → TV(u * ).

As explained earlier, the fact TV(u [k] ) → TV(u * ) and the strong L 1 loc convergence of (u [k] ) k≥0 towards u * also imply Du

[k] * Du * .
Properties of sets appearing in Algorithm 2. An important observation is that, for every k ∈ N, the set E * introduced at Line 5 of Algorithm 2 is a minimizer of the prescribed curvature problem associated to λ [k] η [k] with 1/λ [k] the value of (31) 1 . Moreover, if (u [k] ) is a subsequence as in Proposition 3.6, we have that Φu [k] → Φu * and hence η [k] → η * strongly in L 2 (R 2 ) with η * = (-1/λ)Φ * (Φuy), which also yields λ [k] → 1 2 . This can be exploited by using Lemma 1.11 to obtain "uniform" properties of sets appearing in Algorithm 2, which ultimately yields the stronger convergence results presented below. 1 The value of ( 31) is nonzero unless η [k] = 0, which occurs if and only if y = 0 (in which case u * = 0 is the unique solution of (P λ (y))). To see this, one simply needs to consider a Lebesgue point of η [k] at which the function is nonzero, and to show the objective is strictly positive for su ciently small balls. 2 To see this, it is su cient to note that

´E η 1 P (E) - ´E η 2 P (E) ≤ c2 η1 -η2 L 2 (R 2 )
for every E with 0 < |E| < +∞ and P (E) < +∞.

Additional properties of sequences produced by Algorithm 2. If (u [k] ) k≥0 is a subsequence as in Proposition 3.6 then (λ [k] η [k] ) k≥0 converges strongly in L 2 (R 2 ), and its terms all belong to ∂TV(0). Lemma 1.11 hence yields the existence of R > 0 such that Supp(u [k] ) ⊂ B(0, R) for every k ∈ N. We consequently obtain the following strict BV convergence result.

Proposition 3.7

Let (u [k] ) k≥0 be a sequence produced by Algorithm 2. Then, up to the extraction of a subsequence, we have that (u [k] ) k≥0 converges strictly in BV(R 2 ) to a solution u * of (P λ (y)).

This in turn yields a convergence result for the level sets of the iterates, which is given below (we remind the reader that the inner limit of a sequence of sets is de ned in De nition 1.9).

Corollary 3.8 Let (u [k] ) n≥0 be a subsequence as in Proposition 3.7. Then for almost every t ∈ R, we have

lim n→+∞ |U (t) k U (t) * | = 0 and ∂U (t) * ⊆ lim inf n→+∞ ∂U (t)
k .

Proof : The strong convergence of (u [k] ) k≥0 towards u * in L 1 (R 2 ) and Proposition D.1 yield:

lim n→+∞ |U (t) k U (t) * | = 0 for almost every t ∈ R .
We now x such t ∈ R and let x ∈ ∂U (t) * . We want to show that x ∈ lim inf

n→+∞ ∂U (t) k , which is equivalent to lim sup k→+∞ dist x, ∂U (t) k = 0 .
By contradiction, if the last identity does not hold, we have the existence of r > 0 and of ϕ such that ∀k ∈ N, B(x, r) ∩ ∂U (t) ϕ(k) = ∅ . Hence for all k, we either have

B(x, r) ⊂ U (t) ϕ(k) or B(x, r) ⊂ R 2 \ U (t) ϕ(k) . If B(x, r) ⊂ U (t)
ϕ(k) for a given k, using that η * ∈ ∂TV(u * ) and hence that U (t) * satis es point 4 of Lemma 1.11 we obtain

U (t) ϕ(k) U (t) * ≥ U (t) ϕ(k) \ U (t) * ≥ B(x, r) \ U (t) * ≥ C |B(x, r)| .
We can in the same way show k are obtained as a nite number of unions and intersections of the sets in E [k] and their complements. Although the latter sets have the desired property, it is not stable by such operations, as Figure 13 shows.

U (t) ϕ(k) U (t) * ≥ C |B(x, r)| if B(x, r) ⊂ R 2 \ U (t) ϕ(k) ,

S 1.4.1. Presentation

Several works [START_REF] Bredies | Inverse problems in spaces of measures[END_REF], Boyd et al., 2017, Denoyelle et al., 2019] have advocated for the use of a special nal update, which helps identify the sparse structure of the sought-after signal. Loosely speaking, it would amount in our case to running, at the end of iteration k -1, the gradient ow of the mapping

(a, E) → G( [k] , a, E) (35)
initialized with (a [k] , E [k] ), so as to nd a set of parameters at which the objective is smaller. Formally, this would correspond 1 to nding a curve t → (a i (t), E i (t)) i=1,...,N [k] such that for all t:

∀i ∈ {1, ..., N [k] },        a i (t) = -P (E i (t)) - [k] i ˆEi (t) η(t) , V i (t) = -a i (t) H E i (t) - [k] i η(t) , (36)
1 The formulas in (36) can be formally obtained by using the notion of shape derivative, see [Henrot and Pierre, 2018, Chapter 5].

where V i (t) denotes the normal velocity of the boundary of E i (t) and

η(t) = - 1 λ Φ * (Φu(t) -y) , u(t) = N [k] i=1 a i (t) [k] i 1 E i (t) .
The study of this gradient ow (existence, uniqueness) seems to be challenging, and we did not investigate these questions.

For our purpose, it is enough to introduce the sliding step as a local descent on the mapping de ned in (35) initialized with (a [k] , E [k] ). We however need to ensure this step does not increase the value of the objective, and we hence ask to nd (a, E) = (a i , E i ) 1≤i≤N [k] such that E i is simple for all i and

G N [k] [k] , a, E ≤ G N [k] [k] , a [k] , E [k] . ( 37
)
The resulting algorithm is Algorithm 3. The introduction of Line 15, which is discussed in the next paragraph, ensures that all convergence guarantees derived for Algorithm 2 remain valid.

Algorithm 3: modi ed Frank-Wolfe algorithm applied to (P λ (y)) (with sliding) Data: measurement operator Φ, observations y, regularization parameter λ Result: k] , which is optimal

function u * 1 u [0] ← 0 2 N [0] ← 0 3 while true do 4 η [k] ← -1 λ Φ * Φu [k] -y 5 (E * , * ) ← Argmax E⊂R 2 ∈{-1,1} ´E η [k] /P (E) s.t. E simple, 0 < |E| < +∞ 6 if ´E * η [k] ≤ P (E * ) then 7 output u * ← u [
8 else 9 E [k+1] ← (E [k] 1 , ..., E [k] N [k] , E * ) 10 [k+1] ← ( [k] 1 , ..., [k] N [k] , * ) 11 a [k+1] ← Argmin a∈(R + ) N [k] +1 1 2λ Φ [k+1] E [k+1] a -y 2 H + N [k] +1 i=1 P (E [k+1] i ) a i 12
remove atoms with zero amplitude

13

N [k+1] ← number of atoms in E [k+1] 14 perform a local descent on (a, E) → G N [k+1] ( [k+1] , a, E) initialized with (a [k+1] , E [k+1] )
15 repeat the operations of Lines 11-13

16 u [k+1] ← N [k+1] i=1 a [k+1] i 1 E [k+1] i 17 end 1.4.2

. In uence on convergence

The version of the sliding step we use was introduced in [START_REF] Denoyelle | The Sliding Frank-Wolfe Algorithm and its Application to Super-Resolution Microscopy[END_REF]. In practice, it allows to considerably improve the convergence speed of the algorithm. It also produces sparser solutions: if the solution is a combination of a few indicator functions, removing the sliding step typically produces iterates made of a much larger number of atoms. In fact, it seems that the majority of these atoms correct the crude approximations of the support of the solution made during the rst iterations.

In [START_REF] Denoyelle | The Sliding Frank-Wolfe Algorithm and its Application to Super-Resolution Microscopy[END_REF], the introduction of the sliding step allowed the authors to derive improved convergence guarantees (i.e. nite time convergence) in the context of sparse spikes recovery. Their proof relies on the fact that, at the end of every sliding step, a "critical point" of the objective is reached. However, the existence issues mentioned in Section 1.4.1 make the adaptation of these results to our setting di cult.

A reasonable de nition of "critical point" in our context could be as follows.

De nition 3.9 (Critical point of

G N ) Let N ∈ N * , ∈ {-1, 1} N , a ∈ (R * + )
N and E 1 , ..., E N be simple subsets of R 2 with positive and nite measure. We say that ( i , a i , E i ) i=1,...,N is a critical point of the mapping

G N if ∀i ∈ {1, ..., N },    P (E i ) = i ˆEi η , H E i = i η , (38)
where

η def. = - 1 λ Φ * (Φu -y) and u def. = N i=1 a i i 1 E i .
With this de nition, reaching a critical point at the end of the sliding step implies all sets in E [k] have the same distributional curvature (up to a sign). On the contrary, without sliding, sets are never modi ed and their curvature depends on the iteration during which they were introduced.

We stress that if, for a given iteration, a critical point is reached at the end of the sliding step, then Line 15 can be skipped, since the rst equality in (38) ensure a [k+1] minimizes

a → G N [k+1] ( [k+1] , a, E [k+1] ) .
If a critical point is reached, the fact the sets in E [k] have the same distributional curvature can be exploited to obtain "uniform" density estimates for the level sets of u [k] , in the spirit of [Maggi, 2012, Corollary 17.18]. It is then natural to wonder whether this could be used to prove lim sup n→+∞

∂U (t) n ⊆ ∂U (t) * , (39)
and hence the convergence of (∂U

(t) k ) k≥0 towards ∂U (t) * .
A major obstacle towards this result is that, although Lemma 1.11 provides a uniform upper bound on the perimeter of sets appearing in the algorithm, to our knowledge, it does not seem possible to derive such a bound for the perimeter of the level sets, which prevents one from using the potential weak-* convergence of D1 U (t)

k towards D1 U (t) * .
Link with conic particle gradient descent. Assuming for simplicity that [k] = (1, ..., 1) and, for every i

= j, that H 1 (∂ * E i ∩ ∂ * E j ) = 0, we obtain F   N [k] i=1 a i 1 E i   = G N [k] ( [k] , a, E) ,
and hence see that the sliding step amounts to performing a local descent on the mapping

(a, E) → F   N [k] i=1 a i 1 E i   .
In [Chizat, 2021, Chizat, 2022], a seemingly related task is considered. Given Θ a parameter space, the author investigates the minimization of

J : M + (Θ) → R µ → 1 2λ Φµ -y 2 + µ(Θ)
using conic particle gradient algorithms, which consist in performing gradient-based optimization on

(a, x) → J N i=1 h(a i )δ θ i , (40)
with h : R + → R + a smooth increasing bijection. The main nding of these works is that, for h(a) = a2 and for a speci c choice of metric on R + × Θ (plus a few assumptions), the gradient ow of (40) converges to a minimizer of J. An interesting avenue for future research could be to study how this analysis translates to our setting. On the practical side, one could investigate whether performing multiplicative (rather than additive) updates on the amplitudes a improves the convergence speed of the local descent, as suggested in [Chizat, 2021].

P C

To implement Algorithm 3, one needs two oracles to carry out the operations of Lines 5 and 141 .

In this section, we focus on the rst task (the second is the subject of Section 3.1). Our aim is hence to design a numerical method to approximately solve sup

E⊂R 2 J(E) def. = 1 P (E) ˆE η s.t. E simple, 0 < |E| < +∞ , (41)
which is called a generalized Cheeger problem2 . In section 2.1, we explain its connection with the classical Cheeger problem and provide properties of its solutions, which are called generalized Cheeger sets. We also review existing numerical methods for solving it, and motivate our choice to introduce a slightly di erent procedure. Our approach consists in choosing a (possibly large) integer n ≥ 3, and looking for an optimal set among simple polygons with at most n sides. The existence of such sets is proved in Section 2.2. To approximate this optimal polygon, we optimize the vertices of an initial polygon using a rst order algorithm. The initialization is performed by solving a relaxation of (41) on a xed grid. The details of this numerical method are presented in Sections 2.3 and 2.4. Finally, in Section 2.5, we study how these optimal polygons compare to solutions of (41) in a speci c setting, where η is assumed to be radial.

A C

Related problems. We rst stress that ( 41) is reminiscent of the Cheeger problem, which, given an open bounded set Ω ⊂ R d with Lipschitz boundary, consists in nding a subset E of Ω minimizing P (E)/|E| (see the surveys [Parini, 2011, Leonardi, 2015]). Many generalizations of the Cheeger problem have been considered. Let us mention [START_REF] Ionescu | Generalized Cheeger sets related to landslides[END_REF], in which the usual perimeter and volume are replaced by weighted versions (the integral on ∂ * E and E of some nonnegative weight functions), and [START_REF] Caselles | Anisotropic Cheeger Sets and Applications[END_REF][START_REF] Brasco | The fractional Cheeger problem[END_REF], in which the case of anisotropic and fractional perimeters are respectively investigated.

Our setting. The problem we consider slightly di ers from the ones mentioned above. First, the emphasis is not on the domain Ω anymore (( 41) is solved on the whole plane), but on the weight function η. To put it another way, maximizing J on the whole plane or inside some domain Ω is equivalent as soon as η 2 concentrates1 in a compact subset of Ω. Moreover, contrary to [START_REF] Ionescu | Generalized Cheeger sets related to landslides[END_REF], the function η is not assumed to be nonnegative, and is taken in

L 2 (R 2 ) (instead of L ∞ (Ω) in this last work).
Existence of solutions. Two strategies can be considered to prove the existence of maximizers in (41). The rst is the one hinted in Section 1.2. It consists in proving the existence of maximizers for sup

u∈L 2 (R 2 ) ˆR2 ηu s.t. TV(u) ≤ 1 , (42)
and then showing at least one solution is the indicator function of a simple set (using the form of the extreme points of {TV ≤ 1} and the Krein-Milman theorem, or a result like [Carlier and Comte, 2007, Theorem 2]). The second strategy relies on purely geometric arguments, in the spirit of [Maggi, 2012, Section 12]. One can indeed prove a uniform lower bound on the perimeter of minimizing sequences, and obtain local convergence (up to extraction) towards some limit set, whose optimality can subsequently be shown.

Generic uniqueness.

As explained in [Buttazzo et al., 2007, Proposition 4.1], the solution of ( 41) is in some sense generically unique. Indeed, denoting V (η) the value of (42), we obtain that V is a continuous convex function on L 2 (R 2 ). A theorem of Mazur [Phelps, 1993, Theorem 1.20] then shows V is Gâteaux di erentiable on a dense G δ subset of L 2 (R 2 ). Since ∂V (η) is precisely the solution set of (41), generic uniqueness follows.

Existing numerical methods. In [START_REF] Carlier | Approximation of maximal Cheeger sets by projection[END_REF] a method is proposed to approximate the so-called maximal Cheeger set. It covers the case of a weighted area term and of weighted and anisotropic perimeters. The proposed approach is to solve a discretized version of (42), where the maximization is performed over the set of piecewise constant functions on some xed grid. A method which is very similar to ours is introduced in [START_REF] Ionescu | Boundary variation method for the generalized Cheeger problem[END_REF]. It consists in iteratively deforming an initial domain by computing a shape gradient at each iteration. We discovered the existence of this work during the writing of this manuscript.

E C

The aim of this section is to prove the existence of polygonal generalized Cheeger sets, i.e. maximizers of J among simple polygons with a given number of sides. In fact, we prove a slightly stronger result, namely the existence of maximizers of a relaxed energy which coincides with J on simple polygons, and the existence of a simple polygon among these maximizers.

Existence proofs in polygonal shape optimization problems. To prove the existence of optimizers for shape optimization problems in the class of polygons, a typical strategy is to use a compactness result for the Hausdor convergence of open sets. One constructs a minimizing sequence whose terms are open polygons included in a common compact set, and obtain its convergence (up to extraction) in the Hausdor sense to a "generalized polygon" (see for example [Henrot, 2006, Theorem 3.3.1] for the minimization of the rst Dirichlet eigenvalue of the Laplacian, or [Bucur and Fragalà, 2016, Proposition 9] for the minimization of the Cheeger constant). At the time we considered the problem we are discussing here, we were not aware of this literature, and followed a slightly di erent (but essentially similar) strategy. As explained below, our approach amounts to considering the objective as a function of the vertices of the input polygon.

Notations. In the following, we x an integer n ≥ 3 and denote

X n = x ∈ R n×2 [x 1 , x 2 ], ..., [x n , x 1 ] is simple .
We recall that a polygonal curve is said to be simple if non-adjacent sides do not intersect.

If x ∈ X n , then 1 ∪ n i=1 [x i , x i+1
] is a Jordan curve. It hence divides the plane in two regions, one of which is bounded. We denote this region by E x (it is hence a simple polygon). With this notation the set P n of simple polygons with at most n sides is given by

P n = {E x , x ∈ X n } . 1 If i > n we de ne xi def. = x i mod n , i.e. xn+1 = x1.
Relaxed perimeter and weighted area. We rst begin by de ning relaxed versions of the perimeter and the (weighted) area. To be able to deal with polygons with a number of vertices smaller than n, which will be useful in the following, we de ne for all m ≥ 2 and x ∈ R m×2 the following quantities:

P (x) def. = m i=1 x i+1 -x i and A(x) def. = ˆR2 η χ x ,
where χ x (y) denotes the index (or winding number) of any parametrization of the oriented polygonal curve [x 1 , x 2 ], ..., [x m , x 1 ] around y ∈ R 2 (see for instance [Rudin, 1986, Theorem 10.10]).

In particular, for every x ∈ X m (i.e. for every x ∈ R m×2 de ning a simple polygon), we have χ x = ±1 Ex and

P (x) = P (E x ) and |A(x)| = ˆEx η ,
and hence, as soon as P (x) > 0:

J(E x ) = |A(x)| P (x) .
This naturally leads us to de ne

Y m def. = x ∈ R m×2 P (x) > 0
and to de ne, abusing notation, J(

= |A(x)|/P (x) for every x ∈ Y m .

Properties of the index. As mentioned above, if E x is simple then χ x = ±1 Ex . In the general case χ x is constant on each connected component of

R 2 \ Γ x with Γ x def. = ∪ m i=1 [x i , x i+1 ].
It takes values in {-m, ..., m} and is equal to zero on the only unbounded connected component of R 2 \ Γ x . We also have ∂ supp(χ x ) ⊂ Γ x . Moreover χ x has bounded variation and, for H 1 -almost every y ∈ Γ x , there exists u + Γ (y), u - Γ (y) in {-m, ..., m} such that

Dχ x = (u + Γx -u - Γx ) ν Γx H 1 Γ x .
Positivity of the problem value. Now, we de ne α def.

= sup {J(x), x ∈ Y n }. If η = 0, then the existence of maximizers is trivial. Otherwise, there exists a Lebesgue point x 0 of η at which η is non-zero. Now the family of regular n-gons inscribed in any circle centered at x 0 has bounded eccentricity. Hence, if x n,r de nes a regular n-gon inscribed in a circle of radius r centered at x 0 , the Lebesgue di erentiation theorem 1 ensures that lim

r→0 + ´Exn,r η |E xn,r | > 0 ,
and the fact that α > 0 follows.

1 For details on the Lebesgue di erentiation theorem where balls are replaced by a family of sets with bounded eccentricity, see [Stein, 1993, Chapter 2, Paragraph 3.1].

Lemma 3.10

Let C > 0. There exists R > 0 and c > 0 such that, for every x ∈ Y n J(x) ≥ C =⇒ P (x) ≥ c and x i ≤ R for all i.

Proof : The proof is similar to that of Lemma 1.11.

Upper bound on the perimeter: the integrability of η 2 yields that for every > 0 there exists

R 1 > 0 such that ˆR2 \B(0,R1) η 2 ≤ 2 . ( 43
)
Let > 0 and R 1 > 0 such that (43) holds. We have

P (x) ≤ 1 C |A(x)| ≤ 1 C ˆR2 ∩B(0,R) η χ x + ˆR2 \B(0,R) η χ x ≤ 1 C η L 2 χ x L ∞ |B(0, R)| + χ x L 2 ≤ 1 C η L 2 n |B(0, R)| + 1 √ c 2 |Dχ x |(R 2 ) ≤ 1 C η L 2 n |B(0, R)| + 2n √ c 2 P (x) . Now, taking def. = C 4c 2 n and c = 2n C η L 2 |B(0, R)| ,
we nally get that P (x) ≤ c .

Inclusion in a ball: we take = 1 4c2n and x R 2 > 0 such that ´R2 \B(0,R2) η 2 ≤ 2 . Let us show that supp(χ x ) ∩ B(0, R 2 ) = ∅ .

By contradiction, if supp(χ x ) ∩ B(0, R 2 ) = ∅, we would have:

P (x) ≤ 1 C |A(x)| = 1 C ˆR2 \B(0,R2) η χ x ≤ ˆR2 \B(0,R2) η 2 χ x L 2 ≤ c 2 |Dχ x |(R 2 ) ≤ 2n c 2 P (x) .
Dividing by P (x) > 0 yields a contradiction. Now since

∂ supp(χ x ) ⊂ Γ x , we have diam(supp(χ x )) ≤ P (x) ≤ c which shows supp(χ x ) ⊂ B(0, R) with R def. = c + R 2 .
This in turn implies that x i ≤ R for all i.

Lower bound on the perimeter: the integrability of η 2 shows that, for every > 0, there exists δ > 0 such that

∀E ⊂ R 2 , |E| ≤ δ =⇒ ˆE η 2 ≤ 2 . Taking def. = C/(2 c 2 ), we obtain that if |supp(χ x )| ≤ δ P (x) ≤ 1 C |A(x)| = 1 C ˆsupp(χx) η ≤ 1 C ˆsupp(χx) η 2 |supp(χ x )| ≤ c 2 C P (supp(χ x )) ≤ c 2 C P (x) ,
the last inequality holding because ∂ supp(χ x ) ⊂ Γ x . We get a contradiction since P (x) is positive.

Applying Lemma 3.10 with e.g. C = α/2, and de ning

Y n def. = x ∈ R n×2 P (x) ≥ c and x i ≤ R for all i ,
we see that any maximizer of J over Y n (if it exists) is also a maximizer of J over Y n , and conversely.

Lemma 3.11 Let x ∈ R n×2 . Then, for every a ∈ R 2 , denoting (c 1 c 2 ) the 2 × 2 matrix whose columns are c 1 , c 2 ∈ R 2 , we have:

A(x) = n i=1 sign(det(x i -a x i+1 -a)) ˆax i x i+1 η = n i=1 det(x i -a x i+1 -a) ˆT1 η((x i -a x i+1 -a) y) dy ,
where ax i x i+1 denotes the triangle with vertices a, x i , x i+1 and

T 1 def. = (α, β) ∈ (R + ) 2 α + β ≤ 1
is the unit triangle.

Proof : Let us show that for all a ∈ R 2 we have

χ x = n i=1 sign(det(x i -a x i+1 -a)) 1 axixi+1 (44)
almost everywhere. We have that y ∈ R 2 is in the (open) triangle ax i x i+1 if and only if the ray issued from y directed by y -

a intersects ]x i , x i+1 [. Moreover, if y is in this triangle, then det(x i -a x i+1 -a)) is positive if and only if the triangle yx i x i+1 is oriented counterclockwise. This shows that, if R 2 \ ∪ n i=1 [x i , x i+1 ]
does not belong to any of the segments [a, x i ], evaluating the right hand side of (44) at y amounts to computing the winding number χ x (y) by applying the ray-crossing algorithm described in [START_REF] Hormann | The point in polygon problem for arbitrary polygons[END_REF]. This in particular means that (44) holds almost everywhere, and the result follows.

From Lemma 3.11, we get that A is continuous on R n×2 . This is also the case of P . Now Y n is compact and included in Y n , hence the existence of maximizers of J over Y n , which in turn implies the existence of maximizers of J over Y n . Now, one can see that in each case we have P (x) = P (y) + P (z) and χ x = χ y + χ z almost everywhere, which in turn gives that A(x) = A(y) + A(z). We hence get that P (y) = 0 or P (z) = 0, and hence J(x) = J(y) or J(x) = J(z), or that P (y) > 0 and P (z) > 0, which yields Hence J(x) is smaller than a convex combination of J(y) and J(z), which gives that it is smaller than J(y) or J(z). This shows that y or z is suitable.

We can now prove our nal result, i.e. that there exists

x * ∈ X n such that ∀x ∈ Y n , J(x * ) ≥ J(x) .
Indeed, repeatedly applying the above lemma starting with a maximizer x * of J over Y n , we either have that there exists m with 3 ≤ m ≤ n and x * ∈ X m such that J(x * ) = J(x * ), or that there exists y ∈ Y 2 such that J(x * ) ≤ J(y), which is impossible since in that case J(y) = 0 and J(x * ) = α > 0. We hence have

x * ∈ X m such that ∀x ∈ Y n , J(x * ) = J(x * ) ≥ J(x) .
We can nally build x * ∈ X n such that J(x * ) = J(x * ) by adding dummy vertices to x * , which allows to conclude.

F

In this section, we detail the initialization of our rst order method. We proceed almost exactly as in [START_REF] Carlier | Approximation of maximal Cheeger sets by projection[END_REF], i.e. we solve a discrete version of (42), in which the minimization is performed over the set of piecewise constant functions on a xed grid. The only di erence with this last work is that we do not look for some maximal solution of (42), and hence do not need the strictly concave penalization introduced in [START_REF] Buttazzo | On the selection of maximal Cheeger sets[END_REF] 1 . Discrete problem. Since every solution of (42) has its support included in some ball 2 , we can solve (42) in [-R, R] 2 (with Dirichlet boundary conditions) for a su ciently large R > 0. Let N be a positive integer and h def.

= 2R/N . We denote E h the set of N by N matrices. For every matrix u = (u i,j

) (i,j)∈[1,N ] 2 ∈ E h we de ne ∂ h x u i,j def. = u i+1,j -u i,j , ∂ h y u i,j def. = u i,j+1 -u i,j , (45)
for all (i, j) ∈ [0, N ] 2 , with the convention u i,j = 0 if either i or j is in {0, N + 1}. The isotropic discrete total variation is then de ned by

TV h (u) def. = h N i=0 N j=0 ||∇ h u i,j || 2 = h ∇ h u 2,1 with ∇ h u i,j def. = ∂ h x u i,j , ∂ h y u i,j .
We then solve the following discretized version of (42) for increasingly small values of h > 0

min u∈E h h 2 η h , u s.t. TV h (u) ≤ 1 , (46)
where

η h = 1 h 2 ´Ch i,j η (i,j)∈[1,N ] 2 and (C h i,j ) (i,j)∈[1,N ] 2 is a partition of [-R, R] 2 composed of squares of equal size, i.e. C h i,j def. = [-R + (i -1)h, -R + ih] × [-R + (j -1)h, -R + jh] .
For convenience reasons, we also use the above expression to de ne C h i,j if i or j belongs to {0, N + 1}.

Optimization algorithm. To solve ( 46), we propose to use the primal-dual algorithm introduced in [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF]: we take (τ, σ) such that τ σ D 2 2 < 1 with D def.

= h∇ h and de ne

     φ n+1 = prox σ • 2,∞ (φ n + σ Dū n ) , u n+1 = (u n -τ D * φ n+1 ) -τ h 2 ηh , ūn+1 = 2 u n+1 -u n , (47)
where prox σ • 2,∞ is given by:

prox σ • 2,∞ (φ) = φ -σ proj { • 2,1 ≤1} φ σ .
The projection onto the (2, 1)-unit ball can be computed e ciently (see [Condat, 2016]).

Convergence as h → 0. The following proposition shows that, when the grid becomes ner, solutions of (46) converge to a solution of (42). Its proof is almost the same as the one of [Carlier et al., 2009, Theorem 4.1]. Since the latter however gives a slightly di erent result about the minimization of a quadratic objective (linear in our case) on the total variation unit ball, we decided to include it for the sake of completeness.

Proposition 3.13

Let u h be the piecewise constant function on

(C h i,j ) (i,j)∈[1,N ] 2 , extended to 0 outside [-R, R
] 2 associated to a solution of (46). Then there exists a (not relabeled) subsequence converging strongly in L 1 (R 2 ) and weakly in L 2 (R 2 ) to a solution u * of (42) when h → 0. Moreover, we have Du h * Du.

Proof : First, let us stress that for any function v that is piecewise constant on

(C i,j ) (i,j)∈[1,N ] 2
and that is equal to

0 outside [-R, R] 2 , we have TV(v) = h ∇ h v 1,1
where by abuse of notation ∇ h v is given by (45) with v i,j the value of v in C i,j . Hence TV h (u h ) ≤ 1 for all h implies that TV(u h ) (and hence u h L 2 ) is uniformly bounded in h. There hence exists a (not relabeled) subsequence that converges strongly in L 1 loc (R 2 ) and weakly in L 2 (R 2 ) to a function u, with moreover Du h * Du.

Let us now take

φ = (φ (1) , φ (2) ) ∈ C ∞ c (R 2 , R 2 ) such that ||φ|| ∞ ≤ 1.
The weak-* convergence of the gradients give us that

ˆR2 φ • dDu = lim h→0 ˆR2 φ • dDu h = lim h→0 N i=0 N j=0 ´Ch i,j ∩C h i+1,j φ (1) dH 1 ´Ch i,j ∩C h i,j+1 φ (2) dH 1 • ∇ h u h i,j .
One can moreover show there exists C > 0 such that for h small enough and all (i, j) we have:

ˆCh i,j ∩C h i+1,j φ (1) dH 1 -h φ (1) (x h i+1,j+1 ) ≤ Ch 2 , ˆCh i,j ∩C h i,j+1 φ (2) dH 1 -h φ (2) (x h i+1,j+1 ) ≤ Ch 2 , with x i,j def.
= (-R + i h, -R + j h). We use the above inequalities and the fact φ(x) ≤ 1 for all x to obtain the existence of C > 0 such that for h small enough and for all (i, j) we have:

´Ch i,j ∩C h i+1,j φ (1) dH 1 ´Ch i,j ∩C h i,j+1 φ (2) dH 1 2 ≤ h √ 1 + C h .
This nally yields

N i=0 N j=0 ´Ch i,j ∩C h i+1,j φ (1) dH 1 ´Ch i,j ∩C h i,j+1 φ (2) dH 1 • ∇ h u h i,j ≤ N i=0 N j=0 h √ 1 + C h ∇ h u h i,j = √ 1 + C h TV h (u h ) , which gives ˆR2 φ • dDu ≤ lim sup h→0 √ 1 + C h TV h (u h ) ≤ 1 .
We now have to show that

∀v ∈ L 2 (R 2 ), TV(v) ≤ 1 =⇒ ˆR2 η u ≤ ˆR2 η v . Let v ∈ C ∞ ([-R, R]
2 ) be such that TV(v) ≤ 1. We de ne

v h def. = v i + 1 2 h, j + 1 2 h (i,j)∈[1,N ] 2 .
One can then show that lim h→0 TV h (v h ) = TV(v) = 1 , so that for every δ > 0 we have TV h v h 1+δ ≤ 1 for h small enough. Now this yields

ˆ[-R,R] 2 η u = lim h→0 ˆ[-R,R] 2 η u h ≤ lim h→0 ˆ[-R,R] 2 η v h 1 + δ = ˆ[-R,R] 2 η v 1 + δ .
Since this holds for all δ > 0 we get that

ˆ[-R,R] 2 η u ≤ ˆ[-R,R] 2 η v . (48) Finally, if v ∈ L 2 (R 2 ) is such that v = 0 outside [-R, R] 2
and TV(v) ≤ 1, by standard approximation results (see [Ambrosio et al., 2000, remark 3.22]) we also have that (48) holds, and hence u solves (42). Finally, since u solves (42), its support is included in [-R, R] 2 , which shows the strong L 1 loc (R 2 ) convergence of (u h ) towards u * in fact implies its strong L 1 (R 2 ) convergence.

Extraction of a level set. Since we are interested in nding a simple set E that approximately solves (41), and now have a good way of approximating solutions of (42), we make use of the following result, which is a direct consequence of Proposition 1.3.

Proposition 3.14 Let u be a solution of (42). Then the level sets of u are such that for all t ∈ R * with |U (t) | > 0, the set U (t) solves (41).

If we have v h converging strongly in L 1 (R 2 ) to a solution v * of (42), then for almost every t ∈ R we have that lim

h→0 V (t) h V (t) * = 0 .
The above results hence show we can construct a sequence of sets (E k ) k≥0 such that |E k E * | converges to 0, with E * a solution of (41). However, this convergence only implies that lim sup

k→∞ J(E k ) ≤ J(E * ) ,
and given that E k is a union of squares this inequality is likely to be strict, with the perimeter of E k not converging to the perimeter of E * . From the last paragraph of Section 1.3.1, we know we have to design a numerical method that allows to nd a set at which the value of J is arbitrarily close to J(E * ). This hence motivates the introduction of the re nement step described in the next subsection.

Obtaining a simple polygon. As a nal remark, we note that, even for k large enough, E k could be non-simple. However, using the notations of (59) and ??, since for every set of nite perimeter E, we have that J(E) is a convex combination of the

J(int(γ + i )) i∈I , J(int(γ - i,j )) i∈I,j∈J i
, there is a simple set F in the decomposition of E such that J(F ) ≥ J(E). In practice, such a set can be found by extracting all the contours of the binary image 1 E , and nding the one with highest objective value. This procedure guarantees that the output of the xed grid step is a simple polygon. We stress that in all our experiments, v h is close to being (proportional to) the indicator of a simple set for h large enough, so that its non-trivial level sets are all simple.

O

As explained at the beginning of Section 2, our approach for approximating polygonal Cheeger sets essentially consists in optimizing the vertices of an initial polygon with a rst order method. This is in spirit very similar to so-called shape gradient algorithms (see e.g. [START_REF] Allaire | Chapter 1 -Shape and topology optimization[END_REF]).

We iteratively construct a sequence of polygons by nding at each step a displacement of steepest ascent for J, along which the vertices of the previous polygon are moved. This displacement is found by exploiting Lemmas 3.15 and 3.16 below.

Fixed sign assumption. Given E x 0 the polygon produced by the xed grid initizalization, we only aim at locally optimizing the objective J around x 0 . We can therefore assume the sign of the weighted area term A remains constant and equal to def.

= sign(A(x 0 )) during the optimization. In the following, in order to simplify the exposition, we only consider the case = 1, and consequently seek to maximize

J(x) = A(x) P (x) .
Lemma 3.15 If x ∈ X n we have:

P (E x+h ) = P (E x ) - n i=1 h i , τ + i + τ - i + o ( h ) , (49)
where τ

+ i def. = x i+1 -x i x i+1 -x i and τ - i def. = x i-1 -x i x i-1 -x i .
Proof : If h is small enough we have:

P (E x+h ) = n i=1 x i+1 -x i + h i+1 -h i = n i=1 x i+1 -x i + h i+1 -h i 2 = n i=1 x i+1 -x i 1 + x i+1 -x i , h i+1 -h i x i+1 -x i 2 + o ( h ) = P (E x ) + n i=1 τ + i , h i+1 -h i + o ( h ) ,
and the result follows by re-arranging the terms in the sum.

Lemma 3.16 Let x ∈ X n and η ∈ C 0 (R 2 ). Then we have

ˆEx+h η = ˆEx η + n i=1 h i , w - i ν - i + w + i ν + i + o ( h ) (50)
where ν + i , ν - i are respectively the outward unit normal to E x on [x i , x i+1 ], [x i-1 , x i ] and

w + i def. = x i+1 -x i ˆ1 0 η((1 -t)x i + tx i+1 )(1 -t)dt , w - i def. = x i-1 -x i ˆ1 0 η((1 -t)x i + tx i-1 )(1 -t)dt .
Proof : Our proof relies on the following identity (see Lemma 3.11 for a proof of a closely related formula):

ˆEx η = sign n i=1 det(x i x i+1 ) n i=1 ω(x i , x i+1 ) ,
where ω(a 1 , a 2 ) def.

= det(a 1 a 2 ) ´T1 η((a 1 a 2 ) y) dy and

T 1 def. = (α, β) ∈ (R + ) 2 α + β ≤ 1 is the unit triangle. Assuming η ∈ C 1 (R 2
) and denoting adj(A) the adjugate of a matrix A we get:

ω(a 1 + h 1 , a 2 + h 2 ) = ω(a 1 , a 2 ) + det(a 1 a 2 ) ˆT1 ∇η((a 1 a 2 ) y) • ((h 1 h 2 )y) dy + tr adj(a 1 a 2 ) T (h 1 h 2 ) ˆT1 η((a 1 a 2 ) y) dy + o( h ) = ω(a 1 , a 2 ) + sign(det(a 1 a 2 )) ˆOa1a2 ∇η(y) • ((h 1 h 2 ) (a 1 a 2 ) -1 ) y) dy + tr adj(a 1 a 2 ) T (h 1 h 2 ) |det(a 1 a 2 )| ˆOa1a2 η(y) dy + o( h ) .
Denoting g(y)

def.

= (h 1 h 2 )(a 1 , a 2 ) -1 y, we obtain:

ω(a 1 + h 1 , a 2 + h 2 ) = ω(a 1 , a 2 ) + sign(det(a 1 a 2 )) ˆOa1a2 [∇η • g + η divg] + o( h ) = ω(a 1 , a 2 ) + sign(det(a 1 a 2 )) ˆ∂(Oa1a2) η (g • ν Oa1a2 ) dH 1 + o( h ) ,
where we used Gauss-Green theorem to obtain the last equality. Now if h is small enough then

n i=1 det(x i + h i x i+1 + h i+1 ) and n i=1 det(x i x i+1 )
have the same sign, so that, de ning

g i : y → ((h i h i+1 )(x i x i+1 ) -1 y) we get d ˆE• η (x) . h = n i=1 sign(det(x i x i+1 )) ω i , with def. = sign n i=1 det(x i x i+1 ) and ω i def. = ˆ∂ * (Oxixi+1) η (g i • ν Oxixi+1 ) dH 1 .
Then one can decompose each integral in the sum and show the integrals over [0, x i ] cancel out each other, which allows to obtain

d ˆE• η (x) . h = n i=1 ˆ[xi,xi+1] η (g i • ν i ) dH 1 . But now if y ∈ [x i , x i+1 ] then (x i x i+1 ) -1 y = 1 x i+1 -x i y -x i+1 y -x i ,
and the result follows by re-arranging the terms in the sum. One can then use an approximation argument as in [Maggi, 2012, Proposition 17.8] to show it also holds when η is only continuous.

Remark 3.17

The results given in Lemmas 3.15 and 3.16 are in fact closely linked to the notion of shape derivative1 . If E x is a simple polygon with n vertices, there exists a constant C > 0 such that, given any displacement h of the vertices of E x , one can extend

h to θ ∈ W 1,∞ (R 2 , R 2 ) with θ 1,∞ ≤ C h and such that E x+h = (Id + θ)(E x ) and, for every i ∈ {1, ..., n} ∀t ∈ [0, 1], θ((1 -t)x i + tx i+1 ) = (1 -t)θ(x i ) + tθ(x i+1 ) = (1 -t)h i + th i+1 . ( 51
)
The expressions of the shape derivatives of the perimeter and area (see [Henrot and Pierre, 2018, Chapter 5] for more details, including precise assumptions on the regularity of E and η) yield

P ((Id + θ)(E)) = P (E) + ˆ∂ * E H E • θ dH 1 + o( θ 1,∞ ) , ˆ(Id+θ)(E) η = ˆE η + ˆ∂ * E η θ • ν E dH 1 + o( θ 1,∞ ) ,
where H E denotes the curvature vector of E. To recover the formulas given in the above lemmas, one simply needs to use (51) and the fact

H Ex = - n i=1 (τ + i + τ - i )δ x i with τ + i def. = x i+1 -x i x i+1 -x i and τ - i def. = x i-1 -x i x i-1 -x i .
Considering Lemmas 3.15 and 3.16, we construct our sequence of polygons as follows:

given x t ∈ X n and a step size α t , we de ne the next iterate by x t+1 = x t + α t θ t with θ t j def.

= 1 P (E x t ) θ t area,j -´Ex t η P (E x t ) θ t per,j , θ t area,j def.

= w t- j ν t- j + w t+ j ν t+ j , θ t per,j def.

= -(τ t+ j + τ t- j ) ,

where, for all j, ν t+ j and ν t- j are the outward unit normals of the two sides stemming from x t j , and

τ t+ j def. = x i+1 -x i x i+1 -x i , w t+ j def. = x i+1 -x i ˆ1 0 η((1 -t)x i + tx i+1 )(1 -t)dt , τ t- j def. = x i -x i-1 x i -x i-1 , w t+ j def. = x i-1 -x i ˆ1 0 η((1 -t)x i + tx i-1 )(1 -t)dt .
Lemmas 3.15 and 3.16 in fact show that the displacement θ t we apply to the vertices of

E x t is such that lim α→0 + J(E x t +αθ ) -J(E x t ) α = θ t , θ , (53)
i.e. that it is the displacement of steepest ascent for J at E x t .

Numerical computation of θ t . To compute the integral of η over E x t , we integrate η on each triangle of a su ciently ne triangulation of E x t (this triangulation must be updated at each iteration, and sometimes re-computed from scratch to avoid the presence of ill-shaped triangles).

The integral of η on a triangle and w t+ j , w t- j are computed using standard numerical integration schemes for triangles and line segments. If |T| denotes the number of triangles in the triangulation of E x t , |S T | (resp. |S L |) the number of points used in the numerical integration scheme for triangles (resp. line segments), the complexity of each iteration is of order

O (m (|T| |S T | + n |S L |)).
Comments. Two potential concerns about the above procedure are whether the iterates remain simple polygons (i.e. x t ∈ X n for all t) and whether they converge to a global maximizer of J over P n . We could not prove that the iterates remain simple polygons along the process, but since the initial polygon can be taken arbitrarily close to a simple set solving (41) (in terms of the Lebesgue measure of the symmetric di erence), we do not expect nor observe in practice any change of topology during the optimization. Moreover, even if J could have non-optimal critical points1 , the above initialization allows us to start our local descent with a polygon that hopefully lies in the basin of attraction of a global maximizer.

Critical points. If the sequence of polygons de ned above converges to a simple polygon E x , then E x is such that ∀j ∈ {1, ..., n}, w

+ j ν + j + w - j ν - j = ´Ex η P (E x ) (-τ + i -τ - i ) . (54)
This can be seen as a discrete version of the following rst order optimality condition for solutions of (41):

η = ´E η P (E) H E on ∂ * E . (55)
Note that ( 55) is similar to the optimality condition for the classical Cheeger problem (i.e. with η = 1 and the additional constraint E ⊆ Ω), namely H E = P (E)/ |E| in the free boundary of E (see [Parini, 2011, Proposition 2.4]).

T

In this section, we study how generalized Cheeger sets and their polygonal approximations compare in a particular setting, where the weight function η is radial. To be more precise, we assume that there exists a strictly decreasing nonnegative function η : [0, +∞) → R + such that η(x) = η( x ) for almost every x ∈ R 21 .

We rst prove that, in this setting, (generalized) Cheeger sets associated to η are all disks centered at the origin. Then, we prove that the maximizers of J over P n are regular and inscribed in a circle centered at the origin when n ∈ {3, 4}, and conjecture that this result holds for every n ≥ 3.

To enforce the uniqueness of the generalized Cheeger set, we may also invoke the following assumption, which is for example satis ed by η : r → exp(-r 2 /(2σ 2 )) for all σ > 0.

Assumption 1. The function η is C1 . Moreover, de ning f : r → r η(r), we have f (r) r→+∞ ----→ 0. Finally, there exists r 0 > 0 such that f (r) > 0 for all r < r 0 and f (r) < 0 for all r > r 0 .

Description of generalized Cheeger sets

To prove the above-mentioned result, we rely on Steiner symmetrization. Let us brie y recall its de nition and main property. If E is a set of nite perimeter with nite measure, ν ∈ S 1 and z ∈ R, we denote

E ν,z def. = {t ∈ R | z ν + t ν ⊥ ∈ E} .
The Steiner symmetrization of E with respect to the line through the origin and directed by ν, denoted E s ν , is then de ned by

E s ν def. = {x ∈ R 2 | | x, ν ⊥ | ≤ |E ν, x,ν |/2} .
The fundamental property of Steiner symmetrization is that it preserves volume and does not increase perimeter (see [Maggi, 2012, section 14.1] for more details).

To prove our result, we rst prove two useful lemmas.

Lemma 3.18 Let f : R → R + be even and strictly decreasing on R + . Then for every measurable set A such that |A| < +∞ we have:

ˆA f ≤ ˆAs f ,
where A s def. = (-|A|/2, |A|/2). Moreover, equality holds if and only if |A A s | = 0.

Proof : We have ˆA f = ˆ+∞ 0 |{f 1 A ≥ t}| dt = ˆ+∞ 0 |{f ≥ t} ∩ A| dt .
Since f is even, for all t > 0, there exists α such that {f ≥ t} = [-α, α], so that we have

|{f ≥ t} ∩ A| = |[-α, α] ∩ A| ≤ min(2α, |A|) = |[-α, α] ∩ [-|A|/2, |A|/2]| = |{f ≥ t} ∩ A s | . Hence ˆA f ≤ ˆ+∞ 0 {f ≥ t} ∩ A s dt = ˆAs f . Now if |A A s | > 0, since |A| = |A s | then |A \ A s | = |A s \ A| > 0 and we have ˆAs f = ˆA∩A s f + ˆAs \A f > ˆA∩A s f + f (|A|/2) |A s \ A| ≥ ˆA∩A s f + ˆA\A s f = ˆA f ,
which proves the second part of the result.

Lemma 3.19 Let E ⊂ R 2 be a set of nite perimeter with 0 < |E| < +∞. Then for any ν ∈ S 1 we have ´Es

ν η P (E s ν ) ≥ ´E η P (E) ,
with equality if and only if |E E s ν | = 0.

Proof : From [Maggi, 2012, theorem 14.4] we know that we have P (E s ν ) ≤ P (E). We now perform a change of coordinates in order to have

E s ν = {(x 1 , x 2 ) ∈ R 2 | |x 2 | ≤ |E x1 |/2} with E x1 def. = {x 2 ∈ R | (x 1 , x 2 ) ∈ E} . Now, we have ˆE η = ˆ+∞ -∞ ˆ+∞ -∞ η(x 1 , x 2 ) 1 E (x 1 , x 2 ) dx 2 dx 1 = ˆ+∞ -∞ ˆEx 1 η(x 1 , •) dx 1 , with E x1 = {x 2 ∈ R | (x 1 , x 2 ) ∈ E}.
For almost every x 1 ∈ R we have that E x1 is measurable, has nite measure, and that η(x 1 , •) is nonnegative, even and strictly decreasing on R + . We can hence apply Lemma 3.18 and get that

ˆE η ≥ ˆ+∞ -∞ ˆ(Ex 1 ) s η(x 1 , •) d x1 = ˆEs ν η . (56) Moreover, if |E E s ν | > 0, then since |E E s ν | = ˆ+∞ 0 ˆ+∞ 0 |1 E (x 1 , x 2 ) -1 E s ν (x 1 , x 2 )| dx 2 dx 1 = ˆ+∞ 0 ˆ+∞ 0 1 Ex 1 (x 2 ) -1 (Ex 1 ) s (x 2 ) dx 2 dx 1 = ˆ+∞ 0 E x1 E s x1 dx 1 , we get that L 1 ({x 1 ∈ R | |E x1 (E x1 ) s | > 0})
> 0 and hence that ( 56) is strict.

Using the above lemmas, we may now state the following result.

Proposition 3.20

The generalized Cheeger sets associated to η are disks centered at the origin. Moreover, under Assumption 1, the optimal set is unique and its radius is [Maggi, 2012, Section 14.2], we get that E * is a convex set which is invariant by re ection with respect to any line through the origin, and hence that E * is a disk centered at the origin.

R * = argmax R>0 1 R ˆR 0 rη(r) dr . Proof : If E ⊂ R 2 is
Let us de ne

G(R) def. = 2πJ(B(0, R)) = 1 R ˆR 0 rη(r) dr .
Under Assumption 1 G can be extended to a C 1 function on R + with G(0) = 0, G (0) = η(0)/2 and (denoting f : r → rη(r) as in Assumption 1):

∀R > 0, G (R) = Rf (R) -´R 0 f (r) dr R 2 .
Now de ning F : R → Rf (R)-´R 0 f (r) dr we see that F (R) = Rf (R) for all R > 0, so that F and f have the same variations. Using Assumption 1, we hence have that F is strictly increasing on [0, R 0 ] and strictly decreasing on [R 0 , +∞). Moreover, F (0) = 0 and, since f (R)

R→+∞

-----→ 0, we also get F (R) < 0 for R large enough. Therefore, there exists R * > 0 such that F (and hence G ) is strictly positive on (0, R * ) and strictly negative on (R * , +∞). We hence obtain that R * is the unique maximizer (and the unique critical point) of G. In practice, instead of solving (41), we look for an element of P n (a simple polygon with at most n sides) maximizing J, for some given integer n ≥ 3. It is hence natural to investigate the proximity of these optimal polygons with the "continuous" Cheeger sets. Solving classical geometric variational problems restricted to the set of n-gons is involved, as the Steiner symmetrization procedure might increase the number of sides [Polya and Szegö, 1951, Sec. 7.4]. However, using a trick from Pólya and Szegö, one may prove: Proposition 3.21

Let n ∈ {3, 4}. Then all the maximizers of J over P n are regular and inscribed in a circle centered at the origin.

Proof :

Triangles (n = 3): let E * be a maximizer of J among triangles. Then the Steiner symmetrization of E * with respect to any of its heights through the origin (see Figure 14) is still a triangle, and Lemma 3.19 ensures it has a higher energy except if this operation leaves E * unchanged. As a consequence, E * must be symmetric with respect to all its heights through the origin. This shows E * is equilateral and inscribed in a circle centered at the origin.

Quadrilaterals (n = 4): we notice that if E is a simple quadrilateral, then its Steiner symmetrization with respect to any line perpendicular to one of its diagonals (see Figure 15) is still a simple quadrilateral. We can then proceed exactly as for triangles to prove any maximizer E * of J over P 4 is symmetric with respect to every line through the origin and perpendicular to one of its diagonals. This shows E * is a rhombus centered at the origin. We can now symmetrize with respect to any line through the origin perpendicular to one of its sides to nally obtain that E * must be a square centered at the origin.

Our proof does not extend to n ≥ 5, but the following conjecture is natural:

Conjecture 1

The result stated in Proposition 3.21 holds for all n ≥ 3.

For n ∈ {3, 4} or if Conjecture 1 holds, it remains to compare the radius of optimal polygons with the one of "continuous" Cheeger sets. Let us x a regular n-gon inscribed in the circle of radius R centered at 0, and denote it B n (0, R). Abusing notation, we de ne J(R)

def.

= J(B(0, R)) and J n (R) def.

= J(B n (0, R)), and can now state the following result.

J n -J|| ∞ = O 1 n 2 .
Moreover, under Assumption 1 and assuming f is C 2 with f (r 0 ) < 0, then for n large enough we have that J n has a unique maximizer R * n and

|R * n -R * | = O 1 n . Proof : Recalling that P (B n (0, R)) = 2πR sin(π/n) π/n
we have:

|J(R) -J n (R)| = 1 2πR ˆB(0,R) η - π/n 2πR sin(π/n) ˆBn(0,R) η = 1 - π/n sin(π/n) 1 2πR ˆB(0,R) η + π/n 2πR sin(π/n) ˆB(0,R)\Bn(0,R) η ≤ 1 - π/n sin(π/n) J(R) + π/n 2πR sin(π/n) η ∞ |B(0, R) \ B n (0, R)| .
Now, we de ne

R n (θ) def. = R cos (π/n) cos ((θ mod 2π/n) -π/n) ,
so that in polar coordinates an equation of ∂B n (0, R) is given by r(θ) = R n (θ). We obtain:

|B(0, R) \ B n (0, R)| = ˆ2π 0 ˆR Rn(θ) r dr dθ ≤ 2π sup θ∈[0,2π] ˆR Rn(θ) r dr = π sup θ∈[0,2π] R 2 -R 2 n (θ) = π(R 2 -R 2 cos 2 (π/n)) .
We hence obtain

J n -J ∞ ≤ J ∞ 1 - π/n sin(π/n) + R η ∞ π/n sin(π/n) (1 -cos(π/n))
and we can nally conclude regarding the rst statement. Now, de ning α n (s)

def.

= cos(π/n)/cos(πs/n), we have:

J n (R) = π/n 2πR sin(π/n) ˆ2π 0 ˆRn(θ) 0 rη(r) dr dθ = π/n 2πR sin(π/n) n ˆ2π/n 0 ˆR cos(π/n) cos(θ-π/n) 0 rη(r) dr dθ = π/n 2πR sin(π/n) 2n ˆπ/n 0 ˆR cos(π/n) cos(θ) 0 rη(r) dr dθ = π/n R sin(π/n) ˆ1 0 ˆRαn(s) 0 rη(r) dr ds = π/n sin(π/n) 1 R ˆR 0 r ˆ1 0 (α n (s)) 2 η(rα n (s)) ds dr .
We can now de ne f n : r → r ´1 0 (α n (s)) 2 η(rα n (s)) ds and proceed as in Proposition 3.20:

showing f n is strictly positive on (0, r 1 ) and strictly negative on (r 1 , +∞) for some r 1 > 0 is su cient to prove J n has a unique maximizer (its unique critical point). Now, we have:

∀r ∈ R + , f n (r) = ˆ1 0 (α n (s)) 2 (η(rα n (s)) + rα n (s)η (rα n (s))) ds . The image of [0, 1] by s → rα n (s) is [rcos(π/n), r]. Since r → η(r) + rη (r) = (rη) (r)
is strictly positive on (0, r 0 ) and strictly negative on (r 0 , +∞), we get that f n is strictly positive on (0, r 0 ) and strictly negative on (r 0 /cos(π/n), +∞). It hence remains to investigate the sign of f n on the interval [r 0 , r 0 /cos(π/n)]. But since f is C 2 and f (r 0 ) < 0 there exists > 0 such that f (r) < 0 for every r ∈ (r 0 -, r 0 + ). For n large enough, we hence have:

[r 0 cos(π/n), r 0 /cos(π/n)] ⊂ (r 0 -, r 0 + ) .

This implies that for all r ∈ [r 0 , r 0 /cos(π/n)] we have rα n (s) ∈ (r 0 -, r 0 + ), and hence that f n (r) < 0. This nally shows there exists r 1 > 0 such that f n is strictly positive on (0, r 1 ) and strictly negative on (r 1 , +∞), and the result follows as in Proposition 3.20. Now R * and R * n are respectively the unique solutions of F (0, R) = 0 and F (π/n, R) = 0 with

F (t, R) def. = ˆR 0 f t (r) dr -Rf t (R) , f t (r) def. = r ˆ1 0 α(t, s) 2 η(rα(t, s))ds , α(t, s) def.

= cos(t) cos(ts)

.

Moreover, F is C 1 in a neighborhood of (0, R * ) and ∂ R F (0, R * ) = -R * f (R * ).
From Assumption 1 we know the unique R > 0 such that f (R) = 0 is R 0 , which (from the proof of Proposition 3.20) is di erent from R * . We can hence apply the implicit function theorem and obtain that

|R * n -R * | = |∂ t F (0, R * )| |∂ R F (0, R * )| 1 n + o 1 n ,
which yields the result.

N

Before assessing the practical performance of Algorithm 3, we brie y explain how we implement the sliding step (Line 14).

L

This section is dedicated to the implementation of the sliding step (Line 14 of Algorithm 3). Our approach is exactly the same as in Section 2.4, that is we optimize the vertices of a set of polygons and their corresponding amplitudes using a rst order algorithm. To be more precise, we perform a gradient descent on the mapping

(a, x) → G N ( , a, E x ) , (57)
where E x = (E x 1 , ..., E x N ). Given a step size α t , a vector a t ∈ R N and x t 1 , ..., x t N in X n , we de ne

u t def. = N i=1 a t i 1 E x t i , η t def. = - 1 λ Φ * (Φu t -y) ,
and perform the following update:

a t+1 i def. = a t i -α t h t i , h t i def. = P E x t i -i ˆEx t i η t , x t+1 i,j def. = x t i,j -α t θ t i,j , θ t i,j def. = a t i -(τ t+ i,j + τ t- i,j ) + i θ t data,i,j , θ t data,i,j def.
= Φu ty, w t+ i,j ν t+ i,j + Φu ty, w t- i,j ν t- i,j , where τ t± i,j , ν t± i,j are de ned as in Section 2.4 and

w t+ i,j def. = x t i,j+1 -x t i,j ˆ1 0 ϕ((1 -t) x t i,j + t x t i,j+1 ) (1 -t) dt , w t- i,j def. = x t i,j-1 -x t i,j ˆ1 0 ϕ((1 -t) x t i,j + t x t i,j-1 ) (1 -t) dt .
Using the notations of Section 2.4, the complexity of each iteration is of order

O (N m (|T| |S T | + n |S L |)) .
Comments. We rst stress that the above update is similar to the evolution formally described in (36). Now, unlike the local optimization we perform to approximate Cheeger sets, the sliding step may tend to induce topology changes (see Section 3.3 for an example). This is of course linked to the possible appearance of singularities mentioned in Section 1.4.1. Typically, a simple set may tend to split in two simple sets over the course of the descent. This is a major di erence (and challenge) compared to the sliding steps used in sparse spikes recovery (where the optimization is carried out over the space of Radon measures) [START_REF] Bredies | Inverse problems in spaces of measures[END_REF], Boyd et al., 2017, Denoyelle et al., 2019]. This phenomenon is closely linked to topological properties of the faces of the total (gradient) variation unit ball: its extreme points do not form a closed set for any reasonable topology (e.g. the weak L 2 (R 2 ) topology), nor do its faces of dimension d ≤ k for any k ∈ N. As a result, when moving continuously on the set of faces of dimension d = k, it is possible to "stumble upon" a point which only belongs to a face of dimension d > k.

Handling (or not) topology changes. Our current implementation does not allow to handle these topology changes in a consistent way, and nding a way to deal with them "o -the-grid" is an interesting avenue for future research. It is important to note that not allowing topological changes during the sliding step is not an issue, since all convergence guarantees hold as soon as the output of the sliding step decreases the energy more than the standard update. One can hence stop the local descent at any point before any change of topology occurs, which avoids having to treat them. Still, in order to yield iterates that are as sparse as possible (and probably to decrease the objective as quickly as possible), it seems preferable to allow topological changes.

R

Here, we investigate the practical performance of Algorithm 3. We focus on the case where Φ is a sampled Gaussian convolution operator, i.e.

∀x ∈ R 2 , ϕ(x) = exp - ||x -x i || 2 2σ 2 m i=1
for a given σ > 0 and a sampling grid (x i ) m i=1 . The noise is drawn from a multivariate Gaussian with zero mean and isotropic covariance matrix τ 2 I m . We take λ of the order of 2 log(m) τ 2 .

Numerically certifying that a given function is an approximate solution of (P λ (y)) is di cult. However, as the sampling grid becomes ner, Φ tends to the convolution with the Gaussian kernel, which is injective. Relying on a Γ-convergence argument, one may expect that if u 0 is a piecewise constant image and w is some small additive noise, the solutions of (P λ (y)) with y = Φu 0 + w are all close to u 0 , modulo the regularization e ects of the total variation.

Comparison with other algorithms. We assess the performance of our algorithm by comparing its output to that of a primal dual algorithm minimizing a discretized version of (P λ (y)) on a pixel grid, where the total variation term is replaced by the discrete isotropic total variation or the dual-based discretization of [START_REF] Hintermüller | Functionalanalytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF], Condat, 2017]. To minimize discretization artifacts, we arti cially introduce a downsampling in the forward operator, so that the reconstruction is performed on a grid four times larger than the sampling one.

Experiments. Our rst experiment consists in recovering a function u 0 that is a linear combination of three indicator functions (see Figures 16 and17). During each of the three iterations required to obtain a good approximation of u 0 , a new atom is added to its support. One can see the sliding step is crucial: the large atom on the left, added during the second iteration, is of the xed grid method using the isotropic and Condat's total variation signi cantly re ned during the sliding step of the third iteration, when enough atoms have been introduced. The second experiment (see Figure 18) consists in recovering the indicator function of a set with a hole (which can also be seen as the sum of two indicator functions of simple sets). The support of u 0 and its gradient are accurately estimated. Still, the typical e ects of total variation regularization are noticeable: corners are slightly rounded, and there is a "loss of contrast" in the eye of the pacman. The third experiment (Figure 19) also showcases the rounding of corners, and highlights the in uence of the regularization parameter: as λ decreases, the curvature of the edge set increases. Finally, we provide in Fig.

Figure 20 the results of an experiment on a more challenging task, which consists in reconstructing a natural grayscale image.

Choice of parameters. The number of observations in the rst experiment is 60 × 60, 75 × 75 in the second and third ones, and 64 × 64 in the last one. In all experiments, we solved (46) on a grid of size 80 × 80. In both local descent steps (for approximating Cheeger sets and for the sliding step), the simple polygons have a number of vertices of order 30 times the length of their boundary (100 for the last experiment), and the maximum area of triangles in their inner mesh is 10 -2 (the domain being a square of side 1). The inner triangulation of a simple polygon is obtained by using Richard Shewchuk's Triangle library. The boundary of the polygons are resampled every 30 iterations. Line integrals are computed using the Gauss-Patterson scheme of order 3 (15 points) and integrals on triangles using the Hammer-Marlowe-Stroud scheme of order 5 (7 points).

T

Here, we illustrate the changes of topology that may occur during the sliding step (Line 14 of Algorithm 3). All relevant plots are given in Figure 21. The unknown function (see (a)) is the sum of two indicator functions:

u 0 = 1 B((-1,0),0.6) + 1 B((1,0),0.6) ,
and observations are shown in (b). The Cheeger set computed at Line 5 of the rst iteration covers the two disks (see (c)).

In this setting, our implementation of the sliding step converges to a function similar to (f)1 , and we obtain a valid update that decreases the objective more than the standard Frank-Wolfe update. The next iteration of the algorithm will then consist in adding a new atom to the k] (k = 1, 4) produced by Algorithm 3, outputs of the xed grid method using the isotropic and Condat's total variation approximation, with negative amplitude, so as to compensate for the presence of the small bottleneck. However, it seems natural that the support of (f) should split into two disjoint simple sets, which is not possible with our current implementation. To investigate what would happen in this case, we manually split the two sets (see (g)) and let them evolve independently. The support of the approximation converges to the union of the two disks, which produces an update that decreases the objective even more than (f).

C

In this part, we proposed an iterative algorithm for solving (P λ (y)), which is based on the conditional gradient (or Frank-Wolfe) algorithm. We studied the convergence of its iterates, and presented a natural analog of the "sliding step" introduced in [START_REF] Denoyelle | The Sliding Frank-Wolfe Algorithm and its Application to Super-Resolution Microscopy[END_REF]. Unlike in the latter work, we left the theoretical analysis of this step open, and in particular could not prove that a so-called critical point of the objective can be reached at every iteration.

We then discussed the practical implementation of this theoretical algorithm, and presented two oracles to carry out its main operations. The rst allows to approximate generalized Cheeger sets, and builds up on the method proposed in [START_REF] Carlier | Approximation of maximal Cheeger sets by projection[END_REF]. It consists in optimizing the vertices of a polygon with a shape gradient-like algorithm, initialized with the output of the latter method. To support our approach, we also proved the existence of generalized Cheeger sets in the class of polygons with a given number of sides, and investigated their proximity with their continuous counterpart in a radial setting. The second oracle we proposed allows to implement the sliding step mentioned above. We discussed the topological changes that can appear during the resulting evolution, and argued that handling them in a consistent way is natural but challenging.

Finally, we assessed the performance of the propose numerical method on a few recovery examples, providing comparisons with xed-grid methods. Our approach is particularly suited to reconstruct simple images (i.e. which are the superposition of a few simple shapes), whereas on more complex natural images, existing methods perform better. 110 ensure a notion of pre-certi cate is relevant, a natural way would be to nd out if, as proved in [START_REF] Duval | Exact Support Recovery for Sparse Spikes Deconvolution[END_REF] in the sparse spikes setting, the fact it is a true certi cate is necessary for (noiseless) exact support recovery.

Error bounds. The noise robustness analysis we presented in Part 2 is not quantitative. It does not provide bounds on the distance between solutions of (P 0 (y 0 )) and (P λ (y 0 + w)) in terms of the noise level w H . In the sparse spikes setting, such bounds were derived by analyzing the behaviour of dual certi cates in the so-called far and near regions (see for instance [START_REF] Azaïs | Spike detection from inaccurate samplings[END_REF]). In our case, it seems that the near region analysis is closely linked to the coercivity constant of the second order shape derivative at optimal sets. What seems more challenging is to quantify how the relevant geometric functional behaves in the far region, i.e. for sets that are not smooth normal deformations of a minimizer. We argue that this could be related to the derivation of a quantitative inequality for our functional (we brie y discussed the case of the isoperimetric inequality in Section 2.1.1). To prove such a result, one would need to adapt the selection principle of [START_REF] Cicalese | A Selection Principle for the Sharp Quantitative Isoperimetric Inequality[END_REF], and then use the analysis provided in Section 2.2 to deal with smooth deformations of a minimizer.

Stability for the total variation ow. We argue that the stability analysis of Sections 2.2 and 3.1 could possibly be used to study the L 2 gradient ow of the total variation (see for instance [START_REF] Bellettini | The Total Variation Flow in RN[END_REF]). A relevant question concerning this evolution could be: if the ow is initialized with an M -simple function, is the solution also M -simple for su ciently small times? If so, is the length of the chains preserved? In the particular case where the ow is initialized with the indicator function of a calibrable set, it has been shown in the above reference that the solution remains proportional to the initialization. Answering the more general question formulated above is, to the best of our knowledge, an open problem. Given an initial data u 0 ∈ L 2 (R 2 ), an idea to apply the results of Sections 2.2 and 3.1 would be to consider the minimizing movement scheme

min u∈L 2 (R 2 ) TV(u) + 1 2τ u -u 0 L 2 , (58)
and to prove a stability result in the spirit of Theorem 2.27 when τ → 0 + . The main di erence with our setting is that, in (58), no advantage can be taken from the regularity of the measurement operator, and whether dual solutions converge smoothly to the minimal norm limit certi cate is unclear. Finally, let us mention that, although it is even less clear whether our results could be applied in this setting, one could also consider the same question for the Wasserstein gradient ow of the total variation [START_REF] Carlier | On the total variation Wasserstein gradient ow and the TV-JKO scheme[END_REF].

Model bias. In Part 2, we argued that the class of sparse objects associated to the total variation are M -simple functions. Regularizing an inverse problem with the total variation will hence often produce M -simple reconstructions. In practical applications, the sought-after image might be approximately, but not exactly piecewise constant. An interesting question would hence be to investigate what class of functions can be well approximated with s-sparse M -simple functions, and how this approximation behaves when s grows. This is related to the notion of model bias, i.e. to the error one makes by assuming the unknown signal has a speci c structure. Such an analysis has been conducted for sparse vectors (see for instance [START_REF] Foucart | A Mathematical Introduction to Compressive Sensing[END_REF], Section 2.1]). A related topic is also the study of nonlinear wavelet approximation (see [Mallat, 2008, Section 9.2]), and more speci cally of how the approximation error decreases when the number of coe cients grows.

N

In Part 3, we proposed an algorithm for solving (P λ (y)), which is based on the conditional gradient (or Frank Wolfe) algorithm. It iteratively constructs a sequence of iterates that provably converges to a solution of (P λ (y)). Each iterate is a linear combination of the atoms promoted by the total variation (i.e. indicator functions of simple sets). Then, we discussed the implementation of this theoretical algorithm, and proposed to numerically represent the above-mentioned atoms by indicator functions of simple polygons. We nally assessed the performance of our approach on a few recovery examples, providing comparisons with traditional grid-based methods.

Convergence in nite time. In [Denoyelle et al., 2019, Theorem 3], a nite time convergence result is proved for the so-called sliding Frank-Wolfe algorithm. It holds under the assumption that the dual certi cate associated to the noisy problem is non-degenerate. A crucial ingredient of the proof is that a critical point of the objective is reached at the end of each sliding step. As explained above, we cannot guarantee this in our setting. Still, assuming this is true, it would be interesting to investigate whether the proof of [START_REF] Denoyelle | The Sliding Frank-Wolfe Algorithm and its Application to Super-Resolution Microscopy[END_REF] can be adapted to Algorithm 3.

Topology changes. In Section 1.4, we presented the so-called sliding step (in the spirit of [START_REF] Denoyelle | The Sliding Frank-Wolfe Algorithm and its Application to Super-Resolution Microscopy[END_REF]). We formally described the associated gradient ow, but its study (existence, uniqueness) seems challenging. In Section 3.1 we introduced a possible way to mimic this evolution numerically. We explained that this could lead to the apparition of singularities, and that we could not handle them consistently with our approach. Even if one can circumvent this issue by stopping the evolution before the appearance of a singularity (which does not break any convergence guarantee), it seems preferable to allow them, in order to yield iterates that are as sparse as possible. This topic is related to the more global issue of handling topology changes with purely Lagrangian methods. In the shape optimization community, Eulerian methods (such as the level set method) are preferred to achieve this, as they can treat complicated deformations of a shape in a very robust way [START_REF] Allaire | Chapter 1 -Shape and topology optimization[END_REF]. In a recent work [Lévy, 2022], Bruno Lévy proposed a new Lagrangian free-surface mesh representation for uid simulation. He also mentioned these ideas could possibly be used for shape optimization. It would be interesting to investigate whether these tools can allow for a robust implementation of the sliding step.

Approximation of maximal elements of exposed faces. We mentioned above that the construction of pre-certi cates is highly important to derive identi ability and noise robustness guarantees. Given a pre-certi cate, one can numerically ensure it is a true certi cate by solving the associated Cheeger problem, and showing its value is (approximately) one. If this is the case, it is desirable to nd which functions are associated to this certi cate, i.e. to obtain a complete description of the face it exposes. One can always nd an element of this face by solving the associated Cheeger problem, but how to approximate the complete set of maximizers seems challenging. We argue that a possible way to achieve this would be to adapt the ideas of [START_REF] Buttazzo | On the selection of maximal Cheeger sets[END_REF], Carlier et al., 2009]. These works propose to introduce a small strictly convex penalization of the objective, in order to select some "maximal" solution of the Cheeger problem. Their setting however slightly di ers from ours (the weight functions they deal with are non-negative), and nding out whether their techniques can be adapted to our case would be interesting.

Practical applications. To conclude this part, let us discuss how the work presented here can be relevant for practical applications. The experiments conducted in Section 3.2 suggest that traditional grid-based methods perform better on complex natural images, while our approach is particularly well-suited for reconstructing simple piecewise constant ones. It is hence natural to wonder in which situations one would be interested in the recovery of such signals, and in obtaining a continuous domain representation of them. We mentioned in Section 1.1.1 that this could potentially be the case in cell imaging (see Figure 2). However, it seems natural that the choice of appropriate measurement operators and the relevance of this approach should extensively be discussed with practitioners. A. S

In this section, we collect a few useful de nitions and properties related to sets of nite perimeter. We refer the reader to e.g. [Maggi, 2012] for further details, and to [START_REF] Ambrosio | Connected components of sets of nite perimeter and applications to image processing[END_REF] for more information concerning M -connected components and related notions. In particular, it does not have isolated points. In the following, we always consider such a representative and consequently obtain Supp(D1 E ) = ∂ * E = ∂E .

Distributional curvature. If E ⊂ R 2 is a set of nite perimeter, then the distributional curvature vector of E is

H E : C ∞ c (R 2 , R 2 ) → R de ned by ∀T ∈ C ∞ c (R 2 , R 2 ), H E , T = ˆ∂ * E div E T dH 1 ,
where div E T denotes the tangential divergence of T on E given by

div E T = div T -(DT ν E ) • ν E ,
and DT denotes the di erential of T . A set E is said to have locally integrable distributional curvature if there exists a function H E ∈ L 1 loc (∂ * E; H 1 ) such that

H E = H E ν E H 1 ∂ * E .
For instance, if E is an open set with C 2 boundary, it has a locally summable distributional curvature which is given by the (classical) scalar mean curvature.

M -connected components, holes, exterior. If E ⊂ R 2 has nite perimeter it can be decomposed (up to Lebesgue negligible sets) into an at most countable union of pairwise disjoint indecomposable sets of positive measure, i.e.

B. S :

In this section, we give some de nitions and facts about smooth sets. Although these notions are standard (most of them are for example discussed in [START_REF] Delfour | Shapes and Geometries: Metrics, Analysis, Di erential Calculus, and Optimization[END_REF]), we recall them here to keep the exposition (mostly) self-contained.

Smooth sets. We say that a set is smooth if it is locally the hypograph of a smooth function.

If r > 0, x ∈ R 2 and ν ∈ S 1 , we denote by C(x 0 , r, ν) the square with axis ν and side r centered at x 0 , de ned as follows:

C(x 0 , r, ν)

def. = x 0 + x ∈ R 2 | x, ν | < r, |x -x, ν ν| < r . ( 60
)
With this notation, we have C(0, r, e 2 ) = (-r, r) 2 . We also denote by R ν the rotation which maps e 2 to ν (and hence (-r, r) 2 to C(0, r, ν)). Given a mapping u : (-r, r) → R we de ne graph(u)

def.

= {(z, u(z)), z ∈ (-r, r)} , hypograph(u)

def.

= {(z, t) | z ∈ (-r, r), -r < t < u(z)} .

De nition B.1

Let k ≥ 1. A set E ⊂ R 2 is said to be of class C k if for every x ∈ ∂E there exists r > 0, ν ∈ S 1 and u ∈ C k ([-r, r]) such that R -1 ν (∂Ex) ∩ (-r, r) 2 = graph(u) , R -1 ν (int Ex) ∩ (-r, r) 2 = hypograph(u) .

As stated in [Delfour and Zolesio, 2011, Theorem 5.2], if E is bounded, the compactness of ∂E allows to take r in De nition B.1 independent of x, and the family (u x ) x∈∂E uniformly bounded in C k with equicontinuous k-th derivative.

Convergence of smooth sets. Now, we turn to the de nition of C k convergence for smooth sets.

De nition B.2

Let E be a bounded set of class C k with k ≥ 1. We say that (E n ) n≥0 converges to E in C k if there exists r > 0 and n 0 ∈ N such that

• for every n ≥ n 0 we have ∂E n ⊂ x∈∂E C(x, r, ν E (x))

• for every n ≥ n 0 and x ∈ ∂E there exists u n,x ∈ C k ([-r, r]) such that: x) (int E nx) ∩ (-r, r) 2 = hypograph(u n,x )

R -1 ν E (x) (∂E n -x) ∩ (-r, r) 2 = graph(u n,x ) R -1 ν E (
• denoting (u x ) x∈∂E functions satisfying

R -1 ν E (x) (∂E -x) ∩ (-r, r) 2 = graph(u x ) R -1 ν E (x) (int E -x) ∩ (-r, r) 2 = hypograph(u x )
we have lim In Section 2.2, we make use of the following result, which states that, if a sequence of sets converges in C k , the boundary of its terms can eventually be written as a normal deformation of the boundary of the limit set. A related statement can be found for instance in [Acerbi et al., 2013, Theorem 4.2]. The authors of this last work refer to [Almgren, 1975], which, according to them, essentially contains the result. We provide a detailed proof below for the sake of completeness.

Proposition B.3

If (E n ) n≥0 converges to a bounded set E in C k with k ≥ 2, then for n large enough there exists a mapping ϕ n ∈ C k-1 (∂E) such that ∂E n = (Id + ϕ n ν E )(∂E), and ϕ n C k-1 (∂E) → 0.

Proof : In all the following, we denote by ν (instead of ν E ) the outward unit normal to E. If u ∈ C k (Ω) and δ > 0, we also use the notation

B C k (u, δ) def. = v ∈ C k (Ω) u -v C k (Ω) < δ .
We divide the proof in two main steps.

Step 1: we rst prove that, for every x ∈ ∂E, there exists n 0 ∈ N and two open neighborhoods U, V of x such that, for every n ≥ n 0 , there exists ϕ n : ∂E ∩ U → R of class C k-1 with ∂E n ∩ V = (Id + ϕ n ν)(∂E ∩ U ) and ϕ n C k-1 (∂E∩U ) → 0.

By the C k convergence of (E n ) n≥0 towards E, we can x x ∈ ∂E, r > 0 and n 0 ∈ N such that, in C(x, r, ν E (x)), the set E is the hypograph of u ∈ C k ([-r, r]) and E n is the hypograph of u n ∈ C k ([-r, r]) for every n ≥ n 0 , with u nu C k ([-r,r]) → 0. By an abuse of notation, the normal ν is given in local coordinates by ∀z ∈ (-r, r), ν(z, u(z)) def.

= 1 1 + u (z) 2 (-u (z), 1) .

Finally, de ning U def.

= (-δ 2 , δ 2 ) × (-r, r), U def.

= x + R ν E (x) U , V def.

= x + R ν E (x) V , ϕ v : y → ψ v ((R -1 ν E (x) (yx)) 1 ) , and taking n 0 large enough to have u n ∈ B C k (u, δ 1 ) for every n ≥ n 0 , we obtain

∂E n ∩ V = (Id + ϕ n ν)(∂E ∩ U )
and ϕ n C k-1 (∂E∩U ) → 0 with ϕ n def.

= ϕ un .

Step 2: from now on, we make the dependencies on x explicit. We know from step 1 that for every x ∈ ∂E there exists n 0 (x) ∈ N, a mapping ϕ n,x : ∂E ∩ U x → R, and two open neighborhoods U x , V x of x such that, for every n ≥ n 0 (x)

∂E n ∩ V x = (Id + ϕ n,x ν)(∂E ∩ U x ) .
Moreover, ϕ n,x C k-1 (∂E∩Ux) → 0.

Let us rst show that if y ∈ ∂E ∩ U x1 ∩ U x2 and n ≥ max(n 0 (x 1 ), n 0 (x 2 )) then Now, we cover ∂E by x∈∂E U x ∩ V x and obtain the existence of a nite set I such that E ⊂ i∈I U xi ∩ V xi . For every n ≥ max i∈I n 0 (x i ) we de ne ϕ n globally using the results above, and obtain

ϕ n C k-1 (∂E) ≤ max i∈I ϕ n,xi C k-1 (∂E∩Ux i ) → 0 .
By construction, we also have (Id + ϕ n ν)(∂E) = ∂E n , which concludes the proof.

Normal deformations of smooth sets. Considering the result of Proposition B.3, one could wonder if, given a su ciently smooth ϕ : ∂E → R in a neighborhood of 0, there is a unique smooth set F satisfying ∂F = (Id + ϕ ν E )(∂E) .

This is indeed the case, as shown below.

Lemma B.4

If E is a bounded set of class C k (k ≥ 2), then there exists C > 0 such that, for every ϕ in C k-1 (∂E), the mapping ϕ ν E can be extended to ξ ϕ ∈ C k-1 (R 2 , R 2 ) with

ξ ϕ C k-1 (R 2 ,R 2 ) ≤ C ϕ C k-1 (∂E)
.

Proof : We use the arguments of [Dambrine and Lamboley, 2019, Section 4]. Since E is bounded and of class C k , there exists an open bounded set Ω ⊃ ∂E of class C 1 , such that the projection π ∂E on ∂E belongs to C k-1 (Ω). For every ϕ ∈ C k-1 (∂E), one can then de ne ξ ϕ = (ϕ ν E ) • π ∂E to obtain an extension of ϕ ν E which belongs to C k-1 (Ω). Using Faà di Bruno's formula, one obtains the existence of C > 0 such that

∀ϕ ∈ C k-1 (∂E), ξ ϕ C k-1 (Ω,R 2 ) ≤ C ϕ C k-1 (∂E) .
One can then use an extension theorem (see for instance [Stein, 1970, Chapter 6, Theorem 5]) to obtain the result with Ω = R 2 , up to a modi cation of C.

Moreover, there exists an extension ξ ϕ of ϕ ν E such that E ϕ = (Id + ξ ϕ )(E) and

ξ ϕ C k-1 (R 2 ,R 2 ) < 1 .
In particular, E ϕ is C k-1 -di eomorphic to E.

Proof : Let us rst prove the existence of E ϕ . Let C be as in Lemma B.4. If c < 1/C, then for every ϕ such that ϕ C k-1 (∂E) ≤ c, we have ξ ϕ C k-1 (R 2 , R 2 ) < 1. As a result, we obtain that f ϕ def.

= Id + ξ ϕ is a C k-1 di eomorphism. We hence have that E ϕ If E and c are as in Proposition B.5, for every ≤ c, we say that a set F is a C k-1 -normal deformation of E of size at most if there exists ϕ such that ϕ C k-1 (∂E) ≤ and F = E ϕ . Likewise, if E is such that E c satis es the assumptions of Proposition B.5, we use the same terminology for F if F c = (E c ) ϕ .

C. C

The aim of this section is to prove the following result, which, loosely speaking, states that sets of nite perimeter and nite measure are characterized by their measure-theoretic boundary. In all the following d ≥ 2 is a xed integer and H s denotes the s-dimensional Hausdor measure on R d . Although this result is certainly well-known, we were not able to nd it clearly stated or used, and hence decided to include a proof for the sake of completeness. To prove it, we rely on a tool introduced in [Ambrosio et al., 2001, Section 7]: the topographic function (see Figure 23 for an illustration). To be more precise, if E and F are as in Proposition C.6 and u and v denote their topographic function, we prove that u = v almost everywhere, which, since 1 E = u mod 2 and 1 F = v mod 2, yields the result. This is done by induction, i.e. by proving that {u = k} = {v = k} (mod H d )

for every k ∈ N. Before diving into the proof, we state a few useful lemmas.

Lemma C.7

Let E, F ⊆ R d be of nite perimeter. If E is indecomposable and H d-1 (∂ M F ∩ EM ) = 0 then we have E ⊆ F (mod H d ) or E ⊆ F c (mod H d ).

Proof : Since supp(D1 F ) = ∂ M F (mod H d-1 ) we have |D1 F |( EM ) = 0. The indecomposability of E hence yields, using a result of [START_REF] Dolzmann | Microstructures with nite surface energy: The two-well problem[END_REF], that 1 F is a.e. equal to 0 on E, or a.e. equal to 1 on E, hence the result.

In order to initialize our proof by induction, we need to prove the equality of the zero level set of the topographic functions. From the construction of [Ambrosio et al., 2001, Theorem 6] it is clear that, if u is the topographic function of a set of nite perimeter E with nite measure, then ext(E) = {u = 0}. We can hence conclude using the following lemma. The rst term converges to 0. The second converges to k assuming x ∈ {u ≤ k} (1/2) (which is not restrictive since H d-1 (∂ M {u ≤ k} \ {u ≤ k} (1/2) ) = 0). To show that the last term converges to 0, we write The rst term converges to 0 assuming x ∈ ∂ * {u ≤ k}. Again, this is not restrictive since

H d-1 (∂ M {u ≤ k} \ ∂ * {u ≤ k}) = 0 .
By [Ambrosio et al., 2000, Lemma 3.75] the second term is bounded for H d-1 -a.e. Proof : We already know that {u = 0} = {v = 0}. Let us prove by induction that {u = k} = {v = k} for every k ∈ N. We x k ≥ 1 and assume that {u = l} = {v = l} for every l ≤ k.

Let A ∈ CC M ({u = k + 1}). Since A is indecomposable, u is constant on A and

|Du| = H d-1 ∂ M E = H d-1 ∂ M F = |Dv| ,
we obtain that v is a.e. constant on A. The construction of [Ambrosio et al., 2001, Theorem 6] also yields

H d-1 (∂ M {v ≤ k} ∩ ∂ M A) = H d-1 (∂ M {u ≤ k} ∩ ∂ M A) > 0 .
Using Lemma C.9, we obtain that v = k + 1 a.e. on A, and hence A ⊆ {v = k + 1} (mod H d ).

Since this holds for every A ∈ CC M ({u = k + 1)}) we get {u = k + 1} ⊆ {v = k + 1} (mod H d ).

Exchanging E and F yields {v = k + 1} ⊆ {u = k + 1} (mod H d ), and we can nally conclude.

D. C

The aim of this short section is to prove the following elementary result. (t) . This implies the pointwise convergence of f almost everywhere, but only up to the extraction of a subsequence. Lifting this requirement is in fact possible, as shown below.

Proof : Let t ≥ 0. We rst prove that |{u n ≥ t} \ {u ≥ t}| → 0. For every h > 0 we have: for almost every t ≤ 0 follows from the same arguments.

u n -u L 1 ≥ ˆt t

E. P L 3.2

We now give the proof of Lemma 3.2, whose content is recalled below. = {u ∈ L 2 (R 2 ) | TV(u) ≤ t ≤ t * } are:

• (0, 0) ,

• ( t * 1 E /P (E), t * ) with ∈ {-1, 1}, E ⊂ R 2 simple and 0 < |E| < +∞ .

Proof : Let (u, t) ∈ ext(C).

1. If t = 0 then u = 0.

2. Otherwise we have (u, t) = (1t/t * )(0, 0) + (t/t * )(t * u/t, t * ) , which yields t = t * . Considering that, for any v such that TV(v) ≤ t * , we have:

(0, t * ) = 1 2 (-v, t * ) + 1 2 (v, t * ) ,
we also obtain TV(u) > 0. Now, since (u, t * ) = (1 -TV(u)/t * )(0, t * ) + (TV(u)/t * )(t * u/TV(u), t * ) , we get TV(u) = t * . Finally, let α ∈ (0, 1) and u 1 , u 2 ∈ L 2 (R 2 ) with TV(u 1 ) ≤ t * and TV(u 2 ) ≤ t * . If u = (1α)u 1 + αu 2 then writing (u, t * ) = (1α)(u 1 , t * ) + α(u 2 , t * ) we obtain u = u 1 = u 2 , which shows u ∈ ext({u ∈ L 2 (R 2 ) | TV(u) ≤ t * }). This shows that u = t * 1 E /P (E) for some simple set E such that 0 < |E| < +∞ and ∈ {-1, 1}.

3B

  Exact support recovery 3.1 Stability of the dimension of the faces of the total variation unit ball . . . . . . . 3.2 Main result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 A su cient condition for strict stability . . . . . . . . . . . . . . . . . . . . . . . 3.4 Computing the minimal norm certi cate . . . . . . . . . . . . . . . . . . . . . . -Wolfe algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Proposed algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Convergence results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Sliding step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3 Numerical results 3.1 Local minimization of the objective . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Recovery examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Topology changes during the sliding step . . . . . . . . . . . . . . . . . . . . . . Smooth sets: de nition and convergence C Characterization of sets from their boundary D Convergence of the level sets E Proof of Lemma 3.2

N

  Functions spaces. Let m, n ∈ N * and Ω ⊂ R m be an open set. For any k ∈ N, we denote byC k (Ω, R n ) (C k (Ω) if n = 1) the set of R n -valued functions of class C k on Ω, and by C k b

  problems in imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Variational regularization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Grid-free sparse recovery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Total variation regularization 2.1 Presentation and preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Noise robustness results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 -

 2 Figure 1 -Examples of microtubule images.

Figure 5 - 2 N

 52 Figure 5 -Comparison of grid-based and grid-free approaches for representing simple images.

C1

  Piecewise constant functions 1.1 Exposed faces of the total variation unit ball . . . . . . . . . . . . . . . . . . . . 1.2 Chains associated to a piecewise constant function . . . . . . . . . . . . . . . . 2 The prescribed curvature problem 2.1 Generalities and rst convergence result . . . . . . . . . . . . . . . . . . . . . . 2.2 Stability result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Exact support recovery 3.1 Stability of the dimension of the faces of the total variation unit ball . . . . . . . 3.2 Main result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 A su cient condition for strict stability . . . . . . . . . . . . . . . . . . . . . . . 3.4 Computing the minimal norm certi cate . . . . . . . . . . . . . . . . . . . . . .

Figure 8 -

 8 Figure 8 -Illustration of the result stated in Theorem 2.27. Here u 0 is equal to 0 in A and B. The values taken by u λ,w in A λ,w and B λ,w are close (but not necessarily equal) to 0.

Figure 9 -

 9 Figure 9 -Graphs of R 0 → α(R 0 ) for the Gaussian and ideal low-pass lters, for several values of the variance and cut-o frequency.

Figure 10 -

 10 Figure 10 -Plots of Φ * f E and Φ * g E de ned in (27) for E = 1 B(0,1) and Φ the convolution with the Gaussian kernel with variance σ = 0.2.

Figure 11 -

 11 Figure 11 -Graph of f v de ned in (29) (left: global graph, right: zoom around 1).

Figure 12 -

 12 Figure 12 -Graph of ∂ηv ∂r (1) as a function of σ.

Figure 14 -

 14 Figure 14 -Steiner symmetrization of triangles

  Figure 15 -Steiner symmetrization of quadrilaterals

Figure 16 -

 16 Figure 16 -From left to right: observations, unknown function, output of Algorithm 3, outputs of the xed grid method using the isotropic and Condat's total variation

Figure 19 -

 19 Figure 19 -Left: unknown function, middle: observations, right: output of Algorithm 3 for di erent values of λ

Figure 21 -

 21 Figure 21 -Topology change experiment. (a): unknown signal, (b): observations, (c): weighted Cheeger set, (d,e,f,g): sliding step iterations (with splitting), (h): nal function.

AA

  Sets of nite perimeter B Smooth sets: de nition and convergence C Characterization of sets from their boundary D Convergence of the level sets E Proof of Lemma 3.2

  Reduced boundary. The reduced boundary ∂ * E of a set of nite perimeter E is de ned as the set of points x ∈ Supp (|D1 E |) at whichν E (x) B(x, r)) |D1 E | (B(x, r))exists and satis es |ν E (x)| = 1.Choice of representative. From[Giusti, 1984, Proposition 3.1], we know that if E has nite perimeter, there exists a Lebesgue representative of E with the property that ∀x ∈ ∂E, 0 < |E ∩ B(x, r)| < |B(x, r)| .

Figure 22 -

 22 Figure 22 -Illustration of De nition B.1.

  xu x C k ([-r,r]) = 0

  ϕ n,x1 (y) = ϕ n,x2 (y) .We assume without loss of generality that|ϕ n,x1 (y)| ≤ |ϕ n,x2 (y)|. From step 1, we know that ϕ n,x2 (y) is the unique t ∈ (-γ x2 , γ x2 ) such that y + t ν(y) ∈ ∂E n ∩ C(x 2 , δ 2 , ν(x 2 )).But by de nition y + ϕ n,x1 (y) ν(y) ∈ ∂E n , and since |ϕ n,x1 (y)| ≤ |ϕ n,x2 (y)| we also have that ϕ n,x1 (y) ∈ (-γ x2 , γ x2 ) and y + ϕ n,x1 (y) ν(y) ∈ C(x 2 , δ 2 , ν(x 2 )). We therefore obtain that ϕ n,x1 (y) = ϕ n,x2 (y).

Proposition B. 5

 5 Let E be a bounded open set of class C k (with k ≥ 2). There exists c > 0 such that, for every ϕ ∈ C k-1 (∂E) with ϕ C k-1 (∂E) ≤ c, there is a unique bounded open set of class C k-1 , denoted E ϕ , satisfying ∂E ϕ = (Id + ϕ ν E )(∂E) .

=

  f ϕ (E) is a bounded open set of class C k-1 and satis es ∂E ϕ = f ϕ (∂E) = (Id + ϕ ν E )(∂E).

  Let us now turn to uniqueness. If F and G are two bounded open sets of class C k-1 satisfying (63), in particular we have that ∂F = ∂G. Since F and G are su ciently smooth and bounded, they have nite perimeter and ∂ M F = ∂F = ∂G = ∂ M G. Using Proposition C.6 below (whose proof is the object of Appendix C), we obtain that F = G up to a Lebesgue negligible set. Using that F and G are su ciently smooth open sets, this nally yields F = G.

Figure 23 -

 23 Figure 23 -The topographic function of a set E (gray area).

  with B ± (x, r, ν(x)) = {y ∈ B(x, r) | ± yx, ν > 0}. Since |Du| = H d-1 ∂ M E we also have |a -b| = 1. Let us show that a ≤ k and b ≥ k + 1. We have: B -(x,r,ν(x)) -1 {u≤k} ) ≤ B -(x,r,ν(x)) |u -a| + 2 k |{u ≤ k} ∩ B(x, r)| |B(x, r)| + B(x,r)u (1 B -(x,r,ν(x)) -1 {u≤k} ) .

B

  (x,r) u (1 B -(x,r,ν(x)) -1 {u≤k} ) ≤ |B -(x, r, ν(x)) \ {u ≤ k}| 1 d |B(x, r)| ˆB(x,r) |u| d/(d-1) d-1 d = |B -(x, r, ν(x)) \ {u ≤ k}| |B(x, r)|

  x. We hence obtain a ≤ k. By using the same arguments, one can show b ≥ k + 1 and we nally get a = k and b = k + 1.Proposition C.10Let E, F ⊂ R d be sets of nite perimeter with nite measure and u, v their topographic functions.If ∂ M E = ∂ M F (mod H d-1) then u = v almost everywhere.

Lemma 3. 2

 2 If t * > 0, the extreme points of C def.

  A function u has nite total variation if and only if its distributional gradient, denoted Du, is a nite Radon measure. In this case one has TV(u) = |Du|(R 2 ). Let us now provide more intuition on this functional by considering two classical examples. If u has integrable gradient one has TV(u) = ´R2 |∇u|. If E is a set of class C 1 (see Appendix B for a precise de nition) by the Gauss-Green theorem ´E

  and hence get the inequality for all k. Using that |U Figure 13 -The sets E and F have the weak regularity property (point 4 of Lemma 1.11), but E ∪ F does not, since the density ratio of x with respect to E ∪ F converges to 1. satis es point 4 of Lemma 1.11 (a property that one would crucially need to mimic the proof of Corollary 3.8). The level sets U

		E	x	F
	Convergence of the boundaries of the level sets ? Considering Corollary 3.8, it is natural
	to wonder whether	lim sup n→+∞	∂U k ⊆ ∂U (t) * (t)
	holds, which would yield the convergence of (∂U k ) k≥0 towards ∂U (t) * . However, it is unclear (t) whether U (t)
			(t)
	get a contradiction.			(t) k	U * | → 0, we (t)

k

  a set of nite perimeter such that 0 < |E| < +∞ and ν ∈ S 1 we have:

	´E η P (E)	≤	´Es ν P (E s η ν )	,
	with equality if and only if |E E s ν | = 0. Hence if E * solves (41), arguing as in

  Proposition D.1If (Ω, Σ, |.|) is a measure space and u nu L 1 (Ω) → 0 then for almost every t ∈ R we have:Before giving its proof, let us stress that this property is certainly well-known, but we could not nd it clearly stated or proved. In fact, using Fubini's theorem, one could directly obtain and hence prove the L 1 convergence of f : t → U

			lim n→+∞	U (t) n	U (t) → 0 .	(65)
	lim n→+∞	ˆ+∞ 0	U (t) n	U (t) dt = lim n→+∞ ˆR2	|u n -u| = 0 ,
						(t) n	U

  -h} \ {u ≥ t} = ∅, we obtain that the right-hand side in (66) converges to 0 when h goes to 0, and therefore that:lim n→+∞ |{u n ≥ t} \ {u ≥ t}| = 0 .Let us now prove that |{u ≥ t} \ {u n ≥ t}| → 0. For every h > 0 it holds that:u nu L 1 ≥ Using that |B \ A| ≤ |C \ A| + |B \ C| for every triple of sets (A, B, C) we get: |{u ≥ t} \ {u n ≥ t}| ≤ |{u ≥ t + h} \ {u n ≥ t}| + |{u ≥ t} \ {u ≥ t + h}| ≤ 1 h u nu L 1 + |{u ≥ t} \ {u ≥ t + h}| , t} \ {u n ≥ t}| ≤ |{u ≥ t} \ {u ≥ t + h}| ,(67)Since {u ≥ t} \ {u ≥ t + h} decreases as h decreases andh≥0 {u ≥ t} \ {u ≥ t + h} = {u = t},we obtain that, for almost every t ≥ 0, the right-hand side in (67) converges to 0 and therefore that:lim n→+∞ |{u ≥ t} \ {u n ≥ t}| = 0 .We hence nally have |{u n ≥ t} {u ≥ t}| → 0 for almost every t ≥ 0. The fact|{u n ≤ t} {u ≤ t}| → 0

	and hence
	lim sup
	which yields
	lim sup
		ˆt+h
		t ˆt+h	|{u n ≥ s} {u ≥ s}| ds
	≥	t ˆt+h	|{u ≥ s} \ {u n ≥ s}| ds
	≥

-h |{u n ≥ s} {u ≥ s}| ds ≥ ˆt t-h |{u n ≥ s} \ {u ≥ s}| ds ≥ ˆt t-h |{u n ≥ t} \ {u ≥ t -h}| ds = h |{u n ≥ t} \ {u ≥ t -h}| .

Using that |A \ B| ≤ |A \ C| + |C \ B| for every triple of sets (A, B, C) we get:

|{u n ≥ t} \ {u ≥ t}| ≤ |{u n ≥ t} \ {u ≥ t -h}| + |{u ≥ t -h} \ {u ≥ t}| .

We hence obtain:

|{u n ≥ t} \ {u ≥ t}| ≤ 1 h u nu L 1 + |{u ≥ t -h} \ {u ≥ t}| , n→+∞ |{u n ≥ t} \ {u ≥ t}| ≤ |{u ≥ t -h} \ {u ≥ t}| .

(66)

Since {u ≥ t -h} \ {u ≥ t} decreases as h decreases and h≥0 {u ≥ t t |{u ≥ t + h} \ {u n ≥ t}| ds = h |{u ≥ t + h} \ {u n ≥ t}| . n→+∞ |{u ≥

Our presentation is greatly inspired by[Peyré,

2021, Chapter 8]. We refer the reader to this book draft and the companion website for more details on imaging. Let us also mention the reference book[START_REF] Scherzer | Variational Methods in Imaging[END_REF].2 In many situations of interest, the assumption that Φ is linear is unrealistic. Still, the class of linear inverse problems encompass a large variety of practically relevant applications. We stress that the techniques used to study nonlinear inverse problems signi cantly di er from those we present in this manuscript.

This is the so-called synthesis approach. Another one is the analysis approach (see for instance[Elad et al.,

2007]and[März et al., 2022, Section 1.2.2]).

See Section

2.1.1 for a precise de nition of the total variation and simple sets.

The name BLASSO stands for Beurling LASSO, and was coined in[De Castro and Gamboa,

2012]. We provide more information on this problem in Section 1.3.2 below.

This inequality is for example proved in[Ambrosio et al.,

2000, Theorem 3.47] by using (2) and the coarea formula (see e.g.[Maggi, 2012, Chapter 13]). In this reference u is assumed to be integrable, but the constant m appearing therein is also zero for any p-integrable function with 1 ≤ p < +∞.

This notion is known under the name of inseparable set in the literature on submodular functions (see e.g.[Bach,

2013]).

This is for instance proved in[Chambolle et al.,

2016, proof of Proposition 8].

We stress that this set only depends on η0, and hence on y0, rather than on u0.

One could also wonder how to numerically solve (P0(y0)). As this is less relevant for practical applications, we choose to not cover this question here.

To see this, one simply needs to use Proposition 1.3 with Lemma 1.11.

As underlined by the authors, this consistency is however a bit weak, as it does not provide convergence rates or error bounds. This property is for example also satis ed by the isotropic discretization.

We also stress, as in[Duval,

2022, Section 2.6.2], that the total variation is a particular case of submodular function. The faces of the unit ball de ned by such functions have been studied in the monographs[Fujishige, 2005, Bach, 2013].

Meaning that |{u = ti}| > 0 for all i ∈ {1, ..., m} and, for almost every x ∈ R

, we have u(x) ∈ {ti} m i=1 .

This mapping allows to study the behaviour of the objective in a neighborhood of E with respect to general deformations, while jE is only related to normal deformations.

Here E k denotes the set E de ned in (16) associated to F k .

We recall that Line 11 amounts to solve a LASSO-type problem for which e cient solvers exist.

Indeed, (31) and (41) are equivalent in the sense that (E, ) solves (31) if and only if E solves (41) and is the opposite of the sign of ´E η.

To be more precise, we mean that there exists a compact set K ⊂ Ω such that ´R2 \K η

is su ciently small.

We stress that, although we do not follow this approach, nding some maximal solution of (42) and adding several atoms at each iteration could potentially be interesting. This would lead to what could be called a "polyatomic" Frank-Wolfe algorithm (see[Jarret et al.,

2022] and the related approach of[START_REF] Flinth | On the linear convergence rates of exchange and continuous methods for total variation minimization[END_REF]). 2 Indeed, if u solves (42), then there exists α such that αη ∈ ∂TV(u), and the result follows from Proposition 1.3 and Lemma 1.11

We already introduced this tool in Part

2, where we only considered normal deformations (in contrast with the general deformations we use here).

Here critical point is to be understood in the sense that the limit appearing in (53) is equal to zero for every θ.

In particular, this implies that η is essentially bounded.

This only occurs when λ is small enough. For higher values of λ, the output is similar to (d) or (e).

Existence of a maximizing simple n-gon. Let us now show there exists a maximizer which belongs to X n , i.e. that de nes a simple polygon. To do so, we rely on the following lemma. Loosely speaking, it states that starting from a non-simple polygon, one can always construct another polygon with strictly fewer vertices achieving a higher objective value.

Lemma 3.12

Let m ≥ 3 and x ∈ Y m \ X m . Then there exists m with 2 ≤ m < m and y ∈ Y m such that J(x) ≤ J(y) .

Proof :

x 1 ] is not simple. If there exists i with x i = x i+1 then y = (x 1 , ..., x i , x i+2 , ..., x m ) is suitable, and likewise if x 1 = x m then y = (x 1 , ..., x m-1 )

is suitable. Otherwise we distinguish the following cases:

If there exists i < j with x i = x j : we de ne y = (x 1 , ..., x i , x j+1 , ..., x m ) ∈ R m-(j-i) , z = (x i , x i+1 , ..., x j-1 ) ∈ R j-i .

We notice that 2 ≤ ji < m and 2 ≤ m -(ji) < m.

If there exists i < j with x i ∈]x j , x j+1 [: we necessarily have (i, j) = (1, m). We de ne y = (x 1 , ..., x i , x j+1 , ..., x m ) ∈ R m-(j-i) , z = (x i , x i+1 , ..., x j ) ∈ R j-i+1 .

We again have 2 ≤ m -(ji) < m, and since (i, j) = (1, m), we have ji < m -1 which shows that 2 ≤ ji + 1 < m.

If there exists i < j with x j ∈]x i , x i+1 [: we necessarily have j > i + 1. We de ne

We again have 2 ≤ ji < m, and since

If there exists i < j with

[ and in both cases we fall back on the previously treated cases. The same holds if (i, j) = (1, m). Otherwise, we de ne

u [1] (before sliding)

Supp(Du [1] ), Supp(Du 0 )

u [2] (before sliding)

Supp(Du [2] ), Supp(Du 0 )

u [3] (before sliding)

Supp(Du [3] ), Supp(Du 0 )

Figure 17 -Unfolding of Algorithm 3 for the rst experiment (u [k] denotes the k-th iterate)

Figure 18 -From left to right: unknown function, observations, outputs of the xed grid method using the isotropic and Condat's total variation, output of Algorithm 3, gradients support (red: output of Algorithm 3, black: unknown) P 4

C

We end this manuscript by presenting natural avenues for future research. We rst discuss those related to theoretical recovery guarantees, before turning to numerical methods.

R

Let us start by recalling our contributions to the theoretical analysis of total variation regulatization. In Part 2, we described the structure of the exposed faces of the total variation unit ball, and introduced the class of M -simple functions, which are the sparse objects naturally associated to this regularizer. We also considered solutions of the prescribed curvature problem, and studied their stability under variations of the associated curvature functional. Finally, we introduced a condition under which the jump set of an M -simple function can be exactly recovered in a low noise regime.

Su cient identi ability conditions. In our opinion, an important question which remains mostly open is that of identi ability. Given a function u 0 and a measurement operator Φ, being able to provide su cient conditions ensuring u 0 is the unique solution of (P 0 (y 0 )) is highly desirable. As already mentioned, the only related work we are aware of is [START_REF] Bredies | A perfect reconstruction property for PDE-constrained total-variation minimization with application in Quantitative Susceptibility Mapping[END_REF].

A major di erence with the sparse spikes setting is that even the identi ability of a single atom, i.e. of u 0 = 1 E with E simple, seems di cult to study. Having a clear understanding of this problem could, subsequently, allow to answer the following question: if two atoms E, F are separately identi able and a, b ∈ R * , how does the relative position of E and F impact the identi ability of a 1 E + b 1 F ?

Pre-certi cates. Considering the literature on sparse spikes recovery, the construction of pre-certi cates (i.e. good candidates for being a dual certi cate, possibly of minimial norm) is highly important for studying both identi ability and noise robustness. We have discussed in Section 3.4.2 possible ways to construct such objects in a simpli ed setting. This investigation should however be furthered. In particular, the pre-certi cate we described only uses the rst order optimality condition "in average", and other constraints could hence be considered. To

We notice that if (z, t) ∈ (-r, r)×R and v : (-r, r) → R, then (z, u(z))+t ν(z, u(z)) ∈ graph(v) if and only if

= r/2, then (61) holds and f (v, z, t) is hence well de ned. We claim that f is of class

To see this, one should rst notice that the two mappings

. Moreover, we have f (u, z, 0) = 0 for every z ∈ (-r , r ), and in particular f (u, 0, 0) = 0. Using that u (0) = 0, we also get ∂ t f (u, 0, 0) = 1 = 0. From the implicit function theorem, we hence obtain the existence of δ 1 > 0, of δ 2 ∈ (0, r ), of an open neighborhood W ⊂ (-r , r ) of 0, and of a mapping ψ :

Since f (u, z, 0) = 0, we obtain ψ(u, z) = 0 for every z ∈ (-δ 2 , δ 2 ). Now, using the fact ψ is of class C k-1 , we obtain that all its derivatives up to order k -1 are uniformly continuous on every compact subset of B C k-1 (u, δ 1 ) × (-δ 2 , δ 2 ). By the Ascoli-Arzelá theorem C k ([-r, r]) is compactly embedded into C k-1 ([-r, r]), and we hence have that, for every δ 2 < δ 2 , the family of functions

Now, we de ne

We have that F is of class C k-1 and DF (0, 0) = Id. There hence exists two open neighborhoods

We take δ 2 < δ 2 small enough to have (-δ 2 , δ 2 ) ⊂ U 1 , and

By construction, the set on the left-hand side is included in the one on the right-hand side. Now, if y ∈ graph(v) ∩ V , then denoting (z, t) def.

= F -1 (y), we have that z ∈ (-δ 2 , δ 2 ). Moreover, we also have (v, z, t) ∈ B C k-1 (u, δ 1 ) × (-δ 2 , δ 2 ) × W and f (v, z, t), which shows t = ψ v (z), and hence (62).

Lemma C.8

Assume E and F are two sets of nite measure and nite perimeter in R d . Then, if

we have ext(E) = ext(F ) (mod H d ).

Proof : We have:

which is necessarily ext(F ) since ext(E) has in nite measure. We hence obtain

Exchanging E and F and applying the same argument yields the result.

Before proving the main result, we need a last lemma, which is Lemma C.9. In the following, if u has bounded variation, we denote by J u its approximate jump set and S u its approximate discontinuity set (see [Ambrosio et al., 2000, Sections 3.6 and 3.7]).

Lemma C.9

Let E ⊆ R d be a set of nite perimeter, u its topographic function and k ∈ u(R 2 ). Then for H d-1 -almost every x ∈ ∂ M {u ≤ k} we have x ∈ J u and, denoting ν(x) the outward unit normal to {u ≤ k} at x, the following holds:

u and assume the approximate limit ũ(x) of u at x satis es ũ(x) ≤ k. We hence have:

Since the last term converges to 0 as r → 0 + , we obtain that x ∈ {u ≤ k} M . Likewise if ũ(x) > k we obtain that x ∈ ({u ≤ k} c ) M . This shows ∂ M {u ≤ k} ⊆ S u , and we nally obtain (64). Now, let x ∈ ∂ M {u ≤ k} ∩ J u (which, from our assumptions, is non-empty). We denote by ν the outward unit normal to {u ≤ k}, and have the existence of a = b such that

Conversely, let (u, t), (u 1 , t 1 ) and (u 2 , t 2 ) be in C and α ∈ (0, 1) with

1. If (u, t) = (0, 0), since t 1 , t 2 are nonnegative we obtain t 1 = t 2 = 0 and therefore u 1 = u 2 = 0. Hence (0, 0) ∈ ext(C).

2. If (u, t) = ( t * 1 E /P (E), t * ) with ∈ {-1, 1} and E a simple set such that 0 < |E| < +∞, then since t 1 ≤ t * and t 2 ≤ t * we get t 1 = t 2 = t * . Moreover, using that

we also obtain u 1 = u 2 = u. We therefore have ( t * 1 E /P (E), t * ) ∈ ext(C).

M

Problèmes inverses, variation totale, parcimonie.

R

On s'intéresse dans cette thèse à une famille de problèmes inverses, qui consistent à reconstruire une image à partir de mesures linéaires possiblement bruitées. On cherche à analyser les méthodes de reconstruction variationnelles utilisant un régulariseur spéci que, la variation totale (du gradient). Cette fonctionnelle est utilisée en imagerie depuis les travaux de Rudin Osher et Fatemi, menés en 1992

A

This thesis is devoted to the recovery of piecewise constant images from noisy linear measurements. We aim at analyzing variational reconstruction methods based on total (gradient) variation regularization. The total variation functional has been extensively used in imaging since the 90's. Its minimization is known to produce piecewise constant images, which hence have some kind of sparsity (they can be decomposed as the superposition of a few simple shapes). However, the performance of this regularizer has to our knowledge not extensively been studied from a sparse recovery viewpoint. This thesis aims at bridging this gap. We rst focus on noise robustness results. We assume that sought-after image is sparse, and study the structure of reconstructions in a low noise regime: are they sparse, made of the same number of shapes, and are these shapes close to those appearing in the unknown image? We then turn to numerical methods for total variation regularization. Existing techniques rely on the introduction of a xed spatial discretization, which often yield reconstruction artifacts such as anisotropy or blur. We propose an algorithm which does not su er from this grid bias, and produces a sparse representation of the reconstructed image.

K

Inverse problems, total variation, sparsity.