Stable Calculation of Krawtchouk Functions from Triplet Relations

Abstract

Deployment of the recurrence relation or difference equation to generate discrete classical orthogonal polynomials is vulnerable to error propagation. This issue is addressed for the case of Krawtchouk functions, i.e., the orthonormal basis derived from the Krawtchouk polynomials. An algorithm is proposed for stable determination of these functions. This is achieved by defining proper initial points for the start of the recursions, balancing the order of the direction in which recursions are executed and adaptively restricting the range over which equations are applied. The adaptation is controlled by a user-specified deviation from unit norm. The theoretical background is given, the algorithmic concept is explained and the effect of controlled accuracy is demonstrated by examples.

1. Introduction

Polynomial transformations are well-defined mathematical concepts, with the classical (continuous and discrete) orthogonal polynomials standing out as having a common basis and sharing multiple features. In the digital age, the discrete polynomials are especially popular and those with a finite support (discrete Chebyshev and Krawtchouk) have been used in various signal processing tasks, but dominantly for image processing. Examples of applications of Krawtchouk polynomials can be found in, e.g., [1,2,3,4,5,6,7].

The recurrence relation and difference equation associated with these polynomials seem to constitute a simple recipe for their numerical evaluation. However, the recursive character of these generating schemes makes them vulnerable to error accumulation, which, obviously, becomes increasingly pronounced at higher polynomial degrees and/or larger supports. For many applications, i.e., cases of small support and or low polynomial degree, this issue is of no consequence, and the attention is more focused on computational simplification rather than accuracy. In this paper, we focus on the fundamental issue of accuracy, not so much on the issue of computational parsimony; in fact, extra computational burden is introduced in order to control the accuracy in the function generation process.

For the normalized discrete Chebyshev polynomials several mitigation schemes have been proposed; the latest one a computation scheme characterized by truncation of the iterations under the control of a user-defined deviation from unit norm [8]. The solutions defined for the discrete Chebyshev polynomials do not carry over since they make use of the (anti-)symmetry of the polynomials with respect to the mid of the domain, i.e., a property that does not generally hold for Krawtchouk polynomials. Nevertheless, a first mitigation was proposed in [9] building upon the notion that the direction of execution of the relations (i.e., decreasing/increasing index and polynomial degree) should be carefully chosen in order to reduce the error propagation. However, not all symmetries were exploited, nor does the scheme allow direct control over the accuracy like in [8].

This paper describes a method for generating the Krawtchouk functions where the deviation of generated functions from unit norm is user-controlled. The paper is organized as follows. Section 2 recaps the knowledge on the Krawtchouk polynomials relevant for the proposed algorithm and Section 3 introduces the Krawtchouk functions and the pertinent properties. It also introduces the well-known concepts of center of energy and effective width, notions underlying the proposed scheme. Having put all bits and pieces in the same notational framework, we are able to present the novel algorithm and its performance. The main contribution of this paper is turning the three-term relations into a practical numerical evaluation scheme. Section 4 explains the algorithmic concept and Section 5 illustrates the behavior of the algorithm in a numerical sense. We close with the conclusions.

2. Krawtchouk Polynomials

The Krawtchouk polynomials form an orthogonal set of polynomials on a discrete interval $[0,\cdots ,N]$. We denote the polynomials as $Q_n(k)$, where n is the polynomial degree, $0\le n\le N$. They are orthogonal with respect to a weight function $v(k)$ given by

$$v(k)=\left(\genfrac{}{}{0pt}{}{N}{k}\right)p^k{q}^{N-k}, \mathrm{with} p,q>0, p+q=1.$$

(1)

For the weight function there holds that the ratio of two consecutive weights is the ratio of two first-degree polynomials:

$$r(k)=\frac{v(k)}{v(k-1)}=\frac{-pk+p(N+1)}{qk}.$$

(2)

2.1. Rodrigues’ Formula

The polynomials are expressed by

$$Q_n(k)=\frac{{(-1)}^{n}k!(N-k)!}{n!p^k{q}^{N-k}}{\mathsf{\Delta }}^n\left[\frac{p^k{q}^{N-k}}{(k-n)!(N-k)!}\right],$$

(3)

for $n=0,1,\cdots ,N$ where $\mathsf{\Delta }$ denotes the forward difference operator.

