Swann Marx

Edouard Pauwels

Tillmann Weisser

Didier Henrion

Jean Bernard Lasserre

Semi-algebraic approximation using Christoffel-Darboux kernel

Keywords: approximation theory, convex optimization, moments, positive polynomials, orthogonal polynomials 42C05, 47B32, 41A30

We provide a new method to approximate a (possibly discontinuous) function using Christoffel-Darboux kernels. Our knowledge about the unknown multivariate function is in terms of finitely many moments of the Young measure supported on the graph of the function. Such an input is available when approximating weak (or measure-valued) solution of optimal control problems, entropy solutions to non-linear hyperbolic PDEs, or using numerical integration from finitely many evaluations of the function. While most of the existing methods construct a piecewise polynomial approximation, we construct a semi-algebraic approximation whose estimation and evaluation can be performed efficiently. An appealing feature of this method is that it deals with non-smoothness implicitly so that a single scheme can be used to treat smooth or non-smooth functions without any prior knowledge. On the theoretical side, we prove pointwise convergence almost everywhere as well as convergence in the Lebesgue one norm under broad assumptions. Using more restrictive assumptions, we obtain explicit convergence rates. We illustrate our approach on various examples from control and approximation. In particular we observe empirically that our method does not suffer from the the Gibbs phenomenon when approximating discontinuous functions.

Introduction

In this paper we address the following generic inverse problem. Let f : X → Y x := (x 1 , x 2 , . . . , x p-1 ) → y be a bounded measurable function from a given compact set X ⊂ R p-1 to a given compact set Y ⊂ R, with p ≥ 2. We assume that X is equal to the closure of its interior.

Given d ∈ N, consider a vector of polynomials of total degree at most d b(x, y) : (x, y) ∈ R p → (b 1 (x, y) b 2 (x, y)

• • • b n d (x, y)) ∈ R n d ,
where n d := p + d d which is understood as the binomial coefficient. For example, b may be a vector whose entries are the polynomials of the canonical monomial basis or any orthonormal polynomial basis, e.g. Chebyshev or Legendre. Associated to b and f , let

X b(x, f (x))b(x, f (x)) dx
be the moment matrix of degree 2d, where the integral is understood entry-wise.

Problem 1 (Graph inference from moment matrix). Given the moment matrix of degree 2d, compute an approximation f d of the function f , with convergence guarantees when degree d tends to infinity.

Motivation

Inverse Problem 1 is encountered in several interesting situations. In the weak (or measurevalued) formulation of Optimal Control problems (OCP) [START_REF] Vinter | Convex duality and nonlinear optimal control[END_REF][START_REF] Hernández-Hernández | The linear programming approach to deterministic optimal control problems[END_REF][START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF], Markov Decision Processes [START_REF] Hernández-Lerma | of Handbook of Markov Decision Processes[END_REF], option pricing in finance [START_REF] Lasserre | Pricing a class of exotic options via moments and SDP relaxations[END_REF], stochastic control and optimal stopping [START_REF] Stockbridge | Discussion of dynamic programming and linear programming approaches to stochastic control and optimal stopping in continuous time[END_REF], and some nonlinear partial differential equations (PDEs) [START_REF] Diperna | Measure-valued solutions to conservation laws[END_REF], non-linear non-convex problems are formulated as linear programming (LP) problems on occupation measures. Instead of the classical solution, the object of interest is a measure supported on the graph of the solution. Numerically, we optimize over finitely many moments of this measure. Following the notation introduced above for stating Problem 1, and letting z := (x, y) ∈ R p , the moment matrix of degree 2d reads

M µ,d := b(z)b(z) dµ(z)
corresponding to the measure dµ(z) := I X (x) dx δ f (x) (dy) [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF] supported on the graph {(x, f (x)) : x ∈ X} ⊂ X × Y of function f , where I X denotes the indicator function of X which takes value 1 on X and 0 otherwise, and δ f (x) denotes the Dirac measure at f (x). For instance, in OCPs an optimal occupation measure µ is supported on the graphs of optimal state-control trajectories. In order to recover a particular state resp. control trajectory it suffices to consider the moments of the marginal of the occupation measure µ with respect to time-state resp. time-control. Then with our notation, x is time and y is a state resp. control coordinate. Similarly, for the measure-valued formulation of non-linear first-order scalar hyperbolic PDEs, an occupation measure is supported on the graph of the unique optimal entropy solution. Then with our notation, x is time and space and y is the solution. From the knowledge of the moments of the occupation measure, we want to approximate the solution. The measure µ can be disintegrated into its marginal on X and its conditional on Y given x in X. The latter is also called a parametrized measure or a Young measure, see e.g. [START_REF] Fattorini | Infinite dimensional optimization and control theory[END_REF].

This weak formulation has been used in a number of different contexts to prove existence and sometimes uniqueness of solutions. It turns out that it can also be used for effective computation as it fits perfectly the LP-based methodology described in [START_REF] Helmes | Computing moments of the exit time distribution for Markov processes by Linear Programming[END_REF] and the Moment-SOS (polynomial sums of squares) methodology described in e.g. [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF][START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF][START_REF] Lasserre | The Moment-SOS Hierarchy[END_REF]. In the latter methodology one may thus approximate the optimal solution µ of the measure-valued formulation by solving a hierarchy of semidefinite relaxations of the problem, whose size increases with d; see e.g. [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF] for OCPs and [START_REF] Marx | A moment approach for entropy solutions to nonlinear hyperbolic PDEs[END_REF] for non-linear PDEs. An optimal solution of each semidefinite relaxation is a finite matrix of pseudo moments (of degree at most 2d) which approximate those of µ. This approach allows to approximate values for the corresponding variational problems but it does not provide any information about the underlying minimizing solutions beyond moments of measures supported on these solutions. Therefore an important and challenging practical issue consists of recovering from these moments an approximation of the trajectories of the OCP or PDEs. This is precisely an instance of Problem 1.

Another potential target application of our method is the optimal transportation problem, see e.g. [34, Chapter 1] and references therein. In its original formulation by Monge, it is a highly nonlinear nonconvex optimization problem. Its relaxation by Kantorovich is a linear optimization problem on measures, and hence on moments if the data are semialgebraic. Under convexity assumptions, this linear problem has a unique measure solution called optimal transportation plan, supported on the graph of a function called the optimal transportation map. Our method can be used to approximate separately each coordinate of the transportation map by only considering the submatrix of moments associated with a suitable marginal, extracted from an optimal solution of the semidefinite relaxation (which considers all pseudo-moments up to a given degree). In view of the form of µ, one still obtains the required convergence guarantees (e.g. pointwise), under appropriate assumptions described in the paper.

More generally, moment information about the unknown function f in the format of Problem 1 is available when applying the Moment-SOS hierarchy [START_REF] Lasserre | The Moment-SOS Hierarchy[END_REF] to solve Generalized Moment Problems where the involved Borel measures are Young measures. The necessary moment information is also given when considering empirical measures [START_REF] Lasserre | Sorting out typicality with the inverse moment matrix sos polynomial[END_REF][START_REF] Lasserre | The empirical Christoffel function with applications in data analysis[END_REF][START_REF] Pauwels | Data analysis from empirical moments and the Christoffel function[END_REF] if input data points lie on the graph of an unknown function f (e.g., as is the case in interpolation). On the other hand, in some other applications like image processing or shape reconstruction, moment information is available only for the measure f (x) dx, i.e. y → b(x, y) is linear. Finally, our method and results would apply to measures ν supported on the graph of f , different from µ in [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF]. Provided that ν and µ are mutually absolutely continuous with bounded densities, we would recover similar convergence results modulo constants.

Contribution

We address Problem 1 by providing an algorithm to approximate a (possibly discontinuous) unknown function f from the moment matrix of the measure (1) supported on its graph. The approximation converges to f (in a suitable sense described later on) when the number of moments tends to infinity.

Proposed approximation scheme

As is well-known in approximation theory, the sequence of Christoffel-Darboux polynomials associated with a measure is an appropriate tool to approximate accurately the support of the measure, hence the graph of f in our case. Christoffel-Darboux kernels and functions are closely related to orthogonal polynomials [START_REF] Szegö | Orthogonal polynomials[END_REF], [START_REF] Dunkl | Orthogonal polynomials of several variables[END_REF] and approximation theory [START_REF] Nevai | Géza Freud, orthogonal polynomials and Christoffel functions. a case study[END_REF], [START_REF] Marchi | Multivariate Christoffel functions and hyperinterpolation[END_REF]. Their asymptotic behavior (i.e., when the degree goes to infinity) provides useful and even crucial information on the support and the density of the underlying measure. A quantitative analysis is provided in [START_REF] Máté | Bernstein's inequality in L p for 0 < p < 1 and (C, 1) bounds for orthogonal polynomials[END_REF], [START_REF] Gustafsson | Bergman polynomials on an archipelago: estimates, zeros and shape reconstruction[END_REF] for single dimensional problem and in [START_REF] Kroó | Christoffel functions and universality in the bulk for multivariate orthogonal polynomials[END_REF] in a multivariate setting. Even more recently, in [START_REF] Lasserre | Sorting out typicality with the inverse moment matrix sos polynomial[END_REF], [START_REF] Lasserre | The empirical Christoffel function with applications in data analysis[END_REF] and [START_REF] Pauwels | Data analysis from empirical moments and the Christoffel function[END_REF], Christoffel-Darboux polynomials have been used to approximate the support of Borel measures in a multivariate setting in the context of Machine Learning and Data Science.

