Numerical aspects of finite Hausdorff moment problem by maximum entropy approach

Applied Mathematics and Computation 118 (2001) 133±149 www.elsevier.com/locate/amc

Numerical aspects of ®nite Hausdor€ moment problem by maximum entropy approach Aldo Tagliani Faculty of Economics, Trento University, 38100 Trento, Italy

Abstract Some numerical aspects of the ®nite Hausdor€ moment problem in the framework of maximum entropy have been reviewed, as stability of Lagrange's multipliers computation, entropy decreasing, averages estimate. Spectral properties of Hankel matrices have been considered, determinant and conditioning number have been estimated. Finally entropy estimate has been provided. Ó 2001 Elsevier Science Inc. All rights reserved. Keywords: Hausdor€ moment problem; Moments space; Maximum entropy; Hankel matrix; Conditioning number

1. Introduction The ®nite Hausdor€ moment problem consists in giving a sequence flk g; k ˆ 0; . . . ; M of power moments and looking for a positive real valued function f …x† : ‰0; 1Š ! R‡ satisfying the constraint of moments Z 1 xk f …x† dx; k ˆ 0; . . . ; M …1:1† lk ˆ 0

and l0 ˆ 1, without loss of generality. Clearly, there exists then an in®nite variety of functions whose M ‡ 1 moments coincide and a unique reconstruction of f …x† is impossible. In this case the ®nite Hausdor€ moment problem is a typical underdetermined inverse problem. The classical maximum E-mail address: [email protected] (A. Tagliani). 0096-3003/01/$ - see front matter Ó 2001 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 9 9 ) 0 0 2 1 0 - 6

134

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

entropy method [1] treats the function f …x† as a probability density and its approximation fM …x† is determined maximizing the entropy H ‰f Š, where Z 1 f …x† ln f …x† dx …1:2† H ‰f Š ˆ ÿ 0

subject to constraint (1.1). This is a standard constrained optimization problem, leading to the following explicit representation of the maximum entropy (ME) approximate [2]: ! M X j kj x ; …1:3† fM …x† ˆ exp ÿ jˆ0

to be supplemented by the condition that its ®rst M ‡ 1 moments be given by lk , k ˆ 0; . . . ; M Z 1 xk fM …x† dx; k ˆ 0; . . . ; M; …1:4† lk ˆ 0

where k0 ; . . . ; kM are the Lagrange's multipliers. We obtain also the maximum entropy value attained with the maximum entropy density and given by [2]: H ‰fM Š ˆ

M X jˆ0

kj lj :

…1:5†

Eq. (1.4) constitutes a nonlinear system of M ‡ 1 equations for M ‡ 1 unknown multipliers k0 ; . . . ; kM . When the multipliers obtained from (1.4) are substituted into (1.3) we get explicitly the approximation fM …x† of the true probability density f …x†. Several analytical and computational aspects of the ®nite Hausdor€ moment problem in the approach of ME have been investigated in the literature, as 1. solution's existence [3], 2. convergence [4], 3. stability Refs. [5±7]. In this paper we reconsider the stability problem in recovering fM …x†. Such a problem is a consequence of the ill-posedness of the in®nite-dimensional Hausdor€ moment problem [6]. Thus any accurate computational approach of k's multipliers results to be severely ill-conditioned. Nevertheless, as well documented in the literature, the spread of applications of ME technique pushes us to believe that interesting quantities can be accurately computed in spite of ill-conditioning of Lagrange's multipliers. In other R 1 terms, for instance, we are interested in an accurate evaluation of averages 0 F …x†fM …x† dx, where F …x† is some known function of physical interest, rather than an accurate evaluation of fM …x† itself.

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

135

In general, given …l0 ; . . . ; lM † to which corresponds the true fM …x†, we obtain a computed fMcalc …x† having lj ‡ Dlj ; j ˆ 0; . . . ; M, jDlj j  1. So we investigate how the computed fMcalc …x† a€ects the evaluation of some physical interesting quantities, as averages or entropy, when the true fM …x† is substituted by fMcalc …x†. In a past paper [7] we considered the same problem, varying uniquely lM and ®xing …l0 ; . . . ; lMÿ1 †: pessimistic conclusions about entropy and fM …x† estimate have been reached. In the present paper all the moments lj , j ˆ 0; . . . ; M will vary and di€erent conclusions will be obtained. Some background is in order. 1. The moment space DM  R‡ whose points are the n-ple l ˆ …l1 ; . . . ; lM † is the convex hull of the curve fx; x2 ; . . . ; xM g, x 2 ‰0; 1Š [8]. The existence conditions of fM …x† involve DM , more precisely if (i) the point l ˆ …l1 ; . . . ; lM † is outside DM , the corresponding ®nite Hausdor€ moment problem does not admit any solution, (ii) l 2 oDM (oDM is the boundary of DM ) the only distribution having …l1 ; . . . ; lM † as ®rst moments is a (uniquely determined) convex combination of Dirac's delta, (iii) l is in the interior of DM , there are in®nitely many distributions, one of them being fM …x† [3]. 2. The previous existence conditions can be algebrically formulated involving de®nite positive Hankel matrices, we bound ourselves to quote only one set of them [8]: 2 3 l0


