Fast Computation of Hahn Polynomials for High Order Moments

Discrete Hahn polynomials (DHPs) and their moments are considered to be one of the efficient orthogonal moments and they are applied in various scientific areas such as image processing and feature extraction. Commonly, DHPs are used as object representation; however, they suffer from the problem of numerical instability when the moment order becomes large. In this paper, an efficient method for computation of Hahn orthogonal basis is proposed and applied to high orders. This paper developed a new mathematical model for computing the initial value of the DHP and for different values of DHP parameters ($\alpha$ and $\beta$). In addition, the proposed method is composed of two recurrence algorithms with an adaptive threshold to stabilize the generation of the DHP coefficients. It is compared with state-of-the-art algorithms in terms of computational cost and the maximum size that can be correctly generated. The experimental results show that the proposed algorithm performs better in both parameters for wide ranges of parameter values of ($\alpha$ and $\beta$) and polynomial sizes.

In order to surmount these shortcomings, researches oriented towards discrete orthogonal polynomials (DOPs).
In addition, DOPs are used to solve linear functional differential equations [21].
The robustness of DOPs is essentially based on some important properties such as energy compaction, efficient data processing, numerical stability, robust data analysis, extraction features from the signal, and localization [4,22].
However, the remarkable properties of the most DOMs can be applied only to medium size images, they are not applicable for large size images or high moment orders [23]. This limitation is determined by the DOP overflow, fluctuation of the polynomial values, as well as the high computational cost. Thus, new recurrence algorithms for generating the higher orders are still developed, e.g. for Chebyshev [17] and Krawtchouk [18] polynomials. Recently, researchers have discussed other DOPs such as Charlier polynomials [24] and Hahn polynomials [23].
The recursive algorithms reduce the complexity of calculating the coefficients of DOPs and the propagation of errors [25,26]. We can use either a single recursive formula with respect to degree n or a double recursive formula also with respect to spacial or time coordinate x. The problem of numerical instability is solved by calculating the DOP coefficients with respect to the variable n. However, this calculation is not efficient, when the size of 1D or 2D signals becomes large. For instance, the coefficients of Chebyshev polynomials have numerical instabilities since the squared norm of the scaled Chebyshev polynomials assumes small values. Mukundan [10] propoesed the recurrence algorithm in the x-direction to resolve this issue. After that, many researches began to work on this problem such as [25].
In general, an attention has been paid to the computation cost, which is considered to be an important point that subjects to ill-conditioning, therefore it is taking a substantial consideration in different researches [27,28]. For Meixner moment coefficients, this drawback is resolved via fast and efficient calculation in [29]. For Chebyshev moments, a fast and stable method is proposed by Abdulhussain et al. [17] for higher polynomial degree by a combination of the three-term recurrence relations in the nand x-directions. Daoui et al. [23] proposed a new method using a modified Gram-Schmidt orthogonalization process. This method reduces the numerical error propagation during recursive computations. However, it is relatively slow.
Motivated by this problem for the discrete orthogonal Hahn moments, this study introduces a new algorithm to tackle this issue through composing two recurrence algorithms (n-and x-recurrence relations) and an adaptive threshold to stabilize the generation of the DHP coefficients. The present paper is organized as follows: in Section 2, the preliminaries and the existed three-term recurrence algorithms are presented. In Section 3, the proposed recurrence algorithm is presented. In Section 4, the experimental analysis is performed to evaluate the proposed recurrence algorithm. Finally, conclusions are drawn in Section 5.

2/23
The mathematical definitions and fundamentals of DHP and discrete Hahn moments (DHMs) are introduced in this section.

The mathematical definition of DHP
DHPs of the nth degree H α,β n (x) are defined as the solution of the difference relation [25,30], which is given by where φ(x) and ψ(x) are first and second order functions, respectively, and λ n is a constant. ∆H α,β n (x) and ∇H α,β n (x) represent the forward and backward difference, respectively. The values of φ(x), ψ(x), and λ n are given by [25] φ where α and β are the DHP parameters (α and β > −1 or also α and β < −N ). The values of ∆H α,β n (x) and ∇H α,β n (x) are defined as follows From Eqs. (5) and (6), ∆∇H α,β n (x) can be written where n represent the polynomial degree, x is the signal index, and N is the polynomial size (the number of samples). The solution of the Eq. (1) is where 3 F 2 (·) represents the generalized hypergeometric series which is given by 3/23 and (·) k represents the rising factorial also known as Pochhammer symbol. It is given by (a) k = a(a + 1)(a + 2) · · · (a + k − 1) .

