Expression Templates in C++: Lazy Evaluation That Removes Temporaries in Vector and Matrix Math

By 0xmote · Apr 25, 2026 · 23 min read · Updated Sep 19, 2026 Advanced Tutorial

Key takeaways

How expression templates defer evaluation so that a + b + c runs as one loop with no temporaries, implemented for vectors and matrices, with common errors and a performance comparison.

What Are Expression Templates and Why Do We Need Them?

Problem Scenario: Temporary Objects in Vector Operations

The Problem: In mathematical libraries, a vector operation like result = a + b + c + d creates three temporary objects. A new vector is allocated and copied for each + operation.

class Vector {
public:
    Vector(size_t n) : data(n) {}
    
    Vector operator+(const Vector& other) const {
        Vector result(data.size());
        for (size_t i = 0; i < data.size(); ++i) {
            result.data[i] = data[i] + other.data[i];
        }
        return result;  // Temporary object
    }
    
private:
    std::vector<double> data;
};
// result = a + b + c + d;
// 1. temp1 = a + b      (Temporary object 1)
// 2. temp2 = temp1 + c  (Temporary object 2)
// 3. result = temp2 + d (Temporary object 3)

Issues:

Three memory allocations
Three loops (one for each +)
Reduced cache efficiency Solution: Expression Templates enable lazy evaluation. Instead of calculating a + b + c + d immediately, the operation is stored as an expression tree and computed all at once during assignment.

// Expression Template
Vector result = a + b + c + d;
// 1. expr = Add(Add(Add(a, b), c), d)  (Expression tree, no computation yet)
// 2. result = expr                      (Computed all at once during assignment)

Advantages:

One memory allocation (only for result)
One loop (computed in a single pass)
Improved cache efficiency

flowchart TD
    subgraph normal[Standard Operations]
        n1["a + b → temp1 (allocation)"]
        n2["temp1 + c → temp2 (allocation)"]
        n3["temp2 + d → result (allocation)"]
    end
    subgraph expr[Expression Template]
        e1["a + b + c + d → Expression Tree"]
        e2["result = Expression (1 allocation)"]
        e3["Computed in one loop"]
    end
    n1 --> n2 --> n3
    e1 --> e2 --> e3

The technique traces back to Todd Veldhuizen’s work on the Blitz++ array library in the mid-1990s, and it is the core mechanism behind why Eigen, Armadillo, and Blaze can offer an API that reads like plain mathematical notation (C = A * B + D) while generating code that is competitive with hand-written, loop-fused Fortran or C. The trick is entirely a compile-time one: no runtime interpreter walks the expression tree, no virtual dispatch happens, and no extra data structure survives past compilation. The “tree” that a + b + c + d builds is a nested C++ type — something like VecAdd<VecAdd<Vector, Vector>, Vector> — and the compiler inlines the whole chain into a single loop body as if you had hand-written the fused arithmetic yourself. This is why expression templates are considered a zero-overhead abstraction: the abstraction (writing natural-looking arithmetic) costs nothing extra at runtime compared to writing the fused loop by hand, provided the compiler can see through the template layers (which, with inlining enabled, it reliably does for this pattern). It is worth being upfront about the trade-off before diving into the implementation: expression templates buy you this runtime efficiency by spending compile-time complexity. Every distinct combination of operators and operand types instantiates a new template specialization, which means longer build times, denser and harder-to-read compiler error messages (a single misplaced type in a five-term expression can produce pages of template-instantiation backtrace), and code that is genuinely more difficult to step through in a debugger. This guide walks through the full mechanism so you understand what libraries like Eigen are doing under the hood — not necessarily so you reimplement it from scratch in application code, since the standard advice (covered later in Production Patterns) is to reach for an existing, battle-tested library whenever one exists for your domain.

Basic Structure

Minimal Expression Template

#include <iostream>
#include <vector>
// Base class for expressions
template<typename E>
class VecExpr {
public:
    double operator[](size_t i) const {
        return static_cast<const E&>(*this)[i];
    }
    
