Expression templates

Expression templates are a C++ template metaprogramming technique that builds structures representing a computation at compile time, which structures are evaluated only as needed to produce efficient code for the entire computation.^[1] Expression templates thus allow programmers to bypass the normal order of evaluation of the C++ language and achieve optimizations such as loop fusion.

Expression templates were invented independently by Todd Veldhuizen and David Vandevoorde;^[2]:341[3] it was Veldhuizen who gave them their name.^[2]:341 They are a popular technique for the implementation of linear algebra software.^[1]

Motivation and example

Consider a library representing vectors and operations on them. One common mathematical operation is to add two vectors u and v, element-wise, to produce a new vector. The obvious C++ implementation of this operation would be an overloaded operator+ that returns a new vector object:

class Vec {
    std::vector<double> elems;

  public:
    Vec(size_t n) : elems(n) {}

    double &operator[](size_t i)      { return elems[i]; }
    double operator[](size_t i) const { return elems[i]; }
    size_t size()               const { return elems.size(); }
};

Vec operator+(Vec const &u, Vec const &v) {
    assert(u.size() == v.size());
    Vec sum(u.size());
    for (size_t i = 0; i < u.size(); i++) {
         sum[i] = u[i] + v[i];
    }
    return sum;
}

Users of this class can now write Vec x = a + b; where a and b are both instances of Vec.

A problem with this approach is that more complicated expressions such as Vec x = a + b + c are implemented inefficiently. The implementation first produces a temporary vector to hold a + b, then produces another vector with the elements of c added in. Even with return value optimization this will allocate memory at least twice and require two loops.

Delayed evaluation solves this problem, and can be implemented in C++ by letting operator+ return an object of a custom type, say VecSum, that represents the unevaluated sum of two vectors, or a vector with a VecSum, etc. Larger expressions then effectively build expression trees that are evaluated only when assigned to an actual Vec variable. But this requires traversing such trees to do the evaluation, which is in itself costly.^[4]

Expression templates implement delayed evaluation using expression trees that only exist at compile time. Each assignment to a Vec, such as Vec x = a + b + c, generates a new Vec constructor if needed by template instantiation. This constructor operates on three Vec; it allocates the necessary memory and then performs the computation. Thus only one memory allocation is performed.

An example implementation of expression templates looks like the following. A base class VecExpression represents any vector-valued expression. It is templated on the actual expression type E to be implemented, per the curiously recurring template pattern.

 1 template <typename E>
 2 class VecExpression {
 3   public:
 4     double operator[](size_t i) const { return static_cast<E const&>(*this)[i];     }
 5     size_t size()               const { return static_cast<E const&>(*this).size(); }
 6 
 7     // The following overload conversions to E, the template argument type;
 8     // e.g., for VecExpression<VecSum>, this is a conversion to VecSum.
 9           E& operator()()       { return static_cast<      E&>(*this); }
10     const E& operator()() const { return static_cast<const E&>(*this); }
11 };

The Vec class still stores the coordinates of a fully evaluated vector expression, and becomes a subclass of VecExpression.

class Vec : public VecExpression<Vec> {
    std::vector<double> elems;

  public:
    double operator[](size_t i) const { return elems[i]; }
    double &operator[](size_t i)      { return elems[i]; }
    size_t size() const               { return elems.size(); }

    Vec(size_t n) : elems(n) {}

    // A Vec can be constructed from any VecExpression, forcing its evaluation.
    template <typename E>
    Vec(VecExpression<E> const& vec) : elems(vec.size()) {
        for (size_t i = 0; i != vec.size(); ++i) {
            elems[i] = vec[i];
        }
    }
};

The sum of two vectors is represented by a new type, VecSum, that is templated on the types of the left- and right-hand sides of the sum so that it can be applied to arbitrary pairs of vector expressions. An overloaded operator+ serves as syntactic sugar for the VecSum constructor.

 1 template <typename E1, typename E2>
 2 class VecSum : public VecExpression<VecSum<E1, E2> > {
 3     E1 const& _u;
 4     E2 const& _v;
 5 
 6 public:
 7     VecSum(E1 const& u, E2 const& v) : _u(u), _v(v) {
 8         assert(u.size() == v.size());
 9     }
10 
11     double operator[](size_t i) const { return _u[i] + _v[i]; }
12     size_t size()               const { return _v.size(); }
13 };
14   
15 template <typename E1, typename E2>
16 VecSum<E1,E2> const
17 operator+(E1 const& u, E2 const& v) {
18    return VecSum<E1, E2>(u, v);
19 }

With the above definitions, the expression a + b + c is of type

VecSum<VecSum<Vec, Vec>, Vec>

so Vec x = a + b + c invokes the templated Vec constructor with this type as its E template argument. Inside this constructor, the loop body

elems[i] = v[i];

is effectively expanded (following the recursive definitions of operator+ and operator[] on this type) to

elems[i] = a.elems[i] + b.elems[i] + c.elems[i];

with no temporary vectors needed and only one pass through each memory block.

Applications

Expression templates have been found especially useful by the authors of libraries for linear algebra, i.e., for dealing with vectors and matrices of numbers. Among libraries employing expression template are POOMA,^[5] Boost uBLAS,^[6] Eigen,^[7] Blitz++^[8] and Stan Math Library.^[9]

Outside of vector math, the Spirit parser framework uses expression templates to represent formal grammars and compile these into parsers.

References