Eurovis

/ S Lenti

M Krone

J Schmidt

D Benvenuti

G Fiordeponti

H Cheng

T Catarci

J D Fekete

G Santucci

M Angelini

L Battle

Toward an Interaction-Driven Framework for Modeling Big Data Visualization Systems

Keywords: CCS Concepts, Information systems → Information systems applications;, Human-centered computing → Human computer interaction (HCI);

Designing big data visualization applications is challenging due to their complex yet isolated development. One of the most common issues is an increase in latency that can be experienced while interacting with the system. There exists a variety of optimization techniques to handle this issue in specific scenarios, but we lack models for integrating them in a holistic way, hindering the integration of complementary functionality and hampering consistent evaluation across systems. In response, we present a framework for modeling the big data visualization pipeline which builds a bridge between the Visualization, Human-Computer Interaction, and Database communities by integrating their individual contributions within a single, easily interpretable pipeline. With this framework, visualization applications can become aware of the full end-to-end context, making it easier to determine which subset of optimizations best suits the current context.

Introduction

When users interact with large data in a visualization application, its response time is crucial in keeping users efficient and engaged, particularly for exploratory data analysis (EDA). Even latencies as low as 500 ms can be detrimental [LH14, WK09]. Latencies can be introduced in each of the three layers of the visualization pipeline: (i) data layer; (ii) interaction layer and (iii) rendering layer. For example, for large data and complex queries, even commercial database systems (DBMS) may fail to meet the performance needs of interactive visual analysis systems [BEA * 20], introducing latencies in the data layer. In practice, most visualization applications are designed for small data or resort to ad-hoc mechanisms to try to meet the latency constraints on larger data (e.g., Spotfire [Spo95]) by optimizing only one of the three layers. Although we have good models for optimizing the data and rendering layers (e.g., BIRCH [ZRL96], DEVise [LRB * 97] etc.), we lack one for the interaction layer. Furthermore, no formal framework exists to connect known performance models in the visualization and HCI communities with optimization strategies at the data management level. To fill this gap between these communities we present a framework for modeling how user interactions affect each component of the data visualization pipeline. The main idea of the framework is to introduce a translation layer and to model interactions using augmented statecharts [Har87]. Users' interactions are captured and then translated from low-level user actions into highlevel database queries and rendering information that are annotated in the statechart. The annotated statechart can be explored by the designers to understand which interactions suffer from excessive latency and how to fix it. By applying this framework, it becomes possible to: (i) have a complete model of how users interact with the visualization application that can be annotated with information useful for optimization (e.g., latency constraints, probabilities); (ii) know in advance the SQL queries that could be triggered from the current interaction state and (iii) derive through them the optimized computations needed to (re)render the results.

Framework

The proposed framework, visible in Figure 1, is built on the classic visualization pipeline, modeled into three different layers: (i) a data layer, which collects the relevant data through selection and applies optional binning and aggregation operations; (ii) a rendering layer that uses the binned and aggregated data to generate visualization components and computes the extent of each rendered column to produce axes and scales and (iii) an interaction layer in charge of issuing queries either from dynamic query widgets or direct manipulation of the visualization components. Big data exploration is generally made possible through the cooperation of these layers, but each application is designed to address separate parts of the pipeline. Without a connecting framework, it becomes difficult to understand how a new DBMS optimization strategy applies to new and existing visualization interface designs. As a result, optimizations developed in the database community run the risk of being point solutions at best, and irrelevant at worst in solving the performance issues observed by the visualization community. One of the key ideas behind this framework is the addition of a translation layer to integrate critical information from the other three. It (i) receives query fragments (from the user's activity) from the interaction layer, (ii) combines and transforms them, adding meaningful information about the context in which the interaction happened and (iii) translates them into SQL queries and rendering computations. By specifying interactions as state machines with latency information and leveraging the translation layer, our unified framework enables more versatile management of query load driven by the user's actions and provides a centralized layer of computation that allows us to exploit, rather than modify, the existing interaction and data management layers. We envision this interaction-driven framework to be supported, based on its functionalities, by three different components: (i) Augmented Statechart, (ii) SQL Translator and (iii) Rendering Translator. In the following, we focus on the ongoing work for the first two components.

(i) We define an atomic interaction as a statechart containing the minimum states and transitions needed to represent an interaction, which is typically one state and three transitions for stateful actions or two transitions only for stateless ones. Examples of atomic interactions are: hover, click, zoom, drag, pan and brush. We can combine atomic interactions to model more complex interactions (e.g., "hover" and "click" to model the interaction with a button).

(ii)A user interaction is captured by the interaction layer and sent to the translation layer which extracts the SQL queries to be executed. This feature is implemented through the Label Module (LM), which gives the result back by labeling the statechart. The generations of SQL prepared statement is implemented by the Query Module (QM) which manages a Query Storage (QS), i.e. a data structure containing SQL Fragments. When the statechart is initialized, LM assigns a specific set of query parameters to each of its transitions that contribute to the SQL Statements. Those parameters are referred to as Query Metadata, organized in a set of query labels. When a transition leads to a state containing a Query Metadata, that information is acquired by QM: for each Query Label, a SQL Prepared Statement is created by taking from QS all the different fragments and filled with dimensional values and parameters.

Usage-Scenario

We started experimenting with the proposed framework on a demanding yet popular visualization application, crossfilter [MHH19a]. We modeled the interactions by combining the JSON definitions of atomic interactions and running the statechart using the XState JavaScript library [XSt22], which is used to build the main components of the translation layer. When a user's action is captured by a listener, we pass the information to Xstate which updates the statechart accordingly and labels it with query fragments. Finally, we created a dashboard to show in real-time: (i) each state reachable in n steps from the current one and (ii) each query fragment that can be generated in those steps. A snapshot of this dashboard can be seen in Figure 2 in which we captured its state while the user is idling while brushing on the "Distance in Miles" histogram. The interaction layer collects information about the current state in the statechart, the last transition triggered by the user, the list of two steps paths branching from the current state (e.g., we can reach the "inactive" state if the "MOUSELEAVE" event occurs) and metadata about the context in which the interaction is happening (e.g., the active view and the values of the current brush). The translation layer is updated in real-time and generates the Query Labels (see section 2), visible in the "Active Meta" field of the dashboard.

Conclusions

The proposed interaction-driven framework can enable a higher awareness, both at design and execution time, of the existing interrelations among the four layers considered (data, interaction, rendering, and the newly introduced translation). It can accommodate also different optimization strategies like all the ones listed in ), where the key challenge lies in identifying which queries to actually materialize; User Modeling (e.g., Fore-Cache [START_REF] Battle | Dynamic Prefetching of Data Tiles for Interactive Visualization[END_REF]) where the key challenge is to develop an accurate model of user interaction behavior used to predict future interactions.

Figure 1 :

 1 Figure1: The interaction-driven framework for big data exploration (the translation layer is highlighted in orange).

