Pierre Roux

Yuen-Lam Voronin

Sriram Sankaranarayanan

Validating Numerical Semidefinite Programming Solvers for Polynomial Invariants

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Program analysis techniques using abstract interpretation, especially numerical domain program analysis, rely fundamentally on the ability to reason about constraints expressed in a suitable logic that stems from the abstract domain. Typical reasoning tasks include the problem of checking satisfiability of an assertion in the logic used in emptiness and inclusion checks, and characterizing elements of the cone of consequences of an assertion used to compute the transfer function and join operations [START_REF] Cousot | Abstract Interpretation: A Unified Lattice Model for Static Analysis of Programs by Construction or Approximation of Fixpoints[END_REF]. The process of using basic solver primitives has led to many constraint-based approaches to synthesizing and verifying invariants for programs [START_REF] Adjé | Coupling Policy Iteration with Semidefinite Relaxation to Compute Accurate Numerical Invariants in Static Analysis[END_REF][START_REF] Cousot | Proving Program Invariance and Termination by Parametric Abstraction, Lagrangian Relaxation and Semidefinite Programming[END_REF][START_REF] Gaubert | Static analysis by policy iteration on relational domains[END_REF][START_REF] Gawlitza | Precise fixpoint computation through strategy iteration[END_REF][START_REF] Martin | Improving strategies via SMT solving[END_REF][START_REF] Sankaranarayanan | Scalable analysis of linear systems using mathematical programming[END_REF]. Initial approaches that focused on linear systems [START_REF] Gaubert | Static analysis by policy iteration on relational domains[END_REF][START_REF] Gawlitza | Precise fixpoint computation through strategy iteration[END_REF][START_REF] Sankaranarayanan | Scalable analysis of linear systems using mathematical programming[END_REF] have been generalized to address nonlinear (polynomial) systems [START_REF] Adjé | Coupling Policy Iteration with Semidefinite Relaxation to Compute Accurate Numerical Invariants in Static Analysis[END_REF][START_REF] Xavier Allamigeon | A scalable algebraic method to infer quadratic invariants of switched systems[END_REF][START_REF] Martin | Computing relaxed abstract semantics w.r.t. quadratic zones precisely[END_REF][START_REF] Oulamara | Abstract interpretation with higher-dimensional ellipsoids and conic extrapolation[END_REF]. Other extensions to hybrid systems and stochastic systems have also been proposed [START_REF] Chakarov | Deductive Proofs of Almost Sure Persistence and Recurrence Properties[END_REF][START_REF] Dang | Template-Based Unbounded Time Verification of Affine Hybrid Automata[END_REF][START_REF] Prajna | Safety verification using barrier certificates[END_REF].

However, extensions to polynomial systems necessarily face the challenge of reasoning about polynomial inequality constraints. While the problem of checking satisfiability of these constraints is well-studied, precise solutions to this problem are as yet intractable for large problems. Likewise, computing the cone of consequences precisely is also prohibitively expensive in practice, requiring quantifier elimination. A more tractable alternative uses a convex relaxation from the given polynomial system to a semidefinite programming problem (SDP) using the sum of squares (SOS) relaxation [START_REF] Bernard | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF][START_REF] Naum | Class of global minimum bounds on polynomial functions[END_REF]. Such a relaxation guarantees soundness when used in invariant checking/synthesis tasks, and has been shown to have nice theoretical guarantees. However, in practice, the approach requires us to use SDP solvers. It is well-known that precise solutions of SDPs is a hard and open research question. For instance, there are SDPs which have a feasible solution but no rational feasible solutions. Therefore, most solvers seek an approximate solution. At the same time, it is well known (but not well documented) that numerical SDP solvers are also hard to use in practice. The presence of "numerical issues" leads to unreliable answers from the SDP solvers, that in turn lead to unsound results when employed in program analysis tasks.

In this paper, we characterize numerical issues into many types. At one end of the spectrum, we have issues that arise from floating-point errors and approximate answers, since numerical solvers seldom reach a true optimal solution. At the other end, certain problems are not well posed, depending on the nature of the solution technique used. One common reason involves the failure of strict feasibility. Using actual examples from the literature, we show how the answers from popular numerical SDP solvers can be wrong and potentially mislead even a careful user who pays due attention to the various errors reported by the SDP solver.

Finally, we address some of the numerical problems raised. We first present a sound verification procedure that can check the answer from the solver and help us decide whether the answer is qualitatively correct. Next, we provide a padding procedure that helps reformulate a given problem into a stricter version so that if an approximate, floating-point solver can find a reliable answer to the stricter version, then we conclude feasibility of the original version. We integrate our framework into a polynomial invariant synthesis/verification task, showing how our ideas can successfully address numerical issues arising from the solver.

Motivating Examples

In this section, we illustrate through two examples the scenario where numerical SDPs give seemingly sensible solutions to simple invariant generation problems, and yet the generated invariants are not sound.

Consider the program in Fig. 1. Does there exist an inductive invariant 3 in the form (x 1 , x 2 ) ∈ R 2 p(x 1 , x 2 ) ≥ 0 for some polynomial p? A tractable sufficient condition that guarantees this can be formulated using the SOS optimization approach (see Section 3), resulting in an SDP instance that can be solved by numerical solvers. The widely used SDPT3 [START_REF] Reha H Tütüncü | Solving semidefinite-quadratic-linear programs using SDPT3[END_REF] solver reports a solution. Although all the DIMACS errors [START_REF] Schmieta | Reporting solution quality for the DIMACS library of mixed semidefinite-quadratic-linear programs[END_REF] are less than 10 -8 , not raising any suspicion, we found traces of the program that violate this purported invariant (see Fig. 1).

As another example, we consider a program from ADJ É et al. [START_REF] Adjé | Property-based Polynomial Invariant Generation Using Sums-of-Squares Optimization[END_REF] and the "invariant" they offer, generated with numerical solvers (Fig. 2). Note that the purported invariant is