lM   l0 l1 6 .. 7: .. …1:6† ; . . . ; Dl2M ˆ 4 ... Dl0 ˆ l0 ; Dl2 ˆ . 5 . l1 l2 lM


l2M Then the ®nite Hausdor€ moment problem solvability is related to Hankel determinants jDl2j j, j ˆ 0; 1; . . . ; M positivity. 3. Once …l1 ; . . . ; lMÿ1 † 2 DMÿ1 are assigned then the extremes Mth moment lM consistent with the ®rst M ÿ 1 components …l1 ; . . . ; lMÿ1 † are indicated by ‡ M and lÿ M and lM , respectively. Such values are obtained as lM belongs to oD M provide the M-width of D [8]: ÿ ÿ2M‡2 l‡ M ÿ lM 6 2

…1:7†

(such a relation can be obtained more naturally resorting to canonical mo‡ ÿ ments pk ˆ …lk ÿ lÿ k †=…lk ÿ lk †, k ˆ 1; . . . ; M [9]). The further relation will be used [8]: l2M ÿ lÿ 2M ˆ

jDl2…Mÿ1† j : jDl2M j

…1:8†

4. Varying only one moment li , i ˆ 0; . . . ; M, while the remaining ones are held ®xed, so that kj ˆ kj …li †, j ˆ 0; . . . ; M, di€erentiation of (1.4) provides

136

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

2

Dl2M

3 dk0 =dli 6 7 ..
4 5 ˆ ÿei‡1 ; . dkM =dli

…1:9†

where ei‡1 is the canonical unit vector 2 RM‡1 . 5. From (1.4) (with k ˆ 0) k0 can be obtained as function of k1 ; . . . ; kM : " Z ! # M 1 X 1 j exp ÿ kj x dx : …1:10† k0 ˆ ln l0 0 jˆ1 k0 is a convex function of k1 ; . . . ; kM [2] as well as H ‰fM Š provided by (1.5). For numerical as well as theoretical purposes, the Lagrange's multipliers are calculated minimizing the potential function " Z ! # M M 1 X X 1 j exp ÿ kj x dx ‡ kj l j …1:11† H ‰fM Š ˆ l0 ln l0 0 jˆ1 jˆ1 as setting …oH ‰fM Š†=okj ˆ 0, j ˆ 1; . . . ; M then (1.4) is satis®ed.

2. Stability of fM …x† Let fM …x; kj † be the ME solution corresponding to the set …l0 ; . . . ; lM † and fM …x; kj ‡ Dkj † the ME solution corresponding to …l0 ‡ Dl0 ; . . . ; lM ‡ DlM †, jDlj j  1; 8j. Let us consider the relative error on fM …x† generated by the variation Dkj on kj given by a variation of Dlj on lj , j ˆ 0; . . . ; M: ‰fM …x; Dlj †Š ˆ:

fM …x; kj ‡ Dkj † ÿ fM …x; kj † : fM …x; kj †

…2:1†

Taking into account (1.3), (1.7) and (1.9) and from Taylor expansion we have ‰fM …x; Dlj †Š ˆ exp ˆ exp

ÿ

M X

Dkj xj

ÿ1

jˆ0

ÿ

M X jˆ0

ˆ exp

!

ÿ

M X iˆ0

M X Dkj x Dl Dli i iˆ0

!

j

M X Dkj j Dli x Dli jˆ0

ÿ1 ! ÿ1

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

l0 B .. B . B B B liÿ1 B BX B M Dli 1 B ˆ exp B B iˆ0 jDl2M j li‡1 B B B . B .. B @ l M 0 1 x Dl0 l0 l1 1 ' .. .. jDl2M j .. . . . Dl l l 0

M

M

l1






.. .

.. .

li






x






li‡2






.. .

.. .

lM‡1











xM






lM

.. .

.. .






l2M

M‡1

137

1 lM C .. C C . C C liÿ1‡M C C C M C x C ÿ 1 C li‡1‡M C C C C .. C . A l2M : …2:2†

Adopting a di€erent hypothesis on Dlj then from (2.2) we have (i) Dlj =lj ˆ Dl ˆ constant 8j ‰fM …x; Dlj †Š ˆ

