The overall benefit associated with the sustained virologic response (SVR) has been extensively acknowledged over the past 2 decades with Interferon (IFN)-based regimens and seemed not debatable in the Hepatology community. It includes not only a reduction in the rate of occurrence of hepatocellular carcinoma (HCC), in retrospective (1-2) as well as prospective studies [START_REF] Trinchet | 12 CirVir Group. Complications and competing risks of death in compensated viral cirrhosis (ANRS CO12 CirVir prospective cohort)[END_REF], but also of extra-hepatic manifestations (4). This evidence has been markedly challenged in 2016 by the Barcelona group at the International Liver Congress (5) supported by other groups [START_REF] Conti | Early occurrence and recurrence of hepatocellular carcinoma in HCV-related cirrhosis treated with direct-acting antivirals[END_REF] suggesting that the oral direct-acting antivirals (DAA) which were massively given to "priority" HCV-infected patients, namely those with significant fibrosis, was associated with an unexpected high rate and aggressiveness of tumor recurrence in patients with apparent complete response after local treatment of HCC: among 58 patients with history of HCC treated by DAA, 3 patients died and 16 (27.6 %) had a radiological tumor recurrence after a short median follow-up of 5.7 months ( 5). This unexpectedly high rate and pattern of HCC coinciding with HCV clearance could be explained by a rapid retablishment of normal immunologic homeostasis which may facilitate the emergence of loco-regional or metastatic clones.

The quantitative and qualitative alterations of the immune repertoire which are associated with chronic HCV infection are only slowly reversible [START_REF] Alanio | Bystander hyperactivation of preimmune CD8+ T cells in chronic HCV patients[END_REF] and the restoration of the cytokines profile is only partial; the imbalance of liver immunity with over-and under-expression of some cytokines (8) could also participate in the risk of HCC.

An original Japanese study published in this issue of Liver International adds a path to the large body of evidence (see below) that DAA do not promote or enhance the HCC risk.

Authors nicely investigated the influence of HCV eradication by DAA therapy on HCC development [START_REF] Hidenori | The Impact of HCV Eradication by Directacting Antivirals on the Transition of Precancerous Hepatic Nodules to HCC: a Prospective Observational Study[END_REF]. They prospectively evaluated in 401 DAA-treated patients who did not have a history of HCC, the prevalence and the changes of non-hypervascular hypointense nodules (NHHNs) by gadolinium-ethoxybenzyl-diethylenetriamine pentaacetic acidenhanced magnetic resonance imaging (EOB-MRI) prior to the start of DAA therapy and during a periodic follow up after therapy. After DAA therapy patients underwent routine ultrasound surveillance for HCC every 3 to 6 months and EOB-MRI examination was repeated every 6 to 9 months when NHHNs were detected at baseline and every 6 to 12 months in addition to ultrasound in those without NHHNs at baseline. The progression of NHHNs detected at baseline to typical HCC, as indicated by hypervascularization, and the incidence of newly emergent NHHNs was analyzed. NHHNs were detected in 35 patients (8.7%) with a median size of 8.5 mm and SVR was achieved in 383 patients (95.5%). By comparing SVR patients after a median follow-up of 12.7 months post-SVR with propensity score-matched patients with persistent HCV infection, there was no difference in : 1. the incidence of hypervascularization of NHHNs to typical HCC among patients who had NHHNs at baseline (5 SVR patients and one non SVR patient); and 2. new emergence of NHHNs, 8 of the 350 patients without NHHs at baseline (2.3%) developped such nodules and no hypervascular typical HCC were detected after SVR. This well-done study, which strengths included the prospective analysis and the biopsy-proven diagnosis of HCC since all patients underwent hepatic resection, does not support the hypothesis that, even in high-risk patients with NHHNs during a 2-year observation period after SVR, the eradication of HCV achieved by IFN-free DAA therapy could enhance (or suppress!) HCC development.

In addition to these important results, mMethodological limitations have been pointed out (10) for some studies suggesting a higher risk of de novo HCC or HCC recurrence as well as a higher aggressiveness. Moreover, several retrospective or prospective studies [START_REF]ANRS collaborative study group on hepatocellular carcinoma (ANRS CO22 HEPATHER, CO12 CirVir and CO23 CUPILT cohorts). Lack of evidence of an effect of direct-acting antivirals on the recurrence of hepatocellular carcinoma: Data from three ANRS cohorts[END_REF][START_REF] Nahon | ANRS CO12 CirVir group. Incidence of Hepatocellular Carcinoma After Direct Antiviral Therapy for HCV in Patients With Cirrhosis Included in Surveillance Programs[END_REF][START_REF] Finkelmeyer | Risk of de novo Hepatocellular Carcinoma after HCV Treatment with Direct-Acting Antivirals[END_REF][START_REF] Calvaruso | Rete Sicilia Selezione Terapia-HCV (RESIST-HCV). Incidence of Hepatocellular Carcinoma in Patients With HCV-Associated Cirrhosis Treated With Direct-Acting Antiviral Agents[END_REF] did not confirm such a risk. A systematic review, meta-analysis and meta-regression of 26 studies (IFN-based regimen = 17, oral DAAs = 9) did not confirm a significant risk of both HCC occurrence and recurrence associated with DAA [START_REF] Waziry | Hepatocellular carcinoma risk following direct-acting antiviral HCV therapy. A systematic review, meta-analyses, and meta-regression[END_REF]. Similar results have been reported from the very large US veterans cohort retrospectively analyzed from 1999 to 2015: among the 62,354 HCV patients without detectable HCC, 35,871 (58%) were treated by INF-only (11,988 -33.4%-achieved SVR), 4, 535 (7.2%) were treated by DAAs + INF (2,763 -60.9%-achieved SVR) and 21,948 (35%) were given DAAs only (19,909 -90.7%-achieved SVR)(16); 3,271 incident HCC were diagnosed more than 6 months after initiation of antiviral treatment (mean follow-up: 6.1 yrs). Survival free of HCC by cirrhosis and SVR status was not different according to therapy with the following adjusted hazard ratio (aHR) for HCC: 0.29 (95% CI: 0.23-0.37) in DAA-only, 0.48 (95% CI: 0.32-0.73) in DAA+interferon and 0.32 (95% CI: 0.28-0.37) in interferon-only.

An Italian prospective study including 3 917 patients with a F3-F4 fibrosis score treated by DAA between January 2015 and June 2016 (SVR12 in 97.2 % of F3, in 92.7 % of F4 Child-Pugh A and 80% of F4 Child-Pugh B) also denied an increased HCC risk associated with DAAs : the HCC incidence of 1.18 % person/year (CI95% : 0.92-1.49) in the whole cohort of cirrhotic patients was not different than that of an historical control group [START_REF] Romano | Newly diagnosed hepatocellular carcinoma in patients with advanced hepatitis C treated with DAAs: A prospective population study[END_REF]. Incidence of HCC was related to the severity of the underlying liver disease (0.46% (95% CI 0.12-0.17) in F3, 1.49% (95% CI 1.03-2.08) in CTP-A and 3.61% (95% CI 1.86-6.31) in CTP-B during the first year and 0%, 0.2% and 0.69%, respectively in the second year.

By using weighting methodologies (IPTW, Cox models) given the hererogeneities of treated and untreated patients, at least in countries where treatment priorization was implemented, we established in the ANRS/AFEF Hepather cohort CO-22 first a decreased risk of de novo HCC and second a significant decrease in global, hepatic and extra-hepatic mortality in the former as compared to the latter [START_REF] Carrat | Clinical outcomes in patients with chronic hepatitis C following direct-acting antiviral therapy: a prospective cohort study[END_REF], results which are clearly reassuring: they underline the need of weighting methodologies to avoid misinterpretations.

Thus, if we consider that the DAA-associated risk of HCC is not an issue anymore, the issue of an impact of DAA on the tumor growth and aggressiveness previously suggested [START_REF] Reig | Unexpected early tumor recurrence in patients with HCV-related HCC undergoing interferon-free therapy[END_REF][START_REF] Conti | Early occurrence and recurrence of hepatocellular carcinoma in HCV-related cirrhosis treated with direct-acting antivirals[END_REF][START_REF] Romano | Newly diagnosed hepatocellular carcinoma in patients with advanced hepatitis C treated with DAAs: A prospective population study[END_REF] is still a matter of debate. The influence of DAA therapy on the timing (early or late), frequency and aggressiveness of HCC recurrence after a so-called "curative" treatment clearly needs further investigation.

Very few studies have evaluated the kinetics of tumor growth. We have performed a blinded re-evaluation of the tumor kinetics of the 40 patients of the Hepather cohort who had recurrence of HCC after a so-called "curative" treatment [START_REF]ANRS collaborative study group on hepatocellular carcinoma (ANRS CO22 HEPATHER, CO12 CirVir and CO23 CUPILT cohorts). Lack of evidence of an effect of direct-acting antivirals on the recurrence of hepatocellular carcinoma: Data from three ANRS cohorts[END_REF]: it is noteworthy that 7 patients (17.5%) were not cured and we did not observe any difference regarding the tumor kinetics between those who were given DAA and those who did not (unpublished data).

A residual unmet need is to define in HCV-infected patients undergoing a curative or a palliative treatment of HCC what is the better timing of the antiviral treatment.

1 5 1

 5 Editorial for LIVint-18-00714.R1 « The Impact of HCV Eradication by Direct-acting Antivirals on the Transition of Precancerous Hepatic Nodules to HCC: a Prospective Observational Study » Direct acting antivirals and the risk of hepatocellular carcinoma: the end of a polemic? Stanislas Pol 1,2,3 , Anaïs Vallet-Pichard 2 , Hélène Fontaine 2 , Fabrice Carrat 4,Université Paris Descartes ; 2 Assistance Publique -Hôpitaux de Paris, Hôpital Cochin, Liver Department, Paris; 3 Centre de Recherche Translationnelle and INSERM U-1223, Institut Pasteur, Paris; 4 Sorbonne Université, Institut National de la santé et de la Recherche Médicale, Institut Pierre Louis d'Epidémiologie et de Santé Publique, Paris ; 5 Assistance Publique -Hôpitaux de Paris, Hôpital Saint-Antoine, Unité de Santé Publique, Paris,


Marianne Souaiby

AND DIDIER HENRION Aneel Tanwani

ENSEMBLE APPROXIMATIONS FOR CONSTRAINED DYNAMICAL SYSTEMS USING LIOUVILLE EQUATION

For a class of state-constrained dynamical systems described by evolution variational inequalities, we study the time evolution of a probability measure which describes the distribution of the state over a set. In contrast to smooth ordinary differential equations, where the evolution of this probability measure is described by the Liouville equations, the flow map associated with the nonsmooth differential inclusion is not necessarily invertible and one cannot directly derive a continuity equation to describe the evolution of the distribution of states. Instead, we consider Lipschitz approximation of our original nonsmooth system and construct a sequence of measures obtained from Liouville equations corresponding to these approximations. This sequence of measures converges in weak-star topology to the measure describing the evolution of the distribution of states for the original nonsmooth system. This allows us to approximate numerically the evolution of moments (up to some finite order) for our original nonsmooth system, using a solver that uses finite order moment approximations of the Liouville equation. Our approach is illustrated with the help of two academic examples.

Introduction

In the theory of dynamical systems, studying the evolution of state trajectories, both qualitatively and quantitatively, is a common occurrence. For ordinary differential equations, with a fixed initial condition described by a point in the finite-dimensional vector space, the tools for analyzing the behavior of trajectories are widely available. However, for many applications, it is of interest to consider the evolution of dynamical systems when the initial condition is described by a distribution over some set in the state space. This article explores this latter direction for a particular class of nonsmooth dynamical systems.

If we consider a probability measure to describe the distribution of the initial conditions of a dynamical system, then the time evolution of this initial probability measure with respect to underlying dynamics is the object of our interest. For an autonomous dynamical system described by an ordinary differential equation (ODE) with Lipschitz continuous vector field, the time evolution of this measure is described by a linear partial differential equation (PDE) called the Liouville equation or the continuity equation, see e.g. [START_REF] Villani | Topics in optimal transportation[END_REF]Section 5.4]. The solution to the Liouville equation, that is the probability measure describing the distribution at time t, is the pushforward or image measure of the initial probability measure through the flow map at time t. Lipschitz continuity of the vector field ensures that the flow map of the ODE is invertible, which in turn ensures that the pushforward measure is the unique solution to the Liouville equation. This approach of associating the continuity equation with finite dimensional ODEs has found relevance in numerical optimal control [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF] as well as in several control-theoretic problems [START_REF] Bartsch | A theoretical investigation of Brockett's ensemble optimal control problems[END_REF][START_REF] Brockett | Optimal control of the Liouville equation[END_REF][START_REF] Brockett | Notes on the control of the liouville equation[END_REF].

