Cross product between two lists of vectors

Question

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.

Answer 1

How about this

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

Answer 2

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.