A quicker than outer

Question

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?

Answer 1

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:

• https://mathematica.stackexchange.com/a/21863/12

Answer 2

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]].

Answer 3

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} *)

Answer 4

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

Answer 5

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

0.950001

Answer 6

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 *)