Figure 2 :

 2 Figure 2: A snapshot of the dashboard obtained applying the Interaction-Driven Framework to a crossfilter interface. The user is currently brushing on the "Distance in Miles" histogram.

  Battle and Scheidegger work [BS20]. As examples, we cite Multi-Query Optimization (e.g., SeeDB [VRM * 15], Zenvisage [SKL * 16], and Voyager [WMA * 16]) producing many simultaneous queries; , Materialized Views (e.g., Falcon [MHH19b], imMens [LJH13], and Nanocubes [LKS13]

© 2022 The Author(s) Eurographics Proceedings © 2022 The Eurographics Association.


Matteo Tacchi
email: matteo.tacchi@epfl.ch

Jean Bernard Lasserre

Didier Henrion

Stokes, Gibbs and volume computation of semi-algebraic sets

Keywords: Numerical methods for multivariate integration, real algebraic geometry, convex optimization, Stokes' theorem, Gibbs phenomenon

We consider the problem of computing the Lebesgue volume of compact basic semialgebraic sets. In full generality, it can be approximated as closely as desired by a converging hierarchy of upper bounds obtained by applying the Moment-SOS (sums of squares) methodology to a certain infinite-dimensional linear program (LP). At each step one solves a semidefinite relaxation of the LP which involves pseudo-moments up to a certain degree. Its dual computes a polynomial of same degree which approximates from above the discontinuous indicator function of the set, hence with a typical Gibbs phenomenon which results in a slow convergence of the associated numerical scheme. Drastic improvements have been observed by introducing in the initial LP additional linear moment constraints obtained from a certain application of Stokes' theorem for integration on the set. However and so far there was no rationale to explain this behavior. We provide a refined version of this extended LP formulation. When the set is the smooth super-level set of a single polynomial, we show that the dual of this refined LP has an optimal solution which is a continuous function. Therefore in this dual one now approximates a continuous function by a polynomial, hence with no Gibbs phenomenon, which explains and improves the already observed drastic acceleration of the convergence of the hierarchy. Interestingly, the technique of proof involves recent results on Poisson's partial differential equation (PDE).

Introduction

Consider the problem of computing the Lebesgue volume λ(K) of a compact basic semi-algebraic set K ⊂ R n . For simplicity of exposition we will restrict to the case where K is the smooth superlevel set { x : g(x) ≥ 0 } ⊂ R n of a single polynomial g.

If K is a convex body then several procedures are available; see e.g. exact deterministic methods for convex polytopes [START_REF] Büeler | Exact Volume Computation for Polytopes: A Practical Study[END_REF], or non deterministic Hit-and-Run methods [START_REF] Smith | The hit-and-run sampler: A globally reaching Markov chain sampler for generating arbitrary multivariate distribution[END_REF][START_REF] Zabinsky | Hit-and-Run methods[END_REF] and the more recent [START_REF] Cousins | A practical volume algorithm[END_REF][START_REF] Dadush | Near-optimal deterministic algorithms for volume computation via M-ellipsoids[END_REF]. Even approximating λ(K) by deterministic methods is still a hard problem as explained in e.g. [START_REF] Dadush | Near-optimal deterministic algorithms for volume computation via M-ellipsoids[END_REF] and references therein. In full generality with no specific assumption on K such as convexity, the only general method available is Monte-Carlo, that is, one samples N points according to Lebesgue measure λ normalized on a simple set B (e.g. a box or an ellipsoid) that contains K. If ρ N is the proportion of points that fall into K then the random variable ρ N λ(B) provides a good estimator of λ(K) with convergence guarantees as N increases. However this estimator is non deterministic and neither provides a lower bound nor an upper bound on λ(K).

When K is a compact basic semi-algebraic set, a deterministic numerical scheme described in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] provides a sequence (τ k ) k∈N ⊂ R of upper bounds that converges to λ(K) as k increases. Briefly,

λ(K) = inf p∈R[x] p dλ : p ≥ 1 K on B (1)
τ k = inf p∈R[x] k p dλ : p ≥ 1 K on B , (2)
with x → 1 K (x) = 1 if x ∈ K and 0 otherwise. One can notice that minimizing sequences for (1) and ( 2) also minimize the L 1 (B, λ)-norm p -1 K 1 (with convergence to 0 in the case (1)).

As the upper bound τ k > λ(K) is obtained by restricting the search in [START_REF] Cousins | A practical volume algorithm[END_REF] to polynomials of degree at most k, the infimum is attained and an optimal solution can be obtained by solving a semidefinite program. Of course, the size of the resulting semidefinite program increases with the degree k: this is the so-called Moment-SOS hierarchy; for more details the interested reader is referred to [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF]. Also focusing on compact semi-algebraic sets, [START_REF] Lairez | Computing the volume of compact semi-algebraic sets[END_REF] proposes a symbolic method to compute the volume of K with absolute precision 2 -p , in time O(p log(p) 3+ε ) for any ε > 0 as p → ∞. This is in sharp contrast with the approach considered here, which consists in approximating problem [START_REF] Büeler | Exact Volume Computation for Polytopes: A Practical Study[END_REF] with the sequence of problems (2) indexed by k. Indeed,

• In [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF], for any k ∈ N, τ k is guaranteed to be a converging upper bound for λ(K), i.e.

τ k -λ(K) > 0, while [START_REF] Lairez | Computing the volume of compact semi-algebraic sets[END_REF] also guarantees convergence of the approximant to λ(K) but gives no information on the sign of the difference between the two quantities.

• [START_REF] Lairez | Computing the volume of compact semi-algebraic sets[END_REF] uses symbolic computations that can achieve arbitrary precision, while [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] uses numerical computations based on semidefinite programming, limited to floating-point arithmetic precision.

• The approach of [START_REF] Lairez | Computing the volume of compact semi-algebraic sets[END_REF] can be used to approximate other quantities than the volume, namely real periods of algebraic surfaces. The approach of [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] was extended to approximate sets relevant in systems control, such as regions of attraction or maximal positively invariant sets (see e.g. [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF]). In this context, the present contribution can help improving the Moment-SOS hierarchy for assessing the stability of polynomial differential systems.

When solving problem [START_REF] Cousins | A practical volume algorithm[END_REF], clearly a Gibbs phenomenon1 takes place as one tries to approximate on B and from above, the discontinuous function 1 K by a polynomial of degree at most k. This makes the convergence of the upper bounds τ k very slow (even for modest dimension problems). A trick was used in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] to accelerate this convergence but at the price of loosing monotonicity of the resulting sequence.

In fact [START_REF] Büeler | Exact Volume Computation for Polytopes: A Practical Study[END_REF] is a dual of the following infinite-dimensional Linear program (LP) on measures

sup µ { µ(K) : µ ≤ λ ; µ ∈ M(K) + } (3)
(where M(K) + is the space of finite Borel measures on K). Its optimal value is also λ(K) and is attained at the unique optimal solution µ := λ K = 1 K λ (the restriction of λ to K).

A simple but key observation. As one knows the unique optimal solution µ = λ K of (3), any constraint satisfied by µ (in particular, linear constraints) can be included as a constraint on µ in (3) without changing the optimal value and the optimal solution. While these constraints provide additional restrictions in (3), they translate into additional degrees of freedom in the dual (hence a relaxed version of (1)), and therefore better approximations when passing to the finite-dimensional relaxed version of [START_REF] Cousins | A practical volume algorithm[END_REF]. A first set of such linear constraints experimented in [START_REF] Lasserre | Computing Gaussian and exponential measures of semi-algebraic sets[END_REF] and later in [START_REF] Lasserre | Computing the Hausdorff boundary measure of semi-algebraic sets[END_REF], resulted in drastic improvements but with no clear rationale behind such improvements.

Contribution. The main message and result of this paper is that there is an appropriate set of additional linear constraints on µ in (3) such that the resulting dual (a relaxed version of ( 1)) has an explicit continuous optimal solution with value λ(K). These additional linear contraints (called Stokes constraints) come from an appropriate modelling of Stokes' theorem for integration over K, a refined version of that in [START_REF] Lasserre | Computing Gaussian and exponential measures of semi-algebraic sets[END_REF]. Therefore the optimal continuous solution can be approximated efficiently by polynomials with no Gibbs phenomenon, by the hierarchy of semidefinite relaxations defined in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] (adapted to these new linear constraints). Interestingly, the technique of proof and the construction of the optimal solution invoke results from the field of elliptic partial differential equations (PDE), namely a recent extension of standard Schauder estimates from Dirichlet problems to Neumann formulations.

