Laurent Gosse
email: l.gosse@ba.iac.cnr.it

Olof Runborg
email: olofr@nada.kth.se

Existence, uniqueness and a constructive solution algorithm for a class of finite Markov moment problems

We consider a class of finite Markov moment problems with arbitrary number of positive and negative branches. We show criteria for the existence and uniqueness of solutions, and we characterize in detail the nonunique solution families. Moreover, we present a constructive algorithm to solve the moment problems numerically and prove that the algorithm computes the right solution.

Introduction

We aim at inverting a moment system often associated with the prestigious name of Markov. The original form of the problem is the following. Given a finite set of moments m k for k = 1, . . . , K, find a bounded measurable density function f satisfying

m k = R x k-1 f (x)dx, 0 ≤ f ≤ 1, k = 1, . . . , K. (1)
Condition for the existence of solutions f (x) to this problem is classical [START_REF] Akhiezer | The Classical Moment Problem and Some Related Questions in Analysis[END_REF][START_REF] Akhiezer | Some questions in the theory of moments[END_REF]. In general solutions are not unique, unless more conditions are given, e.g. based on entropy minimization [START_REF] Brenier | Équations de moment et conditions d'entropie pour des modèles cinétiques (French) [Moment equations and entropy conditions for kinetic models[END_REF][START_REF] Brenier | A kinetic formulation for multibranch entropy solutions of scalar conservation laws[END_REF] or L ∞ -minimization [START_REF] Norris | Optimal solutions to the L ∞ moment problem with lattice bounds[END_REF][START_REF] Lewis | Superresolution in the Markov moment problem[END_REF]. A typical result is that the unique solution for even K is piecewise constant, taking values in {0, 1}. More precisely, if K = 2n then f is of the form

f (x) = n j=1 χ [yi,xi] (x) (2)
where χ I (x) is the characteristic function for the interval I and

y 1 < x 1 < y 2 < x 2 < • • • < y n < x n . (3)
See Theorem 3 below in Section 4 and consult e.g. [START_REF] Curto | The truncated complex K-moment problem[END_REF][START_REF] Diaconis | The Markov moment problem and de Finetti's theorem[END_REF][START_REF] Krein | The Markov moment problem and extremal problems[END_REF][START_REF] Simon | The classical moment problem as a self-adjoint finite difference operator[END_REF][START_REF] Talenti | Recovering a function from a finite number of moments[END_REF] for general background on moment problems.

A reduced form of the finite moment problem is to search for solutions to (1) which are precisely of the form [START_REF] Akhiezer | Some questions in the theory of moments[END_REF][START_REF] Brenier | Équations de moment et conditions d'entropie pour des modèles cinétiques (French) [Moment equations and entropy conditions for kinetic models[END_REF]. One then obtains an algebraic problem for the branch values,

m k = 1 k n j=1
x k jy k j , k = 1, . . . , K = 2n.

Finding {x j } and {y j } from {m k } is an ill-conditioned problem when the branch values of the solution come close to each other; the Jacobian of the problem is a Vandermonde matrix and iterative numerical resolution routines require extremely good starting guesses when the matrix degenerates. For less than four moments a direct method based on solving polynomial equations was presented in [START_REF] Runborg | Some new results in multiphase geometrical optics[END_REF]. Routines based on the Simplex algorithm have been proposed in [START_REF] Norris | Optimal solutions to the L ∞ moment problem with lattice bounds[END_REF].

Another algorithm was presented by Koborov, Sklyar and Fardigola in [START_REF] Korobov | Time-optimality and the power moment problem[END_REF][START_REF] Sklyar | The Markov power moment problem in problems of controllability and frequency extinguishing for the wave equation on a half-axis[END_REF] in the slightly modified setting where f takes values in {-1, 1} instead of {0, 1}. It consists of solving a sequence of high degree polynomial equations, constructed through a rather intricate process with unclear stability properties. In [START_REF] Gosse | Finite moment problems and applications to multiphase computations in geometric optics[END_REF] we showed that this algorithm can be drastically simplified and adapted to [START_REF] Brenier | A kinetic formulation for multibranch entropy solutions of scalar conservation laws[END_REF]. Later, in [START_REF] Gosse | Resolution of the finite Markov moment problem[END_REF], we also gave a direct proof that the simplified algorithm indeed computes the correct solution, relying on the classical Newton's identities and Toeplitz matrix theory.

The moment problem has many applications in for instance probability and statistics [START_REF] Gamboa | Large-deviations for random power moment problems[END_REF][START_REF] Diaconis | Application of the method of moments in probability and statistics[END_REF], but also in areas like wave modulation [START_REF] Czarkowski | Solving the optimal PWM problem for single-phase inverters[END_REF][START_REF] Sezan | Incorporation of a-priori moment information into signal recovery and synthesis problems[END_REF] and "shape from moments" inverse problems [START_REF] Golub | A stable numerical method for inverting shape from moments[END_REF]. Our own motivation comes from a quite different field, namely multiphase geometrical optics [START_REF] Brenier | Équations de moment et conditions d'entropie pour des modèles cinétiques (French) [Moment equations and entropy conditions for kinetic models[END_REF][START_REF] Brenier | A kinetic formulation for multibranch entropy solutions of scalar conservation laws[END_REF][START_REF] Gosse | Using K-branch entropy solutions for multivalued geometric optics computations[END_REF][START_REF] Gosse | Two Moment Systems for Computing Multiphase Semiclassical Limits of the Schrödinger Equation[END_REF][START_REF] Gosse | Finite moment problems and applications to multiphase computations in geometric optics[END_REF][START_REF] Runborg | Some new results in multiphase geometrical optics[END_REF]. In this application one needs to solve a system of nonlinear hyperbolic conservation laws. To evaluate the flux function in the partial differential equations (PDEs) a system like (4) must be solved. In a finite difference method this means that the system must be inverted once for every point in the computational grid, repeatedly in every timestep. It is thus important that the inversion can be done fast and accurately; this difficulty has been a bottleneck in computations. In [START_REF] Gosse | Finite moment problems and applications to multiphase computations in geometric optics[END_REF] we used the simplified algorithm mentioned above for numerical implementation inside a shock-capturing finite difference solver. It is our aim here to develop better algorithms and understanding to open the way for the processing of intricate wave-fields with large K, and thus complement the seminal paper [START_REF] Brenier | A kinetic formulation for multibranch entropy solutions of scalar conservation laws[END_REF] where the multiphase geometrical optics PDEs were first proposed.

In this paper we are concerned with a generalization of (4). In the geometrical optics application, the number of moments K is typically not even and one can have a variable number of positive (x k ) and negative (y k ) branches. We thus consider the following problem

m k = nx j=1 x k j - ny j=1 y k j , k = 1, . . . , K, (5)
where n x + n y = K but where n x and n y are not necessarily equal. We study existence and uniqueness of solutions to this problem (Theorem 2). In particular we are interested in how and when uniqueness is lost. For these cases we characterize the family of solutions that exists. The reason is to understand what happens numerically close to degenerate solutions, which is an important feature in the application we have in mind: In the exact solution to the multiphase geometrical optics PDEs the moment problem is typically degenerate for large domains; the numerical approximation is almost degenerate. We also give constructive algorithms to solve [START_REF] Curto | The truncated complex K-moment problem[END_REF] and prove that they generate the right solution (Theorem 1). In a future paper we will study the numerical stability of these algorithms. Experimentally we note, for instance, that to compute the next moment, Algorithm 3 is much more stable than Algorithm 1. The difficulty lies in understanding perturbations around degenerate solutions, which is where the algorithms are most unstable. For this the insights of this paper will be of importance.

Remark 1

The problem (5) can be cast in the form of (1) if one demands that the density function f (x) is of the form

f (x) = nx j=1 sgn(x j ) [H(x) -H(x -|x j |)] - ny j=1 sgn(y j ) [H(x) -H(x -|y j |)] ,
(6) and we rescale the moments m k → km k . For the case n x = n y = n and K = 2n with interlaced branch values (3) this reduces to (2). This paper is organized as follows. In Section 2 we present the algorithms for solving [START_REF] Curto | The truncated complex K-moment problem[END_REF]. Notation and various ways of describing a solution is subsequently introduced in Section 3. Next we derive conditions for existence and uniqueness of solutions in Section 4 and also discuss various properties of the solution, in particular when it is not unique. A theorem proving the correctness of the algorithms is proved in Section 5. Finally, in Section 6, we give additional properties of the elements of our algorithms, and use these to relate our results back to the classical Markov theory.

Algorithms

In this section we detail the algorithms that we propose for solving [START_REF] Curto | The truncated complex K-moment problem[END_REF]. The solution that we obtain is what we call the minimal degree solution, meaning that when the solution is not unique as many branch values as possible are zero. See Section 4 for a precise definition. The algorithms goes as follows; they may fail in case there is no solution to [START_REF] Curto | The truncated complex K-moment problem[END_REF].

Algorithm 1 (Computing {x j } and {y j })

1. Construct the sequence {a k } as follows. Set a 0 = 1 and a k = 0 for k < 0.

For 1 ≤ k ≤ K, let the elements be given as the solution to

     1 -m 1 2 . . . . . . . . . -m K-1 . . . -m 1 K           a 1 a 2 . . . a K      =      m 1 m 2 . . . m K      . ( 7
)
2. Construct the matrix A 1 ∈ R nx×nx as

A 1 =     
Ã0 v = x Ã1 v, (8)
to get the {x j } values of the minimal degree solution to (5).

5.

To compute the {y j } values, the same process is used with m k replaced by -m k and the roles of n x and n y interchanged.

An alternative to Algorithm 1 is as follows:

Algorithm 2 (Computing {x j } and {y j })

1. Construct the matrices Ã0 and Ã1 as in steps 1-3 in Algorithm 1.

2. Denote the first column vector in Ã0 by ã0 by and solve

