summary:The computation of the greatest common divisor (GCD) has many applications in several disciplines including computer graphics, image deblurring problem or computing multiple roots of inexact polynomials. In this paper, Sylvester and Bézout matrices are considered for this purpose. The computation is divided into three stages. A rank revealing method is shortly mentioned in the first one and then the algorithms for calculation of an approximation of GCD are formulated. In the final stage the coefficients are improved using Gauss-Newton method. Numerical results show the efficiency of proposed last two stages

Eckstein, Jiří

Zítko, Jan

English

Institute of Mathematics AS CR, v. v. i.

PANM 17Jiří Eckstein; Jan ZítkoComparison of algorithms for calculation of the greatest common divisor of several polynomialsIn: Jan Chleboun and Petr Přikryl and Karel Segeth and Jakub Šístek and Tomáš Vejchodský (eds.): Programs andAlgorithms of Numerical Mathematics, Proceedings of Seminar. Dolní Maxov, June 8-13, 2014. Institute ofMathematics AS CR, Prague, 2015. pp. 64--70.Persistent URL: http://dml.cz/dmlcz/702664Terms of use:© Institute of Mathematics AS CR, 2015Institute of Mathematics of the Academy of Sciences of the Czech Republic provides access to digitized documentsstrictly for personal use. Each copy of any part of this document must contain these Terms of use.This paper has been digitized, optimized for electronic delivery and stamped with digitalsignature within the project DML-CZ: The Czech Digital Mathematics Libraryhttp://project.dml.czPrograms and Algorithms of Numerical Matematics 17J. Chleboun, P. Přikryl, K. Segeth, J. Š́ıstek, T. Vejchodský (Eds.)Institute of Mathematics AS CR, Prague 2015COMPARISON OF ALGORITHMS FOR CALCULATIONOF THE GREATEST COMMON DIVISOROF SEVERAL POLYNOMIALSJǐŕı Eckstein, Jan ŹıtkoDepartment of Numerical Mathematics,Faculty of Mathematics and Physics, Charles UniversitySokolovská 83, Prague 8, Czech Republicjiri.eckstein@gmail.com, jan zitko@centrum.czAbstractThe computation of the greatest common divisor (GCD) has many applications inseveral disciplines including computer graphics, image deblurring problem or comput-ing multiple roots of inexact polynomials. In this paper, Sylvester and Bézout matricesare considered for this purpose. The computation is divided into three stages. A rankrevealing method is shortly mentioned in the first one and then the algorithms for cal-culation of an approximation of GCD are formulated. In the final stage the coefficientsare improved using Gauss-Newton method. Numerical results show the efficiency ofproposed last two stages.1. IntroductionSylvester matrices (see [1, 3, 5, 6, 10, 11, 14, 15, 16, 17]) or Bézout matrices(see [3, 7, 8, 12]) can be used for the calculation of GCD. We start with Sylvestermatrix. The coefficients of GCD of two polynomials f1 and f2 can be obtainedfrom a Sylvester subresultant Sk(f1, f2) which is formed from the Sylvester matrixS(f1, f2) by deleting the last k − 1 rows, the last k − 1 columns of the coefficientsof f1 and the last k − 1 columns of the coefficients of f2. If ni =deg(fi) for i = 1, 2,n1 ≥ n2, and if for a positive integer d ≤ n2 the subresultant Sd(f1, f2) is the firstrank deficient matrix in the sequenceSn2(f1, f2), Sn2−1(f1, f2), . . . , S1(f1, f2), (1)then d =deg(GCD(f1, f2)). There are two well-known procedures for calculation ofthe rank (or rank deficiency) of a matrix. For small dimension the usage of SVD issufficient (see for example [2] or [9]). The numerical rank revealing algorithm withmany robust examples is in detail described in the papers [11, 17]. The whole processis in details, together with the calculation of GCD and rank determination, explained64in the papers [17], [11], [14] and therefore, in the following, we are considering thecalculation for m polynomials.We now consider m real polynomials f1, f2,. . . fm. Let ni =deg(fi). Denoteg =GCD(f1, f2, . . . , fm). It is assumed that d=deg(g) > 0. The objective is to findpolynomials w1, w2, . . . , wm of degrees n1 − d, n2 − d,. . . , nm − d respectively, suchthat fi = wi g for all considered i, which can be expressed in the form (see [4, 16, 17])Cd(wi)~g = ~fi, for i ∈ {1, 2, . . . , m}, (2)where Cd(wi) is the Cauchy matrix for the polynomial wi with d + 1 columns, i.e.,Cd(wi) ∈ R(ni+1)×(d+1). The symbol ~g denotes the vector of coefficients of g and thesymbols ~fi and ~wi have an analogous meaning. The system (2) can be rewritten inthe formF (x) = b, where x =~g~w1~w2···~wm, b =1~f1~f2···~fm, F (x) =(~r)T~gCd(w1)~gCd(w2)~g···Cd(wm)~g(3)and ~r is a scaling vector (see [16]). Let us remark that ~r,~g ∈ Rd+1, and ~wj ∈ Rnj−d+1.The system (3) represents( m∑j=1nj)+m+1 equations with( m∑j=1nj)+m+1−(m−1) dunknowns and the least square solution (see [9]) is applied. According to the wellknown theory (see [2]) we havegrad[12‖F (x)− b ‖2]= (J(x))T [F (x)− b] , (4)where J(x) is the Jacobian of F and can be easily calculated as a Gateaux derivativeof F . The problem of location of minimum leads to the solution of the system(J(x))T [F (x)− b] = 0. (5)Let us mention the result formulated in [16]: for every scaling vector ~r satisfying~rT~g 6= 0, if GCD(w1, w2, . . . , wm) = 1, then the Jacobi matrix has a full columnrank and therefore F (x) = b. However all these investigations depend on the basicquestion how to find the rank d. This is well known for m = 2 and it is shortlyanalysed in Section 2. In the next section the algorithm using Sylvester matricesfor m ≥ 3 is discussed. In Section 3, the calculation of the greatest common divisorof several univariate polynomials through Bézout-like matrices is considered. Bothstrategies are numerically tested in the last section.652. Calculation of GCD through Sylvester matricesAt the beginning consider the polynomials f1 and f2 of degrees n1 and n2 re-spectively, where n1 ≥ n2. According to the previous section we determine aninteger d such that Sd(f1, f2) is the first rank deficient matrix in (1) and denote theright singular vector of the matrix Sd(f1, f2) = [Cn2−d(f1), Cn1−d(f2)] correspondingto the smallest singular value σmin (Sd(f1, f2)), which is theoretically equal to zero,by [(~w2)T ,−(~w1)T ]T . We have denoted g =GCD(f1, f2). The coefficients of ~g arecalculated as the least square solution of the equationCd(w2)~g = ~f2 or Cd(w1)~g = ~f1. (6)One of these equations (usually the first one) is solved and the second one is usedfor improvement of the result if it is necessary.However, for three or more polynomials it is impossible to apply an analogoustechnique for finding the degree of GCD(f1, f2, f3). A consecutive process is usuallyapplied, which can be formally written for three polynomials in the formd = deg(GCD(f3,GCD(f1, f2))).Numerically, the determination of GCD is usually based on some minimisationmethod which is formally written by (4), (5) and the realization means an infiniteiterative process where only finite number of iterations is implemented. Moreover, ifthe calculation is performed in floating point environment the result is inexact andtherefore an approximation is obtained as a result of the above mentioned minimiza-tion process. This approximation to GCD will be in this paper entitled approximategreatest common divisor - AGCD. This concept is studied and discussed in many pa-pers (see for example [6, 13, 16]). The concept AGCD is mentioned in context withSTLN algorithm (see [10, 15, 13, 6]). In this paper AGCD is the result of the leastsquare procedures which is realized by the Gauss-Newton method. Exact coefficientsare assumed. Let us consider the system (5). By analogy to [16] and [17] we nowpresent the algorithm for several polynomials. The numerical process will be evidentfrom the following algorithm.Algorithm 2.1 (AGCD for m polynomials.)Input: Real polynomials f1, f2,. . . , fm of degrees n1, n2,. . . , nm respectively, vec-tor b defined by (3) and a given tolerance θ. It is assumed thatn1 ≥ n2 ≥ · · · ≥ nm.Output: Polynomial g =AGCD(f1, f2, . . . , fm)beging := fnmfor j = m,m− 1, . . . , 2 doCalculate g =AGCD(g, fnj−1).66endfor j = 1, m dowj(x) := fj(x)/g(x)endPut d :=deg(g);form the vector (x)T = [(~g)T , ( ~w1)T , ( ~w1)T , . . . , ( ~wm)T ]T for the initialapproximation of Gauss-Newton iteration.repeatx+ = x−~r 0 0 0Cd(w1) Cn1−d(g) ·· 0 · · ·· · · · 0Cd(wm) 0 · · Cnm−d(g)† ~r T~gCd(w1)~g − (~f1)··Cd(wm)~g − ( ~fm)x = x+until ‖F (x)− b‖ < θOnce ‖F (x)− b‖ < θ, we extract coefficients of the polynomial g(x)from the vector x.We now have g(x) =GCD (f1, f2, . . . , fm).end of algorithmThe matrix is a block matrix, the non-zero blocks are Cauchy matrices. It con-tains only zero-blocks except for the first column and the diagonal blocks.3. Calculation of GCD using Bézout matricesWe now present a different approach to computing the GCD of several real uni-variate polynomials using Bézoutian matrices (see [3], [8]). The size of this kind ofmatrix depends purely on the degree of one of the polynomials. It will be possibleto determine the degree of GCD of a whole set of polynomials at once. Moreover,its coefficients will be computed at the same time. Let p and q be two polynomials,p(x) = a0xk + a1xk−1 + . . . ak−1x+ ak,q(x) = b0xk + b1xk−1 + . . . bk−1x+ bkof degrees at most k > 0. If deg(p) >deg(q) then some of the first coefficients of qequal zero.The Bézout matrix associated to p and q (see [12]) isB(p, q) =c1,1 · · · c1,k......ck,1 · · · ck,k,67where the coefficients ci,j are defined by the relationp(x)q(y)− p(y)q(x)x− y=k∑i,j=1ci,jxi−1yj−1.In the following, the procedure for computing the AGCD of m polynomials ispresented. To achieve this, a set of polynomials f1, . . . , fm satisfyingk := n1 = deg(f1) > deg(fi), i = 2, 3, . . . , mwill be assumed. In contrast with Sylvester matrices, all the Bézout matrices B(f1, fi),i = 2, 3, . . . , m are square and of the same dimension. Therefore the matrixBf1(f2, . . . , fm) =B(f1, f2)B(f1, f3)...B(f1, fm)can be constructed. Analogously to computation with Sylvester matrices, the degreeof the AGCD equals k − rank (Bf1(f2, . . . , fm)). Its coefficients can be computed bydetermining the linear combinations of column vectors, as described in the algorithmbelow (for details see [8]). The numerical realization of GCD will be again calledAGCD.Algorithm 3.1 (AGCD for m polynomials.)Input: Real polynomials f1, f2,. . . , fm of degrees n1, n2,. . . , nm respectively. It isassumed that k := n1 > max{n2, . . . , nm}.Output: Polynomial g =AGCD(f1, f2, . . . , fm).beginDetermine the d = k − rank(Bf1(f2, . . . , fm)).Let t1, . . . , tk be column vectors of Bf1(f2, . . . , fm) = [t1, . . . , tk].Construct T2 = [tk, tk−1, . . . , td+1] and T1 = [td, td−1, . . . , t1].Calculate QR decomposition of T2, i.e. T2 = QR, where Q ∈ Rk×k is orthogonaland R ∈ Rk×(k−d) is an upper triangular matrix.We set c := (R)−1k−d,k−d and compute wd+1i = c(QTT1)k−d,i, for i = d, . . . , 1.Setting hi := wd+1d−i+1, i = 1, . . . , d and h0 := 1,we finally have g(x) = h0xd + h1xd−1 + . . .+ hd−1x+ hd = GCD(f1, . . . , fm).end of algorithm684. Numerical experimentTo compare the two presented algorithms, let us now have the following polyno-mials:f0 = (x− 0.9)5(x− 0.8)5(x− 0.7)5(x+ 0.3)5(x+ 0.5)5(x+ 0.7)5,f1 = (x− 2)5(x− 0.9)5(x− 0.8)5(x+ 0.5)5(x+ 2)5,f2 = (x− 3)5(x− 0.8)5(x+ 0.5)5(x+ 2)5 andf3 = (x− 0.8)4(x+ 0.5)4.It is easily seen, that GCD(f0, . . . , f3) = f3. Accuracy of these computations is shownin Table 1. The errors made in determining the coefficients are about two orders ofmagnitude smaller in case of Algorithm 2.1 than in the case of Algorithm 3.1.Coefficients Error in coefficientsGCD Algorithm 2.1 Algorithm 3.11.0000-1.2000-1.06001.33200.5361-0.5328-0.16960.07680.02560.0000e+00-9.1038e-15-6.8834e-152.6645e-151.2212e-15-2.6645e-15-2.0539e-15-7.2164e-16-3.8164e-160.0000e+001.8050e-12-7.6916e-13-2.6426e-123.3151e-131.3358e-121.1419e-13-2.3732e-13-5.7697e-14Table 1: Comparison of computational error in AGCD coefficients produced by Al-gorithm 2.1 and Algorithm 3.1.AcknowledgementsThis work was supported by the grant prvouk p47. The authors thank for thissupport.References[1] Barnett, S.: Polynomials and linear control systems. Marcel Dekker, INC., NewYork and Basel, 1983.[2] Björk, Å.: Numerical method for least square problems. SIAM, Philadelphia, 1996.[3] Bini, D. and Pan, V.Y.: Polynomial and matrix computation, vol. 1 fundamentalalgorithms. Birkhuser, 1994.[4] Corless, R.M., Gianni, P.M., Trager, B.M., and Watt, S.M.: The singular valuedecomposition for polynomial systems. In: Proc. ISSAC 95, pp. 195–200. ACMPress, 1995.69[5] Eliaš, J.: Problémy spojené s výpočtem nejvěťśıho společného dělitele. Bachelorthesis, Charles University, Faculty of Mathematics and Physics, 2009.[6] Eliaš, J.: Approximate polynomial greatest common divisor. Master thesis,Charles University, Faculty of Mathematics and Physics, 2012.[7] Diaz-Toca, G.M. and Gonzales-Vega, L.: Barnett’s theorems about greatest com-mon divisor of several univariate polynomials through Bézout-like matrices. J.Symbolic Computation 34, (2002), 59–81doi: 10.1006/jsco. 2002.0542[8] Diaz-Toca, G.M. and Gonzales-Vega, L. : Computing greatest common divisorsand square free decompositions through matrix method: The parametric andapproximate cases. Linear Algebra Appl. 412 (2006), 222–246.[9] Golub, G.H. and Van Loan, C. F.: Matrix computations. 3rd Ed. John HopkinsUniversity Press, Baltimore, USA, 1996.[10] Kaltofen, E., Yang, Z., and Zhi, L.: Structured low rank approximation ofSylvester matrix. Preprint, 2005.[11] Li, T.Y. and Zeng, Z.: A rank-revealing method with updating, downdatingand applications. SIAM J. Matrix Anal. Appl. 26 No. 4 (2005), 918–946.[12] Pták, V.: Explicit expressions for Bézoutians. Linear Algebra Appl. 59 (1984),43–54.[13] Sun, D. and Zhi, L.: Structured low rank approximation of a Bézout matrix.In:MM Research preprints, pp. 207–218. KLMM, AMSS, Academia Sinica, Bei-jing 2006.[14] Winkler, J. R. and Źıtko, J.: Some questions associated with the calculation ofthe GCD of two univariate polynomials. In: Winter School and SNA’07, pp. 130–137. Ostrava, 2007.[15] Winkler, J. R. and Allan, J.D.: Structured total least norm and approximateGCDs of inexact polynomials. Journal of Comp.and Appl. Math. 215 (2006),1–13.[16] Zeng, Z.: The approximate GCD od inexact polynomials, Part I: univariatealgorithm. Preprint, 2004.[17] Źıtko, J. and Eliaš, J.: Application of the rank revealing algorithm for thecalculation of the GCD. In: Winter School and SNA’12, pp. 175–180. Liberec,2012.70

Comparison of algorithms for calculation of the greatest common divisor of several polynomials

http://dml.cz/bitstream/handle/10338.dmlcz/702664/PANM_17-2014-1_12.pdf

Comparison of algorithms for calculation of the greatest common divisor of several polynomials

Abstract

Similar works

Full text

Available Versions

Institute of Mathematics AS CR, v. v. i.