Jean B Lasserre

THE MOMENT-SOS HIERARCHY

Keywords: 1991 Mathematics Subject Classification. 90C26 90C22 90C27 65K05 14P10 44A60 K-Moment problem, positive polynomials, global optimization, semidefinite relaxations

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Consider the optimization problem:

(1.1)

P : f * = inf x { f (x) : x ∈ Ω },
where f is a polynomial and Ω ⊂ R n is a basic semi-algebraic set, that is, (1.2) Ω := { x ∈ R n : g j (x) ≥ 0, j = 1, . . . , m }, for some polynomials g j , j = 1, . . . , m. Problem P is a particular case of Non Linear Programming (NLP) where the data (f, g j , j = 1, . . . , m) are algebraic, and therefore the whole arsenal of methods of NLP can be used for solving P. So what is so specific about P in (1.1)? The answer depends on the meaning of f * in (1.1).

If one is interested in a local minimum only then efficient NLP methods can be used for solving P. In such methods, the fact that f and g j 's are polynomials does not help much, that is, this algebraic feature of P is not really exploited. On the other hand if f * in (1.1) is understood as the global minimum of P then the picture is totally different. Why? First, to eliminate any ambiguity on the meaning of f * in (1.1), rewrite (1.1) as:

(1.3) P : f * = sup { λ : f (x) -λ ≥ 0, ∀x ∈ Ω } because then indeed f * is necessarily the global minimum of P.

In full generality, most problems (1.3) are very difficult to solve (they are labelled NP-hard in the computational complexity terminology) because:

Given λ ∈ R, checking whether "f (x) -λ ≥ 0 for all x ∈ Ω" is difficult.

Indeed, by nature this positivity constraint is global and therefore cannot be handled by standard NLP optimization algorithms which use only local information around a current iterate x ∈ Ω. Therefore to compute f * in (1.3) one needs an efficient tool to handle the positivity constraint "f (x)λ ≥ 0 for all x ∈ Ω". Fortunately if the data are algebraic then:

(1) Powerful positivity certificates from Real Algebraic Geometry (Positivstellensätze in german) are available. (2) Some of these positivity certificates have an efficient practical implementation via Linear Programming (LP) or Semidefinite Programming (SDP). In particular and importantly, testing whether a given polynomial is a sum of squares (SOS) simply reduces to solving a single SDP (which can be done in time polynomial in the input size of the polynomial, up to arbitrary fixed precision).

After the pioneers works of Shor [START_REF] Shor | Nondifferentiable Optimization and Polynomial Problems[END_REF] and Nesterov [START_REF] Nesterov | Squared functional systems and optimization problems[END_REF], Lasserre [START_REF] Lasserre | Optimisation globale et théorie des moments[END_REF][START_REF] Lasserre | Global Optimization with polynomials and the problem of moments[END_REF] and Parrilo [START_REF] Parrilo | Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization[END_REF][START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF] have been the first to provide a systematic use of these two key ingredients in Control and Optimization, with convergence guarantees. It is also worth mentioning another closely related pioneer work, namely the celebrated SDP-relaxation of Goemans & Williamson [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF] which provides a 0.878 approximation guarantee for MAXCUT, a famous problem in nonconvex combinatorial optimization (and probably the simplest one). In fact it is perhaps the first famous example of such a successful application of the powerful SDP convex optimization technique to provide guaranteed good approximations to a notoriously difficult non-convex optimization problem. It turns out that this SDP relaxation is the first relaxation in the Moment-SOS hierarchy (a.k.a. Lasserre hierarchy) when applied to the MAXCUT problem. Since then, this spectacular success story of SDP relaxations has been at the origin of a flourishing research activity in combinatorial optimization and computational complexity. In particular, the study of LP-and SDP-relaxations in hardness of approximation is at the core of a central topic in combinatorial optimization and computational complexity, namely proving/disproving Khot's famous Unique Games Conjecture 1 (UGC) in Theoretical Computer Science. Finally, another "definition" of the global optimum f * of P reads:

(1.4)

f * = inf µ { Ω f dµ : µ(Ω) = 1 }
1 For this conjecture and its theoretical and practical implications, S. Khot was awarded the prestigious Nevanlinna prize at the last ICM 2014 in Seoul [START_REF] Khot | Hardness of approximation[END_REF] .

where the 'inf" is over all probability measures on Ω. Equivalently, writing f as α f α x α in the basis of monomials (where

x α = x α 1 1 • • • x αn n ): (1.5) f * = inf y { α f α y α : y ∈ M (Ω); y 0 = 1 }, where M (Ω) = {y = (y α ) α∈N n : ∃ µ s.t. y α = Ω x α dµ, ∀α ∈ N n }, a convex cone. In fact (1.
3) is the LP dual of (1.4). In other words standard LP duality between the two formulations (1.4) and (1.3) illustrates the duality between the "Ω-moment problem" and "polynomials positive on Ω". Problem (1.4) is a very particular instance (and even the simplest instance) of the more general Generalized Problem of Moments (GPM):

(1.6) inf µ 1 ,...,µp { p j=1 Ω j f j dµ j : p j=1 f ij dµ j ≥ b i , i = 1, . . . , s },
for some functions f ij : R n j → R, i = 1, . . . , s, and sets Ω j ⊂ R n j , j = 1, . . . , p. The GPM is an infinite-dimensional LP with dual:

(1.7) sup

λ 1 ,...,λs≥0 { s i=1 λ i b i : f j - s i=1 λ i f ij ≥ 0 on Ω j , j : 1, . . . , p}.
Therefore it should be of no surprise that the Moment-SOS hierarchy, initially developed for global optimization, also applies to solving the GPM. This is particularly interesting as the list of important applications of the GPM is almost endless; see e.g. Landau [START_REF] Landau | Moments in Mathematics[END_REF].

The MOMENT-SOS hierarchy in optimization

2.1. Notation, definitions and preliminaries. Let R[x] denote the ring of polynomials in the variables x = (x 1 , . . . , x n ) and let R[x] d be the vector space of polynomials of degree at most d (whose dimension is s(d

) := n+d n ). For every d ∈ N, let N n d := {α ∈ N n : |α| (= i α i ) ≤ d}, and let v d (x) = (x α ), α ∈ N n , be the vector of monomials of the canonical basis (x α ) of R[x] d . Given a closed set X ⊆ R n , let P(X ) ⊂ R[x] (resp. P d (X ) ⊂ R[x] d )
be the convex cone of polynomials (resp. polynomials of degree at most 2d) that are nonnegative on d) in the canonical basis of monomials (x α ) α∈N n . For real symmetric matrices, let B, C := trace (BC) while the notation B 0 stands for B is positive semidefinite (psd) whereas B ≻ 0 stands for B is positive definite (pd).

