Didier Henrion

Jean B Lasserre

Martin Mevissen

Mean squared error minimization for inverse moment problems *

Keywords: Moment problems, density estimation, inverse problems, semidefinite programming

We consider the problem of approximating the unknown density u ∈ L 2 (Ω, λ) of a measure µ on Ω ⊂ R n , absolutely continuous with respect to some given reference measure λ, from the only knowledge of finitely many moments of µ. Given d ∈ N and moments of order d, we provide a polynomial p d which minimizes the mean square error (u -p) 2 dλ over all polynomials p of degree at most d. If there is no additional requirement, p d is obtained as solution of a linear system. In addition, if p d is expressed in the basis of polynomials that are orthonormal with respect to λ, its vector of coefficients is just the vector of given moments and no computation is needed. Moreover p d → u in L 2 (Ω, λ) as d → ∞. In general nonnegativity of p d is not guaranteed even though u is nonnegative. However, with this additional nonnegativity requirement one obtains analogous results but computing p d ≥ 0 that minimizes (u -p) 2 dλ now requires solving an appropriate semidefinite program. We have tested the approach on some applications arising from the reconstruction of geometrical objects and the approximation of solutions of nonlinear differential equations. In all cases our results are significantly better than those obtained with the maximum entropy technique for estimating u.

Introduction

Estimating the density u of an unknown measure µ is a well-known problem in statistical analysis, physics or engineering. In a statistical context, one is usually given observations in the form of a sample of independent or dependent identically distributed random variables obtained from the unknown measure µ. And so there has been extensive research on estimating the density based on these observations. For instance, in one of the most popular approaches, the kernel density estimation [START_REF] Parzen | On estimation of a probability density function and mode[END_REF], the density u is estimated via a linear combination of kernel functions -each of them being identified with exactly one observation. The crucial step in this method is to choose an appropriate bandwidth for which minimizing the integrated or the mean-integrated squared error between u and its estimate is most common. Another very popular approach uses wavelets [START_REF] Kerkyacharian | Density estimation by kernel and wavelet methods: optimality of Besov spaces[END_REF][START_REF] Donoho | Density estimation by wavelet thresholding[END_REF][START_REF] Vannucci | Nonparametric density estimation using wavelets[END_REF], an example of approximating a density by a truncated orthonormal series. The coefficients in the truncated wavelet expansion are moment estimates derived from the given identically distributed observations. Again, the approximation accuracy is often measured by the mean-integrated squared error and depends on the number of observations and the degree of the truncation. This approach provides a global density estimate satisfying both good local and periodic approximation properties. For further details the interested reader is referred to [START_REF] Eloyan | Smooth density estimation with moment constraints using mixture distributions[END_REF][START_REF] Izenman | Recent developments in nonparametric density estimation[END_REF][START_REF] Vannucci | Nonparametric density estimation using wavelets[END_REF] and the many references therein.

In another context -arising in challenging fields such as image recognition, solving nonlinear differential equations, spectral estimation or speech processing -no direct observation is available, but rather finitely many moments of the unknown measure µ are given. Then the issue is to reconstruct or approximate the density u based on the only knowledge of finitely many moments, (say up to order d ∈ N), an inverse problem from moments. A simple method due to [START_REF] Teague | Image analysis via the general theory of moments[END_REF] approximates the density u by a polynomial p of degree at most d, so that the moments of the measure pdλ matches those of µ, up to order d. However, and in contrast with more sophisticated approaches, the resulting polynomial approximation p is not guaranteed to be a density (even though u is) as it may takes negative values on the domain of integration. One classical approach to the moment problem is the Padé approximation [START_REF] Baker | Padé approximants[END_REF] which is based on approximating the measure by a (finite) linear combination of Dirac measures. The Dirac measures and their weights in the decomposition are determined by solving a nonlinear system of equations. In the maximum entropy estimation (another classical approach) one selects the best approximation of u by maximizing some functional entropy, the most popular being the Boltzmann-Shannon entropy. In general some type of weak convergence takes place as the degree increases as detailed in [START_REF] Borwein | On the convergence of moment problems[END_REF]. Alternatively the norm of the approximate density is chosen as an objective function [START_REF] Bertero | Linear inverse problems with discrete data. I: General formulation and singular system analysis[END_REF][START_REF] Talenti | Recovering a function from a finite number of moments[END_REF][START_REF] Jones | On extending the orthogonality property of minimum norm solutions in Hilbert space to general methods for linear inverse problems[END_REF][START_REF] Goodrich | L 2 spectral estimation[END_REF], which allows to show a stronger convergence in norm. In [START_REF] Mead | Maximum entropy in the problem of moments[END_REF], maximum entropy and Padé approximates have been compared on some numerical experiments. Finally, piecewise polynomial spline based approaches have also been proposed in [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF].

Motivation. Our main motivation to study the (inverse) moment problem arises in the context of the so-called generalized problem of moments (GPM). The abstract GPM is a infinite-dimensional linear program on some space of Borel measures on R n and its applications seem endless, see e.g. [START_REF] Landau | Moments in Mathematics[END_REF][START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] and the many references therein. For instance, to cite a few applications, the GPM framework can be used to help solve a weak formulation of some ordinary or partial differential equations, as well as some calculus of variations and optimal control problems. The solution u of the original problem (or an appropriate translate) is interpreted as a density with respect to the Lebesgue measure λ on some domain and one computes (or approximates) finitely many moments of the measure dµ := udλ by solving an appropriate finite-dimensional optimization problem. But then to recover an approximate solution of the original problem one has to solve an inverse problem from moments. This approach is particularly attractive when the data of the original problem consist of polynomials and basic semi-algebraic sets. In this case one may define a hierarchy (as the number of moments increases) of so-called semidefinite programs to compute approximations of increasing quality.

Contribution

In this paper we consider the following inverse problem from moments: Let µ be a finite Borel measure absolutely continuous with respect to some reference measure λ on a box Ω of R n and whose density u is assumed to be in L 2 (Ω, λ), with no continuity assumption as in previous works. The ultimate goal is to compute an approximation u d of u, based on the only knowledge of finitely many moments (say up to order d) of µ. In addition, for consistency, it would be highly desirable to also obtain some "convergence"

u d → u as d → ∞.
(a) Firstly, we approximate the density u by a polynomial u * d of degree d which minimizes the mean squared error Ω (u -p) 2 dλ (or equivalently the L 2 (Ω, λ)-norm u -p 2

2 ) over all polynomials p of degree at most d. We show that an unconstrained L 2 -norm minimizer u * d exists, is unique, and coincides with the simple polynomial approximation due to [START_REF] Teague | Image analysis via the general theory of moments[END_REF]; it can be determined by solving a system of linear equations. It turns out that u * d matches all moments up to degree d, and it is even easier to compute if it is expressed in the basis of polynomials that are orthonormal with respect to λ. No inversion is needed and the coefficients of u * d in such a basis are just the given moments. Moreover we show that u * d → u in L 2 (Ω, λ) as d → ∞, which is the best we can hope for in general since there is no continuity assumption on u; in particular u * d k → u almost-everywhere and almost-uniformly on Ω for some subsequence (d k ), k ∈ N. Even though both proofs are rather straightforward, to the best of our knowledge it has not been pointed out before that not only this mean squared error estimate u * d is much easier to compute than the corresponding maximum entropy estimate, but it also converges to u as d → ∞ in a much stronger sense. For the univariate case, in references [START_REF] Provost | Moment-Based Density Approximations[END_REF] and [START_REF] Athanassoulis | The truncated Hausdorff moment problem solved by using kernel density functions[END_REF] the authors address the problem of approximating a continuous density on a compact interval by polynomials or kernel density functions that match a fixed number of moments. In this case, convergence in supremum norm is obtained when the number of moments increases. An extension to the noncompact (Stieltjes) case is carried out in [START_REF] Gavriliadis | The truncated Stieltjes moment problem solved by using kernel density functions[END_REF]. Notice that in [START_REF] Provost | Moment-Based Density Approximations[END_REF] it was already observed that the resulting polynomial approximation also minimizes the mean square error and its coefficients solve a linear system of equations. In [START_REF] Bertero | Linear inverse problems with discrete data. I: General formulation and singular system analysis[END_REF][START_REF] Talenti | Recovering a function from a finite number of moments[END_REF] the minimum-norm solution (and not the minimum distance solution) is shown to be unique solution of a system of linear equations. In [START_REF] Jones | On extending the orthogonality property of minimum norm solutions in Hilbert space to general methods for linear inverse problems[END_REF] the minimal distance solution is considered but it is obtained as the solution of a constrained optimization problem and requires an initial guess for the density estimate.

(b) However, as already mentioned and unlike the maximum entropy estimate, the above unconstrained L 2 -norm minimizer u * d may not be a density as it may take negative values on Ω. Of course, the nonnegative function u * d+ := max[0, u * d ] also converges to u in L 2 but it is not a polynomial anymore. So we next propose to obtain a nonnegative polynomial approximation u * d by minimizing the same L 2 -norm criterion but now under the additional constraint that the candidate polynomial approximations should be nonnegative on Ω. In principle such a constraint is difficult to handle which probably explains why it has been ignored in previous works. Fortunately, if Ω is a compact basic semi-algebraic set one is able to enforce this positivity constraint by using Putinar's Positivstellensatz [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF] which provides a nice positivity certificate for polynomials strictly positive on Ω. Importantly, the resulting optimization problem is convex and even more, a semidefinite programming (SDP) problem which can be solved efficiently by public domain solvers based on interior-point algorithms. Moreover, again as in the unconstrained case, we prove the convergence u * d → u in L 2 (Ω, λ) as d → ∞ (and so almost-everywhere and almost-uniform convergence on Ω as well for some subsequence (u * d k ), k ∈ N) which is far stronger than the weak convergence obtained for the maximum entropy estimate. Notice, in [START_REF] Teague | Image analysis via the general theory of moments[END_REF][START_REF] Provost | Moment-Based Density Approximations[END_REF] methods for obtaining some non-negative estimates are discussed, however these estimates do not satisfy the same properties in terms of mean-square error minimization and convergence as in the unconstrained case. In the kernel density element method [START_REF] Athanassoulis | The truncated Hausdorff moment problem solved by using kernel density functions[END_REF][START_REF] Gavriliadis | The truncated Stieltjes moment problem solved by using kernel density functions[END_REF] a nonnegative density estimate for the univariate case is obtained by solving a constrained convex quadratic optimization problem. However, requiring each coefficient in the representation to be nonnegative as presented there seems more restrictive than the nonnegative polynomial approximation proposed in this paper.