Outline. In Section 2 we recall the primal-dual linear formulation of the volume problem, and we explain why the dual value is not attained, which results in a Gibbs phenomenon. In Section 3 we revisit the acceleration strategy based on Stokes' theorem, with the aim of introducing in Section 4 a more general acceleration strategy and a new primal-dual linear formulation of the volume problem. Our main result, attainment of the dual value in this new formulation, is stated as Theorem 4.2 at the end of Section 4. The drastic improvement in the convergence to λ(K) is illustrated on various simple examples.

Linear reformulation of the volume problem

Consider a compact basic semi-algebraic set

K := {x ∈ R n : g(x) ≥ 0} with g ∈ R[x].
We suppose that K ⊂ B where B is a compact basic semi-algebraic set for which we know the moments B x k dx of the Lebesgue measure λ B , where

x k := x k1 1 x k2 2 • • • x kn n denotes a multivariate monomial of degree k ∈ N n . We assume that Ω := {x ∈ R n : g(x) > 0} is a nonempty open set with closure 2 Ω = K,
and that its boundary

∂Ω = ∂K = K \ Ω is C 1 in
the sense that it is locally the graph of a continuously differentiable function. We want to compute the Lebesgue volume of K, i.e., the mass of the Lebesgue measure λ K :

λ(K) := K dx = R n dλ K (x).
If X ⊂ R n is a compact set, denote by M(X) the space of signed Borel measures on X, which identifies with the topological dual of C 0 (X), the space of continuous functions on X. Denote by M(X) + the convex cone of non-negative Borel measures on X, and by C 0 (X) + the convex cone of non-negative continuous functions on X.

In [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] a sequence of upper bounds converging to λ(K) is obtained by applying the Moment-SOS hierarchy [START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF] (a family of finite-dimensional convex relaxations) to approximate as closely as desired the (primal) infinite-dimensional LP on measures:

max µ µ(K) (4)
s.t. µ ∈ M(K) + λ B -µ ∈ M(B) +
whose optimal value is λ(K), attained for µ := λ K . The LP (4) has an infinite-dimensional LP dual on continuous functions which reads:

inf w B w dλ (5)
s.t. w ∈ C 0 (B) + w| K -1 ∈ C 0 (K) + .
Observe that (5) consists of approximating the discontinuous indicator function 1 K (equal to one on K and zero elsewhere) from above by continuous functions w, in minimizing the L 1 (B)-norm w -1 K 1 . Clearly the infimum λ(K) is not attained. Since K is generated by a polynomial g, and measures on compact sets are uniquely determined by their moments, one may apply the Moment-SOS hierarchy [START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF] for solving [START_REF] Dahl | Semidefinite optimization using MOSEK[END_REF]. The moment relaxation of (4) consists of replacing µ by finitely many of its moments y, say up to degree d ∈ N. Then the cone of moments is relaxed by a linear slice of the semidefinite cone constructed from so-called moment and localizing matrices indexed by d, as defined in e.g. [START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF], and which defines a semidefinite program. Therefore the dual of this semidefinite program (i.e., the dual SOS-hierarchy) is a strengthening of ( 5) where (i) continuous functions w are replaced with polynomials of increasing degree d, and (ii) nonnegativity constraints are replaced with Putinar's SOS-based certificates of positivity [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF] which translate to semidefinite constraints on the coefficients of polynomials; again the interested reader is referred to [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF][START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF] for more details.

For each fixed degree d, a valid upper bound on λ(K) is computed by solving a primal-dual pair of convex semidefinite programming problems (not described here). As proved in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] by combining Stone-Weierstrass' theorem and Putinar's Positivstellensatz [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF],

(i) there is no duality gap between each primal semidefinite relaxation of the hierarchy and its dual, and

(ii) the resulting sequence of upper bounds converges to λ(K) as d increases.

The main drawback of this numerical scheme is its typical slow convergence, observed already for very simple univariate examples, see e.g. [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF]Figs. 4.1 and 4.5]. The best available theoretical convergence speed estimates are also pessimistic, with an asymptoptic rate of log log d [START_REF] Korda | Convergence rates of moment-sum-of-squares hierarchies for volume approximation of semialgebraic sets[END_REF]. Slow convergence is mostly due to the so-called Gibbs phenomenon which is well-known in numerical analysis [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]Chapter 9]. Indeed, as already mentioned, solving (5) numerically amounts to approximating the discontinuous function 1 K from above with polynomials of increasing degree, which generates oscillations and overshoots and slows down the convergence, see e.g. [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF]Figs. 4.2,4.4,4.6,4.7,4.10,4.12].

Example 1. Let K := [0, 1/2] ⊂ B := [-1, 1].
In Figure 1 is displayed the degree-10 and degree-20 polynomials w obtained by solving the dual of SOS strengthenings of problem (4). We can clearly see bumps, typical of a Gibbs phenomenon at points of discontinuity. An idea to bypass this limitation consists of adding certain linear constraints to the finitedimensional semidefinite relaxations, to make their optimal values larger and so closer to the optimal value λ(K). Such linear constraints must be chosen appropriately:

(i) they must be redundant for the infinite-dimensional moment LP on measures (4), and (ii) become active for its finite-dimensional relaxations. This is the heuristic proposed in [START_REF] Lasserre | Computing Gaussian and exponential measures of semi-algebraic sets[END_REF] to accelerate the Moment-SOS hierarchy for evaluating transcendental integrals on semi-algebraic sets. These additional linear constraints on the moments y of µ are obtained from an application of Stokes' theorem for integration on K, a classical result in differential geometry. It has been also observed experimentally that this heuristic accelerates significantly the convergence of the hierarchy in other applied contexts, e.g. in chance-constrained optimization problems [START_REF] Weisser | Computing Approximations and Generalized Solutions using Moments and Positive Polynomials[END_REF].

Introducing Stokes constraints

In this section we explain the heuristic introduced in [START_REF] Lasserre | Computing Gaussian and exponential measures of semi-algebraic sets[END_REF] to accelerate convergence of the Moment-SOS hierarchy by adding linear constraints on the moments of µ . These linear constraints are obtained from a certain application of Stokes' theorem for integration on K.

Stokes' Theorem and its variants