X . A polynomial f ∈ R[x] d is written x → f (x) = α∈N n f α x α , with vector of coefficients f = (f α ) ∈ R s(
The Riesz functional. Given a sequence y = (y α ) α∈N n , the Riesz functional is the linear mapping L y : R[x] → R defined by: (2.1)

f (= α f α x α ) → L y (f ) = α∈N n f α y α .
Moment matrix. The moment matrix associated with a sequence y = (y α ), α ∈ N n , is the real symmetric matrix M d (y) with rows and columns indexed by N n d , and whose entry (α, β) is just y α+β , for every α, β ∈ N n d . Alternatively, let v d (x) ∈ R s(d) be the vector (x α ), α ∈ N n d , and define the matrices

(B o,α ) ⊂ S s(d) by (2.2) v d (x) v d (x) T = α∈N n 2d B o,α x α , ∀x ∈ R n . Then M d (y) = α∈N n 2d y α B o,α . If y has a representing measure µ then M d (y) 0 because f , M d (y)f = f 2 dµ ≥ 0, for all f ∈ R[x] d .
A measure whose all moments are finite, is moment determinate if there is no other measure with same moments. The support of a Borel measure µ on R n (denoted supp(µ)) is the smallest closed set Ω such that µ(R n \ Ω) = 0.

Localizing matrix. With y as above and g ∈ R[x] (with g(x) = γ g γ x γ ), the localizing matrix associated with y and g is the real symmetric matrix M d (g y) with rows and columns indexed by N n d , and whose entry (α, β) is just γ g γ y α+β+γ , for every α, β ∈ N n d . Alternatively, let B g,α ∈ S s(d) be defined by:

(2.3) g(x) v d (x) v d (x) T = α∈N n 2d+deg g B g,α x α , ∀x ∈ R n . Then M d (g y) = α∈N n 2d+degg y α B g,α .
If y has a representing measure µ whose support is contained in the set {x :

g(x) ≥ 0} then M d (g y) 0 for all d because f , M d (g y)f = f 2 gdµ ≥ 0, for all f ∈ R[x] d .
SOS polynomials and quadratic modules.

A polynomial f ∈ R[x] is a Sum-of-Squares (SOS) if there exist (f k ) k=1,...,s ⊂ R[x], such that f (x) = s k=1 f k (x) 2 , for all x ∈ R n . Denote by Σ[x] (resp. Σ[x] d
) the set of SOS polynomials (resp. SOS polynomials of degree at most 2d). Of course every SOS polynomial is nonnegative whereas the converse is not true. In addition, checking whether a given polynomial f is nonnegative on R n is difficult whereas checking whether f is SOS is much easier and can be done efficiently. Indeed let f ∈ R[x] 2d (for f to be SOS its degree must be even),

x → f (x) = α∈N n 2d f α x α . Then f is SOS if and only if there exists a real symmetric matrix X T = X of size s(d) = n+d n , such that: (2.4) X 0; f α = X, B o,α , ∀α ∈ N n 2d
, and this can be checked by solving an SDP.

Next, let x → g 0 (x) := 1 for all x ∈ R n . With a family (g 1 , . . . , g m )

⊂ R[x] is associated the quadratic module Q(g) (= Q(g 1 , . . . , g m )) ⊂ R[x]: (2.5) Q(g) :=    m j=0 σ j g j : σ j ∈ Σ[x], j = 0, . . . , m    ,
and its truncated version

(2.6) Q k (g) :=    m j=0 σ j g j : σ j ∈ Σ[x] k-d j , j = 0, . . . , m    ,
where d j = ⌈deg(g j )/2⌉, j = 0, . . . , m.

Definition 1. The quadratic module Q(g) associated with Ω in (1.2) is said to be Archimedean if there exists M > 0 such that the quadratic polynomial x → M -x 2 belongs to Q(g) (i.e., belongs to Q k (g) for some k).

If Q(g) is Archimedean then necessarily Ω is compact but the reverse is not rue. The Archimedean condition (which depends on the representation of Ω) can be seen as an algebraic certificate that Ω is compact. For more details on the above notions of moment and localizing matrix, quadratic module, as well as their use in potential applications, the interested reader is referred to Lasserre [START_REF] Lasserre | Positive Polynomials and Their Applications[END_REF], Laurent [START_REF] Laurent | Sums of squares, moment matrices and optimization over polynomials[END_REF], Schmüdgen [START_REF] Schmüdgen | The Moment problem[END_REF].

Two certificates of positivity (Positivstellensätze).

Below we describe two particular certificates of positivity which are important because they provide the theoretical justification behind the so-called SDP-and LPrelaxations for global optimization.

Theorem 2.1 (Putinar [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF]). Let Ω ⊂ R n be as in (1.2) and assume that

Q(g) is Archimedean. (a) If a polynomial f ∈ R[x] is (strictly) positive on Ω then f ∈ Q(g). (b) A sequence y = (y α ) α∈N n ⊂ R has a representing Borel measure on Ω if and only if L y (f 2 g j ) ≥ 0 for all f ∈ R[x]
, and all j = 0, . . . , m. Equivalently, if and only if M d (y g j ) 0 for all j = 0, . . . , m, d ∈ N.

There exists another certificate of positivity which does not use SOS.

Theorem 2.2 (Krivine-Vasilescu [START_REF] Krivine | [END_REF][START_REF] Krivine | Quelques propriétés des préordres dans les anneaux commutatifs[END_REF][START_REF] Vasilescu | Spectral measures and moment problems, in Spectral Theory and Applications[END_REF]). Let Ω ⊂ R n as in (1.2) be compact and such that (possibly after scaling) 0 ≤ g j (x) ≤ 1 for all x ∈ Ω, j = 1, . . . , m. Assume also that [1,

g 1 , . . . , g m ] generates R[x]. (a) If a polynomial f ∈ R[x] is (strictly) positive on Ω then (2.7) f (x) = α,β∈N n c α,β m j=1 g j (x) α j (1 -g j (x)) β j ,
g j (x) α j (1 -g j (x)) β j   ≥ 0 for all α, β ∈ N m .
The two facets (a) and (b) of Theorem 2.1 and Theorem 2.2 illustrate the duality between polynomials positive on Ω (in (a)) and the Ω-moment problem (in (b)). In addition to their mathematical interest, both Theorem 2.1(a) and Theorem 2.2(a) have another distinguishing feature. They both have a practical implementation. Testing whether

f ∈ R[x] d is in Q(g) k is
just solving a single SDP, whereas testing whether f van be written as in (2.7) with m i=1 α i + β i ≤ k, is just solving a single Linear Program (LP). 2.3. The Moment-SOS hierarchy. The Moment-SOS hierarchy is a numerical scheme based on Putinar's theorem. In a nutshell it consists of replacing the intractable positivity constraint "f (x) ≥ 0 for all x ∈ Ω" with Putinar's positivity certificate f ∈ Q d (g) of Theorem 2.1(a), i.e., with a fixed degree bound on the SOS weights (σ j ) in (2.6). By duality, it consists of replacing the intractable constraint y ∈ M (Ω) with the necessary conditions M d (g j y) 0, j = 0, . . . , m, of Theorem 2.1(b) for a fixed d. This results in solving an SDP which provides a lower bound on the global minimum. By allowing the degree bound d to increase, one obtains a hierarchy of SDPs (of increasing size) which provides a monotone non-decreasing sequence of lower bounds. A similar strategy based on Krivine-Stengle-Vasilescu positivity certificate (2.7) is also possible and yields a hierarchy of LP (instead of SDPs). However even though one would prefer to solve LPs rather than SDPs, the latter Moment-LP hierarchy has several serious drawbacks (some explained in e.g. [START_REF] Lasserre | An Introduction to Polynomial and Semi-Algebraic Optimization[END_REF][START_REF] Lasserre | Semidefinite programming vs. LP relaxations for polynomial programming[END_REF]), and therefore we only describe the Moment-SOS hierarchy.

Recall problem P in (1.1) or equivalently in (1.3) and (1.4), where Ω ⊂ R n is the basic semi-algebraic set defined in (1.2).

The Moment-SOS hierarchy. Consider the sequence of semidefinite programs

(Q d ) d∈N with d ≥ d := max[deg(f ), max j deg(g j )]: (2.8) Q d : ρ d = inf y { L y (f ) : y 0 = 1; M d (g j y) 0, 0 ≤ j ≤ m } (where y = (y α ) α∈N n 2d ) 2
, with associated sequence of their SDP duals:

(2.9) By duality, the sequence (Q * d ) d∈N forms a hierarchy of SDP-strenghtenings of (1.3). Indeed in (2.9) one has replaced the intractable positivity constraint of (1.3) by the (stronger) Putinar's positivity certificate with degree bound 2d -2d j on the SOS weights σ j 's.

Q * d : ρ * d = sup λ,σ j { λ : f -λ = m j=0 σ j g j ; σ j ∈ Σ[x] d-d j , 0 ≤ j ≤ m}
Theorem 2.3 ( [START_REF] Lasserre | Optimisation globale et théorie des moments[END_REF][START_REF] Lasserre | Global Optimization with polynomials and the problem of moments[END_REF]). Let Ω in (1.2) be compact and assume that its associated quadratic module Q(g) is Archimedean. Then:

(i) As d → ∞, the monotone non-decreasing sequence (ρ d ) d∈N (resp. (ρ * d ) d∈N
) of optimal values of the hierarchy (2.8) (resp. (2.9)) converges to the global optimum f * of P.

(ii) Moreover, let y d = (y d α ) α∈N n 2d be an optimal solution of Q d in (2.8), and let s = max j d j (recall that

d j = ⌈(deg g j )/2⌉). If (2.10) rank M d (y d ) = rank M d-s (y d ) (=: t)
then ρ d = f * and there are t global minimizers x * j ∈ Ω, j = 1, . . . , t, that can be "extracted" from y d by a linear algebra routine.

The sequence of SDP-relaxations (Q d ), d ≥ d, and the rank test (2.10) to extract global minimizers, are implemented in the GloptiPoly software [START_REF] Henrion | Gloptipoly 3: moments, optimization and semidefinite programming[END_REF].

Finite convergence and a global optimality certificate. After being introduced in [START_REF] Lasserre | Optimisation globale et théorie des moments[END_REF], in many numerical experiments it was observed that typically, finite convergence takes place, that is, f * = ρ d for some (usually small) d. In fact there is a rationale behind this empirical observation.

Theorem 2.4 (Nie [START_REF] Nie | Optimality Conditions and Finite Convergence of Lasserre's Hierarchy[END_REF]). Let P be as in (1.3) where Ω in (1.2) is compact and its associated quadratic module is Archimedean. Suppose that at each global minimizer x * ∈ Ω:

• The gradients (∇g j (x * )) j=1,...,m are linearly independent. (This implies existence of nonnegative Lagrange-KKT multipliers λ * j , j ≤ m, such that ∇f (x * ) -m j=1 λ * j ∇g j (x * ) = 0 and λ * j g j (x * ) = 0 for all j ≤ m.)

• Strict complementarity holds, that is, g j (x * ) = 0 ⇒ λ * j > 0. • Second-order sufficiency condition holds, i.e., u, ∇ 2 x (f (x * ) - m j=1 λ * j g j (x * )) u > 0, for all 0 = u ∈ ∇(f (x * ) -m j=1 λ * j g j (x * )) ⊥ . Then f -f * ∈ Q(g), i.e., there exists d * and SOS multipliers σ * j ∈ Σ[x] d * -d j , j = 0, . . . , m, such that: (2.11) f (x) -f * = σ * 0 (x) + m j=1 σ * j (x) g j (x).
With (2.11), Theorem 2.4 provides a certificate of global optimality in polynomial optimization, and to the best of our knowledge, the first at this level of generality. Next, observe that x * ∈ Ω is a global unconstrained minimizer of the extended Lagrangian polynomial f -f * -n j=1 σ * j g j , and therefore Theorem 2.4 is the analogue for non-convex polynomial optimization of the Karush-Kuhn-Tucker (KKT) optimality conditions in the convex case. Indeed in the convex case, any local minimizer is global and is also a global unconstrained minimizer of the Lagrangian f -f * -m j=1 λ * j g j . Also interestingly, whenever the SOS weight σ * j in (2.11) is non trivial, it testifies that the constraint g j (x) ≥ 0 is important for P even if it is not active at x * (meaning that if g j ≥ 0 is deleted from P then the new global optimum decreases strictly). The multiplier λ * j plays the same role in the KKT-optimality conditions only in the convex case. See [START_REF] Lasserre | An Introduction to Polynomial and Semi-Algebraic Optimization[END_REF] for a detailed discussion. Finite convergence of the Moment-SOS-hierarchies (2.8) and (2.9) is an immediate consequence of Theorem 2.4. Indeed by (2.11)

(f * , σ * 0 , . . . , σ * m ) is a feasible solution of Q * d * with value f * ≤ ρ * d ≤ f * (hence ρ * d = ρ d = f * )
. Genericity: Importantly, as proved in Nie [START_REF] Nie | Optimality Conditions and Finite Convergence of Lasserre's Hierarchy[END_REF], the conditions in Theorem 2.4 are generic. By this we mean the following: Consider the class P(t, m) of optimization problems P with data (f, g 1 , . . . , g m ) of degree bounded by t, and with nonempty compact feasible set Ω. Such a problem P is a "point" in the space R (m+1)s(t) of coordinates of (f, g 1 , . . . , g m ). Then the "good" problems P are points in a Zariski open set. Moreover, generically the rank test (2.10) is also satisfied at an optimal solution of (2.8) (for some d); for more details see Nie [START_REF] Nie | Certifying Convergence of Lasserre's Hierarchy via Flat Truncation[END_REF]. Computational complexity: Each relaxation Q d in (2.8) is a semidefinite program with s(2d) = n+2d n variables (y α ), and a psd constraint M d (y) 0 of size s(d). Therefore solving Q d in its canonical form (2.8) is quite expensive in terms of computational burden, especially when using interiorpoint methods. Therefore its brute force application is limited to small to medium size problems. Exploiting sparsity: Fortunately many large scale problems exhibit a structured sparsity pattern (e.g., each polynomial g j is concerned with a few variables only, and the objective function f is a sum i f i where each f i is also concerned with a few variables only). Then Waki et al. [START_REF] Waki | Sums of squares and semidefinite program relaxations for polynomial optimization problems with structured sparsity[END_REF] have proposed a sparsity-adapted hierarchy of SDP-relaxations which can handle problems P with thousands variables. In addition, if the sparsity pattern satisfies a certain condition then convergence of this sparsity-adapted hierarchy is also guaranteed like in the dense case [START_REF] Lasserre | Convergent SDP-relaxations in polynomial optimization with sparsity[END_REF]. Successful applications of this strategy can be found in e.g. Camps and Sznaier [START_REF] Camps | The interplay between Big-data and sparsity in systems identification[END_REF] in Control (systems identification) and in Molzahn and Hiskens [START_REF] Molzahn | Sparsity-Exploiting Moment-Based Relaxations of the Optimal Power Flow Problem[END_REF] for solving (large scale) Optimum Power Flow problems (OPF is an important problem encountered in the management of energy networks).

2.4. Discussion. We claim that the Moment-SOS hierarchy and its rationale Theorem 2.4, unify convex, non-convex (continuous), and discrete (polynomial) Optimization. Indeed in the description of P we do not pay attention to what particular class of problems P belongs to. This is in sharp contrast to the usual common practice in (local) optimization where several classes of problems have their own tailored favorite class of algorithms. For instance, problems are not treated the same if equality constraints appear, and/or if boolean (or discrete variables) are present, etc. Here a boolean variable x i is modeled by the quadratic equality constraint x 2 i = x i . So it is reasonable to speculate that this lack of specialization could be a handicap for the moment-SOS hierarchy.

But this is not so. For instance for the sub-class of convex 3 problems P where f and (-g j ) j=1,...,m are SOS-convex 4 polynomials, finite convergence takes place at the first step of the hierarchy. In other words, the SOS hierarchy somehow "recognizes" this class of easy problems [START_REF] Lasserre | An Introduction to Polynomial and Semi-Algebraic Optimization[END_REF]. In the same time, for a large class of 0/1 combinatorial optimization problems on graphs, the Moment-SOS hierarchy has been shown to provide the tightest upper bounds when compared to the class of lift-and-project methods, and has now become a central tool to analyze hardness of approximations in combinatorial optimization. For more details the interested reader is referred to e.g. Lasserre [START_REF] Lasserre | Semidefinite programming vs. LP relaxations for polynomial programming[END_REF], Laurent [START_REF] Laurent | A comparison of the Sherali-Adams, Lovász-Schrijver and Lasserre relaxations for 0-1 programming[END_REF], Barak [START_REF] Barak | Sum-of-Squares Proofs and the Quest toward Optimal Algorithms[END_REF], Khot [START_REF] Khot | Innapproximability of NP-complete Problems, Discrete Fourier Analysis, and Geometry[END_REF][START_REF] Khot | Hardness of approximation[END_REF] and the many references therein.

The Moment-SOS hierarchy outside optimization

3.1.

A general framework for the Moment-SOS hierarchy. Let Ω i ⊂ R n i be a finite family of compact sets, M (Ω i ) (resp. C (Ω i )) be the space of finite Borel signed measures (resp. continuous functions) on Ω i , i = 0, 1, . . . , s, and let T be a continuous linear mapping with adjoint T * :

T : M (Ω 1 ) × • • • × M (Ω s ) → M (Ω 0 ) C (Ω 1 ) × • • • × C (Ω s ) ← C (Ω 0 ) : T *
Let φ := (φ 1 , . . . , φ s ) and let φ i ≥ 0 stand for φ i is a positive measure. Then consider the general framework:

(3.1) ρ = inf φ≥0 { s i=1 f i , φ i : T(φ) = λ; s i=1 f ij , φ i ≥ b j , j ∈ J},
where J is a finite or countable set, b = (b j ) is given, λ ∈ M (Ω 0 ) is a given measure, (f ij ) j∈J , i = 1, . . . , s, are given polynomials, and •, • is the duality bracket between C (Ω i ) and M (Ω i ) ( h, φ i = Ω i hdφ i ), i = 1, . . . , s.

3 Convex problems P where f and (-gj)j=1,...,m are convex, are considered "easy" and can be solved efficiently.

4 A polynomial f ∈ R[x] is SOS-convex if its Hessian ∇ 2 f is a SOS matrix-polynomial, i.e., ∇f 2 (x) = L(x)L(x) T for some matrix-polynomial L ∈ R[x] n×p .
As we will see, this general framework is quite rich as it encompasses a lot of important applications in many different fields. In fact Problem (3.1) is equivalent to the Generalized Problem of Moments (GPM):

(3.2) ρ = inf φ≥0 { s i=1 f i , φ i : T * p k , φ = p k , λ , k = 0, 1, . . . s i=1 f ij , φ i ≥ b j , j ∈ J},
where the family (p k ) k=0,... is dense in C (Ω 0 ) (e.g. a basis of R[x 1 , . . . , x n 0 ]).

The Moment-SOS hierarchy can also be applied to help solve the Generalized Problem of Moments (GPM) (3.2) or its dual :

(3.3) ρ * = sup (θ j ≥0,γ) { k γ k p k , λ + θ, b : s.t. f i - k γ k (T * p k ) i - j∈J θ j f ij ≥ 0 on Ω i for all i },
where the unknown γ = (γ k ) k∈N is a finite sequence.

3.2.

A hierarchy of SDP-relaxations. Let (3.4) (3.5)

Ω i := { x ∈ R n i : g i,ℓ (x) ≥ 0, i = 1, . . . , m i }, i = 1, . . . ,
ρ d = inf (y i ) { s i=1 L y i (f i ) : s i=1 L y i ((T * p k ) i ) = p k , λ , p k ∈ Γ d s i=1 L y i (f ij ) ≥ b j , j ∈ J d M d (y i ), M d-d ℓ (g iℓ y i ) 0, ℓ ≤ m i ; i ≤ s},
where J d ⊂ J is finite d∈N J d = J. Its dual SDP-hierarchy reads:

(3.6)

ρ * d = sup (θ j ≥0,γ k ) { p k ∈Γ d γ k p k , λ + θ, b : s.t. f i - p k ∈Γ d γ k (T * p k ) i - j∈J θ j f ij = m i ℓ=0 σ i,ℓ g i,ℓ σ i,ℓ ∈ Σ[x 1 , . . . , x n i ] d-d i,ℓ ; i = 1, . . . , s},
As each Ω i is compact, for technical reasons and with no loss of generality, in the sequel we may and will assume that for every i = 1, . . . , s, g i,0 (x) = M i -x 2 , where M i > 0 is sufficiently large. Theorem 3.1. Assume that ρ > -∞ and that for every i = 1, . . . , s, f i0 = 1. Then for every d ≥ d, (3.5) has an optimal solution, and lim d→∞ ρ d = ρ.

Examples in Probability and Computational Geometry.

Bounds on measures with moment conditions. Let Z be a random vector with values in a compact semi-algebraic set Ω 1 ⊂ R n . Its distribution λ on Ω 1 is unknown but some of its moments x α dλ = b α , α ∈ Γ ⊂ N n , are known (b 0 = 1). Given a basic semi-algebraic set Ω 2 ⊂ Ω 1 we want to compute (or approximate as closely as desired) the best upper bound on Prob(Z ∈ Ω 2 ). This problem reduces to solving the GPM:

(3.7) ρ = sup φ 1 ,φ 2 ≥0 { 1, φ 2 : x α , φ 1 + x α , φ 2 = b α , α ∈ Γ; φ i ∈ M (Ω i ), i = 1, 2 },
With Ω 1 and Ω 2 as in (3.4) one may compute upper bounds on ρ by ing the Moment-SOS hierarchy (3.5) adapted to problem (3.7). Under the assumptions of Theorem 3.1, the resulting sequence (ρ d ) d∈N converges to ρ as d → ∞; for more details the interested reader is referred to [START_REF] Lasserre | Bounds on measures satisfying moment conditions[END_REF].

Lebesgue & Gaussian measures of semi-algebraic sets. Let Ω 2 ⊂ R n be compact. The goal is to compute (or approximate as closely as desired) the Lebesgue measure λ(Ω 2 ) of Ω 2 . Then take Ω 1 ⊃ Ω 2 be a simple set, e.g. an ellipsoid or a box (in fact any set such that one knows all moments (b α ) α∈N n of the Lebesgue measure on Ω 1 ). Then:

(3.8) λ(Ω 2 ) = sup φ 1 ,φ 2 ≥0 { 1, φ 2 : x α , φ 1 + x α , φ 2 = b α , α ∈ N n ; φ i ∈ M (Ω i ), i = 1,
ρ * d = inf p∈R[x] 2d { Ω 1 p dλ : p ≥ 1 on Ω 2 ; p ≥ 0 on Ω 1 },
i.e., as → ∞ one tries to approximate the discontinuous function x → 1 Ω 2 (x) by polynomials of increasing degrees. Fortunately there are several ways to accelerate the convergence, e.g. as in [START_REF] Henrion | Approximate volume and integration of basic semi-algebraic sets[END_REF] (but loosing the monotonicity) or in [START_REF] Lasserre | Computing Gaussian & exponential measures of semi-algebraic sets[END_REF] (preserving monotonicity) by including in (3.5) additional constraints on y 2 coming from an application of Stokes' theorem.

For the Gaussian measure λ we need and may take Ω 1 = R n and Ω 2 is not necessarily compact. Although both Ω 1 and Ω 2 are allowed to be 5 The Gibbs' phenomenon appears at a jump discontinuity when one approximates a piecewise C 1 function with a continuous function, e.g., by its Fourier series.

non-compact, the Moment-SOS hierarchy (3.5) still converges, i.e., ρ d → λ(Ω 2 ) as d → ∞. This is because the moments of λ satisfy the generalized Carleman's condition

(3.9) ∞ k=1 R n x 2k i dλ -1/2k
= +∞, i = 1, . . . , n, which imposes implicit constraints on y 1 and y 2 in (3.5), strong enough to guarantee ρ d → λ(Ω 2 ) as d → ∞. For more details see [START_REF] Lasserre | Computing Gaussian & exponential measures of semi-algebraic sets[END_REF]. This deterministic approach is computationally demanding and should be seen as complementary to brute force Monte-Carlo methods that provide only an estimate (but can handle larger size problems).

3.4.

In signal processing and interpolation. In this application, a signal is identified with an atomic signed measure φ supported on few atoms (x k ) k=1,...,s ⊂ Ω, i.e., φ = s k=1 θ k δ x k , for some weights (θ k ) k=1,...,s .

Super-Resolution. The goal of Super-Resolution is to reconstruct the unknown measure φ (the signal) from a few measurements only, when those measurements are the moments (b α ) α∈N n t of φ, up to order t (fixed). One way to proceed is to solve the infinite-dimensional program:

(3.10) ρ = inf φ { φ T V : x α dφ = b α , α ∈ N n t },
where the inf is over the finite signed Borel measures on Ω, and φ T V = |φ|(Ω) (with |φ| being the total variation of φ). Equivalently:

(3.11) ρ = inf φ + ,φ -≥0 { 1, φ + + φ -: x α , φ + -φ -= b α , α ∈ N n t },
which is an instance of the GPM with dual:

(3.12)

ρ * = sup p∈R[x]t { α∈N n t p α b α : p ∞ ≤ 1 },
where p ∞ = sup{|p(x)| : x ∈ Ω}. In this case, the Moment-SOS hierarchy (3.5) with d ≥ d := ⌈t/2⌉, reads:

(3.13) .ρ d = inf y + ,y - {y + 0 + y - 0 : y + α -y - α = b α , α ∈ N n t M d (y ± ) 0; M d (g ℓ y ± ) 0, ℓ = 1, . . . , m},
where Ω = {x : g ℓ (x) ≥ 0, ℓ = 1, . . . , m}.

In the case where Ω is the torus T ⊂ C, Candès and Fernandez-Granda [START_REF] Candès | Towards a Mathematical Theory of Super-Resolution[END_REF] showed that if δ > 2/f c (where δ is the minimal distance between the atoms of φ, and f c is the number of measurements) then (3.10) has a unique solution and one may recover φ exactly by solving the single semidefinite program (3.10) with d = ⌈t/2⌉. The dual (3.12) has an optimal solution p * (a trigonometric polynomial) and the support of φ + (resp. φ -) consists of the atoms z ∈ T of φ such that p * (z) = 1 (resp. p * (z) = -1). In addition, this procedure is more robust to noise in the measurements than Prony's method; on the other hand, the latter requires less measurements and no separation condition on the atoms.

In the general multivariate case treated in [START_REF] De Castro | Exact solutions to Super Resolution on semi-algebraic domains in higher dimensions[END_REF] one now needs to solve the Moment-SOS hierarchy (3.11) for d = d, . . . (instead of a single SDP in the univariate case). However since the moment constraints of (3.11) are finitely many, exact recovery (i.e. finite convergence of the Moment-SOS hierarchy (3.13)) is possible (usually with a few measurements only). This is indeed what has been observed in all numerical experiments of [START_REF] De Castro | Exact solutions to Super Resolution on semi-algebraic domains in higher dimensions[END_REF], and in all cases with significantly less measurements than the theoretical bound (of a tensorized version of the univariate case).

In fact, the rank condition (2.10) is always satisfied at an optimal solution (y + , y -) at some step d of the hierarchy (3.13), and so the atoms of φ + and φ -are extracted via a simple linear algebra routine (as for global optimization). Nie's genericity result [START_REF] Nie | Certifying Convergence of Lasserre's Hierarchy via Flat Truncation[END_REF] should provide a rationale which explains why the rank condition (2.10) is satisfied in all examples.

Sparse interpolation.

Here the goal is to recover an unknown (black-box) polynomial p ∈ R[x] t through a few evaluations of p only. In [START_REF] Josz | Sparse polynomial interpolation: compressed sensing, super resolution, or Prony?[END_REF] we have shown that this problem is in fact a particular case of Super-Resolution (and even discrete Super-Resolution) on the torus T n ⊂ C n . Indeed let z 0 ∈ T n be fixed, arbitrary. Then with β ∈ N n , notice that

p(z β 0 ) = α∈N n d p α (z β 1 01 • • • z βn 0n ) α = α∈N n d p α (z α 1 01 • • • z αn 0n ) β = T n z β d   α∈N n d p α δ z α 0   = T n z β dφ.
In other words, one may identify the polynomial p with an atomic signed Borel measure φ on T n supported on finitely many atoms (z α 0 ) α∈N n t with associated weights (p α ) α∈N n t . Therefore, if the evaluations of the black-box polynomial p are done at a few "powers" (z β 0 ), β ∈ N n , of an arbitrary point z 0 ∈ T n , then the sparse interpolation problem is equivalent to recovering an unknown atomic signed Borel measure φ on T n from knowledge of a few moments, that is, the Super-Resolution problem that we have just described above. Hence one may recover p by solving the Moment-SOS hierarchy (3.13) for which finite convergence usually occurs fast. For more details see [START_REF] Josz | Sparse polynomial interpolation: compressed sensing, super resolution, or Prony?[END_REF].

J * = inf u(t) T 0 L(x(t), u(t)) dt : ẋ(t) = f (x(t), u(t)), t ∈ (0, T ) x(t) ∈ X, u(t) ∈ U, ∀t ∈ (0, T ) x(0) = x 0 ; x(T ) ∈ X T ,
where L, f are polynomials, X, X T ⊂ R n and U ⊂ R p are compact basic semi-algebraic sets. In full generality the OCP problem (3.14) is difficult to solve, especially when state constraints x(t) ∈ X are present. Given an admissible state-control trajectory (t, x(t), u(t)), its associated occupation measure φ 1 up to time T (resp. φ 2 at time T ) are defined by:

φ 1 (A × B × C) := [0,T ]∩C 1 (A,B) ((x(t), u(t)) dt; φ 2 (D) = 1 D (x(T )), for all A ∈ B(X), B ∈ B(U), C ∈ B([0, T ]), D ∈ B(X T ). Then for every differentiable function h : X × [0, T ] → R h(T, x(T )) -h(0, x 0 ) = T 0 ( ∂h(x(t), u(t)) ∂t + ∂h(x(t), u(t)) ∂x f (x(t), u(t))) dt,
or, equivalently, with S := [0, T ] × X × U:

X T h(T, x) dφ 2 (x) = h(0, x 0 )+ S ( ∂h(x, u) ∂t + ∂h(x, u) ∂x f (x, u)) dφ 1 (t, x, u).
Then the weak formulation of the OCP (3.14) is the infinite-dimensional linear program:

(3.15)

ρ = inf φ 1 ,φ 2 ≥0 { S L(x, u) dφ 1 : s.t. X T h(T, •) dφ 2 - S ( ∂h ∂t + ∂h ∂x f ) dφ 1 = h(0, x 0 ) ∀h ∈ R[t, x] }.
It turns out that under some conditions the optimal values of (3.14) and (3.15) are equal, i.e., J * = ρ. Next, if one replaces "for all h ∈ R[t, x, u]" with "for all t k x α u β ", (t, α, β) ∈ N 1+n+p ", then (3.15) is an instance of the GPM (3.2). Therefore one may apply the Moment-SOS hierarchy (3.5).

Under the conditions of Theorem 3.1 one obtains the asymptotic convergence ρ d → ρ = J * as d → ∞. For more details see [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI-relaxations[END_REF] and the many references therein.

Robust control. In some applications (e.g. in robust control) one is often interested in optimizing over sets of the form:

G := {x ∈ Ω 1 : f (x, u) ≥ 0, ∀u ∈ Ω 2 },
where Ω 2 ⊂ R p , and Ω 1 ⊂ R n is a simple set, in fact a compact set such that one knows all moments of the Lebesgue measure λ on Ω 1 .

The set G is difficult to handle because of the universal quantifier. Therefore one is often satisfied with an inner approximation G d ⊂ G, and if possible, with (i) a simple form and (ii) some theoretical approximation guarantees. We propose to approximate G from inside by sets of (simple) form

G d = {x ∈ Ω 1 : p d (x) ≥ 0} where p d ∈ R[x] 2d .
To obtain such an inner approximation G d ⊂ G, define F :

Ω 1 → R, x → F (x) := min u {f (x, u) : u ∈ Ω 2 }. Then with d ∈ N, fixed, solve: (3.16) inf p∈R[x] 2d Ω 1 (F -p) dλ : f (x, u) -p(x) ≥ 0, ∀(x, u) ∈ Ω 1 × Ω 2 }.
Any feasible solution p d of (3. [START_REF] Josz | Sparse polynomial interpolation: compressed sensing, super resolution, or Prony?[END_REF]) is such that

G d = {x : p d (x) ≥ 0} ⊂ G.
In (3.16) Ω 1 (F -p) dλ = F -p 1 (with • 1 being the L 1 (Ω 1 )-norm), and

inf p Ω 1 (F -p) dλ = Ω 1 F dλ =cte + inf p Ω 1 -p dλ = cte -sup p Ω 1 p dλ
and so in (3.16) it is equivalent to maximize Ω 1 pdλ. Again the Moment-SOS hierarchy can be applied. This time one replaces the difficult positivity constraint f (x, u) -p(x) ≥ 0 for all (x, u) ∈ Ω 1 × Ω 2 with a certificate of positivity, with a degree bound on the SOS weights. That is, if

Ω 1 = {x : g 1,ℓ (x) ≥ 0, ℓ = 1, . . . ,
ρ d = sup p∈R[x] 2d Ω 1 p dλ : f (x, u) -p(x) = σ 0 (x, u) + m 1 ℓ=1 σ 1,ℓ (x, u) g i,ℓ (x) + m 2 ℓ=1 σ 2,ℓ (x, u) g i,ℓ (u)
σ i,ℓ ∈ Σ[x, u] d-d i,ℓ , ℓ = 1, . . . , m i , i = 1, 2.
Theorem 3.2 ( [START_REF] Lasserre | Tractable approximations of sets defined with quantifiers[END_REF]). Assume that Ω 1 × Ω 2 is compact and its associated quadratic module is Archimedean. Let p d be an optimal solution of (3.17 Therefore one obtains a nested sequence of inner approximations (G d ) d∈N ⊂ G, with the desirable property that λ(G \ G d ) vanishes as d increases. For more details the interested reader is referred to [START_REF] Lasserre | Tractable approximations of sets defined with quantifiers[END_REF].

Example 1. In some robust control problems one would like to approximate as closely as desired a non-convex set G = {x ∈ Ω 1 : λ min (A(x)) 0} for some real symmetric r × r matrix-polynomial A(x), and where x → λ min (A(x)) denotes its smallest eigenvalue. If one rewrites

G = {x ∈ Ω 1 : u T A(x)u ≥ 0, ∀u ∈ Ω 2 }; Ω 2 = {u ∈ R r : u = 1},
one is faced with the problem we have just described. In applying the above methodology the polynomial p d in Theorem 3.2 approximates λ min (A(x)) from below in Ω 1 , and p d (•) -λ min (A(•)) 1 → 0 as d increases. For more details see [START_REF] Henrion | Convergent relaxations of polynomial matrix inequalities and static output feedback[END_REF].

There are many other applications of the Moment-SOS hierarchy in Control, e.g. in Systems Identification [START_REF] Cerone | Set-membership error-in-variables identification through convex relaxation techniques[END_REF][START_REF] Camps | The interplay between Big-data and sparsity in systems identification[END_REF], Robotics [START_REF] Posa | Stability analysis and control of rigid-body systems with impacts and friction[END_REF], for computing Lyapunov functions [START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF], largest regions of attraction [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF], to cite a few.

3.6. Some inverse optimization problems. In particular:

Inverse Polynomial Optimization. Here we are given a polynomial optimization problem P : f * = min{f (x) : x ∈ Ω} with f ∈ R[x] d , and we are interested in the following issue: Let y ∈ Ω be given, e.g. y is the current iterate of a local minimization algorithm applied to P. Find (3.18) g * = arg min

g∈R[x] d { f -g 1 : g(x) -g(y) ≥ 0, ∀x ∈ Ω },
where

h 1 = α |h α | is the ℓ 1 -norm of coefficients of h ∈ R[x] d .
In other words, one searches for a polynomial g * ∈ R[x] d as close as possible to f and such that y ∈ Ω is a global minimizer of g * on Ω. Indeed if f -g * 1 is small enough then y ∈ Ω could be considered a satisfying solution of P. Therefore given a fixed small ǫ > 0, the test f -g * 1 < ǫ could be a new stopping criterion for a local optimization algorithm, with a strong theoretical justification. Again the Moment-SOS hierarchy can be applied to solve (3.18) as positivity certificates are perfect tools to handle the positivity constraint "g(x)g(y) ≥ 0 for all x ∈ Ω". Namely with Ω as in (1.2), solve:

(3.19) ρ t = min g∈R[x] d { f -g 1 : g(x) -g(y) := m j=0 σ j (x) g j (x), ∀x },
where g 0 (x) = 1 for all x, and σ j ∈ Σ[x] t-d j , j = 0, . . . , m. Other norms are possible but for the sparsity inducing ℓ 1 -norm • 1 , it turns out that an optimal solution g * of (3.19) has a canonical simple form. For more details the interested reader is referred to [START_REF] Lasserre | Inverse polynomial optimization[END_REF].

Inverse Optimal Control. With the OCP (3.14) in §3.5, we now consider the following issue: Given a database of admissible trajectories (x(t; x τ ), u(t, x τ )), t ∈ [τ, T ], starting in initial state x τ ∈ X at time τ ∈ [0, T ], does there exist a Lagrangian (x, u) → L(x, u) such that all these trajectories are optimal for the OCP problem (3.14)? This problem has important applications, e.g., in Humanoid Robotics to explain human locomotion [START_REF]Geometric and Numerical Foundations of Movements[END_REF].

Again the Moment-SOS hierarchy can be applied because a weak version of the Hamilton-Jacobi-Bellman (HJB) optimality conditions is the perfect tool to state whether some given trajectory is ǫ-optimal for the OCP (3.14). Indeed given ǫ > 0 and an admissible trajectory (t, x * (t), u * (t)), let ϕ : [0, T ] × X → R, and L : X × U → R, be such that: (ii) The above positivity constraint (3.20) on [0, T ] × X × U is replaced with a positivity certificate with degree bound on the SOS weights.

(3.20) ϕ(T, x) ≤ 0, ∀ x ∈ X; ∂ϕ(t, x) ∂t + ∂ϕ(t, x) ∂x f (x, u) + L(x, u) ≥ 0, for all (t, x, u) ∈ [0, T ] × X × U, and: ϕ(T, x * (T )) > -ǫ, (3.21
(iii) (3.21) is stated for every trajectory (x(t; x τ ), u(t, x τ )), t ∈ [τ, T ], in the database. Using a discretization {t 1 , . . . , t N } of the interval [0, T ], the positivity constraints (3.21) then become a set of linear constraints on the coefficients of the unknown polynomials ϕ and L.

(iv) ǫ in (3.21) is now taken as a variable and one minimizes a criterion of the form L 1 + γ ǫ, where γ > 0 is chosen to balance between the sparsityinducing norm L 1 of the Lagrangian and the error ǫ in the weak version of the optimality conditions (3.20)- (3.21). A detailed discussion and related results can be found in [START_REF] Pauwels | Linear Conic Optimization for Inverse Optimal Control[END_REF].

Optimal design in statistics.

In designing experiments one models the responses z 1 , . . . , z N of a random experiment whose inputs are represented by a vector t = (t i ) ∈ R n with respect to known regression functions Φ = (ϕ 1 , . . . , ϕ p ), namely: z i = p j=1 θ j ϕ j (t i ) + ε i , i = 1, . . . , N , where θ 1 , . . . , θ p are unknown parameters that the experimenter wants to estimate, ε i is some noise and the (t i )'s are chosen by the experimenter in a design space X ⊆ R n . Assume that the inputs t i , i = 1, . . . , N , are chosen within a set of distinct points x 1 , . . . , x ℓ ∈ X , ℓ ≤ N, and let n k denote the number of times the particular point x k occurs among t 1 , . . . , t N . A design ξ is then defined by: N Φ(x i ) Φ(x i ) T is called the information matrix of ξ. Optimal design is concerned with finding a set of points in X that optimizes a certain statistical criterion φ(M(ξ)), which must be real-valued, positively homogeneous, non constant, upper semi-continuous, isotonic w.r.t. Loewner ordering, and concave. For instance in D-optimal design one maximizes φ(M(ξ)) := log det(M(ξ)) over all ξ of the form (3.22). This is a difficult problem and so far most methods have used a discretization of the design space X .

The Moment-SOS hierarchy that we describe below does not rely an any discretization and works for an arbitrary compact basic semi-algebraic design space X as defined in (1.2). Instead we look for an atomic measure on X (with finite support) and we proceed in two steps: • In the first step one solves the hierarchy of convex optimization problems indexed by δ = 0, 1, . . .. where d is fixed by the number of basis functions ϕ j considered (here the monomials (x α ) α∈N n d ). (Note that (3.23) is not an SDP because the criterion is not linear in y, but it is still a tractable convex problem.) This provides us with an optimal solution y * (δ). In practice one chooses δ = 0.

• In a second step we extract an atomic measure µ from the "moments" y * (δ), e.g. via Nie's method [START_REF] Nie | The A-Truncated K-Moment Problem[END_REF] which consists of solving the SDP: M d+r (y) 0; M d+r-d j (g j y) 0}, where f r is a (randomly chosen) polynomial strictly positive on X . If (y * α (δ)) α∈N n 2d has a representing measure then it has an atomic representing measure, and generically the rank condition (2.10) will be satisfied. Extraction of atoms is obtained via a linear algebra routine. We have tested this two-steps method on several non-trivial numerical experiments (in particular with highly non-convex design spaces X ) and in all cases we were able to obtain a design. For more details the interested reader is referred to [START_REF] De Castro | Approximate Optimal Designs for Multivariate Polynomial Regression[END_REF].

Other applications & extensions. In this partial overview, by lack of space we have not described some impressive success stories of the Moment-SOS hierarchy, e.g. in coding [START_REF] Bachoc | New upper bounds for kissing numbers from semidefinite programming[END_REF], packing problems in discrete geometry [START_REF] De Laat | A semidefinite programming hierarchy for packing problems in discrete geometry[END_REF][START_REF] Schürmann | Computational approaches to lattice packing and covering problems[END_REF]. Finally, there is also a non-commutative version [START_REF] Pironio | Convergent relaxations of polynomial optimization problems with noncommuting variables[END_REF] of the Moment-SOS hierarchy based on non-commutative positivity certificates [START_REF] Helton | A Positivstellensatz for Noncommutative Polynomials[END_REF] and with important applications in quantum information [START_REF] Navascués | A convergent hierarchy of semidefinite programs characterizing the set of quantum correlations[END_REF].

Conclusion

The list of important applications of the GPM is almost endless and we have tried to convince the reader that the Moment-SOS hierarchy is one promising powerful tool for solving the GPM with already some success stories. However much remains to be done as its brute force application does not scale well to the problem size. One possible research direction is to exploit symmetries and/or sparsity in large scale problems. Another one is to determine alternative positivity certificates which are less expensive in terms of computational burden to avoid the size explosion of SOS-based positivity certificates.

(

  where d j =⌉(degg j )/2⌉). By standard weak duality in optimization ρ * d ≤ ρ d for every d ≥ d. The sequence (Q d ) d∈N forms a hierarchy of SDP-relaxations of P because ρ d ≤ f * and ρ d ≤ ρ d+1 for all d ≥ d. Indeed for each d ≥ d, the constraints of Q d consider only necessary conditions for y to be the moment sequence (up to order 2d) of a probability measure on Ω (cf. Theorem 2.1(b)) and therefore Q d is a relaxation of (1.5).

  s, for some polynomials (g i,ℓ ) ⊂ R[x 1 , . . . , x n i ], ℓ = 1, . . . , m i . Let d i,ℓ = ⌈deg(g i,ℓ )/2⌉ and d := max i,j,ℓ [deg(f i ), deg(f ij ), deg(g i,ℓ )]. To solve (3.2), define the "moment" sequences y i = (y i,α ), α ∈ N n i , i = 1, . . . , s, and with d ∈ N, define Γ d := {p k : deg(T * p k ) i ≤ 2d, i = 1, . . . , s}. Consider the hierarchy of semidefinite programs indexed by d ≤ d ∈ N:

3. 5 .

 5 In Control & Optimal Control. Consider the Optimal Control Problem (OCP) associated with a controlled dynamical system:(3.14)

  ). If λ({x ∈ Ω 1 : F (x) = 0}) = 0 then lim d→∞ F -p d 1 = 0 and lim d→∞ λ(G\G d ) = 0.

( 3 . 1 M

 31 [START_REF] Lasserre | Global Optimization with polynomials and the problem of moments[END_REF])ρ δ = sup y {log det(M d (y)) : y 0 = d+δ (y) 0; M d+δ-d j (g j y) 0},

  (3.24) ρ r = supy {L y (f r ) : y α = y * α (δ), ∀α ∈ N n 2d

  m 1 } and Ω 2 = {u : g 2,ℓ (u) ≥ 0, ℓ = 1, . . . , m 2 }, then with d i,ℓ := ⌈(deg(σ i,ℓ )/2⌉, one solves(3.17)

  for all t ∈ [0, T ]. Then the trajectory (t, x * (t), u * (t)) is an ǫ-optimal solution of the OCP (3.14) with x 0 = x * (0) and Lagrangian L. Therefore to apply the Moment-SOS hierarchy:(i) The unknown functions ϕ and L are approximated by polynomials in R[t, x] 2d and R[x, u] 2d , where d is the parameter in the Moment-SOS hierarchy(3.6).

) ∂ϕ(t, x * (t)) ∂t + ∂ϕ(t, x * (t)) ∂x f (x * (t), u * (t)) + L(x * (t), u * (t)) < ǫ,

In Theoretical Computer Science, y is called a sequence of "pseudo-moments".

Research supported by the European Research Council (ERC) through ERC-Advanced Grant # 666981 for the TAMING project.

1


Praveer Singh
email: praveer.singh@enpc.fr

Nikos Komodakis

CLOUD-GAN: CLOUD REMOVAL FOR SENTINEL-2 IMAGERY USING A CYCLIC CONSISTENT GENERATIVE ADVERSARIAL NETWORKS

Keywords: Cloud Removal, Generative Adversarial Networks, Deep Learning, Sentinel-2 Imagery

Cloud cover is a serious impediment in land surface analysis from Remote Sensing images either causing complete obstruction (thick clouds) with loss of information or blurry effects when being semi-transparent (thin clouds). While thick clouds require complete pixel replacement, thin cloud removal is fairly challenging as the atmospheric and landcover information is inter-twined. In this paper, we address this problem and propose a Cloud-GAN to learn the mapping between cloudy images and cloud-free images. The adversarial loss in the proposed method constrains the distribution of generated images to be close enough to the underlying distribution of the non-cloudy images. An additional cycle consistency loss is used to further restrain the generator to predict cloud-free images only of the same scene as reflected in the cloudy images. Our method not only rejects the necessity of any paired (cloud/cloud-free) training dataset but also avoids the need of any additional (expensive) spectral source of information such as Synthetic Aperture Radar imagery which is cloud penetrable. Lastly, we demonstrate the efficacy of our technique by training on an openly available and fairly new Sentinel-2 Imagery dataset consisting of real clouds. We also show significant improvement in PSNR values after removing clouds on synthetic images thus validating the competency of our methodology.

INTRODUCTION

Remote Sensing (RS) imagery is pivotal for varied challenging tasks such as recognizing footprints of buildings [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the inria aerial image labeling benchmark[END_REF], detecting changes in temporarily apart scenes [START_REF] Lyu | Learning a transferable change rule from a recurrent neural network for land cover change detection[END_REF] or semantic segmentation in aerial scenes [START_REF] Audebert | How useful is region-based classification of remote sensing images in a deep learning framework?[END_REF]. Such images are often plagued by films of clouds that partially or completely obstruct the scene. This can be quite annoying for RS experts, especially while observing a city like Paris which witnesses cloudy weather for a major part of the year. Thus, it clearly necessitates the requirement for an automatic technique that detects and removes the cloudy regions in a scene and replaces them with a neat in-painting of the underlying scene. Predicting a scene beneath a cloud is an under-constrained problem and unless we have some prior information, it is largely quite complex to replace clouds with correct underlying details. A way out has been by using multi-temporal images of the same region as done by [START_REF] Xu | Cloud removal based on sparse representation via multitemporal dictionary learning[END_REF] through a Multi-Temporal Dictionary Learning. [START_REF] Huang | Cloud removal from optical satellite imagery with sar imagery using sparse representation[END_REF] used Synthetic Aperture Radar (SAR) Imagery owing to the fact that it can easily penetrate through the clouds. However, SAR imagery is difficult to interpret and have a lower spatial resolution compared to RGB imagery. Additionally, [START_REF] Shen | An effective thin cloud removal procedure for visible remote sensing images[END_REF][START_REF] Xu | Thin cloud removal based on signal transmission principles and spectral mixture analysis," (a) Cloudy Cloud-Free (b) Cloudy Cloud-Free Fig. 6: Failure Cases: Over-smoothening or completely fail to produce images for overly clouded images[END_REF] studied the thin cloud removal problem in the literature. However, most of these approaches were based on conventional hand crafted methods and are limited in terms of performance.

Generative Adversarial Networks (GANs) [START_REF] Isola | Imageto-image translation with conditional adversarial networks[END_REF] have gained immense popularity owing to their remarkable capability in modeling the mapping function between input and output images belonging to target domains. Using an adversarial loss, the GAN's can be trained to produce fake images which are indistinguishable from the real images of target domain. [START_REF] Enomoto | Filmy cloud removal on satellite imagery with multispectral conditional generative adversarial nets[END_REF] used McGANs to predict cloud-free RGB images as well as cloud masks from the input cloudy image. The authors trained the model using pair of cloud-free image and synthetically produced cloudy images (by adding Perlin noise to the original RGB images). Additionally, they also utilize Near-Infrared (NIR) imagery, which is closer to visible range and possess partial cloud penetration capabilities. However, such kind of synthetic clouds produced using Perlin noise, are not realistic and significantly different from actual clouds seen in visible light images. Nevertheless, composing a dataset of real clouds and their cloud-free counterparts is quite a herculean task. We overcome this hurdle, by improvising upon a novel technique, which together with the adversarial loss from traditional GANs, employs a more recent cycle consistent loss [START_REF] Zhu | Unpaired image-to-image translation using cycle-consistent adversarial networks[END_REF], to convert thin cloudy images to cloud-free RGB images. Having a cycle consistency loss constrains the problem, such that if an image is transformed from input domain to target and then back to the input domain, it should look alike to the original image. An additional advantage of our method is that it absolves us from the requirement of an explicit paired cloudy/cloud-free dataset. Moreover, our methodology doesn't require any sort of cloud-penetration sources of imagery such as SAR or NIR. We simply utilize visible range imagery from a fairly new open source dataset (Sentinel-2) to report impressive results (figure 1 and3) clearly showcasing the efficacy of our results. We also report quantitative results on synthetically generated cloudy images, showcasing a significant improvement in PSNR values in Figure 5.

PROPOSED FRAMEWORK AND METHODOLOGY

Our proposed framework broadly consists of modeling a mapping function that can translate cloudy images into cloud-free visible range images. Assuming cloudy images belong to domain X and cloudfree images to domain Y, we define two mapping func-tions G : X -→ Y and F : Y -→ X which are modeled using two generator networks, Generator X2Y and Generator Y 2X , respectively, as illustrated in the figure 2. Training dataset is composed of {x i } N i=1 and {y i } M j=1 samples where x i ∈ X and y i ∈ Y and p data (x) and p data (y) are their data distributions respectively. Also, we have two Discriminator networks, Discriminator X (D X ) and Discriminator Y (D Y ) that distinguish between real data (x and y) and adversaries generated by the two generators i.e. G(x) and F(y).

Hence, in a way the generators and the discriminators compete with each other until they reach to a nash equilibrium. At this instant, the generator is effectively able to match the distribution of generated images (G(x) & G(y)) to the distribution of targeted images (p data (y)&p data (x) respectively) and thus, fool the discriminator. We specifically utilize Least Square GAN's (LSGAN's), which have shown [START_REF] Mao | Least squares generative adversarial networks[END_REF] to generate higher quality images with a much more stable learning process compared to regular GANs. The adversarial objective function is formulated as:

min D Y L LSGAN (D Y , X, Y ) = E y∼p data (y) [(D(y) -1) 2 ]+ E x∼p data (x) [D(G(x)) 2 ] min G L LSGAN (G, X, Y ) = E x∼p data (x) [(D(G(x)) -1) 2 ]
(1) Similarly, for mapping function F and discriminator D X , we have objective given by: min

D X L LSGAN (D X , Y, X), min F L LSGAN (F, Y, X) (2)
Adversarial objective alone is quite under-constrained as the same input can be mapped to any random permutation of target images. Thus, we additionally use a cycle consistency loss which restrains the generator to map a given input x i to a desired output y i . In Figure 2 (a), the loss constrains the generated Cyclic X to match Input X i.e. x -→ G(x) -→ F (G(x)) ≈ x, called as forward consistency. Similarly, for backward consistency, y -→ F (y) -→ G(F (y)) ≈ y. We formulate the cycle consistency objective as:

min G,F L cyc (G, F ) = E x∼p data (x) [ (F (G(x)) -x 1 ]+ E y∼p data (y) [ (G(F (y)) -y 1 ] (3)
Combining the Generator objectives from equations 1 and 2 and cyclic-consistency loss, our final Generator objective is:

min G,F L Gen (G, F ) = L LSGAN (G, X, Y )+ L LSGAN (F, Y, X) + λL cyc (G, F ) (4)
where λ is a regularizing factor that weights the cyclic term with respect to Generative term of the adversarial objective function.

TRAINING AND IMPLEMENTATION DETAILS

Dataset

Our dataset is composed of high resolution Level-1C Sentinel-2 imagery ranging between the year 2015 till 2017. Sentinel-2 is a multi-spectral dataset, with each spectral band is stored as a separate image [START_REF] Onera | Medusa toolbox[END_REF]. For our experiments, we choose images only from visible bands i.e. Blue (B2), Green (B3), Red (B4) all of which have 10 meters of spatial resolution.

Most of our cloud-free images are selected with 0-5% cover while for cloudy images we chose range anywhere between 10 to 100. All the images are downloaded over the Paris region as it is easy to get quite a range of cloudy images. We choose 20 cloudy and 13 cloudless images for training. We then extract 512 × 512 patches from these images. After filtering of unwanted ones, a total of 1677 patches for each cloud and cloud-free dataset were extracted while for testing we had 837 patches. For computational efficiency in training, we resize them to 256 × 256.

Network Architectures

We imbibe the architecture and the naming convention similar to what have been used by [START_REF] Zhu | Unpaired image-to-image translation using cycle-consistent adversarial networks[END_REF]. Generator architecture, uses 6 blocks for 128 × 128 training images and 9 blocks for 256 × 256 or higher resolution images. Additionally, a reflection pad is imbibed to avoid artifacts. The Discriminator architecture consists of a 70 × 70 PatchGAN [START_REF] Isola | Imageto-image translation with conditional adversarial networks[END_REF], classifying 70 × 70 patches as real or fake data. Thus, it can effectively be applied to any input size image and has lesser number of parameters.

Training

Initialization of weights was done through a Gaussian distribution with mean 0 and standard deviation 0.02. Optimization was carried out using ADAM [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], with a batch size of 1 and λ = 10 for all experiments. We perform training from scratch using a learning rate of 0.0002 up-to 200 epochs. The learning rate was kept constant for the first 100 epochs after which it linearly decays to zero until the last epoch. Also, as illustrated in [START_REF] Zhu | Unpaired image-to-image translation using cycle-consistent adversarial networks[END_REF] model oscillations are avoided by using a history of generated images (50) rather than only one.

RESULTS AND DISCUSSIONS

We present the results obtained using Cloud-GAN in Figure 3 (for real clouds) and 4 (for synthetic clouds). Without using any corresponding Cloud-Free pair for a cloudy image, our Cloud-GAN efficiently removes thin clouds spread throughout a scene, as shown in row III, IV, V in Figure 3. More interestingly, it effectively detects small cloudy patches and replaces them with the underlying ground details, as depicted in row I and II in Figure 3. Cloud-GAN interestingly is able to retain finer details like patches of urban settlements, river, fields (row IV) while getting rid of the cloudy film. In some cases e.g., row II, the generated image from our method is more natural and visually more pleasing than the original image which is a byproduct of our method.

We cannot report any quantitative results on the real dataset since we lack paired cloudy-cloud-free images. However, to compensate for that, we report results on 5 synthetic scenes, composed by addition of Perlin noise to cloud free images. We provide the corresponding PSNR results in Figure 5. We observe that even though we trained our model on real dataset, our model substantially outperforms on all synthetic scenes by a significant margin. Figure 4 shows two of these test scenes. Note that we do not provide comparison with [START_REF] Enomoto | Filmy cloud removal on satellite imagery with multispectral conditional generative adversarial nets[END_REF], due to the unavailability of their code and dataset.

We additionally show some special instances of thick cloud (Figure 6) where our model fails to yield credible results. In Figure 6 row I, we see that the model contends by generating an over-smoothed image when the clouds are too opaque. In Figure 6 row II, the model fails completely to produce an image as the clouds have occupied most of the visible area. One of the reasons can be that the network finds no closest sample in the target dataset and hence predicts a spatially smooth region under the cloud or some random noise. The only way to solve this is by the addition of an extra source such as SAR images which can penetrate through these clouds and give us details of the underlying ground details.

CONCLUSIONS

We have proposed a novel technique to remove thin clouds from Sentinel-2 imagery. Our cloud removal technique is employed without using any explicit dataset of paired Cloudy/Cloud-Free images as performed in the past. Our model nullifies the requirement of creating and training on synthetic dataset which is not truly realistic. Our input sources are purely visible range images without any prerequisite for SAR or other cloud-penetration sources. However, presence of thick clouds do necessitate need of an additional high wavelength imagery to gain some knowledge about the underlying ground information which we leave as a topic of future research. Lastly, we report PSNR results on some synthetic test scenes where we see significant improvement in performance thus validating the efficacy of our results.

Fig. 1 :

 1 Fig. 1: Cloud-GAN can effectively remove clouds from thin cloudy satellite imagery without supervision using ground truth

  (a) Forward cycle consistency loss x -→ G(x) -→ F (G(x)) ≈ x (b) Backward cycle consistency loss y -→ F (y) -→ G(F (y)) ≈ y

Fig. 2 :

 2 Fig. 2: Network Architecture: Generator X2Y and Generator Y 2X represent mapping function G : X -→ Y & F : Y -→ X respectively. Discriminator X and Discriminator Y represent D X and D Y respectively.

Fig. 3 :

 3 Fig. 3: Qualitative Results on Real Cloud Dataset. Cloudy images are shown in Row I and Cloud-GAN generated cloud-free images in Row II.

Fig. 4 : 4 SceneFig. 5 :

 445 Fig. 4: Sample images from Synthetic dataset where (a) Scene-5 (b) Scene-4


Praveer Singh
email: praveer.singh@enpc.fr

Nikos Komodakis

EFFECTIVE BUILDING EXTRACTION BY LEARNING TO DETECT AND CORRECT ERRONEOUS LABELS IN SEGMENTATION MASK

Keywords: Deep learning, high-resolution imagery, semantic segmentation, structured prediction, building footprint extraction

Semantic segmentation is pivotal for remote sensing image analysis. Although existing segmentation techniques perform well on similar landscape images, their generalization capability on an entirely different landscape is extremely poor. One of the primary reasons is that they partially or wholly, neglect the underlying relationship that exist in the joint space of input and output variables. Thus, effectively they lack to impose structure in their output predictions which is necessary for successful segmentation. In this paper, we address this problem and propose a novel solution by modeling the joint distribution of input-output variable which in turn enforces some structure in the initial segmentation mask. To this end, we first detect erroneous labels, in the form of Error maps, in the initial building masks. These Error maps are then used to correct the corresponding erroneous labels through a replacement technique. We evaluate our methodology on the benchmark Inria Aerial Image Labeling dataset, which is a large scale high resolution dataset for building footprint segmentation. In contrast to previous methods, our predicted segmentation masks are much closer to ground truth, owning to the fact that they are able to effectively correct both the large errors as well as the blobby effects. We lastly perform on par with other state-of-the-arts, validating the efficacy of our technique.

INTRODUCTION

Recently, we have witnessed an explosion of petabytes of high-resolution remote sensing datasets such as Sentinel 1-5 [START_REF]Sentinel[END_REF], SpaceNet [START_REF]Digital Globe CosmiQWorks and NVIDIA[END_REF] and Inria Aerial Image Labeling dataset [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the inria aerial image labeling benchmark[END_REF]. Till now, these datasets have been manually annotated by experts. However, with the availability of such an enormous amount of high-resolution data, it is truly herculean to label all of them by hand. This necessitates the automatic segmentation of these remote sensing images to quickly and effectively detect varied points of interest such as roads, buildings, forests in an image for tasks such as urban scene planning, green cover monitory and other emergency relief operations such as floods, forest fires and cyclones. Deep learning methods have recently shown significant improvements on automatic semantic labeling tasks for classical datasets such as Vaihingen or Potsdam. [START_REF] Audebert | Fusion of heterogeneous data in convolutional networks for urban semantic labeling[END_REF][START_REF] Audebert | Joint Learning from Earth Observation and OpenStreetMap Data to Get Faster Better Semantic Maps[END_REF] fused semantic maps from multiple sources through a residual correction technique. Whereas, [START_REF] Maggiori | Recurrent neural networks to correct satellite image classification maps[END_REF] corrected shifts in OSM maps through an iterative refinement technique using Recurrent Neural networks. In addition, [START_REF] Marmanis | Classification with an edge: Improving semantic image segmentation with boundary detection[END_REF] used boundary detections to improve the semantic segmentation and report impressive performance on Vaihingen dataset.

Interestingly, [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the inria aerial image labeling benchmark[END_REF] showcased the drawback of models trained over such classical datasets, by highlighting their lack of generalization capability to other cities that are captured under different conditions. They, henceforth build a new dataset covering a much larger surface of the earth including both densely populated urban landscape as well as sparse alpine regions under varied illumination conditions. Subsequently, they split their train and test data such that they come from different cities. Finally, they also addressed problem of blobby-predictions i.e., curvy edges of building masks. By fusing feature maps from different levels of convolution network, both [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the inria aerial image labeling benchmark[END_REF] and [START_REF] Maggiori | High-resolution aerial image labeling with convolutional neural networks[END_REF] combine low level edge information with high level semantic rich information. For the same task, [START_REF] Bischke | Multi-task learning for segmentation of building footprints with deep neural networks[END_REF] uses a multi-task loss in order to approximate the distance transform and the semantic maps.

Lately, most of the aforementioned techniques have partially or fully neglected the idea of enforcing structure in the output label space. This mainly results in these blobby ef- fects due to naive up-sampling during the segmentation masks prediction. On the contrary, a few of those that do enforce structure, either do it through a novice human assumption about the structure of the output label space (in the form of hand engineered CRF pairwise potentials as in [START_REF] Maggiori | Recurrent neural networks to correct satellite image classification maps[END_REF]). Or they rely on semantic maps from other sources to refine the initial predicted labels [START_REF] Audebert | Joint Learning from Earth Observation and OpenStreetMap Data to Get Faster Better Semantic Maps[END_REF][START_REF] Audebert | Fusion of heterogeneous data in convolutional networks for urban semantic labeling[END_REF]. However, these refinements through residual correction can only correct small errors (such as on the boundaries) while leaving major segmentation errors intact.

Built on these observations, we argue that we need to learn a joint structure of input and output variables that can effectively predict the replacements for major segmentation errors and thus, rectify them. Inspired from [START_REF] Gidaris | Detect, replace, refine: Deep structured prediction for pixel wise labeling[END_REF], we henceforth, propose a joint input-output model which segments the building footprints in high resolution imagery. While [START_REF] Gidaris | Detect, replace, refine: Deep structured prediction for pixel wise labeling[END_REF] addresses a continuous variable assignment task (disparity estimation), we solve here a binary classification problem. Our proposed method is buildup of two steps. Firstly, our model predicts an error map E based upon the input image X and initial segmentation Y . We then update the labels in regions of high error probabilities with a new label prediction which in turn relies on X, Y and E for its decision. Through the error map E, we learn a joint distribution between input and output variables which further helps in enforcing structure in the final label prediction. As shown in Figure 2, our model refines the blobby effect by learning the underlying structure of the building footprints which is enforced using the error maps. Finally, this leads our method to perform on par with existing state-of-the-arts on Inria Image Labeling dataset [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the inria aerial image labeling benchmark[END_REF].

PROPOSED FRAMEWORK AND METHODOLOGY

Let's assume our initial input image to be X = x C×H×W i=1 , where C, H and W are the channel, height and width of X respectively. Similarly, Y = y 1×H×W i=1 , be the initial segmentation map. Our technique aims to model a joint relationship between input X and output variables (Y ) to rectify and produce a much more refined version of Y . This can be formulated as Y = F (X, Y ) where Y denote the updated segmentation mask after replacing erroneous labels with new labels.

As shown in Figure 1, our proposed methodology comprises of two major steps. First, we predict errors (E) occurring in the initial building segmentation mask Y with the help of input image X. Next, we utilize these error maps to correct those erroneous labels in Y , whose error probability is large enough (in E) by replacing them with predicted labels. In the following subsections, we explain these two steps in detail.

Error Detection

The error detection component (F e ) computes the probability map (E) to detect the erroneous labels in the initial segmentation mask Y , stated as:

E = F e (X, Y ). (1)
In other words, F e learns from the joint input output space of X & Y to predict an error probability score map where each pixel has a value between 0 and 1. Specifically, it predicts whether or not a particular label y i is erroneous, if so, y i gets replaced with a correct label in the next step. F e can be easily formulated as a deep neural network which requires as such no explicit auxiliary loss and can be learnt under a single umbrella of one loss between corrected segmentation and ground truth segmentation.

Erroneous Label Replacement

The updated label Y ' is a convex combination of the initial segmentation mask Y and updates from the replacement component denoted by F r . It is given as:

Y = E F r (X, Y, E) + (1 -E) Y, (2)
where represents element-wise product. The error map E generated from F e acts as a gateway to restrict F r so as to just focus on the erroneous labels of Y and replace them with the predicted ones. Similar to F e , F r can also be modeled using any deep learning architecture. If we restrict the E probability maps to 0 and 1, the forward pass of Replacement happens as:

Y = Y, if F e (X, Y ) = 0 F r (X, Y, E), if F e (X, Y ) = 1. (3)
This shows that only erroneous labels are being replaced while non-erroneous labels remain intact. For any end loss function L between corrected segmentation mask Y and ground truth labels, its gradient becomes:

dL dF r (.) = 0, if F e (X, Y ) = 0 dL dY , if F e (X, Y ) = 1. (4)
In a way, the gradients update F r only for those regions where erroneous labels are found in Y , thus restricting F r to pay attention and predict replacements only for these particular regions. Additionally, passing Error maps to the F r component, helps it to rely on the correct labels to predict replacements for the new erroneous labels. Altogether, F r improves and makes correction of these erroneous labels, by jointly reflecting upon the Error maps, the input X and the initial segmentation Y .

TRAINING AND IMPLEMENTATION DETAILS

Dataset

We evaluated our method on the Inria Aerial Image Labeling Dataset [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the inria aerial image labeling benchmark[END_REF] which consists of Airborne imagery of urban settlement over the United States and Austria. The entire dataset consists of two classes namely, building and not building. All the images are of size 5000 × 5000 and have a resolution of 30 cm with RGB bands. First 5 images from each class were chosen for validation, while the rest were used for training.

For testing, we upload our test results to the project webpage.

Network Architecture

We initially train a Fully Convolutional Network [START_REF] Shelhamer | Fully convolutional networks for semantic segmentation[END_REF] (FCN) adapted to Resnet-50 [START_REF] He | Deep residual learning for image recognition[END_REF] architecture to generate our initial segmentation mask Y . This model, which we name as Resnet50FCN, is trained to reconstruct ground truth segmentation masks with a given input image X and ground truth labels. We henceforth treat Resnet50FCN as our baseline.

For detection, F e is implemented by using 5 convolutional layers (except last one, each is followed by batch-norm and Relu). The last conv. layer is followed by a soft-max, thus yielding us E between 0 and 1. To follow the input image size, we add an up-sampling layer on top of the Error map E.

For replacement, F r is implemented through compression block (compresses to 1/64 of the input resolution) and decompression block (decompresses to 1/4 of input resolution). These are essentially residual blocks with parameterized skip connection between symmetric layers in decompression and compression blocks respectively.

For additional implementation details of detection and replacement modules, we refer the reader to Section 3.2 of [START_REF] Gidaris | Detect, replace, refine: Deep structured prediction for pixel wise labeling[END_REF].

Training

Using an L1 loss between ground truth and predicted output Y , we optimize using an adam solver [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], with β 1 = 0.9 and β 2 = 0.99. The learning rate starting from 10 -3 is decreases to 3×10 -4 at 12, 10 -4 at 18, 3×10 -5 at 24 and finally 10 5 at 28 epochs. We continue training until 32 epochs. Each epoch consists of 500 batch iterations and each batch consists of 16 training samples where each sample is of size 1024 × 1024. compare our method (ResnetDR) with the baseline (Resnet 50FCN) and other previous best performing results namely, MLP [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the inria aerial image labeling benchmark[END_REF] and Multi-Task Loss [START_REF] Bischke | Multi-task learning for segmentation of building footprints with deep neural networks[END_REF]. In both tables, we report the Intersection over Union (IoU) and Accuracy score (Acc.) of correct pixels in the segmentation mask 1 . As shown in Table 1, we outperformed the previous best method [START_REF] Bischke | Multi-task learning for segmentation of building footprints with deep neural networks[END_REF] on the validation set by a margin of 14.90% on IoU. While on our own baseline Resnet50FCN, we improve by a margin of 1.44%. On the test set too, we outperformed [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the inria aerial image labeling benchmark[END_REF] by 14% and our own baseline by 1.22%.

EXPERIMENTAL RESULTS AND DISCUSSIONS

Quantitative Results

We show the quantitative results in

Qualitative Results

For qualitative analysis, we report our results in Figure 3 where green patch represents predictions for each model, red the ground truth while yellow represents the overlay of ground truth and predictions. We observe that while a major improvement is seen from FCN-MLP to Resnet50FCN, it still is not able to perfectly correct the blobby effects. ResnetDR improves upon these blobby effects by not only refining the boundaries of the segmentations but also regularizing them to better reflect structure of building footprints. For e.g., in Figure 3, in the Chicago case, we note that for the refinery (zoom-in shown in d & e), ResnetDR yields a more structured output in the form of parallel prediction edges. Similarly, for Vienna and Tyrol, predictions shrinking-in / protruding-out from the roofs of the houses for Resnet50FCN are constrained in case of ResnetDR to follow the edges of the roof. All these cases prove that our technique learns the underlying structure of building footprints in the form of regularized and well-defined shapes with straight edges. Simultaneously, it learns to predict error maps that somewhat refine these initial segmentation masks in order to enforce strict rules that govern building footprint structures. 1 Due to the unavailability of the results on test-set by [START_REF] Bischke | Multi-task learning for segmentation of building footprints with deep neural networks[END_REF], we haven't reported them in Table 2

CONCLUSIONS

We present a novel technique for structured semantic segmentation of high resolution satellite imagery. To this end, we propose to learn the joint space of input-output variables. Subsequently, our method enforces structure in the form of predicting a much more regularized building footprint and hence, resolves to a large extent the problem of blobby effects as reported in the past methods. We compare our technique with the state-of-the-art methods on the Inria aerial image labeling dataset, where we perform on par with others.

Fig. 1 :

 1 Fig. 1: The Network Architecture.

Fig. 2 :

 2 Fig.2: Segmentation masks for different techniques along with the final ground truth (GT). WHile Resnet50FCN drastically improves over FCN-MLP[START_REF] Maggiori | Can semantic labeling methods generalize to any city? the inria aerial image labeling benchmark[END_REF], our method (ResnetDR) further corrects the fuzzy or "blobby" effects of Resnet50FCN (caused due to naive up-sampling) by simply enforcing the underlying structure of building footprint shapes.

Fig. 3 :

 3 Fig.3: Qualitative Results. Red segments are ground truth, green are predictions while yellow represents overlap of ground truth and predictions. The blue circles highlight regions with significant changes. We also showcase zoom-in for one of these regions as comparison.

Table 1 and

 1 Table 2 on both the validation and test dataset respectively, where we

	Method	Austin Chicago Kitsap Co. West Tyrol Vienna Overall
	FCN-MLP [1]	IoU Acc. 94.20 61.20	61.30 93.43	51.50 98.92	57.95 96.66	72.13 91.87	64.67 94.42
	Multi-Task [9]	IOU 76.76 Acc. 93.21	67.06 99.25	73.30 97.84	66.91 91.71	76.68 96.61	73.00 95.73
	Resnet50FCN	IoU Acc. 97.02 88.12	81.21 93.04	83.62 99.3	88.15 98.28	87.07 94.65	86.46 96.46
	ResnetDR	IOU 89.42 Acc. 97.37	83.56 94.03	84.57 99.34	89.07 98.42	88.30 95.12	87.90 96.87

Table 1 :

 1 Evaluation results on validation set

	Method	Bellin. Blooming. Inns.	S.Fo Tyrol-E Overall
	FCN-MLP [1]	IoU Acc. 95.37 56.11	50.40 95.27	61.03 61.38 95.37 87.00	62.51 96.61	59.31 93.93
	Resnet50FCN	IoU Acc. 95.90 63.34	63.20 96.61	76.07 74.91 97.18 91.67	77.58 98.01	72.07 95.78
	ResnetDR	IOU Acc. 95.99 64.27	65.85 96.also52 97.30 92.01 77.10 75.86	78.68 98.11	73.30 95.99

Table 2 :

 2 Evaluation results on test set.