Mean squared error minimization
for inverse moment problems^{1}^{1}1D. 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 LAASCNRS, supported by a fellowship within the Postdoctoral Programme
of the German Academic Exchange Service.
Abstract
We consider the problem of approximating the unknown density of a measure on , absolutely continuous with respect to some given reference measure , from the only knowledge of finitely many moments of . Given and moments of order , we provide a polynomial which minimizes the mean square error over all polynomials of degree at most . If there is no additional requirement, is obtained as solution of a linear system. In addition, if 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 in as . In general nonnegativity of is not guaranteed even though is nonnegative. However, with this additional nonnegativity requirement one obtains analogous results but computing that minimizes 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 .
Keywords: Moment problems; density estimation; inverse problems; semidefinite programming.
1 Introduction
Estimating the density of an unknown measure is a wellknown 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 [25], the density 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 meanintegrated squared error between and its estimate is most common. Another very popular approach uses wavelets [16, 6, 30], 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 meanintegrated 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 [7, 11, 30] 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 based on the only knowledge of finitely many moments, (say up to order ), an inverse problem from moments. A simple method due to [29] approximates the density by a polynomial of degree at most , so that the moments of the measure matches those of , up to order . However, and in contrast with more sophisticated approaches, the resulting polynomial approximation is not guaranteed to be a density (even though is) as it may takes negative values on the domain of integration. One classical approach to the moment problem is the Padé approximation [3] 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 by maximizing some functional entropy, the most popular being the BoltzmannShannon entropy. In general some type of weak convergence takes place as the degree increases as detailed in [5]. Alternatively the norm of the approximate density is chosen as an objective function [4, 28, 15, 9], which allows to show a stronger convergence in norm. In [21], maximum entropy and Padé approximates have been compared on some numerical experiments. Finally, piecewise polynomial spline based approaches have also been proposed in [14].
Motivation. Our main motivation to study the (inverse) moment problem arises in the context of the socalled generalized problem of moments (GPM). The abstract GPM is a infinitedimensional linear program on some space of Borel measures on and its applications seem endless, see e.g. [17, 18] 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 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 by solving an appropriate finitedimensional 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 semialgebraic sets. In this case one may define a hierarchy (as the number of moments increases) of socalled 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 and whose density is assumed to be in , with no continuity assumption as in previous works. The ultimate goal is to compute an approximation of , based on the only knowledge of finitely many moments (say up to order ) of . In addition, for consistency, it would be highly desirable to also obtain some “convergence” as .
(a) Firstly, we approximate the density by a polynomial of degree which minimizes the mean squared error (or equivalently the norm ) over all polynomials of degree at most . We show that an unconstrained norm minimizer exists, is unique, and coincides with the simple polynomial approximation due to [29]; it can be determined by solving a system of linear equations. It turns out that matches all moments up to degree , 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 in such a basis are just the given moments. Moreover we show that in as , which is the best we can hope for in general since there is no continuity assumption on ; in particular almosteverywhere and almostuniformly on for some subsequence , . 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 is much easier to compute than the corresponding maximum entropy estimate, but it also converges to as in a much stronger sense. For the univariate case, in references [27] and [2] 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 [8]. Notice that in [27] 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 [4, 28] the minimumnorm solution (and not the minimum distance solution) is shown to be unique solution of a system of linear equations. In [15] 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 norm minimizer may not be a density as it may take negative values on . Of course, the nonnegative function also converges to in but it is not a polynomial anymore. So we next propose to obtain a nonnegative polynomial approximation by minimizing the same 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 semialgebraic set one is able to enforce this positivity constraint by using Putinar’s Positivstellensatz [26] 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 interiorpoint algorithms. Moreover, again as in the unconstrained case, we prove the convergence in as (and so almosteverywhere and almostuniform convergence on as well for some subsequence , ) which is far stronger than the weak convergence obtained for the maximum entropy estimate. Notice, in [29, 27] methods for obtaining some nonnegative estimates are discussed, however these estimates do not satisfy the same properties in terms of meansquare error minimization and convergence as in the unconstrained case. In the kernel density element method [2, 8] 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 [27]. The resulting 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 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 by a polynomial of degree at most via unconstrained and constrained norm minimization, respectively; in both cases we also prove the convergence in (and almostuniform convergence on as well for some subsequence) as 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 wellknown Gibbs effect.
2 Notation, definitions and preliminaries
2.1 Notation and definitions
Let (resp. ) denote the ring of real polynomials in the variables (resp. polynomials of degree at most ), whereas (resp. ) denotes its subset of sums of squares (SOS) polynomials (resp. SOS of degree at most ).
With and a given reference measure on , let be the space of functions on whose square is integrable and let be the convex cone of nonnegative elements. Let (resp. ) be the space of continuous functions (resp. continuous nonnegative functions) on . Let be the space of polynomials nonnegative on .
For every the notation stands for the monomial and for every , let whose cardinal is . A polynomial is written
and can be identified with its vector of coefficients in the canonical basis , . Denote by the space of real symmetric matrices, and by the cone of positive semidefinite elements of . For any the notation stands for positive semidefinite. A real sequence , , has a representing measure if there exists some finite Borel measure on such that
Linear functional
Given a real sequence define the Riesz linear functional by:
Moment matrix
Given , the moment matrix of order associated with a sequence , , is the real symmetric matrix with rows and columns indexed by , and whose entry is , for every . If has a representing measure then because
Localizing matrix
With as above and (with ), the localizing matrix of order associated with and is the real symmetric matrix with rows and columns indexed by , and whose entry is , for every . If has a representing measure whose support is contained in the set then because
2.2 Problem statement
We consider the following setting. For compact, let and be finite Borel measures supported on . Assume that the moments of are known and is absolutely continuous with respect to () with RadonNikodým derivative (or density) , with respect to . The density is unknown but we know finitely many moments of , that is,
(1) 
for some .
The issue is to find an estimate for , such that
(2) 
2.3 Maximum entropy estimation
We briefly describe the maximum entropy method due to [12, 13, 5] as a reference for later comparison with the mean squared error approach.
If one chooses the BoltzmannShannon entropy , the resulting estimate with maximumentropy is an optimal solution of the optimization problem
It turns out that an optimal solution is of the form
for some vector . Hence, computing an optimal solution reduces to solving the finitedimensional convex optimization problem
(3) 
where is the given moment information on the unknown density . If , , is a sequence of optimal solutions to (3), then following weak convergence occurs:
(4) 
for all bounded measurable functions continuous almost everywhere. For more details the interested reader is referred to [5].
Since the estimate 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 secondorder optimization algorithms, computing the gradient or Hessian at a current iterate requires evaluating integrals of the form
which is a difficult task in general, except perhaps in small dimension or .
3 The mean squared error approach
In this section we assume that the unknown density is an element of , and we introduce our mean squared error, or norm approach, for density approximation.
3.1 Density approximation as an unconstrained problem
We now restrict to be a polynomial, i.e.
for some vector of coefficients . We first show how to obtain a polynomial estimate of satisfying (2) by solving an unconstrained optimization problem.
Let denote the sequence of moments of on , i.e., , , with
and let denote the moment matrix of order of . This matrix is easily computed since the moments of are known.
Consider the unconstrained optimization problem
(5) 
Proposition 1
Proof: Observe that for every with vector of coefficients ,
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 wellknown optimality conditions for unconstrained, convex quadratic programs and the fact that is nonsingular because for all . Indeed, if for some then necessarily the polynomial with coefficient vector vanishes on the open set , which implies that , in contradiction with .
Finally, let be the vector of coefficients associated with the monomial , . from we deduce
which is the desired result.
Thus the polynomial minimizing the norm distance to coincides with the polynomial approximation due to [29] defined to be a polynomial which satisfies all conditions (2). Note that this is not the case anymore if one uses an norm distance with .
Next, we obtain the following convergence result for the sequence of minimizers of problem (5), .
Proposition 2
Let be compact with nonempty interior and let be finite with for some open set . Let , , be the sequence of minimizers of problem (5). Then as . In particular there is a subsequence , , such that , almost everywhere and almost uniformly on , as .
Proof: Since is compact, is dense in . Hence as there exists a sequence , , with as . Observe that if then as solves problem (5), it holds that for all , which combined with yields the desired result. The last statement follows from [1, Theorem 2.5.2 and 2.5.3].
Note that computing the norm minimizer is equivalent to solving a system of linear equations, whereas computing the maximum entropy estimate requires solving the potentially hard convex optimization problem (3). Moreover, the convergence in Proposition 2 (and so almosteverywhere and almostuniform convergence for a subsequence) is much stronger than the weak convergence (4). On the other hand, unlike the maximum entropy estimate, the approximation 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 of may be a strict subset of . In the case where is not known or its geometry is complicated, one chooses a set 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 when the degree is fixed. In the maximum entropy method of §2.3 no enclosing frame is chosen. In the approach, choosing 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 of the unknown density 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 norm approach always yields a polynomial estimate .
Remark 3
If is a slightly perturbed version of , the resulting numerical error in and side effects caused by ill conditioning may be reduced by considering the regularized problem
(8) 
where is the Euclidean norm of the coefficient vector of , and (fixed) is a regularization parameter approximately of the same order as the noise in . The coefficient vector of an optimal solution of (8) is then given by
(9) 
The effect of small perturbations in on the pointwise approximation accuracy of for 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.
3.2 Density approximation as a constrained optimization problem
As we have just mentioned, the minimizer of problem (5) is not guaranteed to yield a nonnegative approximation even if on . As we next see, the function also converges to for the norm but

