7
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ Why not just MapThread[Cross, {{{a1, a2, a3}, {b1, b2, b3}}, {{c1, c2, c3}, {d1, d2, d3}}}]? $\endgroup$ – J. M. is away Jun 28 '16 at 16:10
  • 2
    $\begingroup$ You can just do this: MapThread[Cross, {list1, list2}, 3]. But, hm, it's actually slower! $\endgroup$ – march Jun 28 '16 at 16:12
  • 1
    $\begingroup$ 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!) $\endgroup$ – march Jun 28 '16 at 17:09
10
$\begingroup$

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
$\endgroup$
  • 1
    $\begingroup$ 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 :) $\endgroup$ – João Moutinho Jun 28 '16 at 17:37
  • 2
    $\begingroup$ @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. $\endgroup$ – Simon Woods Jun 28 '16 at 17:39
  • 1
    $\begingroup$ 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}). $\endgroup$ – J. M. is away Jun 28 '16 at 23:18
7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ It would be good to compare speed using packed numerical arrays as well, e.g. RandomReal[{-1, 1}, {32, 32, 32, 3}]; $\endgroup$ – Simon Woods Jun 28 '16 at 16:30
  • $\begingroup$ @SimonWoods. Done. $\endgroup$ – march Jun 28 '16 at 16:58
  • $\begingroup$ Interesting results... on my system my version is significantly faster. Do you get similar numbers with the more accurate RepeatedTiming ? $\endgroup$ – Simon Woods Jun 28 '16 at 17:04
  • $\begingroup$ @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. $\endgroup$ – march Jun 28 '16 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.