The Mathematica Compiler

Rob Knapp/Wolfram Research

The Mathematica compiler compiles Mathematica programs into a set of pseudo code instructions that can be interpreted with a virtual machine (“WVM”) built into Mathematica or further processed to generate C code.  This talk will focus on using Compile to get the fastest possible speed.

Even if you have never used the Mathematica function Compile, it has probably been used “under the hood” for a variety of things you have done.

It is used in plotting functions

It is used for numerical computation functions like NDSolve, FindRoot, etc.

It is used for evaluating Table, Map and some other list constructors when possible.

It is used extensively in the part of Mathematica system code written in the Mathematica language.

In all of these cases, the compiler is used to speed up just a part of the computation rather than being used for the entirety.

The most common difficulty I see with many people’s use of the compiler is treating it like a C or Java compiler where an entire program has to be compiled.  In Mathematica the compiler is well integrated with the rest of the system via basic number types and PackedArrays so it is practical and very often more efficient to compile small numerically intensive parts of a computation and leave the rest to ordinary Mathematica evaluation.  An additional advantage of this approach is that you can concentrate on making small pieces of a computation fast, something that is much easier and more likely to produce results than trying to work with the whole computation at once.

Simple example: Longest consecutive subseqences.

This defines a function that returns the longest consecutive subsequence of 0s or 1s in a binary list.

mcs = Function[{bseq},Module[{c = 0, maxc =1, p = -1},
Table[If[elem== p, If[maxc < ++c, maxc = c], c = 1];p=elem; maxc, {elem, bseq}]]]

A good example for testing are the binary digits of π

Short[test = First[RealDigits[N[Pi, 10^6], 2]]]

It takes quite a long time to compute the longest consecutive subsequence

AbsoluteTiming[trial = mcs[test];]

This shows how the longest consecutive subsequence length increases as you go through the sequence

ListPlot[Transpose[{Log[10,N[Range[Length[trial]]]], trial}], PlotRange->{{0,6},All}]

Now lets compile the function, telling the compiler that there is an argument that should be a list of machine integers.

AbsoluteTiming[vmcs = Compile[{{bseq, _Integer, 1}}, Module[{c = 0, maxc =1, p = -1},
Table[If[elem== p, If[maxc < ++c, maxc = c], c = 1];p=elem; maxc, {elem, bseq}]]]]

This returns a CompiledFunction object.  When given arguments matching the types given to Compile, it will evaluate using the WVM

AbsoluteTiming[ctrial = vmcs[test];]

The evaluation was 30-40 times faster and the result is the same

trial == ctrial

The compilation can be taken a step further to generate C code that is then connected to Mathematica through a LibraryFunction.

AbsoluteTiming[cmcs = Compile[{{bseq, _Integer, 1}}, Module[{c = 0, maxc =1, p = -1},
Table[If[elem== p, If[maxc < ++c, maxc = c], c = 1];p=elem; maxc, {elem, bseq}]], CompilationTarget->"C"]]

At the cost of extra compilation time (usually most of this is the C compiler optimization) we get another factor of 5 speedup.

AbsoluteTiming[cmcs[test];]

Even faster yet is to remove checks on integer overflow (certainly not needed here)

cmcs = Compile[{{bseq, _Integer, 1}}, Module[{c = 0, maxc =1, p = -1},
Table[If[elem== p, If[maxc < ++c, maxc = c], c = 1];p=elem; maxc, {elem, bseq}]], CompilationTarget->"C", RuntimeOptions->"Speed"];

AbsoluteTiming[cmcs[test];]

CompiledFunctionTools

The CompiledFunctionTools application contains functions that make it much easier to understand what a CompiledFunction is doing.

Needs["CompiledFunctionTools`"]

The FullForm of a CompiledFunction contains sufficient information for Mathematica to use the WVM if it is saved and restored in another session.  It is not meant to be comprehensible or constructable except through Compile:

FullForm[vmcs]

The CompiledFunctionTools provides a function CompilePrint that can be very useful in seeing what a CompiledFunction actually does:

CompilePrint[vmcs]

Compiler Types

The output of CompilePrint is easier to understand if you are aware of the types the compiler supports:

True | False (B) Boolean values

_Integer (I) machine integer values

_Real (R) machine real values.  This is the default used if the type is not specified.

_Complex (C) machine complex values, i.e. a point in the complex plane with 2 machine reals

Null (V) “Void” or non-values.  This cannot be used as an argument type.

Tensor (T)  list values.

Tensors are arrays with rectangular dimensions and numbers (Integer, Real, or Complex) of all the same types. Outside the scope of the compiler, tensors are exported as PackedArrays.  The rank is the length of the dimensions.

