Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to implement efficiently a transfer-matrix like algorithm. On each iteration, I have two vectors $x=\{x_1,\dots,x_n\}$, $y=\{y_1,\dots,y_n\}$ with real numbers and I need to compute the vector $\{\min(x_1,y_1),\dots,\min(x_n,y_n)\}$. I tried four approaches for computing it:

  1. Uncompiled MapThread[Min,{listX,listY}] call
  2. Compiled MapThread[Min,{listX,listY}] call
  3. Uncompiled Random`Private`MapThreadMin[{listX,ListY}] call
  4. Compiled Random`Private`MapThreadMin[{listX,ListY}] call

(Code see below). The resulting timings were: 4.5s (for 1), 3.5s (for 2), 1.5s (for 3) and 4 reverted to uncompiled evaluation, giving 6.3s.

So my questions are:

  1. Is the uncompiled Random`Private`MapThreadMin[{listX, ListY}] call the fastest way to evaluate the element-wise minimum of two lists, or does anybody have a better idea?
  2. Why does the example using Random`Private`MapThreadMin[{listX, ListY}] fail to compile?

My code examples are:

it1[wd_, len_] := 
    Module[{pot1, fval}, 
    pot1 = RandomVariate[NormalDistribution[], {len, wd}];
    fval = ConstantArray[0., wd];
    Do[fval = MapThread[Min, {RotateLeft[fval], fval}] + pot1[[k]];, {k, 1, len}];
    Return[fval]];

it2 := Compile[{{wd, _Integer}, {len, _Integer}},
    Module[{pot1, fval}, 
    pot1 = RandomVariate[NormalDistribution[], {len, wd}];
    fval = ConstantArray[0., wd];
    Do[fval = MapThread[Min, {RotateLeft[fval], fval}] + pot1[[k]];, {k, 1, len}];
    Return[fval]]];

it3[wd_, len_] := 
    Module[{pot1, fval}, 
    pot1 = RandomVariate[NormalDistribution[], {len, wd}];
    fval = ConstantArray[0., wd];
    Do[fval = Random`Private`MapThreadMin[ {RotateLeft[fval], fval}] + pot1[[k]];, {k, 1, len}];
    Return[fval]];

it4 := Compile[{{wd, _Integer}, {len, _Integer}},
    Module[{pot1, fval}, 
    pot1 = RandomVariate[NormalDistribution[], {len, wd}];
    fval = ConstantArray[0., wd];
    Do[fval = Random`Private`MapThreadMin[ {RotateLeft[fval], fval}] + pot1[[k]];, {k, 1, len}];
    Return[fval]]];

And to obtain the timing values, I used

Table[it1[20, 10] // First, {10000}]; // AbsoluteTiming
Table[it2[20, 10] // First, {10000}]; // AbsoluteTiming
Table[it3[20, 10] // First, {10000}]; // AbsoluteTiming
Table[it4[20, 10] // First, {10000}]; // AbsoluteTiming

Thank you in advance!

share|improve this question
    
"Why does the example using Random`Private`MapThreadMin[{listX, ListY}] fail to compile?" - it's not in the list here. –  J. M. Sep 28 '12 at 13:14
1  
For your consideration: minNew = Compile[{{x, _Real}, {y, _Real}}, Min[x, y], RuntimeAttributes -> {Listable}]. Now, try giving minNew[] two lists as arguments... –  J. M. Sep 28 '12 at 13:22
    
Thanks for the suggestion! I tried your minNew function, it works faster than the pure MapThread call but about 20% slower than Random`Private`MapThreadMin, unfortunately... –  Alex D. Sep 28 '12 at 13:57

1 Answer 1

up vote 6 down vote accepted

I don't seem to have Random`Private`MapThreadMin in version 7.

For Integer data you may wish to try:

a = RandomInteger[1*^6 {-1, 1}, 5*^6];
b = RandomInteger[1*^6 {-1, 1}, 5*^6];

Timing[r1 = MapThread[Min, {a, b}];]
Timing[r2 = a # + b (1 - #) &@UnitStep[b - a];]
r1 === r2

{2.028, Null}

{0.171, Null}

True

For machine-size Real data this compiles nicely:

cf = Compile[{{a, _Real, 1}, {b, _Real, 1}}, a # + b (1 - #) & @ UnitStep[b - a]];

a = N@a;
b = N@b;

Timing[r1 = MapThread[Min, {a, b}];]
Timing[r2 = cf[a, b];]
r1 === r2

{2.137, Null}

{0.109, Null}

True

share|improve this answer
    
+1, I was going to suggest something similar. –  Leonid Shifrin Sep 28 '12 at 13:22
    
@Leonid good, that means this might actually work. :^) –  Mr.Wizard Sep 28 '12 at 13:30
1  
Indeed, this works perfectly, thanks very much! I wrote a compiled version of your suggestion and it runs about 20% faster than my fastest version above (option 3, the Random`Private`MapThreadMin code. I learned of the existence of this function in Mathematica 8 from link). Thanks again! –  Alex D. Sep 28 '12 at 13:56
    
@Alex thanks for the Accept, but I encourage all users to wait 24 hours before Accepting an answer. Who knows what methods may be posted if people are not discouraged from reading your question! –  Mr.Wizard Sep 28 '12 at 14:15
1  
@xzczd You could Fold the operation quite quickly. If your have many lists (more lists than the length of each list, for example) then you would do better to Transpose, e.g. Min /@ Transpose[lists]. –  Mr.Wizard Oct 22 '13 at 14:57

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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