Cholesky_decompose error


I’m very new to Stan - using MATLABStan. I keep getting the following error:

output-1.csv: Rejecting initial value:
output-1.csv: Error evaluating the log probability at the initial value.
output-1.csv: Exception: cholesky_decompose: Matrix m is not positive definite (in ‘C:/Users/rmw61/Documents/GES/GUCalibration/bc_KOH_1x.stan’ at line 126)

Is it possible to view the covariance matrix so I can see what is causing the problem? I’ve tried setting init = 0 and have tried changing the parameter specifications but I keep getting the same error. I’ve had the same code working on a different problem which makes me think it is something to do with the parameter specifications for the new data.

You would have to use print(...) statements.

If your matrix is a free parameter, you should make sure it is positive definite – random initial proposals may not be. If init=0 makes all parameters zero, including the matrix you want to cholesky-decompose you will probably get this error since it is not positive definite (maybe try to set its initial value to the identity matrix instead).
If that doesn’t solve it or help, maybe some of the code would help.

Hi, thanks for responding. I’m not quite sure which bit of the code is most useful? Here is the Stan file -

data {
int<lower=0> n; // number of field data
int<lower=0> m; // number of computer simulation
int<lower=0> p; // number of observable inputs x
int<lower=0> q; // number of calibration parameters t
vector[n] y; // field observations
vector[m] eta; // output of computer simulations
vector[n] xf; // observable inputs corresponding to y
// (xc, tc): design points corresponding to eta
vector[m] xc;
matrix[m, q] tc;
// inputs for prior distributions !!NEW!!
row_vector[q] mu_theta; // mean values for theta priors
row_vector[q] sd_theta; // standard deviations for theta priors
real shape_eta; // shape parameter for prior of emulator
real rate_eta; // rate parameter for prior of emulator
real shape_b; // shape parameter for prior of model bias
real rate_b; // rate parameter for prior of model bias
real shape_e; // shape parameter for prior of random error
real rate_e; // rate parameter for prior of random error
row_vector[p+q] a_eta; // a values for beta_eta priors
row_vector[p+q] b_eta; // b values for beta_eta priors
real a_bias; // a value for beta_b prior
real b_bias; // b value for beta_b prior

transformed data {
int<lower=1> N;
vector[n+m] z; // z = [y, eta]
vector[n+m] mu; // mean vector

N = n + m;
// set mean vector to zero
for (i in 1:N) {
mu[i] = 0;
z = append_row(y, eta);

parameters {
// theta: calibration parameters
// rho_eta: reparameterization of beta_eta
// rho_b: reparameterization of beta_b
// lambda_eta: precision parameter for eta
// lambda_b: precision parameter for bias term
// lambda_e: precision parameter of observation error

row_vector<lower=0,upper=1>[p+q] rho_eta;
real<lower=0,upper=1> rho_b;
real<lower=0> lambda_eta;
real<lower=0> lambda_e;
real<lower=0> lambda_b;
row_vector<lower=0,upper=1>[q] theta;


transformed parameters {
// beta_delta: correlation parameter for bias term
// beta_e: correlation parameter of observation error
row_vector[p+q] beta_eta;
real beta_delta;
beta_eta = -4.0 * log2(rho_eta); // !!NEW!! log2 instead of log
beta_delta = -4.0 * log2(rho_b); // !!NEW!! log2 instead of log

model {
// declare variables
matrix[N, (p+q)] xt;
matrix[N, N] sigma_eta; // simulator covariance
matrix[n, n] sigma_delta; // bias term covariance
matrix[N, N] sigma_z; // covariance matrix
matrix[N, N] L; // cholesky decomposition of covariance matrix
real temp_delta;
row_vector[p+q] temp_eta;

// xt = [[xt,theta],[xc,tc]]
xt[1:n, 1] = xf;
xt[(n+1):N, 1] = xc;
xt[1:n, (p+1):(p+q)] = rep_matrix(theta, n);
xt[(n+1):N, (p+1):(p+q)] = tc;

// diagonal elements of sigma_eta
sigma_eta = diag_matrix(rep_vector((1 / lambda_eta), N));

// off-diagonal elements of sigma_eta
for (i in 1:(N-1)) {
for (j in (i+1):N) {
temp_eta = xt[i] - xt[j];
sigma_eta[i, j] = beta_eta .* temp_eta * temp_eta’;
sigma_eta[i, j] = exp(-sigma_eta[i, j]) / lambda_eta;
sigma_eta[j, i] = sigma_eta[i, j];

// diagonal elements of sigma_delta
sigma_delta = diag_matrix(rep_vector((1 / lambda_b), n));

// off-diagonal elements of sigma_delta
for (i in 1:(n-1)) {
for (j in (i+1):n) {
temp_delta = xf[i] - xf[j];
sigma_delta[i, j] = beta_delta* temp_delta * temp_delta’;
sigma_delta[i, j] = exp(-sigma_delta[i, j]) / lambda_b;
sigma_delta[j, i] = sigma_delta[i, j];

// computation of covariance matrix sigma_z
sigma_z = sigma_eta;
sigma_z[1:n, 1:n] = sigma_eta[1:n, 1:n] + sigma_delta;

// add observation errors
for (i in 1:n) {
sigma_z[i, i] = sigma_z[i, i] + (1.0 / lambda_e);

// Specify distributions for thetas and hyperparameters !!NEW!!

theta[1:q] ~ normal(mu_theta, sd_theta);
rho_eta[1:(p+q)] ~ beta(a_eta, b_eta);
rho_b ~ beta(a_bias, b_bias);
lambda_b ~ gamma(shape_b, rate_b);
lambda_eta ~ gamma(shape_eta, rate_eta);
lambda_e ~ gamma(shape_e, rate_e);

L = cholesky_decompose(sigma_z); // cholesky decomposition
z ~ multi_normal_cholesky(mu, L);

This is called in MATLAB by:

data = struct (‘xf’,xxf,‘xc’, xxc,‘x’,xx,‘y’,yyf,‘eta’,yyc,‘tc’,tt,… % normalized data
‘n’,n,‘m’,m,‘p’,p,‘q’,q,… % dimensions of data
‘mu_theta’,mu_theta,‘sd_theta’,sd_theta,… % parameters for priors calibration parameters
‘shape_eta’, shape_eta,‘rate_eta’,rate_eta,… % parameters for pripor hyperparameters emulator
‘shape_b’,shape_b,‘rate_b’,rate_b,… % parameters for pripor hyperparameters model bias
‘shape_e’,shape_e,‘rate_e’,rate_e); % parameters for pripor hyperparameters random error

calibration = stan (‘file’,name_stan,‘data’,data,‘iter’,iter,‘chains’,chains,‘verbose’,true,‘refresh’,(iter/100),…

Any advice gratefully received, thankyou :)

The problem may be that from random initial parameters used to set up sigma_eta you get a non-positive definite matrix often, which will cause it to fail to initialize after a number of tries.
If you expect the sigma_eta to be positive definite and you have matlab code where you can verify it is for predefined initial values, you can try to looping over the elements of the matrix in Stan to see if you get what you expected. Its twice the work to implement the model in your favorite language and in Stan, but sometimes it’s easier than debugging through printing in Stan directly, at least for me it is.
Hope this helps.

One last hint, if you enclose your code bbetween \texttt{```Stan}, and \texttt{```}, it’ll show syntax highlight and make it more readable.

OK, great thanks for the tips I’ll give that a try. And thanks for the tip about the code - I’ll make sure to do that next time.

1 Like