C++ expression templates explained with a C++20 example

Expression templates make a + b + c return a small object that describes the sum instead of computing it, and the assignment to a vector then runs one loop that computes a[i] + b[i] + c[i] for each element, with no temporary vectors. In C++20 a working version fits in about eighty lines: a class template for the nodes of the expression, operators constrained by a concept, and a vector type that knows how to evaluate a node. For real numeric work, use a library that has already solved the hard parts, Eigen for linear algebra or a ranges pipeline for a one-off fused loop, and write your own only for a numeric type of your own that appears in many expressions.

This guide is written for C++20, which adds the piece older tutorials work around with SFINAE or a CRTP base: a concept that says which types the operators accept. The technique needs nothing newer than C++98, and C++11’s moves changed the case for it, as the counts below show. Every sample was compiled with MSVC 19.44 (/std:c++20 /EHsc /W4 /permissive-) and Clang 19.1 (-std=c++20 -Wall -Wextra, on Microsoft’s standard library) on one Windows x64 machine, without optimization unless the text says otherwise, and the C++23 ranges sample with /std:c++23preview and -std=c++23; the counting programs were run on both, and the counts and errors shown are what those builds produced there.

What changed: the operators are kept to expression types with enable_if or a CRTP base class with a concept.

since C++20

The problem: a + b + c makes temporaries and extra loops

Start with a vector whose operator+ computes its result at once. It counts two things: the element buffers it allocates and the loops it runs over the elements. It is move-only, so every buffer is made in one of its two constructors and none slips past the count. This block, the lazy namespace below and the main under “Count the temporaries” make one program, and this block carries the includes for all of it.

#include <concepts>
#include <cstddef>
#include <cstdio>
#include <functional>
#include <initializer_list>
#include <type_traits>
#include <utility>
#include <vector>

struct Counts {
    int buffers = 0;  // element buffers allocated
    int passes = 0;   // loops over the elements
};
inline Counts counts;

namespace eager {

class Vec {
public:
    explicit Vec(std::size_t n) : data_(n) { ++counts.buffers; }
    Vec(std::initializer_list<double> xs) : data_(xs) {
        ++counts.buffers;
    }
    Vec(Vec&&) = default;             // move-only: declaring the moves
    Vec& operator=(Vec&&) = default;  // deletes the copies

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

private:
    std::vector<double> data_;
};

Vec operator+(const Vec& a, const Vec& b) {
    Vec r(a.size());
    ++counts.passes;
    for (std::size_t i = 0; i < r.size(); ++i) r[i] = a[i] + b[i];
    return r;
}

}  // namespace eager

Vec s = a + b + c; parses as (a + b) + c. The first operator+ allocates a buffer and loops over a and b. The second allocates another and loops over that temporary and c, and the temporary is destroyed at the end of the statement. Since C++17 the second result initializes s directly (guaranteed copy elision); whether the local r inside operator+ is elided too (NRVO) is up to the compiler, and here it doesn’t change the count, because moving a Vec allocates nothing. So: two buffers, two passes. Per element, the eager code reads four doubles and writes two, where one loop would read three and write one, and every further + adds one buffer and one pass. Why declaring the moves deletes the copies is in the rule of zero guide.

The idea: operators build a description, assignment runs it

With expression templates, a + b adds nothing. It returns an object that remembers its two operands and the operation, and whose operator[](i) computes element i when asked. (a + b) + c wraps that object in another one. Only when the result is assigned to a Vec does a loop run, asking the outer node for each element, which asks the inner node, which reads a[i] and b[i]. The type of the expression is its parse tree, and with the code of the next section the compiler confirms it:

lazy::Vec a{1, 2, 3}, b{10, 20, 30}, c{100, 200, 300};

using lazy::BinOp, lazy::Vec;
using Sum = BinOp<BinOp<Vec, Vec, std::plus<>>, Vec, std::plus<>>;
static_assert(std::same_as<decltype(a + b + c), Sum>);