(x1, x2) ∈ x1, x2 x 2 1 + x 2 2 ≤ 1.5 2 while (1) { // Find Inv. p(x1, x2) ≥ 0 x1 = x1 * x2; x2 = -x1; } p(x1, x2) : 1 + 2.46x 2 1 + 2.46x 2 2 -5 × 10 -7 x 4 1 -2.46x 2 1 x 2 2 -5 × 10 -7 x 4 2 -2 0 2 -2 0 2
Fig. 1: (Left) An example program, and "loop invariant" p(x 1 , x 2 ) ≥ 0 synthesized using numerical solvers. (Right) The claimed "invariant" and dashed lines showing violations. (Right) Purported invariant at loop head synthesized using SDP solvers [START_REF] Adjé | Property-based Polynomial Invariant Generation Using Sums-of-Squares Optimization[END_REF].

(x1, x2) ∈ [0.9, 1.1] × [0, 0.2] while (1) { pre_x1 = x1; pre_x2 = x2; if (x1ˆ2 + x2ˆ2 <= 1) { x1 = pre_x1ˆ2 + pre_x2ˆ3; x2 = pre_x1ˆ3 + pre_x2ˆ2; } else { x1 = 0.5 * pre_x1ˆ3 + 0.4 * pre_x2ˆ2; x2 = -0.6 * pre_x1ˆ2 + 0.3 * pre_x2ˆ2; } } 2.510902467 + 0.0050x1 + 0.0148x2 -3.0998x 2 1 + 0.8037x 3 2 + 3.0297x 3 1 -2.5924x 2 2 -1.5266x1x2 + 1.9133x 2 1 x2 + 1.8122x1x 2 2 -1.6042x 4 1 -0.0512x 3 1 x2 + 4.4430x 2 1 x 2 2 + 1.8926x1x 3 2 -0.5464x 4 2 + 0.2084x 5 1 -0.5866x 4 1 x2 -2.2410x 3 1 x 2 2 -1.5714x 2 1 x 3 2 + 0.0890x1x 4 2 + 0.9656x 5 2 -0.0098x 6 1 + 0.0320x 5 1 x2 + 0.0232x 4 1 x 2 2 -0.2660x 3 1 x 3 2 -0.7746x 2 1 x 4 2 -0.9200x1x 5 2 -0.6411x 6 2 ≥ 0
-1 0 1 2 3 4 -1 0 1 -1 0 1 -1 0 1 Fig. 3: (Left)
The candidate invariant from Fig. 2 with arrows showing concrete transitions. The arrows leaving it are counterexamples to its inductiveness. (Right) The invariant of degree 8 whose soundness is proved using the approach in this paper.

indeed not inductive: one can find points in it whose image after one iteration of the loop body exits the invariant (Fig. 3). Fig. 3 also depicts an actual invariant, proved using the method in this paper.

3 Sum of Squares (SOS) and Semidefinite Programs (SDP)

In this section, we provide the background for sum of squares (SOS) optimization from the perspective of proving entailments for program analysis (invariant synthesis/verification) tasks. We then trace the steps along which SOS optimization problems are relaxed to semidefinite programs (SDP). Let R denote the set of real numbers and x : (x 1 , . . . , x n ) denote a vector of realvalued variables. The ring of multivariate real polynomials over x is denoted by R[x]. The degree of any polynomial p(x) ∈ R[x] is denoted by deg(p). A template polynomial is of the form p(c, x) : s j=1 c j p j (x), where p 1 , . . . , p s are basis polynomials and c : (c 1 , . . . , c s ) is a placeholder for parameters serving as scalar multiples of the basis polynomials. A generic template polynomial of degree d > 0 is formed by choosing all monomials of degree up to d as the basis polynomials, and has s = n+d d parameters.

Semi-algebraic Assertions and Entailment Problems

A semi-algebraic assertion ϕ is a finite conjunction of polynomial inequalities:

ϕ : p 1 (x) ≥ 0 ∧ • • • ∧ p m (x) ≥ 0.
It denotes a corresponding semi-algebraic set ϕ : {x ∈ R n | x |= ϕ}. As such, semialgebraic assertions subsume useful abstract domains such as polyhedra and ellipsoids. They also represent a rich class of constraints with a decidable entailment checking problem [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]. We define two classes of problems involving semi-algebraic sets that are commonly used as primitives.

Definition 1 (Entailment Checking). Given two semi-algebraic assertions ϕ and ψ over x, check if ϕ |= ψ, i.e., for all x ∈ R n , if x |= ϕ, then x |= ψ. Definition 2 (Parametric Entailment). Let c : (c 1 , . . . , c s ) represent parameters. The input to a parametric entailment problem consists of k pairs (ϕ i , p i ) k i=1 , wherein ϕ i is a semi-algebraic assertion and p i (c, x) is a template polynomial. The goal is to compute a value c such that all the entailments hold:

(ϕ 1 |= p 1 (c, x) ≥ 0) ∧ • • • ∧ (ϕ k |= p k (c, x) ≥ 0) .
The entailment checking and its analog of parametric entailment checking are fundamental primitives that we will use for synthesizing and checking invariants of programs. The example below illustrates the application of these primitives. Invariant checking: Consider the program in Fig. 4. We wish to prove that all executions remain inside a safe set S :

(x1, x2) ∈ I : (x1, x2) x 2 1 ≤ 1 ∧ x 2 2 ≤
{(x 1 , x 2 ) | |x 1 | ≤ 2 ∧ |x 2 | ≤ 2}. To prove this, we consider an inductive invariant {(x 1 , x 2 ) | p(x 1 , x 2 ) ≥ 0}, where p(x 1 , x 2 ) : 37 -x 2 2 + x 3 1 -2x 2 1 x 2 + 2x 3 2 -12x 4 1 -10x 2 1 x 2 2 -6x 1 x 3 2 -6x 4 2 . (1)
Fig. 4 shows the invariant region p(x 1 , x 2 ) ≥ 0. To show that it is indeed an invariant that establishes S as a safe set, we check that the following conditions hold:

(a) Initial condition: 1 -x 2 1 ≥ 0 ∧ 1 -x 2 2 ≥ 0 |= p(x 1 , x 2 ) ≥ 0, (b) Consecution (loop) conditions: Let τ 1 (x 1 , x 2 ) : (0.687x 1 + 0.558x 2 -0.0001x 1 x 2 ,
-0.292x 1 + 0.773x 2 ) denote the transition enabled by the condition x 1 ≥ x 2 , and τ 2 (x 1 , x 2 ) : (0.369x 1 + 0.532x 2 -0.0001x 2 1 , -1.27x 1 + 0.12x 2 -0.0001x 1 x 2 ) denote the transition enabled by the condition x 1 < x 2 . We require two conditions, corresponding to the two transitions in the loop:

(i)

x 1 -x 2 ≥ 0 ∧ p(x 1 , x 2 ) ≥ 0 |= p1 ≥ 0, where p1 = p • τ 1 , and (ii) x 2 -x 1 ≥ 0 ∧ p(x 1 , x 2 ) ≥ 0 |= p2 ≥ 0, where p2 = p • τ 2 . (c) Safety conditions: p(x 1 , x 2 ) ≥ 0 |= -2 ≤ x 1 ≤ 2 ∧ -2 ≤ x 2 ≤ 2.
The invariant checking problem is then a series of polynomial entailment checking.

Invariant synthesis:

The invariant synthesis problem requires us to synthesize polynomials that satisfy some entailment conditions, such as p(x 1 , x 2 ) in (1) satisfying the initial, loop and safety conditions. To do so, we parameterize a template polynomial as follows:

p(c, x) : c 1 + c 2 x 1 + c 3 x 2 + c 4 x 2 1 + c 5 x 1 x 2 + c 6 x 2 2 + c 7 x 3 1 + c 8 x 1 x 2 2 + c 9 x 2 1 x 2 + c 10 x 3 2 + c 11 x 4 1 + c 12 x 3 1 x 2 + c 13 x 2 1 x 2 2 + c 14 x 1 x 3 2 + c 15 x 4 2 . (2)
We then search for values of c = (c 1 , . . . , c 15 ) such that the following entailments hold:

(a) Initial condition:

1 -x 2 1 ≥ 0 ∧ 1 -x 2 2 ≥ 0 |= p(c, x) ≥ 0, (b) Loop conditions: The loop condition for the transition τ 1 is (x 1 -x 2 ≥ 0 ∧ p(c, x) ≥ 0 |= p1 (c, x) ≥ 0)
, where p1 = p • τ 1 . However, we note that in this condition, a parametric polynomial inequality appears on the antecedent side. This leads to hard bilinear optimization problems that are beyond the scope of this paper (see [START_REF] Sriram Sankaranarayanan | Constructing invariants for hybrid systems[END_REF] for a discussion of this issue). We impose a stronger condition on p that states that p must be non-decreasing for each loop iteration4 :

(i) x 1 -x 2 ≥ 0 |= p1 (c, x) ≥ p(c, x), and (ii) x 2 -x 1 ≥ 0 |= p2 (c, x) ≥ p(c, x), where p2 = p • τ 2 .
The invariant synthesis problem is thus reduced to a parametric entailment problem.

Solving Entailment Problems

There are numerous approaches to solving semi-algebraic entailment checking and parametric entailment problems. We classify these into four broad classes: (a) quantifier elimination for the theory of polynomial inequalities, (b) interval arithmetic with branchand-bound, (c) linear programming relaxations, and (d) sum of squares relaxations that will be the focus of our exposition.

Exact Approaches: It is well-known that the logical theory of polynomial inequalities admits effective decision procedures and a quantifier elimination procedure, originally discovered by TARSKI and further developed by COLLINS, HONG, WEISPFENNING and others [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF][START_REF] Collins | Quantifier elimination for real closed fields by cylindrical algebraic decomposition[END_REF][START_REF] Collins | Partial Cylindrical Algebraic Decomposition for Quantifier Elimination[END_REF][START_REF] Tarski | A decision method for elementary algebra and geometry[END_REF][START_REF] Volker | Quantifier elimination for real algebra-the quadratic case and beyond[END_REF]. These procedures attempt to solve the entailment problem ϕ |= ψ by checking the unsatisfiability of the assertion ϕ(x) ∧ (¬ψ(x)). This problem is known to be NP-hard in theory, and hard to solve, in practice. Typical sizes of problems that can be tackled involve polynomials with ∼ 5 variables, and degrees ∼ 3 [START_REF] Dolzmann | REDLOG: Computer algebra meets computer logic[END_REF].

Likewise, an exact approach for the parametric entailment problem requires to perform a quantifier elimination of the form:

(∀ x) ϕ 1 (x) ⇒ p 1 (c, x) ≥ 0 ∧ • • • ∧ ϕ k (x) ⇒ p k (c, x) ≥ 0 .
Doing so leads to an assertion that can be expressed purely in terms of c. If this assertion is satisfiable, a solution c = c * can be extracted.

Branch-and-Bound (BnB) Approaches: They work over states x that are a priori restricted to a compact set X. They proceed by subdividing X into finitely many interval (hyper-rectangular) cells. Inside each interval, the entailment is evaluated using interval arithmetic [START_REF] Fränzle | Efficient solving of large non-linear arithmetic constraint systems with complex Boolean structure[END_REF][START_REF] Gao | dreal: An SMT solver for nonlinear theories over the reals[END_REF] or a branch-and-bound scheme using linear programming relaxation of the constraints [START_REF] Amin | Linear Relaxations of Polynomial Positivity for Polynomial Lyapunov Function Synthesis[END_REF]. BnB approaches can be used to check whether an entailment ϕ |= ψ holds by checking the unsatisfiability of the assertion ϕ(x) ∧ (¬ψ(x)). They can conclude soundly that the entailment holds or even find a witness x such that ϕ(x) ∧ (¬ψ(x)) is satisfied. Unfortunately, due to computational limitations, these techniques may also terminate without an answer. Recent work on delta-satisfiability procedures have carefully analyzed this condition to conclude that a "nearby" formula is satisfiable [START_REF] Gao | dreal: An SMT solver for nonlinear theories over the reals[END_REF]. BnB approaches can extend beyond polynomial programs and invariants. Currently, BnB approaches are restricted to solving entailment problems. Their application to solving parametric entailment problems remains an open challenge. Part of the challenge involves the optimal subdivision of X to search for solutions c.

Linear Programming (LP) Relaxations: Linear programming relaxations have been considered for checking polynomial entailment and solving parametric entailment problems (see BEN SASSI et al. for details and further references [START_REF] Amin | Linear Relaxations of Polynomial Positivity for Polynomial Lyapunov Function Synthesis[END_REF]). LP approaches are primarily based on so-called Handelman relaxation and reformulation linearization. Given an entailment problem with p 1 , . . . , p k as antecedents. We generate "valid inequalities" that are consequences of the original antecedents. This is achieved by simply multiplying the antecedents together, enriching the set of possible antecedents. This step is inspired by the Handelman positivstellensatz [START_REF] Handelman | Representing polynomials by positive linear functions on compact convex polyhedra[END_REF]. Next, we introduce fresh variables corresponding to each monomial term and turn our polynomial entailment problem into a linear entailment problem that can be checked using solvers. This step is called reformulation linearization technique (RLT) [START_REF] Hanif | A global optimization algorithm for polynomial programming using a reformulation-linearization technique[END_REF]. BEN SASSI et al. show that the generation of linear constraints can be performed in the Bernstein polynomial basis [START_REF] Nanatovich | Démonstration du théoréme de Weierstrass fondée sur le calcul des probabilités[END_REF][START_REF] Rida | The Bernstein polynomial basis: A centennial retrospective[END_REF], rather than the monomial basis to obtain a larger set of valid inequalities. The LP approach has the main advantage that Simplex solvers can be used with exact arithmetic to completely avoid numerical issues. The recent work of MAR ÉCHAL et al. use this approach and generate machine checkable proofs of polynomial entailments [START_REF] Maréchal | Polyhedral approximation of multivariate polynomials using handelman's theorem[END_REF]. However, LP relaxations yield a "weak" proof system that requires higher degree terms or a BnB decomposition of the domain, to prove "simple consequences" such as [START_REF] Amin | Linear Relaxations of Polynomial Positivity for Polynomial Lyapunov Function Synthesis[END_REF]. Interestingly, the Handelman relaxation and RLT are implicit in the polyhedral abstract domain for computing semialgebraic invariants proposed by BAGNARA et al. [START_REF] Bagnara | Generation of Basic Semi-algebraic Invariants Using Convex Polyhedra[END_REF]. A related approach of diagonal SOS (DSOS) has been proposed by ALI AHMADI and MAJUMDAR [START_REF] Amir | DSOS and SDSOS optimization: LP and SOCPbased alternatives to sum of squares optimization[END_REF]. Their approach is based on the SOS relaxation wherein instead of reducing to a SDP, they reduce to an LP by imposing the stronger condition of diagonal dominance on the associated matrix rather than the semi-definiteness condition that will be described subsequently in this section. A full comparison of DSOS with SOS relaxations for program analysis problems is currently open.

-1 ≤ x ≤ 1 ∧ -1 ≤ y ≤ 1 |= (x 2 + y 2 ≥ 0)
Positivstellensatz/Sum of Squares (SOS) Relaxations: The SOS relaxation [START_REF] Bernard | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF] is an incomplete but efficient way to numerically solve polynomial entailment problems.

Definition 3 (SOS Polynomial

). A polynomial p ∈ R[x] is said to be SOS if there exist polynomials h i ∈ R[x] such that for all x, p(x) = i h 2 i (x).
Although not all nonnegative polynomials are SOS, being SOS is a sufficient condition to be nonnegative.

Example 1. Consider p : 2x 4 1 + 2x 3 1 x 2 -x 2 1 x 2 2 + 5x 4 2 . Since p = h 2 1 + h 2 2
, where h 1 :

1 √ 2 2x 2 1 + x 1 x 2 -3x 2 2 and h 2 : 1 √ 2 3x 1 x 2 + x 2 2 , the polynomial p is nonnegative, i.e., p(x 1 , x 2 ) ≥ 0 holds for all x 1 , x 2 ∈ R 2 .
Now consider a polynomial entailment problem of the form:

p 1 (x) ≥ 0 ∧ • • • ∧ p k (x) ≥ 0 ϕ |= p(x) ≥ 0. ( 3
)
Our goal is to write p as a combination of p 1 , . . . , p k in the following form:

p = σ 0 + σ 1 p 1 + • • • + σ k p k (4)
such that σ 0 , . . . , σ k are SOS polynomials over x. Let K denote the semi-algebraic set

k j=1 p j ≥ 0 , let R(K) denote the cone of consequences of K, i.e., R(K) = {p(x) | ϕ |= p(x)
≥ 0} and M (K) denote all polynomials expressible in the form (4):

