# Elementwise, compilable min function

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 RandomPrivateMapThreadMin[{listX,ListY}] call
4. Compiled RandomPrivateMapThreadMin[{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 RandomPrivateMapThreadMin[{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 RandomPrivateMapThreadMin[{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 = RandomPrivateMapThreadMin[ {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 = RandomPrivateMapThreadMin[ {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


-
"Why does the example using RandomPrivateMapThreadMin[{listX, ListY}] fail to compile?" - it's not in the list here. – J. M. Sep 28 '12 at 13:14
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 RandomPrivateMapThreadMin, unfortunately... – Alex D. Sep 28 '12 at 13:57

I don't seem to have RandomPrivateMapThreadMin 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[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[r2 = cf[a, b];]
r1 === r2


{2.137, Null}

{0.109, Null}

True

-
+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
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 RandomPrivateMapThreadMin 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
@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