2.2. Difference Equation

The Krawtchouk polynomials satisfy the second-order difference equation

$$p(N-k-1){\mathsf{\Delta }}^2{Q}_n(k)+[-k+pN-1+n]\mathsf{\Delta }{Q}_n(k)+n{Q}_n(k)=0.$$

Expanding the forward differences yields the triplet relation

$$p(N-k-1)Q_n(k+2)+[(p-q)(k+1)-pN+n]Q_n(k+1)+[q(k+1)]Q_n(k)=0.$$

(4)

2.3. Leading Polynomial Coefficients

Let the polynomial $Q_n(x)$ be defined as

$$Q_n(x)=\sum _{m=0}^n{a}_{n,m}{x}^m.$$

$a_{n,n}$ and $a_{n,n-1}$ are the two leading polynomial coefficients given by

$$a_{n,n}=\frac{1}{n!} \mathrm{and} a_{n,n-1}=-\frac{(n-1)/2+p(N+1-n)}{(n-1)!}.$$

2.4. Norms

The norm of the polynomials under the given weight function is given by

$$h_n^2=\left(\genfrac{}{}{0pt}{}{N}{n}\right)p^n{q}^n$$

(5)

and we take $h_n$ as a positive number.

2.5. Recurrence Relation

The recurrence relation is given by

$$Q_{n+1}(x)=\{\frac{1}{n+1}x+\frac{-p(N-n)+qn}{n+1}\}{Q}_n(x)-pq\frac{N-n+1}{n+1}{Q}_{n-1}(x).$$

(6)

3. Krawtchouk Functions

We define the Krawtchouk functions as the normalized product of the polynomials and the square-root of the weight function [10,11]:

$${\kappa }_n(k)=Q_n(k)\sqrt{v(k)}/h_n.$$

(7)

The set of Krawtchouk functions form an orthonormal basis in an $N+1$-dimensional Euclidian space.

3.1. Difference Equation

From (7) and (4) with

$$\frac{v(k+1)}{v(k+2)}=\frac{q(k+2)}{p(N-k-1)} \mathrm{and} \frac{v(k+1)}{v(k)}=\frac{p(N-k)}{q(k+1)}$$

there follows

$$\begin{array}{ccc}& & \sqrt{pq(N-k)(k+1)}{\kappa }_n(k+1)+[n-p(N-k)-qk]{\kappa }_n(k)+\hfill \\ & & \sqrt{pq(N-k+1)k}{\kappa }_n(k-1)=0.\hfill \end{array}$$

(8)

For convenience we write this as

$$D_1(k){\kappa }_n(k+1)+[n-D_0(k)]{\kappa }_n(k)+D_{-1}(k){\kappa }_n(k-1)=0.$$

3.2. Recurrence Relation

From

$$\frac{h_n}{h_{n+1}}=\sqrt{\frac{n+1}{N-n}}\frac{1}{\sqrt{pq}}$$

and the recurrence relation of the polynomials (6) there follows that the discrete orthonormalized Krawtchouk functions ${\kappa }_n$ satisfy the recurrence relation

$${\kappa }_{n+1}(x)=(A_{n}x+B_n){\kappa }_n(x)-C_n{\kappa }_{n-1}(x)$$

(9)

with

$$\begin{array}{ccc}\hfill A_n& =& \frac{1}{\sqrt{(n+1)(N-n)pq}},\hfill \\ \hfill B_n& =& -\frac{p(N-n)+qn}{\sqrt{(n+1)(N-n)pq}},\hfill \\ \hfill C_n& =& \sqrt{\frac{N-n+1}{N-n}\frac{n}{n+1}}.\hfill \end{array}$$

We note that restructuring of (8) and (9) reveals that the relation between function values at three consecutive k-values (the difference equation) is essentially the same as that for three consecutive orders n (the recurrence relation), except for a systematic sign difference. In particular, we have

$$\begin{array}{ccc}\hfill 1/A_k& =& D_1(k),\hfill \\ \hfill B_k/A_k& =& D_0(k),\hfill \\ \hfill C_k/A_k& =& D_{-1}(k).\hfill \end{array}$$

We introduce the name triplet relations as a generic name for the Equations (8) and (9).

3.3. Symmetries

There are three symmetries in the Krawtchouk functions

