On the computational aspects of Charlier polynomials

Abstract Charlier polynomials (CHPs) and their moments are commonly used in image processing due to their salient performance in the analysis of signals and their capability in signal representation. The major issue of CHPs is the numerical instability of coefficients for high-order polynomials. In this study, a new recurrence algorithm is proposed to generate CHPs for high-order polynomials. First, sufficient initial values are obtained mathematically. Second, the reduced form of the recurrence algorithm is determined. Finally, a new symmetry relation for CHPs is realized to reduce the number of recurrence times. The symmetry relation is applied to calculate 50% of the polynomial coefficients. The performance of the proposed recurrence algorithm is evaluated in terms of computational cost and reconstruction error. The evaluation involves a comparison with existing recurrence algorithms. Moreover, the maximum size that can be generated using the proposed recurrence algorithm is investigated and compared with those of existing recurrence algorithms. Comparison results; indicate that the proposed algorithm exhibits better performance because it can generate a polynomial 44 times faster than existing recurrence algorithms. In addition, the improvement of the proposed algorithm over the traditional recurrence algorithms in terms of maximum-generated size is between 19.25 and 42.85.

On the computational aspects of Charlier polynomials Alaa M. Abdul-Hadi 1 , Sadiq H. Abdulhussain 1 * and Basheera M. Mahmmod 1 Abstract: Charlier polynomials (CHPs) and their moments are commonly used in image processing due to their salient performance in the analysis of signals and their capability in signal representation. The major issue of CHPs is the numerical instability of coefficients for high-order polynomials. In this study, a new recurrence algorithm is proposed to generate CHPs for high-order polynomials. First, sufficient initial values are obtained mathematically. Second, the reduced form of the recurrence algorithm is determined. Finally, a new symmetry relation for CHPs is realized to reduce the number of recurrence times. The symmetry relation is applied to calculate ,50% of the polynomial coefficients. The performance of the proposed recurrence algorithm is evaluated in terms of computational cost and reconstruction error. The evaluation involves a comparison with existing recurrence algorithms. Moreover, the maximum size that can be generated using the proposed recurrence algorithm is investigated and compared with those of existing recurrence algorithms. Comparison results; indicate that the proposed algorithm exhibits better performance because it can generate a polynomial 44 times faster than existing recurrence algorithms. In addition, the improvement of the proposed

PUBLIC INTEREST STATEMENT
Discrete orthogonal polynomials are utilized in different fields such as image processing, speech enhancement, steganography, and others. There are many types of discrete orthogonal polynomials, Discrete Charlier polynomials is wellknown and highly utilized in different applications. The existing algorithms have limitations in size when generated discrete Charlier polynomials. In this article, we present a new algorithm for that is able to generate discrete Charlier polynomial with high order and large size. Thus, large signal and image sizes can be analyzed. algorithm over the traditional recurrence algorithms in terms of maximumgenerated size is between 19.25 and 42.85.

