Semi-algebraic approximation using Christoffel-Darboux kernel

TL;DR

This paper introduces a semi-algebraic approximation method using Christoffel-Darboux kernels that effectively handles discontinuous functions, with proven convergence properties and practical advantages over traditional piecewise polynomial methods.

Contribution

The paper presents a novel semi-algebraic approximation technique based on Christoffel-Darboux kernels, capable of approximating non-smooth functions without prior smoothness assumptions.

Findings

Converges pointwise almost everywhere

Achieves convergence in Lebesgue norm

Does not suffer from Gibbs phenomenon in discontinuous cases

Abstract

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.

Equations

$$\begin{array}[]{rclcl}f&:&X&\to&Y\\ &&\mathbf{x}:=(x_{1},\>x_{2},\ldots,x_{p-1})&\mapsto&y\end{array}$$

$$\mathbf{b}(\mathbf{x},y):(\mathbf{x},y)\in {\mathbb{R}}^p\mapsto \left(b_1(\mathbf{x},y)b_2(\mathbf{x},y)\cdots b_{n_d}(\mathbf{x},y)\right)\in {\mathbb{R}}^{n_d},$$

$${\int }_X\mathbf{b}(\mathbf{x},f(\mathbf{x}))\mathbf{b}(\mathbf{x},f(\mathbf{x}){)}^{\top }d\mathbf{x}$$

$$\mathbf{z}:=(\mathbf{x},y)\in {\mathbb{R}}^p,$$

$${\mathbf{M}}_{\mu ,d}:=\int \mathbf{b}(\mathbf{z})\mathbf{b}(\mathbf{z}{)}^{\top }d\mu (\mathbf{z})$$

$$d\mu (\mathbf{z}):={\mathbb{I}}_X(\mathbf{x})d\mathbf{x}{\delta }_{f(\mathbf{x})}(dy)$$

$$\{(\mathbf{x},f(\mathbf{x})):\mathbf{x}\in X\}\subset X\times Y$$

$$f_d(\mathbf{x}):={\mathrm{argmin}}_{y\in Y}{q}_d(\mathbf{x},y)$$

$$p_1:{\mathbb{R}}^2$$

$$(\mathbf{x},y)$$

