Matrix representation of the shifting operation and numerical properties of the ERES method for computing the greatest common divisor of sets of many polynomials

Journal of Computational and Applied Mathematics 260 (2014) 54–67

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Matrix representation of the shifting operation and numerical properties of the ERES method for computing the greatest common divisor of sets of many polynomials D. Christou a,∗ , N. Karcanias a , M. Mitrouli b a

Systems and Control Centre, School of Engineering and Mathematical Sciences, City University, Northampton Square, EC1V 0HB, London, United Kingdom b

Department of Mathematics, University of Athens, Panepistemiopolis 15773, Athens, Greece

article

info

Article history: Received 14 May 2012 Received in revised form 30 June 2013 Keywords: Univariate real polynomials Greatest common divisor Shifting operation Numerical stability

abstract The Extended-Row-Equivalence and Shifting (ERES) method is a matrix-based method developed for the computation of the greatest common divisor (GCD) of sets of many polynomials. In this paper we present the formulation of the shifting operation as a matrix product which allows us to study the fundamental theoretical and numerical properties of the ERES method by introducing its complete algebraic representation. Then, we analyse in depth its overall numerical stability in finite precision arithmetic. Numerical examples and comparison with other methods are also presented. © 2013 Elsevier B.V. All rights reserved.

1. Introduction The computation of the greatest common divisor (GCD) of polynomials is an algebraic problem which has been intensely studied for many years. The algorithm associated with Euclid’s division method [1] is the oldest known solution to this problem. The computational methods for computing the GCD of real univariate polynomials can be separated in two main categories: (a) The Euclidean type methods which rely on Euclid’s division algorithm and its variations. (b) The matrix-based methods which are based on the processing of a matrix formed directly from the coefficients of the given polynomials. The matrix-based methods may be further classified to those that: (i) form a matrix for two polynomials and work on pairwise computations iteratively, (ii) form a matrix that corresponds to the whole set of polynomials and process it either directly or iteratively. Early GCD algorithms were developed using Euclidean-based methods, applied to two polynomials [2–4]. The Euclidean algorithm is efficient when the polynomials have integer coefficients, but it becomes inefficient when the polynomials have coefficients from the field of real numbers due to the use of finite precision arithmetic, which introduces numerical errors in the solution. In 1985, Schönhage introduced the notion of Quasi-GCD [5], and in 1989, Noda and Sasaki described a

∗

Corresponding author. E-mail addresses: [email protected], [email protected] (D. Christou), [email protected] (N. Karcanias), [email protected] (M. Mitrouli). 0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.09.021

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67

55

special version of Euclid’s algorithm for computing the GCD of a pair of coprime polynomials with inexact coefficients [6]. This approach is amongst the first attempts to define and compute an approximate GCD of polynomials by means of symbolic–numeric computations. The development of numerical stable GCD algorithms which can deal with polynomials of inexact data has attracted a lot of attention the past thirty years and several methods have been proposed [7–16]. The various techniques, which have been developed for the computation of approximate solutions, are based on methodologies where exact properties of these notions are relaxed and appropriate solutions are sought by using a variety of numerical tests. The use of finite precision arithmetic in computer algebra makes the extension of the Euclidean algorithm to sets of many polynomials a rather difficult task. The iterative application of the Euclidean algorithm to two polynomials at a time often results in a total numerical error which might exceed the machine’s fixed numerical tolerance. Conversely, the developed matrix-based methods tend to be more effective in handling sets of several polynomials and produce solutions of better numerical quality. The use of matrices in the problem of computing the GCD of many polynomials appears early in Barnett’s work [17–19], who developed a technique for computing the degree and the coefficients of the GCD by using companion and Sylvester matrices. In [20], Karcanias (1987) has shown that the GCD is an invariant of the row space of the basis matrix of the set of polynomials, as well as that it is also an invariant under the shifting operation. This has led to the development of the current approach that avoids the use of Euclidean divisions. The Extended-Row-Equivalence and Shifting (ERES) method [21] is an iterative matrix-based method developed for the computation of the greatest common divisor (GCD) of sets of real polynomials in one variable. The method exploits the invariance of the greatest common divisor of a set of many polynomials under elementary row transformations and partial column shifting by transforming a basis matrix, which is formed directly from the coefficients of the polynomials of a given set, into a simpler matrix containing the vector of coefficients of the GCD. The previous study of its theoretical and numerical properties in [20–22] revealed the advantage of the ERES method to handle large sets polynomials and to invoke an efficient termination criterion that allows the computation of approximate solutions when the initial data have numerical inaccuracies [23]. The development of the ERES method as it has been described in [21] is an inherently robust algebraic method which defines a special matrix equivalence. However, the complete algebraic representation of the method remained an open issue due to its iterative nature and the luck of an algebraic expression for the partial column shifting transformation, referred to as the shifting operation. The aim is to establish an algebraic relationship between the initial basis matrix of a given set of several polynomials and the last matrix which occurs after the iterative application of the ERES operations and provides the GCD. Apart from its theoretical value, this algebraic representation is significant for the complete analysis of the overall numerical stability of the ERES method which could not be studied before. The main objectives of this paper are: (a) to present the general algebraic representation of the ERES method and discuss its theoretical and practical use, and (b) to analyse the overall numerical stability of the method in finite-precision arithmetic. The paper is structured as follows: In Section 2, the definition and the most important properties of the ERES operations are presented and we give a brief description of how the ERES method is formulated. In Section 3 the major issue of having a matrix representation for the shifting operation is analysed and discussed. The results from the study of the shifting operation are used in Section 4 in order to introduce the general algebraic representation of the ERES method which eventually establishes the ERES representation of the GCD of a set of many polynomials. This representation forms the basis for the analysis of the overall numerical stability of the method in Section 5. Numerical examples and comparison with other methods are also presented. Notation. In the following, N and R denote the sets (fields) of natural and real numbers, respectively. R[s] denotes the ring of polynomials in one variable over R. Capital letters denote matrices and small underlined letters denote vectors. By p(s) we denote a polynomial in s with real coefficients. The greatest common divisor of a set of polynomials P will be denoted by gcd{P }. The following list includes the basic notations that are used in the document. A ∈ Rm×n

v ∈ Rm At

