17
$\begingroup$

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?

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

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:

$\endgroup$
6
  • $\begingroup$ However, if you use a[x] directly the difference will not be so significant, I guess: SetAttributes[a, Listable]; B = a[x - #] & /@ x; // AbsoluteTiming $\endgroup$
    – Dmitri
    Nov 26 '15 at 21:04
  • 2
    $\begingroup$ @Dmitri By setting a to be Listable, you are essentially preventing vectorization. Arithmetic will once again be done element by element. $\endgroup$
    – Szabolcs
    Nov 26 '15 at 21:15
  • $\begingroup$ Sorry, I have to say RTFM to myself :-) $\endgroup$
    – Dmitri
    Nov 26 '15 at 21:18
  • 3
    $\begingroup$ @Szabolcs This gives 2x further speed up: c = Exp[-Abs[x - # & /@ x]];/ $\endgroup$
    – Kuba
    Nov 27 '15 at 7:36
  • 1
    $\begingroup$ 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. $\endgroup$ Nov 30 '15 at 17:27
15
$\begingroup$

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

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

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} *)
$\endgroup$
4
$\begingroup$

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
$\endgroup$
2
$\begingroup$
Exp@-Abs@Outer[#1 - #2 &, #, #] &[Range[-10, 10, 0.02]]; // 
  AbsoluteTiming // First

0.950001

$\endgroup$
1
$\begingroup$

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 *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.