Introduction
Discrete orthogonal polynomials have attracted considerable attention from researchers in numerous scientific fields, particularly speech and image analyses (S. H. Pee et al., 2017), because of they are robust to noise, do not exhibit redundancy, and do not require a priori information (S. H. Abdulhussain, Ramli, Hussain et al., 2019). Consequently, these polynomials have been extensively used in different applications, such as information hiding (Radeaf et al., 2019), face recognition (Akhmedova & Liao, 2019), edge detection (S. H. Abdulhussain, Ramli, Mahmmod et al., 2017), video content analysis (S. H. Abdulhussain, Rahman Ramli et al., 2019), and speech enhancement . In addition, basis functions of orthogonal polynomials can be used as an approximate solution for differential equations (Mizel, 2008). Hu (Ming-Kuei, 1962) introduced the first type of moments, namely, geometric moments. This type of moments are used as a descriptor for images and considered invariant to image manipulation, such as translation, rotation, and scaling. However, geometric moments exhibit information redundancy, and thus, are non-orthogonal, i.e., a signal cannot be reconstructed from these moments.
To address this issue, the concept of continuous orthogonal moments was introduced by Teague (Teague, 1980). For example, Fourier-Mellin (Sheng & Shen, 1994), Legendre (Chong et al., 2004), and Zernike (Khotanzad & Hong, 1990) moments are continuous orthogonal polynomials. Such polynomials improve signal representation by using moments and reduce the effect of noise. However, continuous orthogonal moments involve coordinate transformation and approximation, result in computational load and numerical errors. Accordingly, discrete orthogonal polynomials were introduced in Mukundan (Mukundan et al., 2001) to address the issues of continuous orthogonal polynomials. Hahn (Yap et al., 2007), Krawtchouk (Yap et al., 2003), Tchebichef (Mukundan et al., 2001), Charlier (Zhu et al., 2010 are examples of discrete orthogonal polynomials. To compute discrete orthogonal moments, discrete orthogonal polynomial coefficients (DOPCs) should be computed accurately.
The computation of DOPCs is computationally expensive and may lead to numerical errors due to the presence of a hypergeometric series and gamma functions. Thus, three-term recurrence (TTR) algorithms are utilized to compute DOPCs. The two types of recurrence relations are generally the TTRs in the n-and x-directions. However, different studies  have shown that DOPC values result in numerical propagation error when used to compute DOPCs for high orders. Thus, new recurrence algorithms are presented to address the problem of numerical propagation errors for Tchebichef and Krawtchouk polynomials.
Studies performed on Charlier polynomials (CHPs) are divided into 1) recurrence relation algorithms and 2) moment computation algorithms. Research involving recurrence algorithms has been introduced into n-and x-direction recurrence algorithms. Both types of recurrence algorithms cannot generate high-order polynomials because of the initial values used and the number of recurrence times performed. The second stream of studies involves the computation of Charlier moments (CHMs) with reduced computational cost. However, these studies are based on either the n-or x-direction recurrence algorithm. To the best of our knowledge, there are no previous study has discussed the ability of the CHP. In addition, the main concerns of researchers are increasing the speed and minimizing the error (Razian & MahvashMohammadi, 2017). Motivated by the aforementioned issues and discussion, the current study aims to present a new recurrence relation to compute CHPs efficiently for high-order polynomials.
In this paper, an empirical study is performed to investigate the maximum order that can be generated using the existing algorithms and analyzing their problems. Due to the importance of the initial value and its effect on the computation polynomial coefficients, a sufficient initial values is proposed. In addition, a reduced form of the recurrence algorithm is presented to reduce the computation time as well as the propagation error. Lastly, we present a new symmetry relation that will reduce the numerical propagation error via the reduction of the recurrence times used.
The rest of this paper is structured as follows. The basic aspects of computing CHPs and the existing TTR algorithms are described in Section 2. The proposed recurrence algorithm is presented in Section 3. An experimental analysis is performed to evaluate the proposed recurrence algorithm in Section 4. Lastly, conclusions are drawn in Section 5.

Computation of Charlier polynomials and Charlier moments
CHPs are defined and CHMs are computed in this section. The existing TTR relation is also introduced.