Sum holds the inner node by value and a reference to c, and the inner node holds references to a and b: sizeof(Sum) is 24 on both compilers here, three pointers’ worth, whatever the length of the vectors. To see the type for an expression of your own, paste the code into C++ Insights or Compiler Explorer, both listed on Toolchain.

A minimal expression template in C++20

namespace lazy {

class Vec;

// Being an expression is opt-in, as std::ranges::enable_view is for
// views: a type with operator[] and size() is not enough.
template <class E>
inline constexpr bool enable_vec_expr = false;

template <>
inline constexpr bool enable_vec_expr<Vec> = true;

template <class E>
concept VecExpr = enable_vec_expr<E> &&
    requires(const E& e, std::size_t i) {
        { e[i] } -> std::convertible_to<double>;
        { e.size() } -> std::same_as<std::size_t>;
    };

// A Vec operand is held by reference, a sub-expression by value.
template <class E>
using operand_t =
    std::conditional_t<std::same_as<E, Vec>, const Vec&, E>;

template <class L, class R, class Op>
struct BinOp {
    operand_t<L> l;
    operand_t<R> r;

    double operator[](std::size_t i) const { return Op{}(l[i], r[i]); }
    std::size_t size() const { return l.size(); }
};

template <class L, class R, class Op>
inline constexpr bool enable_vec_expr<BinOp<L, R, Op>> = true;

template <VecExpr L, VecExpr R>
BinOp<L, R, std::plus<>> operator+(const L& l, const R& r) {
    return {l, r};
}

template <VecExpr L, VecExpr R>
BinOp<L, R, std::multiplies<>> operator*(const L& l, const R& r) {
    return {l, r};
}

class Vec {
public:
    explicit Vec(std::size_t n) : data_(n) { ++counts.buffers; }
    Vec(std::initializer_list<double> xs) : data_(xs) {
        ++counts.buffers;
    }
    Vec(Vec&&) = default;
    Vec& operator=(Vec&&) = default;

    template <VecExpr E>               // Vec s = a + b + c;
    Vec(const E& e) : Vec(e.size()) { assign(e); }

    template <VecExpr E>               // r = a + b + c;
    Vec& operator=(const E& e) {       // reuses r's buffer
        assign(e);
        return *this;
    }

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

private:
    template <class E>
    void assign(const E& e) {          // the one loop
        ++counts.passes;
        for (std::size_t i = 0; i < size(); ++i) data_[i] = e[i];
    }

    std::vector<double> data_;
};

}  // namespace lazy

What each part does:

• BinOp<L, R, Op> is a node: two operands and an operation, here std::plus<> or std::multiplies<>. Its operator[] computes one element and stores nothing.
• operand_t decides how a node keeps an operand. A Vec is held by reference, since copying it would copy the data. A node is held by value: it is a few references, and the inner node of a + b + c is a temporary that dies at the end of the statement, so holding it by reference would dangle the moment the expression is stored.
• operator+ and operator* only build nodes, with aggregate initialization.
• Vec gets a constructor and an assignment operator from any expression. Both end in assign, the only loop, and the assignment writes into the buffer r already has.

To keep the sketch short it leaves out size checks (a real type asserts that the operands’ sizes match), scalars (2.0 * a needs a node for a scalar) and aliasing beyond element-wise operations. a = a + b is safe here because element i reads only index i, but an expression that reads other indices, such as a shift or a matrix product, would read elements the loop has already overwritten.

Constraining the operators with a concept

The operators are templates, so without a constraint they would accept any two types and fail later. VecExpr asks for an operator[] whose result converts to double, a size() that returns std::size_t, and an opt-in through enable_vec_expr. Add a scalar to see it work:

lazy::Vec a{1, 2, 3};
lazy::Vec r = a + 2.0;   // error: no operator+ matches

