# Speeding up cross and dot products for a list of vectors

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.

• Just to be clear on your notation: is $f$ supposed to return a list of two unit vectors in $\mathbb{R}^3$? Jun 2 at 18:34
• Also, do you need all 10^6 iterates of the process? Or do you just need the final result? Jun 2 at 18:53
• I should get back to my own work, but a few ideas: (1) If you're going to be renormalizing at each step, you don't need to calculate list32, since its purpose is to renormalize the result. (2) Look into ParallelTable, which may save you some time. Jun 2 at 19:06
• Yes it returns two vectors. Well Actually I do need all of them, basically it would run for 10^3 units of time for 10^5 different ensembles of spin. The list32 was a nice catch, thanks. About ParallelTable, not entirely sure whether it would help since Listable does multithreading, but will check it out anyway. Jun 2 at 19:30

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 = CompileGetElement[x, 1];
x2 = CompileGetElement[x, 2];
x3 = CompileGetElement[x, 3];

y1 = CompileGetElement[y, 1];
y2 = CompileGetElement[y, 2];
y3 = CompileGetElement[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.

• I actually thought I did not write too bad a code, but now I see it was crap. Thanks, I wasn't sure how to use Compile in whatever I was doing, and didn't expect MapThread to be this slow. I was contemplating doing the entire thing on FORTRAN if I did not get anything here. Thanks as always, I learnt a lot from this reply. Jun 2 at 20:57
• You are welcome. =) Jun 2 at 20:59
• Btw. I would not consider your code "crap". It's just that Table, Map and MapThread are not overly fast for such kind of microcomputations. And that's simply because Mathematica is an interpreted language. But fortunately, it provides means to compile some code parts -- which is certainly less painful than coding FORTRAN and link it somehow to the Mathematica session. Jun 3 at 6:39
• Okay doing that now. Jun 3 at 9:05
• mathematica.stackexchange.com/questions/269022/… Jun 3 at 9:29