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.

  • $\begingroup$ Please see this question and the corresponding answers and comments. $\endgroup$
    – Domen
    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$ Feb 15 at 13:20
  • $\begingroup$ Also have a look into the Mathematica function NBodySimulation. $\endgroup$ 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$ Feb 15 at 14:05
  • $\begingroup$ Btw., you should rather use Outer[Subtract, vect1,vect2,1] (mind the fourth argument). $\endgroup$ Feb 15 at 14:10

3 Answers 3

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}


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

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




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


(*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]]


{3, 10000, 10000}


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.


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.

  • $\begingroup$ Thank you for this add to Karl's answer $\endgroup$ Feb 18 at 16:37
  • $\begingroup$ You're welcome. $\endgroup$ Feb 18 at 17:53

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

  • $\begingroup$ Thank you for your answer but your proposal doesnt give the same results .... $\endgroup$ Feb 15 at 9:47
  • $\begingroup$ @RémyGalli Transpose[Subtract[v1, #] & /@ v2] should give the same as Outer[Subtract, v1, v2] $\endgroup$ Feb 15 at 9:49
  • $\begingroup$ Simply exchange v1 and v2. I changed this in my answer. $\endgroup$ 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$ Feb 15 at 10:15
  • $\begingroup$ @DanielHuber It should be Subtract[#, v2] & /@ v1. $\endgroup$ 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.