4
$\begingroup$

In the mathematica document "tutorial/CompilingWolframLanguageExpressions"

1

It says about how to set return type of a called function you used in Compile.

Clear[com]
com[i_] := Binomial[2 i, i]

test=Compile[{x, {i, _Integer}}, x^com[i], {{com[_], _Integer}}]

From the above example, we can see that because the com function evaluates to integer, so we set {{com[_], _Integer}} to let Compile know the return type of com.

But if you inspect it further

Needs["CompiledFunctionTools`"]
CompilePrint[test]

You can see there is MainEvaluate when calling function com.

So, I don't understand the meaning of this example in the document. If a compiled function has a MainEvaluate process, then I think compiling it is just nonsense, for it won't speed up things, right?

2

Then I came up with another question, it is also mentioned in the same document page. If we compile Sqrt as below

sqrtcom1 = Compile[{{x, _Real}}, Sqrt[x]]

we will run into problems if we evaluate sqrtom1[-1.], it will give errors like

CompiledFunction::cfn: Numerical error encountered at instruction 1; proceeding with uncompiled evaluation. >>

This is because the return type of Sqrt is assumed to be real by default. This can be see from

In[20]:= ToCompiledProcedure[sqrtcom1][[4]]

Out[20]= CompiledResult[Register[Real, 1]]

So, theoretically we could solve this by

sqrtcom2 = Compile[{{x, _Real}}, Sqrt[x], {{Sqrt[_], _Complex}}]

But this is not working!! ToCompiledProcedure[sqrtcom2][[4]] still gives CompiledResult[Register[Real, 1]] and sqrtcom2 still gives errors.

Why is it not working?

$\endgroup$
  • $\begingroup$ I doubt the type inferencing will do anything with built-in functions. Just declare the input as _Complex, e.g. In[53]:= sqrtcom1 = Compile[{{x, _Complex}}, Sqrt[x]]; sqrtcom1[-1.1] Out[54]= 0. + 1.04880884817 I $\endgroup$ – Daniel Lichtblau Dec 22 '15 at 15:42
  • $\begingroup$ @DanielLichtblau I know _Complex works. So just as halirutan says in the link provided by Jason B. " one of the situations where the third argument of Compile is required is when your compiled code needs to make a MainEvaluate call". Can I safely change the word "one of" as "only one"? Is there any other situations? $\endgroup$ – matheorem Dec 22 '15 at 15:50
  • $\begingroup$ I timed the different solutions using RepeatedTiming on the root of -4 and found uncompiled Sqrt fastest, followed by sqrtHal and sqrtHal2 accepting real input, then sqrtcom1 accepting complex input, and sqrtcom3 involving a MainEvaluate call. Curiously, sqrtHal was slightly faster than sqrtHal2, but it became clear when I looked at CompilePrint that sqrtHal2 involves an additional complex addition in x + 0. I, which sqrtHal lacks. $\endgroup$ – obsolesced Jul 7 '16 at 17:51
  • $\begingroup$ Regarding why sqrtcom2 doesn't work, it seems, as you may already suspect, Compile ignores the third argument for supported functions. Only when a function requires a MainEvaluate call does Compile look at the third argument for the function's return type. I think the answer by @JasonB supports this conjecture. $\endgroup$ – obsolesced Jul 7 '16 at 18:10
6
$\begingroup$

Let me counter Daniel Lichtblau's answer

This has zero to do with type inferencing.

by saying, the example in the tutorial you linked is all about type inference. It is not about compiling com to make it faster. It is about helping the compiler to deduce the correct type for the expression.

You have to understand one thing: Highlevel Mathematica language is untyped, which means that it is not known upfront whether Sqrt[x] is an integer, a real, a complex or whether it stays as a general expression. It all depends on the value that x has.

Compiled code is completely different, because all variables will have a type. Either explicitly given by you, or derived/assumed by the compiler.

Therefore, your first question is not about whether or not com inside the compile will be too slow. The question is, can the compiler derive the type of com and therefore assume a correct return type.

Since this example does work correctly even without the explicit type hint, let me give a different example. Here, realf is a function that returns a real number, when the input is an integer. In Mathematica 10.3, this leads to an error message:

