A unified regularization theory: the maximum non-extensive entropy principle

A unified regularization theory: the maximum non-extensive entropy principle

Haroldo F. de Campos Velho^I; Elcio H. Shiguemori^II; Fernando M. Ramos^I; João C. Carvalho^III

^IInstituto Nacional de Pesquisas Espaciais (INPE), Laboratório Associado de Computação e Matemática Aplicada (LAC), São José dos Campos, SP, Brazil

^IIInstituto de Estudos Avançados (IEAv), Comando Geral Brasileiro de Tecnologia Aeroespacial (CTA), São José dos Campos, SP, Brazil

^IIIAgência Nacional de Telecomunicações (ANATEL), Brasília, DF, Brazil

E-mails: [haroldo,elcio,fernando]@lac.inpe.br / [e-mail]

ABSTRACT

Tsallis' non-extensive entropy is used as a regularization operator. The parameter ''q'' (non-extensivity parameter) has a central role in the Tsallis' thermostatiscs formalism. Here, several values of q are investigated in inverse problems, using q < 1 and q > 1. Two standard regularization techniques are recovered for special q-values: (i) q = 2 is the well known Tikhonov regularization; (ii) q = 1 is the standard Boltzmann-Gibbs-Shannon formulation for entropy. The regularization feature is illustrated in an inverse test problem: the estimation of initial condition in heat conduction problem. Two methods are studied for determining the regularization parameter, the maximum curvature for the L-curve, and the Morozov's discrepancy principle. The new regularization of higher order is applied to the retrieval of the atmospheric vertical profile from satellite data.

Mathematical subject classification: 35R30.

Key words: Maximum non-extensive entropy pinciple, unified regularization theory, inverse problems.

1 Introduction

Inverse problems belong the class of ill-posed problems. In these type of problems existence, uniqueness and stability of their solutions cannot be ensured.

The regularization theory appeared in the 60's as a general method for solving inverse problems [1]. In this approach, the non-linear least square problem is associated with a regularization term (a priori or additional information), in order to obtain a well-posed problem. A well-known regularization technique was proposed by Tikhonov [1], where smoothness of the unknown function is searched. Similarly to Tikhonov's regularization, the maximum entropy formalism searchs for global regularity and yields the smoothest reconstructions which are consistent with the available data. The maximum entropy principle was first proposed as a general inference procedure by Jaynes [2] on the basis of Shannon's axiomatic characterization of the amount of information [3]. This principle has successfully been applied to a variety of fields. In the middles of 90's, higher order procedures for entropic regularization have beenintroduced [5, 6] (see also [7] and references inside).

In the macrocospic thermodynamics, the entropy is applied to quantify the irreversibility in systems, and the entropy has a fundamental role in the statistical physics. The Shannon's entropy is a measurement of the missing information. In the inverse problem context, it is a shape entropy.

A non-extensive formulation for the entropy has been proposed by Tsallis [8, 9]. Recently, the non-extensive entropic form (S_q) was used as a new regularization operator [10], using only q = 0.5. The q parameter plays a central role in the Tsallis' thermostatistics, in which q = 1 the Boltzmann-Gibbs-Shannon's entropy is recovered. As mentioned, the non-extensive entropy includes as a particular case the extensive entropy (q = 1): in which the maximum entropy principle can be used as a regularization method. However, another important particular case of the non-extensive entropy occurs when q = 2. In such case, the maximum non-extensive entropy principle to the S₂ regularization operator is equivalent to the standard Tikhonov regularization (or zeroth-order Tikhonov regularization) [1, 11].

This new regularization operator was tested for estimating initial condition in heat conduction problem [11, 12]. In the present study several values for q were used for the non-extensive entropic regularization term. Synthetic data with Gaussian white noise corruption were used to simulate experimental data. Two methods were investigated for determining the regularization parameter: the Morozov's discrepancy principle [13], and the maximum curvature scheme of the curve relating smoothness versus fidelity, inspired in Hansen's geometrical criterion [14].

