cmdstanr loo_moment_match fails for models that have a cholesky_factor

Describe the bug Running loo(moment_match = TRUE) on a fit object that has a parameter of type cholesky_factor_cov[K] yields the error "cholesky_factor_free: y is not lower triangular; y[1,3]=x (found before start of program)" where x may vary from one run to the next; in one case it was 6.65e271, in another case it was 1.26e-312.

To Reproduce

create_model <- function() {
  code <- '
data {
  int<lower=1> K;
  int<lower=1> N;
  array[N] vector[K] X;
  real<lower = K - 1 + 1e-3> nu;
	cov_matrix[K] scale_matrix;
}
transformed data {
  cholesky_factor_cov[K] L_scale_matrix = cholesky_decompose(scale_matrix);
  array[N] vector[K] zero = rep_array(rep_vector(0.0, K), N);
}
parameters {
	cholesky_factor_cov[K] LX;
}
model {
  X ~ multi_normal_cholesky(zero, LX);
	LX ~ inv_wishart_cholesky(nu, L_scale_matrix);
}
generated quantities {
  vector[N] log_lik;
  real lp;
  for (i in 1:N)
    log_lik[i] = multi_normal_cholesky_lpdf(X[i] | rep_vector(0, K), LX);
  lp = inv_wishart_cholesky_lpdf(LX | nu, L_scale_matrix) +  // prior
       sum(log(diagonal(LX))) +                              // Jacobian
       sum(log_lik);
}
  '
  cmdstanr::cmdstan_model(cmdstanr::write_stan_file(code))
}

create_data <- function() {
  N <- 2000
  K <- 5
  if (is.null(df))
    df <- K
  Sigma <- diag(rep(1+1/K, K)) - 1/K
  LX <- t(chol(Sigma))
  scale_matrix <- diag(rep(1, K))
  nu <- K
  set.seed(23478)
  X <- lapply(1:N, function(i){ as.vector(mvtnorm::rmvnorm(1, sigma=Sigma)) })
  N1 <- N %/% 100
  mu <- 1.0
  X[1:N1] <- lapply(1:N1, function(i){ mu + as.vector(mvtnorm::rmvnorm(1, sigma=16*Sigma))})
  list(
    K = K,
    N = N,
    X = X,
    nu = nu,
    scale_matrix = scale_matrix
  )
}

dummy_model <- create_model()
dummy_data <- create_data()
dummy_fit <- dummy_model$sample(data = dummy_data)
dummy_loo <- dummy_fit$loo(moment_match = TRUE)

Expected behavior I expect the loo computation to proceed without error.

Operating system MacOS Ventura 13.4

CmdStanR version number 0.6.0

Additional context I have seen problems like this with other models that had cholesky_factor_cov parameters. This is just the simplest example I could create.

Jul 31 '23 18:07 ksvanhorn

What happens if you run:

dummy_fit <- dummy_model$sample(data = dummy_data, sig_figs = 12)
dummy_loo <- dummy_fit$loo(moment_match = TRUE)

Jul 31 '23 18:07 rok-cesnovar

I get this:

Error: Exception: cholesky_factor_free: y is not lower triangular; y[1,3]=7.87685e-309 (found before start of program)

Jul 31 '23 20:07 ksvanhorn

So I guess sig_figs = 18 also errors, right?

Jul 31 '23 20:07 rok-cesnovar

That's right, I just get a different number (y[1,3]=6.89908e-310). And I have verified that all of the posterior draws for LX[i,j], where j > i, are exactly 0.

Jul 31 '23 22:07 ksvanhorn

I just tested this three times without increasing sig_figs and I don't get any errors. @ksvanhorn can you retry with the latest CmdStan?

Nov 29 '23 19:11 avehtari

loo_moment_match fails for models that have a cholesky_factor_cov parameter.