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?

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. — Szabolcs, Nov 26 '15 at 18:51
Exp[-Abs[x - #]] & /@ x is faster than 0.02 seconds on my machine. — Szabolcs, Nov 26 '15 at 18:54
FWIW, the new DistanceMatrix[] function might be useful here. — J. M.'s missing motivation, Nov 26 '15 at 22:58

score 18 · Accepted Answer · edited Apr 13 '17 at 12:55

18

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

edited Apr 13 '17 at 12:55

Community

1

answered Nov 26 '15 at 19:35

Szabolcs

234,956
30
623
1,263

However, if you use a[x] directly the difference will not be so significant, I guess: SetAttributes[a, Listable]; B = a[x - #] & /@ x; // AbsoluteTiming – Dmitri Nov 26 '15 at 21:04
2

@Dmitri By setting a to be Listable, you are essentially preventing vectorization. Arithmetic will once again be done element by element. – Szabolcs Nov 26 '15 at 21:15
Sorry, I have to say RTFM to myself :-) – Dmitri Nov 26 '15 at 21:18
3

@Szabolcs This gives 2x further speed up: c = Exp[-Abs[x - # & /@ x]];/ – Kuba Nov 27 '15 at 07:36
1

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. – Jacob Akkerboom Nov 30 '15 at 17:27
Sorry for my poor choice of variable names – Jacob Akkerboom Nov 30 '15 at 17:38

ybeltukov · Answer 2 · 2015-11-27T22:04:09.853

18

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

edited Nov 27 '15 at 22:04

answered Nov 27 '15 at 21:47

ybeltukov

43,673
5
108
212

1

Gosh, I remembered Outer being fast, but forgot it was only in special cases. Thanks. – Michael E2 Nov 27 '15 at 21:53

score 5 · Answer 3 · answered Nov 26 '15 at 19:21

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

score 4 · Answer 4 · answered Dec 08 '15 at 09:22

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

eldo · Answer 5 · 2015-11-26T19:34:35.220

2

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

0.950001

edited Nov 26 '15 at 19:34

answered Nov 26 '15 at 19:13

eldo

67,911
5
60
168

score 1 · Answer 6 · answered Jun 13 '17 at 19:04

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

A quicker than outer

6 Answers6

Linked