Initial conditions for numerical weather prediction are obtained from the Earth observation system. Radiances measured by meteorological satellites is a component of this observation system. From these radiances, atmospheric temperature and humididy profiles can be determined. Several methodologies and models have been developed to compute this estimation. A method using the second order entropy approach [7] was proposed as a solution for this inverse problem. Here, the higher order non-extensive entropy will be used to improve the previous inverse solution.

Non-extensive entropy as a new regularization operator

A non-extensive form of entropy has been proposed by Tsallis [8] by the expression

where p_i is a probability, and q is a free parameter. In thermodynamics the parameter k is known as the Boltzmann's constant. Similarly as in the mathematical theory of information, k = 1 is considered in the regularization theory. Tsallis' entropy reduces to the the usual Boltzmann-Gibbs-Shanon formula

in the limit q ® 1.

As for extensive form of entropy, the equiprobability condition produces the maximum for the non-extensive entropy function, and this condition is expressed as

where, in the limit N_p ® ¥, S_q diverges if q < 1, and saturates at 1/(q – 1) if q > 1 [9].

The equiprobability condition leads the regularization function defined by the operator S_q(p), given by Eq. (1), to search the smoothest solution for the of the unknown vector p.

The parameter q has a central role in Tsallis' thermostatiscs, and it is called the non-extensivity parameter - see Eqs. (A.5) and (A.6) in the ^Appendix. Figure 1 shows the functional form for the Tsallis' entropy for several values of q. For q < 5/3, the standard central limit theorem applies, implying that if p_i is written as a sum of M random independent variables, in the limit case M ® ¥, the probability density function for p_i in the distribution space is the normal (Gaussian) distribution [9]. However, for 5/3 < q < 3 the Levy-Gnedenko's central limit theorem applies, resulting for M ® ¥ the Levy distribution as the probability density function for the random variable p_i. The index in such Levy distribution is g = (3 – q)/(q – 1) [9].

The non-extensive approach has been used in many different applications, such as in a certain type of anomalous diffusion process [9], as well as the statistical model for data from turbulent flow [15] and from financial market [16]. According to Plastino and Plastino [17], the first experimental confirmation of Tsallis' non-extensive formalism is the Boghosian's approach of the two dimensional pure electron plasma [18]. Some properties of the thermostatistics formalism are described in the ^Appendix.

Regularization unified

Some comments have been addressed for indicating the physical properties of the S_q operator. The goal of this Section is to describe formally the properties for this operator looking at regularization purposes. Regularization properties for entropy operator emerges from the Jaynes' inference criterium: the maximum entropy principle, where all events have the same propability to occur. Implying that all parameters assume the same value: p_i = 1/N_p. The following Lemma extends this result for non-extensive entropy.

Lemma. The non-extensive function S_q is maximum as p_i = 1/N_p for all i.

Proof. The problem is to find the maximum of the function (1), with the following constrain

since p_i represents a probability. Therefore, it is possible to define an objective function where the constrain can be added to the non-extensive function:

where l is the Langrange multiplier. The Lagrange multiplier, in this case, can be determined when a minimum for the objective function J(p) is found, as following

This result can be used to obtain the value of the p_i's that maximizes the function J(p):

that means if p_i = 1/N_p for all i = 1,...,N_p the non-extensive entropy function is maximum.

The next theorem shows that the extensive entropy and Thikhonov's regularizations are particular cases of the non-extensive entropy.

Theorem. For particular values for non-extensive entropy q = 1 and q = 2 are equivalents to the extensive entropy and Tikhonov regularizations, respectively.

Proof. (i) q = 1: Taking the limit,

(ii) q = 2: Remembering that max S₂ is equivalent to min(–S₂), yields

now, for the maximum (minimum) value holds Ñ_pS₂ = 0, therefore

In conclusion: max S₂(p) = min (the zeroth-order Tikhonov regularization).

Inverse analysis

