Accelerating creation of matrices and finding ways for optimal scaling

Question

I have the following specific task, for which I have written code, in Rcpp, but it doesn't scale well and I would like to improve it. I'll describe first the steps of the procedure

Step 1:

We have matrices of dimensions (K+1)x(K+1), which have non zero values on the positions displayed, and this can be generalized for larger K values

K = 3
n = matrix(c(10,10,10,0,10,10,0,0,10,0,0,0,0,0,0,0),K+1,K+1,byrow=TRUE)
n
      [,1] [,2] [,3] [,4]
[1,]   10   10   10    0
[2,]   10   10    0    0
[3,]   10    0    0    0
[4,]    0    0    0    0


K = 4 
n = matrix(c(10,10,10,10,0,10,10,10,0,0,10,10,0,0,0,10,0,0,0,0,0,0,0,0,0),K+1,K+1,byrow=TRUE)
n
      [,1] [,2] [,3] [,4] [,5]
[1,]   10   10   10   10    0
[2,]   10   10   10    0    0
[3,]   10   10    0    0    0
[4,]   10    0    0    0    0
[5,]    0    0    0    0    0

Step 2:

Now we focus on the non-zero cells. Each cell corresponds to a vector of 0 and 1 values.

For example, for K=3 the n[1,1] cell corresponds to the vector c(1,1,1), where the n[2,1] to the vector c(1,1,0) and the cell n[3,1] to the vector c(1,0,0).

Now if we move to the cells of the second column, the value 1 starts from the second position of the vectors, i.e. n[1,2] has as vector the c(0,1,1), n[2,2] has as vector the c(0,1,0) and lastly in a similar manner n[1,3] has c(0,0,1).

The same happens when K=4, the procedure remains the same the only thing that changes is the vector length, i.e. n[1,1] has as vector c(1,1,1,1).

Hence, each cell corresponds to a particular 0-1 vector. What the value of each cell indicates, is the multiplicity of each corresponding vector. So, if K=3 and n[1,1]=10 then we have 10 times the vector c(1,1,1).

So, I created the multiplicity matrix with the following code, I would stick into the K=3 case for the illustration. First, I created the vector-matrix Vector_Matrix which holds all the possible vectors, as described previously c(1,1,1),c(1,1,0),...,c(0,0,1). Then I create the Multiplicity_Matrix, which is created based on the Vector_Matrix by replicating each rows n[i,j] times.

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


bool tell(int x, int y)
{
  return x >= y;
}

// [[Rcpp::export]]

arma::mat Vector_Matrix(int K){
 
  arma::mat P(K*(K+1)/2,K);
  
  int Y = 0;
  for (int step = K; step > 0; --step) {
    for (int y = 0; y < step; ++y, ++Y) {
      for (int x = 0; x < step; ++x) {
        P(Y,x) = tell(x,y);
      }
    }
  }
  return P;
}

// [[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;
}


// [[Rcpp::export]]

arma::mat Multiplicity_matrix(arma::mat n, int K){
  
  arma::mat Res(NonDiagSum(n,K),K);
  arma::mat P = Vector_Matrix(K);
  int k = 0;
  int ind_r = 0;
  int next = 0;
  for(int i=0; i<K; ++i){
    for(int j=0; j<(K-k); ++j){
      
      if( (n(i,j)!=0) ){
        
        for(int t=0; t<K; ++t){  
          for(int iter=0; iter<n(i,j); ++iter){
            
            Res(iter+next,t) = P(ind_r,t);
          }
        }
        next = next + n(i,j);
        }      
      ind_r = ind_r + 1;     
   }
    k = k + 1;    
  }
  return Res;
}

The results from those algorithms would be

K = 3
Vector_Matrix(K)
      [,1] [,2] [,3]
[1,]    1    1    1
[2,]    0    1    1
[3,]    0    0    1
[4,]    1    1    0
[5,]    0    1    0
[6,]    1    0    0


n = matrix(c(2,2,2,0,2,2,0,0,2,0,0,0,0,0,0,0),K+1,K+1,byrow=TRUE)

Multiplicity_Matrix(n,K)
       [,1] [,2] [,3]
 [1,]    1    1    1
 [2,]    1    1    1
 [3,]    0    1    1
 [4,]    0    1    1
 [5,]    0    0    1
 [6,]    0    0    1
 [7,]    1    1    0
 [8,]    1    1    0
 [9,]    0    1    0