Ã1 c ′ = -ã 0 , c ′ = (c 1 , c 2 , . . . , c ñx ) T . (9)
3. Construct the polynomial

P (z) = c ñx + c ñx-1 z + • • • + c 1 z ñx-1 + z ñx .
The roots of P (z) are the {x j } values of the minimal degree solution to (5) (possibly together with some zeros).

4. To compute the {y j } values, the same process is used with m k replaced by -m k and the roles of n x and n y interchanged.

Remark 2 We note that the values of a k in the definition (7) are independent of K, since the system matrix is triangular. We therefore consider the sequence without reference to K in any other respect than the fact that we are only able to compute elements with k ≤ K when we are given K moments. The largest index of the a k -sequence appearing in the matrix A 1 is n y + n x -1 < K. In the matrices Ã0 , Ã1 it is ñy + ñx = n yn x + 2ñ x ≤ n y + n x = K. Hence all three matrices can be constructed from the first K moments. Some properties of the A 1 matrix are detailed in Section 6.

Sometimes one is not interested in finding the individual {x j } and {y j } branch values but just wants the higher moments, defined as

m k = nx j=1 x k j - ny j=1 y k j , (10)
but now for k > K, given a solution {x j }∪{y j } to [START_REF] Curto | The truncated complex K-moment problem[END_REF]. (That this is well-defined is shown later in Theorem 2.) For this case there is another algorithm, which has empirically proven to be more stable than first computing {x j } and {y j } from Algorithm 1 or 2, and then entering the values into [START_REF] Gamboa | Large-deviations for random power moment problems[END_REF]. We stress that this is precisely what is needed in order to compute K-multivalued solutions of the inviscid Burger's equation in geometrical optics, following the ideas of [START_REF] Brenier | A kinetic formulation for multibranch entropy solutions of scalar conservation laws[END_REF].

Algorithm 3 (Computing m K+1 )

1. Construct the A 1 matrix as in steps 1-2 of Algorithm 1.

Let

a 0 = (a ny+1 , a ny+2 , . . . , a ny+nx ) T ∈ R nx ,
and let c = (c 1 , c 2 , . . . , c nx ) T be one solution to

A 1 c = -a 0 . (11)
3. The next moment is given by

m K+1 = -(K + 1) nx j=1 c j a K+1-j - K j=1 m j a K+1-j .
We recall that Algorithm 1 has been shown to be numerically efficient in the paper [START_REF] Gosse | Finite moment problems and applications to multiphase computations in geometric optics[END_REF]. The justification of these algorithms is given in Section 5 where we show the following theorem:

Theorem 1 If a solution to (5) exists then:

(i) In Algorithm 1, the matrix Ã1 is non-singular. The generalized eigenvalue problem in ( 8) is well-defined and the generalized eigenvalues (counting algebraic multiplicity) are the {x j }-values of the minimal degree solution to (5) plus ñx -D min zeros. (See [START_REF] Norris | Optimal solutions to the L ∞ moment problem with lattice bounds[END_REF] for the definition of D min .)

(ii) In Algorithm 2, c ′ is well defined, P (z) = det(zI -Ã-1 1 Ã0 ) ( 12
)
and the roots of P (z) are the {x j }-values of the minimal degree solution to (5) plus ñx -D min zeros.

(iii) In Algorithm 3, the computed moment satisfies

m K+1 = nx j=1 x K+1 j - ny j=1 y K+1 j ,
for all solutions {x j } ∪ {y j } to [START_REF] Curto | The truncated complex K-moment problem[END_REF].

We postpone the proof of Theorem 1 to Section 5. We just note here that the last point in Algorithms 1 and 2 can easily be explained by the symmetry of the problem. Indeed, the negative of ( 5)

-m k = ny j=1 y k j - nx j=1 x k j , k = 1, . . . , K,
is of the same form as (5) itself, with the roles of n x , {x j } and n y , {y j } interchanged.

Preliminaries

We will use three different ways of describing the solution to [START_REF] Curto | The truncated complex K-moment problem[END_REF]. First we have a set of numbers {x j } nx j=1 and {y j } ny j=1 , solving [START_REF] Curto | The truncated complex K-moment problem[END_REF]. We call those numbers branch values. Second, we have a pair of polynomials (p, q) of degrees at most n x and n y respectively in the z variable. Third, we have a pair of coefficient vectors c = (c 0 , . . . , c nx ) T ∈ R nx+1 and d = (d 0 , . . . , d nx ) T ∈ R ny+1 . These three representations are related as [START_REF] Gosse | Two Moment Systems for Computing Multiphase Semiclassical Limits of the Schrödinger Equation[END_REF] and

p(z) = (1 -x 1 z) • • • (1 -x nx z) = c 0 + c 1 z + • • • + c nx-1 z nx-1 + c nx z nx ,
q(z) = (1 -y 1 z) • • • (1 -y ny z) = d 0 + d 1 z + • • • + d ny-1 z ny-1 + d ny z ny . (14)
It is clear that there is a one-to-one correspondence between these ways of describing the solution, if we disregard the ambiguity in the ordering of the numbers {x j } and {y j }. Generally, we will use the notation Deg(p) to denote the degree of a polynomial p, and, for a given coefficient vector c, we systematically write P c to denote the corresponding polynomial [START_REF] Gosse | Two Moment Systems for Computing Multiphase Semiclassical Limits of the Schrödinger Equation[END_REF].

Definition 1 We call the pair of polynomials (p, q) a (polynomial) solution to (5) if 1. The degrees of p and q are at most n x and n y ,

Deg(p) ≤ n x , Deg(q) ≤ n y , (15)
2. They are normalized to one at the origin,

p(0) = q(0) = 1, (16)
3. Their roots {x j } and {ỹ j } satisfy

m k = Deg(p) j=1 x-k j - Deg(q) j=1 ỹ-k j , k = 1, . . . , K. (17)
We note that the roots cannot be zero because of [START_REF] Korobov | Time-optimality and the power moment problem[END_REF].

Next:

Definition 2 A pair of vectors c = (c 0 , . . . , c nx ) T ∈ R nx+1 and d = (d 0 , . . . , d ny ) T ∈ R ny+1
is said to be a (coefficient) solution to (5) if the corresponding pair (P c , P d ) (13)-( 14) realizes a polynomial solution to [START_REF] Curto | The truncated complex K-moment problem[END_REF].

The number of branch values are always n x and n y respectively. Some of them may be zero, and they do not need to be distinct. The number of non-zero branch values are Deg(p) and Deg(q) respectively. The degree of a solution can then also be defined.

Definition 3

The degree of a solution to ( 5) is the number of non-zero x jvalues. This number is equivalent to Deg(p).

Given any polynomial pair satisfying [START_REF] Korobov | Time-optimality and the power moment problem[END_REF], we say that it generates the moment sequence {m k } if m k is given by ( 17) for all k. In turn, each sequence of moments {m k } generates the corresponding {a k } sequence through [START_REF] Diaconis | Application of the method of moments in probability and statistics[END_REF]. We define the big matrix

A =      a ny+1 a ny . . . a ny-nx+1 a ny+2 a ny+1 . . . a ny-nx+2 . . . . . . . . . . . . a ny+nx a ny+nx-1 . . . a ny      ∈ R nx×(nx+1) .
We let the columns of A be denoted a 0 , . . . , a nx and we note that

A =   | | a 0 • • • a nx | |   =   | A 0 a nx |   =   | a 0 A 1 |   , (18)
Hence, A 0 and A 1 constitutes the first and last n x columns of A respectively. When a 0 ∈ range A 1 and a 0 = 0, let

D min = argmin j>0 a 0 ∈ span{a 1 , . . . , a j }, (19)
and set D min = 0 if a 0 = 0. Moreover, define

D max = D min + n x -rank A 1 . (20)

Existence and uniqueness of solutions

In this section we prove results on the existence and uniqueness of solutions to [START_REF] Curto | The truncated complex K-moment problem[END_REF]. We aim at establishing the following theorem:

Theorem 2

(i) There exists a solution to [START_REF] Curto | The truncated complex K-moment problem[END_REF] if and only if

a 0 ∈ range(A 1 ). ( 21
) (ii) If d is the degree of a solution to (5), then D min ≤ d ≤ D max .
(iii) When ( 21) holds, there is a unique solution (p * , q * ) of minimal degree D min . For this solution, x j = y i for all indices i, j representing non-zero branch values. Moreover, Deg(q * ) ≤ n yn x + rank A 1 with equality if

D min < rank A 1 .
(iv) When ( 21) holds, a polynomial pair (p, q) is a solution if and only if p = p * r and q = q * r where r(z) is a polynomial satisfying r(0) = 1 and Deg(r) ≤ D max -D min .

(v) The minimal degree solution is the only solution to (5) if and only if the matrix A 1 is non-singular.

(vi) Let {x j } and {y j } be a solution to [START_REF] Curto | The truncated complex K-moment problem[END_REF]. Then the higher moments defined in [START_REF] Gamboa | Large-deviations for random power moment problems[END_REF] are well-defined.

Let us proceed with several remarks:

Remark 3 In particular it follows from (i) that there exists a solution as soon as the matrix A 1 is non-singular.

Remark 4 Since (5) is a system of polynomial equations of degree K, one could expect there to be a finite number of solutions, typically K solutions. However, because of the special structure of the equations there is either one unique solution (when A 1 is non-singular) or inifintely many solutions (when A 1 is singular).

Remark 5

The form (p * r, q * r) of solutions can also be stated as follows: All solutions have a core set of values {x j }, j = 1, . . . , Deg(p * ) = D min and {y i }, i = 1, . . . , Deg(q * )corresponding to non-zero branch values of the minimal degree solution, where x j = y i for all those i, j. One can then add an optional set of non-zero branch values {x Dmin+j }, and {y Deg(q * )+j }, for j = 1, . . . , D max -D min such that x Dmin+j = y Deg(q * )+j .

To prove this theorem we first establish some utility results in the next subsection. We then derive different ways of characterizing the solution in Section 4.2, which are subsequently used to prove Theorem 2 in Section 4.3.

