Why expression templates matter ?

Imagine we are using a matrix library, with a naive implementation of +, that simply adds column and row-wise each coefficient. Thus, for a, m1, m2, m3 being for example square matrices of order n, if we want to compute the sum of m1, m2 and m3 and to put the result in a, corresponding to the following code :

a = m1 + m2 + m3

then the computation will look like the following tree :

Evaluation tree of m1 + m2 + m3

Thus, there will be a first for loop to compute m1 + m2 which will result in a matrix we’ll call t, then another one to compute t + m3. For 3×3 matrices, it’s ok. Imagine n is actually 40000. Doing two for loops of 40000 iterations where we would do a single one… quite annoying, isn’t it ? Actually, we would rather want an evaluation tree like the following :

The desired evaluation tree of m1 + m2 + m3

Here comes expression templates 😉

It’ll need a bit of refactoring though. First, let’s say our matrix type is the following (we’ll only deal with float to avoid ambedding an additional typename template parameter representing the number type we use) :

Now, we’ll have to define a Domain Specific Embedded Language for matrix operations. Like in any language people design, we will have a tree representing what’s happening in the code. In our case, it’ll represent the evaluation tree of the matrix operations we’re dealing with. Thus, like in any expression tree, we need an Expression type. Ours will look like this :

Thanks to the definition of Expression, we can embed matrix operations in Expressions, but for the moment we can’t do the reverse way, that is we can’t convert an Expression to a matrix. So let’s write an operator= in the matrix class.

So you see, this is the only moment where we call the operator()(unsigned int, unsigned int) of the Expression type. Thus, this is the only moment when the whole expression is being evaluated. Let’s study the following code.

instance. The coefficients aren’t computed until operator= is called. Indeed, we have the following expression tree.

Expression tree

That is, we know which operations we have to call, on which matrices, but nothing is evaluated. The only place where we call operator() (unsigned int, unsigned int) on a value of type Expression is in matrix<N>::operator=, and this call actually computes each coefficient, one by one, applying the whole computation tree (here, two calls to ‘+’) for each coefficient. This way, we only execute the two for loops once, and it’ll remain the same whatever the number of computations is. Moreover, the only change we had to make to our matrix class was to add an operator= to be able to assign an expression to a matrix.

Yeah indeed. My code is just a beginning. We should have a temporary object created in the operator=, and then putting this one’s data in this’ data.

I consider writing about expression templates / DSELs for a less common task, in this one I’ll take care of avoiding such “bugs”. By the way, if you have any idea for the “domain” of the DSEL, don’t hesitate. And thanks for reading & commenting !

[…] Why expression templates matter ? « Alp Mestanogullari's Blog alpmestan.wordpress.com/2009/11/01/why-expression-templates-matter – view page – cached Imagine we are using a matrix library, with a naive implementation of +, that simply adds column and row-wise each coefficient. Thus, for a, m1, m2, m3 being for example square matrices of order n,… (Read more)Imagine we are using a matrix library, with a naive implementation of +, that simply adds column and row-wise each coefficient. Thus, for a, m1, m2, m3 being for example square matrices of order n, if we want to compute the sum of m1, m2 and m3 and to put the result in a, corresponding to the following code : a = m1 + m2 + m3 (Read less) — From the page […]

Awesome post! I have been using Eigen recently, a linear algebra library, and they stated support for expression templates. I didn’t know what that meant, but I assumed it was along the lines of this. It is nice to see a simple and clear implementation, thanks.

Nice post. But a problem arise when you decide to do operations on generic types: how to choose the right type between matrix and matrix? Traits can bring an elegant solution to this problem but I don’t know if it’s possible to reduce the number of classes (e.g for symetric cases).

You can create a matrix class being built combinating many “blocks” (the storage, etc), well something like a policy-based matrix, or policy/traits-based, something like that. But then we rather enter the generative programming world, totally, whereas here it was just about having an AST for a minimal DSEL for matrixes. Very minimal.