C++
Idiom
since C++98
Advanced

Expression Templates

Encode arithmetic expression trees as nested types so chained operations compile to a single loop with no intermediate heap allocations.

Expression Templatessince C++98

A metaprogramming technique that represents arithmetic sub-expressions as nested template types, deferring element-wise evaluation until assignment to eliminate heap allocations for intermediate results.

Overview

When you write result = a + b + c with a naive vector type, operator+ is called twice. The first call allocates storage, fills it, and returns a temporary t1 = a + b. The second call repeats this for result = t1 + c. For an n-element vector you pay two heap allocations, three full passes, and two destructions — where a hand-written loop would need one pass and zero allocations.

Expression templates solve this by making operator+ return a lightweight expression node — a struct holding const-references to its operands — instead of a materialised result. The type of a + b + c becomes Add<Add<Vec, Vec>, Vec>. No storage is allocated. The single element-wise loop only runs when the expression is assigned to a Vec. A good optimiser sees through the entire call chain and emits code identical to the hand-written loop.

This approach dates to the C++98/03 era. Every serious numerical library uses a variant of it: Eigen, Blaze, Armadillo, and xtensor all build their zero-overhead abstraction on this foundation.

Syntax

The canonical structure uses CRTP (C++98) to give every expression node a common base, enabling a single constrained operator= on the concrete type.

cpp
#include <vector>
#include <cstddef>
#include <initializer_list>

// C++98: CRTP expression base
template<typename Derived>
struct VecBase {
    const Derived& derived() const { return static_cast<const Derived&>(*this); }
    double      operator[](std::size_t i) const { return derived()[i]; }
    std::size_t size()                    const { return derived().size(); }
};

// Expression node: element-wise addition
template<typename L, typename R>
struct Add : VecBase<Add<L, R>> {
    const L& lhs;
    const R& rhs;
    Add(const L& l, const R& r) : lhs(l), rhs(r) {}
    double      operator[](std::size_t i) const { return lhs[i] + rhs[i]; }
    std::size_t size()                    const { return lhs.size(); }
};

// Expression node: element-wise multiplication
template<typename L, typename R>
struct Mul : VecBase<Mul<L, R>> {
    const L& lhs;
    const R& rhs;
    Mul(const L& l, const R& r) : lhs(l), rhs(r) {}
    double      operator[](std::size_t i) const { return lhs[i] * rhs[i]; }
    std::size_t size()                    const { return lhs.size(); }
};

// Expression node: scalar multiplication
template<typename E>
struct Scale : VecBase<Scale<E>> {
    double   scalar;
    const E& expr;
    Scale(double s, const E& e) : scalar(s), expr(e) {}
    double      operator[](std::size_t i) const { return scalar * expr[i]; }
    std::size_t size()                    const { return expr.size(); }
};

// Concrete storage type
class Vec : public VecBase<Vec> {
    std::vector<double> data_;
public:
    explicit Vec(std::size_t n, double fill = 0.0) : data_(n, fill) {}
    Vec(std::initializer_list<double> il) : data_(il) {}

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

    // Assign from any expression — one loop, no temporaries  (C++98)
    template<typename Expr>
    Vec& operator=(const VecBase<Expr>& rhs) {
        const Expr& e = static_cast<const Expr&>(rhs);
        for (std::size_t i = 0; i < data_.size(); ++i)
            data_[i] = e[i];
        return *this;
    }

    // Construct from expression — size deduced from the tree  (C++98)
    template<typename Expr>
    explicit Vec(const VecBase<Expr>& rhs) : data_(rhs.size()) { *this = rhs; }
};

// Operators return expression nodes — not Vec
template<typename L, typename R>
Add<L, R>   operator+(const VecBase<L>& l, const VecBase<R>& r) { return {l.derived(), r.derived()}; }

template<typename L, typename R>
Mul<L, R>   operator*(const VecBase<L>& l, const VecBase<R>& r) { return {l.derived(), r.derived()}; }