vt ρ(A) deg{p(s)} ∥v∥2 ∥A∥2 ∥A∥∞ ,

:= ≈

Matrix A with elements from R arranged in m rows and n columns (m, n ∈ N and m, n ≥ 2). Column vector with m ≥ 2 elements from R. Transpose matrix of A. Transpose vector of v . The rank of a matrix A. The degree of a polynomial p(s). 

µ |vi |2 . √ i=1 The Euclidean norm of A : ∥A∥2 = max eigenvalue of At A. ν The infinity norm of A : ∥A∥∞ = max1≤i≤µ j=1 |aij |. The Euclidean norm of v : ∥v∥2 =

Mathematical operator which denotes equality by definition. Mathematical operator which denotes equality by input. Mathematical operator which denotes approximate equality.

56

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67

2. Definition of the ERES operations and background results 2.1. Background theory Consider the set of real polynomials in one variable (univariate polynomials):





Pm,n = pi (s) ∈ R[s], i = 1, 2, . . . , m with n = max (deg{pi (s)} ≥ 1) .

(1)

i

We represent the polynomials pi (s) with respect to the highest degree n as pi (s) = ai,n sn + ai,n−1 si,n−1 + · · · + ai,1 s + ai,0 ,

ai,n ̸= 0, i = 1, 2, . . . , m.

(2)

Definition 1. For any Pm,n set, we define a vector representative (vr), p (s) and a basis matrix Pm represented as m

, pm ] · en (s) = Pm · en (s), m−1

