Hello,

Code: Select all

Eigen::SparseMatrix<float> A; // large sparse system with 0 and 1 values only
Eigen::MatrixXf B1, B2;
// the operations here takes most of the time so i was wondering if there might be some
// redundant temp copies
auto At = A.transpose();
Eigen::MatrixXf AtA = (At * A).eval(); // if i dont use eval() it wont compile
B2 = At * B1;


on the other hand LDLT decomposition and solve runs pretty fast.

Code: Select all

Eigen::LDLT<Eigen::MatrixXf> solver(AtA);
solver.solveInPlace(B2);


I see that you are storing the result of the A'A into a dense matrix, so:
- How much sparse is A ? (A.nonZeros())
- What about A'A, is it still sparse? ((AtA.array()!=0).count())
- What are the sizes of A? (A.rows(), A.cols())

Hi

A = 326500 x 10880
A-non zero = 326500 x 680
A'A-non zero = 118151144 (so its really dense)

Perhaps try Eigen::SparseMatrix<int> A if they are guaranteed integers. You can avoid using AtA (I think) doing this,

Code: Select all

B2 = (At * A).eval(); // if i dont use eval() it wont compile
B2 = B2 * B1;


In the devel branch, sparse product can be directly evaluated within a dense result, thus saving some computation and memory accesses to deal with the sparsity. For your pb size, AtA takes 210s instead of 270s. Moreover, in your case, using LLT should be enough, and for such problem size, LLT is one order of magnitude faster (5s instead of 50s) because it implements a cache-friendly algorithm.

@ggael were the changes u are talking about included in last 3.2.7 release?
do i need to change the way i create dense matrix out of (At*A) ?
thanx

It is in the devel branch only, and cannot be backported to 3.2. WritingMatrixXf AtA = (At * A) will do the job.

i was able to avoid additional allocations with following syntax:
AtA.noalias() = Eigen::MatrixXf(At * A);

i'm using the fresh clone of bitbucket repo.

To avoid one dense temporary you really have to write:

AtA = At * A;

There is no need for noalias because a dense matrix cannot alias with sparse matrices.

After updating to 3.2.7 i'm getting following error:
SparseSparseProductWithPruning.h:50:7: error: no member named 'reserve' in 'Eigen::Matrix<float, -1, -1, 0, -1, -1>'
res.reserve(estimated_nnz_prod);
~~~ ^
Eigen/src/SparseCore/SparseSparseProductWithPruning.h:91:15: note: in instantiation of function template specialization 'Eigen::internal::sparse_sparse_product_with_pruning_impl<Eigen::SparseMatrix<float, 0, int>, Eigen::SparseMatrix<float, 0, int>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >' requested here
internal::sparse_sparse_product_with_pruning_impl<Lhs,Rhs,ResultType>(lhs, rhs, _res, tolerance);
^
Eigen/src/SparseCore/SparseProduct.h:121:94: note: in instantiation of member function 'Eigen::internal::sparse_sparse_product_with_pruning_selector<Eigen::SparseMatrix<float, 0, int>, Eigen::SparseMatrix<float, 0, int>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, 0, 0, 0>::run' requested here
internal::sparse_sparse_product_with_pruning_selector<_LhsNested, _RhsNested, Dest>::run(lhs(),rhs(),result,m_tolerance);
^
Eigen/src/Core/Assign.h:529:107: note: in instantiation of function template specialization 'Eigen::SparseSparseProduct<Eigen::SparseMatrix<float, 0, int>, const Eigen::SparseMatrix<float, 0, int> &>::evalTo<Eigen::Matrix<float, -1, -1, 0, -1, -1> >' requested here
static EIGEN_STRONG_INLINE Derived& evalTo(ActualDerived& dst, const ActualOtherDerived& other) { other.evalTo(dst); return dst; }
^
Eigen/src/Core/Assign.h:578:65: note: in instantiation of function template specialization 'Eigen::internal::assign_selector<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::SparseSparseProduct<Eigen::SparseMatrix<float, 0, int>, const Eigen::SparseMatrix<float, 0, int> &>, false, false>::evalTo<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::SparseSparseProduct<Eigen::SparseMatrix<float, 0, int>, const Eigen::SparseMatrix<float, 0, int> &> >' requested here
return internal::assign_selector<Derived,OtherDerived,false>::evalTo(derived(), other.derived());
^
Eigen/src/Core/PlainObjectBase.h:483:13: note: in instantiation of function template specialization 'Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >::operator=<Eigen::SparseSparseProduct<Eigen::SparseMatrix<float, 0, int>, const Eigen::SparseMatrix<float, 0, int> &> >' requested here
Base::operator=(other.derived());
^
Eigen/src/Core/Matrix.h:184:20: note: in instantiation of function template specialization 'Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >::operator=<Eigen::SparseSparseProduct<Eigen::SparseMatrix<float, 0, int>, const Eigen::SparseMatrix<float, 0, int> &> >' requested here
return Base::operator=(other);
^
note: in instantiation of function template specialization 'Eigen::Matrix<float, -1, -1, 0, -1, -1>::operator=<Eigen::SparseSparseProduct<Eigen::SparseMatrix<float, 0, int>, const Eigen::SparseMatrix<float, 0, int> &> >' requested here
AtA = At * A;
----------------------------
AtA -- is dense
At/A -- sparse