    size_t size() const {
        return static_cast<const E&>(*this).size();
    }
};
// Addition expression
template<typename LHS, typename RHS>
class VecAdd : public VecExpr<VecAdd<LHS, RHS>> {
public:
    VecAdd(const LHS& l, const RHS& r) : lhs(l), rhs(r) {}
    
    double operator[](size_t i) const {
        return lhs[i] + rhs[i];
    }
    
    size_t size() const { return lhs.size(); }
    
private:
    const LHS& lhs;
    const RHS& rhs;
};
// Vector class
class Vector : public VecExpr<Vector> {
public:
    Vector(size_t n) : data(n) {}
    
    double& operator[](size_t i) { return data[i]; }
    double operator[](size_t i) const { return data[i]; }
    size_t size() const { return data.size(); }
    
    // Assigning an Expression Template
    template<typename Expr>
    Vector& operator=(const VecExpr<Expr>& expr) {
        const Expr& e = static_cast<const Expr&>(expr);
        for (size_t i = 0; i < size(); ++i) {
            data[i] = e[i];  // Lazy evaluation
        }
        return *this;
    }
    
private:
    std::vector<double> data;
};
// Operator
template<typename LHS, typename RHS>
VecAdd<LHS, RHS> operator+(const VecExpr<LHS>& lhs, const VecExpr<RHS>& rhs) {
    return VecAdd<LHS, RHS>(
        static_cast<const LHS&>(lhs),
        static_cast<const RHS&>(rhs)
    );
}
int main() {
    Vector a(3), b(3), c(3);
    a[0] = 1; a[1] = 2; a[2] = 3;
    b[0] = 4; b[1] = 5; b[2] = 6;
    c[0] = 7; c[1] = 8; c[2] = 9;
    
    Vector result(3);
    result = a + b + c;  // Expression tree, computed during assignment
    
    for (size_t i = 0; i < result.size(); ++i) {
        std::cout << result[i] << ' ';
    }
    std::cout << '\n';  // 12 15 18
}

Key Point: a + b + c returns an expression object of type VecAdd<VecAdd<Vector, Vector>, Vector>, which is computed all at once during result = ....

Why CRTP Instead of a Virtual Base Class

The VecExpr<E> base class above uses the Curiously Recurring Template Pattern (CRTP) rather than the more familiar approach of a virtual base class with operator[] declared virtual. This choice is deliberate and central to why expression templates are fast. A virtual function call goes through a vtable lookup at runtime — cheap in absolute terms, but it also blocks the compiler from inlining across the call boundary, because the compiler cannot know at compile time which override will actually run. For a tight numeric loop executed millions of times, losing inlining means losing auto-vectorization too, since the compiler can no longer see that the loop body is a simple, branch-free arithmetic expression. CRTP sidesteps this entirely. static_cast<const E&>(*this) resolves to a concrete, fully known type at compile time — there is no dynamic dispatch, no vtable, and nothing standing between the compiler and the actual arithmetic. The cost of this design is that VecExpr cannot be used polymorphically in the traditional sense (you cannot hold a std::vector<VecExpr*> of mixed expression types), but that limitation is fine here: expression templates are meant to be consumed and discarded within a single expression, never stored or passed around as an abstract interface. This is the general shape of the trade-off CRTP always makes — you give up runtime polymorphism to get compile-time (static) polymorphism, and it pays off exactly in hot loops like this one.

Why `operator[]` Returns `double`, Not `const double&`

Notice that VecAdd::operator[] returns lhs[i] + rhs[i] by value, not by reference. This is not a stylistic choice — it is required by the semantics of a lazy expression. The result of lhs[i] + rhs[i] is a temporary double that exists only for the duration of the expression; there is no persistent storage anywhere in the tree to return a reference to, because no computation has actually happened until an element is asked for. Every operator[] call on an expression node recomputes its result from its children on demand — this is exactly what “lazy” means here. It also explains why expression templates are ill-suited to random-access patterns that read the same element repeatedly: each read re-walks and re-evaluates the whole subtree, so if your access pattern revisits the same indices many times, you can end up doing more arithmetic than the naive eager version, not less. The technique is a net win specifically for the common “compute once per assignment, single forward pass” pattern shown in Vector::operator=.

Implementing Vector Operations

Adding Multiplication and Subtraction

