# Speed improvement for root finding

I'm wondering whether it is possible to speed up the following code that I'm using to find points on the intersection between a line and the zero locus of an equation.

The idea is the following. First, generate two points $$\vec{p}$$ and $$\vec{q}$$ randomly on the 9-dimensional sphere. Each of these points is given as a 5-vector of complex numbers (viewing 10 real coordinates as 5 complex coordinates). Using these two points, you can parametrise a line in 5 complex dimensions as $$\vec{l}(t) = \vec{p} + t\, \vec{q},$$ where $$t$$ parametrises where on the line you are.

I then want to find the intersection of this line with the equation $$Q(\vec{z}) \equiv z_1^4+\ldots+z_5^4 = 0 ,$$ where $$(z_1,\ldots,z_5)$$ are again complex coordinates. (In the full problem, $$Q(z)$$ is actually of arbitrary degree and so I cannot find the roots of $$Q=0$$ analytically.) The intersection is given by $$Q(\vec{p} + t\, \vec{q})=0,$$ which is then a quartic equation for $$t$$. There are generically 4 solutions $$t_i$$ to such an equation. Substituting these values of $$t$$ back into $$\vec{l}(t)$$, we find 4 points (4 sets of coordinates $$\vec{z}_i=\vec{l}(t_i)$$) where the line intersects $$Q=0$$.

I then repeat this using many randomly generated points $$\vec{p}$$ and $$\vec{q}$$ to get a large number of points (usually on the order of 5-10 million).

My problem is that even for 1M points, this takes roughly 90s or so, which is a good chunk of the total runtime of my calculation. At the moment, I have a compiled function that finds the random points on the sphere (by sampling a normal distribution and then scaling the length of the vector to 1), a non-compiled function that does the root finding, and then a final ParallelTable function that repeats this for as many points as I want (and converts the result to a packed array, since I then do some numerical linear algebra with it).

I'd be very interested in any way to speed this up. I come back to this every month or so, and fail at squeezing anymore speed out of it.

I was hoping to get some improvement from compiling the root-finding function as this seems to be the slowest part, but couldn't see a way to do this. I was also wondering if generating all of the random points on the sphere first and then using Listable in some capacity might help, but again I couldn't quite see how to get this to work.

A nagging voice in the back of my head says I should just use C instead, but I've been amazed at how close Mathematica can usually get, so I'm not giving up yet! Thanks for the taking the time to read this!

The code is:

(* define equation we want to solve for Q=0 *)
(* quartic equation so generically 4 roots *)
degree=4;
dim=5;

(* equation of the form z[[1]]^degree+... *)
Q[z_]:=Sum[z[[i]]^degree,{i,1,dim}];