When the vector field is not Lipschitz continuous, then the study of the evolution of the initial distribution is more involved. The first occurrence of continuity equations corresponding to nonsmooth ODEs occurs in [START_REF] Diperna | Ordinary differential equations, transport theory and sobolev spaces[END_REF]. Continuity equations corresponding to one-sided Lipschitz vector fields have been studied in [START_REF] Bouchut | One dimensional transport equation with discontinuous coefficients[END_REF][START_REF] Bouchut | Uniqueness and weak stability for multidimensional transport equations with one-sided Lipschitz coefficients[END_REF]. In [START_REF] Ambrosio | Transport equation and cauchy problem for non-smooth vector fields[END_REF], the authors consider less regular ODEs and study uniqueness of solutions for (Lebesgue) almost-all initial conditions by using the Liouville equation.

The dynamical systems for which we want to study the evolution of probability measures (describing the distribution of states) are the so-called constrained systems described by differential inclusions. In particular, given a closed convex set S ⊂ R n , and a continuously differentiable function f : R n → R n , we describe the evolution of constrained systems via the differential inclusion

(1) ẋ ∈ f (x) -N S (x)
where N S (x) ∈ R n denotes the outward normal cone to the set S at the point x ∈ R n . Since the normal cone takes a zero value in the interior of S, it is clear that the right-hand side of (1) is potentially discontinuous at the boundary of the set S. One can also think of (1) as an evolution variational inequality, described as ẋ(t) -f (x(t)), y -x(t) 0, for all y ∈ S, x(t) ∈ S, t ∈ [0, T ]. Such dynamical systems have been a matter of extensive study in past decades due to their relevance in engineering and physical systems. A recent survey article [START_REF] Brogliato | Dynamical systems coupled with monotone set-valued operators: Formalisms, applications, well-posedness, and stability[END_REF], and a research monograph [START_REF] Adly | A variational approach to nonsmooth dynamics[END_REF], provide an overview of different research oriented directions in the literature pertaining to system (1) and its connections to different classes of nonsmooth mathematical models. Analysis of such systems requires tools from variational analysis, nonsmooth analysis, set-valued analysis [START_REF] Aubin | Set-Valued Analysis[END_REF][START_REF] Mordukhovich | Variational Analysis and Generalized Differentiation[END_REF][START_REF] Rockafellar | Variational Analysis[END_REF]. Results based on stability analysis with computational aspects have been addressed recently by the authors in [START_REF] Souaiby | Cone-copositive Lyapunov functions for complementarity systems: Converse result and polynomial approximation[END_REF].

For a fixed initial condition, x(0) ∈ S, the question of existence and uniqueness of solution to system (1) has already been well-established in the literature, and the origins of such works can be found in [START_REF] Moreau | Preliminary version in: Problème d'évolution associé à un convexe mobile d'un espace Hilbertien[END_REF], see [START_REF] Edmond | BV solutions of nonconvex sweeping process differential inclusion with perturbation[END_REF] for a recent exposition. However, if we consider the initial conditions described by a probability measure, then the evolution of this measure under the dynamics of (1) has received very little attention in the literature. One can study such problems by considering stochastic versions of (1) by adding a diffusion term on the right-hand side. Such systems first came up in the study of variational inequalities arising in stochastic control [START_REF] Bensoussan | Applications des inégalités variationnelles en contrôle stochastique[END_REF], and in the literature, we can find results on existence and uniqueness of solutions in appropriate function space. In [C 95], this is done by considering Yosida approximations of the maximal monotone operator, whereas [START_REF] Bernardin | Multivalued stochastic differential equations: Convergence of a numerical scheme[END_REF] provides a proof based on time-discretization of system (1). These approaches have been generalized for prox-regular set S in [START_REF] Bernicot | Stochastic perturbation of sweeping process and a convergence result for an associated numerical scheme[END_REF], and the case where the drift term contains Young measures [START_REF] Castaing | Some problems in optimal control governed by the sweeping process[END_REF][START_REF] Castaing | A Skorokhod problem governed by a closed convex moving set[END_REF]. One could also, in principle, formulate a partial differential equation with set-valued elements and study the solutions of such equations under appropriate hypothesis, which is the case in [START_REF] Bonnet | Differential inclusions in Wasserstein spaces: The Cauchy-Lipschitz framework[END_REF] but it is not clear how to derive the corresponding set-valued partial differential equation for system (1) and whether the resulting inclusion would satisfy the necessary hypothesis for well-posedness. Different from these approaches, and inspired by the fact that the evolution of a probability measure for single-valued dynamical system is described by Liouville equation, it is natural to ask whether the evolution of a probability measure under the dynamics of system (1) can be studied using Liouville equation. To the best of authors' knowledge, this approach has only been adopted in [START_REF] Marino | Measure sweeping processes[END_REF], where the authors consider system of form (1) without the drift term f (•). Since the right-hand side of (1) is set-valued, it is not immediately clear how the divergence term in the Liouville equation is to be interpreted. In [START_REF] Marino | Measure sweeping processes[END_REF], the authors consider approximations to the solutions of Liouville equation associated with (1), which are similar to time-stepping algorithm. That is, a time-discretization technique is introduced which is based on projecting the density function on to the constraint set with respect to Wasserstein metric.

In this paper, we consider a different route for computing the approximate solution of system (1) in the space of probability measures. Inspired by the concepts presented in [START_REF] Ambrosio | Transport equation and cauchy problem for non-smooth vector fields[END_REF], our basic idea is to consider Lipschitz approximations of system (1). The particular approximations that we work with are the ones obtained by Yosida-Moreau regularization and are parameterized by a positive scalar converging to zero. We can then associate a single-valued Liouville equation to each of these approximants, and establish convergence of the resulting sequence of measures. Unlike [START_REF] Marino | Measure sweeping processes[END_REF], our approach for numerically solving the Liouville equation does not depend upon discretization in time, or space for that matter. Instead, we use functional discretization: we choose a family of test functions (the monomials) on which the evolution measure and the associated moments are then approximated numerically by a hierarchy of semidefinite programs. Furthermore, we also show that the support of the sequence of measures converges (with respect to the Hausdorff distance) to the support of the pushforward measure for the nonsmooth system. These analytical results allow us to get an approximation of the actual solution.

Since the pushforward measure, at each time instant, is an infinite-dimensional object, it can be challenging to approximate it numerically. One possibility -that we do not explore here -could to use Monte-Carlo probabilistic algorithms. Instead, we investigate a purely deterministic approach: in order to get a quantitative measure of the distribution of state at any time instant, which involves building a hierarchy of moments defined by the action of a finite Borel measure on polynomial test functions, and encoding the positivity constraints on moment matrix by using sum-of-squares (SOS) decomposition. This technique, called moment-SOS hierarchy [START_REF] Henrion | The Moment-SOS Hierarchy[END_REF] has been used in several engineering problems, and for our purposes, it allows us to approximate numerically the moments (up to some finite order) associated with the pushforward measure. Also, using the recent developments on approximating the support of a measure with the Christoffel-Darboux kernel [START_REF] Lasserre | The empirical Christoffel function with applications in data analysis[END_REF], we can approximate the support of the pushforward measure, and hence the trajectories corresponding to a certain initial distribution.

The remainder of the article is organized as follows: In Section 2, we formalize the problem and introduce the basic mathematical elements necessary for doing so. In Section 3, we construct Lipschitz approximations of our initial dynamical system. In Sections 4 and 5, we study certain properties of the sequence of measures associated with approximations constructed in Section 3. Numerical aspects for approximating the moments, and support, of the probability measure describing the evolution of system dynamics are also discussed in Sections 4 and 5. We illustrate our results with the help of two academic examples in Section 6. Some concluding remarks with possible future directions appear in Section 7, followed by an Appendix which collects some additional tools used in the development of our results.

ż(t) = g(t, z(t)), z(0) = z 0 ,
over a given time interval [0, T ], where g : [0, T ] × R n → R n is a given vector field and z(t) ∈ R n is the state. For each t ∈ [0, T ], let us consider the flow map G t : R n → R n , so that the mapping z 0 → G t (z 0 ) provides the value of state trajectory of (2) at time t, and moreover it satisfies

(3) ∂ t G t (z 0 ) = g(t, G t (z 0 )), G 0 (z 0 ) = z 0 , (t, z 0 ) ∈ [0, T ] × R n .
In this article, we consider the evolution of dynamical systems when the initial condition is defined probabilistically. In particular, we use the notation z(0) ∼ ξ 0 to mean that z(0) is a random variable whose law is a given probability measure, or density function ξ 0 ∈ P(R n ), where P(S) denotes the set of probability measures supported on S.

This model allows to capture an initial spatial distribution of particles. To define the corresponding density function at time t 0, denoted by ξ t ∈ P(R n ), we consider the pushforward or image measure of ξ 0 through the flow map G t (•). That is, let

ξ t := G t ξ 0 , (4)
so that, for every Borel subset B ⊂ R n , it holds that

ξ t (B) = ξ 0 (G -1 t (B)) = ξ 0 ({z ∈ R n : G t (z) ∈ B}).
The evolution of ξ t is described by the following PDE, called the continuity or Liouville equation:

(5)

∂ t ξ t + div(ξ t g) = 0,
with the initial condition:

(6) ξ| t=0 = ξ 0 .

The Liouville equation (5) should be understood in the sense of distributions, i.e. the family of probability measures t → ρ t is a measure-valued solution of (5)-(6) if :

• it is continuous in the sense that for every compactly supported continuous function φ : R n → R, the map

M φ : t → R n φ dξ t is continuous on [0, ∞) with M φ (0) = R n φ dξ 0 .
• for every r > 0 and every v ∈ C 1 ([0, T ] × R n ) such that v(T, •) = 0 and v(t, •) is supported on a closed ball of radius r for every t ∈ [0, T ], one has

T 0 R n (∂ t v(t, z) + ∂ z v(t, z) • g(t, z)) dξ t (z) dt = - R n v(0, z) dξ 0 (z).
The equivalence between the solutions of ODE (2) and PDE (5), is established in the following result, see e.g. [Vil03, Theorem 5.34]:

Theorem 2.1. For each t ∈ [0, T ], let G t : R n → R n be a diffeomorphism so that (3) holds. Given ξ 0 ∈ P(R n ), let ξ t be defined as in (4). Then, ξ t is the unique solution of the Liouville equation (5)-( 6) over the time interval [0, T ].

The importance of the Liouville PDE relies on its linearity in the probability measure ξ t , whereas the Cauchy ODE is nonlinear in the state trajectory z(t). This PDE governs the time evolution of a measure transported by the flow of a nonlinear dynamical system. The nonlinear dynamics is then replaced by a linear equation on measures. It is important to note that, in Theorem 2.1, the equivalence is established under the assumption that G t is a diffeomorphism for each t ∈ [0, T ], which in particular requires that the flow map G t is invertible. ODEs with Lipschitz vector fields have this property, but when the vector field is not Lipschitz continuous in state variable, the backward invertibility assumption may not hold, or the flow map G t may itself not be uniquely defined.

2.2. Ensembles of constrained system. In this paper, we are interested in studying a class of dynamical systems described by the variational inequalities

(7) ż(t) ∈ f (t, z(t)) -N S(t) (z(t)), z(0) ∼ ξ 0 ,
over an interval [0, T ] for some given T > 0, where f : [0, T ] × R n → R n is a given vector field, S : [0, T ] ⇒ R n a compact and convex-valued mapping, and we recall that the normal cone to S at z is defined by

(8) N S (z) := {λ ∈ R n | λ, z -z 0, ∀z ∈ S} .
If we consider a point z ∈ int(S), the interior of S, then N S (z) = 0 and by convention, we let N S (z) := ∅ for all z ∈ S. The formalism of system (7) with inclusion naturally allows us to describe dynamics constrained to evolve in set S. Using the depiction in Figure 2, it is seen that, during the evolution of a trajectory, if z(t) is in interior of S, then N S (z(t)) = 0 and the motion of the trajectory continues according to the differential equation ż(t) = f (t, z(t)). Whenever z(t) is on the boundary, we add a vector from the set -N S (z(t)), which restricts the motion of the state trajectory in tangential direction on the boundary of the constraint set S. The foregoing discussion motivates us to consider the following definition of a solution to (7) originating from a point mass in S(0): An absolutely continuous function z : [0, T ] → R n is called a solution to system (7) if there exists a selection

η(t) ∈ N S(t) (z(t)) such that ż(t) = f (t, z(t)) -η(t),
holds Lebesgue a.e. on [0, T ], and for each t ∈ [0, T ], we have z(t) ∈ S(t). The reader may refer to [START_REF] Brogliato | Dynamical systems coupled with monotone set-valued operators: Formalisms, applications, well-posedness, and stability[END_REF] for different formalisms and methods for describing the selection rule η.

For this paper, we emphasize that, in (7), ξ 0 ∈ P(S(0)) is a probability measure that specifies the distribution of the initial state. For each t ∈ [0, T ], let us denote the flow map by F t : S(t) → S(t), so that z 0 → F t (z 0 ) is the value at time t of the state trajectory of (7) with z(0) = z 0 . Given this random initial condition, the state at each time t can also be interpreted as a random variable in S(t), i.e. z(t)

∼ S = R 2 + N S (z) f (z) -NS(z) z(0) N S (0) f (z) Figure 2. State trajectories in constrained system with S = R 2 + .
ξ t ∈ P(S(t)) defined by ξ t := F t ξ 0 . However, unlike Lipschitz continuous ODEs, the mapping F t is not invertible in general. An example illustrating this fact is given next.

Example 1

(Flow map not invertible). Let f (z) = Az with A = 0 1 -1 0 and S = R 2
+ and let z 0 be a given initial condition, with angle θ 0 . For t θ 0 , we have

z(t) = F t (z 0 ) = e At z 0 = cos(t) sin(t)
-sin(t) cos(t) z 0 . And for t θ 0 , we have

z(t) = [|z 0 | 0] . For example if z 0 = [1 1] , it holds θ 0 = π 4 and then for t θ 0 , we have z(t) = √ 2 2 √ 2 2 - √ 2 2 √ 2 2 1 1 = [ √ 2 0] . The flow map reads z(t) = F t (z 0 ) = e At z 0 if t θ 0 [|z 0 | 0] if t θ 0 .
Indeed, as we can observe, the flow map is not invertible since given a state z(t) for a given time t θ 0 , it is not possible to retrieve the initial condition z 0 .

As a consequence of Example 1, it is seen that the flow map associated with dynamical system (7) is not necessarily invertible, and hence the conditions of Theorem 2.1 are not satisfied in general for such systems. On the other hand, for each t ∈ [0, T ], the forward flow map F t is well-defined and therefore the solution ξ t := F t ξ 0 exists and is uniquely defined. However, it is not possible to write down the evolution equation for ξ t , like Liouville equation for smooth ODEs, due to nonsmooth set-valued dynamics in (7). Recent literature in this direction deals with such problems, either by studying partial differential equations with set-valued mappings [START_REF] Bonnet | Differential inclusions in Wasserstein spaces: The Cauchy-Lipschitz framework[END_REF] or by introducing an approximation based on time discretization [START_REF] Marino | Measure sweeping processes[END_REF]. In this article, our goal is to find alternate methods based on functional discretization with monomial basis to approximate the measure ξ t and propose computational algorithms to calculate such approximations numerically. 2.3. Problem Formulation. We consider the dynamical system (7) with flow map F t : R n → R n . For a given ξ 0 ∈ P(S(0)), since there is no direct derivation of the PDE for characterizing the evolution of ξ t := F t ξ 0 , we compute an approximation of ξ t as follows:

• Construct a sequence of ODEs with Lipschitz continuous right-hand sides which approximate the solution of ODE (7) for a fixed initial condition. This construction is based on a regularization of (7), and results in a sequence parameterized by a scalar λ > 0.

• Exploit the regularity of the approximating ODE to construct a sequence of measures ξ λ t := F λ t ξ 0 . • When λ tends to 0, prove that ξ λ t converges to ξ t := F t ξ 0 in the weak-star topology. In particular, all finite order moments of ξ λ t converge to the moments of ξ t .

• When λ tends to 0, prove the convergence of the support of ξ λ t to the support of ξ t in the Hausdorff metric.

From a computational viewpoint, the by-product of the above results is that, for a fixed λ > 0, one can invoke efficient numerical methods for computing moments associated with the probability measure ξ λ t and the support of ξ λ t . This allows us to compute an approximation of the moments and support of ξ t associated with nonsmooth system (7).

Lipschitz Approximation

The first step in our analysis is to compute an approximation of the solutions of (7) by using Moreau-Yosida regularization. The development carried out here is inspired by [START_REF] Brogliato | Dynamical systems coupled with monotone set-valued operators: Formalisms, applications, well-posedness, and stability[END_REF]. We introduce a sequence of approximate solutions, the so-called Moreau-Yosida approximants {z λ } λ>0 , which are obtained by solving the following ODE parameterized by λ > 0:

(9) żλ (t) = f (t, z λ (t)) - 1 λ (z λ (t) -proj(z λ (t), S(t))), z λ (0) = z 0 ∈ S(0)
over the interval [0, T ], where proj(z, S) is the (unique) Euclidean projection of vector z onto convex set S. It is observed that, for each λ > 0, the right-hand side of ( 9) is (globally) Lipschitz continuous, and therefore, there exists a continuously differentiable trajectory z λ : [0, T ] → R n such that (9) holds for every t ∈ [0, T ].

The relation between the solution of the inclusion (7) and the approximants {z λ } λ>0 holds under the following assumptions:

Assumption 1. There exists a constant L f > 0 such that, for each t ∈ [0, T ], |f (t, z)| L f (1 + |z|), ∀z ∈ R n |f (t, z 1 ) -f (t, z 2 )| L f |z 1 -z 2 |, ∀z 1 , z 2 ∈ R n .
Assumption 2. The mapping S : [0, T ] ⇒ R n is closed and convex-valued for each t ∈ [0, T ], and S(•) varies in a Lipschitz continuous manner with time, that is, there exists a constant L S 0, such that

d H (S(t 1 ), S(t 2 )) L S |t 1 -t 2 |, ∀t 1 , t 2 ∈ [0, T ].
The notation d H (A, B) means the Hausdorff distance between sets A and B, that is,

(10) d H (A, B) := max sup y∈A d(y, B), sup x∈B d(x, A)
where d(x, A) denotes the Euclidean distance between vector x and set A.

Theorem 3.1. Under Assumptions 1-2, consider the sequence of solutions {z λ } λ>0 to parameterized ODE (9) on an interval [0, T ]. Then, as λ → 0, the sequence converges uniformly to a Lipschitz continuous function z : [0, T ] → R n , the unique solution to the differential inclusion (7).

The proof of this theorem is described in the remainder of this section. Note that the right-hand side of (7) is typically seen as a Lipschitz perturbation of a time-varying maximal monotone operator and the existence of solutions for such systems is proved by constructing a sequence of solutions and studying their convergence. Such sequences are either obtained by time-discretization or Moreau-Yosida regularization. The simplest case of static maximal monotone operators is studied in [START_REF] Brezis | Opérateurs Maximaux Monotones et Semi-Groupes des Contractions dans les Espaces de Hilbert[END_REF] but the situation is more complex when the domain of the multivalued mapping is time-dependent. Sweeping processes (defined by taking f ≡ 0 in (7) and deterministic initial condition) form a particular case of such systems and for the most part, the existence of solutions is proved using a discretization algorithm originating from [START_REF] Moreau | Preliminary version in: Problème d'évolution associé à un convexe mobile d'un espace Hilbertien[END_REF]; see [START_REF] Camlibel | Convergence of proximal solutions for evolution inclusions with time-dependent maximal monotone operators[END_REF] for some recent developments. Our contribution in Theorem 3.1 provides a proof based on Moreau-Yosida regularization for differential inclusions with normal cones associated to time-varying sets and Lipschitz perturbation, which does not appear in the literature. One is obviously helped by the machinery developed for the proofs with no perturbations but the intermediate calculations are different and we find it instructive to provide a proof for self-contained exposition. For the proof that follows, certain calculations, leading to the intermediate lemmas used in the proof, have been included in the appendix.