Typically, inverse problems are ill-posed - existence, uniqueness and stability of their solutions cannot be ensured. An inverse solution can be formulated to obtain existence and uniqueness, but this solution can still be unstable under the presence of noise in the experimental data. Hence, it requires some regularization technique, i.e., the incorporation in the inversion procedure of some available information about the true solution. Following the Tikhonov's approach [1], a regularized solution is obtained by choosing the function p* that minimizes the following functional

where is the experimental data, F(p) is the answer computed from the forward model, W[p] denotes the regularization operator, a is the regularization parameter, and || · ||₂ is the L₂ norm.

The regularization parameter a is chosen by two methods: numerically, assuming that a bound d (or the 'statistics') of the measurement error is known, i.e., ||F_computed – || < d — this numerical procedure is based on Morozov's discrepancy principle [13]; graphically, finding out the point of maximum curvature in the curve , a type of L-curve [14, 19].

Optimization algorithm

The optimization problem is iteratively solved by the quasi-newtonian optimizer routine from the NAG Fortran Library [20]. This algorithm is designed to minimize an arbitrary smooth function subject to constraints (simple bound, linear or non-linear constraints), using a sequential programming method.

This routine has been successfully used in several previous works: in geophysics, hydrologic optics, and meteorology.

Backward heat conduction

An example of the use for this new generalized regularization is the identification of the intial condition in heat conduction transfer. The direct (forward) problem consists of a transient heat conduction problem in a slab with adiabatic boundary condition and initially at a temperature denoted by f(x). The mathematical formulation of this problem is given by the following heat equation

where T(x,t) (temperature), f(x) (initial condition), x (spatial variable) and t (time variable) are dimensionless quantities. The set of partial differential equations is solved by using a central finite difference approximation for space variable O(Dx²), and explit Euler method for numerical time integration O(Dt) [21].

This problem has been used for testing different methodologies in inverse problems [12, 11, 22, 23, 24], and it is badly conditioned problem [12].

The numerical experiment with the non-extensive entropy is based on two test functions, the triangular function

and semi-triangular function

The experimental data (measured temperatures at a time t > 0), which intrinsically contains errors, is obtained by adding a random perturbation to the exact solution of the direct problem, such that

where s is the standard deviation of the errors, and µ is a random variable taken from a Gaussian distribution, with zero mean and unitary variance. All tests were carried out using 5% of noise (s = 0.05).

It is important to observe that the spatial grid consists of 101 points (N_x = 101), and the time-integration is performed up to t = 0.01. The residue R(f_a) and the error E(f_a) are defined by

If we effectively want to apply some kind of regularization, which means a > 0 in Eq. (5), then the discrepancy principle - an a-posteriori parameter choice rule - implies that a suitable regularized solution can be obtained. Since the spatial resolution is N_x = 101, the optimum a is reached for R(f_*) ~ N_xs² = 0.2525. Table 1 shows the square diference term R(f_*) obtained for different values of a, and the optimum value is pointed out for each case in bold font.

A set of tables (Tables 1 and 2) presents the least squares (or residual) term R(f_a) and the error E(f_a) between the approximated (or calculated) solution f_a and the exact solution f_exact obtained for two values of q = 0.5 and q = 2.0, from a family of regularization non-extensive entropy functions S_q, for different values of a. The value of a satisfying the discrepancy principle is pointed out for each test functions by bold font. Regularized solutions are presented in Figures 2 and 3.

Figures 3a-3d

The parameter vector was always subjected to the following simple bounds: 1.2 > f_k > –0.2 for the triangular test function, and 1.2 > f_k > 0 for the semi-triangular test function, with k = 1, 2,..., N_x.

Another criterion for finding the regularization parameter was also investigated, and it is based on the maximum curvature in the L-curve [14]. Figures 4a-4b show the L-curve for triangular test function using q = 0.5 and 2.0, respectively. The L-curve for semi-triangular test function is displayed in Figures 5a-5b. The regularization parameter a is chosen at the corner of the L-curve.

The numerical values for the regularization parameters estimated by Morozov's and Hansen's criteria are shown in Table 3.

Reconstructions using a as computed by Hansen's criterion [14] are shown in Figures 6a-6d for the triangular test function, and Figures 7a-7d for the semi-triangular test function. The best reconstruction was obtained using q = 2.5, and the worst for q = 0.5.