template<typename E>
Scale<E>    operator*(double s, const VecBase<E>& e) { return {s, e.derived()}; }
template<typename E>
Scale<E>    operator*(const VecBase<E>& e, double s) { return {s, e.derived()}; }

Examples

Single-pass evaluation

cpp
Vec a{1.0, 2.0, 3.0};
Vec b{4.0, 5.0, 6.0};
Vec c{7.0, 8.0, 9.0};

// Deduced type: Add<Add<Vec,Vec>, Vec>  — zero allocations
Vec result(a.size());
result = a + b + c;          // one loop: result[i] = a[i] + b[i] + c[i]

// Mixed: 2*a[i] + b[i]*c[i] — still one loop
Vec result2(2.0 * a + b * c);

At -O2 the compiler inlines every operator[] call in the tree; the emitted assembly is indistinguishable from a hand-written loop.

Eigen: expression templates in the wild

Eigen has used this pattern since its first release and extends it to matrix products, reductions, and block operations.

cpp
#include <Eigen/Dense>

Eigen::VectorXd a(1000), b(1000), c(1000);
// fills omitted

// Type: Eigen::CwiseBinaryOp<...> — no allocation until assignment
Eigen::VectorXd r = a + b + c;

// .eval() forces immediate materialisation — needed before aliasing ops
auto t = (a + b).eval();

// Avoid a defensive copy for non-aliasing matrix products  (Eigen-specific API)
Eigen::MatrixXd A(100,100), B(100,100), C(100,100);
C.noalias() += A * B;   // skips the default temporary that guards against aliasing

C++20 Simplifications

C++20 concepts (C++20) replace CRTP as the constraint mechanism. The expression node types are unchanged; only the constraint syntax improves.

cpp
// C++20: concept models the VecBase protocol without inheritance
template<typename T>
concept VecExpr = requires(const T& e, std::size_t i) {
    { e[i]     } -> std::convertible_to<double>;
    { e.size() } -> std::convertible_to<std::size_t>;
};

// Abbreviated function template — C++20
VecExpr auto operator+(const VecExpr auto& l, const VecExpr auto& r) {
    return Add{l, r};   // CTAD deduces Add<L, R>  — C++17
}

// Vec::operator= with concept constraint instead of CRTP base parameter
template<VecExpr Expr>
Vec& Vec::operator=(const Expr& rhs) {
    for (std::size_t i = 0; i < data_.size(); ++i)
        data_[i] = rhs[i];
    return *this;
}

When a type mismatch occurs, the compiler reports which VecExpr requirement the offending type fails to satisfy — far more readable than the multi-screen instantiation backtraces produced by unconstrained C++98/11 templates.

Best Practices

Expose .eval() explicitly. Any expression type should provide a method that forces materialisation:

cpp
template<typename Derived>
struct VecBase {
    Vec eval() const;   // defined after Vec is complete
};

Callers use (a + b).eval() to state intent unambiguously — helpful at API boundaries and before storing an expression in a named variable.

Keep the element type generic. Replace double with a template parameter T on every node and Vec. The expression tree pattern is type-agnostic; locking it to double rules out float, int, and complex numbers for no reason.

Assert the CRTP invariant in C++17. Before C++20, add a static_assert in operator= to produce a readable diagnostic instead of a wall of substitution failures:

cpp
template<typename Expr>
Vec& operator=(const VecBase<Expr>& rhs) {
    static_assert(std::is_base_of_v<VecBase<Expr>, Expr>,  // C++17
                  "rhs must derive from VecBase<Expr>");
    ...
}

Common Pitfalls

Dangling references

Every expression node stores const& to its operands. If an operand is a temporary returned by value, it is destroyed before the assignment loop runs:

cpp
// UB: temporaries from factory() are destroyed before Vec's constructor loop
Vec bad(makeVec(1) + makeVec(2));

