# How to speed up Outer[Subtract, vect1,vect2]?

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.

• 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. Feb 15 at 13:20
• Also have a look into the Mathematica function NBodySimulation. Feb 15 at 13:21
• 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). Feb 15 at 14:05
• Btw., you should rather use Outer[Subtract, vect1,vect2,1] (mind the fourth argument). Feb 15 at 14:10

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

• Relevant answer from @ybeltukov: "Outer is highly optimized for several built-in functions (Plus, Times, List)". 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


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.

• Thank you for this add to Karl's answer Feb 18 at 16:37
• You're welcome. 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


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