M (K) = p(x) p = σ 0 + k j=1 σ j p j , σ i SOS . Theorem 1 (Putinar's Positivstellensatz). For all K, M (K) ⊆ R(K). Conversely, if K is compact and M (K) contains a polynomial of the form s(x) = n i=1 x 2 i -L for some constant L > 0, then M (K) = R(K).
Proof. We will prove the "easy" direction that

M (K) ⊆ R(K). Let p ∈ M (K). There exist SOS polynomials σ 0 , . . . , σ k such that p = σ 0 + k i=1 σ i p i . Let x be such that p i (x) ≥ 0 for all i ∈ [1, k]. We have that σ i (x) ≥ 0 since each σ i is a SOS polynomial. Therefore, we conclude that p(x) = σ 0 (x) + k i=1 σ i (x)p i (x) ≥ 0.
For the converse, we refer the reader to Putinar's work [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF].

Thus, a polynomial entailment problem of the form (3) is relaxed to an SOS problem:

find : polynomials σ 0 , . . . , σ k ∈ R d [x] s.t. p = σ 0 + k i=1 σ i p i , σ 0 , . . . , σ k are SOS. (5)
First, we choose a degree limit d > 0 (d must be an even number because all positive polynomials have even maximum degree), and select templates σ 0 (c (0) , x), . . . , σ k (c (k) , x) with unknowns c (0) , . . . , c (k) . We then require that p be equal to a polynomial combination of p 1 , . . . , p k with "multipliers" σ 0 , . . . , σ k as shown above. This yields a set of linear equations involving c (0) , . . . , c (k) and the coefficients of p, obtained by comparing both sides monomial by monomial and setting their coefficients to be the same. Finally, we require σ 0 , . . . , σ k to be SOS. This will be tackled through a reduction to a semidefinite programming (SDP) problem, as will be explained subsequently.

Example 2. Consider the initial condition check for the program in Fig. 4:

p 1 (x 1 , x 2 ) ≥ 0 ∧ p 2 (x 1 , x 2 ) ≥ 0 |= p(x 1 , x 2 ) ≥ 0 with p i (x 1 , x 2 ) = 1 -x 2
i and p given in (1). Our goal here is to find polynomials σ 0 , σ 1 , σ 2 such that p = σ 0 + σ 1 p 1 + σ 2 p 2 . For simplicity, let us write

σ 0 = c 1 + c 2 x 1 + • • • + c 15 x 4 2 , σ 1 = d 1 + • • • + d 15 x 4 2 and σ 2 = e 1 + • • • + e 15 x 4 2
. We obtain equality constraints by equating terms corresponding to the same monomial on both sides:

c 1 + d 1 + e 1 =
37 (comparing constant terms), . . . , -e 15 = 0 (comparing x 6 2 ).

The SOS problem seeks to satisfy these equalities, and additionally make σ 0 , σ 1 , σ 2 SOS. Solving this as an SDP problem (as will be explained below), we obtain: σ1 ≈ 11 -

0.13x1+1.5x2+24x 2 1 -3x1x2+8.2x 2 2 and σ2 ≈ 8.8+0.63x1-1.4x2+6.5x 2 1 +1.6x1x2+18x 2 2 .
SOS formulation of parametric entailment problems. Consider now a parametric entailment problem of the form (ϕ j |= p j (c, x) ≥ 0) for j = 1, . . . , K involving parameters c. Let us write ϕ j :

p j1 (x) ≥ 0 ∧ • • • ∧ p j l (x) ≥ 0.
This is reduced to a sum of squares problem:

find : polynomials σ j,0 , . . . , σ j,j l ∈ R d [x], j ∈ {1, . . . , K} s.t. p j = σ j,0 + k i=1 σ j,i p ji , j ∈ {1, . . . , K} , σ j,0 , . . . , σ j,j l are SOS, j ∈ {1, . . . , K}.

(6)

The unknowns include the coefficients c involved in each p j (c, x) for the original parametric entailment and the coefficients c (j,i) corresponding to SOS multipliers σ j,i .

Next, we provide a reduction from SOS problems to a well known class of optimization problems: semidefinite programs (SDPs). Any polynomial p of degree 2d (a nonnegative polynomial is necessarily of even degree) can be written as a quadratic form in the vector z of all monomials of degree less or equal d:

p(x) = z T Q z, (7)
where z = 1, x 1 , . . . , x n , x 1 x 2 , . . . , x d n T and Q is a constant symmetric matrix.

Example 3. Consider p(x 1 , x 2 ) : 2x 4 1 + 2x 3 1 x 2 -x 2 1 x 2 2 + 5x 4 2 .
To satisfy the equality

p(x 1 , x 2 ) =   x 2 1 x 2 2 x 1 x 2   T 
 q 11 q 12 q 13 q 12 q 22 q 23 q 13 q 23 q 33

    x 2 1 x 2 2 x 1 x 2   = q 11 x 4 1 + 2q 13 x 3 1 x 2 + (q 33 + 2q 12 )x 2 1 x 2 2 + 2q 23 x 1 x 3 2 + q 22 x 4 2 ,
the equalities q 11 = 2, 2q 13 = 2, q 33 + 2q 12 = -1, 2q 23 = 0 and q 22 = 5 must hold.

Two possible examples for the matrix Q are shown below:

Q =   2 1 1 1 5 0 1 0 -3   , Q =   2 -3 1 -3 5 0 1 0 5   .
The polynomial p is then SOS if and only if there exists a positive semidefinite matrix Q satisfying [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]. A matrix Q is said to be positive semidefinite, denoted by Q 0, when for all vectors y, y T Q y ≥ 0. A matrix Q is said to be positive definite, denoted by Q 0, when for all nonzero vectors y, y T Q y > 0.

Example 4. In Example 3, the first matrix Q is not positive semidefinite (for y : [0, 0, 1]

T , y T Q y = -3). However, the second matrix Q is positive semidefinite as it can be written

Q = L T L with L = 1 √ 2 2 -3 1 0 1 3
(then, for all y, y T Qy = (Ly) T (Ly) = Ly 2 2 ≥ 0). This gives the SOS decomposition of Example 1:

p(x 1 , x 2 ) = 1 2 (2x 2 1 + x 1 x 2 -3x 2 2 ) 2 + 1 2 (3x 1 x 2 + x 2 2 ) 2 .
As a result, SOS programming problems can be written as semidefinite optimization problems involving matrices. Let z be a vector of monomials over x chosen so that we may write each polynomial σ i (c (i) , x) as a quadratic form σ i (c (i) , x) = z T C i z. Thus, the SOS programming problems ( 5) and ( 6) can be written down as an SDP problem:

find : c, C 0 , . . . , C k s.t. a T i c + k j=0 tr(A i,j C j ) = b i , i = 1, . . . , m, C j 0, j = 0, . . . , k. (8)
wherein the vector c encodes the parameters c (i) of ( 5) (or c and c (j,i) of ( 6)) and the A i,j and C j are symmetric matrices. Note that the expression tr(XY ) equals n i=1 n j=1 (X) i,j (Y ) i,j for n × n matrices X, Y when X T = X. Example 5. We check whether the entailment (p 1 (x, y) : x -y ≥ 0 |= p(x, y) :

x -y + 2x 2 -2y 2 + x 3 + x 2 y -xy 2 -y 3 ≥ 0) is true using an SOS relaxation: we look for degree 2 SOS polynomials σ 0 , σ 1 ∈ R[x, y] such that p = σ 0 + σ 1 p 1 . In other words, we seek coefficients c (0) : (c 1 , . . . , c 6 ) and c (1) : (c 7 , . . . , c 12 ) such that σ 0 (x, y) : c 1 +c 2 x+c 3 y+c 4 x 2 +c 5 xy+c 6 y 2 and σ 1 (x, y) : c 7 +c 8 x+c 9 y+c 10 x 2 +c 11 xy+c 12 y 2 are SOS and the coefficients of p and σ 0 + σ 1 p 1 coincide, i.e., with z : [1, x, y] T , comparing coeffs of σ0 + σ1 p1 and p :

                       const. term c1 = 0, coeff. of x c2 + c7 = 1, coeff. of y c3 -c7 = -1, coeff. of x 2 c4 + c8 = 2, coeff. of xy c5 -c8 + c9 = 0, coeff. of y 2 c6 -c9 = -2, coeff. of x 3 c10 = 1, coeff. of x 2 y c11 -c10 = 1, coeff. of xy 2 c12 -c11 = -1, coeff. of y 3 -c12 = -1,
comparing coefficients of σ0 and z T C0 z same for σ1 and z T C1 z.

:

         c1 = (C0)1,1, c2 = (C0)1,2 + (C0)2,1, c3 = (C0)1,3 + (C0)3,1, c4 = (C0)2,2, c5 = (C0)2,3 + (C0)3,2, c6 = (C0)3,3, (9)
Each of the m = 22 equality constraints in ( 9) is then encoded as in [START_REF] Amin | Linear Relaxations of Polynomial Positivity for Polynomial Lyapunov Function Synthesis[END_REF]. For instance, the second constraint on σ 0 is encoded by the vector a 12 = [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] T and the matrices A 12,0 = 0 -1 0 -1 0 0 0 0 0 and A 12,1 = 0.

Thus, we have eliminated the formal variables x from the problem and reduced it to finding matrices that satisfy some linear equality constraints, and are positive semidefinite. In fact, moving one step further, we write a single unknown matrix C in the block diagonal form:

C = Diag(c + 1 , c - 1 , . . . , c + s , c - s , C 0 , C 1 , . . . , C k ), encoding c i as c + i -c - i with c + i , c - i ∈ R + .
This allows us to write (8) as:

find : C s.t. tr(A i C) = b i , i = 1, . . . , m, C 0. ( 10
)
Problems that follow this form, or equivalently [START_REF] Amin | Linear Relaxations of Polynomial Positivity for Polynomial Lyapunov Function Synthesis[END_REF], are called semidefinite programming problems (SDP). They form a well known class of convex optimization problems that generalize linear programs and can be solved numerically even for large 5 problem matrices C. Numerical solvers also allow to optimize a linear objective function of the coefficients of C. Finally, we define the notion of strict feasibility.

Definition 4 (Strict Feasibility). The SDP in (10) is said to be strictly feasible when there exists a solution to the problem wherein the matrix C is positive definite.

If every feasible solution C to the problem ( 10) is positive semidefinite but not positive definite (in other words, the matrix has zero eigenvalues, or alternatively is rank deficient), the problem is said to fail strictly feasibility.

Remark 1. For a strictly feasible problem, there exist solutions C such that any C in a neighborhood from C and satisfying the equality constraints tr(A i C) = b i is also a solution. In contrary, problems that are not strictly feasible are also said to have an empty (relative) interior because, for any solution C, there exist C arbitrarily close from C that satisfy the equality constraints but are not solutions. This is illustrated on Fig 5.

{M | M 0} { C } C {M | tr(AiM ) = bi} {M | M 0} C C {M | tr(AiM ) = bi}

Verified SDPs

In Section 3, we laid out a procedure for formulating invariant checking and synthesis as a general SOS feasibility problem, which in turn is an SDP feasibility problem. There are SDPs with rational problem data whose solutions are irrational [START_REF] Tuncel | Polyhedral and semidefinite programming methods in combinatorial optimization[END_REF]. However, under some regularity conditions, SDP problems can theoretically be solved efficiently up to arbitrarily small error tolerance (see e.g., [START_REF] Wolkowicz | Handbook of semidefinite programming[END_REF]). In practice, many numerical solvers are available to solve SDP instances satisfactorily (see [START_REF] Anjos | Introduction to semidefinite, conic and polynomial optimization[END_REF]Part III]). Currently, the default choice for solving SDPs are specialized second order methods (i.e., using second order derivatives) called interior point methods (IPMs) 6 .