Dl
jDl2M j ˆ Dl jDl2M j

…2:3†

and then Z

1 0

2

2 ‰fM …x; Dlj †ŠfM …x† dx ˆ …Dl† :

(ii) Dlj =nj ˆ Dl ˆ constant 8j, 0 < n 6 1 so 0 1 x 1 l l 0 1 Dl n l l 1 2 ‰fM …x; Dlj †Š ˆ ÿ jDl2M j . .. .. . . . . M n lM lM‡1 and then [10] Z 1 ‰fM …x; Dlj †ŠfM …x† dx ˆ Dl ÿ 0

independent on M.

…2:4† that








.. .




xM lM lM‡1 .. . l 2M

…2:5†

138

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

(iii) Dlj ˆ Dl ˆ constant 8j, then we have Z 1 ‰fM …x; Dlj †ŠfM …x† dx ˆ Dl ÿ 0

…2:6†

independent on M. Eqs. (2.3)±(2.6) lead to optimistic conclusions about fM …x† recovering. (iv) DlM ˆ Dl, Dlj ˆ 0; j < M, then from (2.2) we have ‰fM …x; Dlj †Š ˆ

PM …x† Dl; jDl2M j

…2:7†

where

1 l0 PM …x† ˆ . .. l

Mÿ1

x l1 .. .

lM

xM lM .. .







.. .




l2Mÿ1



…2:8†

is the M-degree orthogonal polynomial de®ned by the ®rst 2M moments of fM …x† [8]. Taking into account (1.7), (1.8) and [8], the following relationship holds:  2 Z 1 Z 1 Dl 2  ‰fM …x; Dlj †ŠfM …x† dx ˆ PM2 …x†fM …x† dx jDl2M j 0 0  2 Dl jDl2M j
jDl2…Mÿ1† j ˆ jDl2M j ˆ

…Dl†2 …Dl†2 …Dl†2 > P ÿ l2M ÿ lÿ l‡ 2ÿ4M‡2 2M 2M ÿ l2M 2

ˆ …Dl† 24Mÿ2 :

…2:9†

Eq. (2.9) raises pessimistic conclusions about fM …x† recovering as M assumes high values. Eq. (2.9) can be interpreted if (1.7) is considered. As only lM varies then the corresponding point l ˆ …l1 ; . . . ; lMÿ1 ; lM ‡ DlM † moves towards the moment space boundary so that the density fM …x† tends to degenerate in a set of Dirac's delta.

3. Entropy decreasing Let us consider the entropy H ‰lj Š corresponding to …l0 ; . . . ; lM † and H ‰lj ‡ Dlj Š corresponding to …l0 ‡ Dl0 ; . . . ; lM ‡ DlM †. Taking into account (1.9) we have

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

139

DH ˆ H ‰lj ‡ Dlj Š ÿ H ‰lj Š M X ‰…kj ‡ Dkj †…lj ‡ Dlj † ÿ kj lj Š ˆ jˆ0

ˆ

M X jˆ0

ˆ

M X

…kj Dlj ‡ lj Dkj ‡ Dlj Dkj † "

jˆ0

M X Dkj kj Dlj ‡ …lj ‡ Dlj † Dl Dli i iˆ0

#

0 M M X X Dkj 1 Dl0 kj Dlj ‡ lj Dl ÿ ˆ jDl2M j ... Dli i jˆ0 iˆ0 Dl M 0 l


l 0 M Dl M l0


lM X 0 1 kj Dlj ÿ ˆ .. .. jDl2M j ... .


. jˆ0 Dl l


l M M 2M 0 Dl0


DlM Dl l0


lM 0 1 ÿ .. .. jDl2M j ... .


. Dl lM


l2M M M X Dl0 jDl2M j kj Dlj ÿ ˆ jDl 2M j jˆ0 0 Dl0


DlM


lM 1 Dl0 l0 : ÿ .. .. jDl2M j ... . . . . . Dl l


l "

#

M

M

Dl0 l0 .. . lM







.. .




DlM lM .. . l2M

…3:1†

2M

From (3.1) we have (i) Dlj =lj ˆ Dl ˆ constant 8j DH ˆ Dl

M X jˆ0

kj lj ÿ l0 Dl ÿ

ˆ Dl…H ÿ l0 † ÿ l0 …Dl†

1 …Dl†2 l0 jDl2M j jDl2M j 2

…3:2†

so that DH  Dl:

…3:3†

140

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

(ii) Dlj =nj ˆ Dl ˆ constant 8j, 0 < n 6 1, then from (1.3) and [10] we have 0 1 n


nM


