Chebyshev approximation

1. Chebyshev approximation

Pafnuty Lvovich Chebyshev was a Russian mathematician and considered to be the founding father of Russian mathematics. Chebyshev is known for his fundamental contributions to the fields of probability, statistics, mechanics, and number theory [1].

Chebyshev was the first who came up with the idea of approximating functions in the uniform norm. He asked himself if it is possible to represent a continuous function $f(x)$ on the closed interval $[a, b]$ by a polynomial in such a way that the maximum error at any point $x \in [a, b]$ is controlled? Or in other words, is it possible to construct a polynomial so that the error is minimized [2]?

The key to using polynomials to evaluate functions is not to think of polynomials as being composed of linear combinations of $1, x, x^2, ...$, but as linear combinations of Chebyshev polynomials [3].

Without further ado, let's just dive into the Chebyshev polynomial of the first kind. The recursive expression for generating Chebyshev polynomial of the $n$-th order is [4]:

$$T_0(x) = 1, \quad T_1(x) = x, \quad T_{n}(x) = 2 x T_{n-1}(x) - T_{n-2}(x). \tag{1}$$

This doesn't say much, but if we plot the first five polynomials (up to order 4) in the range from -1 to 1 we can see that the functions resemble the periodic behavior of the Fourier series.

The Chebyshev approximation is the linear superposition of Chebyshev polynomials with coefficients $a_i$ to the polynomial of order $i$. With known coefficients, we can approximate a function as:

$$f(x) \approx \sum_{i=0}^n a_i T_i(x)$$

For practical purposes, one can calculate directly the Chebyshev polynomials along with the calculation of the approximating value. The overhead of calculating the Chebyshev polynomials is small regarding the ordinary polynomial calculation. We define the function approximation, for a coefficient list $a$, in the following Python snippet:

def getChebyApprox(coeff : list, x: float):
    T0 = 1.0
    T1 = x
    retVal = T0*coeff[0] + T1*coeff[1]
    for i in range(2, len(coeff)):
        Ti = 2*x*T1 - T0
        retVal += Ti*coeff[i]
        T0 = T1
        T1 = Ti

    return retVal

The function is meant to be used for x-values in the range $[-1,1]$. For an arbitrary range, one can easily transform the x-axis into a new axis that is of the defined range of $[-1, 1]$. The affine transformation is given in the next paragraph.

1.1. Extending the approximating range

Although the properties of the Chebyshev polynomials are defined on the range [-1, 1], the approximation can be used for any range by transforming the variable $x$. If we want to approximate function in the range [a, b], we can use $u$ instead of $x$:

$$u = x\frac{2}{b - a} - \frac{a + b}{b - a}, \quad x \in [a,b] \quad \Rightarrow \quad u \in [-1, 1].\quad$$

With the relation that transforms the axis back to $x$:

$$x = u\frac{b-a}{2} + \frac{a + b}{2}$$

We can write a simple function that returns the factor and offset for the affine transformations from x to u and vice versa.

def getAffineValues(a : float = -1, b : float = 1):
    facx = 2/(b-a)
    offx = -(a+b)/(b-a) 
    facu = 0.5*(b - a)
    offu = 0.5*(a + b)
    return facx, offx, facu, offu

2. Orthogonality

One of the most important properties of Chebyshev polynomials and Fourier series is that they are both orthogonal functions. Orthogonality is a powerfull property that simplifies the discovery of individual components in a linear superposition of functions. Mathematically orthogonality of the Chebyshev polynomials is expressed as:

