# Does NRoots own an abstract counterpart? If not, can we write one?

We know when solving linear algebra equations, despite its abstract syntax, LinearSolve is much faster compared to Solve:

n = 1000;
m = DiagonalMatrix@RandomReal[{1, 2}, n];
b = RandomReal[1, n];
sol1 = LinearSolve[m, b]; // AbsoluteTiming
x = Table[Unique[], {n}];
sol2 = Solve[m.x == b, x]; // AbsoluteTiming
Last /@ First@sol2 == sol1

{0.124800, Null}
{1.216800, Null}
True


I wonder if NRoots owns such a counterpart, too?

…To be honest, I came across this problem because a friend of mine claimed that a matlab function roots() is one order of magnitude faster than NRoots[]. The following is the comparison between NRoots and roots.

NRoots in Mathematica 9.0.1:

n = 10; m = 10000; list = RandomInteger[{1, 100}, {m, n + 1}];
equations = list.x^Range[0, n];
NRoots[#, x] & /@ equations; // AbsoluteTiming


{7.984318, Null}

roots in Matlab R2014a:

 clear;clc;
n = 10;
m = 10000;
lst=randi(100,m,n+1);
tic;
for i=1:m
roots(lst(i,:));
end
toc


Elapsed time is 0.756003 seconds.

Given the performance difference between LinearSolve and Solve, I believe if the analysis for the structure of polynomial is eliminated, Mathematica will have the same speed.

If such an internal function doesn't exist, can we have a self-made one that eliminates the preprocessing? Anyway, though I almost don't know the syntax of Matlab, the definition of roots() seems to be simple so I think it won't be that hard to implement the same algorithm in Mathematica.

OK, since the target of this question hasn't been completely achieved (NRoots owns 3 methods, currently only "CompanionMatrix" is implemented), let me add something to try to draw more attention. You remember this?:

It's said that Sam Derbyshire spent 4 days for a similar one, while the above image only takes me no more than 3 minutes in my old 2GHz laptop with roots[] that chyaong showed in his answer below:

(* I modified chyaong's roots[] a little for this specific task. *)
modifiedroots[c_List] :=
Block[{a = DiagonalMatrix[ConstantArray[1., Length@c - 1], -1]},
a[[1]] = -c;
Eigenvalues[a]];

s = 1000;
l = ConstantArray[0., {s + 1, s + 1}];

l[[#, #2]] += 1. & @@@
Round[1 + s Rescale@
Function[z, {Im@z, Re@z}, Listable]@
Flatten[modifiedroots /@ Tuples[{-1., 1.}, 18]]]; // AbsoluteTiming

ArrayPlot[UnitStep[80 - l] l, ColorFunction -> "AvocadoColors"]

{161.737000, Null}


Let's ignore the fact that my image is only degree 19.

This is not the end, according to my simple test (I simply add different Methods to the 2nd sample of this post), the "Aberth" method is probably more efficient than the CompanionMatrix inside NRoots, so if it's implemented in an abstract form, the above image will likely to be done in a even shorter time!

• (1) NRoots takes a Method option. The three documented possibilities are "Aberth", "JenkinsTraub", and "CompanionMatrix". (I do not know if there are undocumented possibilities). From what I could see, the third of those possibilities should be similar to what roots() does. – Daniel Lichtblau Feb 5 '15 at 16:16
• (2) I confess I fail to see an analogy to Solve/LinearSolve, unless NRoots plays the role of LinearSolve with Roots being the "more general" function (insofar as it handles exact and symbolic input). I suppose one could go from the coefficient list and directly form a companion matrix and invoke Eigenvalues. Would that be faster or as reliable? I don't know. – Daniel Lichtblau Feb 5 '15 at 16:18
• To the downvoter, I am interested in what's wrong with my question, would you please elaborate. I'm not trying to complain here, I just want improve my question if possible. – xzczd Feb 6 '15 at 12:23
• I don't know exactly (I didn't vote one way or the other) but I can say a bit about things I thought could bear improvement. One issue is that it could have been phrased differently, maybe as "Are there different methods of NRoots that could have impact on performance vs. quality?" As I noted in a comment, the analogy to Solve/LinearSolve is just not clear to me. A bigger issue, in my mind, is that there is no clear example where NRoots is seen to perform relatively badly as compared to roots(). – Daniel Lichtblau Feb 6 '15 at 16:45
• @DanielLichtblau That sounds reasonable. Question edited. (Thanks for the test from Amita!) – xzczd Feb 9 '15 at 9:16

Translated that roots() function from Matlab code to Mathematica, about 4 times faster than NRoots .

Clear["*"];
n=2000;
m=RandomReal[1,{n,10}];
res1=x/.(ToRules@NRoots[FromDigits[#,x]==0.,x]& /@ m);//AbsoluteTiming
(*-------------------------------*)
roots[c_List]:=Block[{a},
a=DiagonalMatrix[ConstantArray[1.,Length@c-2],-1];
a[[1]]=-Rest@c/First@c;
Eigenvalues[a]
];

res2=Join @@ roots /@ m;//AbsoluteTiming

Union@Chop[Sort@res1 - Sort@res2]


Output:

{1.080062, Null}

{0.241012, Null}

{0}

n = 10;
m = 10000;
list = RandomReal[{1, 100}, {m, n + 1}];
equations = list.x^Range[0, n];
res1 = x /. (ToRules@NRoots[#, x] & /@ equations); // AbsoluteTiming
res2 = Join @@ roots /@ Reverse /@ list; // AbsoluteTiming
Sort@res1 - Sort@res2 // Chop // Union


• If you try this with modestly higher degree (so preprocessing time does not dominate) you will see comparable timings. E.g with n = 100; m = RandomReal[1, {n, 100}]; Also the residuals from the NRoots result will be smaller. – Daniel Lichtblau Feb 8 '15 at 20:19
• @DanielLichtblau Yeah, and eliminating the preprocessing is the target of this question. I just edited the question and made this part clearer (I think). – xzczd Feb 9 '15 at 9:18
• Ah. No, I'm not aware of a lower level function that extracts numeric roots from a list of polynomial equations. – Daniel Lichtblau Feb 9 '15 at 16:18
• As for the higher degree case I had mentioned, here is a plausible test. n = 100; m = 100; list = RandomReal[{1, 100}, {m, n + 1}]; equations = list.x^Range[0, n]; res1 = x /. (ToRules@NRoots[#, x] & /@ equations); // AbsoluteTiming res2 = Join @@ roots /@ Reverse /@ list; // AbsoluteTimingl – Daniel Lichtblau Feb 9 '15 at 16:18
• Any plan to implement the other 2 method of NRoots[]? Though I admit this question is motivated by my competitive mood against Matlab, now its target is to get the counterpart of NRoots[]` anyway :) – xzczd Feb 10 '15 at 10:55