Re: Approximate/asymptotic factorization
- To: mathgroup at smc.vnet.net
- Subject: [mg73676] Re: Approximate/asymptotic factorization
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sat, 24 Feb 2007 02:20:55 -0500 (EST)
- References: <200702210635.BAA17453@smc.vnet.net> <ermdo5$i0d$1@smc.vnet.net> <paul-2CDC4B.22092823022007@news.uwa.edu.au>
Paul Abbott wrote: > In article <ermdo5$i0d$1 at smc.vnet.net>, > Daniel Lichtblau <danl at wolfram.com> wrote: > > >>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. > > > OK, I've got that ref and also got > > N. Karmarkar and Y. N. Lakshman "Approximate polynomial greatest > common divisors and nearest singular polynomials", Proceedings of the > 1996 international symposium on Symbolic and algebraic computation, > 35 - 39, 1996 Right. I think theirs is more the remainder sequence approach. I also looked at it a bit yesterday, but decided 4M was simpler to demo. >>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}]] >> ] > > > I note that > > 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}]] > > can be replaced by the more intelligible (to me) > > Join[ > NestList[RotateRight, PadRight[coeffs1, l1 + l2 - 2], l2 - 2], > NestList[RotateRight, PadRight[coeffs2, l1 + l2 - 2], l1 - 2] > ] > > which is a direct implementation of the description at MathWorld. I also thought the code looked a bit bulky but did not really have time to consider how best to recode it. > >>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]; > > > Perhaps Decompose will also be useful here, since both p1 are > polynomials in z^4 (which explains the observed 4-fold symmetry). I had missed that. Yes, deflating might well be helpful in terms of improving stability of the various algorithms. >>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&]; > > > What is S and where is sylmat defined? I assume that > > sylmat = SylvesterMatrix[p1, p2, z] > > S = SingularValueList[sylmat] Clearly I was in too much of a rush. Yes, sylmat is the Sylvester matrix, and I used full SVD to get {U,S,V}. Sorry for the omission in the response. Also newS is a transformation of S to remove small singular values. Actually it was not needed. I first used it with default SingularValueDecomposition, then changed the latter to use Tolerance->eps which automates stripping the small values. Below is the code in question. Again, it is a sort of agglomeration of tactics and some is not necessary. sylmat = SylvesterMatrix[p1, p2, z]; eps = 10^(-3); {U,S,V} = SingularValueDecomposition[N[sylmat], Tolerance->eps]; newS = S /. aa_/;Abs[aa]<eps :> 0; nulls = Select[S,Max[Abs[#]]<eps&]; k = Length[sylmat] - Length[nulls] newmat = Take[newS.Transpose[V], k]; I think newS is the same as S now, since I use Tolerance->eps. So probably this aspect of the code could be made more concise. >>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]; > > > What is newS? > > >>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. > > > Thanks for this -- all very useful and interesting. > > 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 Looks like I had several omissions and cross-purpose code hiccups. Sorry you (PA) had to be the proof reader. I saw in my other MathGroup post yesterday that I made a similar mistake there as well, in the sense that what I wrote did not accurately reflect the final version of code I used (that one was a relatively tame wording issue). I just hope my bug fixing yesterday was more methodical. Daniel Lichtblau Wolfram Research
- References:
- Approximate/asymptotic factorization
- From: Paul Abbott <paul@physics.uwa.edu.au>
- Approximate/asymptotic factorization