Accurate and Efficient Evaluation of Chebyshev Tensor Product Surface

A Chebyshev tensor product surface is widely used in image analysis and numerical approximation. This article illustrates an accurate evaluation for the surface in form of Chebyshev tensor product. This algorithm is based on the application of error-free transformations to improve the traditional Clenshaw Chebyshev tensor product algorithm. Our error analysis shows that the error bound is u + O(u2) × cond(P, x, y) in contrast to classic scheme u × cond(P, x, y), where u is working precision and cond(P, x, y) is a condition number of bivariate polynomial P(x, y), which means that the accuracy of the computed result is similar to that produced by classical approach with twice working precision. Numerical experiments verify that the proposed algorithm is stable and efficient.


Introduction
Chebyshev polynomials have been extended to almost all mathematical and physical discipline, including spectral methods, approximation theory, and representation of potentials [1][2][3][4].Bivariate Chebyshev polynomials have gained attention of the computer vision researchers [5,6].Over the years, researchers have focused on the implementation of Chebyshev tensor product series in image analysis [5][6][7].The Chebyshev tensor product series can be used to approximate an image, which is essentially regarded as a two-dimensional spatial function [8].Two separable univariate Chebyshev polynomials that are discrete and orthogonal can approximate two-dimensional signal.Mukundan et al. [5] introduce a new discrete Chebyshev tensor product based on Chebyshev polynomials.This discrete Chebeshev tensor product functions show the effectiveness as feature descriptors.Rahmalan et al. [6] propose a novel approach based on discrete orthogonal Chebyshev tensor product for an efficient image compression.Recently, Omar et al. [7] propose a novel method for fusing images using Chebyshev tensor product series.All above need an image reconstruction from a finite Chebyshev tensor product surface.Thus, developing fast and reliable algorithms to evaluate the Chebyshev tensor product series are of challenging interest [9].The Clenshaw tensor product algorithm (CTP) [10] is one of algorithms that are used to evaluate Chebyshev tensor product series.
In order to get a high-precision approximation of an image, it is essential to evaluate the series accurately.Particularly, we require higher level of accurate numeric results for ill-conditioned cases.Li et al. 's double-double [11] (doubledouble numbers are represented as an unevaluated sum of a leading double and a trailing double) is a library used to improve the accuracy of numerical computation.However, the algorithm is time-consuming when an input image becomes larger.
Error-free transformation studied by Ogita et al. [12] is another direct possible method to improve the accuracy apart from increasing the working precision.Compensated algorithms to evaluate the univariate polynomials in different basis have been proposed in [12][13][14][15].Inspired by their work, we extend the univariate compensated algorithm to tensor product case using the compensated Clenshaw algorithm [16] for evaluation of Chebyshev series [15].We perform a com-pensated Clenshaw tensor product algorithm (for simplicity we denote it by CompCTP algorithm) to evaluate the polynomials expressed in Chebyshev tensor product form.The proposed algorithm produces the same accuracy as using twice working precision.
Since the image is fundamentally treated as two-dimensional spatial function, we use general two-dimensional function to illustrate our algorithm in the sequel.This paper has the following layout.In Section 2, we introduce some preliminaries and basic algorithms underlying our algorithm.In Section 3, we propose the compensated algorithm to compute surface in form of Chebyshev tensor product.In Section 4, we analyze forward error bound of the algorithm.In Section 5, a series of numerical experiments illustrate the accuracy and efficiency of the proposed algorithm.

Preliminaries and
Error-Free Transformations 2.1.Basic Notations.At the present time, IEEE 64-bit floating arithmetic standard is implemented, which is sufficiently accurate for most scientific applications.Throughout the paper, we presume that all the computations are performed using IEEE-754 [17] standard in double precision so that neither overflow nor underflow occurs.We assume that the computations are produced in a floating-point arithmetic which yields the models where op ∈ {+, −, ×, /} and  is the working precision.For brevity we denote fl( op ) =  o , o ∈ {⊕, ⊖, ⊗, ⊘}.

Accurate Algorithm to Evaluate Chebyshev Tensor Product Surface
In this section, we perform a compensated algorithm (Comp-CTP) to evaluate Chebyshev tensor product series based on EFTs.The technology is to extend compensated Clenshaw algorithm to polynomials expressed in Chebyshev tensor product.In order to extend Clenshaw algorithm to tensor product case, we express the series as Therefore, we write Clenshaw tensor product algorithm (CTP) to evaluate the Chebyshev tensor product series with a nested Clenshaw algorithm (Algorithm 1).
Based on previous analysis, we apparently know that ê = ê2 ⊕ ê3 is the approximation of .So the result of CompCTP algorithm (, ) = P(, ) ⊕ ê is more accurate than the floating-point result P(, ) of Algorithm 1.

