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 *)

(* equation of the form z[[1]]^degree+... *)

(* compiled function to generate a random point on a (2*d-1)-dimensional sphere *)
(* express as a complex point in C^d *)
(* generate a 2*d-vector (a point) in R^(2*d) using rotationally symmetric normal distribution *)

(* normalise vector to 1, giving a point on sphere S^(2*d-1) *)

(* 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]],

(* Function to generate degree # of points by intersecting the line (p + tq) with Q=0 *)
(* t is variable that we solve for *)
(* 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 *)

(* function to generate at least N points that lie on Q = 0 *)

(* {0.0004921`,{5}} *)
(* {0.0013839`,{4,5}} *)

(* timing for at least 1000000 points *)
(* {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*)

(*equation of the form z[[1]]^degree+...*)




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;


findPoints$test[dim, degree, 1000000]//Dimensions//AbsoluteTiming
(* {72.7542245`,{1000000,5}} *)
  • $\begingroup$ Does it improve if you set "CompilationTarget"->"C"? $\endgroup$ Oct 8, 2020 at 20:24
  • 3
    $\begingroup$ Last code takes 9.41695 s on my computer with Parallelize[ , Method -> "CoarsestGrained"] $\endgroup$ Oct 13, 2020 at 21:20
  • 1
    $\begingroup$ Try to use it as follows findPoints$test[d_, deg_, N_] := Parallelize[(), Method -> "CoarsestGrained"] $\endgroup$ Oct 13, 2020 at 22:16
  • 2
    $\begingroup$ 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? $\endgroup$
    – Carl Woll
    Oct 15, 2020 at 0:27
  • 1
    $\begingroup$ 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. $\endgroup$ Oct 15, 2020 at 16:48

2 Answers 2


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

  Table[findPoints$test[dim, degree, 1000000/nk, i] // Dimensions, {i,
     nk}], Method -> "FinestGrained"] // AbsoluteTiming

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

With other options we have same result, for example,

  Table[findPoints$test[dim, degree, 1000000/nk, i] // Dimensions, {i,
     nk}], Method -> Automatic] // AbsoluteTiming

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

This result found out with 50% CPU. With 100% CPU it takes 10 s and so on.



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

Your Answer

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

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