Hi,

I was under the impression that Eigen was able to perform some optimizations
in matrix products when selfadjointView() is used. Specifically, I was
trying to compare the time of these two matrix products:

r1 = full.transpose()*sym*full;
r2 = full.transpose()*sym.selfadjointView<Eigen::Upper>()*full;

With Eigen 3.0.2 and MS VC++ 2008, the second product takes slightly
longer than the first. I'm building with Release configuration and whether
I run with SSE2 on or off doesn't seem to affect this comparison.

Am I remembering incorrectly or perhaps this is a platform-specific issue?
I've posted my test code below. Any suggestions appreciated.

Thanks.

Bill

------------------------------------------------
static void symTest2()
{
const int n=20;
Eigen::MatrixXd sym(n, n), full(n,n);
for(int i=0; i<n; i++) {
double val = i+1;
for(int j=0; j<n; j++) {
sym(i,j) = val;
sym(j,i) = val;
full(i,j) = i*2 + j;
full(j,i) = j*2 + 3*i;
}
}

const int num=1000000;
boost::timer timer;
Eigen::MatrixXd r1, r2;
for(int i=0; i<num; i++) {
r1 = full.transpose()*sym*full;
}
printf("elapsed time for full=%8.2f\n", timer.elapsed());
timer.restart();
for(int i=0; i<num; i++)
r2 = full.transpose()*sym.selfadjointView<Eigen::Upper>()*full;
printf("elapsed time for sym=%8.2f\n", timer.elapsed());
std::cout << "diff=\n" << r1-r2 << std::endl;
}

In this case the .selfadjointView<Eigen::Upper>() only tells Eigen to use only the upper half and consider it as a selfadjoint matrix. This does not reduce the operation count. Eigen cannot know that your result will be selfadjoint. You can tell it by doing:

r2.triangularView<Upper>() = full.transpose()*sym.selfadjointView<Eigen::Upper>()*full;

that means: compute the upper triangular view only. Then you can use the result with:

r2.selfadjointView<Eigen::Upper>()

Thanks very much, Gael.
By following your instructions,

r2.triangularView<Eigen::Upper>() = full.transpose()*sym.selfadjointView<Eigen::Upper>()*full;

executes in about 1/2 the time of the full case, as expected.

But I'm still confused about how to use the result.
After the execution of this product, r2.rows() and r2.cols() are both
equal zero so
r2.selfadjointView<Eigen::Upper>()
is also empty.

In my code, I have this declaration for r2:
Eigen::MatrixXd r2;

If instead, I declare r2
Eigen::MatrixXd r2(n,n)

the contents of r2 after the matrix product are as expected.
But the two forms of matrix
product now execute in approximately the same time.

I've also scanned the Eigen test cases but haven't found one that
is like mine.

I'll appreciate any insights.

Bill

hm, in debug mode you should get an assertion. You have to resize r2 yourself. It is like doing r2.block(....) = something;...

>hm, in debug mode you should get an assertion.
I was running with optimization turned on but after rebuilding
in debug mode, the operation completes without an error-- just that
the LHS matrix is empty.

>You have to resize r2 yourself. It is like doing r2.block(....)
OK, I see. But as I said, when I declare r2 as
Eigen::MatrixXd r2(n,n)
the matrix operation using selfadjointView takes just as long as
with the full matrix. And r2 contains terms in both upper and lower
triangles-- which I guess explains why the computation takes as long
as for the full matrix case.

Bill

Eigen::MatrixXd r2(n,n);
r2.setZero();
r2.triangularView<Eigen::Upper>() = full.transpose()*sym.selfadjointView<Eigen::Upper>()*full;

should really works, with the lower half of r2 being full of 0. In the best case you expect to cut down the computation time by 0.75, not more. Also note that the gain might not be visible for small matrices.

after some experiments I get a factor of 0.75 on 1000^2 matrices, and 0.9 for 20x20 matrices.

Thanks very much for going to the trouble to take a close look at this.
Clearly I was mistaken about how much performance improvement to expect.

Bill