Theorem 3.1 (Stokes' Theorem). Let Ω ⊂ R n be a C 1 open set 3 with closure K. For any (n -1)-differential form ω on K, it holds ∂Ω ω = Ω dω. Corollary 3.2. In particular, for u ∈ C 1 (K) n and ω(x) = u(x) • n Ω (x) dσ(x)
, where the dot is the inner product, σ is the surface or Hausdorff measure on ∂Ω and n Ω is the outward pointing normal to ∂Ω, we obtain the Gauss formula

∂Ω u(x) • n Ω (x) dσ(x) = Ω div u(x) dx. ( 6
)
With the choice u(x) := u(x) e i where u ∈ C 1 (K) and e i is the vector of R n with one at entry i and zeros elsewhere, for i = 1, . . . , n, we obtain the dual Gauss formula

∂Ω u(x) n Ω (x) dσ(x) = Ω grad u(x) dx. (7)
Proof. These are all particular cases of [10, Theorem 6.10.2].

Original Stokes constraints

Associated to a sequence y = (y k ) k∈N n ∈ R N n , introduce the Riesz linear functional L y : R[x] → R which acts on a polynomial p := k p k x k ∈ R[x] by L y (p) := k p k y k .
Thus, if y is the sequence of moments of λ K , i.e. y k := K x k dx for all k ∈ N n , then L y (p) = K p(x)dx and by [START_REF] Henrion | The Moment-SOS Hierarchy -Lectures in Probability[END_REF] with u(x) := x k g(x):

L y (grad(x k g)) = K grad(x k g(x)) dx = ∂K x k g(x) n K (x) dσ(x) = 0,
since by construction g vanishes on ∂K. Thus while in the infinite-dimensional LP (4) one may add the linear constraints

K grad(x k g) dµ = 0 ∀k ∈ N n ,
without changing its optimal value λ(K), on the other hand inclusion of the linear moment constraints

L y (grad(x k g)) = 0 , |k| ≤ 2d + 1 -deg(g) (8)
in the moment relaxation with pseudo-moments y of degree at most d, will decrease the optimal value of the initial relaxation. In practice, it was observed that adding constraints (8) dramatically speeds up the convergence of the Moment-SOS hierarchy, see e.g. [START_REF] Lasserre | Computing Gaussian and exponential measures of semi-algebraic sets[END_REF][START_REF] Weisser | Computing Approximations and Generalized Solutions using Moments and Positive Polynomials[END_REF]. One main goal of this paper is to provide a qualitative mathematical rationale behind this phenomenon.

Infinite-dimensional Stokes constraints

In [START_REF] Tacchi | Exploiting sparsity in semi-algebraic set volume computation[END_REF], Stokes constraints were formulated in the infinite-dimensional setting, and a dual formulation was obtained in the context of the volume problem. Using ( 6) with u = gv (which vanishes on ∂K) and v ∈ C 1 (K) n arbitrary, yields:

K (grad g(x) • v(x) + g(x) div v(x)) dx = ∂K gv n K dσ = 0 ,
which can be written equivalently (in the sense of distributions) as

(grad g)λ K -grad(gλ K ) = 0 .
This allows to rewrite problem (4) as

max µ µ(K) (9)
s.t. µ ∈ M(K) + λ B -µ ∈ M(B) + (grad g)µ -grad(gµ) = 0
without changing its optimal value λ(K) attained at µ = λ K . Using infinite-dimensional convex duality as in e.g. the proof of Theorem 2 in [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF], the dual of LP [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF]

reads inf v,w B w dλ (10)
s.t. v ∈ C 1 (K) n w ∈ C 0 (B) + w| K -div(gv) -1 ∈ C 0 (K) + .
Crucial observation. Notice that w in [START_REF] Hubbard | Vector Calculus, Linear Algebra and Differential Forms (A Unified Approach[END_REF] is not required to approximate 1 K from above anymore. Instead, it should approximate 1 + div(gv) on K and 0 outside K. Hence, provided that 1 + div(gv) = 0 on ∂K, w might be a continuous function for some well-chosen v ∈ C 1 (K) n , and therefore an optimal solution of (10) (i.e., the infimum is a minimum). As a result, the Gibbs phenomenon would disappear and convergence would be faster. The issue is then to determine whether the infimum in (10) is attained or not. And if not, are there other special features of problem [START_REF] Hubbard | Vector Calculus, Linear Algebra and Differential Forms (A Unified Approach[END_REF] that can be exploited to yield more efficient semidefinite relaxations ?

New Stokes constraints and main result

In the previous section, the Stokes constraint

K (v(x) • grad g(x) + g(x) div v(x)) dµ(x) = 0
or equivalently (in the sense of distributions)

(grad g)µ -grad(gµ) = 0 (11)
(with µ ∈ M(K) + being the Lebesgue measure on K) was obtained as a particular case of Stokes' theorem with u = gv in [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF]. Instead, we can use a more general version with u not in factored form, and also use the fact that ∀x ∈ ∂K, 0 = grad g

(x) = -| grad g(x)| n K (x) (here |y| := √ y • y is the n-dimensional Euclidean norm), to obtain K div u(x) dµ(x) = - ∂K u(x) • grad g(x) dν(x) ,
or equivalently (in the sense of distributions)

grad µ = (grad g)ν , (12)
with µ ∈ M(K) + being the Lebesgue measure on K and ν ∈ M(∂K) + being the measure having density 1/| grad g(x)| with respect to the (n -1)-dimensional Haussdorff measure σ on ∂K. The same linear equation was used in [START_REF] Lasserre | Computing the Hausdorff boundary measure of semi-algebraic sets[END_REF] to compute moments of the Hausdorff measure. In fact, equation ( 12) is a generalization of equation [START_REF] Korda | Convergence rates of moment-sum-of-squares hierarchies for volume approximation of semialgebraic sets[END_REF] in the following sense. [START_REF] Lairez | Computing the volume of compact semi-algebraic sets[END_REF], then µ also satisfies [START_REF] Korda | Convergence rates of moment-sum-of-squares hierarchies for volume approximation of semialgebraic sets[END_REF].

Lemma 4.1. If ν ∈ M(∂K) + is such that µ ∈ M(K) + satisfies
Proof. Equation ( 12) means that K div u(x) dµ(x)

+ ∂K u(x) • grad g(x) dν(x) = 0 for all u ∈ C 1 (K) n . In particular if u = gv for some v ∈ C 1 (K) n then (12) reads K (v(x) • grad g(x) + g(x) div v(x)) dµ(x) = 0 ,
which is precisely [START_REF] Korda | Convergence rates of moment-sum-of-squares hierarchies for volume approximation of semialgebraic sets[END_REF].

Hence we can incorporate linear constraints ( 12) on µ and ν, to rewrite problem (4) as

max µ,ν µ(K) (13)
s.t. µ ∈ M(K) + ν ∈ M(∂K) + λ B -µ ∈ M(B) + (grad g)ν -grad µ = 0
without changing its optimal value λ(K) attained at µ = λ K and ν = σ/| grad g|. Notice that LP [START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF] involves two measures µ and ν whereas LP (9) involves only one measure µ.

Next, by convex duality as in e.g. the proof of Theorem 2 in [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF], the dual of (13) reads

inf u,w B w dλ (14)
s.t. u ∈ C 1 (K) n w ∈ C 0 (B) + w| K -div u -1 ∈ C 0 (K) + -(u • grad g)| ∂K ∈ C 0 (∂K) + .
Our main result states that the optimal value of the dual ( 14) is attained at some continuous function

(w, u) ∈ C 0 (B) + × C 1 (K) n .
Therefore, in contrast with problem (5), there is no Gibbs phenomenon at an optimal solution of the (finite-dimensional) semidefinite strengthening associated with [START_REF] Lasserre | Convex optimization and parsimony of L p -balls representation[END_REF].

Let Ω i , i = 1, . . . , N denote the connected components of Ω, and let

m Ωi (g) := 1 λ(Ω i ) Ωi g dλ.
Theorem 4.2. In dual LP [START_REF] Lasserre | Convex optimization and parsimony of L p -balls representation[END_REF] the infimum is a minimum, attained at

w (x) := g(x) N i=1 1 Ωi (x) m Ωi (g) , x ∈ B , and u (x) := grad u(x) ,
where u solves the Poisson PDE

-∆u(x) = 1 -w (x), x ∈ Ω ∂ n u(x) = 0, x ∈ ∂Ω.
Remark 1. The Moment-SOS hierarchy associated to LPs ( 13) and ( 14) yields upper bounds for the volume. Theorem 4.2 is designed for these LPs but it has a straightforward counterpart for lower bound volume computation, obtained by replacing K with B \ Ω in the previous developments, i.e. computing upper bounds of λ(B \ Ω). However, two additional technicalities should then be considered:

• This work only deals with semi-algebraic sets defined by a single polynomial; actually, it immediately generalizes to finite intersections of such semi-algebraic sets, as long as their boundaries do not intersect (i.e. here K should be included in the interior of B): the constraints on boundaries should just be splitted between the boundaries of the intersected sets.

• This work heavily relies on the fact that the boundary of the considered set should be smooth; for this reason, computing lower bounds of the volume implies that one chooses a smooth bounding box B (typically a euclidean ball, ellipsoid or p ball), which rules out simple sets like the hypercube [-1, 1] n .

Upon taking into account these technicalities, Theorem 4.2 still holds, allowing to deterministically compute upper and lower bounds for the volume, with arbitrary precision. Of course in practice, one is limited by the performance of state-of-art SDP solvers.

Proof of main result

Theorem 4.2 is proved in several steps as follows:

• we show that the optimal dual solution satisfies a Poisson PDE;

• we study the Poisson PDE on a union of connected domains;

• we construct an explicit optimum for problem [START_REF] Lasserre | Convex optimization and parsimony of L p -balls representation[END_REF].

Equivalence to a Poisson PDE

u ∈ C 1 (K) n and h ∈ C 0 (K) + solving h = 0 on ∂Ω, (15a)
-div u = 1 -h in Ω, (15b)
u • n Ω = 0 on ∂Ω. (15c)
Proof. Let (u, h) solve [START_REF] Lasserre | Computing Gaussian and exponential measures of semi-algebraic sets[END_REF]. Using (15a), one can define

w(x) = h(x) if x ∈ K 0 if x ∈ B \ K.
Then (u, w) is feasible for ( 14) and one has

B w dλ = Ω h dλ (15b) = Ω (1 + div u) dλ (6)
= λ(Ω)

+ ∂Ω u • n Ω dσ (15c) = λ(Ω)
so that (u, w) is optimal. Conversely, let (u, w) be an optimal solution of problem [START_REF] Lasserre | Convex optimization and parsimony of L p -balls representation[END_REF]. We know that (µ , ν ) = (λ Ω , σ/| grad g|) is optimal for problem [START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF]. Then, the KKT optimality conditions ensure complementary slackness:

Ω (w| Ω -div u -1) dλ = 0, ( 16a
) ∂Ω u • grad g | grad g| dσ = 0. (16b)
Since w| Ω -div u -1 is nonnegative, (16a) yields (15b) with h := w| Ω . Likewise, since -(u • grad g)| ∂Ω is nonnegative, (16b) yields (15c) and thus, using [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF], it holds Ω div u dλ = 0. Eventually, (16a) yields Ω w dλ = λ(Ω) = B w dλ by optimality of w, so that B\Ω w dλ = 0 and, since w is nonnegative, w| B\Ω = 0. Continuity of w finally allows to conclude that w = 0 on ∂Ω, which is exactly (15a).

From Lemma 5.1, existence of an optimum for ( 14) is then equivalent to existence of a solution to [START_REF] Lasserre | Computing Gaussian and exponential measures of semi-algebraic sets[END_REF], which we rephrase as follows, defining f := 1 -h and u = grad u with u ∈ C 2 (K), and where ∆u := div grad u is the Laplacian of u, and

∂ n u := grad u • n Ω . Lemma 5.2. If there exist u ∈ C 2 (K) n and f ∈ C 0 (K) solving -∆u = f in Ω, (17a)
∂ n u = 0 on ∂Ω, (17b)
f ≤ 1 in Ω, (17c)
f = 1 on ∂Ω, (17d)
then problem [START_REF] Lasserre | Convex optimization and parsimony of L p -balls representation[END_REF] has an optimal solution.

This rephrasing is a Poisson PDE (17a) with Neumann boundary condition (17b), whose source term f is a parameter subject to constraints (17c) and (17d).

Remark 2 (Loss of generality). Looking for u under the form u = grad u makes us loose the equivalence. Indeed, while ( 14) and ( 15) are equivalent, existence of a solution to [START_REF] Löfberg | YALMIP : A Toolbox for Modeling and Optimization in MATLAB[END_REF] is only a sufficient condition for existence of an optimum for [START_REF] Lasserre | Convex optimization and parsimony of L p -balls representation[END_REF], since [START_REF] Lasserre | Computing Gaussian and exponential measures of semi-algebraic sets[END_REF] might have only solutions u that are not gradients.

Remark 3 (Invariant set for gradient flow). From a dynamical systems point of view, the constraint in [START_REF] Lasserre | Convex optimization and parsimony of L p -balls representation[END_REF] which states that the inner product of u = grad u with grad g is non-positive on ∂Ω, means that we are looking for a velocity field or control u in the form of the gradient of a potential u such that K is an invariant set for the solutions t

∈ R → x(t) ∈ R n of the Cauchy problem ẋ(t) = -grad u(x(t)), x(0) ∈ B
after what we just have to define h := 1 + ∆u on Ω.

Regular solutions to the Poisson PDE

It remains to prove existence of solutions to problem [START_REF] Löfberg | YALMIP : A Toolbox for Modeling and Optimization in MATLAB[END_REF]. First, notice that PDE (17a) together with its boundary condition (17b) enforces an important constraint on the source term f , namely its mean must vanish:

Ω f dλ = 0. (18)
Indeed, if (f, u) solves ( 17), then

Ω f dλ (17a) = - Ω ∆u dλ (6) = - ∂Ω grad u • n Ω dσ (17b) = 0 .
Moreover, the following holds.

Lemma 5.3 (Existence and regularity on a connected domain). Suppose that Ω is connected. Let the source term f be Lipschitz continuous on K and have zero mean on Ω. Then there exists u ∈ C 2 (K) satisfying (17a) and (17b).

Proof. This is a direct application of [START_REF] Nardi | Schauder estimate for solutions of Poisson's equation with Neumann boundary condition[END_REF]: for α ∈ (0, 1), since K is bounded (let R > 0 be such that K ⊂ {x ∈ R n : x ≤ R}) and f is Lipschitz (let L be its Lipschitz constant on K), one has for x, y ∈ K that

|f (x) -f (y)| ≤ L x -y ≤ L x -y 1-α x -y α ≤ L (2R) 1-α <∞ x -y α , so that f ∈ C 0,α (K) := ϕ ∈ C 0 (K) : sup x,y∈K |ϕ(x) -ϕ(y)| x -y α < ∞ ,
and [START_REF] Nardi | Schauder estimate for solutions of Poisson's equation with Neumann boundary condition[END_REF] yields a solution

u ∈ C 2,α (K) := ϕ ∈ C 2 (K) : sup x,y∈K H(ϕ)(x) -H(ϕ)(y) x -y α < ∞
to the Poisson PDE (17a) with Neumann boundary condition (17b), where H(ϕ) = ∂ 2 ϕ ∂x i ∂x j i,j is the Hessian matrix of ϕ.

Remark 4. Assuming that ∂Ω is C ∞ instead of C 1 is actually without loss of generality since Ω is a semi-algebraic set: as soon as ∂Ω is locally the graph of a C 1 function, it is smooth.

In Lemma 5.3, we assumed that Ω is connected, so that we could apply the results of [START_REF] Nardi | Schauder estimate for solutions of Poisson's equation with Neumann boundary condition[END_REF]. To tackle non-connected sets, we recall that Ω is a semi-algebraic set, hence it has a finite number of connected components Ω 1 , . . . , Ω N . Moreover, the regularity of ∂Ω ensures that the Ω i have disjoint closures, so that we can trivially compute a solution u i on each Ω i and glue the u i together into a solution u := N i=1 1 Ωi u i on the whole Ω. Remark 5. Tackling the non-connected case requires that f has zero mean on each connected component of Ω: ∀i ∈ {1, . . . , N }, Ωi f dλ = 0. Remark 6. Lemma 5.3 automatically enforces -∆u = 1 on ∂Ω, which is crucial for the continuity of the optimization variable w.

Explicit optimum for volume computation with Stokes constraints

Our optimization problem does not feature only the Poisson PDE with Neumann condition: it also includes constraints (17c) and (17d) on the source term. Consequently, a Lipschitz continuous function f on K with zero integral over any connected component of Ω and satisfying (17c) and (17d) remains to be constructed. We keep the notations of Section 5.2 and suggest as candidate

x → f (x) := 1 -g(x) N i=1 1 Ωi (x) m Ωi (g) . ( 19
)
By definition, g = 0 on ∂Ω, so that (17d) automatically holds. Moreover, both g and 1 Ωi are nonnegative on K, so that (17c) also holds.

In terms of regularity, f is continuous and piecewise polynomial, so it is Lipschitz continuous on K.

Eventually, let i ∈ {1, . . . , N } so that Ω i is a connected component of Ω. Then, by definition, ∂Ω i ⊂ ∂Ω, and one has Ωi

f dλ = Ωi 1 -g(x) N i=1 1 Ωi (x) m Ωi (g) dx = λ(Ω i ) - 1 m Ωi (g) Ωi g(x) dx = 0,
by definition of m Ωi (g). This, together with Lemmata 5.2 and 5.3, concludes the proof of Theorem 4.2. Indeed, one can check that for the resulting w (x) = g(x)

N i=1

1 Ω i (x) m Ω i (g) , B w dλ = N i=1 1 m Ωi (g) Ωi g dλ = N i=1 λ(Ω i ) = λ(Ω) = λ(K).

Examples

To illustrate how efficient can be the introduction of Stokes constraints for volume computation, we consider the simple setting where K is a Euclidean ball included in B the unit Euclidean ball, as well as some basic variations around this case, where K is a non-euclidean ball or a union of balls, or B a non-euclidean ball. Indeed drastic improvements on the convergence are observed. All numerical examples were processed on a standard laptop computer under the Matlab environment with the SOS parser of YALMIP [START_REF] Löfberg | YALMIP : A Toolbox for Modeling and Optimization in MATLAB[END_REF], the moment parser GloptiPoly [START_REF] Henrion | GloptiPoly 3: moments, optimization and semidefinite programming[END_REF] and the semidefinite programming solver of MOSEK [START_REF] Dahl | Semidefinite optimization using MOSEK[END_REF]. For an interested reader, the codes used to obtain the results presented in this section are available online: https://homepages.laas. fr/henrion/software/stokesvolume/

Practical implementation

Following the Moment-SOS hierarchy methodology for volume computation as described in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF], in the (finite-dimensional) degree d semidefinite strengthening of dual problem ( 14) with unit euclidean ball as the bounding box B:

• w ∈ R[x] d and u ∈ R[x]
n d are polynomials of degree at most d; • the positivity constraint w ∈ C 0 (B) + is replaced with a Putinar certificate of positivity on B, that is:

w(x) = σ 0 (x) + σ 1 (x)(1 -|x| 2 ) , ∀x ∈ R n ,
where σ 0 (resp. σ 1 ) is an SOS polynomial of degree at most d (resp. d -2);

• the positivity constraint w| K -div u -1 ∈ C 0 (K) + is replaced with a Putinar certificate of positivity on K, that is:

w(x) -div u(x) -1 = ψ 0 (x) + ψ 1 (x) g(x) , ∀x ∈ R n ,
where ψ 0 (resp. ψ 1 ) is an SOS polynomial of degree at most d (resp. d -deg(g));

• the positivity constraint (u • grad g)| ∂K ∈ C 0 (∂K) + is replaced with a Putinar certificate of positivity on ∂K, that is:

-u(x) • grad g(x) = η 0 (x) + η 1 (x) g(x) , ∀x ∈ R n ,
where η 0 is an SOS polynomial of degree at most d and η 1 is a polynomial of degree at most d -deg(g);

• the linear criterion B w dλ translates into a linear criterion on the vector of coefficients of w, as B x α dλ is available in closed-form.

The above identities define linear constraints on the coefficients of all the unknown polynomials. Next, stating that some of these polynomials must be SOS translates into semidefinite constraints on their respective unknown Gram matrices. The resulting optimization problem is a semidefinite program, called the SOS strengthening of problem [START_REF] Lasserre | Convex optimization and parsimony of L p -balls representation[END_REF], and is in lagrangian duality with the so-called moment relaxation of problem [START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF]. Using the strong duality property, we interchangeably use the SOS strengthenings and moment relaxations, as they are equivalent; for more details the interested reader is referred to e.g. [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF].

Bivariate disk

Let us first illustrate Theorem 4.2 for computing the area of the disk The degree d = 16 polynomial approximation w obtained by solving the SOS strengthening of linear problem ( 5) is represented at the left of Figure 2. We can see bumps and ripples typical of a Gibbs phenomenon, since the polynomial should approximate from above the discontinuous indicator function 1 K as closely as possible. A rather loose upper bound of 1.1626 is obtained on the volume λ(K) = π 4 ≈ 0.7854. In comparison, the degree d = 16 polynomial approximation w obtained by solving the SOS strengthening of linear problem ( 14) is represented at the right of Figure 2. As expected from the proof of Theorem 4.2, the poynomial should approximate from above the continuous function g1 K λ(K)/( gλ K ). The resulting polynomial approximation is smoother and yields a much improved upper bound of 0.7870.

K := {x ∈ R 2 : g(x) = 1/4 -(x 1 -1/2) 2 -x 2 2 ≥ 0} included in the unit disk B := {x ∈ R 2 : 1 -x 2 1 -x 2 2 ≥ 0}.

Higher dimensions

In Table 1 we report on the dramatic acceleration brought by Stokes constraints in the case of the Euclidean ball

K := {x ∈ R 3 : g(x) = (3/4) 2 -|x| 2 ≥ 0} of dimension n = 3 included in the unit ball B.
We specify the relative errors on the bounds obtained by solving moment relaxations with and without Stokes constraints, together with the computational times (in seconds), for a relaxation degree d ranging from 4 to 20. We observe that tight bounds are obtained already at low degrees with Stokes constraints, sharply contrasting with the loose bounds obtained without Stokes constraints. However, we see also that the inclusion of Stokes constraints has a computational price.

In Table 2 we report the relative errors on the bounds obtained with and without Stokes constraints, together with the computational times (in seconds), for a relaxation degree equal to d = 10 (left) resp. d = 4 (right) and for dimension n ranging from 1 to 5 (left) resp. from

n
Table 2: Relative errors (%) and computational times (in brackets in seconds) for solving the degree d = 10 (left) and d = 4 (right) moment relaxation approximating the volume of a ball of increasing dimensions n.

Higher dimensional problems can be addressed only if the problem description has some sparsity structure, as explained in [START_REF] Tacchi | Exploiting sparsity in semi-algebraic set volume computation[END_REF]. Also, depending on the geometry of the problem, and for larger values of the relaxation degree, alternative polynomial bases may be preferable numerically than the monomial basis which is used by default in Moment and SOS parsers (see [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF]Fig.4.5]).

Changing the bounding box

Choosing the unit euclidean ball as our bounding box B is the easiest and most standard choice, but one could wonder what happens if we take another set, for example an p ball for p > 2.

Let B n p := {x ∈ R n : x p p = n i=1 |x i | p ≤ 1} denote the unit p ball in dimension n. We now compute the area of the bivariate disk K = {x ∈ R 2 : (3/4) 2 -x 2 1 -x 2 2 ≥ 0} included in the unit p ball B = B 2
p for p = 2, 4, 6, 8, 10. To that end, we use the closed formula for the Lebesgue moments on B 2 p (see Appendix A) :

B 2 p x k dx = 0 ∀k ∈ N n \ (2N) n and B 2 p x 2k dx = 2 (1 + |k|) p B 1 + 2k 1 p , 1 + 2k 2 p
where Γ(x) := ∞ 0 e -t t x-1 dt and B(x, y) := Γ(x)Γ(y) Γ(x + y) are Euler's Gamma and Beta functions, so that in particular one has Γ(1 + x) = x Γ(x).

We then perform the computations using Matlab's beta or gamma commands. Our numerical results are reported in Table 3. Unsurprisingly, as in the case of the euclidean bounding box, Stokes constraints drastically improve the accuracy of the moment relaxations. However, the number p has an influence on both accuracy (decreasing with p) and computational time (global tendency to slightly decrease with p) in the original as well as Stokes-augmented hierarchies. This can be explained by analyzing the influence of p on the SOS strengthenings described in section 6.1: indeed, the only change is that 1 -|x| 2 is replaced with 1 -x p p in the SOS representation of constraint w ∈ C(B) + , so that σ 1 now has degree at most d -p instead of d -2. Consequently, with increasing p, the size of the Gram matrix of σ 1 becomes smaller, slightly reducing the size of the corresponding SDP problem, which can result in a reduction of the computational time (although other factors impact the computational time, hence a non monotonic function of p). Conversely, as the degree of σ 1 is more limited, this brings less freedom in the search for an optimal solution, hence reducing the accuracy of the SOS strengthening for a fixed degree d.

d
In terms of the efficiency of Stokes constraints when the bounding box is an p ball, we highlight the fact that the loss in accuracy is less important in the Stokes-augmented hierarchy than in the standard formulation: Stokes constraints are somewhat more robust to the increase in the degree of the polynomial describing the bounding box. Regarding computational times, when they decrease with p, we observe that this decrease is more important with Stokes constraints than without: here again, Stokes constraints lead to a better behaved hierarchy.

Other sets

So far we only computed the volume of euclidean balls. In order to further explore the influence of the input polynomials on the efficiency of Stokes constraints as well as the Moment-SOS hierarchy in general, we now switch to computing the volume of more sophisticated semi-algebraic sets in two dimensions: first, we proceed from the euclidean ball to generic p balls, as we did with the bounding box; second, we test the limits of the scheme by approximating the volume of a non-convex, non-connected double disk.

4 disk

As discussed in section 6.4, the degree of the polynomials involved in the Moment-SOS hierarchy has a direct influence on its accuracy, as a higher input degree means, for a fixed degree of the hierarchy, less degrees of freedom to optimize over. We now show on one example what the practical influence of the degree over our scheme's accuracy is, by computing the approximate volume of the 4 disk:

K := x ∈ R 2 : (25/72) 4 -x 4 1 -x 4 2 ≥ 0 .
It is clear that one has λ(K) = (25/72) 2 λ(B 2 4 ) with, using formula (21) from Appendix A, λ(B 2 4 ) = Γ(1/4) 2

2

√ π so that λ(K) ≈ 0.4471. We implement the degree 16 SOS strengthenings corresponding to the standard and Stokes-augmented problems, and plot the resulting w in Figure 3. Again, the original SOS strengthening is flawed by a strong Gibbs phenomenon that introduces a large error in the volume approximation (we get a bound of 0.8511, i.e. a relative error of 90%), characterized by wide oscillations on the boundary of the 4 disk K. The Stokes-augmented version gives a tighter bound of 0.4653 (relative error 4%, still more than for the euclidean disk, but much less than without Stokes constraints). Moreover, an interesting feature appears here that was not visible on Figure 2 in the case of the euclidean disk: we observe small oscillations of w around 0 on B \ K. This can be expected as w is a non-zero polynomial, so it cannot vanish on a set of positive Lebesgue measure. These observations confirm our predictions that the lower the degree of the involved polynomials, the more accurate the SDP relaxations. However, we are now going to show that some other parameters should be considered when discussing the accuracy of the Moment-SOS hierarchy for volume computation, such as the geometry of the considered set K.

Disconnected double disk

We finally test our numerical scheme on a non-convex, non-connected semi-algebraic set: As usual, on Figure 4 we observe a Gibbs phenomenon in the standard volume approximation scheme (with a bound of 0.8551 instead of π 8 ≈ 0.3927, i.e. a relative error of 118%), as well as wide oscillations near the boundary of K. As for the Stokes-augmented scheme, again we get a better bound of 0.4671 (relative error 19%, interestingly higher than in all the previous cases with Stokes constraints, but still much more accurate than without Stokes constraints). More striking here, even in the Stokes-augmented SOS strengthening, w is seen clearly oscillating. However, this should not be mistaken for a consequence of the Gibbs phenomenon, as in this new formulation w is proved to approximate a Lipschitz continuous function. As a consequence to the Stone-Weierstrass theorem, those remaining oscillations are bound to ultimately vanish as the degree d goes to infinity, while in the case of the Gibbs phenomenon, the oscillations do not ultimately disappear (only their contribution to w dλ ultimately vanishes).

K := x ∈ R 2 : 1/16 -(x 1 -1/2) 2 -x 2 2 (x 1 + 1/2) 2 + x 2 2 -1/16 ≥ 0 .
A possible explanation for this oscillatory phenomenon is that, despite the regularity of the optimizer in the infinite dimensional problem ( 14), the SOS strengthenings are still very demanding for the polynomial w: indeed, it is requested to be as close to 0 as possible outside K while being sufficiently large in K so that its integral is bigger than λ(K). In the case of a non-connected K, w is thus literally requested to oscillate. To conclude on this example, we highlight the fact that, in addition to the degree of the polynomial g defining K, the geometry of K (typically: its number of connected components and how they are distributed in the bounding box B) plays a key role in the accuracy of the Moment-SOS hierarchy for computing its volume, both in the original and Stokes-augmented versions. Indeed, it is this geometry that is likely to generate (or, on the contrary, prevent) an oscillatory behavior in the approximating polynomial w, when one gets rid of the Gibbs phenomenon by complementing the hierarchy with Stokes constraints. This is particularly visible when comparing our examples in Sections 6.5.1 and 6.5.2, where K is described by degree 4 polynomials, but the schemes are far more accurate in the convex case than in the disconnected case, especially when one adds Stokes constraints.

Conclusion

In this paper we proposed a new primal-dual infinite-dimensional linear formulation of the problem of computing the volume of a smooth semi-algebraic set generated by a single polynomial, generalizing the approach of [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] while still allowing the application of the Moment-SOS hierarchy. The new dual formulation contains redundant linear constraints arising from Stokes's Theorem, generalizing the heuristic of [START_REF] Lasserre | Computing Gaussian and exponential measures of semi-algebraic sets[END_REF]. A striking property of this new formulation is that the dual value is attained, contrary to the original formulation. As a consequence, the corresponding dual SOS hierarchy does not suffer from the Gibbs phenomenon, thereby accelerating the convergence.

Numerical experiments (not reported here) reveal that the values obtained with the new Stokes constraints (with a general vector field) are closely matching the values obtained with the original Stokes constraints of [START_REF] Lasserre | Computing Gaussian and exponential measures of semi-algebraic sets[END_REF] (with the generating polynomial factoring the vector field). It may be then expected that the original and new Stokes constraints are equivalent. However at this stage we have not been able to prove equivalence.

The proof of dual attainment builds upon classical tools from linear PDE analysis, thereby building up a new bridge between infinite-dimensional convex duality and PDE theory, in the context of the Moment-SOS hierarchy. We expect that these ideas can be exploited to prove regularity properties of linear reformulations of other problems in data science, beyond volume approximation. For example, it would be desirable to design Stokes constraints tailored to the infinite-dimensional linear reformulation of the region of attraction problem [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF] or its sparse version [START_REF] Tacchi | Approximating regions of attraction of a sparse polynomial differential system[END_REF].

In terms of practical implementation, while still observed with Stokes constraints, the dependence on the degree of the input polynomials, already discussed in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF], seems to be of less importance. However, the dependence in the geometry of K now seems to prevail, as Stokes constraints add information on this geometry; more precisely, the simpler the geometry, the more efficient the constraints: the smooth and convex case leads to the best increase in accuracy, but the dual attainment still holds even on disconnected smooth sets. Also, experiments carried out in [START_REF] Lasserre | Computing Gaussian and exponential measures of semi-algebraic sets[END_REF][START_REF] Tacchi | Exploiting sparsity in semi-algebraic set volume computation[END_REF] show that even in the non-smooth case, Stokes constraints drastically improve the accuracy of the volume computing scheme.