We discuss in Section 4.1 issues leading to inaccuracy or poor solution quality in SDP solving via IPMs. Then in Section 4.2, we consider solutions to guarantee soundness.

Sources of Solution Inaccuracy in Solving SDPs

Before we discuss ways to ensure soundness of solutions to the invariant checking and synthesis problem generated by SDP solvers, we first focus on a few issues that could possibly make an SDP solution inaccurate, leading to potential unsoundness. We will concentrate on SDP solutions obtained from general IPMs.

How do IPMs Work?

The convergence of a general IPM assumes strict feasibility (Def. 4). Using positive definite matrices as initial points, a general IPM repeatedly solves a perturbed linearization of the Karush-Kuhn-Tucker optimality conditions for a search direction, and moves along that search direction with a fractional step size that maintains the positive definiteness of the iterates. (See e.g. [START_REF] Tuncel | Polyhedral and semidefinite programming methods in combinatorial optimization[END_REF][START_REF] Wolkowicz | Handbook of semidefinite programming[END_REF].)

In the following, we discuss the four potential issues that can cause solution inaccuracy when a general IPM is used for obtaining SDP solutions: (1) inexact termination, (2) failure of strict feasibility, (3) ill-conditioning and (4) floating-point errors.

Inexact termination. The first source of inaccuracy stems from the fact that IPMs usually do not converge in finitely many iterations. Iterations are then stopped when some stopping criterion is met, for instance when the equalities in [START_REF] Borchers | CSDP, a C library for semidefinite programming[END_REF] are -approximated (i.e., |tr(A i C) -b i | ≤ ) or when the number of iterations becomes too big. Thus IPMs only produce approximate solutions. Nonetheless, under strict feasibility assumption, most common IPMs enjoy a convergence result of the following form: for any ε ∈ (0, 1), if an appropriate initial point is chosen, then it takes at most a number of steps polynomial in the problem size and log( 1 ε ) to obtain an ε-approximate solution. (See e.g. [START_REF] Wolkowicz | Handbook of semidefinite programming[END_REF]Ch. 10].) Failure of Strict Feasibility. Strictly feasibility is a desirable property. For instance, as seen above, it guarantees that a SDP can be solved to arbitrary accuracy by an IPM. "Random" SDP problems are strictly feasible with probability one [START_REF] Dür | The Slater condition is generic in linear conic programming[END_REF]Theorem 3.2]. However, strict feasibility can fail systematically for SDP instances arising from applications due to the inherent problem structure. In particular, strict feasibility can fail for entailment problems (Sect. 3.1), as shown in the following example. Example 6. The SDP feasibility problem in Example 5 fails the strict feasibility. Indeed, for any solution (c, C 0 , C 1 ), the equality constraints imply

(C 0 ) 1,1 = 0, hence 7 (C 0 ) 1,2 = (C 0 ) 2,1 = (C 0 ) 1,3 = (C 0 ) 3,1 = 0 which means that C 0 is rank deficient.
While the failure of the strict feasibility in small instances such as Example 6 usually does not cause much numerical issues, significant inaccuracy can often be observed as the number of variables and the degrees of the polynomials increase [START_REF] Waki | Strange behaviors of interiorpoint methods for solving semidefinite programming problems in polynomial optimization[END_REF]. Facial reduction techniques proposed by BORWEIN and WOLKOWICZ can be used for preprocessing SDP instances that are not strictly feasible [START_REF] Borwein | Facial reduction for a cone-convex programming problem[END_REF]. A more efficient version using linear programming reduction, called partial facial reduction, was proposed by PERMENTER and PARRILO [START_REF] Permenter | Partial facial reduction: simplified, equivalent sdps via approximations of the psd cone[END_REF].

Ill Conditioning. The coefficients of the polynomials in the entailment problems can influence the condition number of the linear system that is solved in IPMs, and a large condition number can affect the convergence of IPMs. While most SDP solvers use preconditioning to enhance the numerical stability, it is important to caution against the possible inaccuracy caused simply by the large input coefficients, which can occur even when preconditioning is used.

Example 7. Consider the entailment checking problem instance:

q1(x, y) x+y ≥ 0 ∧ q2(x, y) γ•(x 2 +y-1) ≥ 0 ∧ q3(x, y) x-4y 2 ≥ 0 |= p(x, y) ≥ 0, (11)
where p : (x 2 +y 2 )(q 1 (x, y)+q 2 (x, y)+q 3 (x, y)+8) and γ is a user-specified constant.

For any γ ∈ R, (11) is true, and the corresponding SOS problem has an obvious solution (σ i : x 2 + y 2 for i = 1, . . . , 3) that is independent of γ. Even though theoretically the solution set remains the same for varying γ, we see from Table 1 that a mere change in the value of γ can affect the solution accuracy in some SDP solvers: in this example, SDPT3 appears more robust against ill conditioning than SeDuMi. The large value of γ worsens the conditioning of the linear system solved in each iteration of an IPM and can lead to significant inaccuracy.

Table 1: The relative residual norm of the solutions returned by SDPT3 [START_REF] Reha H Tütüncü | Solving semidefinite-quadratic-linear programs using SDPT3[END_REF] and Se-DuMi [START_REF] Jos | Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[END_REF] for varying values of γ. Floating-Point Errors. For the sake of efficiency, IPMs are implemented using floatingpoint arithmetic. Thus, the precision of the floating-point format used limits the accuracy of the result. The most commonly used floating-point format offers a precision of about 10 -16 for arithmetic operations and SDP solvers usually offer accuracies around 10 -8 [START_REF] Borchers | CSDP, a C library for semidefinite programming[END_REF][START_REF] Yamashita | A high-performance software package for semidefinite programs: SDPA 7[END_REF]. Higher accuracies can be reached using more precise (and expensive) floating-point formats such as done by the SDPA-GMP solver (see [START_REF] Anjos | Introduction to semidefinite, conic and polynomial optimization[END_REF]Ch. 24] and [START_REF] Nakata | A numerical evaluation of highly accurate multiple-precision arithmetic version of semidefinite programming solver: SDPA-GMP, -QD and -DD[END_REF]).

Proving Soundness

Now we describe several different techniques for proving that SDP feasibility problems [START_REF] Borchers | CSDP, a C library for semidefinite programming[END_REF] arising from the SOS formulation of parametric entailment problems admit solutions. These techniques can be separated in two main approaches:

(a) Techniques that attempt to get an actual solution. They are able to solve some empty interior problems but this is often expensive. (b) Techniques that prove the existence of an actual solution, nearby to an approximate one. They require strict feasibility but are much cheaper.

After a quick review of the first approach, we detail the second one, since numerical tests in Sect. 5 indicate that it is the most useful one for proving polynomial invariants.

Deriving Exact Solutions As already mentioned in Section 3.2, the problem (10) is decidable. Unfortunately, even the most recent algorithms [START_REF] Henrion | Exact algorithms for linear matrix inequalities[END_REF] are not meant to be competitive with numerical solvers. Another approach consists in assuming the existence of a rational solution, using a numerical solver and attempting by various means to project its approximate solution to an exact rational solution [START_REF] Harrison | Verifying nonlinear real formulas via sums of squares[END_REF][START_REF] Kaltofen | Exact certification in global polynomial optimization via sums-of-squares of rational functions with rational coefficients[END_REF][START_REF] Monniaux | On the generation of positivstellensatz witnesses in degenerate cases[END_REF][START_REF] Peyrl | Computing sum of squares decompositions with rational coefficients[END_REF][START_REF] Platzer | Real world verification[END_REF]. These methods are truly impressive as they are able to solve some empty interior problems. It is also worth noting that since they provide an exact SOS decomposition, mechanically checking it with a proof assistant like Isabelle/HOL or Coq is particularly simple [START_REF] Harrison | Verifying nonlinear real formulas via sums of squares[END_REF][START_REF] Monniaux | On the generation of positivstellensatz witnesses in degenerate cases[END_REF].

Unfortunately they require heavy computations in rational arithmetic, which incurs the risk of an exponential blow-up of the size of the denominators.

Proving Existence of a Nearby Solution We now assume that ( 10) is strictly feasible, call a numerical solver that returns an approximate solution C and attempt to derive from it a proof that there exists an actual solution C (without actually computing C), based on the following proposition, whose proof is similar to that of [START_REF] Löfberg | Pre-and post-processing sum-of-squares programs in practice[END_REF]Theorem 4]. This suggests the following method to prove that a SOS problem is feasible:

Step 1. Obtain an approximate solution C.

Step 2. Compute (an overapproximation of) := max i∈{1,...,m} |tr(A i C) -b i |.

Step 3. Check that C -s I 0 (which implies s ≤ λ min ( C)).

Step 1 is achieved using a numerical solver and Step 2 is performed using floatingpoint interval arithmetic. The hard step is to provide a sound and efficient way to check C -s I 0. We rely on a check suggested by the following theorem. Let F be a floating-point format with unit roundoff eps and underflow unit eta. For any symmetric floating-point matrix M ∈ F s×s with 2(s + 2)eps < 1, define α : (s+1)eps 1-(2s+2)eps tr(M ) + 4(s + 1) (2(s + 2) + max i M i,i ) eta. -M ij = M ij , for any i = j; -M ii ≤ M ii -α, for any i; and the Cholesky algorithm implemented in floating-point arithmetic succeeds on M , i.e., "concludes" that M is positive semidefinite, Theorem 2 is used to prove that C -s I 0, as follows:

compute M := C -s I using floating-point arithmetic with rounding toward -∞.

It follows that the error ( C -s I) -M will be a diagonal matrix with nonnegative entries. Hence, if M 0 then C -s I 0, as well.

check that M is symmetric and that 2(s + 2)eps < 1; compute M := M -α I with rounding toward -∞ (the closest M to M -α I, the more likely its Cholesky decomposition is to succeed); compute the Cholesky decomposition of M .

If the Cholesky decomposition succeeds (which happens when, e.g., λ min ( C) ≥ s +2α [START_REF] Demmel | On floating point errors in Cholesky[END_REF]), then C -s I 0.

Remark 2. For the IEEE 754 [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF] binary64 format with rounding to nearest 8 , eps = 2 -53 ( 10 -16 ) and eta = 2 -1075 ( 10 -323 ). Thus, the hypothesis 2(s + 2)eps < 1 is always satisfied for practical values of s. Moreover, for typical values (s ≤ 1000 and elements of M of order of magnitude 1), α ≤ 10 -10 . This is negligible in front of s ∼ 10 -8 s (10 -8 being the typical default stopping tolerance), which means that the incompleteness of this positive definiteness check is not an issue in practice.

Steps 2 and 3 can be performed in only O(s 3 ) floating-point operations (cost of the Cholesky decomposition) so the cost of the whole method is dominated by the call to the numerical SDP solver in Step 1.

Remark 3. For ease of exposition, the above technique was presented on the whole matrix C, although it is preferable to apply it on each block C j of C.