We propose a simple scheme to approximate the graph of f based on the knowledge of the moment matrix of degree 2d of µ. To that end we first compute the Christoffel-Darboux polynomial using a spectral decomposition of the moment matrix. The Christoffel-Darboux polynomial is an SOS polynomial q d (x, y) of degree 2d in p variables. For every fixed x we define f d (x) := argmin y∈Y q d (x, y)

which is a semi-algebraic function1 , assuming for the moment for the ease of exposition that the above argmin is uniquely defined. This class of functions is quite large. For example, all polynomials of degree at most d can be expressed using this technique: let r be a polynomial in x of degree d, then q : (x, y) → (r(x) -y) 2 is a degree 2d SOS polynomial whose partial minimization in y yields y = r(x) for all x. However, this class contains many more functions, including non-smooth semi-algebraic functions such as signs or absolute values. In particular this class of functions can be used to describe efficiently discontinuous functions, a typical case encountered in e.g. OCP problems with bang-bang controls and solutions with shocks for nonlinear PDEs.

Example 1 (Sign function as SOS partial minimum). Consider the polynomial

p 1 : R 2 → R (x, y) → 4 -3xy -4y 2 + xy 3 + 2y 4 .
(

) 2
One can easily check that

argmin y∈Y p 1 (x, y) = sign(x) :=    -1 if x < 0 {-1, 1} if x = 0 1 if x > 0 for any x ∈ X := [-1, 1] and Y ⊂ R. Note that since p 1 (x, •) is positive for x ∈ [-1, 1],
it can be squared without changing the argmin in y and hence we obtain a similar representation of the sign function in the form of partial minimization in y of the degree 8 SOS polynomial q(x, y) := p 2 1 (x, y), represented in Figure 1. Example 2 (Absolute value as SOS partial minimum). The reader may check that the argmin in y of the (square of the positive) polynomial 11 -12x 4 y -6x 2 y 2 + 4x 2 y 3 + 3y 4 is equal to |x| for all x ∈ [-1, 1].

To the best of our knowledge, this work is the first contribution where this class of semialgebraic functions is used for graph approximations. The present work shows how these functions may be used to approximate discontinuous functions accurately.

Comparison to existing approximation approaches:

A classical alternative to our approach would be to use L 2 -norm Legendre or Chebyshev approximations of the function f which are also based on moment information. However these approaches only use moment information on the measure f (x) dx, i.e. y → b(x, y) is linear. Moreover, the support of this measure is not the graph of f .

We claim that the full moment information provides useful additional data on the graph of f which we can access using Christoffel-Darboux kernels associated with the measure µ in [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF]. Note that in interpolation applications, we have the possibility to estimate the higher order moments of f from finitely many evaluations of f through Riemann integral approximations or Monte-Carlo approximations for example. However, in situations where we have neither access to higher moments nor pointwise evaluation of f , our method cannot be applied; signal processing applications are a typical case of the latter situations.

In general, approximating a discontinuous function f is a challenge. Indeed, most well-known techniques suffer from the Gibbs phenomenon, i.e. the approximation produces oscillations at each point of discontinuity of f , see e.g. [START_REF] Gottlieb | On the Gibbs phenomenon and its resolution[END_REF] for a good survey on this topic. The main tools usually rely on properties of orthogonal polynomials [START_REF] Szegö | Orthogonal polynomials[END_REF] and the resulting approximations are based on a finite number of Fourier coefficients of the latter functions, i.e., typically first degree moment information on f . Projecting a discontinuous function on a class of infinitely smooth functions is the typical mechanism producing Gibbs phenomenon. In order to get rid of such a curse, additional techniques and prior information is needed. A possible approach is reported in [START_REF] Eckhoff | Accurate and efficient reconstruction of discontinuous functions from truncated series expansions[END_REF] in the univariate case (p = 2 in our notations), where a good approximation of locations of discontinuity points and jump magnitude is obtained by solving an appropriate (univariate) polynomial equation. Recent developments have effectively shown that such approaches may tame the Gibbs phenomenon [START_REF] Batenkov | Algebraic Fourier reconstruction of piecewise smooth functions[END_REF][START_REF] Batenkov | Complete algebraic reconstruction of piecewise-smooth functions from Fourier data[END_REF]. Iterative numerical methods can also be used to identify the points of discontinuities of f (and of its derivatives) so as to construct accurate approximations locally in each identified interval, see e.g. [START_REF] Pachón | Piecewise-smooth chebfuns[END_REF] in the case of Chebyshev polynomials. However such ad-hoc techniques are very specific to the univariate setting.

On the contrary, our approach is not based on projection on a subspace of smooth functions, or identification of points of discontinuities of the function to be approximated. It is based on geometric approximation of the support of a singular measure using semi-algebraic techniques. An important feature of this approach is that the resulting approximant functions are not necessarily smooth, and furthermore, discontinuities (if any) can be treated implicitly only based on the whole moment information. As a result, we obtain a single approximation scheme, which (i) may adapt to the smoothness features of the target function f without requiring prior knowledge of it, and (ii) can be used for multivariate f , both points being important challenges regarding numerical approximation.

Description of our contribution

1. We first provide a numerical scheme which allows to approximate the compact support of a measure which is singular with respect to the Lebesgue measure. We need to adapt the strategy of [START_REF] Lasserre | The empirical Christoffel function with applications in data analysis[END_REF] which covered the absolutely continuous case. This result may be considered of independent interest and will be instrumental to providing convergence guarantees for our approach.

2. Next, given a degree d ∈ N, we provide an approximation f d of the function f and prove that the sequence (f d ) d∈N converges pointwise to f almost everywhere as well as in L 1norm as d goes to infinity (under broad assumptions on f ). Furthermore, if we assume more regularity on f , then we also provide estimates for the rate of convergence. More precisely, we obtain O(d -1/2 ) for multivariate Lipschitz functions and O(d -1/4 ) for univariate functions of bounded variation.

3. Finally, we provide some numerical examples to illustrate the efficiency of the method. We first use our algorithm to approximate discontinuous or non-smooth solutions of OCP or PDE problems based on the Moment-SOS hierarchy. These experiments empirically demonstrate the absence of Gibbs phenomena. We also provide an example where only samples of the measure under consideration are given and show that our algorithm also works well, even for moderate size samples, showing that moment input data could in principle be approximated using numerical integration methods.

Example 3 (Sign function). To give a flavor of what can be obtained numerically, consider the measure (1) supported on the graph {(x, f (x)) :

x ∈ X} ⊂ X × Y of the function x → f (x) := sign(x), with X := Y := [-1, 1].
In Figure 2 (right) the resulting approximation f 2 with a moment matrix of size 6 and degree 4 (i.e. 15 moments) cannot be distinguished from f . On the other hand, on the left, their Chebsyhev interpolants of degrees 4 and 20 (obtained with chebfun [START_REF] Driscoll | Chebfun Guide[END_REF]) illustrate the typical Gibbs phenomenon, namely oscillations near the discontinuity points that cannot be attenuated by increasing the degree. This phenomenon can be reduced or suppressed by identifying the discontinuity points and splitting X into intervals (as in e.g. [START_REF] Pachón | Piecewise-smooth chebfuns[END_REF] and also implemented in chebfun), but this strategy works only in the univariate case. In contrast, our algorithm does not attempt to approximate directly a univariate function with one or several univariate polynomials of increasing degree, but with the argmin of a bivariate polynomial of increasing degree. Moreover, our algorithm works also for multivariate functions, as shown by numerical examples later on.

Organization of the paper

Section 2 introduces the Christoffel-Darboux polynomial, its regularized version and our semialgebraic approximant. In Section 3 the main results of the paper are collected, while their proofs are provided in Section 4. More precisely, we first give some quantitative estimates on the support of the measure µ and then prove the L 1 convergence of our approximant. Section 5 discusses computational issues and with the help of a simple Matlab prototype it illustrates the efficiency of our method on some examples. Finally, Section 6 collects some concluding remarks together with further research lines to be followed.

Notation

