The cost model¶

Some part of an expression is read many times. In an einsum, an operand is read once for every value of the letters it lacks; in a broadcast, once for every position it is broadcast over. For each such part an evaluation has two ways to go. It can recompute the part at every read, or compute it once into a buffer and read that. Which is better depends on how dear the part is and how often it is read. strata decides with a cost model, at compile time, from the types alone, and nothing in your code says so.

The decision never changes what you get, only how fast, and how much memory a temporary uses. That is why it is allowed to be a heuristic.

Every code block on this page is part of examples/11_scratch.cpp or examples/14_cost_model.cpp, which are built and run with the tests. The output shown is theirs.

The question¶

Take exp(a / 10) as an operand of a matrix product with b. Each element of it is read once for every column of b:

// exp(a / 10) is read once for every column of b: computed into a buffer first, 12 calls of exp
// rather than 60. `materializes` says so, for a static_assert where a buffer is expected; a + a,
// read once, never needs one.
const auto product = einsum<"ij,jk->ik">(exp(a * 0.1), b);
static_assert(materializes<decltype(product)>);
static_assert(!materializes<decltype(a + a)>);
std::println("exp(a / 10) b:\n{::.2f}", Tensor{product});

exp(a / 10) b:
[[ 0.76,  5.44, 10.11, 14.79, 19.47],
 [ 0.77,  5.49, 10.21, 14.94, 19.66],
 [ 0.78,  5.55, 10.32, 15.09, 19.86]]

a has 12 elements and b has 5 columns. Recomputing exp at each read is 60 calls; computing it into a buffer first is 12. materializes<E> is the question asked as a static_assert: does evaluating E compute a part of it into a buffer first? And a + a, read once, never needs one.

The rule¶

Per element of the part, with \(r\) reads of each element, the evaluation buffers it when

where \(c_{\text{fill}}\) is what it costs to compute one element into the buffer, \(w\) the cost of writing it there, \(c_{\text{read}}\) the cost of reading it back, and \(c_{\text{inline}}\) what an element costs when it is recomputed. Both \(w\) and \(c_{\text{read}}\) are 1. Both sides are times the size of the part, which cancels, so the rule is per element. The inequality is strict, so a tie does not pay for a buffer, and read once, a part never does.

The costs are in the units of Eigen's cost model: reading a value, adding or multiplying two are 1.

What an element costs¶

The functions have their own prices:

So, for example, \(\exp(m)\) costs \(14 + 1 = 15\) an element, \(m + m\) costs \(1 + 1 + 1 = 3\), and \(\sqrt{m} \cdot 2\) costs \(5 + 1 + 1 + 0 = 7\), the scalar being free. A sum along an axis of 4 costs \(4 \cdot (1 + 1) = 8\).

A dynamic extent is taken as 64, wherever a number is needed. Per element that gives the answer a polynomial in the real extents would give, since both the cost and the reuse grow with them.

How often a part is read¶

Element-wise nodes and reductions read each element of an operand once. An operand of an einsum is read once for every position of the letters it does not have: \(a_{ij}\) in "ij,jk->ik" lacks \(k\), so it is read \(n_k\) times, a dynamic \(n_k\) counting as 64. When parts are nested, the reads multiply on the way down from the root.

The numbers¶

With static extents the counts are small and the decision is close. In a product of a \(3 \times 4\) matrix by a \(4 \times 2\) one, each element of the left operand is read twice:

With a dynamic extent, \(r = 64\) and even a sum is worth it: \(3 + 1 + 64 = 68 < 3 \cdot 64 = 192\). Read once, nothing is, however dear. The tests check all of these figures.

Seeing it¶

A function of your own can name its price with a static constexpr double cost. This one also counts its calls, so the decision can be watched from outside:

// A function of our own that counts its calls, and names what one costs: in Eigen's units, where
// reading a value or adding two is 1.
template <double Cost>
struct Counted {
    static constexpr double cost = Cost;
    static inline std::size_t calls = 0;

    auto operator()(double x) const -> double {
        ++calls;
        return x * 2.0;
    }
};

using Cheap = Counted<1.0>;
using Dear = Counted<30.0>;

// A 3 x 4 matrix times a 4 x 2 one: each of the 12 elements of `a` is read once for each of the
// 2 columns of `b`.
const Tensor<double, Shape<3, 4>> a = strata::full(Shape<3, 4>{}, 1.0);
const Tensor<double, Shape<4, 2>> b = strata::full(Shape<4, 2>{}, 1.0);

const Tensor<double, Shape<3, 2>> cheap = einsum<"ij,jk->ik">(map(Cheap{}, a), b);
const Tensor<double, Shape<3, 2>> dear = einsum<"ij,jk->ik">(map(Dear{}, a), b);
std::println("a function costing 1: {} calls", Cheap::calls);
std::println("a function costing 30: {} calls", Dear::calls);
std::println("the same result: {} and {}", cheap[0, 0], dear[0, 0]);

a function costing 1: 24 calls
a function costing 30: 12 calls
the same result: 8 and 8

The cheap function is recomputed: each of the 12 elements is read twice, 24 calls. The dear one is computed once into a buffer, 12 calls. The result is the same. The same decisions as static_asserts:

// Read twice, a function costing 30 is worth a buffer; one costing 1 is not.
static_assert(!materializes<decltype(einsum<"ij,jk->ik">(map(Cheap{}, a), b))>);
static_assert(materializes<decltype(einsum<"ij,jk->ik">(map(Dear{}, a), b))>);

// Read once, nothing is worth a buffer, however dear.
static_assert(!materializes<decltype(strata::exp(a + a))>);

// A dynamic extent is taken as 64: read that often, even a sum is worth a buffer.
const Tensor<double, Shape<Dyn, Dyn>> x{Shape<Dyn, Dyn>{3, 4}, 1.0};
const Tensor<double, Shape<Dyn, Dyn>> y{Shape<Dyn, Dyn>{4, 2}, 1.0};
static_assert(!materializes<decltype(einsum<"ij,jk->ik">(a + a, b))>);
static_assert(materializes<decltype(einsum<"ij,jk->ik">(x + x, y))>);

Where the buffers' memory comes from¶

By default a buffer is on the heap, aligned. scratch chooses another place: an allocator, a std::pmr memory resource such as a stack arena, or a memory policy. The nearest scratch above a part gives its memory. The 11_scratch program in the examples shows it, including an arena on the stack.

What the model does not know¶

It is a first version, and it says so in its code.

• The prices are Eigen's, to be replaced by measurements of strata's own. Only the functions in the table above are priced; the rest of <cmath> costs 1.
• A buffer costs the same at any size. A buffer larger than the cache is dearer to read in practice, and the model does not see it.
• A dynamic extent is always 64, never the real one. A dynamic shape that happens to be tiny can be buffered when a recompute would have been cheaper.
• What a reducer does once at the end, such as the division of a mean, is not counted.

For that reason what materializes answers may change as the model is tuned. Use it in a static_assert to guard a case you care about, in a test or a benchmark, not as a promise about every expression.

See also¶

• Evaluation and assignment in the API reference.
• einsum: the main source of parts that are read many times.
• Expressions, Assignment and aliasing.