I am attempting to fit a stan model which includes a product term of 2 sampled parameters. In my original formulation of the model I was able to produce a product and estimate a linear effect without sampling that product. I decided to re-formulate the model in matrix form, however this has led to the product term being treated as a parameter and is now sampled. This has resulted in sub-optimal performance regarding estimates as well as computation time. I assumed this is a problem on my end with how I am specifying the product.

This is an excerpt from the transformed parameters block. This is the model in the non-matrix formulation :

```
vector [N] mueta;
matrix [N ,2] xi ;
xi = zi * diag_pre_multiply(sigmaxi,L1)';
mueta = b1[1] +
b1[2]*xi[,1] +
b1[3]*xi[,2] +
b1[4]*(xi[,2].*xi[,1]) ;
```

This is my matrix formulation:

```
vector [N] mueta;
matrix [N ,4] xi ;
xi[,1] = rep_vector(1, N);
xi[,2:3] = zi * diag_pre_multiply(sigmaxi,L1)' ; // Cholesky decomp produces factor scores for xi.
xi[,4] = xi[,2].*xi[,3] ;
mueta = xi * b1 ;
```

Where xi[,1] is an N length vector of 1s (intercept), xi[,2:3] are factor scores, and xi[,4] is a product (interaction) of xi[,2] and xi[,3].

I would like to estimate all 4 effects (b1:b4) in the matrix formulation without sampling xi[,4]. Is there a way to specify that xi[,4] is the product of xi[,2] and xi[,3] and should not be sampled? Maybe this is not possible? If so, the matrix formulation (which should be more computationally efficient) is not realistic to specify for this type of model.