The Euclidean space of real-valued symmetric matrices of size n is denoted by S n . Given a set X in Euclidean space, its diameter is denoted by diam(X) and its volume, or Lebesgue measure, is denoted by vol(X). For k ≥ 1, the Lebesgue space L k (X) consists of functions on X whose k-norms are bounded. Given a positive Borel measure µ, we denote by supp(µ) its support, defined as the smallest closed subset whose complement has measure zero. Throughout the paper, p denotes the dimension of the ambient space. Consistently with the notations introduced in Section 1 for stating Problem 1, we let z = (x y) ∈ R p . We denote by R[z] the algebra of multivariate polynomials of z ∈ R p with real coefficients. For a given degree d ∈ N, the dimension of the vector space of polynomials of degree less than or equal to d is denoted by

n d := p + d d .
2

Christoffel-Darboux approximation

This section describes our main approximant based on the Christoffel-Darboux kernel. We first introduce notations and definitions, describe our regularization scheme for the Christoffel-Darboux kernel and then describe our functional approximant based on moments.

Polynomials and moments

Following the notations introduced for stating Problem 1, any polynomial q ∈ R[z] of degree d can be expressed in the polynomial basis b(z) as q : z → q b(z) with q ∈ R n d denoting its vector of coefficients. Recall that the moment matrix of degree 2d of the measure µ reads

M µ,d = b(z)b(z) dµ(z) ∈ S n d . Since M µ,d is positive semi-definite, it has a spectral decomposition M µ,d = PEP , (3)
where P ∈ R n d ×n d is an orthonormal matrix whose columns are denoted p i , i = 1, 2, . . . , n d , and satisfy p i p i = 1 and p i p j = 0 if i = j, and E ∈ S n(d) is a diagonal matrix whose diagonal entries are eigenvalues e i+1 ≥ e i ≥ 0 of the moment matrix. Each column

p i ∈ R n d is the vector of coefficients of a polynomial p i ∈ R[z], i = 1, . . . n d , so that p i M µ,d p i = e i = p 2 i (z)dµ(z), p i M µ,d p j = 0 = p i (z)p j (z)dµ(z), i = j. (4)

Approximating the support from moments

Let us assume for now that the support of the measure µ has nonempty interior, then M µ,d is positive definite for any d ∈ N, i.e., e i > 0, i = 1, . . . , n d . In this case, one can define the Christoffel-Darboux polynomial

q µ,d (z) := n d i=1 p 2 i (z) e i = b(z) M -1 µ,d b(z). (5)
It is known that sublevel sets of q µ,d can be used to recover the support of µ for large d, see for example [START_REF] Lasserre | The empirical Christoffel function with applications in data analysis[END_REF] for an overview. The goal of this work is to approximate the function f . From a set theoretic perspective, this amounts to approximating the graph of f whose closure is actually the support of the measure µ in [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF]. Hence the sublevel sets of q µ,d are interesting candidates for this goal. However, in the case of the graph of a function, the construction given in ( 5) is not valid anymore since this graph is a singular set so that M µ,d may not be positive definite and invertible. In this singular setting, one should ideally consider the following extended value Christoffel-Darboux polynomial:

q e µ,d (z) := +∞ if ∃i, e i = 0, p i (z) = 0 ei>0 p 2 i (z) ei = b(z) M † µ,d b(z) otherwise, (6)
where † denotes the Moore-Penrose pseudo inverse. This is a natural extension of the Christoffel-Darboux polynomial to the singular case [START_REF] Pauwels | Data analysis from empirical moments and the Christoffel function[END_REF] and amounts to working in the Zariski closure of the graph of f .

Regularization scheme

Spectral filtering: Computing an object such as in ( 6) can be numerically sensitive since it essentially relies on pseudo-inverse which requires an eigenvalue thresholding scheme. This means that small perturbations of the moment matrix may lead to large changes of the output. Furthermore the candidate function takes finite values only on an algebraic set, and this situation is difficult to handle in finite precision arithmetic. One may rewrite the extended value polynomial (6) in the following form

q e µ,d (z) = n d i=1 g(e i )p 2 i (z)
where g : [0, +∞) → [0, +∞] with g(s) = 1 s for any s > 0 and g(0) = +∞. One approach to restore stability is to use regularization techniques which replace the pseudo-inversion expressed through the function g, by spectral filtering expressed through a different spectral function (see for example [START_REF] Vito | Learning sets with separating kernels[END_REF] for an illustration in the context of support estimation). Since the function g is not regular, instead of studying the above defined extended value polynomial, one rather looks at the following polynomial

n d i=1 g β (e i )p 2 i (z) (7)
where g β is a parametrized family of spectral filtering regular functions satisfying, for any β > 0,

g β : [0, +∞) → [0, +∞]. Common examples include
Tikhonov regularization:

g β : s → 1 β + s , Spectral cut-off: g β : s → 1 β I [0,β] (s) + 1 s I (β,+∞) (s),
Ideal low-pass:

g β : s → 1 β I [0,β] (s).
We choose to work with the Tikhonov regularization as it has an intuitive measure space intepretation. We believe that our results can be generalized to different spectral filters.

Tikhonov regularization and measures:

Applying the Tikhonov spectral filter to [START_REF] Vito | Learning sets with separating kernels[END_REF] yields the following polynomial

z → n d i=1 p 2 i (z) e i + β = b(z) (M µ,d + βI n d ) -1 b(z) (8)
where I n d denotes the identity matrix of size n d . In order to use analytic tools, we need to provide an interpretation of the addition of diagonal elements in terms of measures. One therefore has to choose a polynomial basis for which the diagonal matrix is the moment matrix of a reference Borel measure on R p that we will denote µ 0 . We make the following assumption which will be standing throughout the paper.

Assumption 1 (Reference measure and polynomial basis).

• The reference measure µ 0 is absolutely continuous with respect to the Lebesgue measure and it has compact support. • The polynomial basis b is orthonormal with respect to the bilinear form induced by µ 0 , that is b i (z)b j (z)dµ 0 (z) = 1 if i = j and 0 otherwise.

The first part of Assumption 1 ensures that the moment matrix of µ 0 is always positive definite. The second part of Assumption 1 provides the following relation:

M µ,d + βI n d = M µ+βµ0,d . (9)
An easy example of such a measure should be the following: considering a function f whose domain of definition is contained in the unit cube of dimension p -1 and which takes values in [-1, 1], one might pick the uniform measure on the unit cube of dimension p, for which moments are easy to compute.

Most importantly, using the notation in [START_REF] Coste | An introduction to semialgebraic geometry[END_REF], this allows to express the polynomial of interest (8) as follows.

Definition 1 (Regularized Christoffel-Darboux polynomial). The regularized Christoffel-Darboux polynomial is the SOS

q µ+βµ0,d (z) := n d i=1 p 2 i (z) e i + β . (10)
This provides a geometric interpretation of the regularization parameter as a combination of two measures: µ which is supported on the graph of the function of interest and µ 0 which is a reference measure, used for regularization purposes. The supported boundedness hypothesis in Assumption 1 will allow to provide quantitative estimates in further sections and it could in principle be replaced by a fast decreasing tail condition. An important example for measures satisfying Assumption 1 is the restriction of Lebesgue measure to the unit hypercube together with Legendre polynomials which form an orthonormal basis.

Making Assumption 1 is a slight restriction for which a few comments are in order. Firstly, this could be relaxed in various ways to remove the restriction on the polynomial basis, for example:

• Replace the identity matrix by the moment matrix of µ 0 ; • Add assumption on the spectrum of the moment matrix µ 0 .

These would lead to a lot of technical complications and we find our results clearer and easier to state under Assumption 1. Secondly, working numerically with polynomials in the standard monomial basis is problematic in many situations. Better conditioned polynomial bases are often those enjoying orthonormality properties with respect to a certain reference measure, such as e.g. Chebyshev or Legendre polynomials. We would like to argue here that the restrictions induced by Assumption 1 are quite benign since it is already common in practice to work in such polynomial bases for numerical reasons.

Semi-algebraic approximant

Definition 2 (Semi-algebraic approximant). The regularized Christoffel-Darboux semialgebraic approximant f β,d is defined as follows:

x ∈ X → f β,d (x) := min{argmin y∈Y q µ+βµ0,d (x, y)}.

(11)

Remark 1. If X and Y are semi-algebraic, then the set-valued map which associates to each x ∈ X the set argmin y∈Y q µ+βµ0,d (x, y) is semi-algebraic. Recall that a map is semi-algebraic if its graph is semi-algebraic. By the Tarski-Seidenberg Theorem (see for example [5, Theorem 2.6]), any first order formula involving semi-algebraic sets describes a semi-algebraic set. Since minima are described by first order formulas, the argmin of a polynomial on the compact semi-algebraic set Y is a compact semialgebraic subset of Y , which is itself a subset of the real line. Hence the argmin set has a minimal element and the function in [START_REF] Eckhoff | Accurate and efficient reconstruction of discontinuous functions from truncated series expansions[END_REF] is well defined.

