RESEARCH ARTICLE


On Orthogonal Polynomials and Finite Moment Problem



Fawaz Hjouj1, *, Mohamed Soufiane Jouini1
1 Department of Mathematics, Khalifa University, Abu Dhabi, UAE


Article Metrics

CrossRef Citations:
8
Total Statistics:

Full-Text HTML Views: 950
Abstract HTML Views: 490
PDF Downloads: 528
ePub Downloads: 251
Total Views/Downloads: 2219
Unique Statistics:

Full-Text HTML Views: 534
Abstract HTML Views: 316
PDF Downloads: 321
ePub Downloads: 197
Total Views/Downloads: 1368



Creative Commons License
© 2022 Hjouj and Jouini

open-access license: This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International Public License (CC-BY 4.0), a copy of which is available at: https://creativecommons.org/licenses/by/4.0/legalcode. This license permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

* Address correspondence to this author at the Department of Mathematics, Khalifa University, Abu Dhabi, UAE; E-mail: fawaz.hjouj@ku.ac.ae


Abstract

Background:

This paper is an improvement of a previous work on the problem recovering a function or probability density function from a finite number of its geometric moments, [1]. The previous worked solved the problem with the help of the B-Spline theory which is a great approach as long as the resulting linear system is not very large. In this work, two solution algorithms based on the approximate representation of the target probability distribution function via an orthogonal expansion are provided. One primary application of this theory is the reconstruction of the Particle Size Distribution (PSD), occurring in chemical engineering applications. Another application of this theory is the reconstruction of the Radon transform of an image at an unknown angle using the moments of the transform at known angles which leads to the reconstruction of the image form limited data.

Objective:

The aim is to recover a probability density function from a finite number of its geometric moments.

Methods:

The tool is the orthogonal expansion approach. The Shifted-Legendre Polynomials and the Chebyshev Polynomials as bases for the orthogonal expansion are used in this study.

Results:

A high degree of accuracy has been obtained in recovering a function without facing a possible ill-conditioned linear system, which is the case with many typical approaches of solving the problem. In fact, for a normalized template function f on the interval [0, 1], and a reconstructed function ; the reconstruction accuracy is measured in two domains. One is the moment domain, in which the error (difference between the moments of f and the moments of ) is zero. The other measure is the standard difference in the norm -space ||f- || which can be ≈ 10-6 or less.

Conclusion:

This paper discusses the problem of recovering a function from a finite number of its geometric moments for the PSD application. Linear transformations were used, as needed, so that the function is supported on the unit interval [0, 1], or on [0, α] for some choice of α. This transformation forces the sequence of moments to vanish. Then, an orthogonal expansion of the Scaled Shifted Legendre Polynomials, as well as the Chebyshev Polynomials, are developed. The result shows good accuracy in recovering different types of synthetic functions. It is believed that up to fifteen moments, this approach is safe and reliable.

Keywords: Reconstruction, Function, Moments, Orthogonal expansion, Particle size, Distributions image, processing.



1. INTRODUCTION

Consider an unknown nonnegative function f with compact support on the interval (a,b.) The geometric moments of the function f are defined by:

(1)

In fact, the problem of determining f from a given finite sequence Mkf takes different forms due to the related physical and mathematical assumptions. Specifically, if Mfk is known for all orders, then f can be recovered completely, but if only a finite number of these moments are given, then we have the well-known inverse ill –conditioned problem. Furthermore, it is known that the support of the objective function f categorizes the problem into different scenarios [1-4]. The Hausdorff moment problem, where the support of f is a closed interval (a, b), the Stieltjes moment problem, where the support of f on (0, +∞), and the Hamburger moment problem, where f is supported on (−∞, +∞).

Reconstructing a function from a finite number of its moments is an essential tool for many problems in science and technology. Among many applications that involve the moment problems, three major fields are described. Namely, chemical engineering applications, tomography, and image processing. The most known existing methods of solving the moment problem are reviewed in [1]. A general review of these applications is now given.

In the field of chemical engineering, the Particle Size Distribution (PSD) is a crucial process that results in several industrial operations in chemical engineering, especially in dynamic multiphase flows. The size distribution of particles can be essential for several applications. For example, this property may have an impact on the efficiency of processes like filtration, as a small product size may result in a significant increase in processing time. Moreover, PSD generally represents a pivotal result to assess the excellent quality of chemical process output. Methods to retrieve PSD can be divided into three main categories [1-4]. The first approach is based on analytical inversion techniques. In the second approach, PSD is recovered using discrete linear equations sets implemented in the independent model algorithm. The third approach consists of reconstructing PSD from known moments. In the last decade, several research studies in chemical engineering focused on PSD applications using moment-based techniques to determine distributions [4-10]. Moreover, most of the applications implement a limited, finite number of moments for the reconstruction introducing possible inaccuracies. This moment constraint makes chemical engineering applications restricted only to elementary distribution shapes such as Gaussian or lognormal functions.

