Positive definite warnings

[ScottClaessens]

@ErikRingen We have discussed this before, but often the model throws (many) positive definite warnings at the beginning of sampling, as shown below. My understanding is that these warnings are harmless. But users may be wonder about the meaning of these warnings and whether they should be concerned.

library(coevolve)

coev_fit(
  data = authority$data,
  variables = list(
    political_authority = "ordered_logistic",
    religious_authority = "ordered_logistic"
  ),
  id = "language",
  tree = authority$phylogeny,
  dist_mat = authority$distance_matrix,
  # additional arguments for cmdstanr::sample()
  chains = 4,
  parallel_chains = 4,
  seed = 1
)

Note: Distance matrix detected. Gaussian processes over spatial distances have been included for each variable in the model using the 'exp_quad' covariance kernel.
Compiling Stan program...
Running MCMC with 4 parallel chains...

Chain 1 Rejecting initial value:
Chain 1   Error evaluating the log probability at the initial value.
Chain 1 Exception: cholesky_decompose: Matrix m is not positive definite (in 'C:/Users/scla896/AppData/Local/Temp/RtmpeQ4qpu/model-f40261f757c.stan', line 183, column 4 to column 46)
Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 2 Rejecting initial value:
Chain 2   Error evaluating the log probability at the initial value.
Chain 2 Exception: cholesky_decompose: Matrix m is not positive definite (in 'C:/Users/scla896/AppData/Local/Temp/RtmpeQ4qpu/model-f40261f757c.stan', line 183, column 4 to column 46)
Chain 2 Rejecting initial value:
Chain 2   Error evaluating the log probability at the initial value.
Chain 2 Exception: cholesky_decompose: Matrix m is not positive definite (in 'C:/Users/scla896/AppData/Local/Temp/RtmpeQ4qpu/model-f40261f757c.stan', line 183, column 4 to column 46)
Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 4 Iteration:  100 / 2000 [  5%]  (Warmup)
...

The warnings point to this part of the Stan code, suggesting that it is specific to the distance Gaussian Processes:

// distance covariance functions
  for (j in 1:J) {
    matrix[N_tips,N_tips] dist_cov;
    matrix[N_tips,N_tips] L_dist_cov;
    for ( i in 1:(N_tips-1) )
      for ( m in (i+1):N_tips ) {
        dist_cov[i,m] = sigma_dist[j] * exp(-(square(dist_mat[i,m]) / (2.0 * square(rho_dist[j]))));
        dist_cov[m,i] = dist_cov[i,m];
      }
    for ( q in 1:N_tips )
      dist_cov[q,q] = sigma_dist[j] + 0.01;
    L_dist_cov = cholesky_decompose(dist_cov);
    dist_v[,j] = L_dist_cov * dist_z[,j];
  }

We could suppress these warnings. But is there any way we can improve the Stan code here? For example, by declaring the matrices as cov_matrix and cholesky_factor_cov instead?

[ErikRingen]

Good question! Maybe the alternative parameterizations/priors of the GP using lengthscale is better in this sense? This is like the "McElreath" parameterization of GPs but it is not what most software uses (including the built-in Stan functions).

[ScottClaessens]

I agree we could probably reparametrize. I'll take a look at this soon.