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 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. >>
*)
share|improve this question
    
You probably need to do the threading manually. –  Mr.Wizard Feb 14 '12 at 19:46
    
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 Feb 14 '12 at 19:53
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 Feb 14 '12 at 19:58
add comment

3 Answers

up vote 11 down vote accepted

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.

share|improve this answer
    
How did you figure that out?? That Boole is not welcome inside MapThread. –  István Zachar Feb 14 '12 at 21:18
    
@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 Feb 14 '12 at 21:24
4  
+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. Feb 14 '12 at 21:34
    
@OleksandrR. of course, silly me, thanks for pointing that out. –  Andy Ross Feb 14 '12 at 21:42
add comment

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
];
share|improve this answer
    
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 Feb 14 '12 at 20:10
    
@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 Feb 14 '12 at 20:16
    
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 Feb 14 '12 at 20:19
    
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 Feb 14 '12 at 20:20
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 Feb 14 '12 at 21:22
show 9 more comments

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.

share|improve this answer
add comment

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.