[Date Index] [Thread Index] [Author Index]

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>