Both compilers reject the line where it is written. Clang reports that the candidate was ignored because double does not satisfy VecExpr, since enable_vec_expr<double> evaluated to false; MSVC reports that the associated constraints are not satisfied and names lazy::VecExpr<double>. With the operators declared as plain template <class L, class R>, the same line builds a BinOp<Vec, double, std::plus<>> and fails inside BinOp::operator[], where r[i] is applied to a double, at the end of a chain of “in instantiation of” notes that leads back to your line. That is the effect cppreference’s page on constraints and concepts describes: a violated constraint is caught early in instantiation, and the error is easier to follow.

The opt-in is there because the syntax alone accepts too much. std::vector<double> has an operator[] and a size() too: with the enable_vec_expr test removed from the concept, v + w on two vectors, after a using namespace lazy;, compiled on both compilers into a BinOp holding copies of both vectors. The standard library makes views opt in the same way, with std::ranges::enable_view. For writing concepts, and for what they do and don’t check, see C++20 concepts: how to write and use them. Concepts need C++20; GCC 10, Clang 10 and Visual Studio 2019 16.3 are the first releases the C++20 table of Since When lists for them.

Count the temporaries and loops, before and after

This main runs the same two statements on both vector types and prints the counts:

template <class V>
void report(const char* name) {
    V a{1, 2, 3}, b{10, 20, 30}, c{100, 200, 300};
    V r(3);

    counts = {};
    V s = a + b + c;
    std::printf("%s V s = a + b + c;  buffers %d, passes %d, s[2] %g\n",
                name, counts.buffers, counts.passes, s[2]);

    counts = {};
    r = a + b + c;
    std::printf("%s r = a + b + c;    buffers %d, passes %d, r[2] %g\n",
                name, counts.buffers, counts.passes, r[2]);
}

int main() {
    report<eager::Vec>("eager");
    report<lazy::Vec>("lazy ");
}

Both builds printed the same counts, and the same sums (333 in every line):

The counts don’t depend on the length of the vectors or on the optimizer: they count what the source asks for. They say nothing about time. Whether one pass instead of two, and no temporary, makes your code faster depends on the length, on whether the data fits in cache, and on what the optimizer already does with the eager loops. Measure your own sizes and types on a release build with a benchmarking library such as Google Benchmark or nanobench, both on Toolchain.

The auto trap: expressions that outlive their operands

auto deduces the node type, not Vec:

lazy::Vec a{1, 2, 3}, b{10, 20, 30};

auto sum = a + b;               // a BinOp: references to a and b
a[0] = 100;
std::printf("%g\n", sum[0]);    // 110: computed now, from the new a[0]

Both builds printed 110. sum is a description that reads a and b each time it is indexed, so a loop that reads it twice computes it twice. Worse, a node keeps references to its Vec operands, and a temporary operand dies at the end of the full-expression:

// Returns a BinOp that refers to the temporary Vec{1, 1, 1}, which
// is destroyed when the function returns: reading it is UB.
auto shifted(const lazy::Vec& v) { return v + lazy::Vec{1, 1, 1}; }

// Returns a Vec: the sum is computed while the temporary is alive.
lazy::Vec shifted_ok(const lazy::Vec& v) {
    return v + lazy::Vec{1, 1, 1};
}

shifted compiles on both compilers without a warning, and auto bad = a + lazy::Vec{1, 2, 3}; dangles the same way at its semicolon. Naming the result type, as shifted_ok does, evaluates the expression while its operands are alive. Eigen gives its users the same advice on its Common pitfalls page: in short, don’t use auto with Eigen’s expressions unless you are sure of what you are doing, because the variable holds an abstract expression rather than a matrix; where you do use it, call .eval() on the whole expression. Clang can catch both dangling forms: with [[clang::lifetimebound]] on the operators’ parameters, Clang 19 warned on each (-Wreturn-stack-address for shifted, -Wdangling for the variable). MSVC 19.44 doesn’t recognize that attribute and says so (warning C5030), so put it behind a macro if you build with both.

What it costs: compile times, errors and debugging