// Subtraction expression
template<typename LHS, typename RHS>
class VecSub : public VecExpr<VecSub<LHS, RHS>> {
public:
    VecSub(const LHS& l, const RHS& r) : lhs(l), rhs(r) {}
    
    double operator[](size_t i) const {
        return lhs[i] - rhs[i];
    }
    
    size_t size() const { return lhs.size(); }
    
private:
    const LHS& lhs;
    const RHS& rhs;
};
// Scalar multiplication expression
template<typename E>
class VecScale : public VecExpr<VecScale<E>> {
public:
    VecScale(double s, const E& e) : scalar(s), expr(e) {}
    
    double operator[](size_t i) const {
        return scalar * expr[i];
    }
    
    size_t size() const { return expr.size(); }
    
private:
    double scalar;
    const E& expr;
};
// Operators
template<typename LHS, typename RHS>
VecSub<LHS, RHS> operator-(const VecExpr<LHS>& lhs, const VecExpr<RHS>& rhs) {
    return VecSub<LHS, RHS>(
        static_cast<const LHS&>(lhs),
        static_cast<const RHS&>(rhs)
    );
}
template<typename E>
VecScale<E> operator*(double scalar, const VecExpr<E>& expr) {
    return VecScale<E>(scalar, static_cast<const E&>(expr));
}
int main() {
    Vector a(3), b(3), c(3);
    a[0] = 1; a[1] = 2; a[2] = 3;
    b[0] = 4; b[1] = 5; b[2] = 6;
    c[0] = 7; c[1] = 8; c[2] = 9;
    
    Vector result(3);
    result = 2.0 * a + b - c;  // Expression tree
    
    for (size_t i = 0; i < result.size(); ++i) {
        std::cout << result[i] << ' ';
    }
    std::cout << '\n';  // -1 -1 -3
}

The Type Explosion Problem

Look closely at what 2.0 * a + b - c actually resolves to as a C++ type: VecSub<VecAdd<VecScale<Vector>, Vector>, Vector>. Every additional operator nests another template layer, and every distinct combination of operand types and operators is a distinct type that the compiler has to instantiate, parse, and (with optimizations on) inline separately. This is the direct cost side of the trade-off mentioned earlier — a numerical expression with a dozen chained operations, mixing vectors, matrices, and scalars, can produce a type name hundreds of characters long, and if something in that expression fails to compile (a size mismatch handled at the wrong layer, an ambiguous overload), the resulting error message will print that entire nested type, often several times, once per template instantiation involved. This is the single most common complaint developers have about expression-template-heavy libraries like Eigen: the usage code is beautifully simple, but diagnosing a usage mistake can mean scrolling through a wall of template noise to find the one relevant line. In practice, nobody writes out these nested types by hand — you always assign the result of an expression directly into a concrete Vector or Matrix (as result = ... does above), letting operator=’s template parameter Expr absorb whatever type the compiler generated. The one place this bites people is when someone tries to cache an intermediate expression in a local variable with auto, which is covered in Common Errors below — it looks convenient but is one of the most reliable ways to introduce a dangling-reference bug in this style of code.

The Reference-Member Lifetime Constraint

Every expression node — VecAdd, VecSub, VecScale — stores its operands as const LHS& / const E& references, not by value. This is essential for performance: copying a Vector into every node of the tree would defeat the entire point of avoiding allocations. But it comes with a hard constraint that is easy to violate by accident — every operand referenced anywhere in the expression tree must stay alive for the entire lifetime of the expression, which in practice means it must outlive the single statement that builds and immediately consumes the tree. As long as expressions are built and evaluated in one line (result = 2.0 * a + b - c;), the named vectors a, b, and c are guaranteed to be alive for that whole statement, so this is safe. The danger only appears once an expression object escapes that single statement — through a return value or a stored auto variable — which is exactly what Error 1 below demonstrates.

Matrix Operations

Extending the Pattern to Matrices

The same CRTP structure used for Vector applies directly to Matrix — only the indexing changes from a single index to a (row, col) pair.

