Stan::math's quad_form and Eigen matrix multiplication

General
1

I have a question regarding the stan::math implementation of quad_form

What interests me - in the matrix-matrix version - is the return b.transpose() * a * b; return part of what looks to be something like a C++ lambda (it looks like it takes a pointer of the arguments so there is probably some cleaning up happening behind the scenes, I imagine this is useful for handling autodiff).

The question I have how Stan handles the matrix multiplication here. I imagine since you are using Eigen matrices here, then I suppose it is really just handling it how Eigen does it. Does it use expression templates here, or is that just for arithmetic operations? Does it do anything special to cut down on temporaries when there are multiple matrix multiplications at once?

2

The return wrapper make_holder() is used because Eigen operations return expression templates, rather than fully evaluated results. The risk with returning expression templates is that particular objects (matrices/vectors/scalars) might fall out of scope and no longer be available when all of the expression templates are finally evaluated (discussed more in this Eigen FAQ).

The make_holder() wrapper is used to make sure that the objects associated with the returned expression template do not fall out of scope before it is evaluated

3

@andrjohns Thanks, but I’m not sure my focus is on make_holder. I suppose what I’m really interested in is how the expression b.transpose() * a * b is actually evaluated by Eigen (when it is evaluated). If I understand this explanation on lazy evaluation here, there is special treatment for matrix products. It is passed around as an expression, but when you evaluate them, it assigns the result to a temporary and then copied the result over to the original matrix. My question is whether anything special happens when they chain the matrix products. The 2nd and 3rd circumstances on that link suggest that since it is a costly operation then they store the results. This suggests that if you evaluate something like `result = b.transpose() * a * b’ it would get re-written like

tmp1 = b.transpose() * a;
tmp2 = tmp1 * b;
result = tmp2.copy;

So that there are two temporaries. But I suppose I wanted to make sure I was understanding it correctly.

4

There is not a single answer for how this evaluated, since Eigen has multiple pathways for evaluating its expression templates. The given pathway is determined by the expression itself as well as the size of the matrices/vectors under manipulation.

If you want a specific answer on how Eigen handles this particular operation, rather than how it’s used in the Math library, then you’ll need to ask the Eigen developers themselves. The easiest way to do this is to ask a question on StackOverflow with the [eigen] or [eigen3] tags

1 Like
5

Thanks