Moment-based approach for several types of image processing and tomography problems attracted the attention of many researchers, for instance [11-20]. For example, in image processing, we compute the image moments from the given projections for the purpose of shape classification or object position and orientation, recovering a general linear transformation of images, segmentation, texture analysis [12-16], and others. In the field of tomography, including digital rock physics, we deal with the well-known application of reconstructing an image from its moments or the reconstruction of the Radon Transforms from their moments [17-21]. Specifically, we use the connection between the projection moments and the image moments, and the connections among the projection moments from different views. Using these known relations [13-17], we obtain all possible image moments from the given projections. The image can then be recovered.

As said, a full review of the most known existing methods of solving the moment problem is reviewed [1].

This work considers the Hausdorff moment problem and proposes the use of orthogonal polynomial expansion techniques to solve this problem. Without loss of generality, f can be transformed and normalized so that -fxdx=1 on the support (0,1). In this way, the moments of tend to 0 as k.

This work is organized as follows: In the next section, our proposed method is presented. Section 3 is concerned with experiments and discussion. Finally, the conclusion can be found in the last section.

2. MATERIALS AND METHODS

Our tool is the orthogonal expansion. Mathematics Literature provides a rich theory and tools for the classical orthogonal expansion of functions [22-29]. The choice of several options of orthogonal bases that serve and fit the physical, geometry, and settings of our current problem is explored. Undoubtedly, the role of the Fourier series, the Legendre’s and Chebyshev polynomials in many fields of pure and applied sciences are well-established. An orthogonal expansion that allows us to represent the expansion's coefficients in terms of the moments rather than the function itself since it is not known, is required. It is found, for example, that the Fourier series is not relevant to this problem. Expressing the Fourier series coefficients in terms of moments is possible, but the numerical implementation immediately would suffer the overflow issues. On the other hand, the Legendre series and Chebyshev series proved very well stable solutions using a few moments of the target function.

2.1. Using Legendre Polynomials

The Legendre polynomials ϕm(x) form a complete orthogonal system over the interval [-1,1]:

(2)

where δmn denotes the Kronecker delta.

Since we require an approximation on the interval [ 0,1], we first use the mapping (transformation)

(3)

that moves the interval [0, α], α>0 to the interval [−1,1], and it preserves the orthogonality of the polynomials, say, Qα,n on the interval [0,α]. In particular, the orthogonality in equation (2) becomes:

(4)

We have a handy explicit formula for Qα,n, namely:

(5)

The Legendre series for f(x) is given by:

(6)

By orthogonality of Qα,n the expansion coefficient cn is written as:

(7)

Using (5), we write (7) as:

(8)

But so (8) is written as:

(9)

Using (9), equation (6) is written as:

(10)

If we consider the N-partial sum of (10), we obtain:

(11)

Thus, the approximating function Legr(x) is calculated from the given moments. Numerical discussion on (11) is presented in the next section.

2.2. Using Chebyshev Polynomials

The Shifted Chebyshev Polynomial of the second kind [25], say Vnx, is considered and arguments like what is done above are developed. These polynomials are orthogonal on the interval [0,1] concerning the inner product

(12)

with the weight function x-x2. Furthermore, the following is an analytical formula for computing these polynomials,

(13)

We can express our objective function in the following manner:

(14)

Using (12), we solve for the coefficients,

(15)
(16)

Using (13), we write (15) as

(17)

Thus, (14) is written as

(18)

Finally, we consider the partial sum of (18) to obtain:

(19)

Examples and discussion of these results are given in the next section.

3. RESULTS AND DISCUSSION

Now, the following remarks concerning our method are addressed. First, in (3.1), our theory using synthetic functions is tested and validated using MATLAB. In doing so, moments of order up to fifteen or less are employed, and a reliable approximation using either one of the two families of polynomials of this study is shown. Second, in (3.2), the behavior of a series of polynomials with a large number of terms and other properties is considered. This work is not a classical interpolation since interpolation requires a sample of the objective function. In the classical theory of interpolation, the Chebyshev series is believed to be a better approximation if the interpolation nodes are chosen to be the Chebyshev nodes. In the context of this work, the two series (11) and (19) are not comparable in terms of accuracy, given that we want to rely on a few terms of these series. However, they exhibit different behaviors in regard to errors. Overall, thanks to the restriction of using only a few moments, we will not need many terms. It is believed that up to fifteen moments, this approach is safe and reliable.