Padding the SDP Problem. Naturally, all this requires that the least eigenvalue of the solution returned by the numerical solver be larger than s . It could seem that is known only after numerically solving the SDP problem, since it is computed from its result in Step 2. In fact, will be less than the stopping criterion of the solver, which is known in advance. Thus instead of solving [START_REF] Borchers | CSDP, a C library for semidefinite programming[END_REF], we solve the slightly modified problem The simple criterion in Proposition 1 assumes SDP problems translated from SOS problems. On the other hand, the tool VSDP [START_REF] Härter | VSDP: verified semidefinite programming[END_REF][START_REF] Jansson | Rigorous error bounds for the optimal value in semidefinite programming[END_REF] verifies the solutions of general SDP problems using interval arithmetic results. Remark 4. Mechanically checking proofs generated by the three step method of this section is an ongoing project. To this end, Theorem 2 has been verified [START_REF] Roux | Formal proofs of rounding error bounds[END_REF] in Coq.

Experiments

8 Type double in C.

x ∈ x ∈ R n ∧ k j=1 ij (x) ≥ 0 while (1) if (g(x) <= 0) x = τ1(x) else x = τ2(x)
Fig. 6: Benchmarks form.

This section presents an experimental evaluation of the methods described in Sections 3 and 4 on the examples of ADJ É et al. [START_REF] Adjé | Property-based Polynomial Invariant Generation Using Sums-of-Squares Optimization[END_REF]. We first synthesize polynomial invariants for these programs, following [START_REF] Adjé | Property-based Polynomial Invariant Generation Using Sums-of-Squares Optimization[END_REF], then attempt to formally prove their soundness. As seen in Section 2, these formal proofs are particularly worthwhile as synthesizing incorrect invariants is quite easy. Considered programs are of the form in Fig. 6. An invariant p(c, x) ≥ 0 can be provided by any solution 9 of the parametric entailment problem:

     i 1 (x) ≥ 0 ∧ • • • ∧ i k (x) ≥ 0 |= p(c, x) ≥ 0 g(x) ≤ 0 |= (p • τ 1 )(c, x) ≥ p(c, x) g(x) ≥ 0 |= (p • τ 2 )(c, x) ≥ p(c, x).
(See Section 3.) Thus, any solution of the following SOS problem gives an invariant:

find : polynomials σ j ∈ R d-di j [x] (j ∈ {1, . . . , k}), σ k+1 , σ k+2 ∈ R d-dg [x] s.t. p - k j=1 σ j i j is SOS, (p • τ 1 ) -p + σ k+1 g is SOS, (p • τ 2 ) -p -σ k+2 g is SOS, σ 1 , . . . , σ k+2 are SOS, ( 12
)
where d is the degree of p and d i1 , . . . , d i k and d g are the degrees of i 1 , . . . , i k and g respectively (all assumed to be less than d).

Table 2 gives the time needed to synthesize candidate invariants of degree d equal to 4, 6, 8 and 10 by solving the above SOS problem. "Example 4" in this table corresponds to the program of Fig. 2. The candidate invariant obtained for degree d = 6 is given in Fig. 2 and displayed in Fig. 3. The one obtained for d = 8 is also displayed in Fig. 3. "Example 8" corresponds to Fig. 4.

Unfortunately, the problem (12) usually has an empty interior 10 . This means that the candidate invariant obtained from numerical solvers does not precisely satisfy [START_REF] Chakarov | Deductive Proofs of Almost Sure Persistence and Recurrence Properties[END_REF]. In fact, there often exist values x 0 such that i 1 (x 0 ) ≥ 0, . . . , i k (x 0 ) ≥ 0 and p(x 0 ) is a tiny negative value. To fix that, we look for a small 11 c ∈ R such that p + c -k j=1 σ j i j is SOS for SOS polynomials σ j . This is done using the padding technique of Section 4. Times in table 2 include this fixing step.

We now attempt to prove that the fixed candidate invariants p are correct by considering the following entailment checking problem

i 1 (x) ≥ 0 ∧ • • • ∧ i k (x) ≥ 0 |= p(x) ≥ 0 g(x) ≤ 0 ∧ p(x) ≥ 0 |= (p • τ 1 )(x) ≥ 0 g(x) ≥ 0 ∧ p(x) ≥ 0 |= (p • τ 2 )(x) ≥ 0. (13
) 9 To get a "small" invariant, one minimizes the radius of the ball enclosing it [START_REF] Adjé | Property-based Polynomial Invariant Generation Using Sums-of-Squares Optimization[END_REF]. 10 Assignments τ often admit a fixpoint x0 = τ (x0) meaning that the condition (p•τ )-p+σ g ≥ 0 boils down in x0 to σ(x0) g(x0) ≥ 0 implying σ(x0) = 0 when g(x0) < 0. 11 In practice, c < 10 -3 when coefficients of p are of order of magnitude 1.

We first evaluated methods looking for exact solutions with the implementation of MONNIAUX and CORBINEAU [START_REF] Monniaux | On the generation of positivstellensatz witnesses in degenerate cases[END_REF]. Table 3 gives the results. The checking process is split in two parts: init for the initialization property (first entailment of ( 13)) and ind. for the inductiveness property (remaining entailments). As seen in the table, most of the initialization properties are indeed proved but proofs of the inductiveness property fail for all but the smallest example. This can be explained by the size of the corresponding SDP problems. For the initialization property, the largest block is a matrix of size . This is too much 12 to perform heavy computations with exact rational arithmetic.

Although [START_REF] Chakarov | Deductive Proofs of Almost Sure Persistence and Recurrence Properties[END_REF] usually has an empty interior, it is worth noting that this unfortunate property is due to the relaxation and is not intrinsic to the problem. Indeed, the loop body τ of the considered programs are usually strictly contractive, i.e., the image of the invariant {x | p(x) ≥ 0} by τ is included in its interior. When τ is continuous, this means that any polynomial p close enough from p also defines an invariant {x | p(x) ≥ 0}. In fact, the entailment checking problem (13) commonly leads to strictly feasible SDP problems. Thus, the method presented in Section 4.2 can be used to efficiently prove the soundness of a large part of the candidate invariants, as seen in Table 4. The time needed to compute the proofs (Table 4) is comparable to the time needed to synthesize the invariants (Table 2). Indeed, most of this time is spent running SDP solvers.

These results are confirmed by VSDP [START_REF] Härter | VSDP: verified semidefinite programming[END_REF][START_REF] Jansson | Rigorous error bounds for the optimal value in semidefinite programming[END_REF] when we provide it the SDP problems corresponding to [START_REF] Collins | Quantifier elimination for real closed fields by cylindrical algebraic decomposition[END_REF] and the strictly feasible solutions we computed using SDP solvers. This again indicates that these numerical verification methods only induce a very small overhead compared to the time required to run SDP solvers.

Implementation. The SOS to SDP translation described in Section 3, as well as the validation method described in Section 4.2 have been implemented in our OCaml library OSDP. It offers an interface to the SDP solvers Csdp [START_REF] Borchers | CSDP, a C library for semidefinite programming[END_REF], Mosek [START_REF] Mosek Aps | The MOSEK C optimizer API manual Version 7[END_REF], SDPA [START_REF] Yamashita | A high-performance software package for semidefinite programs: SDPA 7[END_REF] and SDPA-GMP [START_REF] Nakata | A numerical evaluation of highly accurate multiple-precision arithmetic version of semidefinite programming solver: SDPA-GMP, -QD and -DD[END_REF] and is available at http://cavale.enseeiht.fr/osdp/. Results from Tables 2 and4 have been obtained thanks to a small static analyzer relying on the library and available, along with all benchmarks, at http://cavale. enseeiht.fr/validatingSDP2016/. All computations were performed with the Mosek solver on a Xeon @ 2.67GHz.

Conclusion

Thus far, we have reviewed the use of SOS relaxations and numerical SDP solvers to solve polynomial problems arising in static analysis of programs. We presented some examples and experiments showing that, although erroneous results are often obtained from numerical solvers, rigorous proofs of soundness are possible. Moving forward, we wish to examine the application of our approach inside theorem provers and applications to hybrid systems, as well. Table 2: Time to synthesize candidate invariants for benchmarks [START_REF] Adjé | Property-based Polynomial Invariant Generation Using Sums-of-Squares Optimization[END_REF]. n is the number of variables and d τ the degree of the polynomial assignments. All times are in seconds, TO means timeout (900s) and MO out of memory (4GB).

Fig. 2 :

 2 Fig. 2: (Left) An example program taken from from ADJ É et al. [1] (Example 4).(Right) Purported invariant at loop head synthesized using SDP solvers[START_REF] Adjé | Property-based Polynomial Invariant Generation Using Sums-of-Squares Optimization[END_REF].

Fig. 4 :

 4 Fig. 4: (Left) An example program. (Right) Invariant p(x 1 , x 2 ) ≥ 0 , along with executions that start inside the initial set I (square).

Fig. 5 :

 5 Fig. 5: The line represents the equality constraints tr(A i C) = b i and the shaded area the matrices C 0. The set of solutions is the intersection of the line and the shaded area. (Left) A strictly feasible SDP problem. (Right) An empty interior problem.

9 SDPT3 2 . 1 ×

 921 10 -8 5.4 × 10 -10 5.1 × 10 -9 2.3 × 10 -8 SeDuMi 5.5 × 10 -9 2.6 × 10 -9 3.3 × 10 -5 0.00023

Proposition 1 .

 1 If (10) results from the SOS programming problem (5) or (6), and C ∈ R s×s satisfies the inequality s max i∈{1,...,m} |tr(A i C) -b i | ≤ λ min ( C) (the smallest eigenvalue of C), then (10) admits an actual solution C.

Theorem 2 .

 2 [START_REF] Rump | Verification of positive definiteness[END_REF] Corollary 2.4]) M 0 if there exists M ∈ F s×s such that the following conditions hold:

  find : C s.t. tr(A i C) = b i , i = 1, . . . , m, C -s I 0,which is an SDP (up to the change of variable C → C + s I).

n+ d 2 n × n+ d 2 n 2 n× n+ d dτ 2 n

 2222 whereas for inductiveness it is of size n+ d dτ

d = 4 d

 4 = 6 d = 8 d = 10 Example 4 (n = 2, dτ = 3) 0.25 0.95 3.31 8.85 Example 5 (n = 3, dτ = 2) 0.48 2.83 22.75 112.37 Example 6 (n = 4, dτ = 2) 2.12 64.07 TO MO Example 7 (n = 2, dτ = 3) 0.25 0.96 3.15 10.45 Example 8 (n = 2, dτ = 2) 0.17 0.34 0.74 1.93

  init ind. init ind. init ind. init ind. Example 4 (n = 2, dτ = 3) 1.43 NS 3.35 TO 19.80 MO 142.33 MO Example 5 (n = 3, dτ = 2) 3.82 TO 142.49 MO TO MO TO MO Example 6 (n = 4, dτ = 2) 32.20 TO TO MO ----Example 7 (n = 2, dτ = 3) 1.48 NS 3.36 TO 18.36 MO 120.40 MO Example 8 (n = 2, dτ = 2) 1.93 12.81 3.78 NS 26.29 TO 193.79 TO

d = 4 d = 6 d

 6 = 8 d = 10 init ind. init ind. init ind. init ind. Example 4 (n = 2, dτ = 3) 0.04 NS 0.06 NS 0.06 0.30 0.07 NS Example 5 (n = 3, dτ = 2) 0.06 0.18 0.09 0.26 0.16 0.80 0.27 2.52 Example 6 (n = 4, dτ = 2) 0.10 0.30 0.27 1.11 ----Example 7 (n = 2, dτ = 3) 0.05 NS 0.05 0.15 0.06 0.25 0.07 NS Example 8 (n = 2, dτ = 2) 0.04 0.07 0.03 NS 0.04 NS 0.05 NS