Proof of Theorem 3.1. The basic idea of the proof is to show that the sequence {z λ } λ>0 satisfies bounds ensuring uniform convergence to a function z(•) solving (7). This development is carried out in four steps.

Step 1: Estimates on the sequence {z λ } λ>0 . As a first step, to obtain bounds on the norm of z λ (.), let us begin by computing bounds on the norm of żλ (•) as stated in the following lemma, whose proof is given in A.

Lemma 3.2. For each λ > 0, it holds

(11) | żλ (t)| 2L f + L f |z λ (t)| + L f max 0 s t |z λ (s)| + L S ,
where L f , L S were introduced in Assumptions 1 and 2 respectively.

Based on Lemma 3.

d dt |z λ (t)| 2 2L f |z λ (t)| 2 + 2L f |z λ (t)| • max 0 s t |z λ (s)| + (4L f + 2L S )|z λ (t)|. Let y λ (t) = |z λ (t)| 2 , so d dt y λ (t) 2L f y λ (t) + 2L f y λ (t) • max 0 s t y λ (s) + (4L f + 2L S ) y λ (t).
Since the right-hand side of this differential inequality results in a nonnegative and nondecreasing function, it follows that y λ (t) y λ (t), for all t ∈ [0, T ], where

y λ satisfies d dt y λ (t) = 2L f y λ (t) + 2L f y λ (t). y λ (t) + (4L f + 2L S ) y λ (t) = 4L f y λ (t) + (4L f + 2L S ) y λ (t). ( 13
)
By using the substitution v(t) = ( y λ (t))

1 2 in (13), it yields v(t) = 2L f v(t) + 2L f + L S . The solution of this differential equation is v(t) = e 2L f t v(0) + (e 2L f t -1) (2L f +L S ) 2L f . Consequently, |z λ (t)| 2 = y λ (t) y λ (t) = v(t) 2 , and we obtain (14) |z λ (t)| e 2L f T |z λ (0)| + (e 2L f T -1) (2L f + L S ) 2L f , so that |z λ (t)| is bounded on the interval [0, T ], independently of λ.
Step 2: Extracting a converging subsequence. Based on the estimates in Step 1, there exists a subsequence of z λ (•) which converges to z(•). More formally, the following statement is obtained. The proof of Lemma 3.3 is a consequence of the Arzelà-Ascoli theorem since the sequence {z λi } i∈N is continuously differentiable and { żλi } i∈N is uniformly bounded. The limit function z(•) is also Lipschitz continuous in this case.

Step 3: Limit is a solution. To finish the proof of Theorem 3.1, we just need to show that the limit z(•) satisfies the differential inclusion (7). This particular step requires a variational inequality, which is stated in the following lemma, and its proof is given in A.

Lemma 3.4. There exists a continuous function ϕ : [0, T ] → R n that satisfies ϕ(t 1 ) ∈ S(t 1 ) and ϕ(s)

+ s t1 f (r, z λ (r)) dr ∈ S(s) for each s ∈ [t 1 , t 2 ], with t 1 , t 2 ∈ [0, T ]. Moreover, it holds that, (15) t2 t1 ϕ(s), ż(s) -f (s, z(s)) ds 1 2 z(t 2 ) - t2 t1 f (r, z(r)) dr 2 -z(t 1 ) 2 .
We now complete the proof of Theorem 3.1 by showing that the limit of the converging subsequence z(

•) satisfies ż(t) ∈ f (t, z(t)) -N S(t) (z(t)) that is, ξ - z(t), ż(t) -f (t, z(t))
0, for any ξ ∈ S(t) and for almost every t 0. This is indeed the case, since for every ξ ∈ S(t), we can take a Lipschitz continuous function ϕ : [t, T ] → R n with ϕ(t) = ξ such that, due to Lemma 3.4, we get

[t,t+ ] ϕ(s), ż(s) -f (s, z(s)) ds 1 2 z(t + ) - t+ t f (r, z(r)) dr 2 -z(t) 2 ,
and by letting ϕ(s) = ξ -(ξ -ϕ(s)), we obtain

[t,t+ ] ξ, ż(s) -f (s, z(s)) ds - [t,t+ ] ξ -ϕ(s), ż(s) -f (s, z(s)) ds 1 2 z(t + ) - t+ t f (r, z(r)) dr + z(t), z(t + ) - t+ t f (r, z(r)) dr -z(t) , which implies ξ, z(t + ) -z(t) - t+ t f (s, z(s)) ds - t+ t ξ -ϕ(s), ż(s) -f (s, z(s)) ds 1 2 z(t + ) - t+ t f (r, z(r)) dr + z(t), z(t + ) - t+ t f (r, z(r)) dr -z(t) .
From this, we get

ξ - 1 2 z(t + ) - t+ t f (r, z(r)) dr + z(t) , z(t + ) -z(t) - t+ t f (s, z(s)) ds t+ t ξ -ϕ(s), ż(s) -f (s, z(s)) ds -max s∈[t,t+ ] |ξ -ϕ(s)|| ż(s) -f (s, z(s))| -max s∈[t,t+ ] |ξ -ϕ(s)|| ż(s)| -L f max s∈[t,t+ ] |ξ -ϕ(s)|(1 + |z(s)|). Since z(•) is Lipschitz continuous, z(•) is bounded on [0, T ]
ξ - 1 2 z(t + ) - t+ t f (r, z(r)) dr + z(t) , z(t + ) -z(t) - t+ t f (s, z(s)) ds -M max s∈[t,t+ ] |ξ -ϕ(s)| -M L f max s∈[t,t+ ] |ξ -ϕ(s)|,
for some constant M > 0. Letting tend to zero, and recalling that ξ = ϕ(t) ∈ S(t),

we get ξ -z(t), ż(t) -f (t, z(t)) 0, for each ξ ∈ S(t),
and hence, z(•) satisfies the differential inclusion (7).

Remark 3.5. In the literature, we can find several proofs of convergence of solutions obtained from Moreau-Yosida regularization to the solution of systems closely related to (7), see for example [BT20, KM96, NT19]. The proof technique adopted here closely follows the outline given in [START_REF] Brogliato | Dynamical systems coupled with monotone set-valued operators: Formalisms, applications, well-posedness, and stability[END_REF], but the difference here is that we add the Lipschitz perturbation f (t, z) on the right-hand side of (7), which modifies certain calculations.

Convergence of Measures

Using the results from the previous section on the convergence of solutions for fixed initial condition, we now study the evolution of probability measures for system (7). As before, let us assume that z(0) is a random variable whose law is a given probability measure ξ 0 ∈ P(S(0)). We recall that the flow map for system (7) is denoted by F t , so that t → z(t) := F t (z 0 ) is the unique solution to (7).

For the Lipschitz approximation given in (9), consider the map F λ t : R n → R n , so that t → z λ (t) := F λ t (z 0 ) defines the unique solution to (9). Since the right-hand side of (9) is Lipschitz continuous for each λ > 0, we can consider a sequence of probability measures ξ λ t ∈ P(S(t)) defined as ξ λ t := F λ t ξ 0 for each t ∈ [0, T ] and λ > 0. From Theorem 2.1, it follows that ξ λ t satisfies the partial differential equation:

(16) ∂ t ξ λ t + div(ξ λ t f λ t ) = 0 in the sense of distributions, with the initial condition ξ t=0 = ξ 0 , and (17)

f λ t (z) := f (t, z) - 1 λ z -proj(z, S(t)) .
On the other hand, we do not know how to derive a meaningful PDE for ξ t . However, in the sequel, we show that the probability measure ξ t can be approximated by ξ λ t as λ → 0. This way, a good numerical approximation of ξ λ t would also provide an approximation of ξ t . 4.1. Weak-star convergence. We first show convergence in the weak-star topology. This allows us to approximate the evolution of the moments of the measure ξ t using the moments of ξ λ t . Given a measure ξ, we denote its support by supp(ξ), defined as the smallest closed set whose complement has zero measure with respect to ξ. Equivalently, it is the smallest closed set for which every point has a neighborhood of positive measure with respect to ξ.

Proposition 4.1. Let v : R n → R be a continuous function, and assume that ξ 0 has bounded support. Then,

(18) lim λ→0 R n v(z) dξ λ t (z) = R n v(z) dξ t (z).
Proof. By definition of the pushforward measure ξ λ t , it holds

(19) R n v(z) dξ λ t (z) = R n v(F λ t (y)) dξ 0 (y)
for all continuous functions v. From Theorem 3.1, for each t ∈ [0, T ], we have lim λ→0 z λ (t) = z(t), which is equivalent to

lim λ→0 F λ t (y) = F t (y), ∀y ∈ S(0).
Since v is any continuous function, this implies

lim λ→0 v(F λ t (y)) = v(F t (y)).
By assumption, v • F λ t is bounded on the bounded set supp(ξ 0 ). This allows us to invoke Lebesgue's dominated convergence theorem to get

(20) lim λ→0 R n v(F λ t (y)) dξ 0 (y) = R n v(F t (y))dξ 0 (y).
Hence, (19) and (20) yield

lim λ→0 R n v(z) dξ λ t (z) = R n v(F t (y))dξ 0 (y).
Using again the change of variables formula, we obtain

lim λ→0 R n v(z) dξ λ t (z) = R n v(z) dξ t (z)
for all continuous functions v on R n . Therefore, the equality in (18) is proved.

Remark 4.2. In the proof of Proposition 4.1, the boundedness of supp(ξ 0 ) was used to invoke dominated convergence theorem. The result of Proposition 4.1 extends in some cases where supp(ξ 0 ) is unbounded. In particular, if it can be shown that there exists a function g : [0, T ] × R n → R 0 such that, for each λ > 0,

F λ t (y) g(t, y), t ∈ [0, T ]
then the convergence in (18) holds for all continuous functions v which satisfy

R n v(g(t, y))dξ 0 (y) < ∞, t ∈ [0, T ].
4.2. Relations describing moments. An immediate consequence of Proposition 4.1 is that we can get a desired approximation of the moments of ξ t by choosing appropriate test functions v. This amounts to computing the moments of ξ λ t . We will now explore numerical techniques which allow us to compute the solution of (16) by computing the desired moments.

Toward this end, we first recall the notion of occupation measure associated with the trajectories of a nonlinear ODE (2). In the following, we denote the indicator function of a set B by I B that is, I B (z) = 1 when z ∈ B and I B (z) = 0 when z ∈ B.

Definition 1. Given an initial condition z 0 ∈ R n , the occupation measure of a trajectory G t (z 0 ) is defined by

