3
$\begingroup$

Suppose I have a list of (complex, in my case) numbers, and I want to construct the coefficients of the monic polynomial of which they are the roots. The obvious way:

poly[n_]:= CoefficientList[Expand[Product[x-list[[i]],{i, 1, n}]], x]

Seems to take quartic time in $n,$ and for $n=5000$ takes about a minute on my MacBook Pro.

The less obvious way:

poly2[n_]:= CoefficientList[CharacteristicPolynomial[DiagonalMatrix[list], x], x]

Is much faster (takes 8 seconds for the same computation) but it's hard to believe that this is optimal. Any ideas?

Edit There is also InterpolatingPolynomial, but that is slower than the CharacteristicPolynomial scheme.

More interesting edit Upon meditating on this for a while, I came up with the following:

cfromr[roots_] := 
 With[{n = Length[roots], 
   f = Compile[{{x, _Complex}}, 
     Product[(x - roots[[i]]), {i, 1, Length[roots]}]]},
  Module[{fvals = Array[f[Exp[(# - 1) 2 Pi I/(n + 1)]] &, n + 1]},
   InverseFourier[fvals]/Sqrt[n + 1]]]

Let's try it:

cfromr[{1}]
{-1., 1.}

cfromr[{1, 1}]
{1., -2., 1.}

cfromr[{1, 2, 3}]
{-6., 11., -6., 1.}

Cool! And of course, it is fast as blazes (orders of magnitude faster than the previous efforts - try it). But, let's just try it for something a little more complicated, and compare it with the boneheaded approach below:

stupidCfromR[roots_] := 
 CoefficientList[Product[x - roots[[i]], {i, 1, Length[roots]}], x]

OK, 

foo = RandomReal[{0, 1}, 100]
{0.615659, 0.749481, 0.877498, 0.179083, 0.580726, 0.297616, 0.436366, \
0.885845, 0.207169, 0.979251, 0.768706, 0.962705, 0.934818, 0.558659, \
0.810448, 0.578337, 0.880773, 0.0389867, 0.561441, 0.636913, \
0.747634, 0.061617, 0.686383, 0.149683, 0.413282, 0.659633, 0.534624, \
0.266132, 0.16876, 0.29776, 0.156488, 0.678895, 0.287585, 0.192998, \
0.527663, 0.352505, 0.718803, 0.380109, 0.44967, 0.46576, 0.693758, \
0.557772, 0.811107, 0.824323, 0.77316, 0.365756, 0.283581, 0.849779, \
0.0787828, 0.180125, 0.0536948, 0.406494, 0.570083, 0.458642, \
0.209918, 0.0254337, 0.340365, 0.280486, 0.0969694, 0.902567, \
0.311759, 0.445944, 0.734187, 0.556372, 0.968549, 0.722562, 0.885184, \
0.642906, 0.264817, 0.541621, 0.793586, 0.259047, 0.0489964, \
0.183852, 0.499949, 0.0692111, 0.714986, 0.793983, 0.525797, \
0.508116, 0.0739616, 0.0726041, 0.542157, 0.143221, 0.980601, \
0.994549, 0.542997, 0.77058, 0.706582, 0.858526, 0.00312537, \
0.0235296, 0.961969, 0.436431, 0.611778, 0.167991, 0.424293, \
0.0989057, 0.256478, 0.732805}

stupidCfromR[foo]
{6.56546*10^-45, -5.01907*10^-42, 1.56046*10^-39, -2.8729*10^-37, 
 3.64875*10^-35, -3.47096*10^-33, 2.60308*10^-31, -1.59376*10^-29, 
 8.17047*10^-28, -3.57511*10^-26, 1.35545*10^-24, -4.50701*10^-23, 
 1.32751*10^-21, -3.49263*10^-20, 8.2664*10^-19, -1.77082*10^-17, 
 3.4517*10^-16, -6.15044*10^-15, 1.00596*10^-13, -1.5158*10^-12, 
 2.11112*10^-11, -2.72564*10^-10, 3.27087*10^-9, -3.65718*10^-8, 
 3.81829*10^-7, -3.72996*10^-6, 0.0000341544, -0.00029365, \
0.00237425, -0.0180782, 0.129804, -0.879936, 5.63813, -34.1814, \
196.26, -1068.19, 5515.59, -27039.4, 125940., -557668., 
 2.34902*10^6, -9.41749*10^6, 3.5953*10^7, -1.30763*10^8, 
 4.53272*10^8, -1.49802*10^9, 4.72179*10^9, -1.41987*10^10, 
 4.07431*10^10, -1.11587*10^11, 2.91749*10^11, -7.28283*10^11, 
 1.73593*10^12, -3.9513*10^12, 8.58888*10^12, -1.78288*10^13, 
 3.53415*10^13, -6.68947*10^13, 1.20892*10^14, -2.08562*10^14, 
 3.43421*10^14, -5.39605*10^14, 8.08845*10^14, -1.15627*10^15, 
 1.57581*10^15, -2.04655*10^15, 2.53174*10^15, -2.98174*10^15, 
 3.34142*10^15, -3.56063*10^15, 3.60542*10^15, -3.46646*10^15, 
 3.1619*10^15, -2.73363*10^15, 2.23776*10^15, -1.73254*10^15, 
 1.26708*10^15, -8.74157*10^14, 5.68038*10^14, -3.47091*10^14, 
 1.9906*10^14, -1.06931*10^14, 5.36788*10^13, -2.5117*10^13, 
 1.09229*10^13, -4.40038*10^12, 1.63609*10^12, -5.5903*10^11, 
 1.74674*10^11, -4.96231*10^10, 1.27307*10^10, -2.92548*10^9, 
  5.96249*10^8, -1.06468*10^8, 
 1.63982*10^7, -2.13439*10^6, 228286., -19264.2, 1202.6, -49.3738, 1}

Looks reasonable. Next,

cfromr[foo]
{5.9802 + 0.0792079 I, -8.1594 + 0.0953864 I, 
10.3109 + 0.107405 I, -11.2502 - 0.144703 I, 
11.7326 + 0.0258995 I, -10.0821 + 0.0336182 I, 
8.24025 + 0.123999 I, -5.02346 + 0.598038 I, 0.450765 + 0.232466 I, 
4.02386 - 0.0386024 I, -10.1325 - 0.354894 I, 
15.4001 + 0.0682209 I, -20.8642 + 0.00860923 I, 
25.4178 + 0.0414025 I, -29.7352 - 0.0478316 I, 
32.0813 + 0.0872942 I, -34.4719 - 0.186356 I, 
34.9878 - 0.140137 I, -34.796 - 0.438696 I, 
34.1201 - 0.266264 I, -32.825 - 0.34352 I, 
31.1676 + 0.0377669 I, -28.6403 - 0.273897 I, 
26.9552 - 0.226578 I, -24.5299 - 0.210857 I, 
21.7318 - 0.166584 I, -18.3817 - 0.406169 I, 
14.7006 - 0.211768 I, -10.6488 + 0.0218902 I, 
5.89532 + 0.0502537 I, -2.43193 + 0.145602 I, -2.29822 + 0.0280126 I,
10.1683 + 0.0739904 I, -41.196 + 0.297882 I, 
204.612 + 0.22332 I, -1077.32 + 0.356418 I, 
5524.19 - 0.0830177 I, -27046.6 - 0.241576 I, 
125945. - 0.265022 I, -557670. + 0.314207 I, 
2.34902*10^6 + 0.110874 I, -9.41749*10^6 + 0.0519993 I, 
3.5953*10^7 - 0.0733092 I, -1.30763*10^8 - 0.0274269 I, 
4.53272*10^8 + 0.136665 I, -1.49802*10^9 + 0.139398 I, 
4.72179*10^9 - 0.00431823 I, -1.41987*10^10 - 0.0238292 I, 
4.07431*10^10 - 0.0859524 I, -1.11587*10^11 + 0.413876 I, 
2.91749*10^11 - 0.124791 I, -7.28283*10^11 + 0.0234331 I, 
1.73593*10^12 - 0.0233547 I, -3.9513*10^12 + 0.137854 I, 
8.58888*10^12 - 0.0113619 I, -1.78288*10^13 + 0.328116 I, 
3.53415*10^13 - 0.310418 I, -6.68947*10^13 + 0.0321282 I, 
1.20892*10^14 - 0.197183 I, -2.08562*10^14 + 0.314149 I, 
3.43421*10^14 - 0.127756 I, -5.39605*10^14 + 0.0133408 I, 
8.08845*10^14 + 0.0230662 I, -1.15627*10^15 - 0.272123 I, 
1.57581*10^15 - 0.192791 I, -2.04655*10^15 - 0.0992356 I, 
2.53174*10^15 + 0.342712 I, -2.98174*10^15 - 0.448979 I, 
3.34142*10^15 - 0.382704 I, -3.56063*10^15 - 0.161659 I, 
3.60542*10^15 - 0.245491 I, -3.46646*10^15 + 0.181787 I, 
3.1619*10^15 + 0.53769 I, -2.73363*10^15 + 0.379394 I, 
2.23776*10^15 + 0.0504824 I, -1.73254*10^15 + 0.3809 I, 
1.26708*10^15 + 0.062272 I, -8.74157*10^14 + 0.230549 I, 
5.68038*10^14 - 0.115016 I, -3.47091*10^14 + 0.317305 I, 
1.9906*10^14 - 0.0884983 I, -1.06931*10^14 - 0.00104238 I, 
5.36788*10^13 - 0.25693 I, -2.5117*10^13 + 0.218549 I, 
1.09229*10^13 + 0.26878 I, -4.40038*10^12 - 0.140603 I, 
1.63609*10^12 + 0.0812383 I, -5.5903*10^11 + 0.14043 I, 
1.74674*10^11 - 0.402052 I, -4.96231*10^10 + 0.272428 I, 
1.27307*10^10 - 0.0000964737 I, -2.92548*10^9 - 0.0669364 I, 
5.96249*10^8 - 0.168353 I, -1.06468*10^8 + 0.42484 I, 
1.63982*10^7 - 0.511921 I, -2.13439*10^6 - 0.268171 I, 
228291. - 0.242538 I, -19266.8 + 0.125755 I, 
1203.94 + 0.0223555 I, -48.8476 + 0.0424496 I, -1.905 - 0.0992187 I}

Notice that this is complete nonsense. The numbers are not real, and their magnitude is completely off. Either I am insane, or Mathematica is. Any suggestions?

Improve this question
$\endgroup$
9
  • $\begingroup$ Why the need to Expand in your first approach? $\endgroup$
    – David G. Stork
    Jun 12, 2017 at 23:49
  • $\begingroup$ Closely related Finding the coefficient of a certain power in a generating function $\endgroup$
    – Artes
    Jun 12, 2017 at 23:56
  • $\begingroup$ @DavidG.Stork How else would you find the coefficients? $\endgroup$
    – Igor Rivin
    Jun 13, 2017 at 0:05
  • $\begingroup$ @DavidG.Stork Ah, I see, the Expand[] is implicit, but the system clearly does it, because removing the Expand[] does not speed things up. $\endgroup$
    – Igor Rivin
    Jun 13, 2017 at 0:07
  • $\begingroup$ @Artes Related, but not that closely. Computing a single coefficient is quite different from computing all of them (the tradeoffs are quite different). $\endgroup$
    – Igor Rivin
    Jun 13, 2017 at 0:10

3 Answers 3

Reset to default
6
$\begingroup$

Multiplying two polynomials is equivalent to using ListConvolve on the coefficients of the polynomial. The only issue is that using ListConvolve will experience a catastrophic loss of precision. This can be mitigated by arbitrarily increasing the precision before using ListConvolve. Here is a sample set of roots:

roots = RandomComplex[1+I, 5000];

Here is essentially the obvious method from your OP:

res1 = CoefficientList[Times @@ (roots-x), x]; //AbsoluteTiming

{52.9876, Null}

Now, for the ListConvolve method. The first step is to use the obvious method in chunks of 100, since it is pretty fast for this size root list:

coeffList = Table[
    CoefficientList[Times @@ (k-x), x],
    {k, Partition[roots, 100]}
]; //AbsoluteTiming

{0.477353, Null}

Now, we use ListConvolve on the coefficient lists, but arbitrarily raising the precision first:

res2 = Fold[
    ListConvolve[SetPrecision[#2,100], SetPrecision[#1,100], {1,-1}, 0]&,
    coeffList
]; //AbsoluteTiming

{27.2074, Null}

Finally, the CharacteristicPolynomial method:

res3 = CoefficientList[
    CharacteristicPolynomial[DiagonalMatrix[roots], x],
    x
]; //AbsoluteTiming

{42.9382, Null}

Now, let's compare the outputs:

res1[[-10;;]]
res2[[-10;;]]//N
res3[[-10;;]]

{-1.59795*10^26 - 1.65294*10^26 I, 5.87868*10^23 + 8.83793*10^21 I, -9.57174*10^20 + 9.32321*10^20 I, 2.99454*10^16 - 2.65645*10^18 I, 3.17046*10^15 + 3.23058*10^15 I, -6.42539*10^12 - 4.82769*10^10 I, 5.18789*10^9 - 5.12975*10^9 I, -23330.2 + 6.21159*10^6 I, -2487.88 - 2497.24 I, 1}

{-1.59795*10^26 - 1.65294*10^26 I, 5.87868*10^23 + 8.83793*10^21 I, -9.57174*10^20 + 9.32321*10^20 I, 2.99454*10^16 - 2.65645*10^18 I, 3.17046*10^15 + 3.23058*10^15 I, -6.42539*10^12 - 4.82769*10^10 I, 5.18789*10^9 - 5.12975*10^9 I, -23330.2 + 6.21159*10^6 I, -2487.88 - 2497.24 I, 1.}

{-1.59795*10^26 - 1.65294*10^26 I, 5.87868*10^23 + 8.83793*10^21 I, -9.57174*10^20 + 9.32321*10^20 I, 2.99454*10^16 - 2.65645*10^18 I, 3.17046*10^15 + 3.23058*10^15 I, -6.42539*10^12 - 4.82769*10^10 I, 5.18789*10^9 - 5.12975*10^9 I, -23330.2 + 6.21159*10^6 I, -2487.88 - 2497.24 I, 1}

The last 10 coefficients are in reasonable agreement, and the first 10 are as well. However, the middle coefficients are not in agreement. For instance:

res1[[3000]]
res2[[3000]]//N
res3[[3000]]

-3.911349964088013*10^1154 + 6.655806866305701*10^1154 I

-7.24028240645896*10^1153 + 1.217768123179758*10^1154 I

-7.240282406458935*10^1153 + 1.217768123179755*10^1154 I

I'm inclined to believe the ListConvolve approach, because I can arbitrarily use even higher precision, and the output doesn't change:

res4 = Fold[
    ListConvolve[SetPrecision[#2,500], SetPrecision[#1,500], {1,-1}, 0]&,
    coeffList
]; //AbsoluteTiming

{31.0709, Null}

res4[[3000]]//N

-7.24028240645896*10^1153 + 1.217768123179758*10^1154 I

Finally, it is possible to tweak the ListConvolve approach by using different chunks. For example, here I use chunks of 500:

coeffList2 = Table[
    CoefficientList[Times @@ (k-x), x],
    {k, Partition[roots, 500]}
]; //AbsoluteTiming

{2.46872, Null}

res5 = Fold[
    ListConvolve[SetPrecision[#2,1500], SetPrecision[#1,1500], {1,-1}, 0]&,
    coeffList2
]; //AbsoluteTiming

{10.4263, Null}

Notice that this combination is more than twice as fast as using chunks of 100, but it requires insanely high precision to get outputs that agree with the 100 chunk approach:

res5[[3000]]//N

-7.24028240645890*10^1153 + 1.217768123179749*10^1154 I

res5[[-10;;]]//N

{-1.59795*10^26 - 1.65294*10^26 I, 5.87868*10^23 + 8.83793*10^21 I, -9.57174*10^20 + 9.32321*10^20 I, 2.99454*10^16 - 2.65645*10^18 I, 3.17046*10^15 + 3.23058*10^15 I, -6.42539*10^12 - 4.82769*10^10 I, 5.18789*10^9 - 5.12975*10^9 I, -23330.2 + 6.21159*10^6 I, -2487.88 - 2497.24 I, 1. + 0. I}

For example, a precision of 1000 is not enough:

res5 = Fold[
    ListConvolve[SetPrecision[#2,1000], SetPrecision[#1,1000], {1,-1}, 0]&,
    coeffList2
]; //AbsoluteTiming

res5[[-10;;]]//N

{7.1049, Null}

{0. + 0. I, 0. + 0. I, 0. + 0. I, 0. + 0. I, 0. + 0. I, 0. + 0. I, 0. + 0. I, 0. + 0. I, 0. + 0. I, 0. + 0. I}

Improve this answer
$\endgroup$
1
$\begingroup$

The coefficients of each term are sums of products of subsets of the roots with alternating sign. e.g (for an even number of roots)

coeffs = Subsets[roots, {#}] & /@ Range[Length[roots]];

plusminus = ConstantArray[{-1, 1}, Quotient[Length[coeffs], 2]] // Flatten

coefficientlist= (Total[Times @@@ #] & /@ coefs)*plusminus
Improve this answer
$\endgroup$
1
  • $\begingroup$ Yes, but this makes the computation EXPONENTIAL time. $\endgroup$
    – Igor Rivin
    Jun 13, 2017 at 0:53
1
$\begingroup$

Some (slight) savings are made by replacing Product by Times in your code

I would expect that Series would be most efficient for such a computation (it is more efficient than CharacteristicPolynomial and should be optimal), e.g.,

With[{p = Length[foo] + 1}, 
   PadLeft[(Times @@ (x - foo) + O[x]^p)[[3]], p]]

The PadLeft is required for zero roots.

Wilkinson's polynomial illustrates that the location of the roots can be very sensitive to perturbations in the coefficients of the polynomial.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Actually, the CharacteristicPolynomial method is faster for large root lists. $\endgroup$
    – Carl Woll
    Jun 14, 2017 at 16:00
  • $\begingroup$ Actually, computing the Series first and then applying CharacteristicPolynomial seems to be faster by almost a factor of 3 (at least for random integer coefficients) than just using CharacteristicPolynomial. However, for random reals, the times are comparable (tested with p=2000). $\endgroup$
    – TheDoctor
    Jun 19, 2017 at 5:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.