Reliable Recurrence Algorithm for High-Order Krawtchouk Polynomials

Krawtchouk polynomials (KPs) and their moments are promising techniques for applications of information theory, coding theory, and signal processing. This is due to the special capabilities of KPs in feature extraction and classification processes. The main challenge in existing KPs recurrence algorithms is that of numerical errors, which occur during the computation of the coefficients in large polynomial sizes, particularly when the KP parameter (p) values deviate away from 0.5 to 0 and 1. To this end, this paper proposes a new recurrence relation in order to compute the coefficients of KPs in high orders. In particular, this paper discusses the development of a new algorithm and presents a new mathematical model for computing the initial value of the KP parameter. In addition, a new diagonal recurrence relation is introduced and used in the proposed algorithm. The diagonal recurrence algorithm was derived from the existing n direction and x direction recurrence algorithms. The diagonal and existing recurrence algorithms were subsequently exploited to compute the KP coefficients. First, the KP coefficients were computed for one partition after dividing the KP plane into four. To compute the KP coefficients in the other partitions, the symmetry relations were exploited. The performance evaluation of the proposed recurrence algorithm was determined through different comparisons which were carried out in state-of-the-art works in terms of reconstruction error, polynomial size, and computation cost. The obtained results indicate that the proposed algorithm is reliable and computes lesser coefficients when compared to the existing algorithms across wide ranges of parameter values of p and polynomial sizes N. The results also show that the improvement ratio of the computed coefficients ranges from 18.64% to 81.55% in comparison to the existing algorithms. Besides this, the proposed algorithm can generate polynomials of an order ∼8.5 times larger than those generated using state-of-the-art algorithms.


Introduction
Digital image processing plays an essential role in several aspects of our daily lives. Image signals are subject to several processes such as transmission [1], enhancement [2], transformation [3], hiding [4], and compression [5,6]. Similarly to image processing, speech signals processing is also essential [7], and involves several stages such as transfer [8], acquisition [9], and coding [10]. Pattern recognition, which is considered an automated process, This is attributed to the fact that pitfalls may happen even when small errors in floating numbers occur. As such, there is an essential need to reduce the number of recurrences, especially when the polynomial size is increased. Furthermore, such a reduction could also lead to a reduction in the propagation error, thereby leading to a more stable computation of polynomial coefficients, as desired. The work in [31] proposes a modified recurrence algorithm in the n direction (RAN) by partitioning the KP array into two partitions. Therefore, only 50% of the coefficients need to be computed. However, the partitioning of the KP array generates a larger polynomial size, which is undesirable. On a similar basis, the work in [42] proposes a recurrence algorithm in the x direction (RAX) by partitioning the KP plane into two partitions. Specifically, the x direction of the recurrence algorithm is used to compute the KPCs. The results show that the RAX algorithm outperforms the RAN algorithm. It is worth noting that the RAN and RAX algorithms use a symmetric property to compute the polynomial coefficients of the second portion. A novel bi-recursive relation algorithm in the n direction (BRRN) was proposed in [43]. In this method, the KP array is divided into four partitions. However, the KPC coefficients are computed for two partitions only, i.e., 50% of the coefficients are computed. Then, a symmetric property is used to compute the KPCs for the remaining partitions. The results indicate that the BRRN algorithm provides higher gain than the RAX algorithm for limited values of parameter p, i.e., the polynomial size. Abdulhussain et al. [44] developed an algorithm and presented new properties of orthogonal polynomials such that the KP plane is divided into four portions and only the KPCs for one portion are computed. For this, the size of the generated polynomials is increased, but it is still limited, especially for parameter p less than 0.25 and greater than 0.75. This is because the initial values or sets become zero as the polynomial size increases. Recently, the work in [25] proposed a recurrence relation algorithm that has the ability to compute KPCs with very large sizes. However, the proposed algorithm is limited to the parameter value of p = 0.5.
The existing algorithms suffer from the following limitations: (1) no initial value is provided; (2) the propagation error is high; and (3) the implementation of these algorithms is limited to a specific value of the parameter p. They also suffer from numerical instabilities, especially when the polynomials orders and sizes become high. Therefore, an advanced and reliable recurrence algorithm for high-order polynomials and large sizes is required. Therefore, a new recurrence algorithm is presented in this paper, which handles the numerical instabilities issue of using high orders of polynomials and large sizes. The proposed algorithm is able to compute the KPCs for all values of the parameter p. In addition to this, this paper presents the development of a new mathematical model for computing the initial value of p. In particular, the initial value is accurately computed for all values of parameter p. Furthermore, a new relation to compute the values of the initial sets is derived. To this end, a diagonal recurrence relation is introduced. The proposed diagonal recurrence algorithm is derived from the existing n and x directions of the recurrence algorithm. The diagonal and the existing recurrence algorithm are exploited to compute the KP coefficients. The KP coefficients are then used for one partition after dividing the KP plane into four partitions. To compute the KP coefficients in other partitions, a symmetric property relation is utilized.
Organization of the paper: This paper is organized as follows. Section 2 presents the mathematical formulations of the orthogonal polynomials and moments. In Section 3, the methodology of the proposed recurrence algorithm is provided. This methodology involves providing a discussion about the initial value selection of parameter p. In addition, this section explains how the Krawtchouk polynomial's coefficients can be computed. In order to characterize the performance of the proposed approaches, Section 4 provides the numerical results. Finally, conclusions are discussed in Section 5.
Notation: In this paper, the operator transpose is denoted by (·) T and ( a b ) denotes the binomial coefficients.