lM 1 l0 l1 2 M X j …Dl† n l l2


lM‡1 1 kj n ÿ l0 Dl ÿ DH ˆ Dl 2jDl2M j . . jˆ0 .. . . . . . .


.. nM l l


l M M‡1 2M M X Pj2 …n† ˆ ÿ Dl…ÿ ln fM …n† ÿ l0 † ÿ …Dl†2 jˆ0

2

ˆ ÿ Dl…ÿ ln fM …n† ÿ l0 † ÿ

…Dl† qM …n†

…3:4†

so that DH  Dl:

…3:5†

Here PM …n† represents the M-degree orthogonal polynomial previously introduced, qM …n† the maximum mass which can be concentrated at the point x ˆ n, under the condition that ®rst M ‡ 1 moments of fM …x† be assigned [10]. Eqs. (3.3) and (3.5) lead to conclude that the entropy estimate is reliable. (iii) DlM 6ˆ 0, and Dlj ˆ 0; j < M. Let us consider …l0 ; . . . ; lMÿ1 †, the corresponding ME density fMÿ1 …x†, its Mth computed moment lM and entropy H ‰fMÿ1 Š. Then, given …l0 ; . . . ; lMÿ1 ; lM † we consider the corresponding ME density fM …x†. Such a density is characterized by kM ˆ 0 and H ‰fM Š ˆ H ‰fMÿ1 Š. From (1.9), with i ˆ M, we have M X dH ‰fM Š dkj ˆ lj ‡ kM ˆ kM ˆ 0; …3:6† dl dlM lM ˆl M jˆ0 M

jDl2…Mÿ1† j d2 H ‰fM Š dkM : ˆ ˆÿ dlM jDl2M j dl2M

…3:7†

From Taylor expansion at lM and taking into account (3.6), (3.7), (1.7) and (1.8) we have ÿDH ˆ ÿ H ‰fM …lM ‡ DlM †Š ‡ H ‰fM …lM †Š ˆ

jDl2…Mÿ1† j 1 …DlM †2 2 jDl2M j

ˆ

1 …DlM † 1 …DlM † > P 24Mÿ3 …DlM †2 : ‡ ÿ † …l ÿ l † 2 …l2M ÿ lÿ 2 2M 2M 2M

2

2

Here Dl2…Mÿ1† and Dl2M are evaluated substituting lM by lM .

…3:8†

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

141

4. Averages estimate Rather than f …x† itself, one is most often interested in averages of some interesting physical bounded functions F …x†. Let us consider two sets of moments …l0 ; . . . ; lM † and …l0 ‡ Dl0 ; . . . ; lM ‡ DlM † to which the densities corand the corresponding respond to fM …x; kRj † and fM …x; kj ‡ Dkj †, respectively, R1 1 averages hFM i ˆ 0 F …x†fM …x; kj † dx, hFMcalc i ˆ 0 F …x†fM …x; kj ‡ Dkj † dx. We have DFM ˆ hFM i ÿ hFMcalc i Z 1 F …x†‰fM …x; kj † ÿ fM …x; kj ‡ Dkj †Š dx ˆ 0

ˆ ÿ

Z

1 0

F …x†‰fM …x; Dlj †ŠfM …x; kj † dx

…4:1†

(i) if Dlj =lj ˆ Dl ˆ constant 8j then from (2.3) ‰fM …x; Dlj †Š ˆ Dl and then Z 1 …4:2† F …x†fM …x† dx 6 KjDlj jDFM j ˆ jDlj
0

K being a constant so that jF …x†j 6 K. Thus the averages can be reliably estimated whatever M. 5. Lagrange's multipliers computation Taking into account (2.3)±(2.6), (3.3), (3.5) and (4.2) we reach the conclusion that an accurate Lagrange multipliers computation is not mandatory. Rather it is preferable to have an accurate estimate of fM …x† so that its moments be affected by small relative or absolute errors. Thus the method we will propose to calculate Lagrange's multipliers should substitute minimization of (1.9) or the solution of the non-linear system (1.4). Let us vary all the moments …l0 ; . . . ; lM †. Taking into account (1.9), the corresponding k's variation is M X Dkj Dl Dli i iˆ0 l0 1 . ˆÿ . jDl2M j . l M

Dkj ˆ




.. .




ljÿ1 .. . ljÿ1‡M

Dl0 .. . DlM

lj‡1 .. .

lj‡1‡M




.. .




lM .. : . l2M

…5:1†

(i) if Dlj =lj ˆ Dl ˆ constant 8j then Dk0 ˆ ÿDl;

Dkj ˆ 0; 1 6 j 6 M

…5:2†

142

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