Remark 2. For clarity of exposition we describe our main results by considering that the partial minimization in y over Y is exact. As will be seen from arguments in the proof, especially in the proof of Lemma 7, approximate minimization up to a factor of the order d p+2 enjoys similar approximation properties. Indeed, in Remark 5, we mention the precision γ d 2 , where γ d is chosen later on with more justifications in [START_REF] Gottlieb | On the Gibbs phenomenon and its resolution[END_REF] to be of order d p+2 .

The two parameters d and β control the behavior of the approximant f β,d . In latter sections, we describe an explicit dependency between d and β which allows to construct a sequence of regularization parameters (β d ) d∈N , and we investigate the asymptotic properties of the sequence of approximants (f β d ,d ) d∈N .

Main results

Our main theoretical contribution is an investigation of the relations between the function f to be approximated and its regularized Christoffel-Darboux approximant f β,d under Assumption 1. In particular we are interested in building an explicit sequence (β d ) d∈R and investigating the convergence

f β d ,d (x) → f (x) for x ∈ X, fixed, as well as the convergence f -f β d ,d L 1 (X) → 0, when d → ∞.

Convergence under continuity assumptions

The following section describes our main result regarding convergence of the approximant f β d ,d in [START_REF] Eckhoff | Accurate and efficient reconstruction of discontinuous functions from truncated series expansions[END_REF] under different regularity assumptions for the function f to be approximated. Let us define

δ 0 := diam(supp(µ + µ 0 )), m := µ(R p ), m 0 := µ 0 (R p ).
Theorem 1. Under Assumption 1 and with the choice

β d = 2 3- √ d in Definition 2, it holds: (i) If the set S ⊂ X of continuity points of f is such that X \ S has Lebesgue measure zero, then f β d ,d (x) → d→∞ f (x)
for almost all x ∈ X, and

f -f β d ,d L 1 (X) → d→∞ 0.
(ii) If f is L-Lipschitz on X for some L > 0, then for any d > 1 and any r > p,

f -f β d ,d L 1 (X) ≤ vol(X) δ 0 √ d -1 (1 + L) + diam(Y ) 8(m + m 0 )(3r) 2r e p 2 d
p p e 2r-p d r-p .

Remark 3.

• As proved recently in [START_REF] Silvestre | Oscillation properties of scalar conservation laws[END_REF], the solutions to scalar conservation laws are continuous almost everywhere, i.e. the Lebesgue measure of the set of discontinuity points reduces to 0.

It is then clear that Item (i) of Theorem 1 can be applied directly to the case of scalar conservation laws.

• Thanks to Egorov's Theorem, pointwise convergence almost everywhere implies almost uniform convergence, that is uniform convergence up to a subset of X whose measure can be taken arbitrarily small. Since we manipulate bounded functions, this in turn implies convergence in L 1 . • For Lipschitz continuous functions, we obtain an O(d -1/2 ) convergence rate in L 1 norm by setting r = p + 1/2. • The convergence rate for Lipschitz functions is slow and we observe in practice a much faster convergence. We conjecture that faster rates can be obtained for special classes of functions.

Theorem 1 is a special case of a more general result described in Theorem 3 and proven later on.

Convergence for univariate functions of bounded variation

Spaces of functions of bounded variations are of interest because many PDE problems are formulated on such spaces, see [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF] for an introduction. Modern construction of such spaces is done by duality through measure theoretic arguments. Our main proof mechanisms rely on pointwise properties of the function f , which are not completely captured by the measure theoretic construction.

We prove L 1 convergence for univariate bounded variation functions. The reason we are limited to the univariate setting is that we can use the classical definition of total variation which is directly connected to pointwise properties of the function of interest. We conjecture that the proposed approximation scheme is also convergent for multivariate functions of bounded variation, but we leave this question for future work. Definition 3. Let f : R → R be a measurable function. The (Jordan) total variation norm of f is given by

V (f ) = sup n∈N sup x0<x1<x2<...<xn n i=1 |f (x i ) -f (x i-1 )|.
f -f β d ,d L 1 (X) ≤ vol(X) 2δ 0 √ d -1 + d -1 4 + diam(Y ) 8(m + m 0 )(3r) 2r e 4 d 4e 2r-2 d r-2 + 4d 1 4 V (f )δ 0 √ d -1 .
We remark that we obtain a convergence rate in O(d -1/4 ). This result is a special case of a more general result described in Theorem 4 and proven later on.

Robustness to small perturbations

In many situations, one only has access to an approximation of the regularized Christoffel-Darboux polynomial. This is for example the case when applying the Moment-SOS hierarchy. At the end, one indeed obtains pseudo-moments of the measure under consideration, i.e. a vector of real numbers close to the actual moments of the measure. The moment matrix is then not M µ+β d µ0 as in [START_REF] Driscoll | Chebfun Guide[END_REF] but a matrix M close to it. The effect of this perturbation on the Christoffel-Darboux polynomial is captured by the following result.

Lemma 1. Assume that the approximate moment matrix M ∈ S n d is positive definite and let

α := I n d -M 1 2 µ+β d µ0 M -1 M 1 2 µ+β d µ0
where we used the matrix operator norm. Then the polynomial q α

d : z → b(z) M -1 b(z) satisfy sup z∈R p 1 - q α d (z) q µ+β d µ0,d (z) ≤ α.
Proof : For any z ∈ R p , we have

q α d (z) q µ+β d µ0,d (z) -1 = 1 q µ+β d µ0,d (z) b(z) (M -1 -M -1 µ+β d µ0 )b(z) = 1 q µ+β d µ0,d (z) (M -1 2 µ+β d µ0 b(z)) (I n d -M 1 2 µ+β d µ0 M -1 M 1 2 µ+β d µ0 )M -1 2 µ+β d µ0 b(z) ≤ I n d -M 1 2 µ+β d µ0 M -1 M 1 2 µ+β d µ0 = α. 2
Remark 4 (Accuracy of the moments). In this paper, we assume the existence of a positive number α which makes the link between an approximated moment matrix and the real one, which exists in general, but estimates for α are not known in general. Such an analysis has been performed for the moments themselves in [START_REF] Nie | On the complexity of Putinar's Positivstellensatz[END_REF], but not for moment matrices. Indeed, the bounds linking moment matrices and their corresponding depend nonlinearly on too many variables to obtain easily bounds on the approximated moment matrix and the real one. This is a topic of further investigation.

More generally, we can consider a robust Christoffel-Darboux function satisfying the following inequality.

Assumption 2. For a given α ∈ [0, 1), let (β d ) d∈N be a sequence of positive numbers and (q α d ) d∈N be a sequence of continuous functions over R p such that for any d ∈ N and any z ∈ R p , we have

1 - q α d (z) q µ+β d µ0,d (z) ≤ α.
Note that Assumption 2 ensures that

(1 -α)q µ+β d µ0,d (z) ≤ q α d (z) ≤ (1 + α)q µ+β d µ0,d (z). ( 12
)
Furthermore, one can always choose q α d (z) = q µ+β d µ0,d (z) and α = 0 which corresponds to the nominal case. The robust approximant is then defined similarly as in Definition 2.

Definition 4 (Robust semi-algebraic approximant). Given a degree d ∈ N, a regularizing parameter β > 0, and a scalar α ∈ [0, 1), our robust approximant f α β d ,d is defined as follows:

x ∈ X → f α β d ,d (x) := min{argmin y∈Y q α d (x, y)}.

The main technical result of this paper is the following from which Theorem 1 directly follows. (i) If the set S ⊂ X of continuity points of f is such that X \ S has Lebesgue measure zero, then

f α β d ,d (x) → d→∞ f (x)
for almost all x ∈ X, and

f -f α β d ,d L 1 (X) → d→∞ 0.
(ii) If f is L-Lipschitz on X for some L > 0, then for any d > 1 and any r > p,

f -f α β d ,d L 1 (X) ≤ vol(X) δ 0 √ d -1 (1 + L) + diam(Y ) 1 + α 1 -α 8(m + m 0 )(3r) 2r e p 2 d
p p e 2r-p d r-p .

Furthermore, we have the following robust convergence result for univariate functions of bounded variation, from which Theorem 2 directly follows.

f -f α β d ,d L 1 (X) ≤ vol(X) 2δ 0 √ d -1 + d -1 4 + diam(Y ) 1 + α 1 -α 8(m + m 0 )(3r) 2r e 4 d 4e 2r-2 d r-2 + 4d 1 4 V (f )δ 0 √ d -1 .
The next section is dedicated to the proof of these theorems.

Proofs

This section is divided into several subsections. Subsection 4.1 gives some quantitative results on the estimation of the support of µ which does not depend on the nature of µ and could be of independent interest. More precisely, we show that the regularized Christoffel-Darboux polynomial takes large values outside the support of µ and smaller values inside. This is expressed by describing properties of certain sublevel sets of the polynomial being close to the support of µ.

In subsection 4.2 we translate these geometric results in functional terms. In subsection 4.3, we prove our main results: the argument of the minimum of the regularized Christoffel-Darboux polynomial is close to the graph of the function f .

Estimation of the support

