## Mathematica Programming » Performance Tuning

Mathematica is what is often called a high-level programming language, because when you are writing programs, you don't need to worry about how the computer is doing its job. This provides notable gains in simplicity and scalability of large programs, but it comes at the cost of the speed and efficiency of low-level computations.

## Basics of Compilation

To get around this Mathematica provides the function `Compile`

which will generate a low-level function that we can then apply with fantastic results. For example, in the section on `Nest`

I provided a random vector generator:

This creates a compiled function that will perform the task at hand, in this case taking a vector of `Real`

numbers and returning a vector. Just so we can see the difference let's create the same function, uncompiled and compare the two:

` ````
randomVector=
Compile[{{mag,_Real}},
mag*(#/Norm[#])&@RandomReal[{-1,1},3]
]
```

` ````
randomVectorU=
Function[mag,
mag*Normalize@RandomReal[{-1,1},3]
];
```

Then we'll compare them using a little timing function:

` ````
compareTiming[f1_,f2_,its_:1]:=
With[{
tf1=First@AbsoluteTiming[Do[f1,its]],tf2=First@AbsoluteTiming[Do[f2,its]]},
(tf2-tf1)/Max@{tf1,tf2}
];
compareTiming~SetAttributes~HoldAll;
compareTiming[
randomVector[1],
randomVectorU[1],
10^6
]
```

` ``0.5244177474084153``

So we get a 52% speed boost from essentially just changing `Function`

to `Compile`

on the simplest function imaginable. If we had a more involved function the speed boost would likely be even larger. But since it's so simple we can easily dissect what is giving the speed boost:

### Analyzing Compile

We'll build from the bottom up:

` ````
randomReal3=Compile[{},RandomReal[{-1,1},3]];
randomReal3U=Function[RandomReal[{-1,1},3]];
compareTiming[randomReal3[], randomReal3U[],10^6]
```

` ``0.4059389866935647``

This alone is 40% faster! Of course when we know how many vectors we'll need ahead of time the fastest thing is to use the syntax `RandomReal[range,{count,3}]`

which can, for example, generate `10^8`

random vectors in a fraction of a second, which is much, much faster than our compiled form:

` ``compareTiming[Do[randomReal3[], 10^7], RandomReal[{-1, 1}, {3, 10^7}]]`

