A contradiction result with hand calculation

Modeling
1

Hi all,

I need to define a matrix of the following form:

Cova_trans = diag_matrix(rep_vector(1, R))- matrix_exp(-deltat[k] * diag_matrix(rep_vector(0.1, R))) * diag_matrix(rep_vector(1, R)) * matrix_exp(-deltat[k] * diag_matrix(rep_vector(0.1, R))');
print(Cova_trans);
xi[k] ~ multi_normal(matrix_exp(-deltat[k] * Gamma) * xi[k-1], Cova_trans);

This is a simpler version but just take this one for my question.

By hand calculation, Cova_trans is of the form qI_3 where (0<q<1). (R=3 is provided by data block).

We know that I_3 is positive definite so is qI_3. To be sure, I use R, function is.positive.definite(matrix, tol=1e-1) and obtain the same result of positive definiteness.

However, when I run the code, Stan gives an error:

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite.  last conditional variance is 0.  (in 'model1d881fc227d2_modelALS' at line 170)


Initialization between (-2, 2) failed after 100 attempts. 
 Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
[1] "error occurred during calling the sampler; sampling not done"

I use print function and see that in the console in R I obtain:

[[0.999999,0,0],[0,0.999999,0],[0,0,0.999999]]
[[0.999997,0,0],[0,0.999997,0],[0,0,0.999997]]
[[0.999986,0,0],[0,0.999986,0],[0,0,0.999986]]
[[0.834701,0,0],[0,0.834701,0],[0,0,0.834701]]
[[0.834701,0,0],[0,0.834701,0],[0,0,0.834701]]
[[0.997521,0,0],[0,0.997521,0],[0,0,0.997521]]
[[0.996302,0,0],[0,0.996302,0],[0,0,0.996302]]
[[0.993262,0,0],[0,0.993262,0],[0,0,0.993262]]

Clearly these matrices are of the form qI_3 and they are positive definite.

So my question is: Why does Stan says that the matrix is not positive definite contradicting with hand calculation. Why does Stan gives error even Stan prints positive definite matrices?

Thank you for reading!
P/S: I am sorry that I cannot post the full code. It is not allowed at this moment.

2

Matrices can be numerically indefinite even if they are positive definite by construction, particularly when there are numerically unstable things like subtraction and matrix exponentiation going on. Usually, it can be made to initialize by specifying a narrower range of initial values in the unconstrained space or else by specifying the initial values of things yourself.

2 Likes
3

Thank you so much @bgoodri!

I have divided deltat[k] by a constant and also give suitable initial values and the model is now running.