In this section we build a polynomial sublevel set that will be instrumental for our proofs of convergence in functional terms. Note that in practice this sublevel set is not computed: we just focus on the argmin of the regularized Christoffel-Darboux polynomial. The contents of this section may be considered of independent interest. For any d ∈ N and r ∈ N such that r > p and for any α ∈ [0, 1) , define

γ d := 1 -α 8(m + m 0 ) e 2r d r (3r) 2r (14)
Figure 3: Illustration of the result of Theorem 5, the dotted curve represents the boundary of the set S d and the considered function f is the absolute value. The theorem states that (i) most of the points of the graph of f will be in S d and that (ii) all points in S d will be close to the graph of f . More precisely, (i) states that the measure of R 1 ∪ R 2 ∪ R 3 will vanish and (ii) the distance r of any point in S d to the graph of f will go to zero for d → ∞, respectively.

and

S d := {z ∈ R p : q α d (z) < γ d }. (15)
We aim at proving that the sublevel set S d is approaching the support of µ as d goes to infinity, with a given convergence rate. It is precisely quantified with the following result which is illustrated in Figure 3.

Theorem 5. For d > 1 it holds (i) µ({z ∈ R p : z ∈ S d }) ≤ 1 + α 1 -α 8(m + m 0 )(3r) 2r e p 2 d
p p e 2r-p d r-p .

(ii) For any z ∈ S d ,

dist(z, supp(µ)) ≤ δ 0 √ d -1 .
The results and techniques that we will use to prove this theorem are adapted from [START_REF] Lasserre | The empirical Christoffel function with applications in data analysis[END_REF] which considers the absolutely continuous setting without regularization. Proof of Theorem 5 (i) Using ( 10) and ( 4), we obtain

R p q µ+β d µ0,d (z)dµ(z) = n d i=1 e i e i + β d ≤ n d ≤ d p e p p e p 2 d (16)
where the last inequality is given in [START_REF] Lasserre | The empirical Christoffel function with applications in data analysis[END_REF]Lemma 6.5]. Using Markov's inequality [36, Page 91] and (16) yields

µ({z ∈ R p : q µ+β d µ0,d (z) ≥ γ d 1 + α }) ≤ R p (1 + α)q µ+β d µ0,d (z)dµ(z) γ d ≤ (1 + α) d p e p p e p 2 d γ d . (17)
Now using ( 12) and ( 15), we have the following implications

z ∈ S d ⇔ q α d (z) ≥ γ d ⇒ q µ+β d µ0,d (z) ≥ γ d 1 + α . Hence µ({z ∈ R p : z ∈ S d }) ≤ µ({z ∈ R p : q µ+β d µ0,d (z) ≥ γ d 1+α }).
Using the expression of γ d in ( 14) and the inequality [START_REF] Hernández-Hernández | The linear programming approach to deterministic optimal control problems[END_REF], one has:

µ({z ∈ R p : z ∈ S d }) ≤ 1 + α 1 -α 8(m + m 0 )(3r) 2r e p 2 d
p p e 2r-p d r-p which concludes the proof of item (i) of Theorem 5. 2

To carry out the proof of Theorem 5 (ii), we begin with a few lemmas. The following result is classical, see e. Lemma 2. Let d ∈ N, z ∈ R p , β > 0, and q be a polynomial of degree at most d. Then

q 2 (z) R p q 2 (z)d(µ + βµ 0 )(z) ≤ q µ+βµ0,d (z).
The following Lemma defines the needle polynomial, introduced first in [START_REF] Kroó | Christoffel functions and universality in the bulk for multivariate orthogonal polynomials[END_REF], and gives a quantitative result crucial for our analysis.

Lemma 3 (Existence of a needle polynomial). Let B δ denote the euclidean ball of radius δ. Then, for all δ ∈ (0, 1) and d ∈ N, d > 0, there exists a polynomial q of degree 2d such that

q(0) = 1, q(z) ∈ [-1, 1] for all z ∈ B 1 , and |q(z)| ≤ 2 1-δd for all z ∈ B 1 \ B δ . (18)
A detailed proof is provided in [START_REF] Lasserre | The empirical Christoffel function with applications in data analysis[END_REF]Lemma 6.3]. Thanks to the latter lemma, we can characterize the behavior of the regularized Christoffel-Darboux polynomial q µ+β d µ0,d outside the support of µ:

Lemma 4. Let d ∈ N, d > 1 and z ∈ R p . Recall that β d := 2 3- √ d . Assume that dist(z, supp(µ)) ≥ δ0 √ d-1 . Then 2 √ d-3 m + m 0 ≤ q µ+β d µ0,d (z). ( 19
)
Proof of Lemma 4: Let d > 1, z ∈ R p , δ = dist(z, supp(µ)). Let d ∈ N and t > 0, arbitrary for the moment. Consider the affine map T : w → w-z δ+δ0 . Let q be the degree 2d polynomial given as in Lemma 3 such that

q(0) = 1, q(w) ∈ [-1, 1] for all w ∈ B 1 , and
|q(w)| ≤ 2 1-δd δ+δ 0 for all w ∈ B 1 \ B δ δ+δ 0 . ( 20
)
Let r = q • T . The polynomial r satisfies

|r(z )| ≤ 1, ∀z ∈ supp(µ + µ 0 ), r(z ) ≤ 2 1-δd δ+δ 0 , ∀z ∈ supp(µ), r(z) = 1. ( 21
)
Using Lemma 2 and the fact that r(z) = 1 we obtain

R p r 2 (w)d(µ + βµ 0 )(w) -1
≤ q µ+βµ0,2d (z) [START_REF] Lasserre | The Moment-SOS Hierarchy[END_REF] and

R p r 2 (w)d(µ + βµ 0 )(w) -1 ≤ q µ+βµ0,2d +1 (z) (23)
From ( 21), we deduce

R p r 2 (w)d(µ + βµ 0 )(w) ≤ 2 2-δ(2d ) δ+δ 0 m + βm 0 (24) R p r 2 (w)d(µ + βµ 0 )(w) ≤ 2 3-δ(2d +1) δ+δ 0 m + βm 0 (25)
Combining ( 22), ( 23) , ( 24) and ( 25), we obtain the following bounds

q µ+βµ0,2d (z) ≥ 2 3-δ(2d ) δ+δ 0 m + βm 0 -1 , q µ+βµ0,2d +1 (z) ≥ 2 3-δ(2d +1) δ+δ 0 m + βm 0 -1 . ( 26
)
Recall that d and β were arbitrary. Now we can choose d = d/2 , β = β d in one of the identities in ( 26) (depending on the parity of d) to obtain

q µ+β d µ0,d (z) ≥ 2 3-δd δ+δ 0 m + β d m 0 -1 ≥ 2 √ d-3 m + m 0 , (27)
where the last inequality follows because the right hand side is strictly increasing as a function of δ and δ ≥ δ0

≥ (2r) (1 -log(3r)) .
Proof of Lemma 5: A simple analysis shows that the minimum is attained at x = 2r log [START_REF] Bareiss | Sylvester's identity and multistep integer-preserving gaussian elimination[END_REF] . The lower bound follows because 2 log(2) ≤ 3. 2 Lemma 6. For any d ∈ N, we have

2 √ d-3 m + m 0 ≥ γ d 1 -α .
Proof of Lemma 6: Using the definition of γ d in ( 14), we have

log γ d 8(m + m 0 ) 1 -α = 2r(1 -log(3r) + log( √ d)) ≤ log(2) √ d, (28)
where the inequality follows from Lemma 5 with x = √ d. This proves the desired result. 2

Proof of Theorem 5 (ii): We prove the result by contraposition. We have the following chain of implications for z

∈ R p , dist(z, supp(µ)) ≥ δ 0 √ d -1 ⇒ 2 √ d-3 m + m 0 ≤ q µ+β d µ0,d (z) ⇒ γ d /(1 -α) ≤ q µ+β d µ0,d (z) ⇒ γ d ≤ q α d (z) ⇒ z ∈ S d , (29)
where the first implication is from Lemma 4, the second implication is due to Lemma 6, the third implication is deduced from [START_REF] Fattorini | Infinite dimensional optimization and control theory[END_REF] and the last implication is from the definition of S d in [START_REF] Gustafsson | Bergman polynomials on an archipelago: estimates, zeros and shape reconstruction[END_REF]. 2

Estimation of functions

We now translate Theorem 5 in functional terms. Considering that z can be written as follows z = (x, y), let us introduce a specific set which will be of interest throughout the proof:

I d := {x ∈ X : inf y∈Y q α d (x, y) ≥ γ d }. ( 30
)
Lemma 7. Suppose that d > 1. Then, we have

I d dx ≤ 1 + α 1 -α 8(m + m 0 )(3r) 2r e p 2 d p p e 2r-p d r-p . ( 31
)
Proof of Lemma 7: For all x ∈ X and all A ⊂ R p measurable, one has