µ(A × B|z 0 ) := A I B (G t (z 0 )) dt for every A, respectively B, contained in the Borel σ-algebra of [0, T ], respectively R n .
A geometric interpretation is that µ measures the time spent by the graph of the trajectory (t, G t (z 0 )) in a given subset A × B of [0, T ] × R n . An analytic interpretation is that integration with respect to µ is equivalent to time-integration along a system trajectory, that is,

[0,T ] v(t, G t (z 0 )) dt = [0,T ] R n v(t, z) µ(dt, dz|z 0 ) for every test function v ∈ C([0, T ] × R n ).
As a consequence of this last interpretation, we observe that the Liouville equation ( 16) can be equivalently written as a linear PDE satisfied by the occupation measures

dµ λ := dt dξ λ t , with µ λ 0 := δ 0 ξ 0 , µ λ T := δ T ξ T , which is (21) ∂ t µ λ + div(µ λ f λ ) + µ λ T = µ λ 0
which again should be understood in the sense of distributions, i.e.

R+ R n (∂ t v(t, z) + ∂ z v(t, z) • f λ (t, z)) dµ λ (t, z) = R+ R n v(t, z)(dµ λ T (t, z) -dµ λ 0 (t, z)), for all continuously differentiable functions v over [0, T ] × R d such that v(t, •) is supported on a compact set of R n for every t ∈ [0, T ].
We compute approximate moments of µ λ by applying the moment-SOS hierarchy. This method consists of minimizing a functional subject to the following constraints:

• The Liouville equation (21) expressed in the sense of distributions, as a linear constraints on the moments of µ λ and µ λ T .

• Necessary linear matrix inequality (LMI) constraints based on the dual of

Putinar's Positivstellensatz.

We will see in the following how to formulate the Liouville equation ( 21) as a linear moment constraint.

Let g be a polynomial vector field defined as

g : (z 1 , z 2 , . . . , z n ) z ∈ R n → (g 1 , g 2 , . . . , g n ) ∈ R n ,
and v be a monomial test function, with a maximum degree d ∈ N, defined as

v : (t, z) → t a z b := t a z b1 1 z b2 2 • • • z bn n , for all a ∈ N and b ∈ N n , with a + b 1 + b 2 + • • • + b n d. Besides, let us denote (22) y a-1,b := T 0 R n t a-1 z b dµ λ (t, z) and (23) y T a,b := T 0 R n t a z b dµ λ T (t, z), (24) y 0 a,b := T 0 R n t a z b dµ λ 0 (t, z).
We next provide an expression for the solution of ( 21). In what follows, we let e i ∈ R n denote the vector whose only non-zero entry is equal to one at position i.

Proposition 4.3. The Liouville equation ( 21) is equivalently expressed as:

(25) y T a,b -y 0 a,b = ay a-1,b + n i=1 T 0 R n b i t a z b-ei g i (z) dµ λ (t, z)
which are linear constraints that link the moments of the initial measure, terminal measure and occupation measure.

Proof. Choosing v(t, z) = t a z b as a monomial test function, the Liouville equation ( 21) is then written as

∂ t µ λ , v + div(µ λ g), v + µ λ T , v = µ λ 0 , v , which implies (26) T 0 R n (∂ t v(t, z) + ∂ z v(t, z) • g(z)) dµ λ (t, z) = T 0 R n v(t, z) (dµ λ T (t, z) -dµ λ 0 (t, z)).
We have ∂ t v(t, z) = at a-1 z b , and

∂ z v(t, z) = (b 1 t a z b1-1 1 z b2 2 • • • z bn n , b 2 t a z b1 1 z b2-1 2 • • • z bn n , . . . , b n t a z b1 1 z b2 2 • • • z bn-1 n ).
Replacing ∂ t v(t, z) and ∂ z v(t, z) by their expressions in (26) yields

T 0 R n (at a-1 z b + n i=1 b i t a z b-ei g i (z)) dµ λ (t, z) = T 0 R n t a z b dµ λ T (t, z) - T 0 R n t a z b dµ λ 0 (t, z)
which is the expected statement by using the notations ( 22), ( 23) and (24).

Numerical computation.

Based on the result of Proposition 4.3, we now describe a numerical method for computing y T a,b . It is assumed that the initial measure dµ 0 is given, which allows us to compute y 0 a,b . We next describe the main steps involved in writing a semidefinite program for calculating y T a,b corresponding to the measure dµ λ . Note that, for each λ > 0, the measure µ λ is supported on a subset of R n+1 . In what follows, we provide some elements of construction for our algorithm for a finite Borel measure µ supported on R p . Given a Borel probability measure µ and α ∈ N p , we let

y α (µ) = R p z α dµ(z),
where we recall that z α := z α1 z α2 . . . z αp p . We consider the set {α ∈ N

p ; α 1 + • • • + α p
d} with graded lexicographic order, and denote it by N p d ; for example, with p = 2, d = 2, N 2 2 = {(0, 0), (1, 0), (0, 1), (2, 0), (1, 1), (0, 2)}. The cardinality of N p d is s(d) := p+d d , which is the number of monomials of degree at most d. The sequence y = (y α (µ)) α∈N p therefore encodes the moments of the measure µ.

The moment matrix of degree d associated with a Borel measure µ, denoted by M d (µ) is a matrix of dimension s(d) × s(d), whose rows and columns are indexed by monomials of degree at most d. For α, β ∈ N p d , the corresponding entry in M d (µ) is defined by (M d (µ)) α,β := y α+β (µ). As an example, once again with p = 2, d = 2, M 2 (µ) ∈ R 6×6 , and the element in second row (α = (1, 0)), third column (β = (0, 1)), corresponds to

R 2 z 1 z 2 dµ(z). To see an alternate representation of M d (µ), let b d (z) := (z α ) α∈N p d ∈ R[z] s(d) d
denote the vector of monomials of degree less than or equal to d, with graded lexicographic order. If the sequence {y α } α∈N p has a representing measure µ, i.e. y α = R p z α dµ(z) for all α ∈ N p , we can use the equivalent definition

M d (µ) := R p b d (z)b d (z) dµ(z)
, where the integral is understood entrywise. We can also define the localizing matrix of degree d with respect to a given q(z) ∈ R[z] by

M d-deg(q)/2 (qµ) := R p q(z)b d (z)b d (z) dµ(z)
where x denotes the smallest integer greater than x.

Assume that X ⊂ R n is a compact basic semialgebraic set i.e. The moment-SOS hierarchy, based on Proposition 4.4, allows us to compute approximate moments of the occupation measure and terminal measures. Recall that the moments of the initial measure are given since the initial measure is given. We fix a degree d ∈ N and we consider the linear constraint (25) linking moments of degree up to d, and subject to the constraints that the localizing matrices of the occupation measure and terminal measure, truncated to moments of degree up to d, are all positive semi-definite. This results in a finite-dimensional feasibility problem describe by linear matrix inequalities. The higher is the relaxation degree d, the better are the approximate moments, in the sense that when d tends to infinity, Proposition 4.4 and linear constraint (25) ensure that we have indeed moments of measures satisfying the Liouville equation. For a finite relaxation degree d, it may however happen that the approximate moments are not genuine moments of occupation measures. This is why the approximate moments are sometimes referred to as pseudo-moments.

X := {z ∈ R n : p k (z) 0, k = 0, . . . , n X } for given p k ∈ R[z], k = 0, . . . , n X . Let p 0 (z) =
The LMI constraints are automatically constructed by the msdp command in Gloptipoly for Matlab [START_REF] Henrion | Gloptipoly 3: moments, optimization and semidefinite programming[END_REF]. For more details about the LMI constraints, we refer the reader to [Hen13, Section 3.3] or the two introductory chapters of [START_REF] Henrion | The Moment-SOS Hierarchy[END_REF].

Convergence of Support of Measures

For several applications, it is important to approximate the support of the measure ξ t , since it provides a probabilistic estimate of the state trajectories at time t ∈ [0, T ]. Once again, our goal is to approximate the support of ξ t by the support of ξ λ t where ξ λ t satisfies (16).

5.1. Hausdorff convergence of support. We first show that supp(ξ λ t ) converges in the Hausdorff distance to supp(ξ t ).

Proposition 5.1. For each t ∈ [0, T ], it holds

(27) lim λ→0 d H (supp(ξ λ t ), supp(ξ t )) = 0.
Proof. First, let A λ t := supp(ξ λ t ) and A t := supp(ξ t ). For proving that lim

λ→0 d H (A λ t , A t ) = 0,
we need to prove the following two limits:

(28) lim For proving (28), we first observe that sup

y λ ∈A λ t d(y λ , A t ) = sup y λ ∈A λ t inf x∈At |y λ -x|,
and hence it needs to be shown that for every y λ ∈ A λ t , there exists x ∈ A t such that |x -y λ | converges to zero as λ converges to zero. Since y λ ∈ A λ t , there exists z 0 ∈ supp(ξ 0 ) such that y λ = F λ t (z 0 ). By choosing x = F t (z 0 ) ∈ A t , it follows from Theorem 3.1 that lim λ→0 F λ t (z 0 ) = F t (z 0 ), or equivalently, |x -y λ | converges to 0 as λ → 0.

For proving (29), we similarly observe that

sup x∈At d(x, A λ t ) = sup x∈At inf y λ ∈A λ t |x -y λ |.
Following the same idea as before, let us take x ∈ A t , then there exists z 0 ∈ supp(ξ 0 ) such that x = F t (z 0 ). By choosing y λ = F λ t (z 0 ) ∈ A λ t , it again follows from Theorem 3.1 that |x -y λ | converges to 0 as λ → 0, and (29) is obtained. 5.2. Approximation of support. Just like the approximation of moments, we can provide some numerical methods to approximate the support of the sequence of measures ξ λ t . By Proposition 5.1, by computing such an approximation for λ > 0 sufficiently small, we get an approximation of the support of the probability measure ξ t for the original system.

The technique we present is based on approximating the support of a measure by the sublevel sets of a polynomial function. In particular, for a finite Borel measure µ, we consider the moment matrix M d (µ) and introduce the mapping,

R n x → Λ ξ,d (x) := b d (x) M d (µ) -1 b d (x) ∈ R,
which we call Christoffel-Darboux polynomial. Thus, the basic idea behind the construction of the support of the measure µ is to use the finite order moments, and show that the sublevel sets of the Christoffel-Darboux polynomial indeed converge to the actual support of µ. This technique has been proposed in [START_REF] Lasserre | The empirical Christoffel function with applications in data analysis[END_REF] for stationary measures under certain hypothesis. Here, we show that by placing certain hypothesis on the initial measure ξ 0 , the approximations ξ λ t obtained by the Liouville equation satisfy the required hypothesis, which allow us to approximate the support of ξ λ t by constructing the corresponding Christoffel-Darboux polynomial. The following statement shows the existence of a sublevel set that approximates the support of the sequence of measures ξ λ t , when λ and t ∈ [0, T ] are fixed. Proposition 5.2. Let ξ 0 be absolutely continuous with respect to the Lebesgue measure and let us suppose that supp(ξ 0 ) is compact. For a fixed λ > 0, and t ∈ [0, T ], consider ξ λ t obtained by solving (16), and M λ d,t (ξ λ t ) the corresponding moment matrix of degree d. For every > 0 (small enough), there exists d ∈ N (large enough) and γ d > 0, such that the sublevel set