// Safe: materialise each operand first
Vec tmp1 = makeVec(1);
Vec tmp2 = makeVec(2);
Vec good(tmp1 + tmp2);

This is the most common production bug with expression templates. Libraries that need to support temporaries store operands by value when the operand is an rvalue — a technique Eigen applies via internal::ref_selector.

Aliasing in place-wise operations

For element-wise nodes each output index reads only its own input index, so v = v + v is safe: v[i] is read before it is overwritten. However, any operation where output index i depends on multiple input indices — matrix-matrix multiply, reductions, convolutions — requires a defensive copy unless the library can prove the expression is alias-free. This is why Eigen defaults to a temporary for C = A * B and requires .noalias() to opt out.

Template instantiation depth

Deep expression trees like a + b + c + d + ... + z produce deeply nested types that can exceed default template instantiation depth limits on some compilers (-ftemplate-depth on GCC/Clang). Real libraries keep trees shallow with N-ary expression nodes or cap expression depth via traits.

C++20 Ranges as a Modern Alternative

C++20 ranges offer lazy pipelines that look superficially like expression templates, but the two mechanisms are complementary rather than interchangeable. Ranges operate on iterables and produce another range; they carry no static knowledge of vector sizes and produce no expression tree — just a chain of adaptor objects.

cpp
// C++20: lazy element-wise addition via ranges
#include <ranges>
#include <vector>
#include <algorithm>

std::vector<double> a{1, 2, 3}, b{4, 5, 6};

// zip pairs (a[0],b[0]), (a[1],b[1]), ... then transforms lazily — no allocation
auto sum_view = std::views::zip(a, b)
    | std::views::transform([](auto [x, y]) { return x + y; });

// Materialise into a vector when needed
std::vector<double> result(sum_view.begin(), sum_view.end());

For general algorithms on arbitrary containers, ranges eliminate the boilerplate of writing custom expression nodes. The pipeline above compiles to a single pass with no heap allocations when the output is pre-sized. Where ranges fall short is in domain-specific linear algebra: they have no concept of vector dimension, cannot fuse a matrix-vector product with a subsequent addition at the type level, and provide no mechanism for size mismatch detection at compile time. Eigen-style expression templates remain the right tool when the shape of the computation must be known and verified statically.

Expression TemplatesC++20 Ranges
Lazy evaluationYesYes
Static size knowledgeYes (via traits)No
Works on arbitrary iterablesNoYes
Custom reduction fusionYesLimited
Standard libraryNoYes (<ranges>)

Constexpr Expression Templates

The expression template pattern is fully compatible with consteval/constexpr evaluation. All that is required is that the storage type uses constexpr-capable allocation — which std::vector satisfies from C++20 onward, making heap-backed Vec usable in constant expressions without any structural changes.

cpp
// C++20: constexpr expression template nodes
#include <vector>
#include <cstddef>

template<typename Derived>
struct VecBase {
    constexpr const Derived& derived() const { return static_cast<const Derived&>(*this); }
    constexpr double         operator[](std::size_t i) const { return derived()[i]; }
    constexpr std::size_t    size()                    const { return derived().size(); }
};

template<typename L, typename R>
struct Add : VecBase<Add<L, R>> {
    const L& lhs; const R& rhs;
    constexpr Add(const L& l, const R& r) : lhs(l), rhs(r) {}
    constexpr double      operator[](std::size_t i) const { return lhs[i] + rhs[i]; }
    constexpr std::size_t size()                    const { return lhs.size(); }
};

template<typename E>
struct Scale : VecBase<Scale<E>> {
    double scalar; const E& expr;
    constexpr Scale(double s, const E& e) : scalar(s), expr(e) {}
    constexpr double      operator[](std::size_t i) const { return scalar * expr[i]; }
    constexpr std::size_t size()                    const { return expr.size(); }
};

