Converting simplex arrays to do matrix multiplication

Hi,
In my Stan model I sample two arrays of simplexes, and then do matrix multiplication with them. The way I do it is to copy the arrays into intermediate matrices, and multiply the intermediates. I wondered if there was a better way?

Relevant part of my model:

parameters {
    simplex[S] weights[G];
    simplex[C] signatures[S];
}
transformed parameters {
    matrix<lower=0>[G, C] probs;  // matrix product of weights * signatures
    {
        matrix[G, S] weights_mat;
        matrix[S, C] signatures_mat;
    
        for (g in 1:G) {
            for (s in 1:S) {
                weights_mat[g, s] = weights[g, s];
            }
        }
        
        for (s in 1:S) {
            for (c in 1:C) {
                signatures_mat[s, c] = signatures[s, c];
            }
        }
        probs = weights_mat * signatures_mat;
    }
}
1 Like

Best I could come up with on the spot. I transposed the matrices so that the copies are columnwise:

matrix[S, G] weights_mat_t;
matrix[C, S] signatures_mat_t;

for (g in 1:G) {
  weights_mat_t[:, g] = weights[g];
}

for (s in 1:S) {
  signatures_mat_t[:, s] = signatures[s];
}

probs = (signatures_mat_t * weights_mat_t)';

Nope. We could always write more conversion algorithms to make it a one-liner. Loops that just copy elements are efficient in Stan, so it’s not an efficiency thing, just convenience. The efficiency derives from using matrix ops. So overall, it should be a win.

I would suggest putting it in a function so it doesn’t clutter up your model.

The easiest way to do this is as follows (assuming the matrix doesn’t have zero rows or zero columns):

matrix vector_array_to_matrix(vector[] x) {
  matrix[size(x), rows(x[1])] y;
  for (m in 1:size(x))
    y[m] = x[m]';
  return y;
}
1 Like

Great, thank you. I didn’t realise you could assign directly to matrix columns like that. Would this be more efficient because it aligns better with the way the data is stored in the underlying Eigen objects, i.e. column-major for matrices, row major for arrays?

Thanks again for taking the time to respond!

Thanks for the clarification, Bob. I suspected that I couldn’t avoid copying, but also thought that my way of doing it with nested loops wasn’t the neatest solution. I’ll adopt your suggestion to write it as a function.

Thanks again!

1 Like

Yup, the underlying objects are just Eigen matrices. Only reason I did things with the column major copies was I couldn’t really figure out any other way to make it better :D. Bob’s right though. The copies there won’t matter compared to the matrix multiplies.