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$ Commented Jun 12, 2017 at 23:49
  • $\begingroup$ Closely related Finding the coefficient of a certain power in a generating function $\endgroup$ Commented Jun 12, 2017 at 23:56
  • $\begingroup$ @DavidG.Stork How else would you find the coefficients? $\endgroup$ Commented 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$ Commented 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$ Commented 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$ Commented 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$ Commented 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$ Commented Jun 19, 2017 at 5:01

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.