$$\begin{array}{ccc}\hfill {\kappa }_n(k)& =& {(-1)}^{k-n}{\kappa }_k(n),\hfill \\ \hfill {\kappa }_n(k)& =& {\kappa }_{N-k}(N-n),\hfill \\ \hfill {\kappa }_n(k)& =& {(-1)}^{k-n}{\kappa }_{N-n}(N-k).\hfill \end{array}$$

Two of these three symmetries are independent, and the third one is implied by the other two. Roughly speaking, this means that the Krawtchouk functions need to be calculated on a (well-chosen) quarter of the whole $k,n$ plane; mirroring, with possibly sign reversal, can be used to obtain the function values at the remaining $k,n$ positions.

3.4. Special Case

There is one special parameter choice, namely $p=0.5$: this setting introduces another axis of (anti-)symmetry expressed as

$${\kappa }_n(k)={(-1)}^n{\kappa }_n(N-k).$$

This is a symmetry also present in the discrete Chebyshev polynomials and can be exploited to define proper starting points for deploying difference or recurrence relations in combination with a norm-based stop criterion, thus mitigating numerical issues [8].

3.5. Limiting Cases

Limiting case are the boundaries of the range allowed for p, where

$$\begin{array}{ccc}\hfill \underset{p\to 0}{lim}{\kappa }_n(k)& =& {\delta }_{k-n},\hfill \\ \hfill \underset{p\to 1}{lim}{\kappa }_n(k)& =& {(-1)}^{k-n}{\delta }_{N-k+n},\hfill \end{array}$$

with $\delta$ the Kronecker delta.

3.6. Center of Energy and Effective Width

The fact that at well-defined parts of the $k,n$ plane, the function $\kappa$ takes on large values while at other parts it attains only small values, is a basic ingredient for our algorithm. For insight and later deployment, we introduce the well-known concepts of center of energy and effective width of a function.

Definition 1.

The center of energy $C_f$ of a function f on a discrete domain is defined as

$$M_f=\frac{{\sum }_{x}x{f}^2(x)}{{\sum }_x{f}^2(x)}.$$

Definition 2.

The effective width of the energy $W_f$ of a function f on a discrete domain is defined as

$$W_f^2=\frac{{\sum }_x{(x-M_f)}^2{f}^2(x)}{{\sum }_x{f}^2(x)}.$$

Theorem 1

([10]). For an orthonormal system $\{g_n\}$, $n=0,1,2,\cdots$, in ${\ell }^2$ defined by a 3-term recurrence relation

$$(x-{\gamma }_n)g_n(x)={\alpha }_n{g}_{n+1}(x)+{\beta }_n{g}_{n-1}(x)$$

with $g_{-1}=0$, the center of energy of the function $g_n$ is given by

$$M_{g_n}={\gamma }_n$$

and for its squared effective width holds

$$W_{g_n}^2={\alpha }_n^2+{\beta }_n^2.$$

In the remainder, we will use the notations $M_n$ and $W_n$ as the center and effective width for the nth order Krawtchouk function. In particular we have:

$$M_n={\gamma }_n=p(N-n)+qn$$

and for its squared effective width holds

$$W_n^2={\alpha }_n^2+{\beta }_n^2=(1+C_n^2)/A_n^2=pq[N+2n(N-n)].$$

In view of the symmetry relations, similar equations evolve when considering energies along the function order (i.e., at constant k). Examples of lines of the center of energy and the areas of energy concentration as defined by the effective width are provided in Figure 1.

3.7. Parameter Symmetry

The parameters p and q ($q=1-p$) are almost interchangeable. To make the dependence on p explicit, we denote the Krawtchouk functions as ${\kappa }_n(k;p)$. We then have

$${\kappa }_n(k;p)={(-1)}^n{\kappa }_n(N-k;1-p).$$

4. Controlled Accuracy Strategy

Calculation of the Krawtchouk functions using the recurrence relation and/or the difference equation runs into numerical problems when executed (far) outside of the effective width. The numerical inaccuracies occur at the position where the function is nearly zero, and an effective strategy is to simply set these values to zero. However, we want to carry this out in a way such that we have control over the attained approximation; in particular, the proposed scheme tries to truncate operations such that the deviation of the norm is controlled.