holds. Under the hypothesis Dlj =lj ˆ Dl ˆ constant 8j a twofold goal is reached 1. entropy and averages are reliably calculated, 2. Eq. (5.2) leads to the conjecture that even kj 's computation is accurate. Thus a viable route for kj 's computation should be addressed to satisfy Dlj =lj ˆ Dl ˆ constant, j ˆ 0; . . . ; M, for instance in the sense of least squares. 6. Spectral properties of Hankel's matrices We state some results concerning spectral properties of Dl2M matrices (1.6) arising in the Hausdor€ ®nite moment problem in the framework of ME. Thus our results concern only particular Hankel matrices. Preliminarily we quote some known results. (i) For any real positive de®nite order n Hankel matrix Mn : Cond2 Mn P 3
2nÿ6

…6:1†

holds [11]. R1 (ii) Let Mn ˆ ‰mi‡jÿ2 Š ˆ 0 xi‡jÿ2 dw…x†, i; j ˆ 1; . . . ; n generated by a positive, non-decreasing, non-constant function of bounded variation on the interval ‰0; 1Š w…x†. Then the asymptotic estimate lim inf…Cond2 Mn †1=n P 16

…6:2†

kmin 6 2ÿ4n‡2

…6:3†

n!1

and hold [12], where kmin denotes the minimum Mn eigenvalue. (iii) Let x…x† be a positive integrable function on ‰0; 1Š, x…x† P n > 0, lj , j P 0 its moments and Mn the corresponding Hankel matrix, then [13] lim …Cond2 Mn †1=n ˆ lim …Cond2 Hn †1=n ˆ e3:525... ;

n!1

n!1

…6:4†

where Hn denotes n order Hilbert matrix. Eq. (6.4) states essentially each Hankel matrix, generated by a strictly positive weight x…x†, conditioned like the Hilbert matrix. Preliminarily we provide an estimate of Hankel matrix determinant. Taking into account (1.9) (with i ˆ M), (1.7) and (1.8) we have ÿ

dkM jDl2…Mÿ1† j 1 1 ˆ ˆ > ‡ P 24Mÿ2 ÿ dlM l2M ÿ l2M l2M ÿ lÿ jDl2M j 2M

so that from jDl2…Mÿ1† j > 24Mÿ2 jDl2M j we obtain recursively 1 ˆ jDl0 j > 22 jDl2 j > 22
26 jDl4 j >


> 22
26 . . . 24Mÿ2 jDl2M j

…6:5†

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

143

and then jDl2M j < 2ÿ2M

2

…6:6†

The maximum value of Hankel determinant jDl2M j over the moment space D2M is obtained resorting to canonical moments pk , k ˆ 1; . . . ; 2M [9]. The moments sequence maximizing jDl2M j is generated by a measure having M ‡ 1 equal weights 1=…M ‡ 1† located at the zeroes of the polynomial x…1 ÿ x†L0M …x†, where LM …x† is the Mth Legendre polynomial on the interval ‰0; 1Š. Such a sequence of canonical moments is given by [9]: 1 p2jÿ1 ˆ ; 2

p2j ˆ

M ÿj‡1 ; 2…M ÿ j† ‡ 1

j ˆ 1; . . . ; M:

Let us express (1.4) in vector form /…k† ˆ l;

…6:7†

where l ˆ …l0 ; . . . ; lM †, k ˆ …k0 ; . . . ; kM †. Then the …l; k† entry of the Jacobian matrix J /…k† is the moment ll‡k of fM …x†. Being J /…k† ˆ Dl2M symmetric de®nite positive then the map /…k† ˆ l is globally invertible with regularity [14], so that k ˆ /ÿ1 …l† is meaningful. Let us vary all the moments lj , j ˆ 0; . . . ; M into lj ‡ Dlj . Next consider k ˆ /ÿ1 …l† and k ‡ Dk ˆ /ÿ1 …l ‡ Dl†. We have Dk ˆ /ÿ1 …l ‡ Dl† ÿ /ÿ1 …l† from which kDkk2 ˆ kJ /ÿ1 …l†Dlk2 6 kJ /ÿ1 …l†k2
kDlk2 :

…6:8†

Being J /ÿ1 …k† symmetric de®nite positive, from (6.8) we obtain kDkk2 6 q…Dlÿ1 2M †
kDlk2 ˆ

1 kDlk2 ; kmin …Dl2M †

…6:9†

where q…
† denotes the spectral radius. Varying uniquely lM into lM ‡ DlM , while lj , j < M are kept constant and taking into account (6.5) we have jDkM j > 24Mÿ2 jDlM j:

…6:10†

Combining (6.9) and (6.10) we obtain jDlM j
24Mÿ2 < jDkM j < kDkk2 6 ˆ

