Structured Matrix Methods Computing the Greatest Common Divisor of Polynomials

Abstract This paper revisits the Bézout, Sylvester, and power-basis matrix representations of the greatest common divisor (GCD) of sets of several polynomials. Furthermore, the present work introduces the application of the QR decomposition with column pivoting to a Bézout matrix achieving the computation of the degree and the coeffcients of the GCD through the range of the Bézout matrix. A comparison in terms of computational complexity and numerical effciency of the Bézout-QR, Sylvester-QR, and subspace-SVD methods for the computation of theGCDof sets of several polynomials with real coeffcients is provided.Useful remarks about the performance of the methods based on computational simulations of sets of several polynomials are also presented.


Introduction
The greatest common divisor (GCD) of a polynomial set is proven to be very important to many applications in applied mathematics and engineering. Several methods have been proposed for the computation of the GCD of sets of polynomials. Most of them are based on the Euclidean algorithm. They are designed to process two polynomials at a time [5] and can be applied iteratively when a set of more than two polynomials is considered [18,21]. Conversely, there exist e cient matrix-based methods which can compute the degree and the coe cients of the GCD by applying speci c transformations to a matrix formed directly from the coe cients of the polynomials of the entire given set [2,16,19] The greatest common divisor has a signi cant role in Control Theory, Network Theory, signal and image processing [14,24] and in several other areas of mathematics. A number of important invariants for Linear Systems rely on the notion of the greatest common divisor of many polynomials. In fact, it is instrumental in de ning system notions such as zeros, decoupling zeros, zeros at in nity or notions of minimality of system representations. Conversely, Systems and Control methods provide concepts and tools, which enable the development of new computational procedures for GCD [17]. A major challenge for the control theoretic applications of the GCD is that frequently we have to deal with a very large number of polynomials. It is this requirement that makes the pairwise type approaches for GCD not suitable for such applications [25]. However, matrix-based methods tend to have better performance and quite good numerical stability, especially in the case of large sets of polynomials, because they use the entire set of polynomials.
Barnett's GCD method [2] is a well known method for computing the GCD of several polynomials through the construction of the companion matrix C A of a properly selected polynomial from the given set and the  , g) or Bez(f (s), g(s), is an n × n symmetric matrix which is constructed from the coe cients of the polynomials as follows: More speci cally, the elements b i,j of the Bézout matrix are given by where l = min(i − , j − ), ur = υr = if r > n, and |ur υs| = us υr − ur υs.
The next two theorems, presented in [9,10], refer to the properties of the GCD of polynomials through Bézout matrices.

Theorem 1. Let f (s) and g(s) two polynomials in one variable as given in De nition 1. The greatest common divisor of the polynomials f (s) and g(s), denoted by
If c , c , ..., cn are the columns of the Bézout matrix B(f , g) with rank n − k, then i) the last n − k columns, i.e. c k+ , . . . , cn, are linearly independent, and ii) every column c i for i = , , . . . , k can be written as a linear combination of c k+ , . . . , cn : with d a non-zero real number. Then, the GCD of the polynomials f and g, denoted by gcd(f , g), is Remark 1. Let f , g be two polynomials of degree n and p, respectively, and let k = max{n, p}. Then deg{gcd(f , g)} = k − rank(B(f , g)) or equivalently rank(B(f , g)) = k − deg{gcd(f , g)} ≤ k. The equality holds when the polynomials are coprime. Otherwise, rank(B(f , g)) < k, which means that the Bézout matrix is rank de cient.
An important issue arising from Theorem 2 is the determination of the coe cients of the GCD of the entire set of polynomials. Next, exploiting the rank de ciency property of the Bézout matrix when a non-trivial GCD exists, we propose the application of the rank revealing QR factorization to a Bézout matrix. Theorem 3. (QR factorization with column pivoting (QRCP) for rank de cient Bézout matrices). Let B ∈ R n×n and rank(B) = r < n, where B is a Bézout matrix as de ned in (1). Then, there always exist a permutation matrix Π of order n and a n × n orthogonal matrix Q [12] such that  (4), we can directly obtain the coecients d i of the gcd(f , g) through the Bézout-QRCP method. (This is fully demonstrated in Example 1 in Section 3). The application of the QRCP method to Bézout matrices simultaneously reveals the rank and an orthogonal base for the range of the Bézout matrix. Thus, by following Theorem 2 the coe cients of the GCD can easily be determined in a more e cient way.
Considering the case of sets of several polynomials, the following de nition of an extended form of the Bézout matrix is given. De nition 2. We consider the set of m + real univariate polynomials: De nition 3. Let u, v , . . . , vm be m + polynomials, with u a polynomial of maximal degree n. Let B i be the Bézout matrix of polynomials u, v i , for i = , . . . , n. Then the generalized Bézout matrix is de ned as follows: Remark 3. Theorems 1, 2, and 3 also hold for the generalized Bézout matrix.