Existing Recurrence Algorithms
The three term recurrence relations are employed because of both the time consumption and the insufficient precision of the hypergeometric series in Eq. (8). In this section, the existing recurrence algorithms and their analysis are briefly presented.

The Three Term Recurrence Relation in the n-direction (TTRRnd)
The DHP of the nth degree at the xth index is given by [25] H α,β n (x; N ) = where the parameters of the n-direction recurrence algorithms are

AB EĤ
The limitation of the n-direction recurrence algorithm arises from the initial valuesĤ α,β

The Three Term Recurrence Relation in the x-direction (TTRRxd)
DHP of nth degree at the xth index is computed as [25] H α,β n (x;

6/23
where the coefficients of the x-direction recurrence algorithm are with initial valuesĤ α,β To reduce the time required to compute the DHPCs, the following symmetry relation [23] is employed Using the symmetry relation, Eq.

Recurrence Relation Based on Gram-Schmidt orthonormalization process (RRGSOP)
Recently, Daoui et al. [23] presented an algorithm based on Gram-Schmidt orthonormalization process (GSOP) and the n-direction recurrence relation to compute DHP. The GSOP is used to overcome the problem of the instability in the DHPCs. The presented algorithm begins with the computation of the initial setsĤ α,β Then, the recurrence relation in the n-direction is employed to compute the coefficients for n > 1. Finally, the GSOP is applied to minimize the numerical errors generated by the n-direction recurrence algorithm. However, the GSOP-based recurrence algorithm satisfies the orthogonality condition, it has three issues. First, the algorithm is not able to correctly generate the coefficients of the DHP when α = β, which is observed from the results in [23]. Second, the algorithm is unable to generate DHP for a wide range of parameters α and β because of the formula used to compute the initial values. Third, the GSOP-based recurrence algorithm has high computational cost due to the nested loops of the employed GSOP to minimize the error for each polynomial degree, which in turn increases the number of operations required to compute the coefficients of the DHP.

Proposed Recurrence Algorithm
In this section, the design of the proposed algorithm is presented. For simplicity DHPĤ α,β n (x; N ) is denoted aŝ H α,β n (x). Fig 3 shows DHP for α = β, which is considered to be more general than the case α = β.

Initial Value and Initial Sets
The initial values are crucial and prerequisite for computation of the DHPC values. The existing algorithms (TTRRnd, TTRRxd, and RRGSOP) compute the initial sets. While TTRRnd and RRGSOP utilize Eqs. (19) and (20), TTRRxd utilizes Eqs. (23) and (24) to determine the initial values. These equations are problematic because of the gamma and binomial functions. Thus, infinity (Inf) or not a number (NaN) are occurred in different environments such as C++, python, and MATLAB, i.e., the coefficients of the initial sets are not correctly computed. In addition, the existing algorithm computes the initial sets using the same formula, which in turn leads to increase the computation time.

Utilization of the TTRRnd and TTRRxd in the proposed algorithm
The proposed recurrence algorithm is designed based on merging the TTRRnd and TTRRxd. For the parts P1 and P2, the coefficients of the DHP are computed using TTRRxd with a maximum degree of 40% of N . For the part P1, the modified TTRRxd (mTTRRxd) is applied as followŝ where M = 0.4N to ensure non-zero initial sets. The parameters ν 1 , ν 2 , and ν 3 of the mTTRRxd are defined by For the part P2, the mTTRRxd is applied backwardly as followŝ
Step 3: Compute the coefficients of DHP in part P1 using mTTRxd given in Eq. (37), while the coefficients of DHP in part P2 are computed using backward mTTRRxd given in (39).
Step 4: Compute the coefficients of the DHP in parts P3 and P4 using backward mTTRxd given in (40)  Step 5: For each value computed at the xth index, if the absolute value of the previously computed coefficient is less than the currently computed coefficients and less than 10 −6 , the process is terminated and moved to the next value of n, i.e., n + 1. We call the zero coefficients part P6.
It should be noted that for the case of α = β, the proposed algorithm computes the coefficients in parts P1 and P3, while the coefficients in part P5 are computed using symmetry relation given in Eq. (25).

Experimental Results
The performance of the proposed algorithm is evaluated against existing algorithms in this section. The evaluation is performed in terms of the maximum size of the generated polynomial, signal reconstruction, and computational cost.

Analysis of Maximum Size
In this experiment, we searched maximum size of DHP generated by the proposed and existing algorithms for different values of parameters α and β. If R is the matrix of valuesĤ α,β n (x) with the size N × N , then R × R T should be the identity matrix 1(N) of the size N × N . So, mean(|1(N ) − R × R T |) can be used as the average orthogonality error.
We searched maximum N satisfying two criteria, average orthogonality error less than 10 −5 and the computing time less than 1 minute.
The maximum size for each algorithm is reported in Table 1. From it we can observe that the proposed algorithm is able to generate DHP with different polynomial parameters without propagation error. The main problem with existing algorithms relies on: 1) the formula used to compute the initial value, and 2) the propagation error generated when the  because of the time limit 1 minute, while the typical time of computation of the other algorithms is less than 2 s. The proposed algorithm outperforms the existing algorithms in terms of the maximum size and degree that can be generated.