Retrieval of atmospheric temperature profile

The vertical structure of temperature and water vapor plays an important role in the meteorological process of the atmosphere. However due a several conditions, there is a lack of observation in several regions of the Earth. In this sense, the retrieval of temperature and humidity profiles from satellite radiance data became important for applications such as weather analyses and data assimilation in numerical weather predictions models.

Interpretation of satellite radiances in terms of meteorological fields requires the inversion of the Radiative Transfer Equation (RTE) where measurements of radiation performed in different frequencies are related to the energy from different atmospheric regions. The degree of indetermination is associated with the spectral resolution and the number of spectral channels. Moreover, usually this solution is very unstable regarding the noises in the measuring process. Also, several methodologies and models have been developed to improve the satellite data processing. Due to the difficulty of obtaining correct RTE solutions, several approaches and methods were developed to extract information from satellite data [25, 26, 27, 28]. The direct problem may be expressed by [29]

where I_l is the spectral radiance, l is the channel frequency; Á is the layer tospace atmospheric transmittance function, the subscript s denotes surface [12]; and B is the Planck function which is a function of the temperature T (or pressure p):

being h the Planck constant, c the light speed, and k_B the Boltzmann constant. For practical purposes, equation (21) is discretized using central finite differences:

with I_iº I_li (0), N_l is the number of channels in the satellite, and N_p is the number of the atmospheric layers considered.

Some previous results have employed a generalization of the standard maximum entropy principle (MaxEnt-0) for solving this inverse problem [7, 30]: the higher order entropy approach. The same strategy can be applied here. Therefore, the non-extensive entropy of order-g is defined as

and

where g = 0, 1, 2,..., and D is a discrete difference operator. The standard MaxEnt-0 can be derived from (23) and (24) imposing g = 0 and q = 1.A small value should be added to the difference operator (say V = 10^–15) to assure a definite quantity for all values of q.

Figures 8-10 present the atmospheric temperature retrieval achieved using radiance data from the High Resolution Radiation Sounder (HIRS-2) of NOAA-14 satellite. HIRS-2 is one of the three sounding instruments of the TIROS Operational Vertical Sounder (TOVS). The results for the non-extensive MaxEnt-0 (for short: NE-MaxEnt-0) did not produce good results - Figure 8. However, results with NE-MaxEnt-1 and NE-MaxEnt-2 are compared to those obtained with MaxEnt-2 (such results have already been analysed against the profile computed by ITPP-5, a TOVS processing package [7] employed by weather service research centers throughout the world), and to in situ radiosonde measurements.

From figures 9 and 10 is hard to identify the better performance among the retrievals. Table 5 compares the RMS of NE-MaxEnt-2 and NE-MaxEnt-1 related to the MaxEnt-2. The layer defined by 500-250 hPa there is no significant diffences among the different regularization operators. But, for the layer 250-50 hPa, the best performance was obtained by the NE-MaxEnt-1, for the 700-500 hPa better result for the MaxEnt-0, and for the region 1000-750 hPa the NE-MaxEnt-2 (q = 0.5) has presented the better inversion.

Conclusion

The implict strategy and the regularization techniques adopted in this work yield good results in reconstructing the initial condition of the heat equation and the atmospheric temperature profile.

For the initial condition estimation in the heat transfer process, the Morozov's discrepancy principle was efficient to estimate the regularization parameter in the analyzed cases. The Hansen's criterion also produced good estimates for the regularization parameter a. According to the Tables 1 and 2, the value of regularization parameter calculated by the Morozov's principle tends to over-estimate the value of a. However, looking at Table 4, it is possible to realize that the parameter a computed by the Hansen's criterion is closer to the optimum. Nevertheless, sometimes is hard to obtain the L-curve. One case specially difficult was found to q = 2.5, for some values of a was not possible to obtain a solution (no convergence). A possible solution to convergence would be to change the deterministic optimizer by a stochastic one. Table 4 shows that an appropriated choice of the non-extentive parameter q can improve the reconstruction. Of course, other schemes for determining the regularization parameter can be used (see [31]).