(c) Our approach is illustrated on some challenging applications. In the first set of problems we are concerned with recovering the shape of geometrical objects whereas in the second set of problems we approximate solutions of nonlinear differential equations. Moreover, we demonstrate the potential of this approach for approximating densities with jump discontinuities, which is harder to achieve than for the smooth, univariate functions discussed in [START_REF] Provost | Moment-Based Density Approximations[END_REF]. The resulting L 2 -approximations clearly outperform the maximum entropy estimates with respect to both running time and pointwise approximation accuracy. Moreover, our approach is able to handle sets Ω more complicated than a box (as long as the moments of the measure udλ are available) as support for the unknown density, whereas such sets are a challenge for computing maximum entropy estimates because integrals of a combination of polynomials and exponentials of polynomials must be computed repeatedly.

Outline of the paper

In Section 2 we introduce the notation and we state the problem to be solved. In Section 3 we present our approach to approximate an unknown density u by a polynomial u * d of degree at most d via unconstrained and constrained L 2 -norm minimization, respectively; in both cases we also prove the convergence u * d → u in L 2 (Ω, λ) (and almost-uniform convergence on Ω as well for some subsequence) as d increases. In Section 4 we illustrate the approach on a number of examples -most notably from recovering geometric objects and approximating solutions of nonlinear differential equations -and highlight its advantages when compared with the maximum entropy estimation. Finally, we discuss methods to improve the stability of our approach by considering orthogonal bases for the functions spaces we use to approximate the density. And we discuss the limits of approximating discontinuous functions by smooth functions in connection with the well-known Gibbs effect.

2 Notation, definitions and preliminaries

Notation and definitions

Let R[x] (resp. R[x] d ) denote the ring of real polynomials in the variables x = (x 1 , . . . , x n ) (resp. polynomials of degree at most d), whereas Σ[x] (resp. Σ[x] d ) denotes its subset of sums of squares (SOS) polynomials (resp. SOS of degree at most 2d).

With Ω ⊂ R n and a given reference measure λ on Ω, let L 2 (Ω, λ) be the space of functions on Ω whose square is λ-integrable and let L 2 + (Ω, λ) ⊂ L 2 (Ω, λ) be the convex cone of nonnegative elements. Let C(Ω) (resp. C + (Ω)) be the space of continuous functions (resp. continuous nonnegative functions) on Ω. Let P (Ω) be the space of polynomials nonnegative on Ω.

For every α ∈ N n the notation x α stands for the monomial

x α 1 1 • • • x αn n and for every d ∈ N, let N n d := {α ∈ N n : j α j ≤ d} whose cardinal is s(d) = n+d d . A polynomial f ∈ R[x] is written x → f (x) = α∈N n f α x α
and f can be identified with its vector of coefficients f = (f α ) in the canonical basis (x α ), α ∈ N n . Denote by S n the space of real n × n symmetric matrices, and by S n + the cone of positive semidefinite elements of S n . For any A ∈ S n + the notation A 0 stands for positive semidefinite. A real sequence y = (y α ), α ∈ N n , has a representing measure if there exists some finite Borel measure µ on R n such that

y α = x α dµ(x), ∀ α ∈ N n .

Linear functional

Given a real sequence y = (y α ) define the Riesz linear functional L y : R[x] → R by:

f (= α f α x α ) → L y (f ) = α f α y α , f ∈ R[x].

Moment matrix

Given d ∈ N, the moment matrix of order d 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 d) .

entry (α, β) is y α+β , for every α, β ∈ N n d . If y has a representing measure µ then M d (y) 0 because f , M d (y)f = f 2 dµ ≥ 0, ∀ f ∈ R s(

Localizing matrix

With y as above and g ∈ R[x] (with g(x) = γ g γ x γ ), the localizing matrix of order d 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 γ g γ y (α+β+γ) , for every α, β ∈ N n d . If y has a representing measure µ whose support is contained in the set {x :

g(x) ≥ 0} then M d (g y) 0 because f , M d (g y)f = f 2 g dµ ≥ 0, ∀ f ∈ R s(d) .

Problem statement

We consider the following setting. For Ω ⊂ R n compact, let µ and λ be σ-finite Borel measures supported on Ω. Assume that the moments of λ are known and µ is absolutely continuous with respect to λ (µ λ) with Radon-Nikodým derivative (or density) u : Ω → R + , with respect to λ. The density u is unknown but we know finitely many moments y = (y α ) of µ, that is,

y α := Ω x α u(x)dλ(x) = Ω x α dµ(x), ∀α ∈ N n d , (1)
for some d ∈ N.

The issue is to find an estimate u d : Ω → R + for u, such that

Ω x α u d (x)dλ(x) = y α , ∀α ∈ N n d . (2)

Maximum entropy estimation

We briefly describe the maximum entropy method due to [START_REF] Jaynes | Information theory and statistical mechanics[END_REF][START_REF] Jaynes | Information theory and statistical mechanics II[END_REF][START_REF] Borwein | On the convergence of moment problems[END_REF] as a reference for later comparison with the mean squared error approach.

If one chooses the Boltzmann-Shannon entropy H(u) := -u log u, the resulting estimate with maximum-entropy is an optimal solution of the optimization problem max

u d Ω H(u d )dλ s.t. Ω x α u d (x)dλ(x) = y α , ∀α ∈ N n d .
It turns out that an optimal solution u * d is of the form d) . Hence, computing an optimal solution u * d reduces to solving the finite-dimensional convex optimization problem max