Y I A (x, y)δ f (x) (dy) = I A (x, f (x)) (32)
and hence

µ(A) = Y I A (x, y)dµ(x, y) = X I A (x, f (x))dx = I A dx, (33)
where

I A := {x ∈ X : (x, f (x)) ∈ A}. One can see from (30) that x ∈ I d implies that (x, f (x)) ∈ S d and hence I d ⊂ I S c
d where S c d denotes the complement of S d given in [START_REF] Gustafsson | Bergman polynomials on an archipelago: estimates, zeros and shape reconstruction[END_REF]. We deduce that

I d dx ≤ I S c d dx = µ(S c d ) ( 34
)
and the result follows from Item (i) of Theorem 5. 2

Remark 5. Letting Ĩd := {x ∈ X : inf y∈Y q α d (x, y) ≥ γ d /2} it can be seen using the exact same arguments that a bound on Ĩd dx holds similarly as in Lemma 7 with a multiplicative factor of 2. This can be used to handle the situation where the argmin in [START_REF] Eckhoff | Accurate and efficient reconstruction of discontinuous functions from truncated series expansions[END_REF] is computed up to a precision of the order γ d /2. See also Remark 2.

Thanks to Lemma 7, it is sufficient to prove the convergence of the approximated function in the set d as in [START_REF] Golub | Matrix Computations. 3rd Edition[END_REF] and let I d be defined by [START_REF] Pachón | Piecewise-smooth chebfuns[END_REF]. Then for any r > p, we have

I c d := X \ I d . Proposition 1. Let d ∈ N, d > 1. Let α ∈ [0, 1) be as in Assumption 2 with β d = 2 3- √ d . Consider f α β d ,
X |f (x) -f α β d ,d (x)|dx ≤ I c d |f (x) -f α β d ,d (x)|dx + diam(Y ) 1 + α 1 -α 8(m + m 0 )(3r) 2r e p 2 d
p p e 2r-p d r-p .

Proof of Proposition 1: Since y takes values in the compact set Y , it is clear that

f α β d ,d ∈ L ∞ (X). We have X |f (x) -f α β d ,d (x)|dx = X\I d |f (x) -f α β d ,d (x)|dx + I d |f (x) -f α β d ,d (x)|dx ≤ I c d |f (x) -f α β d ,d (x)|dx + f -f α β d ,d L ∞ (I d ) I d dx ≤ I c d |f (x) -f α β d ,d (x)|dx + diam(Y ) 1 + α 1 -α 8(m + m 0 )(3r) 2r e p 2 d p p e 2r-p d r-p (35)
where we have used Lemma 7 for the last inequality. 2

Proofs of the main theorems

We are now in position of proving Theorem 3. We start with the Lipschitz case, which is the simplest and conveys most of the ideas.

Proof of Theorem 3 (ii) Using Proposition 1: It remains to bound the term

I c d |f (x) -f α β d ,d (x)|dx. (36)
For any x ∈ X define

u d (x) ∈ argmin u∈X x, f α β d ,d (x) -(u, f (u)) , (37)
with an arbitrary choice in the case where the argmin is not unique. Note that by continuity, the graph of f is closed so that the minimum is attained. Using the definition of I d in [START_REF] Pachón | Piecewise-smooth chebfuns[END_REF], the fact that

x ∈ I c d implies that x, f α β d ,d (x) ∈ S d . (38)
Moreover, Theorem 5 implies that

|f α β d ,d (x) -f (u d (x))| ≤ δ 0 √ d -1 , (39)
x -

u d (x) ≤ δ 0 √ d -1 ,
where • denotes the usual Euclidean distance in R p-1 . Therefore, using Lipschitz continuity of f , we have, for any x ∈ I c d ,

|f α β d ,d (x) -f (x)| ≤ |f α β d ,d (x) -f (u d (x))| + |f (x) -f (u d (x))| ≤ δ 0 √ d -1 (1 + L) . (40)
We deduce that

I c d |f (x) -f α β d ,d (x)|dx ≤ vol(X) δ 0 √ d -1 (1 + L) (41)
which concludes the proof. 2

We now turn to case (i), starting with the pointwise convergence.

Proof of Theorem 3 (i):

We rely on a slighlty different use of Lemma 7. Choose r = p + 2 and let

I := {x ∈ X : ∀d 0 ∈ N, ∃d ∈ N, d ≥ d 0 , x ∈ I d } = ∩ d0∈N ∪ d≥d0 I d . (42)
Lemma 7 ensures that vol(

I d ) = O(1/d 2 ) so that vol (∪ d≥d0 I d ) ≤ d≥d0 vol(I d ) → d0→∞ 0.
We have vol(I) = lim d0→∞ (∪ d≥d0 I d ) = 0. This means that we have the two following properties, for almost every x ∈ X:

• f is continuous at x (by assumption),

• ∃d 0 ∈ N, ∀d ∈ N, d ≥ d 0 , x ∈ I d (because vol(I) = 0).
Fix any such x and for any d ≥ d 0 consider

(u d , v d ) ∈ argmin u∈X,v∈Y x, f α β d ,d (x) -(u, v) , s.t. (u, v) ∈ supp(µ) (43)
with an arbitrary choice when the argmin is not unique. Note that the support of µ is actually the closure of the graph of f so that the minimum is attained. Using the definition of I d in (30), we have that x ∈ I c d for all d ≥ d 0 which implies that

x, f α β d ,d (x) ∈ S d , (44)
and Theorem 5 implies that

|f α β d ,d (x) -v d | ≤ δ 0 √ d -1 , x -u d ≤ δ 0 √ d -1 .
Since (u d , v d ) ∈ supp(µ) and supp(µ) is the closure of the graph of f , there exists h d ∈ X, such that

|f α β d ,d (x) -f (h d )| ≤ 2δ 0 √ d -1 , x -h d ≤ 2δ 0 √ d -1 . (45)
This concludes the proof of pointwise convergence since

|f α β d ,d (x) -f (x)| ≤ |f α β d ,d (x) -f (h d )| + |f (x) -f (h d )|
and both terms tend to 0 as d → ∞, using (45) and continuity of f at x. Convergence in L 1 follows from Egorov's Theorem (see e.g. [33, chapter 18]): for any > 0, there exists S ⊂ X, measurable, of Lebesgue measure smaller than such that f α β d ,d → f uniformly on X \ S . We have

f α β d ,d -f L 1 (X) = X |f α β d ,d (x) -f (x)|dx = S |f α β d ,d (x) -f (x)dx + X\S |f α β d ,d (x) -f (x)|dx ≤ vol(S ) diam(Y ) + vol(X) f α β d ,d -f (x) L ∞ (X\S ) ≤ diam(Y ) + vol(X) f α β d ,d -f (x) L ∞ (X\S ) .
By uniform convergence the second term goes to 0 as d → ∞, this shows that lim sup

d→∞ f α β d ,d -f L 1 (X) ≤ diam(Y ).
Moreover, since > 0 was arbitrary, the limit is 0. 2

We now turn to the proof of convergence in L 1 norm for univariate functions of bounded variation.

Lemma 8. Let f : R → R be such that V (f ) is finite and f vanishes outside a segment I. Let a, b > 0 be positive constants. Let

J = {t ∈ I : ∃u, |t -u| ≤ b, u ∈ I, |f (u) -f (t)| > a} . Then J dx ≤ 2V (f ) b a .
Proof of Lemma 8: This is a packing argument. Let J 0 = J and follow the recursive process,

for k ∈ N, k ≥ 1,          if J k-1 ∩ J = ∅ let t k ∈ J k-1 ∩ J, u k ∈ J, |t k -u k | ≤ b let J k = J k-1 \ [t k -b, t k + b] otherwise stop.
This process must stop after a finite number of iterations. Indeed, the set

I k = I \ ∪ k i=1 [t k - b, t k + b]
consists of a finite union of intervals. At iteration k, either one of these intervals is a subset of [t k+1 -b, t k+1 +b] and then it is removed entirely from I k , or otherwise t k+1 is contained in an interval which contains either [t k+1 -b, t k+1 ] or [t k+1 , t k+1 + b] and the Lebesgue measure of I k+1 is reduced by at least b compared to I k , possibly creating a new interval.

Let K be the last iteration, so that J K ∩ J = ∅. By the iterative process, at each step, the measure of J k is reduced by at most 2b and we have 0 =

J K ∩J dx ≥ J K-1 ∩J dx -2b ≥ . . . ≥ J0∩J dx -2Kb so that J dx ≤ 2Kb.
Finally, since the intervals [t k , u k ], k = 1, . . . , K are disjoint, we have

V (f ) ≥ K i=1 |f (t k ) -f (u k )| ≥ Ka ≥ a 2b J dx,
which proves the desired result. 2

Proof of Theorem 4: For any d > 1 consider

J d = x ∈ X : ∃u ∈ X, |x -u| ≤ 2δ 0 /( √ d -1), |f (u) -f (t)| > d -1/4 .
By Lemma 8, for any d > 1 we have

