16

Suppose I have two arrays of data, which I want to divide:

n = 1000;
a = RandomInteger[{1, 9}, {n, n}];
b = RandomInteger[{1, 9}, {n, n}];
First @ AbsoluteTiming @ (a/b)

0.637064

But suppose that two thirds of the values are actually zero. I'd want the ones dividing by zero to just remain zero.

mask = RandomChoice[{0, 0, 1}, {n, n}];
c = mask a;
d = mask b;
First @ AbsoluteTiming @ Quiet[c/d /. Indeterminate -> 0]

8.568857

Now it has become very slow. How can we speed this up?

Edit: Thanks for all the answers guys. Below are some timings:

timings

wxffles
  • 14,246
  • 1
  • 43
  • 75

5 Answers5

8

Simply turn all the 0/0s into 0/1s:

ans2 = c/(1 - Sign@d + d); // AbsoluteTiming

{0.7630000, Null}

The above answer only works for specific example in the question i.e. it only handles positive lists and identical mask. For more general cases one can use:

n = 1000;
(* a and b contain non-positive elements now *)
origin := RandomInteger[{-9, 9}, {n, n}]
a = origin; b = origin;

(* c and d now use different masks *)    
mask := RandomChoice[{0, 0, 1}, {n, n}];
c = mask a; d = mask b;

ans2 = With[{s = Unitize@d}, Divide[s c, 1 - s + d]]; // AbsoluteTiming

{0.606000, Null}

Notice that Divide is faster than /.

Community
  • 1
xzczd
  • 65,995
  • 9
  • 163
  • 468
  • 1
    The best general answer so far IMHO. – Mr.Wizard Oct 28 '14 at 08:55
  • 1
    Note that you can use Unitize@d in place of Abs@Sign@d. – Mr.Wizard Oct 28 '14 at 08:59
  • 1
    @Mr.Wizard It makes sense that one might know beforehand that 0/0 ought to be 0 - in any case, Indeterminate -> 0 is specified in the question. The question leaves open what should happen to a result like 1/0 or ComplexInfinity. One might want infinity, zero or a different result, such as a maximum value. The general answer here effectively takes the union of the masks, which is a plausible (even common) use-case. But I think one ought to be careful in stating which generalized problem it is a solution to. - +1 long before now. :) – Michael E2 Oct 28 '14 at 12:40
  • Clever and fast, but I'm going to wonder what the heck it is doing when I come back to read the code a few months later. Comments are for the weak. – wxffles Oct 28 '14 at 21:39
  • @wxffles I think the most probable obstacle for one to understand the code is unfamiliarity with Unitize and vectorization :) …OK, I've added a bit explanation, just a bit though. – xzczd Nov 01 '14 at 03:29
8

SparseArray can help, given the size and nature of the mask. It's slightly faster to convert c and d to sparse arrays than to convert a and b.

mask = SparseArray@RandomChoice[{0, 0, 1}, {n, n}];

First@AbsoluteTiming[
 c = mask a;
 d = mask b;
 Quiet[foo2 = Block[{Indeterminate = 0}, SparseArray[c] / SparseArray[d]]]
 ]
(*
  0.151565
*)

Compare:

First@AbsoluteTiming[
  c = mask a;
  d = mask b;
  foo1 = Quiet[c/d /. Indeterminate -> 0]
  ]
(*
  6.492534
*)

foo2 == foo1
(*
  True
*)
Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • @kguler I'm not sure why you deleted your Block trick, which I've used because it's faster than any other trick I've tried. But let me know if there's a problem. My use of SparseArray is a significant difference, but Replace does not work on it. Block and Map were the only tools I had, and `Map was way too slow. – Michael E2 Oct 28 '14 at 04:37
  • 1
    The trick will fail when c and d are using different masks. (Block won't help because the FullForm of the output of ComplexInfinity generated by 1/0 etc. is DirectedInfinity[].) Anyway, it's the fastest method for OP's specific example so far. – xzczd Oct 28 '14 at 05:38
  • 2
    @xzczd Actually it seems very odd to have different masks. In the OP's Q, it's clear that Indeterminate is to go to 0. 1/0 does not become Indeterminate and there is no indication that ComplexInfinity should become 0. (I decided it probably should not.) – Michael E2 Oct 28 '14 at 12:15
  • I like this answer the best. It's a bit wordy, but clear what it is doing. It's fast and flexible (if you want to do more than a simple divide). It's not quite so fast if you want to do a Normal afterwards - about the same speed as the case with no zeroes. – wxffles Oct 28 '14 at 21:34
  • Michael, just saw your comment (somehow i think @name mechanism is broken; i don't get pinged when mentioned in comments). The only reason i deleted my answer was I forgot to refresh the kernel in a long/crowded mma session and got confused with conflicting results:) (+1) – kglr Oct 29 '14 at 03:05
  • @kguler The at-name only works if you've posted a comment already. I put it there just in case you saw this answer. I knew you wouldn't get pinged by it, but I didn't know how to ping you. – Michael E2 Oct 29 '14 at 11:49
7
SetAttributes[f, Listable];
f@__ := 0
f[a_, b_] /; b != 0 := a/b
First@AbsoluteTiming[f[c, d];]

Edit

On the same spirit, this is 50% faster:

SetAttributes[g, Listable];
g[_, 0] := 0
g[a_, b_] := Divide[a, b]
First@AbsoluteTiming[g[c, d];]
Dr. belisarius
  • 115,881
  • 13
  • 203
  • 453
2
c/(d /. (0 | 0. -> Infinity));
Basheer Algohi
  • 19,917
  • 1
  • 31
  • 78
2

For the given example where no values in b are zero before the mask one can use:

Divide[mask*a, b]

Note the use of explicit Divide for optimum performance. Timings:

(* your data *)

(r1 = Quiet[c/d /. Indeterminate -> 0]); // AbsoluteTiming // First
(r2 = Divide[mask*a, b]);                // AbsoluteTiming // First

r1 == r2

4.578262

0.212012


True

For the general problem I cannot beat xzczd's method.

Community
  • 1
Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371