Standalone GQ takes as much time as fitting on windows

Hello!

I have a hierarchical model taking approximately 3 minutes to fit.

//
// This Stan program defines a hierarchical model, with a
// vector of values isotopic 'y' modeled as normally distributed
// with hierarchical slopes between C/N and y
// being dependant upon soil C content
//

data {
  // Integers
  int<lower=0,upper=1> prior_only;
  int<lower=0> N; //nb of data points
  int<lower=0> U; //nb of units
  
  // Response
  vector[N] y;
  
  // Covariates
  vector[U] SOC;
  vector[N] CN;
  array[N] int unit;
  
  // Priors
  real lam_sigma_beta_1; real lam_sigma_beta_2;
  real lam_sigma;
  real sd_gamma_0;
  real sd_gamma;
  real mean_beta_0; real sd_beta_0;
}
parameters {
  // Hyperparameters
  cholesky_factor_corr[2] rho; //correlation among hierarchical effects
  vector<lower=0>[2] sigma_beta; //sd of hierarchical effects
  real<multiplier=sd_gamma_0> gamma_0; //int of the slope-SOC relationship
  real<multiplier=sd_gamma> gamma; //slope of the slope-SOC relationship
  // Parameters
  real<offset=mean_beta_0,multiplier=sd_beta_0> beta_0; //int of the isotope-CN relationship
  matrix[U,2] beta_raw; //slope of the isotope-CN relationship
  real<lower=0> sigma; //residual sd
}
transformed parameters{
  // Hyperparameters
  matrix[U,2] mu_beta = append_col(rep_vector(0,U), gamma_0 + SOC * gamma); //linear predictor of slopes
  matrix[2,2] Sigma_beta = diag_post_multiply(rho, sigma_beta); //cov matrix of hierarchical effects
  // Parameters
  matrix[U,2] beta = mu_beta + beta_raw * Sigma_beta; //final hierarchical parameters
  // Linear predictor
  vector[N] mu = beta_0 + beta[unit,1] + CN .* beta[unit,2];
}
model {
  // Hyperpriors
  rho ~ lkj_corr_cholesky(2);
  sigma_beta[1] ~ exponential(lam_sigma_beta_1);
  sigma_beta[2] ~ exponential(lam_sigma_beta_2);
  gamma_0 ~ normal(0, sd_gamma_0);
  gamma ~ normal(0, sd_gamma);
  // Priors
  beta_0 ~ normal(mean_beta_0, sd_beta_0);
  to_vector(beta_raw) ~ std_normal();
  sigma ~ exponential(lam_sigma);
  // Likelihood
  if(prior_only == 0){
    y ~ normal(mu, sigma);
  }
}
generated quantities{
  real y_rep[N] = normal_rng(mu, sigma);
}

When I had ubuntu OS, this very simple standalone generated quantities took some seconds to compute.

//
// This Stan program defines the generated quantities for a hierarchical model, with a
// vector of values isotopic 'y' modeled as normally distributed
// with hierarchical slopes between C/N and y
// being dependant upon soil C content
//

data {
  // Integers
  int<lower=0,upper=1> prior_only;
  int<lower=0> N; //nb of data points
  int<lower=0> N_pred; //nb of C/N values for predictions
  int<lower=0> U; //nb of units
  int<lower=0> U_pred; //nb of SOC values for predictions
  
  // Response
  vector[N] y;
  
  // Covariates
  vector[U] SOC;
  vector[U_pred] SOC_pred;
  vector[N] CN;
  vector[N_pred] CN_pred;
  int unit[N];
  int unit_pred[N_pred];
  
  // Priors
  real lam_sigma_beta_1; real lam_sigma_beta_2;
  real lam_sigma;
  real sd_gamma_0;
  real sd_gamma;
  real mean_beta_0; real sd_beta_0;
}
parameters {
  // Hyperparameters
  cholesky_factor_corr[2] rho; //correlation among hierarchical effects
  vector<lower=0>[2] sigma_beta; //sd of hierarchical effects
  real<multiplier=sd_gamma_0> gamma_0; //int of the slope-SOC relationship
  real<multiplier=sd_gamma> gamma; //slope of the slope-SOC relationship
  // Parameters
  real<offset=mean_beta_0,multiplier=sd_beta_0> beta_0; //int of the isotope-CN relationship
  matrix[U,2] beta_raw; //slope of the isotope-CN relationship
  real<lower=0> sigma; //residual sd
}
transformed parameters{
  // Hyperparameters
  matrix[U,2] mu_beta = append_col(rep_vector(0,U), gamma_0 + SOC * gamma); //linear predictor of slopes
  matrix[2,2] Sigma_beta = diag_post_multiply(rho, sigma_beta); //cov matrix of hierarchical effects
  // Parameters
  matrix[U,2] beta = mu_beta + beta_raw * Sigma_beta; //final hierarchical parameters
  // Linear predictor
  vector[N] mu = beta_0 + beta[unit,1] + CN .* beta[unit,2];
}
generated quantities{
  vector[U_pred] mu_beta_pred = gamma_0 + SOC_pred * gamma;
  vector[N_pred] mu_pred = beta_0 + beta[unit_pred,1] + CN_pred .* beta[unit_pred,2]; 
}

Now that I switched to windows, it takes almost the same time as model fitting! I use as many cores in the GQ calls as in the fitting call.

Does anyone already noticed that? Is there a way to make GQ faster?

I am using cmstanR 0.5.2 with the last cmdstan.

Thank you very much!
All the best,
Lucas

what are the values for U, N, U_pred, and N_pred?

GQ is just doing math on doubles - each iteration will run the calculations in the transformed parameters and generated quantities block. but if you’re manipulating really big arrays in the transformed parameters block, that could be expensive, and if you’re writing out really big vectors in the generated quantities block, that too could be a problem.