Preliminaries
This section presents the Krawtchouk polynomials and their recurrence relation. To this end, the n-th order of the Krawtchouk polynomials based on the hypergeometric series is given asK The weighted function ω(x, p) and the norm function ρ(n, p) are used to generate the weighted and normalized KP coefficients as given in [31]. To this end, the weighted and normalized KP coefficients can be written as in (2) and (3), respectively: The Pochhammer symbol (·) c , which is known as an ascending or rising factorial function, can be written as [45] (a) c = Γ(a + c) where Γ(.) denotes the Gamma function. To this end, using the weight and norm functions, the weighted and normalized Krawtchouk polynomials of the n-th order for a signal of size N are given as [31] n, x = 0, 1, . . . , N − 1; p ∈ (0, 1), where 2 F 1 describes the hypergeometric series and can be written as

Proposed Recurrence Algorithm
This section describes the methodology of the proposed recurrence algorithm.

Computing the Initial Value
The problem with traditional approaches for computing the initial value in the n = 0 and x = 0 directions is the numerical instability. For example, the traditional methods provide zero values of the initial K p 0 (0)-which is unstable. The initial value can be computed as The expression in (8) makes the initial value (K p 0 (0)) decrease to zero for different values of parameter p, especially for large polynomial sizes N. Figure 1 shows the values of K p 0 (0) for different values of parameter p and size N. The results show that the value of K p 0 (0) starts to fall to zero when p becomes larger than 0.1. Specifically, as the values of parameter p increase, the value of K p 0 (0) falls to zero. For example, for p = 0.15, the value of K p 0 (0) becomes zero for N > 5000 while for p = 0.4, the value of K p 0 (0) becomes zero earlier for N > 2000. This makes it impossible to compute the rest of the KP coefficients' values. Therefore, there is an essential need to find an efficient method for computing the initial value of p, which prevents the initial value (K p 0 (0)) from dropping to zero. To this end, this paper identifies the suitable non-zero values in the KP plane that need to be used as an initial value. As such, there is an essential need to plot the values of coefficient K p 0 (x) for different values of parameter p. Figure 2 shows the plots of the values of coefficient K p 0 (x) for different values of parameter p. Clearly, the results in Figure 2 show that the values of K p 0 (x) start with a small number and gradually reach the peak. Then, the values drop to a very small number. In addition, we observe that using non-small values as an initial set of parameter p seems useful to compute other values of the KP coefficients.   Figure 2 demonstrates that the peak values can be located at x = N p. In this paper, the value of x = N p is denoted by x 0 . To this end, a general formula for computing K p 0 (x 0 ) can be written as Computing K p 0 (x 0 ) using expression (9) may also lead to unstable numerical values of coefficients with errors, especially for high-order polynomials. This is because the binomial coefficients function tends to be very large and close to infinity. To demonstrate this behavior, Figure 3 is provided.  Figure 3 shows that the initial values (K p 0 (x 0 )) are still inaccurate where these values record either NaN or Inf values. This is due to the nature of the polynomial coefficients that are obtained by the expression in (9). Thus, the initial values (K p 0 (x 0 )) seem difficult to be computed for large polynomial sizes (N). To overcome this issue, this paper proposes an efficient and suitable approach that makes the value of K p 0 (x 0 ) commutable. It is worth noting that the values obtained from the polynomial coefficients' formula, especially the Gamma function, should be reduced when the coefficients become large since their argument value is increased. This can be achieved using arg = exp(ln(arg)) = e ln (arg) and the initial values can be computed as where z is given as A proof of expression (10) is presented in Appendix A. Figure 4 shows a plot of the proposed initial values of (K p 0 (x 0 )) in the developed expression in (10) for various values of parameter p as a function of polynomials size N. The results show that the proposed initial values are more computable for wide ranges of parameter p and large polynomial sizes N, as desired. Hence, this signifies the feasibility of the proposed formula for practical implementations compared with state-ofthe-art equations.