// C++20: std::vector is constexpr-capable
class Vec : public VecBase<Vec> {
    std::vector<double> data_;   // heap allocation allowed in constexpr — C++20
public:
    constexpr explicit Vec(std::size_t n, double fill = 0.0) : data_(n, fill) {}
    constexpr Vec(std::initializer_list<double> il) : data_(il) {}
    constexpr double      operator[](std::size_t i) const { return data_[i]; }
    constexpr double&     operator[](std::size_t i)       { return data_[i]; }
    constexpr std::size_t size()                    const { return data_.size(); }

    template<typename Expr>
    constexpr Vec& operator=(const VecBase<Expr>& rhs) {
        const Expr& e = static_cast<const Expr&>(rhs);
        for (std::size_t i = 0; i < data_.size(); ++i) data_[i] = e[i];
        return *this;
    }
    template<typename Expr>
    constexpr explicit Vec(const VecBase<Expr>& rhs) : data_(rhs.size()) { *this = rhs; }
};

template<typename L, typename R>
constexpr Add<L,R> operator+(const VecBase<L>& l, const VecBase<R>& r) {
    return {l.derived(), r.derived()};
}
template<typename E>
constexpr Scale<E> operator*(double s, const VecBase<E>& e) {
    return {s, e.derived()};
}

// Compile-time verification of the entire expression tree
consteval bool check_sum() {
    Vec a{1.0, 2.0, 3.0};
    Vec b{4.0, 5.0, 6.0};
    Vec c{7.0, 8.0, 9.0};
    Vec result(a + b + c);                   // C++20
    return result[0] == 12.0 && result[2] == 18.0;
}
static_assert(check_sum(), "constexpr ET sum incorrect");

Note the key prerequisite: std::vector became constexpr-capable in C++20 (P0784R7). In C++17 and earlier you must use a fixed-size array (std::array) as the backing store, which limits Vec to compile-time-known sizes — still useful for SIMD-width vectors or small matrix libraries, but not general-purpose.

Expression Templates for String Building

Expression templates are not confined to linear algebra. Any domain where an eagerly evaluated operator produces wasteful intermediates is a candidate. String concatenation is a classic example: naïve chaining with std::string::operator+ creates N-1 heap-allocated temporaries for N operands.

cpp
// C++17: string expression template
#include <string>
#include <string_view>
#include <cstddef>

// Expression node: concatenation of two sub-expressions
template<typename L, typename R>
struct StrCat {
    const L& lhs;
    const R& rhs;

    StrCat(const L& l, const R& r) : lhs(l), rhs(r) {}

    // Two-pass strategy: measure total length, then allocate once
    std::size_t total_size() const {
        return lhs.total_size() + rhs.total_size();
    }
    void append_to(std::string& out) const {
        lhs.append_to(out);
        rhs.append_to(out);
    }

    // Implicit conversion materialises in a single allocation  (C++17)
    operator std::string() const {
        std::string result;
        result.reserve(total_size());
        append_to(result);
        return result;
    }
};

// Leaf wrapper for std::string_view
struct StrLeaf {
    std::string_view sv;
    explicit StrLeaf(std::string_view s) : sv(s) {}
    std::size_t total_size()            const { return sv.size(); }
    void        append_to(std::string& out) const { out.append(sv); }
    operator std::string() const { return std::string(sv); }
};

// Operator builds the expression tree
template<typename L, typename R>
StrCat<L, R> operator+(const L& l, const R& r) { return {l, r}; }

// Usage
void example() {
    StrLeaf first{"Hello"}, sep{", "}, last{"World!"};

    // Type: StrCat<StrCat<StrLeaf,StrLeaf>,StrLeaf>
    // one reserve() + three append() calls — zero intermediate strings
    std::string result = first + sep + last;   // C++17

    // Contrast: naïve approach creates two temporaries
    std::string naive = std::string("Hello") + ", " + "World!";
}