A Moments of the Lebesgue measure on an p ball For our numerical experiments, we use closed formulae for the moments of the Lebesgue measure on an p ball. These formulae can be derived in a quite straightforward fashion using [14, Theorem 2.2].

Definition A.1 (Positively homogeneous functions). For d ∈ R, h : R n -→ R + is said to be positively homogeneous of degree d if for all x ∈ R n \ {0}, λ > 0, one has h(λx) = λ d h(x).

Lemma A.2 (Lebesgue moments and positively homogeneous functions, Theorem 2.2 in [START_REF] Lasserre | Convex optimization and parsimony of L p -balls representation[END_REF]). Let h : R n -→ R be a positively homogeneous function of degree d ∈ R \ {0}. Let B h := {x ∈ R n : h(x) ≤ 1} be such that λ(B h ) < ∞. Then, for all k ∈ N n ,

B h x k dx = Γ 1 + n + |k| d -1
R n

x k e -h(x) dx.

From this we deduce the following closed formula for moments of the Lebesgue measure on an p ball: Lemma A.3 (Lebesgue moments on the unit p ball). Let p > 1. The even moments of the Lebesgue measure on B n p are given, for k = (k 1 , . . . , k n ) ∈ (2N) n , by

B n p x k dx = 2 n p n Γ 1 + n + |k| p -1 n i=1 Γ 1 + k i p . ( 21
)
The odd moments are B n p

