# Algorithm that creates in a systematic form, matrix with zero-one elements

I'm working on Rcpp and I wrote the following code

#include <RcppArmadilloExtensions/sample.h>
using namespace Rcpp;

// [[Rcpp::export]]

arma::mat SIZA(int K){

arma::mat P;
P.ones(K,K);
P = trimatu(P);

for(int idx=1; idx<K; ++idx){
arma::mat W;
arma::mat Hlp;

W.ones(K-idx,K-idx);
W = trimatu(W);
W = join_horiz(W,Hlp.zeros(K-idx,idx));

P = join_vert(P,W);
}

return P;
}


Essentially it creates matrices of the following form

SIZA(2)
[,1] [,2]
[1,]    1    1
[2,]    0    1
[3,]    1    0

SIZA(3)
[,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

SIZA(4)
[,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    0    1    1    1
[3,]    0    0    1    1
[4,]    0    0    0    1
[5,]    1    1    1    0
[6,]    0    1    1    0
[7,]    0    0    1    0
[8,]    1    1    0    0
[9,]    0    1    0    0
[10,]    1    0    0    0


What I'm curious and interested in is finding a way to make this chunk of code quicker. But from my knowledge, I do not see an obvious way to achieve that. What I'm thinking is to avoid using the FOR LOOP but it seems to be necessary.

Any ideas would be great!

• Please edit your question's title to reflect what the code actually does instead of using a generic title.
– Null
Commented Mar 22, 2021 at 16:16

First of all I would try to avoid all the matrix operations like trimming, joining, etc.

And instead just find a function that tells what value is on each position

bool tell(int x, int y, int K);


This would of course require more mathematical skill then programming, but you might eventually end up doing just this:

arma::mat P;
int maxY = K * (K-1) / 2
for (int y = 0; y < maxY; ++y) {
for (int x = 0; x < K; ++x) {
P.set(x, y, tell(x, y, K) ? 1 : 0)
}
}


If tell(x,y) has O(1) complexity, this will probably be faster.

Another option that I can think of is similar, but more iterative.

Fill the topmost KxK part of the matrix with the pattern. Then another (K-1)x(K-1) part of the metrix, fill with the (K-1) pattern while filling the last column with zeroes. then another part of size (K-2)x(K-2) fill with pattern of size K-2 while filling the last two columns with zeroes. And so on and so on until trying to fill pattern 0x0 where you terminate. But you'll end up on something similar to the previous

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.set(x, Y, tell(x, y) ? 1 : 0)
}
for (int r = step; r < K; ++r) {
P.set(r, Y, 0)
}
}
}


here the tell function actually becomes quite trivial and third parameter is not even needed

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


Whether the tell function for the first implementation is possible in O(1) and how to do it I leave to you and the mathematician inside of you :)

• Thank you a lot for your time!!! So, overall from what I understand you say that it is less time-consuming to use FOR LOOP with a function of complexity O(1) instead of using build-in functions ?? Commented Mar 23, 2021 at 15:23
• Yes, for loop on its own is not inefficient. Actually many of the built-in functions are likely to perform their own loops inside. But where you use the built-in functions to generate matrices and then copy pieces of them to another matrix is where you get inefficient. Because you do quite more operations then KxM which is the size of the resulting matrix. My algorithm works directly on the resulting matrix, no temporary matrices are created and only KxM assignments are performed. Commented Mar 24, 2021 at 5:27

# Improving performance

The for-loop is not the issue here. If you need to iterate over something, a for-loop is the natural thing to use, and most of the time it will be the fastest way.

The main issue I see with this code is that you have temporary variables holding matrices. Those can be expensive. For example, let's look at this line:

W = trimatu(W);


I don't know if the function trimatu() takes the argument by value (requiring a copy) or by reference, but in any case, the return value is a new matrix object that needs to be allocated, and then it is assigned back to W. Depending on the implementation of arma::mat::operator=(), this might be a copy again. So at least one copy is made, possibly three copies are made here.

It would be more efficient to just declare a matrix of the right size up front, and fill that in with ones and zeroes yourself, like so:

arma::mat P(K * (K + 1) / 2, K);

for (int n = K; n > 0; --n) {
// add an n * n triangular matrix to P
for (...) { // rows
for (...) { // columns
P(row, column) = ...;
}
}
}

return P;


Or something of the form slepic suggested. Don't worry about the nested for-loops in this case, you are still only doing as many operations as there are elements in the matrix you are returning.

# Making the code more readable

Try to avoid single-letter variable names, and give them some descriptive names. Exceptions are very commonly used variable names like i and j for loop indices, x, y and z for coordinates, and so on. Another reason for single-letter vairable names might be if you are implementing equations from a paper and want to keep the variable names identical to the symbols used in the formulas, although in that case I would definitely add comments to the code refering to the source of those equations.

To make the structure of the function more readable, I would either go two ways. The first way is not the most efficient way, but it is observing that the result of SIZA(n) contains the result of SIZA(n - 1), so you could rewrite your code to be recursive:

arma::mat SIZA(int K) {
if (K == 0) {
// end of recursion reached
return {};
}

arma::mat P;

// make P be a triangular matrix of size K * K
...

// Append SIZA(n - 1) to the result
P.join_vert(join_horiz(SIZA(K - 1), arma::mat(..., 1, 0)));

return P;
};


The second option is to take an efficient version using for-loops, but split off parts of it into another function:

static void add_triangular_part(arma::mat &P, int row, int K) {
// Add a K * K triangular matrix starting at the given row to P
...
}

arma::mat SIZA(int K) {
arma::mat P(K * (K + 1) / 2, K);

for (int n = K, row = 0; n > 0; row += n--) {
$$$$

• Thank you a lot for your helpful explanation! Because I'm new could you please tell me what do you mean here Don't worry about the nested for-loops in this case, you are still only doing as many operations as there are elements in the matrix you are returning. What I know is that for example a matrix K is of complexity K^2, so it takes two for loops to fill in the matrix, you mean something similar ? Commented Mar 23, 2021 at 15:24
• What I'm saying is that if you are filling all the elements of the matrix, then assuming calculating the value for each element is trivial, then the complexity for a K by K matrix is indeed O(K^2), regardless of whether you use one for-loop or two nested for-loops, or even no visible for-loops (for example by using std::generate()`). Commented Mar 23, 2021 at 16:44