# Cross product between two lists of vectors

I'm trying to create a function that does the cross product between two lists of the same dimensions, element by element. Each entry of the list is a 3D vector. Something like this:

a = {{a1, a2, a3}, {b1, b2, b3}};
b = {{c1, c2, c3}, {d1, d2, d3}};

listCross[a, b] = {Cross[{a1, a2, a3}, {c1, c2, c3}], Cross[{b1, b2, b3}, {d1, d2, d3}]};


All my lists have 3 dimensions, so something like Dimensions[a] = {32, 32, 32, 3} would be typical, and I have the following code that works:

listCross[list1_, list2_] := Module[{vec11, vec12, vec13, result},
vec11 = -list1[[All, All, All, 3]]list2[[All, All, All, 2]] + list1[[All, All, All, 2]]list2[[All, All, All, 3]];
vec12 = +list1[[All, All, All, 3]]list2[[All, All, All, 1]] - list1[[All, All, All, 1]]list2[[All, All, All, 3]];
vec13 = -list1[[All, All, All, 2]]list2[[All, All, All, 1]] + list1[[All, All, All, 1]]list2[[All, All, All, 2]];
result = Partition[Partition[MapThread[{#1, #2, #3}&, {Flatten[vec11], Flatten[vec12], Flatten[vec13]}], Dimensions[list1][[3]]], Dimensions[list1][[2]]]
];


It exploits the result of Cross[{a1, a2, a3}, {b1, b2, b3}] and generalizes to lists with list sums. In the end it uses Partition and MapThread to join the x, y and z components into a new list with the same dimensions as the input lists.

However this is not as fast as I would like it to be, and I haven't been able to come up with anything better. I also tried doing the cross product in it's matrix-vector product form, creating a block diagonal sparse matrix from the first list that I would multiply with the Flatten of the second list, but I couldn't create the matrix fast enough.

Does anyone have any thoughts on this? I also know nothing about Compile, so no idea if it would be of use here. Any help will be greatly appreciated.

• Why not just MapThread[Cross, {{{a1, a2, a3}, {b1, b2, b3}}, {{c1, c2, c3}, {d1, d2, d3}}}]? Commented Jun 28, 2016 at 16:10
• You can just do this: MapThread[Cross, {list1, list2}, 3]. But, hm, it's actually slower! Commented Jun 28, 2016 at 16:12
• If you really feel the need to accept an answer right now, I suggest accepting @SimonWoods' answer, since his is by far faster than mine and not that much more verbose. But it also might be worth waiting a little before accepting an answer, because there are certainly other people on this site who are very good at performance-tuning, and not having an accepted answer is more likely to attract more good answers. (So feel free to de-accept my answer!) Commented Jun 28, 2016 at 17:09

furious[a_, b_] := Module[{a1, a2, a3, b1, b2, b3, c},
{a1, a2, a3} = Transpose[a, {2, 3, 4, 1}];
{b1, b2, b3} = Transpose[b, {2, 3, 4, 1}];
c = {-a3 b2 + a2 b3, a3 b1 - a1 b3, -a2 b1 + a1 b2};
Transpose[c, {4, 1, 2, 3}]]


Timing results (from march's answer) for version 10.4.1

list1 = RandomReal[{-1, 1}, {32, 32, 32, 3}];
list2 = RandomReal[{-1, 1}, {32, 32, 32, 3}];

l1 = listCrossMarch1[list1, list2]; // RepeatedTiming // First
l2 = listCrossMarch2[list1, list2]; // RepeatedTiming // First
l3 = listCross[list1, list2]; // RepeatedTiming // First
l4 = furious[list1, list2]; // RepeatedTiming // First

2.67

0.6064

0.0386

0.0015

• First of all, thanks a lot for the answers from everyone! Just want to add that running these same tests on my system I get the same timings you're listing here. As suggested by march I'm going to leave this up a couple more hours to see if anyone else has anything to add :) Commented Jun 28, 2016 at 17:37
• @JoãoMoutinho We often suggest leaving a question at least 24 hours before accepting an answer - this gives users in all time zones an opportunity to contribute. Commented Jun 28, 2016 at 17:39
• I suspect you can squeeze out further gains if you can somehow exploit LeviCivitaTensor[] (recall that Cross[{a, b, c}, {u, v, w}] == {u, v, w}.LeviCivitaTensor[3].{a, b, c}). Commented Jun 28, 2016 at 23:18

Interestingly enough, MapThreading Cross works but is much slower:

Using sample lists:

list1 = Array[c, {20, 20, 20, 3}];
list2 = Array[d, {20, 20, 20, 3}];


We can perform this operation in the following two ways, using MapThread:

listCrossMarch1[list1_, list2_] := MapThread[Cross, {list1, list2}, 3]
listCrossMarch2[list1_, list2_] := MapThread[{#1[[2]] #2[[3]] - #1[[3]] #2[[2]], #1[[3]] #2[[1]] - #1[[1]] #2[[3]], #1[[1]] #2[[2]] - #1[[2]] #2[[1]]} &, {list1, list2}, 3]


As we can see below, there is a lot of overhead associated with Cross apparently, since that version is much slower than coding the cross-product explicitly. The MapThread version with the explicit cross-product (rather than Cross) is almost as fast as the OP's version and much cleaner to write down. Simon Woods' answer is the fastest (and also very clean).

l1 = listCrossMarch1[list1, list2]; // AbsoluteTiming // First
l2 = listCrossMarch2[list1, list2]; // AbsoluteTiming // First
l3 = listCross[list1, list2]; // AbsoluteTiming // First
l4 = furious[list1, list2]; // AbsoluteTiming // First
l1 === l2 === l3 === l4
(* 1.436369 *)
(* 0.190366 *)
(* 0.120962 *)
(* 0.084740 *)
(* True *)


As suggested by Simon Woods, here are Timings with packed arrays of real numbers. Using

list1 = RandomReal[{-1, 1}, {20, 20, 20, 3}];
list2 = RandomReal[{-1, 1}, {20, 20, 20, 3}];


we do

l1 = listCrossMarch1[list1, list2]; // AbsoluteTiming // First
l2 = listCrossMarch2[list1, list2]; // AbsoluteTiming // First
l3 = listCross[list1, list2]; // AbsoluteTiming // First
l4 = furious[list1, list2]; // AbsoluteTiming // First
(* 0.429275 *)
(* 0.094540 *)
(* 0.008337 *)
(* 0.028592 *)


The OP's messy version is significantly the fastest here! Simon Woods' answer still does a great job, and Cross still has lots of overhead.

• It would be good to compare speed using packed numerical arrays as well, e.g. RandomReal[{-1, 1}, {32, 32, 32, 3}]; Commented Jun 28, 2016 at 16:30
• @SimonWoods. Done. Commented Jun 28, 2016 at 16:58
• Interesting results... on my system my version is significantly faster. Do you get similar numbers with the more accurate RepeatedTiming ? Commented Jun 28, 2016 at 17:04
• @SimonWoods. Unfortunately, I have V10.0, so I don't have RepeatedTiming. So my timing measurements are always suspect. Perhaps you could do the timings instead? It wouldn't hurt. Commented Jun 28, 2016 at 17:07