- Email: [email protected]
Computing the polynomial remainder sequence via Bézout matrices
Journal of Computational and Applied Mathematics 250 (2013) 244–255
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Computing the polynomial remainder sequence via Bézout matrices Skander Belhaj ∗ Manouba University, ISAMM, Campus Universitaire de la Manouba, 2010 Tunis, Tunisia University of Tunis El Manar, ENIT-LAMSIN, BP 37, 1002 Tunis, Tunisia
article
abstract
info
Article history: Received 22 December 2011 Received in revised form 15 February 2013
In this paper, we address the task of computing the polynomial remainder sequence appearing in the Euclidean algorithm applied to two polynomials u(x) and v(x) of degree n and m, respectively, m < n via renewed approach based on the approximate block diagonalization of a Bézout matrix B (u, v). More specifically, this algorithm consists of using the efficient Schur complementation method given in terms of a Bézout matrix B(u, v) associated with the input polynomials. We also compare the proposed approach to the approximate block diagonalization of a Hankel matrix H (u, v) in which the successive transformation matrices are upper triangular Toeplitz matrices. All algorithms have been implemented in Matlab. Numerical experiments performed with a wide variety of test problems show the effectiveness of this algorithm in terms of efficiency, stability and robustness. © 2013 Elsevier B.V. All rights reserved.
MSC: 65Fxx Keywords: Hankel matrix Bézout matrix Block diagonalization Euclidean algorithm Polynomial remainder sequence Schur complementation
1. Introduction Let u(x) =
n k=0
uk x k
and
v(x) =
m
v k xk
(1)
k=0
be two polynomials with deg (u(x)) = n, deg (v(x)) = m and m < n. The classical Euclidean algorithm applied to u(x) and v(x) returns a sequence of quotients qk (x) and remainders rk (x), such that r−1 (x) = u(x), r0 (x) = v(x) rk−2 (x) = rk−1 (x) qk (x) − rk (x) ,
k = 1, . . . , K
(2)
where −rk (x) is the polynomial remainder of the division of rk−2 (x) and rk−1 (x) and rK (x) is the greatest common divisor (GCD) of u(x) and v(x). The question of computing the quotients and the remainders is beginning to be studied by Collins and Brown (see [1–3]) via efficient Euclidean-based methods. Naturally, the classical algorithm turns out to be unstable when applied to approximate polynomials. One might think of designing stabilized numerical versions of the Euclidean algorithm to the approximate case (see [4,5]). Several papers have been devoted to describe the natural relation between the Euclidean algorithm and the block LU factorization of the Hankel and Bézout matrices. Sure enough, Bini and Gemignani (see [6–9])
∗
Correspondence to: Manouba University, ISAMM, Campus Universitaire de la Manouba, 2010 Tunis, Tunisia. E-mail address: [email protected].
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.02.029
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
245
compute the coefficients of the polynomials in (2) via renewed approach based on a block LU factorization (see [10,11]) of the Hankel matrix H (u, v), associated with u(x) and v(x). More specifically, At H (u, v) A = D
(3)
where A is an upper triangular matrix with entries equal to one and D is a block diagonal matrix and each block is a lower Hankel triangular matrix. Lately, Ben Atti and Diaz-Toca [12] calculate the coefficients of the polynomials in (2) via a different block LU factorization of the Hankel matrices [13] from the classical one. Naturally, when the input polynomials (1) are introduced via a perturbation (add k × 10−13 to all entries of ‘‘exact’’ input polynomials with k is taken randomly in (−1, 1)), the beautiful relation between the Euclidean algorithm and the block LU factorization of the Hankel and Bézout matrices may be lost. In addition, the polynomial remainder sequence (prs), is well known to be an ill-conditioned problem. In fact, small perturbations of the input polynomials may generate big variations in the sequence (even if computed exactly). Moreover, the floating point computation of the polynomial remainder sequence, performed by means of subsequent division, is a numerically unstable computation. In fact, just after the first step, the errors introduced by the floating point arithmetic provide a perturbation to the input of the subsequent steps which form again a prs which is an ill-conditioned problem. Thus these errors are largely amplified and spoil the final result. This is why we try to generalize this to the computation of approximate versions of these polynomials and factorizations. For more details on this research topic we refer the reader to [14–22] and to the references therein. It is shown that an n × n Hankel matrix admits an approximate block diagonalization in which the successive transformation matrices are upper triangular Toeplitz matrices [23,24]. More explicitly, lHϵm1
Θ1,2
Θ2,1 Dϵ = Aϵ HAϵ = . ..
lHϵm2
t
Θn,1
..
. ···
··· .. . .. . Θn,n−1
Θ1,n .. . Θn−1,n
(4)
lHϵml
where every block lHϵmj , for m1 + m2 + · · · + ml = n, l = 1, 2, . . . is an approximate lower Hankel matrix, all entries of Θj,k ,
for j, k = 1, 2, . . . , n and j ̸= k are very close to zero and Aϵ is an approximate upper triangular matrix. This algorithm is revised in order to obtain the approximate quotients polynomials and the approximate remainders polynomials appearing in the Euclidean Algorithm applied to two polynomials u(x) and v(x) of degree n and m, respectively, m < n [25]. The main contribution of this paper is to ensure the natural relation between the Euclidean algorithm and the approximate block LU factorization diagonalization of Hankel and Bézout matrices associated with u(x) and v(x) and to propose an efficient Schur-based algorithm of the corresponding Bézout matrix B(u, v) (see the works by Bini and Gemignani [6–9]) which provide the approximate quotients and the approximate remainders during the process execution. The paper is organized as follows: Section 2 gives some notations and definitions. Some theoretical results associated with the approximate block diagonalization of Hankel matrices are recalled in Section 3. Section 4 presents an efficient algorithm for an approximate block diagonalization of n × n Bézout matrices applied to two polynomials u(x) and v(x) of degree n and m, respectively, m < n and relates this approach to the Euclidean Algorithm in order to obtain the approximate quotients and the approximate remainders appearing during the process execution. In Section 5 we illustrate our method with examples given by an implementation of the procedures using Matlab and we propose a statistical study performed with a wide variety of test problems to show the effectiveness of this algorithm in terms of speed, stability and robustness. Finally, a summary and future research are given in Section 6 to complete the paper. 2. Notations and definitions We begin the section introducing the following notations:
• H (S ) ∈ Rn×n denotes the Hankel matrix associated with a list S of length (2n − 1). This means that the first row is given by the first n terms of S and the last column is given by the last n terms of S.
• lH (S ) ∈ Rn×n denotes the Hankel triangular matrix (with respect to the antidiagonal) associated with a list S of length (2n − 1) such that the last column is defined by S. • uH (S ) ∈ Rn×n denotes the Hankel triangular matrix (with respect to the antidiagonal) associated with a list S of length (2n − 1) such that the first column is defined by S. • H (S ; m; n) ∈ Rm×n denotes the Hankel matrix associated with a list S of length (m + n − 1). The first row is given by the first n terms of S and the last column is given by the last m terms of S.
• T (S ; m; n) ∈ Rm×n denotes the Toeplitz matrix associated with a list S of length (m + n − 1). The first row is given by the first m terms of S and the last column is given by the last n terms of S.
• uT (S ) ∈ Rn×n denotes the upper Toeplitz triangular matrix associated with a list S such that the first row is defined by S.
• lT (S ) ∈ Rn×n denotes the lowerToeplitz triangular matrix associated with a list S such that the last row is defined by S. p • Let p ∈ N. Let Σp ∈ Rp×p , Σp = εjk j, k=1 , where all entries of Σp are zero except that εj+k,j = εk for j+k, k = 1, 2, . . . , p.
246
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
p are zero except that εj+k,j for j + k, k = 1, 2, . . . , p. p ∈ Rp×p , Σ = εjk p , where all entries of Σ • Let p ∈ N. Let Σ j, k=1 • Jp = lH (1, 0, . . . , 0), p ∈ N. p−1
• • • •
Given P ∈ Rn×m , P = Jm P t Jn . T (S ; m; n) = Jm H (S ; m; n). uT (S ) = Jl lH (S ), where l is the length of L. Let a ∈ R. Let µ > 0, V (a, µ) = (a − µ; a + µ) is a neighborhood of a.
2.1. Hankel matrices k k Let u(x) = k=0 uk x and v(x) = k=0 vk x be two polynomials of degree n and m, respectively, where m < n. The power series expansion R (x) of the function v(x)/u(x) at the infinity
n
R ( x) =
∞
m
hk x −k
(5)
k=0
defines the n × n Hankel matrix, H = H (u, v), associated with u(x) and v(x),
.
.. .
··· ··· .. .
hn
hn +1
···
h1 h2
H = ..
h2 h3
hn hn+1 h2n−1
.. . .
(6)
In addition, every nonsingular Hankel matrix can be viewed as a Hankel matrix associated with the given polynomials. Proposition 2.1. For any nonsingular n × n Hankel matrix H there exists two coprime polynomials u(x) and v(x) of degree n and m, respectively, where m < n, such that H = H (u, v). The polynomial u(x) and v(x) are related to H by the following equalities: H (u0 , . . . , un−1 )t = −un (hn+1 , . . . , h2n ) ,
(vn−1 , . . . , v0 )t = lT (h1 , . . . , hn ) (un , . . . , u1 )t where h2n is any number and lT (h1 , . . . , hn ) is the lower triangular Toeplitz matrix defined by the list [h1 , . . . , hn ]. Proof (See [7]).
Remark 2.1. Suppose that H has the following structure: H = lH (hn , . . . , h2n−1 ) . Thus, Proposition 2.1 concludes, for any polynomial v(x),
0
(vn−1 , . . . , v0 )t = ...
hn
..
. ···
t t (un , . . . , u1 ) = (0, . . . , 0, hn un ) . 0
If v(x) = 1 then v0 = 1 and un = 1/hn . Thus, H represents a Hankel matrix associated with u(x) and v(x) = 1 : H = H (u, 1). Suppose that we know all the entries of a Hankel matrix H and all the coefficients of u(x). Thus, the algorithm which provides all the coefficients of v(x) can be formulated as follows: Algorithm 2.1 (Construct v(x) from u(x) and a Hankel Matrix H). Given a Hankel matrix H and u(x) = m k polynomial of degree n, this algorithm provides all the coefficients of v(x) = k=0 vk x . 1. 2. 3. 4. 5.
Set H = [hn−m+2 , . . . , hn−1 ] and U = [u0 , . . . , un−1 , un ]. Set n = length (U ), U0 = u0 and U = [u1 , . . . , un−1 , un ]. Set n = length (U ). Compute U = Jn ∗ U t . Compute V = (Jn ∗ lT (H (1, 1 : n)) ∗ U )t .
n
k=0
uk x k a
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
247
2.2. Bézout matrices Let u(x) =
n
k =0
b (x, y) =
uk xk and v(x) =
m
k=0
u (x) v (y) − u (y) v (x) x−y
vk xk be two polynomials of degree n and m, respectively, where m < n. Let
.
The rational function b (x, y) is easily proved to be a polynomial in two variables and can therefore be written as b (x, y) = n−1 j k b x y . The matrix B x , y = b ( ) j,k 0≤j,k≤n−1 is called the Bézout matrix or Bezoutian associated with u(x) and v(x). j, k=0 j,k We next state more different representations of the Bézout matrix. Proposition 2.2. The Bezoutian matrix can be written as: B (u, v) = uH (u1 , . . . , un ) uT (v0 , . . . , vn−1 ) − uH (v1 , . . . , vn ) uT (u0 , . . . , un−1 ) where if m < n the coefficients vm+1 , . . . , vn are assumed to be zero. Proof (See [26]).
As a consequence of the previous result, we have: Corollary 2.1. The entries of B (u, v) can be recursively calculated as follows: bj,k+1 = bj+1,k + uj vk − uk vj setting bn+1,k = b0,k = 0 for all j, k and assuming vm+1 = · · · = vn = 0 as above. From Proposition 2.2, we observe that the matrix B (u, 1) is an upper Hankel triangular matrix B (u, 1) = uH (u1 , . . . , un ). Moreover, we see that B (u, 1)−1 = H (u, 1). Suppose now that we know all the entries of a Bézout matrix B(u, v) and all the coefficients of u(x). We finish the section proposing the algorithm which provides all the coefficients of v(x) as follows: Algorithm 2.2 (Construct v(x) from u(x) and a Bézout matrix B(u, v)). Given a Bézout matrix B(u, v) and u(x) = m k a polynomial of degree n, this algorithm provides all the coefficients of v(x) = k=0 vk x . 1. 2. 3. 4.
n
k =0
uk x k
Set U = [u0 , . . . , un ] and n = length (U ). Compute B(u, 1) = uH (u1 , . . . , un ). Compute H = H (u, 1) ∗ Jn B(u, v)Jn ∗ H (u, 1). Call Algorithm 2.1 to conclude V = [v0 , . . . , vm ].
3. Approximate block diagonalization of a Hankel matrix We introduce the approximate block diagonalization of a Hankel matrix presented in [23]. Lemma 3.1. Let n ∈ N ∗ . Let h = H (h1 , . . . , h2n−1 ) be a square Hankel matrix of order n. Suppose that hj = εj with εj ∈ V (0, µ) for j = 1, 2, . . . , p − 1 and hp ̸∈ V (0, µ). Then h has the form h = H ε1 , . . . , εp−1 , hp , . . . , h2n−1 =
h11 h21
h12 h22
,
(7)
where h11 = lH hp , . . . , h2p−1 + Jp Σp ,
h22 = H h2p+1 , . . . , h2n−1 ,
h12 = H hp+1 , . . . , hn+p−1 ; p; n − p ,
h21 = ht12 .
We can successively construct from h the following two matrices:
• A square lower Hankel triangular matrix H of order (2n − p), 0 0 H13 0 h11 h12 , H = lH hp , . . . , h2n−1 = H31
ht12
where H31 = H13 = lH (hp , . . . , hn−1 ). • A square upper triangular Toeplitz matrix T , t11 T = J2n−p H =uT hp , . . . , h2n−1 = 0 0
t12 t22 0
t13 t23 t33
,
where t11 = t33 = uT hp , . . . , hn−1 , t22 = Jp h11 , t13 = Jn−p h22 , t12 = Jn−p ht12 , and t23 = Jp h12 .
(8)
h22
(9)
248
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
Lemma 3.2. Let T be an upper triangular Toeplitz matrix, nonsingular with non-zero diagonal. Then T −1 = uT µ1 , . . . , µ2n−p and has the following block decomposition,
T −1 −1 T = 0 0
−1 T −1 12
11
T
−1 −1 t11 T −1 13 = 0 T −1 23
22
0
T
33
P
Q P , −1 t11
−1 t22 0
0
(10)
where P = T (µ2 , . . . , µn ; p; n − p) ,
P = Jn−p P t Jp ,
t22 P + t23 t11 = 0(p,n−p) ,
−1 h11 P + Mt11 = 0(p,n−p) ,
−1 t11 P + t12 t22 = 0(n−p,p) ,
−1 t11 Q + t12 P + t13 t11 = 0(n−p,n−p) .
−1
Theorem 3.1. Let h = H ε1 , . . . , εp−1 , hp , . . . , h2n−1 be an approximate Hankel matrix (approximation of an ‘‘exact’’ Hankel
matrix, i.e. H 0, . . . , 0, hp , . . . , h2n−1 ) where εj ∈ V (0, µ) for j = 1, 2, . . . , p − 1 with hp ̸∈ V (0, µ) and
t = uT hp , . . . , hn+p−1 ,
t −1 = uT (µ1 , . . . , µn ) .
Then h′ = t −1
t
ht −1 =
ϵ′
h′
11 t ϵ′
h′22
,
(11)
where −1 h′11 = lH µ1 , . . . , µp + t22
t
−1 Jp Σp t22 ,
h′22 = −H µp+2 , . . . , µ2n−p + P t Jp Σp P ,
−1 t ϵ ′ = t22 J p Σp P .
Remark 3.1. Theorem 3.1 gives an approximate reduction for real Hankel matrix. Thus if we have εj = 0, Theorem 3.1 provides the exact case. To iterate this result, h′22 must be a Hankel matrix, but −H µp+2 , . . . , µ2n−p is a Hankel matrix and P t Jp Σp P is only a symmetric matrix. If we choose εj very close to zero, we can conclude that Σp ≈ 0p , and the process converges to the exact case. To solve this problem, we can choose h′22 = −H µp+2 , . . . , µ2n−p + Θ where Θ is a Hankel matrix built with the first column and the last row of P t Jp Σp P. Remark 3.2. The key of our approach is based on the inversion of triangular Toeplitz matrix. Thus, we use superfast techniques described in [27]. Then, we give the fast algorithm to compute the block factorization as follows: Algorithm 3.1 (Approximate Block Diagonalization for Hankel Matrix). Given a list S = [ε1 , . . . , εp−1 , hp , . . . , h2n−1 ] where εj ∈ V (0, µ) for j = 1, 2, . . . , p − 1 with hp ̸∈ V (0, µ), this algorithm computes the approximate block diagonalization of Hankel matrix via upper triangular Toeplitz matrix. 1. Define a Hankel matrix h = H ε1 , . . . , εp−1 , hp , . . . , h2n−1 .
2. Define an upper triangular Toeplitz matrix t = uT hp , . . . , hn+p−1 .
3. 4. 5. 6.
Compute t −1 . t Compute h′ = t −1 ht −1 . Set h′11 = h′ (1 : p, 1 : p) and h′22 = h′ (p + 1 : n, p + 1 : n). Recursively apply Algorithm 3.1 to h = H h′22 (1 : n − p, 1) h′22 (n − p, 1 : n − p) ,
obtaining (4). 4. Block diagonalization of Hankel and Bézout matrices: connection to the Euclidean algorithm The topic of this section is the relationship between the Euclidean algorithm applied to a pair of polynomials u(x) and v(x) and the LU factorization of the associated Bézout or Hankel matrix. The main result essentially states that the entries of the approximate triangular matrices of the block triangular factorization of H (u, v) or Jn B(u, v)Jn yield the coefficients of all the polynomials generated by the Euclidean scheme.
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
249
4.1. The Euclidean algorithm Recall that the classical Euclidean algorithm for the computation of the GCD goes as follows. Let u(x) and v(x) be two polynomials in R [x] with deg (u(x)) = n, deg (v(x)) = m and m < n. The classical Euclidean algorithm applied to u(x) and v(x) returns a sequence of quotient qk (x) and remainder rk (x), such that r0 (x) = v(x) r−1 (x) = u(x), (12) rk−2 (x) = rk−1 (x) qk (x) − rk (x) , k = 1, . . . , K where −rk (x) is the remainder of the division of rk−2 (x) and rk−1 (x). The algorithm stops when a remainder rK +1 (x) = 0 is found; then rK (x) is the desired greatest common divisor (GCD) of u(x) and v(x). Observe that, since the degrees of the rk (x)’s are strictly decreasing, the algorithm must terminate in a finite number of steps. One might think of adapting the Euclidean algorithm to the approximate case by simply converting all ‘‘small’’ coefficients into 0 after each polynomial division (where ‘‘small’’ coefficients are coefficients whose absolute value is smaller than a fixed tolerance). Unfortunately, this naive method is highly unstable, as was already observed in [4]. It is possible, however, to design stabilized versions of the Euclidean algorithm. A first step in this direction is found in [5], where the authors introduce a normalization rule on remainders and calculate the approximate GCD as follows. Algorithm 4.1 (Computing the Polynomial Remainder Sequence and the Polynomial Quotients Via Euclid). Given u(x) and v(x) two polynomials in R [x] of degree n and m, respectively, where m < n, this algorithm computes the approximate quotients and the approximate remainders with a fixed tolerance less than a small positive number ϵ, 0 < ϵ ≪ 1. 1. Compute qk as the quotient of rk−1 divided by rk . 2. Set rk−1 = qk rk + max {1, ∥qk ∥∞ } · rk+1 . 3. Recursively apply Algorithm 4.1 for k = 0, 1, . . . , K and stops when a remainder |rK +1 | < ϵ . Remark 4.1. It must be pointed out that any variant of the approximate Euclidean algorithm (see for example [28,29]) can be used. We propose to implement Algorithm 3.1 to make a comparison with our method (see Section 5). 4.2. Approximate reduction of a Hankel matrix associated with two polynomials We recall the correlation between the approximate Euclidean algorithm applied to two polynomials u(x) and v(x) and the approximate block diagonalization of a Hankel matrix [25]. Corollary 4.1. Let H (u, v) = H (ε1 , . . . , εn−m−1 , hn−m , . . . , h2n−1 ) be an approximate Hankel as defined m matrix n k k in Theorem 3.1 associated with two coprime polynomials in R [x] , u(x) = k=0 uk x and v(x) = k=0 vk x , deg (u(x)) = n, deg (v(x)) = m and m < n. Let t = uT (hn−m , . . . , h2n−m−1 ) ,
t −1 = uT (µ1 , . . . , µn ) .
Then
t
−1 t
J B (q, 1) J −1 H (u, v) t = n−m ′ t n−m ϵ
ϵ′ , H (v, r )
(13)
where −1 J B (q, 1) J = Jn−m B (q, 1) Jn−m + Jn−m t22
H (v, r ) = H (v, r ) + P t Jn−m Σn−m P ,
t
−1 Jn−m Σn−m t22 Jn−m ,
−1 t ϵ ′ = t22 Jn−m Σn−m P ,
q(x) and r (x) are the polynomial quotient and remainder of the division u (x) /v (x). Let us devise the algorithm for the computation of the approximate polynomials quotients and the approximate polynomials remainders appearing in the Euclidean algorithm applied to two polynomials u(x) and v(x) of degree n and m, respectively, m < n based on the above corollary. Algorithm 4.2 (Computing Remainder Sequence and the Polynomial Quotients Via H (u, v)). Given u(x) = n m the Polynomial k k k=0 uk x and v(x) = k=0 vk x two polynomials of degree n and m, respectively, where m < n, this algorithm computes the approximate polynomials quotients and the approximate polynomials remainders appearing in the Euclidean Algorithm 1. 2. 3. 4. 5. 6. 7.
Construct a Hankel matrix H = H (0, . . . , 0, hn−m , . . . , h2n−1 ). Define an upper triangular matrix t = uT (hn−m , . . . , h2n−m−1 ). Toeplitz t Compute t −1 and h′ = t −1 ht −1 . Set h′11 = h′ (1 : n − m, 1 : n − m) and h′22 = h′ (n − m + 1 : n, n − m + 1 : n). Recover the coefficients of the quotient polynomial from step 3. Recover the coefficients of the remainder polynomial from Algorithm 2.1 by using steps 1. and 5. Recursively apply Algorithm 4.2 to h = H h′22 (1 : m, 1) h′22 (m, 1 : m) , obtaining all the coefficients appearing in the Euclidean Algorithm.
250
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
4.3. Approximate reduction of a Bézout matrix associated with two polynomials We explore the relationship between the approximate Euclidean algorithm applied to two polynomials u(x) and v(x) and the approximate block diagonalization of a Bézout matrix (see the works by Bini and Gemignani [6–9]). The following proposition explains how the polynomials generated by the approximate Euclidean algorithm can be obtained through Schur complements of the approximate permuted Bezoutian Jn B (u, v) Jn =
Pm
i
X
Xt Y
,
Pmi = Pmi + J Σ . where Pmi = (Jn B(u, v)Jn )mi and k k Proposition 4.1. Let u(x) = k=0 uk x and v(x) = k=0 vk x be two coprime polynomials with deg (u(x)) = n, deg (v(x)) = m and m < n. Let H be the associated Hankel matrix as in Theorem 3.1 and let 0 = m0 < m1 < · · · < ml = n be integers Pmi ̸∈ V (0, µ) for i = 1, . . . , l and such that det Hmi ̸∈ V (0, µ) for i = 0, . . . , l and det Hk ∈ V (0, µ) otherwise. Then det det Pmi ∈ V (0, µ) for j ̸= mi , i = 1, . . . , l. Moreover, if Smi is the Schur complement of Pmi and Smi is the Schur complement of Pmi in Jn B(u, v)Jn , the following relation holds:
n
m
−1 t Smi = Y − XPm X = Jn B(ri , ri+1 )Jn i
−1 −1 and Smi = Smi + X t Pm J Σ X t Pm i i
t
,
where {ri (x)}1,...,l is the polynomial remainder sequence obtained by applying the Euclidean algorithm to u(x) and v(x). Proof. We may inductively apply the same proof in [7].
Corollary 4.2. Under the hypothesis of Proposition 4.1. Then Jn B(u, v)Jn =
=
I −1 XPm i
0 I
T
0 I
−1 XPm T i
Θ Sm
Pmi Θt
−1 t Pm X i I
I 0
i
(i)
Jn B(u , v (i) )Jn + J Σ
Θ Sm
t Θ
i
Tt 0
−1 t T t Pm X i I
(14)
where
−1 −1 t Smi = Smi + X t Pm J Σ X t Pm , i i −1 t (i) −1 = T −1 Θ , Jn , Θ = −J Σ XPmi , Θ T = lT un , . . . , un−mi +1 B(u , 1) and the polynomials u(i) (x) , v (i) (x) verify H (u, v)mi = H u(i) , v (i) . −1 t Smi = Y − XPm X = Jn B(ri , ri+1 )Jn , i
Proof. The corollary follows from [7] and Proposition 4.1.
From Corollary 4.2, we present a parallel algorithm for an approximate block diagonalization of Jn B(u, v)Jn . Let us compare this method with Algorithm 4.2. We start computing the first diagonal block (a lower Hankel matrix). We denote by DBe´ zout the block diagonal matrix of Jn B(u, v)Jn . Let D11 be the first block. Then m1 = n − m, h11 = H(u, v)n−m = lH (un−m , . . . , u2n−2m−1 ) and so, by Remark 2.1, we can choose v (1) (x) = 1 such that H (u, v)n−m = H u(1) , 1 . Thus,
Pmi = Jn B(u(1) , v (1) )Jn + J Σ = Jn B(u(1) , 1)Jn + J Σ 1 = Jn H (u(1) , 1)−1 Jn + J Σ = Jn h− 11 Jn + J Σ = Jn B(q1 , 1)Jn + J Σ .
So, D11 = Jn B(q1 , 1)Jn + J Σ . Successively, we obtain an approximate upper matrix ABe´ zout and DBe´ zout such that AtBe´ zout Jn B(u, v)Jn ABe´ zout = DBe´ zout . Moreover the successive matrices B(ri , ri+1 ) directly provide the approximate remainders. Thus, we obtain the same blocks Dii for both Bézout and Hankel matrices, except the last one. Therefore, the coefficient vector of the (i + 1)-th approximate remainder ri+1 can be expressed as
−1 t −1 −1 ri+1 = Jn Y − XPm X Jn + Jn X t Pm J Σ X t Pm i i i
t
Jn en−mi ,
that is, ri+1 can be computed by solving an mi × mi system with Hankel-like matrix Pmi . Thus, we give an algorithm which computes the approximate polynomials quotients and the approximate polynomials remainders appearing in the Euclidean algorithm applied to two polynomials u(x) and v(x) of degree n and m, respectively, m < n based on the above corollary.
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
251
Algorithm 4.3 (Computing the Polynomial Remainder Sequence and the Polynomial Quotients via Jn B (u, v) Jn ). Given u(x) = n m k k k=0 uk x and v(x) = k=0 vk x two polynomials of degree n and m, respectively, where m < n, this algorithm computes the approximate polynomials quotients and the approximate polynomials remainders appearing in the Euclidean Algorithm. 1. Construct the Bézout matrix B (u, v). 2. Reduce Jn B (u, v) Jn as (14).
t
−1 t −1 −1 −1 X + X t Pm J Σ X t Pm . Pm and Sm1 = Y − XPm 3. Compute Pm1 = Jn B(u(1) , v (1) )Jn + J Σ , 1 1 1 1
4. Recover the coefficients of the first polynomial quotient from Algorithm 2.2 applied to Pm1 = Jn B(u(1) , 1)Jn + J Σ . 5. Recover the coefficients of the first polynomial remainder from Algorithm 2.2 applied to Sm1
−1 −1 X t Pm J Σ X t Pm 1 1
t
= Jn B(r1 , r2 )Jn +
with r1 = v (1) .
6. Recursively apply Algorithm 4.3 to Jn B(u, v)Jn = Smi , obtaining all the coefficients appearing in the Euclidean Algorithm. Remark 4.2. Observe that the assumption in Proposition 4.1 that u and v are coprime is no loss of generality. Indeed, if rL = GCD(u, v) were a polynomial of positive degree, then the Euclidean remainder sequence {ri+1 } can be recovered from the sequence { ri+1 } generated by the Euclidean algorithm applied to u(x)/rL (x) and v(x)/rL (x). 5. Numerical examples The algorithms presented in this paper have been implemented in Matlab 7.0.4.287 (R2007A) and run on an Intel(R) Core(TM)2 CPU T5600 laptop with a 1.83 GHz processor and 2046 Mb of RAM and applied to a wide variety of randomly generated polynomials u(x) and v(x) of degree n and m, respectively, where m < n. In the examples, we introduce input polynomials via a perturbation (add k × 10−13 to all entries of ‘‘exact’’ input polynomials with k is taken randomly in (−1, 1)) of another ‘‘exact’’ input polynomials, whose sequence of quotients and remainders is ‘‘exactly’’ known.
• deg: the degree of the computed approximate quotients and/or the approximate remainders. This result is particularly interesting when the degree of the approximate quotient and/or remainder sequence of the input polynomials is very sensitive to the choice of ϵ . • fails: when the algorithm does not recognize the ‘‘exact’’ polynomial quotient degree or the ‘‘exact’’ polynomial remainder degree. • The relative accuracy of the Hankel method and the Bézout method are, respectively: ErrorRHankel
=
H (u, v) − H (u, v)
1
∥H (u, v)∥1
,
ErrorRBe´ zout
=
B(u, v) − B(u, v) ∥B(u, v)∥1
1
.
(15)
• The absolute accuracy of the difference between the block diagonal matrices via Hankel and Bézout methods, respectively: ErrorAHankel = DHankel − DHankel 1 ,
DBe´ zout − DBe´ zout 1 . ErrorABe´ zout =
(16)
5.1. Parallelism steps The first example in this section is taken from [12]. Let u(x) = 6x9 + 24x8 + 44x7 + 162x6 + 60x5 + 273x4 + 32x3 + 193x2 − 70x − 10 and v(x) = 2x7 + 6x6 + 6x5 + 40x4 − 28x3 + 65x2 − 19x − 2. Thus, 0 12 36 36 J9 B (u, v) J9 = 240 −168 390 −114 −12
12 84 180 384 792 −282 1446 −468 −48
36 180 324 936 932 150 2006 −744 −68
36 384 963 1544 4992 −2722 8628 −2726 −264
240 792 932 4992 −1960 6756 16 −984 −60
−168 −282 150
−2722 6756
−8908 9041
−2447 −146
390 1446 2006 8628 16 9041 5037 −2714 −146
−114 −468 −744 −2726 −984 −2447 −2714 539 264
−12 −48 −68 −264 −60 . −146 −344 264
−50
We present a table (for every steps) which shows the parallelism between the approximate block factorization for H (u, v) and J9 B (u, v) J9 and the Euclidean Algorithm.
252
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
Step 1: r−1 (x) = r0 (x) q1 (x) − r1 (x), with qE1 (x) = 2.9999x2 + 3.0000x + 3.9999 2 qH 1 (x) = 3.0000x + 3.0000x + 3.9999
qB1 (x) = 2.9999x2 + 3.0000x + 3.9999 r1E (x) = −2.0000x4 − 5.9999x3 + 3.9999x2 − 11.9999x + 2.0000 r1H (x) = −2.0000x4 − 5.9999x3 + 3.9999x2 − 12.0000x + 1.9999 r1B (x) = −2.0000x4 − 5.9999x3 + 3.9999x2 − 11.9999x + 2.0000 Step 2: r0 (x) = r1 (x) q2 (x) − r2 (x), with qE2 (x) = −0.9999x3 + 3.7334e−11x2 − 4.9999x + 0.9999 3 2 qH 2 (x) = −0.9999x − 8.2589e−11x − 4.9999x + 0.9999
qB2 (x) = −0.9999x3 − 3.5625e−08x2 − 4.9999x + 0.9999 r2E (x) = r2H (x) = r2B (x) = −1.0000x2 − 2.9999x + 3.9999 Step 3: r1 (x) = r2 (x) q3 (x) − r3 (x), with qE3 (x) = 1.9999x2 − 3.5625e−8x + 3.9999 2 qH 3 (x) = 1.9999x + 6.4096e−7x + 3.9999
qB3 (x) = 1.9999x2 − 4.9382e−5x + 3.9999 r3E (x) = r3H (x) = r3B (x) = 13.9999 5.2. Tolerance–sensitive degree Let u(x) and v(x) two polynomials generated randomly of degree 54 and 40, respectively. When applying the classical Euclidean algorithm, there are four steps. The ‘‘exact’’ polynomial quotients degree {deg (q1 ) , deg (q2 ) , deg (q3 ) , deg (q4 )} and the ‘‘exact’’ polynomial remainders degree {deg (r1 ) , deg (r2 ) , deg (r3 ) , deg (r4 )} are {14, 13, 12, 14} and {27, 15, 1, 0}, respectively. Thus, we show the approximate quotients and remainders degrees founded for different values of ϵ . Step 1
ϵ
10−2 , . . . , 10−9 10−10 , . . . , 10−15 10−16 , . . . , 10−20 Step 2
ϵ
10 , . . . , 10 10−9 , . . . , 10−11 10−12 , . . . , 10−20 Step 3 −2
−8
ϵ
10 , . . . , 10 10−7 , . . . , 10−9 10−10 , . . . , 10−20 Step 4 −2
−6
ϵ
10 , . . . , 10 10−7 , 10−8 10−9 , . . . , 10−20 −2
−6
deg qB1 14 14 14
deg qH 1 14 14 14
deg qE1 14 14 14
deg r1B 27 27 27
deg r1H 27 27 27
deg r1E 27 fails fails
deg qB2 13 13 13
deg qH 2 13 13 fails
deg qE2 13 13 fails
deg r2B 15 15 15
deg r2H 15 fails fails
deg r2E 15 fails fails
deg qB3 12 12 12
deg qH 3 12 12 fails
deg qE3 12 12 fails
deg r3B 1 1 1
deg r3H 1 1 fails
deg r3E 1 fails fails
deg qB4 14 14 14
deg qH 4 14 1 fails
deg qE4 14 1 fails
deg r4B 0 0 0
deg r4H 0 0 fails
deg r4E 0 fails fails
The results show that the Bézout-based method is able to separate the approximate quotients and the approximate remainders degrees while the tolerance decreases in comparison with Euclid and/or Hankel-based methods. 5.3. Roots comparison In this section we compare the roots of the approximate quotients and remainders as defined in Section 5.2 with a tolerance ϵ = 10−5 . So, we describe the roots of the approximate quotients and remainders for the three methods in Fig. 1.
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
253
Fig. 1. Roots comparison of the approximate quotients and remainders as defined in Section 5.2 with a tolerance ϵ = 10−5 .
Thus, the example shows that the roots will not be distinguishable from each other. 5.4. Efficiency We measure the absolute error (16) with a fixed tolerance ε = 10−5 . We check this for several Hankel matrix H (u, v) and Bézout matrix Jn B(u, v)Jn associated with two polynomials u and v perturbed 50 times. Thereby, Fig. 2 represents the effectiveness of the block diagonalization using H (u, v) and Jn B(u, v)Jn of size 2448.
254
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
Fig. 2. ErrorAHankel (left figure) and ErrorABe´ zout (right figure) for n = 2448.
Table 1 The relative error in floating arithmetic with 15 digits. n
ErrorRHankel
ErrorRBe´ zout
1024 1224 1424 1624 1824 2048 2448
5.784074308862707e−012 7.377238376056904e−007 1.229383772472314e−007 1.402838670935590e−006 1.552356261630665e−007 1.482326667465912e−005 3.863794639419495e−005
1.086190389822859e−008 1.503623817980172e−007 5.902898664239044e−012 4.631154488009333e−007 8.760591568775624e−008 1.042891806603727e−009 1.162851904025527e−008
Fig. 3. The relative error (∗: Hankel, +:Bézout) for n = 1424.
5.5. Stability Now, we compare the approach based on H (u, v) and Jn B(u, v)Jn for their numerical stability. The main objective of this part is to improve the numerical stability of the approximate block diagonalization using Bézout matrices because the Bézout matrix coefficients are bounded by the Bézout Φ = 2nµν (with ∥u∥∞ ≤ µ and ∥v∥∞ ≤ ν ) and the Hankel matrix coefficients H (u, v) may increase exponentially. Also, if u and v are integer coefficients, the intermediate calculations steps can be preserved in integers numbers. See [6] for more details. Thus, we give the relative error (15), in floating arithmetic with 15 digits for H (u, v) and Jn B(u, v)Jn associated with two ‘‘exact’’ polynomials u and v . Thus, Table 1 shows the effectiveness of Algorithm 4.3 in terms of stability, specially when n increases. Moreover, Fig. 3 represents the stability behavior for the coefficients of 1424 × 1424 Hankel and Bézout matrices constructed with two polynomials.
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
255
6. Concluding remarks In this paper, we introduced a fast algorithm for computing the polynomials generated by the Euclidean Algorithm via an approximate block LU factorization of a Bézout matrix associated with the input polynomials. Numerical experiments performed with a wide variety of test problems, showed the effectiveness of this algorithm in terms of efficiency, stability and robustness. Acknowledgments I thank Prof. D. A. Bini for suggesting me this problem and guiding me throughout the course of this investigation and I owe a special debt of gratitude to Prof. H. Lombardi and Prof. G. M. Diaz-Toca for the helpful discussions. I also thank the anonymous referee for those helpful comments. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]
G.E. Collins, Subresultants and reduced polynomial remainder sequences, J. Assoc. Comput. Mach. 14 (1967) 128–142. W.S. Brown, On Euclid’s Algorithm and the computation of polynomial greatest common divisors, J. Assoc. Comput. Mach. 18 (1971) 478–504. W.S. Brown, J.F. Traub, On Euclid’s Algorithm and the theory of subresultants, J. Assoc. Comput. Mach. 18 (1971) 505–514. A. Schönhage, Quasi-GCD Computations, J. Complexity 1 (1) (1985) 118–137. M.-T. Noda, T. Sasaki, Approximate GCD and its application to ill-conditioned algebraic equations, J. Comput. Appl. Math. 38 (1–3) (1991) 335–351. D. Bini, V. Pan, Polynomial and Matrix Computations, Vol. 1, in: Fundamental Algorithms, Birkhäuser, 1994. D. Bini, L. Gemignani, Fast parallel computation of the polynomial remainder sequence via Bézout and Hankel matrices, SIAM J. Comput. 24 (1995) 63–77. L. Gemignani, GCD of polynomials and Bézout matrices, in: Proceedings of the 1997 International Symposium on Symbolic and Algebraic Computation (Kihei, Maui, HI), ACM Press, New York, 1997, pp. 271–277. D.A. Bini, L. Gemignani, Fast fraction-free triangularization of Bezoutians with applications to subresultant chain computation, Linear Algebra Appl. 284 (1–3) (1998) 19–39. W.B. Gragg, A. Lindquist, On partial realization problem, Linear Algebra Appl. 50 (1983) 277–319. G. Heining, K. Rost, Algebraic Methods for Toeplitz-Like Matrices and Operators, Birkhäuser Verlag, Basel, 1984. G.M. Diaz-Toca, N. Ben Atti, Block LU factorization of Hankel and Bézout matrices and Euclidean algorithm, Int. J. Comput. Math. 86 (2009) 135–149. N. Ben Atti, G.M. Diaz-Toca, Block diagonalization and LU-equivalence of Hankel matrices, Linear Algebra Appl. 412 (2006) 247–269. M. Mitrouli, N. Karcanias, C. Koukouvinos, Numerical performance of the matrix pencil algorithm computing the greatest common divisor of polynomials and comparison with other matrix-based methodologies, J. Comput. Appl. Math. 76 (1–2) (1996) 89–112. N. Karmarkar, Y.N. Lakshman, Approximate polynomial greatest common divisor and nearest singular polynomials, in: Proceedings of the International Symposium on Symbolic and Algebraic Computation, ACM Press, New York, 1996, pp. 35–39. R.M. Corless, S.M. Watt, L. Zhi, QR factoring to compute the GCD of univariate inexact polynomials, IEEE Trans. Signal Process. 52 (12) (2004) 3394–3402. C.J. Zarowski, X. Ma, F.W. Fairman, QR-factorization method for computing the greatest common divisor of polynomials with inexact coefficients, IEEE Trans. Signal Process. 48 (11) (2000) 3042–3051. V.Y. Pan, Computation of approximate polynomial GCDs and an extension, Inform. and Comput. 167 (2) (2001) 71–85. Z. Zeng, The approximate GCD of inexact polynomials, Part 1: a univariate algorithm, 2004, Preprint. B. Li, Z. Liu, L. Zhi, A structured rank-revealing method for Sylvester matrix, J. Comput. Appl. Math. 213 (1) (2008) 212–223. J.R. Winkler, J.D. Allan, Structured total least norm and approximate GCDs of inexact polynomials, J. Comput. Appl. Math. 215 (1) (2008) 1–13. J.R. Winkler, X. Lao, The calculation of the degree of an approximate greatest common divisor of two polynomials, J. Comput. Appl. Math. 235 (6) (2011) 1587–1603. S. Belhaj, A fast method to block-diagonalize a Hankel matrix, Numer. Algorithms 47 (2008) 15–34. S. Belhaj, Computing the block factorization of complex Hankel matrices, Computing 87 (3–4) (2010) 169–186. S. Belhaj, Block factorization of a Hankel matrix and Euclidean Algorithm, Math. Model. Nat. Phenom. 5 (7) (2010) 38–44. P. Lancaster, M. Tismenetsky, The Theory of Matrices, second ed., in: Computer Science and Applied Mathematics, Academic Press Inc., Orlando, FL, 1985. Fu-Rong Lin, Wai-Ki Ching, M.K. Ng, Fast inversion of triangular Toeplitz matrices, Theoret. Comput. Sci. 315 (2004) 511–523. V. Hribernig, Sensitivity of algebraic algorithms, Ph.D. Thesis at the Technical University of Vienna, 1994. V. Hribernig, H.J. Stetter, Detection and validation of clusters of polynomial zeros, J. Symbolic Comput. 24 (6) (1997) 667–681.
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Computing the polynomial remainder sequence via Bézout matrices Skander Belhaj ∗ Manouba University, ISAMM, Campus Universitaire de la Manouba, 2010 Tunis, Tunisia University of Tunis El Manar, ENIT-LAMSIN, BP 37, 1002 Tunis, Tunisia
article
abstract
info
Article history: Received 22 December 2011 Received in revised form 15 February 2013
In this paper, we address the task of computing the polynomial remainder sequence appearing in the Euclidean algorithm applied to two polynomials u(x) and v(x) of degree n and m, respectively, m < n via renewed approach based on the approximate block diagonalization of a Bézout matrix B (u, v). More specifically, this algorithm consists of using the efficient Schur complementation method given in terms of a Bézout matrix B(u, v) associated with the input polynomials. We also compare the proposed approach to the approximate block diagonalization of a Hankel matrix H (u, v) in which the successive transformation matrices are upper triangular Toeplitz matrices. All algorithms have been implemented in Matlab. Numerical experiments performed with a wide variety of test problems show the effectiveness of this algorithm in terms of efficiency, stability and robustness. © 2013 Elsevier B.V. All rights reserved.
MSC: 65Fxx Keywords: Hankel matrix Bézout matrix Block diagonalization Euclidean algorithm Polynomial remainder sequence Schur complementation
1. Introduction Let u(x) =
n k=0
uk x k
and
v(x) =
m
v k xk
(1)
k=0
be two polynomials with deg (u(x)) = n, deg (v(x)) = m and m < n. The classical Euclidean algorithm applied to u(x) and v(x) returns a sequence of quotients qk (x) and remainders rk (x), such that r−1 (x) = u(x), r0 (x) = v(x) rk−2 (x) = rk−1 (x) qk (x) − rk (x) ,
k = 1, . . . , K
(2)
where −rk (x) is the polynomial remainder of the division of rk−2 (x) and rk−1 (x) and rK (x) is the greatest common divisor (GCD) of u(x) and v(x). The question of computing the quotients and the remainders is beginning to be studied by Collins and Brown (see [1–3]) via efficient Euclidean-based methods. Naturally, the classical algorithm turns out to be unstable when applied to approximate polynomials. One might think of designing stabilized numerical versions of the Euclidean algorithm to the approximate case (see [4,5]). Several papers have been devoted to describe the natural relation between the Euclidean algorithm and the block LU factorization of the Hankel and Bézout matrices. Sure enough, Bini and Gemignani (see [6–9])
∗
Correspondence to: Manouba University, ISAMM, Campus Universitaire de la Manouba, 2010 Tunis, Tunisia. E-mail address: [email protected].
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.02.029
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
245
compute the coefficients of the polynomials in (2) via renewed approach based on a block LU factorization (see [10,11]) of the Hankel matrix H (u, v), associated with u(x) and v(x). More specifically, At H (u, v) A = D
(3)
where A is an upper triangular matrix with entries equal to one and D is a block diagonal matrix and each block is a lower Hankel triangular matrix. Lately, Ben Atti and Diaz-Toca [12] calculate the coefficients of the polynomials in (2) via a different block LU factorization of the Hankel matrices [13] from the classical one. Naturally, when the input polynomials (1) are introduced via a perturbation (add k × 10−13 to all entries of ‘‘exact’’ input polynomials with k is taken randomly in (−1, 1)), the beautiful relation between the Euclidean algorithm and the block LU factorization of the Hankel and Bézout matrices may be lost. In addition, the polynomial remainder sequence (prs), is well known to be an ill-conditioned problem. In fact, small perturbations of the input polynomials may generate big variations in the sequence (even if computed exactly). Moreover, the floating point computation of the polynomial remainder sequence, performed by means of subsequent division, is a numerically unstable computation. In fact, just after the first step, the errors introduced by the floating point arithmetic provide a perturbation to the input of the subsequent steps which form again a prs which is an ill-conditioned problem. Thus these errors are largely amplified and spoil the final result. This is why we try to generalize this to the computation of approximate versions of these polynomials and factorizations. For more details on this research topic we refer the reader to [14–22] and to the references therein. It is shown that an n × n Hankel matrix admits an approximate block diagonalization in which the successive transformation matrices are upper triangular Toeplitz matrices [23,24]. More explicitly, lHϵm1
Θ1,2
Θ2,1 Dϵ = Aϵ HAϵ = . ..
lHϵm2
t
Θn,1
..
. ···
··· .. . .. . Θn,n−1
Θ1,n .. . Θn−1,n
(4)
lHϵml
where every block lHϵmj , for m1 + m2 + · · · + ml = n, l = 1, 2, . . . is an approximate lower Hankel matrix, all entries of Θj,k ,
for j, k = 1, 2, . . . , n and j ̸= k are very close to zero and Aϵ is an approximate upper triangular matrix. This algorithm is revised in order to obtain the approximate quotients polynomials and the approximate remainders polynomials appearing in the Euclidean Algorithm applied to two polynomials u(x) and v(x) of degree n and m, respectively, m < n [25]. The main contribution of this paper is to ensure the natural relation between the Euclidean algorithm and the approximate block LU factorization diagonalization of Hankel and Bézout matrices associated with u(x) and v(x) and to propose an efficient Schur-based algorithm of the corresponding Bézout matrix B(u, v) (see the works by Bini and Gemignani [6–9]) which provide the approximate quotients and the approximate remainders during the process execution. The paper is organized as follows: Section 2 gives some notations and definitions. Some theoretical results associated with the approximate block diagonalization of Hankel matrices are recalled in Section 3. Section 4 presents an efficient algorithm for an approximate block diagonalization of n × n Bézout matrices applied to two polynomials u(x) and v(x) of degree n and m, respectively, m < n and relates this approach to the Euclidean Algorithm in order to obtain the approximate quotients and the approximate remainders appearing during the process execution. In Section 5 we illustrate our method with examples given by an implementation of the procedures using Matlab and we propose a statistical study performed with a wide variety of test problems to show the effectiveness of this algorithm in terms of speed, stability and robustness. Finally, a summary and future research are given in Section 6 to complete the paper. 2. Notations and definitions We begin the section introducing the following notations:
• H (S ) ∈ Rn×n denotes the Hankel matrix associated with a list S of length (2n − 1). This means that the first row is given by the first n terms of S and the last column is given by the last n terms of S.
• lH (S ) ∈ Rn×n denotes the Hankel triangular matrix (with respect to the antidiagonal) associated with a list S of length (2n − 1) such that the last column is defined by S. • uH (S ) ∈ Rn×n denotes the Hankel triangular matrix (with respect to the antidiagonal) associated with a list S of length (2n − 1) such that the first column is defined by S. • H (S ; m; n) ∈ Rm×n denotes the Hankel matrix associated with a list S of length (m + n − 1). The first row is given by the first n terms of S and the last column is given by the last m terms of S.
• T (S ; m; n) ∈ Rm×n denotes the Toeplitz matrix associated with a list S of length (m + n − 1). The first row is given by the first m terms of S and the last column is given by the last n terms of S.
• uT (S ) ∈ Rn×n denotes the upper Toeplitz triangular matrix associated with a list S such that the first row is defined by S.
• lT (S ) ∈ Rn×n denotes the lowerToeplitz triangular matrix associated with a list S such that the last row is defined by S. p • Let p ∈ N. Let Σp ∈ Rp×p , Σp = εjk j, k=1 , where all entries of Σp are zero except that εj+k,j = εk for j+k, k = 1, 2, . . . , p.
246
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
p are zero except that εj+k,j for j + k, k = 1, 2, . . . , p. p ∈ Rp×p , Σ = εjk p , where all entries of Σ • Let p ∈ N. Let Σ j, k=1 • Jp = lH (1, 0, . . . , 0), p ∈ N. p−1
• • • •
Given P ∈ Rn×m , P = Jm P t Jn . T (S ; m; n) = Jm H (S ; m; n). uT (S ) = Jl lH (S ), where l is the length of L. Let a ∈ R. Let µ > 0, V (a, µ) = (a − µ; a + µ) is a neighborhood of a.
2.1. Hankel matrices k k Let u(x) = k=0 uk x and v(x) = k=0 vk x be two polynomials of degree n and m, respectively, where m < n. The power series expansion R (x) of the function v(x)/u(x) at the infinity
n
R ( x) =
∞
m
hk x −k
(5)
k=0
defines the n × n Hankel matrix, H = H (u, v), associated with u(x) and v(x),
.
.. .
··· ··· .. .
hn
hn +1
···
h1 h2
H = ..
h2 h3
hn hn+1 h2n−1
.. . .
(6)
In addition, every nonsingular Hankel matrix can be viewed as a Hankel matrix associated with the given polynomials. Proposition 2.1. For any nonsingular n × n Hankel matrix H there exists two coprime polynomials u(x) and v(x) of degree n and m, respectively, where m < n, such that H = H (u, v). The polynomial u(x) and v(x) are related to H by the following equalities: H (u0 , . . . , un−1 )t = −un (hn+1 , . . . , h2n ) ,
(vn−1 , . . . , v0 )t = lT (h1 , . . . , hn ) (un , . . . , u1 )t where h2n is any number and lT (h1 , . . . , hn ) is the lower triangular Toeplitz matrix defined by the list [h1 , . . . , hn ]. Proof (See [7]).
Remark 2.1. Suppose that H has the following structure: H = lH (hn , . . . , h2n−1 ) . Thus, Proposition 2.1 concludes, for any polynomial v(x),
0
(vn−1 , . . . , v0 )t = ...
hn
..
. ···
t t (un , . . . , u1 ) = (0, . . . , 0, hn un ) . 0
If v(x) = 1 then v0 = 1 and un = 1/hn . Thus, H represents a Hankel matrix associated with u(x) and v(x) = 1 : H = H (u, 1). Suppose that we know all the entries of a Hankel matrix H and all the coefficients of u(x). Thus, the algorithm which provides all the coefficients of v(x) can be formulated as follows: Algorithm 2.1 (Construct v(x) from u(x) and a Hankel Matrix H). Given a Hankel matrix H and u(x) = m k polynomial of degree n, this algorithm provides all the coefficients of v(x) = k=0 vk x . 1. 2. 3. 4. 5.
Set H = [hn−m+2 , . . . , hn−1 ] and U = [u0 , . . . , un−1 , un ]. Set n = length (U ), U0 = u0 and U = [u1 , . . . , un−1 , un ]. Set n = length (U ). Compute U = Jn ∗ U t . Compute V = (Jn ∗ lT (H (1, 1 : n)) ∗ U )t .
n
k=0
uk x k a
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
247
2.2. Bézout matrices Let u(x) =
n
k =0
b (x, y) =
uk xk and v(x) =
m
k=0
u (x) v (y) − u (y) v (x) x−y
vk xk be two polynomials of degree n and m, respectively, where m < n. Let
.
The rational function b (x, y) is easily proved to be a polynomial in two variables and can therefore be written as b (x, y) = n−1 j k b x y . The matrix B x , y = b ( ) j,k 0≤j,k≤n−1 is called the Bézout matrix or Bezoutian associated with u(x) and v(x). j, k=0 j,k We next state more different representations of the Bézout matrix. Proposition 2.2. The Bezoutian matrix can be written as: B (u, v) = uH (u1 , . . . , un ) uT (v0 , . . . , vn−1 ) − uH (v1 , . . . , vn ) uT (u0 , . . . , un−1 ) where if m < n the coefficients vm+1 , . . . , vn are assumed to be zero. Proof (See [26]).
As a consequence of the previous result, we have: Corollary 2.1. The entries of B (u, v) can be recursively calculated as follows: bj,k+1 = bj+1,k + uj vk − uk vj setting bn+1,k = b0,k = 0 for all j, k and assuming vm+1 = · · · = vn = 0 as above. From Proposition 2.2, we observe that the matrix B (u, 1) is an upper Hankel triangular matrix B (u, 1) = uH (u1 , . . . , un ). Moreover, we see that B (u, 1)−1 = H (u, 1). Suppose now that we know all the entries of a Bézout matrix B(u, v) and all the coefficients of u(x). We finish the section proposing the algorithm which provides all the coefficients of v(x) as follows: Algorithm 2.2 (Construct v(x) from u(x) and a Bézout matrix B(u, v)). Given a Bézout matrix B(u, v) and u(x) = m k a polynomial of degree n, this algorithm provides all the coefficients of v(x) = k=0 vk x . 1. 2. 3. 4.
n
k =0
uk x k
Set U = [u0 , . . . , un ] and n = length (U ). Compute B(u, 1) = uH (u1 , . . . , un ). Compute H = H (u, 1) ∗ Jn B(u, v)Jn ∗ H (u, 1). Call Algorithm 2.1 to conclude V = [v0 , . . . , vm ].
3. Approximate block diagonalization of a Hankel matrix We introduce the approximate block diagonalization of a Hankel matrix presented in [23]. Lemma 3.1. Let n ∈ N ∗ . Let h = H (h1 , . . . , h2n−1 ) be a square Hankel matrix of order n. Suppose that hj = εj with εj ∈ V (0, µ) for j = 1, 2, . . . , p − 1 and hp ̸∈ V (0, µ). Then h has the form h = H ε1 , . . . , εp−1 , hp , . . . , h2n−1 =
h11 h21
h12 h22
,
(7)
where h11 = lH hp , . . . , h2p−1 + Jp Σp ,
h22 = H h2p+1 , . . . , h2n−1 ,
h12 = H hp+1 , . . . , hn+p−1 ; p; n − p ,
h21 = ht12 .
We can successively construct from h the following two matrices:
• A square lower Hankel triangular matrix H of order (2n − p), 0 0 H13 0 h11 h12 , H = lH hp , . . . , h2n−1 = H31
ht12
where H31 = H13 = lH (hp , . . . , hn−1 ). • A square upper triangular Toeplitz matrix T , t11 T = J2n−p H =uT hp , . . . , h2n−1 = 0 0
t12 t22 0
t13 t23 t33
,
where t11 = t33 = uT hp , . . . , hn−1 , t22 = Jp h11 , t13 = Jn−p h22 , t12 = Jn−p ht12 , and t23 = Jp h12 .
(8)
h22
(9)
248
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
Lemma 3.2. Let T be an upper triangular Toeplitz matrix, nonsingular with non-zero diagonal. Then T −1 = uT µ1 , . . . , µ2n−p and has the following block decomposition,
T −1 −1 T = 0 0
−1 T −1 12
11
T
−1 −1 t11 T −1 13 = 0 T −1 23
22
0
T
33
P
Q P , −1 t11
−1 t22 0
0
(10)
where P = T (µ2 , . . . , µn ; p; n − p) ,
P = Jn−p P t Jp ,
t22 P + t23 t11 = 0(p,n−p) ,
−1 h11 P + Mt11 = 0(p,n−p) ,
−1 t11 P + t12 t22 = 0(n−p,p) ,
−1 t11 Q + t12 P + t13 t11 = 0(n−p,n−p) .
−1
Theorem 3.1. Let h = H ε1 , . . . , εp−1 , hp , . . . , h2n−1 be an approximate Hankel matrix (approximation of an ‘‘exact’’ Hankel
matrix, i.e. H 0, . . . , 0, hp , . . . , h2n−1 ) where εj ∈ V (0, µ) for j = 1, 2, . . . , p − 1 with hp ̸∈ V (0, µ) and
t = uT hp , . . . , hn+p−1 ,
t −1 = uT (µ1 , . . . , µn ) .
Then h′ = t −1
t
ht −1 =
ϵ′
h′
11 t ϵ′
h′22
,
(11)
where −1 h′11 = lH µ1 , . . . , µp + t22
t
−1 Jp Σp t22 ,
h′22 = −H µp+2 , . . . , µ2n−p + P t Jp Σp P ,
−1 t ϵ ′ = t22 J p Σp P .
Remark 3.1. Theorem 3.1 gives an approximate reduction for real Hankel matrix. Thus if we have εj = 0, Theorem 3.1 provides the exact case. To iterate this result, h′22 must be a Hankel matrix, but −H µp+2 , . . . , µ2n−p is a Hankel matrix and P t Jp Σp P is only a symmetric matrix. If we choose εj very close to zero, we can conclude that Σp ≈ 0p , and the process converges to the exact case. To solve this problem, we can choose h′22 = −H µp+2 , . . . , µ2n−p + Θ where Θ is a Hankel matrix built with the first column and the last row of P t Jp Σp P. Remark 3.2. The key of our approach is based on the inversion of triangular Toeplitz matrix. Thus, we use superfast techniques described in [27]. Then, we give the fast algorithm to compute the block factorization as follows: Algorithm 3.1 (Approximate Block Diagonalization for Hankel Matrix). Given a list S = [ε1 , . . . , εp−1 , hp , . . . , h2n−1 ] where εj ∈ V (0, µ) for j = 1, 2, . . . , p − 1 with hp ̸∈ V (0, µ), this algorithm computes the approximate block diagonalization of Hankel matrix via upper triangular Toeplitz matrix. 1. Define a Hankel matrix h = H ε1 , . . . , εp−1 , hp , . . . , h2n−1 .
2. Define an upper triangular Toeplitz matrix t = uT hp , . . . , hn+p−1 .
3. 4. 5. 6.
Compute t −1 . t Compute h′ = t −1 ht −1 . Set h′11 = h′ (1 : p, 1 : p) and h′22 = h′ (p + 1 : n, p + 1 : n). Recursively apply Algorithm 3.1 to h = H h′22 (1 : n − p, 1) h′22 (n − p, 1 : n − p) ,
obtaining (4). 4. Block diagonalization of Hankel and Bézout matrices: connection to the Euclidean algorithm The topic of this section is the relationship between the Euclidean algorithm applied to a pair of polynomials u(x) and v(x) and the LU factorization of the associated Bézout or Hankel matrix. The main result essentially states that the entries of the approximate triangular matrices of the block triangular factorization of H (u, v) or Jn B(u, v)Jn yield the coefficients of all the polynomials generated by the Euclidean scheme.
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
249
4.1. The Euclidean algorithm Recall that the classical Euclidean algorithm for the computation of the GCD goes as follows. Let u(x) and v(x) be two polynomials in R [x] with deg (u(x)) = n, deg (v(x)) = m and m < n. The classical Euclidean algorithm applied to u(x) and v(x) returns a sequence of quotient qk (x) and remainder rk (x), such that r0 (x) = v(x) r−1 (x) = u(x), (12) rk−2 (x) = rk−1 (x) qk (x) − rk (x) , k = 1, . . . , K where −rk (x) is the remainder of the division of rk−2 (x) and rk−1 (x). The algorithm stops when a remainder rK +1 (x) = 0 is found; then rK (x) is the desired greatest common divisor (GCD) of u(x) and v(x). Observe that, since the degrees of the rk (x)’s are strictly decreasing, the algorithm must terminate in a finite number of steps. One might think of adapting the Euclidean algorithm to the approximate case by simply converting all ‘‘small’’ coefficients into 0 after each polynomial division (where ‘‘small’’ coefficients are coefficients whose absolute value is smaller than a fixed tolerance). Unfortunately, this naive method is highly unstable, as was already observed in [4]. It is possible, however, to design stabilized versions of the Euclidean algorithm. A first step in this direction is found in [5], where the authors introduce a normalization rule on remainders and calculate the approximate GCD as follows. Algorithm 4.1 (Computing the Polynomial Remainder Sequence and the Polynomial Quotients Via Euclid). Given u(x) and v(x) two polynomials in R [x] of degree n and m, respectively, where m < n, this algorithm computes the approximate quotients and the approximate remainders with a fixed tolerance less than a small positive number ϵ, 0 < ϵ ≪ 1. 1. Compute qk as the quotient of rk−1 divided by rk . 2. Set rk−1 = qk rk + max {1, ∥qk ∥∞ } · rk+1 . 3. Recursively apply Algorithm 4.1 for k = 0, 1, . . . , K and stops when a remainder |rK +1 | < ϵ . Remark 4.1. It must be pointed out that any variant of the approximate Euclidean algorithm (see for example [28,29]) can be used. We propose to implement Algorithm 3.1 to make a comparison with our method (see Section 5). 4.2. Approximate reduction of a Hankel matrix associated with two polynomials We recall the correlation between the approximate Euclidean algorithm applied to two polynomials u(x) and v(x) and the approximate block diagonalization of a Hankel matrix [25]. Corollary 4.1. Let H (u, v) = H (ε1 , . . . , εn−m−1 , hn−m , . . . , h2n−1 ) be an approximate Hankel as defined m matrix n k k in Theorem 3.1 associated with two coprime polynomials in R [x] , u(x) = k=0 uk x and v(x) = k=0 vk x , deg (u(x)) = n, deg (v(x)) = m and m < n. Let t = uT (hn−m , . . . , h2n−m−1 ) ,
t −1 = uT (µ1 , . . . , µn ) .
Then
t
−1 t
J B (q, 1) J −1 H (u, v) t = n−m ′ t n−m ϵ
ϵ′ , H (v, r )
(13)
where −1 J B (q, 1) J = Jn−m B (q, 1) Jn−m + Jn−m t22
H (v, r ) = H (v, r ) + P t Jn−m Σn−m P ,
t
−1 Jn−m Σn−m t22 Jn−m ,
−1 t ϵ ′ = t22 Jn−m Σn−m P ,
q(x) and r (x) are the polynomial quotient and remainder of the division u (x) /v (x). Let us devise the algorithm for the computation of the approximate polynomials quotients and the approximate polynomials remainders appearing in the Euclidean algorithm applied to two polynomials u(x) and v(x) of degree n and m, respectively, m < n based on the above corollary. Algorithm 4.2 (Computing Remainder Sequence and the Polynomial Quotients Via H (u, v)). Given u(x) = n m the Polynomial k k k=0 uk x and v(x) = k=0 vk x two polynomials of degree n and m, respectively, where m < n, this algorithm computes the approximate polynomials quotients and the approximate polynomials remainders appearing in the Euclidean Algorithm 1. 2. 3. 4. 5. 6. 7.
Construct a Hankel matrix H = H (0, . . . , 0, hn−m , . . . , h2n−1 ). Define an upper triangular matrix t = uT (hn−m , . . . , h2n−m−1 ). Toeplitz t Compute t −1 and h′ = t −1 ht −1 . Set h′11 = h′ (1 : n − m, 1 : n − m) and h′22 = h′ (n − m + 1 : n, n − m + 1 : n). Recover the coefficients of the quotient polynomial from step 3. Recover the coefficients of the remainder polynomial from Algorithm 2.1 by using steps 1. and 5. Recursively apply Algorithm 4.2 to h = H h′22 (1 : m, 1) h′22 (m, 1 : m) , obtaining all the coefficients appearing in the Euclidean Algorithm.
250
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
4.3. Approximate reduction of a Bézout matrix associated with two polynomials We explore the relationship between the approximate Euclidean algorithm applied to two polynomials u(x) and v(x) and the approximate block diagonalization of a Bézout matrix (see the works by Bini and Gemignani [6–9]). The following proposition explains how the polynomials generated by the approximate Euclidean algorithm can be obtained through Schur complements of the approximate permuted Bezoutian Jn B (u, v) Jn =
Pm
i
X
Xt Y
,
Pmi = Pmi + J Σ . where Pmi = (Jn B(u, v)Jn )mi and k k Proposition 4.1. Let u(x) = k=0 uk x and v(x) = k=0 vk x be two coprime polynomials with deg (u(x)) = n, deg (v(x)) = m and m < n. Let H be the associated Hankel matrix as in Theorem 3.1 and let 0 = m0 < m1 < · · · < ml = n be integers Pmi ̸∈ V (0, µ) for i = 1, . . . , l and such that det Hmi ̸∈ V (0, µ) for i = 0, . . . , l and det Hk ∈ V (0, µ) otherwise. Then det det Pmi ∈ V (0, µ) for j ̸= mi , i = 1, . . . , l. Moreover, if Smi is the Schur complement of Pmi and Smi is the Schur complement of Pmi in Jn B(u, v)Jn , the following relation holds:
n
m
−1 t Smi = Y − XPm X = Jn B(ri , ri+1 )Jn i
−1 −1 and Smi = Smi + X t Pm J Σ X t Pm i i
t
,
where {ri (x)}1,...,l is the polynomial remainder sequence obtained by applying the Euclidean algorithm to u(x) and v(x). Proof. We may inductively apply the same proof in [7].
Corollary 4.2. Under the hypothesis of Proposition 4.1. Then Jn B(u, v)Jn =
=
I −1 XPm i
0 I
T
0 I
−1 XPm T i
Θ Sm
Pmi Θt
−1 t Pm X i I
I 0
i
(i)
Jn B(u , v (i) )Jn + J Σ
Θ Sm
t Θ
i
Tt 0
−1 t T t Pm X i I
(14)
where
−1 −1 t Smi = Smi + X t Pm J Σ X t Pm , i i −1 t (i) −1 = T −1 Θ , Jn , Θ = −J Σ XPmi , Θ T = lT un , . . . , un−mi +1 B(u , 1) and the polynomials u(i) (x) , v (i) (x) verify H (u, v)mi = H u(i) , v (i) . −1 t Smi = Y − XPm X = Jn B(ri , ri+1 )Jn , i
Proof. The corollary follows from [7] and Proposition 4.1.
From Corollary 4.2, we present a parallel algorithm for an approximate block diagonalization of Jn B(u, v)Jn . Let us compare this method with Algorithm 4.2. We start computing the first diagonal block (a lower Hankel matrix). We denote by DBe´ zout the block diagonal matrix of Jn B(u, v)Jn . Let D11 be the first block. Then m1 = n − m, h11 = H(u, v)n−m = lH (un−m , . . . , u2n−2m−1 ) and so, by Remark 2.1, we can choose v (1) (x) = 1 such that H (u, v)n−m = H u(1) , 1 . Thus,
Pmi = Jn B(u(1) , v (1) )Jn + J Σ = Jn B(u(1) , 1)Jn + J Σ 1 = Jn H (u(1) , 1)−1 Jn + J Σ = Jn h− 11 Jn + J Σ = Jn B(q1 , 1)Jn + J Σ .
So, D11 = Jn B(q1 , 1)Jn + J Σ . Successively, we obtain an approximate upper matrix ABe´ zout and DBe´ zout such that AtBe´ zout Jn B(u, v)Jn ABe´ zout = DBe´ zout . Moreover the successive matrices B(ri , ri+1 ) directly provide the approximate remainders. Thus, we obtain the same blocks Dii for both Bézout and Hankel matrices, except the last one. Therefore, the coefficient vector of the (i + 1)-th approximate remainder ri+1 can be expressed as
−1 t −1 −1 ri+1 = Jn Y − XPm X Jn + Jn X t Pm J Σ X t Pm i i i
t
Jn en−mi ,
that is, ri+1 can be computed by solving an mi × mi system with Hankel-like matrix Pmi . Thus, we give an algorithm which computes the approximate polynomials quotients and the approximate polynomials remainders appearing in the Euclidean algorithm applied to two polynomials u(x) and v(x) of degree n and m, respectively, m < n based on the above corollary.
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
251
Algorithm 4.3 (Computing the Polynomial Remainder Sequence and the Polynomial Quotients via Jn B (u, v) Jn ). Given u(x) = n m k k k=0 uk x and v(x) = k=0 vk x two polynomials of degree n and m, respectively, where m < n, this algorithm computes the approximate polynomials quotients and the approximate polynomials remainders appearing in the Euclidean Algorithm. 1. Construct the Bézout matrix B (u, v). 2. Reduce Jn B (u, v) Jn as (14).
t
−1 t −1 −1 −1 X + X t Pm J Σ X t Pm . Pm and Sm1 = Y − XPm 3. Compute Pm1 = Jn B(u(1) , v (1) )Jn + J Σ , 1 1 1 1
4. Recover the coefficients of the first polynomial quotient from Algorithm 2.2 applied to Pm1 = Jn B(u(1) , 1)Jn + J Σ . 5. Recover the coefficients of the first polynomial remainder from Algorithm 2.2 applied to Sm1
−1 −1 X t Pm J Σ X t Pm 1 1
t
= Jn B(r1 , r2 )Jn +
with r1 = v (1) .
6. Recursively apply Algorithm 4.3 to Jn B(u, v)Jn = Smi , obtaining all the coefficients appearing in the Euclidean Algorithm. Remark 4.2. Observe that the assumption in Proposition 4.1 that u and v are coprime is no loss of generality. Indeed, if rL = GCD(u, v) were a polynomial of positive degree, then the Euclidean remainder sequence {ri+1 } can be recovered from the sequence { ri+1 } generated by the Euclidean algorithm applied to u(x)/rL (x) and v(x)/rL (x). 5. Numerical examples The algorithms presented in this paper have been implemented in Matlab 7.0.4.287 (R2007A) and run on an Intel(R) Core(TM)2 CPU T5600 laptop with a 1.83 GHz processor and 2046 Mb of RAM and applied to a wide variety of randomly generated polynomials u(x) and v(x) of degree n and m, respectively, where m < n. In the examples, we introduce input polynomials via a perturbation (add k × 10−13 to all entries of ‘‘exact’’ input polynomials with k is taken randomly in (−1, 1)) of another ‘‘exact’’ input polynomials, whose sequence of quotients and remainders is ‘‘exactly’’ known.
• deg: the degree of the computed approximate quotients and/or the approximate remainders. This result is particularly interesting when the degree of the approximate quotient and/or remainder sequence of the input polynomials is very sensitive to the choice of ϵ . • fails: when the algorithm does not recognize the ‘‘exact’’ polynomial quotient degree or the ‘‘exact’’ polynomial remainder degree. • The relative accuracy of the Hankel method and the Bézout method are, respectively: ErrorRHankel
=
H (u, v) − H (u, v)
1
∥H (u, v)∥1
,
ErrorRBe´ zout
=
B(u, v) − B(u, v) ∥B(u, v)∥1
1
.
(15)
• The absolute accuracy of the difference between the block diagonal matrices via Hankel and Bézout methods, respectively: ErrorAHankel = DHankel − DHankel 1 ,
DBe´ zout − DBe´ zout 1 . ErrorABe´ zout =
(16)
5.1. Parallelism steps The first example in this section is taken from [12]. Let u(x) = 6x9 + 24x8 + 44x7 + 162x6 + 60x5 + 273x4 + 32x3 + 193x2 − 70x − 10 and v(x) = 2x7 + 6x6 + 6x5 + 40x4 − 28x3 + 65x2 − 19x − 2. Thus, 0 12 36 36 J9 B (u, v) J9 = 240 −168 390 −114 −12
12 84 180 384 792 −282 1446 −468 −48
36 180 324 936 932 150 2006 −744 −68
36 384 963 1544 4992 −2722 8628 −2726 −264
240 792 932 4992 −1960 6756 16 −984 −60
−168 −282 150
−2722 6756
−8908 9041
−2447 −146
390 1446 2006 8628 16 9041 5037 −2714 −146
−114 −468 −744 −2726 −984 −2447 −2714 539 264
−12 −48 −68 −264 −60 . −146 −344 264
−50
We present a table (for every steps) which shows the parallelism between the approximate block factorization for H (u, v) and J9 B (u, v) J9 and the Euclidean Algorithm.
252
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
Step 1: r−1 (x) = r0 (x) q1 (x) − r1 (x), with qE1 (x) = 2.9999x2 + 3.0000x + 3.9999 2 qH 1 (x) = 3.0000x + 3.0000x + 3.9999
qB1 (x) = 2.9999x2 + 3.0000x + 3.9999 r1E (x) = −2.0000x4 − 5.9999x3 + 3.9999x2 − 11.9999x + 2.0000 r1H (x) = −2.0000x4 − 5.9999x3 + 3.9999x2 − 12.0000x + 1.9999 r1B (x) = −2.0000x4 − 5.9999x3 + 3.9999x2 − 11.9999x + 2.0000 Step 2: r0 (x) = r1 (x) q2 (x) − r2 (x), with qE2 (x) = −0.9999x3 + 3.7334e−11x2 − 4.9999x + 0.9999 3 2 qH 2 (x) = −0.9999x − 8.2589e−11x − 4.9999x + 0.9999
qB2 (x) = −0.9999x3 − 3.5625e−08x2 − 4.9999x + 0.9999 r2E (x) = r2H (x) = r2B (x) = −1.0000x2 − 2.9999x + 3.9999 Step 3: r1 (x) = r2 (x) q3 (x) − r3 (x), with qE3 (x) = 1.9999x2 − 3.5625e−8x + 3.9999 2 qH 3 (x) = 1.9999x + 6.4096e−7x + 3.9999
qB3 (x) = 1.9999x2 − 4.9382e−5x + 3.9999 r3E (x) = r3H (x) = r3B (x) = 13.9999 5.2. Tolerance–sensitive degree Let u(x) and v(x) two polynomials generated randomly of degree 54 and 40, respectively. When applying the classical Euclidean algorithm, there are four steps. The ‘‘exact’’ polynomial quotients degree {deg (q1 ) , deg (q2 ) , deg (q3 ) , deg (q4 )} and the ‘‘exact’’ polynomial remainders degree {deg (r1 ) , deg (r2 ) , deg (r3 ) , deg (r4 )} are {14, 13, 12, 14} and {27, 15, 1, 0}, respectively. Thus, we show the approximate quotients and remainders degrees founded for different values of ϵ . Step 1
ϵ
10−2 , . . . , 10−9 10−10 , . . . , 10−15 10−16 , . . . , 10−20 Step 2
ϵ
10 , . . . , 10 10−9 , . . . , 10−11 10−12 , . . . , 10−20 Step 3 −2
−8
ϵ
10 , . . . , 10 10−7 , . . . , 10−9 10−10 , . . . , 10−20 Step 4 −2
−6
ϵ
10 , . . . , 10 10−7 , 10−8 10−9 , . . . , 10−20 −2
−6
deg qB1 14 14 14
deg qH 1 14 14 14
deg qE1 14 14 14
deg r1B 27 27 27
deg r1H 27 27 27
deg r1E 27 fails fails
deg qB2 13 13 13
deg qH 2 13 13 fails
deg qE2 13 13 fails
deg r2B 15 15 15
deg r2H 15 fails fails
deg r2E 15 fails fails
deg qB3 12 12 12
deg qH 3 12 12 fails
deg qE3 12 12 fails
deg r3B 1 1 1
deg r3H 1 1 fails
deg r3E 1 fails fails
deg qB4 14 14 14
deg qH 4 14 1 fails
deg qE4 14 1 fails
deg r4B 0 0 0
deg r4H 0 0 fails
deg r4E 0 fails fails
The results show that the Bézout-based method is able to separate the approximate quotients and the approximate remainders degrees while the tolerance decreases in comparison with Euclid and/or Hankel-based methods. 5.3. Roots comparison In this section we compare the roots of the approximate quotients and remainders as defined in Section 5.2 with a tolerance ϵ = 10−5 . So, we describe the roots of the approximate quotients and remainders for the three methods in Fig. 1.
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
253
Fig. 1. Roots comparison of the approximate quotients and remainders as defined in Section 5.2 with a tolerance ϵ = 10−5 .
Thus, the example shows that the roots will not be distinguishable from each other. 5.4. Efficiency We measure the absolute error (16) with a fixed tolerance ε = 10−5 . We check this for several Hankel matrix H (u, v) and Bézout matrix Jn B(u, v)Jn associated with two polynomials u and v perturbed 50 times. Thereby, Fig. 2 represents the effectiveness of the block diagonalization using H (u, v) and Jn B(u, v)Jn of size 2448.
254
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
Fig. 2. ErrorAHankel (left figure) and ErrorABe´ zout (right figure) for n = 2448.
Table 1 The relative error in floating arithmetic with 15 digits. n
ErrorRHankel
ErrorRBe´ zout
1024 1224 1424 1624 1824 2048 2448
5.784074308862707e−012 7.377238376056904e−007 1.229383772472314e−007 1.402838670935590e−006 1.552356261630665e−007 1.482326667465912e−005 3.863794639419495e−005
1.086190389822859e−008 1.503623817980172e−007 5.902898664239044e−012 4.631154488009333e−007 8.760591568775624e−008 1.042891806603727e−009 1.162851904025527e−008
Fig. 3. The relative error (∗: Hankel, +:Bézout) for n = 1424.
5.5. Stability Now, we compare the approach based on H (u, v) and Jn B(u, v)Jn for their numerical stability. The main objective of this part is to improve the numerical stability of the approximate block diagonalization using Bézout matrices because the Bézout matrix coefficients are bounded by the Bézout Φ = 2nµν (with ∥u∥∞ ≤ µ and ∥v∥∞ ≤ ν ) and the Hankel matrix coefficients H (u, v) may increase exponentially. Also, if u and v are integer coefficients, the intermediate calculations steps can be preserved in integers numbers. See [6] for more details. Thus, we give the relative error (15), in floating arithmetic with 15 digits for H (u, v) and Jn B(u, v)Jn associated with two ‘‘exact’’ polynomials u and v . Thus, Table 1 shows the effectiveness of Algorithm 4.3 in terms of stability, specially when n increases. Moreover, Fig. 3 represents the stability behavior for the coefficients of 1424 × 1424 Hankel and Bézout matrices constructed with two polynomials.
S. Belhaj / Journal of Computational and Applied Mathematics 250 (2013) 244–255
255
6. Concluding remarks In this paper, we introduced a fast algorithm for computing the polynomials generated by the Euclidean Algorithm via an approximate block LU factorization of a Bézout matrix associated with the input polynomials. Numerical experiments performed with a wide variety of test problems, showed the effectiveness of this algorithm in terms of efficiency, stability and robustness. Acknowledgments I thank Prof. D. A. Bini for suggesting me this problem and guiding me throughout the course of this investigation and I owe a special debt of gratitude to Prof. H. Lombardi and Prof. G. M. Diaz-Toca for the helpful discussions. I also thank the anonymous referee for those helpful comments. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]
G.E. Collins, Subresultants and reduced polynomial remainder sequences, J. Assoc. Comput. Mach. 14 (1967) 128–142. W.S. Brown, On Euclid’s Algorithm and the computation of polynomial greatest common divisors, J. Assoc. Comput. Mach. 18 (1971) 478–504. W.S. Brown, J.F. Traub, On Euclid’s Algorithm and the theory of subresultants, J. Assoc. Comput. Mach. 18 (1971) 505–514. A. Schönhage, Quasi-GCD Computations, J. Complexity 1 (1) (1985) 118–137. M.-T. Noda, T. Sasaki, Approximate GCD and its application to ill-conditioned algebraic equations, J. Comput. Appl. Math. 38 (1–3) (1991) 335–351. D. Bini, V. Pan, Polynomial and Matrix Computations, Vol. 1, in: Fundamental Algorithms, Birkhäuser, 1994. D. Bini, L. Gemignani, Fast parallel computation of the polynomial remainder sequence via Bézout and Hankel matrices, SIAM J. Comput. 24 (1995) 63–77. L. Gemignani, GCD of polynomials and Bézout matrices, in: Proceedings of the 1997 International Symposium on Symbolic and Algebraic Computation (Kihei, Maui, HI), ACM Press, New York, 1997, pp. 271–277. D.A. Bini, L. Gemignani, Fast fraction-free triangularization of Bezoutians with applications to subresultant chain computation, Linear Algebra Appl. 284 (1–3) (1998) 19–39. W.B. Gragg, A. Lindquist, On partial realization problem, Linear Algebra Appl. 50 (1983) 277–319. G. Heining, K. Rost, Algebraic Methods for Toeplitz-Like Matrices and Operators, Birkhäuser Verlag, Basel, 1984. G.M. Diaz-Toca, N. Ben Atti, Block LU factorization of Hankel and Bézout matrices and Euclidean algorithm, Int. J. Comput. Math. 86 (2009) 135–149. N. Ben Atti, G.M. Diaz-Toca, Block diagonalization and LU-equivalence of Hankel matrices, Linear Algebra Appl. 412 (2006) 247–269. M. Mitrouli, N. Karcanias, C. Koukouvinos, Numerical performance of the matrix pencil algorithm computing the greatest common divisor of polynomials and comparison with other matrix-based methodologies, J. Comput. Appl. Math. 76 (1–2) (1996) 89–112. N. Karmarkar, Y.N. Lakshman, Approximate polynomial greatest common divisor and nearest singular polynomials, in: Proceedings of the International Symposium on Symbolic and Algebraic Computation, ACM Press, New York, 1996, pp. 35–39. R.M. Corless, S.M. Watt, L. Zhi, QR factoring to compute the GCD of univariate inexact polynomials, IEEE Trans. Signal Process. 52 (12) (2004) 3394–3402. C.J. Zarowski, X. Ma, F.W. Fairman, QR-factorization method for computing the greatest common divisor of polynomials with inexact coefficients, IEEE Trans. Signal Process. 48 (11) (2000) 3042–3051. V.Y. Pan, Computation of approximate polynomial GCDs and an extension, Inform. and Comput. 167 (2) (2001) 71–85. Z. Zeng, The approximate GCD of inexact polynomials, Part 1: a univariate algorithm, 2004, Preprint. B. Li, Z. Liu, L. Zhi, A structured rank-revealing method for Sylvester matrix, J. Comput. Appl. Math. 213 (1) (2008) 212–223. J.R. Winkler, J.D. Allan, Structured total least norm and approximate GCDs of inexact polynomials, J. Comput. Appl. Math. 215 (1) (2008) 1–13. J.R. Winkler, X. Lao, The calculation of the degree of an approximate greatest common divisor of two polynomials, J. Comput. Appl. Math. 235 (6) (2011) 1587–1603. S. Belhaj, A fast method to block-diagonalize a Hankel matrix, Numer. Algorithms 47 (2008) 15–34. S. Belhaj, Computing the block factorization of complex Hankel matrices, Computing 87 (3–4) (2010) 169–186. S. Belhaj, Block factorization of a Hankel matrix and Euclidean Algorithm, Math. Model. Nat. Phenom. 5 (7) (2010) 38–44. P. Lancaster, M. Tismenetsky, The Theory of Matrices, second ed., in: Computer Science and Applied Mathematics, Academic Press Inc., Orlando, FL, 1985. Fu-Rong Lin, Wai-Ki Ching, M.K. Ng, Fast inversion of triangular Toeplitz matrices, Theoret. Comput. Sci. 315 (2004) 511–523. V. Hribernig, Sensitivity of algebraic algorithms, Ph.D. Thesis at the Technical University of Vienna, 1994. V. Hribernig, H.J. Stetter, Detection and validation of clusters of polynomial zeros, J. Symbolic Comput. 24 (6) (1997) 667–681.