The two-pass strategy — walk the tree to measure total length, then walk it again to fill a pre-reserved buffer — is exactly what absl::StrCat and fmt::format implement internally. The ET encoding makes the two passes explicit in the type system. For very long chains (hundreds of operands) the single allocation can reduce total memory traffic by an order of magnitude compared to quadratic naive concatenation.

Best Practice: prefer std::string::reserve + append for small concatenations. The StrCat ET pattern pays off only when the chain length is not known at the call site or when the operands are themselves computed lazily. For three or four literal strings, absl::StrCat or C++20's std::format is simpler and equivalent in performance.

Aliasing: The Store-Before-Read Problem

Aliasing occurs when the output of an expression overlaps with one of its inputs. For element-wise operations this is harmless — each result[i] depends only on a[i] and b[i], so overwriting result[i] cannot affect a later result[j]. Matrix multiplication breaks this invariant: computing C[i][j] requires the entire i-th row of A and j-th column of B, so if C aliases either operand the stored value is corrupted mid-loop.

cpp
// C++98 / C++17: illustrating the aliasing hazard
Eigen::MatrixXd A(3,3), B(3,3);
A << 1,2,3, 4,5,6, 7,8,9;
B = A;

// UNSAFE without a temporary: A[0][j] is overwritten before A[1][j] is read
// Eigen inserts a defensive copy by default for matrix products
A = A * B;           // Eigen detects aliasing, uses a temporary automatically

// .noalias() asserts no aliasing — Eigen skips the defensive copy
// ONLY use when C is provably disjoint from A and B
Eigen::MatrixXd C(3,3);
C.noalias() = A * B;   // safe: C does not alias A or B

Eigen detects aliasing via internal::evaluator_traits::is_aliasing, which compares the data pointers of the destination and each leaf in the expression tree at runtime. The .noalias() call attaches a NoAlias<MatrixXd> wrapper to the assignment, which is an expression-tree annotation that disables the defensive copy path entirely — it does not change evaluation semantics, only the copy guard.

cpp
// Conceptual sketch: .noalias() as an expression annotation  (C++11)
template<typename Derived>
struct MatBase {
    NoAlias<Derived> noalias() { return NoAlias<Derived>{derived()}; }
};

template<typename Dest, typename Src>
void assign(Dest& dst, const NoAlias<Src>& src) {
    // skip aliasing check — evaluate directly into dst
    eval_to(dst, src.expr);
}

For the small toy Vec type above, a simple compile-time check can catch the most obvious aliasing:

cpp
// C++17: static_assert that detects same-object self-assignment in a numeric ET
template<typename Expr>
Vec& operator=(const VecBase<Expr>& rhs) {
    // If Expr IS Vec, check that &rhs != this  (runtime, but triggered at assignment)
    if constexpr (std::is_same_v<Expr, Vec>) {
        if (static_cast<const Vec*>(&rhs) == this)
            return *this;  // self-assignment is a no-op for element-wise ops
    }
    for (std::size_t i = 0; i < data_.size(); ++i)
        data_[i] = static_cast<const Expr&>(rhs)[i];
    return *this;
}

Common Pitfall: never use .noalias() speculatively. If the assertion is wrong and aliasing actually exists, the result is a silent data corruption bug that is extremely difficult to reproduce because it depends on matrix size, loop order, and SIMD vectorisation.

Type Explosion and Compile-Time Cost

Every additional term in an expression chain adds one level of template instantiation. The type of a + b + c + d + e + f is Add<Add<Add<Add<Add<Vec,Vec>,Vec>,Vec>,Vec>,Vec> — six levels deep. You can inspect this with decltype:

cpp
// C++11: visualise the nested type
Vec a{1}, b{2}, c{3}, d{4}, e{5}, f{6};
using T = decltype(a + b + c + d + e + f);
// T = Add<Add<Add<Add<Add<Vec,Vec>,Vec>,Vec>,Vec>,Vec>
static_assert(!std::is_same_v<T, Vec>, "not yet evaluated");