Utility results

We start with a useful lemma on Taylor coefficients for a product of functions: Lemma 1 Suppose f , g and h are analytic functions in a neighborhood of zero satsifying f (z) = g(z)h(z). Let f have the Taylor expansion

f (z) = ∞ k=0 f k z k ,
and let {g k }, and {h k } be the corresponding coefficients for g(x) and h(x) respectively. Then

f k = k j=0 g j h k-j . ( 22
)
Proof: Since the functions are analytic the coefficients are given as

f k = 1 k! d k dz k f (z) z=0 = 1 k! d k dz k g(z)h(z) z=0 = 1 k! k j=0 c jk g (j) (0)h (k-j) (0),
where c jk = k!/j!(kj)! are the binomial coefficients. But g (j) (0) = j!g j and h (k-j) (0) = (kj)!h k-j and therefore ( 22) follows.

Remark 6 The sum [START_REF] Sezan | Incorporation of a-priori moment information into signal recovery and synthesis problems[END_REF] is in fact precisely an elementwise description of multiplication of a lower triangular k × k Toeplitz matrix by a vector. In the notation of [START_REF] Gosse | Resolution of the finite Markov moment problem[END_REF], it would read f = T (g)h.

As was already known by Markov, the exponential transform of the moment sequence plays an important role in the analysis of these problems, see e.g. [START_REF] Akhiezer | The Classical Moment Problem and Some Related Questions in Analysis[END_REF][START_REF] Akhiezer | Some questions in the theory of moments[END_REF]. We show here that {a k } is a version of the exponential transform of {m k }.

Lemma 2 Suppose {m k } is generated by the polynomials p(z) and q(z) and {a k } is generated by {m k }. Let m(z) be defined as

m(z) = m 1 z + 1 2 m 2 z 2 + 1 3 m 3 z 3 + • • • . ( 23
)
Then if [START_REF] Korobov | Time-optimality and the power moment problem[END_REF] holds,

e m(z) = q(z) p(z) = a 0 + a 1 z + a 2 z 2 + • • • , (24)
written as its Taylor expansion around z = 0.

Proof: Let us first show that m(z) is a well-defined analytic function at zero. We have

m(z) = ∞ k=0 m k z k k = ∞ k=0 nx j=1 x k j z k k - ∞ k=0 ny j=1 y k j z k k = -nx j=1 log(1 -x j z) + ny j=1 log(1 -y j z).
The last step is allowed when |z| < 1/ max ij (|x j |, |y i |), which is true for small enough z since p(0) = 0. This also shows that the function is analytic at zero. Moreover,

e m(z) = ny j=1 (1 -y j z) nx j=1 (1 -x j z) = q(z) p(z) .
Finally, setting a(z) := exp(m(z)) and differentiating gives

za ′ (z) = zm ′ (z)a(z),
where all three functions are analytic at zero. Let a(z) have the Taylor coefficients

{ã k }. Then za ′ (z) = ã1 z + 2ã 2 z 2 + 3ã 3 z 3 • • • and clearly zm ′ (z) = m 1 z + m 2 z 2 + • • • . By Lemma 1, for k ≥ 1, kã k = k j=1 m j ãk-j .
Since ã0 = q(0)/p(0) = 1, we see that a k and ãk satisfy the same non-singular linear system of equations ( 7), and therefore a k = ãk , showing [START_REF] Sklyar | The Markov power moment problem in problems of controllability and frequency extinguishing for the wave equation on a half-axis[END_REF]. We now have the following basic characterization of a solution.

Lemma 3 Suppose p(z) and q(z) are two polynomials satisfying [START_REF] Gosse | Resolution of the finite Markov moment problem[END_REF][START_REF] Korobov | Time-optimality and the power moment problem[END_REF]. They form a polynomial solution to [START_REF] Curto | The truncated complex K-moment problem[END_REF] if and only if their quotient has the Taylor expansion around z = 0

q(z) p(z) = a 0 + a 1 z + • • • + a K z K + O z K+1 , (25)
where {a k } is generated by {m k }. Moreover, if (p, q) is a solution then (p, q) is also a solution if and only if the pair satisfies [START_REF] Gosse | Resolution of the finite Markov moment problem[END_REF][START_REF] Korobov | Time-optimality and the power moment problem[END_REF] and p/q = p/q where these fractions are defined.

Proof: Let { mk } be generated by p and q and suppose (25) holds. Then, as in the of proof of Lemma 2 for 1

≤ k ≤ K ka k = k j=1 mj a k-j .
Since {m k } satisfy the linear system (7), we have after subtraction,

m n -mn = - n-1 k=1 (m k -mk )a n-k , m 1 = m1 ,
for n = 2, . . . , K. By induction mk = m k for 1 ≤ k ≤ K, showing that (p, q) solves ( 5). On the other hand, if (p, q) is a solution, then (25) must hold by [START_REF] Sklyar | The Markov power moment problem in problems of controllability and frequency extinguishing for the wave equation on a half-axis[END_REF] in Lemma 2.

For the last statement, the "if" part is obvious since both pairs then satisfy [START_REF] Talenti | Recovering a function from a finite number of moments[END_REF]. To show the "only if" part, suppose both (p, q) and (p, q) are solutions. By definition they satisfy [START_REF] Gosse | Resolution of the finite Markov moment problem[END_REF][START_REF] Korobov | Time-optimality and the power moment problem[END_REF], and by [START_REF] Talenti | Recovering a function from a finite number of moments[END_REF],

q(z) p(z) - q(z) p(z) = q(z)p(z) -p(z)q(z) p(z)p(z) = O(z K+1 ).
Since p(0)p(0) = 1 we must have that (q(z)p(z)p(z)q(z))/z K+1 is bounded as z → 0. But since the degree of qppq is at most K = n x + n y this is only possible if it is identically zero. Hence q(z)p(z) = p(z)q(z) which concludes the proof.

Characterization of the solution

In this section we show three Propositions that characterize solutions to [START_REF] Curto | The truncated complex K-moment problem[END_REF] in terms of polynomials, coefficient vectors and the column vectors of the Amatrix in [START_REF] Lewis | Superresolution in the Markov moment problem[END_REF]. We start by expressing the uniqueness properties of the solution in terms of its polyomial representation.

Proposition 1 Suppose the pairs (p, q) and (p, q) are both polynomial solutions to [START_REF] Curto | The truncated complex K-moment problem[END_REF]. Then,

(i) Deg(p) -Deg(q) = Deg(p) -Deg(q).
(ii) If Deg(p) ≤ Deg(p), and if there is no polynomial r(z) such that p = pr, then there is another solution (p, q) with Deg(p) < Deg(p). In particular, if Deg(p) = Deg(p) but p = p, there is such a lower degree solution. Proof:

(i) The statement follows directly from Lemma 3, since qp = pq implies that Deg(q) + Deg(p) = Deg(p) + Deg(q).

(ii) We let p(z) = r p (z)p(z) + s p (z), q(z) = r q (z)q(z) + s q (z), be the unique polynomial decomposition of (p, q) such that r p , r q , s p , s q are polynomials, Deg(s p ) < Deg(p) and Deg(s q ) < Deg(q). Since pq = pq by Lemma 3, we get pq(r qr p ) = qs pps q .

Unless r q = r p the degree of the left hand side is at least Deg(p) + Deg(q), while the degree of the right hand side is at most max (Deg(q) + Deg(s p ), Deg(p) + Deg(s q )) < Deg(q) + Deg(p).

Hence, r q = r p and qs p = ps q . Since q, p ≡ 0 it follows that either s p and s q are both zero or both non-zero. Suppose s p ≡ 0 and s q ≡ 0. Write s p (z) = z mp sp (z) and s q (z) = z mq sq (z) where sp (0) = 0 and sq (0) = 0. Since z mp sp (z)q(z) = z mq sq (z)p(z)

and also q(0) = p(0) = 1, the lowest degree term in the left and right hand side polynomials are z mp and z mq respectively, and therefore m p = m q . Consequently, sp (z)q(z) = sq (z)p(z), and sp (0) = sq (0). We can then take p(z) = sp (z)/s p (0) and q(z) = sq (z)/s q (0). They satisfy p(z)q(z) = q(z)p(z), p(0) = q(0) = 1, while Deg(p) = Deg(s p ) ≤ Deg(s p ) < Deg(p) and similarly Deg(q) < Deg(q) ≤ n y . Hence (p, q) is a polynomial solution by Lemma 3. It has degree strictly less than (p, q), which shows the first statement in (ii). If Deg(p) = Deg(p) and p = p then there is no r(z) satisfying the requirements, showing the second statement in (ii).

(iii) We finally let r(z) be any polynomial with Deg(r) ≤ Deg(p) -Deg(p) and r(0) = 1. We then set p = pr and q = qr. These polynomials trivially satisfy ( 16) and [START_REF] Talenti | Recovering a function from a finite number of moments[END_REF]. Since Deg(p) = Deg(r) + Deg(p) ≤ Deg(p) ≤ n x and Deg(q) = Deg(r) + Deg(q) ≤ Deg(p) -Deg(p) + Deg(q) = Deg(q) ≤ n y , they also satisfy [START_REF] Gosse | Resolution of the finite Markov moment problem[END_REF] and thus are a polynomial solution by Lemma 3. In particular we can take r(z) of degree m.

A solution to (5) can also be characterized in terms of the coefficient vectors. We have the following Proposition.

d k = min(k,nx) j=0 c j a k-j , k = 0, . . . , n y . (26)
Proof: Suppose first that c is in the null-space of A, c 0 = 1 and {d k } is given by (26). Extend the coefficient sequences by setting c k = 0 for k > n x and d k = 0 for k > n y . Since c is in the null-space of A, we get k j=0 c j a k-j = 0 when n y + 1 ≤ k ≤ n x + n x = K, and in conclusion

