CODEWITHSUNDEEP
sparse-matrixeigen
yannick
yannick

Reputation: 397

remove empty rows of an Eigen::SparseMatrix

I have built a sparse matrix mat from a list of triplets

Eigen::SparseMatrix<double, Eigen::RowMajor> mat(Nbins,Ndata);
mat.setFromTriplets(tripletList.begin(), tripletList.end());

Now I would like to create a new matrix ret, which only contains the rows of the previous matrix which are not empty. I do it as follows

Eigen::SparseMatrix<double, Eigen::RowMajor> ret(Nbins,Ndata);
unsigned Nrow=0;
for (unsigned i=0; i<Nbins; ++i) {
  auto mrow = mat.row(i);
  if (mrow.sum()>0) {
    ret.row(Nrow++) = mrow;
  }
}
ret.conservativeResize(Nrow,Ndata);

However, doing it this way is slow and inefficient. Slow because quick profiling suggests it spends most of its time on ret.row(Nrow++) = mrow;. Inefficient because we are also copying all the data twice.

Is there a better solution? I feel one has to fiddle with the inner vectors but I get confused by them and I don't know how user-proof it is to play with them.

EDIT: In my application, matrices are row major, and I want to remove empty rows. mat is not needed, just ret. All coefficients are positive hence the way I check for nonzero rows. The triplets are sorted but column-major. There are no duplicate triplets.

Upvotes: 2

Views: 1066

Answers (1)

yannick
yannick

Reputation: 397

Found it! Instead of writing a hand-made setFromTriplets, I went with a modification of the tripletList. The interface of Eigen::Triplet makes it very easy.

//get which rows are empty
std::vector<bool> has_value(Nbins,false);
for (auto tr : tripletList) has_value[tr.row()] = true; 

//create map from old to new indices
std::map<unsigned,unsigned> row_map;
unsigned new_idx=0;
for (unsigned old_idx=0; old_idx<Nbins; old_idx++) 
   if(has_value[old_idx])
      row_map[old_idx]=new_idx++;

//make new triplet list, dropping empty rows
std::vector<Eigen::Triplet<double> > newTripletList;
newTripletList.reserve(Ndata);
for (auto tr : tripletList) 
   newTripletList.push_back(
      Eigen::Triplet<double>(row_map[tr.row()],tr.col(),tr.value()));

//form new matrix and return
Eigen::SparseMatrix<double, Eigen::RowMajor> ret(new_idx,Ndata);
ret.setFromTriplets(newTripletList.begin(), newTripletList.end());

Upvotes: 1

Related Questions