` ``-0.9477078222470973``

Ignoring the correct way to do this, though, and considering out surprising result, this is likely an artifact of type-checking in the compiled function that has to be performed at each call in the uncompiled one.

When we add our normalization now:

` ````
randomVectorNoMag=
Compile[{},
(#/Norm[#]&)@RandomReal[{-1,1},3]
];
randomVectorNoMagU=
Function[
Normalize@RandomReal[{-1,1},3]
];
compareTiming[randomVectorNoMag[], randomVectorNoMagU[], 10^6]
```

` ``0.4923576135986387``

Now we're a full 50% faster than we were (that this is greater than our original 40 is likely just a lack of consistency between runs).

Note that I'm in the `CompiledFunction`

I'm not using `Normalize`

. This is because `Normalize`

cannot be directly compiled.

We can see this by looking at the intermediate representation generated when we use it:

` ````
<<CompiledFunctionTools`
CompilePrint@
Compile[{},
Normalize@RandomReal[{-1,1}, 3]
]//Snippet[#, -5]&
```

` ````
1 R0 = I1
2 R1 = I2
3 T(R1)1 = RandomReals[ R0, R1, T(I1)0]]
4 T(R1)2 = MainEvaluate[ Hold[Normalize][ T(R1)1]]
5 Return
```

The call to `"MainEvaluate"`

means that the compiled function couldn't handle this and thus it returns the value to the main evaluator to process.

### Drastic Examples

The previous case was simple so that it could easily be analyzed, but for more complicated cases we can much bigger gains. Consider the primary gradient function we used in the section on `FixedPoint`

:

` ````
gradient=
Compile[{{pos,_Real,1}},
{
Piecewise[{
{Sqrt@-(#+2),#<-3},
{(-(#+2)^3),-3<#<0},
{-((#-2)^3),0<=#<3},
{-Sqrt@(#-2),3<=#}
},#]&@pos[[1]],
Piecewise[{
{Sqrt@-(#),#<-1},
{(-(#)^3),-1<=#<=1},
{-Sqrt@(#),1<=#}
},#]&@pos[[2]],
Piecewise[{
{Sqrt@-(#),#<-1},
{(-(#)^3),-1<=#<=1},
{-Sqrt@(#),1<=#}
},#]&@pos[[3]]
}
];
```

And then consider its uncompiled cousin:

` ````
gradientU=
Function[pos,
{
Piecewise[{
{Sqrt@-(#+2),#<-3},
{(-(#+2)^3),-3<#<0},
{-((#-2)^3),0<=#<3},
{-Sqrt@(#-2),3<=#}
},#]&@pos[[1]],
Piecewise[{
{Sqrt@-(#),#<-1},
{(-(#)^3),-1<=#<=1},
{-Sqrt@(#),1<=#}
},#]&@pos[[2]],
Piecewise[{
{Sqrt@-(#),#<-1},
{(-(#)^3),-1<=#<=1},
{-Sqrt@(#),1<=#}
},#]&@pos[[3]]
}
];
With[{grid=RandomReal[{-5,5},{10^5,3}]},
compareTiming[
gradient@g,
gradientU@g,
{g,grid}]
]
```

` ``0.8961919038258733``

So the simple application of `Compile`

provides a 90% speed boost.

As another example, one can implement numerical integration of a volume over a surface via:

` ````
compDotCross=
Compile[{{coords, _Real, 2}, {abc,_Integer,1}},
coords[[ abc[[1]] ]].
Cross[coords[[ abc[[2]] ]], coords[[ abc[[3]] ]]]
];
vol[s_?RegionQ]:=
With[
{coords = MeshCoordinates[s]},
(MeshCells[s, 2]/.Polygon[g_List] :> compDotCross[coords, g])
//Total
]/6
```

then get an interesting surface:

` ````
surf=
RegionUnion@@
With[
{
graphics=
FirstCase[
ChemicalData["Caffeine", "SpaceFillingMoleculePlot"],
_GraphicsComplex,
None,
∞
]
},
Cases[
graphics[[2]],
Sphere[s_, r_]:>Ball[graphics[[1, s]], r]
]
]//BoundaryDiscretizeRegion[#, MaxCellMeasure->{1->20}]&
```

And compute its volume:

` ``vol@surf//RepeatedTiming`

` ````
{0.0989207000000000003`2.,1.6728430521905467`*^8}
Volume@surf//RepeatedTiming
```

` ``{0.0010948070175438598`2.,1.6728430521905503`*^8}`

Our little compiled implementation is much slower than the built-in implementation, of course, but it's still incredibly fast for the level of mesh refinement.

Personally, I this it for computing atomic orbital volumes derived from orbital contour plots.

Lets then to a version without the compiled dot-cross:

` ````
uncompDotCross=
Function[{coords, abc},
coords[[ abc[[1]] ]].
Cross[coords[[ abc[[2]] ]], coords[[ abc[[3]] ]]]
];
volU[s_?RegionQ]:=
With[
{coords = MeshCoordinates[s]},
(MeshCells[s, 2]/.Polygon[g_List] :> uncompDotCross[coords, g])
//Total
]/6
compareTiming[vol@surf, volU@surf]
```

` ``0.9539375762845481``

Ours is a whopping 95% faster. As the surface refinement increases this can be the difference between a function being nearly instantaneous versus a major speed bump in the research process.

We can also compare the version where we compile to C (i.e. we use `CompilationTarget->"C"`

` ````
compDotCross=
Compile[{{coords, _Real, 2}, {abc,_Integer,1}},
coords[[ abc[[1]] ]].
Cross[coords[[ abc[[2]] ]], coords[[ abc[[3]] ]]],
CompilationTarget->"C"
];
volC[s_?RegionQ]:=
With[
{coords = MeshCoordinates[s]},
With[{compiledDotCross=compDotCross},
(MeshCells[s, 2]/.Polygon[g_List] :> compiledDotCross[coords, g])
//Total
]/6
]
compareTiming[volC@surf, vol@surf]
```

` ``0.13463908089141324``

And simply compiling to C made the function 13% faster. More performance can be eked but by using the difference `RuntimeOptions`

and things.

### Using Compile

It's often unclear whether a given function can be compiled or not, but happily Mathematica provides all of the possible compilable functions in `Compile`CompilerFunctions`

. We can use this to make a `compileableQ`

type function:

` ````
compilableQ[f_]:=
MemberQ[Compile`CompilerFunctions[], f]
```

For example, we can see that, say, `RandomVariate`

is compilable:

` ``compilableQ@RandomVariate`

` ``True`

Alternatively we can figure out which `Random⋆`

functions are compilable:

` ````
compilableSelect[p_]:=
Select[Compile`CompilerFunctions[], StringMatchQ[ToString@#,p]&]
compilableSelect@"Random*"
```

` ``{RandomChoice,RandomSample,RandomInteger,RandomVariate,RandomComplex,RandomReal,Random}`

Once we have that, `Compile`

has a very simple call structure that is similar to `Module`

:

` ````
Compile[
{
{var1,type1,arraySize1},
{var2,type2,arraySize2},
...
},
expression,
ops
]
```

This will generate a `CompiledFunction`

which is just a function in either Mathematica byte code or C byte code, as specified by the `CompilationTarget`

option. In most cases Mathematica byte code is more than fast enough, but there are cases where the C code is a significantly better choice.