Re: Approximate/asymptotic factorization

*To*: mathgroup at smc.vnet.net*Subject*: [mg73630] Re: [mg73509] Approximate/asymptotic factorization*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Fri, 23 Feb 2007 04:40:18 -0500 (EST)*References*: <200702210635.BAA17453@smc.vnet.net>

Paul Abbott wrote: > Consider the monomial H[4,4] > > 256*(23625 + 126000*z^4 - 7200*z^8 + 3840*z^12 + 256*z^16) > > where H[m,n] is a Generalized Hermite Polynomial (Noumi & Yamada 1998), > a rational solution to the Painlevï¿½ equation P_IV. A plot of the roots > of H[4,4] displays an interesting quasi-rectangular structure. > > Then consider the monomial Q[4,4] > > 1528899609315374375625 + 6795109374734997225000*z^4 - > 560866170613047390000*z^8 + 153399294526645440000*z^12 + > 2734590598399296000*z^16 - 167891551278796800*z^20 + > 2948686820352000*z^24 - 40649991782400*z^28 + 277762867200*z^32 - > 920125440*z^36 + 1048576*z^40 > > where Q[m,n] is a Generalized Okamoto Polynomial, another rational > solution to P_IV. A plot of the roots displays an even more interesting > structure, with a quasi-rectangular core and 4 triangular "lobes" > (related, in some way, to the the Yablonskiiï¿½Vorobï¿½ev Polynomials, Q[n], > which are Rational Solutions of P_II) > > If you overlay plots of the roots of H[4,4] and Q[4,4] what is most > striking is that the "common" roots are very similar. In other words, > H[4,4] is "approximately" a factor of Q[4,4]. > > Is there a "standard" way of quantifying this approximate factorization? > > This approximate factorization holds for any m and n and appears to > improve as m and n increases. Perhaps there is a limit in which the > factorization holds asymptotically? > > Cheers, > Paul > > _______________________________________________________________________ > Paul Abbott Phone: 61 8 6488 2734 > School of Physics, M013 Fax: +61 8 6488 1014 > The University of Western Australia (CRICOS Provider No 00126G) > AUSTRALIA http://physics.uwa.edu.au/~paul I'm not sure how to tackle the asymptotics you have in mind. For particular examples this is very much the province of univariate approximate gcd computation. There has been a good deal of work in this area since around 1995, with some work even earlier. A good place to start would be the "4M" article of Corless, Gianni, Trager, and Watt, "The singular value decomposition for polynomial systems", ISSAC 1995. The idea of their approach is to use the SVD of the Sylvester matrix in order to determine the degree of a "good" approx gcd, then derive that gcd also using the SVD. For your example I show some quick-and-dirty code that, I think, does more or less what they have in mind. For the Sylvester matrix I cribbed and slightly modified code from MathWorld: http://mathworld.wolfram.com/SylvesterMatrix.html SylvesterMatrix[poly1_, poly2_, var_] := Module[ {coeffs1, coeffs2, l1, l2}, coeffs1 = Reverse[CoefficientList[poly1, var]]; coeffs2 = Reverse[CoefficientList[poly2, var]]; l1 = Length[coeffs1]; l2 = Length[coeffs2]; Join[Table[ Join[Table[0,{j}], coeffs1, Table[0,{l2-2-j}]], {j,0,l2-2}], Table[Join[Table[0,{j}], coeffs2, Table[0,{l1-2-j}]], {j,0,l1-2}]] ] h44 = 256*(23625 + 126000*z^4 - 7200*z^8 + 3840*z^12 + 256*z^16); q44 = 1528899609315374375625 + 6795109374734997225000*z^4 - 560866170613047390000*z^8 + 153399294526645440000*z^12 + 2734590598399296000*z^16 - 167891551278796800*z^20 + 2948686820352000*z^24 - 40649991782400*z^28 + 277762867200*z^32 - 920125440*z^36 + 1048576*z^40; I'll normalize so largest coefficient is unity. normPoly[poly_, var_] := Expand[poly / Max[Abs[CoefficientList[poly,var]]]] p1 = normPoly[h44,z]; p2 = normPoly[q44,z]; I'll use an epsilon of 10^(-3) as tolerance for what singular values to remove. eps = 10^(-3); We count the number of singular values that were zero. nulls = Select[S,Max[Abs[#]]<eps&]; We now get the degree of our approximate gcd as below (for the epsilon we chose, it is 12). k = Length[sylmat] - Length[nulls] To form that gcd we'll row reduce a certain matrix and use the final row. newmat = Take[newS.Transpose[V], k]; coeffrow = Last[RowReduce[newmat]]; gcdapprox = coeffrow . z^Reverse[Range[0,Length[coeffrow]-1]] So here is what I get. In[150]:= InputForm[gcdapprox] Out[150]//InputForm= 6.980499621431014 + 0.000585041745605589*z + 0.0035109293315991083*z^2 + 0.002602929711965681*z^3 + 33.10928622044708*z^4 + 0.0019204205746833373*z^5 + 0.01152459048089943*z^6 + 0.008544087788492293*z^7 - 3.49084586854123*z^8 + 0.00011025782928825875*z^9 + 0.0006624721373616564*z^10 + 0.0004911544753497421*z^11 + z^12 I'll note that an epsilon of around 10^(-1) would have given a degree 16 approx gcd (not equal to a scalar multiple of p1, though). One can get an idea of how close are the various roots by looking at them as below. s1 = Sort[z /. NSolve[p1==0, z]]; s2 = Sort[z /. NSolve[p2==0, z]]; sg = Sort[z /. NSolve[gcdapprox==0, z]]; Playing tricks with Intersection and a comparison criterion might automate this process. Note that there are other approaches to approximate univariate gcds using QR decomposition, remainder sequences, approximate factoring, and probably some others I am forgetting. Daniel Lichtblau Wolfram Research

**References**:**Approximate/asymptotic factorization***From:*Paul Abbott <paul@physics.uwa.edu.au>