Is it possible to disable symmetry checks?

Is it possible to turn off the symmetry check?
Commonly, BLAS calls only look at the lower or upper triangle of the input matrix, treating the other triangle as the mirror image.

I ask because of errors like:

Chain 11: Rejecting initial value:
Chain 11:   Error evaluating the log probability at the initial value.
Chain 11: Exception: Exception: cholesky_decompose: m is not symmetric. m[1,2] = -4.37401, but m[2,1] = -4.37401  (in 'modela8982e43d111_model6_marginalize_interp' at line 178)

For numerical reasons, the matrix will often not be precisely symmetric, although within rounding, as the exception itself makes clear: obviously -4.37401 == -4.37401.

It’s also worth pointing out that I am calculating the symmetric matrix via:

non_missing_block - quad_form(inverse_spd(missing_block), association_block)

instead of

non_missing_block - quad_form_sym(inverse_spd(missing_block), association_block)

because quad_form_sym would complain of similar numerical differences, eg

Exception: Exception: quad_form_sym: A is not symmetric. A[72,78] = 1.87279e+07, but A[78,72] = 1.87279e+07

Floating point errors mean that matrices that ought to be precisely symmetric will often fail the symmetry checks.
For example:

SAMPLING FOR MODEL 'model6_marginalize_interp' NOW (CHAIN 13).
Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[37,68] = 1.03544e+07, but A[68,37] = 1.03544e+07  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[1,67] = -13273.6, but A[67,1] = -13273.6  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[8,68] = 275203, but A[68,8] = 275203  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[1,68] = -97994.3, but A[68,1] = -97994.3  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[1,73] = -54295.9, but A[73,1] = -54295.9  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[56,74] = -4.43735e+07, but A[74,56] = -4.43735e+07  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13 failed on all 100 initial values, causing the chain to abort. Every single one of the error messages looked something like A[56,74] = -4.43735e+07, but A[74,56] = -4.43735e+07 – clearly, the error is smaller than the number of significant figures the print method keeps.

I can of course replace quad_form_sym with quad_form. This seemingly solves the problem.
After another operation, I do a cholesky_decompose, which does sometimes throw an exception, but rarely enough that I’ve not seen it cause a chain to crash.
That is the approach I have generally been forced to take. An alternative is trying to force symmetry in the input.

Is there a way to disable the symmetry check? Or, at least decrease its sensitivity.

EDIT:
For what it is worth, this was the offending bit of code:

quad_form_sym(inverse_spd(missing_block), association_block)

I replaced that with this, which I assume is more stable than just using quad_form, and ought to guarantee symmetry of the output:

crossprod(mdivide_left_tri_low(cholesky_decompose(missing_block), association_block))

which does not seem to throw any errors. I am guessing the inverse_spd produced large magnitudes, which may have tripped an absolute tolerance? Or large magnitudes reduce stability?
Either way – code works, but it would be nice to have more functions with an api like mdivide_left_tri_low, that cast triangular (or symmetric) views by focusing on only a lower/upper half of a matrix.

EDIT:
Runtime was/is a major issue with that model.
So I was shocked to see that switching to the crossprod(mdivide_left_tri_low(cholesky(...))) version decreased the 1000 transition-10-leapfrog time roughly 3-fold (~ 800 seconds to ~250 seconds)!
I am pleased.
Seems like trying multiple variations can be worthwhile. And that I should probably favor “div” methods, even over specialized inverses (like inverse_spd) in general?
EDIT:
Performance difference is closer to 5x based on sampling rate.
Maybe I should be glad I couldn’t disable a symmetry check, otherwise I probably wouldn’t have tried this work-around.

3 Likes

If you want to know why the difference occur it may be that the more numerically stable version is more accurate and results in a higher acceptance rate for a trajectory which in turn allows adaptation to push for a higher stepsize in the integrator. Pretty easy to look for this in diagnostics and it’s generally relevant for speeding up sampling

@chriselrod, Thank you so much for posting the solution to this! It helped me fix my problem in far less time than if I were to do this on my own. I really appreciate it that you came back to edit the post after finding a way to make it work.

1 Like