1 kDlk2 kmin …Dl2M †

1 jDlM j kmin …Dl2M †

…6:11†

from which the estimate of kmin …Dl2M †, according to [12]: kmin …Dl2M † 6 2ÿ4M‡2 :

…6:12†

So we obtain 8M an estimate of conditioning number of …M ‡ 1†-order Hankel matrix Dl2M :

144

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

Cond2 …Dl2M † ˆ q…Dl2M †
q…Dlÿ1 2M † ˆ

kmax …Dl2M † > kDl2M k2 24Mÿ2 > 24Mÿ2
maxfli‡jÿ2 g i;j kmin …Dl2M †

ˆ 24Mÿ2 l0 ˆ 24Mÿ2

…6:13†

to be compared to asymptotic estimate (6.2). 6.1. A special case: Hilbert matrices As lj ˆ 1=…j ‡ 1†, j ˆ 0; . . . ; M then ME density fM …x†  1 8M, as can be easily proved, so that Dl2M ˆ HM‡1 , where HM‡1 ˆ ‰hi;j Š ˆ 1=…i ‡ j ‡ 1†, i; j ˆ 0; . . . ; M denotes …M ‡ 1†-order Hilbert matrix. From (6.13): Cond2 …HM‡1 † > 24Mÿ2

…6:14†

holds. Eq. (6.14) can be improved, resorting to the following relationship characterizing Hilbert matrices [9]: jHM j 1 …2M†!…2M ‡ 1†! ˆ ˆ 4 jHM‡1 j l2M ÿ lÿ …M!† 2M

…6:15†

so that from (6.11), (6.13) and Stirling formula Cond2 …HM‡1 † >

…2M†!…2M ‡ 1†! 4

…M!†

'

24M …2M ‡ 1† : pM

…6:16†

8M be compared to asymptotic estimate Cond2 …HM‡1 † ' 24M =pM given in [12]. Eqs. (6.14) and (6.16) can be further improved as follows. Let CM ˆ ‰ci;j Š, i; j ˆ 0; . . . ; M ÿ 1 be the inverse of the lower triangular Cholesky factor of the T CM ˆ HMÿ1 , where superscript T stands for transpose and Hilbert matrix HM , CM ÿ1 ÿ1 HM ˆ ‰hi;j Š, i; j ˆ 0; . . . ; M ÿ 1. The entries ci;j of the lower triangular matrix are the coecients of shifted Lagrange polynomials Li …x† on ‰0; 1Š normalized so that kLi …x†k2 ˆ 2i ‡ 1, with positive leading coecient [6]: Li …x† ˆ

i X

ci;j xj ;

i ˆ 0; 1; . . .

…6:17†

jˆ0

with cM;j ˆ …2M ‡ 1†

1=2

   M M ‡j …ÿ1† ; j j j

From (1.9), with i ˆ M, we obtain

j ˆ 0; . . . ; M:

…6:18†

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

2 6 4

145

3

dk0 =dlM 7 .. ÿ1 eM‡1 5 ˆ ÿHM‡1 . dkM =dlM

…6:19†

and then

2 M  X dkj ÿ1 ÿ1 ˆ ‰HM‡1 eM‡1 ŠT
‰HM‡1 eM‡1 Š: dl M jˆ0

…6:20†

From (6.11), and (6.18)±(6.20): q…HM‡1 † kHM‡1 k2 maxi;j fhi;j g ˆ > ÿ1 ÿ1 † q…HM‡1 † kmin …HM‡1 † q…HM‡1 " 2 #1=2 M  X 1 jDkj2 dkj > ˆ ˆ dlM kmin …HM‡1 † kDlM k jˆ0 " #1=2 " #1=2 M  M 2 X X ˆ hÿ1 ˆ c2M;M c2M;j j;M

Cond2 …HM‡1 † ˆ

jˆ0

 ˆ …2M ‡ 1†

jˆ0

2M

" X M 

M

M

jˆ0

j

from which the estimate 

2M Cond2 …HM‡1 † > …2M ‡ 1† M

" X M  jˆ0

M j



M ‡j

2 #1=2 …6:21†

j 

M ‡j j

2 #1=2 :

…6:22†

Eq. (6.22) improves the previous estimates (6.14) and (6.16). 7. Hankel matrices preconditioning Eq. (6.4) has been obtained under the hypothesis moment vector l ˆ …l0 ; . . . ; l2M † characterizing Hankel matrix entries Dl2M be far from moment space boundary (equivalently x…x† P n > 0). Then (6.4) suggests that Hankel matrices are conditioned essentially the same as Hilbert matrices. In [13] it has been stated that Hilbert matrix is a good preconditioner for Hankel matrix, more precisely T †6 Cond2 …CM‡1 Dl2M CM‡1