Table 3 :

 3 Checking the candidate invariants with the implementation of MONNIAUX and CORBINEAU[START_REF] Monniaux | On the generation of positivstellensatz witnesses in degenerate cases[END_REF]. All times are in seconds, NS means that no proof is found, TO means timeout (900s) and MO out of memory (4GB).

Table 4 :

 4 Checking the candidate invariants with the method of Section 4 (computing strictly feasible SDP solutions and verifying them). All times are in seconds. As seen in Section 2, counter-examples are easily found for Ex. 4, d = 4 and 6 and Ex. 7, d = 4. No such counter-examples were found for the other unproved cases and it remains unknown whether they are actually inductive or not.

		d = 4	d = 6	d = 8	d = 10
		init ind. init ind. init ind. init ind.
	Example 4 (n = 2, dτ = 3)	0.05 NS 0.07 NS 0.19 3.03 0.17 NS
	Example 5 (n = 3, dτ = 2)	0.08 0.33 0.23 2.20 0.74 14.55 2.50 92.15
	Example 6 (n = 4, dτ = 2)	0.22 1.52 1.26 38.94 ----
	Example 7 (n = 2, dτ = 3)	0.05 NS 0.07 0.85 0.19 3.32 0.17 NS
	Example 8 (n = 2, dτ = 2)	0.05 0.13 0.07 NS 0.09 NS 0.15 NS

Table 5 :

 5 Rechecking the proofs of Table4with VSDP[START_REF] Härter | VSDP: verified semidefinite programming[END_REF][START_REF] Jansson | Rigorous error bounds for the optimal value in semidefinite programming[END_REF] (verifying given strictly feasible SDP solutions). All times are in seconds.

In the remainder of this paper, the word "invariant" is used for inductive invariant.

Control theorists call (opposite of) such functions Lyapunov functions[START_REF] Féron | From control systems to control software[END_REF].

Typically, matrices Cj can be of size n × n for n up to a few hundreds.

There also exist first order methods handling larger problems but with less accurate solutions.

If a matrix is PSD and one of its diagonal entry (e.g. the (1, 1) entry) equals 0, then the entire row and column that contain that diagonal entry (e.g., the first row and column) equal 0.

Pierre Roux, Yuen-Lam Voronin, and Sriram Sankaranarayanan

For n = 2, dτ = 3 and d = 8, n+ d

Acknowledgments:

The authors would like to thank Didier Henrion, Pierre-Loïc Garoche and Assalé Adjé for interesting discussions on this subject.


Michael Röckner

Francesco Russo

Stochastic

Keywords: Fokker-Planck equation stochastic partial differential equations, infinite volume, porous media type equation, multiplicative noise, stochastic Fokker-Planck type equation. 2010 AMS-classification: 35R60, 60H15, 82C31

The purpose of the present note consists of first showing a uniqueness result for a stochastic Fokker-Planck equation under very general assumptions. In particular, the second order coefficients may be just measurable and degenerate.

We also provide a proof for uniqueness of a stochastic porous media equation in a fairly large space.

Introduction

We consider real functions e 0 , . . . , e N , • • • : R → R fulfilling Assumption 2.3 below.

In particular they are H -1 -multipliers, see Definition 2.1.

Let T > 0 and (Ω, F, P ), be a fixed probability space. Let (F t , t ∈ [0, T ]) be a filtration fulfilling the usual conditions and we suppose F = F T . Let µ(t, ξ), t ∈ [0, T ], ξ ∈ R, be a random field of the type

µ(t, ξ) = ∞ i=1 e i (ξ)W i t + e 0 (ξ)t, t ∈ [0, T ], ξ ∈ R,
where W i , i ≥ 1 are independent continuous (F t )-Brownian motions on (Ω, F, P ), which are fixed from now on until the end of the paper. For technical reasons we will sometimes set W 0 t ≡ t. We now consider a random field a and a deterministic function ψ as follows.

Assumption 1.1. a : [0, T ] × R × Ω → R + is a bounded progressively measurable random field. Assumption 1.2. ψ is monotone increasing, Lipschitz such that ψ(0) = 0. Let x 0 ∈ S ′ (R). We will consider the following two types of equations (1.1) and (1.2). The first one is a (linear) stochastic Fokker-Planck equation, the second one a stochastic porous media type equation, i.e.

   ∂ t z(t, ξ) = ∂ 2 ξξ ((az)(t, ξ)) + z(t, ξ)∂ t µ(t, ξ), z(0, • ) = x 0 , (1.1) and ∂ t X(t, ξ) = 1 2 ∂ 2 ξξ (ψ(X(t, ξ)) + X(t, ξ)∂ t µ(t, ξ), X(0, dξ) = x 0 .

(1.2)

They are both to be understood in the sense of (Schwartz) distributions. Their precise sense will be given in Remark 3.2 a) and in Definition 4.1. The stochastic multiplication above is of Itô type. In this paper we confine ourselves to the case of the underlying space being R 1 .

Fokker-Planck equations have been investigated until now in the deterministic framework, i.e. when e i = 0, i ≥ 1. There is a huge literature about existence and uniqueness in this case, see e.g. [START_REF] Bogachev | Fokker-Planck-Kolmogorov equations[END_REF] and references therein. More particularly, concerning uniqueness, in addition we draw the attention to Proposition 3.4 [START_REF] Blanchard | Probabilistic representation for solutions of a porous media type equation[END_REF] and Theorem 3.1 of [START_REF] Belaribi | Uniqueness for Fokker-Planck equations with measurable coefficients and applications to the fast diffusion equation[END_REF]. As far as we know this is the first time that a Fokker-Planck equation as (1.1) is considered in the literature, in particular for uniqueness, except for the unpublished work by the same authors [START_REF] Barbu | A stochastic Fokker-Planck equation and double probabilistic representation for the stochastic porous media type equation[END_REF]. We point out that we can allow degenerate coefficients in the second order term.

Concerning porous media equations, both in the deterministic and stochastic cases, there is a huge number of contributions, especially in finite volume. As far as the infinite volume case is concerned, in the deterministic situation a good framework is the classical Benilan-Crandall approach of the seventies; in the stochastic case some recent significant contributions have been made, see [START_REF] Ren | Stochastic generalized porous media and fast diffusion equations[END_REF][START_REF] Röckner | Non-monotone stochastic generalized porous media equations[END_REF][START_REF] Barbu | Stochastic porous media equations in R d[END_REF] and in particular [START_REF] Barbu | Stochastic Porous Media Equation[END_REF] and references therein. As mentioned, this paper draws however the attention on uniqueness for equations in the sense of distributions, within a large solutions class.

For instance, in the deterministic case a typical result in that sense is the paper [START_REF] Brezis | Uniqueness of solutions of the initial-value problem for u t -∆ϕ(u) = 0[END_REF] of Brezis and Crandall, which establishes uniqueness in the sense of distributions

in the class (L 1 ∩ L ∞ )([0, T ] × R d ).
Here we consider the equation (1.2) in the sense of distributions and we investigate uniqueness in the class of progressively measurable random fields The paper is organized as follows. After, this introduction and Section 2 devoted to preliminaries, in Section 3 the uniqueness Theorem 3.1 for an SPDE of Fokker-Planck type is formulated and proved. This, in turn, is an important ingredient for the probabilistic representation of a solution to a stochastic porous media type equation, see [START_REF] Barbu | A stochastic Fokker-Planck equation and double probabilistic representation for the stochastic porous media type equation[END_REF]. In the final section 4, using the same ideas as in Section 3 we prove a uniqueness result for (1.2), see Theorem 4.2.

X : Ω × [0, T ] × R such that [0,T ]×R X 2 (s, ξ)dsdξ < ∞ a.

Preliminaries

First we introduce some basic recurrent notations. M(R) denotes the space of signed Borel measures with finite total variation. We recall that S(R) is the space of the Schwartz fast decreasing test functions with its usual topology. S ′ (R) is its dual, i.e. the space of Schwartz tempered distributions. On S ′ (R), the map (I -∆)

s 2 , s ∈ R, is well-defined, via Fourier transform. For s ∈ R, H s (R) denotes the classical Sobolev space consisting of all functions f ∈ S ′ (R) such that (I -∆) s 2 f ∈ L 2 (R). We introduce the norm f H s := (I -∆) s 2 f L 2 ,
where • L p is the classical L p (R)-norm for 1 ≤ p ≤ ∞. In the sequel, we will often simply denote H -1 (R), by H -1 and L 2 (R) by L 2 . Furthermore, W r,p denotes the classical Sobolev space of order r ∈ N in L p (R) for 1 ≤ p ≤ ∞. Definition 2.1. Given a function e belonging to L 1 loc (R) ∩ S ′ (R), we say that it is an H -1 -multiplier, if the map ϕ → ϕe is continuous from S(R) to H -1 with respect to the H -1 -topology on both spaces. C(e) denotes the norm of this operator and we will call it multiplier norm. We remark that ϕe is always a well-defined Schwartz tempered distribution, whenever ϕ is a fast decreasing test function.

Remark 2.2. Let e : R → R. If e ∈ W 1,∞ (for instance if e ∈ W 2,1 ), then e is a H -1 (R)-multiplier.
Indeed, by duality arguments, to show this, it is enough to show the existence of a constant C(e) such that

eg H 1 C(e) g H 1 , ∀ g ∈ S(R).
(2.1)

Now (2.1) follows easily by the derivation product rules with for instance

C(e) = √ 2 e 2 ∞ + e ′ 2 ∞ 1 2 .
Here we fix some conventions concerning measurability. Any topological space E is naturally equipped with its Borel σ-algebra B(E). For instance B(R) (resp.

B([0, T ]) denotes the Borel σ-algebra of R (resp. [0, T ]).
In the whole paper, the following assumption on µ will be in force.

Assumption 2.3.

1. Each e i , i ≥ 0, belongs to the Sobolev space W 1,∞ .

2.

∞ i=1

(e i ) ′ 2 ∞ + e i 2 ∞ < ∞.
With respect to the random field µ, we introduce a notation for the Itô type stochastic integral below.

Let Z = (Z(s, ξ), s ∈ [0, T ], ξ ∈ R) be a random field on (Ω, F, (F t ), P ) such that T 0 R |Z(s, ξ)|dξ 2 ds < ∞ a.
s. and it is an L 1 (R)-valued, (F s )-progressively measurable process. Then, provided, Assumption 2.3 holds, the stochastic integral

[0,t]×R Z(s, ξ)µ(ds, ξ) := ∞ i=0 t 0 R e i (ξ)Z(s, ξ)dξ dW i s , t ≥ 0, is well-defined. More generally, if s → Z(s, •) is a measurable map [0, T ] × Ω → M(R), such that T 0 Z(s, •) 2 var ds < ∞, then the stochastic integral [0,t]×R Z(s, dξ)µ(ds, ξ) := ∞ i=0 t 0 R e i (ξ)Z(s, dξ) dW i s , t ≥ 0, is well-defined.
3 On the uniqueness of a Fokker-Planck type SPDE

The theorem below plays the analogous role as Theorem 3.8 in [START_REF] Blanchard | Probabilistic representation for solutions of a porous media type equation[END_REF] or Theorem 3.1 in [START_REF] Belaribi | Uniqueness for Fokker-Planck equations with measurable coefficients and applications to the fast diffusion equation[END_REF]. We recall that our Fokker-Planck SPDE has possibly degenerate measurable coefficients.

Theorem 3.1. We suppose that Assumptions 1.1 and 2.3 hold. Let z 1 , z 2 be two measurable random fields belonging ω a.s. to

C([0, T ], S ′ (R)) such that z 1 , z 2 :]0, T ]× Ω → M(R).
We moreover suppose the following.