template<typename E>
class MatExpr {
public:
    double operator()(size_t r, size_t c) const {
        return static_cast<const E&>(*this)(r, c);
    }
    size_t rows() const { return static_cast<const E&>(*this).rows(); }
    size_t cols() const { return static_cast<const E&>(*this).cols(); }
};
template<typename LHS, typename RHS>
class MatAdd : public MatExpr<MatAdd<LHS, RHS>> {
public:
    MatAdd(const LHS& l, const RHS& r) : lhs(l), rhs(r) {}
    double operator()(size_t r, size_t c) const { return lhs(r, c) + rhs(r, c); }
    size_t rows() const { return lhs.rows(); }
    size_t cols() const { return lhs.cols(); }
private:
    const LHS& lhs;
    const RHS& rhs;
};
class Matrix : public MatExpr<Matrix> {
public:
    Matrix(size_t r, size_t c) : rows_(r), cols_(c), data(r * c) {}
    double& operator()(size_t r, size_t c) { return data[r * cols_ + c]; }
    double operator()(size_t r, size_t c) const { return data[r * cols_ + c]; }
    size_t rows() const { return rows_; }
    size_t cols() const { return cols_; }
    template<typename Expr>
    Matrix& operator=(const MatExpr<Expr>& expr) {
        const Expr& e = static_cast<const Expr&>(expr);
        for (size_t r = 0; r < rows_; ++r)
            for (size_t c = 0; c < cols_; ++c)
                (*this)(r, c) = e(r, c);
        return *this;
    }
private:
    size_t rows_, cols_;
    std::vector<double> data;
};
template<typename LHS, typename RHS>
MatAdd<LHS, RHS> operator+(const MatExpr<LHS>& lhs, const MatExpr<RHS>& rhs) {
    return MatAdd<LHS, RHS>(static_cast<const LHS&>(lhs), static_cast<const RHS&>(rhs));
}

Why Matrix Multiplication Doesn’t Fit the Same Template

Element-wise operators (+, -, scalar *) map one output element to a fixed, small set of input elements — cheap to evaluate lazily per access. Matrix multiplication requires summing over an entire row/column per output element (O(n) work per element instead of O(1)), so naively lazy-evaluating A * B * C inside operator() would redo that summation on every single access — actually making it slower than eager evaluation. Production expression-template libraries (Eigen, Blaze) special-case matrix multiplication to evaluate eagerly into a temporary rather than folding it into the lazy expression tree.

Common Errors and Solutions

Error 1: Dangling References to Temporaries

VecAdd<Vector, Vector> makeExpr() {
    Vector a(3), b(3);
    return a + b;  // ❌ VecAdd stores const Vector& to locals that are about to be destroyed
}
Vector result = makeExpr();  // undefined behavior: reads freed memory

Fix: never return an expression-template object from a function — only ever assign the expression directly into a concrete Vector/Matrix in the same statement/scope where the operands are alive. This same bug has a quieter, more common cousin that has nothing to do with function returns:

Vector a(3), b(3);
auto expr = a + b;   // ❌ looks harmless — expr is just a local variable
// ... some other code runs here, maybe reassigns a or b, or a/b go out of scope ...
Vector result = expr;  // reads whatever a and b's references now point to (or freed memory)

auto deduces expr to the expression-template type itself (VecAdd<Vector, Vector>), not to a Vector. Because that type holds references, not copies, expr is only valid for as long as a and b are unchanged and alive — and unlike the direct result = a + b; idiom, nothing about this code visually signals that constraint. The rule of thumb that avoids both variants of this bug: never let an expression-template value outlive the statement that produced it. If you need to store or reuse a computed result, assign it into a concrete Vector/Matrix immediately, forcing evaluation, rather than into auto.

Error 2: Forgetting the CRTP `static_cast`

template<typename E>
class VecExpr {
public:
    double operator[](size_t i) const {
        return (*this)[i];  // ❌ infinite recursion — calls VecExpr::operator[] again,
                             //     not the derived class's
    }
};

Fix: always static_cast<const E&>(*this) to reach the derived implementation, as shown in the base template.

Error 3: Mismatched Sizes with No Runtime Check

Vector a(3), b(5);
Vector result(3);
result = a + b;  // ❌ no bounds check — reads past the end of b's data silently

Fix: add a size assertion in the expression’s constructor or operator[] for debug builds:

VecAdd(const LHS& l, const RHS& r) : lhs(l), rhs(r) {
    assert(l.size() == r.size() && "Vector size mismatch");
}

Error 4: Overusing Expression Templates for Simple Cases

For a one-off a + b with no chaining, an expression template adds template-instantiation compile-time cost for no runtime benefit — the technique pays off specifically when multiple operations chain (2.0 * a + b - c) and the intermediate-allocation savings compound.

Production Patterns

SFINAE/Concepts Guard Against Mixed Expression Types

Constrain the operators so they only combine VecExpr-derived types, preventing accidental instantiation with unrelated types that happen to have a matching operator[].

template<typename LHS, typename RHS>
    requires std::derived_from<LHS, VecExpr<LHS>> && std::derived_from<RHS, VecExpr<RHS>>
VecAdd<LHS, RHS> operator+(const LHS& lhs, const RHS& rhs) {
    return VecAdd<LHS, RHS>(lhs, rhs);
}

Reference Real Libraries Instead of Reinventing This

In production numerical code, reach for Eigen or Blaze rather than hand-rolling expression templates — they already handle the dangling-reference pitfalls, SIMD vectorization, and matrix-multiplication special-casing this guide covers. Understanding the technique (this guide’s goal) is what lets you read and debug their internals, or apply the same pattern to a different domain (e.g. a query-builder DSL or lazy string-formatting library) where no existing library fits.

Combining with SIMD

Once an expression tree is fully built and about to be evaluated, the innermost loop (for i in 0..size) is exactly the kind of tight, branch-free loop that benefits from auto-vectorization or explicit SIMD intrinsics — expression templates and SIMD are complementary, not competing, techniques.

The Aliasing Trap: When `result` Also Appears on the Right-Hand Side

There is one more production-grade pitfall that the minimal implementation in this guide does not defend against: aliasing, where the assignment target also appears somewhere inside the expression being assigned to it. Consider result = result + a; — with the naive element-wise loop in Vector::operator=, each iteration does data[i] = e[i], and e[i] for this expression reads result[i] + a[i]. Because the loop writes result[i] and reads result[i] in the same iteration (before that particular index is overwritten), this specific case happens to be safe. But shift to something like reversing a vector into itself, or certain matrix expressions where computing output element (r, c) requires reading input elements at different indices than (r, c) — a matrix transpose-in-place is the classic example — and the lazy, read-as-you-go evaluation model silently corrupts data, because earlier loop iterations overwrite values that later iterations still need to read. This is precisely why real libraries like Eigen ship both a default (conservative) assignment path that detects potential aliasing and evaluates into a temporary first, and an explicit .noalias() escape hatch (result.noalias() = a + b;) that lets the caller assert “I know these don’t alias, skip the safety check and evaluate directly.” Aliasing detection is genuinely hard to get right in general — Eigen’s own documentation has historically flagged specific expression shapes it cannot safely detect — which is another concrete reason to prefer a mature library over a hand-rolled one for anything beyond a learning exercise: the aliasing logic alone represents years of accumulated edge-case handling that is easy to underestimate when first implementing this pattern.

Complete Example: Mathematical Library