$$\int_{-1}^{1} T_i(x) T_j(x) \frac{dx}{\sqrt{1-x^2}} = \left\{ \begin{aligned} &0 \quad &i \neq j, \\ &\pi \quad &i = j = 0, \\ &\pi/2 \quad &i = j \neq 0. \end{aligned} \right.$$

Probably you can't see how important the above integral is yet, but let's see how can we use it to find the coefficients used to approximate an arbitrary function $y(x)$. If we replace the approximating sum in the orthogonality expression we get:

$$\int_{-1}^{1} y(x) T_j(x) \frac{dx}{\sqrt{1-x^2}} \approx \sum_{i=0}^n a_i \int_{-1}^{1} T_i(x) T_j(x) \frac{dx}{\sqrt{1-x^2}}.$$

It is now obvious that the integral extracts exactly specific components ($a_k$) from the approximating sum since it non-zero only when $k=m$. If we rewrite the above expression we get formula for the coefficients of the approximation:

$$\begin{aligned} a_0 &= \frac{1}{\pi} \int_{-1}^{1} y(x) T_0(x) \frac{dx}{\sqrt{1-x^2}} \\ a_i &= \frac{2}{\pi} \int_{-1}^{1} y(x) T_i(x) \frac{dx}{\sqrt{1-x^2}}. \end{aligned}$$

The same property is used for calculating coefficients of the Fourier transform. We can use simple numerical integration for practical purposes. Implemented in Python, the function that returns the Chebyshev coefficients for a callable or sampled function $f$, can be written as:

import math

def getChebyCoeffIntegral(xList : list, f, order: int, a: float = -1, b: float = 1):
    facx, offx, _, _ = getAffineValues(a, b)
    polyCnt = order + 1
    coeff = [0]*polyCnt
    for i in range(polyCnt):
        integral = 0
        ai = [0]*polyCnt
        ai[i] = 1
        for j in range(1, len(xList)):
            x = 0.5*(xList[j-1] + xList[j])
            u = x*facx + offx
            du = (xList[j] - xList[j-1])*facx            
            Ti = getChebyApprox(ai, u)
            # Sampled function
            if type(f) == list:
                y = 0.5*(f[j-1] + f[j])
            # Callable function
            else: 
                y = f(x)
            integral += Ti*y*du/math.sqrt(1 - u**2)
        coeff[i] = integral/math.pi
        if i > 0:
            coeff[i] *= 2

    return coeff

2.1. Chebyshev nodes

Chebyshev nodes are the roots of the Chebyshev polynomials [5]. Node values are given with the following formula, for a polynomial of order $n$:

$$x_k = \cos\left(\frac{2k-1}{2n}\pi\right), \quad k = 0, 1, ..., n-1.$$

For the order of n=4, we can represent the Chebyshev nodes on the existing figure as above:

2.2. Discrete orthogonality

Chebyshev polynomials have another usefull property, they are orthogonal at the Chebyshev nodes ($x_k$). The discrete orthogonality is expressed as:

$$\sum_{k=0}^{n-1}T_i(x_k)T_j(x_k) = \left\{ \begin{aligned} &0 \quad &i \neq j, \\ &n \quad &i = j = 0, \\ &n/2 \quad &i = j \neq 0. \end{aligned} \right.$$

The discrete orthogonality can be also used to approximate the coefficients $a_k$ by calculating the function value only at the Chebyshev nodes. It is an approximation of the integral form given in the previous section. If we substitute the function $y(x_k)$ instead of $T_i(x_k)$, we can derive the expression for approximately calculating the Chebyshev components.

$$\sum_{k=0}^{n-1}y(x_k)T_j(x_k) \approx \sum_{i=0}^{n}a_i\sum_{k=0}^{n-1}T_i(x_k)T_j(x_k).$$

By rearranging the above expression we get a discrete formula for calculating the Chebyshev components:

$$\begin{aligned} a_0 = \frac{1}{n}\sum_{k=0}^{n-1}y(x_k)T_0(x_k) \\ a_i = \frac{2}{n}\sum_{k=0}^{n-1}y(x_k)T_i(x_k) \end{aligned}$$

Let's write the function for getting the Chebyshev coefficients using the discrete orthogonality, for a callable function $f$:

import math

def getChebyCoeffNodes(f, order: int, a: float = -1, b: float = 1):
    polyCnt = order + 1
    coeff = [0]*polyCnt
    _, _, facu, offu = getAffineValues(a, b)
    ukList = [math.cos(math.pi*(2*k + 1)/(2*order)) for k in range(order)]
    for i in range(polyCnt):
        sum = 0
        ai = [0]*polyCnt
        ai[i] = 1 
        for k in range(order):
            uk = ukList[k]
            xk = uk*facu + offu
            Ti = getChebyApprox(ai, uk)
            y = f(xk)
            sum += Ti*y
        coeff[i] = sum/order
        if i > 0:
            coeff[i] *= 2

    return coeff

3. Chebyshev to ordinary polynomial coefficients

Transforming Chebyshev polynomials to ordinary polynomials can be easily performed. This transformation doesn't have much practical use, since the number of instructions on a machine will remain practically the same. The ordinary polynomial can be useful for representing the Chebyshev approximation in an ordinary form. The derivation of the method is cumbersome because it involves matrix arithmetic, but it can be stated with the following matrix expression:

$$\begin{bmatrix} 1 & 0 & \!\!\!\!-1 & 0 & 1 & 0 & \!\!\!\!-1 \\ 0 & 1 & 0 & \!\!\!\!-3 & 0 & 5 & 0 \\ 0 & 0 & 1 & 0 & \!\!\!\!-4 & 0 & 9 \\ 0 & 0 & 0 & 1 & 0 & \!\!\!\!-5 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & \!\!\!\!-6 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ a_4 \\ a_5 \\ a_6 \end{bmatrix} = \begin{bmatrix} b_0 \\ b_1 \\ b_2/2 \\ b_2/4 \\ b_4/8 \\ b_5/16 \\ b_6/32 \end{bmatrix}$$

In the above expression Chebyshev coefficients are named with $a_i$ and the ordinary coefficients are named $b_i$. The approximation can be written both ways:

$$f(x) = \sum_{i=0}^n a_iT_i(x) = \sum_{i=0}^n b_ix^i.$$

The conversion matrix is upper triangular of size $n \times n$, and has non-zero coefficients $c_{i,j}$ given by the algorithm:

$$\begin{aligned} c_{i,j} &= 1, \quad &i=j \\ c_{1,j} &= 0, \quad &j=2,4,6,... \\ c_{1,j} &= -c_{1,i-2}, \quad &j=3,5,7,... \\ c_{2,j} &= -\mathrm{sgn}(c_{2,j-2})\left( \lvert c_{2,j-2} \rvert + 2 \right), \quad &j=4,6,8,... \\ c_{i,j} &= -\mathrm{sgn}(c_{i,j-2})\left( \lvert c_{i,j-2} \rvert + \lvert c_{i-1,j-1} \rvert \right), \quad &i \geq 3, \quad j>i, \quad j=5,7,9,... \\ \end{aligned}$$

This algorithm can be easily written in Python by exploiting numpy for matrix manipulation:

import numpy as np

def convertChebyshevCoeff(aCoeff: list):
    n = len(aCoeff)
    cMatrix = np.zeros((n, n), dtype=int)
    for i in range(n):
        cMatrix[i, i] = 1
    for i in range(2, n, 2):
        cMatrix[0, i] = -cMatrix[0, i-2]
    for i in range(3, n, 2):
        cMatrix[1, i] = -np.sign(cMatrix[1, i-2])*(2 + np.abs(cMatrix[1, i-2]))
    for i in range(2, n):
        for j in range(i+2, n, 2):
            cMatrix[i, j] = -np.sign(cMatrix[i, j-2])*(np.abs(cMatrix[i, j-2]) + np.abs(cMatrix[i-1, j-1]))
    cNorm = [1]*n
    for i in range(2, n):
        cNorm[i] = 2**(i-1)

    return np.matmul(cMatrix, aCoeff)*cNorm

The conversion technique is taken from the paper [6]:

Maryam Shams Solary (2018) Converting a Chebyshev polynomial to an ordinary polynomial with series of powers form, Journal of Interdisciplinary Mathematics, 21:7-8, 1519-1532, DOI: 10.1080/09720502.2017.1303949

4. Practical use of the Chebyshev approximation

The use of Chebyshev approximation for function calculation has its pros and cons. These depend mainly on the type of the function, and the possibility of it being easily approximated by a polynomial.

Pros:

• Speed: Chebyshev approximation is a fast method for machine computation since it uses a relatively small number of two operations: multiplications and adding. The speeds of computation up to a polynomial order of 6, are comparable to the speed of a lookup table with linear interpolation.
• Memory: Memory requirements are extremely small - only for the approximation coefficients (in opposite to a lookup table)
• Error: Approximation error is evenly spread in the range [-1, 1], so that a Chebyshev approximation avoids the Runge's penomenon. Runge's phenomenon is the oscillation at the edges of an interval that occurs when using polynomial interpolation with polynomials of a high degree over a set of equidistant interpolation points.
• Almost perfect: Chebyshev approximation is close to the best polynomial approximation to a continuous function under the maximum norm, also called the "minimax" criterion.

Cons:

• Not suitable for all functions: There are functions that can't be represented well enough with a polynomial of low order. These functions manifest the non-polynomial property by a small decrease in the Chebyshev polynomial coefficients with order increment. A rule of thumb is that a decrement of an order of magnitude is a good indicator that a function can be approximated with a Chebyshev polynomial.

  For example, the first function is more suitable for approximating than the second function:

  f(x)	range	$a_0$	$a_1$	$a_2$	$a_3$	$a_4$	$a_5$
  $\exp(x)$	[0,1]	1.7534	0.85039	0.10521	0.0087221	0.00054231	-4.6926e-16
  $\frac{1}{1+x^2}$	[-3,3]	0.30404	0	-0.29876	0	0.12222	0

5. Examples

Let's examine the approximations of two functions given above, with the implementations in Python as described:

$$\begin{aligned} f_1(x) &= \exp(x), \quad x \in [0,1], \\ f_2(x) &= \frac{1}{1+x^2}, \quad x \in [-3, 3]. \end{aligned}$$

5.1. Function with a good polynomial property

The function $\exp(x)$ can be nicely approximated with a polynomial. Here we investigate the two methods for calculating the coefficients, from integral and from nodes.

a = 0
b = 1
sampleSize = 100000
dx = (b-a)/(sampleSize - 1)

xList = [i*dx - a for i in range(sampleSize)]
f1List = [math.exp(x) for x in xList]

facx, offx, _, _ = getAffineValues(a, b)

a1ListInt = getChebyCoeffIntegral(xList, math.exp, 5, a, b)
a1ListNds = getChebyCoeffNodes(math.exp, 5, a, b)

f1ApproxIntList = [getChebyApprox(a1ListInt, x*facx + offx) for x in xList]
f1ApproxNdsList = [getChebyApprox(a1ListNds, x*facx + offx) for x in xList]

Coefficients from integral calulation 
    a = 1.7511, 0.8483, 0.10068, 0.0066296, -0.0039847, -0.0020655
Coefficients from nodes calulation 
    a = 1.7534, 0.85039, 0.10521, 0.0087221, 0.00054231, -4.6926e-16

We see that the approximation in both cases gives good results. For more information about the quality, an error of the approximation is more suitable to look at. From figure 4, we can conclude that the nodes calculation gives much better results for continuous functions.

The approximation error for the integral calculation of coefficients can be increased by increasing the number of sampling points for the calculation. We see in figure 5 that if we increase the number of samples for the numerical integration, the coefficients become better for approximating the function. We can also notice that the coefficients decrease more rapidly for better approximations.

Coefficients from integral calculation (100000 samples)
    a = 1.7511, 0.8483, 0.10068, 0.0066296, -0.0039847, -0.0020655
Coefficients from integral calculation (500000 samples)
    a = 1.7524, 0.84946, 0.10318, 0.0077863, -0.0014816, -0.00090867
Coefficients from integral calculation (2500000 samples)
    a = 1.7529, 0.84997, 0.1043, 0.0083036, -0.00036216, -0.00039138

5.2. Function with a bad polynomial property

The function $\frac{1}{1+x^2}$ is not suitable for polynomial approximation. We can use the same code as for the previous function. When we use the polynomial of order 5, the approximation has a big error that can be seen in figure 6. One solution is to increase the polynomial order. For order 20 we get a good enough approximation but the order is so high that it doesn't make sense anymore to use an approximation at all.

def funcF2(x) -> float:
    return 1/(1 + x**2)

a = -3
b = 3
sampleSize = 100000
dx = (b-a)/(sampleSize - 1)

xList = [i*dx + a for i in range(sampleSize)]
f2List = [funcF2(x) for x in xList]

facx, offx, _, _ = getAffineValues(a, b)

a2List5Nds = getChebyCoeffNodes(funcF2, 5, a, b)
f2Approx5NdsList = [getChebyApprox(a2List5Nds, x*facx + offx) for x in xList]

a2List20Nds = getChebyCoeffNodes(funcF2, 20, a, b)
f2Approx20NdsList = [getChebyApprox(a2List20Nds, x*facx + offx) for x in xList]

Coefficients from nodes calulation (order 5)
    a = 0.3411, 2.2204e-17, -0.38935, -6.1062e-17, 0.26955, 4.6818e-17
Coefficients from nodes calulation (order 20) 
    a = 0.31623, 3.4694e-17, -0.32855, -1.6653e-17, 0.17068, -1.6653e-17, -0.088659, -1.6653e-17, 0.046045, 5.1348e-17, -0.023895, -5.5511e-18, 0.012365, -6.5226e-17, -0.006331, 2.9837e-17, 0.0031105, 8.4655e-17, -0.0012725, -8.1532e-17, 1.673e-17