My code computes many quadratic forms of the form

Code: Select all

VectorXd x; MatrixXd M;
double n = x.transpose() * M.selfadjointView<Eigen::Lower>() * x;


and my profiling tests are hinting that a temporary is being created, i.e.

Code: Select all

VectorXd temp = M..selfadjointView<Eigen::Lower>() * x;
n = x.transpose() * temp;


This seems like the expected behavior, and a sterling example of using noalias -- the problem is that I can't call noalias on the double and I'd prefer not to create an Eigen object (Vector1d?).

Is there a way to force a noalias on the RHS or do I need to so something like

Code: Select all

Vector1d t;
n = t.noalias() = x.transpose() * M.selfadjointView<Eigen::Lower>() * x;


Thanks!

-Mike

what are the typical sizes of your matrices ? Unless they are very small, this temporary is really needed such that efficient matrix*vector product can take place. If your matrices are very small you can try to enforce the use of a lazy product:

x.dot(M.selfadjointView<Eigen::Lower>().lazyProduct(x))

you might also try:

x.dot(x.transpose().lazyProduct(M.selfadjointView<Eigen::Lower>()))

(the best choice depends on the storage order).

Moreover, if your matrices are small it might be a good idea to compute the full M such that vectorization can take place, especially if it used multiple times.

Please report your findings on it !!

The matrices are small for these tests, but eventually will be large when the code is applied to real problems. I guess I had assumed that it would be easy enough to vectorize the quadratic form a la

Code: Select all

lhs = 0

for(int i = 0; i < n; ++i)
{
    lhs += x(i) * ( M.row(i) * x);
}


depending on the storage order, etc.

Now I do have auxiliary vectors instantiated in the code for other purposes, but they could be used here. If the optimizations require that the matrix-vector product be explicitly computed would it be best to simply use the existing vectors,

Code: Select all

aux = M * x;
n = x.dot(aux);


to avoid using the heap?

Thanks for the help!

-Mike

yes for large enough matrices it is considerably faster to evaluate into a temporary, but you have to use noalias to really avoid heap allocation:

aux.noalias() = M * x;

Right, forgot the noalias on the temp. Thanks again for the help and the great code.

-Mike