2M ‡ 1 ; n

…7:1†

T T CM ˆ HMÿ1 and CM‡1 Dl2M CM‡1 is where CM‡1 , previously de®ned, satis®es CM symmetric de®nite positive. We provide a qualitative result about conditioning

146

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

of Hankel matrices, involved in the Hausdor€ moment problem by the ME approach, dropping the restriction x…x†  fM …x† P n > 0. Let us consider modi®ed moments l ˆ …l0 ; . . . ; lM † ˆ CM‡1 lT , where l ˆ …l0 ; . . . ; lM † are the given power moments. Then, according to the ME technique, given the partial information l , the approximate density fM …x† is " # M X kj Lj …x† …7:2† fM …x† ˆ exp ÿ jˆ0

constrained by Z 1 Lk …x†fM …x† dx; k ˆ 0; . . . ; M: lk ˆ 0

…7:3†

Fixed all lj , j < M, while lM varies continuously, di€erentiating (7.3) we have 2 3 dk0 =dlM 6 7 .. T …7:4† CM‡1 Dl2M CM‡1 4 5 ˆ ÿeM‡1 . dkM =dlM

from which, taking into account CM‡1 is lower triangular, ÿ

ÿ1 j
jDl2…Mÿ1† j dkM jHM‡1 jHM j jDl2…Mÿ1† j
ˆ ˆ  ÿ1 dlM jHM‡1 j jDl2M j jHM j
jDl2M j

…7:5†

and from (6.15), (6.16), (1.7) and (1.8): ÿ

4 jDl2…Mÿ1† j dkM …M!† pM > 4M
24Mÿ2 ˆ  2 …2M ‡ 1† dlM …2M†!…2M ‡ 1†! jDl2M j p : ˆ 4…2 ‡ …1=M††

…7:6†

Resorting to modi®ed moments l ˆ CM‡1 lT is equivalent both using Legendre polynomials Lj …x† in place of monomials xj and the preconditioned T in place of Dl2M . It is well known that the expression of matrix CM‡1 Dl2M CM‡1 a polynomial in the monomial basis is far worse conditioned than its expansion in Legendre polynomials [15] and the use of orthogonal expansions in the numerical solution of the ®nite moment problem is a standard procedure [6]. Let us express (7.3) in vector form l ˆ /…k†, analogously (6.7). By a similar procedure we have k ˆ /ÿ1 …l †, and calling J /…k† the Jacobian matrix we have, similarly (6.8) and (6.9), jDkM j < kDkk2 6 kJ /ÿ1 …l †k2
jDlM j ˆ

1
jDlM j: T kmin …CM‡1 Dl2M CM‡1 †

…7:7†

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

Taking into account (7.6) we have 1 p > : T kmin …CM‡1 Dl2M CM‡1 † 4…2 ‡ …1=M††

147

…7:8†

T symmetric de®nite positive, to evaluate Being CM‡1 Dl2M CM‡1 T T …C Dl C †. Following Cond2 …CM‡1 Dl2M CM‡1 † it needs to estimate kmax PM M‡1 2M M‡1 [13], for each M-degree polynomial p…x† ˆ jˆ0 cj Lj …x† and letting c ˆ …c0 ; . . . ; cM † we have Z 1 T p2 …x†fM …x† dx ˆ cCM‡1 Dl2M CM‡1 cT ; …7:9† 0

Z

1 0

2

p …x†fM …x† dx 6 kp

2

2 …x†k1

Z

1 0

fM …x† dx

2

ˆ kp2 …x†k1 6 …2i ‡ 1†ccT

…7:10†

so that T cCM‡1 Dl2M CM‡1 cT 6 2M ‡ 1 ccT and then

…7:11†

T † 6 2M ‡ 1: kmax …CM‡1 Dl2M CM‡1

…7:12†

From (7.8) and (7.12) we can conclude simply T † > 1: Cond2 …CM‡1 Dl2M CM‡1

…7:13†

Comparing (6.13) and (7.13) we can argue that the Hilbert matrix HM‡1 should be a good preconditioner for the class of Hankel's matrices. 8. Entropy estimate We provide an estimate of H ‰fM Š given by (1.5) without explicitly calculating fM …x†. Let us consider the vector lH ˆ f1=…j ‡ 1†g, j ˆ 0; . . . ; M whose entries are the Hilbert matrix entries and l ˆ …l0 ; . . . ; lM † the vector of given moments. As l ˆ lH the Hankel matrix and Hilbert matrix coincide. Then the l corresponding ME density fMH …x†  1, 8M, so that kj ˆ 0, j ˆ 0; . . . ; M and lH H ‰fM Š ˆ 0, 8M. Let us consider the vector n ˆ …n0 ; . . . ; nM † belonging to the segment joining lH and l. Once given n we consider the corresponding ME density fMn …x† from which the next moments nj , j > M can be obtained. From (1.5) and (1.9) we have M oH ‰fM Š X okj ˆ lj ‡ kk ˆ kk ÿ d0k ; ol olk k jˆ0