The Fundamental Computation of the Initial Values
Typically, for any orthogonal polynomial, the computation of coefficients requires the evaluation of a significant number of fundamental initials. Thus, based on the first initial value K p 0 (x 0 ), computed using the proposed formula, the KP coefficients are obtained Figure 5). Therefore, this section shows how the aforementioned coefficients values are computed.
First, K p 0 (x 1 ) is computed using the proposed derived formula, which provides the two terms relation between the K p 0 (x 0 ) and K p 0 (x 1 ). To this end, this relation/ratio between the coefficients can be formulated as where x 1 = x 0 + 1 = N p + 1. Thus, the expression in (11) can be further simplified to: Then, K p 1 (x 0 ) and K p 1 (x 1 ) can be computed using a two-term recurrence relation with K p 0 (x 0 ) and K p 0 (x 0 ), respectively. To derive the recurrence relation of the proposed approach, the following formulas are used: From (13) and (14), K p 1 (x) can be simplified to: Using the expression in (15)

The Computation of the Initial Sets
In this section, the computation of the initial sets is discussed. These initial sets are shown in Figure 6. Figure 6 shows parts of the KP coefficients, which are covered in this section. The initial sets are defined in the ranges of x = x 0 , x 1 and n = 2, 3, . . . , x. To compute the values of the initial set, the recurrence in the n direction is used. To this end, the formulation of recurrence is given as After computing the initial sets, the values in the ranges x = x 0 , x 1 and n = 0, 1, . . . , x are used as the initials to compute the rest of the KP coefficients values.

Computation of the Coefficients Values for KP
In this section, the rest of the coefficient values are computed. These coefficients are shown in Figure 7. As depicted in Figure 7, there are two main parts, which are located at the left (Part 1) and the right (Part 2) sides of the initial sets. In addition, the coefficients are located at the right side of the initial sets and can be divided into three sub-parts. These parts are Part 2-1, Part 2-2, and Part 2-3. The detailed description of the computation of each part is presented in the following subsections.
The values of KP coefficients become unstable as the index of x goes towards n. This is because the values of the coefficients tend to be less than 10 −7 . To overcome this issue, the condition of a threshold value is used to stop the recurrence for each value of index n. The proposed condition is given by

Computation of the Coefficients Located at Part 2-1
In this section, the values of KP coefficients in Part 2-1, given in Figure 7, were computed using a forward x recurrence relation as given in (21): The aforementioned recurrence relation, which is used to compute the values in Part 2-1, is subject to the following condition:

Computation of the Coefficients Located at Part 2-2
This section presents two new recurrence relations' approaches to compute the KP coefficient values diagonally. This diagonal calculation is given in Part 2-2 Figure 7. The values in the diagonal of Figure 7 are then used to compute the coefficients' values in Part 2-3 in Figure 7. This is because the recurrence relation in the n direction cannot be used to compute the coefficients' values. Consequently, some values in Part 2 become zero, which results from the condition used to prevent the occurrence of unstable values.
This paper derives the recurrence relations provided in Figure 8. From Figure 8a, it can be seen that the elements computed for x 0 and x 1 can be used to compute the coefficients along the main diagonal n = x and n = x − 1. Furthermore, to compute the coefficients' values of KP K p n (x + 1), the coefficient value K p n+1 (x) is computed using the n direction recurrence algorithm. The similarity across the main diagonal (n = x) is exploited for simplicity where K p n (x + 1) = K p n+1 (x). To this end, the KP coefficients along n = x − 1 are computed as To compute the values at the main diagonal where n = x, a new recurrence relation approach is developed. This is achieved by combining both n and x directions recurrences. Suppose that the values at (n, x + 1), and (n − 1, x + 1) are known (see circulated values I and K in Figure 8d). Then, the value at n + 1, x + 1 (see circulated values L in Figure 8d) can be computed using the n direction recurrence relation as The value at (n − 1, x + 1) can be computed using the x direction recurrence relation as Substituting Equation (25) in (24) yields the following general expression of the recurrence relation: This recurrence relation is termed as the four-term recurrence relation in the n − −x direction. This new development approach is used to compute the KP coefficients in the range x = x 1 + 1, x 1 + 2, . . . , N/2 + 1 and n = x − 1, and x = x 1 + 1, x 1 + 2, . . . , N/2 and n = x as shown in Figure 9:

Computation of the Rest of the KP Coefficients
This subsection provides the computation of the rest of the KP coefficients. To this end, the rest of the coefficients can be computed using a similarity relation of the KP. The coefficients in the ranges x = 0, 1, . . . , N/2 − 1 and n = x + 1, x + 2, . . . , N − x − 1 are given as K p n (x) = K p x (n).
The coefficients in the ranges x = 0, 1, . . . , N − 1 and n = N − x, N − x + 1, . . . , N − 1 are computed using the following expression: In addition, to calculate the KP coefficients for p > 0.5, firstly the value of p is set to 1 − p. Then, the KP coefficients are computed using the proposed methodology. Finally, the following formula is applied for all coefficients [44]:

Summary of the Proposed Algorithm
In this subsection, a summary of the proposed algorithm is presented. To this end, a flow chart of the proposed recurrence is shown in Figure 10. In addition to this, a pseudocode is presented (see Algorithm 1) for more clarification. In addition, 3D plots of the KP coefficients are given in this subsection.  Flag=True; p ← p − 1 4: end if 5: x 0 ← N p, x 1 ← x 0 + 1 6: Compute K p 0 (x 0 ) using (10) 7: Compute K p 0 (x 1 ) using (12) 8: Compute K p 1 (x 0 ) and K p 1 (x 1 ) using (16) and (17)  9: Compute initial set 10: for x = x 0 : x 1 do 11: for n = 2 : x do 12: Compute K p n (x) using (18) Figures 11 and 12 show a 3D plot of the KP coefficients, which are generated using the proposed recurrence algorithm with N = 2000 and different values of the p parameter ranging between <0.5 and >0.5, respectively.

Numerical Results and Analyses
This section presents the results obtained using the proposed recurrence algorithm. In addition, a comprehensive comparison is conducted with the existing recurrence algorithms. The comparison is carried out in terms of the energy compaction, reconstruction error, and computation cost.

Energy Compaction Analysis
The order of moments n impacts the process of signal reconstruction, energy compaction, and information retrieval. The order of the KP moments is given by n = 0, 1, 2, . . . , N − 1. The energy compaction is utilized to check the impact of using KP to transform a large fraction of the signal energy into relatively few coefficients of moments. To find the impact of using the KP parameter (p) on the energy compaction property, the procedure given by [46] is employed. The stationary Markov sequence with length N and zero mean is analyzed. A matrix L with covariance coefficients (ρ) is defined as [27]: The matrix L is then transformed to the Krawtchouk domain. As such, the coefficients in the main diagonal of the transformed matrix (S) are computed. The matrix S represents the variance σ 2 l , which can be computed as where R denotes the KP matrix and (·) T refers to the matrix transpose operation. In addition, the normalized restriction error (J m ) can be computed using: where σ 2 q represents the order of σ 2 l sorted in descending order. In the experiment, the normalized restriction error is performed by considering different covariance coefficients (ρ) and different values of the parameters MNP. Figure 13 shows the normalized restriction error for different values of parameters (p) with the covariance coefficient ρ = 0.93. It can be observed from Figure 13 that when the value of p is equal to 1 − p, the normalized restriction error becomes equal. For example, the normalized restriction error for p = 0.05 is equal to p = 1 − 0.05 = 0.95. In addition, the energy compaction is influenced by the KP parameter p. For instance, as the parameter p increases from 0.05 to 0.45, the performance of the KP in terms of energy compaction is changed. Furthermore, the energy compaction at parameter p = 0.45 shows better performance for parameter p = 0.05 because the normalized restriction error (J m ) reaches zero values. However, small values of parameter p shows better performance in terms of feature extraction, as proven in [44]. Thus, it can be concluded that KP provides further performance improvement as the parameter p reaches 0.5. Furthermore, a more accurate result can be achieved when the parameter p is deviates from 0.5 [44].