The definition of CHP
The definition of CHP C p n x ð Þ for the nth order is obtained from (Zhu et al., 2010) as follows: where p represents the parameter of the CHP, and 2 F 0 is the hypergeometric series which is given by (Koekoek et al., 2010): where ðaÞ b represents the ascending factorial, which is also known as the Pochhammer symbol, and defined as follows: On the basis of Equations (1) and (2), CHP can be expressed as: CHP fulfills the orthogonality conditions such that : where ω C ðx; pÞ and ρ C ðn; pÞ are the weight function and squared norm of CHP C p n x ð Þ, respectively. The weight function is defined as follows (Jahid et al., 2019): and the squared norm ρ C ðn; pÞ is expressed as follows (Jahid et al., 2019): The computation of CHP coefficients (CHPCs) by using Equation (4) produces numerical instability. Therefore, a weighted and normalized CHP (WNCHP) is required (S. H. . The nth order of WNCHP is expressed as follows:

Definition of Charlier moments
Moments, also known as transform coefficients, are defined as scalar quantities used to represent signals without redundancy (Abdulhussain et al., 2019a;. To represent a 1D signal f ðxÞ in the moment domain of CHP, CHMs are computed as follows: where ϕ n represents the CHMs of the function f ðxÞ, Ord is the maximum order used to represent the signal. The signal f ðxÞ can be reconstructed from CHMs as follows: For a 2D signal f ðx; yÞ with a size of N 1 Â N 2 , 2D CHMs ϕ nm is computed as follows: n ¼ 0; 1; . . . ; Ord 1 and m ¼ 0; 1; . . . ; Ord 2 where Ord 1 and Ord 2 represent the maximum order utilized for signal characterization. To reconstruct the 2D signal f ðx; yÞ from CHMs, the process is performed as follows: Notably, the reconstructed signal is f ¼ f when all the moments are used in the reconstruction process.

Existing three-term recurrence algorithms
The computation of CHPCs by using Equation (8) is considered computationally expensive because of the utilization of the hypergeometric and gamma functions. Consequently, the TTR algorithm is used to compute CHPCs . Existing TTR algorithms are performed in the n-and x-directions. These TTR algorithms are presented and analyzed in detail in the following subsections.

Recurrence Algorithm in the n-direction
The TTR algorithm of WNCHP in the n-direction is presented as follows (Zhu et al., 2010): where n ¼ 1; 2; . . . ; N À 1; and x ¼ 0; 1; . . . ; N À 1, such that and initial values of TTR arê where x ¼ 0; 1; . . . ; N À 1. The TTR algorithm in the n-direction computes the WNCHP coefficients in the nth order by using the WNCHP values of the ðn À 1Þth and ðn À 2Þth orders. For the sake of simplicity the TTR algorithm in the n-direction is denoted by TTRn.
In practice, this TTRn algorithm is inapplicable to large signal sizes and high-order polynomials. Figure 1 shows the WNCHP for different size (N) and parameter (p) values. The coefficient values and distribution of WNCHP is affected by the value of the CHP parameter (p). When N ¼ 150, WNCHP is generated without any numerical instability; however, stable WNCHP coefficients cannot be generated when N ¼ 180 because of the numerical instabilities and accumulated errors resulting from the nature of the recursive computation of polynomial coefficients . In addition, numerical errors occur in TTR because the initial values are zero (as shown in Figure 1) due to the nature of the formula used. To solve this problem, a TTR algorithm in the xdirection is proposed.
In practice, the TTRx algorithm cannot handle signals with high orders and large sizes. Figure 2 shows the WNCHPs generated using the TTRx algorithm for different size (N) and parameter (p) values. When N ¼ 150, WNCHP is completely generated without numerical errors. However, the TTRx algorithm cannot generate firm WNCHP coefficients when N ¼ 180 because of numerical instability (see Figure 2(d)), and the zero values of the computed initials (see Figure 2(e) and 2(f)).

Proposed recurrence relation
The proposed algorithm for computing CHPCs is presented in this section. In the following subsection, initial value analysis, CHP identity, and the proposed TTR algorithm based on the initial value and CHP identity are discussed. Then, a summary of the proposed algorithm is provided for further clarification. For more clarification, the nomenclatures and the abbreviation used in this paper are listed in Tables 1 and 2, respectively.

Initial values
The computation of the initial value is essential and affects the other values of CHPCs.
Both TTR algorithms (i.e., TTRn and TTRx) rely on two sets of computed initial values. For example,Ĉ p 0 x ð Þ andĈ p 1 x ð Þ are the two initial sets used in TTRn. In practice, computing the initial values by using Equations (15) and (16) for TTRn or Equations (19) and (20) for TTRx is infeasible because of zero and incorrectly computed values (Figures 1(d) and 2(d)). Figure 1(d) shows that the values within the range n ¼ 0 and x ¼ 172; . . . ; 174 are generated with zerovalue coefficients. In addition, the values within the range n ¼ 0 and x ¼ 175; 176; . . . ; 179 are incorrectly generated due to the nature of the formula used. The incorrectly generated coefficient values are indicated as "Not a Number" in several testing environments, such as MATLAB, C++, and Python. In Figure 2(d), the values within the range x ¼ 0 and n ¼ Figure 1. WNCHP using TTRn for the following signal length and parameter values: From Equation (21), the factorial swiftly increases as the values of index x increase. In addition, the term p x increases as x increases. In practice, these terms are computed incorrectly. Note that, the same previous problem of TTRn is taken place for TTRx.
To overcome this problem, one initial value can be computed and a two-term recurrence can be used. This process can be performed by first computingĈ p 0 0 ð Þand then estimating the values of C p n 0 ð Þ recursively within the range n ¼ 1; 2; . . . ; N À 1 as follows: In practice, however, the initial valueĈ p 0 0 ð Þ cannot be computed for the parameter value p ¼ 1480. Thus, the best initial value that can be used to start computing CHPCs is a coefficient with a recognizable value. To achieve a recognizable initial value and considering that all the CHPCs rely onĈ p n 0 ð Þ, the coefficient ofĈ p n 0 ð Þ is inspected. Figure 3 shows the values ofĈ p n 0 ð Þ for different values of parameter p when N ¼ 300.
In 3, the peak value ofĈ p n 0 ð Þ for different parameter p values is located at n ¼ p À 1 and n ¼ p.
For example, the peak ofĈ Nevertheless, this initial value cannot be estimated for large values of parameter p due to exponentiation and factorial functions. The key solution to this problem is using the logarithmic function as follows: where logΓ is the logarithmic gamma function (Hart, 1978). Figure 4 presents the result ofĈ   To compute the remaining CHPCs forĈ p n 0 ð Þ, a two-term recurrence relation can be obtained and used. For the range n ¼ p À 2; p À 3; . . . ; 0, the two-term recurrence relation is expressed as follows: For the range n ¼ p þ 1; p þ 2; . . . ; N À 1 is presented as follows: The two-term recurrence relation can be used to compute the CHPCs ofĈ n ¼ 0; 1; . . . ; N À 1

CHP identity and recurrence relation
In this section, we present the WNCHP identity that will reduce the numerical error of CHPCs based on TTR. This identity is used to compute CHPCs for one portion and then find CHPCs for the other portion. This identity is symmetric to the primary diagonal (n ¼ x).

The TTR algorithm
From the presented WNCHP identity and the computed initial values, CHPCs in P1 are computed using the modified TTRx algorithm, as follows:

Summary of the proposed algorithm
This section summarizes the steps of the proposed TTR algorithm to further understand the computation of WNCHP. The steps are summarized as follows.
Step 2-The value ofĈ p n 0 ð Þ are computed as follows.
The aforementioned steps are illustrated in Figure 6.

Experimental analysis
This section presents the performance evaluation of the proposed method in terms of computational cost, energy compaction, signal reconstruction ability, and maximum size of the generated polynomial.

Analysis of computational cost
+ The performance of the proposed algorithm in terms of computational cost in obtaining CHPCs should be demonstrated. The computational cost of the proposed algorithm is compared with those of existing algorithms, namely, TTRn and TTRx. In addition, a recent algorithm is included in Figure 6. Illustration of the steps of the proposed algorithm.
The computation time experiment is performed on TTRn, TTRx, and the proposed algorithm using different CHP sizes and parameter p ¼ N=2. Figure 7 shows the average computation time results for 10 runs. The time required to generate CHPCs using the TTRn and TTRx algorithms is comparable. For example, the time required to generate CHPCs using the TTRn and TTRx algorithms at N ¼ 2000 is ,0.619 s and ,0.607 s, respectively. In addition, the time required to generate CHPCs remarkably increases polynomial size increases when using the TTRn and TTRx algorithms because of 1) the number of computed coefficients, i.e., N Â N, and 2) the mathematical expression used in the computation of CHPCs. However, the computation time exhibits higher reduction than those of existing algorithms. For example, at N ¼ 2000, the computation time is ,0.014 s, i.e., the execution time improvement percentage (ETIP) is ,97.7%. Notably, ETIP is computed using the following formula (Hosny, 2008): where Time P and Time E are the computation times of the proposed and existing TTR algorithms, respectively. To further illustrate the advantage of the proposed algorithm, Table 3 lists the ETIP values of the proposed polynomial compared with those of existing recurrence algorithms (TTRn and TTRx). The proposed algorithm reduces the time required to generate CHP by ,97.70%. To further explain the obtained results, it is well-known that the expensive operation that influence the computation cost of the recurrence algorithm are multiplication, division, and recurrence times. Due to the employment of the identity relation in the proposed algorithm, only 50% of the coefficients are computed. Accordingly, 50% of the expensive operations have been excluded. For instance, at N ¼ 1000, there are 10,000 coefficients should be computed using the existing TTR algorithms. Each coefficient is computed by performing the TTR relation and each recurrence relation has parameters need to be computed. These parameters involve multiplication and division operation, which increase the execution time. On the other hand, for the proposed algorithm, only 5000 coefficients are computed, i.e., 50% of recurrence relations are excluded. Thus, the proposed algorithm reduces the number of operations required and it is considered effective in terms of computational complexity for a wide range of CHP sizes.
To emphasis the ability of the proposed algorithm in terms of the computational cost, the execution time speedup ratio (ETSR) is calculated between the proposed algorithm and TTRn, TTRx, and GSOP. Table 4 lists the ETSR between the proposed algorithm and existing algorithms. The ETSR is computed as follows: It can be inferred from Table 4 that the proposed algorithm outperforms the existing recurrence algorithms in terms of speedup ratio. For instance, when N is 500, the average ETSR is 44x over TTRn and TTRx. The obtained speedup ratio is due to the reduction of the recurrence times. In addition, the average speedup ratio of the proposed algorithm over the GSOP is 17906x. The GSOP algorithm is considered computationally overload because of the utilization of the nested loops in its implementation.