• Compile time. Every distinct expression shape is a new type. a + b, a + b + c and a * b + c each instantiate a new BinOp and new specializations of Vec’s constructor, assignment and assign. A code base with many formulas instantiates many of them, and that work happens in every translation unit that uses them.
• Errors. Diagnostics spell out the whole tree: the unconstrained error above named lazy::BinOp<lazy::Vec, double, std::plus<void>>, and a real formula nests far deeper. Constraints move the error to the call, which is the main reason to write them.
• Debugging. Stepping into s[i] goes through one operator[] per node, and a debugger shows a node’s references, not values. To look at an intermediate result, evaluate it into a Vec first.
• New failure modes. Lifetimes and aliasing, from the sections above, are now the type’s problem; the eager version had neither.

Do moves make them unnecessary?

Partly. C++11 lets operator+ see that its left operand is a temporary and reuse that buffer:

namespace eager {

// When the left operand is a temporary, add into its buffer.
Vec operator+(Vec&& a, const Vec& b) {
    ++counts.passes;
    for (std::size_t i = 0; i < a.size(); ++i) a[i] += b[i];
    return a;  // C++20: an rvalue reference parameter moves implicitly
}

}  // namespace eager

Add it after the first block, and the eager rows of the table become one buffer and two passes for both statements, on both compilers: a + b allocates the one buffer, and (a + b) + c adds c into it and moves it out. Moves remove the allocations; they don’t remove the passes, because each operator finishes its loop before the next one starts. Whether the second pass costs anything you can measure depends on your sizes, which is the benchmark’s question again. return a; needs no std::move: since C++20 a return moves an rvalue reference parameter implicitly (P1825R0). Code that has to build under the C++17 rules writes return std::move(a);, because those rules treat a as an lvalue there, and Vec can’t be copied. The move semantics guide has the rules for each standard.

Where you already use them: Eigen, std::valarray and others

• Eigen. Its arithmetic returns expression objects, which is why its documentation warns about auto. For linear algebra, use it rather than growing a sketch like this one: matrix products and aliasing are problems it has already worked through, and its documentation has a page on aliasing.
• std::valarray. The standard lets its operators return proxy objects so that an expression such as v1 = a * v2 + v3 can run as one loop, and cppreference’s std::valarray notes say libstdc++ and libc++ use expression templates for it. Microsoft’s STL doesn’t: there, decltype(va + vb) is std::valarray<double>, as a static_assert confirmed on both compilers here. So auto with a valarray expression can deduce a proxy on one library and a valarray on another.
• xtensor. Its README describes “an extensible expression system enabling lazy broadcasting”, whose values are computed on access or when the expression is assigned to an xarray.
• Ranges views. A view describes a computation and runs it when iterated, the same idea for sequences. C++23’s std::views::zip_transform gives this article’s fused loop with no type of your own:

#include <ranges>
#include <vector>

std::vector<double> a{1, 2, 3}, b{10, 20, 30}, c{100, 200, 300};

auto add3 = [](double x, double y, double z) { return x + y + z; };
auto r = std::views::zip_transform(add3, a, b, c)
       | std::ranges::to<std::vector>();     // C++23: one loop

It needs C++23 mode: here it compiled with /std:c++23preview and /std:c++latest on MSVC 19.44 and with -std=c++23 on Clang 19, and under C++20 both compilers rejected zip_transform and std::ranges::to (the latter has a row in the C++23 table of Since When). With a counting operator new in a release build, both builds made one allocation, the one for r. For a fused loop you need once, this or a plain for loop is a better choice than expression templates of your own.

If you do need to build an embedded language of your own, Boost.YAP is a library for exactly that; its manual opens with An Expression Template Primer, which defines expression templates as “templates that you write to capture expressions so that they can be transformed and/or evaluated lazily”. Write your own when you own a numeric type that appears in many expressions, the temporaries or passes show up when you profile, and no library fits. More generic code, from function templates to concepts, is in Templates.