Accelerating Mathematica with C++ libraries

I love Mathematica. It's a great tool for prototyping ideas and quickly figuring out which ones have merit and which ones don't. However, that convenience comes at a price: a lot of code written in Mathematica runs considerably slower than in a compiled language like C++. But how much slower is "considerably slower"?

For example, in a recent project, I was using Mathematica to process finite element meshes and I needed to calculate element qualities. So, I implemented a common quality metric ($q=\frac{6}{\sqrt{2}}\frac{V}{{\ell }_{rms}^3}$$q = \frac{6}{\sqrt{2}} \frac{V}{\ell_{rms}^3}$) in Mathematica:

quality[x_] := Module[{e1, e2, e3, e4, e5, e6, v, Lrms},
  e1 = x[[2]] - x[[1]];
  e2 = x[[3]] - x[[1]];
  e3 = x[[4]] - x[[1]];
  e4 = x[[4]] - x[[2]];
  e5 = x[[4]] - x[[3]];
  e6 = x[[3]] - x[[2]];
  v = Det[{e1, e2, e3}]/6.0;
  Lrms = Sqrt[(e1.e1 + e2.e2 + e3.e3 + e4.e4 + e5.e5 + e6.e6)/ 6.0];
  6 Sqrt[2] v/(Lrms*Lrms*Lrms)
]

It took Mathematica about 16 seconds to compute qualities of the elements in a modestly sized mesh of ~1000000 tetrahedra. This seemed unusually slow, so I decided to implement the same calculation in C++ and measure the performance relative to Mathematica

xxxxxxxxxx
double volume(const mat<4,3> & tet) {
  auto e1 = tet[1] - tet[0];
  auto e2 = tet[2] - tet[0];
  auto e3 = tet[3] - tet[0];
  return dot(cross(e1, e2), e3) / 6.0;
}

double quality(const mat<4,3> & tet) {
  double L_rms = sqrt((norm_squared(tet[1] - tet[0]) + 
                       norm_squared(tet[2] - tet[1]) + 
                       norm_squared(tet[0] - tet[2]) + 
                       norm_squared(tet[0] - tet[3]) + 
                       norm_squared(tet[1] - tet[3]) + 
                       norm_squared(tet[2] - tet[3])) / 6.0); 
  double V = volume(tet);
  return 6.0 * sqrt(2.0) *  V / (L_rms * L_rms * L_rms);
}

The timings for a 1000000 element mesh comparison are:

So, depending on the number of threads, the C++ implementation is 3-4 orders of magnitude faster than the native Mathematica one. Using Mathematica's built in Compile[] function helps performance, but not very much. In my experience Compile[] is also pretty frustrating to use, since the rules of what you can and can't do in a compiled function aren't documented, so the compiled function frequently ends up resorting to the native (i.e. slow) Mathematica implementation.

Luckily, there are a few ways to write C++ extensions for Mathematica.

WSTP (MathLink) and LibraryLink

These options let you write C++ libraries and use them from Mathematica in theory, but they're unpleasant to use in practice. They involve a lot of boilerplate in the C++ implementation, the examples aren't helpful and they don't seem to integrate nicely into realistic C++ projects that have dependencies and build systems. So, despite my trying to use these tools for years, I could never manage to find a workflow that actually made it worth the time investment.

Wolfram seems to understand that these tools are unpleasant to use, so they created tools to try and make them easier to wield (e.g. LibraryLink Utilities, "LLU"), but even these convenience tools are themselves clumsy.

CCompilerDriver package

Another option is to use the existing CCompilerDriver package from within Mathematica. This lets you pass strings containing C++ source code to CreateLibrary[] to compile it to a binary file that can be linked into Mathematica. To demonstrate how easy it is to use, their own documentation shows an example of writing a function that adds 1 to its input:

add1src = "
#include \"WolframLibrary.h\"

DLLEXPORT mint WolframLibrary_getVersion(){
  return WolframLibraryVersion;
}

DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) \
{
  return 0;
}

DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData \
libData) {
  return;
}

DLLEXPORT int constantzero(WolframLibraryData libData, mint Argc, \
MArgument *Args, MArgument Res){
   MArgument_setInteger(Res, 0);
   return LIBRARY_NO_ERROR;
}

DLLEXPORT int add1(WolframLibraryData libData,
      mint Argc, MArgument *Args, MArgument Res) {
  mint I0;
  mint I1;
  I0 = MArgument_getInteger(Args[0]);
  I1 = I0 + 1;
  MArgument_setInteger(Res, I1);
  return LIBRARY_NO_ERROR;
}

";
  
add1lib = CreateLibrary[add1src, "add1"]
add1 = LibraryFunctionLoad[add1lib, "add1", {Integer}, Integer]

As you can see, defining a simple function with this approach also requires a lot of setup and boilerplate!

Further still, that C++ code is a string inside Mathematica, so you have to escape quotations, you get no syntax highlighting, no autocompletion, no intellisense, no vim/emacs bindings, etc. It's completely impractical to use this with an actual C++ project with multiple sourcefiles and external dependencies.

wll-interface

Thankfully, there's a great github repo here: https://github.com/njpipeorgan/wll-interface/tree/master

It provides an elegant, minimal wrapper around the LibraryLink interface. Instead of 40 lines of boilerplate, the add1 function from earlier with wll-interface is:

#include "wll_interface.h"
double add1(double x){ return x + 1; }
DEFINE_WLL_FUNCTION(add1)

which is about as minimal as one could possibly hope for. It also defines easy-to-use container classes for the relevant Mathematica data types (e.g. tensors, sparse arrays) that feel like using STL containers. Perhaps best of all, it integrates seamlessly into existing C++ projects. With just a couple of lines of CMake, wll-interface can be added to an existing project:

include(FetchContent)
FetchContent_Declare(
  wll_interface
  GIT_REPOSITORY https://github.com/samuelpmish/wll_interface.git
  GIT_TAG main
)
FetchContent_MakeAvailable(wll_interface)

...

# needs to be a shared library for Mathematica to load it
add_library(my_library SHARED ...)
target_link_libraries(my_library PUBLIC wll_interface)

From there, we can build the shared library and load its symbols into Mathematica with LibraryFunctionLoad[]. Here are the timings of the wll-interface approach, compared to the previous results:

So, there's definitely still some overhead when going through LibraryLink, but a speedup of 1000x for practically no extra effort is worth it to me!

Summary

wll-interface is an incredible library, I'm really grateful for njpipeorgan's work developing and maintaining it. It makes writing C++ extensions for Mathematica really simple, which addresses the biggest drawback associated with working in Mathematica. This way, the C++ library benefits not only the users of that library, but also the prototyping of future features in Mathematica.