# 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:= 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= {2.99032, Null}


I think it spends too much time. A modification

In:= B = Partition[Map[a, Flatten[Outer[#1 - #2 &, x, x]]], 1001]; // AbsoluteTiming
Out= {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. Nov 26 '15 at 18:51
• Exp[-Abs[x - #]] & /@ x is faster than 0.02 seconds on my machine. Nov 26 '15 at 18:54
• FWIW, the new DistanceMatrix[] function might be useful here. Nov 26 '15 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 Nov 26 '15 at 21:04
• @Dmitri By setting a to be Listable, you are essentially preventing vectorization. Arithmetic will once again be done element by element. Nov 26 '15 at 21:15
• Sorry, I have to say RTFM to myself :-) Nov 26 '15 at 21:18
• @Szabolcs This gives 2x further speed up: c = Exp[-Abs[x - # & /@ x]];/
– Kuba
Nov 27 '15 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. Nov 30 '15 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. Nov 27 '15 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 *)