d k = k j=0 c j a k-j , k = 0, . . . , K. (27)
Upon noting that {c k } ∞ k=0 and {d k } ∞ k=0 are the Taylor coefficients of P c and P d , and since P c (0

) = c 0 = 1, P d (0) = d 0 = a 0 c 0 = 1, Lemma 1 shows that P d (z) = P c (z) a 0 + a 1 z + • • • + a K z K + O z K+1 , (28)
and by Lemma 3 we have that (P c , P d ) is a solution to [START_REF] Curto | The truncated complex K-moment problem[END_REF]. Conversely, if (P c , P d ) is a solution, then c 0 = P c (0) = 1 and by Lemma 1 we get that (27) holds. For k = n y + 1, . . . , K this also implies that c is in the null-space of A.

The final Proposition of this section relates the degree of the solution to the column vectors of A and the linear spaces they span.

Proposition 3 Let V j = span{a 1 , . . . , a j } and V 0 j = span{a 0 , . . . , a j }. Set V 0 = V 0 -1 = ∅. Then (i) There is a solution if and only if a 0 ∈ V nx = Range(A 1 ).

(ii) There is a solution of degree j ≥ 0 if and only if

a 0 ∈ V j , and a j ∈ V 0 j-1 . ( 29
) (iii) When a 0 ∈ V nx then a 0 ∈ V d , V 0 d = V d , if and only if d ≥ D max .
(iv) When a 0 ∈ V nx the vectors a 1 , . . . , a Dmin , (when D min > 0) a Dmax+1 , . . . , a nx , (when D max < n x ), are all linearly independent. Moreover,

a j ∈ V Dmin , V j = V Dmin , j = D min , . . . , D max .
Proof:

(i) By Proposition 2 there exists a solution to (5) if and only if there is a coefficient vector c = (1, c ′ ) T in the null-space of A, i.e.

Ac = A 1 c + a 0 = 0.
But such a vector c exists if and only if a 0 is in the range of A 1 . This shows (i).

(ii) Again by Proposition 2 there is a solution of degree j if and only if there is a vector c = (c 0 , c 1 , . . . , c j , 0, . . . , 0) T such that

0 = Ac = c 0 a 0 + c 1 a 1 + • • • + c j a j , (30)
with c j = 0 and c 0 = 1. For j = 0 this is clearly equivalent to a 0 = 0 or a 0 ∈ V 0 = V 0 -1 . For j > 0 the existence of c j -coefficients satisfying (30) is equivalent to the left condition in (29). Moreover, if a j = V 0 j-1 = span{a 0 , . . . , a j-1 }, then we must have c j = 0 to satisfy (30), and c cannot represent a solution of degree j. On the other hand, if c j = 0 and a

j = c ′ 0 a 0 + • • • + c ′ j-1 a j-1 for some non-zero coefficients c ′ k , then a 0 + c ′′ 1 a 1 + • • • + c ′′ j-1 a j-1 + a j = 0, with c ′′ k = (1 + c ′ 0 )c k -c ′ k
, represents a solution of degree j. This shows (ii).

(iii) The statement is obvious in case D min = 0. If D min > 0 there are scalars such that

a 0 = v 1 a 1 + • • • + v Dmin a Dmin , (31)
by [START_REF] Norris | Optimal solutions to the L ∞ moment problem with lattice bounds[END_REF]. Hence, a 0 ∈ V Dmin and since the V j spaces are nested, V j ⊂ V j+1 , we have a 0 ∈ V d for d ≥ D min . Moreover, the minimal property of D min ensures that v Dmin = 0 in (31), so that a 0 ∈ V d when d < D min .

(iv) To show that when D min > 0 the vectors a 1 , . . . , a Dmin are linearly independent, we use (31) and note that P c (z) with c = (1, -v 1 , . . . , -v Dmin , 0, . . . , 0) T is a polynomial solution to [START_REF] Curto | The truncated complex K-moment problem[END_REF]. Suppose now that the there are non-zero coefficients c ′ j such that

c ′ 1 a 1 + • • • + c ′ Dmin a Dmin = 0. Then P c ′ with c ′ = (1, c ′ 1 -v 1 , .
. . , c Dminv Dmin , 0, . . . , 0) T is another polynomial solution to [START_REF] Curto | The truncated complex K-moment problem[END_REF]. Moreover, by the minimality property of D min we must have c Dmin -v Dmin = 0 and therefore Deg(P c ) = Deg(P c ′ ) = D min . But by (ii) in Proposition 1 this implies that there is yet another solution P c ′′ of degree strictly less than D min . Hence, there are coefficients c ′′ j such that [START_REF] Norris | Optimal solutions to the L ∞ moment problem with lattice bounds[END_REF]. The vectors must therefore be linearly independent.

a 0 + c ′′ 1 a 1 + • • • + c ′′ d a d = 0, with d < D min , contradicting
Suppose D * ≥ D min is the highest degree of an existing solution. Since P c (z) is a solution of degree D min we get from (iii) in Proposition 1 that there are solutions of all intermediate degrees D min , . . . , D * . Hence, from (ii), a j ∈ V 0 j-1 for j = D min , . . . , D * and from (iii) a j ∈ V j-1 for j = D min +1, . . . , D * . Noting that if a j+1 ∈ V j then V j = V j+1 we can conclude inductively that V Dmin = • • • = V D * and a j ∈ V Dmin for j = D min , . . . , D * . We now have three different cases:

Suppose there are non-zero coefficients α k such that

α D * +1 a D * +1 + • • • + α nx a nx = 0,
and let k * be the highest index of all non-zero coefficients,

α k * = 0. Then a k * ∈ V 0 k * -1
and there is a solution of degree k * by (ii), a contradiction to the definition of D * . Hence, the vectors in (32) must be linearly independent and

D * = n x -dim V nx = D min + n x -rank A 1 = D max ,
showing (iv) for this case.

3. If D * < n x and D min > 0 we have

V nx = span{a 1 , . . . , a Dmin , a D * +1 , . . . , a nx }. ( 33
)
Suppose there are non-zero coefficients α k such that

α 1 a 1 + • • • + α Dmin a Dmin + • • • + α D * +1 a D * +1 + • • • + α nx a nx = 0.
Since a 1 , . . . , a Dmin are linearly independent at least one α k with k > D * must be non-zero. By the same argument as above in case two we then get a contradiction and the vectors in (33) must be linearly independent. Hence,

D * = D min + n x -dim V nx = D min + n x -rank A 1 = D max ,
showing this final case.

Proof of Theorem 2

To prove Theorem 2 we essentially have to combine the results from Propositions 1 and 3. The statement (i) is given directly by (i) in the latter. For the remaining points we have:

(ii) From (ii) in Proposition 3 we see that a 0 ∈ V d and a d ∈ V 0 d-1 . It follows from (iii) in Proposition 3 that d ≥ D min . On the other hand, if D max < n x and d > D max it says that V 0 d-1 = V d-1 . Hence, a d ∈ V d-1
which contradicts the linear independence of a Dmax , . . . , a nx established in point (iv) of Proposition 3.

(iii) We note that by [START_REF] Norris | Optimal solutions to the L ∞ moment problem with lattice bounds[END_REF] there are scalars v 1 , . . . , v Dmin such

a 0 = v 1 a 1 + • • • + v Dmin a Dmin . ( 34
)
Hence, a 0 ∈ V Dmin and since v Dmin = 0, we also have a Dmin ∈ V 0 Dmin . By (ii) in Proposition 3 there is thus a solution of degree D min which we denote (p * , q * ). Since a 1 , . . . , a Dmin are linearly independent by (iii) in Proposition 3, the coefficients in (34) are unique and therefore also the D min -degree solution is unique. Moreover, suppose that x j = y i = x * = 0 for some i, j. Then p * and q * would have a common factor (1zx * ), and by Lemma 3 also p(z) := p * (z)/(1-zx * ) and q(z) := q * (z)/(1-zx * ) would be a solution. But this is impossible since Deg(p) < Deg(p * ) = D min . By (iv), shown below, a solution is given by (p * r, q * r) where r(0) = 1 and Deg(r) = D max -D min . Hence n y ≥ Deg(q * r) = Deg(q * ) + n xrank A 1 . Suppose finally that D min < rank A 1 and that Deg(q * ) < n y -

n x + rank A 1 . Let Deg(r) = D max + 1 -D min . Then (p * r, q * r) is still a solution by Lemma 3 since (p * , q * ) is a solution, Deg(p * r) = D max + 1 = n x + D min + 1 -rank A 1 ≤ n x and Deg(q * r) < n y -n x + rank A 1 + D max + 1 -D min = n y + 1.
This contradicts (ii) and therefore Deg(q * ) = n y -n x +rank A 1 , concluding the proof of (iii).

(iv) We first note that there exists a solution of degree D max by Proposition 3 since if D max > D min we have a 0 ∈ V 0 Dmax-1 and a Dmax ∈ V Dmin = V Dmax-1 = V 0 Dmax-1 . Hence, (iii) in Proposition 1 shows that any polynomial pair of the stated type is a solution. On the other hand, if the polynomial solution is not of this type, then (ii) in Proposition 1 says there is a solution of degree strictly less than D min , contradicting (ii) above.

(v) We suppose first that A 1 is non-singular. Then rank A 1 = n x so that D min = D max and the uniqueness is given by (iii) above. If, on the contrary, A 1 is singular then D max > D min and since we can then pick infinitely many polynomials r(z) in (iv), we have infinitely many solutions.

(vi) This is a consequence of (iv). The solution can be represented by (p * r, q * r) for some polynomial r(z) with r(0) = 1. Let 1/x j for j = 1, . . . , D min and 1/y j for j = 1, . . . , Deg(q * ) be the roots of p * (z) and q * (z) respectively. Let 1/z j for j = 1, . . . , Deg(r) be the roots of r(z). Then