(30) S λ d,t := {z ∈ R p | b d (z) M λ d,t (ξ λ t ) -1 b d (z) γ d } satisfies (31) d H (S λ d,t , supp(ξ λ t )) , as d → +∞.
Proof. For each λ > 0 and t ∈ [0, T ], if we show that

• The set supp(ξ λ t ) is compact and has nonempty interior. • It holds that ξ λ t is absolutely continuous with respect to the Lebesgue measure.

then, the statement follows by applying [LP19, Theorem 3.11] to the measure ξ λ t . The aforementioned properties basically follow from the fact that, for a fixed t ∈ [0, T ] and λ > 0, the mapping F λ t : R n → R n is a homeomorphism obtained from the solution of an ODE with Lipschitz continuous right-hand side (9). Let L λ denote the (uniform with respect to time) Lipschitz constant for the mapping on the right-hand side of (9). One can readily show that for a pair of initial conditions y 0 , z 0 and y t := F λ t (y 0 ),

z t := F λ t (z 0 ), it holds that |z 0 -y 0 | exp (-L λ t) |z t -y t | |z 0 -y 0 | exp (L λ t).
Using this estimate, and recalling that ξ λ 0 = ξ 0 , it readily follows that supp(ξ λ t ) is compact and has nonempty interior under the given hypothesis on ξ 0 . Absolute continuity of ξ λ t with respect to Lebesgue measure holds if ξ λ t is absolutely continuous with respect to ξ 0 . The later indeed holds because for every

F λ t implies that (32) ξ 0 (A) = 0 ⇒ ξ λ t (A) = ξ 0 ((F λ t ) -1 (A)) = 0,
whence the desired result follows.

Illustrative examples

In this section, we give two examples that illustrate the computation of the pseudo-moments associated with ξ λ t of the regularized system (9) in the case where f : R 2 → R 2 , by applying the methods discussed in this article.

Example 2. Consider the constrained system (7) of Example 1 where f (z) = Az with A = 0 1 -1 0 and S = R 2 + . Let us write the regularized system (9) in polar coordinates (r, θ) as follows:

ṙ(t) = 0, θλ (t) = -1 -1 λ (θ λ (t) -proj(θ λ (t), S(t))).
or equivalently:

(33) ṙ(t) = 0, θλ (t) = -1 -1 λ (θ λ (t) -max(θ λ (t), 0)).
Let d = 4 be the degree of relaxation, and let us choose different values of the regularization parameter λ ∈ {0.05, 0.1, 0.5}. We introduce the initial measure as a Dirac measure with respect to time product a uniform measure in [0, 1] × [0, 1 2 ] with respect to the state.

We calculate the moments of the initial measure to replace it directly in Liouville constraint (25), where the variables z 1 and z 2 in (25) are respectively r and θ. For all (a, b 1 , b 2 ) ∈ N 3 , with a + b 1 + b 2 d, the moment of the initial measure is then given as

T 0 R n t a z b1 1 z b2 2 dµ 0 (t, z) = T 0 R n t a z b1 1 z b2 2 δ 0 (dt)λ [0,1] (dz 1 )λ [0, 1 2 ] (dz 2 ) = 0 a 1 0 z b1 1 dz 1 1 2 0 2z b2 2 dz 2 = 0 a 1 b 1 + 1 1 b1+1 -0 b1+1 2 b 2 + 1 1 2 b2+1 -0 b2+1 .
Then we apply the moment-SOS hierarchy [START_REF] Henrion | The Moment-SOS Hierarchy[END_REF] which allows us to approximate numerically the moments of the unknown occupation measure and terminal measure. For different values of the terminal time T ∈ {0, 0.1, 0.2, . . . , 1}, this gives us:

• The evolution of the moment r(t) 2 dµ λ T as a function of time, which we observe numerically is a constant for different values of the regularization parameter λ.

• The evolution of the moment θ(t) 2 dµ λ T as a function of time for different values of the regularization parameter λ, which is illustrated on Figure 3.

Example 3. We now consider an example with moving set S, which is described as follows:

(34) S(t) := {x ∈ R 2 | g i (t, x) 0, i = 1, 2, 3, 4}
where, for each i ∈ {1, 2, 3, 4}, we take g i to be an affine function with time-varying coefficients:

g i (t, x) = a i (t)(x -c(t)) -b i with a i (t), c(t) ∈ R 2 , b i ∈ R 2 >0 for each t ∈ [0, T ]
. Furthermore, the coefficients a i (•) satisfy the following constraints a 2 (t) = -a 1 (t), a 3 (t)a 1 (t) = 0, a 4 (t) = -a 3 (t),

a i (t) a i (t) = 1, i ∈ {1, 2, 3, 4},
for each t ∈ [0, T ], so that the resulting constraint is a rectangle moving in plane with center at c(t). An illustration of the constraint set, for a fixed t, appears in Figure 4.

For the sake of simulations, we take

a 1 (t) = 1 -t 2 1 + t 2 , 2t 1 + t 2 , a 3 (t) = -2t 1 + t 2 , 1 -t 2 1 + t 2 b 1 = b 2 = b 3 = b 4 = 1 4 c(t) = 2t 2 -1, 4t 3 -3t .
The dynamical system under consideration in this example is

(35) ẋ(t) ∈ -N S(t) (x(t))
and the Moreau-Yosida regularization of this system is, For each t ∈ [0, T ], S(t) has an affine representation. We let S j (t) := {x ∈ R 2 | g j (t, x) 0}, and the index of active constraints is defined as I act = {j ∈ {1, 2, 3, 4} | g j (t, x) > 0}. Using the fact that the vectors a i (t) have norm equal to 1 for each t 0 and are mutually orthogonal, we can write the projection vector as follows:

(36) żλ = f λ (t, z λ ) := - 1 λ (z λ (t) -proj(z λ (t), S(t))). c(t) a 1 (t) a 3 (t) a 4 (t) a 2 (t)
proj(z, S(t)) = z - j∈Iact (a j (t) (z -c(t)) -b j (t))a j (t).
Using the approach proposed in this paper, we consider the evolution of the uniformly distributed initial condition over the set S(0) via the Liouville equation associated with (36). For numerical computation of the pseudo-moments associated to the occupation measure, we use GloptiPoly. The implementation requires us to provide data in the form of polynomials, and since the vector field in (36) is in the rational form, we have to define a new (scaled) occupation measure as explained next. For the approach discussed in Section 4.2, we solve the following equation

[0,T ] R 2 ∂ t v(t, z) + ∂ z v(t, z) • h λ (t, z) (1 + t 2 ) 2 dµ λ (t, z) = [0,T ] R 2 v(t, z)(dµ λ T (t, z) -dµ λ 0 (t, z)),
for a given test function v(t, z), with the measure µ λ (t, z) as the unknown, and h λ (t, z) = (1 + t 2 ) 2 f λ (t, z) being a polynomial function. To rewrite the foregoing term purely in terms of polynomial functions, instead of rationals, we now consider a new measure of the form

dν λ (t, x) = dµ λ (t, x) (1 + t 2 ) 2 = dt dξ λ (x|t) (1 + t 2 ) 2 ,
so that the modified transport equation takes the following form:

[0,T ] R 2 ∂ t v(t, z)(1 + t 2 ) 2 + ∂ z v(t, z) • h λ (t, z) dν λ (t, z) = [0,T ] R 2 v(t, z)(dµ λ T (t, z) -dµ λ 0 (t, z))
with ν λ , µ λ T and µ λ 0 as the unknowns. By choosing v(t, z) = t a z b , all the entries are in polynomial form. It is therefore possible to compute the moments associated to the measure dµ λ of the following form:

(37) [0,T ] R 2 t a x b dµ λ T (t, x) = [0,T ] R 2 (t a x b )(1 + t 2 ) 2 dν λ T (t, x),
where the term on the right-hand sides are computed by our solver.

For illustration of our approach, our objective is to compute the empirical moments associated with the particles which evolve according to the equation (36), and then compare them with the moments which are obtained by the distribution satisfying the Liouville equation with polynomial coefficients associated to system (36). The simulation results are plotted in Figure 6 where we plot the value of (37) by taking a = 0, b = (1 0) and b = (0 1) for T ∈ [0, 1] and different values of λ. Along with the ensemble of the trajectories of the differential inclusion (35) with different initial conditions, we plot the empirical moments using the time and space discretization of the regularized ordinary differential equation (36), and the pseudo-moments of the terminal measure computed with GloptiPoly and MOSEK for certain relaxation order. Solving a moment LMI relaxation of order 4 takes a few seconds on a standard desktop PC. We observe experimentally that increasing the relaxation order does not really improve the accuracy of the results: the optimization problem becomes larger and worse conditioned. Similarly, decreasing the regularization parameter generate large coefficients in the vector field and hence in the coefficients of the LMI problem, and as a result, the optimization problem becomes worse conditioned. On the other hand, the accuracy of the empirical moments computed by time and space discretization of the ODE (36) improves as we decrease the regularization parameter λ.

Conclusions

In this article, we studied the time evolution of nonsmooth constrained dynamical systems when the initial condition is described by a probability measure. Unlike conventional ODEs, it is not obvious how to describe the time evolution of the image measure by the flow as a Liouville PDE. To circumvent this issue, we propose an approximation technique based on constructing Lipschitz approximations for the original nonsmooth system, and then using the Liouville equation for the approximate Lipschitz dynamics. Numerical methods for computing the approximation of solutions of Liouville equation then allow us to compute the moments and support of the probability measures associated to the original system. While evolution of probability measure for a class of constrained systems has been studied in [START_REF] Marino | Measure sweeping processes[END_REF], here we adopted a different approach for computing the approximations of differential equations associated with the evolution of probability measure. To seek generality in the class of systems studied in our work, one could, just like [START_REF] Marino | Measure sweeping processes[END_REF], consider a congestion constraint formulated as a uniform bound on the density of the measure. With the method proposed in this article, we can readily deal with such constraints, which would translate into a linear semidefinite constraint on the moments.

To seek improvements in the approach adopted in this paper, it is observed that the proposed Lipschitz approximations are difficult to simulate numerically. In particular, for the illustrated examples, we implemented the projection map onto a set by splitting the Liouville equation in different parts, where each of them corresponds to the region where the approximating ODE is continuous. One could use some recent work on approximating ODEs with twice differentiable right-hand side [START_REF] Colombo | Existence and uniqueness of solutions for an integral perturbation of moreau's sweeping process[END_REF] to see if the resulting implementation is easier to simulate for a broader class of constraint sets.

