Speeding the calculation of a complex polynomial from its roots

Question

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?

Answer 1

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}

Answer 2

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

Answer 3

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.