. Representation of the GCD through Sylvester matrices
Let R[s] be the ring of real polynomials in one variable. We consider two polynomials a(s), b(s) ∈ R[s], with degrees deg{a(s)} = n and deg{b(s)} = p, respectively, where p ≤ n. a(s) = an s n + a n− s n− + . . . + a s + a , an ≠ b(s) = bp s p The resultant matrix or Sylvester matrix S ∈ R (n+p)×(n+p) of the two polynomials a and b is de ned by an a n− a n− .  In order to reduce the computational complexity of the QR factorization applied to S, we construct the modi ed Sylvester matrix, S * , through row interchanges as it is described in [26]. The form of S * is the following: an a n− . . . an−p a n−p− a n−p− .
The following result provides a matrix representation of the standard factorization of the GCD of a set of several polynomials based on Sylvester and Toeplitz-like matrices. De nition 4. We consider the set of polynomials P ≡ P m+ ,n as de ned in (9): i) We can de ne a p × (n + p) matrix associated with a(s): an a n− a n− · · · a a · · · an a n− · · · · · · a a . . .
. . . · · · an a n− · · · · · · a a        and n × (n + p) matrices associated with each b i (s), i = , , . . . , m: A generalized Sylvester matrix or generalized resultant for the set P is de ned by: ii) The matrix S P is the basis matrix of the set of polynomials which is also referred to as the Sylvester Resultant set of the given set P, [11,29].
Theorem 4 also holds for the generalized Sylvester matrix S P which has a special structure. By reordering its blocks through row-interchanging, we can construct the modi ed generalized Sylvester matrix S * P which has n same blocks [26]. Thus, it can be handled more e ciently in respect of the required number of oating-point operations during the implementation of the QR factorization.
an a n− . . . an−p a n−p− a n−p− . . . a a a . . . an . . . a n−p+ an−p a n−p− . . . a a a a . . . .

Representation of the GCD through a power-basis matrix
We consider again a set of real polynomials in one variable as de ned in (9). The polynomials b i (s) are presented with respect to the maximum degree n of the leading polynomial a(s) as where Then, for any set P m+ ,n , a vector representative p(s) and an associated matrix P m+ ∈ R (m+ )×(n+ ) are de ned by . . , m. The next theorem [7] provides a representation of the GCD by using a matrix which is constructed directly from the coe cients of the polynomials according to a power-basis vector. Theorem 5. Let P m+ ,n a set of m real univariate polynomials of maximum degree n ∈ N and P m+ ∈ R (m+ )×(n+ ) the matrix which is constructed according to the power-basis vector e n (s). The application of elementary row operations and shifting to P m+ results in a matrix G ∈ R (m+ )×(n+ ) with rank(G) = , which satis es the equation: represent the applied elementary row operations and the application of the shifting operation, respectively. The last non-zero row of G provides the coe cients of the GCD of the set P m+ ,n .

Remark 4.
It is important to stress that the shifting matrix S does not a ect the numerical GCD solution, but it plays a signi cant role when a symbolic form of the GCD is under consideration. Therefore, it is always computed in symbolic-rational form (for more details the interested reader may refer to [7]).
The advantages of the representation of the GCD through a power-basis matrix are the following: i) Implementation using symbolic-rational computations.
The relation (18) is particularly useful when symbolic-rational computations are involved, because it provides a compact way to represent the GCD solution as a product of matrices with the least possible dimensions allowed by the input data (see [7] for more details).
Recent experimental results showed that (18) provides a direct way to compute the GCD of sets of multivariate polynomials, when the procedure that constructs the initial matrix is appropriately adjusted. Considering the case of sets of m + polynomials in two variables (s, t) (bivariate polynomials) with real coe cients, the a basis matrix P m+ can be formed according to the bivariate power-basis: where n, r are the maximum powers of the variables s, t, respectively. The dimension of the corresponding basis vector e n,r (s, t) is equal to (n + ) · (r + ). Similar base vectors can be formed for polynomials in several variables.
If we consider the column vectors e n (s) = [s n , . . . , s, ] T and e r (t) = [t r , . . . , t, ] T , then the matrix P m+ is structured according to the bivariate power-basis vector: If we consider the power-basis vectors e n,r (s, t) and e n·r+n+r (s), we can see that both vectors can yield the same power-basis matrix P m+ for a given set of polynomials. Therefore, we may assume a one-toone correspondence between the two power-basis vectors, and if we apply Theorem 5 to a set of bivariate polynomials of degrees (n, r), we can treat them as polynomials of maximum degree n r + n + r and the column dimension of the matrix P m+ will be equal to (n + )(r + ). An analytical example for the case of bivariate polynomials is presented in the appendix.