#include <iostream>
#include <vector>
#include <cassert>
template<typename E>
class VecExpr {
public:
    double operator[](size_t i) const { return static_cast<const E&>(*this)[i]; }
    size_t size() const { return static_cast<const E&>(*this).size(); }
};
template<typename LHS, typename RHS>
class VecAdd : public VecExpr<VecAdd<LHS, RHS>> {
public:
    VecAdd(const LHS& l, const RHS& r) : lhs(l), rhs(r) { assert(l.size() == r.size()); }
    double operator[](size_t i) const { return lhs[i] + rhs[i]; }
    size_t size() const { return lhs.size(); }
private:
    const LHS& lhs;
    const RHS& rhs;
};
template<typename E>
class VecScale : public VecExpr<VecScale<E>> {
public:
    VecScale(double s, const E& e) : scalar(s), expr(e) {}
    double operator[](size_t i) const { return scalar * expr[i]; }
    size_t size() const { return expr.size(); }
private:
    double scalar;
    const E& expr;
};
class Vector : public VecExpr<Vector> {
public:
    explicit Vector(size_t n) : data(n) {}
    double& operator[](size_t i) { return data[i]; }
    double operator[](size_t i) const { return data[i]; }
    size_t size() const { return data.size(); }
    template<typename Expr>
    Vector& operator=(const VecExpr<Expr>& expr) {
        const Expr& e = static_cast<const Expr&>(expr);
        for (size_t i = 0; i < size(); ++i) data[i] = e[i];
        return *this;
    }
private:
    std::vector<double> data;
};
template<typename LHS, typename RHS>
VecAdd<LHS, RHS> operator+(const VecExpr<LHS>& lhs, const VecExpr<RHS>& rhs) {
    return VecAdd<LHS, RHS>(static_cast<const LHS&>(lhs), static_cast<const RHS&>(rhs));
}
template<typename E>
VecScale<E> operator*(double scalar, const VecExpr<E>& expr) {
    return VecScale<E>(scalar, static_cast<const E&>(expr));
}
int main() {
    Vector a(3), b(3), c(3);
    a[0] = 1; a[1] = 2; a[2] = 3;
    b[0] = 4; b[1] = 5; b[2] = 6;
    c[0] = 7; c[1] = 8; c[2] = 9;
    Vector result(3);
    result = 2.0 * a + b + c;  // one fused pass, no temporaries
    for (size_t i = 0; i < result.size(); ++i) std::cout << result[i] << ' ';
    std::cout << '\n';  // 16 20 24
}

Performance Comparison

Temporaries Avoided

// Naive operator overloading: each + allocates a new Vector
Vector result = a + b + c;
// 1) operator+(a, b) allocates temp1, computes temp1 = a + b
// 2) operator+(temp1, c) allocates temp2, computes temp2 = temp1 + c
// 3) result = temp2 (copy or move)
// Total: 2 heap allocations + 2 full passes over the data
// Expression templates: builds a tree, evaluates once in operator=
Vector result2 = a + b + c;
// VecAdd<VecAdd<Vector,Vector>,Vector> — no allocation until operator= runs
// Total: 0 extra heap allocations, 1 pass over the data

Benchmark Shape (Illustrative)

Vector size	Naive (ms)	Expression templates (ms)	Speedup
1,000	0.01	0.005	~2x
100,000	1.2	0.4	~3x
10,000,000	180	55	~3.3x

The gap widens with vector size because the naive version’s extra full passes over larger data cost proportionally more, while expression templates always do exactly one pass regardless of chain length.

When the Gains Disappear

For very short chains (a + b alone, no further chaining) or very small vectors, the constant overhead of extra template instantiations can erase the benefit — profile before assuming expression templates are a win in a specific hot path, the same way you would for any other optimization.

Why the Speedup Isn’t Purely About Allocation Count

It is tempting to attribute the entire performance gain to “fewer malloc calls,” but allocation count is only part of the story, and for small vectors it may not even be the dominant part. The second factor — memory traffic — often matters more at scale. The naive three-temporary version of a + b + c + d writes and then re-reads each intermediate result from main memory (or, if you’re lucky, cache) three separate times: once per +. The expression-template version reads each source element once and writes the final output element once, with the whole a[i] + b[i] + c[i] + d[i] computation happening in registers in between. On modern hardware where memory bandwidth, not raw ALU throughput, is frequently the bottleneck for simple element-wise arithmetic, cutting the number of full array passes from three down to one is often the larger contributor to the speedup — which is also why the benchmark table above shows the advantage growing with vector size: at 1,000 elements everything comfortably fits in L1/L2 cache and the gap is modest, but at 10,000,000 elements the naive version is paying real DRAM bandwidth costs on each of its extra passes that the fused version simply never incurs. This also explains why profiling matters more than intuition here: the actual speedup you see depends heavily on your specific hardware’s cache hierarchy, the vector size relative to cache capacity, and whether the compiler successfully auto-vectorizes the fused loop. Numbers measured on one machine and one compiler version are not a promise about another — treat the table above as illustrative of the shape of the effect (allocation-bound and memory-bandwidth-bound workloads benefit more as data grows), not as a number to quote for your own workload.