Analysis of Reconstruction Error
In this section, the proposed recurrence algorithm is evaluated by carrying out reconstruction error analysis (REA). This REA was conducted for the proposed and the existing works. The REA was performed using an image formed from 16 well-known images as shown in Figure 14. In addition, the comparison was performed on an image with a large size, i.e., (4096 × 4096). Different values of parameter p were considered in the analysis. These values are p = 0.1, 0.2, 0.3, and 0.4. First, the WNKP (R) is generated using the proposed and existing algorithms. Then, the KMs (ψ) of the image are computed. Then, the image is reconstructed from the computed moments using a limited number of moments. Finally, the normalized mean square error (NMSE) is calculated between the input image and the reconstructed version of the image. Hence, the NMSE is given as [44] NMSE(I, I Rec ) = ∑ x,y (I(x, y) − I Rec (x, y)) 2 ∑ x,y (I(x, y)) 2 , where parameters I and I Rec denote the original image and the reconstructed image, respectively. The NMSE and the reconstructed image for p = 0.1 and p = 0.2 are shown in Figure 15 and Figure 16, respectively. The first row depicts the reconstructed images by utilizing the FRK [44], while the second row represents the reconstructed images using the proposed algorithm. The FRK algorithm is unable to reconstruct the image of the low-order moments. In addition, it is unable to reconstruct the image with high orders. However, the proposed algorithm is capable of fully reconstructing the image of different order values. In addition, the NMSE is minimized in the proposed algorithm using a moment order of 680, which is ∼ 16%. Figure 15 shows that the proposed algorithm achieves an NMSE of 0.72 while the FRK algorithm records a value of 0.84. This implies that the proposed algorithm outperforms the FRK algorithm. Moreover, the NMSE reaches zero when the proposed algorithm is used, while it records 0.64 when the FRK is used. The limitation with the FRK algorithm is due to the initial set that was computed, which makes the KPCs values tend towards zero, and thus, the NMSE is increased. Figure 17 provides a plot of experiments using p = 0.3. The results show that the proposed algorithm provides better NMSE than the FRK starting from a moment order of 680. In addition, the performance improvement increases when the moment order increases until it reaches the full order of (4096). At the full order, the NMSE using the proposed algorithm reaches 0, while it reaches 0.18 using the FRK algorithm. Figure 15. The NMSE performance comparing the proposed algorithm and the RAK algorithm [44] with p = 0.1. Figure 16. The NMSE performance comparing the proposed algorithm and the RAK algorithm [44] with p = 0.2. Figure 17. The NMSE performance comparing the proposed algorithm and the RAK algorithm [44] with p = 0.3. Figure 18 shows a new performance evaluation using p = 0.4. The results show that the proposed algorithm has the ability to accurately generate the KP coefficients while for the FRK, the KP coefficients remain inaccurate. This is attributed to the zero initial value obtained in the FRK algorithm. It is also worth noting that the proposed algorithm is able to generate KP coefficients for large polynomial sizes and at a high polynomial order, which was experimentally found to be greater than 8192. Figure 18. The NMSE performance comparing the proposed algorithm and the RAK algorithm [44] with p = 0.4. Finally, this paper provides a comparison of a maximum polynomial size between the proposed and existing algorithms. The maximum size is obtained for different values of the parameter p. For each recurrence algorithm, the polynomial (R) is computed first with a size and order of N × N. Then, the following criterion is applied: where Th is the threshold value over which 10 −5 is chosen. I(N) is the identity matrix with a size of N × N. It is worth noting that the maximum value of N is set to 20,480. Table 1 lists the obtained maximum size for each algorithm for different values of parameter p. This table demonstrates that the proposed algorithm is capable of generating an orthogonal polynomial with a large size and for different values of parameter p. The proposed algorithm outperforms the existing recurrence algorithms for all values of the parameters p considered. For example, at p = 0.05, the proposed algorithm generates a size of 82× larger than the RAN algorithm, 243× larger than the RAX algorithm, and 16× larger than the FRK algorithm. Thus, the proposed recurrence algorithm can be used to process large signals sizes quickly and accurately.