The new regularization technique used in this work, namely the maximum non-extensive entropy, worked very well for the backwards heat equation for all parameter q tested (zeroth-order). The choice q = 0.5 made by Rebollo-Neira et al. [10] was linked to a previous result of physical relevance related to the relaxation of two-dimensional turbulence [18]. There is no reason to restrict the regularization operator S_q at q = 0.5. Actually, the worse reconstructions for the triangular test function were obtained using S_0.5!

Concerning to the atmospheric temperature retrieval, some improvement was obtained using the higher order of the non-extensive entropic approach. We note that the most important layer for numerical weather prediction lies into the levels 1000-700 hPa. The results suggest that we need to combine different techniques, considering different atmospheric layers, in order to have a better inverse solution.

Acknowledgment. Authors acknowledge the CNPq and FAPESP, Brazilian agencies for research support. H.F. de Campos Velho also thanks to the Prof. Dr. João Batista da Silva (UFPA, Brazil) and Dra. Valeria Cristina Barbosa (LNCC, Brazil) for the usefull discussions and support.

[20] E04UCF routine, NAG Fortran Library Mark 17, Oxford, UK (1995).

Received: 15/XII/06. Accepted:15/XII/06.

#697/06.

A1: Non-extensive entropy:

A2: q-expectation of an observable:

Properties

1. If q ® 1:

2. Non-extensive entropy is positive: S_q > 0.

3. Non-extensivity

4. Max S_q under constrain O_q = (canonical ensemble):

where the

_i is the energy of state i, O_q = U_q is the non-extensive form to the internal energy, and the normalization factor Z_q (partition function), for 1 < q < 3, is given by

For q = 1 yields