Another potential direction of research that comes out from this work is the possibility of using the proposed tools for optimal control problems. As was done for ODEs [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF], it is possible to use the formalism of Liouville equation for optimal control problems. The optimal control for the class of nonsmooth systems studied in this paper is a challenging problem, and it has been addressed recently in [CCMN21, CHHM16, VBP20]. It would be interesting to see if the methods proposed in this paper provide a numerically constructive solution to such challenging problems.

Appendix A. Proofs of Lemmas in Theorem 3.1 A.1. Proof of Lemma 3.2. For each λ > 0, the dynamics for z λ in (9) yield

| żλ (t)| = |f (t, z λ (t)) - 1 λ (z λ (t) -proj(z λ (t), S(t)))| |f (t, z λ (t))| + 1 λ |z λ (t) -proj(z λ (t), S(t))|. (38)
For the first term on the right-hand side of (38), we have that

(39) |f (t, z λ (t))| L f (1 + |z λ (t)|).
For the second term on the right-hand side of (38), we introduce the function We take limits with respect to λ → 0. Since z λ (•) converges pointwise to z(•), we have ϕ λ (s), żλ (s) -f (s, z λ (s)) -→ ϕ(s), ż(s) -f (s, z(s)) for each s ∈ [t 1 , t 2 ], and z λ (t 2 )-t2 t1 f (r, z λ (r)) dr 2 -→ z(t 2 )-t2 t1 f (r, z(r)) dr 2 , and z λ (t 1 ) 2 -→ z(t 1 ) 2 . Therefore, this yields to t2 t1 ϕ(s), ż(s) -f (s, z(s)) ds 1 2 z(t 2 ) -t2 t1 f (r, z(r)) dr 2 -z(t 1 ) 2 , and Lemma 3.4 is then proved.