p (s) = [p1 (s), . . . , pm (s)]t = [p , . . . , p m

1

where Pm ∈ Rm×(n+1) , en (s) = [1, s, . . . , sn−1 , sn ]t and p ∈ Rn+1 for all i = 1, . . . , m. i

The matrix Pm is formed directly from the coefficients of the polynomials of the set Pm,n and it has the least possible dimensions.

= 0 and pc ̸= 0, then c = w(Pm,n ) is called the order of Pm,n and s is an elementary divisor of the GCD. The set Pm,n is considered to be a c-order set and will be called proper if c = 0, and nonproper if c ≥ 1. Definition 2. If c is the integer for which p = · · · = p c

1

c −1

For a nonproper set Pm,n with w(Pm,n ) = c, we can always consider the corresponding proper one Pm,n−c by dismissing the c leading zero columns. Then gcd{Pm,n } = sc · gcd{Pm,n−c }. In the following without loss of generality we assume that Pm,n is proper. Definition 3 (ERES Operations). Given a set Pm,n of many polynomials with a basis matrix Pm the following operations are defined [20]: (a) Elementary row operations with scalars from R on Pm . (b) Addition or elimination of zero rows on Pm . (c) If at = [0, . . . , 0, al , . . . , an+1 ] ∈ Rn+1 , al ̸= 0 is a row of Pm then we define as the shifting operation shf : shf(at ) = [al , . . . , an+1 , 0, . . . , 0] ∈ Rn+1 . By shf(Pm,n ), we shall denote the set obtained from Pm,n by applying shifting to every row of Pm (matrix shifting). Type (a)–(c) operations are referred to as Extended-Row-Equivalence and Shifting (ERES) operations. The ERES operations without applying the shifting operation are referred to as ERE operations. The following theorem describes the properties characterizing the GCD of any given Pm,n . Theorem 1 ([20]). For any set Pm,n , with a basis matrix Pm , ρ(Pm ) = r and gcd{Pm,n } = φ(s) we have the following properties: (i) If RP is the row space of Pm , then φ(s) is an invariant of RP (e.g. φ(s) remains invariant after the execution of elementary row operations on Pm ). Furthermore if r = dim(RP ) = n + 1, then φ(s) = 1. (ii) If w(Pm,n ) = c ≥ 1 and Pm∗,n = shf(Pm,n ), then

  φ(s) = gcd{Pm,n } = sc · gcd Pm∗,n . (iii) If Pm,n is proper, then φ(s) is invariant under the combined ERES set of operations. 2.2. The formulation of the ERES method and the computation of the GCD of a set of polynomials ERES operations preserve the GCD of any Pm,n and thus can be easily applied iteratively in order to obtain a modified basis matrix with much simpler structure [20]. The ERES method in its simplest form consists of three basic procedures: 1. Construction of the proper basis matrix Pm for the set Pm,n . 2. Application of elementary row operations to the processed matrix, which practically involves row reordering, triangularization, and elimination of zero rows (ERE operations). 3. Application of the shifting operation to the nonzero rows of the processed matrix.

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67

57

After successive applications of the ERES operations to the initial basis matrix, the maximal degree of the given set of polynomials is reduced, and after a finite number of steps the resulting matrix has rank 1. At this point, the process is terminated and, considering that all the arithmetic operations are performed accurately (e.g. by using symbolic–rational operations), any row of the last matrix specifies the coefficients of the GCD of the set. The iterative application of the processes of triangularization and shifting forms the core of the ERES method and we shall refer to it as the main procedure of the method. The main problem in the formulation of an algebraic expression, which will establish the connection between the initial basis matrix Pm and the final matrix that contains the coefficients of the GCD, requires appropriate matrix representations of the ERE operations and the shifting operation, respectively. The ERE row operations, i.e. triangularization, deletion of zero rows and reordering of rows, can be represented [24,25] by a matrix R ∈ Rr ×m , r < m, which converts the initial rectangular matrix Pm into an upper trapezoidal form. Conversely, the matrix representation of the shifting operation is not straightforward. The problem of the matrix representation of the shifting operation for real matrices has remained open until now. Solving this problem is crucial for establishing a general matrix representation of the ERES method for all the performed iterations. Therefore, in the following section we aim to find the simplest possible algebraic relation between a matrix and its shifted form. 3. The shifting operation for real matrices The shifting operation is a special matrix transformation which is not very common in the literature of algebra. In Definition 3 the shifting operation is defined for real vectors as a permutation of consecutive zero elements. Specifically, having a real vector a = [0, . . . , 0, ak+1 , . . . , an ] ∈ Rn ,

ak+1 ̸= 0

   k elements

the shifting operation is defined as shf : shf(a) = [ak+1 , . . . , an , 0, . . . , 0] ∈ Rn . This definition can be extended to the case of real matrices. Definition 4. Given a matrix A = [a1 , a2 , . . . , am ]t ∈ Rm×n , where ai ∈ Rn for i = 1, 2, . . . , m are the row-vectors of A, the shifting operation for matrices is defined as the application of vector-shifting to every row of A. This transformation will be referred to as matrix-shifting and the shifted form of A will be denoted by shf(A) , A∗ = [shf(a1 ), shf(a2 ), . . . , shf(am )]t ∈ Rm×n . It is important to notice that the shifting operation, as defined here, permutes the elements of a vector without changing their values, and this is a basic requirement for the shifting operation in the study of the numerical properties of the ERES method. The vector-shifting can be represented by the multiplication: shf(a) = a · Jk,n where Jk,n is an appropriate n × n permutation matrix which is a square binary matrix that has exactly one entry 1 in each row and each column, and zeros elsewhere. Each such matrix represents a specific permutation of k elements and for the vector-shifting it has the form:

 Jk,n =

On−k In − k

Ik Ok



∈ Rn×n

(3)

where Ii denotes the i × i identity matrix and Oi denotes the i × i zero matrix for i = k, n − k. Although it is rather simple to represent the vector-shifting transformation with a simple vector–matrix multiplication, it is not obvious how to represent the matrix-shifting transformation, because in general the application of vector-shifting to the rows of a matrix alters the structure of the columns in a non-uniform way. The problem of representing the matrixshifting by using an appropriate matrix–matrix multiplication is challenging for the study of the theoretical and numerical properties of the ERES method. We are particularly interested in finding an algebraic relationship between a real matrix and its shifted form in the class of upper trapezoidal matrices. This type of matrices occurs after applying Gaussian elimination, or other triangularization method, and they have the following generic form: a11 0



A=  ..

.

0

a12 a22

..

. ···

··· ··· .. .

a1m a2m

.. .

··· ··· .. .

a1n a2n 

0

amm

···

amn



m×n ..  ∈R , .

m < n.

(4)

58

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67

Then, the shifted form of A, which is obtained by the matrix-shifting transformation as defined in Definition 4, is a11  a22



A∗ =   .. .

am m

a12

···

a1m

···

··· .. .

a2m

···

a2n . ..

···

am n

.. .

..

. 0

···

Definition 5. (a) If A = a1 , a2 , . . . , am every i = 1, 2, . . . , m, such that



A=

m 

t

a1n 0 



m×n ..  ∈R , .

m < n.

(5)

0

 t ∈ Rm×n , then we can define the matrices Ai = 0, . . . , ai , . . . , 0 ∈ Rm×n , for

Ai

(6)

i=1

where ai is the ith row of A and 0 ∈ Rn is a n-dimensional zero vector. (b) We define the permutation matrices Ji ∈ Rn×n for i = 1, 2, . . . , m, so that every Ji gives the appropriate shifting to each Ai respectively. Therefore, shf(A) =

m 

Ai Ji .

(7)

i=1

Since a11 ̸= 0, we note that J1 = In , where In is the n × n identity matrix. Remark 1. If A has full rank, then, since it is defined as an upper trapezoidal matrix with aii ̸= 0 for all i = 1, . . . , m, it is 1 1 right-invertible. Let us denote its right inverse by A− ∈ Rn×m . Hence, A A− = Im . r r The following theorem establishes the connection between a nonsingular upper trapezoidal matrix and its shifted form. Theorem 2. If A ∈ Rm×n , 2 ≤ m < n, is a non-singular upper trapezoidal matrix with rank ρ(A) = m and shf(A) ∈ Rm×n is the matrix obtained from A by applying shifting to its rows, then there exists a square matrix S ∈ Rn×n such that: shf(A) = A · S .

(8)

The matrix S will be referred to as the shifting matrix of A. Proof. Let A∗ = shf(A). We shall use the notation described in Definition 5 and we will follow the next method to determine the shifting matrix S ∈ Rn×n . 1. Apply to the original matrix A the n × n(m + 1) block matrix: S (1) = J1



···

1 A− r

Jm



(9)

such that: A(1) = A · S (1) . 2. Multiply the matrix A(1) by the n(m + 1) × mn block matrix: In On

On In

.. .

··· ··· .. .

On On

On (A1 − A) J1

On (A2 − A) J2

··· ···

In

  

S (2) =  

.. .



.. .

     

(10)

(Am − A) Jm

where On denotes the n × n zero matrix and In the n × n identity matrix. Hence, A(2) = A(1) · S (2) . 3. Multiply the matrix A(2) by the mn × n block matrix:

  In

S

(3)

  =  ...  In

and hence, A(3) , A∗ = A(2) · S (3) .

(11)

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67

59

The final matrix S = S (1) · S (2) · S (3) has the form: S=

m  

1 −1 In − A− r A + A r A i Ji



(12)

i=1

and satisfies the equation: A∗ = A · S.



1 In the proof of Theorem 2 the right inverse matrix A− of A is not unique when m < n. Conversely, the pseudo-inverse r matrix AĎ ∈ Rn×m of A can be uniquely determined by calculating the singular value decomposition of A [25], such that

A A Ď = Im . Therefore, an alternative expression of the representation (12) of the shifting matrix S can be given, if we use the pseudoinverse matrix of A. This is S=

m  

In − AĎ A + AĎ Ai Ji



(13)

i=1

and it is more appropriate for the numerical computation of the shifting matrix. Example 1. Consider the following randomly selected matrix A and its shifted form A∗ : 1 0 0

2 1 0

 A=

3 2 2

4 3 4

5 0 , 6



1 1 2

2 2 4

 ∗

A =

3 3 6

4 0 0

5 0 . 0



According to (13), the corresponding shifting matrix is:



0

  2 −   3  1 S=  3   1   3 

0

− −

36 43 25 43 9

36

80

67 

43 161

43 319

129 113

258 17

43  487  

− −

86 34

258 145

43 9

129 9

86

86

− −

258 95 258 23

258  25  



258  179 



−

86 ∗

and it can be easily verified that A = A · S.



258   37 86



Obviously, the shifting operation alters the structure of a matrix. Therefore, even if the original matrix has full rank, the corresponding shifted matrix may not have full rank. For instance, in the previous Example 1 we have ρ(A) = 3 and ρ(A∗ ) = 2. However, in the case where both A and A∗ have full rank we obtain the following result. Corollary 1. If A ∈ Rm×n , m < n, is a nonsingular upper trapezoidal matrix with rank ρ(A) = m and A∗ ∈ Rm×n is the shifted matrix of A with rank ρ(A∗ ) = m, then there exists an invertible matrix S ∈ Rn×n with rank ρ(S ) = n, such that A ∗ = A · S ⇔ A = A ∗ · S −1

(14)

where S −1 denotes the inverse of S. The previous corollary can be proven by following the same steps as in the proof of Theorem 2. We only have to change appropriately the set of permutation matrices Ji , i = 1, 2, . . . , m to achieve the proper shifting, and compute the inverse or pseudo-inverse of A∗ . Therefore, we conclude that the shifting of a nonsingular upper trapezoidal matrix is a reversible process, unless the shifted matrix is rank deficient. Remark 2. The results from Theorem 2 and Corollary 1 can also be applied to a square upper triangular matrix provided that this matrix is invertible. The relations (8) and (14) have a key role in the general algebraic representation of the ERES method, since (a) in every iteration of the main procedure there is always a nonsingular upper trapezoidal matrix which is formed after the application of the ERE operations, and (b) the process stops when a rank-deficient matrix occurs and shifting cannot be applied.

60

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67

4. The general algebraic representation of the ERES method ERES is an iterative matrix-based method where, until now, only the ERE operations (i.e. triangularization and reordering of rows) could be represented by a matrix R ∈ Rm×m . With the introduction of the representation of the shifting operation as a matrix product in Theorem 2 it is now possible to form an algebraic expression representing all the required transformations until a rank-1 matrix is reached. This representation provides a link between the initial basis matrix and the last matrix which gives the coefficients of the GCD. This algebraic relationship is described by the following theorem. Theorem 3. Given a set Pm,n of m real univariate polynomials of maximum degree n ∈ N and its basis matrix Pm ∈ Rm×(n+1) , the application of the ERES operations to Pm results in a matrix G ∈ Rm×(n+1) with rank ρ(G) = 1, which satisfies the equation: G = R · Pm · S

(15) (n+1)×(n+1)

m×m

where R ∈ R and S ∈ R represent the applied row transformations (ERE operations) and the application of the shifting operation, respectively. The GCD of Pm,n is then represented by gcd{Pm,n } = e1 · G · en (s)

(16)

where e1 = [1, 0, . . . , 0] ∈ R and en (s) = [1, s, s , . . . , s ] . m

2

n t

Proof. Given a set Pm,n of m > 2 polynomials and its basis matrix Pm , let P (1) := Pm be the initial matrix and P (k) is the processed matrix at the beginning of the kth iteration, k ∈ N. The superscript ‘‘(k)’’, k = 1, 2, 3, . . . , will be used in all matrices to indicate the number of iteration of the main procedure. The ERES operations are performed in the following order: 1. Construct the permutation matrix J (1) ∈ Rm×m which reorders the rows of the initial matrix such that the first row corresponds to the polynomial with the lowest degree in the set. 2. Apply the elementary row transformations (ERE operations) by using an appropriate lower triangular matrix L(1) ∈ Rm×m . 3. Delete (or reorder) the zero rows by using an appropriate permutation matrix Z (1) ∈ Rr1 ×m , r1 ≤ m. 4. Apply the shifting operation by using an appropriate square matrix S (1) ∈ R(n+1)×(n+1) . Therefore, after performing the above transformations, the resulting matrix is P (2) = Z (1) · L(1) · J (1) · P (1) · S (1) . If we set R

(1)

=Z

(1)

·L

(1)

·J

(1)

(17)

, then it follows

P (2) = R(1) · P (1) · S (1) .

(18)

Eq. (18) represents the first complete iteration of the main procedure of the ERES method. The whole process terminates when a matrix with rank equal to 1 appears. This can be practically achieved in less than n + 1 iterations. Therefore, after the kth iteration we have P (k+1) = R(k) · P (k) · S (k) ,

k = 1, 2, . . .

(19)

and, if the final number of iterations is ℓ ∈ N, then P (ℓ+1) = R(ℓ) · · · R(1) · Pm · S (1) · · · S (ℓ) ⇔ P (ℓ+1) = R˜ · Pm · S

(20)

where we denote by R˜ =

ℓ 

R(k)

and

k=1

S=

ℓ 

S (k) .

(21)

k=1

Obviously, the matrices P (k) do not necessarily have the same dimensions as the initial matrix Pm due to the frequent deletion of the produced zero rows during the iterative main procedure of the method. However, for theoretical purposes we may preserve the original dimensions of the basis matrix. This can be easily achieved if we change the permutation matrix Z (k) so as to move the zero rows of P (k) to the bottom of the matrix instead of deleting them. Therefore, in this case P (k) can have the same dimensions as Pm , but we actually continue to work with a rk × (n + 1) submatrix of P (k) with a decreasing row dimension rk < m. Since the last matrix P (ℓ+1) has rank equal to 1, every row gives the coefficients of the GCD. But if we triangularize P (ℓ+1) once more by using an appropriate matrix L(ℓ+1) , then the final matrix G ∈ Rm×(n+1) contains the coefficients of the GCD in its first row and it has zeros elsewhere. Hence, G = L(ℓ+1) · P (ℓ+1)

and R = L(ℓ+1) · R˜

(22)

and the main result in (15) derives from the combination of (20)–(22). Then, if e1 = [1, 0, . . . , 0] ∈ R [1, s, s2 , . . . , sn ]t , the GCD is given by

m

gcd{Pm,n } = e1 · R · Pm · S · en (s) = e1 · G · en (s).



and en (s) = (23)

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67

61

The next example demonstrates the application of the ERES method to a set of three polynomials. Example 2. Consider the set of polynomials:

  p1 (s) = s3 − 2s2 − 11s + 12  = p2 (s) = 2s3 − 11s2 + 13s − 4 ,   p3 (s) = s2 − s − 12

Pm,n

m = 3, n = 3

with gcd{Pm,n } = s − 4. The initial basis matrix is:

 Pm =

−11

12 −4 −12

13

−2 −11

−1

1

1 2 0

1

 

 3 ×4

∈R

,

s en (s) =  2  . s s

(24)

3

The iterative main procedure of the ERES method will start with the matrix Pm . After two iterations of the main procedure, the final matrix will have rank equal to 1 and its first row gives the vector of coefficients of the GCD. The matrix R, which represents all the necessary row transformations, and the matrix S, which represents all the shifting transformations, have the form:

 529687 5

9

1

 14  10 −

28 9

4

 R=

−

7 1



  0

7 0

0

 397579  528432   397579 and S =   523412   397579 

− −

120713

551

69 

397579 85273

23387 2204

397579 56487

23387 8816

4369  276  

503332

397579 623527

23387 35264

397579

397579

23387

4369  . 1104  



(25)

4369   4416

4369

The computation of the shifting matrices during the iterations corresponds to (13). The final matrix is

 −4 G=

0 0

1 0 0

0 0 0

0 0 0

 ∈ R3×4

(26)

and it can be verified that: G = R · Pm · S . If en (s) = [1, s, s2 , s3 ]t and e1 = [1, 0, 0, 0], then the GCD of the set Pm,n can be expressed as: gcd{Pm,n } = e1 · R · Pm · S · en (s) = s − 4. 

(27)

5. Analysis of the numerical stability of the ERES method Based on the relation (15) we can now develop a detailed analysis of the numerical stability of the method for all the performed iterations of its main procedure. 5.1. The basics of the ERES algorithm The development of an effective numerical algorithm for the ERES method, requires [21]: (a) to find a robust numerical procedure for the application of the ERE operations, (b) to develop a proper termination criterion for the algorithm, and finally (c) to find a reliable way to extract the coefficients of the GCD from the last rank-1 matrix. These requirements are the most essential parts of the ERES algorithm (Fig. 1) and in the following we will consider them in the context of a numerical implementation in a floating-point computational environment. Having a basis matrix of a set of polynomials, the ERES algorithm involves row addition or row multiplication, row reordering, elimination of zero rows and shifting. We refer to this process as the main procedure of the ERES algorithm. The most reliable and stable numerical method for applying elementary row operations is the Gaussian elimination with partial pivoting (GEPP) [24,25]. Hence, an upper triangular or trapezoidal form of the basis matrix is computed. The shifting operation is merely a permutation of the leading consecutive zero elements in a row. The algorithm’s termination criterion relies on the proper detection of the final unity rank matrix which is based on the numerical computation of the singular values of an associated normalized matrix obtained at the end of each iteration of

62

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67

Fig. 1. The ERES algorithm.

the main procedure. We shall refer to it as the Rank-1 procedure of the ERES algorithm. This property is detected numerically according to the inequality [21]:

|σ1 −

√ µ| ≤ εt and σi ≤ εt , i = 2, 3, . . . , µ

(28)

where σi are the singular values of a normalized matrix with µ < m rows, which is formed at the end of the main procedure, and σ1 corresponds to the maximum singular value. The tolerance εt is referred to as the termination accuracy of the ERES algorithm in finite precision computations. When a matrix P with εt -rank = 1 is achieved, the vector of the coefficients of the GCD is given by the first row of the matrix: G1 = σ1 · u · w t

(29)

where u and w are the first columns of the orthogonal matrices U and W of the singular value decomposition P = U · Σ · W t , respectively [21]. 5.2. Estimation of the total numerical error of the ERES method (k)

Let us denote by P (k) = [pij ] ∈ Rrk ×(n+1) the matrix to be processed at the kth iteration of the main procedure of the ERES algorithm with rk ≤ m and k ∈ N. The superscript ‘‘(k)’’ will be used in all matrices to indicate the number of iteration of the main procedure, and P (1) := Pm . If we denote by d > 0 the degree of the GCD, the constant application of the shifting operation gradually zeros the last n − d columns of Pm , and therefore we will denote by nk the number of the first consecutive nonzero columns of P (k) . Clearly, nk ≤ n + 1, and hence we may consider P (k) as an rk × nk matrix. In the following, we shall consider the numerical error for each individual step of the main procedure of the ERES algorithm and we will conclude with the total numerical error of the method. We measure this error by using the infinity matrix norm ∥ · ∥∞ which is invariable under the row-reordering and shifting transformations. Starting the main procedure, we must reorder the rows of P (k) such that the first row corresponds to the least degree polynomial. This is a basic step of the ERES method, but, in order to prevent a further change during the process of GEPP, it (k) is required to scale the elements of the matrix, such that the first element p11 of the first column can be larger in magnitude than the other elements of the first column. This normalization can be achieved if we multiply P (k) by an appropriate diagonal matrix N1 ∈ Rrk ×rk with ∥N1 ∥∞ = 1, [22], such that P˜ (k) = N1 P (k) + E1

(30)

where E1 is the error matrix and

∥E1 ∥∞ ≤ nk u ∥P (k) ∥∞

(31) 1 −t b 2

−53

=2 in binary 64 bits arithmetic, known as ‘‘double precision’’). However, where u is the machine precision (u = this normalization is not always necessary to be performed and in practice this error is negligible compared to the error produced from the GEPP process. A further normalization of P (k) can be added so that every element of P (k) is bounded by unity. This normalization can be written as [24]: (k)  P (k) = N (k) P (k) + EN

where N N

(k)

(k)

(32)

is a diagonal rk × rk matrix of the form:

  (k) −1 (k) −1 1 = diag ∥p(1jk) ∥− 2 , ∥p2j ∥2 , . . . , ∥prj ∥2

for j = 1, 2, . . . , nk

(k)

and EN is the error matrix with

∥EN(k) ∥∞ ≤ nk u ∥P (k) ∥∞ . Consequently, ∥ P (k) ∥∞ ≤ nk .

(33)

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67

63

The backward error analysis of the Gaussian elimination with partial pivoting [26] shows that the computed upper and lower triangular matrices L(k) and U (k) satisfy: (k)

L(k) · U (k) =  P (k) + EG (k)

n2k

∥E G ∥∞ ≤

(34)

(k)

ρk u ∥P ∥∞

(35)

(k)

where EG is the error matrix and generally, ∥ P (k) ∥∞ ≤ ∥P (k) ∥∞ . The term ρk denotes the growth factor [24], which is n k −1 ρk ≤ 2 for GEPP, but in practice this bound is not attainable [26]. The upper triangular matrix U (k) will eventually give the next matrix P (k+1) after the deletion of its zero rows and shifting: −1 (k)  U (k) = L(k) P (k) + EG





(36)

and P (k+1) = Z (k) · U (k) · S (k) .

(37)

Therefore, at the end of the kth iteration of the main procedure we have: −1 P (k+1) = Z (k) · L(k)



(k)

N (k) N1 · J (k) · P (k) + EN







 + EG(k) S (k) .

(38)

Then, we may set −1 R(k) = Z (k) · L(k) · N (k) · N1 · J (k)

(39)

and the error matrix for the kth iteration is −1 (k) (k) E (k) = Z (k) · L(k) EN + EG S (k) .





(40)

The reordering of rows, the deletion of zero rows, and the shifting are error-free transformations, since they do not alter the values of the data. Especially for the shifting operation, there is no need to compute the shifting matrices S (k) . We only use them in order to connect the matrices P (k) which are generated in every iteration of the main procedure. Thus, for practical reasons we may set ∥S (k) ∥∞ = 1. According to the form of the matrices in (40) it is ∥Z (k) ∥∞ = 1, and for normalized matrices we have ∥L(k) ∥∞ ≤ nk [26]. Hence, if we combine the relations (33), (35) and (40) we conclude with the result in the next lemma, which describes the numerical error E (k) produced in every iteration of the main procedure of the ERES algorithm. −1

Lemma 1. The matrix P (k+1) , k ∈ N, which is produced after the numerical processing of the matrix P (k) ∈ Rrk ×nk during the main procedure of the ERES algorithm, satisfies the equation: P (k+1) = R(k) · P (k) · S (k) + E (k)

(41)

  ∥E (k) ∥∞ ≤ n3k ρk + n2k u ∥P (k) ∥∞

(42)

with

where R(k) denotes the matrix for the combined ERE operations and S (k) denotes the matrix for the shifting transformation during the kth iteration of the main procedure. If we denote by ℓ the total number of iterations of the main procedure of the ERES algorithm, then the total numerical error E for all the performed iterations is E=

ℓ 

E (k)

(43)

k=1

with

 ∥E ∥∞ ≤

ℓ  

n3k

ρk +

n2k

 

u ∥Pm ∥∞ .

(44)

k=1

However, the error in (44) depends on the column dimension nk of P (k) which decreases in a non-uniform way. But, since nk ≤ n + 1 for all k = 1, 2, . . . , ℓ, we can establish a higher (theoretical) bound which characterizes the overall numerical stability of the ERES algorithm in finite precision arithmetic.

64

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67

Theorem 4 (Numerical Stability of the ERES Method). Given a set of m real univariate polynomials of maximum degree n ∈ N, Pm ∈ Rm×(n+1) the basis matrix, and εt > 0 a small tolerance, the iterative application of the ERES operations to P (1) := Pm using numerical finite precision computations results in a matrix P (ℓ+1) ∈ Rrℓ+1 ×(n+1) with numerical ε -rank = 1 which satisfies the equation: P (ℓ+1) = R · Pm · S + E

Rrℓ+1 ×m

(45) (n+1)×(n+1)

where R ∈ and S ∈ R represent the combined row transformations (ERE operations) and the application of the shifting operation, respectively. The matrix E provides the total numerical error, such that

  ∥E ∥∞ ≤ ℓ (n + 1)3 ρ + O(n2 ) u ∥Pm ∥∞

(46)

where ℓ denotes the total number of applications of the ERES operations to the basis matrix Pm . In (45) the matrix R is defined as the matrix product of all R(k) as given in (39), and S is defined as the matrix product of all S for k = 1, 2, . . . , ℓ. In (46) the term ρ denotes the growth factor which corresponds to the first Gaussian elimination for the basis matrix Pm , and O(n2 ) denotes a polynomial function of n with maximum degree 2. The proof of the above theorem follows from Lemma 1 and the preceding results. The singular value decomposition (SVD) during the Rank-1 procedure is applied to a copy of the matrix P (k) only when it is required [21]. Thus, the processing of P (k) in the Rank-1 procedure does not numerically affect the data during the iterations the main procedure [23]. The preliminary stage in the SVD algorithm is the bidiagonal reduction of P (k) and in most bidiagonal reduction methods the error is expressed in the following form [24,25]: (k)

P (k) + δ P (k) = U B V t (k)

(47) (k)

∥δ P ∥2 ≤ f (rk , nk ) u ∥P ∥2

(48)

where B is bidiagonal, U and V are orthogonal, and f (rk , nk ) is a modestly growing function of the dimensions of P (k) [24,25], where rk ≤ m and nk ≤ n + 1. The error (48) is not accumulated during the iterations of the main procedure of the ERES algorithm, but a small error of the form (48) for rk = rℓ+1 and nk = d + 1 must be taken into account for the final solution. Therefore, the total numerical error mainly comes from the application of the processes of normalization and Gaussian elimination during the iterative main procedure which is given in (46). 5.3. Remarks on the numerical performance of the ERES algorithm The main advantages of the ERES algorithm are: (a) the processing of all the polynomials of a given set simultaneously by creating an initial basis matrix of the least possible dimensions, and (b) its ability to constantly reduce the dimensions of the initial matrix and hence reduce the amount of data during the processing, which results in fast data processing and lower usage of computer memory. Tests on sets of more than two polynomials have shown that the numerical ERES algorithm behaves well producing sufficiently accurate results [22,27,28]. Example 3. Consider the set of polynomials:

Pm,n

  p1 (s) = 0.87 s4 − 31.14 s3 + 108.21 s2 − 55.38 s − 32.01       p2 (s) = −0.65 s4 + 22.76 s3 − 63.74 s2 + 12.87 s − 32.98 = 3 2 p (s) = −0.16 s + 6.49 s − 46.67 s + 86.33      3  p4 (s) = 0.30 s3 − 9.96 s2 + 10.20 s + 52.38

with m = 4, n = 4 and exact GCD g (s) = s2 − 35 s + 97. In this example, the ERES method, the Matrix Pencil method (MPGCD) [29], the subspace method (SubGCD) [30], the quasi-GCD method (QuasiGCD) [5], and the QRGCD method [9] are used in order to compute the GCD of the given set in a typical 16-digits arithmetic system (machine precision u ≈ 2.2 10−16 ). We evaluate the results by measuring the relative error between the exact and the computed solution, which is given by Rel =

∥v−g ∥2 , ∥g ∥2

where v, g are the coefficient vectors

of the provided solution v(s) and the exact GCD g (s) respectively. The required time of processing in milliseconds is also provided and these results are presented in Table 1. The solution given by ERES has the smallest relative error found and it was obtained in relative short time.  Considering the previous Example 3, it is important to stress that QuasiGCD and QRGCD algorithms1 are designed to work with two polynomials. In the case of many polynomials they work iteratively with two polynomials at a time, but they tend 1 QuasiGCD and QRGCD are parts of the SNAP package in MAPLE 16 (Maplesoft, Waterloo Maple Inc.).

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67

65

Table 1 Comparison of GCD methods (Example 3). Method

Rel. error

Time (ms)

MPGCD QuasiGCD SubGCD QRGCD ERES

1.4367948 10−11 9.1428733 10−13 1.0736840 10−14 1.0123822 10−15 8.7271766 10−16

78 62 15 93 31

to fail, especially in the case of large sets of polynomials. Conversely, ERES, MPGCD, and SubGCD are matrix-based algorithms which are designed to work with all the polynomials simultaneously, either in a direct or iterative way, and produce better results. In Table 2 we present a sample of the results given by the above algorithms for the computation of the GCD of randomly selected sets of many polynomials. Generally, the ERES and SubGCD algorithms succeeded in producing a solution with very small relative error close enough to machine precision (u ≈ 2.2 10−16 ). QRGCD algorithm also succeeded in producing solutions with the same degree as the exact GCD, but for higher values of its tolerance (eps > 10−16 ) and consequently larger relative error. The other two methods, MPGCD and QuasiGCD, failed or produced a solution with smaller degree than the exact GCD in the most cases of large sets of polynomials. More tests and comparison with other methods can be found in [28]. Regarding the total numerical error of the ERES method, it is obvious from (46) that it depends on how many iterations ℓ of the main procedure of the algorithm are performed. Practically, ℓ is much less than the maximum polynomial degree n and it strongly depends on the linear dependence of the coefficient vectors of the polynomials, i.e. the rank of the initial basis matrix. For instance, for large sets of polynomials of high degree if ρ(Pm ) ≪ n, then we expect a number of iterations close to n. Conversely, if ρ(Pm ) = n + 1 − d, where d denotes the degree of the GCD, then the GCD can be computed in just two iterations. The following example demonstrates this case in a general form. Example 4. Consider an arbitrary set of polynomials Pm,n , m = 8, n = 4 and d = 1. Let Pm ∈ Rm×(n+1) be the initial basis matrix in arbitrary form, where its non-zero elements are symbolized with ‘‘∗’’ and assume that ρ(Pm ) = n + 1 − d = 4. Then, by using the ERES method we get: Iteration 1:

 ∗ ∗ ∗  ∗  ∗  ∗  ∗ ∗

∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗

∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗

∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗

  ∗ ∗ ∗ 0  ∗  Gaussian 0 0  ∗ elimination  0  ∗  −→ 0 ∗   0 ∗ ∗ 0

∗ ∗ ∗ ∗

0

0 0

0 ∗ Gaussian 0 0 elimination  0 0 −→ ∗ 0

∗ ∗ 0 0 0 0 0 0

∗ ∗ ∗ 0 0 0 0 0

∗ ∗ ∗ ∗ 0 0 0 0

 ∗ ∗  ∗  Zero row deletion ∗  ∗ ∗ and shifting ∗ 0  −→ ∗ 0 

∗ ∗ ∗ ∗

∗ ∗ ∗ 0

∗ ∗ 0 0

 ∗ 0 0 0

0 0

Iteration 2:



∗ ∗ ∗ ∗

∗ ∗ ∗

∗ ∗





∗ ∗ 0 0

0

∗ ∗ 0

0 0

∗ ∗





0 ∗ 0 Shifting ∗ 0 −→ ∗

∗

∗

∗ ∗ ∗ ∗

0 0 0 0

0 0 0 0



0 0 . 0 0

The last matrix has exactly d + 1 consecutive non-zero columns and, provided that it has rank 1, every row gives the coefficients of the GCD with degree d = 1.  However, there are cases where the iterative nature of the ERES method acts as disadvantage, especially when the number of iterations is high and inexact data are present. GEPP is a quite stable numerical process, but although pivoting keeps the multipliers bounded by unity, the elements in the reduced matrices still can grow arbitrarily [24,26] during the iterations. Therefore, in fixed-precision arithmetic the error in (44) may become prohibitively large. This problem motivated the search for a new kind of implementation for the ERES method which improves its performance and reliability and made it suitable for computing an approximate GCD for sets of polynomials with inaccurate data. A major step towards this direction is the usage of different types of arithmetic, such as rational and numeric (variable floating-point) arithmetic. The benefits from the mixture of rational and numeric computations, known as hybrid computations, are significant and thus hybrid computations are widespread nowadays. The development of the hybrid implementation of the ERES algorithm optimized for the approximate GCD problem has been described in [23], and it is referred to as the Hybrid ERES algorithm (H-ERES). More information about the hybridization of the ERES method and comparison with other methods can be found in [23,28].

66

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67 Table 2 Comparison of GCD methods for random sets of polynomials. Set

Method

Tolerance

Rel. error

Time (ms)

(A) m = 10 n = 10 d=3

ERES SubGCD MPGCD QRGCD QuasiGCD

10−16 10−16 10−8 10−15 10−15

5.0398803 10−16 4.3369269 10−15 4.0010207 10−10 1.3419339 10−15 4.0030939 10−6

125 94 202 109 94

(B) m = 20 n = 40 d = 10

ERES SubGCD MPGCD QRGCD QuasiGCD

10−16 10−16 FAIL 10−10 FAIL

4.5948046 10−16 1.2802986 10−13 FAIL 2.1016595 10−10 FAIL

8549 4181 FAIL 4024 FAIL

(C) m = 50 n = 20 d=5

ERES SubGCD MPGCD QRGCD QuasiGCD

10−16 10−16 10−5 10−12 FAIL

2.4916347 10−16 2.2638888 10−14 9.9727246 10−1 7.3988778 10−13 FAIL

343 484 827 750 FAIL

m = number of polynomials in the set, n = maximum degree of polynomials, d = degree of the exact GCD.

6. Conclusions In this paper, the fundamental theoretical and numerical properties of the ERES method for computing the greatest common divisor of sets of many polynomials were presented and analysed. The general algebraic representation of the method is presented in Theorem 3 and requires the shifting operation to be written as a matrix product just like the elementary row operations. In Theorem 2 it is now proven that the shifting operation applied to upper trapezoidal matrices with full rank can be represented as a simple matrix product. The main problem in this study was to analyse the overall numerical stability of the ERES method. In Theorem 4, a total numerical error bound is now established for all the iterations of the ERES algorithm which indicates that, under certain conditions, the method is numerically stable for large sets of polynomials. These results are demonstrated through numerical examples. The accuracy of the solutions given by the ERES algorithm reveals that ERES is an efficient method for the computation of the greatest common divisor of sets of many polynomials compared to other GCD methods. References [1] T.L. Heath, The Thirteen Books of Euclid’s Elements, second ed., Dover Publications, 1956. [2] W. Blankiship, A new version of Euclid’s algorithm, American Mathematics Monthly 70 (1963) 742–744. [3] W.S. Brown, On Euclid’s algorithm and the computation of polynomials greatest common divisors, Journal of the Association for Computer Machinery 18 (4) (1971) 478–504. [4] I.S. Pace, S. Barnett, Comparison of algorithms for calculation of g.c.d of polynomials, International Journal of Control 4 (2) (1973) 211–216. [5] A. Schönhage, Quasi-GCD computations, Journal of Complexity 1 (1985) 118–137. [6] M.T. Noda, T. Sasaki, Approximate GCD and its applications to ill-conditioned algebraic equations, Journal of Computational and Applied Mathematics 38 (1991) 335–351. [7] N. Karcanias, C. Giannakopoulos, M. Hubbard, Almost zeros of a set of polynomials of R[s], International Journal of Control 38 (1983) 1213–1238. [8] P. Chin, R.M. Corless, G.F. Corliss, Optimization strategies for the approximate GCD problem, in: Proc. ISSAC’98, Rostock, Germany, 1998, pp. 228–235. [9] R.M. Corless, S.M. Watt, L. Zhi, QR factoring to compute the GCD of univariate approximate polynomials, IEEE Transactions on Signal Processing 52 (12) (2004) 3394–3402. [10] I.Z. Emiris, A. Galligo, H. Lombardi, Certified approximate univariate GCDs, Journal of Pure and Applied Algebra 117 & 118 (1997) 229–251. [11] N. Karmarkar, Y.N. Lakshman, On approximate GCDs of univariate polynomials, Journal of Symbolic Computation 26 (6) (1998) 653–666. [12] V.Y. Pan, Computation of approximate polynomial GCDs and an extension, Information and Computation 167 (2001) 71–85. [13] D. Rupprecht, An algorithm for computing certified approximate GCD of n univariate polynomials, Journal of Pure and Applied Algebra 139 (1999) 255–284. [14] Z. Zeng, A method computing multiple roots of inexact polynomials, Mathematics of Computation 74 (2005) 869–903. [15] D.A. Bini, P. Boito, Structured matrix-based methods for polynomial ϵ -gcd: analysis and comparison, in: Proc. ISSAC’07, Waterloo, Ontario, Canada, 2007, pp. 9–16. [16] J.R. Winkler, X. Lao, The calculation of the degree of an approximate greatest common divisor of two polynomials, Journal of Computational and Applied Mathematics 235 (2011) 1587–1603. [17] S. Barnett, Greatest common divisor of several polynomials, Proceedings of the Cambridge Philosophical Society (1971) 263–268. [18] S. Barnett, Greatest common divisor from generalized Sylvester matrices, Proceedings of the Cambridge Philosophical Society 8 (1980) 271–279. [19] L. Gonzales-Vega, An elementary proof of Barnett’s theorem about the greatest common divisor of several univariate polynomials, Linear Algebra and its Applications 247 (1996) 185–202. [20] N. Karcanias, Invariance properties and characterisation of the greatest common divisor of a set of polynomials, International Journal of Control 46 (1987) 1751–1760. [21] M. Mitrouli, N. Karcanias, Computation of the GCD of polynomials using Gaussian transformation and shifting, International Journal of Control 58 (1993) 211–228. [22] M. Mitrouli, N. Karcanias, C. Koukouvinos, Further numerical aspects of the ERES algorithm for the computation of the greatest common divisor of polynomials and comparison with other existing methodologies, Utilitas Mathematica 50 (1996) 65–84. [23] D. Christou, N. Karcanias, M. Mitrouli, The ERES method for computing the approximate GCD of several polynomials, Applied Numerical Mathematics 60 (2010) 94–114.

D. Christou et al. / Journal of Computational and Applied Mathematics 260 (2014) 54–67 [24] [25] [26] [27]

67

B.N. Datta, Numerical Linear Algebra and Applications, second ed., SIAM, Philadelphia, USA, 2010. G.H. Golub, C.F. Van Loan, Matrix Computations, second ed., The John Hopkins University Press, Baltimore, London, 1989. J.H. Wilkinson, Rounding Errors in Algebraic Processes, Her Majesty’s Stationary Office, London, 1963. 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, Journal of Computational and Applied Mathematics 76 (1996) 89–112. [28] D. Christou, N. Karcanias, M. Mitrouli, D. Triantafyllou, Numerical and symbolical methods for the GCD of several polynomials, in: Numerical Linear Algebra in Signals, Systems and Control, in: Lecture Notes in Electrical Engineering, vol. 80, Springer, 2011, pp. 123–144. [29] N. Karcanias, M. Mitrouli, A matrix pencil based numerical method for the computation of the GCD of polynomials, IEEE Transactions on Automatic Control 39 (1994) 977–981. [30] W. Qui, Y. Hua, K. Abed-Meraim, A subspace method for the computation of the GCD of polynomials, Automatica 33 (4) (1997) 255–284.