[10,]    0    1    0
[11,]    1    0    0
[12,]    1    0    0

where we see that each row of the Vector_Matrix has been replicated two times, because of the cells n[i,j]=2.

Step 3:

Now as a third step, where the output of that step is our main goal, we have a Bernoulli trial on each non-zero cell of the Multiplicity_Matrix and we save the result matrix.

// [[Rcpp::export]]

arma::mat Bern_Matrix(arma::mat n, int K, double p){
  
  int N = NonDiagSum(n,K);
  arma::mat P(N,K);
  P = Multiplicity_Matrix(n, K);
  for(int i=0; i<N; ++i){
    for(int j=0; j<K; ++j){
        P(i,j) = P(i,j)*R::rbinom(1,p);
    }
  }  
  return P;
    }

So, if we run the Bern_matrix we have

K = 3
n = matrix(c(2,2,2,0,2,2,0,0,2,0,0,0,0,0,0,0),K+1,K+1,byrow=TRUE)

Bern_Matrix(n, K, 0.5)
       [,1] [,2] [,3]
 [1,]    1    1    1
 [2,]    1    0    0
 [3,]    0    1    1
 [4,]    0    1    0
 [5,]    0    0    1
 [6,]    0    0    1
 [7,]    1    0    0
 [8,]    0    1    0
 [9,]    0    1    0
[10,]    0    1    0
[11,]    1    0    0
[12,]    1    0    0

Where the Bern_Matrix can be regarded as the Multiplicity_Matrix but with the non-zero values, i.e. 1, being changed to 0 with probability 0.5

Conclusion:

What I want to achieve is to make those algorithms even faster if possible, and also to find a way that will make the algorithms scale well as the number of elements inside the cells of the matrix nincreases.

Side Note: I tried to find a way that will avoid the creation of all those matrices and will go directly from the matrix n to the Bern_Matrix without calculating the Vector_Matrix and Multiplicity_Matrix, the code for that is the following

// [[Rcpp::export]]


arma::mat Direct_Matrix(arma::mat n, int K, double p){
  
  int N = NonDiagSum(n,K);
  arma::mat H(N,K);
  H.zeros();
  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){
        H.submat(iter+next,j,iter+next,K-1-i) = as<arma::rowvec>(Rcpp::rbinom( K-j-i, 1, p ));
      }
      next += n(i,j);
    }
  }
  return H;
}

Sorry for the long question, and if something it is needed to be added, I will happily add it and clarify it.

Answer 1

Your Direct_Matrix() function is already doing things quite well; it mainly avoids going multiple times over the same memory. However, there might be some more improvements possible. First of all:

arma::mat H(N,K);
H.zeros();

This zeroes all elements, but we are actually going to overwrite about half of the elements anyway. It might be better to avoid zeroing the matrix, and do this only inside the loop for those elements we know should be zero.

The triple nested for-loop is fine, it looks bad but it just does as many iterations as the number of rows in the result. But inside it we have:

H.submat(iter+next,j,iter+next,K-1-i) = as<arma::rowvec>(Rcpp::rbinom( K-j-i, 1, p ));

Especially if Rcpp::binom() cannot be inlined, this might be a bit inefficient, as it has to pass around a whole temporary vector. I would try to write it as:

for (int k = 0; k < i; ++k)
    H(j, k) = Rccp::rbinom(1, p);
for (int k = i; k < K; ++k)
    H(j, k) = 0;

Of course, now we have much more function calls inside the loop. I'd try to use std::bernoulli_distribution() next instead of Rcpp::rbinom():

std::random_device rd;
std::mt19937 gen(rd());
std::bernoulli_distribution distrib(p);
...
for (int k = 0; k < i; ++k)
    H(j, k) = distrib(gen);
for (int k = i; k < K; ++k)
    H(j, k) = 0;

If p is exactly 0.5 though, you could go for another approach: fill an unsigned integer with a random number, and then you use the individual bits to give you the value for each element in the matrix. Especially if K <= 64, this can be done very efficiently on 64-bits machines. I would then even merge the two loops:

std::random_device rd;
std::mt19937_64 gen(rd());
...
auto bits = gen(); // Generate 64 random bits
bits &= (1UL << i) - 1; // Zero all but the first i bits

for (int k = 0; k < K; ++k) {
    H(j, k) = bits & 1;
    bits >>= 1;
}