Error Analysis of CompCTP Algorithm
In this section, we carry out an error bound of ComCTP algorithm for Chebyshev tensor product surface.Firstly, we consider the error bound of the CTP algorithm.According to Theorem 1, we deduct the Theorem 3.  (
Finally, we show the forward error bound of CompCTP algorithm using previous analysis.

Experimental Results
In this section, we perform a series of numerical experiments.The programs are written in Matlab code and run with the software of MATLAB R2015b.We consider the polynomials with floating-point numbers coefficients and floatingpoint element (, ) expressed in Chebyshev tensor product.Considering the efficiency, we use CTP with quad-double arithmetic (QDCTP) [22] to accurately evaluate polynomial instead of symbolic toolbox.
Generally, we can get an accurate result using Algorithm 1 as the Chebyshev tensor product series is well-conditioned.But, we need a more accurate algorithm when the problem is ill-conditioned.We consider a bivariate polynomial in area [0, 1] × [0, 1] proposed by [14]  (, ) = ( − 0.75) 3 ( − 0.2) 3 ( − 0.75) and convert it to a Chebyshev tensor product form.We extend the univariate conversion algorithm [23] to change a bivariate polynomial in power form to a Chebyshev tensor product form.We transform the coefficients via the Matlab symbolic toolbox.Since the coefficients of polynomials in Chebyshev tensor product which are evaluated by us are floating-point numbers, they need to be rounded to the nearest floating-point elements.
We use CTP, CompCTP, and QDCTP (CTP algorithm along with quad-double arithmetic) [22] to compute the value of Chebyshev tensor product polynomial (, ) at 400 grid points near the multiple root (0.75, 0.2). Figure 1 shows the surface generated by different algorithms.It is obvious that the compensated algorithm can approximate the expected smooth surface as accurate as that using CTP algorithm with quad-double arithmetic, when the results are rounded to the working precision.Observe that the surface of CTP is a folding interface and varies slightly in the direction .The reason is that we firstly compute α() = Clenshaw((, : ), ) at the th step in the loop of the CTP algorithm and then compute α(0) = Clenshaw(α, ) leading to the more obviously influences of the round-off errors along .Thus, the CompCTP algorithm can compute a desired result.
Figure 2 performs an absolute forward error using the algorithms CTP and CompCTP for 400 points.It is clear that the CompCTP reduces the error better than CTP algorithm.Besides, we observe that the relative errors of CompCTP algorithm are smaller than the working precision  even near the point (0.75, 0.2).Therefore, the experiments verify our estimation of relation (37).
Next we consider the relative error bounds for illconditioned polynomials.We produce a series of illconditioned polynomials in Chebyshev tensor product and evaluate them using CTP, CompCTP, and DDCTP (CTP with double-double precision in Algorithms 3-7).We choose degree  ×  (where ,  is parameter from (6)) with condition numbers changing from 10 3 to 10 36 .The algorithm generating the ill-conditioned polynomials is shown in Algorithms 3-7 (GenPolyCTP), which is similar to algorithm GenPoly in [24].Considering the size of the problem and computational efficiency, we choose  ×  = 6 × 7. We plot the results in Figure 3. Obviously, we observe that the CompCTP illustrates an expected result, where the CompCTP is more stable and accurate than CTP.The relative errors are equal to or smaller than  when cond(, , ) ≤ 1/.While cond(, , ) ∈ [1/, 1/ 2 ], the relative errors increase almost linearly.Meanwhile, we notice that the CompCTP shows high accuracy because the rounding errors are recorded to approximate the real errors even under ill-condition.Besides, we also compute the ill-conditioned polynomials using DDCTP based on the Bailey's quad-double [2,25] arithmetic.Comparing with DDCTP, we can find that the results of CompCTP algorithm are as accurate as those computed using double-double arithmetic, which are illustrated by our numerical experiment.
Finally, we focus on the computational complexity of all the algorithms.

Algorithm 2 :
Compensated Clenshaw algorithm for evaluation of Chebyshev tensor product surface.
() acts as the theoretical error of ∑  =0     () at th step; combining Theorem 1 we get =0     ()  () be a polynomial in Chebyshev tensor product.If the condition number for polynomial evaluation of the (, ) at entry (, ) is defined by