Each level is a distinct template instantiation. For a chain of length N the compiler must instantiate O(N) types, each with its own operator[] and size() specialisations. In a large numerical codebase with many expression variables this can dominate compilation time.

Three practical mitigations:

(a) Break chains with .eval() every N terms. Materialising after every four or five terms caps instantiation depth and gives the compiler smaller inlining units. The cost is one extra allocation per break point, which is usually negligible.

cpp
// C++98: manual chain break
Vec t1(a + b + c + d);   // depth 3
Vec result(t1 + e + f);  // depth 2

(b) Use an N-ary node that flattens the tree. Rather than left-associative binary nodes, an AddN<Pack...> variadic template holds all operands in a std::tuple and iterates over them in a single operator[]:

cpp
// C++17: N-ary addition node — depth stays at 1 regardless of operand count
template<typename... Ops>
struct AddN : VecBase<AddN<Ops...>> {
    std::tuple<const Ops*...> ops;
    explicit AddN(const Ops&... o) : ops{&o...} {}

    double operator[](std::size_t i) const {
        return std::apply([i](auto*... p){ return ((*p)[i] + ...); }, ops);  // C++17 fold
    }
    std::size_t size() const { return std::get<0>(ops)->size(); }
};

(c) Raise the compiler's instantiation depth limit. As a stopgap while refactoring, pass -ftemplate-depth=1000 (GCC/Clang) or /D_VARIADIC_MAX=15 (MSVC) to allow deeper recursion. This does not fix the root cause and increases compiler memory usage, but it unblocks compilation during development.

cpp
# GCC / Clang: raise template instantiation depth limit
g++ -ftemplate-depth=1000 -O2 my_numerical_code.cpp

Best Practice: profile compile times with -ftime-report (GCC) or -ftime-trace (Clang) before committing to a deep chain design. Type explosion is invisible at runtime but can make incremental rebuilds painfully slow in production codebases.

The auto Storage Pitfall

Storing an expression node in an auto variable is the most common user-facing bug with expression templates. The type deduced by auto is the expression node type — not Vec — so the stored object holds references to its operands. If those operands go out of scope or are modified before the expression is evaluated, the result is undefined behaviour.

cpp
// C++11: subtle bug — expr holds references into a and b
Vec a{1.0, 2.0, 3.0};
Vec b{4.0, 5.0, 6.0};

auto expr = a + b;   // type: Add<Vec,Vec> — stores const& to a and b

a[0] = 99.0;         // modifies a — expr[0] will now return 99.0 + 4.0 = 103.0
double v = expr[0];  // reads through the dangling/stale reference — surprising!

// Worse: if a or b is a local that goes out of scope:
auto make_bad = []() {
    Vec x{1, 2, 3};
    return x + x;   // returns Add<Vec,Vec> holding refs to a local — UB after return
};
auto bad_expr = make_bad();
double undef = bad_expr[0];   // UB: x has been destroyed

The fix is straightforward: either force evaluation immediately with .eval(), or spell out the concrete type.

cpp
// Fix 1: materialise immediately — type is Vec, no dangling reference  (C++98)
Vec expr_val = a + b;   // triggers operator=(VecBase<Expr>), stores a copy

// Fix 2: use .eval() at the point of capture  (C++11)
auto expr_safe = (a + b).eval();   // type: Vec — owns its data

// Fix 3: use explicit type annotation — communicates intent clearly  (C++11)
Vec expr_named(a + b);

This pitfall is especially subtle because expr[i] compiles without warning and produces correct results during the initial testing phase when a and b are still in scope and unmodified. The bug only manifests after refactoring moves the construction into a factory function or stores the expression across a scope boundary.

Common Pitfall: treat auto expr = <expression template>; as a code smell in any codebase that uses expression templates. Linters and static analysis tools (clang-tidy's cppcoreguidelines-avoid-non-const-global-variables, custom checks) can flag this pattern. A comment convention like // ET: must not outlive operands is better than nothing if you must store expression nodes.

See Also