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

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

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 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.

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.