Analysis of polynomial order and energy compaction
The distributions of signal moments typically differ in each transform (Abdulhussain et al., 2019a). In addition, the arrangement of moment sequence indices n is important in retrieving significant signal information. Thus, the moment energy distribution for discrete Charlier transform (DCHT) should be known prior to the analysis of signal reconstruction. Such analysis can be performed by following the mathematical procedure presented in (Jain, 1989;Zhu et al., 2010), which is based on a stationary Markov sequence with zero mean and length N. Suppose that the matrix M is the covariance matrix with different covariance coefficients (r). This matrix is given as follows (Jain, 1989): After transforming the covariance matrix (M) into the Charlier transform domain, the diagonal of the transformed matrix S is used to represent the variance σ 2 l of the transform coefficients. To obtain matrix S, the following function is applied: where R refers to the WNCHP matrix, and ðÁÞ T represents the transpose operation of a matrix. Different values of parameter p and the coefficient covariance r are used in this experiment. The results of the experiment are listed in Table 5 by using p ¼ N=4; N=2; . . . ; and 3N=4 and r ¼ 0:85; and 0:95.
The maximum value of the variance (σ 2 l ) of DCHT is located at index l ¼ 0, and the variance decreases as the index value l increases. Consequently, the DCHT moment order, which can be used to reconstruct signal information, is given as follows: n ¼ 0; 1; 2; . . . ; N À 1.
Energy compaction is considered an important property in checking the performance of discrete transform. It is used to evaluate the tendency of CHP to transform a large fraction of signal energy into relatively few coefficients of moments. To check the effect of parameter value p on the energy compaction performance of CHP, the normalized restriction error (NRE) (Jain, 1989) is utilized as follows (Jain, 1989): where J m is the NRE, and σ 2 q is the arrangement of σ 2 l in descending order. NRE is computed using different values of parameter p and covariance coefficients r. Figure 8 presents NRE (J m ) as a function of the number of retained samples q. Two values of covariance coefficients (r) are used. EC is influenced by the CHP parameter p. For example, when parameter p increases from N=6 to N=3, the performance of CHP in terms of energy compaction increases because the NRE of p ¼ N=3 reaches the minimum value before p ¼ N=6. When the value of p increases to more than N=3, NRE is smaller than that of p ¼ N=3. Therefore, CHP at p ¼ N=3 uses less moments than the other values of p to reconstruct a signal efficiently.

