5
$\begingroup$

My Nbody code use a lot of Outer[Subtract, vect1,vect2]. A lot means ten thousands to one hundred thousands times !! This is at the heart of my Nbody code, and I need to speed up this instruction. Each vector contains for example N=10000 values, and it scales like N^2.

By compiling this instruction , I obtain a factor of 4 faster (4 CPU !!). Not bad , but I dont want to buy 128 CPU !.

I dont understand for example why Outer[Times,vect1,vect2] is about 100 times faster than Outer[Subtract, vect1,vect2].It would be great to get this speed-up !

Is there a way to speed up Outer[Subtract, vect1,vect2]? Using C++, or using GPU , or else ??

Thank you for getting some advices.

Improve this question
$\endgroup$
6
  • $\begingroup$ Please see this question and the corresponding answers and comments. $\endgroup$
    – Domen
    Commented Feb 15 at 10:26
  • 1
    $\begingroup$ This is an XY-problem. You are using an all-pairs algorithm which is inherintly $O(N^2)$. You should read about the Barnes-Hut method for $N$-body problems; it approximates the correct result in $O(N \log(N))$ time and space. $\endgroup$
    – Henrik Schumacher
    Commented Feb 15 at 13:20
  • $\begingroup$ Also have a look into the Mathematica function NBodySimulation. $\endgroup$
    – Henrik Schumacher
    Commented Feb 15 at 13:21
  • 1
    $\begingroup$ Even if you don't use Barnes-Hut: Using Outer leads to $O(N^2)$ memory consumption. And since the memory hierachy of this size rather slow, you should rather go for a (compiled!) double loop that needs only $O(N)$ memory to store the results (positions, velocities, forces). $\endgroup$
    – Henrik Schumacher
    Commented Feb 15 at 14:05
  • $\begingroup$ Btw., you should rather use Outer[Subtract, vect1,vect2,1] (mind the fourth argument). $\endgroup$
    – Henrik Schumacher
    Commented Feb 15 at 14:10

3 Answers 3

Reset to default
16
$\begingroup$ 
x = RandomReal[1., 10000];
y = RandomReal[1., 10000];
a = Outer[Subtract, x, y]; // RepeatedTiming
b = Outer[Plus, x, -y]; // RepeatedTiming 
a == b

{13.667, Null}

{0.270628, Null}

True

Improve this answer
$\endgroup$
1
  • 7
    $\begingroup$ Relevant answer from @ybeltukov: "Outer is highly optimized for several built-in functions (Plus, Times, List)". $\endgroup$
    – Domen
    Commented Feb 15 at 10:27
5
$\begingroup$

In case you want to do that for vector-valued data (because you mentiond N-body problems): Outer seems to be best optimized for scalar outputs. So working to "structure of arrays" format might be useful.

n = 10^4;
P = RandomReal[{-1, 1}, {n, 3}];
Q = RandomReal[{-1, 1}, {n, 3}];

RTimes = Outer[Times, P, Q, 1]; // AbsoluteTiming // First
RSubtract = Outer[Subtract, P, Q, 1]; // AbsoluteTiming // First
RMinus = Outer[Plus, P, -Q, 1]; // AbsoluteTiming // First

2.14531

15.3862

1.84677

Here the version that I mean. (I am using Karl's very good idea here.)

(*Convert to structure of arrays format.*)
PT = Transpose[P];
QT = Transpose[Q];

(*Get result in structure of arrays format.*)
RT = {
  Outer[Plus, PT[[1]], -QT[[1]]], 
  Outer[Plus, PT[[2]], -QT[[2]]], 
  Outer[Plus, PT[[3]], -QT[[3]]]
}; // AbsoluteTiming // First

Dimensions[RT]

(*Convert to array of structures format (just for comparison; don't do that in your code!)*)
R = Transpose[RT, {3, 1, 2}];

Max[Abs[RMinus - R]]

0.350867

{3, 10000, 10000}

0.

Nonetheless, you should rather use the Barnes-Hut method or the fast multipole method for approximating the interactions instead. There are a bazillion of code bases out there. Most of them will be written in Fortran, C, or C++, though. One cannot code them efficiently in an interpreted language like Mathematica. So I suggest you look for some C++ library and link it with LibraryLink.

Edit

I also compared to a parallelized implementation in C++ and even a blocked version (for avoiding cache misses). Both take about 0.143006 seconds. Single-threaded they do it in 0.325479, which shows that this problem is likely memory bound. (I know from other experiments that two of my 8 CPU cores can fully satiate the memory bandwitdh of my machine.) So I don't think that you can get out much better performance here.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you for this add to Karl's answer $\endgroup$
    – Rémy Galli
    Commented Feb 18 at 16:37
  • $\begingroup$ You're welcome. $\endgroup$
    – Henrik Schumacher
    Commented Feb 18 at 17:53
2
$\begingroup$

I can not explain what Outer with Subtract does. But if you do not use Outer you can get a similar performance as with Time by writing:

n = 10^4;
v1 = RandomReal[{-1, 1}, n];
v2 = RandomReal[{-1, 1}, n];

Outer[Times, v1, v2]; // AbsoluteTiming
Outer[Subtract, v1, v2]; // AbsoluteTiming
(v2 - #) & /@ v1; // AbsoluteTiming

enter image description here

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thank you for your answer but your proposal doesnt give the same results .... $\endgroup$
    – Rémy Galli
    Commented Feb 15 at 9:47
  • $\begingroup$ @RémyGalli Transpose[Subtract[v1, #] & /@ v2] should give the same as Outer[Subtract, v1, v2] $\endgroup$
    – Sjoerd Smit
    Commented Feb 15 at 9:49
  • $\begingroup$ Simply exchange v1 and v2. I changed this in my answer. $\endgroup$
    – Daniel Huber
    Commented Feb 15 at 10:00
  • $\begingroup$ @Sjoerd Smit Thanks a lot. I gain a factor of 5. Any other idea ? By compiling for example ? $\endgroup$
    – Rémy Galli
    Commented Feb 15 at 10:15
  • $\begingroup$ @DanielHuber It should be Subtract[#, v2] & /@ v1. $\endgroup$
    – Sjoerd Smit
    Commented Feb 15 at 10:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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