# Accurate and Efficient Evaluation of Chebyshev Tensor Product Surface

Authors: Keshan He; Peibing Du; Hao Jiang; Chongwen Duan; Hongxia Wang; Lizhi Cheng
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

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 isu+Ou2×condP,x,y in contrast to classic scheme u×cond(P,x,y), where u is working precision and condP,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.


---

## Body

## 1. Introduction

Chebyshev polynomials have been extended to almost all mathematical and physical discipline, including spectral methods, approximation theory, and representation of potentials [1–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–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] (double-double 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–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 compensated 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 Section2, 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.


## 2. 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(1)flxopy=xopy1+ε1=xopy1+ε2,ε1,ε2≤u,where op∈+,-,×,/ and u is the working precision. For brevity we denote flxopy=xoy, o∈⊕,⊖,⊗,⊘. Besides, we denote γn≔nu/(1-nu)=nu+O(u2) [18] and use a^ and fl(a) as the computed element of a.Finally, we recall the Chebyshev polynomials and analogous Chebyshev polynomials [15, 19]. The forms of Chebyshev polynomials and analogous Chebyshev polynomials definite in the interval [-1,1] with three term recurrence are shown in (2) and (3), respectively.(2)T0x=1,T1x=x,⋮Tn+1x=2xTnx-Tn-1x,(3)T~0x=1,T~1x=x,⋮T~n+1x=2xT~nx+T~n-1x.


### 2.2. Error-Free Transformations

Rounding errors are an unavoidable consequence of working in finite precision arithmetic [18]. Error-free transformations (EFTs) are a technology of the floating-point operation +,-,×, which transforms any pair of floating-point numbers (a,b) into a new pair (x,y) with x=flaopb and aopb=x+y to obtain an accurate result. Two algorithms of EFTs are Donald et al.’s  TwoSum [20] (compensated summation of two floating-point numbers) and Dekker’s TwoProd algorithm [21] (compensated product of two floating-point numbers).The Clenshaw algorithm [16] is a recursive method to compute a linear combination of Chebyshev series p(x)=∑j=0ncjTj(x). Reviewing work [15], the forward error bound of Clenshaw algorithm satisfies (4).Theorem 1 (see [15]).
Letp(x)=∑j=0ncjTj(x) be a polynomial at point x and Clenshaw (p,x) denote the numerical result of Clenshaw algorithm; then(4)Clenshawp,x-px≤γ3n-1∑j=0ncjT~jx.Combining EFTs with Clenshaw algorithm, [15] proposes a compensated Clenshaw algorithm (CompClenshaw in Algorithms 3–7) to evaluate univariate finite Chebyshev series. The algorithm shows a smaller forward error bound than Clenshaw algorithm.Theorem 2 (see [15]).
Letp(x)=∑j=0ncjTj(x) be a finite Chebyshev series. The forward error bound of compensated Clenshaw algorithm (CompClenshaw)verifies(5)CompClenshawp,x-px≤upx+γ3n-12∑j=0ncjT~jx.Since the elementγ3n-12 is Ou2, comparing to (4), CompClenshaw is more stable to get an accurate result.


## 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(1)flxopy=xopy1+ε1=xopy1+ε2,ε1,ε2≤u,where op∈+,-,×,/ and u is the working precision. For brevity we denote flxopy=xoy, o∈⊕,⊖,⊗,⊘. Besides, we denote γn≔nu/(1-nu)=nu+O(u2) [18] and use a^ and fl(a) as the computed element of a.Finally, we recall the Chebyshev polynomials and analogous Chebyshev polynomials [15, 19]. The forms of Chebyshev polynomials and analogous Chebyshev polynomials definite in the interval [-1,1] with three term recurrence are shown in (2) and (3), respectively.(2)T0x=1,T1x=x,⋮Tn+1x=2xTnx-Tn-1x,(3)T~0x=1,T~1x=x,⋮T~n+1x=2xT~nx+T~n-1x.


## 2.2. Error-Free Transformations

Rounding errors are an unavoidable consequence of working in finite precision arithmetic [18]. Error-free transformations (EFTs) are a technology of the floating-point operation +,-,×, which transforms any pair of floating-point numbers (a,b) into a new pair (x,y) with x=flaopb and aopb=x+y to obtain an accurate result. Two algorithms of EFTs are Donald et al.’s  TwoSum [20] (compensated summation of two floating-point numbers) and Dekker’s TwoProd algorithm [21] (compensated product of two floating-point numbers).The Clenshaw algorithm [16] is a recursive method to compute a linear combination of Chebyshev series p(x)=∑j=0ncjTj(x). Reviewing work [15], the forward error bound of Clenshaw algorithm satisfies (4).Theorem 1 (see [15]).
Letp(x)=∑j=0ncjTj(x) be a polynomial at point x and Clenshaw (p,x) denote the numerical result of Clenshaw algorithm; then(4)Clenshawp,x-px≤γ3n-1∑j=0ncjT~jx.Combining EFTs with Clenshaw algorithm, [15] proposes a compensated Clenshaw algorithm (CompClenshaw in Algorithms 3–7) to evaluate univariate finite Chebyshev series. The algorithm shows a smaller forward error bound than Clenshaw algorithm.Theorem 2 (see [15]).
Letp(x)=∑j=0ncjTj(x) be a finite Chebyshev series. The forward error bound of compensated Clenshaw algorithm (CompClenshaw)verifies(5)CompClenshawp,x-px≤upx+γ3n-12∑j=0ncjT~jx.Since the elementγ3n-12 is Ou2, comparing to (4), CompClenshaw is more stable to get an accurate result.


## 3. Accurate Algorithm to Evaluate Chebyshev Tensor Product Surface

In this section, we perform a compensated algorithm(CompCTP) 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(6)Px,y=∑i=0m∑j=0naijTjyTix≔∑i=0mαi,yTix.Therefore, we write Clenshaw tensor product algorithm (CTP) to evaluate the Chebyshev tensor product series with a nested Clenshaw algorithm (Algorithm 1).Algorithm 1:Clenshaw algorithm for evaluation of Chebyshev tensor product surface.
functionres=CTPP,x,y% assuming P is the coefficients matrix of polynomial P(x,y)
α
=
z
e
r
o
s
(
1
,
m
)
for  i=0:m  do
α
^
(
i
)
=
C
l
e
n
s
h
a
w
(
P
(
i
,
:
)
,
y
)
end for
α
^
(
0
)
=
C
l
e
n
s
h
a
w
(
α
^
,
x
)
C
T
P
(
P
,
x
,
y
)
=
P
^
(
x
,
y
)
=
r
e
s
=
α
^
(
0
)
%  P^(x,y) is numerical result of P(x,y) using Clenshaw algorithmSubstitutingCompClenshaw  for  Clenshaw, we put forward a compensated Clenshaw tensor product algorithm to improve CTP algorithm. We call the compensated CTP algorithm as CompCTP (Algorithm 2).Algorithm 2:Compensated Clenshaw algorithm for evaluation of Chebyshev tensor product surface.
function[rescomp,resorg,err]=CompCTP(P,x,y)
α
=
z
e
r
o
s
(
1
,
m
)
for  i=0:m  do
[
α
^
(
i
)
,
e
^
1
(
i
)
]
=
C
o
m
p
C
l
e
n
s
h
a
w
(
P
(
i
,
:
)
,
y
)
end for
[
α
^
(
0
)
,
e
^
2
]
=
C
o
m
p
C
l
e
n
s
h
a
w
(
α
^
,
x
)
e
^
3
=
C
l
e
n
s
h
a
w
(
e
^
1
,
x
)
e
r
r
=
e
^
2
⊕
e
^
3
r
e
s
o
r
g
=
P
^
(
x
,
y
)
=
α
^
(
0
)
CompCTP(P,x,y)=rescomp=resorg⊕errAlgorithm 3:Generate polynomial P(x,y)=∑i=0m∑j=0naijTi(x)Ti(y) in the form of Chebyshev tensor product.
[
p
,
v
a
c
t
,
c
a
c
t
]
=
G
e
n
P
o
l
y
C
T
P
(
m
,
n
,
x
,
y
,
v
,
c
e
x
p
)
%v—the expected evaluation of the polynomial;
%cexp—the expected condition number;
%x,y—the coordinates;
%vact—the actual evaluation of the polynomial;
%cact—the actual condition number;
a
i
,
j
=
0, for i=0:m, j=0:n;
n
u
m
=
(
m
+
1
)
∗
(
n
+
1
);
d
2
=
c
e
i
l
(
n
u
m
/
2
);
v
b
=
l
o
g
2
(
c
e
x
p
∗
a
b
s
(
v
)
);
I
=
J
=
z
e
r
o
s
(
1
,
n
u
m
);
p
e
r
m
=
r
a
n
d
p
e
r
m
(
n
u
m
);
Changeperm to coordinates (I,J),
where(I(i),J(i)) is a random coefficient of the surface in accord with perm(i) in some sort order;
a
I
(
1
)
,
J
(
1
)
=
(
2
∗
r
a
n
d
-
1
)
/
T
I
(
1
)
(
x
)
T
J
(
1
)
(
y
);
a
I
(
2
)
,
J
(
2
)
=
(
2
∗
r
a
n
d
-
1
)
∗
2
v
b
/
T
I
(
2
)
(
x
)
T
J
(
2
)
(
y
);
for  i=3:d2  do
a
I
(
i
)
,
J
(
i
)
=
(
2
∗
r
a
n
d
-
1
)
∗
2
v
b
∗
r
a
n
d
/
T
I
(
i
)
(
x
)
T
J
(
i
)
(
y
);
end for
l
o
g
2
v
=
l
o
g
2
(
a
b
s
(
v
)
);
m
=
[
(
2
∗
r
a
n
d
(
1
,
n
u
m
-
d
2
)
-
1
)
∗
2
l
i
n
s
p
a
c
e
(
v
b
,
l
o
g
2
v
,
n
u
m
-
d
2
)
];
for  i=d2+1:num-1  do
a
I
(
i
)
,
J
(
i
)
=
(
m
(
i
-
d
2
)
-
Q
D
C
T
P
s
(
F
,
x
,
y
)
)
/
T
I
(
i
)
(
x
)
T
J
(
i
)
(
y
);
end for
a
I
(
n
u
m
)
,
J
(
n
u
m
)
=
(
v
-
Q
D
C
T
P
(
F
,
x
,
y
)
)
/
T
I
(
n
u
m
)
(
x
)
T
J
(
n
u
m
)
(
y
);
v
a
c
t
=
Q
D
C
T
P
(
P
,
x
,
y
);
c
a
c
t
=
Q
D
C
T
P
(
a
b
s
(
P
)
,
x
,
y
)
/
a
b
s
(
v
a
c
t
);Algorithm 4:[16] Clenshaw algorithm to evaluate finite Chebyshev series.
functionres=Clenshaw(p,x)
bn+2=bn+1=0
for  j=n:-1:1  do
b
j
=
2
x
b
j
+
1
-
b
j
+
2
+
c
j
end for
b
0
=
x
b
1
-
b
2
+
c
0
Clenshaw(p,x)=res=b0Algorithm 5:[15] compensated Clenshaw algorithm to evaluate Chebyshev series accurately.
functionres=CompClenshaw(p,x)
b^n+2=b^n+1=0
ε
b
n
+
2
=
ε
b
n
+
1
=
0
for  j=n:-1:1  do
[
s
,
π
j
]
=
T
w
o
P
r
o
d
(
b
^
j
+
1
,
2
x
)
[
v
,
σ
j
]
=
T
w
o
S
u
m
(
s
,
-
b
^
j
+
2
)
[
b
^
j
,
β
j
]
=
T
w
o
S
u
m
(
v
,
c
j
)
ω
^
j
=
π
j
⊕
σ
j
⊕
β
j
ε
b
^
j
=
2
x
⊗
ε
b
^
j
+
1
⊖
ε
b
^
j
+
2
⊕
ω
^
j
end for
[
s
,
π
0
]
=
T
w
o
P
r
o
d
(
b
^
1
,
x
)
[
v
,
σ
0
]
=
T
w
o
S
u
m
(
s
,
-
b
^
2
)
[
b
^
0
,
β
0
]
=
T
w
o
S
u
m
(
v
,
c
0
)
ω
^
0
=
π
0
⊕
σ
0
⊕
β
0
ε
b
^
0
=
x
⊗
ε
b
^
1
⊖
ε
b
^
2
⊕
ω
^
0
CompClenshaw(p,x)=res=b^0⊕εb0Algorithm 6:[15] Clenshaw algorithm in double-double format.
functionDDClenshaw(ph,pl,x)
b=zeros(2,N+2)
for  j=N:-1:2  do
[
r
h
1
,
r
l
1
]
=
p
r
o
d
_
d
d
_
d
(
b
(
1
,
j
+
1
)
,
b
(
2
,
j
+
1
)
,
2
∗
x
)
[
r
h
2
,
r
l
2
]
=
a
d
d
_
d
d
_
d
d
(
r
h
1
,
r
l
1
,
-
b
(
1
,
j
+
2
)
,
-
b
(
2
,
j
+
2
)
)
b
1
,
j
,
b
2
,
j
=
a
d
d
_
d
d
_
d
d
(
r
h
2
,
r
l
2
,
p
h
(
j
)
,
p
l
(
j
)
)
end for
[
r
h
1
,
r
l
1
]
=
p
r
o
d
_
d
d
_
d
(
b
(
1,2
)
,
b
(
2,2
)
,
2
∗
x
)
[
r
h
2
,
r
l
2
]
=
a
d
d
_
d
d
_
d
d
(
r
h
1
,
r
l
1
,
-
b
(
1,3
)
,
-
b
(
2,3
)
)
[
b
(
1
,
j
)
,
b
(
2
,
j
)
]
=
a
d
d
_
d
d
_
d
d
(
r
h
2
,
r
l
2
,
p
h
(
1
)
,
p
l
(
1
)
)
[rh,rl]=[b(1,j),b(2,j)]Algorithm 7:Clenshaw Chebyshev tensor product in double-double format.
function[rh,rl]=DDCTP(P,x,y)
for  i=1:m  do
[
f
1
(
i
)
,
f
2
(
i
)
]
=
D
D
C
l
e
n
s
h
a
w
(
A
(
i
,
:
)
,
a
0
,
y
) % here a0=zeros(n+1,1)
end for
[rh,rl]=DDClenshaw(f1,f2,x)According to Algorithm2, combining (7)∑i=0mα^iTix=P^x,y+e2,αi=α^i+e1i,0≤i≤m,e3=∑i=0me1iTix,0≤i≤m, we have(8)P^x,y+e2=∑i=0mαi-e1iTix=∑i=0m∑j=0naijTixTjx-∑i=0me1iTix=∑i=0m∑j=0naijTixTjx-e3; that is(9)Px,y=P^x,y+e,where e=e2+e3, e1(i) is the theoretical error produced by α^(i)=Clenshaw(P(i,:),y) at ith step, and e2 is the theoretical error generated by α^(0)=Clenshaw(α^(i),x).Based on previous analysis, we apparently know thate^=e^2⊕e^3 is the approximation of e. So the result of CompCTP algorithm P¯(x,y)=P^(x,y)⊕e^ is more accurate than the floating-point result P^(x,y) of Algorithm 1.


## 4. Error Analysis ofCompCTP Algorithm

In this section, we carry out an error bound ofComCTP 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.Theorem 3.
Let us consider a Chebyshev tensor product surfaceP(x,y)=∑i=0m∑j=0naijTi(x)Tj(y) and suppose that 3(m+n)u<1, where u is the working precision. Then the value P^(x,y) computed in floating-point arithmetic through Algorithm 1 satisfies(10)P^x,y-Px,y≤γ3m+n-2∑i=0m∑j=0naijT~ixT~jy.Proof.
SupposeP^y(x)≔∑i=0mp^i(y)Ti(x), where p^i(y) is the numerical result of Clenshaw algorithm by ∑j=0naijTj(y) at ith step. Using Theorem 1, we obtain (11)Px,y-P^x,y≤Px,y-P^yx+P^yx-P^x,y≤∑i=0mTix∑j=0naijTjy-p^iy+∑i=0mp^iyTix-P^x,y≤γ3n-1∑i=0m∑j=0naijT~ixT~jy+γ3m-1∑i=0mp^iyT~ix≤γ3n-1∑i=0m∑j=0naijT~ixT~jy+γ3m-11+γ3n-1∑i=0m∑j=0naijT~ixT~jy≤γ3m+n-2∑i=0m∑j=0naijT~ixT~jy.In order to analyze the error bound ofCompCTP algorithm, we need a lemma (Lemma 5). Firstly, we review a lemma in [15].Lemma 4 (see [15]).
Givenp(x)=∑j=0ncjTj(x) and πj,σj,βj,j=n-1:-1:0, is the round-off error of EFTs in CompClenshaw (Algorithm 5), one obtains(12)∑j=0n-1πj+σj+βjT~jx≤3nu1+γ3n∑j=0ncjT~jx.Lemma 5.
Letεb0 and εb^0 be the error of theoretical and numerical error of CompClenshaw, respectively. One can get(13)εb0-εb^0≤γ3n-1γ3n+1∑j=0ncjT~jx.Proof.
Obviously,εb^0 is the Clenshaw algorithm with coefficient ω^j in CompClenshaw (Algorithm 5). Using (4) we obtain(14)εb0-εb^0≤γ3n-1∑j=0nω^jT~jx=γ3n-1∑j=0nπj⊕σj⊕βjT~jx≤γ3n-1∑j=0nπj+σj+βj1+γ2T~jx≤γ3n-11+γ2∑j=0nπj+σj+βjT~jx.Using Lemma 4 to (14), we have (15)εb0-εb^0≤3nu1+γ3n-1γ3n-11+γ2∑j=0ncjT~jx. We obtain 3nuγ3n-1(1+γ3n-1)(1+γ2)≤γ3n-13nu(1+γ3n+1)≤γ3n-1γ3n+1. Then we deduct relation (13).Finally, we show the forward error bound ofCompCTP algorithm using previous analysis.Theorem 6.
LetP(x,y)=∑i=0m∑j=0naijTi(x)Tj(y) be a Chebyshev tensor product surface. The CompCTP algorithm satisfies(16)CompCTPP,x,y-Px,y≤uPx,y+3γ3m+12+γ3n+12∑i=0m∑j=0naijT~ixT~jy.Proof.
According to Algorithm2, we get the numerical result of CompCTP algorithm as (17)P¯x,y=P^x,y⊕e^=P^x,y+e^1+δ,δ≤u. Using relation (9) and (18)Px,y=P^x,y+e,e=e2+e3,e3=∑i=0me1iTix,we obtain(19)P¯x,y=Px,y1+δ+e^-e1+δ.Next, we consider the bound of e^-e. For e^=e^2⊕e^3=(e^2+e^3)(1+δ),δ<u, and e=e2+e3=(e2+e3)(1+δ)-δ(e2+e3), we get(20)e^-e≤ue2+e3+1+ue2-e^2+e3-e^3.Then, according to e=e2+e3=Px,y-P^x,y, we have(21)e2+e3≤γ3m+n∑i=0m∑j=0naijT~ixT~jy.Let us consider e2-e^2. Because e^2 is the error from CompClenshaw algorithm, according to Lemma 5, we obtain(22)e2-e^2≤γ3m-1γ3m+1∑i=0mα^iT~ix.Using the relation(23)α^i≤1+γ3n-1∑j=0naijT~jy,0≤i≤m and (22), we have(24)e2-e^2≤γ3m-1γ3m+11+γ3n-1∑i=0m∑j=0naijT~ixT~jy.Finally, we focus on the bound of e3-e^3. We assume e3mid=∑i=0me^1(i)Ti(x), using triangle inequality(25)e3-e^3≤e3-e3mid+e3mid-e^3;then apparently we have(26)e3-e3mid=∑i=0me1i-e^1iT~ix≤γ3n-1γ3n+1∑i=0m∑j=0naijT~ixT~jy.Actually(27)e3mid-e^3≤γ3m-1∑i=0me^1iT~ix;then(28)e^1i≤e1i+e1i-e^1i,0≤i≤m,where e1(i) acts as the theoretical error of ∑j=0naijTj(y) at ith step; combining Theorem 1 we get(29)e1i≤γ3n-1∑j=0naijT~jy.Combining (13), (28), and (29) we obtain(30)e^1i≤γ3n-1∑j=0naijT~jy+γ3n-1γ3n+1∑j=0naijT~jy.Synthetically, combining (27) and (30) we derive(31)e3mid-e^3≤γ3m-1∑i=0me^1iT~ix≤γ3m-1∑i=0me1i+e1i-e^1iT~ix≤γ3m-1γ3n-1+γ3n-1γ3n+1∑i=0m∑j=0naijT~ixT~jy.According to (25), (26), and (31), we obtain(32)e3-e^3≤γ3n-1γ3n+1+γ3m-1γ3n-1+γ3n-1γ3n+1∑i=0m∑j=0naijT~ixT~jy,using (20), (21), (24), and (32), the error bound yields(33)e-e^≤αm,n∑i=0m∑j=0naijT~ixT~jy,where α(m,n) is (34)αm,n=uγ3m+n+1+uγ3m-1γ3m+11+γ3n-1+γ3n-1γ3n+11+γ3m-1+γ3m-1γ3n-1. By the properties of the element γn [18], we deduct some inequations as follows: (35)1+uuγ3m+n≤12γ3m+12+γ3n+12;1+u21+γ3m-1γ3n-1γ3n+1≤21+u2γ3n-1γ3n+1≤2γ3n+12;1+u2γ3n-1γ3m-1≤γ3mγ3n≤12γ3m+12+γ3n+12;1+u21+γ3n-1γ3m+1γ3m-1≤2γ3m+12.  So we have (1+u)α(m,n)≤3(γ3m+12+γ3n+12) and combine (19) with (33); then we obtain relation (16).LetP(x,y)=∑i=0m∑j=0naijTi(x)Tj(y) be a polynomial in Chebyshev tensor product. If the condition number for polynomial evaluation of the P(x,y) at entry (x,y) is defined by(36)condP,x,y=P~x,yPx,y=∑i=0m∑j=0naijT~ixT~jy∑i=0m∑j=0naijTixTjy, we show the relative error bound of the CTP and CompCTP algorithm in Corollary 7.Corollary 7.
LetP(x,y)=∑i=0m∑j=0naijTixTj(y) be a polynomial in Chebyshev tensor product. The relative error bound of the CTP algorithm and CompCTP algorithm yields(37)CTPP,x,y-Px,yPx,y≤γ3m+n-2condP,x,y,CompCTPP,x,y-Px,yPx,y≤u+3γ3m+12+γ3n+12condP,x,y.As Corollary7 illustrates that if 3γ3m+12+γ3n+12condP,x,y<u, the error of the result evaluated by CompCTP is bounded by working precision u. Besides that, the computed result generated by the CompCTP scheme is as accurate as if evaluated in double-double precision.


## 5. 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 floating-point elementx,y 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 Algorithm1 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] (38)Px,y=x-0.753x-0.23y-0.753y-0.23and 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 useCTP, CompCTP, and QDCTP (CTP algorithm along with quad-double arithmetic) [22] to compute the value of Chebyshev tensor product polynomial P(x,y) 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 x. The reason is that we firstly compute α^(i)=Clenshaw(P(i,:),y) at the ith step in the loop of the CTP algorithm and then compute α^(0)=Clenshaw(α^,x) leading to the more obviously influences of the round-off errors along x. Thus, the CompCTP algorithm can compute a desired result.Figure 1
The surface usingCTP, QDCTP, and CompCTP algorithm to evaluate Chebyshev tensor product series.Figure2 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 u even near the point (0.75,0.2). Therefore, the experiments verify our estimation of relation (37).Figure 2
Accuracy of the absolute error ofCTP algorithm (a) and CompCTP algorithm (b).
(a)
(b)Next we consider the relative error bounds for ill-conditioned polynomials. We produce a series of ill-conditioned polynomials in Chebyshev tensor product and evaluate them usingCTP, CompCTP, and DDCTP (CTP with double-double precision in Algorithms 3–7). We choose degree m×n (where m,n is parameter from (6)) with condition numbers changing from 103 to 1036. 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 m×n=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 u when cond(P,x,y)≤1/u. While cond(P,x,y)∈1/u,1/u2, 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’squad-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.Figure 3
The relative forward error by usingCTP, CompCTP, and DDCTP to evaluate ill-conditioned polynomials.Finally, we focus on the computational complexity of all the algorithms.Clenshaw: 3n+4CompClenshaw: 36n+38DDClenshaw: 52n+53CTP: (3n+4)(m+1)+3m+4CompCTP: (36n+38)(m+1)+39m+43DDCTP: (52n+53)(m+1)+(52m+53)Considering the previous comparisons of the accuracy, we can confirm thatCompCTP algorithm is as accurate as computation with DDCTP algorithm (shown in Algorithms 3–7). However, CompCTP only requires on the average about 69.2% of flops. So our algorithm is more efficient than CTP with double-double arithmetic.


## 6. Conclusion

We present an accurate and efficient algorithm for evaluation of Chebyshev tensor product surface, which is based on the Clenshaw algorithm and error-free transformations. The error analysis shows thatCompCTP algorithm can get the same accuracy as that computed by the traditional CTP algorithm with twice working precision. Besides, this compensated algorithm can run more efficiently than CTP algorithm with double-double precision. Experiments illustrate that our algorithm is stable and accurate even in some ill-condition cases.



---

*Source: 1729287-2017-09-27.xml*