Error in call to cholesky decomposition; require symmetric matrix

[Cody Ross]

All,

I am working on large multilevel model in which I have a variance-covariance matrix, SigmaTheta[k], to be estimated in each cluster; if necessary I can come up with a mimimum working example to post, but basically:
I specify weak priors on Sigma[k] using:
 SigmaTheta[k] <- diag_matrix(sdTheta[k]) * RhoTheta[k] * diag_matrix(sdTheta[k]);

Then I run a cholesky decomposition:
   L[k]<- cholesky_decompose( SigmaTheta[k] );

  Theta[i] <- muStateTheta[Cluster[i]]  + L[Cluster[i]]  * ATheta[i];               
 
I get the following error:

Error in call to cholesky decomposition; require symmetric matrix, but found; x[0,5]=-8.69387:0; x[5,0]=-8.69387:0
error occurred during calling the sampler; sampling not done

On the second run I get:
 x[0,2]=0.0800156:0; x[2,0]=0.0800156:0

These outputs look symmetric to me... Is there a way to adjust the sensitively of the mechanism which checks the parameter declarations?


Running in Rstan,
R version 3.0.0 (2013-04-03)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252  
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                         
[5] LC_TIME=English_United States.1252   

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods 
[8] base    

other attached packages:
[1] rstan_1.3.0   inline_0.3.11 Rcpp_0.10.3 

loaded via a namespace (and not attached):
[1] codetools_0.2-8 tools_3.0.0   
>

[Ben Goodrich]

The threshold can't be changed (unless you dig into the source) but permits a discrepancy of about 10^-8, which if anything is too lenient. Can you post the whole model? There must be some reason why the intermediate calculations are losing too much precision.

Probably the solution is to do

L[k] <- diag_matrix(sdTheta[k]) * cholesky_decompose(RhoTheta[k]);

Ben

[Cody Ross]

Thanks Ben,

I'll give this a try right now, and post the full model if I can't correct the problem.

[Cody Ross]

Ben,

That worked beautifully. But why?

[Ben Goodrich]

Not totally sure why your original way of doing it was causing you to lose so much precision, but evidently it did at least for one proposed set of parameters.

As to why it is the case that the Cholesky factor of a covariance matrix is equal to the product of a diagonal matrix of standard deviations and a Cholesky factor of a correlation matrix, that is just the way it is. So, basically the solution worked by doing a Cholesky decomposition of the correlation matrix (which still had good precision and thus passed the symmetry check) and then pre-multiplying to get the Cholesky factor that would result (at least if there were infinite precision) if you first multiplied the correlation matrix from both sides by the diagonal matrix of standard deviations and then took the Cholesky factor of the result.

Ben