• [1] A.N. Tikhonov and V.I. Arsenin, Solutions of Ill-posed Problems. John Wiley & Sons, 1977.
• [2] E.T. Jaynes, Information theory and statistical mechanics. Phys. Rev., 106 (1057), 620.
• [3] C.D. Rodgers, Retrieval of the atmospheric temperature and composition from remote measurements of thermal radiation. Rev. Geophys. Space Phys., 14 (1976), 609-624.
• [4] C.R. Smith and W.T. Grandy (Eds.). Maximum-Entropy and Bayesian Methods in Inverse Problems, in Fundamental Theories of Physics, Reidel, Dordrecht (1985).
• [5] F.M. Ramos and H.F. de Campos Velho, Reconstruction of geoelectric conductivity distributions using a minimum first-order entropy technique, 2nd International Conference on Inverse Problems on Engineering, Le Croisic, France, vol. 2, 199 (1996).
• [6] H.F. de Campos Velho and F.M. Ramos, Numerical inversion of two-dimensional geoelectric conductivity distributions from eletromagnetic ground data. Braz. J. Geophys., 15 (1997), 133.
• [7] F.M. Ramos, H.F. de Campos Velho, J.C. Carvalho and N.J. Ferreira, Novel Approaches on Entropic Regularization. Inverse Problems, 15(5) (1999), 1139-1148.
• [8] C. Tsallis, Possible generalization of Boltzmann-Gibbs statistics. J. Statistical Physics, 52 (1988), 479.
• [9] C. Tsallis, Nonextensive statistics: theoretical, experimental and computational evidences and connections. Braz. J. Phys., 29 (1999), 1.
• [10] L. Rebollo-Neira, J. Fernandez-Rubio and A. Plastino, A non-extensive maximum entropy based regularization method for bad conditioned inverse problems. Phys. A, 261 (1998), 555.
• [11] W.B. Muniz, F.M. Ramos and H.F. de Campos Velho, Entropy- and Tikhonov-based regularization techniques applied to the backwards heat equation. Comp. Math. Appl., 40 (2000), 1071.
• [12] W.B. Muniz, H.F. de Campos Velho and F.M. Ramos, A comparison of some inverse methods for estimating the initial condition of the heat equation. J. Comp. Appl. Math., 103 (1999), 145.
• [13] V.A. Morosov, Methods for Solving Incorrectly Posed Problems Springer Verlag (1984).
• [14] P.C. Hansen, Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev., 34 (1992), 561.
• [15] F.M. Ramos, R.R. Rosa, C. Rodrigues Neto, M.J.A. Bolzan, L.D.A. Sa and H.F. de Campos Velho, Nonextensive statistics and three-dimensional fully developed turbulence. Phys. A, 295 (2001a), 250.
• [16] F.M. Ramos, C. Rodrigues Neto, R.R. Rosa, L.D. Abreu Sa and M.J.A. Bolzan, Generalized thermostatistical description of intermittency and nonextensivity in turbulence and financial markets. Nonlinear Anal.-Theor., 23 (2001b), 3521.
• [17] A. Plastino and A.R. Plastino, Tsallis's entropy and Jaynes' information theory formalism. Braz. J. Phys., 29 (1999), 50.
• [18] B.M.R. Boghosian, Thermodynamic description of two-dimensional Euler turbulence using Tsallis' statistics. Phys. Rev., E 53 (1996), 4754.
• [19] T. Reginska, A regularization parameter in discrete ill-posed problem. SIAM J. Sci. Comput., 17 (1996), 740.
• [21] J.D. Hoffmann, Numerical Methods for Engineers and Scientists. McGraw-Hill, 1993.
• [22] E. Issamoto, F.T. Miki, J.I. da Luz, J.D. da Silva, P.B. de Oliveira and H.F. de Campos Velho, An inverse initial condition problem in heat conduction: a neural network approach, Braz. Cong. on Mechanical Eng. (COBEM), Campinas (SP), Brazil, Proc. in CD-ROM - paper code AAAGHA, page 238 in the Abstract Book (1999).
• [23] F.T. Miki, E. Issamoto, J.I. da Luz, P.B. de Oliveira, H. F. de Campos Velho and J.D. da Silva, A neural network approach in a backward heat conduction problem, Braz. Conference on Neural Networks, Proc. in CD-ROM, paper code 0008, pages 19-24, Săo José dos Campos (SP), Brazil (1999).
• [24] E.H. Shiguemori, H.F. de Campos Velho and J.D.S. da Silva, Estimation of initial condition in heat conduction by neural network, in this Proceedings - ICIPE-2002, May 26-31, Angra dos Reis (RJ), Brazil.
• [25] C.D. Rodgers, Retrieval of the atmospheric temperature and composition from remote measurements of thermal radiation. Rev. Geophys. Space Phys., 14 (1976), 609-624.
• [26] S. Twomey, Introduction to the matehmatics of inversion in remote sensing and interative measurements. Amsterdam, Elsevier Scientific, 1977.
• [27] W.L. Smith, H.M. Woolf and A.J. Schriener, Simultaneous retrieval os surface and atmospheric parameters: a physical analytically direct approach. Adv. In Rem. Sens., 7 (1985).
• [28] M.T. Chahine, Inverse Problem in Radiative Transfer: determination of atmospheric parameters. Jour. Atmos. Sci., 27 (1970), 960.
• [29] K.N. Liou, An introduction to atmospheric radiation. Academic Press, Orlando, 1982.
• [30] J.C. Carvalho, F.M. Ramos, N.J. Ferreira and H.F. de Campos Velho, Retrieval of Vertical Temperature Profiles in the Atmosphere Inverse Problems in Engineering (3ICIPE), Proceedings in CD-ROM, under paper code HT02 (see also in the internet: http://www.me.ua.edu/3icipe/fin3prog.htm) - Proc. Book: pp. 235-238, Port Ludlow, Washington, USA, June 13-18, UEF-ASME (2000).
• [31] M. Bertero and P. Boccacci, Introduction to Inverse Problems in Imaging, Institute of Physics, (1999).

Appendix

Publication Dates

• Publication in this collection
  19 Mar 2007
• Date of issue
  2006

History

• Accepted
  15 Dec 2006
• Received
  15 Dec 2006