Stokes, Gibbs and volume computation of semi-algebraic sets

We consider the problem of computing the Lebesgue volume of compact basic semi-algebraic 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 [1], or non deterministic Hit-and-Run methods [20,25] and the more recent [2,3].Even approximating λ(K) by deterministic methods is still a hard problem as explained in e.g.[3] 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 [9] provides a sequence (τ k ) k∈N ⊂ R of upper bounds that converges to λ(K) as k increases.Briefly, 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 (2) 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 [9].Also focusing on compact semi-algebraic sets, [12] 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 (1) with the sequence of problems (2) indexed by k.Indeed, • In [9], for any k ∈ N, τ k is guaranteed to be a converging upper bound for λ(K), i.e. τ k − λ(K) > 0, while [12] also guarantees convergence of the approximant to λ(K) but gives no information on the sign of the difference between the two quantities.
• [12] uses symbolic computations that can achieve arbitrary precision, while [9] uses numerical computations based on semidefinite programming, limited to floating-point arithmetic precision.
• The approach of [12] can be used to approximate other quantities than the volume, namely real periods of algebraic surfaces.The approach of [9] was extended to approximate sets relevant in systems control, such as regions of attraction or maximal positively invariant sets (see e.g.[6]).In this context, the present contribution can help improving the Moment-SOS hierarchy for assessing the stability of polynomial differential systems.
When solving problem (2), 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 [9] to accelerate this convergence but at the price of loosing monotonicity of the resulting sequence.
In fact ( 1) is a dual of the following infinite-dimensional Linear program (LP) on measures (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 (2).A first set of such linear constraints experimented in [15] and later in [16], 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 [15].Therefore the optimal continuous solution can be approximated efficiently by polynomials with no Gibbs phenomenon, by the hierarchy of semidefinite relaxations defined in [9] (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 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 and that its boundary 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 : 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 [9] a sequence of upper bounds converging to λ(K) is obtained by applying the Moment-SOS hierarchy [13] (a family of finite-dimensional convex relaxations) to approximate as closely as desired the (primal) infinite-dimensional LP on measures: whose optimal value is λ(K), attained for µ := λ K .The LP (4) has an infinite-dimensional LP dual on continuous functions which reads: 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 [13] for solving (4).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.[13], 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 [19] which translate to semidefinite constraints on the coefficients of polynomials; again the interested reader is referred to [9,13] 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 [9] by combining Stone-Weierstrass' theorem and Putinar's Positivstellensatz [19], (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.[9,Figs. 4.1 and 4.5].The best available theoretical convergence speed estimates are also pessimistic, with an asymptoptic rate of log log d [11].Slow convergence is mostly due to the so-called Gibbs phenomenon which is well-known in numerical analysis [23,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.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 [15] 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 [24].

Introducing Stokes constraints
In this section we explain the heuristic introduced in [15] 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
, 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 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 Proof.These are all particular cases of [10, Theorem 6.10.2].

Original Stokes constraints
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 (7) with u(x) := x k g(x): since by construction g vanishes on ∂K.Thus while in the infinite-dimensional LP (4) one may add the linear constraints without changing its optimal value λ(K), on the other hand inclusion of the linear moment constraints 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.[15,24].One main goal of this paper is to provide a qualitative mathematical rationale behind this phenomenon.

Infinite-dimensional Stokes constraints
In [21], 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: which can be written equivalently (in the sense of distributions) as This allows to rewrite problem (4) as without changing its optimal value λ(K) attained at µ = λ K .Using infinite-dimensional convex duality as in e.g. the proof of Theorem 2 in [6], the dual of LP (9) Crucial observation.Notice that w in (10) 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 (10) that can be exploited to yield more efficient semidefinite relaxations ?

New Stokes constraints and main result
In the previous section, the Stokes constraint (with µ ∈ M(K) + being the Lebesgue measure on K) was obtained as a particular case of Stokes' theorem with u = gv in (6).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 or equivalently (in the sense of distributions) 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 [16] to compute moments of the Hausdorff measure.In fact, equation ( 12) is a generalization of equation (11) in the following sense.(12), then µ also satisfies (11).
Next, by convex duality as in e.g. the proof of Theorem 2 in [6], the dual of (13) reads Our main result states that the optimal value of the dual ( 14) is attained at some continuous function Therefore, in contrast with problem (5), there is no Gibbs phenomenon at an optimal solution of the (finite-dimensional) semidefinite strengthening associated with (14).
Let Ω i , i = 1, . . ., N denote the connected components of Ω, and let Theorem 4.2.In dual LP (14) the infimum is a minimum, attained at where u solves the Poisson PDE 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 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 (14).
From Lemma 5.1, existence of an optimum for ( 14) is then equivalent to existence of a solution to (15), 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 then problem (14) 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 (17) is only a sufficient condition for existence of an optimum for (14), since (15) 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 (14) 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 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 (17).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: Indeed, if (f, u) solves ( 17), then 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 [18]: 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 and [18] yields a solution 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 [18].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 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 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

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 [17], the moment parser GloptiPoly [8] and the semidefinite programming solver of MOSEK [4].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 [9], in the (finite-dimensional) degree d semidefinite strengthening of dual problem ( 14) with unit euclidean ball as the bounding box B: 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: 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: 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: 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 (14), and is in lagrangian duality with the so-called moment relaxation of problem (13).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.[9].

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.

Higher dimensions
In Table 1 we report on the dramatic acceleration brought by Stokes constraints in the case of the Euclidean ball 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 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 [21].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 [9,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.
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.
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: 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).
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 [9] 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 [15].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 [15] (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 [6] or its sparse version [22].
In terms of practical implementation, while still observed with Stokes constraints, the dependence on the degree of the input polynomials, already discussed in [9], 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 [15,21] 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).

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

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

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

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