m k = Dmin j=1 x k j + Deg(r) j=1 z k j - Deg(q * ) j=1 y k j - Deg(r) j=1 z k j = Dmin j=1 x k j - Deg(q * ) j=1 y k j ,
which is independent of r(z) and uniquely determined because (p * , q * ) is unique.

Proof of Theorem 1

We can now use the results in Section 4 to prove Theorem 1.

(i-ii) To show the statements about Algorithms 1 and 2 we consider the reduced problem

m k = ñx j=1 xk j - ñy j=1 ỹk j , k = 1, . . . , K, (35)
where ñx = rank A 1 ≤ n x , ñy = n yn x + ñx ≤ n y and K = ñx + ñy ≤ K.

The moments m k in the left hand side are the same as in [START_REF] Curto | The truncated complex K-moment problem[END_REF]. First, we consider the minimal solution (p * , q * ) of [START_REF] Curto | The truncated complex K-moment problem[END_REF]. By (iv) in Proposition 3 we must have Deg(p * ) = D min ≤ rank A 1 = ñx . Moreover, by (iii) in Theorem 2, Deg(q * ) ≤ n yn x + rank A 1 = ñy .

It follows from Lemma 3 that (p * , q * ) is also a solution to (35). Second, let (p * , q * ) be the minimal degree solution to (35). Then by (iv) in Theorem 2 there is a polynomial r(z) with r(0) = 1 such that p * = p * r and q * = q * r. But then (p * , q * ) is also a solution to (5) by Lemma 3. By the uniqueness of the minimal degree solution of (5) it follows that r ≡ 1 and p * = p * ,q * = q * . Suppose now that there is another polynomial r(z) with r(0) = 1, Deg(r) > 0 such that (p * r, q * r) is a solution to (35). Then Deg(p * r) = D min + Deg(r) ≤ ñx = rank A 1 . Hence, D min < rank A 1 and therefore by (iii) in Theorem 2 we have Deg(q * ) = n y -n x +rank A 1 = ñy .

Thus, Deg(q * r) > ñy which is impossible if (p * r, q * r) is a solution. Hence, (p * , q * ) is the unique solution to (35) and therefore Ã1 is non-singular by (v) in Theorem 2.

Since Ã1 is invertible, the generalized eigenvalue problem (8) and c ′ are well-defined. Moreover, we can construct Ã-1 1 Ã0 . By (9), Ã-1

1 Ã0 =         -c 1 1 0 • • • 0 -c 2 0 1 • • • 0 . . . . . . . . . . . . . . . -c ñx-1 0 0 . . . 1 -c ñx 0 0 • • • 0        
, which is a companion matrix. It is well-known that for those matrices the elements in the first column are the coefficients of its characteristic polynomial. This is shown as follows: let M ij be the minor of V := zI -Ã-1 1 Ã0 , i.e. the determinant of the matrix obtained by removing row i and column j. Then, the determinant can be expanded by minors, for any j,

det(V ) = ñx i=1 (-1) i+j v ij M ij , V = {v ij }. Taking j = 1, we get M i,1 = det(diag(z, . . . , z, -1, . . . , -1)) with i -1 occurrences of -1, so that M i,1 = z ñx-i (-1) i-1 . Therefore, det(V ) = (-1) 2 (c 1 + z)M 1,1 + ñx i=2 (-1) i+1 c i M i,1 = c 1 z ñx-1 + z ñx + ñx i=2 c i z ñx-i = P (z),
which is exactly [START_REF] Gosse | Using K-branch entropy solutions for multivalued geometric optics computations[END_REF]. This shows that the results of Algorithms 1 and 2 are identical, since the generalized eigenvalues in (8) are exactly the roots of P (z).

It remains to show what the roots are. Let à = [ã 0 Ã1 ] be the A-matrix related to (35). Clearly, c = (1, c ′T ) T is in the null-space of à and hence P c (z) is the unique solution to (35). But for z = 0,

P (z) = c ñx + c ñx-1 z + • • • + c 1 z ñx-1 + z ñx = z ñx c ñx z ñx + c ñx-1 z ñx-1 + • • • + c 1 z + 1 = z ñx P c (1/z) = z ñx (1 -x 1 /z)(1 -x 2 /z) • • • (1 -x Dmin /z) = z ñx-Dmin (z -x 1 )(z -x 2 ) • • • (z -x Dmin ),
which extends to z = 0 by continuity. This concludes the proof of points (i) and (ii).

(iii) Let (p, q) be a polynomial solution to (5) and c the corresponding coefficient solution. From Lemma 2 we have

q(z) = p(z)e m(z) ,
where m(z) is defined in [START_REF] Simon | The classical moment problem as a self-adjoint finite difference operator[END_REF]. For the (K + 1)-th Taylor coefficient of the left and right hand side we have by Lemma 2 and Lemma 1,

0 = nx j=0 c j a K+1-j ⇒ a K+1 = - nx j=1 a K+1-j c j , (36)
since the k-th Taylor coefficient of q and p is zero for k > n x and k > n y respectively. Finally, the last row of (7) extended to size K + 1 gives

m K+1 = (K + 1)a K+1 - K j=1 m j a K+1-j .
Together the last two equations show point (iii).

6 Properties of A 1 and Markov's Theorem

We now look more in detail on the structure of the A 1 matrix. In particular we look at the implications of A 1 R being positive definite. Then we get an explicit simplified formula for the matrix and our results also shed some light on the relationship of our results to the classical Markov theorem on the existence and uniqueness of solutions to the finite moment problem (1) discussed in the introduction. For this we need to define the matrix

R =   1 . . . 1   ,
and note that left (right) multiplication by R reverses the order of rows (columns) of a matrix. In our notation we can then formulate Markov's theorem as follows Theorem 3 (Markov) Suppose K = 2n is even and n = n x = n y . There is a unique piecewise continuous function f (x) satisfying

m k = k R x k-1 f (x)dx, 0 ≤ f ≤ 1, k = 1, . . . , K, (37)
if A 1 R is symmetric positive definite and the matrix

a 0 A 1 a K+1 a T 0 (38)
is singular. This f is of the form in [START_REF] Akhiezer | Some questions in the theory of moments[END_REF][START_REF] Brenier | Équations de moment et conditions d'entropie pour des modèles cinétiques (French) [Moment equations and entropy conditions for kinetic models[END_REF].

Remark 7

The theorem does not rule out other forms of f (x) a priori, and without the second condition in (38) such solutions are indeed possible. It only considers the case n x = n y , i.e. problem (4), and says nothing about the possibility of other solution types, e.g. when the {x j } and {y j } are not interlaced as in [START_REF] Brenier | Équations de moment et conditions d'entropie pour des modèles cinétiques (French) [Moment equations and entropy conditions for kinetic models[END_REF].

We start by introducing some new notation that will be used throughout this section. If {x j } and {y j } is a solution of ( 5) and (p, q) is the corresponding polynomial solution as defined in [START_REF] Gosse | Two Moment Systems for Computing Multiphase Semiclassical Limits of the Schrödinger Equation[END_REF][START_REF] Gosse | Finite moment problems and applications to multiphase computations in geometric optics[END_REF], we can introduce the new polynomials p r (z) = z nx p(1/z) and q r (z) = z ny q(1/z) to describe the solution. Defining them by continuity at z = 0, we have

p r (z) = (z -x 1 ) • • • (z -x nx ), q r (z) = (z -y 1 ) • • • (z -y ny ). (39)
Furthermore, we assume that the number of distinct roots of p r (x j -branch values) is ñ. We also order the roots such that we can write

p r (z) = (z -x 1 ) 1+η1 (z -x 2 ) 1+η2 • • • (z -x ñ) 1+ηñ
, where 1 + η j is the multiplicity of the root x j , so that

n x = Deg(p r ) = ñ + ñ ℓ=1 η ℓ .
We start the analyis with a Lemma giving explicit expressions for the a k values.

Lemma 4 For k ≥ 0,

a ny-nx+1+k = ñ j=1 1 η! lim z→xj d ηj dz ηj (z -x j ) 1+ηj z k q r (z) p r (z) . (40)
Proof: This result follows from an application of the residue theorem in complex analysis as follows. Let C r be the circle in the complex plane with radius r. Since the roots of p(z) are non-zero, the function q/p is analytic within and on C ε if ε is taken small enough, and the Cauchy integral formula gives

a k = 1 k! d k dz k q(z) p(z) z=0 , k ≥ 0, 0, k < 0, = 1 2πi Cε q(z) p(z)z k+1 dz. Setting f (z) := q r (z) p r (z) = z ny-nx q(1/z) p(1/z) . (41)
and changing variable z → 1/z we get

a ny-nx+1+k = 1 2πi Cε q(z) p(z)z ny-nx+k+2 dz = 1 2πi Cε f (1/z) z k+2 dz = 1 2πi C 1/ε z k f (z)dz.
Hence, a ny-nx+1+k is given by the sum of the residues of z k f (z) (assuming we take small enough ε). By (41) and the restriction k ≥ 0 we see that its poles are located at the x j -values and they have multiplicities 1 + η j at x j . Then (40) follows from the residue formula for a pole of a function g(z) at z * with multiplicity η + 1,

Res(g, z * ) = 1 η! lim z→z * d η dz η (z -z * ) 1+η g(z).
When the branch values {x j } are distinct the expression for the a k elements simplifies. They can then be expressed as sums of the powers of {x j } in a way similar to the moments m k , but with weights different from one. We can also give a more concise description of the matrices A 0 and A 1 , which can be factorized into a product of Vandermonde and diagonal matrices. More precisely, we let V be the Vandermonde matrix

V =        1 1 • • • 1 x 1 x 2 • • • x nx x 2 1 x 2 2 • • • x 2 nx . . . • • • • • • . . . x nx-1 1 x nx-1 2 • • • x nx-1 nx       
, and introduce the diagonal matrices,

W =    w 1 . . . w nx   , X =    x 1 . . . x nx   ,
where w j are the weights defined as

