sparse sparse product A^T*A optim in Eigen lib

1.1k views Asked by At

In the case of multiple of same matrix matA, like

matA.transpose()*matA, 

You don't have to compute all result product, because the result matrix is symmetric(so only if the m>n), in my specific case is always symmetric! square.

So its enough the compute only for. ex. lower triangular part and rest only copy..... because the results of the multiple 2nd and 3rd row, resp.col, is the same like 3rd and 2nd.....And etc....

So my question is , exist way how to tell Eigen, to compute only lower part. and optionally save to only lower trinaguler part the product?

    DATA = SparseMatrix<double>((SparseMatrix<double>(matA.transpose()) * matA).pruned()).toDense();
2

There are 2 answers

5
Avi Ginsburg On BEST ANSWER

According to the documentation, you can evaluate the lower triangle of a matrix with:

m1.triangularView<Eigen::Lower>() = m2 + m3;

or in your case:

m1.triangularView<Eigen::Lower>() = matA.transpose()*matA;

(where it says "Writing to a specific triangular part: (only the referenced triangular part is evaluated)"). Otherwise, in the line you've written Eigen will calculate the entire sparse matrix matA.transpose()*matA.

Regarding saving the resulting m1 matrix, it is the same as saving whatever type of matrix it is (Eigen::MatrixXt or Eigen::SparseMatrix<t>). If m1 is sparse, then it will be only half the size of a straightforward matA.transpose()*matA. If m1 is dense, then it will be the full square matrix.

0
Charlie S On

https://eigen.tuxfamily.org/dox/classEigen_1_1SparseSelfAdjointView.html

The symmetric rank update is defined as:

B = B + alpha * A * A^T

where alpha is a scalar. In your case, you are doing A^T * A, so you should pass the transposed matrix instead. The resulting matrix will only store the upper or lower portion of the matrix, whichever you prefer. For example:

SparseMatrix<double> B;
B.selfadjointView<Lower>().rankUpdate(A.transpose());