xtensor features offer alternative to NumPy in C++

Python, and especially NumPy, has been a cornerstone in scientific computing for quite a while now. Initially set out to reimplement the most valuable features of Matlab, NumPy might already be more popular in a couple of science fields. While NumPy is quite fast, it remains bound to the Python interpreter, which is working with an incredibly dynamic language — slowing it down. On the opposite, C++ doesn’t have this problem: You only wait once, and that is for the initial compilation to take place!

And boy, C++ has caught up, in terms of language expressiveness. C++11 and C++14 delivered heaps of language goodness. In this blog post, I will introduce you to xtensor, a new framework which aims to replicate NumPy’s features in modern C++. And it uses a lot of template metaprogramming C++14 goodness under the hood, making it incredibly easy to use.

While there exists a number of well-established linear algebra libraries for C++, like Eigen or Armadillo, they are all limited to two or 3 matrix dimensions. And of course there are various implementations of BLAS for C++, but all of them are pretty barebones quite inconvenient to use, requiring working with bare pointers and matrix strides. xtensor is a clear departure from that, rethinking the linear algebra library from the ground up. For this reason, there is also no BLAS library currently in use at xtensor! This might seem counter-intuitive, but we plan to create our own SIMD accelerated routines.

Anyway, to get to the meat of the article, let’s demonstrate a few features of xtensor:

xt::xarray a = {{1,2,3}, {4,5,6}, {7,8,9}};
auto b = xt::arange(0, 1000000000);
xt::xtensor c = xt::view(b, range(_, 3));
auto d = a * c + 3;
std::cout << d << std::endl;

What is happening here!? In the first line, we instantiate a container of the `xarray` type which holds a bunch of doubles. An xarray is a container with a dynamic number of dimensions. This is an important distinction to the xtensor. An xtensor has a fixed number of dimensions. The xarray is filled with values from an initiializer_list (a C++11 feature), and the shape is deduced from the initializer list as well (3, 3).

On the second line we see an `xgenerator`. Everybody is probably familiar with the NumPy function `arange`. Here we instantiate an `arange` from 0 to 1000000000. That’s A LOT of memory we’d be wasting if we’d do that in NumPy. In xtensor, `xgenerators` are only evaluated upon access. That means, the `arange` function costs almost nothing in storage and generates the required values on the fly.

In the next line we instantiate an xtensor. Note the second template parameter (`1`) specifying the number of dimensions this xtensor has. This is necessary as the xtensor has a stack allocated shape, which usually is more performant than the dynamic shape of the xarray. We also create a view into the xgenerator! The view takes the values from `_` to `3`. The underscore is another new C++ feature – a universal placeholder. In this case, it works like the empty colon in Python. The view above is the same as taking the first three elements, and equivalent to the following NumPy: `c = b[:3]`.

Now that we have two xexpressions (a xtensor and a xarray), we can use them to do some mathematics. A simple multiplication with an addition is performed in line 3. You might have noticed, that this expression contains 3 different shapes: (3, 3), (3,) and a scalar. Again, xtensor got inspired by NumPy’s broadcasting rules and applies the same.

But unlike NumPy, the mathematical calculation is only executed upon access. So you could multiply incredibly large arrays without clogging your memory, even if you only need 3 specific entries – as they will only get calculated upon access.

Well, the last line is quite standard if you know a bit of C++, just printing the `xexpression` to the command line. If you do so, you’ll notice that the developers spend quite some time getting the output right. Just like NumPy, it’s properly indented, and line breaks are inserted etc.

This was just a quick introduction to xtensor to get you excited. In the meantime, you should definitly have a look at the xtensor to NumPy cheatsheet to find out which features of NumPy are already available in xtensor (and more are added approximately every 2 weeks): https://xtensor.readthedocs.io/en/latest/numpy.html. In an upcoming post I’ll show you how to implement a custom `xfunction`, and a custom `xreducer`. So stay tuned.

©ODSC2017