realf[i_] := 1.5*i;
f1 = Compile[{{i, _Integer}}, realf[i]]

Calling f1[3] will give you a warning, saying that realf will be the reason that the uncompiled version of f1 is used. If you check the f1 with CompilePrint you will find the line

I1 = MainEvaluate[ Hold[realf][ I0]]

The important part is that I1 means integer register. So because Compile has no information about realf, it assumes this call will be of type integer which is wrong. If we change the definition to

f2 = Compile[{{i, _Integer}}, realf[i], {{realf[_], _Real}}]

and check the compiled code again, we see that now a real register is used for the result of realf.

R0 = MainEvaluate[ Hold[realf][ I0]]

Therfore, f2[3] will run without message since the types are consistent. Nevertheless, realf will be an external call that is evaluated by the kernel.

What Daniel's answer is showing you is, that in the specific example of com being defined as Binomial, you can expand the call and indeed compile all instructions to gain a lot of speed.

As for your second example,

Compile[{{x, _Real}}, Sqrt[x]]

there is one additional thing to note: You call Sqrt[x] where x is the input of type Real. Therefore, the compiler deduces that you want the square-root that works on reals. It seems, not even the type hint at the end of Compile prevents this, but there is a simpler solution:

sqrtHal = Compile[{{x, _Real}}, 
  Module[{xx = 0. I},
   xx = x;
   Sqrt[xx]
  ]
]

Look what we did: we created another variable xx and by giving an initial complex value, we force the type-system to assume xx to be complex. The rest works like in most programming languages. When we assign xx=x a type-conversion takes place and xx is still complex where the imaginary part is zero and the real part is x. Furthermore, the correct complex Sqrt function is selected and therefore

sqrtHal[-1]
(* 0. + 1. I *)

works without problems.

As soon as you have understood the reasons behind this, you easily find another solution:

sqrtHal2 = Compile[{{x, _Real}},
  Sqrt[x + 0. I]
 ]