Implementation of the GCD computation through matrix factorizations
The computation of the degree of the GCD of polynomials is an issue of great importance in GCD computations. The computation of the coe cients of the GCD depend on the correct determination of the degree of the GCD either in the Sylvester, or the Bézout case. In both cases QR factorization can be used for the e cient computation of the Sylvester matrix S P and the Bézout matrix B. However, several other approaches have been proposed for the computation of the degree of the GCD for the speci c types of matrices.
Winkler et al. [28] proposed two methods for the calculation of the degree of the GCD of polynomials. The rst method is based on the two subspaces that are enabled from the partitioned structure of the Sylvester matrix. The degree of the GCD is computed through monitoring the change of the angle of these two subspaces, by deliting rows and columns of the Sylvester matrix. The second method is based on considering the change in the error between two estimates of the GCD of the polynomials as a function of its degree. The properties and the application of the Sylvester matrix to GCD computations is also studied in [27].
Kaltofen et al. [15] presented a di erent method based on structured total least norm (STLN) algorithms applied to Sylvester matrices for computing the GCD of polynomials and Bini and Boito [4] presented a fast algorithm for computing the GCD of two polynomials, which uses Bézout and Sylvester matrices. More particularly, this algorithm is based on the reduction of their displacement structure to Cauchy-like one.

. Computational complexity
Next,we analyze the numerical complexity of QR-based methods applied to the generalized Sylvester or the Bézout matrices. These factorizations compute an orthogonal matrix Q and an upper triangular matrix R which can produce the coe cients of the GCD of the polynomials. The computational complexity is measured in ops, where 1-op corresponds to the computational time required for one multiplication and one addition of oating-point numbers.
. . Computation of the GCD through QR factorization a) Sylvester-QR factorization Theorem 6 ([2, 3]). Let S P be the generalized Sylvester matrix of m + polynomials. If S P = QR is the QR factorization of S P , the the last non zero row of R gives the coe cients of the GCD of the polynomials.
The complexity of the QR factorization applied to S P is high. More precisely, for m + polynomials of maximum degree n and second maximum degree p the generalized Sylvester matrix is of size (mn + p) × (n + p) and thus, the required complexity is ops. If m n p the complexity is O( n ). The modi ed Sylvester matrix S * P has n same blocks and, if we properly exploit the special structure of S * P (see (16)) in QR factorization, the required complexity is decreased to ops. If m n p, the required ops are about O( n log n), [26].

. . Computation of the GCD through QR factorization with column pivoting (Bézout-QRCP method)
Since the mn × n Bézout B matrix is always rank de cient when a non-trivial GCD exists, it is more e cient to extract the coe cients h i appeared in (4) using Remark 2, which indicates that the coe cients of the GCD of the polynomials can be derived from the QRCP factorization of the Bézout matrix. The complexity of the QRCP factorization is O mn r − r (mn + n) + r ops [8], where r is the rank of B, which is less than the ops required by the classical QR factorization. The appropriate correspondence of the columns of the original and the permuted matrix, which reveal the GCD coe cients (Remark 2), is symbolically implemented. In the case where the rank de ciency of B is high, the Bézout-QRCP method becomes more e cient.

. . Computation of the GCD through Singular Value Decomposition (Subspace method)
A similar approach with Sylvester-like matrices is the method presented in [23]. Given a set of univariate polynomials P m+ ,n , the rst two steps of the developed algorithm in [23] involves the construction of an (m + )(n + ) × ( n + ) generalized Sylvester matrix¹Ŝ P from the input polynomials (as demonstrated in Example 2) and the computation of the left null space of the transposedŜ T P via singular value decomposition. If we denote by U ∈ R ( n+ )×k the basis matrix for the computed left null space ofŜ T P and C is the ( n + ) × ( n + − k) Toeplitz matrix of a polynomial of degree k with arbitrary coe cients, then the GCD vector is actually the unique (up to a scalar) solution of the system U T C = . The degree of the GCD is k = dim{range{U }}. The computational cost of this method, which generally is rather high, is dominated by the singular value decomposition of the generalized Sylvester matrixŜ P , which requires O( m n + m n ) ops. If m n the required ops are about O( n ). Table 1 summarizes the required computational complexity for each of the aforementioned methods.