vol(J d ) ≤ 4δ 0 d 1/4 V (f ) √ d -1 . (46)
Choose any d > 1 and any x such that x ∈ I d and x ∈ J d . Consider

(u d , v d ) ∈ argmin u∈X,v∈Y x, f α β d ,d (x) -(u, v) s.t. (u, v) ∈ supp(µ) (47)
with an arbitrary choice when the argmin is not unique. Note that the support of µ is actually the closure of the graph of f so that the minimum is attained. Using the definition of

I d in (30), x ∈ I c d implies that x, f α β d ,d (x) ∈ S d (48)
and Theorem 5 implies that

|f α β d ,d (x) -v d | ≤ δ 0 √ d -1 , |x -u d | ≤ δ 0 √ d -1 .
Since (u d , v d ) ∈ supp(µ) and supp(µ) is the closure of the graph of f , there exists h d ∈ X, such that

|f α β d ,d (x) -f (h d )| ≤ 2δ 0 √ d -1 , |x -h d | ≤ 2δ 0 √ d -1 . ( 49
) Now since x ∈ J d and |x -h d | ≤ 2δ0 √ d-1 , we have |f (x) -f (h d )| ≤ d -1/4 . This entails |f α β d ,d (x) -f (x)| ≤ |f α β d ,d (x) -f (h d )| + |f (x) -f (h d )| ≤ 2δ 0 √ d -1 + d -1/4 .
The latter expression does not depend on x which was arbitrarily chosen outside of I d and J d . We deduce that

f -f α β d ,d L 1 (X) ≤ vol(X) 2δ 0 √ d -1 + d -1/4 + diam(Y ) (vol(I d ) + vol(J d ))
and the result follows by invoking Lemma 7 and using Inequality (46). 2

5 Numerical examples

Computational tractability

Working with a large class of approximation functions may pose computational difficulties. An advantage of our Christoffel-Darboux semi-algebraic approximant is that it can be computed efficiently.

If the input moments are exactly known and given in rational form, the Christoffel-Darboux polynomial to be partially minimized is obtained through formal inversion of the moment matrix. This operation has efficient implementations, namely polynomial time algorithms over rational entries, an example is given in [START_REF] Bareiss | Sylvester's identity and multistep integer-preserving gaussian elimination[END_REF].

In most of the applications, the moments are however known only approximately, and the Christoffel-Darboux polynomial is constructed via the numerical eigenvalue decomposition of the approximate moment matrix. Since the moment matrix is symmetric, its eigenvalue decomposition can be computed efficiently with numerically stable algorithms in floating point arithmetic [START_REF] Golub | Matrix Computations. 3rd Edition[END_REF].

In addition, the computational overhead of evaluating the semi-algebraic approximant at a given point x is that of minimizing a univariate polynomial over the segment [-1, 1]. The Lipschitz constant of a univariate polynomial with coefficients p = (p 0 , . . . , p 2d ) over [-1, 1] is at most p 1 and hence grid search finds an -accurate solution to this problem using 2 p 1 evaluations. In our case the entries of p are polynomials in x which are deduced from moment data so that for a fixed d estimating and evaluating our semi-algebraic approximant up to a fixed arbitrary precision with rational inputs (moment matrix and x) can be done in polynomial time. Note also that our analysis shows that a level of precision of the order d -p-2 is sufficient so that the cost of the overall procedure has a complexity which is polynomial in the bit size of the moment matrix, the target evaluation point x, as well as in d, the degree bound.

Prototype code

We provide a simple Matlab prototype to validate our algorithm. All the examples described below are reproducible, and the Matlab scripts can be found at homepages.laas.fr/henrion/software/momgraph

The calling syntax of the main routine is

Sign function

Consider the measure supported on the graph of the sign function f (x) = sign(x) whose moments in the monomial basis on X := [-1, 1] are

1 -1 x a1 f (x) a2 dx = (-1) a2 0 -1 x a1 dx+(1) a2 1 0 x a1 dx = (-1) a2 (0 a1+1 -(-1) a1+1 ) + 1 -0 a1+1 a 1 + 1 for (a 1 , a 2 ) ∈ N 2 .
Here is an example of the use of momgraph to recover the sign function from the (floating point approximations) of the (exact) moments: [START_REF] Eckhoff | Accurate and efficient reconstruction of discontinuous functions from truncated series expansions[END_REF].

>> M %
On Figure 4 we see that the resulting degree 4 semi-algebraic approximation cannot be distinguished from the graph of the sign function. Our semi-algebraic approximant is x → argmin y∈[-1,1] q(x, y) with q the Christoffel-Darboux polynomial constructed as the sum of squares of the polynomials returned by the momgraph function: We see in particular that the first polynomial is (x, y) → 1 -y 2 with a large scaling factor. This polynomial vanishes on the graphs of the functions y → -1 and y → 1. The other polynomials are instrumental to determining which one of the two graphs corresponds to a given value of x.

Discontinuous functions

Let us revisit the discontinuous univariate examples of [START_REF] Eckhoff | Accurate and efficient reconstruction of discontinuous functions from truncated series expansions[END_REF]. Since in this case we do not have the analytic moments of the measure supported on the graph of the function f to input to our algorithm, we use the empirical moment matrix computed by uniform sampling, i.e.

M = 1 N N k=1 b(x k , f (x k ))b(x k , f (x k )) (50)
for N sufficiently large, i.e. 10 3 , and b the monomial basis vector. Degree 10 semi-algebraic approximations are reported on Figure 5 for three benchmarks [START_REF] Eckhoff | Accurate and efficient reconstruction of discontinuous functions from truncated series expansions[END_REF]Examples 65,66, 67] of discontinuous functions f , appropriately scaled in X = Y = [-1, 1]. We observe that the second rightmost discontinuity in the middle example is not detected. Increasing the degree of the approximations does not fix the issue, and we believe that it is due to the poor resolution of the monomial basis. It would be more appropriate to use here a complex exponential basis (i.e. Fourier coefficients) or an orthogonal basis (e.g. Chebyshev or Legendre polynomials).

Interpolation

Suppose now that we have access only to the values {f (x k )} k=1,...,N of the function to be approximated at given sampling points {x k } k=1,...,N , for N small. Our algorithm takes as input the empirical moment matrix (50). On Figure 6, we revisit [11, Example 65] to study the effect of the number of samples N on the quality of the approximation, for a uniform distribution of samples. We see that with 20 samples the function is already well approximated.

Recovering trajectories for optimal control

In [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF], the moment-SOS hierarchy is applied to solve numerically non-linear optimal control of ODEs with polynomial data and semi-algebraic state and control constraints. Non-linear optimal control is formulated as a linear problem on moments of occupation measures supported on optimal trajectories. Let us show how numerical approximations of these moments obtained by semidefinite programming can be input to our algorithm to approximate optimal state and control trajectories.

Let us revisit the state-constrained double integrator problem of [23, Section 5.1] to approximate the time optimal trajectories. After a scaling of time, state and control, we use the Matlab interface GloptiPoly 3 and the conic solver MOSEK to compute the pseudo-moments of the occupation measure of degree up to 8. This can be achieved in less than 2 seconds on a standard desktop computer. From this output, we construct the 45-by-45 moment matrices of the control and state marginals, conditioned w.r.t. time. Using our notations, the independent variable x is time, while the dependent variable y is respectively the control, the first state and the second state. For this example, the analytic trajectories are available for comparison. We see on Figure 7 that the state trajectory approximations are tight, whereas the control trajectory approximation misses partly the central region corresponding to the saturation of the second state. Indeed, since it is obtained by solving numerically the degree 8 semidefinite relaxation of the moment-SOS hierarchy, the approximated moment matrix differs from the exact moment matrix, and this has an impact on the quality of the Christoffel-Darboux approximation. For this example, we can construct analytically the exact moment matrix of the control trajectory and observe that indeed its Christoffel-Darboux semi-algebraic approximation of degree 8 identifies well the optimal control trajectory switching times, see Figure 8.

Bivariate examples

Consider the indicator function

f (x) := I {x∈R 2 :x 2 1 +x 2 2 ≤1/4} (x)
of a centered disk of radius 1/2. We compute the emprical moments obtained by sampling 100 2 points on a uniform grid of X := [-1, 1] 2 . With this input, our algorithm computes the degree 8 semi-algebraic approximation reported on Figure 9, to be compared with the Chebyshev polynomial approximation obtained from 100 2 points by the chebfun2 command, showing the typical Gibbs phenomenon.

We perform the same computations for the piecewise constant function f (x) := I {x∈R 2 :x 2 1 +x 2 2 ≤1/4} (x) -

1 2 I {x∈R 2 :(x1+ 1 2 ) 2 +(x2+ 1
2 ) 2 ≤1/4} (x)

Discontinuous solutions of non-linear PDEs