3.1. Demonstration

In testing our formulas (11) and (19), the support of the testing functions will be the interval (0, 1). Observe that in (11); f0 or f1 need not to be zero, while in (19), it is assumed that f(0) = f(1) = 0. Measures of the closeness of the reconstructed function to the original function can be:

(20)
(21)

As a first example, we consider the function:

(22)

normalized on [0,1], shown in Fig. (1a) (the green graph). The Legr(x), from (11), with moments of orders N=2,4,6,8,and 10 are computed. We display these graphs together with the graph of f in (Fig. 1a-e). The corresponding error values elk are calculated in Table 1, with a plot of elk in (Fig. 1f).

Table 1. Calculations of figure 1: ( (a-e) The function f(x)=sin(2πx)+1 on [0,1] recovered using (11) with different moments: N=2,4,6,8,10, along with the corresponding errors (20). (f) plot of elk).
Figure 1 Moment
N
elk
a 2 0.9584
b 4 0.2084
c 6 0.0662
d 8 0.0092
e 10 0.0004
Fig. (1). (a-e) The function f(x)=sin(2πx)+1 on (0,1) recovered using (11) with different moments: N=2,4,6,8,10, along with the corresponding errors (20). (f) plot of elk .

different moments: N=2,4,6,8,10, along with the corresponding errors (20). (f) plot of elk

In the second example, consider the function:

(23)

normalized on 0,1 . The Legr(x), from (11), with moments N=6,8,10,12, and 15 are computed. We display these graphs together with the graph of f (the green graph); (Fig. 2a-e). The corresponding error values elk are calculated in Table 2 with a plot of elk in Fig. (2f).

Fig. (2). a-e)_The function f(x)=exp [ -80(x-0.3)^2 ] + 1/2 exp [ -80 (x-0.6)^2 ] on (0,1) recovered using (19) with different moments: N=6, 8, 10, 12, 15 along with the corresponding errors (21). (f) plot of elk.

Table 2. Calculation of figure 2 ((a-e)_The function f(x)=exp [ -80(x-0.3)^2 ] + 1/2 exp [ -80 (x-0.6)^2 ] on [0,1] recovered using (19) with different moments: N=6, 8, 10, 12, 15 along with the corresponding errors (21). (f) plot of elk).
Figure 2 Moment
N
elk
a 6 1.1954
b 8 1.1054
c 10 0.1682
d 12 0.0396
e 15 0.0037

The assumption that f(0) = f(1) = 0 is natural in some applications. For example, in image processing and tomography, we can assume that the image is supported in the interior of the unit circle. Similarly, we consider the value of PSD out of the support (0, 1) to be small enough and can be set as zero. For the rest of our experiments and discussion, we will be under this assumption.

recovered using (19) with different moments: N=6, 8, 10, 12, 15 along with the corresponding errors (21). (f) plot of elk

Working with both families of polynomials, we test Legr(x) and Cheb(x) for the function:

fx=sinπxe10x, normalized on (0,1).

Results are shown in Table 3 and Fig. (3). Both are approximated very well using either (11) or (19) within the first twelve moments.

Fig. (3). (a-e) The function f(x)=sin(πx) e^10x normailzed on [0,1] recovered using (11), red graphs, and (19), blue graphs, with different moments: N=3,5,7,10,12 along with the corresponding errors (20) and (21). (f) plots of [el]_k and [ec]_k.

Table 3. Calculation of figure 3 ((a-e) The function f(x)=sin(πx) e^10x normailzed on (0,1) recovered using (11), red graphs, and (19), blue graphs, with different moments: N=3,5,7,10,12 along with the corresponding errors (20) and (21). (f) plots of elk and eck).
Figure 4 Moment
N
elk eck
a 3 0.0389 0.0067
b 5 0.0220 0.0060
c 7 0.0046 0.0034
d 10 0.0001 0.0085
e 12 0.0000 0.0015

Fig. (4). (a-b) The function fx=sinπx normailzed on [0,1] , with Legr(x) from (11) and Chebx from (19) using N=25
(c-d) The function fx=sinπx normailzed on [0,1] , with Legr(x) from (11) and Chebx from (19) using N=30