x k dx = 0 for k ∈ N n \ (2N) n .

Proof. Let k ∈ N n . Since x → x p p is positively homogeneous of degree p, we can use Lemma A.2:

B n p x k dx = Γ 1 + n + |k| p -1 R n
x k e -x p p dx. (a)

Then, Fubini's theorem ensures that R n

x k e -x p p dx = n i=1 +∞ -∞

x ki i e -|xi| p dx i .

Then, if one of the k i is odd, then the corresponding x i → x ki i e -|xi| p is an odd function, so that its integral over R is 0. Thus, this yields that B n p x k dx = 0. Finally, we compute for k ∈ 2N x k dx

= Γ 1 + 1 + k p 2 k + 1 = 2 k + 1 1 + k p Γ 1 + k p (c)
where we used the factorial property Γ(1 + x) = x Γ(x). ( 21) is then obtained by combining equations (a), (b) and (c).

Figure 1 :

 1 Figure 1: Gibbs effect occurring when approximating from above with a polynomial of degree 10 (left red curve) and 20 (right red curve) the indicator function of an interval (black curve).

Lemma 5 . 1 .

 51 Problem (14) has an optimal solution iff there exist

Figure 2 :

 2 Figure 2: Degree 16 polynomial approximations of the disk's area obtained without Stokes constraints (left) and with Stokes constraints (right).