d λ (t) = inf y∈S(t) |y -z λ (t)|, so that d λ (t) = d S(t) (z λ (t)). It is seen that d λ (t) = |z λ (t) -proj(z λ (t), S(t))|. So 1 λ |z λ (t) -proj(z λ (t), S(t))| = 1 λ d λ (t

  Lemma 3.3. There exists a subsequence {z λi } i∈N which converges uniformly to a Lipschitz continuous function z(•) on [0, T ].

  1 and let one of the inequalities p k (z) 0 be of the form R -n i=1 z 2 i 0 where R is a sufficiently large positive constant. Proposition 4.4. (Putinar's Theorem) The sequence of moments y has a representing measure supported on X if and only if M d-degp k /2 (p k µ), k = 0, . . . , n X are positive semidefinite for all d ∈ N.

Figure 3 .

 3 Figure 3. First order moment of the second state (vertical axis) of the occupation measure of the regularized system, as a function of time (horizontal axis), for different values of the regularization parameter (top curve λ = 0.5, middle curve λ = 0.1, bottom curve λ = 0.05)

Figure 4 .

 4 Figure 4. Box constraint formed by the intersection of four halfspaces where the coefficients a i and the center c change with time.

Figure 5 .

 5 Figure5. From left to right -thin black lines: box S(t) for t ∈ {0, 0.6, 0.8, 1}; light gray: ensemble of points initially uniformly supported on S(0) and transported by the flow; dark gray: 9 distinguished points transported by the flow. We observe that the points quickly concentrate on a corner of the box.

Figure 6 .

 6 Figure6. For two values of regularization parameter λ, we represent in dark gray the empirical first moments of the trajectory of (36) (obtained by sampling initial conditions) and in black the first pseudo-moments of the terminal measure computed with Glopti-Poly and MOSEK for the relaxation of order 4 (moment-SOS of degree 8). In light gray we represent for reference the ensemble trajectories for a uniform distribution of the initial condition in the box S(0).

ϕfff

  (s) := lim λ→0 ϕ λ (s), we observe that ϕ(•) is continuous and ϕ(s) ∈ S(s) := S(s) -s t1 f (r, z(r)) dr. Let zλ (s) := proj(z λ (s), S(s)); then s → zλ (s) is a continuous mapping. Since ϕ λ (s) + s t1 f (r, z λ (r)) dr ∈ S(s) and λ is positive, it follows from the definition of the projections thatϕ λ (s) + s t1 f (r, z λ (r)) dr -zλ (s), żλ (s) -f (s, z λ (s)) (r, z λ (r)) dr -zλ (s), zλ (s) -z λ (s) 0.Thenϕ λ (s), żλ (s) -f (s, z λ (s)) zλ (s) -s t1 f (r, z λ (r)) dr, żλ (s) -f (s, z λ (s)) , which implies that, t2 t1 ϕ λ (s), żλ (s) -f (s, z λ (s)) ds (r, z λ (r)) dr, żλ (s) -f (s, z λ (s)) ds.Since at the points where z λ (.) is differentiable, we have zλ (s), żλ (s) -f (s, z λ (s)) = zλ (s) -z λ (s), żλ (s) -f (s, z λ (s))+ z λ (s), żλ (s) -f (s, z λ (s)) = 1 λ |z λ (s) -z λ (s)| 2 0 + z λ (s), żλ (s) -f (s, z λ (s)) , it follows that, zλ (s), żλ (s) -f (s, z λ (s)) z λ (s), żλ (s) -f (s, z λ (s)) ,and,t2 t1 ϕ λ (s), żλ (s) -f (s, z λ (s)) ds t2 t1 z λ (s) -s t1 f (r, z λ (r)) dr, żλ (s) -f (s, z λ (s)) ds. (r, z λ (r)) dr, żλ (s) -f (s, z λ (s)) ds λ (t 2 ) -t2 t1 f (r, z λ (r)) dr 2 -z λ (t 1 ) 2 ,hence, we obtain that t2 t1 ϕ λ (s), żλ (s) -f (s, z λ (s)) ds1 2 z λ (t 2 ) -t2 t1f (r, z λ (r)) dr 2 -z λ (t 1 ) 2 .

  Evolution of probability measure ξ t w.r.t. time and space.

	evolution of ξ	z
	t 0	t 1	t 2	t
	2. Preliminaries and Problem formulation
	2.1. Evolution of ensembles. Let us consider the time-varying ODE
	(2)

t Figure 1.

  and differentiable almost everywhere on [0, T ]. Hence, for almost every t ∈ [0, T ], where z(•) is differentiable, dividing the last inequality by , we get

  ). To obtain a bound on d λ (t), we compute the derivative of d 2For the first term in the limit, we use thatd 2 S(t+ ) (z λ (t + )) -d 2 S(t) (z λ (t + )) d H (S(t + ), S(t)) d S(t+ ) (z λ (t + )) + d S(t) (z λ (t + )) | |L S d S(t+ ) (z λ (t + )) + d S(t) (z λ (t + )) . (z λ (t) + żλ (t)) -d 2 S(t) (z λ (t)) + d S(t) (z λ (t + )) -d S(t) (z λ (t) + żλ (t)) d S(t) (z λ (t + )) + d S(t) (z λ (t) + żλ (t)) . Since z λ (.) is differentiable, z λ (t + ) = z λ (t) + żλ (t) + O( ) and hence d S(t) (z λ (t + )) -d S(t) (z λ (t) + żλ (t)) = O( ). This implies that (t) f (t, z λ (t)) + 2L S d λ (t). -(t-s)/λ f (s, z λ (s)) + L S ds.Or, d λ (0) = |z 0 -proj(z 0 , S(0))| = 0 since z 0 ∈ S(0) and we have that f satisfies (39), then it follows-(t-s)/λ (L f + L f |z λ (s)| + L S ) ds.And therefore, substituting (39) and (43) in (38), we get| żλ (t)| L f + L f |z λ (t)| + 1 λThe bound of | żλ (t)| is then expressed as| żλ (t)| 2L f + L f |z λ (t)| + L f max 0 s t |z λ (s)| + L S .A.2. Proof of Lemma 3.4. For fixed t 1 , t 2 ∈ [0, T ], letS λ (s) := S(s)with s ∈ [t 1 , t 2 ]. Under the assumptions imposed on S(•) and the mapping f (•, •), the set-valued mapping S λ (•) is Hausdorff continuous. Thus, by a theorem on continuous selections[START_REF] Michael | Continuous selections. I[END_REF], there exists a continuous function ϕ

	By substitution of (41) and (42) in equation (40), we obtain
	2d λ (t) ḋλ (t) = λ (t) + 2d λ Dividing by 2d λ (t), we get d dt d 2 λ (t) -2 λ d 2
						d dt	d λ (t) -	1 λ	d λ (t) + f (t, z λ (t)) + L S ,
	which implies that,
										t
	λ (t): S(t) (z λ (t)) S(t+ ) (z λ (t + )) -d 2 S(t) (z λ (t)) d 2 = lim λ (t) = d 2 d dt →0 d 2 = lim d λ (t) e -t/λ d λ (0) + d dt 0 e (43) (40) 1 λ d λ (t) t 1 λ 0 e We have
	L f λ	0	t	e -(t-s)/λ ds =	L f λ	e -t/λ λe s/λ t 0	=	L f λ	e -t/λ λe t/λ -λ
							= L f 1 -e -t/λ	L f .
	(41) Similarly,
	L S For the second term in the limit, we first notice that λ
	Besides, we have	d 2 S(t) (z λ (t + )) -d 2 S(t) (z λ (t))
	= d 2 S(t) (z λ (t + )) -d 2 S(t) (z λ (t)) = d 2 S(t) (z λ (t) + żλ (t)) L f λ t 0 e -(t-s)/λ |z λ (s)| ds t L f λ 0 e -(t-s)/λ ds L f S(t) d 2 L f max 0 s t |z λ (s)|.	. max 0 s t	|z λ (s)|
										-d 2 S(t) (z λ (t)).
	And,									s
				lim →0	1 d 2 S(t) (z λ (t + )) -d 2 S(t) (z λ (t)) t1 f (r, z λ (r)) dr,
				= ∇d 2 S(t) (z λ (t)), żλ (t)
	(42)			= 2 z = -2 λ	d 2

→0 d 2 S(t+ ) (z λ (t + )) -d 2 S(t) (z λ (t + )) + d 2 S(t) (z λ (t + )) -d 2 S(t) (z λ (t)) . λ (t) -proj(z λ (t), S(t)), żλ (t) λ (t) + 2 z λ (t) -proj(z λ (t), S(t)), f (t, z λ (t)) . t 0 e -(t-s)/λ (L f + L f |z λ (s)| + L S ) ds L f + L f |z λ (t)| + L f λ t 0 e -(t-s)/λ ds+ L f λ t 0 e -(t-s)/λ |z λ (s)| ds + L S λ t 0 e -(t-s)/λ ds. t 0 e -(t-s)/λ ds L S . λ : [t 1 , t 2 ] →

R n such that ϕ λ (s) ∈ S λ (s) for each s ∈ [t 1 , t 2 ] and by construction, we have ϕ λ (t 1 ) ∈ S λ (t 1 ) = S(t 1 ). By taking λ → 0, we see that z λ → z; and letting

This work was supported by the ANR project ConVan with grant number ANR-17-CE40-0019-01. 1


Adrien Pavao

Zhengying Liu

Isabelle Guyon

Filtering participants improves generalization in competitions and benchmarks

We address the problem of selecting a winning algorithm in a challenge or benchmark. While evaluations of algorithms carried out by third party organizers eliminate the inventor-evaluator bias, little attention has been paid to the risk of over-fitting the winner's selection by the organizers. In this paper, we carry out an empirical evaluation using the results of several challenges and benchmarks, evidencing this phenomenon. We show that a heuristic commonly used by organizers consisting of pre-filtering participants using a trial run, reduces over-fitting. We formalize this method and derive a semi-empirical formula to determine the optimal number of top k participants to retain from the trial run.

Introduction

Competitions and benchmarks in Artificial Intelligence (AI), Machine Learning (ML), and other scientific domains, also called challenges or contests, have been playing an increasing role in research and development. Hundreds of competitions are organized each year on platforms such as Kaggle, AICrowd, CodaLab, and others. Competitions are now part of most major conferences and their protocols are peer reviewed, to increase the rigor of evaluations and check such things as sufficient test set sizes and possible bias in data. The code of the winners is generally open-sourced, to maximize their impact. One example often cited is the start of the "deep learning era" usually traced back to the success of AlexNet in the ImageNet benchmark [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF].

The ambition of competitions is generally to recommend algorithms that could perform well on new tasks resembling that of the competition. Thus competitions are a problem in which the organizers perform a learning task: from the tasks of the challenge, they select an algorithm that should perform well on new future tasks. Organizer over-fitting occurs when the number of participants is large and rankings are noisy, increasing the chance of poorly selecting a winner. Competition organizers face a sad paradox: the larger the number of participants, the more "successful" their competition, but also the greater the risk to overfit the particular competition setting.

A heuristic often employed in sports, chess, and other types of competition is to use eliminatory trial runs to filter participants for the final competition phase. In this paper, we investigate this strategy in the context of 2-phase machine learning competitions. To alleviate over-fitting the challenge setting, we propose to use the development phase to filter participants before entering the final phase.

The algorithm selection problem can be traced back at least as far as 1976 [START_REF] Rice | The algorithm selection problem[END_REF]. Its relationship with the problem of meta-learning and the risk of over-fitting at the meta-level has been pointed out [START_REF] Liu | Automated Deep Learning : Principles and Practice[END_REF]. However, to the best of our knowledge, this is the first paper that considers judging competitions or conducting a benchmark as an algorithm selection problem, and considers the risk of organizer over-fitting. This paper addresses the problem of determining a single winner from rankings already obtained in two-phase competitions. We do not address here the process of generating a sound and reproducible ranking. However, many related problems have been investigated in the literature: the test set size needed to get good error rate estimations [START_REF] Guyon | What size test set gives good error rate estimates[END_REF], score distributions for stochastic algorithms [START_REF] Reimers | Reporting score distributions makes a difference: Performance study of LSTM-networks for sequence tagging[END_REF], the selection of the worst run to reduce chance in competitions [START_REF] Baz | Advances in metadl: Aaai 2021 challenge and workshop[END_REF], or the problem of fusing scores from multiple "judges" (multiple tasks and/or multiple metrics) [START_REF] Pavao | Judging competitions and benchmarks: a candidate election approach[END_REF].

A question related to the problem we address is that of "out of domain generalization", which has been addressed from a variety of standpoints (although not in the context of judging competitions). For example, recent work has been warning about the danger that meta-learning algorithms may overfit the set of meta-training tasks [START_REF] Shen | Towards out-of-distribution generalization: A survey[END_REF]. Propositions to reduce over-fitting at the "meta-level" include algorithm-specific regularization methods like dropout [START_REF] Tseng | Regularizing meta-learning via gradient dropout[END_REF], Bayesian mechanisms [START_REF] Yoon | Bayesian model-agnostic meta-learning[END_REF], and meta-augmentation of datasets [START_REF] Rajendran | Meta-learning requires meta-augmentation[END_REF] (which could easily be applied to challenges). Learning-theoretic bounds on the meta-generalization gap have also been proposed [START_REF] Amit | Meta-learning by adjusting priors based on extended pac-bayes theory[END_REF], introducing prior knowledge through setting an experience-dependent prior for novel task. This relates to the algorithm proposed in this paper: we use the development phase as a prior to select the winner in the final phase.

Proposed algorithm: top-k

In the sequel, we consider multiple phase competitions or benchmarks, in which participants are ranked multiple times. We assume that such rankings are drawn i.i.d. from a distribution of rankings. This abstracts from the need of describing more precisely how such rankings are obtained, which may widely vary from one competition to another. We revisit the i.i.d. assumption in Section 5.

In a two-phase competition, the development phase is the "trial run" during which participant are filtered, and the subsequent final phase is used to determine the winner. The proposed algorithm (Algorithm 1) retains the top k participants in the development phase, and then selects the winner among them on the basis of the best rank in the final phase.

Clearly, keeping only a few participants in the final phase may save time and resources and may be logistically necessary if the final phase is a live competition. The question addressed in the remainder of the paper is the following: does the top-k algorithm yield as good or better (meta-)generalization than the "vanilla" method of simply selecting the winner in the final phase, without pre-filtering? Algorithm 1: top-k method Let: n ∈ N, n ≥ 1, be the number of participants, k ∈ N, 1 ≤ k ≤ n, be a chosen quota, rankD(i) be the rank of participant i in the development phase, rankF (i) be the rank of participant i in the final phase, i = 1 • • • n.

Select the winner as: i * = argmin rankD(i)≤k rankF (i)

Empirical results

Although, from the participants' point of view, their algorithms are being tested in the final phase, from the organizers' point of view the final phase is used for (meta-)training1 . Hence we need a "post-challenge phase" for (meta-)testing, to evaluate meta-generalization, if we do not have direct access to the distribution of rankings (which is only possible for synthetic data). We plot meta-learning curves as a function of k, the number of participants pre-selected in the top-k of the development phase:

-Meta-training error: score in the final phase of the declared winner.

-Meta-test error: score in the post-challenge phase of the declared winner.

-Generalization gap: difference between meta-test and meta-training error.

Experiments on real data

We ran the top-k algorithm on 4 real meta-datasets coming from previously run challenges or benchmarks: AutoDL, AutoML, OpenML and Statlog, provided by [START_REF] Sun-Hosoya | Activmetal: Algorithm recommendation with active meta learning[END_REF]. We simulate phases by averaging the scores on subsets of tasks. The learning curves shown in Figure 1 were obtained by averaging over 10000 datasplit trials. We first observe that the meta-training error monotonically decreases as a function of k. This is to be expected since the meta-training error is the minimum score in the final phase of the top-k in the development phase: as k increases, one should get smaller and smaller scores by taking the minimum over a larger subset of values. The meta-training error basically varies with 1 k [START_REF] Gordon | On the minimum of several random variables[END_REF]. We also observe that the meta-test error goes through a minimum (desired behavior). The last point represents the vanilla method consisting in selecting the candidate with the smallest score in the final phase. We remark that the first and last point have similar values. This is not surprising since the first point is analogous to the vanilla method, but using the development phase to select the winner instead of the final phase.

Experiments on synthetic data

To gain more insight into the problem, and understand under which condition the top-k method allows us to meta-generalize better than simply selecting the winner in the final phase, we generated synthetic data, as follows: Consider a competition with n ∈ N participants. Assume that there is an ideal true (but unknown) ranking of participants g (g is for "generalization"). A synthetic empirical ranking obtained in a challenge phase is generated from g by repeated permutations of pairs of neighbors: a position i is drawn at random from {1, • • • , n -1} and the participants i and i + 1 are exchanged. We repeat this operation s times. The smaller s, the more the empirical rankings will be correlated to the true ranking g (and to one another). We generate in this way 3 fake participant rankings, one for each phase: D, F , and P . We can then perform similar experiments as with real data and compute meta-training error and meta-test error of the top-k method.

We show in Figure 1 the results of experiments on synthetic data next to those on real data. On each row, we took care of matching the distance d between phase rankings (as measured by the Kendall τ ). We are pleased to see that qualitatively the real and synthetic data curves remarkably resemble one another: monotonic decrease of meta-training curve; meta-test curve going through an optimum; first and last point of meta-test curve nearly identical (performance of vanilla method). We observe that the optimum value of k is relatively small, even for large numbers of participants (last row). It increases with d, the distance between rankings. While s (the number of swaps applied to the ideal ranking, in synthetic data) can go to infinity, the expected value of d is bounded by half the Kendall τ distance between a ranking and its reverse, i.e. n(n-1)

4

. Noting that the random process to create artificial rankings is essentially a random walk of the winner, we propose a semi-empirical formula for the optimum of k2 :

k * ≃ 1 + d n
We validate the formula using simulations3 .

Discussion and conclusion

The remarkable resemblance between real and synthetic curves when the Kendall τ distance d between phase rankings is matched, allowed us to propose a semiempirical formula to predict the optimal value of the number of challenge participants k selected during the "trial run": k * ≃ 1 + d n . However, this requires knowing d a priori or estimating it. A more conservative choice of k may anyway be preferable in practice, because when rankings are "noisy" (d large), there is a lot of variance (not visible in our experimental results, which are averaged over 10000 trials). So, instead of using k * , we advocate eliminating participants who do not outperform the baseline methods provided by the organizers with the "starting kit" (which may include well performing methods from previous challenges). This should be more acceptable to the participants than setting a hard threshold on the number of entrants of the final phase, and will at least eliminate the least serious participants who just submit the "starting kit". The validity of our paradigm rests on the assumptions that rankings in various phases are drawn i.i.d. from a distribution of rankings, and without ties. This assumption may be violated if, for instance, the tasks of the development phase and final phase are of different nature or difficulty, or if the participants overfit the tasks of the development phase (e.g. by making many submissions and seeking feed-back from a leaderboard). This last problem does not seem to be severe in practice [START_REF] Roelofs | A meta-analysis of overfitting in machine learning[END_REF].

Further work include: Extending to 3-best selection (instead of winner selection); Handling ranking with ties; Combining top-k with meta-cross-validation; Other winner selection methods; Considering more than two phases; Finalising the formal analysis of the synthetic problem provided in supplemental material.

Fig. 1 :

 1 Fig. 1: Learning curves with top-k method (averaged over 10000 trials): Meta-training error, Meta-test error, Generalization gap. Rank error is the rank r of the selected winner in F and P, normalized by (r -1)/(n -1). The blue vertical bar marks the optimum, s is the number of random swaps, d is the average Kendall τ distance [6] between phase rankings, and the grey vertical bar marks k * ≃ 1 + d n , an estimation of the optimal value of k.

The development phase during which the participants "practice" (i.e. perform a form of training) is used as a prior by the organizers.

Although regular random walk calculations are simple, the boundary conditions imposed in this problem increase the difficulty of obtaining a simple exact formula.

The code to reproduce all experiments and supplemental material can be found at https://github.com/didayolo/metagen-esann-2022.