Multivariate Data with many DVs - Covariance Matrix is not positive definite

[Matt]

Hi,
I'm not sure where this should go, so I'll stick it here but let me know if you want it somewhere else.

I was running into problems with a more complex model when I started adding in more dependent variables. So I simplified it down to try and isolate where the problem might be. It seems as though stan has a limit of approximately 15 dependent variables for a multi-normal IV node, and at 15 DVs it starts to struggle. I get the following error for the seed set at 666 and 15 DVs:

Error : Error in function stan::prob::multi_normal_log(d): Covariance matrix is not positive definite. Covariance matrix(0,0) is 2.4779795485011187.
error occurred during calling the sampler; sampling not done

And sometimes it starts eating massive amounts of memory at around 20 DVs (compared to JAGS which handles a similar specification with much less memory) - I think it also doesn't release the memory from the run when it doesn't complete, but this may be an RStudio - Rstan interface issue not related to stan. Possibly also related to this last post here?

Example code and system specs below. Let me know if you need more info.

albMRegMod = '
  data {
    int<lower=1> nData;
    int<lower=1> nVars;
    real x3[nData];
    matrix[nData, nVars] y;
    matrix[nVars, nVars] Imat;
    vector[nVars] b3p;
  }
  parameters {
    vector[nVars] b0;
    vector[nVars] b3;
    cov_matrix[nVars] Omega3;
    real<lower=0, upper=100> sd[nVars];
  }
  transformed parameters {
    vector[nVars] mu[nData];
    for ( np in 1:nVars){
      for ( i in 1:nData){
        mu[i, np] <- b0[np] + b3[np] * x3[i];
      }
    }
  }

  model {
//Priors
    b0 ~ normal( 0 , 10);
    b3 ~ multi_normal(b3p , Omega3);
    Omega3 ~ inv_wishart(nVars+1, Imat);

//Likelihood
    for(np in 1:nVars){
      for(i in 1:nData){
        y[i, np] ~ normal(mu[i, np], sd[np]);
      }
    }
  }
'
seed = 666
nVars = 15 #12 works, 13:14 give lots of warning messages, 15+ gives errors
nData = 300
b0p = rep(0, nVars)
b3p = rep(0, nVars)
Imat <- diag(1, nVars)
y <- matrix(rnorm(nVars*nData), nrow=nData, ncol=nVars)
x3 <- rnorm(nData)
dataList = list( y = y, x3 = x3,
                 Imat = Imat, b0p = b0p, b3p = b3p,
                 nVars = nVars, nData = nData)

f1 <- stan(model_code = albMRegMod, data = dataList, chains = 1, iter=500,
           seed = seed, chain_id = 1,
           pars = c("b0", "b3"))
 

R version 2.15.1 running in RStudio
Stan version 1.1.1 (downloaded today)
Updated other R packages today also
OS Linux 64 bit openSuse KDE
8 Gb RAM, 4-core Intel, hyperthreaded to 8-cores
Intel onboard graphics

[Ben Goodrich]

On Wednesday, February 6, 2013 4:01:17 AM UTC-5, Matt wrote:

I was running into problems with a more complex model when I started adding in more dependent variables. So I simplified it down to try and isolate where the problem might be. It seems as though stan has a limit of approximately 15 dependent variables for a multi-normal IV node, and at 15 DVs it starts to struggle.

If a parameter is declared as a cov_matrix, then it is positive definite with probability one. But numerically it can fail the check for positive definiteness nonetheless.
 

I get the following error for the seed set at 666 and 15 DVs:

Error : Error in function stan::prob::multi_normal_log(d): Covariance matrix is not positive definite. Covariance matrix(0,0) is 2.4779795485011187.
error occurred during calling the sampler; sampling not done

That shouldn't be a fatal exception; it should just reject the proposal. Marcus did a bunch of stuff with multi_normal for v1.1.1 so maybe the code is overreacting to semi-definiteness now?
 

And sometimes it starts eating massive amounts of memory at around 20 DVs (compared to JAGS which handles a similar specification with much less memory) - I think it also doesn't release the memory from the run when it doesn't complete, but this may be an RStudio - Rstan interface issue not related to stan. Possibly also related to this last post here?

I think JAGS uses a parameterization of the multivariate normal with the precision matrix, in which case you don't have to invert anything and thereby save RAM. As of v1.1.1, Stan has that too. It is called multi_normal_prec() I think. You would just call that on your cov_matrix; just remember in the output that if you really want to analyze a covariance, you need to invert that matrix.

Ultimately, the problem is driven by the fact that this prior

Omega3 ~ inv_wishart(nVars+1, Imat);

puts too much mass in the corners of the parameter space where the matrix is barely positive definite. Increasing the degrees of freedom would help, as would using a lkj_corr prior on the correlation matrix with shape parameter > 1 and some other prior for standard deviations.

Or see the chapter called Optimizing Stan Code, which has reparameterizations of the multivariate normal and inverse Wishart that could avoid some numerical problems in addition to being faster.

Ben

[Matt]

Thanks Ben. The re-parameterisation with multi_normal_cholesky works well from the optimizing section, but I haven't yet figured out how to successfully use the inverse Wishart re-parameterisation (the last attempt exhausted my RAM). Will keep tweaking.

Any plans for a multi_student_t_cholesky?

[Ben Goodrich]

On Thursday, February 7, 2013 4:42:37 AM UTC-5, Matt wrote:

Thanks Ben. The re-parameterisation with multi_normal_cholesky works well from the optimizing section, but I haven't yet figured out how to successfully use the inverse Wishart re-parameterisation (the last attempt exhausted my RAM). Will keep tweaking.

Reparameterizing to avoid calling inv_wishart() --- along the lines of page 106 in the manual --- should use way less RAM. Maybe something went wrong?

Any plans for a multi_student_t_cholesky?

Maybe Marcus is planning something. The log density would be easy enough to implement in Stan code. However, generating a Cholesky factor of a correlation matrix without first generating the correlation matrix is tedious in Stan code, although there is C++ code in transform.hpp that you could try to translate.

Ben

[Marcus Brubaker]

With the new varis that we did last week, there is a lot to be done to make wishart, inv_wishart and multi_student_t much faster.  I just did some small things with multi_student_t now that I will push once I confirm the unit tests are OK but there's more to be done.

Cheers,

Marcus

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.