Analysis of the Energy Compaction and Reconstruction Error
Discrete transforms are dissimilar because of their moments distribution [24]. The sequence of moment indices is essential for the reconstruction of signal information. Therefore, the DHP distribution of the moment energy needs to be investigated before the signal reconstruction analysis. The procedure presented by Jian [33] has been followed to find the distribution of moments. The procedure can be summarized as follows 14/23 1. A covariance matrix CM with length N and zero mean is given by [33] 2. The covariance matrix CM is transformed into the discrete Hahn moment domain as follows where T is the transformed matrix that is utilized to describe the transform coefficients σ 2 l , and R is the DHP matrix generated with a size of N and degree N .

The diagonal coefficients of the matrix T are considered
The aforementioned procedure is carried out using two values of covariance coefficients, ρ = 0.85 and ρ = 0.95, different values of DHP parameters, and length N = 16. The results are reported in Table 2. It can inferred from Table 2 that the variance values σ 2 l of the DHP are arranged such that the maximum value is located at l = 0 and the variance values reduced as the variance index increased. Thus, the DHP moment order, which is utilized to reconstruct signal information, is given by: n = 0, 1, . . . , N − 1.  One of the important properties of the orthogonal-polynomial-based discrete transform is the energy compaction property. This property is employed to measure the tendency of the DHP to reconstruct a large amount of the signal information using a small number of moments coefficients. To examine the impact of the DHP parameters α and β on 15/23 the energy compaction, the restriction error, J , is used as follows [33] J where σ 2 a represents σ 2 l ordered descendingly.  For more evaluation of the proposed algorithm, the normalized mean square error (NMSE) is computed between the image and its reconstructed version. The formula for NMSE, E, is given in [24] E(I, I r ) = x,y (I(x, y) − I r (x, y)) 2 x,y (I(x, y)) 2 (45) The image "Fruits" shown in Fig 8 is    hand, the proposed algorithm is able to compute DHPCs and to reconstruct the image correctly for different values of DHP parameters and image sizes.

Computational Cost Analysis
In this section, the proposed algorithm is evaluated in terms of computational cost. The execution time is performed for the proposed algorithm and compared to that of the existing algorithms. The execution time experiment is carried out using different DHP sizes. The experiment is performed 10 times and the average time for each algorithm is reported as shows in Fig 11. It can be observed from Fig 11 that the execution time of the proposed algorithm for α = β is less than that of the the proposed algorithm for α = β in Fig 11a. In addition, the proposed algorithm shows less computation time than that of the TTRRnd because of the proposed algorithm reduces the formula used for computation as depicted in Fig 11a. Compared to TTRRxd, the proposed algorithm shows higher execution time than that of the TTRRxd.
On the other hand, the execution time required to generate DHP using RRGSOP is higher than that of the proposed algorithm. The average improvement ratio computational cost for the proposed algorithm is ∼0.76, 2.52, and 289.59 over TTRRxd, TTRRnd, and RRGSOP, respectively.
The experiment was carried out using MATLAB environment on MSI-GT60 laptop with a memory of 16GB and core i7-4700MQ CPU.

Conclusion
In this paper, a new recurrence algorithm for DHP is introduced. The proposed algorithm uses the logarithmic gamma function for computation of the initial values so that it is able to compute initial values for different DHP parameters 19/23