k ˆ 0; . . . ; M;

…8:1†

148

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

where d0k denotes d-Kronecker. As l ˆ lH then from (8.1) we have oH ‰fM Š ˆ ÿd0k ; olk

k ˆ 0; . . . ; M:

From (8.1) and taking into account (1.9) o2 H ‰fM Š okj ˆ oli olk lˆn olk lˆn n0


nkÿ1 . . 1 .. . . ˆÿ . jDl2M jlˆn . . n


n M

kÿ1‡M

…8:2† the Hessian is obtained

.. . ek‡1 .. .

nk‡1 .. .




.. .

nM .. .

nk‡1‡M






n2M

From Taylor expansion at l ˆ n we obtain  M  X 1 oH ‰fM Š lj ÿ H ‰fM Š ˆ H ‰fM Šjn ‡ 1‡j olj n jˆ0 1 T ÿ …l ÿ lH †Dlÿ1 2M …l ÿ lH † 2 n 1 T ˆ ÿ …l ÿ lH †Dlÿ1 2M …l ÿ lH † < 0: 2 n

:

…8:3†

…8:4†

Being Dlÿ1 2M symmetric de®nite positive T …l ÿ lH †…Dlÿ1 1 1 2M †jn …l ÿ lH † 6 6 kl ÿ lH k2 kmax …Dl2M †jn kmin …Dl2M †jn

…8:5†

and from (8.4) kmin …Dl2M †jn 6

kl ÿ lH k2 6 kmax …Dl2M †jn : ÿ2H ‰fM Š

From the Gerschgorin theorem   M X 1 ˆ k maxfl; lH gk1 max lj ; kmax …Dl2M †jn 6 1‡j jˆ0

…8:6†

…8:7†

and from (8.6), (8.7) and ME principle H ‰f Š 6 H ‰fM Š 6 ÿ

kl ÿ lH k2 2k maxflj ; lH gk1

holds for each M, providing an estimate of H ‰f Š.

…8:8†

A. Tagliani / Appl. Math. Comput. 118 (2001) 133±149

149

References [1] E.T. Jaynes, Where do we stand on maximum entropy, in: R.D. Levine, M. Tribus (Eds.), The Maximum Entropy Formalism, MIT Press, Cambridge, MA, 1978, pp. 15±118. [2] H.K. Kesavan, J.N. Kapur, Entropy Optimization Principles with Applications, Academic Press, New York, 1992. [3] I. Csiszar, I-divergence geometry of probability distributions and minimization problems, Ann. Probab. 1 (1975) 146±158. [4] J.M. Borwein, A.S. Lewis, Convergence of best entropy estimates, SIAM J. Optim. 1 (1991) 191±205. [5] G. Inglese, A note about minimum relative entropy solutions of ®nite moment problems, Numer. Funct. Anal. Optim. 16 (1995) 1143±1153. [6] G. Talenti, Recovering a function from a ®nite number of moments, Inverse Problems 3 (1987) 501±517. [7] A. Tagliani, Hausdor€ moment problem and maximum entropy: a uni®ed approach, Appl. Math. Comput. 105 (1999) 291±305. [8] S. Karlin, L.S. Shapley, Geometry of Moment Spaces, AMS Memoirs, vol. 12, AMS, Providence, RI, 1953. [9] H. Dette, W.J. Studden, The Theory of Canonicals Moments with Applications in Statistics Probability and Analysis, Wiley/Interscience, New York, 1997. [10] N.I. Achieser, The Classical Moment Problem and Some Related Questions in Analysis, Hafner, New York, 1965. [11] E.E. Tyrtyshnikov, How bad are Hankel matrices?, Numer. Math. 67 (1994) 261±269. [12] J.M. Taylor, The condition of Gram matrices and related problems, Proc. Roy. Soc. Edinb. 80 (1978) 45±56. [13] D. Fasino, The spectral properties of Hankel matrices and numerical solution of ®nite moment problem, J. Comput. Appl. Math. 65 (1995) 145±155. [14] W.H. Fleming, Functions of Several Variables, Addison-Wesley, Reading, MA, 1965. [15] W. Gautschi, The condition of orthogonal polynomials, Math. Comput. 26 (1972) 923±924.