different moments: N=3,5,7,10,12 along with the corresponding errors (20) and (21). (f) plots of elk and eck

3.2. The Case of Large N

Although our approach never needs a large N, the completeness of this discussion needs to address several concerns regarding the approximation by polynomials. The Weierstrass approximation theorem [25] states that for every continuous function f(x) defined on an interval (a,b), there exists a set of Polynomials pix, iN that approximates f(x) with uniform convergence over (a,b) as N→∞ . However, this theorem does not provide a general method of finding such a set of polynomials. A set of polynomials may even diverge as N increases. This typically occurs in an oscillating pattern that magnifies near the end of the interpolation points. This is known as Runge's phenomenon [26]. In fact, oscillation at the edges of an interval occurs when polynomials of a high degree over a set of equispaced interpolation points are used. The phenomenon is similar to the Gibbs phenomenon in Fourier series approximations [27]. Runge's phenomenon comes from: The magnitude of the n-th order derivative of this function grows quickly for large degrees, and the equidistance between points leads to a larger Lebesgue constant for large N.

In the classical interpolation practice, this issue using different approaches is mitigated, such as change of interpolation points, using the S-Runge algorithm without resampling, use of constrained minimization, use of piecewise polynomials, Least-squares fitting, and others. Perhaps, the most common of these is the standard Chebyshev points. Although our work here is not a classical interpolation since interpolations require the knowledge of a sample of f(x) on (a,b). Both families of polynomials used here are no exception, and they, too, suffer the same behavior for large N. However, the Legendre expansion exhibits different rates of error at different values in (a,b) [28, 29], while in the Chebyshev case the error is more uniformly distributed over the interval (a, b).

For visual display, we repeat the test of the function x=sinπx normailzed on (0,1), with Legr(x) and Chebx using N=25,30, Results are shown in Fig. (4).

CONCLUSION

This paper discusses the problem of recovering a function from a finite number of its geometric moments for the PSD application. Linear transformations were used, as needed so that the function is supported on the unit interval

[0, 1], or on [0, α] for some choice of α. This transformation forces the sequence of moments to vanish. Then, an orthogonal expansion of the scaled shifted Legendre polynomials, as well as the Chebyshev polynomials, developed. The result shows good accuracy in recovering different types of synthetic functions. It is believed that up to fifteen moments, this approach is safe and reliable.

LIST OF ABBREVIATIONS

PSD = Particle Size Distribution
f = Function

CONSENT FOR PUBLICATION

Not applicable.

AVAILABILITY OF DATA AND MATERIALS

The data supporting the findings of the article is available within the article.

FUNDING

None.

CONFLICT OF INTEREST

The authors declare no conflict of interest, financial or otherwise.

ACKNOWLEDGEMENTS

Declared none.

REFERENCES