w j = q r (x j ) p ′ r (x j ) . (42)
(Note that p r has only simple roots when {x j } are distinct, so p ′ r (x j ) = 0.) Then we can show Proposition 4 If {x j } are distinct, then for k ≥ 0,

a ny-nx+1+k = nx j=1 w j x k j , (43)
and

A 1 R = V W V T , A 0 R = V W XV T . (44)
Proof: When {x j } are distinct η j = 0 for all j and the expression (40) for the x j -residue simplifies,

lim z→xj (z -x j )z k q r (z) p r (z) = x k j q r (x j ) p ′ r (x j )
. This shows (43). For (44) we set b k = a ny-nx+1+k . Then

A 1-r R =      b r b r+1 . . . b r+nx b r+1 b r+2 . . . b r+nx+1 . . . . . . . . . . . . b r+nx b r+nx+1 . . . b r+2nx      ∈ R nx×nx , r = 0, 1.
From (43) we then have, for k ≥ 0,

     b k b k+1 . . . b k+nx      = nx j=1 w j      x k j x k+1 j . . . x k+nx j      = nx j=1 w j x k j      1 x j . . . x nx j      = V      w 1 x k 1 w 2 x k 2 . . . w nx x k nx      = V W      x k 1 x k 2 . . . x k nx      . Consequently, A 1-r R = V W      x r 1 x r+1 1 . . . x r+nx 1 x r 2 x r+1 2 . . . x r+nx 2 . . . . . . . . . . . . x r nx x r+1 nx . . . x r+nx nx      = V W X r V T ,
which concludes the proof. We now consider the implications of a positive definite A 1 R. It turns out that this is a necessary and sufficient condition to guarantee both distinct {x j } values and positive weights. We get Theorem 4 The matrix A 1 R is symmetric positive definite if and only if {x j } are distinct and the weights are strictly positive, w j > 0 for j = 1, . . . , n x .

Proof: We use the same notation as in Lemma 4 and set

S j (z) = 1 η j ! (z -x j ) 1+ηj q r (z) p r (z) .
We note that S j (z) is smooth and regular close to z = x j . Then by Lemma 4, for k ≥ 0,

a ny-nx+1+k = ñ j=1 lim z→xj d ηj dz ηj z k S j (z).
Next, we let v = (v 1 , . . . , v nx ) T be an arbitrary vector in R nx and recall that P v (z) is the corresponding n x -1 degree polynomial

P v (z) = v 1 + v 2 z + • • • + v nx z nx-1 . Then v T A 1 Rv = nx j=1 nx k=1 v j v k a ny-nx+j+k-1 = nx j=1 nx k=1 ñ ℓ=1 lim z→x ℓ d η ℓ dz η ℓ z j+k-2 S ℓ (z)v j v k = ñ ℓ=1 lim z→x ℓ d η ℓ dz η ℓ S ℓ (z) nx j=1 nx k=1 z j+k-2 v j v k = ñ ℓ=1 lim z→x ℓ d η ℓ dz η ℓ S ℓ (z)P v (z) 2 . ( 45
) If ñ + ñ j=1 ⌊η j /2⌋ ≤ n x -1, (46)
we can take

P v (z) = (z -x 1 ) 1+η1 (z -x 2 ) 1+η2 • • • (z -x ñ) 1+ηñ , ηj = ⌊η j /2⌋. Since 2(1 + ηℓ ) = 2 + 2⌊η ℓ /2⌋ ≥ 2 + 2(η ℓ /2 -1) > η ℓ and d ℓ dz ℓ f (z)(z -z * ) k z=z * = 0, 0 ≤ ℓ < k,
for all smooth enough f (z), we get v T A 1 Rv = 0, which contradicts the positivity of

A 1 R. Hence, ñ + ñ j=1 ⌊η j /2⌋ > n x -1 = ñ + ñ ℓ=1 η ℓ -1.
Since for any integer n > 0 we have ⌊n/2⌋ ≤ n -1 it follows that all η ℓ = 0 and ñ = n x . Hence, if A 1 R is positive definite, then {x j } are distinct.

To show the theorem it is now enough to show that, when {x j } are distinct, A 1 R is positive if and only if the weights are positive. From (45) we then have

v T A 1 Rv = nx ℓ=1 S ℓ (x ℓ )P v (x ℓ ) 2 = nx ℓ=1 w ℓ P v (x ℓ ) 2 .
Clearly, when all w ℓ > 0, this expression is positive for v = 0, and A 1 R is positive definite. To show the converse, we take P v (z) to be the Lagrange basis polynomials L j (z) of degree n x -1 defined as

L j (x i ) = 1, i = j, 0, i = j. . If A 1 R is positive then 0 < v T A 1 Rv = nx ℓ=1 w ℓ L j (x ℓ ) 2 = w j .
This can be done for each j, which concludes the proof.

We can now relate our conclusions with those in Markov's Theorem 3. We consider all solutions to (5), instead of those given by the integral relation (37) with a piecewise continuous function f (x). The extra condition (38) is then automatically satisfied, and we note that the positivity of A 1 R guarantees a unique solution also in our space of density functions [START_REF] Czarkowski | Solving the optimal PWM problem for single-phase inverters[END_REF]. We view this as a corollary of Theorems 2 and 4.

Corollary 1 If there exists a solution to [START_REF] Curto | The truncated complex K-moment problem[END_REF], then the matrix in (38) is singular. When n x = n y there is a unique solution to (5) of the form (3) if and only if A 1 R is symmetric positive definite.

Proof: We start by proving the singularity of (38). By (ii) in Proposition 2 a coefficient solution c = (c 0 , . . . , c nx ) T = (c 0 , cT ) T satisfies Ac = 0. Since A = (a 0 A 1 ) it remains to prove that c 0 a K+1 + a T 0 c = 0. This was already proved in (36).

Next, we prove the "if" part of the second statement. If A 1 R is symmetric positive definite it is non-singular and by (i), (iii) and (v) in Theorem 2 the minimal degree solution exists and is unique and x j = y i for all i, j. (If x j = 0 for some j, then there is no zero y i -value since Deg(q * ) = n by point (iii).) By Theorem 4 the corresponding branch values {x j } are distinct. It remains to show that, upon some reordering, the {x j } and {y j } are interlaced as in [START_REF] Brenier | Équations de moment et conditions d'entropie pour des modèles cinétiques (French) [Moment equations and entropy conditions for kinetic models[END_REF].

Order the x j -values in an increasing sequence and let m k be the number of y j -values such that y j < x k . Clearly, m k is increasing and 0 ≤ m k ≤ n y . Moreover, sgn(q r (x k )) = (-1) ny-m k and since lim z→∞ p ′ r (z) > 0, we also have sgn(p ′ r (x k )) = (-1) nx-k . Hence, by also using the fact that n y = n x , sgn(w k ) = (-1) ny-m k +nx-k = (-1) m k +k .

We conclude that m k + k is even, which implies that m k is in fact strictly increasing. Then, for k = 1, . . . , n x -1, we have m k+1 ≥ m k + 1 and

n x ≥ m nx ≥ m k + n x -k ⇒ m k ≤ k.
Similarly, m k ≥ m 1 + k -1 ≥ k -1, so k -1 ≤ m k ≤ k, and therefore 2k -

1 ≤ m k + k ≤ 2k.
Finally, since m k + k is even we must have m k = k, which implies that the values are interlaced. We now consider the "only if" part. If there is a solution of the form (3), then the {x j }-values are obviously distinct and m k = k. By Proposition 4 the weights are then given by (42) and they are positive since, as above, sgn(w k ) = (-1) m k +k = 1. It follows from Theorem 4 that A 1 R is positive definite.

Outlook

Several interesting issues may be worth mentioning:

1. Computational complexity in a finite difference implementation: one can consult the article [START_REF] Gosse | Finite moment problems and applications to multiphase computations in geometric optics[END_REF] where practical implementation issues and several examples of increasing complexity have been addressed in the context of geometric optics problems. In particular, comparisons with Lagrangian (ray-tracing) solutions are shown.

2. Extension to higher dimensions: nothing seems to be done in this direction at the time being; see however the last sections of [START_REF] Putinar | A renormalized Riesz potential and applications[END_REF] and the routines based on complex variables in [START_REF] Golub | A stable numerical method for inverting shape from moments[END_REF][START_REF] Elad | Shape from moments-An estimation theory perspective[END_REF] for "shape from moments".

3. A very special case of the trigonometric moment problem can be solved by means of a slight variation of the algorithms presented here, in [START_REF] Gosse | Finite moment problems and applications to multiphase computations in geometric optics[END_REF] and in Section IV.A of [START_REF] Elad | Shape from moments-An estimation theory perspective[END_REF]. That is to say, one tries to invert the following set of equations:

Let us state that in case the n + 1 real frequencies λ j are known, the set of complex amplitudes µ j are found by solving a Vandermonde system:

     1 • • • 1 exp(iλ 0 ) • • • exp(iλ n ) . . . . . . exp(inλ 0 ) • • • exp(inλ n )           µ 0 µ 1 . . . µ n      =      m 0 m 1 . . . m n     
.

The frequencies can be found through a byproduct of [START_REF] Elad | Shape from moments-An estimation theory perspective[END_REF][START_REF] Gosse | Finite moment problems and applications to multiphase computations in geometric optics[END_REF] as we state now: let us suppose n is odd (i.e. the number of equations is even), we form the two matrices,

A 1 =    m 0 • • • m n-1 2 . . . . . . m n-1 2 • • • m n-1    , A 2 =    m 1 • • • m n+1 2 . . . . . . m n+1 2 • • • m n    ,
and then the frequencies can be obtained through a generalized eigenvalue problem, A 1 v j = λ j A 2 v j , j = 0, ..., n. This kind of algorithm can be used to check the accuracy of the classical FFT and will be studied in a forthcoming article.

(

  iii) If Deg(p) ≤ Deg(p), any polynomial pair (pr, qr) is a solution if r(z) is a polynomial satisfying r(0) = 1 and Deg(r) ≤ Deg(p) -Deg(p). In particular, if Deg(p) ≤ m ≤ Deg(p) there is a solution (p, q) with Deg(p) = m.