In [START_REF] Marx | A moment approach for entropy solutions to nonlinear hyperbolic PDEs[END_REF], the moment-SOS hierarchy is applied to solve numerically a class of non-linear PDEs for which we known that classical (i.e. differentiable) solutions do not exist. The advantage of optimizing over occupation measures is that they can be supported on graphs of weak (i.e. possibly discontinuous) solutions. Let us show how approximate moments of these measures computed by semidefinite programming can be processed by our algorithm so as to recover these discontinuous solutions. We focus on the Burgers equation and choose the initial data (a function of one space coordinate, at time zero) in a way that at a given time a shock appears, i.e. the solution becomes a discontinuous function of the space coordinate. Once the shock appeared, it propagates through, i.e. the discontinuity remains but its location varies. In Figures 12,13, and 14 we show the graphs obtained from the moment relaxations proposed in [START_REF] Marx | A moment approach for entropy solutions to nonlinear hyperbolic PDEs[END_REF]. In all cases we use the 969 triviate moments of degree 16 of the occupation measure (supported on time, space, and solution) to recover the graph of the approximated solution. For comparison we also sketch the analytic solution with red lines.

For the graphs in Figures 12 and13, the approximated moments match the analytic moments up to an error of the order of 10 -8 . Our semi-algebraic approximations are almost identical to the analytic solution.

For the graph in Figure 14 the approximated moments are noticeably incorrect, i.e., the error is of order 10 -4 . Nevertheless, our semi-algebraic approximation is able to reproduce the graph of the solution quite accurately. In particular the propagation of the shock is retrieved from the moment data. However, the approximation is erroneous when the solution passes over from its continuous to its discontinuous part.

Conclusion

In this paper, we describe a new technique to estimate discontinuous functions from moment data, based on Christoffel-Darboux kernels. Instead of using polynomial or piecewise polynomial approximants, we use a class of semi-algebraic approximants, namely arguments of minima of polynomials. This is another occurrence of a lifting technique: instead of using only moments depending linearly on the function so as to recover directly the function, we use also moments depending non-linearly on the function so as to approximate the support of a measure concentrated on the graph of the function. We provide functional analytic and geometric convergence proofs. Finally, some numerical examples illustrate the efficiency of our algorithm. We believe that this work opens the way to many other further research lines:

• When applying the Moment-SOS hierarchy, the moments are numerical approximations of the real ones. It would be interesting to provide a sensibility analysis of the application of our algorithm for the real moments and the approximated ones. We believe that such an analysis can be performed, since promising results were achieved recently in [START_REF] Klep | Minimizer extraction in polynomial optimization is robust[END_REF] for the case of zero dimensional manifolds, i.e. unions of finitely many points.

• It could also be interesting to investigate in a more quantitative way why the Gibbs phenomenon might be avoided or at least attenuated with the technique we provide. We believe that it is mainly due to the semi-algebraic point of view we are following.

• We could also check whether our algorithm works as well when considering only the knowledge of Fourier coefficients, namely moments depending linearly on the function that we want to approximate. In many problems, this is the only measurement that we might have. This is therefore a partial moment information, and we may want to complement it with estimates of higher degree moments. This makes the problem challenging.

Acknowledgments

We are grateful to Quentin Vila for his technical input on discontinuous solutions of PDEs, and to Milan Korda, Victor Magron and Matteo Tacchi for interesting discussions. This work was partly funded by the ERC Advanced Grant Taming and was also conducted in the framework of the regional programme "Atlanstic 2020, Research, Education and Innovation in Pays de la Loire, supported by the French Region Pays de la Loire and the European Regional Development

Figure 1 :

 1 Figure 1: SOS polynomial q(x, y) whose argmin in y is the sign of x on [-1, 1].

Figure 2 :

 2 Figure 2: On the left, Chebyshev interpolant of degrees 4 (gray) and 20 (black) of the step function (red), featuring the typical Gibbs phenomenon. On the right, the proposed semi-algebraic approximation of degree 4 (black) of the same step function (red). The approximation cannot be distinguished from the step function.

Theorem 2 .

 2 Under Assumption 1, let X and Y be intervals of the real line, and assume that V (f ) < +∞. With the choice β d = 2 3- √ d in Definition 2, for any r > 2 and for any d > 1, it holds

Theorem 3 .

 3 Under Assumptions 1 and 2, and with the choice α ∈ [0, 1) and β d = 2 3- √ d in Definition 4, it holds:

Theorem 4 .

 4 Under Assumptions 1 and 2, let X and Y be intervals of the real line and assume that V (f ) < +∞. With the choice α ∈ [0, 1) and β d = 2 3- √ d in Definition 4, for any r > 2 and for any d > 1, it holds

  g. [25, Remark 3.6.] and [20, Equation (1.1.)].

√ d- 1 . 2 Lemma 5 .

 125 This proves the desired result. Let us give two additional simple technical lemmas. For any r > 0, min x>0 {log(2)x -(2r) log(x)} = (2r) 1 -log 2r log(2)

1 (

 1 [Y,P] = momgraph(M,X)It takes as an input an approximate moment matrixM = X b(x, f (x))b(x, f (x)) dx(in Matlab floating point format) for b the monomial basis vector (in grevlex ordering), and a collection of points{x 1 , x 2 , . . . , x N } ⊂ X ⊂ [-1, 1] p-in Matlab floating point format) with p > 1.It outputs an approximation{y 1 , y 2 , . . . , y N } ⊂ Y := [-1, 1](in Matlab floating point format) of the values {f (x 1 ), f (x 2 ), . . . , f (x N )}, as well as a matrix P of coefficients (in Matlab floating point format) of the vector of polynomials p whose sum of squares yields the Christoffel-Darboux polynomial. Our implementation is straightforward, not optimized for efficiency. The regularization parameter β is set to the default value of 10 -8 , and the Christoffel-Darboux polynomial is computed from the eigenvalue decomposition (Matlab's command eig) of the approximate moment matrix M + βI.

Figure 4 :

 4 Figure 4: Graph of the sign function (red) and its degree 2 (left) and degree 4 (right) semi-algebraic approximations (black).

Figure 5 :

 5 Figure 5: Degree 10 semi-algebraic approximations (black) for the discontinuous univariate functions (red) of Examples 65 (left), 66 (middle) and 67 (right) of reference [11].

  >> mpol x y; b = mmon([x y],2); % GloptiPoly monomial vector of degree 2 >> P*b 6-by-1 polynomial vector (1,1):7071.0678-7071.0678y^2 (2,1):0.86713+9.4x^2-9.7305xy+0.86713y^2 (3,1):2.4315x-1.3013y (4,1):-0.53517+1.2757x^2+1.137xy-0.53517y^2 (5,1):0.29635x+0.55374y (6,1):0.29443+0.10761x^2+0.15643xy+0.29443y^2

Figure 6 :

 6 Figure 6: Degree 10 semi-algebraic approximations (black) of a discontinuous function (red) computed from empirical moments evaluated at 10 (upper left), 20 (upper right), 30 (lower left) and 40 (lower right) uniformly distributed samples.

Figure 7 :

 7 Figure 7: Minimum time double integrator with state constraint: control (left), first state (middle) and second state (right) trajectories (red) and their degree 8 semi-algebraic approximations (black) constructed from the pseudo-moments of the occupation measure.

Figure 8 :

 8 Figure 8: Minimum time double integrator with state constraint: control trajectories (red) and their degree 8 semi-algebraic approximations (black) constructed from the pseudo-moments (left, same as left of Figure 7) and from the analytic moments (right) of the occupation measure.

Figure 9 :

 9 Figure 9: Degree 4 (left) semi-algebraic approximation, and Chebyshev polynomial approximation (right) of the indicator function of a disk.

Figure 10 :

 10 Figure 10: Degree 8 (left) and degree 16 (right) semi-algebraic approximations of the superposition of signed indicator functions of two disks.

Figure 11 :

 11 Figure 11: Contour plots of the absolute error between the two disk indicator function and its degree 8 (left) and degree 16 (right) semi-algebraic approximations of Figure 10.

Figure 12 :

 12 Figure 12: Graph of the solution (a function of time and space) recovered from approximate moments for the Burgers PDE: Discontinuous initial data. The shock propagates linearly with time.

Figure 13 :

 13 Figure 13: Graph of the solution (a function of time and space) recovered from approximate moments for the Burgers PDE: Initial condition chosen to produce a shock at final time.

Figure 14 :

 14 Figure 14: Graph of the solution (a function of time and space) recovered from approximate moments for the Burgers PDE: Initial function chosen such that the shock occurs at t = 1 2 .

A semi-algebraic function is a function whose graph is semi-algebraic, i.e. defined with finitely many polynomial inequalities.

Fund. E. Pauwels and J.B. Lasserre are also partially supported by the AI Interdisciplinary Institute ANITI funding through the french program "Investing for the Future PIA3", under the Grant agreement number ANR-19-PI3A-0004.