Specified in arguments by giving the rank, e.g. {bseq, _Integer, 1} means a tensor of integers with a single dimension, in other words a vector.

The rank must be known at compile time

Mathematica “external” evaluation

When the compiler gets an expression it cannot translate into instructions for the virtual machine, it is handled by doing a Mathematica evaluation of the expression called from the virtual machine.

As an example, consider a function that bins values hit when iterating a discrete map of the unit interval to itself.

binnest = Compile[{{r, _Real}, {n0, _Integer}, {n, _Integer}, {bins, _Integer}}, 
Module[{x = 1./Pi, mb = bins - 1., rounded},
Do[x = f[r,x],{n0}];
rounded = Table[x = f[r,x];Ceiling[mb*x],{n}];
Tally[rounded]]];

There is no way that Compile can know about f at this point, since it is undefined, so to resolve it, the interpreter must do a Mathematica evaluation each time the value is needed at run time.  You can see this easily using CompilePrint

CompilePrint[binnest]

Something else you can see is that the expected type of f[r, x] is real.  The type was inferred from the types of the arguments and is really only a guess.

It is important that the type be right from the MainEvaluate or the CompiledFunction will not be able to run.  If you turn on the “CompilerWarnings” messge group, the compiler will issue a message whenever it has to try to guess the type:

On["CompilerWarnings"];

binnest = Compile[{{r, _Real}, {n0, _Integer}, {n, _Integer}, {bins, _Integer}}, 
Module[{x = 1./Pi, mb = bins - 1., rounded},
Do[x = f[r,x],{n0}];
rounded = Table[x = f[r,x];Ceiling[mb*x],{n}];
Tally[rounded]]];

You can give Compile a third argument specifying types for things that are will need to be resolved by evaluation.  
Here we use {_f, _Real} to mean anything matching the pattern _f (i.e. with Head f) is expected to be of real type.

binnest = Compile[{{r, _Real}, {n0, _Integer}, {n, _Integer}, {bins, _Integer}}, 
Module[{x = 1./Pi, mb = bins - 1., rounded},
Do[x = f[r,x],{n0}];
rounded = Table[x = f[r,x];Ceiling[mb*x],{n}];
Tally[rounded]], {{_f, _Real}}];

The advantage of having an evaluation here is that the function can be changed without recompiling and as long as the function returns the type you gave, it will just work.

This uses the discrete logistic map.

AbsoluteTiming[Block[{f = Function[{r, x}, r x (1. - x)]}, binest[3.5, 10^5, 10^5, 1000]]]

This uses the tent map

AbsoluteTiming[Block[{f = Function[{r, x}, r (1. - Abs[1. - 2. x])]}, binest[.75, 10^5, 10^5, 1000]]]

While there is some advantage to being able to change the function so easily, that comes at a price.  If the function is known at compile time, then the computation can be done much faster.  Using the option CompilationOptions->{“InlineExternalDefinitions”->True} tries to get current definitions and uses them for the compilation.  With this, it is not necessary to specify the type for f because the compiler has it explicilty.

logistic = Block[{f = Function[{r, x}, r x (1. - x)]},Compile[{{r, _Real}, {n0, _Integer}, {n, _Integer}, {bins, _Integer}}, 
Module[{x = 1./Pi, mb = bins - 1., rounded},
Do[x = f[r,x],{n0}];
rounded = Table[x = f[r,x];Ceiling[mb*x],{n}];
Tally[rounded]], 
CompilationOptions->{"InlineExternalDefinitions"->True}]];

Now its more than a lilttle bit faster!

AbsoluteTiming[logistic[3.5, 10^5, 10^5, 1000]]

Conclusion: If you can avoid having Mathematica evaluations in your CompiledFunction, it will typically run a lot faster.

RuntimeAtttributes->Listable

Now suppose you want to evaluate the logistic function at several values of the parameter r, say to make a bifurcation diagram.  One way is to use Map

rtest = Table[r, {r, 3.,4.,.001}];

