Ordering-invariant latent factor model

I assume your problem is entirely about rotation and not prediction. The easiest idea that comes to my mind is to apply some standard rotation (e.g. varimax) to the draws after-the-fact. My guess is that you wouldn’t be able to do this within Stan on a per-draw basis, but I could be wrong.

Edit 2: I realize now that your use of cholesky_factor_cov was to (much more easily) achieve what I describe below.

An alternative (edit: at least for a 2-component solution) would be to constrain one cell in W to 0 so that it lines up with the same variable across re-orderings, which is a very specific rotation. That should keep the magnitude of the loadings identical. I don’t understand enough about Cholesky matrices to know whether this would work with that data structure, but here is what it would look like with a standard matrix.

data{
  int i; // which cell should be constrained to 0
  int j[K-1]; // which cells should be estimated
}
parameters{
  vector[K-1] W1_raw;
  matrix[K, D-1] W_raw;
}
transformed parameters{
  matrix[K, D] W;
  W[i,1] = 0;
  W[j,1] = W1_raw;
  W[,2:D] = W_raw;
}

Maybe you could achieve the same thing with a very narrow prior on that one cell, thereby avoiding all of this reconstruction. But I’m not sure.