본문으로 건너뛰기 The Complete Guide to C++ Expression Templates | Lazy Eva...

The Complete Guide to C++ Expression Templates | Lazy Evaluation and Mathematical Library Optimization

The Complete Guide to C++ Expression Templates | Lazy Evaluation and Mathematical Library Optimization

이 글의 핵심

A complete guide to Expression Templates. Learn about lazy evaluation, eliminating temporary objects, optimizing vector operations, and implementing Eigen-style libraries.

originalId: cpp-expression-template

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

Table of Contents

  1. Basic Structure
  2. Implementing Vector Operations
  3. Matrix Operations
  4. Common Errors and Solutions
  5. Production Patterns
  6. Complete Example: Mathematical Library
  7. Performance Comparison

1. 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 = ....

2. 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
}

3. 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.

4. 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.

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.

5. 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.

6. 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
}

7. 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 sizeNaive (ms)Expression templates (ms)Speedup
1,0000.010.005~2x
100,0001.20.4~3x
10,000,00018055~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.

Other articles related to this topic.



Keywords Covered in This Article (Related Search Terms)

This article covers C++, expression-template, template, optimization, lazy, eigen.