Does A * A.transpose() is implemented to do half of the operations instead of going trough all matrix indexes? If not, can someone point me the best approach to implement it using Eigen?

computing only one half is worth it for large enough matrices only, you can do:

M.setZero();
M.selfadjointView<Lower>().rankUpdate(A);

that is enough for products, cholesky, etc. If you want a complete M:

M.triangularView<StrictlyUpper>() = M.transpose();

I have a very large Nx6 matrix and solving this system is being the bottleneck of my app. I think it worth the effort.

I'm not sure how to use your piece of code in my case. Currently I'm doing

Matrix<float, Dynamic, 6> A;
Matrix<float, Dynamic, 1> b;
Matrix<float, 6, Dynamic> At;

At = A.transpose();

x = (At * A).llt().solve(At * B);

I think that optimizing At * A would greatly increase the performance of this code, am I right?

No need to evaluate At, simply do:

Matrix<float, 6, 6> AtA;
AtA.setZero();
AtA.selfadjointView<Lower>().rankUpdate(A);
Matrix<float, 6, 1> x = AtA.selfadjointView<Lower>().llt().solve(A.transpose()*B);

make sure you compile with optimizations (-O2), Eigen is very slow otherwise.

Tried your approach and didn't worked. There is no error messages, but the system is not solved.
Im on Windows 7 and using MSVC2010 environment, compiling as 32bit. Any thoughts?

what do you mean by "the system is not solved" ?


(AtA.selfadjointView<Lower>() * x - A.transpose()*B).norm()/B.norm() should be small, and same for (A*x-B).norm()/B.norm();

System is not solved means that it don't returns the same result as
(At * A).llt().solve(At * b). The returning x is completely wrong using the selfAdjointView expression. To be honest, the returning value looks a lot with memory garbage.

You are telling me that it will just works if
(AtA.selfadjointView<Lower>() * x - A.transpose()*B).norm()/B.norm() or (A*x-B).norm()/B.norm() are small?

I've made a simple test here. The expression:

AtA.setZero();
AtA.selfadjointView<Eigen::Lower>().rankUpdate(A);
AtA.triangularView<Eigen::StrictlyUpper>() = AtA.transpose();

Should return exactly the same as
AtA = At * A;

Am I right?

If that's the case, there is some mistake here, because they are not equal in my test.

sorry, an update operaton U U^T, not U^T U, so you need to do:

AtA.selfadjointView<Eigen::Lower>().rankUpdate(A.transpose());

The previous code should give you an assert when compiling without -DNDEBUG.