$\endgroup$
  • $\begingroup$ Thank you so much. I think your sqrtHal is equivalent to Compile[{{x, _Complex}}, right? $\endgroup$ – matheorem Dec 22 '15 at 16:23
  • 2
    $\begingroup$ @matheorem Your second question: You forget that usually our algorithms have more than a few lines. Therefore, it is easily possible that you can speed up large parts of your algorithm, but there is an instruction that cannot be compiled. If the overall runtime is still much better, even if we have some MainEvaluate calls in it, then you would still use the compiled function, right? The problem here on stackexchange is that you somehow get a wrong view, because so many people come here and try to speed up one-liners by wrapping Compile around it. This of course rarely helps. $\endgroup$ – halirutan Dec 22 '15 at 16:38
  • 1
    $\begingroup$ My remark was to the effect that the type of com was not relevant to the requirement of invoking MainEvaluate. Moreover I saw nothing to indicate that Compile was getting that particular type inference incorrect. $\endgroup$ – Daniel Lichtblau Dec 22 '15 at 16:52
  • 1
    $\begingroup$ @DanielLichtblau Yes, this is the reason I had to use a different example. Nevertheless, exactly th com example is used in the compile tutorial to explain type inference. $\endgroup$ – halirutan Dec 22 '15 at 16:58
  • 1
    $\begingroup$ For a small matrix it may still save an iota to use LinearSolve. But also note that now you are in the sort of situation @Halirutan indicated in a comment, to the effect that Compile might be advantageous for the other stuff. Explicit loops can be faster when used in Compile, for example, since iterator variables and bonds checks are done in low level code rather than via the interpreter (the main evaluator). $\endgroup$ – Daniel Lichtblau Dec 22 '15 at 17:06
5
$\begingroup$

This has zero to do with type inferencing. There is no Compile support for Binomial. This variant should show the same MainEvaluate invocation:

test2 = Compile[{x, {i, _Integer}}, x^Binomial[2*i, i]]

Let's see if there are evaluator calls:

test2 // InputForm

(* Out[40]//InputForm=
CompiledFunction[{10, 10.3, 5468}, {_Real, _Integer}, 
 {{3, 0, 0}, {2, 0, 0}, {3, 0, 1}}, {{2, {2, 0, 1}}}, {0, 4, 2, 0, 0}, 
 {{15, 1, 0, 2}, {47, Binomial, 2, 0, 2, 2, 0, 0, 2, 0, 3}, 
  {41, 263, 3, 0, 0, 2, 0, 3, 3, 0, 1}, {1}}, 
 Function[{x, i}, x^Binomial[2*i, i]], Evaluate] *)

Same story except slightly more efficient because it is a direct call to Binomial rather than to com. One can evade this, and get a modest speed improvement, by resorting to the FunctionExpand form of Binomial.

test3 = compile[{x, {i, _Integer}}, x^FunctionExpand[com[i]]] /. 
   compile -> Compile;
test3 // InputForm

(* Out[45]//InputForm=
CompiledFunction[{10, 10.3, 5468}, {_Real, _Integer}, 
 {{3, 0, 0}, {2, 0, 0}, {3, 0, 2}}, {{3.141592653589793, {3, 0, 1}}, 
  {2, {2, 0, 1}}, {1, {2, 0, 4}}, {0.5, {3, 0, 3}}}, {0, 6, 7, 0, 0}, 
 {{15, 1, 0, 2}, {41, 263, 2, 0, 1, 2, 0, 2, 2, 0, 3}, 
  {40, 59, 3, 0, 1, 3, 0, 2}, {10, 0, 4}, {13, 3, 4, 5}, 
  {40, 76, 3, 0, 5, 3, 0, 4}, {12, 4, 0, 2}, {40, 76, 2, 0, 2, 2, 0, 5}, 
  {10, 5, 5}, {40, 60, 3, 0, 5, 3, 0, 6}, {10, 3, 5}, {16, 5, 2, 4, 6, 5}, 
  {41, 263, 3, 0, 0, 3, 0, 5, 3, 0, 2}, {1}}, 
 Function[{x, i}, x^((2^(2*i)*Gamma[1/2 + i])/(Sqrt[Pi]*Gamma[1 + i]))], 
 Evaluate] *)

We get modest speed gains as we proceed.

AbsoluteTiming[Do[test[4., 3], {10^6}]]
AbsoluteTiming[Do[test2[4., 3], {10^6}]]
AbsoluteTiming[Do[test3[4., 3], {10^6}]]

(* Out[46]= {1.410492, Null}

Out[47]= {0.904184, Null}

Out[48]= {0.513575, Null} *)

One might get further speedup using Listability and maybe parallelism.

$\endgroup$
  • $\begingroup$ BTW, you could use Evaluate for Compile[{x, {i, _Integer}}, Evaluate[x^FunctionExpand[com[i]]]] $\endgroup$ – matheorem Dec 23 '15 at 1:51
2
$\begingroup$

I can't pretend to fully understand this, there are others here like Leonid or Halirutan that talked about this. I notice on this page in the section about "Type Propagation" it says

This shows that for this input, the compiler chooses the result of the square root to also be a real number. This is less general but faster.

So maybe it treats Sqrt specially? This is a workaround,

sqrt[x_] := Sqrt[x];
sqrtcom3 = Compile[{{x, _Real}},
   sqrt[x]
   , {{sqrt[_], _Complex, 0}}];
sqrtcom3[-4.]
(* 0. + 2. I *)
$\endgroup$
  • $\begingroup$ Worth mentioning this method involves a MainEvaluate call, which can be seen in CompilePrint, the OP wants to avoid. $\endgroup$ – obsolesced Jul 7 '16 at 17:07
  • $\begingroup$ Also, I'm curious why the complex type specification is {sqrt[_], _Complex, 0} instead of just {sqrt[_],_Complex}. Both produce equivalent answers, but I don't see any mention of a 3rd argument in the documentation. Is it supposed to specify a default value? $\endgroup$ – obsolesced Jul 7 '16 at 17:57
  • $\begingroup$ Ok I realized the 3rd argument of the type specification specifies the rank of the array of that type, as it does for the input type specification. Unfortunately this wasn't evident from the docs. $\endgroup$ – obsolesced Jul 8 '16 at 2:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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