Function does not compile with Greater in it

13

6

I have to speed up a matrix-calculation and would like to use Compile for it, though it fails and produces an unknown error message:

compiledFunc = 
  Compile[{{w, _Real, 2}, {v, _Integer, 1}, {hb, _Real, 
     1}, {vb, _Real, 1}}, Module[{h, hs, vr, hr},
    h = 1./(1. + Exp[-(w.v + hb)]);
    hs = Boole@Thread[h > RandomReal[{0, 1}, {Length@hb}]];
    vr = 1./(1. + Exp[-((Transpose@w).hs + vb)]);
    hr = 1./(1. + Exp[-(w.vr + hb)]);
    w + (0.01*(Outer[Times, h, v] - Outer[Times, hr, vr]))
    ], Parallelization -> True];

(*
  ==> Compile::argcompten: The comparison, Greater, is invalid for tensor arguments. >>
*)

Any idea how I can overcome this and/or speed it up?

Update:

According to the answers, the problem is that Thread is not compilable. Substituting it with MapThread indeed allows compilation. But then running the code gives another error. (I thought that since it is about the same piece of code, I won't start a new post.)

cFunc = Compile[{{w, _Real, 2}, {v, _Integer, 1}, {hb, _Real, 
     1}, {vb, _Real, 1}}, Module[{h, hs, vr, hr},
    h = 1./(1. + Exp[-(w.v + hb)]);
    hs = Boole@
      MapThread[Greater, {h, RandomReal[{0, 1}, {Length@hb}]}];
    vr = 1./(1. + Exp[-((Transpose@w).hs + vb)]);
    hr = 1./(1. + Exp[-(w.vr + hb)]);
    w + (0.01*(Outer[Times, h, v] - Outer[Times, hr, vr]))],
   Parallelization -> True];

cFunc[w, v, hb, vb]

(*
   CompiledFunction::cfse: Compiled expression 1.7+3.2 I should be a machine-size real number. >>

   CompiledFunction::cfex: Could not complete external evaluation at instruction 2; proceeding with uncompiled evaluation. >>
*)

István Zachar

Posted 2012-02-14T19:36:02.090

Reputation: 44 135

You probably need to do the threading manually. – Mr.Wizard – 2012-02-14T19:46:49.180

MMA is assuming that h is a tensor. How do you determine that a tensor is larger than some real number? You need to provide a metric in order for Greater to work. – Matariki – 2012-02-14T19:53:51.663

1@Matariki The purpose of the Thread is to thread the > over elements, so what you say is not a problem. However, Thread does not get compiled, and that is actually the problem. – acl – 2012-02-14T19:58:11.333

Answers

11

The following seems to work based on what I assumed would be valid input.

compiledFunc = 
  Compile[{{w, _Real, 2}, {v, _Integer, 1}, {hb, _Real, 
     1}, {vb, _Real, 1}}, 
   Module[{h, hs, vr, hr, rr}, h = 1./(1. + Exp[-(w.v + hb)]);
    rr = RandomReal[{0, 1}, Length@hb];
    rr = h - rr;
    hs = UnitStep[rr] Unitize[rr];
    vr = 1./(1. + Exp[-((Transpose@w).hs + vb)]);
    hr = 1./(1. + Exp[-(w.vr + hb)]);
    w + (0.01*(Outer[Times, h, v] - Outer[Times, hr, vr]))], 
   Parallelization -> True];

Edit:

Since you are comparing to random reals we can safely assume none of the elements of hb will ever equal the pseudo-random values and remove the Unitize[rr] from the code. The above is for a more general comparison where equality can occur. Here is the cleaned up but less general version.

compiledFunc = 
  Compile[{{w, _Real, 2}, {v, _Integer, 1}, {hb, _Real, 
     1}, {vb, _Real, 1}}, 
   Module[{h, hs, vr, hr}, h = 1./(1. + Exp[-(w.v + hb)]);
    hs = UnitStep[h-RandomReal[{0, 1}, Length@hb]];
    vr = 1./(1. + Exp[-((Transpose@w).hs + vb)]);
    hr = 1./(1. + Exp[-(w.vr + hb)]);
    w + (0.01*(Outer[Times, h, v] - Outer[Times, hr, vr]))], 
   Parallelization -> True];

As for the second part of the question, apparently Compile doesn't like Boole being applied outside of MapThread. OleksandrR explained why in the comments. Compile cannot work with the boolean tensor that is created by MapThread. You can replace the line with

MapThread[Boole[#1 > #2] &, {h, RandomReal[{0, 1}, Length@hb]}]

and it should work.

Andy Ross

Posted 2012-02-14T19:36:02.090

Reputation: 18 640

when there is some error from using map in compile, with adding IntegerPart@ to the list that will be map, the error will not occur again. but why? – jack cilba – 2016-01-11T21:42:15.120

How did you figure that out?? That Boole is not welcome inside MapThread. – István Zachar – 2012-02-14T21:18:36.767

@IstvánZachar I used <<"CompiledFunctionTools" and then CompilePrint[compiledFunc] and noticed a nasty call to MainEvaluate around the thing. Bringing the Boole inside fixed the issue. As to why this is the case, I'm not sure. – Andy Ross – 2012-02-14T21:24:02.657

6+1. Just wanted to add (as I think it is not necessarily obvious) that the reason why the MapThread[Greater, _] formation doesn't compile is that boolean tensors are not supported in the Mathematica VM. By placing Boole inside MapThread, boolean scalars generated by Greater get coerced into a real tensor, which is supported. – Oleksandr R. – 2012-02-14T21:34:39.210

@OleksandrR. of course, silly me, thanks for pointing that out. – Andy Ross – 2012-02-14T21:42:25.280

10

The problem is not because of Greater, but Thread, which is not compilable.

MemberQ[Compile`CompilerFunctions[], Thread]
Out[1]= False

However, MapThread is! So you just need to replace the offending line involving Thread with:

hs = MapThread[Greater, {h, RandomReal[{0, 1}, {Length@hb}]}];

and it compiles fine. Here's the complete function:

compiledFunc = 
    Compile[{{w, _Real, 2}, {v, _Integer, 1}, {hb, _Real, 1}, {vb, _Real, 1}}, 
        Module[{h, hs, vr, hr},
            h = 1./(1. + Exp[-(w.v + hb)]);
            hs = MapThread[Greater, {h, RandomReal[{0, 1}, {Length@hb}]}];
            vr = 1./(1. + Exp[-((Transpose@w).hs + vb)]);
            hr = 1./(1. + Exp[-(w.vr + hb)]);
            w + (0.01*(Outer[Times, h, v] - Outer[Times, hr, vr]))
        ], 
    Parallelization -> True
];

rm -rf

Posted 2012-02-14T19:36:02.090

Reputation: 85 395

when I tried this, a) it didn't compile, b) I was asked to click on something that took me to a page that says, inter alia: "If you have reached this page from Mathematica, it means that Mathematica has discovered that one of these checks has failed." :) – acl – 2012-02-14T20:10:15.807

@acl Hmm... I have no clue. I tried compiling again in a fresh kernel and it had no problems. I've included the full function with the line replaced. Could you please try it in a new kernel? I'm using 8.0 for Mac OS X x86 (64-bit) (October 5, 2011) – rm -rf – 2012-02-14T20:16:18.200

I diligently checked that post you have linked, before posting my question, but since I thought that Greater is the weak link, I did not check on Thread. The error message is thus rather useless in this case, one wonders why it does not point to Thread instead of the function that is to be threaded. – István Zachar – 2012-02-14T20:19:49.453

Of course the crash is not reproducible, it only happened randomly when I ran this! However, to see what I meant about compilation, try SetSystemOptions[ "CompileOptions" -> "CompileReportExternal" \[Rule] True] before compiling, or CompilePrint[compiledFunc] after you have compiled it. – acl – 2012-02-14T20:20:02.237

@IstvánZachar It didn't say anything because the check failed at Greater. If it had been something that was valid, then it would've definitely moved on and warned/complained about Thread. Try using Boole@Thread[{1, 2, 3} == {4, 5, 6}] instead. I get the error Compile::extscalar: "Thread[{1,2,3}=={4,5,6}] cannot be compiled and will be evaluated externally." – rm -rf – 2012-02-14T20:27:12.417

@R.M You only get that because you have just set the option CompileReportExternal as I described in my comment above. You are not normally warned about this. – acl – 2012-02-14T20:29:55.750

@acl Oh, right. I forgot to mention that in my comment (too late to edit, so please note this, Istvan). Btw, I don't know what you (acl) mean is incorrect here... Hints? – rm -rf – 2012-02-14T20:32:49.710

@R.M maybe "incorrect" is not the right word, as your explanation is just fine (as is everybody else's). What I mean is that your compiledFunc emits the extscalar message, and if you use CompilePrint you'll see that it calls the main evaluator to do most of its work (ie, it's not actually compiled). – acl – 2012-02-14T20:35:46.977

@R.M It's not MapThread that is the problem. – acl – 2012-02-14T20:43:17.590

I have withdrawn my acceptance, since it now gives another error message, so I don't consider it "case closed". – István Zachar – 2012-02-14T20:57:39.127

1@IstvánZachar My mistake. I forgot to include Boole in my MapThread function. Please try Boole@(#1 > #2) & as the function instead of Greater. It works for me (no errors with compile report turned on). If it works for you too, then I'll edit it in. Otherwise, I have no clue what else is incorrect... – rm -rf – 2012-02-14T21:22:02.813

I see Andy's edit covers that part. – rm -rf – 2012-02-14T21:25:01.440

It seems to work with the Boole inside MapThread. Any idea how to give a shared accept? Both you and Andy gave answers that helped me to understand what is going on. – István Zachar – 2012-02-14T21:27:08.407

@IstvánZachar if you like the MapThread approach better than UnitStep I suggest you accept RM's answer. There isn't a huge speed difference as far as I can tell so its a matter of preference. – Andy Ross – 2012-02-14T21:44:19.283

9

Thread isn't compilable.

There are at least two ways of being warned that something is not being compiled: The first is to set SetSystemOptions[ "CompileOptions" -> "CompileReportExternal"->True] (this gives Compile::extscalar messages upon compilation), and examination of the compiled function after the fact with Needs["CompiledFunctionTools`"] followed by CompilePrint[compiledFunc]. This will show calls to the main evaluator.

To get this to compile, I find threading by hand the most readable approach:

compiledFunc = 
 Compile[{{w, _Real, 2}, {v, _Integer, 1}, {hb, _Real, 1}, {vb, _Real,
     1}}, Module[{h, hs, vr, hr}, h = 1./(1. + Exp[-(w.v + hb)]);
   hs = Table[0., {i, Length@hb}];
   Do[
    hs[[i]] = Boole[h[[i]] > RandomReal[{0, 1}]],
    {i, Length@hb}
    ];
   vr = 1./(1. + Exp[-((Transpose@w).hs + vb)]);
   hr = 1./(1. + Exp[-(w.vr + hb)]);
   w + (0.01*(Outer[Times, h, v] - Outer[Times, hr, vr]))], 
  Parallelization -> True]

As for the rest of the code, I am sure it can be further optimised.

acl

Posted 2012-02-14T19:36:02.090

Reputation: 19 146