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=Developer`ToPackedArray[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}} *)
"CompilationTarget"->"C"? $\endgroup$Parallelize[ , Method -> "CoarsestGrained"]$\endgroup$findPoints$test[d_, deg_, N_] := Parallelize[(), Method -> "CoarsestGrained"]$\endgroup$