Symmetric covariance matrix causes exception in multi_normal_lpdf for being non-symmetric

Modeling
1

Dear Stan Forum,

I encountered a weird issue about a symmetric matrix reported not being symmetric for the covariance matrix of multi_normal_lpdf.

The exact error message is the following:

“Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: Exception: multi_normal_lpdf: Covariance matrix is not symmetric. Covariance matrix[1,17] = 0.000174503, but Covariance matrix[17,1] = 0.000255977 (in ‘/Users/…/model.stan’, line 21, column 4 to column 142) (in ‘/Users/…/model.stan’, line 83, column 4 to line 84, column 86)”

Line 21 mentioned in the message calls the multi_normal_lpdf function, which relies on the covariance matrix QtMOmegaMtQ:

multi_normal_lpdf(Qtrating_benchmarked_i | x*QtMmu_C, QtMOmegaMtQ)

The covariance matrix was defined in line 19:

matrix[num_calibration-d,num_calibration-d] QtMOmegaMtQ = (QtM * QtM') *(x^2 * sigma_epsilon^2+sigma_eta^2);

The above covariance matrix is symmetric because it is of the form “matrix * matrix’ * scalar”.

Lines 83 and 84 mentioned in the message are about the 1-dimensional integration over a function that calls multi_normal_lpdf.

I attached the minimal code for replicating this error message. To run the code, first compile model.stan and then run bash sampling.sh

Any help is greatly appreciated! Thank you very much!

model.stan (4.1 KB)

(Uploading .json data file and .sh script file seems not allowed. You can find them in here: Dropbox - test - Simplify your life)

2

Try a print statement to check the contents of QtMOmegaMtQ. Also wonder if numeric error has creeped in? Maybe do those multiplications as additions on the log scale then exponentiate the result? Alternatively you could loop over the lower tri and assign the partner from the upper tri to enforce symmetric before the multi_normal, but you’d probably want to figure out why things aren’t symmetric before doing that.

1 Like
3

Thanks very much, @mike-lawrence!

As you mentioned, some numerical errors indeed have slipped in: numbers that should be the same theoretically (on the order of exp(-17)) turned out different numerically. For the suggestion of turning multiplications into additions through the log-then-exponentiation procedure, I’m aware of the log_sum_exp function in Stan, which applies to a vector of numbers. Are you aware of a similar Stan function that applies to matrix multiplication, or do I need to do the procedure step by step, i.e., looping over the calculation of each element of the resulting matrix?

4

I don’t know for certain that this will fix the problem (I suspect it will), but you could check if tcrossprod is more stable. See here: 5.5 Dot Products and Specialized Products | Stan Functions Reference

2 Likes
5

tcrossprod made the error message gone! Thank you, @jsocolar!

1 Like