Proposition 2

 2 The pair c = (c 0 , . . . , c nx ) T ∈ R nx+1 and d = (d 0 , . . . , d ny ) T ∈ R ny+1 is a coefficient solution to (5) if and only if(i) c 0 = 1,(ii) c is in the null-space of A, (iii)

1 .

 1 If D * = n x then V Dmin = V nx and by (20) we get D * = rank A 1 -D min +D max = dim V nx -D min +D max = dim V Dmin -D min +D max = D max since either a 1 , . . . , a Dmin are linearly independent or D min = 0 and V Dmin = ∅. This shows (iv) for D * = n x . 2. If D * < n x and D min = 0 then V Dmin = V D * = ∅ and V nx = span{a D * +1 , . . . , a nx }.

µ

  j exp(ikλ j ) = m k , k = 0, ..., n.

  Compute the rank of A 1 . Let ñx = rank A 1 and ñy = n yn x + ñx .

				a ny	a ny-1	. . . a ny-nx+1	
				a ny+1 . . .	a ny . . .	. . . a ny-nx+2 . . . . . .	   	.
			a ny+nx-1 a ny+nx-2 . . .	a ny
	3. Construct the matrices Ã0 , Ã1 ∈ R ñx×ñx as
				a ñy+1	a ñy	. . . a ñy-ñx+2	
	Ã0 =	   	a ñy+2 . . .	a ñy+1 . . .	. . . a ñy-ñx+3 . . . . . .	   	,
				a ñy+ñx a ñy+ñx-1 . . .	a ñy+1
				a ñy	a ñy-1	. . .	a ñy-ñx+1	
	Ã1 =	   	a ñy+1 . . .	a ñy . . .	. . . a ñy-ñx+1+1 . . . . . .	    .
		a ñy+ñx-1 a ñy+ñx-2 . . .	a ñy
	4. Solve the generalized eigenvalue problem


Guy-Michel Cloarec
email: cloarecg@eeng.dcu.ie.

Senior Member, IEEE John Ringwood
email: ring-wood@eeng.dcu.ie.

Ing Dip

Dr G-M Cloarec

Incorporation of Statistical Methods in Multi-step Neural Network Prediction Models

Keywords: Multi-step Ahead Prediction, Network Committees, Topology Optimization, Statistical Methods, Bootstrap, Principal Components Analysis

This paper addresses the problems associated with multi-step ahead prediction neural networks models. We will see how some concepts from the statistical theory field can be applied in various ways to improve the modelling. The generalization and error autocorrelation problems will be addressed using topological and methodological approach among which network committees, statistical bootstrap and principal component analysis will play a key role. These methods will be applied to the sunspot time series e

I. INTRODUCTION

T HE determination of the optimal structure in neural networks modelling is a critical issue can cannot be adressed analytically.

When some a priori knowledge or understanding of the problem is available or when the complexity of the function to map is low, it is possible to find a topology which is close to optimal. When the complexity of the problem increases, we have to use all the available methods, including techniques from the statistical theory field [l], [2].

Time series modelling is a class of problem where the goal is to approximate an hypothetical function describing the behaviour of a dynamical system on the basis of observations.

The modelling task is then complicated by the ignorance of the optimal input space and all the model uncertainties tend to downgrade the prediction accuracy. Moreover, in applications such as predictive control that require multi-step ahead predictions, and therefore the use of past predicted values in the model, the optimality of the topology and of the network parameters are essential to ensur reliability.

It is well known that modelling methods are problem dependant and that any modelling methodology cannot pretend to be completly universal.

However, it is possible to draw some general lines on the modelling methodology in the case of multi-step ahead prediction of time series with neural networks.

The purpose of this paper will therefore be to demonstrate how the use of statistical methods can improve the performance of neural networks models on critical times series problem. The benchmark time series expressed on Figure 1 [3] examined here will be the sunspot series well known for its pseudochaotic dynamics and probably non-stationary behaviour.

The paper is organized as follows : Section II expresses the modelling problems associated with sunspot prediction and gives a review of the previous studies. Section III will aim at explaining the rational of the different techniques and methods used proposed our the modelling. The results of our modelling and the comparisons with other studies will be done in section IV while section V will conclude this study.

II. SUNSPOT SERIES FORECASTING

The sunspot series is composed of the yearly mean of an arbitrarily defined sunspot number RI that expresses the number of spots or group of spots on the surface of the sun [3]. Although such phenomena have been observed since 1700, there is still no physical explanation.

The resulting time series is thought to be chaotic with a pseudoperiod of 12 years with some dynamical system behaviour. In order to enable comparisons with previous studies, the data set obtained from the Sunspot Index Data Center [3] that covers the yearly means from 1700 to 1996, was split into three subsets expressed on Figure 1. The first one called training set covers values from 1700 to 1920. The test set covers the 1921-1955 period and the validation set . Since the number of samples available is critical for statistical coherence, used methods that enabled to use larger sets that will be refered as the extended truining set covering the 1700-1955 period and the extended validation set covering the 1956-1996 period.

The modelling performed uses a feed-forward neural network that performs the mapping of the input space, the twelve previous sunspot numbers, to the output space i.e. the next sunspot number.

III. A METHODS IN MULTI-STEP NEURAL NETWORKS PREDICTION MODELING

A. Input Space Orthogonulixation and Principal Component Analysis

The determination of the optimal topology is one of the most critical issues in neural network modelling [2], [9]. Only a minimal structure will enable good prediction and generalization properties.

Despite attempts to estimate the problem complexity using Vapnik-Chervonenkis dimensions bounds [15], there is no theoretical result which gives the characteristics of the optimal neural network. The curse of dimensionality [2] expresses the problem of an exponentially growing network as the input space grows. The first step in the modelling is therefore to optimize the input space. This is performed by the orthogonalization of the input space using the principal components.

In this methodology new inputs are selected using linear systems model selection methods the We proportion of variimplemented linear models using principal components as inputs and statistical decision methods [ll] to determine the best principal components in respect to prediction error. The rational for using orthogonalized inputs was to feed the network with the minimal information, or in other words to reduce the hypervolume of the input space and therefore the network size. Moreover, the use of uncorrelated inpu ts improves convergence of the gradient based learning algorithms .

the The orthogonalization was performed on data series made of the past 12 sunspot numbers [5]. The three first principal components, the main extracted features, can be associated respectively with a two or three years delayed signal and a nine years delayed signal, with the third component corresponding to the so-called super-period [3] that could also be interpreted as a signal expressing the nonstationarity nature of the system (see Figure 2). To determine the optimal topology, one should use the PCA on a trained network which is sufficiently large to perform off-line neuron pruning. Altough the term suficiently large can be estimated using Vapnik-Chervonenkis dimensions [15], we used the results of the input space orthogonalization and of a linear modelling exercice (the complexity of the function is proportional to the number of parameters required in the linear model). Initially, a 6-6-l MLP neural network was trained to map the transformed and reduced input space to the output nals of the input and hidden layers space. were an

The output sigalyzed using the principal component analysis method and possible redundant neurons were detected.

The process was reiterated 10 times on different training subsets in order to achieve a good statistical inference on the topology determination and neuron pruning.

The analysis indicated that a 4-3-1 structure would reduce correlations between layers and make the networks close to optimally regularized.

One of the problems main in neural network modelling is the overfitting problem [9] that leads to poor generalization. Assuming that generalization is possible, the problem is then to determine a reliable generalization error estimate that will be used

C. Network Committees

The concept of network committees or ensembles of networks [a], [22] comes from the field of statistical theory and is a natural extension of the estimates accuracy improvement through statistical data analysis. The commitees are collections of finite numbers of neural networks trained on different sample sets, with different initial conditions or learning criteria. Their used is known to improve the generalization performance by decreasing the effect of spatial and temporal dependencies.

Once trained, the network committee predictions can be used to improve the decision making process during the model structural determination.

The predictions can also be combined to produce a better prediction.

The network commitees were used during the network topology determination in conjunction with the prinicpal component analysis. By enabling statistical inference, the model selection process was made on the basis of confidence limits.

The committees were also used to perform improved prediction. While the committees are usually combined using linear combiners [22], [2], we choosed to use a neural network, designated m the combining network. The overall structure of the predictor is expressed on Figure 4 This network was designed to combine the multi-step ahead predictions of the committees to give a better prediction. Multistep ahead predictors are known to be very unstable and very sensitive to initial errors so that a neural network must be used to overcome the error propagation problem. The networks were trained using different initial conditions on bootstrap ensembles originated from the same limited training set. It is important to obtain a population of networks with spread characteristics so that the performance ana1ysis and the prediction combination are improved.

The need for spread characteristics is also the rational for the use of the early stop training methods which gives emphasis on generalization rather than structural reduction. A drawback of spread characteristics [22] is that a portion of the committee downgrades the combination per- formance and so must be eliminated through a statistical process.

The PCA method was used to determine the initial topology of the committee of networks as well as the combination network. It has been shown in [20] that the more the networks are regularised, the better are the generalization estimates with resulting improvements with the early stop training method. The networks were trained, neurons were pruned using PCA results and the weights were subsequently put through an elimination process. The resulting corn .mittee of networks was regularized but with tics different enough to be combined efficiently.

character-

D. Bias/Variance Dilemma and Trajectory Learning

The Bias/Variance Dilemma [l], [2] refers to the trade off between the prediction error bias and variance that the training algorithm must perform to give satisfactory generalization properties. This is usually expressed by the MSE criterion that can be decomposed into bias and variance components.

The characteristics of the problem to solve may require that control of that trade off be exercised.

When committees are used, the trained networks exhibit different, characteristics with regard to prediction error bias and variance. It is therefore possible to use that distribution in the combination process to influence the trade off. Moreover, it is possible to refine or to regularize further these networks with bias or variance oriented criteria.