Analysis of reconstruction error and maximum size
In this section, experiments are performed to evaluate the proposed algorithm against existing algorithms in terms of reconstruction error and maximum size can be generated. Different image sizes and CHP parameter values are used in the experiment. A cameraman image is selected as the test image for this experiment. The test image is resized several times to obtain I. For each image size,, WNCHP R is generated using TTRn, TTRx, and the proposed algorithm. Then, the image is transformed into the moment domain M. Thereafter, the image is transformed into the spatial domain I r . Lastly, the normalized mean square error (NMSE) between the input and reconstructed images is computed as follows (S. H. : where E n represents the NMSE between images I and I r . The image is resized three times (N ¼ 144; 160; and 256), and the three values of the CHP parameter are (p ¼ bN=3c; bN=2c; and ; b2N=3c). The experiment results of NMSE for TTRn, TTRx, and the proposed algorithm are presented in Figure 9.
As shown in Figure 9, the image can be fully reconstructed using TTRn, TTRx, and the proposed algorithm for image size N Â N ¼ 144 Â 144 and all the values of parameter p. However, for image size N Â N ¼ 160 Â 160, the image can be fully reconstructed using TTRn for parameter values bN=3c and bN=2c but fails to recover any image information for any order when the parameter value is b2N=3c. In addition, the reconstruction of the image by using TTRx can be fully performed for parameter values bN=3c and b2N=3c, and the reconstructed image exhibits an error when the parameter value is bN=2c for a moment value greater than 40 because of the zero values of the polynomial coefficients for high-order polynomials (see Figure 2). For image size N Â N ¼ 256 Â 256, image reconstruction by using TTRx fails to reconstruct the image for any parameter value due to the propagation error in the computation of CHP. Meanwhile, TTRn successfully reconstructs the image only for CHP parameter values bN=3c and b2N=3c but fails to reconstruct any image information for the parameter value of bN=2c. By contrast, the proposed algorithm can efficiently reconstruct the image for all the image sizes and CHP parameter values. The maximum size that can be generated is further investigated using the proposed polynomial. To find the maximum size that can generate CHP without propagation error for each algorithm, the following steps are performed for each algorithm: (1) Small values are set accordingly for polynomial order N and parameter p.
(2) The cameraman image I is resized to N Â N.
(3) CHP R with an order of N is generated.
(5) Image I r ¼ R T Â M Â R is reconstructed.
(6) The NMSE between I and I r is found.
(7) If the NMSE less than 10 À3 , then the size is increased; otherwise, the previous value of N is considered the maximum size that the algorithm can generate. Note for small image size the threshold is set to 18 Â 10 À3 Table 6 lists the results obtained from the aforementioned procedure. The results include the maximum size for TTRn, TTRx, and the proposed algorithms corresponding to different values of the CHP parameter p. The results show that the proposed algorithm is better than existing algorithms in terms of the size that can be generated for CHP with different values of parameter p. Improvement in the proposed algorithm and the best performance in existing algorithms (TTRn and TTRx) are also reported to for better illustration. The improvement values show that the minimum improvement is 19.25 at p ¼ 5N=6, which increases to 42.85 at p ¼ N=2.

Conclusion
This study introduces a new identity of WNCHP and the reduced form to compute the initial values of WNCHP. Then, they are used to present a new recurrence algorithm for the efficient computation of WNCHP with high order. The analyses of the recurrence relation and numerical instability of CHP are explained. The identity relation of CHP to the main diagonal is beneficial for reducing the number of recurrence times. The CHP plane is divided into two portions, CHPCs are computed for one portion, and then the identity relation is utilized to find the CHPCs of the second portion. The analysis clearly indicates that the proposed recurrence algorithm has significantly less computation load than existing recurrence algorithms. In addition to its ability to reduce numerical propagation errors, the proposed algorithm also accelerates the computational speed of CHPCs. Although the proposed recurrence algorithm achieves remarkable results, its performance can still be improved to generate higher-order polynomials by examining the causes of the zero values in the coefficients, particularly for parameter p values that are less than N=2.