$$\mathrm{argmin}_{y\in Y}p_{1}(\mathbf{x},y)=\mathrm{sign}(\mathbf{x}):=\left\{\begin{array}[]{ll}-1&\text{if}\quad\mathbf{x}<0\\ \{-1,1\}&\text{if}\quad\mathbf{x}=0\\ 1&\text{if}\quad\mathbf{x}>0\\ \end{array}\right.$$

$${\mathbf{M}}_{\mu ,d}=\int \mathbf{b}(\mathbf{z})\mathbf{b}(\mathbf{z}{)}^{\top }d\mu (\mathbf{z})\in {\mathbb{S}}^{n_d}.$$

$${\mathbf{M}}_{\mu ,d}={\mathbf{PEP}}^{\top },$$

$$\begin{array}[]{llll}\mathbf{p}_{i}^{\top}\mathbf{M}_{\mu,d}\,\mathbf{p}_{i}&=&e_{i}&=\int p_{i}^{2}(\mathbf{z})d\mu(\mathbf{z}),\\ \mathbf{p}_{i}^{\top}\mathbf{M}_{\mu,d}\,\mathbf{p}_{j}&=&0&=\int p_{i}(\mathbf{z})p_{j}(\mathbf{z})d\mu(\mathbf{z}),\quad\,i\neq j.\end{array}$$

$$q_{\mu ,d}(\mathbf{z}):=\sum _{i=1}^{n_d}\frac{p_i^2(\mathbf{z})}{e_i}=\mathbf{b}(\mathbf{z}{)}^{\top }{\mathbf{M}}_{\mu ,d}^{-1}\mathbf{b}(\mathbf{z}).$$

$$q^{e}_{\mu,d}(\mathbf{z}):=\left\{\begin{array}[]{ll}+\infty&\text{ if }\exists i,\,e_{i}=0,\,p_{i}(\mathbf{z})\neq 0\\ \sum_{e_{i}>0}\frac{p^{2}_{i}(\mathbf{z})}{e_{i}}=\mathbf{b}(\mathbf{z})^{\top}\mathbf{M}_{\mu,d}^{\dagger}\mathbf{b}(\mathbf{z})&\text{ otherwise},\end{array}\right.$$

$$q_{\mu ,d}^e(\mathbf{z})=\sum _{i=1}^{n_d}g(e_i)p_i^2(\mathbf{z})$$

$$\sum _{i=1}^{n_d}{g}_{\beta }(e_i)p_i^2(\mathbf{z})$$

$$g_{\beta }:s\mapsto \frac{1}{\beta +s},$$

$$g_{\beta }:s\mapsto \frac{1}{\beta }{\mathbb{I}}_{[0,\beta ]}(s)+\frac{1}{s}{\mathbb{I}}_{(\beta ,+\infty )}(s),$$

$$g_{\beta }:s\mapsto \frac{1}{\beta }{\mathbb{I}}_{[0,\beta ]}(s).$$

$$\mathbf{z}\mapsto \sum _{i=1}^{n_d}\frac{p_i^2(\mathbf{z})}{e_i+\beta }=\mathbf{b}(\mathbf{z}{)}^{\top }({\mathbf{M}}_{\mu ,d}+\beta {\mathbf{I}}_{n_d}{)}^{-1}\mathbf{b}(\mathbf{z})$$

$${\mathbf{M}}_{\mu ,d}+\beta {\mathbf{I}}_{n_d}={\mathbf{M}}_{\mu +\beta {\mu }_0,d}.$$

$$q_{\mu +\beta {\mu }_0,d}(\mathbf{z}):=\sum _{i=1}^{n_d}\frac{p_i^2(\mathbf{z})}{e_i+\beta }.$$

$$\mathbf{x}\in X\mapsto f_{\beta ,d}(\mathbf{x}):=\min \{{\mathrm{argmin}}_{y\in Y}{q}_{\mu +\beta {\mu }_0,d}(\mathbf{x},y)\}.$$

$${\mathrm{argmin}}_{y\in Y}{q}_{\mu +\beta {\mu }_0,d}(\mathbf{x},y)$$

$${\delta }_0:=\mathrm{diam}(\mathrm{supp}(\mu +{\mu }_0)),m:=\mu ({\mathbb{R}}^p),m_0:={\mu }_0({\mathbb{R}}^p).$$

$$f_{{\beta }_d,d}(\mathbf{x})\underset{d\to \infty }{\to }f(\mathbf{x})$$

$$\Vert f-f_{{\beta }_d,d}{\Vert }_{{\mathcal{L}}^1(X)}\underset{d\to \infty }{\to }0.$$

$$\Vert f-f_{{\beta }_d,d}{\Vert }_{{\mathcal{L}}^1(X)}$$

$$\le$$

$$V(f)=\underset{n\in \mathbb{N}}{\sup }\underset{x_0<x_1<x_2<\dots <x_n}{\sup }\sum _{i=1}^n|f(x_i)-f(x_{i-1})|.$$

Full text

Keywords: approximation theory, convex optimization, moments, positive polynomials, orthogonal polynomials.

2020 MSC: 42C05, 47B32, 41A30

1. 1 Introduction

2. 1.1 Motivation

3. 1.2 Contribution

4. 1.3 Organization of the paper

5. 2 Christoffel-Darboux approximation

6. 2.1 Polynomials and moments

7. 2.2 Approximating the support from moments

8. 2.3 Regularization scheme

9. 2.4 Semi-algebraic approximant

10. 3 Main results

11. 3.1 Convergence under continuity assumptions

12. 3.2 Convergence for univariate functions of bounded variation

13. 3.3 Robustness to small perturbations

14. 4 Proofs

15. 4.1 Estimation of the support

16. 4.2 Estimation of functions

17. 4.3 Proofs of the main theorems

18. 5 Numerical examples

19. 5.1 Computational tractability

20. 5.2 Prototype code

21. 5.3 Sign function

22. 5.4 Discontinuous functions

23. 5.5 Interpolation

24. 5.6 Recovering trajectories for optimal control

25. 5.7 Bivariate examples

26. 5.8 Discontinuous solutions of non-linear PDEs

27. 6 Conclusion

28. 7 Acknowledgments

1 Introduction

In this paper we address the following generic inverse problem. Let

[TABLE]

be a bounded measurable function from a given compact set $X\subset\mathbb{R}^{p-1}$ to a given compact set $Y\subset\mathbb{R}$, with $p\geq 2$. We assume that $X$ is equal to the closure of its interior.

Given $d\in\mathbb{N}$, consider a vector of polynomials of total degree at most $d$

[TABLE]

where $n_{d}:=\begin{pmatrix}p+d\\ d\end{pmatrix}$ which is understood as the binomial coefficient. For example, $\mathbf{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 $\mathbf{b}$ and $f$, let

[TABLE]

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.

1.1 Motivation

Inverse Problem 1 is encountered in several interesting situations. In the weak (or measure-valued) formulation of Optimal Control problems (OCP) [39, 17, 23], Markov Decision Processes [18], option pricing in finance [24], stochastic control and optimal stopping [37], and some non-linear partial differential equations (PDEs) [8], 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

[TABLE]

the moment matrix of degree $2d$ reads

[TABLE]

corresponding to the measure

[TABLE]

supported on the graph

[TABLE]

of function $f$, where $\mathbb{I}_{X}$ denotes the indicator function of $X$ which takes value $1$ on $X$ and [math] otherwise, and $\delta_{f(\mathbf{x})}$ denotes the Dirac measure at $f(\mathbf{x})$.

For instance, in OCPs an optimal occupation measure $\mu$ 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 $\mu$ with respect to time-state resp. time-control. Then with our notation, $\mathbf{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, $\mathbf{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 $\mu$ can be disintegrated into its marginal on $X$ and its conditional on $Y$ given $\mathbf{x}$ in $X$. The latter is also called a parametrized measure or a Young measure, see e.g. [12].

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 [16] and the Moment-SOS (polynomial sums of squares) methodology described in e.g. [23, 21, 22]. In the latter methodology one may thus approximate the optimal solution $\mu$ of the measure-valued formulation by solving a hierarchy of semidefinite relaxations of the problem, whose size increases with $d$; see e.g. [23] for OCPs and [26] 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 $\mu$. 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 $\mu$, 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 [22] to solve Generalized Moment Problems where the involved Borel measures are Young measures. The necessary moment information is also given when considering empirical measures [31, 25, 32] 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(\mathbf{x})\>d\mathbf{x}$, i.e. $y\mapsto\mathbf{b}(\mathbf{x},y)$ is linear. Finally, our method and results would apply to measures $\nu$ supported on the graph of $f$, different from $\mu$ in (1). Provided that $\nu$ and $\mu$ are mutually absolutely continuous with bounded densities, we would recover similar convergence results modulo constants.

1.2 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 [38], [10] and approximation theory [28], [6]. 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 [27], [15] for single dimensional problem and in [20] in a multivariate setting. Even more recently, in [31], [25] and [32], 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 $\mu$. 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}(\mathbf{x},y)$ of degree $2d$ in $p$ variables. For every fixed $\mathbf{x}$ we define

[TABLE]

which is a semi-algebraic function111A semi-algebraic function is a function whose graph is semi-algebraic, i.e. defined with finitely many polynomial inequalities., 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 $\mathbf{x}$ of degree $d$, then $q:(\mathbf{x},y)\mapsto(r(\mathbf{x})-y)^{2}$ is a degree $2d$ SOS polynomial whose partial minimization in $y$ yields $y=r(\mathbf{x})$ for all $\mathbf{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 non-linear PDEs.

Example 1** (Sign function as SOS partial minimum).**

Consider the polynomial

[TABLE]

One can easily check that

[TABLE]

for any $\mathbf{x}\in X:=[-1,1]$ and $Y\subset\mathbb{R}$. Note that since $p_{1}(\mathbf{x},\cdot)$ is positive for $\mathbf{x}\in[-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(\mathbf{x},y):=p_{1}^{2}(\mathbf{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-12\mathbf{x}^{4}y-6\mathbf{x}^{2}y^{2}+4\mathbf{x}^{2}y^{3}+3y^{4}$ is equal to $|\mathbf{x}|$ for all $\mathbf{x}\in[-1,1]$.

To the best of our knowledge, this work is the first contribution where this class of semi-algebraic 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 $\mathscr{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(\mathbf{x})\>d\mathbf{x}$, i.e. $y\mapsto\mathbf{b}(\mathbf{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 $\mu$ in (1). 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. [14] for a good survey on this topic. The main tools usually rely on properties of orthogonal polynomials [38] 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 [11] 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 [3, 4]. 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. [30] 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

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 [25] 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\in\mathbb{N}$, we provide an approximation $f_{d}$ of the function $f$ and prove that the sequence $(f_{d})_{d\in\mathbb{N}}$ converges pointwise to $f$ almost everywhere as well as in $\mathscr{L}^{1}$-norm 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 $\{(\mathbf{x},f(\mathbf{x})):\mathbf{x}\in X\}\subset X\times Y$ of the function $\mathbf{x}\mapsto f(\mathbf{x}):=\mathrm{sign}(\mathbf{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 [9]) 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. [30] 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.

1.3 Organization of the paper

Section 2 introduces the Christoffel-Darboux polynomial, its regularized version and our semi-algebraic 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 $\mu$ and then prove the $\mathscr{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 $\mathbb{S}^{n}$. Given a set $X$ in Euclidean space, its diameter is denoted by $\mathrm{diam}(X)$ and its volume, or Lebesgue measure, is denoted by $\mathrm{vol}(X)$. For $k\geq 1$, the Lebesgue space $\mathscr{L}^{k}(X)$ consists of functions on $X$ whose $k$-norms are bounded. Given a positive Borel measure $\mu$, we denote by $\mathrm{supp}(\mu)$ 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 $\mathbf{z}=(\mathbf{x}\>\>y)\in\mathbb{R}^{p}$. We denote by $\mathbb{R}[\mathbf{z}]$ the algebra of multivariate polynomials of $\mathbf{z}\in\mathbb{R}^{p}$ with real coefficients. For a given degree $d\in\mathbb{N}$, the dimension of the vector space of polynomials of degree less than or equal to $d$ is denoted by $n_{d}:=\begin{pmatrix}p+d\\ d\end{pmatrix}$.

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.

2.1 Polynomials and moments

Following the notations introduced for stating Problem 1, any polynomial $q\in\mathbb{R}[\mathbf{z}]$ of degree $d$ can be expressed in the polynomial basis $\mathbf{b}(\mathbf{z})$ as $q\colon\mathbf{z}\mapsto\mathbf{q}^{\top}\mathbf{b}(\mathbf{z})$ with $\mathbf{q}\in\mathbb{R}^{n_{d}}$ denoting its vector of coefficients. Recall that the moment matrix of degree $2d$ of the measure $\mu$ reads

[TABLE]

Since $\mathbf{M}_{\mu,d}$ is positive semi-definite, it has a spectral decomposition

[TABLE]

where $\mathbf{P}\in\mathbb{R}^{n_{d}\times n_{d}}$ is an orthonormal matrix whose columns are denoted $\mathbf{p}_{i}$, $i=1,2,\ldots,n_{d}$, and satisfy $\mathbf{p}_{i}^{\top}\mathbf{p}_{i}=1$ and $\mathbf{p}_{i}^{\top}\mathbf{p}_{j}=0$ if $i\neq j$, and $\mathbf{E}\in\mathbb{S}^{n(d)}$ is a diagonal matrix whose diagonal entries are eigenvalues $e_{i+1}\geq e_{i}\geq 0$ of the moment matrix. Each column $\mathbf{p}_{i}\in\mathbb{R}^{n_{d}}$ is the vector of coefficients of a polynomial $p_{i}\in\mathbb{R}[\mathbf{z}]$, $i=1,\ldots n_{d}$, so that

[TABLE]

2.2 Approximating the support from moments

Let us assume for now that the support of the measure $\mu$ has nonempty interior, then $\mathbf{M}_{\mu,d}$ is positive definite for any $d\in\mathbb{N}$, i.e., $e_{i}>0$, $i=1,\ldots,n_{d}$. In this case, one can define the Christoffel-Darboux polynomial

[TABLE]

It is known that sublevel sets of $q_{\mu,d}$ can be used to recover the support of $\mu$ for large $d$, see for example [25] 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 $\mu$ in (1). Hence the sublevel sets of $q_{\mu,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 $\mathbf{M}_{\mu,d}$ may not be positive definite and invertible. In this singular setting, one should ideally consider the following extended value Christoffel-Darboux polynomial:

[TABLE]

where $\dagger$ denotes the Moore-Penrose pseudo inverse. This is a natural extension of the Christoffel-Darboux polynomial to the singular case [32] and amounts to working in the Zariski closure of the graph of $f$.

2.3 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

[TABLE]

where $g\colon[0,+\infty)\mapsto[0,+\infty]$ with $g(s)=\frac{1}{s}$ for any $s>0$ and $g(0)=+\infty$. 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 [7] 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

[TABLE]

where $g_{\beta}$ is a parametrized family of spectral filtering regular functions satisfying, for any $\beta>0$, $g_{\beta}\colon[0,+\infty)\mapsto[0,+\infty]$. Common examples include

[TABLE]

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 (7) yields the following polynomial

[TABLE]

where $\mathbf{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 $\mathbb{R}^{p}$ that we will denote $\mu_{0}$. We make the following assumption which will be standing throughout the paper.

Assumption 1** **(Reference measure and polynomial basis).

• •

The reference measure $\mu_{0}$ is absolutely continuous with respect to the Lebesgue measure and it has compact support.

• •

The polynomial basis $\mathbf{b}$ is orthonormal with respect to the bilinear form induced by $\mu_{0}$, that is $\int b_{i}(\mathbf{z})b_{j}(\mathbf{z})d\mu_{0}(\mathbf{z})=1$ if $i=j$ and [math] otherwise.

The first part of Assumption 1 ensures that the moment matrix of $\mu_{0}$ is always positive definite. The second part of Assumption 1 provides the following relation:

[TABLE]

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 (5), 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

[TABLE]

This provides a geometric interpretation of the regularization parameter as a combination of two measures: $\mu$ which is supported on the graph of the function of interest and $\mu_{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 $\mu_{0}$;

• •

Add assumption on the spectrum of the moment matrix $\mu_{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.

2.4 Semi-algebraic approximant

Definition 2** **(Semi-algebraic approximant).

The regularized Christoffel-Darboux semi-algebraic approximant $f_{\beta,d}$ is defined as follows:

[TABLE]

Remark 1**.**

If $X$ and $Y$ are semi-algebraic, then the set-valued map which associates to each $\mathbf{x}\in X$ the set

[TABLE]

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 semi-algebraic subset of $Y$, which is itself a subset of the real line. Hence the argmin set has a minimal element and the function in (11) 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 31, approximate minimization up to a factor of the order $d^{p+2}$ enjoys similar approximation properties. Indeed, in Remark 5, we mention the precision $\frac{\gamma_{d}}{2}$, where $\gamma_{d}$ is chosen later on with more justifications in (14) to be of order $d^{p+2}$.

The two parameters $d$ and $\beta$ control the behavior of the approximant $f_{\beta,d}$. In latter sections, we describe an explicit dependency between $d$ and $\beta$ which allows to construct a sequence of regularization parameters $(\beta_{d})_{d\in\mathbb{N}}$, and we investigate the asymptotic properties of the sequence of approximants $\left(f_{\beta_{d},d}\right)_{d\in\mathbb{N}}$.

3 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_{\beta,d}$ under Assumption 1. In particular we are interested in building an explicit sequence $\left(\beta_{d}\right)_{d\in\mathbb{R}}$ and investigating the convergence $f_{\beta_{d},d}(\mathbf{x})\to f(\mathbf{x})$ for $\mathbf{x}\in X$, fixed, as well as the convergence $\|f-f_{\beta_{d},d}\|_{\mathscr{L}^{1}(X)}\to 0$, when $d\to\infty$.

3.1 Convergence under continuity assumptions

The following section describes our main result regarding convergence of the approximant $f_{\beta_{d},d}$ in (11) under different regularity assumptions for the function $f$ to be approximated. Let us define

[TABLE]

Theorem 1**.**

Under Assumption 1 and with the choice $\beta_{d}=2^{3-\sqrt{d}}$ in Definition 2, it holds:

• (i)

If the set $S\subset X$ of continuity points of $f$ is such that $X\setminus S$ has Lebesgue measure zero, then

[TABLE]

for almost all $\mathbf{x}\in X$, and

[TABLE]

• (ii)

If $f$ is $L$-Lipschitz on $X$ for some $L>0$, then for any $d>1$ and any $r>p$,

[TABLE]

Remark 3**.**

• •

As proved recently in **[35]**, the solutions to scalar conservation laws are continuous almost everywhere, i.e. the Lebesgue measure of the set of discontinuity points reduces to [math]. 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 $\mathscr{L}^{1}$.

• •

For Lipschitz continuous functions, we obtain an $O(d^{-1/2})$ convergence rate in $\mathscr{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.

3.2 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 [1] 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 $\mathscr{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\colon\mathbb{R}\mapsto\mathbb{R}$ be a measurable function. The (Jordan) total variation norm of $f$ is given by

[TABLE]

Theorem 2**.**

Under Assumption 1, let $X$ and $Y$ be intervals of the real line, and assume that $V(f)<+\infty$. With the choice $\beta_{d}=2^{3-\sqrt{d}}$ in Definition 2, for any $r>2$ and for any $d>1$, it holds

[TABLE]

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.

3.3 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 $\mathbf{M}_{\mu+\beta_{d}\mu_{0}}$ as in (9) but a matrix $\mathbf{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 $\mathbf{M}\in\mathbb{S}^{n_{d}}$ is positive definite and let

[TABLE]

where we used the matrix operator norm. Then the polynomial $q^{\alpha}_{d}\colon\mathbf{z}\mapsto\mathbf{b}(\mathbf{z})^{\top}\mathbf{M}^{-1}\mathbf{b}(\mathbf{z})$ satisfy

[TABLE]

Proof: For any $\mathbf{z}\in\mathbb{R}^{p}$, we have

[TABLE]

$\Box$

Remark 4** (Accuracy of the moments).**

In this paper, we assume the existence of a positive number $\alpha$ which makes the link between an approximated moment matrix and the real one, which exists in general, but estimates for $\alpha$ are not known in general. Such an analysis has been performed for the moments themselves in [29], 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 $\alpha\in[0,1)$, let $(\beta_{d})_{d\in\mathbb{N}}$ be a sequence of positive numbers and $(q^{\alpha}_{d})_{d\in\mathbb{N}}$ be a sequence of continuous functions over $\mathbb{R}^{p}$ such that for any $d\in\mathbb{N}$ and any $\mathbf{z}\in\mathbb{R}^{p}$, we have

[TABLE]

Note that Assumption 2 ensures that

[TABLE]

Furthermore, one can always choose $q^{\alpha}_{d}(\mathbf{z})=q_{\mu+\beta_{d}\mu_{0},d}(\mathbf{z})$ and $\alpha=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\in\mathbb{N}$, a regularizing parameter $\beta>0$, and a scalar $\alpha\in[0,1)$, our robust approximant $f^{\alpha}_{\beta_{d},d}$ is defined as follows:

[TABLE]

The main technical result of this paper is the following from which Theorem 1 directly follows.

Theorem 3**.**

Under Assumptions 1 and 2, and with the choice $\alpha\in[0,1)$ and $\beta_{d}=2^{3-\sqrt{d}}$ in Definition 4, it holds:

• (i)

If the set $S\subset X$ of continuity points of $f$ is such that $X\setminus S$ has Lebesgue measure zero, then

[TABLE]

for almost all $\mathbf{x}\in X$, and

[TABLE]

• (ii)

If $f$ is $L$-Lipschitz on $X$ for some $L>0$, then for any $d>1$ and any $r>p$,

[TABLE]

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

Theorem 4**.**

Under Assumptions 1 and 2, let $X$ and $Y$ be intervals of the real line and assume that $V(f)<+\infty$. With the choice $\alpha\in[0,1)$ and $\beta_{d}=2^{3-\sqrt{d}}$ in Definition 4, for any $r>2$ and for any $d>1$, it holds

[TABLE]

The next section is dedicated to the proof of these theorems.

4 Proofs

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

4.1 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\in\mathbb{N}$ and $r\in\mathbb{N}$ such that $r>p$ and for any $\alpha\in[0,1)$ , define

[TABLE]

and

[TABLE]

We aim at proving that the sublevel set $S_{d}$ is approaching the support of $\mu$ 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)

[TABLE]

• (ii)

For any $\mathbf{z}\in S_{d}$,

[TABLE]

The results and techniques that we will use to prove this theorem are adapted from [25] which considers the absolutely continuous setting without regularization.

Proof of Theorem 5 (i) Using (10) and (4), we obtain

[TABLE]

where the last inequality is given in [25, Lemma 6.5]. Using Markov’s inequality [36, Page 91] and (16) yields

[TABLE]

Now using (12) and (15), we have the following implications

[TABLE]

Hence $\mu(\{\mathbf{z}\in\mathbb{R}^{p}:\mathbf{z}\not\in S_{d}\})\leq\mu(\{\mathbf{z}\in\mathbb{R}^{p}:q_{\mu+\beta_{d}\mu_{0},d}(\mathbf{z})\geq\frac{\gamma_{d}}{1+\alpha}\})$. Using the expression of $\gamma_{d}$ in (14) and the inequality (17), one has:

[TABLE]

which concludes the proof of item $(i)$ of Theorem 5.   $\Box$

To carry out the proof of Theorem 5 (ii), we begin with a few lemmas. The following result is classical, see e.g. [25, Remark 3.6.] and [20, Equation (1.1.)].

Lemma 2**.**

Let $d\in\mathbb{N}$, $\mathbf{z}\in\mathbb{R}^{p}$, $\beta>0$, and $q$ be a polynomial of degree at most $d$. Then

[TABLE]

The following Lemma defines the needle polynomial, introduced first in [20], and gives a quantitative result crucial for our analysis.

Lemma 3** (Existence of a needle polynomial).**

Let $B_{\delta}$ denote the euclidean ball of radius $\delta$. Then, for all $\delta\in(0,1)$ and $d\in\mathbb{N}$, $d>0$, there exists a polynomial $q$ of degree $2d$ such that

[TABLE]

A detailed proof is provided in [25, Lemma 6.3]. Thanks to the latter lemma, we can characterize the behavior of the regularized Christoffel-Darboux polynomial $q_{\mu+\beta_{d}\mu_{0},d}$ outside the support of $\mu$:

Lemma 4**.**

Let $d\in\mathbb{N}$, $d>1$ and $\mathbf{z}\in\mathbb{R}^{p}$. Recall that $\beta_{d}:=2^{3-\sqrt{d}}$. Assume that $\mathrm{dist}(\mathbf{z},\mathrm{supp}(\mu))\geq\frac{\delta_{0}}{\sqrt{d}-1}$. Then

[TABLE]

Proof of Lemma 4: Let $d>1$, $\mathbf{z}\in\mathbb{R}^{p}$, $\delta=\mathrm{dist}(\mathbf{z},\mathrm{supp}(\mu))$. Let $d^{\prime}\in\mathbb{N}$ and $t>0$, arbitrary for the moment. Consider the affine map $T:\mathbf{w}\mapsto\frac{\mathbf{w}-\mathbf{z}}{\delta+\delta_{0}}$. Let $q$ be the degree $2d^{\prime}$ polynomial given as in Lemma 3 such that

[TABLE]

Let $r=q\circ T$. The polynomial $r$ satisfies

[TABLE]

Using Lemma 2 and the fact that $r(\mathbf{z})=1$ we obtain

[TABLE]

and

[TABLE]

From (21), we deduce

[TABLE]

[TABLE]

Combining (22), (23) , (24) and (25), we obtain the following bounds

[TABLE]

Recall that $d^{\prime}$ and $\beta$ were arbitrary. Now we can choose $d^{\prime}=\lfloor d/2\rfloor$, $\beta=\beta_{d}$ in one of the identities in (26) (depending on the parity of $d$) to obtain

[TABLE]

where the last inequality follows because the right hand side is strictly increasing as a function of $\delta$ and $\delta\geq\frac{\delta_{0}}{\sqrt{d}-1}$. This proves the desired result.   $\Box$

Let us give two additional simple technical lemmas.

Lemma 5**.**

For any $r>0$,

[TABLE]

Proof of Lemma 5: A simple analysis shows that the minimum is attained at $x=\frac{2r}{\log(2)}$. The lower bound follows because $\frac{2}{\log(2)}\leq 3$.   $\Box$

Lemma 6**.**

For any $d\in\mathbb{N}$, we have

[TABLE]

Proof of Lemma 6: Using the definition of $\gamma_{d}$ in (14), we have

[TABLE]

where the inequality follows from Lemma 5 with $x=\sqrt{d}$. This proves the desired result.   $\Box$

Proof of Theorem 5 (ii): We prove the result by contraposition. We have the following chain of implications for $\mathbf{z}\in\mathbb{R}^{p}$,

[TABLE]

where the first implication is from Lemma 4, the second implication is due to Lemma 6, the third implication is deduced from (12) and the last implication is from the definition of $S_{d}$ in (15).   $\Box$

4.2 Estimation of functions

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

[TABLE]

Lemma 7**.**

Suppose that $d>1$. Then, we have

[TABLE]

Proof of Lemma 31: For all $\mathbf{x}\in X$ and all $A\subset\mathbb{R}^{p}$ measurable, one has

[TABLE]

and hence

[TABLE]

where $I_{A}:=\{\mathbf{x}\in X:(\mathbf{x},f(\mathbf{x}))\in A\}$. One can see from (30) that $\mathbf{x}\in I_{d}$ implies that $(\mathbf{x},f(\mathbf{x}))\not\in S_{d}$ and hence $I_{d}\subset I_{S_{d}^{c}}$ where $S_{d}^{c}$ denotes the complement of $S_{d}$ given in (15). We deduce that

[TABLE]

and the result follows from Item (i) of Theorem 5.   $\Box$

Remark 5**.**

Letting $\tilde{I}_{d}:=\{\mathbf{x}\in X:\inf_{y\in Y}\,q^{\alpha}_{d}(\mathbf{x},y)\geq\gamma_{d}/2\}$ it can be seen using the exact same arguments that a bound on $\int_{\tilde{I}_{d}}d\mathbf{x}$ holds similarly as in Lemma 31 with a multiplicative factor of $2$. This can be used to handle the situation where the argmin in (11) is computed up to a precision of the order $\gamma_{d}/2$. See also Remark 2.

Thanks to Lemma 31, it is sufficient to prove the convergence of the approximated function in the set $I_{d}^{c}:=X\setminus I_{d}$.

Proposition 1**.**

Let $d\in\mathbb{N}$, $d>1$. Let $\alpha\in[0,1)$ be as in Assumption 2 with $\beta_{d}=2^{3-\sqrt{d}}$. Consider $f^{\alpha}_{\beta_{d},d}$ as in (13) and let $I_{d}$ be defined by (30). Then for any $r>p$, we have

[TABLE]

Proof of Proposition 1: Since $y$ takes values in the compact set $Y$, it is clear that $f_{\beta_{d},d}^{\alpha}\in\mathscr{L}^{\infty}(X)$. We have

[TABLE]

where we have used Lemma 31 for the last inequality.   $\Box$

4.3 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

[TABLE]

For any $\mathbf{x}\in X$ define

[TABLE]

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 (30), the fact that $\mathbf{x}\in I_{d}^{c}$ implies that

[TABLE]

Moreover, Theorem 5 implies that

[TABLE]

where $\|\cdot\|$ denotes the usual Euclidean distance in $\mathbb{R}^{p-1}$. Therefore, using Lipschitz continuity of $f$, we have, for any $\mathbf{x}\in I_{d}^{c}$,

[TABLE]

We deduce that

[TABLE]

which concludes the proof.   $\Box$

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 31. Choose $r=p+2$ and let

[TABLE]

Lemma 31 ensures that $\mathrm{vol}(I_{d})=O(1/d^{2})$ so that

[TABLE]

We have $\mathrm{vol}(I)=\lim_{d_{0}\to\infty}\left(\cup_{d\geq d_{0}}I_{d}\right)=0$. This means that we have the two following properties, for almost every $\mathbf{x}\in X$:

• •

$f$ is continuous at $\mathbf{x}$ (by assumption),

• •

$\exists d_{0}\in\mathbb{N}$, $\forall d\in\mathbb{N}$, $d\geq d_{0}$, $\mathbf{x}\not\in I_{d}$ (because $\mathrm{vol}(I)=0)$.

Fix any such $\mathbf{x}$ and for any $d\geq d_{0}$ consider

[TABLE]

with an arbitrary choice when the argmin is not unique. Note that the support of $\mu$ 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 $\mathbf{x}\in I_{d}^{c}$ for all $d\geq d_{0}$ which implies that

[TABLE]

and Theorem 5 implies that

[TABLE]

Since $(\mathbf{u}_{d},v_{d})\in\mathrm{supp}(\mu)$ and $\mathrm{supp}(\mu)$ is the closure of the graph of $f$, there exists $\mathbf{h}_{d}\in X$, such that

[TABLE]

This concludes the proof of pointwise convergence since

[TABLE]

and both terms tend to [math] as $d\to\infty$, using (45) and continuity of $f$ at $\mathbf{x}$.

Convergence in $\mathscr{L}^{1}$ follows from Egorov’s Theorem (see e.g. [33, chapter 18]): for any $\epsilon>0$, there exists $S_{\epsilon}\subset X$, measurable, of Lebesgue measure smaller than $\epsilon$ such that $f_{\beta_{d},d}^{\alpha}\to f$ uniformly on $X\setminus S_{\epsilon}$. We have

[TABLE]

By uniform convergence the second term goes to [math] as $d\to\infty$, this shows that

[TABLE]

Moreover, since $\epsilon>0$ was arbitrary, the limit is [math].   $\Box$

We now turn to the proof of convergence in $\mathscr{L}^{1}$ norm for univariate functions of bounded variation.

Lemma 8**.**

Let $f\colon\mathbb{R}\mapsto\mathbb{R}$ be such that $V(f)$ is finite and $f$ vanishes outside a segment $I$. Let $a,b>0$ be positive constants. Let

[TABLE]

Then

[TABLE]

Proof of Lemma 8: This is a packing argument. Let $J_{0}=J$ and follow the recursive process, for $k\in\mathbb{N}$, $k\geq 1$,

[TABLE]

This process must stop after a finite number of iterations. Indeed, the set $I_{k}=I\setminus\cup_{i=1}^{k}[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}\cap J=\emptyset$. By the iterative process, at each step, the measure of $J_{k}$ is reduced by at most $2b$ and we have

[TABLE]

so that

[TABLE]

Finally, since the intervals $[t_{k},u_{k}]$, $k=1,\ldots,K$ are disjoint, we have

[TABLE]

which proves the desired result.   $\Box$

Proof of Theorem 4: For any $d>1$ consider

[TABLE]

By Lemma 8, for any $d>1$ we have

[TABLE]

Choose any $d>1$ and any $x$ such that $x\not\in I_{d}$ and $x\not\in J_{d}$. Consider

[TABLE]

with an arbitrary choice when the argmin is not unique.

Note that the support of $\mu$ is actually the closure of the graph of $f$ so that the minimum is attained. Using the definition of $I_{d}$ in (30), $x\in I_{d}^{c}$ implies that

[TABLE]

and Theorem 5 implies that

[TABLE]

Since $(u_{d},v_{d})\in\mathrm{supp}(\mu)$ and $\mathrm{supp}(\mu)$ is the closure of the graph of $f$, there exists $h_{d}\in X$, such that

[TABLE]

Now since $x\not\in J_{d}$ and $|x-h_{d}|\leq\frac{2\delta_{0}}{\sqrt{d}-1}$, we have $|f(x)-f(h_{d})|\leq d^{-1/4}$. This entails

[TABLE]

The latter expression does not depend on $x$ which was arbitrarily chosen outside of $I_{d}$ and $J_{d}$. We deduce that

[TABLE]

and the result follows by invoking Lemma 31 and using Inequality (46).   $\Box$

5 Numerical examples

5.1 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 [2].

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 [13].

In addition, the computational overhead of evaluating the semi-algebraic approximant at a given point $\mathbf{x}$ is that of minimizing a univariate polynomial over the segment $[-1,1]$. The Lipschitz constant of a univariate polynomial with coefficients $\mathbf{p}=(p_{0},\ldots,p_{2d})$ over $[-1,1]$ is at most $\|\mathbf{p}\|_{1}$ and hence grid search finds an $\epsilon$-accurate solution to this problem using $\frac{2\|\mathbf{p}\|_{1}}{\epsilon}$ evaluations. In our case the entries of $\mathbf{p}$ are polynomials in $\mathbf{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 $\mathbf{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 $\mathbf{x}$, as well as in $d$, the degree bound.

5.2 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

[Y,P] = momgraph(M,X)

It takes as an input an approximate moment matrix

[TABLE]

(in Matlab floating point format) for $\mathbf{b}$ the monomial basis vector (in grevlex ordering), and a collection of points

[TABLE]

(in Matlab floating point format) with $p>1$. It outputs an approximation

[TABLE]

(in Matlab floating point format) of the values $\left\{f(\mathbf{x}_{1}),\,f(\mathbf{x}_{2}),\ldots,f(\mathbf{x}_{N})\right\}$, as well as a matrix $\mathbf{P}$ of coefficients (in Matlab floating point format) of the vector of polynomials $\mathbf{p}$ whose sum of squares yields the Christoffel-Darboux polynomial.

Our implementation is straightforward, not optimized for efficiency. The regularization parameter $\beta$ 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 $\mathbf{M}+\beta\mathbf{I}$.

5.3 Sign function

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

[TABLE]

for $(a_{1},\>a_{2})\in\mathbb{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:

M % moment matrix of degree 4 for the sign function

M =

2.0000         0         0    0.6667    1.0000    2.0000
     0    0.6667    1.0000         0         0         0
     0    1.0000    2.0000         0         0         0
0.6667         0         0    0.4000    0.5000    0.6667
1.0000         0         0    0.5000    0.6667    1.0000
2.0000         0         0    0.6667    1.0000    2.0000

X = linspace(-1,1,1e3)’; % samples for evaluation [Y,P] = momgraph(M,X); % Christoffel-Darboux approximation plot(X,graph(X),’-r’,’linewidth’,6); hold on; % exact graph plot(X,Y,’-k’,’linewidth’,3); xlabel(’x’); ylabel(’y’); % approximate graph

This code corresponds to a degree 4 approximation from a moment matrix of size 6 with 15 moments. A degree 2 approximation can be obtained from its 3 by 3 submatrix

M(1:3,1:3)

ans =

2.0000         0         0
     0    0.6667    1.0000
     0    1.0000    2.0000

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 $\mathbf{x}\mapsto\mathrm{argmin}_{y\in[-1,1]}q(\mathbf{x},y)$ with $q$ the Christoffel-Darboux polynomial constructed as the sum of squares of the polynomials returned by the momgraph function:

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

We see in particular that the first polynomial is $(\mathbf{x},y)\mapsto 1-y^{2}$ with a large scaling factor. This polynomial vanishes on the graphs of the functions $y\mapsto-1$ and $y\mapsto 1$. The other polynomials are instrumental to determining which one of the two graphs corresponds to a given value of $\mathbf{x}$.

5.4 Discontinuous functions

Let us revisit the discontinuous univariate examples of [11]. 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.

[TABLE]

for $N$ sufficiently large, i.e. $10^{3}$, and $\mathbf{b}$ the monomial basis vector. Degree 10 semi-algebraic approximations are reported on Figure 5 for three benchmarks [11, 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).

5.5 Interpolation

Suppose now that we have access only to the values $\{f(\mathbf{x}_{k})\}_{k=1,\ldots,N}$ of the function to be approximated at given sampling points $\{\mathbf{x}_{k}\}_{k=1,\ldots,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.

5.6 Recovering trajectories for optimal control

In [23], 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 $\mathbf{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.

5.7 Bivariate examples

Consider the indicator function

[TABLE]

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

[TABLE]

obtained as a superposition of signed indicator functions of two disks. Its degree 8 and 16 semi-algebraic approximations are reported on Figure 10. The absolute pointwise error between the approximations and the original function is displayed on Figure 11.

5.8 Discontinuous solutions of non-linear PDEs

In [26], 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 [26]. 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 and 13, 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.

6 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 [19] 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.

7 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 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.