i) z 1 -z 2 ∈ L 2 ([0, T ] × R) a.s. ii) t → (z 1 -z 2 )(t, •) is an (F t )-progressively measurable S ′ (R)-valued process.
iii) z 1 , z 2 are solutions to (1.1). such that

T 0 z i (s, •) 2 var ds < ∞ a.s.
Then

z 1 ≡ z 2 .
Remark 3.2. a) By a solution z of equation (1.1) we mean the following: for

every ϕ ∈ S(R), ∀t ∈ [0, T ], R ϕ(ξ)z(t, dξ) = x 0 , ϕ + t 0 ds R a(s, ξ)ϕ ′′ (ξ)z(s, dξ) + [0,t]×R ϕ(ξ)z(s, dξ)µ(ds, ξ) a.s. b) Let z = z 1 -z 2 . Since z is ω a.s. in L 2 ([0, T ]; L 2 (R)∩M(R)) ⊂ L 2 ([0, T ]; H -1 (R)), t 0 z(s, •)µ(ds, • ) belongs ω a.s. to C([0, T ]; H -1 (R))
and so also to C([0, T ]; H -2 (R)) ω a.s. On the other hand t 0 (az) ′′ (s, •)ds can be seen as a Bochner integral in H -2 (R). In particular any solutions z 1 , z 2 to (1.1) are such that z = z 1z 2 admits a modification whose paths belong (a.s.) to C([0, T ];

H -2 (R)) ∩ L 2 ([0, T ]; L 2 (R)). Since z i , i = 1, 2, are continuous with values in S ′ (R), their difference is indistinguishable with the mentioned modification. Consequently for ω a.s. z(t, •) ∈ C([0, T ]; H -2 (R)) and outside a P -null set N 0 , we have (in S ′ (R) and H -2 (R) ) z(t, •) = t 0 (az) ′′ (s, •)ds + t 0 z(s, •)µ(ds, •).
(3.1) c) By assumption i), possibly enlarging the P -null set N 0 we get the following.

For ω / ∈ N 0 , for almost all t ∈]0, T ],

t 0 (az)(s, •)ds ′′ ∈ H -1 (R) and so t 0 (az)(s, • )ds ∈ H 1 dt a.e. Proof of Theorem 3.1. Let z = z 1 -z 2 .
We fix the null set N 0 and so ω will always lie outside N 0 introduced in Remark 3.2 c). Let φ be a mollifier with compact support and φ ε = 1 ε φ( • ε ) be a generalized sequence of mollifiers converging to the Dirac delta function. We set

g ε (t) = z ε (t) 2 H -1 = R z ε (t, ξ)((I -∆) -1 z ε )(t, ξ)dξ, where z ε (t, ξ) = R φ ε (ξ -y)z(t, dy). Since t → z(t, • ) is continuous in H -2 (R), t → z ε (t, • ) is continuous in L 2 (
R) and so also in H -1 (R). We look at the equation fulfilled by z ε . The identity (3.1) produces the following equality in L 2 (R) and so in H -1 (R):

z ε (t, • ) = t 0 [(a(s, • )z(s, • )) ⋆ φ ε ] ′′ -(a(s, • )z(s, • )) ⋆ φ ε ds (3.2) + t 0 (a(s, • )z(s, • )) ⋆ φ ε ds + ∞ i=0 t 0 (e i z)(s, • ) ⋆ φ ε dW i s .
We apply (I -∆) -1 and we get

(I -∆) -1 z ε (t, • ) = - t 0 (a(s, • )z(s, • )) ⋆ φ ε ds (3.3) + t 0 (I -∆) -1 [(a(s, • )z(s, • )) ⋆ φ ε ] ds + ∞ i=0 t 0 (I -∆) -1 (e i z)(s, • ) ⋆ φ ε dW i s .
We apply Itô's formula to g ε . For a general introduction to infinite dimensional Hilbert space valued stochastic calculus, see [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF], [START_REF] Prévôt | A concise course on stochastic partial differential equations[END_REF] or [START_REF] Liu | Stochastic partial differential equations: an introduction[END_REF]. Taking into account, (3.2), (3.3) and that f, g H -1 = f, (I -∆) -1 g L 2 , we now obtain

g ε (t) = 2 t 0 z ε (s, • ), dz ε (s, • ) H -1 (3.4) + ∞ i=1 t 0 (e i z(s, • )) ⋆ φ ε , (e i z(s, • )) ⋆ φ ε H -1 ds = -2 t 0 z ε (s, • ), (a(s, • )z(s, • )) ⋆ φ ε L 2 ds + 2 t 0 z ε (s, • ), (I -∆) -1 ((a(s, • )z(s, • )) ⋆ φ ε ) L 2 ds + ∞ i=1 t 0 (e i z(s, • )) ⋆ φ ε , (e i z(s, • )) ⋆ φ ε H -1 ds + 2 t 0 z ε (s, • ), (e 0 z)(s, •) ⋆ φ ε H -1 ds + M ε t
where

M ε t = 2 ∞ i=1 t 0 z ε (s, • ), (e i z)(s, • ) ⋆ φ ε H -1 dW i s . (3.5)
Below we will justify that (3.5) is well-defined. We summarize (3.4) into

g ε (t) = g ε (t) + M ε t , t ∈ [0, T ].
We remark that

∞ i=1 T 0 z(s, • ), (e i z)(s, • ) H -1 2 ds = ∞ i=1 T 0 (e i z)(s, • ), (I -∆) -1 z(s, • ) L 2 2 ds (3.6) ≤ ∞ i=1 T 0 (e i z)(s, • ) 2 L 2 z(s, •) 2 H -2 ds ≤ ∞ i=1 e i 2 ∞ sup s∈[0,T ] z(s, •) 2 H -2 T 0 z(s, • ) 2 L 2 ds, because z : [0, T ] → H -2 is a.s. continuous by Remark 3.2 b).
Consequently,

M t = ∞ i=1 t 0 z(s, • ), (e i z)(s, • ) H -1 dW i s ,
is a well-defined local martingale.

It is also not difficult to show that for ε > 0,

∞ i=1 T 0 z ε (s, • ), (e i z)(s, • ) ⋆ φ ε 2 H -1 ds < ∞, a.s.
and so M ε defined in (3.5) is also a well-defined local martingale.

By assumption we have of course (ω

∈ N 0 ) [0,T ]×R (z ε (s, ξ) -z(s, ξ)) 2 dsdξ -→ ε→0 0, (3.7)
[0,T ]×R

((az) ⋆ φ ε -az) 2 (s, ξ)dsdξ -→ ε→0 0, (3.8)
[0,T ]×R

((e i z)(s, •) ⋆ φ ε -(e i z)(s, •)) 2 (ξ)dsdξ -→ ε→0 0, (3.9)
for every i ≥ 0, because z, az, e i z ∈ L 2 ([0, T ] × R), i ≥ 0. By usual estimates on convolutions, there is a universal constant C such that

[0,T ]×R ((e i z)(s, •)⋆φ ε ) 2 (ξ)dsdξ ≤ [0,T ]×R z 2 (s, ξ)(e i ) 2 (ξ)dsdξ ≤ C e i 2 ∞ z 2 L 2 ([0,T ]×R) . (3.10)
By Lebesgue dominated convergence theorem, using (3.9), it follows that (for ω / ∈

N 0 ), ∞ i=0 T 0 (e i z)(s, •) ⋆ φ ε -e i z(s, • ) 2 L 2 ds → ε→0 0. (3.11)
Using (3.7) and (3.11), it is not difficult to show that (for ω / ∈ N 0 )

∞ i=0 T 0 z ε (s, • ), (e i z)(s, •) ⋆ φ ε H -1 -z(s, • ), (e i z)(s, • ) H -1 2 ds (3.12) converges to zero. Now (for ω / ∈ N 0 ), ∞ i=0 T 0 (e i z)(s, • ) ⋆ φ ε 2 H -1 -(e i z)(s, • ) 2 H -1 ds (3.13) ≤ √ 2 ∞ i=0 T 0 (e i z)(s, • ) ⋆ φ ε -(e i z)(s, • ) 2 H -1 ds T 0 (e i z)(s, • ) ⋆ φ ε 2 H -1 + (e i z)(s, • ) 2 H -1 ds ≤ √ 2 ∞ i=0 T 0 (e i z)(s, • ) ⋆ φ ε -(e i z)(s, • ) 2 H -1 ds ∞ i=0 T 0 (e i z)(s, • ) ⋆ φ ε 2 H -1 + (e i z)(s, • )
g(t) = -2 t 0 z(s, • ), a(s, • )z(s, • ) L 2 ds (3.14) + 2 t 0 z(s, • ), (I -∆) -1 (a(s, • )z(s, • )) L 2 ds + 2 t 0 z(s, • ), (e 0 z)(s, • ) H -1 ds + ∞ i=1 t 0 (e i z)(s, • ), (e i z)(s, • ) H -1 ds.
The convergence of the second term of the right-hand side of (3.4) to the second term of the right-hand side of (3.14) holds again due to (3.7) and (3.8), cutting the difference in two pieces and using Cauchy-Schwarz's inequality. On the other hand the convergence of (3.12) to zero implies that M ε → M ucp, so that the ucp limit of g ε (t) + M ε t is equal to g(t) + M t . So, after a possible modification of the P -null set N 0 , setting g(t

) := z(t, • ) 2 H -1 , for ω / ∈ N 0 , we have g(t) + 2 t 0 z(s, • ), a(s, • )z(s, • ) L 2 ds (3.15) = 2 t 0 (I -∆) -1 z(s, • ), a(s, • )z(s, • ) L 2 ds + 2 t 0 ds z(s, • ), (e 0 z)(s, • ) H -1 ds + ∞ i=1 t 0 (e i z)(s, • ), (e i z)(s, • ) H -1 ds + M t .
By the inequality

2bc b 2 a ∞ + c 2 a ∞ , b, c ∈ R, it follows that 2 t 0 < (I -∆) -1 z(s, • ), (az)(s, • ) > L 2 ds a ∞ t 0 (I -∆) -1 z(s, • ) 2 L 2 ds + 1 a ∞ t 0 < (az)(s, • ), (az)(s, • ) > L 2 ds a ∞ t 0 z(s, • ) 2
H -2 ds + 1 a ∞ a ∞ t 0 z(s, • ), az(s, • ) L 2 ds. Since • H -2 • H -1 , (3.15) gives now (for ω / ∈ N 0 ), g(t) + t 0 z(s, • ), (az)(s, • ) L 2 ds M t + ∞ i=1 t 0 (e i z)(s, • ), (e i z)(s, • ) H -1 ds + 2 t 0 z(s, • ), (e 0 z)(s, • ) H -1 ds + a ∞ t 0 z(s, • ) 2 H -1 ds.
Since e i , i ≥ 0 are H -1 -multipliers with norm C(e i ), (for ω / ∈ N 0 )

g(t) + t 0 z(s, • ), (az)(s, • ) L 2 ds (3.16) M t + C t 0 z(s, • ) 2 H -1 ds = M t + C t 0 g(s)ds, ∀ t ∈ [0, T ],
where

C = ∞ i=1 C(e i ) 2 + 2C(e 0 ) + a ∞ .
We proceed now via localization which is possible because t → t 0 z(s, • ) 2 ds and t → z(t, •) H -2 are continuous P a.s. Let (ς ℓ ) be the sequence of stopping times