. . Remarks upon the computational complexity of the methods
As it is shown in Table 1, for a given polynomial set P m+ ,n , the complexity of the Sylvester and the modied Sylvester QR is a function of the maximum polynomial degrees n and p. The required ops of the modi ed Sylvester QR are signi cantly less than that of the classical Sylvester QR. When p << n the complexity of the modi ed Sylvester QR can be further reduced. If n p, the Sylvester QR method requires O n m ops, while for its modi ed version O n log n + + n m log n ops are needed. Since log n becomes signi cantly less than n as the degree n increases, the modi ed Sylvester QR becomes more and more ecient than the classical one for higher degrees n. When m n p, the complexity of the modi ed method is decreased by one order comparing to the complexity of the classical method. The modi ed Sylvester QR has remarkable performance when applied to sets of many polynomials with high degrees.
The Bézout QRCP exploits the rank de ciency n−r of the matrices, which is equal to the degree of the GCD of the polynomials. Thus, the higher the GCD degree is (i.e higher rank de ciency of the Bézout matrix) the more e cient the method becomes. If the rank of the Bézout matrix r is signi cantly less than the maximum degree n of the polynomials, then the complexity of the Bézout QRCP method is one order less comparing to the complexity of the classical Bézout QR. When m n p, the complexity of the Bézout QRCP method is less than that of the modi ed Sylvester QR for r < log n + . The subspace method requires signi cantly more oating-point operations than the other methods. This disadvantage is o set by the fact that this method gives more accurate values (i.e. lower relative errors) for the coe cients of the computed GCD comparing to the other methods.

. Demonstrative Examples
The following example demonstrates the steps of the current Bézout-QRCP method for computing the GCD of set of many polynomials. Example 1. We consider the pair of real univariate polynomials of degree 5: 1 We must note that the constructed matrixŜ P , although it is de ned as a generalized Sylvester matrix in [23], it is slightly di erent than the matrix de ned by (15), because the second maximum degree p of the polynomials is not considered at all.  The exact GCD is s − s + . The Bézout matrix of the given polynomials in the set P , is where bc i , i = , , . . . , are the columns of the initial Bézout matrix B ∈ R × . The following factorization is achieved by applying the QR factorization with column pivoting (QRCP) to where and After applying the QRCP factorization, the permuted Bézout matrix Bperm = B · Π is = bc bc bc bc bc = bc bc bc bc bc The lowest right × part of R is considered to be zero and, thus, QRCP indicates that r = rank(B) = and deg{gcd(P , )} = − = . From Theorem 2 we know that the last 3 columns of the initial Bézout matrix B in (25), i.e. bc , bc , and bc , are linear independent. Therefore, the rst two columns of B, bc and bc , can be written as a linear combination of bc , bc and bc . Thus, from (4) in Theorem 2 we have: Let d s + d s + d be the GCD of the polynomials. The coe cients h ( ) and h ( ) give the coe cients d and d , respectively, and the constant term d is 1.
Using QRCP, the coe cients h ( ) and h ( ) of the GCD are derived from the correspondence of the columns of B and Bperm. According to Theorem 3, the columns q , q and q of Q generate the range of Bperm. From (8) we have: Since the columns bc and bc of the initial Bézout matrix B correspond to bc and bc of the permuted Bézout matrix Bperm, respectively, it is necessary to express the columns bc and bc as linear combinations of the columns bc , bc and bc . Since each column q i , i = , , is given by an analytic formula as the solution of the lower triangular system, formed from the rst, the second, and the fourth equation of (33), we symbolically substitute in the third and the fth equation of (33) and we obtain: Therefore, we conclude that and we obtain the quadratic polynomial: . s − . s + If we convert it to a monic polynomial, dividing by . , we nally compute the GCD of the polynomials in P , . That is In the following, we consider the computation of the GCD using the Bézout-QR, Sylvester-QR, power-basis, and the subspace-SVD methods. Example 2. Let us consider the next set of three univariate polynomials: of degree . Their exact GCD is s − s + .

i) GCD from Bézout matrices using QRCP decomposition
The generalized Bézout matrix of the given polynomials in the set P , is and and c , c , c are the columns of B.
We apply the QRCP factorization to B, such that where q , q , q , q , and q are the columns of Q. The lowest right × part of R is zero and thus, QRCP indicates that r = rank(B) = . The degree of the GCD is deg{gcd(P , )} = − r = . Theorem 2 denotes that the last column b of the initial Bézout matrix B in (38) is linear independent and the other columns b and b are multiples of b . Working similarly with Example 1 we conclude that: ii) GCD from Sylvester matrices using QR decomposition.
The generalized Sylvester matrix of the polynomials in P , , as de ned by (15), is: Applying the QR factorization (without column pivoting), we have The last non zero row of R gives the coe cients of the GCD. If the fourth row of R is divided by the element R , , then gcd(P , ) = s − . s + . .