x → u * d (x) := exp   |α|≤d u α x α   for some vector u d = (u α ) ∈ R s(
u d ∈R s(d)    y, u d - Ω exp   |α|≤d u α x α   dλ(x)    (3)
where y = (y α ) is the given moment information on the unknown density u. If (u * d ), d ∈ N, is a sequence of optimal solutions to (3), then following weak convergence occurs:

lim d→∞ Ω ψ(x)u * d (x) dλ(x) = Ω ψ(x) u(x) dλ(x), (4)
for all bounded measurable functions ψ : Ω → R continuous almost everywhere. For more details the interested reader is referred to [START_REF] Borwein | On the convergence of moment problems[END_REF].

Since the estimate u * d is an exponential of a polynomial, it is guaranteed to be nonnegative on Ω and so it is a density. However, even though the problem is convex it remains hard to solve because in first or second-order optimization algorithms, computing the gradient or Hessian at a current iterate u d = (u α ) requires evaluating integrals of the form

Ω x α exp   |α|≤d u α x α   dλ(x), α ∈ N n d
which is a difficult task in general, except perhaps in small dimension n = 1 or 2.

The mean squared error approach

In this section we assume that the unknown density u is an element of L 2 (Ω, λ), and we introduce our mean squared error, or L 2 -norm approach, for density approximation.

Density approximation as an unconstrained problem

We now restrict u d to be a polynomial, i.e.

x → u d (x) := |α|≤d u α x α for some vector of coefficients u = (u α ) ∈ R s(d) . We first show how to obtain a polynomial estimate u * d ∈ R[x] d of u satisfying (2)
by solving an unconstrained optimization problem. Let z denote the sequence of moments of λ on Ω, i.e., z = (z α ), α ∈ N n , with

z α = Ω x α dλ, ∀ α ∈ N n ,
and let M d (z) denote the moment matrix of order d of λ. This matrix is easily computed since the moments of λ are known.

Consider the unconstrained optimization problem min

u d ∈R[x] d u -u d 2 2 = Ω (u -u d ) 2 dλ . (5) Proposition 1 Let Ω ⊂ R n have a nonempty interior and let λ(O) > 0 for some open set O ⊂ Ω. A polynomial u * d ∈ R[x]
d is an optimal solution of problem ( 5) if and only if its vector of coefficients u * d ∈ R s(d) is an optimal solution of the unconstrained quadratic optimization problem min

u d ∈R s(d) { u T d M d (z) u d -2 u T d y }. (6)
Then u * d := M -1 d (z)y is the unique solution of ( 6), and

u * d ∈ R[x] d satisfies: Ω x α u * d dλ = y α = Ω x α u dλ, ∀ α ∈ N n d . (7)
Proof: Observe that for every

u d ∈ R[x] d with vector of coefficients u d ∈ R s(d) , Ω (u -u d ) 2 dλ = Ω u 2 d dλ -2 Ω u d u dλ u d dµ + Ω u 2 dλ = u T d M d (z)u d -2 u T d y + Ω u 2 dλ.
The third term on the right handside being constant, it does not affect the optimization and can be ignored. Thus, the first claim follows.

The second claims follows from the well-known optimality conditions for unconstrained, convex quadratic programs and the fact that M d (z) is nonsingular because M d (z) 0 for all d ∈ N. Indeed, if q T M d (z)q = 0 for some 0 = q ∈ R s(d) then necessarily the polynomial q ∈ R[x] d with coefficient vector q vanishes on the open set O, which implies that q = 0, in contradiction with q = 0.

Finally, let e α ∈ R s(d) be the vector of coefficients associated with the monomial

x α , α ∈ N n d . from M d (z)u * d = y we deduce y α = e T α M d (z)u * d = Ω x α u * d dλ
which is the desired result.

Thus the polynomial u * d ∈ R[x] d minimizing the L 2 -norm distance to u coincides with the polynomial approximation due to [START_REF] Teague | Image analysis via the general theory of moments[END_REF] defined to be a polynomial which satisfies all conditions [START_REF] Athanassoulis | The truncated Hausdorff moment problem solved by using kernel density functions[END_REF]. Note that this is not the case anymore if one uses an L p -norm distance with p > 2.

Next, we obtain the following convergence result for the sequence of minimizers of problem [START_REF] Borwein | On the convergence of moment problems[END_REF]

Proof: Since Ω is compact, R[x] is dense in L 2 (Ω, λ). Hence as u ∈ L 2 (Ω, λ) there exists a sequence (v k ) ⊂ R[x], k ∈ N, with u -v k 2 → 0 as k → ∞. Observe that if d k = deg v k then as u * d k solves problem (5), it holds that u -v k 2 ≥ u -u * d k 2
for all k, which combined with u -v k 2 → 0 yields the desired result. The last statement follows from [1, Theorem 2.5.2 and 2.5.3].

Note that computing the L 2 norm minimizer u * d is equivalent to solving a system of linear equations, whereas computing the maximum entropy estimate requires solving the potentially hard convex optimization problem [START_REF] Baker | Padé approximants[END_REF]. Moreover, the L 2 -convergence u * du 2 → 0 in Proposition 2 (and so almost-everywhere and almost-uniform convergence for a subsequence) is much stronger than the weak convergence (4). On the other hand, unlike the maximum entropy estimate, the L 2 -approximation u * d ∈ R[x] d is not guaranteed to be nonnegative on Ω, hence it is not necessarily a density. Methods to overcome this shortcoming are discussed in §3.2.

Remark 1 In general the support K := supp µ of µ may be a strict subset of Ω = supp λ. In the case where K is not known or its geometry is complicated, one chooses a set Ω ⊃ K with a simple geometry so that moments of λ are computed easily. As demonstrated on numerical examples in Section 4, choosing an enclosing frame Ω as tight as possible is crucial for reducing the pointwise approximation error of u * d when the degree d is fixed. In the maximum entropy method of §2.3 no enclosing frame Ω ⊃ K is chosen. In the L 2 -approach, choosing Ω ⊃ K is a degree of freedom that sometimes can be exploited for a fine tuning of the approximation accuracy.

Remark 2 From the beginning we have considered a setting where an exact truncated moment vector y of the unknown density u ≥ 0 is given. However, usually one only has an approximate moment vector ỹ and in fact, it may even happen that ỹ is not the moment vector of a (nonnegative) measure. In the latter case, the maximum of the convex problem (3) is unbounded, whereas the L 2 -norm approach always yields a polynomial estimate

u * d ∈ R[x] d .
Remark 3 If ỹ is a slightly perturbed version of y, the resulting numerical error in u d and side effects caused by ill conditioning may be reduced by considering the regularized problem min

u d Ω (u -u d ) 2 dλ + u d 2 (8)
where u d 2 is the Euclidean norm of the coefficient vector

u d ∈ R s(d) of u d ∈ R[x] d
, and > 0 (fixed) is a regularization parameter approximately of the same order as the noise in ỹ. The coefficient vector of an optimal solution u * d ( ) of ( 8) is then given by

u * d ( ) = (M d + I) -1 ỹ. ( 9
)
The effect of small perturbations in y on the pointwise approximation accuracy of u * d for u is demonstrated on some numerical examples in Section 4. However, a more detailed analysis of the sensitivity of our approach for noise or errors in the given moment information is beyond the scope of this paper.

Density approximation as a constrained optimization problem

As we have just mentioned, the minimizer u * d of problem ( 5) is not guaranteed to yield a nonnegative approximation even if u ≥ 0 on Ω. As we next see, the function x → u * d+ (x) := max[0, u * d (x)] also converges to u for the L 2 -norm but

• it does not satisfy the moments constraints (i.e., does not match the moments of dµ = udλ up to order d);

• it is not a polynomial anymore (but a piecewise polynomial).

In the sequel, we address the second point by approximating the density with a polynomial nonnegative on Ω, which for practical purposes, is easier to manipulate than a piecewise polynomial. We do not address explicitly the first point, which is not as crucial in our opinion. Note however that at the price of increasing its degree and adding linear constraints, the resulting polynomial approximation may match an (a priori) fixed number of moments.

Adding a polynomial nonnegativity constraint to problem (5) yields the constrained optimization problem: min

u d ∈R[x] d { u -u d 2 2 : u d ≥ 0 on Ω} (10)
which, despite convexity, is untractable in general. We consider two alternative optimization problems to enforce nonnegativity of the approximation; the first one considers necessary conditions of positivity whereas the second one considers sufficient conditions for positivity.

Necessary conditions of positivity

Consider the optimization problem: min

u d ∈R[x] d { u -u d 2 2 : M d (u d z) 0 }, (11)
where z = (z α ) is the moment sequence of λ and M d (u d z) is the localizing matrix associated with u d and z. Observe that problem ( 11) is a convex optimization problem because the objective function is convex quadratic and the feasible set is defined by linear matrix inequalities (LMIs).

The rationale behind the semidefiniteness constraint M d (u d z) 0 in [START_REF] Izenman | Recent developments in nonparametric density estimation[END_REF] follows from a result in [START_REF] Lasserre | A new look at nonnegativity on closed sets and polynomial optimization[END_REF] which states that if supp λ = Ω and M d (u d z) ≥ 0 for all d then u d ≥ 0 on Ω.

Lemma 1 If Ω is compact then P (Ω) is dense in L 2 (Ω, λ) + with respect to the L 2 -norm.

Proof: As Ω is compact the polynomials are dense in L 2 (Ω, λ). Hence there exists a sequence (u

d ) ⊂ R[x], d ∈ N, such that u -u d 2 → 0 as d → ∞. But then the sequence (u + d ), d ∈ N, with u + d (x) := max[0, u d (x)
] for all x ∈ Ω, also converges for the L 2 -norm. Indeed,

Ω (u -u d ) 2 dλ = Ω∩{x:u d (x)<0} (u -u d ) 2 dλ + Ω∩{x:u d (x)≥0} (u -u d ) 2 dλ ≥ Ω∩{x:u d (x)<0} u 2 dλ + Ω∩{x:u d (x)≥0} (u -u d ) 2 dλ = Ω (u -u + d ) 2 dλ = u -u + d 2 2 . So let k ∈ N and d k ∈ N be such that u -u d k 2 < k -1 and so u -u + d k 2 < k -1 . As u + d k is continuous and Ω is compact, by the Stone-Weierstrass theorem there exists a sequence (v d k ) ⊂ R[x], ∈ N, that converges to u + d k for the supremum norm. Hence sup x∈Ω |u + d k -v d k | < k -1
for all ≥ k (for some k ). Therefore, the polynomial

w d k := v d k k + k -1 is positive on Ω and u -w d k k 2 ≤ u -u + d k 2 <k -1 + u + d k -w d k k 2 ≤ k -1 + u + d k -v d k k 2 <k -1 λ(Ω) 1/2 + v d k k -w d k k 2 =k -1 λ(Ω) 1/2 ≤ k -1 + 2k -1 λ(Ω) 1/2 .
Therefore we have found a sequence (w

d k k ) ⊂ P (Ω), k ∈ N, such that u -w d k k 2 → 0 as k → ∞.
Proposition 3 Let Ω be compact with nonempty interior. Then problem [START_REF] Izenman | Recent developments in nonparametric density estimation[END_REF] has an optimal solution u * d for every d ∈ N, and u -u * d 2 → 0 as d → ∞.

Proof: Fix d and consider a minimizing sequence (u ) ⊂ R[x] d with u -u 2 2 monotonically decreasing and converging to a given value as → ∞. We have u 2 ≤ u 2 + u -u 2 ≤ u 2 + u -u 0 2 for all ∈ N. Therefore as • 2 defines a norm on the finite dimensional space R[x] d (Ω has nonempty interior) the whole sequence (u ) is contained in the ball {v : v 2 ≤ u 2 + u -u 0 2 }. As the feasible set is closed, problem [START_REF] Izenman | Recent developments in nonparametric density estimation[END_REF] has an optimal solution.

Let (u * d ) ⊂ R[x]
be a sequence of optimal solutions of problem [START_REF] Izenman | Recent developments in nonparametric density estimation[END_REF]. By Lemma 1, P (Ω) is dense in L 2 (Ω) + . Thus there exists a sequence

(v k ) ⊂ R[x] ≥0 , k ∈ N, with u -v k 2 → 0 as k → ∞. As v k ≥ 0 on Ω then necessarily M d (v k z) 0 for all d and all k. In particular M d (v k d z) 0 where k d := max{k : deg v k ≤ d}. Therefore v k d ∈ R[x] d is a feasible solution of problem (11) which yields u -v k d 2 2 ≥ u -u *

Via sufficient conditions of positivity

Let b d (x) := (1, x 1 , . . . , x n , x 2 1 , x 1 x 2 , . . . , x d n ) T denote the standard monomial basis of R[x] d .
Let Ω be a basic compact semi-algebraic set defined by Ω = {x ∈ R n | g j (x) ≥ 0, j = 1, . . . , m} for some polynomials g j ∈ R[x], j = 1, . . . , m with d j = (deg g j )/2 . Then, with d ≥ max j d j , consider the optimization problem: min

u d ∈R[x] d u -u d 2 2 s.t. u d (x) = b d (x) T A 0 b d (x) + m j=1 b d-d j (x) T A j b d-d j (x)g j (x), A 0 ∈ S s(d) + , A j ∈ S s(d-d j ) + , j = 1, . . . , m. (12)
Since the equality constraints are linear in u d and the entries of A j , j = 0, . . . , m, the feasible set of ( 12) is a convex LMI set. Moreover, the objective function is convex quadratic.

Note that whereas the semidefinite constraint M d (u d z) 0 in [START_REF] Izenman | Recent developments in nonparametric density estimation[END_REF] was a relaxation of the nonnegativity constraint u d ≥ 0 on Ω, in (12) a feasible solution u d is necessarily nonnegative on Ω because the LMI constraint on u d is a Putinar certificate of positivity on Ω. However, in problem [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] we need to introduce m + 1 auxiliary matrix variables A j , whereas ( 11) is an optimization problem in the original coefficient vector u d and does not require such a lifting. Thus, problem [START_REF] Izenman | Recent developments in nonparametric density estimation[END_REF] is computationally easier to handle than problem [START_REF] Jaynes | Information theory and statistical mechanics[END_REF], although both are convex SDPs, which are substantially harder to solve than problem [START_REF] Donoho | Density estimation by wavelet thresholding[END_REF].

Proposition 4

Let Ω be compact with nonempty interior. For every d ≥ max j d j , problem ( 12) has an optimal solution u * d and u -u *

* d 1 ≤ u * d 2 ≤ u 2 + u -u * 0 2 . As u * d is feasible, u * d (x) = b d (x) T A d 0 b d (x) + m j=1 b d-d j (x) T A d j b d-d j (x)g i (x),
for some real symmetric matrices A d j 0, j = 0, . . . , m. Rewriting this as

u * d (x) = A d 0 , b d (x)b d (x) T + m j=1 A d j , b d-d j (x)b d-d j (x) T ,
and integrating with respect to λ yields

A d 0 , M d (z) + m j=1 A d j , M d-d j (g j z) = Ω u * d dλ ≤ u * d 1 ≤ a, with a := u * d 1 ≤ u 2 + u -u * 0 2 . Hence, for every d, A 0 , M d (z) ≤ a and A d j , M d-d j (g j z) ≤ a, j = 1, . . . , m.
As M d (z) 0, M d-d j (g j z) 0, and A d j 0, j = 1, . . . , m, we conclude that all matrices A d j are bounded. Therefore the minimizing sequence (u * d , (A d j )) belongs to a closed bounded set and as the mapping v → u -v 2 2 is continuous, an optimal solution exists.

From Lemma 1 there exists

(v k ) ⊂ R[x] ≥0 , k ∈ N, such that u -v k 2 2 → 0 as k → ∞. Using properties of norms, u -v k 2 -k -1 ≤ u -(v k + k -1 ) 2 ≤ u -v k 2 + k -1
, and so u -(v k + k -1 ) 2 → 0 as k → ∞. Moreover, as v k + k -1 is strictly positive on Ω, by Putinar's Positivstellensatz [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF], there exists d k such that

v k (x) = b d k (x) T A 0 b d k (x) + m j=1 b d k -d j (x) T A j b d k -d j (x)g j (x), ∀x,
for some real matrices A j 0, j = 0, . . . , m. So letting Remark 4 Since Ω is compact, Proposition 3 and 4 imply that minimizers of the two constrained L 2 norm minimization problems ( 11) and ( 12) converge almost uniformly to the unknown density u as in the unconstrained case.

d + k = max[deg v k , d k ], the polynomial v k + k -1 is a feasible solution of problem (12) whenever d ≥ d + k and with value u -(v k + k -1 ) 2 2 ≥ θ d + k . Hence u -v d + k 2 2 → 0 as k → ∞,
Remark 5 Our approach can handle quite general sets Ω and K as support and frame for the unknown measure in the unconstrained and constrained cases, the only restriction being that (i) Ω and K are basic compact semi-algebraic sets, and (ii) one can compute all moments of λ on Ω. In contrast, to solve problem (3) by local minimization algorithms using gradient and possibly Hessian information, integrals of the type

Ω x α exp β u β x β dλ(x)
must be evaluated. Such evaluations may be difficult as soon as n ≥ 3. In particular in higher dimensions, cubature formulas for approximating such integrals are difficult to obtain if Ω is not a box or a simplex.

Numerical experiments

In this section, we demonstrate the potential of our method on a range of examples. We measure the approximation error between a density u and its estimate u d by the average error

¯ d := Ω |u(x) -u d (x)|dλ(x)
and the maximum pointwise error

ˆ d := max x∈Ω |u(x) -u d (x)|.
In some examples we also consider ¯ o d and ˆ o d , the respective errors on particular segments of the interior of Ω. We compare the performance of our approach to the maximum entropy estimation from §2.3. Both methods are encoded in Matlab. We implemented the mean squared error (L 2 ) minimization approach for the standard monomial basis, which results in solving a linear system in the unconstrained case. In the constrained case we apply SeDuMi to solve the resulting SDP problem. In both cases, the numerical stability can be improved by using an orthogonal basis (such as Legendre or Chebychev polynomials) for R[x]. In order to solve the unconstrained, concave optimization problem (3), we apply the Matlab Optimization Toolbox command fminunc as a black-box solver. The observed performance of the maximum entropy estimation (MEE) method may be improved when applying more specialized software. Example 1 First we consider the problem of retrieving the density u(x) = x 1 + x 2 on [0, 1] 2 given its moments

y α = Ω x α u(x)d x = 1 (α 1 +1)(α 2 +2) + 1 (α 1 +2)(α 2 +1
) , which was considered as a test case in [START_REF] Mnatsakanov | Moment-recovered approximations of multivariate distributions: the Laplace transform inversion[END_REF]. This example is a priori very favorable for our technique since the desired u is a polynomial itself. When solving problem (5) for d ∈ {3, 5, 10} we do obtain the correct solution u * d = (0, 1, 1, 0, . . .) T in less than 0.1 secs. Thus, the pointwise approximation is much better than the one achieved in [START_REF] Mnatsakanov | Moment-recovered approximations of multivariate distributions: the Laplace transform inversion[END_REF]. Moreover, u * d ≥ 0 without adding LMI constraints. The question arises of how the polynomial approximation behaves if the moment vector contains some noise, i.e. if it does not exactly coincide with the moment vector of the desired u. We solve problem (5) for ỹ := y + where the maximal relative componentwise perturbation between y and ỹ is less than 3%. The pointwise error between u * 10 (y) and the solution for the perturbed moment vector u * 10 (ỹ) is sufficiently small, as pictured in Figure 1.

Example 2 Next we consider recovering the function u : [-1, 1] → R with u(x) = |x| as a first example of a nondifferentiable function. In a first step we solve problem (5) for the exact moment vector with entries y k := In a second step we consider a perturbed moment vector ỹ as input. When solving problem (5), we observe that both errors ˆ d , ¯ d on the interior of [-1, 1] -in particular at the nondifferentiable point x = 0 -decrease for increasing d, whereas the pointwise error increases at the boundary of the domain. Although providing good approximations for u on the entire interior of [-1, 1], the estimates u * 30 and u * 50 take negative values at the boundary {-1, 1}. This is circumvented by solving problem [START_REF] Jaynes | Information theory and statistical mechanics[END_REF]. The new estimates are globally nonnegative while their approximation accuracy is only slightly worse than in the unconstrained case, as reported in Table 2 and pictured in Figure 2 (right).

1 -1 |x|x k dx = 1+(-1) k k+2 , k = 0, . . . ,
Example 3 Consider the functions u 1 , u 2 : [-1, 1] → R, u 1 (x) = |x| 1 2 , u 2 (x) = | 1 4 -x 2 | 1 2
, which are both not locally Lipschitz. Applying mean squared error minimization for d ∈ {20, 50} to the exact moment vectors yields accurate polynomial approximations for both functions on their entire domains, even at the boundary and at the points where the functions are not locally Lipschitz, cf. Figure 3.

Recovering geometric objects

One of the main applications of density estimation from moments is the shape reconstruction of geometric objects in image analysis. There has been extensive research on this topic, c.f. [START_REF] Teague | Image analysis via the general theory of moments[END_REF][START_REF] Liao | On image analysis by moments[END_REF] and the references therein. The reconstruction of geometric objects is a particular case of density estimation when u = I K , i.e. the desired density u is the indicator function of the geometric object K ⊂ Ω ⊂ R n . Its given moments y α = Ω x α I K (x)dλ(x) = K x α dλ(x) do not depend on the frame Ω. However, Ω does enter when computing the matrix M d (z) in problem [START_REF] Donoho | Density estimation by wavelet thresholding[END_REF]. As indicated in Remark 1 and demonstrated below, the choice of the enclosing frame Ω for K is crucial for the pointwise approximation accuracy for a fixed order d. Since u has a special structure, we derive an estimate K d for K by choosing a superlevel set of u * d as proposed in [START_REF] Teague | Image analysis via the general theory of moments[END_REF]:

K d := {x ∈ R n | u * d (x) ≥ 1 2 }.
Example 4 A first object we recover is shaped like the letter E as in [START_REF] Teague | Image analysis via the general theory of moments[END_REF]. We determine the mean squared error minimizer for the same moment vector y with d ∈ {3, 5, 8, 10}, but for three different frames Ω. As pictured in Figure 4, we are able to reconstruct the estimates K d derived in [START_REF] Teague | Image analysis via the general theory of moments[END_REF] when Ω is chosen tight. Moreover, we observe that this good approximation of the severely nonconvex set K and its discontinuous indicator function for a small number of moments d depends heavily on the choice for Ω. The wider the frame, the worse gets the approximation accuracy of the truncated estimate. Applying maximum entropy estimation for d ∈ {3, 5, 8} yields density estimates of comparable accuracy than mean squared error minimization, cf. Figure 5. However, the computational effort is much larger: for d ∈ {3, 5, 8}, the computational times for maximum entropy estimation are 82, 902 and 4108 seconds, respectively, whereas the mean squared error minimizer can be determined in less than one second for these values of d.

Example 5 Secondly we consider recovering an F-shaped set, a less symmetric example than the E-shaped set of Example 4. As previously, we observe that the approximation accuracy of K d relies heavily on the frame Ω. Even though the set K has a complicated geometry, K 10 approximates K accurately if Ω is chosen sufficiently tight, cf. Figure 6.

Example 6 Consider approximating

K := {x ∈ R 2 | x 1 (x 2 1 -3x 2 2 ) + (x 2 1 + x 2
2 ) 2 ≥ 0}, a nonconvex region enclosed by a trefoil curve. Since this curve is of genus zero, its moment vector y can be determined exactly. Again, we need to choose an appropriate frame Ω ⊃ K. The results for Ω = B(0, 1) and d ∈ {3, 5, 8, 10} are pictured in Figure 7.

Approximating solutions of differential equations

Another important class of moment problems arises from the numerical analysis of ordinary and partial differential equations. Solutions of certain nonlinear differential equations can be understood as densities of measures associated with the particular differential equation. We obtain an approximate solution from the moment vector of this measure by solving an SDP problem. Approaches for deriving moment vectors associated with solutions of nonlinear differential equations have been introduced in [START_REF] Mevissen | Moment and SDP relaxation techniques for smooth approximations of problems involving nonlinear differential equation[END_REF][START_REF] Henrion | The generalized problem of moments for nonlinear differential equations[END_REF] and are omitted here. We assume that the moment vector of a measure whose density is a solution of the respective differential equation is given in the following examples.

Example 7 Given the moment vectors of a solution (u, v) of the following reaction- diffusion equation [START_REF] Mimura | Asymptotic behaviors of a parabolic system related to a planktonic prey and predator model[END_REF]:

1 20 u + 1 9 (35 + 16u -u 2 ) u -u v = 0, 4v -1 + 2 5 v v + u v = 0, u (0) = u (5) = v (0) = v (5) = 0, 0 ≤ u, v ≤ 14 on [0, 5],
we apply both maximum entropy estimation and mean squared error minimization to approximate the desired solution. For the numerical results, see Table 3 and Figure 8. We observe that the mean squared error minimizers provide accurate pointwise approximations for (u, v) on the entire domain, whereas the maximum entropy estimates provides a fairly accurate pointwise approximation on some segment of the domain only. Moreover, the mean squared error minimizers are obtained extremely fast as solutions of linear systems compared with the maximum entropy estimates. Thus, in this example, mean squared error minimization is clearly superior to maximum entropy when estimating densities from moments.

Example 8 Given the moment vector of the nontrivial, positive solution of the Allen-Cahn bifurcation PDE: we apply both maximum entropy and mean squared error minization. The numerical results are reported in Table 4. The approximation accuracy for both methods is comparable, with mean squared error minimization being slightly more precise. However, applying maximum entropy estimation is limited for this problem as the cases d > 6 are numerically too heavy to be solved in reasonable time. Mean squared error minimization yields increasingly better estimates for the desired solution within seconds, as pictured in Figure 9.

u xx + u yy + 22u(1 -u 2 ) = 0 on [0, 1] 2 , u = 0 on ∂[0, 1] 2 , 0 ≤ u ≤ 1 on [0, 1] 2 , (13)
Example 9 Given the moment vector of the viscosity solution of the classical eikonal

d method ¯ d ¯ o d time (sec.) 3 MEE 5
(x) = I [0.5,1] (x). PDE: u 2 x + u 2 y -1 = 0 on [0, 1] 2 , u = 0 on ∂[0, 1] 2 , (14)
we apply both maximum entropy and mean squared error minization. Both methods provide better pointwise approximates for increasing degrees d, cf. Table 5 and Figure 10. Again, the advantage of using mean squared error minimization is its negligible computation time.

Approximating indicator functions

In all examples discussed so far, the polynomial approximates u * d of solution u obtained by solving the unconstrained problem (5) have been nonnegative on Ω. It was therefore not necessary to solve the constrained problems [START_REF] Izenman | Recent developments in nonparametric density estimation[END_REF] or [START_REF] Jaynes | Information theory and statistical mechanics[END_REF].

Example 10 As a first simple example where solving one of the constrained problems is required, consider u : [0, 1] → R, u(x) = I [0.5,1] (x) and its given vector of moments. We solve problems ( 5) and ( 12) for d ∈ {10, 50, 100}. As illustrated in Figure 11, solving problem [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] provides globally nonnegative estimates for u. This comes at the price of decreasing approximation accuracy compared to solving the unconstrained problem, since the feasible set of problem ( 12) is a subset of the feasible set of (5), cf. Table 6.

Orthogonal bases and the Gibbs effect

Density estimation via L 2 norm minimization over the truncated space of polynomials is a special case of approximating a function belonging to a complicated functional space by a linear combination of basis elements of an easier, well-understood functional space. This general setting has been considered for both real and complex functional spaces. In the following we discuss the particular case when the basis of the easier functional space is orthogonal, and also effects resulting from truncating the infinite dimensional series expansion of the unknown function.

The complex case

The classic, complex analogue of the real, inverse moment problem discussed in the previous sections is the problem of approximating a periodic function u by a trigonometric polynomial. An orthonormal basis for the space of trigonometric polynomials is given by (e -ikx ) k∈Z . It provides the Fourier series expansion: The Fourier approximation for a periodic function is therefore the trigonometric analogue of the real polynomial approximation for a function obtained by mean squared error minimization in our approach. As in the real case we have almost uniform convergence. However, u d does not converge to u uniformly if u piecewise continuously differentiable with jump discountinuities, due to an effect known as the Gibbs phenomenon. It states that the truncated approximation shows a near constant overshoot and undershoot near a jump discontinuity. This overshoot does not vanish, it only moves closer to the jump for increasing d.

Since the truncated Fourier series is the trigonometric analogue of unconstrained mean squared norm minimizer of problem [START_REF] Borwein | On the convergence of moment problems[END_REF], the question arises of how the trigonometric estimate behaves for a periodic functions when adding nonnegativity constraints as in problems [START_REF] Izenman | Recent developments in nonparametric density estimation[END_REF] and ( 12) which aim at preventing the estimate from over-or undershooting near jump discontinuities.

In order to derive a tractable SDP problem, the nonnegativity constraints for the trigonometric polynomial need to be relaxed or tightened to LMI constraints. The difference with the real case will be the moment and localizing matrices being of Toeplitz type in contrast with the Hankel type matrices in the constraints of problem [START_REF] Izenman | Recent developments in nonparametric density estimation[END_REF].

The real case

In our discussions in the previous sections we approximated an unknown density by a linear combination of elements of the monomial basis of the space of real polynomials. This choice of a basis for R[x] has several theoretical and practical shortcomings. For once, it is not an orthogonal basis with respect to the Lebesgue measure, and moreover the moment matrix M d (z) in problem [START_REF] Borwein | On the convergence of moment problems[END_REF], whose inverse needs to be computed to determine u * d , is severely ill-conditioned [START_REF] Bertero | Linear inverse problems with discrete data. I: General formulation and singular system analysis[END_REF][START_REF] Teague | Image analysis via the general theory of moments[END_REF]. As already pointed out in [START_REF] Teague | Image analysis via the general theory of moments[END_REF], when choosing an orthogonal basis such as Legendre or Chebychev polynomials for R[x], the L 2 norm minimizer u * d can be determined according to a closed-form formula. Thus, we do not even have to solve a linear system of equations. This is essentially the real analogue of the closed form formula for the coefficients in the Fourier series approximation of periodic functions.

In §5.1 we discussed the Gibbs phenomenon. For reasons outlined there, we expect to observe this effect in the real case as well, when approximating an L 2 integrable function with jump discontinuities. In fact, the function from Example 10 illustrates that.

As shown in Figure 11, we observe an over-and undershoot on both sides of the discontinuity. The amplitude of this overshoot does not decrease for increasing d, but the overshoot moves closer to the jump. When adding the global nonnegativity constraint for the polynomial estimate, the undershoot at the left side disappears. However, it is compensated for by a weaker, overall pointwise approximation accuracy of the estimate. This observation is a first, partial answer to the question raised in the complex case.

Conclusion

We introduced an approach for estimating the density of a measure given a finite number of its moments only, in the multivarite case and with no continuity assumption on the density. As an estimate we choose the polynomial minimizing the L 2 norm distance, or mean squared error to the unknown density. We have shown that this estimate is easy to determine by solving a linear system of equations. Moreover, it converges almost uniformly towards the desired density when the degree increases. By minimizing the mean squared error subject to additional linear matrix inequality constraints, which translates to solving a semidefinite program, we obtain a density estimate guaranteed to be nonnegative on the support of the measure. Also in the constrained case, we have shown almost uniform convergence of the nonnegative mean squared error minimizer towards the unknown density, when the degree increases. Mean squared error minimization is often superior to maximum entropy estimation in terms of approximation accuracy and computation time, as demonstrated for a number of examples. Moreover, it is able to handle general, basic, compact semialgebraic sets as support for the unknown measure, which present a challenge to maximum entropy estimation in higher dimension.

  , d ∈ N. Proposition 2 Let Ω be compact with nonempty interior and let λ be finite with λ(O) > 0 for some open set O ⊂ Ω. Let (u * d ), d ∈ N, be the sequence of minimizers of problem (5). Then u -u * d 2 → 0 as d → ∞. In particular there is a subsequence (d k ), k ∈ N, such that u * d k → u, λ-almost everywhere and λ-almost uniformly on Ω, as k → ∞.

d 2 2

 2 for all d. Combining with u -v k d 2 → 0 yields the desired result.

d 2 2

 2 → 0 as d → ∞. Proof: With d ≥ max j d j fixed, consider a minimizing sequence (u * d ) ⊂ R[x] d for problem (12). As in the proof of Proposition 3 one has u

2 2

 2 and by monotonicity of the sequence ( u -u * d ), d ∈ N, the result follows.

Figure 1 :

 1 Figure 1: Degree 10 approximation of the density from exact moment vector (left) vs. perturbed moment vector (right).

Figure 2 :

 2 Figure 2: Mean squared error minimizers for u(x) = |x| from exact (left) and perturbed (right) moments. Blue: u, green: u 20 , red: u 30 , cyan: u 50 .

Figure 3 :

 3 Figure 3: Mean squared error minimizers for u 1 (left) and u 2 (right) from exact moments. Blue: u, green: u 20 , cyan: u 50 .

Figure 4 :

 4 Figure 4: Recovering the letter E with mean squared error minimization; tight frames (top) vs. loose frames (bottom); original (left) vs. estimates for degrees d ∈ {3, 5, 8, 10} (from left to right).

Figure 5 :

 5 Figure 5: Recovering the letter E with maximum entropy estimation; original (left) vs. estimates for degrees d ∈ {3, 5, 8} (from left to right).

Figure 6 :

 6 Figure 6: Recovering the letter F with mean squared error minimization; tight frames (top) vs. loose frames (bottom); original (left) vs. estimates for degrees d ∈ {3, 5, 8, 10} (from left to right)..

Figure 7 :

 7 Figure 7: Trefoil region K (left) and ((K d \ K) ∪ (K \ K d )) for d ∈ {3, 5, 8, 10} (from left to right).

Figure 8 :

 8 Figure 8: Estimates for u (left) and v (right) of the reaction-diffusion equation via maximum entropy estimation (above) and mean squared error minization (below). Blue: u, green: u * 10 , red: u * 20 , cyan: u * 30 , magenta: u * 50 .

Figure 9 :

 9 Figure 9: Estimates for solution of the Allen-Cahn bifuraction PDE by maximum entropy estimation (above) and mean squared error minization (below). Exact solution (left) vs estimates for d ∈ {4, 6} and d ∈ {4, 6, 10, 12}, respectively.

Figure 10 :

 10 Figure 10: Estimates for solution (left) of eikonal PDE by maximum entropy estimation (center) for d = 6 and mean squared error minimization (right) for d = 10.

Figure 11 :

 11 Figure 11: Unconstrained and constrained mean squared error minimization estimates. Blue: u, green: u * 10 unconstrained, red: u * 10 constrained, cyan: u * 100 unconstrained, magenta: u * 100 constrained.

u(x) = k∈Z y k e ikx , where y k = 1 2ππ

 1 -π u(x)e -ikx dx in the univariate case, or u(x, y) = k,l∈Z y k,l e ikx e ily , where y k,l = (x, y)e -ikx e -ily dx dy in the bivariate case. One obtains a trigonometric polynomial approximation for u when truncating this series expansion at some d ∈ N,u d (x) = d k=-d y k e ikx .It is a well-known fact in Fourier theory that 2π 0 (u -u d ) 2 dx → 0 for d → ∞.

Table 1 :

 1 d corresponding to the density u. The resulting estimates u * d for d ∈ {20, 30, 50} provide a highly accurate pointwise approximation of u on the entire domain, as reported in Table1and pictured in Figure2(left). Estimating the density u(x) = |x| from an exact moment vector.

		d problem	¯ d	ˆ d
		20	(5)	0.0031 0.0296
		30	(5)	0.0022 0.0265
		50	(5)	0.0021 0.0251
	d problem	¯ d	ˆ d	¯ o d	ˆ o d
	10	(5)	0.0252 0.2810 0.0207 0.1358
	20	(5)	0.0244 0.7934 0.0142 0.0952
	30	(5)	0.0237 1.0956 0.0112 0.1024
	30	(12)	0.0176 0.4705 0.0106 0.0937
	50	(5)	0.0236 1.4591 0.0088 0.1023
	50	(12)	0.0206 0.6632 0.0118 0.0979

Table 2 :

 2 Estimating the density u(x) = |x| from a perturbed moment vector.

Table 3 :

 3 Approximation accuracy of maximum entropy estimation (MEE) and mean squared error (L 2 ) optima for the reaction-diffusion equation.

	d function method ¯ d (M )	ˆ d	¯ o d	ˆ o d	time (sec.)
	10	u	MEE	2.0536 4.1063 2.1679 4.1063	38
	30	u	MEE	1.7066 3.4775 1.7974 3.4775	2
	50	u	MEE	1.6978 3.3858 1.7792 3.3858	820
	10	u	L 2		0.4611 1.4506 0.4619 1.4506	0.1
	30	u	L 2		0.3942 1.2431 0.4101 1.2431	0.1
	50	u	L 2		0.3941 1.2173 0.4004 1.2173	0.1
	10	v	MEE	0.5782 2.0982 0.5205 1.0806	221
	30	v	MEE	0.5978 9.6060 0.4592 1.2356	645
	50	v	MEE	0.6853 23.5993 0.4117 1.2973	3306
	10	v	L 2		0.3429 5.4767 0.1765 0.5617	0.1
	30	v	L 2		0.3286 12.4501 0.1024 0.6253	0.1
	50	v	L 2		0.3744 14.2223 0.1454 0.5907	0.1
			d method	¯ d	¯ o d	time (sec.)
			4	MEE	1.3e-2 1.1e-2	566
			6	MEE	9.7e-3 8.9e-3	2489
			4	L 2		1.8e-3 1.4e-3	0.1
			6	L 2		4.6e-4 2.8e-4	0.1
			10	L 2		5.5e-4 2.1e-4	0.1
			12	L 2		5.7e-4 2.1e-4	0.1

Table 4 :

 4 Approximation accuracy of maximum entropy estimation (MEE) and mean squared error (L 2 ) optima for the Allen-Cahn bifurcation PDE.

Table 5 :

 5 Approximation accuracy of maximum entropy estimation (MEE) and mean squared error (L 2 ) minization for the eikonal PDE.

			.4e-2 5.1e-2	20
	6	MEE	1.9e-2 1.9e-2	2192
	3	L 2	2.8e-2 2.6e-2	0.1
	6	L 2	2.8e-2 2.6e-2	0.1
	10	L 2	1.4e-2 1.4e-2	0.2

Table 6 :

 6 Mean squared error minimization for u

	d problem ¯ d	ˆ d	time (sec.)
	10	(5)	0.08 0.50	0.1
	10	(12)	0.11 0.56	0.9
	50	(5)	0.05 0.50	0.1
	50	(12)	0.08 0.54	2.7
	100	(5)	0.05 0.50	0.3
	100	(12)	0.07 0.55	37.5

* D. Henrion acknowledges support by project number 103/10/0628 of the Grant Agency of the Czech Republic. The major part of this work was carried out during M. Mevissen's stay at LAAS-CNRS, supported by a fellowship within the Postdoctoral Programme of the German Academic Exchange Service.


Mohamed Amine Kafi
email: kafi@mail.cerist.dz

Yacine Challal

Djamel Djenouri

Abdelmadjid Bouabdallah

Khelladinadjib Lyes

Badache

A study of Wireless Sensor Network Architectures and Projects for Traffic Light Monitoring

Keywords: WSN, ITS, traffic management

Vehicular traffic is increasing around the world, especially in urban areas. This increase results in a huge traffic congestion, which has dramatic consequences on economy, human health, and environment. Traditional methods used for traffic management, surveillance and control become inefficient in terms of performance, cost, maintenance, and support, with the increased traffic. Wireless Sensor Networks (WSN) is an emergent technology with an effective potential to overcome these difficulties, and will have a great added value to intelligent transportation systems (ITS) overall. In this survey, we review traffic light projects and solutions. We discuss their architectural and engineering challenges, and shed some light on the future trends as well.

Introduction

Urban and metropolitan areas worldwide witness tremendous increase of vehicular traffic. This increase results in huge traffic congestion that induces dramatic consequences on economy, human health, and environment. The classical methods for traffic management, surveillance and control become inefficient in terms of performance, cost, maintenance, and support.

Conventional Intelligent transportation systems (ITS) detect vehicles in predefined positions, and use wired technologies (for co mmunication and power supply), which increases their construction and maintenance cost. Different vehicles detection techniques are used by conventional ITS such as aboveground sensors (video image processing, micro wave radar, laser radar, passive infrared, ultrasonic, and passive acoustic array) but their accuracy is environ ment dependent. Other detection techniques use underground sensors (inductive loops, micro-loop probes, and pneumatic road tubes). These techniques cause disruption of traffic upon installation and repair. These detection techniques are bulky, power-hungry, and expensive to install, maintain and repair. These characteristics subvert the scalability of ITS and affect their majo r objectives [START_REF] Khanafer | WSN Architectures for Intelligent Transportation Systems[END_REF].

Develop ments in embedded Systems and wireless technology have enabled the creation of small sized and cheap wireless sensors. Wireless sensor networks (WSNs) are co mposed of spatially distributed nodes that communicate and sense the surrounding physical world. Each node contains sensors, a processor, a memory, a radio, and an energy source. WSN have many advantages in terms of wireless enabled communication, low cost and flexibility. This technology has a great potential to overcome existing difficulties of ITS.

With WSN, different types of motes can be used to sense, process and send sensed data in purpose of mastering co mplex situations in an optimized manner and enabling online adaptive traffic control systems. Data of interest include position, traffic condition, local weather, images, acceleration, etc. Some possible sensors for ITS include magneto resistive, light, pressure, video, etc.

In [START_REF] Lee | A survey of urban vehicular sensing platforms[END_REF], a survey of only vehicular sensor networks (VSN) p latforms was presented, while [START_REF] Sharma | Applications of Wireless Sensor Network in Intelligent Traffic System: A Review[END_REF] discusses many u rban applications. In our survey we review WSN based (on road or on vehicles) traffic light projects and solutions. We discuss their architectures and highlight some open issues.

The remainder of this paper will be organized as fo llo ws. In section 2, we present some applicat ions of ITS. In section 3, we d iscuss some requirements of WSN usage in ITS systems. We present in section 4 relevant architectures that are used in ITS t raffic light projects based on WSN technology. In section 5, we present and analyze some projects of traffic light management relying on WSN. In section 6 a b rief discussion is presented. Finally, we end up this paper with co nclusions and open issues in section 7.

ITS applications

Intelligent Transportation Systems (ITS) aim to imp rove traffic quality by enhancing safety and reducing traveling time and fuel consumption through gathering, organizing, analy zing, using, and sharing traffic informat ion.

In tradit ional ITS heavy and expensive detection systems are used to collect traffic data. In [START_REF] Elena | A Summary of Vehicle Detection and Surveillance Technologies used in Intelligent Transportation Systems, the Vehicle Detector Clearinghouse[END_REF], a summary of trad itional Vehicle Detection and Surveillance Technologies is presented. Co mmunicating this informat ion to decision maker system is done in wired costly manner. These two drawbacks are deleted by the evolution of embedded systems and wireless technologies which give birth to WSN networks based ITS characterized by their miniature size, ease of installation and certainly their wireless communicat ion that offer scalability with less cost. ITS applications are nu merous and can be categorized depending on many factors. [5] Authors' present a large variety of ITS possible applicat ions in both developed and developing countries . Follo wing a possible classification is presented. Drivers' safety by transmitting accidents and weather information and dispatch ambulances and fire trucks.

Traffic Management to minimize congestion and optimize road capacity utilizat ion. Th is category includes traffic optimization by applying automated traffic assistance decisions and real-time traffic light handling in both isolated and interconnected intersections. Smart cities. Th is category includes: navigation orientations before and during a trip to minimize cost, pollution prevention, parking locations finding and optimization, road maintenance and repair, public transport management, electronic t ransactions to allow fees payment on roads and parking.

Requirements and challenges

We distinguish five main fundamental requirements that should be fulfilled by any enabling technology: reliability, security, interoperability, real-time co mmunication, and mu ltimodal sensing.

Data reliability is of a high importance for urban traffic monitoring and management, but this becomes challenging when using WSN due to the lo w coverage and limited capacity of sensor motes, added to the lossy nature of wireless commun ication channels. Therefore, ensuring reliab le end -to-end co mmunication must be taken into account by any proposed WSN architecture and co mmunication protocols. Many applications need real-time information collection to take in-time decisions. Traffic light monitoring is a typical examp le, where dynamic crossroad management depends on real-time traffic distribution.

The rapid develop ment of wireless technologies and different types of sensors imposes on the proposed solutions to take into account hardware heterogeneity to ensure interoperability and scalability.

Co mmunicat ion security must be handled efficiently by the proposed solutions. Indeed, wireless nature of communicat ions in WSN is more prone to security attacks that wired technologies. Wireless communicat ion can be easily eavesdropped or denied through jamming. This constitutes a priv ileged target for hackers and terroris m with potentially high collateral damage. Moreover, motes may be deployed without physical protection or surveillance. This increases the risk of co mpro mising the motes by intruders and get access to sensitive data or credentials to pursue attacks through privilege escalation. Motes may also store sensitive informat ion relat ing to users privacy. Therefore, traffic management solutions based on WSN must take into consideration the vulnerable nature of WSN and their components. Specifically, the commun ication protocols must be secure enough to prevent attacks and sufficiently reliable and fau lt tolerant to guarantee service continuity in case of failures.

Traffic control is a mu lti-object ive optimizat ion problem that requires the optimization of various weighted variables as: gas emissions , number of stops, noise levels due to traffic, etc.This diversity in optimization variab les requires a diversity of sensors that is actually one of the fundamental design objectives of WSN. We believe that the mult imodal nature of WSN is a strong point th at allows feeding the traffic control mult i-objective problem with the required informat ion in real-time to take the best traffic management decisions on time.

These challenges should be overcome by any WSN based traffic management solution. So me requirements are already fulfilled by WSN technology which make it a favorite candidate for effective ITS systems.

WSNs architectures for Traffic light management

Network architecture for ITS applications using WSN technology changes from an applicat ion to anoth er, depending on the needs and the cost. Information exchange can be performed either through ad -hoc communication, or using infrastructure, or hybrid. We also distinguish two types of sensors: on -road sensors and on-vehicle sensors. The combination of sensor types and communication paradig ms g ives birth to various wireless sensor network architectures for ITS applicat ions as illustrated in Figure 1. Decisions can be taken even locally or cooperatively. In Figure 2.b, a co mmunicat ion example is presented. (iii) Hybrid ad-hoc sensor network, wh ich is more robust and combines the two previous deployments. Both on -road and in-vehicle sensors, ( Figure 2.C) exchange traffic information to cooperatively take correct and real-t ime decisions.

Hybrid ad-hoc and infrastructure-enabled monitoring:

In this paradig m, static on-road and mobile onvehicle sensors and devices commun icate using available infrastructures (wifi, W imax, cellular, BS) o r using multi-hop ad-hoc communications in the absence of infrastructure as shown in Figure 4. Th is kind of architecture is the most efficient and takes advantage of all the availab le architectures.

Traffic Light Control Projects

The goal of t raffic light signals in an intersection is to schedule conflicting traffic to move in a secure and manner and exploit ing the road capacity as possible. For managing a traffic light controlled intersection, there are many parameters to be defined [START_REF] Papageorgiou | Review of road traffic control strategies[END_REF], following the most important ones.

Signal cycle is the repetition of the signal co mbinations and its duration is a "cycle t ime". A stage (or phase) is a part of the signal cycle, during which one set of streams can move securely. Constant lost times (a few seconds) between stages are necessary to avoid streams interference. Split is the green duration of each stage that should be optimized according to the demands. Offset is the time difference between green phases of successive intersections that may give rise to a "green wave", where traffic flows over several intersections without need to stop on mult iple intersections.

Traffic light control system can manage one intersection or a network of intersections giving so birth to: isolated strategies and coordinated strategies. Also, the system control can use "fixed-t ime" strategies derived off-line based on historical demands. It uses optimizat ion codes. The system can be "trafficresponsive" using real-time measurements. Optimizing isolated intersection behavior using local traffic informat ion achieves locally optimal performance which often provides a subopt imal system performance. To achieve g lobal optimal perfo rmance, neighbouring intersections must exchange their traffic informat ion. TRANSYT [START_REF] Robertson | TRANSYT method for area traffic control[END_REF] is the most fixed-t ime coordinated strategy used. It is also used as a reference to test real-time methods' improvements. Optimization algorith ms introduce small changes to the control parameters (cycle time, split …) to enhance system performance. Many optimizing algorithms were p roposed, including fuzzy logic algorith ms [START_REF] Lee | Fuzzy controller for intersection group[END_REF], neural network -fu zzy [START_REF] Srinivasan | Neural networks for real-time traffic signal control[END_REF], mu lti-objective genetic algorith ms [START_REF] Chen | A Dynamic Optimization Method for Traffic Signal Timings Based on Genetic Algorithm[END_REF][START_REF] Chen | Real-coded genetic algorithm for signal timing optimization of a single intersection[END_REF][START_REF] Nakamura | Genetic algorithm approach for adaptive offset optimization for the fluctuation of traffic flow[END_REF], Markov Process [START_REF] Li | The model analysis of vehicles situation and distribution in intersections based on Markov process[END_REF], etc.

Many parameters for quantifying the effectiveness of such algorithms were introduced such as the control delay or average wait ing time, wh ich is the delay that a vehicle spends due to the presence o f the signal control, the maximu m intersection throughput which is the nu mber of vehicles passing through the intersection, etc. Minimizing the delay at intersections can be achieved by selecting a shortest cycle length to produce less "red" signal time length and shorter queues but without influencing the delay by phase changing overhead.

Several adaptive traffic control systems have been implemented. The most important ones are Split, Cycle and Offset Optimization Technique (SCOOT) [START_REF] Robertson | Optimizing Networks of Traffic Signals in Real T ime-The SCOOT Method[END_REF] and Sydney Coordinated Adaptive Traffic System (SCATS) [START_REF] Sims | The Sydney Coordinated Adaptive T raffic (SCAT) System Philosophy and Benefits[END_REF]. SCOOT is the traffic adaptive version of TRA NSYT. It is based on inductive loop detectors placed on every link to an intersection (just downstream fro m the previous junction). It is designed to optimize signal timing by continuously measuring traffic demand. It adapts to changing traffic using a series of optimization routines. A few seconds before every phase change, the SCOOT split optimizer calculates whether it is better to advance or delay the scheduled chang e by up to 4 s, or to leave it unaltered. In a similar manner, the cycle time of a group of junctions may be incremented or reduced by a few seconds every few minutes.

SCAT contains three types of controllers: a central controller, regional controllers, a nd local controllers. The local controller is installed at each intersection. It collects data fro m detectors, processes it, and executes taken traffic control decisions. The regional controller controls several local controllers. It analyzes local controller information and implements the corresponding signals. The central controller monitors the entire system and provides management support, backups, and fault analysis.

On-road only sensors

Hierarchical architecture for urban traffic system

In this work, Curiac and Vo losencu [START_REF] Daniel-Ioan | Urban Traffic Control System Architecture Based on Wireless Sensor -Actuator Networks[END_REF] propose using a hierarchical three-level strategy to control urban traffic in cities. The lowest level controls single intersection traffic lights by using a wireless sensoractuator network (WSAN) cluster. The second level (base stations), or ZTC unit, executes the zonal traffic coordination. The highest level is Traffic Coord ination unit (TCU) that optimizes the traffic parameters fo r the city; it is imp lemented as a knowledge-based decisional system on a secured computer. This hierarchical control p rovides flexibility for changing the control application (adding or removing monitored entit ies, integrating other applications, etc).

Each WSAN cluster contains several sensors and one actuator (cluster head) -also termed ITU (intersection traffic unit) -that replaces the intersection traffic light controller by aggregating sensor measurement (vehicles queue length, vehicle time to cros s the intersection, etc) and performs the local control algorith m (optimizing traffic parameters for a specific crossroad) and forwards data to other clusters (intersections) or to base station with the role of ZTC. ZTC unit coordinates the ITUs to optimize the traffic by redirect ing the traffic flow, implementing "green waves", etc.

TCU distributes the tasks to the ZTC units and coordinates interactions between different city utilities systems (Disaster Management, Transportation, Public Health, etc). But authors don't give details on the algorithms used in this system. They focus only on the architectures and large lines on the possible relation between elements. Also the exp loitation of TCU decisions is not so detailed. Ho w does the system inform clients or drivers to use these collected information (through displays or internet browsers…). The security mechanisms involved to meet requirements is not highlighted enough. Also the kind of sensors used is not discussed sufficiently. But the architecture us ed is appropriate if co mponents are well crafted.

Sensys

Sensys Networks' [18] wireless vehicle detection system (W VDS) was designed to work seamlessly within SCOOT's architecture, which was orig inally deployed with inductive loops. Each W VDS sensor is placed at the middle of a lane, t ransmitting detection informat ion to a repeater located at short distance away, typically at the nearest pole to the traffic control cabinet.

The system consists of two parts: the wireless sensor network and the access point . Traffic informat ion is generated at the sensor nodes and then transferred to the access point. The communicat ion protocol PEDAMACS (Power Efficient and Delay A ware Mediu m Access Protocol for Sensor Net works) [START_REF] Ergen | PEDAMACS: power efficient and delay aware medium access protocol for sen sor networks[END_REF] is used to increase the lifet ime of the s ystem up to several years, while still guaranteeing the timely arrival of info rmation fro m the sensor nodes to the access point. The traffic management center collects the informat ion fro m each access point to analyze traffic conditions and take appropriat e actions such as adjusting the traffic light durations. It permits archiving and doing Statistics, using its statistics (SNAPS) proprietary software or the Sensys System Manager (SSM ) appliance, by connecting to the AP via Ethernet or cellu lar modem. SNAPS/SSM co llects W VDS parameters (RF signal strengths, link quality, battery levels) to generate d iagnostic alarm alerts for maintaining agency to insure high levels of detector availability.

But authors don't give information on security mechanis ms used to protect communication. A lso the algorith ms used to imp lement networked intersection are not presented. The commun ication protocol used is efficient to enhance the system life time and reduce overhead.

Tubaishat et al.

Tubaishat et al. [START_REF] Tubaishat | Adaptive Traffic Light Control with Wireless Sensor Networks[END_REF] propose a system composed of wireless sensors (on the lanes going in and out the intersection), localized traffic flow model policy, and the coordination of the traffic lights agents. Sensors send vehicles' informat ion (number, speed, queue length, etc) to the neares t Intersection Control Agent (ICA) wh ich determines the intersection flow model and dynamically control the traffic lights in real-t ime by maximizing vehicles flow and reducing wait ing time wh ile maintaining fairness among the other traffic lights. Each traffic light is responsible for controlling traffic on three lanes (the right lane turns right only, center lane goes straight or left and the left lane goes left only). Every lane contains sensors monitoring the traffic before entering and after leaving th e intersection. Multip le intersection agents exchange informat ion to control intersection network.

Authors use the same sensors and commun ication protocol as Sensys [18]. But authors don't give informat ion on the flow policy used or the security requirements.

TSCA & TSTMA

Youcef et al. [START_REF] Khalil | Intelligent Traffic Light Flow Control System Using Wireless Sensors Networks[END_REF] present an adaptive traffic light control system based on Wireless Sensor Network (WSN) fo r controlling the traffic flow sequences on single and mult iple intersections. The system consists of two parts. The first one is the WSN -or TSNs (traffic sensor nodes) -installed on roadside. It collects, periodically, the traffic data (vehicle speed, and length of the vehicles), and sends it to control bo x (basestation) that run control algorith ms.

The second part is the controller, wh ich contains the traffic system communication algorithm (TSCA), and the traffic signals time manipulat ion algorith m (TSTMA). TSCA finds and controls the communicat ion routes between sensors and the BS by a simp le scheduling algorith m that minim izes the collecting data time, as well as the interfacing with the traffic control bo x.

The TSTMA sets time durations for the traffic signals according to the vehicles' nu mber on each segment and the control box applies it on the traffic signals. The selection process of phases is performed every cycle, with dynamic order of phases based on the queues (lanes) that hold maximu m lengths, so that the average queue length and the average waiting time are min imal. For mu ltiple intersections management, TSTMA performs the efficient dynamic control to perform the green wave synchronization. A testbed was also developed and deployed for real measurements.

Authors present the algorithm fo r signal sequences in both isolated and interconnected intersections. But they don't talk about security of commun ication nor kind of sensor used (magnetic…). Also the possibility of connecting the system with internet is not enough highlighted to explain the exp loitation with other ITS systems or to be displayed by drivers.

Adaptive sequenced Traffic Light Control

Zhou et al. [START_REF] Zhou | Adaptive T raffic Light Control in Wireless Sensor Network-based Intelligent Transportation System[END_REF][START_REF] Zhou | Adaptive Traffic Light Control of Multiple Intersections in WSN-Based IT S[END_REF] p ropose an adaptive traffic light control algorith m to adjust both the sequence and length of traffic lights depending on the real time traffic detected instead of just adjusting the green ligh t length in a fixed sequence of traffic lights in the purpose of maximizing intersection throughput and minimizing average wait ing time. The algorithm contains three steps: vehicle detection, green light sequence determination and light length determinatio n. It considers a number of traffic factors (traffic volume, waiting time, vehicle density, etc) to determine green light sequence and the optimal green light length.

In this system the intersection has four d irections, each one with two lanes (going forward and turning left) and on each lane t wo sensors are used (arrival and departure). A total of sixteen sensor nodes are placed on the eight lanes to detect the traffic flow.

Another version of this algorith m is also proposed for mu ltip le intersections in the purpose to schedule the sequences and lengths of traffic lights of intersections cooperatively, to decide next green light establishment and duration. The objective is to increase the network throughput, decrease the average delay and average number o f stops. Light length determination algorith m can determine the duration of next g reen light in each intersection using local traffic volu me and traffic condition fro m neighbor intersections. The proposed algorithm is implemented on transportation testbed iSensNet [24].

Authors propose traffic light management algorithms for both isolated and interconnected intersections. But they don't give information about the commu nicat ion scheme to be used or security protocols. Also the sensors used are not mentioned. The integration of this traffic light system with other ITS possible systems or internet is not discussed.

Combined on-road and on-vehicles sensors

Adaptive Traffic Lights Using Car-to-Car Communication

Gradinescu et al. [START_REF] Gradinescu | Adaptive Traffic Lights Using Carto-Car Communication[END_REF] propose an adaptive traffic light system using vehicles and fixed intersection controller nodes, equipped with short range wireless communication. Veh icles periodically transmit informat ion (position, identificat ion number, speed, direction, state and a timestamp ) about themselves and other cars they know about. The traffic light controller listens to this information to take appropriate decisions.

Neighbor intersection controllers may have wired co mmunications (exchanging additional information).

The upstream signal controller forwards in formation to the downstream one (helping it taking timing decision in advance). At each cycle, the timing plan generation process establishes a plan for the following cycle to minimize control delay (the vehicles' delay caused by the signal control). A fter determining the cycle length, the green splits for each phase are allocated to produce equal degrees of saturation on each link. If the green phase for a lane has finished, but cars keep coming while there is no demand on the conflict ing lanes, the green phase is extended until an acceptable maximu m pedestrians wait ing time.

Authors present the communication scheme between vehicles and the controller but t heir algorith m to optimize the traffic is not so detailed. They don't present how to manage mult iple intersections. Also they don't present any security mechanis m. Ho w does the system react to vehicles don't involving in the system? The authors don't also precise how to join this system to internet or other ITS systems, or drivers displays.

WITS (A Wireless Sensor Network for Intelligent Transportation System)

Chen et al. [START_REF] Chen | WIT S: A Wireless Sensor Network for Intelligent Transportation System[END_REF] present a system for gathering and transferring traffic informat ion. It contains three types of nodes, the vehicle unit (in vehicle), the roadside unit (along both sides of road); and the intersection unit (on the intersection).

The vehicle unit measures the vehicle parameters, such as location (through broadcasting from road side units), etc and transfers them to the roadside units. The roadside unit gathers the information of the vehicles around (each unit collects vehicles' informat ion in one d irection) and transfers it to the intersection unit. The intersection unit receives and analy zes the info rmation fro m other units, and passes them to the strategy sub-system.

Intersection unit needs the number of vehicles in every lane that will reach the intersection before the signal phase ends. Roadside units collect the vehicle message and transform them into the data required by the intersection unit. Roadside unit collects the vehicle informat ion in its scope, aggregates them with the data from its upstream node, and passes them to the downstream node. The strategy sub-system calculates an optimized scheme to control and/or guide the execution sub -system.

Authors suppose that all vehicles dispose the vehicle unit, which is not always the case, and base their system on this assumption. The system use road side unit just to collect and d efuse informat ion but if used also to sense vehicles, system reliability will be enhanced and the above vehicle unit assumption will not be needed. The authors propose a system for travel minimizing, but they don't propose an algorithm fo r interconnected intersections and their algorith m don' take into account other parameters to min imize such as cost and fuel consumption. They don't also give details on security mechanis ms and how drivers will consume in formation (through displays…) or how to connect the s ystem with internet.

RFID & sensor-based Traffic Management

Salama et al. [START_REF] Salama | Intelligent Cross Road T raffic Management System (ICRT MS)[END_REF] design a system for managing and controlling traffic lights based on distributed long range photoelectric sensors in distances prior to and after the traffic lights. The algorithm is based on the total calculated relat ive weight of each road. The system accepts informat ion about any emergency case through RFID. Authors don't give enough details on the algorithm and the possibility to monitor mu ltip le intersections. The position of mu ltiple sensors before the intersection is not justified because vehicles can be detected just with the first sensor and the following can be used as relays. Authors don't also involve security mechanis ms or how to extend the system to work with other ITS systems or how to join it to internet. Using RFID to handle emergency is interesting so that the light signal be emergency enabled.

On-vehicles only sensors

Virtual traffic lights (VTL)

Ferreira et al. [START_REF] Ferreira | Self-Organized Traffic Control[END_REF] replace traffic lights fro m roadside infrastructures to dynamic in-vehicle virtual signs protocol through vehicle-to-vehicle co mmunicat ions. This solution optimizes the throughput of the entire road network (rather than only materialized traffic lighted intersections).

The vehicles act as mobile traffic state sensors resolving intersection conflicts, while elected vehicles act as temporary road junction infrastructures and broadcast traffic light messages to drivers through invehicle displays. This is done by vehicle-to-vehicle co mmunications (using DSRC). It is supposed that all vehicles have: DSRC devices, the same digital road map, and a global positioning system GPS device (to guarantee global time and position synchronization). Each vehicle contains a database of intersections supporting virtual traffic light VTL. When approaching intersection, the vehicle requests for a running VTL to be obeyed.

If it does not exist and crossing conflicts is perceived, a VTL must be created between approaching vehicles. The vehicles approaching the intersection elect one of the m to beco me the leader of the VTL. This leader will be p resented with a red light, and stopped at the intersection while leading it. Du ring this operation, other vehicles are passive, receiv ing traffic light messages and presenting them to the driver. When the green light is in the leader's lane, a new leader must be elected to maintain the VTL.

If there are vehicles stopped before a red light, the leader selects one to become the new leader. If there are no stopped vehicles, then a new leader will be elected by the same process whenever necessary. If the VTL is no longer required, the cycle is interrupted and vehicles proceed without stopping.

Authors don't exp lain the algorith m used in this paradig m and if the system support mu ltip le intersection data or not. Also supposing that all vehicles dispose this system is so strong because applying the algorith m in this manner don't only give poor results as the case of supposition of previous works in this paper but engender catastrophic results because a veh icle suppose having green signal but another one suppose having right priority.

Discussion

Mobile systems, such as vehicular sensor networks, sense the environment with better granularity and at higher scale, co mpared to static sensor networks (particu larly over large areas), and instrument a larger geographical area with a less number of sensors. But with traffic management applications where critical decisions must be taken in real-t ime, such as traffic light monitoring, supposing all vehicles dispose onboard sensors or displays is not yet the case. Therefore, the taken decisions may be wrong and engender catastrophic results in life and materials. So, using VSNs or hybrid on-road and VSNs for requesting real time traffic informat ion and avoiding traffic jam is very efficient but not for taking vital decisions in real time manner (as traffic light monitoring in [START_REF] Ferreira | Self-Organized Traffic Control[END_REF]).

For traffic light signaling systems using WSNs, many parameters may be used to enhance the efficiency of the whole system, added to the collected real-t ime data concerning vehicles (speed, queue length…).

For examp le, a traffic light network of a whole city can use weather and climate (air pollution) values to choose sequences time length.

To benefit co mp letely fro m WSN and VSN capabilities, the conception of a whole system comb ining the two paradig ms is very interesting. A complete ITS system using the two technologies permit the monitoring and management of traffic vehicles and furthermore the traffic infrastructure. In literature and through the previous study a complex but co mplete scheme co mbin ing on road and on vehicles sensors using different kinds of commun ication technologies will clearly permit the overall mon itoring. It is explicit that many d rawbacks need mo re investigation and research to permit fu ll co mmercialization, which is already starting. Applying security co mmun ications, privacy protocols, middle -ware p latforms different wireless technologies will gracefu lly assist WSN based ITS systems.

Conclusion

Road traffic is becoming an important problem in all countries through the world, especially in industrialized countries. Th is state obliges thinking to build ITS systems with high dynamicity and lo w congestion and incidents. WSNs technology helps a lot in designing flexible and cheaper ITS systems because of its easy installation, extension and maintenance. In last years, many projects using WSNs for ITS systems were established. In this paper, we surveyed some existing works on ITS traffic lights, showing their architectural aspect and some weaknesses observed related to each type of architectures to help in designing WSNs based ITS systems.

Figure

  Figure 1WSN Architectures for IT S Figure 2Ad-Hoc architectures 4.1. Ad -hoc paradigm: In this paradig m sensors do not have a specific backbone but they exchange and forward collected data in an ad-hoc manner. Node deployment can be classified into: (i) On road sensor network, where all sensors are imp lanted shallowly inside the carriage or on poles next to the road. In this case, sensors are static. Sensors communicate in a mult i-hop way (without using any infrastructure) as shown in Figure 2.a. (ii) Vehicular sensor network, where all sensors are included in vehicles. In V2V (vehicle -to-vehicle) co mmunicat ions, mobile nodes directly co mmunicate to each other without any need of infrastructure.Decisions can be taken even locally or cooperatively. In Figure2.b, a co mmunicat ion example is presented. (iii) Hybrid ad-hoc sensor network, wh ich is more robust and combines the two previous deployments. Both on -road and in-vehicle sensors, ( Figure2.C) exchange traffic information to cooperatively take correct and real-t ime decisions.

Figure 3

 3 Figure 3 Infrastructure-based architectures Figure 4Hybrid Adhoc and infrastructure 4.2 Infrastructure-enabled monitoring:In this kind of communications all the above sub classes can be found but in the presence of base stations (BSs) that can be also relayed to each other or to servers and Internet through either wireless or wired links. Th is includes wifi, W imax, cellu lar, DSRC, and sensors. DSRC (Dedicated Short-Range communication) [6] is a Wireless Access Protocol used in Vehicular Environments. It operates in the 5.9 GHz frequency range, supports vehicle speed up to 120 mph (192 Km/h), t ransmission range of 300m (up to 1000 m), and data rate of 6 Mbps (up to 27 Mbps).(i) On-road sensors with BS, as in Figure3.a.(ii) Vehicu lar sensors with BS, as in Figure3.b. (iii) Hybrid on-road and on-vehicles sensors with BS, as shown in Figure3.c. The BS makes decisions and notifies actuators. It may also communicate to take g lobal decisions.

  The system can open a comp lete path for such emergency cases from the next traffic until reaching the destination. Reporting the emergency is done by installation o f RFID (Radio Frequency Identificat ion) tag in the concerned vehicles (police cars, ambu lance, fire trucks, etc) and it is detected by RFID reader located besides road. The emergency vehicle driver can be equipped with a mob ile device that sends a specific radio signal to the central control system to open the road.