AbsoluteTiming[data1 = Map[logistic[#, 10^5, 10^5, 1000]&, rtest];]

Another way is to allow it to work for multiple values of r by giving it the Listable runtime attribute.

logistic = Block[{f = Function[{r, x}, r x (1. - x)]},Compile[{{r, _Real}, {n0, _Integer}, {n, _Integer}, {bins, _Integer}}, 
Module[{x = 1./Pi, mb = bins - 1., rounded},
Do[x = f[r,x],{n0}];
rounded = Table[x = f[r,x];Ceiling[mb*x],{n}];
Tally[rounded]],CompilationOptions->{"InlineExternalDefinitions"->True}, RuntimeAttributes->Listable]];

AbsoluteTiming[data2 = logistic[rtest, 10^5, 10^5, 1000];]

This works very much like the Listable attribute you can give to a symbol

SetAttributes[g, Listable];
g[r_, x_] := If[x < 1/2, 2 r x, 2r(1- x)]

g[{0,1/4,1/2,1}, x]

g[r,{0,1/4,1/2,1}]

Parallelization

The Parallelization option directs how to handle Listable evaluation.  If Parallelization->True, then the separate instances of the function will be evaluated in $ProcessorCount separate threads in parallel.

logistic = Block[{f = Function[{r, x}, r x (1. - x)]},Compile[{{r, _Real}, {n0, _Integer}, {n, _Integer}, {bins, _Integer}}, 
Module[{x = 1./Pi, mb = bins - 1., rounded},
Do[x = f[r,x],{n0}];
rounded = Table[x = f[r,x];Ceiling[mb*x],{n}];
Tally[rounded]], 
CompilationOptions->{"InlineExternalDefinitions"->True}, RuntimeAttributes->Listable, Parallelization->True]];

$ProcessorCount

AbsoluteTiming[data2 = logistic[rtest, 10^5, 10^5, 1000];]

logistic = Block[{f = Function[{r, x}, r x (1. - x)]},Compile[{{r, _Real}, {n0, _Integer}, {n, _Integer}, {bins, _Integer}}, 
Module[{x = 1./Pi, mb = bins - 1., rounded},
Do[x = f[r,x],{n0}];
rounded = Table[x = f[r,x];Ceiling[mb*x],{n}];
Tally[rounded]], 
CompilationOptions->{"InlineExternalDefinitions"->True}, RuntimeAttributes->Listable, Parallelization->True, CompilationTarget->"C"]];

CompilationTarget

The default compilation target is instructions for the (WVM) virtual machine.  However, if you have a function that you think will be used many times and speed is important, you can compile to C code that is then connected to Mathematica as LibraryFunction.

Compilation to C code typically takes MUCH more time than for the virtual machine.  Most of the extra cost is typically in the optimization phase of the C compiler.

AbsoluteTiming[logistic = Block[{f = Function[{r, x}, r x (1. - x)]},Compile[{{r, _Real}, {n0, _Integer}, {n, _Integer}, {bins, _Integer}}, 
Module[{x = 1./Pi, mb = bins - 1., rounded},
Do[x = f[r,x],{n0}];
rounded = Table[x = f[r,x];Ceiling[mb*x],{n}];
Tally[rounded]], 
CompilationOptions->{"InlineExternalDefinitions"->True}, RuntimeAttributes->Listable, Parallelization->True, CompilationTarget->"C"]];]

This runs the compiled C code in parallel:

AbsoluteTiming[data2 = logistic[rtest, 10^5, 10^5, 1000];]

Here is the C code:

Needs["CCodeGenerator`"];
CCodeStringGenerate[logistic,"logistic", CodeTarget->"WolframLibrary"]

RuntimeOptions

In some cases, the C code generated can be much faster if less error checking is done and more of the code can be inlined (as opposed to using function calls).  The easiest way to turn off error checking is to choose RuntimeOptions->”Speed”.  RuntimeOptions->”Quality” does the most error checking and the result is most likely to be the same as what you would get with ordinary Mathematica evaluation.

In the logistic function example, changing this option does not make any difference, however, in other examples it can make a profound difference.

Now that we have a function that runs quite fast, it is possible to make a very nice looking bifurcation diagram.

BifurcationDiagram[f_, {rmin_, rmax_, nr_}, nb_, n_] := BifurcationDiagram[f, {rmin, rmax, nr}, nb, {n,n}];BifurcationDiagram[f_, {rmin_, rmax_, nr_}, nb_, {ns_, nt_}] := 
Module[{mapper, data, r, S},
mapper = Compile[{{r, _Real}, {n0, _Integer}, {n, _Integer}, {bins, _Integer}}, 
Module[{x = 1./Pi, mb = bins - 1., rounded},
Do[x = f[r,x],{n0}];
rounded = Table[x = f[r,x];Floor[mb*x],{n}];
Tally[rounded]], 
 RuntimeAttributes->Listable, Parallelization->True(*,CompilationTarget->"C"*)];
data = mapper[Table[r, Evaluate[{r, N[rmin], N[rmax], N[(rmax - rmin)/(nr- 1)]}]], ns, nt, nb];
S = SparseArray[Table[{bins, freq} = Transpose[d]; SparseArray[Partition[nb - bins, 1]->Log[N[freq]],{nb}],{d,data}], {nr, nb}];
ArrayPlot[Transpose[S], ColorFunction->Function[If[# == 0, Black, Hue[#]]], ImageSize->{nr, nb}]]

BifurcationDiagram[Function[{r, x}, r x (1. - x)], {3,4,512}, 512, {10^4, 10^5}]

BifurcationDiagram[Function[{r, x}, r (1. - Abs[1. - 2. x])], {.5,1.,512}, 512, {10^4, 10^5}]

BifurcationDiagram[Function[{r, x},  Sin[π r x(1-x)]], {1,4,512}, 512, {10^4, 10^5}]

CompiledFunction calls

You may have noticed that I did not use the option CompilationTarget->”C” in the definition of BifurcationDiagram above.  This is because unless you wanted to do a very large plot, the compilation time would take as much or more time than the rest of the computation.  CompilationTarget->”C” is only the best choice if you are going to use a function many times or the computation does take a very long time.

When Compile encounters another CompiledFunction or LibraryFunction, it is able to use special instructions to call those functions very efficiently, particularly in the case of LibraryFunction.  This means that if the function given to BifurcationDiagram is a CompiledFunction created with CompilationTarget->”C”, the mapper function will call this function efficiently.

AbsoluteTiming[clogistic = Compile[{{r, _Real}, {x, _Real}}, r x(1.-x), CompilationTarget->"C"];]

calllogistic = Block[{f = clogistic},Compile[{{r, _Real}, {n0, _Integer}, {n, _Integer}, {bins, _Integer}}, 
Module[{x = 1./Pi, mb = bins - 1., rounded},
Do[x = f[r,x],{n0}];
rounded = Table[x = f[r,x];Ceiling[mb*x],{n}];
Tally[rounded]], 
CompilationOptions->{"InlineExternalDefinitions"->True}, RuntimeAttributes->Listable, Parallelization->True]];

AbsoluteTiming[data2 = logistic[rtest, 10^5, 10^5, 1000];]

The way BifurcationDiagram is defined, the function will be used directly in the mapper function.

AbsoluteTiming[BifurcationDiagram[Function[{r, x}, r x (1. - x)], {3,4,512}, 512, {10^4, 10^5}];]

AbsoluteTiming[BifurcationDiagram[clogistic, {3,4,512}, 512, {10^4, 10^5}];]

Compilation on demand

Sometimes if you want to write a general code, you don’t know ahead of time exactly what types will be used.  To get around this, you can define functions that will compile for appropriate types as needed.

Clear[escape];
escape[z_] := Module[{i = 0, x = 0}, While[Abs[x] < 2 && i < 1000,x = x*x + z; i++]; i] /; NumberQ[z]

Timing[Plot[escape[x],{x,0.25,.26}, PlotRange->All]]

zero[Real] = 0.;
zero[Complex] = 0. + 0. I;

Clear[cescape];cescape[type_] := (cescape[type] = With[{zt = zero[type]},Compile[{{z, Blank[type]}}, Module[{i = 0, x =zt}, While[Abs[x] < 2 && i < 1000,x = x*x + z; i++]; i], CompilationOptions->{"InlineExternalDefinitions"->True}, RuntimeAttributes->Listable, Parallelization->True, RuntimeOptions->"Speed",CompilationTarget->"C"]])

Timing[Plot[cescape[Real][x],{x,0.25,.26}, PlotRange->All]]

ArrayOf[type_][v_] := ArrayQ[v, _, MatchQ[#, Blank[type]]&];

Clear[escape];
escape[z_Real] := cescape[Real][z];
escape[z_Complex] := cescape[Complex][z];
escape[z_?(ArrayOf[Real])] := cescape[Real][z];
escape[z_?(ArrayOf[Complex])] := cescape[Complex][z];
escape[z_] := Module[{i = 0, x = 0}, While[Abs[x] < 2 && i < 1000,x = x*x + z; i++]; i] /; NumberQ[z];
escape[z_] := Map[escape, z, {ArrayDepth[z]}] /; ArrayQ[z, _, NumberQ];

Timing[Plot[escape[x],{x,.25,.26}]]

AbsoluteTiming[data = escape[Table[x + (0. + 1. I) y,{y,-1.25,1.25, .005}, {x,-2, .5, .005}]];]

ArrayPlot[data, AspectRatio->1, ImageSize->Dimensions[data]]

Some ideas to take with you

Be aware of types

Use CompiledFunctionTools and PrintCompile to see what a CompiledFunction is doing

Speed things up.  Here are things to try:

Try to avoid Mathematica evaluations

Check that types are as consistent as possible

If you want to run many instances think about using RuntimeAttributes->Listable possibly with Parallelization->True

If you expect to use the function many times, try CompilationTarget->”C”

For lots more information, go to the The Mathematica Compiler User Guide

Here is what I think is the most important thing of all:  For most computations, you’ll have far better success in using the compiler to speed up parts of the computation rather than trying to use Compile to generate one function that does it all.