8
$\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.

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

2 Answers 2

Reset to default
11
$\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
Improve this answer
$\endgroup$
3
  • 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$ Commented Jun 28, 2016 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$ Commented Jun 28, 2016 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$ Commented Jun 28, 2016 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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ It would be good to compare speed using packed numerical arrays as well, e.g. RandomReal[{-1, 1}, {32, 32, 32, 3}]; $\endgroup$ Commented Jun 28, 2016 at 16:30
  • $\begingroup$ @SimonWoods. Done. $\endgroup$ Commented Jun 28, 2016 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$ Commented Jun 28, 2016 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$ Commented Jun 28, 2016 at 17:07

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.