Hi,
I'm working on a non-linear least squares solver which requires the computation of the product of a sparse matrix and its transpose. The Jacobian for this problem is fairly sparse (on average, 7% of the elements are filled in each row). Specifically, I'm computing J_transpose * J. The matrix J is 9971 x 584 in this case.

While profiling my code, I noticed that the majority of the time (75% of execution time) was spent in computing this sparse matrix product. This was somewhat unexpected - I thought the Cholesky decomposition or the construction of the Jacobian would be the bottleneck. I'm using the latest stable release of Eigen (3.1.3) and Visual Studio 2012.

Am I doing something wrong? I've tried this by making lhs dense, but that didn't make any difference. It seems like the product itself takes a while to compute. Is there any way of speeding up this sparse matrix product?

Thanks a lot!

The relevant code snippet:

Code: Select all

      
Eigen::SparseMatrix<Float> lhs(this->numFreeVariables, this->numFreeVariables);
Eigen::SparseMatrix<Float, Eigen::RowMajor, int> jacobian;

...

lhs = (this->jacobian.transpose() * this->jacobian).pruned(1e-6);         



The right hand side "this->jacobian" has to be evaluated into a temporary column major matrix to be able to perform the sparse-sparse product. It might be that this transpose operation is the bottleneck. You can check by explicitly computing it:

Code: Select all

Eigen::SparseMatrix<Float> tmp(this->jacobian);
lhs = this->jacobian.transpose() * tmp;


You might also compare with and without the .pruned(). The underlying algorithms are not the same. It would also be interesting to know what are the number of non zeros into the resulting lhs as well as the running times you are observing. I also assume you tested in release mode with the compiler optimizations on.

Finally, if your goal is to do least-square solving, you might try the SparseQR class in the devel branch which can directly solve jacobian * x = b without forming the normal equation.

Thanks very much, Gael. I'll get back to you with details on the number of non-zeros inserted into the lhs.

The transpose operation isn't terribly expensive (~4% of the run time). Both pruned and unpruned versions have similar run times. I didn't realize that there was a sparse QR in the devel branch. I'll give that a shot today.