The trade off between prediction error bias and variance is a critical issue when multi-step ahead forecasting is involved. It is well known that when past predictions are used to estimate the future values, the resulting prediction error is autocorrelated and the predictor may become unstable.

A way to improve the predictor robustness and which enables a longer prediction horizon to be achieved is to emphasize the reduction of the error bias compared to the variance. However, the determination of such a trade off is difficult to estimate.

One method of performing the optimal trade off is to traine the combination network with predicted inputs. The recursive nature of the learning process requires non conventional training algorithms gathered under the trajectory learning acronym [9].

The trajectory learning algorithm is difficult to set up and because of the autocorrelated error is prone to instability. We propose here a new method that has the advantage of performing trajectory learning but still uses classical training algorithms.

The key idea is to use multi step ahead predictions obtained by network committees trained with single step ahead prediction error criterion and optimize for error bias and variance, and then to combine them so that the overall multi step ahead prediction is improved. The combination network which has the multi step ahead predictions as inputs is trained using a classical learning algorithm. In other words, we use the statistical properties of the network commitees to improve the prediction and then increase the reliable horizon of prediction To be efficient, the method requires good single-step ahead predictors so that their multi-step ahead predictions enable an efficient recombination.

This was obtain by assessing the generalization properties of the networks. The robustness of the predictors was estimated for horizons of prediction varying between one and five years.

The horizon of prediction is determined on the basis of the prediction errors obtained by the committees of networks optimized for single-step ahead prediction and pruned in regard of the bias or variance criteria. Figure 7 expresses the prediction error obtained for horizon of prediction varying from one to seven years and obtained on the extended training data set.

When the horizon of prediction has been defined, it is possible to implement the weight elimination process. In order to increase the spread, i.e. the variations of the network population, a simple elimination without retraining algorithm supervised by prediction error bias and variance is performed on the network committee. The network committee is divided in two groups composed respectively of the networks that minimize the error bias and variance. One by one, the weights are eliminated from the network and the multi-step prediction is analysed with respect to error and variance criteria. If a small increase or a decrease of the prediction error criterion is observed, the weight is pruned.

This simple algorithm, although far from being optimal, enables a drastic reduction of the number of parameters. Figure 5 shows a network on which a weight elimination otpimized for minimizing the error variance has been per-

2516

formed. The width of the weigths is proportional to their amplitude rank between the two layers. The size of the neurons is propotional to the rankink of the absolute sum of weights leaveing the neurons (neurons are ranked within layers). When color graphs are enabled, the color of the neurons depends on the rank of the absolute sum of the weights which fires the neurons. The use of such representations enables the viszlalization the network behaviour, to gives information on the critica paths within a network, as well as indications for supervised regularization.

The modeling procedure is outlined in Figure 9.

IV. RESULTS

Sunspot series forecasting is a benchmark problem for nonlinear modelling [4], [16] and neural network modelling in particular [5], [6], [7]. One difficulty is the very nature of the problem : an astronomical phenomena for which only a limited number of observations (or reconstructions) are available, which is particularly small with regard to the length pseudoperiods observed, and the lack of physical understanding of the choatic sun behaviour.

[4] performed its modeling in 1980 and following studies [51, ET1 evaluated their methods using the same limited sets. For example, both the test and validation sets were based on two pseudoperiods which, considering the erratic behaviour of the series, do not guarantee the significance of the performance estimates.

With more data available and by using methods (such as statistical bootstrap) that enable to use larger data sets and better performance estimates, we performed the modelling using the data from 1700 to 1955 as the training set and 1956 to 1997 as the validation set.

However, to enable comparisons with previous studies, the modelling was first performed using the same data set. The results results are gathered on Table I, with the double figures indicating mean and standard deviation respectively. When performed on the extended training set, the singlestep ahead predictions obtained on the extended validation sets are more representative of the modelling performance. The predictions and the actual sunspot numbers, along with the 95 % confidence limits, are expressed on Figure 6.

Obtaining satisfactory single-step prediction was only the first step of the multi-step ahead modelling process. Using the committee, we have performed mutli-step prediction with various horizon of prediction.

The statistics Figure 7 shows that beyond four-step ahead horizon, the prediction of the individual network committee members are not reliable.

A combination network, that aimed at combining the four-step ahead prediction of a selected subcommittee of networks was now trained. The subcommittee members were selected in relation to their four-step ahead prediction error bias and variance. The final committee of networks was composed of 16 networks optimal for error variance and 8 networks optimal for error bias.

The prediction normalized MSE obtained on the extended training sets was 0.3119 while we obtained 0.3361 on the extended validation set for H=4. The four-steps ahead prediction is shown in Figure 8.

It is difficult to assess the method based on the sunspot series application since the non-stationarities make the evaluation not necessarily representative on small data sets.

V. CONCLUSIONS

A general methodology for multi-step neural network modelling has been presented with an emphasis on the various concepts and techniques related to the field of statistical theory.

The transformation of the input space has been performed to optimize the neural network topology and an early stop training method using the bootstrap en- sembles, has been presented. Furthermore, network committees combined with a weight elimination method were used to deal with the bias/variance dilemma and to implement a new approach for trajectory learning. These proved to be satisfactory for the sunspot prediction. This modelling methodology was applied to the sunspot series, given its benchmark status, but the very limited data sets and the non-stationarity of the system made modelling and model estimation more critical. This was partly overcame by the statistical methods which enabled us to . use larger effective sets for both training and validation.

  Fig. 1. Sunspot Numbers Variations.Modelling of the sunspot series must observe several constrains. The first is the very limited number of data with regard to the 12 years pseudoperiod.Most of the statistical analysis and identification methods rely on statistical coherence of the data set, i.e. data sets on which it is possible to carry out generalization assessments.Another problem associated the limited number of data and the lack of understanding of the physical phenomena is the question of the time invariance of the sun considered as a dynamical system. The reduced observation time scale, with regard to astronomical time scales and the important variance of the pseudoperiodic signals both in amplitude and frequency (the cycle varies between 6 and 14 years), increases the uncertainty on the nature of the system studied[7]. The small data set leads us to consider the full set of reconstructed

  [lo], [If], by the normalization of the transformed input space 1121, [13] and by the selection of appropriate neural networks topology based on heuristic understanding of the problem [14]. Principal Component Analysis (PCA) is a multivariate analysis method that aims at determining correlation relationships between different variables [12], [13]. In neural network time series modelling context, these variables can be either the tapped delayed inputs or the output signals from each network layer [ 141. The PCA method can be used to transform the data so they have maximum variations and are orthogonal. The crosscorrelation criteria are then used to evaluate the level of information contained within the transformed signals (the principal components scores), so that a dimensionality reduction can be performed. The orthogonalized signals can be interpreted as extracted feutures in multivariate data analysis, where the dimensionality reduction can be called compression, for example in the communications field [12], [13]. The first six principal components were selected, although only the first three explained 95 % of the data variance. The sunspot series and the first three components scores are expressed on Figure 2. The choice of the principal ante components is critical since is not a causality criterion.

Fig. 2 .

 2 Fig. 2. Sunspot Series and Principal Components.

  for model selection and optimization. Classical methods use single-sample statistics (AIC, FPE, MDL criteria) [16], [ll] that estimate the model performance on a single set (usually the training set), or the split-sample method that estimates the error on a test set which is independent of the training set. Other methods such as cross-validation and bootstrap [17], [US], use statistical theory principles to obtain better estimates of the error without using test sets (which is an advantage with limited data sets problems). In k-folded cross validation, the training set is split into k distinct equal size sets and k-l models are estimated using training sets composed of k-l subsets where the one left out each time is different. The bootstrap method is used to decrease the statistical dependencies within subsets. The method consists to draw with replacement samples at random to form a given size training subsets [20]. For cross validation and bootstrap methods, the optimality hides in so called 7 ratio problem [19] which relates to the choice of the number of subsets and the size of the subsets. A greater number of subsets improves the statistical inference, but may result with unacceptable small bootstrap set sizes, which do not adequately represent the coherent characteristics of the data. Small subsets may emphasize idiosyncrasies which increase the variance of the estimate error. We arbitrarily defined 36 samples (3 pseudoperiods) to be the minimal subset size although previous studies [4], [5] used smaller sets (23) for test and validation sets. The statistical bootstrap will give us statistics that will be used for inference estimation which is useful for model selection and generalization error estimation. The training was performed on a 4-3-l network using (successively) 3 bootstraps of 120 samples (y = 0.5) and at each epoch the generalization error is estimated using 15 bootstraps of 36 samples each. Because of the hazards due to random initial conditions, the training was performed on an ensemble of 20 networks. The backpropagation algorithm was used during the first 10 epochs and was subsequently replaced by the Levenberg-Marquardt algorithm supervised by an early stop criterion. Remedies to poor generalization are numerous and are known as regularization methods [9]. The early stop method, the oldest form of regularization, aims at the overfitting problem [20] and gives good results when good generalization estimates are available. It has been shown in [20] that the use of statistical bootstrap estimates of the generalization property are a viable means for implementing an early stop rule.

Figure 3

 3 Figure 3 shows the variations of the model response MSE on the extended training and validation sets as well as the bootstrap generalization estimates when the network was trained to convergence.The early stop criterion was implemented on the basis of the minima of the generalization

  Fig. 4. Multi-step Prediction Neural Network Structure.

  Fig. 5. Weight Elimination Optimized for Variance Criterion.

Fig. 6 .

 6 Fig. 6. One Step Ahead Forecast on the Validation Set.

Fig. 7 .

 7 Fig. 7. Prediction Error Statistics on Extended Sets.

  Fig. 8. Four-step Ahead Forecast on the Validation Set.

TABLE I COMPARISONS OF THE NORMALIZED ERRORFOR SINGLE-STEP AHEAD FORECASTING.

 I

 , vol. 3, pp 193-209, 1990. [6] , vol. 8, pp 190-196, 1995.