it does not satisfy the moments constraints (i.e., does not match the moments of up to order );

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:
(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:
(11) 
where is the moment sequence of and is the localizing matrix associated with and . 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 in (11) follows from a result in [19] which states that if and for all then on .
Lemma 1
If is compact then is dense in with respect to the norm.
Proof: As is compact the polynomials are dense in . Hence there exists a sequence , , such that as . But then the sequence , , with for all , also converges for the norm. Indeed,
So let and be such that and so . As is continuous and is compact, by the StoneWeierstrass theorem there exists a sequence , , that converges to for the supremum norm. Hence for all (for some ). Therefore, the polynomial is positive on and
Therefore we have found a sequence , , such that as .
Proposition 3
Let be compact with nonempty interior. Then problem (11) has an optimal solution for every , and as .
Proof: Fix and consider a minimizing sequence with monotonically decreasing and converging to a given value as . We have for all . Therefore as defines a norm on the finite dimensional space ( has nonempty interior) the whole sequence is contained in the ball . As the feasible set is closed, problem (11) has an optimal solution.
Let be a sequence of optimal solutions of problem (11). By Lemma 1, is dense in . Thus there exists a sequence , , with as . As on then necessarily for all and all . In particular where . Therefore is a feasible solution of problem (11) which yields for all . Combining with yields the desired result.
Via sufficient conditions of positivity
Let denote the standard monomial basis of . Let be a basic compact semialgebraic set defined by for some polynomials , with . Then, with , consider the optimization problem:
(12) 
Since the equality constraints are linear in and the entries of , , the feasible set of (12) is a convex LMI set. Moreover, the objective function is convex quadratic.
Note that whereas the semidefinite constraint in (11) was a relaxation of the nonnegativity constraint on , in (12) a feasible solution is necessarily nonnegative on because the LMI constraint on is a Putinar certificate of positivity on . However, in problem (12) we need to introduce auxiliary matrix variables , whereas (11) is an optimization problem in the original coefficient vector and does not require such a lifting. Thus, problem (11) is computationally easier to handle than problem (12), although both are convex SDPs, which are substantially harder to solve than problem (6).
Proposition 4
Let be compact with nonempty interior. For every , problem (12) has an optimal solution and as .
Proof: With fixed, consider a minimizing sequence for problem (12). As in the proof of Proposition 3 one has . As is feasible,
for some real symmetric matrices , . Rewriting this as
and integrating with respect to yields
with . Hence, for every ,
As , , and , , we conclude that all matrices are bounded. Therefore the minimizing sequence belongs to a closed bounded set and as the mapping is continuous, an optimal solution exists.
From Lemma 1 there exists , , such that as . Using properties of norms, , and so as . Moreover, as is strictly positive on , by Putinar’s Positivstellensatz [26], there exists such that
for some real matrices , . So letting , the polynomial is a feasible solution of problem (12) whenever and with value . Hence as , and by monotonicity of the sequence , , the result follows.
Remark 4
Remark 5
Our approach can handle quite general sets and as support and frame for the unknown measure in the unconstrained and constrained cases, the only restriction being that (i) and are basic compact semialgebraic 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
must be evaluated. Such evaluations may be difficult as soon as . In particular in higher dimensions, cubature formulas for approximating such integrals are difficult to obtain if is not a box or a simplex.
4 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 and its estimate by the average error
and the maximum pointwise error
In some examples we also consider and , 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 () 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 . In order to solve the unconstrained, concave optimization problem (3), we apply the Matlab Optimization Toolbox command fminunc as a blackbox 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
on given its moments , which was considered as a test case in [24]. This example is a priori very favorable for our technique since the desired is a polynomial itself. When solving problem (5) for we do obtain the correct solution in less than secs. Thus, the pointwise approximation is much better than the one achieved in [24]. Moreover, 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 . We solve problem (5) for where the maximal relative componentwise perturbation between and is less than . The pointwise error between and the solution for the perturbed moment vector is sufficiently small, as pictured in Figure 1.