An efficient way for removing rows with zero elements iteratively while constructing the matrix

Question

Based on a previous post, Accelerating creation of matrices and finding ways for optimal scaling, we managed to accelerate the way that I construct a matrix in Rcpp (inputs for example are at the end of the post). The code used is the following:

#include <RcppArmadillo.h>
    
// [[Rcpp::depends(RcppArmadillo)]]


// [[Rcpp::export]]
int NonDiagSum(arma::mat n, int K){
  int sum = 0; 
  for(int i=0; i<(K+1); ++i){
    sum += n(i,K-i);
  } 
  return accu(n) - sum;
}

std::random_device rd;
std::mt19937 gen(rd());

// [[Rcpp::export]]

arma::mat Generate(arma::mat n, int K, double p){
  
  std::bernoulli_distribution distrib(p);
  
  int N = NonDiagSum(n,K);
  NumericMatrix H(N,K);

  int next = 0;
  
  for(int j=0; j<K; ++j){
    for(int i=0; i<(K-j); ++i){
      for(int iter=0; iter<n(i,j); ++iter){
        
        for(int k=j; k<(K-i); ++k){
          H(iter+next,k) = distrib(gen);   
        }
      
      }
      next += n(i,j);
    }
  }
  return H;
}

The main part of the code is the function Generate(), which produces a matrix with N rows and K columns and each cell can take the value 0 or 1.

Because 0 and 1 are generated in a probabilistic way it is expected to have rows which are entirely of zero elements. My goal, is to remove those rows, in an efficient way.

I think the naive way to approach this problem, is to create an additional function that will calculate the row sum of the matrix produced by the Generate() function, and then create a second function that will remove the rows for which the sum is equal to zero, i.e. row_sum[i]==0. Those function are the following (I assume that those function have the least amount of complexity and cannot be further improved??):

// [[Rcpp::export]]

arma::vec RowSum(arma::mat H, int K){
  
  int N=size(H)[0];
  
  arma::vec s(N);
  s.zeros();
  
  for(int i=0; i<N; ++i){
    for(int k=0; k<K; ++k){
      s[i] += H(i,k);
    }
  }
  return s;
}


// [[Rcpp::export]]


arma::mat Extract(arma::mat H, int K){
  
  arma::vec s = RowSum( H,  K);
  
  int idx = 0;
  
  for(int i=0; i<size(s)[0]; ++i){
    if(s[i]==0){
      H.shed_row(i-idx);
      idx += 1;
    }
  }    
  return H;
}

And I could use the Extract() function inside the Generate() function as

// [[Rcpp::export]]

arma::mat Generate(arma::mat n, int K, double p){

  ... //Same as before

  return Extract(H,K); //Use the Extract() function on the output matrix H
}

which will solve my problem, but it will burden the calculations enormously.

However, I believe that we do not trully need the RowSum() function as those operations are implemented already inside the Generate() function. Hence, the Generate() function which calculates also the row sums is the following:

// [[Rcpp::export]]


arma::mat Generate(arma::mat n, int K, double p){
  
  std::bernoulli_distribution distrib(p);
  
  int N = NonDiagSum(n,K);
  NumericMatrix H(N,K);
  
  NumericVector row_sum(N);
 
  int next = 0;

  for(int j=0; j<K; ++j){
    for(int i=0; i<(K-j); ++i){
      for(int iter=0; iter<n(i,j); ++iter){
        
        
        for(int k=j; k<(K-i); ++k){
          H(iter+next,k) = distrib(gen);
          row_sum(iter+next) += H(iter+next,k);
        }       
      }
      next += n(i,j);
    }
  }
  return H;
}

Then ideally, there are I think two approaches to remove the rows which have zero sum. The first one is to find a way for the removing zero rows after the creation of each row, i.e. after the implamentation of the code part (in Generate() function)

for(int k=j; k<(K-i); ++k){
     H(iter+next,k) = distrib(gen);
     row_sum(iter+next) += H(iter+next,k);
  }  

which personally I cannot find a way at the momment to do that.

The second approach is the following, inside the Generate() function instead of defining H as a NumericMatrix we could define it as arma::mat in order to use the row removing function H.shed_row() function. Based on that the new Generate() code could be

// [[Rcpp::export]]


arma::mat Generate(arma::mat n, int K, double p){
  
  
  ...//Same as before, we construct the matrix H
  
  //Here we remove the rows with zero row_sum
  int idx=0;
  for(int i=0; i<N; ++i){
    if(s[i]==0){
      H.shed_row(i-idx);
      idx += 1;
    }
  }
 
  return H;
}

I assume that the optimal way incorporates keeping the H matrix as NumericMatrix because it avoids the comands arma::mat H and H.zeros() and that the part that we remove the rows should be right after the implementation of the code part

for(int k=j; k<(K-i); ++k){
     H(iter+next,k) = distrib(gen);
     row_sum(iter+next) += H(iter+next,k);
  }   

Sorry for the length of the post, I just wanted to give all the details. If any clarification or something more needed I'm willing to help. Any suggestion would be really helpful!

Potential inputs that can be used for examples are the following:

p = 0.8
K=4
n = matrix(c(0,242,0,272,9222,0,10,0,123,0,0,0,0,0,0,0,131,0,0,0,0,0,0,0,0),K+1,K+1,byrow=TRUE)
n
     [,1] [,2] [,3] [,4] [,5]
[1,]    0  242    0  272 9222
[2,]    0   10    0  123    0
[3,]    0    0    0    0    0
[4,]    0  131    0    0    0
[5,]    0    0    0    0    0

```

Answer 1

Code review is normally meant to review complete and working code, not for getting answers to questions like "how should I ...?" However, since you presented code for all the alternatives you've come up with, I'll review it.

Avoid doing things you have to fix up later

First creating rows which are all zeroes and then removing them is of course more work than not creating those rows at all. So if there is a cheap way to predict that a row will be filled with all zeroes, then we can use that to avoid the work. And there is: for a series of \$n\$ boolean values, each with probability \$p\$ to be 1, the chance of all of them being zero is \$(1 - p)^n\$. So naively thought, you could just check this first before adding a row:

std::bernouilly_distribution empty_row_distrib(std::pow(1 - p, K - i - j));

for (int iter=0; iter<n(i,j); ++iter) {
    if (empty_row_distrib(gen)) {
        /* skip row */
    } else {
        /* generate row */
    }
}

But there are some problems with the above code. First, if we get to the part where we do generate a row, we still might produce a row of all zeroes. So we'd have to check for all zeroes anyway, and to ensure we still have the correct distribution in the end, we'd have to do another attempt at generating the row, until we get one that is not all zeroes. Also, we had to generate yet another random number and compare it. So there is an overhead, and I think this approach is only going to be faster if p is quite small.

Otherwise, I would recommend the following:

Overwrite rows instead of removing them

After adding a row, check if it is all zeroes. If so, don't increment the row index, so the next iteration will just overwrite the same row. When you are at the end, just shrink the matrix using arma::mat::resize() to the number of actual rows written to.

This is not so different from calling shed_row() on individual rows, but this way we just combine everything into a single operation.