In case of $p=q=0.5$, an effective strategy is to run the recurrence relation over two values of k in the mid. Next, from the middle the difference equation in combination with the symmetry can be deployed to determine, for a certain order, the function at decreasing and increasing values of k until the sum of squares (norm) is sufficiently close to unity. The function at values of k which have not been iterated yet are set to zero. This is akin to the strategy for the discrete Chebyshev polynomials [8].

This strategy does not apply for other values of p, as the middle does not have this symmetry and may not be a good starting point for the recursion. The concept that is applied is based on the earlier introduced concept of center of energy. Essentially, it is balancing act: we know the final center of energy, and we can determine the center of energy associated with the samples determined so far. The range of k values is extended in the direction opposite to where the current center of gravity is relative to the final one.

In more detail, the procedure works as follows. For each function order n, the initialization starts at the two points of k closest to the center of energy $M_n$. We now extend the range of k-values either at the lower or higher border of the existing range, thus increasing the range over which samples has been calculated. The side at which the range is increased is selected based on the center of gravity of the function values so far calculated: if the center of gravity is lower than the ultimate center of gravity, the border is increased at high end range of already addressed values of k; otherwise, the range of k-values is extended at the lower boundary. Initialization at each n is carried out using the recurrence relation with function values at the lower orders $n-1$ and $n-2$. The whole procedure starts at $n=0$ by the determination of the square-root of the binomial distribution for the two points of k closest to $M_0$—i.e., around the peak of the binomial distribution defined by N and p.

5. Results

In this section, we will first look into the generation of the zeroth-order Krawtchouk function specifically considering the control (user-specified and attained deviation from norm) and number of samples that need to be calculated. Next, we consider the more general case and experimentally assess the maximal deviations from the norm over the entire basis as well as maximum deviation from orthogonality.

5.1. Zeroth-Order Function and Number of Non-Zero Elements

To gain insight into the effect of the control mechanism, we consider how the truncation affects the calculation of the zeroth-order function. We show two different segment sizes represented by $N=200$ and 2000; we refrain from looking at small sizes because the effects there are minimal (see Introduction): the algorithm is intended for large segments and/or high function orders. We take five different values of p between 0 and 1/2; the effects for the upper half of the p range is implied due to the symmetry in p and q.

In Figure 2, the attained error is plotted as a function of the specified error in the norm of the zeroth-order function. It shows that the attained deviation indeed does not exceed the specified error independent of parameter settings N and p.

The fraction of samples of the zeroth-order Krawtchouk function that are not calculated by the triplet relation but simply set to zero is shown in Figure 3. It is clear that the segment size and p-value have a large influence: increasing N or taking p further from 0.5 increases the relative number of samples actually skipped in the calculation. Especially at larger N, the amount of calculations not executed is very high: for $N=2000$ and an $ϵ$ close to machine accuracy (${10}^{-13}$) less than 20% of the samples need to be calculated; for p close to the range boundaries the savings are even larger.

5.2. Deviations from Norm and Orthogonality

By definition, the algorithm holds the deviation from norm in check. The maximum of the attained deviation over all function orders at a certain parameter settings is always slightly less than the specified $ϵ$. This is shown in Figure 4.

However, there is no explicit control over the orthogonality. It is therefore of interest to check if this is well-behaved in some sense. We experimented with this in the same way as before: we took various parameter settings (two values for N, and five values for p) and determined the maximum deviation from orthogonality for a range of stopping criteria. The results are shown in Figure 5. Not unexpectedly, it was found that the maximum scales with $\sqrt{ϵ}$: the dashed line in the figure. The maximum deviation is found to be always less than $\sqrt{ϵ}$: more detailed inspection showed that the maximum non-orthogonality is close to $\sqrt{ϵ}/3$.

6. Conclusions

Calculation of the Krawtchouk functions using triplet relations (recurrence relation and difference equation) is prone to error accumulation, resulting in large deviations from the actual data for already relatively small support and function order. An algorithm is proposed to solve this issue. It is controlled by a parameter specifying the allowed deviation from unit norm. The theoretical background is given, the algorithmic concept is explained and the effect of controlled accuracy is demonstrated by examples. The experiments demonstrate that not only the norm is controlled, but the maximum absolute deviation from orthonormality is held in check as well.