[1] F. Hjouj, and M.S. Jouini, "Review and improvement of the finite moment problem", Open Chem. Eng. J., vol. 14, no. 1, pp. 17-24, 2020.
[2] A.K. Roy, and S.K. Sharma, "A simple analysis of the extinction spectrum of a size distribution of Mie particles. J. Opt. A: Pure Appl", Opt., vol. 7, pp. 675-684, 2005.
[3] K. Bandyopadhyay, A.K. Bhattacharya, P. Biswas, and D.A. Drabold, "Maximum entropy and the problem of moments: A stable algorithm", Phys. Rev. E Stat. Nonlin. Soft Matter Phys., vol. 71, no. 5, p. 057701, 2005.
[4] V. John, I. Angelov, A.A. Öncül, and D. Thévenin, "Techniques for the reconstruction of a distribution from a finite number of its moments", Chem. Eng. Sci., vol. 62, no. 11, pp. 2890-2904, 2007.
[5] R.M. Mnatsakanov, "Hausdorff moment problem: Reconstruction of probability density functions", Stat. Probab. Lett., vol. 78, no. 13, pp. 1869-1877, 2008.
[6] H.M. Hulburt, and S. Katz, "Some problems in particle technology", Chem. Eng. Sci., vol. 19, no. 8, pp. 555-574, 1964.
[7] G.A. Athanassoulis, and P.N. Gavriliadis, "The truncated Hausdorff moment problem solved by using kernel density functions", Probab. Eng. Mech., vol. 17, no. 3, pp. 273-291, 2002.
[8] R.V. Abramov, "An improved algorithm for the multidimensional moment-constrained maximum entropy problem", J. Comput. Phys., vol. 226, no. 1, pp. 621-644, 2007.
[9] R. Abramov, "A practical computational framework for the multidimensional moment-constrained maximum entropy principle", J. Comput. Phys., vol. 211, no. 1, pp. 198-209, 2006.
[10] G. Inglese, "A note about the discretization of finite moment problems", Inverse Probl., vol. 10, no. 2, pp. 401-414, 1994.
[11] S.R. Deans, The Radon Transform and Some of Its applications. John Wiley& Sons, Inc: New York, USA, 1983.
[12] F. Hjouj, and M.S. Jouini, "On the Radon transform and linear transformations of images", ACM International Conference Proceeding Series, 2019pp. 26-31 New York, NY, United States
[13] F. Hjouj, and M.S. Jouini, "On image registration using the radon transform: review-and-improvement", DMIP '21: 2021 4th International Conference on Digital Medicine and Image Processing, 2021pp. 17-23 New York, NY, United States
[14] H. Arjah, M. Hjouj, and F. Hjouj, "Low dose brain CT, comparative study with brain post processing algorithm", 2nd International Conference on Digital Medicine and Image Processing, DMIP 2019, 2019pp. 1-7 New York, NY, United States
[15] T.J. Wang, and T.W. Sze, "The image moment method for the limited range CT image reconstruction and pattern recognition", Pattern Recognit., vol. 34, no. 11, pp. 2145-2154, 2001.
[16] B. Xiao, J.F. Ma, and J-T. Cui, "Combined blur, translation, scale and rotation invariant image recognition by Radon and pseudo-Fourier–Mellin transforms", Pattern Recognit., vol. 45, no. 1, pp. 314-321, 2012.
[17] F. Hjouj, "Towards tomography with random orientation", ACM International Conference Proceeding Series, 2019pp. 49-53 New York, NY, United States.
[18] A. Andreychenko, L. Mikeev, and V. Wolf, "Model reconstruction for moment-based stochastic chemical kinetics", ACM Trans Model Comput Simul, vol. 25, pp. 12-12, 2015.
[19] M.S. Jouini, F. Bouchaala, M.K. Riahi, M. Sassi, H. Abderrahmane, and F. Hjouj, "Multifractal Analysis of Reservoir Rock Samples Using 3D X-Ray Micro Computed Tomography Images", IEEE Access, vol. 10, pp. 67898-67909, 2022.
[20] M.S. Jouini, F. Bouchaala, E. Ibrahim, and F. Hjouj, "Permeability and porosity upscaling method using machine learning and digital rock physics", Conference Proceedings, 83rd EAGE Annual Conference & Exhibition, vol. 2022, 2022pp. 1-5 Madrid, Spain.
[21] M. Tembely, A. AlSumaiti, M. Jouini, and K. Rahimov, "The effect of heat transfer and polymer concentration on non-newtonian fluid from pore-scale simulation of rock X-ray Micro-CT", Polymers, vol. 9, no. 12, p. 509, 2017.
[22] E. Kreyszig, Advanced Engineering Mathematics, 9th ed John Willey & Sons Ltd: USA, 2006.
[23] F.W.J. Olver, D.W. Lozier, R.F. Boisvert, and C.W. Clark, NIST Handbook of Mathematical Functions. Cambridge University Press: Cambridge, United Kingdom, 2010.
[24] W. Bosch, "On the computation of derivatives of legendre functions", Phys. Chem. Earth PT. A, vol. 25, pp. 655-659, 2000.
[25] J.C. Mason, "Hands comb", In: Chebyshev Polynomials. Chapman and Hall/CRC press company: Boca Raton, 2003.
[26] J.P. Boyd, "Defeating Gibbs Phenomenon in Fourier and Chebyshev spectral methods for solving differential equations", In: A. Jerri, Ed., Gibbs Phenomenon.. Sampling Publishing: Potsdam, New York, 2007.
[27] T.H. Fay, and P.H. Kloppers, "The Gibbs’ phenomenon for fourier-bessel series", Int. J. Math. Educ. Sci. Technol., vol. 34, no. 2, pp. 199-217, 2003.
[28] I. Babuška, and H. Hakula, "Pointwise error estimate of the legendre expansion: The known and unknown features", Comput. Methods Appl. Mech. Eng., vol. 345, pp. 748-773, 2019.
[29] R.L. Burden, J.D. Faires, and A.M. Burden, Numerical Analysis, 10th ed Cengage Learning: Boston, USA, 2016.