# A quicker than outer

I need to get a matrix $\{a(x_i-x_j)\}$, where $x_i$ form a partition of an interval, $a(x)$ is a given function. I use

In[67]:= a[x_?NumericQ] := N[Exp[-Abs[x]]];
x = Table[-10 + 0.02 (j - 1), {j, 1, 1001}];
A = Outer[a[#1 - #2] &, x, x]; // AbsoluteTiming

Out[69]= {2.99032, Null}


I think it spends too much time. A modification

In[209]:= B = Partition[Map[a, Flatten[Outer[#1 - #2 &, x, x]]], 1001]; // AbsoluteTiming
Out[209]= {2.88966, Null}


does not help too much as well. Can I do this faster?

• The ideas from Leonid's answer here should be applicable to your problem as well. Try to vectorize, i.e. do arithmetic on vectors instead of scalars. Commented Nov 26, 2015 at 18:51
• Exp[-Abs[x - #]] & /@ x is faster than 0.02 seconds on my machine. Commented Nov 26, 2015 at 18:54
• FWIW, the new DistanceMatrix[] function might be useful here. Commented Nov 26, 2015 at 22:58

Vectorization will help a lot:

a[x_?NumericQ] := N[Exp[-Abs[x]]];
x = Table[-10 + 0.02 (j - 1), {j, 1, 1001}];
A = Outer[a[#1 - #2] &, x, x]; // AbsoluteTiming
(* {2.11988, Null} *)

B = Exp[-Abs[x - #]] & /@ x; // AbsoluteTiming
(* {0.016182, Null} *)

A == B
(* True *)


Notice that I am doing arithmetic on vectors the size of x instead of scalars. This is much faster than element-wise computation.

The idea is from Leonid's answer here:

• However, if you use a[x] directly the difference will not be so significant, I guess: SetAttributes[a, Listable]; B = a[x - #] & /@ x; // AbsoluteTiming Commented Nov 26, 2015 at 21:04
• @Dmitri By setting a to be Listable, you are essentially preventing vectorization. Arithmetic will once again be done element by element. Commented Nov 26, 2015 at 21:15
• Sorry, I have to say RTFM to myself :-) Commented Nov 26, 2015 at 21:18
• @Szabolcs This gives 2x further speed up: c = Exp[-Abs[x - # & /@ x]];/
– Kuba
Commented Nov 27, 2015 at 7:36
• I'm not really used to the term vectorization in this context, but should it really exclude using listable functions? My timings indicate that (res2 = Thread[Unevaluated@Subtract[a, b]]) // AbsoluteTiming // First is 4 times faster than a-b for random reals, which appears to "do arithmetic element by element". Interestingly, the version with Thread seems to parallelise automatically. Commented Nov 30, 2015 at 17:27

Outer is highly optimized for several built-in functions (Plus, Times, List). Therefore

Exp@-Abs@Outer[Plus, #, -#] &@Range[-10, 10, 0.02]; // RepeatedTiming
(* {0.025, Null} *)


gives ~50x speedup over Outer[#1 - #2&, #, #] and ~15x speedup over Outer[Subtract, #, #]. Also is a bit faster then Kuba's Exp[-Abs[x - # & /@ x]].

• Gosh, I remembered Outer being fast, but forgot it was only in special cases. Thanks. Commented Nov 27, 2015 at 21:53

Yes, two things help. The first is that Subtract is going to execute faster than #1 - #2 &, and the other is that all the operations involved in a are Listable, so getting rid of that _?NumericQ restriction speeds things up greatly. On my computer, this amounts to an order of magnitude speedup:

With[{x = Table[-10 + 0.02 (j - 1), {j, 1, 1001}]},
Outer[a[#1 - #2] &, x, x]]; // AbsoluteTiming
(* {2.29455, Null} *)

With[{x = Table[-10 + 0.02 (j - 1), {j, 1, 1001}]},
a[Outer[Subtract, x, x]]]; // AbsoluteTiming
(* {0.213449, Null} *)


It appears Outer can be reasonably compiled to C.

cfu = Compile[
{{x, _Real, 1}}
,
Outer[Exp[-Abs[# - #2]] &, x, x]
,
CompilationTarget -> "C"
]


Timings

A = cfu@x; // RepeatedTiming
B = Exp[-Abs[x - #]] & /@ x; // RepeatedTiming
A === B

{0.014,Null}
{0.020,Null}
True

Exp@-Abs@Outer[#1 - #2 &, #, #] &[Range[-10, 10, 0.02]]; //
AbsoluteTiming // First


0.950001

Borrowing from @Leonid Shifrin's approach, we can create a slightly better version than Outer gaining some improvement if not a lot

with built-in Outer and a as Listable:

SetAttributes[a, Listable];
a[x_?NumericQ] := N[Exp[-Abs[x]]];
x = Table[-10 + 0.02 (j - 1), {j, 1, 1001}];
(b = a[Outer[#1 - #2 &, x, x]];) // RepeatedTiming
(* 2.62 seconds *)


defining our own version of Outer

ourOuter[f_, arg1_, arg2_] := Module[{auxf},
Attributes[auxf] = {Listable};
auxf[x_][y_] := f[x, y];
Through[auxf[arg1][arg2]]
];

ClearAll[a];
SetAttributes[a, Listable];
a[x_] := N[Exp[-Abs[x]]];

(c = With[{x = Table[-10 + 0.02 (j - 1), {j, 1, 1001}]},
a[ownOuter[Subtract, x, x]]
];)//RepeatedTiming
(* 1.7 seconds *)

b == c
(* True *)