# Accelerating creation of matrices and finding ways for optimal scaling

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)]]

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.

• I am guessing that repeated calls of Rcpp::rbinom are the bottleneck. What is the current performance? and target? Jun 30 at 13:39
• @minem For example, for the matrix 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) when I use the Direct_Matrix which avoids mainy computations it takes on average 8.769453e-09 second, whereas for n = matrix(c(650,800,900,0,500,750,0,0,800,0,0,0,0,0,0,0),K+1,K+1,byrow=TRUE) it takes on average 4.474615e-07 seconds. I would like to push it down as much as possible. Jun 30 at 13:47
• @minem Imagine that I will have to call the Direct_matrix in my MCMC algorithm on each iteration. That is ok when the values on the cells of n are small but when they increase it takes extremely much more time. Jun 30 at 13:49

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

• This is a really helpful answer, I did manage to reduce the computational time with those suggestions!! So, for the H.zeros(); because each row will have zeros on the first cells and last cells I think I would have to use 3 loops, one for the first zeros, for the ones in the middle and for the zeros in the end. What I did is instead of using arma::mat H(N,K)  and H.zeros(); I used NumericMatrix H(N,K) which is much faster and initializes the matrix with zeros. Jul 1 at 11:45
• Also, as you pointed out I changed the .submat which requires the use of a non optimal binomial version with a For Loop that goes exactly to the places where we want to conduct the trial Jul 1 at 11:47
• And lastly, the most significant time reduction was given by std::bernoulli_distribution distrib(p) and I was wondering why?? I mean a random draw from a Bernoulli distribution seems already easy and fast to implement what additionally they do that makes it even faster? Jul 1 at 11:48
• I didn't delve into the internals of Rccp, but could it be that Rccp::rbinom() returns dynamically allocated memory? That means a lot of new and delete behind the scenes. Jul 1 at 19:31