ς ℓ := inf{t ∈ [0, T ]| t 0 ds z(s, • ) 2 L 2 ≥ ℓ, z(t, •) 2 H -2
ℓ}.

(3.17)

If { } = ∅ we convene that ς ℓ = +∞. Clearly, the stopped processes M ς ℓ are (square integrable) martingales starting at zero. We evaluate (3.16) at t∧ς ℓ . Taking expectation we get

E(g(t ∧ ς ℓ )) E(M ς ℓ ∧t ) =0 +CE t∧ς ℓ 0 g(s)ds C t 0 dsE(g(s ∧ ς ℓ )). By Gronwall's lemma it follows that E(g(t ∧ ς ℓ )) = 0 ∀ ℓ ∈ N ⋆ . Since g is a.s.
continuous and lim ℓ→∞ t ∧ ς ℓ = T a.s., for every t ∈ [0, T ], by Fatou's lemma we get

E(g(t)) = E lim inf ℓ→∞ g(t ∧ ς ℓ ) lim inf ℓ→∞ E g(t ∧ ς ℓ ) = 0,
and the result follows.

4 Uniqueness for the porous media equation with noise

We first discuss first in which sense the SPDE (1.2) has to be understood.

Definition 4.1. A random field X = (X(t, ξ, ω), t ∈ [0, T ], ξ ∈ R, ω ∈ Ω) is said to be a solution to (1.
2) if P a.s. we have the following.

• X ∈ C([0, T ]; S ′ (R)) ∩ L 2 ([0, T ]; L 1 loc (R)).
• X is an S ′ (R) -valued (F t )-progressively measurable process.

• For any test function ϕ ∈ S(R) with compact support, t ∈]0, T ], we have

R X(t, ξ)ϕ(ξ)dξ = R ϕ(ξ)x 0 (dξ) + 1 2 t 0 ds R ψ(X(s, ξ, •))ϕ ′′ (ξ)dξ + [0,t]×R X(s, ξ)ϕ(ξ)µ(ds, ξ)dξ a.s.
We can state now the uniqueness theorem for the stochastic porous media equation.

Theorem 4.2. Suppose that Assumptions 1.2 and 2.3 hold. Then equation (1.2) admits at most one solution among the random fields

X :]0, T ] × R × Ω → R such that [0,T ]×R X 2 (s, ξ)dsdξ < ∞ a.s. (4.2)
Remark 4.3. Let X be a solution of (1.2) verifying (4.

2).

i) There is a P -null set N 0 , so that for ω ∈ N 0 , X(t,

• ) ∈ L 2 (R) for almost all t ∈ [0, T ].
ii) Condition (4.2) also implies that

T 0 X(s, • ) 2 H -1 ds < ∞ a.s.
iii) Since ψ is Lipschitz and ψ(0) = 0, (4.2) implies that

T 0 ψ(X(r, •)) 2 L 2 dr < ∞ a.s. So, t 0 dsψ(X(s, •)) is a Bochner integral with values in L 2 (R). iv) Consequently, t → ∆ t 0 ψ(X(s, • ))ds is continuous from [0, T ] to H -2 and
so also in S ′ (R); since e i , i ≥ 0, are H -1 -multipliers verifying Assumption 2.3, by Kolmogorov's lemma t → t 0 X(s, •)µ(ds, •) admits a version which belongs to C [0, T ]; H -1 (R) . Since x 0 ∈ S ′ (R) and X ∈ C([0, T ]; S ′ (R)) a.s., it follows that for ω not belonging to a null set, we have

X(t, • ) = x 0 + ∆   t 0 ψ(X(s, • ))ds   + t 0 X(s, •)µ(ds, •), t ∈ [0, T ],
as an identity in S ′ (R).

v) If x 0 ∈ H -1 , then X ∈ C [0, T ]; H -2 , for ω ∈ N 0 , N 0 a P -null set. vi) If x 0 ∈ H -s for some s ≥ 2, then X ∈ C ([0, T ]; H -s ), for ω ∈ N 0 , N 0 a P -null set.
vii) We consider a sequence of mollifiers (φ ε ) converging to the Dirac measure.

Then

X ε (t, • ) = X(t, •) ⋆ φ ε belongs a.s. to C [0, T ]; L 2 (R) . Remark 4.4. Since ψ is Lipschitz, there is α > 0 such that (ψ(r) -ψ(r)) (r -r) α (ψ(r) -ψ(r)) 2 .
Remark 4.5.

1. We note that condition 2. in Assumption 2.3 is more general than (3.1) of [START_REF] Barbu | Stochastic porous media equations in R d[END_REF], which can be reformulated here as follows. (c) (e i ) is an orthonormal system of H -1 and

∞ i=1 (e i ) ′ 2 ∞ + e i 2 ∞ + e i 2 H -1 < ∞.
2. An easy adaptation of Theorem 3.4 of [START_REF] Barbu | Stochastic porous media equations in R d[END_REF], in order to take into account e 0 , constitutes an existence result for (1.2): it says the following. Besides Assumptions 4.6 and 1.2, let us suppose moreover that x 0 ∈ L 2 or ψ is non-degenerate (i.e. ψ(x)

x ≥ 0, ∀x = 0). Then, there is a random field X such that

E [0,T ]×R X 2 (s, ξ)dsdξ < ∞, with t → t 0 ψ(X(s, •))ds ∈ C([0, T ]; H 1 (R)) a.s.
3. So, under the assumptions of item 2. above, the solution X is unique among those fulfilling (4.2).

Proof. Let (φ ε , ε > 0) be a sequence of mollifiers as in Remark 4.3 vii). Let X 1 , X 2 be two solutions of (1.2). For i = 1, 2, we set

(X i ) ε (t, •) = X i (t, •) ⋆ φ ε . We set X = X 1 -X 2 and X ε = (X 1 ) ε -(X 2 ) ε which a.s. belongs to C [0, T ]; L 2 (R) ⊂ C [0, T ]; H -1 . We set g ε (t) := X ε (t, • ) 2 H -1 = R (I -∆) -1 X ε (t, • ) (ξ)X ε (t, ξ)dξ.
Itô's formula gives

g ε (t) = 2 t 0 < X ε (s, • ), X ε (ds, • ) > H -1 + ∞ i=1 t 0 (e i X)(s, •) ⋆ φ ε 2 H -1 ds. (4.3)
On the other hand we have

X ε (t, • ) = t 0 ∆ ψ(X 1 (s, • )) -ψ(X 2 (s, • )) ⋆ φ ε ds + t 0 φ ε ⋆ (Xµ(ds, • )), (4.4
) where the notation of the latter integral is self-explanatory. So

(I -∆) -1 X ε (t, • ) = - t 0 ψ(X 1 (s, • )) -ψ(X 2 (s, • )) ⋆ φ ε ds (4.5) + t 0 (I -∆) -1 ψ(X 1 (s, • )) -ψ(X 2 (s, • )) ⋆ φ ε ds + ∞ i=0 t 0 (I -∆) -1 (e i X(s, • )) ⋆ φ ε dW i s .
We define

M t = ∞ i=1 t 0 < (I -∆) -1 X(s, • ), e i X(s, • ) > L 2 dW i s .
We observe that M is well-defined and it is a local martingale. Indeed, by Remark 4.3 v), X ∈ C([0, T ]; H -2 ). So by similar arguments as in (3.6),

∞ i=1 t 0 < (I -∆) -1 X(s, • ), e i X(s, • ) > 2 L 2 ds ≤ sup s∈[0,T ] X(s, •) 2 H -2 ∞ i=1 e i 2 ∞ T 0 X(s, • ) 2 L 2 ds < ∞. (4.6)
Using (4.3), (4.4) and (4.5) we get

g ε (t) = ∞ i=1 t 0 (e i X)(s, •) ⋆ φ ε 2 H -1 ds (4.7) -2 t 0 < X ε (s, • ), ψ(X 1 (s, • )) -ψ(X 2 (s, • )) ⋆ φ ε > L 2 ds + 2 t 0 < X ε (s, • ), (I -∆) -1 ψ(X 1 (s, • )) -ψ(X 2 (s, • )) ⋆ φ ε > L 2 ds + 2 t 0 < X ε (s, • ), (I -∆) -1 e 0 X(s, • ) ⋆ φ ε > L 2 ds + M ε t ,
where M ε is the local martingale defined by

M ε t = ∞ i=1 t 0 < X ε (s, • ), (I -∆) -1 e i X(s, • ) ⋆ φ ε > L 2 dW i s ,
which is again well-defined by similar arguments as in the proof of (4.6). Taking into account (4.2) and the Lipschitz property for ψ, we can take the limit when ε → 0 in (4.7) and for g(t) := X(t, • ) The proof is then completed by localization as in (3.17) at the end of Section 3.

ACKNOWLEDGMENTS

2 H - 1

 21 ds converges to zero, because of (3.10) and (3.11) and Assumption 2.3. Taking into account (3.7), (3.8), (3.12) and (3.13) we obtain (for ω / ∈ N 0 ), that lim ε→0 g ε (t) = g(t), t ∈ [0, T ], where

Assumption 4 . 6 .

 46 (a) e i ∈ W 1,∞ for every i ≥ 0.(b) e i , i ≥ 0, belong to H 1 .

2 H - 1 0 X 2 t 0 Xi=1 t 0 Xα t 0 ( 2 L 2 ds + α t 0 ψ X 1

 2102000201 , to obtaing(t) + 2 t (s, • ), ψ X 1 (s, • )ψ X 2 (s, • ) L 2 ds (∆) -1 X(s, • ), ψ X 1 (s, • )ψ X 2 (s, • ) L 2 ds + (s, • ), e 0 X(s, • ) H -1 ds + M t . The convergence M ε → M when ε → 0 is ucp, since ∞ ε (s, • ), [e i X(s, • )] ⋆ φ ε H -1 -X(s, • ), e i X(s, • ) H -1 2 ds -→ ε→0 0,which follows by similar arguments as in the proof of (3.12), using (4.2). 2ab a 2 α + b 2 α, for a, b ∈ R, α being the constant appearing at Remark 4.4, the second term of the right-hand side of equality (4.8) is bounded by 1 I -∆) -1 X(s, • ) (s, • )ψ X 2 (s, • )

Financial

  support through the SFB 701 at Bielefeld University and NSF-Grant 0606615 is gratefully acknowledged. The second named author benefited partially from the support of the "FMJH Program Gaspard Monge in optimization and operation research" (Project 2014-1607H). The authors are grateful to Viorel Barbu for stimulating discussions. The authors are grateful to the Referee whose comments have stimulated them to drastically improve the first version of the paper.

  (s, • )ψ X 2 (s, • ) , X(s, • ) L 2 ds. • ), ψ X 1 (s, • )ψ X 2 (s, • ) L 2 ds 2

	t ψ X 1 This together with (4.8) gives 2 H -1 ds + 0
		t
	g(t) +	0	X(s, t 0	X(s, • ) 2 H -1 ds
		∞	t
	+	i=1	0	(e i X)(s, •)
				∞ i=1	C(e i ) 2 +	1 α	)	0	t	g(s)ds.

t 0 X(s, • ), e 0 X(s, • ) H -1 ds + 1 α 2 H -1 ds + M t , t ∈ [0, T ] dP -a.s.

Since e i , i ∈ N, are H -1 -multipliers and taking into account Assumption 2.3, we get g(t) M t + (2C 0 +