iii) GCD from Sylvester matrices using singular value decomposition.
A di erent approach [23] can be followed if we compute the left null space of the generalized Sylvester matrix S P by using the singular value decompositionŜ T P = U Σ V T , wherê , . , . , . , .
when we divide every element by . .

iv) GCD from a power-basis matrix.
The power-basis matrix of the polynomials in the set P , is for e (s) = [s , s, ] T . If we apply Theorem 5 using symbolic-rational computations, we obtain the matrices and the rows of the matrix provide the coe cients of the GCD up to a scalar multiple.

Numerical examples
In this section we compare the Bézout-QRCP method, which is proposed in this work, the Bézout-QR developed in [6], the Sylvester-QR for the modi ed case, and the subspace-SVD method [23]. All the aforementioned methods are based on numerically stable procedures such as QR decomposition (with or without column pivoting) and singular value decomposition [12] which are widely used in numerical applications and are included in high-performance computational software packages, such as Matlab and Maple. We run several computational simulations in Matlab on an AMD-A6 Dual-Core 3.6 GHz -8Gb machine using many di erent sets of polynomials with randomly selected coe cients in order to test the numerical behavior and performance of the described methods by measuring the processing time ( Fig. 1) and the relative error ( Fig. 2) when an exact GCD is known. The results obtained showed that, in general, all methods can provide reliable results within an acceptable range of accuracy.
Furthermore, those results were obtained after normalizing the rows of the input matrix (i.e. the elements each row are divided by the norm of the row) and were slightly better than those obtained without normalization. This is a scaling technique which is widely used and generally improves the accuracy of the results. In the current computational simulations we used the Euclidean norm. Two operations, scaling by the geometric mean and relative scaling are also suggested in [28]. These operations are applied to the coe cients of the initial polynomials before the factorization of the Sylvester or the Bézout matrix. Example 3. In Table 2 and Table 3 we summarize the results obtained regarding the numerical relative error for the computed GCD of the polynomial sets in Example 1 and Example 2, respectively.

Algorithm
Tolerance Rel. Error The tolerance indicates the di erent levels of precision (numerical accuracy) where a number is considered to be zero. For the particular sets of polynomials a tolerance between − and − was selected with an exception for the subspace-SVD method which couldn't determine the correct degree of the GCD when the tolerance was smaller than − .

Conclusions
In this paper, we proposed the application of the QR factorization with column pivoting to a Bézout matrix in order to compute the coe cients of the GCD of sets of several polynomials in a more e cient way. We also presented an overview of the most frequently applied structured matrix-based representations: i) the Sylvester QR, ii) the Bézout QR, and iii) the subspace SVD method. All these methods start with a structured matrix constructed directly from the coe cients of the given polynomial set. The size of the block banded generalized Sylvester matrix is larger than that of the generalized Bézout, which consists of blocks of symmetric matrices. The subspace method involves a Sylvester-like matrix of biggest dimensions. Taking into account the structure of each of the above matrices, we compared the methods theoretically with respect to their computational complexity.
From this study, we can recommend the application of the most appropriate method according to the given polynomial set. Speci cally, for a large set of polynomials with high degree, the modi ed Sylvester QR method is more preferable. As the number of the polynomials in the set decreases the Bézout QRCP becomes more e cient. Additionally, the results obtained from the simulations considering an exact GCD showed that both methods produced accurate results in reasonable time limits. However, the subspace-SVD method appears to be more stable for large sets of polynomials, but the processing time increases at a high rate as the number of polynomials increases.
The study of the approximate GCD case is also a topic of great interest. A thorough comparison among the existing methods and possible extension of the QRCP method to the approximate case is under consideration. Furthermore, a proper framework for the algebraic and geometric properties of the GCD of sets of many polynomials in a multidimensional space is currently under study in order to de ne and evaluate exact or approximate multivariate GCDs given by the QRCP method. This is a challenging problem for further research, because several real-time applications, such as image and signal processing, rely on GCD methods where multivariate polynomials (especially in two variables) are used.
If we apply Theorem 5 using symbolic-rational computations, we nally obtain the matrices R ∈ R × and S ∈ R × as given in (54) and (55), respectively, where such that The last non-zero row provides the coe cients of the GCD. Actually, any non-zero row of G provides the GCD up to a scalar multiple. The GCD vector is g = [ , , , , , , , − , ] and gcd(P) = g · e , (s, t) = s t − t +