Figure 3 :

 3 Figure 3: Degree 16 polynomial approximations of the area of the 4 disk obtained without Stokes constraints (left) and with Stokes constraints (right).

Figure 4 :

 4 Figure 4: Degree 16 polynomial approximations of the area of the double disk obtained without Stokes constraints (left) and with Stokes constraints (right).

Table 1 :

 1 Relative errors (%) and computational times (in brackets in seconds) for solving moment relaxations of increasing degrees d approximating the volume of ball of dimension n = 3.6 to 10 (right). When d = 10 and n = 5 the semidefinite relaxation features 6006 pseudomoments without Stokes constraints, and 12194 pseudo-moments with Stokes constraints. We see that introducing Stokes constraints incurs a computational cost, to be compromised with the expected quality of the bounds.

			d without Stokes with Stokes
		3	4	88% (0.03s)	18% (0.04s)
		3	8	57% (0.16s)	1.0% (0.44s)
		3 12	47% (1.97s)	0.0% (4.63s)
		3 16	43% (23.9s)	0.0% (30.1s)
		3 20	41% (142s)	0.0% (206s)
	n d without Stokes with Stokes	n d without Stokes	with Stokes
	1 10	17% (0.05s)	0.0% (0.03s)	6 4	190% (0.25s)	45.1% (1.03s)
	2 10	35% (0.09s)	0.2% (0.25s)	7 4	203% (0.32s)	60.0% (4.88s)
	3 10	56% (0.52s)	0.3% (1.19s)	8 4	221% (0.42s)	78.6% (8.45s)
	4 10	72% (9.74s)	0.4% (22.8s)	9 4	245% (1.15s)	102% (45.1s)
	5 10	79% (150s)	0.6% (669s)	10 4	278% (3.10s)	131%

