R Session Aborted when Replicating a Stan Factor Analysis model


#1

Update: the issue was solved if I request the Stan package and compile the model without running the C++ toolchain configuration line

Sys.setenv(LOCAL_CPPFLAGS = '-march=native')

Not sure why it is the case. When would it be appropriate to do the configuration though?

----------------- Original post below-----------
Hi,

I was trying to replicate a factor analysis model done by Rick Farouni with mostly the original code, but the RStudio session would abort, saying R encountered a fatal error and the session was terminated, as soon as the console started saying “SAMPLING FOR MODEL ‘fa’ NOW (Chain 1).”

My RStan version is 2.18.2. RStudio version 1.1.463. R version 3.5.2. Windows PC 8GB RAM 2.8Hz i5 CPU. I tried deleting the .rds file before running the model, but that did not help.

My Makevars file has

CXXFLAGS=-O3 -Wno-unused-variable -Wno-unused-function

My code is below. Could it be simply an issue of running out of RAM? I wonder if the optimization of the model itself, as previously discussed here and on the Google group, has anything to do with this at all. Any suggestion on how to work around it would be much appreciated!

library("rstan")
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
Sys.setenv(LOCAL_CPPFLAGS = '-march=native')

#simulate the data
library("MASS")
set.seed(42)
D <-3
P <- 10 
N <-300

mu_theta <-rep(0,D) # the mean of eta
mu_epsilon<-rep(0,P) # the mean of epsilon
Phi<-diag(rep(1,D))
Psi <- diag(c(0.2079, 0.19, 0.1525, 0.20, 0.36, 0.1875, 0.1875, 1.00, 0.27, 0.27))
l1 <- c(0.99, 0.00, 0.25, 0.00, 0.80, 0.00, 0.50, 0.00, 0.00, 0.00)
l2 <- c(0.00, 0.90, 0.25, 0.40, 0.00, 0.50, 0.00, 0.00, -0.30, -0.30)
l3<-  c(0.00, 0.00, 0.85, 0.80, 0.00, 0.75, 0.75, 0.00, 0.80, 0.80)
L <-cbind(l1,l2,l3) # the loading matrix

Theta <-mvrnorm(N, mu_theta, Phi) # sample factor scores
Epsilon <-mvrnorm(N, mu_epsilon, Psi) # sample error vector
Y<-Theta%*%t(L)+Epsilon# generate observable data

stan model

data {
  int<lower=1> N;                // sample size
  int<lower=1> P;                // number of items
  matrix[N,P] Y;                 // data matrix of order [N,P]
  int<lower=1> D;              // number of latent dimensions 
}
transformed data {
  int<lower=1> M;
  vector[P] mu;
  M  <- D*(P-D)+ D*(D-1)/2;  // number of non-zero loadings
  mu <- rep_vector(0.0,P);
}
parameters {    
  vector[M] L_t;   // lower diagonal elements of L
  vector<lower=0>[D] L_d;   // lower diagonal elements of L
  vector<lower=0>[P] psi;         // vector of variances
  real<lower=0>   mu_psi;
  real<lower=0>  sigma_psi;
  real   mu_lt;
  real<lower=0>  sigma_lt;
}
transformed parameters{
  cholesky_factor_cov[P,D] L;  //lower triangular factor loadings Matrix 
  cov_matrix[P] Q;   //Covariance mat
{
  int idx1;
  int idx2;
  real zero;
  idx1 <- 0;
  idx2 <- 0; 
  zero <- 0;
  for(i in 1:P){
for(j in (i+1):D){
  idx1 <- idx1 + 1;
  L[i,j] <- zero; //constrain the upper triangular elements to zero 
}
  }
  for (j in 1:D) {
  L[j,j] <- L_d[j];
for (i in (j+1):P) {
  idx2 <- idx2 + 1;
  L[i,j] <- L_t[idx2];
} 
  }
} 
Q<-L*L'+diag_matrix(psi); 
}
model {
// the hyperpriors 
   mu_psi ~ cauchy(0, 1);
   sigma_psi ~ cauchy(0,1);
   mu_lt ~ cauchy(0, 1);
   sigma_lt ~ cauchy(0,1);
// the priors 
  L_d ~ cauchy(0,3);
  L_t ~ cauchy(mu_lt,sigma_lt);
  psi ~ cauchy(mu_psi,sigma_psi);
//The likelihood
for( j in 1:N)
Y[j] ~ multi_normal(mu,Q); 
}

#fit model to the data

fa.data <-list(P=P,N=N,Y=Y,D=D)
# a function to generate intial values that are slightly jittered for each chain.
init_fun = function() {
  init.values<-list(L_t=rep(0,24)+runif(1,-.1,.1),
                    L_d=rep(.5,D)+runif(1,-.1,.1),
                    psi=rep(.2,P)+runif(1,-.1,.1),
                    sigma_psi=0.15+runif(1,-.1,.1),
                    mu_psi=0.2++runif(1,-.1,.1),
                    sigma_lt=0.5+runif(1,-.1,.1),
                    mu_lt=0.0+runif(1,-.1,.1))
  return(init.values); 
} 

#compile the model
library("parallel")
Nchains <- 1
Niter <- 300
fa.model<- stan("fa.stan", 
                data = fa.data,
                chains = 1, iter = Niter, init = init_fun, seed = 42,
                pars=c("L","psi","sigma_psi","mu_psi","sigma_lt","mu_lt"))
#this would crash on my laptop during poterior sampling

#2

Glad you got that worked out – kinda weird error haha.

The -march=native flag enables some processor-specific optimizations. Maybe somehow the compiler tried to target a different CPU than the one you actually have and so things are blowing up when you try to run the sampler. I’m not sure.

You’re not probably missing too much by skipping it. -O3 is probly good enough.


#3

Thanks for your input. I’m new to Stan but would love to better appreciate its optimization capability!