Speeding up cross and dot products for a list of vectors

Question

I am working on a code where I perform dot and cross product operations on a large list of vectors multiple times. I am using MapThread to achieve this but I feel the speed of operation is not up to the mark.

The operation I want to achieve is this, $f(\mathbf{S}_1,\mathbf{S}_2)=\frac{1}{a^2+b^2+c^2(1-(\mathbf{S}_1.\mathbf{S}_2)^2)+2 a b \mathbf{S}_1.\mathbf{S}_2}(a \mathbf{S}_1+b \mathbf{S}_2-c (\mathbf{S}_1 \times \mathbf{S}_2), a \mathbf{S}_2+b \mathbf{S}_1+c (\mathbf{S}_1 \times \mathbf{S}_2))$.

$\mathbf{S}_1$ and $\mathbf{S}_2$ are three dimensional vectors on the unit sphere and $a,b,c$ are arbitrary. What I need to compute is $f(\mathbf{\tilde{S}}_{2n-1},\mathbf{\tilde{S}}_{2n}), \hspace{0.2in} n=1...L/2$ where $\mathbf{\tilde{S}}$ is a list of normalized $3D$ vectors of length $L$.

Here is the sample code I have now,

N1=1000;
r1 = ArcCos[RandomReal[{-1, 1}, N1]]; 
r2 = RandomReal[{0, 2*Pi}, N1]; a=2;b=3;c=4;
spinsinit = Transpose[{Sin[r1]*Cos[r2], Sin[r1]*Sin[r2], Cos[r1]}]; spins = spinsinit; \[Tau] = 1; 
qq1 = Table[list1 = spins[[2*Range[N1/2] - 1]]; list2 = spins[[2*Range[N1/2]]]; list3 = MapThread[Dot, {list1, list2}]; list31 = MapThread[Cross, {list1, list2}]; list32=Sqrt[a^2+b^2+c^2 (1-list3^2)+2 a b list3
];     list4 = (a list1 + b list2 + c list31)/(list32);   list5 = (a list2 + b list1 - c list31)/(list32);spins = Normalize /@ Flatten[Transpose[Join[{list4, list5}]], 1],{i,1,100}];

In context to this code, I have two questions.

1. It seems I have to renormalize the vectors after some iterations because somehow numerical errors creep in and the results blow up, is there anyway more efficient way to tackle this phenomena?

2. I don't think this code is fully optimized, the sample code takes $\sim 1-2$ seconds to run, and I need to repeat this operation for around $10^6$ times at least, it just does not seem feasible. So any improvements will be greatly appreciated.

Answer 1

Almost the fastest way would be to use Compile to generate a CompiledFunction.

randomSpherePoint[n_] := Module[{z, \[Phi], r},
   \[Phi] = RandomReal[{0, 2 Pi}, n];
   z = RandomReal[{-1, 1}, n];
   r = Sqrt[1. - z z];
   Transpose[{r Cos[\[Phi]], r Sin[\[Phi]], z}]
   ];

cf = Compile[{{x, _Real, 1}, {y, _Real, 1}, {a, _Real}, {b, _Real}, {c, _Real}, {iters, _Integer}},
   Block[{x1, x2, x3, y1, y2, y3, u1, u2, u3, v1, v2, v3, ufactor, 
     vfactor},
    
    x1 = Compile`GetElement[x, 1];
    x2 = Compile`GetElement[x, 2];
    x3 = Compile`GetElement[x, 3];
    
    y1 = Compile`GetElement[y, 1];
    y2 = Compile`GetElement[y, 2];
    y3 = Compile`GetElement[y, 3];
    
    u1 = u2 = u3 = v1 = v2 = v3 = 0.;
    
    Table[
     u1 = a x1 + b y1 + c (-x3 y2 + x2 y3);
     u2 = a x2 + b y2 + c (x3 y1 - x1 y3);
     u3 = a x3 + c (-x2 y1 + x1 y2) + b y3;
     v1 = b x1 + a y1 - c (-x3 y2 + x2 y3);
     v2 = b x2 + a y2 - c (x3 y1 - x1 y3);
     v3 = b x3 - c (-x2 y1 + x1 y2) + a y3;
     ufactor = 1./Sqrt[u1 u1 + u2 u2 + u3 u3];
     vfactor = 1./Sqrt[v1 v1 + v2 v2 + v3 v3];
     
     x1 = u1 ufactor;
     x2 = u2 ufactor;
     x3 = u3 ufactor;
     
     y1 = v1 vfactor;
     y2 = v2 vfactor;
     y3 = v3 vfactor;
     
     {{x1, x2, x3}, {y1, y2, y3}}
     
     , {i, 1, iters}]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

Note that cf is a CompiledFunction with RuntimeAttributes -> {Listable} and Parallelization -> True. So it automatically threads and parallelizes. It has a bit different data layout for inputs and outputs. (See the Flatten-code below to see how to convert the result into your format.)

Let's run an experiment. On my 8 core the timings and the errors are as follows:

a = 2.;
b = 3.;
c = 4.;
N1 = 10000;
iters = 100;

x = randomSpherePoint[N1/2];
y = randomSpherePoint[N1/2];

spins = Riffle[x, y];
result1 = Table[
     list1 = spins[[1 ;; N1 ;; 2]]; list2 = spins[[2 ;; N1 ;; 2]]; 
     list3 = MapThread[Dot, {list1, list2}]; 
     list31 = MapThread[Cross, {list1, list2}]; 
     list32 = Sqrt[a^2 + b^2 + c^2 (1 - list3^2) + 2 a b list3]; 
     list4 = (a list1 + b list2 + c list31)/(list32); 
     list5 = (a list2 + b list1 - c list31)/(list32); 
     spins = Normalize /@ Flatten[Transpose[Join[{list4, list5}]], 1]
     , {i, 1, iters}]; // AbsoluteTiming // First

result2 = cf[x, y, a, b, c, iters]; // AbsoluteTiming // First

Max[Abs[result1 - Flatten[result2, {{2}, {1, 3}, {4}}]]]

8.67334

0.004071

3.98293*10^-15

That's a speed-up by a factor of about 2000.

Final remark:

One could certainly squeeze out a bit more performance by writing the C code oneselves and creating a LibraryLink function. But that is also a bit more fiddly.