Table 3 :

 3 Relative errors (%) and computational times (in brackets in seconds) for solving the degree d = 8 (left) and d = 16 (right) moment relaxations approximating the volume of ball of dimension n = 2 embedded in an p ball bounding box for p = 2, 4, 6, 8, 10.

	p without Stokes with Stokes	d	p without Stokes with Stokes
	8 2	39% (0.13s)	0.6% (0.14s)	16 2	27% (0.47s)	0.0% (0.70s)
	8 4	52% (0.13s)	1.2% (0.14s)	16 4	33% (0.41s)	0.0% (0.52s)
	8 6	57% (0.13s)	1.7% (0.13s)	16 6	36% (0.33s)	0.0% (0.50s)
	8 8	59% (0.13s)	2.0% (0.14s)	16 8	37% (0.30s)	0.0% (0.53s)
	8 10	61% (0.13s)	2.2% (0.14s)	16 10	38% (0.35s)	0.0% (0.49s)

The Gibbs phenomenon appears at a jump discontinuity when one numerically approximates a piecewise C 1 function with a polynomial function, e.g. by its Fourier series; see e.g.[START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF] Chapter 9].

Whereas K is a rather standard notation for a compact semi-algebraic set in polynomial optimization, the addition of the notation Ω is motivated by the conventions used for open sets in the differential geometry and PDE analysis literature. To account for both uses, we denote by Ω the interior of K.

An open set is said to be C 1 if its boundary is locally the graph of a C 1 function (up to reordering coordinates and changing orientation, see e.g.[START_REF] Evans | Partial Differential Equations[END_REF] Section C.1.]).

Acknowledgments

This work benefited from discussions with Swann Marx. The work of M. Tacchi was funded by the French company Réseau de Transport d'Électricité, as well as the Swiss National Science Foundation under the "NCCR Automation" grant n • 51NF40_180545.

The work of J.B. Lasserre was partly funded by the AI Interdisciplinary Institute ANITI through the French "Investing for the Future PIA3" program under the Grant agreement n • ANR-19-PI3A-0004.