(* compiled function to generate a random point on a (2*d-1)-dimensional sphere *)
(* express as a complex point in C^d *)
genPoint$S=Compile[{{d,_Integer}}, Module[{x}, (* generate a 2*d-vector (a point) in R^(2*d) using rotationally symmetric normal distribution *) x=RandomVariate[NormalDistribution[],2d]; (* normalise vector to 1, giving a point on sphere S^(2*d-1) *) x=x/Norm[x]; (* convert real 10-vector to complex 5-vector (view as point in C^dim) *) Part[x,1;;d]+ I Part[x,d+1;;2 d]], "RuntimeOptions"->"Speed",CompilationOptions->{"InlineExternalDefinitions"->True},RuntimeAttributes->Listable,Parallelization->True]; (* Function to generate degree # of points by intersecting the line (p + tq) with Q=0 *) (* t is variable that we solve for *) genPoint[d_]:=Module[{t,line}, (* define line as (p + tq) where p and q are random points on S^(2*dim - 1) written as complex 5-vectors *) line=genPoint$$S[d]+t genPoint$$S[d]; (* solve for t in Q(p+tq)=0 - find degree=4 solutions as Q is quartic *) (* substitute solutions back into (p + tq) to find points *) (* get 4 x dim array as output - 4 sets of points, each specified by a complex 5-vector *) line/.{NRoots[Q[line]==0,t,Method->"JenkinsTraub"]//ToRules}] (* function to generate at least N points that lie on Q = 0 *) findPoints[d_,deg_,N_]:=Module[{output}, output=DeveloperToPackedArray[Flatten[ParallelTable[genPoint[d],{i,1,Ceiling[N/deg]}],1]]; output] genPoint$S[dim]//Dimensions//AbsoluteTiming
(* {0.0004921,{5}} *)
genPoint[dim]//Dimensions//AbsoluteTiming
(* {0.0013839,{4,5}} *)

(* timing for at least 1000000 points *)
findPoints[dim,degree,1000000]//Dimensions//AbsoluteTiming
(* {90.7323466,{1000000,5}} *)


Edit: I've managed to find a ~20% improvement by batching the computation of the points on the sphere. It also produces a packed array by default. The new code is simply:

(*define equation we want to solve for Q=0*)
(*quartic equation so generically 4 roots*)
degree=4;
dim=5;

(*equation of the form z[[1]]^degree+...*)
Q[z_]:=Sum[z[[i]]^degree,{i,1,dim}];

findPoints$test[d_,deg_,N_]:=( x=RandomVariate[NormalDistribution[],{Ceiling[N/deg],2d}]; y=RandomVariate[NormalDistribution[],{Ceiling[N/deg],2d}]; x=x/(Norm/@x); y=y/(Norm/@y); cx=Part[x, All, 1 ;; d] + I Part[x, All, 1 + d ;; 2 d]; cy=Part[y, All, 1 ;; d] + I Part[y, All, 1 + d ;; 2 d]; lines=cx + t cy; Flatten[(#/.{NRoots[Q[#]==0,t,Method->"JenkinsTraub"]//ToRules})&/@lines,{{1,2},{3}}]) findPoints$test[dim, degree, 1000000]//Dimensions//AbsoluteTiming
(* {72.7542245,{1000000,5}} *)

• Does it improve if you set "CompilationTarget"->"C"? Oct 8, 2020 at 20:24
• Last code takes 9.41695 s on my computer with Parallelize[ , Method -> "CoarsestGrained"] Oct 13, 2020 at 21:20
• Try to use it as follows findPoints$test[d_, deg_, N_] := Parallelize[(), Method -> "CoarsestGrained"] Oct 13, 2020 at 22:16 • You said in a comment to the deleted answer that you were particularly interested in a Q that is degree 5 in the z. Of course degree 5 polynomials don't have a solution in radicals. Degree 4 polynomials do, and it is possible to vectorize root finding of degree 4 polynomials, achieving a large speed up. Of course this approach won't work at all for degree 5 polynomials. So, is there any benefit to providing an answer using this appoach? Oct 15, 2020 at 0:27 • Hm. Numerical root finding typically profits from good initial guesses. One could exploit this as follows: Rewrite the random point generation as a Markov chain with decently small increments. Then you can use the solutions of the old point pair as initial guesses for the new point pair and perform a few Newton iterations. If Newton's method converges: Fine. Otherwise, use NRoots as fallback. If the Markov chain is sufficiently long, you still get the product of sphere uniformuly covered. Oct 15, 2020 at 16:48 ## 2 Answers This code devoted for testing to demonstrate effect of parallel job. Let take number of kernels used equal nk, then we can distribute computation as follows SeedRandom[1234](*we set randomizer for testing only!*); degree = 4; dim = 5; LaunchKernels[]; nk =$KernelCount;

Q[z_] := Sum[z[[i]]^degree, {i, 1, dim}];

findPoints$test[d_, deg_, N_, i_] := (x = RandomVariate[NormalDistribution[], {Ceiling[N/deg], 2 d}]; y = RandomVariate[NormalDistribution[], {Ceiling[N/deg], 2 d}]; x = x/(Norm /@ x); y = y/(Norm /@ y); cx = Part[x, All, 1 ;; d] + I Part[x, All, 1 + d ;; 2 d]; cy = Part[y, All, 1 ;; d] + I Part[y, All, 1 + d ;; 2 d]; lines = cx + t cy; root[i] = Flatten[(# /. {NRoots[Q[#] == 0, t, Method -> "JenkinsTraub"] // ToRules}) & /@ lines, {{1, 2}, {3}}]);  In the case with nk=4 we have Parallelize[ Table[findPoints$test[dim, degree, 1000000/nk, i] // Dimensions, {i,
nk}], Method -> "FinestGrained"] // AbsoluteTiming

Out[]= {19.8538, {{250000, 5}, {250000, 5}, {250000, 5}, {250000,
5}}}


With other options we have same result, for example,

Parallelize[
Table[findPoints$test[dim, degree, 1000000/nk, i] // Dimensions, {i, nk}], Method -> Automatic] // AbsoluteTiming Out[]= {20.8837, {{250000, 5}, {250000, 5}, {250000, 5}, {250000, 5}}}  This result found out with 50% CPU. With 100% CPU it takes 10 s and so on. Edit Tried an anylytical solution as shown below. Sorry, but didn't see a programming error. So it turned out that it is much slower. I show it anyway. Learn, how anylytical solutions look like ... p = {p1, p2, p3, p4, p5}; q = {q1, q2, q3, q4, q5}; Q[z_] := Total@(z^4) sol[{p1_, p2_, p3_, p4_, p5_}, {q1_, q2_, q3_, q4_, q5_}] = t /. Solve[Q[p + t q] == 0, t, Quartics -> False]  Solution writen in vector form is Clear[p, q, p1, p2, p3, p4, p5, q1, q2, q3, q4, q5] sol4 = {Root[p .p^3 + (4 p^3 .q) #1 + (6 p^2 .q^2) #1^2 + (4 p.q^3) #1^3 + (q .q^3) #1^4 &, 1], Root[p .p^3 + (4 p^3 .q) #1 + (6 p^2 .q^2) #1^2 + (4 p.q^3) #1^3 + (q .q^3) #1^4 &, 2], Root[p .p^3 + (4 p^3 .q) #1 + (6 p^2 .q^2) #1^2 + (4 p.q^3) #1^3 + (q .q^3) #1^4 &, 3], Root[p .p^3 + (4 p^3 .q) #1 + (6 p^2 .q^2) #1^2 + (4 p.q^3) #1^3 + (q .q^3) #1^4 &, 4] };  Analysing gives Roots as a Binomial sum of p and q vector. rsum[d_] := p.p^(d - 1) + Sum[(Binomial[d, j] p^(d - j) . q^j ) #1^j, {j, 1, d - 1}] + q.q^(d - 1) #1^d sol[d_] := Table[Root[Evaluate[rsum[d]] &, k], {k, 1, d}] sol[4] == sol4 (* True *)  Your function to generate points genPoint$S =
Compile[{{d, _Integer}},
Module[{x},(*generate a 2*d-vector (a point) in R^(2*
d) using rotationally symmetric normal distribution*)
x = RandomVariate[NormalDistribution[], 2 d];
(*normalise vector to 1,giving a point on sphere S^(2*d-1)*)
x = x/Norm[x];
(*convert real 10-vector to complex 5-vector (view as point in  C^
dim)*)Part[x, 1 ;; d] + I Part[x, d + 1 ;; 2 d]],
"RuntimeOptions" -> "Speed",
CompilationOptions -> {"InlineExternalDefinitions" -> True},
RuntimeAttributes -> Listable, Parallelization -> True]

{pp, qq} = {genPoint$$S[5], genPoint$$S[5]}

ss[p_, q_] = Evaluate[sol[6]];
pp + # qq & /@ ss[pp, qq] // TableForm


Here dimension 6

(ss[p_, q_] = Evaluate[sol[6]];
ParallelTable[({pp, qq} = {genPoint$$S[5], genPoint$$S[5]};
pp + # qq & /@ ss[pp, qq]), {10000}];) // AbsoluteTiming
`