How to select rows from two matrices and multiply them

Question

I am trying to loop through the rows in two matrices and multiply them but each attempt produces the first error below, which is clear, and a second error on dimension mismatch from the last two attempts. I have tried to change the format of the extracted rows to a form that is acceptable to stan but I don't know how to coerce them to the relevant format. How do I select the rows and multiply them please?

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
No matches for:
row_vector * row_vector

Expression is ill formed.

---

16:       for(i in 1:Ns) {
17: //        yh[i, n] = xnew[n]* beta[i];
18:         yh[i, n] = xnew[n]* beta[i, ];
                                        ^
19: //        yh[i, n] = xnew[n]* row(beta, i);

---

My rstan code

library("rstan") 
txt <- 
  'data {
    int<lower=0> N;   
    int<lower=0> K;     
    int Ns;
    matrix[N, K] xnew;   
    matrix[Ns, K] beta;
  }
  parameters {
  }
  model {
  }
  generated quantities {
    matrix[Ns, N] yh;

    for(n in 1:N) {
      for(i in 1:Ns) {
//        yh[i, n] = xnew[n]* beta[i];
          yh[i, n] = xnew[n]* beta[i, ];
//        yh[i, n] = xnew[n]* row(beta, i);
//        yh[i, n] = xnew[n,]* row(beta, i);
//        yh[i, n] = to_vector(row(xnew, n))* to_matrix(row(beta, i));
//        yh[i, n] = to_vector(row(xnew, n))* row(beta, i);         
      }  
    }
  }

'
stan_model(model_code=txt) 

For clarity this is what i am trying to do in base R

set.seed(1)
Ns=10; N=2; K=3
beta = matrix(rnorm(Ns*K), ncol=K)
xnew = matrix(rnorm(N*K), ncol=K)

yh=matrix(nr=Ns, nc=N)
for(n in 1:N) {
  for(i in 1:Ns) {
    p = as.numeric(xnew[n, , drop=FALSE] %*% beta[i,])
    yh[i, n] = p
  }  
}

#tcrossprod(beta, xnew)

EDIT:

This seems to do the trick:

yh[i, n] = dot_product(row(xnew, n), row(beta, i));

But is there a way to calculate this without looping through each row? (I can't see anything at https://mc-stan.org/docs/2_18/functions-reference/dot-products-and-specialized-products.html)

Answer 1

Well I made a song and dance over this ... but you can just use matrix multiplication: yh = beta* xnew'; EDIT: following Bob's advice in the comment below I moved the matrix transpose to the transformed data block.

So the full code:

txt <- 
  "data {
    int<lower=0> N;   
    int<lower=0> K;     
    int Ns;
    matrix[N, K] xnew;   
    matrix[Ns, K] beta;
  }
  transformed data{
    matrix[K, N] xnew_t = xnew';
  }
  parameters {
  }
  model {
  }
  generated quantities {
    matrix[Ns, N] yh;
    yh = beta* xnew_t;
  }

"

fit <- stan(model_code=txt,  data = list(beta=beta, xnew=xnew, Ns=10, N=2, K=3 ),    
            chains = 1, seed = 1, iter=1, algorithm = "Fixed_param")

ex_samp = as.matrix(fit)
all.equal(yh, matrix(ex_samp[-length(ex_samp)], nc=2))

The relevant docs: https://mc-stan.org/docs/2_18/functions-reference/matrix-arithmetic-operators.html