Computation of the Cost Analysis
The computation cost is considered an important factor to evaluate the performance of the proposed recurrence algorithm [44]. Thus, the computation cost of the proposed algorithm is compared with the existing works. These algorithms the are recurrence algorithm in the x direction (RAX), the recurrence algorithm in the n direction (RAN), and FRK. The computation cost is performed using the number of computed coefficients. Table 2 shows the ratio of computed coefficients (RCCs) using the proposed algorithm for different polynomials' sizes. Full moments' orders are considered. From Table 2, it can be observed that the RCC for small values of p and 1 − p achieves a low percentage. This percentage increases as the value of p increases towards 0.5. The average RCC for p = 0.05 and p = 0.95 is 9.22% while for p = 0.5, the average RCC is 20.35%. This is because the effective coefficient is shaping a circle for p = 0.5, and they are shaping a rotated ellipse as the parameter p deviates towards 0 or 1, which allows the number of effective coefficients to be reduced. Consequently, the percentage of the computed coefficients is reduced.  Table 3 demonstrates the performance improvement ratio between the proposed and the existing algorithms, namely RAN, RAX, and FRK. The RAN and RAX algorithms compute 50% of the KPCs for all values of parameter p because the nx plane is divided into two portions. On the other hand, the FRK algorithm computes 25% of the KPCs as it divides the nx plane into four partitions. The improvement ratio between the proposed and existing algorithms is obtained as Improvement Ratio = 1 − RCC of the proposed algorithm RCC of the existing algorithm .
According to Table 3, the proposed polynomial shows an improvement ratio of 18.64% at p = 0.5 when compared to the FRK algorithm and 59.32% when compared to the RAN and the RAX algorithm. The improvement ratio increases as the value of the parameter p deviates towards 0 and 1. For example, the results show that at parameter p = 0.25, the proposed method achieves an improvement ratio of 29.18% compared to the FRK algorithm, and 64.59% compared to other algorithms. In addition, a maximum improvement ratio of 63.11% is achieved by the proposed algorithm when it compares with the FRK algorithm while it is 81.55% when it compares with the RAN and RAX algorithms.

Conclusions
This paper described that state-of-the-art algorithms suffer from high error and that the implementation of these algorithms is limited to a specific value of parameter p. In addition to this, the state-of-the-art algorithms do not provide any computation of the initial value. Furthermore, to date, no algorithm is proposed for computing the KP coefficients with high-order polynomial and large polynomial sizes. To address these limitations, a new recurrence relation was proposed in this paper. To this end, a new initial value was presented and derived. In addition to this, a new diagonal recurrence relation was introduced. The proposed algorithm divided the KP plane into four triangles and only the coefficients in the upper triangle are computed. The coefficients in the upper triangle were divided into fundamental initial sets, initial sets, Part-1 and Part-2, respectively. The n-recurrence, x-recurrence, backwards x-recurrence, and diagonal recurrence relations were used to compute the values of the coefficients in the upper triangle. In addition, the identities were used to compute the values of the coefficients in the other triangles. The proposed algorithm was evaluated and compared with the existing recurrence algorithms. The comparison was carried out based on the reconstruction error, energy compaction, and computation cost. The experimental results showed that the proposed algorithm achieves a remarkable improvement over the existing algorithms in terms of the maximum size generated and the number of coefficients computed. Although the proposed algorithm outperforms state-of-the-art algorithms, the computational complexity is still high and can be further reduced. This can be achieved by implementing the proposed algorithm in a multi-processing environment (parallelizing) rather than in sequential form, as considered in this paper. Our future work is also directed towards implementing the proposed algorithm for KP with other orthogonal polynomials. This is expected to result in a new OP that has the potential of using orthogonal polynomials as well as the properties of KP coefficients.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Proof of Equation (10)
This section presents a proof through Equation (10): By taking the exponential and natural log (e log ), then (A1) can be simplified to: